spyrit.core.recon.Tikhonov

class spyrit.core.recon.Tikhonov(meas_op, sigma: tensor, approx=False)[source]

Bases: Module

Tikhonov regularization (aka as ridge regression).

It estimates the signal \(x\in\mathbb{R}^{N}\) the from linear measurements \(y = Ax\in\mathbb{R}^{M}\) corrupted by noise by solving

\[\| y - Ax \|^2_{\Gamma{-1}} + \|x\|^2_{\Sigma^{-1}},\]

where \(\Gamma\) is covariance of the noise, and \(\Sigma\) is the signal covariance. In the case \(M\le N\), the solution can be computed as

\[\hat{x} = \Sigma A^\top (A \Sigma A^\top + \Gamma)^{-1} y,\]

where we assume that both covariance matrices are positive definite. The class is constructed from \(A\) and \(\Sigma\), while \(\Gamma\) is passed as an argument to forward(). Passing \(\Gamma\) to forward() is useful in the presence of signal- dependent noise.

Note

  • \(x\) can be a 1d signal or a vectorized image/volume. This can be specified by setting the meas_shape attribute of the measurement operator.

  • The above formulation assumes that the signal \(x\) has zero mean.

Args:

  • meas_op : Measurement operator (see meas).

Its measurement operator has shape \((M, N)\), with \(M\) the number of measurements and \(N\) the number of pixels in the image.

  • sigma : Signal covariance prior, of shape \((N, N)\).

  • diagonal_approximation : A boolean indicating whether to set

the non-diagonal elements of \(A \Sigma A^T\) to zero. Default is False. If True, this speeds up the computation of the inverse \((A \Sigma A^T + \Sigma_\alpha)^{-1}\).

Attributes:

  • meas_op : Measurement operator initialized as meas_op.

  • diagonal_approximation : Indicates if the diagonal approximation

is used.

  • img_shape : Shape of the image, initialized as meas_op.img_shape.

  • sigma_meas : Measurement covariance prior initialized as

\(A \Sigma A^T\). If diagonal_approximation is True, the non-diagonal elements are set to zero.

  • sigma_A_T : Covariance of the missing measurements initialized

as \(\Sigma A^T\).

  • noise_scale : Hidden parameter to use to scale the noise

regularization. It is used in the computation of the inverse: \((A \Sigma A^T + noisescale \times \Sigma_\alpha)^{-1}\). Default is 1.

Example:
>>> B, H, M, N = 85, 17, 32, 64
>>> sigma = torch.rand(N, N)
>>> gamma = torch.rand(M, M)
>>> A = torch.rand([M,N])
>>> meas  = Linear(A, meas_shape=(1,N))
>>> recon = Tikhonov(meas, sigma)
>>> y = torch.rand(B,H,M)
>>> x = recon(y, gamma)
>>> print(y.shape)
>>> print(x.shape)
torch.Size([85, 17, 32])
torch.Size([85, 17, 64])

Methods

divide(y, gamma)

Computes the division \(y \cdot (\sigma_lpha imes noisescale + (A \Sigma A^T))^{-1}\).

forward(y, gamma)

Reconstructs the signal from measurements and noise covariance.