Functions `parallelproj.functions`¶

class parallelproj.functions.C1AffineObjective(loss: C1Function, op: LinearOperator, s: Array | None = None)[source]¶

Bases: C1Function

Composes a prediction-space C1Function with an affine forward model.

Turns \(g(\bar{y})\) into \(f(x) = g(A x + s)\) using the chain rule:

\[\nabla_x f(x) = A^H \nabla_{\bar{y}} g(A x + s).\]

Any scaling is carried exclusively by the beta attribute of loss. When \(s\) is None the pure linear model \(\bar{y} = A x\) is used, avoiding the addition entirely.

Parameters:

loss (C1Function) – A C1Function operating on the prediction space.
op (LinearOperator) – The linear part of the forward model \(A\).
s (Array or None, optional) – Additive contamination term \(s\) (e.g. scatter/randoms). None (default) selects the pure linear model \(\bar{y} = A x\).

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

Listmode MLEM, OSEM, and SVRG

class parallelproj.functions.C1Function(beta: float = 1.0)[source]¶

Bases: ABC

Abstract base class for continuously differentiable (C1) scalar functions with an optional scalar scale factor \(\beta\).

The public interface (__call__(), gradient(), call_and_gradient()) evaluates \(\beta f(x)\). Subclasses only implement the unscaled private methods _call() and _gradient(); \(\beta\) is applied automatically.

Parameters:: beta (float, optional) – Multiplicative scale factor \(\beta\) applied to the function value and all derivatives. Defaults to 1.0.

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float[source]¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array][source]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array[source]¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

Listmode MLEM, OSEM, and SVRG

class parallelproj.functions.C1FunctionWithConjProx(beta: float = 1.0)[source]¶

Bases: C1Function, FunctionWithConjProx

Abstract base class for C1 functions that also admit a closed-form proximal operator of their convex conjugate.

Combines C1Function (gradient, call) with FunctionWithConjProx (prox_convex_conj()). Subclasses must implement _call(), _gradient(), and _prox_convex_conj().

MRO: C1FunctionWithConjProx -> C1Function -> FunctionWithConjProx -> ABC

The beta property and __call__ are resolved from C1Function. The prox_convex_conj() public method is resolved from FunctionWithConjProx.

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

Parameters:: beta (float)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

prox(x: Array, sigma: float | Array) → Array¶

Proximal operator of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)}(x) = x - \sigma \, \text{prox}_{(\beta f)^* / \sigma}(x / \sigma)\]

Subclasses with a cheaper closed-form direct prox may override this.

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array¶

Proximal operator of the convex conjugate of \(\beta f\).

Uses the identity

\[\text{prox}_{\sigma (\beta f)^*}(y) = \beta \, \text{prox}_{(\sigma / \beta)\, f^*}(y / \beta)\]

to reduce to the unscaled _prox_convex_conj().

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.C2AffineObjective(loss: C2Function, op: LinearOperator, s: Array | None = None)[source]¶

Bases: C2Function, C1AffineObjective

Composes a prediction-space C2Function with an affine forward model.

Extends C1AffineObjective with second-order information. For \(f(x) = g(A x + s)\) the Hessian-vector product is:

\[H_f(x)\, v = A^H \bigl(\operatorname{diag}(H_g(\bar{y})) \odot (A v)\bigr)\]

where \(\bar{y} = A x + s\) (or \(\bar{y} = A x\) when \(s\) is None) and \(\odot\) denotes elementwise multiplication. Any scaling is carried exclusively by the beta attribute of loss.

Parameters:

loss (C2Function) – A C2Function operating on the prediction space.
op (LinearOperator) – The linear part of the forward model \(A\).
s (Array or None, optional) – Additive contamination term \(s\) (e.g. scatter/randoms). None (default) selects the pure linear model \(\bar{y} = A x\).

Examples

