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.

This class plays the same role as HadamSplit2d, but replaces the 2D Hadamard transform used there by a 2D S-matrix transform. Considering the acquisition of \(M\) square DMD patterns of size \(h\), it computes

\[y =\mathcal{N}\left(\mathcal{S}\left(S_{1d} X S_{1d}^T\right)\right),\]

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, \(S_{1d} \in \{0,1\}^{h\times h}\) is the (1D) Walsh-ordered S-matrix (see spyrit.misc.walsh_hadamard.walsh_S_matrix()), and \(X \in \mathbb{R}^{h\times h}\) is the (2D) image.

Unlike the Hadamard matrix used by HadamSplit2d, the S-matrix \(S_{1d}\) only takes values in \(\{0, 1\}\). It is therefore already directly realizable as a set of DMD patterns, and there is no need to split it into positive and negative components as LinearSplit-based classes do. For this reason, HadamSmatrix2d inherits directly from Linear (not LinearSplit), and any preprocessing needed downstream reduces to a scaling of the raw measurements (no “unsplitting” step is required).

The S-matrix is built from a Hadamard matrix of order \(h+1\), so \(h+1\) (not \(h\)) must be a power of two.

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 \(H = S_{1d}\otimes S_{1d}\), 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 of the 1D S-matrix \(S_{1d}\) – i.e. of the acquisition matrix itself, 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.

In short: scramble acts on the acquisition matrix \(S_{1d}\) (columns), order acts on the sequence of measurements in \(y\) (rows of \(H\), built from \(S_{1d}\) after scrambling has already been applied, if any). The two options are independent and can be combined.

If scramble is True, the columns of \(S_{1d}\) are randomly permuted (with a fixed seed for reproducibility). This is useful e.g. to decorrelate the acquisition order from the natural Walsh ordering. Because permuting columns destroys the symmetry of \(S_{1d}\), its inverse is no longer given directly by iwalsh_S_matrix(); instead, it is obtained exactly (with no extra matrix inversion) using spyrit.core.torch.reindex() (see T1d and the note below), and the adjoint operator is computed using \(S_{1d}^T\) explicitly rather than relying on \(S_{1d}\) being symmetric.

Note

\(S_{1d}\) is not orthogonal: unlike the Hadamard matrix, \(S_{1d}^{-1} \neq S_{1d}^T / h\). The adjoint (transpose) of the measurement operator and its pseudo-inverse are therefore genuinely different operators here, computed respectively from \(S_{1d}^T\) and from \(T_{1d}\) (see T1d).

Note

Let column_perm be the random permutation used to scramble the columns of \(S_{1d}\) : \(S_{1d}^{\text{scrambled}} = S_{1d} P\), where \(P\) is the corresponding (orthogonal) permutation matrix. Since \(P\) is orthogonal, \((S_{1d}P)^{-1} = P^{-1}S_{1d}^{-1} = P^TS_{1d}^{-1}\). Building \(S_{1d}P\) gathers the columns of \(S_{1d}\) directly by column_perm; by the structure of permutation matrices, gathering the rows of \(S_{1d}^{-1}\) with that very same (non-inverted) column_perm array already computes \(P^TS_{1d}^{-1}\) – no separate inverse-permutation array is needed at the indexing level, since transposing \(P\) is already accounted for by switching which axis (rows vs. columns) is gathered, not by inverting the indices themselves. Concretely, this means the same column_perm array is passed to reindex() for both matrices, but with different arguments: reindex(S1d, column_perm, axis=”cols”, inverse_permutation=True) for \(S_{1d}\), and reindex(T1d, column_perm, axis=”rows”, inverse_permutation=False) for \(T_{1d}\) – using column_perm’s argsort (the literal “inverse indices”) in either of these two calls would silently compute the wrong matrix.

Args:

h (int): Image size \(h\). \(h+1\) 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 \(H\) (i.e. the sequence of measurements in \(y\)). Not to be confused with scramble, which permutes the columns of \(S_{1d}\) (the acquisition matrix itself).

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 \(S_{1d}\) (the acquisition matrix) are randomly permuted. Defaults to False. Not to be confused with order, which permutes the rows of the system matrix \(H\) (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:

S1d (torch.tensor): 1D S-matrix of shape \((h,h)\), values in \(\{0,1\}\). Its columns are permuted if scramble is True (see column_perm).

T1d (torch.tensor): Inverse of S1d, of shape \((h,h)\). When scramble is True, obtained exactly by permuting the rows of the unscrambled inverse with column_perm (see __init__()), since the closed-form inverse from spyrit.misc.walsh_hadamard.iwalsh_S_matrix() only applies directly to the unscrambled matrix.

column_perm (torch.tensor, optional): The random column permutation applied to \(S_{1d}\). Only set as an attribute when scramble is True.

H (torch.tensor): The 2D measurement matrix given by \(S_{1d}\otimes S_{1d}\), subsampled to self.M rows. Computed on the fly (not stored) to avoid materializing a \(h^2 \times h^2\) matrix.

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 \(S_{1d}\) 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

>>> h = 63  # h + 1 = 64 = 2**6
>>> meas_op = HadamSmatrix2d(h, 2000)
>>> print(meas_op.S1d.shape)
torch.Size([63, 63])
>>> print(meas_op.M)
2000
>>> x = torch.rand(4, h, h)
>>> y = meas_op(x)
>>> print(y.shape)
torch.Size([4, 2000])

Example 2: With full sampling (M = \(h^2\), the default), fast_pinv() exactly recovers the image, since no measurement is discarded and T1d is the exact matrix inverse of S1d.

>>> h = 15  # h + 1 = 16 = 2**4
>>> 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)
>>> 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 the columns of S1d, breaking its symmetry, but full-sampling recovery via fast_pinv() remains exact (see the note above on why the adjoint and pseudo-inverse are handled correctly regardless).

>>> meas_op_scrambled = HadamSmatrix2d(h, scramble=True, seed=42)
>>> print(torch.allclose(meas_op.S1d, meas_op.S1d.T))       # unscrambled: symmetric
True
>>> print(torch.allclose(meas_op_scrambled.S1d, meas_op_scrambled.S1d.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 S1d itself.

>>> order = torch.rand(h, h)
>>> meas_op_ordered = HadamSmatrix2d(h, M=100, order=order)
>>> print(torch.equal(meas_op_ordered.S1d, meas_op.S1d))
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 (only multiplications with the 1D S-matrix `self.S1d` are required).
`fast_pinv`(m[, vectorize])	Apply the pseudo-inverse of the measurement matrix.
`forward`(x)	Simulate noisy measurements \(y = \mathcal{N}(\mathcal{S}(S_{1d} X S_{1d}^T))\).
`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.