Exercises

$\log p(\mathbf{x}) = -\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu}) + \text{const}.$

$\nabla_\mathbf{x}\log p(\mathbf{x}) = -\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu}).$

The score points from $\mathbf{x}$ back toward the mean $\boldsymbol{\mu}$, scaled by the inverse covariance. $\blacksquare$

With the linear schedule, $\bar{\alpha}_t$ decreases from $\sim 1$ to $\sim 0$. Evaluating numerically: $\bar{\alpha}_{t^*} = 0.5$ occurs at $t^* \approx 520$.

Before this point, the signal dominates; after, noise dominates. This is the crossover point of the diffusion process.

The first half of the diffusion process ($t < 520$) operates in the high-SNR regime where the score network sees mostly signal; the second half ($t > 520$) operates in the low-SNR regime where the score must extrapolate from noise. The cosine schedule shifts this crossover to later steps, spending more time in the informative regime.

$\mathbf{x}_t = \sqrt{\bar{\alpha}_t}\mathbf{x}_0 + \sqrt{1-\bar{\alpha}_t}\boldsymbol{\epsilon}$, so $\mathbb{E}[\mathbf{x}_t] = \sqrt{\bar{\alpha}_t}\boldsymbol{\mu}$, $\text{Var}(\mathbf{x}_t) = \bar{\alpha}_t\sigma_0^2 + (1-\bar{\alpha}_t)$, $\text{Cov}(\mathbf{x}_0, \mathbf{x}_t) = \sqrt{\bar{\alpha}_t}\sigma_0^2$.

$\mathbb{E}[\mathbf{x}_0 \mid \mathbf{x}_t] = \boldsymbol{\mu} + \frac{\sqrt{\bar{\alpha}_t}\sigma_0^2}{\bar{\alpha}_t\sigma_0^2 + 1-\bar{\alpha}_t}(\mathbf{x}_t - \sqrt{\bar{\alpha}_t}\boldsymbol{\mu}).$

The marginal $p_t(\mathbf{x}_t) = \mathcal{N}(\sqrt{\bar{\alpha}_t}\boldsymbol{\mu}, (\bar{\alpha}_t\sigma_0^2 + 1-\bar{\alpha}_t)\mathbf{I})$, so $\nabla\log p_t = -(\bar{\alpha}_t\sigma_0^2+1-\bar{\alpha}_t)^{-1}(\mathbf{x}_t-\sqrt{\bar{\alpha}_t}\boldsymbol{\mu})$. Substituting into Tweedie's formula yields the same expression. $\checkmark$

The DPS gradient step has effective step size $\lambda_{\text{eff}} = \zeta / \sigma^2_{n}$. This plays the same role as $1/\rho$ in PnP-ADMM, where $\rho$ is the penalty parameter.

In PnP-ADMM: $\mathbf{x}^{k+1} = D_\sigma(\mathbf{z}^k)$ (denoiser step), $\mathbf{z}^{k+1} = \arg\min \frac{1}{2\sigma^2_{n}}\|\mathbf{A}\mathbf{z}-\mathbf{y}\|^2 + \frac{\rho}{2}\|\mathbf{z}-\mathbf{x}^{k+1}\|^2$ (data step).

In DPS: the score network provides the denoiser, and the guidance gradient provides the data step. The ratio $\zeta/\sigma^2_{n}$ corresponds to $1/\rho$ — both control the balance between prior and data fidelity.

1. Higher reconstruction quality ($+1$--$3$ dB PSNR)
2. Posterior sampling enables uncertainty quantification
3. Stronger prior captures complex scene statistics

1. Higher computational cost ($2$--$10\times$ more NFEs)
2. Requires more training data for the score network
3. Approximate posterior — no convergence guarantees

$\log p(\mathbf{y} \mid \hat{\mathbf{x}}_0) = -\frac{1}{2\sigma^2_{n}}\|\mathbf{y} - f(\hat{\mathbf{x}}_0)\|^2 + \text{const}.$