>>> import numpy as np
>>> from parallelproj.operators import MatrixOperator
>>> from parallelproj.functions import NegPoissonLogL, C2AffineObjective
>>> # 4-element image space, 6-element sinogram space
>>> A = np.random.rand(6, 4)
>>> op = MatrixOperator(A)
>>> s = 0.1 * np.ones(6)          # scatter/randoms contamination
>>> data = np.ones(6)             # measured counts
>>> x = np.ones(4)                # image estimate
>>> v = np.ones(4)                # direction vector
>>>
>>> loss = NegPoissonLogL(data)
>>> loss.beta = 0.5
>>> aff_obj = C2AffineObjective(loss, op, s)
>>>
>>> fx = aff_obj(x)                           # scalar value, scaled by beta=0.5
>>> grad = aff_obj.gradient(x)                          # shape (4,), scaled by beta=0.5
>>> hv   = aff_obj.hessian_diag_vec_prod(x, v)          # shape (4,), scaled by beta=0.5

Pure linear forward model with virtual bins (zero rows in \(A\)) handled via NegPoissonLogLSafe and a user-supplied mask:

>>> from parallelproj.functions import NegPoissonLogLSafe
>>> A2 = np.zeros((6, 4))
>>> A2[:4, :] = np.random.rand(4, 4)   # last 2 rows are virtual (all zero)
>>> op2   = MatrixOperator(A2)
>>> data2 = np.array([2., 1., 3., 0., 0., 0.])  # virtual bins measure 0
>>> mask  = np.array([True, True, True, True, False, False])
>>>
>>> loss2 = NegPoissonLogLSafe(data2, mask)
>>> aff_obj2 = C2AffineObjective(loss2, op2)            # no contamination s
>>>
>>> fx2   = aff_obj2(x)                                 # scalar function value, no nan
>>> grad2 = aff_obj2.gradient(x)                        # shape (4,), no nan
>>> hv2   = aff_obj2.hessian_diag_vec_prod(x, v)        # shape (4,), no nan

A regularised objective combining a data fidelity term and a roughness penalty via SumC2Function:

\[h(x) = \underbrace{f_{\text{PL}}(Ax + s)}_{\text{data fidelity}} + \underbrace{\beta_{\text{reg}} \cdot \tfrac{1}{2}\|Dx\|_2^2}_{\text{roughness penalty}}\]

where \(D\) is a finite forward difference operator.

>>> from parallelproj.operators import FiniteForwardDifference
>>> from parallelproj.functions import HalfSquaredL2Deviation
>>> beta_reg = 0.1
>>> D   = FiniteForwardDifference(x.shape)            # finite differences in image space
>>> reg = C2AffineObjective(HalfSquaredL2Deviation(beta=beta_reg), D)
>>> data_fidelity = C2AffineObjective(NegPoissonLogL(data), op, s)
>>>
>>> obj_fun = data_fidelity + reg                     # SumC2Function
>>> obj_val = obj_fun(x)                              # scalar function value
>>> grad = obj_fun.gradient(x)                        # shape (4,)
>>> hv   = obj_fun.hessian_diag_vec_prod(x, v)        # shape (4,)

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

Listmode MLEM, OSEM, and SVRG

hessian_diag_vec_prod(x: Array, v: Array) → Array¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

Convergence comparison: SGD vs SVRG with quadratic regularization

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

class parallelproj.functions.C2Function(beta: float = 1.0)[source]¶

Bases: C1Function

Abstract base class for twice continuously differentiable (C2) scalar functions with an optional scalar scale factor \(\beta\).

Extends C1Function with curvature information. The public hessian_diag_vec_prod() returns \(\beta \operatorname{diag}(H_f(x))\, v\). Subclasses implement the unscaled _hessian_diag_vec_prod().

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

Parameters:: beta (float)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

Listmode MLEM, OSEM, and SVRG

hessian_diag_vec_prod(x: Array, v: Array) → Array[source]¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

Convergence comparison: SGD vs SVRG with quadratic regularization

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

class parallelproj.functions.C2FunctionWithConjProx(beta: float = 1.0)[source]¶

