spyrit.core.recon.TikhonovMeasurementPriorDiag.forward

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

Computes the Tikhonov regularization with prior in the measurement domain.

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: A batch of measurement vectors \(m\)

• x_0: A batch of prior images \(x_0\)

• var: A batch of measurement noise variances \(\Sigma_\alpha\)

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

Shape:

• x: \((*, M)\)

• x_0: \((*, N)\)

• var \((*, M)\)

• Output: \((*, N)\)

Example:

>>> B, H, M = 85, 32, 512
>>> sigma = torch.rand([H**2, H**2])
>>> recon_op = TikhonovMeasurementPriorDiag(sigma, M)
>>> Ord = torch.ones((H,H))
>> meas = HadamSplit(M, H, Ord)
>>> y = torch.rand([B,M], dtype=torch.float)
>>> x_0 = torch.zeros((B, H**2), dtype=torch.float)
>>> var = torch.zeros((B, M), dtype=torch.float)
>>> x = recon_op(y, x_0, var, meas)
torch.Size([85, 1024])