Exercises

$\mathcal{L}(\theta) = \frac{1}{2}\|\mathbf{y} - \mathbf{A}f_\theta(\mathbf{z})\|^2.$$Only the network weights$\theta$are optimised via gradient descent. The input$\mathbf{z}$is sampled once from$\mathcal{N}(\mathbf{0}, \mathbf{I})$and fixed. The reconstruction is$\hat{\mathbf{x}} = f_{\theta^*}(\mathbf{z})$.$\square$

$\mathbb{E}[\|f_\theta(\mathbf{y}_1) - \mathbf{y}_2\|^2] = \mathbb{E}[\|f_\theta(\mathbf{y}_1) - \mathbf{x} - \mathbf{w}_2\|^2]KATEXPLACEHOLDER0END= \mathbb{E}[\|f_\theta(\mathbf{y}_1) - \mathbf{x}\|^2] - 2\mathbb{E}[\langle f_\theta(\mathbf{y}_1) - \mathbf{x}, \mathbf{w}_2\rangle] + \mathbb{E}[\|\mathbf{w}_2\|^2].$$

Since $\mathbf{w}_2$ is independent of $\mathbf{y}_1$ (and hence of $f_\theta(\mathbf{y}_1)$ and $\mathbf{x}$) with $\mathbb{E}[\mathbf{w}_2] = \mathbf{0}$: $$\mathbb{E}[\langle f_\theta(\mathbf{y}_1) - \mathbf{x}, \mathbf{w}_2\rangle] = 0.$$

The N2N loss equals the supervised loss plus a constant $\mathbb{E}[\|\mathbf{w}_2\|^2]$. The minimiser over $\theta$ is the same. $\square$

$[f(\mathbf{y})]_i = \sum_j W_{ij} y_j$. Therefore $\frac{\partial [f]_i}{\partial y_i} = W_{ii}$ and $\operatorname{div}(f) = \sum_i W_{ii} = \operatorname{tr}(\mathbf{W})$.

$\text{SURE} = \frac{1}{N}\|(\mathbf{I} - \mathbf{W})\mathbf{y}\|^2 - \sigma^2 + \frac{2\sigma^2}{N}\operatorname{tr}(\mathbf{W}).$$Minimising over$\mathbf{W}$yields the Wiener filter$\mathbf{W}^* = \mathbf{C}_x(\mathbf{C}_x + \sigma^2\mathbf{I})^{-1}$.$\square$

$\mathcal{L}_{\text{DC}} = \mathbb{E}_\mathbf{y}\bigl[\|\mathbf{A}f_\theta(\mathbf{y}) - \mathbf{y}\|^2\bigr].$$This ensures the reconstruction is consistent with the measurements (constrains the range of$\mathbf{A}^H$).

$\mathcal{L}_{\text{EI}} = \frac{1}{K}\sum_{k=1}^K \mathbb{E}_\mathbf{y}\bigl[\|f_\theta(\mathbf{A}T_k f_\theta(\mathbf{y})) - T_k f_\theta(\mathbf{y})\|^2\bigr].$$This requires the network to produce consistent reconstructions under transformations (constrains the null space of$\mathbf{A}$).

$\mathcal{L} = \mathcal{L}_{\text{DC}} + \lambda\,\mathcal{L}_{\text{EI}}$ where $\lambda$ balances the two objectives. $\square$

1. Dynamic range: Natural images: 8-bit (0--255). RF reflectivity: 40+ dB dynamic range, often log-scaled. Mitigation: Normalise to log-magnitude before feeding to the foundation model.

2. Complex values: Natural images are real-valued (RGB). RF signals are complex (amplitude + phase). Mitigation: Use 2-channel (real/imaginary) representation, or separate magnitude/phase processing.

3. Texture statistics: Natural images have smooth textures and sharp edges. RF images have speckle (multiplicative noise), sidelobes, and grating lobes. Mitigation: Fine-tune on simulated RF data, or use LoRA adaptation with a small RF dataset. $\square$

With 1.5M parameters and 16K pixels, the network is $\sim 90\times$ over-parameterised. It has sufficient capacity to fit any function, including pure noise.