Bases: C2Function, C1FunctionWithConjProx

Abstract base class for C2 functions that also admit a closed-form proximal operator of their convex conjugate.

Combines C2Function (Hessian diagonal) with C1FunctionWithConjProx (gradient, call, dual prox). Subclasses must implement _call(), _gradient(), _hessian_diag_vec_prod(), and _prox_convex_conj().

MRO: C2FunctionWithConjProx -> C2Function -> C1FunctionWithConjProx -> C1Function -> FunctionWithConjProx -> ABC

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

Parameters:: beta (float)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

hessian_diag_vec_prod(x: Array, v: Array) → Array¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

prox(x: Array, sigma: float | Array) → Array¶

Proximal operator of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)}(x) = x - \sigma \, \text{prox}_{(\beta f)^* / \sigma}(x / \sigma)\]

Subclasses with a cheaper closed-form direct prox may override this.

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array¶

Proximal operator of the convex conjugate of \(\beta f\).

Uses the identity

\[\text{prox}_{\sigma (\beta f)^*}(y) = \beta \, \text{prox}_{(\sigma / \beta)\, f^*}(y / \beta)\]

to reduce to the unscaled _prox_convex_conj().

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.FunctionWithConjProx(beta: float = 1.0)[source]¶

Bases: ABC

Abstract base class for functions with a closed-form proximal operator of their convex conjugate.

This class is a standalone root – it does not require the function to be differentiable. Non-smooth functions (e.g. total variation, indicator functions) that admit a closed-form \(\text{prox}_{\sigma f^*}\) should inherit directly from this class.

Differentiable functions that additionally have a closed-form dual prox should use C1FunctionWithConjProx or C2FunctionWithConjProx instead.

The public prox_convex_conj() handles the \(\beta\) scaling automatically. Subclasses implement only the unscaled private methods _call() and _prox_convex_conj().

A default prox() is provided via Moreau’s identity for convenience; subclasses may override it with a more efficient closed-form if available.

Parameters:: beta (float, optional) – Multiplicative scale factor \(\beta\). Defaults to 1.0.

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float[source]¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value \(\beta f(x)\).
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

prox(x: Array, sigma: float | Array) → Array[source]¶

Proximal operator of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)}(x) = x - \sigma \, \text{prox}_{(\beta f)^* / \sigma}(x / \sigma)\]

Subclasses with a cheaper closed-form direct prox may override this.

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array[source]¶

Proximal operator of the convex conjugate of \(\beta f\).

Uses the identity

\[\text{prox}_{\sigma (\beta f)^*}(y) = \beta \, \text{prox}_{(\sigma / \beta)\, f^*}(y / \beta)\]

to reduce to the unscaled _prox_convex_conj().

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.FunctionWithProx(beta: float = 1.0)[source]¶

Bases: ABC

Abstract base class for functions with a closed-form proximal operator.

This class is a standalone root – it does not require the function to be differentiable. Functions (e.g. indicator functions, L1 norm) that admit a closed-form \(\text{prox}_{\sigma f}\) should inherit directly from this class.

The public prox() handles the \(\beta\) scaling automatically. Subclasses implement only the unscaled private methods _call() and _prox().

A default prox_convex_conj() is provided via Moreau’s identity for convenience; subclasses may override it with a more efficient closed-form if available.

Parameters:: beta (float, optional) – Multiplicative scale factor \(\beta\). Defaults to 1.0.

Examples

PDHG and SPDHG for PET reconstruction with a directional TV prior

__call__(x: Array) → float[source]¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value \(\beta f(x)\).
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

prox(x: Array, sigma: float | Array) → Array[source]¶

Proximal operator of \(\beta f\).

\[\text{prox}_{\sigma (\beta f)}(x) = \text{prox}_{(\sigma \beta) f}(x)\]

so the effective step size passed to the unscaled _prox() is \(\sigma \beta\).

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array[source]¶

Proximal operator of the convex conjugate of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)^*}(y) = y - \sigma \, \text{prox}_{(\beta f) / \sigma}(y / \sigma)\]

