spyrit.core.meas.HadamSmatrix2d

class spyrit.core.meas.HadamSmatrix2d(h: int, M: int = None, order: tensor = None, fast: bool = True, reshape_output: bool = False, scramble: bool = False, seed: int = 0, *, noise_model=Identity(), dtype: dtype = torch.float32, device: device = device(type='cpu'))[source]

Bases: Linear

Simulate 2D S-matrix acquisitions.

Considering the acquisition of \(M\) DMD patterns of size \(h\times h\), the class computes

\[y = \mathcal{N}\left(\mathcal{S}\left(Y\right)\right), \quad Y = \frac{\mathrm{sum}(X)\cdot J_h - H_h\,X\,H_h^T}{2},\]

where \(\mathcal{N} \colon\, \mathbb{R}^M \to \mathbb{R}^M\) represents a noise operator (e.g., Gaussian), \(\mathcal{S} \colon\, \mathbb{R}^{h\times h} \to \mathbb{R}^M\) is a subsampling operator, \(J_h\) is the \(h\times h\) all-ones matrix, and \(H_h \in \{-1,+1\}^{h\times h}\) is the same 1D Walsh-ordered Hadamard matrix used by HadamSplit2d (see HadamSplit2d.H1d).

Equivalently,

\[y = \mathcal{N}\left(\mathcal{S}(A\,\mathrm{vec}(X))\right), \quad A = \max(0,-(H_h\otimes H_h)),\]

where \(A \in \{0,1\}^{h^2 \times h^2}\) is the acquisition matrix – see A.

Note

By definition, the first row and column of \(H_h \in \{-1,+1\}^{h\times h}\) are all-ones. Therefore, \(A = \max(0, -(H_h\otimes H_h))\) has an all-zero first row and an all-zero first column, which correspond to the trivial measurement and the unrecoverable top-left pixel.

Warning

Because \(H_h\)’s first row and column are both all-ones (the usual Walsh-ordering convention), \(A\) has an all-zero row and an all-zero column. scramble, when used, only permutes columns \(1,\dots,h-1\) of \(H_h\), leaving column 0 fixed (see column_perm), so both of the following hold regardless of scramble:

\(A\)’s all-zero row means the very first measurement is always exactly 0 (at natural-order index 0, the position order keeps first by default), regardless of \(X\) – the “all DMD micromirrors off” pattern.
\(A\)’s all-zero column means the top-left pixel, \(X[0,0]\), has no effect whatsoever on any measurement, and is therefore fundamentally unrecoverable. By convention, input images are expected to have \(X[0,0]=0\) (see fast_measure()); fast_pinv() always reconstructs it as exactly 0.

Note

order and scramble both rely on permutations, but they act on different things and serve different purposes – do not confuse them:

order permutes the rows of the (subsampled) 2D system matrix \(A\), i.e. it reorders the \(M\) measurements in the output vector \(y\) (and selects which \(M\) of the \(h^2\) possible measurements are kept, if \(M<h^2\)). It does not change what any individual measurement pattern looks like, only the sequence in which the measurements appear (e.g. by decreasing variance/significance, via order).
scramble permutes the columns \(1,\dots,h-1\) of \(H_h\) – i.e. of the (1D) matrix used to build the acquisition matrix, before any row reordering/subsampling happens. This changes which pixels of the image each individual measurement pattern probes, not the order in which measurements are returned.

The two options are independent and can be combined.

If scramble is True, columns \(1,\dots,h-1\) of \(H_h\) are randomly permuted (with a fixed seed for reproducibility; column 0 is always left fixed, see the warning above). This is useful e.g. to decorrelate the acquisition order from the natural Walsh ordering.

Note

\(A\) is not orthogonal, so the adjoint (transpose) of the measurement operator and its pseudo-inverse are genuinely different operators here. Both, however, admit a closed form expressed only in terms of \(H_h\) (and its transpose), so both remain fast and separable across rows and columns – no \(h^2\times h^2\) matrix and no explicit matrix inversion is ever needed, even when scramble is True (see fast_adjoint() and fast_pinv() for the exact formulas).

Args:

h (int): Image size \(h\), which must be a power of 2.

M (int, optional): Number of measurements. Defaults to \(h^2\) (no subsampling).

order (torch.tensor, optional): Order matrix \(O\) that defines the measurements to keep. The first component of \(y\) will correspond to the index where order is the highest. Permutes the rows of the system matrix \(A\) (i.e. the sequence of measurements in \(y\)). Not to be confused with scramble, which permutes the columns of \(H_h\).

fast (bool, optional): Whether to use the fast, separable computation of the 2D S-transform. If False, it uses (memory-heavy) matrix-vector products with the full measurement matrix. Defaults to True.

reshape_output (bool, optional): Whether to reshape the output of the adjoint and pseudo-inverse methods to images. If False, outputs are vectors.

scramble (bool, optional): If True, the columns of \(H_h\) are randomly permuted. Defaults to False. Not to be confused with order, which permutes the rows of the system matrix \(A\) (i.e. the sequence of measurements in \(y\)), applied after scrambling.

seed (int, optional): Seed used to generate the random column permutation when scramble is True, ensuring reproducibility. Defaults to 0.

noise_model (see spyrit.core.noise): Noise model \(\mathcal{N}\). Defaults to torch.nn.Identity().

dtype (torch.dtype, optional): Data type of the measurement matrix. Defaults to torch.float32.