Empirically, DIP overfitting begins after $\sim 2{,}000$--$5{,}000$ iterations. The spectral bias delays overfitting: low-frequency signal ($\sim 50\%$ of energy) is fit in $\sim 500$ iterations; mid-frequency details in $\sim 2{,}000$; high-frequency noise in $\sim 5{,}000$+. The optimal stopping point depends on SNR.

A Deep Decoder with 50K parameters ($\sim 3\times$ the pixel count) cannot memorise noise due to under-parameterisation. It converges without overfitting but may miss fine details. The U-Net DIP achieves higher peak PSNR (at the optimal stopping point) but requires careful early stopping. $\square$

The divergence is $\operatorname{div}(f_\lambda) = |\{i : |y_i| > \lambda\}| \triangleq K(\lambda)$.

$$\text{SURE}(\lambda) = \frac{1}{N}\sum_i \min(y_i^2, \lambda^2) - \sigma^2 + \frac{2\sigma^2}{N}K(\lambda).$$

For Gaussian signal + Gaussian noise with known sparsity fraction $s = K_{\text{true}}/N$, the optimal threshold scales as $\lambda^* \approx \sigma\sqrt{2\log(N/K_{\text{true}})}$ (the universal threshold).

In practice, minimise SURE numerically: evaluate $\text{SURE}(\lambda)$ for a grid of $\lambda$ values and choose the minimiser. This requires no knowledge of $s$. $\square$

A circular shift by $\Delta$ pixels: $[T_\Delta\mathbf{x}]_n = x_{(n-\Delta) \bmod N}$. In Fourier domain: $[\mathbf{F}T_\Delta\mathbf{x}]_k = e^{-j2\pi k\Delta/N}\hat{x}_k$.

$\mathbf{A}T_\Delta\mathbf{x} = \mathbf{P}_\Omega\,\operatorname{diag}(e^{-j2\pi k\Delta/N})\,\hat{\mathbf{x}}$.

The phase rotation does not change which frequencies are measured, but it modulates them differently. Different shifts create a system of equations at each measured frequency.

The EI constraint $f_\theta(\mathbf{A}T_\Delta\hat{\mathbf{x}}) = T_\Delta\hat{\mathbf{x}}$ for multiple shifts $\{\Delta_1, \ldots, \Delta_K\}$ forces the network to produce consistent reconstructions under different phase modulations. This implicitly constrains the unmeasured frequencies through the network's learned inductive bias. $\square$

Define $h_\theta(\mathbf{y}) = \mathbf{A}f_\theta(\mathbf{y}) : \mathbb{R}^M \to \mathbb{R}^M$. Applying standard SURE: $$\text{GSURE} = \frac{1}{M}\|\mathbf{y} - h_\theta(\mathbf{y})\|^2 - \sigma^2 + \frac{2\sigma^2}{M}\operatorname{div}(h_\theta).$$

$\operatorname{div}(h_\theta) = \operatorname{tr}(\mathbf{A}\,\mathbf{J}_{f_\theta})$ where $\mathbf{J}_{f_\theta} \in \mathbb{R}^{N \times M}$ is the Jacobian of $f_\theta$. MC estimate: $\widehat{\operatorname{div}} = \mathbf{b}^\top \mathbf{A}\,\mathbf{J}_{f_\theta}\,\mathbf{b}$ with $\mathbf{b} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_M)$.

GSURE estimates $\mathbb{E}[\|\mathbf{A}(f_\theta(\mathbf{y}) - \mathbf{x})\|^2]/M$, which depends only on the range-space error. If $\mathbf{v} \in \text{null}(\mathbf{A})$, then $\mathbf{A}(\hat{\mathbf{x}} + \mathbf{v}) = \mathbf{A}\hat{\mathbf{x}}$, so GSURE is identical for $\hat{\mathbf{x}}$ and $\hat{\mathbf{x}} + \mathbf{v}$. Additional regularisation (TV, EI, prior) is needed. $\square$

Noise2Void masks pixel $i$ and predicts it from $\{y_j : j \neq i\}$. The loss is $|f_\theta(\mathbf{y}_{\setminus i})_i - y_i|^2$. The optimal predictor minimises this by predicting $\mathbb{E}[y_i \mid \{y_j : j \neq i\}]$.

If noise is independent: $\mathbb{E}[y_i \mid \{y_j\}_{j \neq i}] = \mathbb{E}[x_i \mid \{y_j\}_{j \neq i}]$ (the network learns to denoise).