Subclasses with a cheaper closed-form dual prox may override this.

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.HalfSquaredL2Deviation(data: Array | None = None, weights: Array | None = None, beta: float = 1.0)[source]¶

Bases: C2FunctionWithConjProx

Half squared L2 deviation from reference data, with optional weights.

Implements

\[f(x) = \frac{1}{2} \sum_i w_i (x_i - d_i)^2\]

where \(w_i = 1\) when no weights are supplied (reducing to the standard \(\tfrac{1}{2}\|x - d\|_2^2\)).

Gradient:

\[\nabla f(x) = w \odot (x - d)\]

Diagonal Hessian-vector product:

\[\operatorname{diag}(H_f(x)) \odot v = w \odot v\]

The \(\tfrac{1}{2}\) prefactor is chosen so that the gradient contains no factor of 2, keeping expressions clean when \(\beta = 1\).

Parameters:

data (Array or None, optional) – Reference array \(d\). None (default) is equivalent to \(d = 0\) but avoids the subtraction entirely.
weights (Array or None, optional) – Non-negative weight array \(w\) of the same shape as data (or x when data is None). None (default) is equivalent to unit weights but avoids the multiplication entirely.
beta (float, optional) – Multiplicative scale factor \(\beta\). Defaults to 1.0.

Examples

Convergence comparison: SGD vs SVRG with quadratic regularization

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

hessian_diag_vec_prod(x: Array, v: Array) → Array¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

prox(x: Array, sigma: float | Array) → Array¶

Proximal operator of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)}(x) = x - \sigma \, \text{prox}_{(\beta f)^* / \sigma}(x / \sigma)\]

Subclasses with a cheaper closed-form direct prox may override this.

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array¶

Proximal operator of the convex conjugate of \(\beta f\).

Uses the identity

\[\text{prox}_{\sigma (\beta f)^*}(y) = \beta \, \text{prox}_{(\sigma / \beta)\, f^*}(y / \beta)\]

to reduce to the unscaled _prox_convex_conj().

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.MixedL21Norm(beta: float = 1.0)[source]¶

Bases: FunctionWithConjProx

Mixed L2-L1 norm (isotropic TV semi-norm on a gradient field).

For an array \(g\) whose first axis enumerates the gradient directions (as produced by FiniteForwardDifference), the norm is

\[f(g) = \sum_{\mathbf{i}} \|g_{:, \mathbf{i}}\|_2\]

where the sum runs over all spatial multi-indices \(\mathbf{i}\) and the L2 norm is taken along axis 0.

The convex conjugate is the indicator of the mixed \(L_{\infty,2}\) unit ball:

\[\begin{split}f^*(p) = \begin{cases} 0 & \|p_{:, \mathbf{i}}\|_2 \leq 1 \; \forall \mathbf{i} \\ +\infty & \text{otherwise} \end{cases}\end{split}\]

so its proximal operator is a pointwise projection onto the L2 unit ball along axis 0 (independent of \(\sigma\)):

\[\left(\text{prox}_{\sigma f^*}(y)\right)_{:, \mathbf{i}} = \frac{y_{:, \mathbf{i}}}{\max\!\left(\|y_{:, \mathbf{i}}\|_2,\, 1\right)}\]

The direct prox() (block soft-thresholding) is available via Moreau’s identity.

Parameters:: beta (float, optional) – Multiplicative scale factor \(\beta\) (regularization weight). Defaults to 1.0.

Examples

PDHG and SPDHG for PET reconstruction with a directional TV prior

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value \(\beta f(x)\).
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

prox(x: Array, sigma: float | Array) → Array¶

Proximal operator of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)}(x) = x - \sigma \, \text{prox}_{(\beta f)^* / \sigma}(x / \sigma)\]

Subclasses with a cheaper closed-form direct prox may override this.

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array¶

Proximal operator of the convex conjugate of \(\beta f\).

Uses the identity