$\nabla_{\mathbf{x}_t}\log p(\mathbf{y} \mid \hat{\mathbf{x}}_0) = \frac{1}{\sigma^2_{n}}\frac{\partial\hat{\mathbf{x}}_0}{\partial\mathbf{x}_t}^T\,\mathbf{J}_f^T(\hat{\mathbf{x}}_0)\,(\mathbf{y} - f(\hat{\mathbf{x}}_0)),$$where$\mathbf{J}_f(\hat{\mathbf{x}}0) = \partial f/\partial\mathbf{x}\big|{\hat{\mathbf{x}}_0}$is the Jacobian of the forward model. In practice, this is computed via automatic differentiation through both$f$ and the Tweedie estimate.

For $f(\mathbf{x}) = \mathbf{A}\mathbf{x}$, the Jacobian is $\mathbf{J}_f = \mathbf{A}$, recovering the linear DPS gradient. $\blacksquare$

$r = m = n/4$ (range space dimension), $n - r = 3n/4$ (null space dimension).

Only $25\%$ of the reconstruction content is determined by the measurements. The remaining $75\%$ must be inferred by the diffusion prior. This makes the quality of the prior critical: a weak prior produces artefacts in the null-space components, while a strong prior fills in plausible content.

For RF imaging at 75% undersampling, the prior must capture scene-specific statistics (point scatterers, clutter) — a natural-image prior would hallucinate inappropriate textures.

$\mathbf{A}\hat{\mathbf{x}}_0^{\text{DDNM}} = \mathbf{A}\hat{\mathbf{x}}_0 - \mathbf{A}\mathbf{A}^\dagger(\mathbf{A}\hat{\mathbf{x}}_0 - \mathbf{y}).$$

For full row rank: $\mathbf{A}\mathbf{A}^\dagger = \mathbf{I}_m$, so

$$= \mathbf{A}\hat{\mathbf{x}}_0 - (\mathbf{A}\hat{\mathbf{x}}_0 - \mathbf{y}) = \mathbf{y}. \quad\blacksquare$$

The correction $\mathbf{A}^\dagger(\mathbf{A}\hat{\mathbf{x}}_0 - \mathbf{y})$ lies in the row space of $\mathbf{A}$. Therefore, the null-space component of $\hat{\mathbf{x}}_0$ is preserved: $(\mathbf{I} - \mathbf{A}^\dagger\mathbf{A})\hat{\mathbf{x}}_0^{\text{DDNM}} = (\mathbf{I} - \mathbf{A}^\dagger\mathbf{A})\hat{\mathbf{x}}_0$.

• DPS: $3C_{\text{net}} + O(mn)$ (forward, backward, matrix-vector)
• DDRM: $C_{\text{net}} + O(rn)$ (forward, projection via precomputed SVD, where $r = \text{rank}(\mathbf{A})$)

DDRM is $\sim 3\times$ faster per step.

The SVD of $\mathbf{A}$ costs $O(\min(m,n)^2\max(m,n))$. This amortises over all $S$ steps, adding $O(\min(m,n)^2\max(m,n)/S)$ per step.

DDRM is overall more efficient when $S \gg \min(m,n)^2\max(m,n) / (2C_{\text{net}})$. For typical imaging problems with moderate $m, n$ ($< 10^4$) and $S \geq 100$, DDRM wins. For large-scale RF imaging ($m, n > 10^5$), the SVD is prohibitive and DPS is preferred.

$\mathbf{z}_k \to (\mathbf{A}^{H}\mathbf{A})^{-1}\mathbf{A}^{H}\mathbf{y} = \mathbf{A}^\dagger\mathbf{y}$ (pseudoinverse, pure data fidelity, no prior).

$\mathbf{z}_k \to \rho^{-1}(\rho^{-1}\mathbf{A}^{H}\mathbf{A} + \mathbf{I})^{-1}(\rho^{-1}\mathbf{A}^{H}\mathbf{y} + \hat{\mathbf{x}}_k) \to \hat{\mathbf{x}}_k$ (pure prior, no data fidelity).

The parameter $\rho$ interpolates between data fidelity and prior. In DiffPIR, $\rho$ is scheduled to decrease over iterations: early iterations rely on the prior (high $\rho$, low noise), later iterations enforce data fidelity (low $\rho$). $\blacksquare$