If noise is correlated: $\mathbb{E}[y_i \mid \{y_j\}_{j \neq i}]$ includes a noise prediction component, so the network learns to partially reproduce noise rather than fully remove it.

After matched filtering in RF imaging, the noise is convolved with the point spread function, creating correlations over the mainlobe width ($\sim \lambda/(2\sin\theta_{\max})$). A N2V network can predict the masked pixel's noise from its correlated neighbours, resulting in noisy reconstructions. $\square$

Supervised gradient: $\nabla_\theta \|f_\theta(\mathbf{y}_1) - \mathbf{x}\|^2 = 2\mathbf{J}_{f_\theta}^\top(f_\theta(\mathbf{y}_1) - \mathbf{x})$.

N2N gradient: $\nabla_\theta \|f_\theta(\mathbf{y}_1) - \mathbf{y}_2\|^2 = 2\mathbf{J}_{f_\theta}^\top(f_\theta(\mathbf{y}_1) - \mathbf{x} - \mathbf{w}_2)$.

The N2N gradient has an additional term $-2\mathbf{J}_{f_\theta}^\top\mathbf{w}_2$.

$\text{Var}[\nabla\mathcal{L}_{\text{N2N}}] = \text{Var}[\nabla\mathcal{L}_{\text{sup}}] + 4\sigma^2\mathbb{E}[\|\mathbf{J}_{f_\theta}\|_F^2]$.

The extra variance is proportional to $\sigma^2$ and the Jacobian norm.

Higher gradient variance means N2N requires: (1) smaller learning rate, (2) larger batch size, or (3) more training iterations to achieve the same convergence. Typically, N2N needs $\sim 2$--$5\times$ more iterations than supervised training, but the per-iteration cost is the same. $\square$

$\mathcal{L}(\theta) = \frac{1}{2}\|\mathbf{y} - \mathbf{A}f_\theta(\mathbf{z})\|^2 + \lambda\sum_i \|\nabla [f_\theta(\mathbf{z})]_i\|_2.$$

DIP's spectral bias learns low frequencies first. TV penalises high-frequency gradients, reinforcing this tendency. The combined effect:

• Early iterations: DIP learns low-frequency structure; TV is inactive (low gradient values).
• Mid iterations: DIP starts fitting mid-frequency details; TV selectively preserves edges while suppressing oscillations.
• Late iterations: TV prevents the network from fitting high-frequency noise, extending the useful training window.

TV significantly reduces the need for early stopping: the PSNR curve has a plateau (from $\sim 3{,}000$ to $\sim 8{,}000$ iterations) rather than a sharp peak ($\sim 3{,}200$ without TV). However, for very long runs ($>10{,}000$ iterations), the network can still overfit in the TV-null directions (constant regions). In practice, DIP+TV is much more robust to the stopping criterion than vanilla DIP. $\square$

For shift $\Delta$: $\mathbf{A}T_\Delta\mathbf{x} = \mathbf{P}_\Omega\operatorname{diag}(e^{-j2\pi k\Delta/N})\hat{\mathbf{x}}$. The EI constraint $f_\theta(\mathbf{A}T_\Delta\hat{\mathbf{x}}) = T_\Delta\hat{\mathbf{x}}$ requires correct reconstruction from these modulated measurements.

With $N$ shifts $\Delta = 0, 1, \ldots, N-1$, the EI constraints create a system of $N \times M$ equations. For each measured frequency $k \in \Omega$, the $N$ shifts provide $N$ modulated observations with phases $e^{-j2\pi k\Delta/N}$ for $\Delta = 0, \ldots, N-1$.

The reconstruction must be consistent with all these modulated views, which constrains not just the measured coefficients but also the unmeasured ones through the network.

The system has $NM$ equations for $N$ unknowns (the Fourier coefficients). Since $M \geq 1$, this is over-determined, and the unique solution (for a perfect network) is the true signal. The condition for exact recovery is that shifts generate enough "virtual diversity" to span all frequencies through the reconstruction network. $\square$

For $Y \sim \text{Poisson}(\lambda)$: $\mathbb{E}[\lambda \cdot g(Y)] = \mathbb{E}[Y \cdot g(Y)]$ for any $g$ with $\mathbb{E}[|Y \cdot g(Y)|] < \infty$.