\[\text{prox}_{\sigma (\beta f)^*}(y) = \beta \, \text{prox}_{(\sigma / \beta)\, f^*}(y / \beta)\]

to reduce to the unscaled _prox_convex_conj().

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.NegPoissonLogL(data: Array)[source]¶

Bases: C2FunctionWithConjProx

Negative Poisson log-likelihood as a function of expected counts.

Implements

\[f(\bar{y}) = \sum_i \bar{y}_i - y_i \log \bar{y}_i\]

and its gradient w.r.t. the expected counts \(\bar{y}\):

\[\nabla_{\bar{y}} f = 1 - \frac{y}{\bar{y}}.\]

This class operates directly on the predicted counts \(\bar{y}\), independently of how they were computed. Use C2AffineObjective to compose with a forward model \(\bar{y}(x) = A x + s\).

Parameters:: data (Array) – Measured data \(y\).

Examples

Convergence comparison: MLEM vs OSEM vs SVRG

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

hessian_diag_vec_prod(x: Array, v: Array) → Array¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

prox(x: Array, sigma: float | Array) → Array¶

Proximal operator of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)}(x) = x - \sigma \, \text{prox}_{(\beta f)^* / \sigma}(x / \sigma)\]

Subclasses with a cheaper closed-form direct prox may override this.

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array¶

Proximal operator of the convex conjugate of \(\beta f\).

Uses the identity

\[\text{prox}_{\sigma (\beta f)^*}(y) = \beta \, \text{prox}_{(\sigma / \beta)\, f^*}(y / \beta)\]

to reduce to the unscaled _prox_convex_conj().

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.NegPoissonLogLListmode(lm_op: LinearOperator, sensitivity_image: Array, contamination_list: Array, contamination_sinogram_sum: float)[source]¶

Bases: C2Function

Negative Poisson log-likelihood for listmode (event-by-event) data.

Implements the listmode equivalent of NegPoissonLogL with an affine forward model \(\bar{y}_e = (A_{\text{LM}}\, x)_e + s_e\). The function value is

\[f(x) = \langle \text{sens},\, x \rangle + c_{\text{sino}} - \sum_{e=1}^{N_{\text{ev}}} \log\bigl((A_{\text{LM}}\,x)_e + s_e\bigr)\]

where

\(\text{sens} = A_{\text{full}}^T \mathbf{1}\) is the sensitivity image (backprojection of all-ones from the full sinogram),
\(c_{\text{sino}} = \sum_i s_i^{\text{sino}}\) is the scalar sum of the contamination over all sinogram bins, and
the sum runs over all \(N_{\text{ev}}\) detected events.

This is mathematically equivalent to the sinogram NegPoissonLogL because \(\sum_e \log \bar{y}_{j_e} = \sum_i y_i \log \bar{y}_i\) (each event from bin \(i\) contributes \(\log \bar{y}_i\) and there are \(y_i\) such events).

The gradient with respect to the image \(x\) is

\[\nabla_x f(x) = \text{sens} - A_{\text{LM}}^T \!\left(\frac{1}{A_{\text{LM}}\,x + s_{\text{LM}}}\right)\]

and the diagonal Hessian-vector product is

\[\operatorname{diag}(H_f(x)) \odot v = A_{\text{LM}}^T \!\left( \frac{A_{\text{LM}}\, v}{(A_{\text{LM}}\,x + s_{\text{LM}})^2} \right).\]

Note

Unlike NegPoissonLogL, this class operates directly in image space (it takes \(x\) as input, not the predicted counts \(\bar{y}\)). It therefore cannot be composed with C2AffineObjective; the forward model is built in internally.

Parameters:

lm_op (LinearOperator) – Listmode forward projector \(A_{\text{LM}}\) mapping an image (shape: n_voxels) to per-event predicted counts (shape: n_events).
sensitivity_image (Array) – Sensitivity image \(A_{\text{full}}^T \mathbf{1}\) (same shape as the image \(x\)).
contamination_list (Array) – Per-event additive contamination \(s_{\text{LM}}\) (same shape as the listmode output, i.e. n_events).
contamination_sinogram_sum (float) – Scalar sum of the contamination over all sinogram bins, \(c_{\text{sino}} = \sum_i s_i^{\text{sino}}\).