$\log\text{SNR}(t) = \log\bar{\alpha}_t - \log(1-\bar{\alpha}_t)$. For the linear schedule, $\bar{\alpha}_t$ decreases roughly exponentially, making $\log\text{SNR}$ approximately linear in $t$.

If $\log\text{SNR}$ is linear in $t$, then uniform steps $\Delta t = T/S$ give uniform steps $\Delta\log\text{SNR} \approx (\log\text{SNR}(T) - \log\text{SNR}(0))/S$.

Each DDIM step covers the same "amount of denoising" in the log-SNR sense. This is why uniform time subsampling works well for the linear schedule. For non-linear schedules (cosine), non-uniform subsampling may be preferred.

The guidance gradient has magnitude $\propto \zeta\|\mathbf{A}\hat{\mathbf{x}}_0 - \mathbf{y}\|/\sigma^2_{n}$. As $\sigma^2_{n} \to 0$, this diverges unless $\mathbf{A}\hat{\mathbf{x}}_0 = \mathbf{y}$.

In the noiseless limit, DPS requires the Tweedie estimate to be exactly measurement-consistent at every step, which is unrealistic. In practice, $\zeta$ must be decreased as $\sigma^2_{n}$ decreases, keeping $\zeta/\sigma^2_{n}$ bounded. A common heuristic: $\zeta = \zeta_0 \cdot \sigma^2_{n}$, making the effective step size $\zeta_0$ independent of noise level.

$\mathbf{y} \mid \mathbf{x}_t$: $\mathbf{x}_0 \sim \mathcal{N}(\hat{\mathbf{x}}_0, r_t^2\mathbf{I})$, $\mathbf{y} \mid \mathbf{x}_0 \sim \mathcal{N}(\mathbf{A}\mathbf{x}_0, \sigma^2_{n}\mathbf{I})$. By the convolution property: $\mathbf{y} \mid \mathbf{x}_t \sim \mathcal{N}(\mathbf{A}\hat{\mathbf{x}}_0, \sigma^2_{n}\mathbf{I} + r_t^2\mathbf{A}\mathbf{A}^{H})$.

$\nabla_{\hat{\mathbf{x}}_0}\log p(\mathbf{y} \mid \mathbf{x}_t) = \mathbf{A}^{H}(\sigma^2_{n}\mathbf{I} + r_t^2\mathbf{A}\mathbf{A}^{H})^{-1}(\mathbf{y} - \mathbf{A}\hat{\mathbf{x}}_0).$$

Apply the chain rule through the Tweedie estimate: $$\nabla_{\mathbf{x}_t}\log p(\mathbf{y} \mid \mathbf{x}_t) = \frac{\partial\hat{\mathbf{x}}_0}{\partial\mathbf{x}_t}^T\mathbf{A}^{H}(\sigma^2_{n}\mathbf{I} + r_t^2\mathbf{A}\mathbf{A}^{H})^{-1}(\mathbf{y} - \mathbf{A}\hat{\mathbf{x}}_0).$$

At large $t$: $r_t^2$ is large, the covariance is dominated by $r_t^2\mathbf{A}\mathbf{A}^{H}$, and guidance is weak (prior dominates). At small $t$: $r_t^2 \to 0$, recovering the DPS gradient. $\blacksquare$

$\mathbf{A}\hat{\mathbf{x}} = \mathbf{U}_r\boldsymbol{\Sigma}_r\mathbf{V}_r^H(\mathbf{V}_r\boldsymbol{\Sigma}_r^{-1}\mathbf{U}_r^H\mathbf{y} + \mathbf{V}_\perp\boldsymbol{\eta}).$$

Using $\mathbf{V}_r^H\mathbf{V}_r = \mathbf{I}_r$ and $\mathbf{V}_r^H\mathbf{V}_\perp = \mathbf{0}$:

$$= \mathbf{U}_r\boldsymbol{\Sigma}_r\boldsymbol{\Sigma}_r^{-1}\mathbf{U}_r^H\mathbf{y} = \mathbf{U}_r\mathbf{U}_r^H\mathbf{y}.$$

