spyrit.core.inverse.Tikhonov

class spyrit.core.inverse.Tikhonov(meas_op: Linear, sigma: tensor, approx=False, reshape_output: bool = True)[source]

Bases: Module

Tikhonov regularization.

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).

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

approx : 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}\).

reshape_output : A boolean indicating whether to reshape the output to the shape of the image. Default is True.

Attributes:

meas_op : Measurement operator initialized as meas_op.

sigma : Signal covariance prior initialized as sigma.

approx : Indicates if the diagonal approximation is used.

reshape_output : Indicates if the output is reshaped.

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

sigma_meas : Measurement covariance prior initialized as \(A \Sigma A^T\). If approx is True, the non-diagonal elements are set to zero. It is pre-computed at initialization to speed up future computations.

sigma_A_T : Covariance of the missing measurements initialized as \(\Sigma A^T\). It is computed at initialization to speed up future computations.

Example:

>>> from spyrit.core.meas import Linear
>>> 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)
torch.Size([85, 17, 32])
>>> print(x.shape)
torch.Size([85, 17, 1, 64])

Methods

`divide`(y, gamma)	Computes the division \(y \cdot (\Sigma \alpha + (A \Sigma A^T))^{-1}\).
`forward`(y, gamma)	Reconstructs the signal from measurements and noise covariance.