Examples

Listmode MLEM, OSEM, and SVRG

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

Listmode MLEM, OSEM, and SVRG

hessian_diag_vec_prod(x: Array, v: Array) → Array¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

class parallelproj.functions.NegPoissonLogLSafe(data: Array, mask: Array, beta: float = 1.0)[source]¶

Bases: C2FunctionWithConjProx

Negative Poisson log-likelihood, safe for bins with zero expectation.

Identical to NegPoissonLogL but correctly handles virtual bins where the predicted counts \(\bar{y}_i = 0\). Naively evaluating \(0 \cdot \log 0\), \(0 / 0\), or \(0 / 0^2\) yields nan; this class avoids those cases via a user-supplied boolean mask \(m_i\) (True = active bin, False = virtual bin).

Note

The mask cannot be derived from the measured data \(y\) alone. With Poisson noise, any active bin can observe \(y_i = 0\) while still having \(\bar{y}_i > 0\), which is well-defined and requires no special treatment. The mask must reflect the structure of the forward model: a bin is virtual when no image voxel contributes to it (i.e. the corresponding row of \(A\) is zero) and the contamination \(s_i = 0\). This information is geometry-dependent and must be provided by the user. If \(s_i > 0\) for all bins, no virtual bins exist and NegPoissonLogL can be used directly.

Implements

\[f(\bar{y}) = \sum_i \bar{y}_i - y_i \log \bar{y}_i\]

with gradient

\[\nabla_{\bar{y}} f = 1 - \frac{y}{\bar{y}}\]

and diagonal Hessian-vector product

\[\operatorname{diag}(H_f(\bar{y})) \odot v = \frac{y}{\bar{y}^2} \odot v.\]

For virtual bins (\(m_i = \text{False}\), \(\bar{y}_i = 0\)) the mathematically correct limiting values are substituted:

Expression	Limit
\(\bar{y}_i - y_i \log \bar{y}_i\)	\(0\)
\(1 - y_i / \bar{y}_i\)	\(1\)
\(y_i / \bar{y}_i^2\)	\(0\)

Parameters:

data (Array) – Measured data \(y\).
mask (Array) – Boolean array of the same shape as data. True for active bins (\(\bar{y}_i > 0\) is guaranteed), False for virtual bins (\(\bar{y}_i = 0\) by construction).
beta (float, optional) – Multiplicative scale factor \(\beta\). Defaults to 1.0.

Examples

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

hessian_diag_vec_prod(x: Array, v: Array) → Array¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

prox(x: Array, sigma: float | Array) → Array¶

Proximal operator of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)}(x) = x - \sigma \, \text{prox}_{(\beta f)^* / \sigma}(x / \sigma)\]

Subclasses with a cheaper closed-form direct prox may override this.

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array¶

Proximal operator of the convex conjugate of \(\beta f\).

Uses the identity

\[\text{prox}_{\sigma (\beta f)^*}(y) = \beta \, \text{prox}_{(\sigma / \beta)\, f^*}(y / \beta)\]

to reduce to the unscaled _prox_convex_conj().

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.NonNegativeIndicator(beta: float = 1.0)[source]¶

Bases: FunctionWithProx

Indicator function of the non-negative orthant.

\[\begin{split}f(x) = \begin{cases} 0 & x \geq 0 \\ +\infty & \text{otherwise} \end{cases}\end{split}\]

The proximal operator is the projection onto the non-negative orthant:

\[\text{prox}_{\sigma f}(x) = \max(x, 0)\]

independent of \(\sigma\) and \(\beta\).

The dual prox (via Moreau’s identity) is the projection onto the non-positive orthant:

\[\text{prox}_{\sigma f^*}(y) = \min(y, 0)\]