Since $\mathbf{y} \in \text{range}(\mathbf{A}) = \text{span}(\mathbf{U}_r)$, we have $\mathbf{U}_r\mathbf{U}_r^H\mathbf{y} = \mathbf{y}$. $\blacksquare$

DPM-Solver++ with $S = 25$: reduces time from 60s to $60 \times 25/200 = 7.5$ seconds. Still above 5s.

DPM-Solver++ ($S = 20$) + gradient checkpointing (saves memory, costs $\sim 30\%$ more compute per step): $60 \times 20/200 \times 1.0 = 6$ seconds.

Alternatively: DPM-Solver++ ($S = 20$) without backpropagation (replace DPS guidance with DDNM-style projection): $60 \times 20/200 \times (1/3) = 2$ seconds. This sacrifices the gradient-based guidance but is sufficient for a structured $\mathbf{A}$.

For $\mathbf{A}$ with known SVD: DDNM + DPM-Solver++ ($S = 20$) achieves $\sim 2$ seconds with exact measurement consistency. For general $\mathbf{A}$: DPS + DPM-Solver++ ($S = 25$) achieves $\sim 5$ seconds with approximate consistency.

The DDIM trajectory in continuous time follows: $d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2\nabla\log p_t(\mathbf{x})]dt$.

With guidance: $d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2(\nabla\log p_t(\mathbf{x}) + \frac{\zeta}{\sigma^2_{n}}\nabla\log p(\mathbf{y}|\hat{\mathbf{x}}_0))]dt$.

As $t \to 0$: $p_t \to p_0 = p$ (data distribution), $\hat{\mathbf{x}}_0 \to \mathbf{x}$ (Tweedie estimate converges to the observation). The ODE becomes a gradient flow on: $\nabla[-\log p(\mathbf{x}) - \frac{\zeta}{2\sigma^2_{n}}\|\mathbf{y}-\mathbf{A}\mathbf{x}\|^2]$, which is the MAP objective with regularisation weight $\zeta$. $\blacksquare$

Define $\tilde{\mathbf{x}} = [\text{Re}(\mathbf{c})^T, \text{Im}(\mathbf{c})^T]^T \in \mathbb{R}^{2n}$, $\tilde{\mathbf{y}} = [\text{Re}(\mathbf{y})^T, \text{Im}(\mathbf{y})^T]^T \in \mathbb{R}^{2m}$, and the real-valued sensing matrix: $$\tilde{\mathbf{A}} = \begin{bmatrix} \mathbf{A}_{R} & -\mathbf{A}_{I} \\ \mathbf{A}_{I} & \mathbf{A}_{R} \end{bmatrix} \in \mathbb{R}^{2m \times 2n}.$$

The complex measurement $\mathbf{y} = \mathbf{A}\mathbf{c} + \mathbf{w}$ is equivalent to $\tilde{\mathbf{y}} = \tilde{\mathbf{A}}\tilde{\mathbf{x}} + \tilde{\mathbf{w}}$.

DPS on the 2-channel representation with $\tilde{\mathbf{A}}$ is mathematically identical to complex-valued DPS. The score network operates on $\tilde{\mathbf{x}} \in \mathbb{R}^{2n}$ and the guidance gradient is computed via the real-valued chain rule. $\blacksquare$

$\mathcal{L}(\theta) = \underbrace{\mathbb{E}_{t,\mathbf{x}_0,\boldsymbol{\epsilon}}\|\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_\theta(\mathbf{x}_t, t)\|^2}_{\mathcal{L}_{\text{DSM}}} + \lambda\underbrace{\mathbb{E}_{t,\mathbf{x}_0,\boldsymbol{\epsilon}}\|\mathbf{y} - \mathbf{A}\hat{\mathbf{x}}_0(\mathbf{x}_t;\theta)\|^2}_{\mathcal{L}_{\text{physics}}}$$where$\hat{\mathbf{x}}_0(\mathbf{x}_t;\theta) = (\mathbf{x}_t - \sqrt{1-\bar{\alpha}t}\boldsymbol{\epsilon}\theta(\mathbf{x}_t,t))/\sqrt{\bar{\alpha}_t}$.