Equivalently: $\mathbb{E}[\lambda \cdot (g(Y+1) - g(Y))] = \mathbb{E}[(Y - \lambda) \cdot g(Y)]$.

For denoiser $f_\theta$ applied to Poisson data: $$\text{PURE}(f_\theta) = \frac{1}{N}\sum_i \bigl[f_\theta(\mathbf{y})_i - y_i\bigr]^2 + \frac{2}{N}\sum_i y_i\bigl[f_\theta(\mathbf{y} + \mathbf{e}_i)_i - f_\theta(\mathbf{y})_i\bigr] - \frac{1}{N}\|\mathbf{y}\|_1$$ where $\mathbf{e}_i$ is the $i$-th canonical vector.

The divergence term $y_i[f_\theta(\mathbf{y} + \mathbf{e}_i)_i - f_\theta(\mathbf{y})_i]$ is a finite-difference approximation that replaces the Gaussian divergence.

PURE requires $N$ forward passes (one per pixel perturbation), which is much more expensive than MC-SURE (one extra backward pass). Efficient approximations use random subsets of pixels. $\square$

Network $f_\theta: \mathbb{R}^d \to \mathbb{R}^{2N}$ outputs 2 channels: $[\operatorname{Re}(\hat{x}); \operatorname{Im}(\hat{x})]$. The complex image is $\hat{\mathbf{x}} = f_\theta^{(1)}(\mathbf{z}) + jf_\theta^{(2)}(\mathbf{z})$.

For multiplicative speckle with Goodman model: $$\mathcal{L}(\theta) = \sum_i \left[\frac{|y_i|}{|\hat{x}_i|^2} + 2\log|\hat{x}_i|\right] + \mu\,\|\mathbf{P}_\Omega\mathbf{F}\hat{\mathbf{x}} - \mathbf{y}_{\text{raw}}\|^2$$ combining the negative log-likelihood with data fidelity on raw measurements.

Add a phase consistency loss: $\|\angle(\mathbf{A}\hat{\mathbf{x}}) - \angle(\mathbf{y})\|^2$ on the measured frequencies. The complex DIP naturally preserves phase through the 2-channel representation.

Real DIP on magnitude: PSNR $\approx 26$ dB, no phase. Complex DIP: PSNR $\approx 27.5$ dB, phase error $< 15°$ on strong scatterers. The improvement comes from jointly optimising magnitude and phase. $\square$

$\mathcal{L} = \underbrace{\|\mathbf{P}_{\Omega_2}\mathbf{F}f_\theta(\mathbf{y}_1) - \mathbf{y}_2\|^2}_{\text{measurement splitting}} + \lambda_{\text{EI}}\underbrace{\sum_k \|f_\theta(\mathbf{A}T_{g_k}f_\theta(\mathbf{y}_1)) - T_{g_k}f_\theta(\mathbf{y}_1)\|^2}_{\text{equivariance}} + \lambda_{\text{DC}}\underbrace{\|\mathbf{A}f_\theta(\mathbf{y}_1) - \mathbf{y}\|^2}_{\text{data consistency}}.$$

• Measurement splitting constrains frequencies in $\Omega_2$ (cross-validates on held-out measurements).
• Data consistency constrains frequencies in $\Omega = \Omega_1 \cup \Omega_2$.
• Equivariance constrains frequencies in $\Omega^c$ (null space) via symmetry-induced virtual measurements.

Together, they provide supervision for all frequency components.

The combination is more robust than either method alone: measurement splitting handles non-symmetric scenes, while EI handles scenes where the measurement split is too sparse. $\square$

Encode $\mathbf{A}$ via its truncated SVD: $\mathbf{A} \approx \mathbf{U}_r\boldsymbol{\Sigma}_r\mathbf{V}_r^H$. Feed the singular values $\{\sigma_1, \ldots, \sigma_r\}$ and a low-dimensional representation of $\mathbf{U}_r, \mathbf{V}_r$ to the network as side information.

The network architecture: $f_\theta(\mathbf{y}, \text{SVD}(\mathbf{A}))$ uses cross-attention or FiLM conditioning to modulate features based on the operator.

For shift-invariant operators, encode $\mathbf{A}$ via its point spread function (PSF). The PSF is a compact 2D/3D function that captures the operator's spatial characteristics. This is natural for RF imaging where the PSF depends on array geometry and frequency.