Parameters:: beta (float, optional) – Multiplicative scale factor \(\beta\). Defaults to 1.0. For an indicator function the value is 0 or \(+\infty\) regardless of \(\beta > 0\), but \(\beta\) affects the effective step size passed to prox().

Examples

PDHG and SPDHG for PET reconstruction with a directional TV prior

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value \(\beta f(x)\).
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

prox(x: Array, sigma: float | Array) → Array¶

Proximal operator of \(\beta f\).

\[\text{prox}_{\sigma (\beta f)}(x) = \text{prox}_{(\sigma \beta) f}(x)\]

so the effective step size passed to the unscaled _prox() is \(\sigma \beta\).

Parameters:

x (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)}(x)\).

Return type:

Array

Examples

prox_convex_conj(y: Array, sigma: float | Array) → Array¶

Proximal operator of the convex conjugate of \(\beta f\) via Moreau’s identity.

\[\text{prox}_{\sigma (\beta f)^*}(y) = y - \sigma \, \text{prox}_{(\beta f) / \sigma}(y / \sigma)\]

Subclasses with a cheaper closed-form dual prox may override this.

Parameters:

y (Array) – Input array.
sigma (float or Array) – Step-size parameter \(\sigma > 0\).

Returns:

\(\text{prox}_{\sigma (\beta f)^*}(y)\).

Return type:

Array

Examples

class parallelproj.functions.SumC1Function(functions: Sequence[C1Function])[source]¶

Bases: C1Function

Sum of an arbitrary number of C1Function objects.

Represents

\[h(x) = \sum_{k} f_k(x)\]

where each \(f_k\) may itself carry its own \(\beta_k\). The gradients add accordingly:

\[\nabla h(x) = \sum_{k} \nabla f_k(x).\]

Instances are most conveniently created via the + operator on any two C1Function objects, but can also be constructed directly to sum more than two terms at once.

Parameters:: functions (Sequence[C1Function]) – The functions \(f_k\) to sum. Must contain at least one element.

Examples

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

class parallelproj.functions.SumC2Function(functions: Sequence[C2Function])[source]¶

Bases: C2Function, SumC1Function

Sum of an arbitrary number of C2Function objects.

Extends SumC1Function with second-order information. The diagonal Hessian-vector product is:

\[\operatorname{diag}\!\left(H_h(x)\right) \odot v = \sum_{k} \operatorname{diag}\!\left(H_{f_k}(x)\right) \odot v.\]

Parameters:: functions (Sequence[C2Function]) – The functions \(f_k\) to sum. Must contain at least one element.

Examples

Convergence comparison: SGD vs SVRG with quadratic regularization

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

__call__(x: Array) → float¶

Evaluate \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the function.
Returns:: Scaled scalar function value.
Return type:: float

Examples

property beta: float¶: Multiplicative scale factor \(\beta\).

call_and_gradient(x: Array) → tuple[float, Array]¶

Evaluate \(\beta f(x)\) and its gradient simultaneously.

Parameters:: x (Array) – Point at which to evaluate.
Returns:: Scaled function value and gradient.
Return type:: tuple[float, Array]

Examples

gradient(x: Array) → Array¶

Gradient of \(\beta f(x)\).

Parameters:: x (Array) – Point at which to evaluate the gradient.
Returns:: Array of the same shape as x containing \(\beta \nabla f(x)\).
Return type:: Array

Examples

hessian_diag_vec_prod(x: Array, v: Array) → Array¶

Scaled diagonal Hessian-vector product: \(\beta \operatorname{diag}(H_f(x))\, v\).

Parameters:

x (Array) – Point at which to evaluate the Hessian diagonal.
v (Array) – Vector to multiply with the Hessian diagonal, same shape as x.

Returns:

Array of the same shape as x containing \(\beta \operatorname{diag}(H_f(x)) \odot v\) (elementwise).

Return type:

Array

Examples

Functions parallelproj.functions¶

Functions `parallelproj.functions`¶