device (torch.device, optional): Device of the measurement matrix. Defaults to torch.device(“cpu”).

Attributes:

H1d (torch.tensor): 1D Walsh-ordered Hadamard matrix of shape \((h,h)\), values in \(\{-1,+1\}\) (the same matrix as HadamSplit2d.H1d). Columns \(1,\dots,h-1\) are permuted if scramble is True (see column_perm); column 0 is always left fixed.

column_perm (torch.tensor, optional): The random column permutation applied to \(H_h\); column_perm[0] == 0 always. Only set as an attribute when scramble is True.

h (int): Image size \(h\).

A (torch.tensor): The 2D measurement matrix given by \(\max(0, -(H_h\otimes H_h))\), subsampled to self.M rows. Computed on the fly (not stored) to avoid materializing a \(h^2 \times h^2\) matrix.

H (torch.tensor): Alias for A, kept for compatibility with code that expects a measurement operator to expose H (see A’s docstring).

M (int): Number of measurements \(M\).

N (int): Number of pixels in the image, equal to \(h^2\).

meas_shape (torch.Size): Shape of the measurement patterns. Equal to \((h, h)\).

meas_dims (torch.Size): Dimensions of the image the acquisition matrix applies to. Equal to (-2, -1).

order (torch.tensor): Order matrix \(O\). Only affects the row order (sequence) of the measurements in \(y\); unrelated to scramble, which affects the columns of \(H_h\) instead.

indices (torch.tensor): Indices used to reorder the measurement vector (derived from order). Used by the method reindex().

Example 1:

Basic construction and simulated (subsampled) measurements. The very first measurement is always (up to floating-point error) 0, since it corresponds to the trivial all-zero row of \(A\) (see the warning above).

>>> h = 64  # h must be a power of two
>>> meas_op = HadamSmatrix2d(h, 2000)
>>> print(meas_op.H1d.shape)
torch.Size([64, 64])
>>> print(meas_op.M)
2000
>>> x = torch.rand(4, h, h)
>>> y = meas_op(x)
>>> print(y.shape)
torch.Size([4, 2000])
>>> print(torch.allclose(y[:, 0], torch.zeros(4), atol=1e-2))
True

Example 2:

With full sampling (M = \(h^2\), the default) and the top-left pixel set to 0 (see the warning above), fast_pinv() recovers the image exactly.

>>> h = 16
>>> meas_op = HadamSmatrix2d(h)  # M defaults to h**2 (full sampling)
>>> print(meas_op.M == meas_op.N)
True
>>> x = torch.rand(2, h, h)
>>> x[:, 0, 0] = 0  # convention: top-left pixel carries no information
>>> y = meas_op.measure(x)
>>> x_hat = meas_op.fast_pinv(y, vectorize=False)
>>> print(torch.allclose(x, x_hat, atol=1e-4))
True

Example 3:

With subsampling (M < \(h^2\)), fast_pinv() only approximates the image (see the note in fast_pinv()), unlike the exact recovery obtained above with full sampling.

>>> meas_op_sub = HadamSmatrix2d(h, M=100)
>>> y_sub = meas_op_sub.measure(x)
>>> x_hat_sub = meas_op_sub.fast_pinv(y_sub, vectorize=False)
>>> print(torch.allclose(x, x_hat_sub, atol=1e-4))
False

Example 4:

scramble permutes columns \(1,\dots,h-1\) of H1d (never column 0), breaking its symmetry – but full-sampling recovery via fast_pinv() remains exact, including at the top-left pixel, exactly as in Example 2.

>>> meas_op_scrambled = HadamSmatrix2d(h, scramble=True, seed=42)
>>> print(torch.allclose(meas_op.H1d, meas_op.H1d.T))       # unscrambled: symmetric
True
>>> print(torch.allclose(meas_op_scrambled.H1d, meas_op_scrambled.H1d.T))  # scrambled: not symmetric
False
>>> y2 = meas_op_scrambled.measure(x)
>>> x_hat2 = meas_op_scrambled.fast_pinv(y2, vectorize=False)
>>> print(torch.allclose(x, x_hat2, atol=1e-4))
True

Example 5:

order and scramble act independently (see note above): changing order selects/reorders which measurements are kept, but does not affect H1d itself.

>>> order = torch.rand(h, h)
>>> meas_op_ordered = HadamSmatrix2d(h, M=100, order=order)
>>> print(torch.equal(meas_op_ordered.H1d, meas_op.H1d))
True

Methods

`adjoint`(m[, unvectorize])	Apply the adjoint (transpose) of the measurement matrix.
`fast_adjoint`(m[, unvectorize])	Apply the adjoint of the measurement matrix using the separability of the 2D S-transform.
`fast_measure`(x)	Simulate noiseless measurements using the separability of the 2D S-transform.
`fast_pinv`(m[, vectorize])	Apply the pseudo-inverse of the measurement matrix.
`forward`(x)	Simulate noisy measurements \(y = \mathcal{N}(\mathcal{S}(\mathrm{vec}^{-1}(A\,\mathrm{vec}(X))))\).
`measure`(x)	Simulate noiseless measurements.
`reindex`(x[, axis, inverse_permutation])	Sorts a tensor along a specified axis using `self.indices`.
`set_matrix_to_inverse`(matrix_name)
`unvectorize`(input)	Unflatten the measured dimensions.
`vectorize`(input)	Flatten the measured dimensions.