Train an operator encoder $\psi(\mathbf{A})$ that maps any linear operator to a fixed-dimensional embedding. During RAM pretraining, the encoder learns to extract operator-relevant features (rank, condition number, spatial frequency coverage).

For new operators not seen during training: (1) compute the conditioning vector, (2) run the foundation model for an initial reconstruction, (3) fine-tune using SURE or EI losses for 50--100 iterations. This combines the foundation model's broad prior with operator-specific adaptation. $\square$

Test set: 100 simulated RF scenes (50 sparse point-scatterers, 50 extended). Forward model: partial Fourier with compression ratio $M/N = 0.25$. SNR: 10, 20, 30 dB.

• Low SNR (10 dB), sparse: DIP (strong architectural bias for sparse signals)
• Low SNR (10 dB), extended: EI (symmetries provide null-space info)
• High SNR (30 dB), sparse: SURE+PnP (denoiser quality dominates)
• High SNR (30 dB), extended: N2N (strongest supervision among self-supervised)
• Domain shift at test time: FM with LoRA (robust to distribution changes)

Primary: PSNR (magnitude), SSIM, phase error (degrees). Secondary: computational time, memory, hyperparameter sensitivity. $\square$

Unroll $K$ ADMM iterations: $\hat{\mathbf{x}}_K = \text{ADMM}_K(\mathbf{y}; f_\theta)$ where the denoiser $f_\theta$ is applied at each proximal step. The end-to-end map is $g_\theta(\mathbf{y}) = \text{ADMM}_K(\mathbf{y}; f_\theta)$.

Apply GSURE to the end-to-end reconstruction: $$\text{GSURE}(g_\theta) = \frac{1}{M}\|\mathbf{y} - \mathbf{A}g_\theta(\mathbf{y})\|^2 - \sigma^2 + \frac{2\sigma^2}{M}\operatorname{div}_\mathbf{y}(\mathbf{A}g_\theta).$$

The divergence $\operatorname{div}(\mathbf{A}g_\theta)$ is computed via backpropagation through the entire unrolled ADMM (MC estimate with one probe vector).

Two nested optimisation loops interact:

• Outer loop: gradient descent on $\theta$ to minimise GSURE.
• Inner loop: $K$ ADMM iterations for reconstruction.

For convergence: (1) the ADMM must converge for each fixed $\theta$ (requires the denoiser to be firm non-expansive), (2) the GSURE gradient must be unbiased (requires Gaussian noise). In practice, $K = 5$--$10$ iterations and the denoiser is regularised for non-expansiveness via spectral normalisation. $\square$

Each Tx-Rx pair provides $\mathbf{y}_{i,j} = \mathbf{A}_{i,j}\mathbf{c} + \mathbf{w}_{i,j}$ where $\mathbf{A}_{i,j}$ depends on the steering vectors $\mathbf{a}(\phi_i, \theta_i)$ and $\hat{\mathbf{a}}(\hat{\phi}_j, \hat{\theta}_j)$. The combined system has $IJK$ measurements for $Q$ voxels.

For scenes invariant under rotations by angle $\alpha$: the rotation $T_\alpha$ permutes the voxel grid, and the EI loss enforces $f_\theta(\mathbf{A}T_\alpha f_\theta(\mathbf{y})) = T_\alpha f_\theta(\mathbf{y})$.

For the rotation to mix range and null spaces, the array must not have the same rotational symmetry as the scene (otherwise the measurements are invariant under rotation, providing no information).

If the array has $P$-fold symmetry (e.g., UPA with $P = 4$), the measurements of a rotated scene can be computed from the measurements of the original scene by permuting Tx-Rx indices. This provides "free" data augmentation that does not require re-measurement.

Full recovery requires the group action to be transitive on the null space. For a UPA with $N_t \times N_r$ elements at wavelength $\lambda$ and scene diameter $D$, the null space has dimension $Q - \text{rank}(\mathbf{A})$. The number of independent constraints from $K$ group elements is $\sim K \cdot \text{rank}(\mathbf{A}) \cdot (1 - \text{overlap})$ where overlap measures redundancy between transformed measurements. Recovery requires this to exceed $Q - \text{rank}(\mathbf{A})$. $\square$