spyrit.core.inverse.TikhonovMeasurementPriorDiag.forward

TikhonovMeasurementPriorDiag.forward(x: tensor, x_0: tensor, var: tensor) → tensor

Computes the Tikhonov regularization with prior in the measurement domain.

This method, unlike the forward_no_prior() method, allows for a non-zero mean image prior \(x_0\). We approximate the solution as:

\[\begin{split}\hat{x} = x_0 + F^{-1} \begin{bmatrix} m_1 \\ m_2\end{bmatrix}\end{split}\]

with \(m_1 = D_1(D_1 + \Sigma_\alpha)^{-1} (m - GF x_0)\) and \(m_2 = \Sigma_1 \Sigma_{21}^{-1} m_1\), where \(\Sigma = \begin{bmatrix} \Sigma_1 & \Sigma_{21}^\top \\ \Sigma_{21} & \Sigma_2\end{bmatrix}\) and \(D_1 =\textrm{Diag}(\Sigma_1)\). Assuming the noise covariance \(\Sigma_\alpha\) is diagonal, the matrix inversion involved in the computation of \(m_1\) is straightforward.

This is an approximation to the exact solution

\[\begin{split}\hat{x} &= x_0 + F^{-1}\begin{bmatrix}\Sigma_1 \\ \Sigma_{21} \end{bmatrix} [\Sigma_1 + \Sigma_\alpha]^{-1} (m - GF x_0)\end{split}\]

See Lemma B.0.5 of the PhD dissertation of A. Lorente Mur (2021): https://theses.hal.science/tel-03670825v1/file/these.pdf

Args:

x (torch.tensor): A batch of measurement vectors \(m\) with shape \((*, M)\).

x_0 (torch.tensor): A batch of prior images \(x_0\) with shape \((*, \sqrt{N}, \sqrt{N})\).

var (torch.tensor): A batch of measurement noise variances \(\Sigma_\alpha\) with shape \((*, M)\).

meas_op (torch.tensor): A measurement operator that provides \(GF\) and \(F^{-1}\).

Output:

(torch.tensor): A batch of images with shape \((*, \sqrt{N}, \sqrt{N})\).

Example:

>>> from spyrit.core.meas import HadamSplit2d
>>> from spyrit.core.inverse import TikhonovMeasurementPriorDiag
>>> import torch

>>> acqu = HadamSplit2d(32, 400)
>>> sigma = torch.rand([32*32, 32*32])
>>> recon_op = TikhonovMeasurementPriorDiag(acqu, sigma)
>>> y = torch.rand([10, 3, 400])
>>> x0 = torch.rand([10, 3, 32, 32])
>>> var = torch.rand([10, 3, 400])
>>> x = recon_op(y, x0, var)
>>> print(x.shape)
torch.Size([10, 3, 32, 32])