spyrit.core.meas.HadamSmatrix2d.fast_pinv

HadamSmatrix2d.fast_pinv(m: tensor, vectorize: bool = False) → tensor

Apply the pseudo-inverse of the measurement matrix.

Args:

m (torch.tensor): Measurement \(m\) of length self.M.

vectorize (bool): Whether to apply vectorize() after computation of the pseudo-inverse.

Returns:

torch.tensor: Vectorized (or image-shaped) signal \(x\) of length self.N.

Note

\(A\) is given by \(A=\begin{pmatrix}0&0\\0&S\end{pmatrix}\), where \(S\) is the classical (order-\((h^2-1)\)) S-matrix derived from the genuine order-\(N\) Hadamard matrix \(H_N=H_h\otimes H_h\) introduced below (the Kronecker product of two Hadamard matrices is itself Hadamard). Being block-diagonal, \(A\)’s pseudo-inverse is therefore also block-diagonal:

\[\begin{split}A^+ = \begin{pmatrix} 0 & 0 \\ 0 & S^{-1} \end{pmatrix}.\end{split}\]

How it works. Writing \(H_N = H_h\otimes H_h\) (order \(N=h^2\), the full, non-separable Hadamard-like matrix such that \(A = (J_N - H_N)/2\)), \(A^+\) can be computed without ever forming \(S^{-1}\) explicitly:

\[A^+ y = -\frac{2}{N}\, H_N^T y,\]

applied to a measurement vector \(y\) (padded with a leading 0). This is exact at every entry except the leading one (always 0 for \(y\) in the range of \(A\), but not automatically 0 in this formula – see below).

How it is implemented. Since \(H_N = H_h\otimes H_h\), applying \(H_N^T\) to a vector is the same as applying \(H_h^T\) to both axes of its \(h\times h\) reshaping: for a (zero-padded) measurement image \(Y\),

\[X = -\frac{2}{N}\, H_h^T\, Y\, H_h,\]

which only requires multiplications with the 1D matrix \(H_h\) (applied to the rows and columns of \(Y\)) – no matrix inversion, and no \(h^2\times h^2\) matrix is ever materialized. If the number of measurements is smaller than the number of pixels, the measurement vector is zero-padded before inversion, exactly as done in HadamSplit2d.fast_pinv().

Note

Applying fast_pinv() to the (noiseless, full-sampling) measurement of an image returns that same image back exactly, at every pixel except the top-left one, \(X[0,0]\) (see the warning in the class docstring): that single pixel is fundamentally unrecoverable (it never affects any measurement), so this method arbitrarily sets it to 0 – the block-diagonal \(A^+\) above shows this convention is exact, not an approximation: \(X[0,0]\) is genuinely a free parameter of the formula above, which the top-left 0 block of \(A^+\) fixes at 0 (the minimum-norm, Moore-Penrose choice).

Note

When self.M equals self.N, this reconstruction is exact (up to \(X[0,0]\)). This holds whether or not scramble is used.