$\nabla_\theta\mathcal{L}_{\text{physics}} = -\frac{2\lambda}{\sqrt{\bar{\alpha}_t}}\mathbb{E}\left[(\mathbf{y}-\mathbf{A}\hat{\mathbf{x}}_0)^T\mathbf{A}\frac{\sqrt{1-\bar{\alpha}_t}}{\sqrt{\bar{\alpha}_t}}\nabla_\theta\boldsymbol{\epsilon}_\theta\right].$$

This requires backpropagation through the network (same as DPS guidance, but during training rather than inference).

1. Sample $(\mathbf{x}_0, \mathbf{y})$ from the simulation database
2. Sample $t, \boldsymbol{\epsilon}$; compute $\mathbf{x}_t$
3. Forward pass: $\boldsymbol{\epsilon}_\theta(\mathbf{x}_t, t)$
4. Compute $\hat{\mathbf{x}}_0$ via Tweedie
5. Compute $\mathcal{L} = \|\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_\theta\|^2 + \lambda\|\mathbf{y} - \mathbf{A}\hat{\mathbf{x}}_0\|^2$
6. Backpropagate and update $\theta$

Define $\mathbf{z}(t) = e^{-\int_T^t f(s)ds}\mathbf{x}(t)$. Then $d\mathbf{z}/dt = e^{-\int_T^t f(s)ds}\frac{g(t)^2}{2\sigma_t}\boldsymbol{\epsilon}_\theta$.

Approximate $\boldsymbol{\epsilon}_\theta(\mathbf{x}(s), s) \approx \boldsymbol{\epsilon}_\theta(\mathbf{x}(t), t)$ for $s \in [t-\Delta t, t]$. Integrating:

$$\mathbf{z}(t-\Delta t) - \mathbf{z}(t) = \boldsymbol{\epsilon}_\theta(\mathbf{x}(t), t)\int_{t}^{t-\Delta t}e^{-\int_T^s f(\tau)d\tau}\frac{g(s)^2}{2\sigma_s}ds.$$

For the VP-SDE with $f(t) = -\frac{1}{2}\beta(t)$ and $g(t) = \sqrt{\beta(t)}$, evaluating the integral and transforming back to $\mathbf{x}$ coordinates recovers the DDIM update formula. The approximation error is $O(\Delta t^2)$ (first-order), which explains why DDIM needs $\sim 50$--$200$ steps for good quality. $\blacksquare$

$\log p(\mathbf{y} \mid \hat{\mathbf{x}}_0) = \sum_{i=1}^m\left[y_i\log(\mathbf{A}_{i}^{T}\hat{\mathbf{x}}_0) - \mathbf{A}_{i}^{T}\hat{\mathbf{x}}_0 - \log(y_i!)\right].$$

$\nabla_{\hat{\mathbf{x}}_0}\log p(\mathbf{y} \mid \hat{\mathbf{x}}_0) = \sum_{i=1}^m\left(\frac{y_i}{\mathbf{A}_{i}^{T}\hat{\mathbf{x}}_0} - 1\right)\mathbf{A}_{i} = \mathbf{A}^{T}\left(\frac{\mathbf{y}}{\mathbf{A}\hat{\mathbf{x}}_0} - \mathbf{1}\right),$$

where the division is element-wise.

$\nabla_{\mathbf{x}_t}\log p(\mathbf{y} \mid \hat{\mathbf{x}}_0) = \frac{\partial\hat{\mathbf{x}}_0}{\partial\mathbf{x}_t}^T\mathbf{A}^{T}\left(\frac{\mathbf{y}}{\mathbf{A}\hat{\mathbf{x}}_0} - \mathbf{1}\right).$$

1. Numerical instability: When $\mathbf{A}_{i}^{T}\hat{\mathbf{x}}_0 \to 0$, the gradient diverges. Requires clipping or smoothing.
2. Non-negativity: $\mathbf{x}_0$ must be non-negative for the Poisson model to be valid. The diffusion model may generate negative values, requiring projection.
3. No closed-form proximal: Unlike the Gaussian case, there is no closed-form solution for the Poisson data-fidelity step in DiffPIR-style methods.