Exercises

For each component: $\min_{x_i} \tfrac{1}{2}(x_i - v_i)^2 + \lambda|x_i|$. Setting the subgradient to zero: $x_i - v_i + \lambda\,\partial|x_i| \ni 0$. Solution: $x_i = \operatorname{sign}(v_i)\max(|v_i| - \lambda, 0)$.

Laplace prior: $p(x_i) \propto e^{-|x_i|/b}$, so $-\log p(\mathbf{x}) = \tfrac{1}{b}\|\mathbf{x}\|_1$. MAP estimator: $\operatorname{prox}_{\sigma^2/b\,\|\cdot\|_1}(\mathbf{v})$. With $\lambda = \sigma^2/b$, this is soft-thresholding. $\square$

$\mathbf{x}^{(k+1)} = \mathcal{D}_\sigma\!\bigl(\mathbf{x}^{(k)} - \alpha\mathbf{A}^{H}(\mathbf{A}\mathbf{x}^{(k)} - \mathbf{y})\bigr).$$

The proximal operator at step size $\alpha$ solves $\min_\mathbf{x} \tfrac{1}{2}\|\mathbf{x} - \mathbf{v}\|^2 + \alpha\lambda R(\mathbf{x})$, which is MAP denoising at variance $\alpha\lambda$. Hence $\sigma^2 = \alpha\lambda$, or $\sigma = \sqrt{\alpha\lambda}$. $\square$

Case 1: $|a|, |b| > \tau$, same sign. $|\mathcal{S}_\tau(a) - \mathcal{S}_\tau(b)| = |a - b|$. ✓

Case 1b: Opposite signs. $|\mathcal{S}_\tau(a) - \mathcal{S}_\tau(b)| \leq |a| - \tau + |b| - \tau < |a| + |b| \leq |a - b|$. ✓

Case 2: $|a| \leq \tau$. Then $\mathcal{S}_\tau(a) = 0$, so $|\mathcal{S}_\tau(a) - \mathcal{S}_\tau(b)| = |\mathcal{S}_\tau(b)| \leq |b| - \tau \leq |b - a|$ (since $|a| \leq \tau$). ✓

$\|\mathcal{S}_\tau(\mathbf{a}) - \mathcal{S}_\tau(\mathbf{b})\|^2 = \sum_i |\mathcal{S}_\tau(a_i) - \mathcal{S}_\tau(b_i)|^2 \leq \sum_i |a_i - b_i|^2 = \|\mathbf{a} - \mathbf{b}\|^2$. $\square$

$(\mathbf{A}^{H}\mathbf{A} + \rho\mathbf{I})^{-1} = \mathbf{F}^H(\mathbf{D}_\Omega + \rho\mathbf{I})^{-1}\mathbf{F}$.

$\mathbf{c}^{(k+1)} = \mathbf{F}^H\!\left[ \frac{\mathbf{D}_\Omega\mathbf{F}\mathbf{A}^{H}\mathbf{y} + \rho\,\mathbf{F}(\mathbf{z}^{(k)} - \mathbf{u}^{(k)})} {\mathbf{D}_\Omega + \rho}\right].$$One FFT for$\mathbf{F}(\mathbf{z}^{(k)} - \mathbf{u}^{(k)})$, element-wise ops, one IFFT. Total:$2 \times O(N\log N)$.$\square$

DRUNet concatenates $\sigma$ as a constant channel, allowing one network to denoise at any noise level. For PnP, this maps directly to the ADMM penalty parameter: $\sigma = \sqrt{\lambda/\rho}$.

If $\sigma$ is too large (stronger denoising than needed), the denoiser over-smooths and destroys fine details. The reconstruction converges to an excessively blurry image.

If $\sigma$ is too small, the denoiser is too weak to suppress artefacts. The iterates may oscillate or produce a noisy reconstruction with residual measurement artefacts. $\square$

$\nabla R_\text{RED} = \frac{1}{2}\nabla\!\left[\|\mathbf{x}\|^2 - \mathbf{x}^T\mathcal{D}_\sigma(\mathbf{x})\right] = \mathbf{x} - \frac{1}{2}\!\left[\mathcal{D}_\sigma(\mathbf{x}) + \mathbf{J}_\mathcal{D}^T(\mathbf{x})\,\mathbf{x}\right].$$

If $\mathbf{J}_\mathcal{D} = \mathbf{J}_\mathcal{D}^T$ and local homogeneity $\mathbf{J}_\mathcal{D}\mathbf{x} \approx \mathcal{D}(\mathbf{x})$ holds: $\nabla R_\text{RED} = \mathbf{x} - \mathcal{D}_\sigma(\mathbf{x})$.

The exact gradient involves the vector-Jacobian product $\mathbf{J}_\mathcal{D}^T\mathbf{x}$, expensive via autodiff. Using the approximation $\nabla R \approx \mathbf{x} - \mathcal{D}_\sigma(\mathbf{x})$ introduces error $\approx \tfrac{1}{2}\|(\mathbf{J}_\mathcal{D}^T - \mathbf{J}_\mathcal{D})\mathbf{x}\|$. $\square$

$\gamma = (1/49)^{1/29} \approx 0.875$. $\sigma_k = 0.196 \times 0.875^{k-1}$. Decreases rapidly early; slow refinement at the end.

$\sigma_k = 0.004 + \tfrac{0.192}{2}(1 + \cos(\pi(k-1)/29))$. Decreases slowly at both endpoints, rapidly in the middle.

For RF imaging with structured artefacts, the cosine schedule is often preferred: it spends more iterations at intermediate noise levels where the denoiser effectively suppresses artefact-like noise. The geometric schedule converges faster but spends fewer iterations on the mid-range transition. $\square$

$\mathbf{p} = \operatorname{prox}_f(\mathbf{a})$ satisfies $\mathbf{a} - \mathbf{p} \in \partial f(\mathbf{p})$. $\mathbf{q} = \operatorname{prox}_f(\mathbf{b})$ satisfies $\mathbf{b} - \mathbf{q} \in \partial f(\mathbf{q})$.

Since $f$ is convex, $\partial f$ is monotone: $\langle (\mathbf{a} - \mathbf{p}) - (\mathbf{b} - \mathbf{q}),\; \mathbf{p} - \mathbf{q}\rangle \geq 0$.

Let $\mathbf{d} = \mathbf{p} - \mathbf{q}$, $\mathbf{e} = (\mathbf{a} - \mathbf{p}) - (\mathbf{b} - \mathbf{q})$. Monotonicity: $\langle \mathbf{e}, \mathbf{d}\rangle \geq 0$. $\mathbf{a} - \mathbf{b} = \mathbf{d} + \mathbf{e}$, so $\|\mathbf{a} - \mathbf{b}\|^2 = \|\mathbf{d}\|^2 + 2\langle\mathbf{d},\mathbf{e}\rangle + \|\mathbf{e}\|^2 \geq \|\mathbf{d}\|^2 + \|\mathbf{e}\|^2$, which is the desired inequality. $\square$

$[\mathbf{z}_1]_j = \operatorname{ReLU}([\mathbf{U}_1\mathbf{x} + \mathbf{b}_1]_j)$. The affine argument is convex. ReLU is convex and non-decreasing. Composition preserves convexity.

Assume each $[\mathbf{z}_{\ell-1}]_i$ is convex in $\mathbf{x}$. Since $[\mathbf{W}_\ell^{(z)}]_{ji} \geq 0$, the sum $\sum_i [\mathbf{W}_\ell^{(z)}]_{ji}[\mathbf{z}_{\ell-1}]_i$ is a non-negative weighted sum of convex functions, hence convex. Adding the affine term $[\mathbf{U}_\ell\mathbf{x} + \mathbf{b}_\ell]_j$ preserves convexity. Applying ReLU (convex, non-decreasing) preserves convexity.

$\Phi_\theta(\mathbf{x}) = \mathbf{w}^T\mathbf{z}_D + c$ with $\mathbf{w} \geq 0$ is a non-negative weighted sum of convex functions plus a constant: convex. $\square$

At $\mathbf{x}^*$: $\mathbf{x}^* = \mathbf{x}^* - \alpha\bigl[\mathbf{A}^{H}(\mathbf{A}\mathbf{x}^* - \mathbf{y}) + \lambda(\mathbf{x}^* - \mathcal{D}_\sigma(\mathbf{x}^*))\bigr]$. Subtracting $\mathbf{x}^*$ from both sides and dividing by $\alpha > 0$: $\mathbf{A}^{H}(\mathbf{A}\mathbf{x}^* - \mathbf{y}) + \lambda(\mathbf{x}^* - \mathcal{D}_\sigma(\mathbf{x}^*)) = \mathbf{0}$. $\square$

The $\mathbf{c}$-update uses complex arithmetic: $\mathbf{c}^{(k+1)} = (\mathbf{A}^{H}\mathbf{A} + \rho\mathbf{I})^{-1}(\mathbf{A}^{H}\mathbf{y} + \rho(\mathbf{z}^{(k)} - \mathbf{u}^{(k)}))$, which is a complex linear solve (same structure as real ADMM).

Stack: $\tilde{\mathbf{v}} = [\operatorname{Re}(\mathbf{v}); \operatorname{Im}(\mathbf{v})] \in \mathbb{R}^{2Q}$. Apply DRUNet as a 2-channel denoiser. Reconstruct: $\mathbf{z} = \tilde{z}_{1:Q} + j\tilde{z}_{Q+1:2Q}$.

$\|\mathcal{D}(\tilde{\mathbf{a}}) - \mathcal{D}(\tilde{\mathbf{b}})\|_{\mathbb{R}^{2Q}} \leq \|\tilde{\mathbf{a}} - \tilde{\mathbf{b}}\|_{\mathbb{R}^{2Q}} = \|\mathbf{a} - \mathbf{b}\|_{\mathbb{C}^Q}$. So non-expansiveness on $\mathbb{R}^{2Q}$ implies non-expansiveness on $\mathbb{C}^Q$. $\square$

Layer 1: $\mathbf{h}_1(\mathbf{x}) = \operatorname{ReLU}(\mathbf{W}_1\mathbf{x} + \mathbf{b}_1)$. $\|\mathbf{h}_1(\mathbf{x}) - \mathbf{h}_1(\mathbf{y})\| \leq \|\mathbf{W}_1(\mathbf{x}-\mathbf{y})\| \leq \|\mathbf{W}_1\|\|\mathbf{x}-\mathbf{y}\|$. By composition: $L \leq \|\mathbf{W}_2\|\|\mathbf{W}_1\|$.

Initialise $\tilde{\mathbf{u}} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. Each training step: $\tilde{\mathbf{v}} = \mathbf{W}^T\tilde{\mathbf{u}}/\|\mathbf{W}^T\tilde{\mathbf{u}}\|$, then $\tilde{\mathbf{u}} = \mathbf{W}\tilde{\mathbf{v}}/\|\mathbf{W}\tilde{\mathbf{v}}\|$. Spectral norm estimate: $\hat{\sigma}_1 = \tilde{\mathbf{u}}^T\mathbf{W}\tilde{\mathbf{v}}$.

Weight clipping ($|\mathbf{W}_{ij}| \leq c$) constrains individual entries, often destroying gradient signal. Spectral normalisation $\bar{\mathbf{W}} = \mathbf{W}/\hat{\sigma}_1$ divides by the largest singular value, preserving the weight's directional structure and allowing large entries in directions that do not expand the output norm. $\square$

$\|\nabla F(\mathbf{x}) - \nabla F(\mathbf{y})\| \leq (\|\mathbf{A}\|^2 + \lambda)\|\mathbf{x} - \mathbf{y}\|$ (assuming $\nabla R_\text{RED}$ is 1-Lipschitz, i.e., $\mathcal{D}_\sigma$ is 0-Lipschitz, which holds if $\mathcal{D}_\sigma$ is a constant). In general: $\beta = \|\mathbf{A}\|^2 + \lambda(1 + L_\mathcal{D})$ where $L_\mathcal{D}$ is the Lipschitz constant of $\mathcal{D}_\sigma$.

Required step size: $\alpha \leq 1/\beta$. Convergence rate for strongly convex $F$ with $\mu$-strong convexity: $F(\mathbf{x}^{(k)}) - F^* \leq (1 - \mu\alpha)^k(F(\mathbf{x}^{(0)}) - F^*)$. $\square$

Let $\mathbf{A} = \mathbf{I}$, $\mathbf{y} = \mathbf{0}$, $c = 1.5$, $\alpha = 0.1$.

$\mathbf{x}^{(k+1)} = 1.5(\mathbf{x}^{(k)} - 0.1\mathbf{x}^{(k)}) = 1.35\,\mathbf{x}^{(k)}$. The spectral radius is $1.35 > 1$, so $\|\mathbf{x}^{(k)}\| = 1.35^k\|\mathbf{x}^{(0)}\| \to \infty$. $\square$

OAMP's state evolution shows that the OAMP denoiser input is asymptotically Gaussian with a known variance. The BG MMSE denoiser computes $\mathbb{E}[\mathbf{x}_i | v_i]$ exactly, giving the minimum MSE component-wise. Convergence of state evolution to a fixed point corresponds to convergence of OAMP to the MMSE estimate.

A general deep denoiser (DRUNet) is trained on natural images and does not compute the BG MMSE denoiser. Even if the denoiser is optimal for natural images, it sub-optimally handles BG signals because the BG distribution differs from the natural image prior.

If the PnP denoiser is specifically the BG MMSE denoiser $\mathcal{D}_\sigma(\mathbf{v}) = \mathbb{E}_\text{BG}[\mathbf{x}|\mathbf{v}]$ at the correct variance $\sigma^2$, then PnP-ADMM converges to the same fixed point as OAMP and achieves the MMSE under state evolution. This is achieved by using a matched denoiser — but then PnP reduces to OAMP with state evolution tracking. $\square$

$\Phi_\theta(x_1, x_2) = \text{softplus}(\mathbf{u}^T\mathbf{x} + b) + \text{softplus}(-\mathbf{u}^T\mathbf{x} + b)$ where $\mathbf{u} = (1, -1)^T$ and $b = 0$. This gives $\Phi_\theta(\mathbf{x}) = 2\log(1 + e^{x_1 - x_2}/2)$.

$\text{softplus}(a) = \log(1 + e^a)$ is convex. A linear function of $\mathbf{x}$ is affine (convex). The composition of a convex function with an affine function is convex. The sum of convex functions is convex.

As $|x_1 - x_2| \to \infty$: $\Phi_\theta(\mathbf{x}) \to |x_1 - x_2|$ (exact). For small $|x_1 - x_2| \leq \epsilon$: $\Phi_\theta(\mathbf{x}) \approx \log(2) + (x_1 - x_2)^2/(4\log(2)) \geq 0$. Maximum error: $|\Phi_\theta - |x_1 - x_2|| \leq \log(2) \approx 0.693$, achieved at $x_1 = x_2$. $\square$

$\mathcal{D}_\sigma(\mathbf{x}) = \mathbf{x} + \sigma^2\nabla_\mathbf{x}\log p_\sigma(\mathbf{x}).KATEXPLACEHOLDER0END\mathbf{x} - \mathcal{D}_\sigma(\mathbf{x}) = -\sigma^2\nabla_\mathbf{x}\log p_\sigma(\mathbf{x}).$$

The RED update: $$\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} - \alpha[\underbrace{\nabla f(\mathbf{x}^{(k)})}_\text{data fidelity} + \lambda(\mathbf{x}^{(k)} - \mathcal{D}_\sigma(\mathbf{x}^{(k)}))].$$ Substituting: $$= \mathbf{x}^{(k)} - \alpha\nabla f(\mathbf{x}^{(k)}) + \alpha\lambda\sigma^2\nabla_\mathbf{x}\log p_\sigma(\mathbf{x}^{(k)}).$$ This is a data-fidelity gradient step plus a score-function ascent step in $\log p_\sigma$ (the smoothed image log-density).

The equivalence uses $p_\sigma(\mathbf{x}) \approx p(\mathbf{x})$ (image distribution evaluated at a clean point). This is accurate when $\sigma$ is small relative to the image structure, i.e., in late iterations when the estimate is close to the truth. $\square$

$\rho = \lambda/\sigma_1^2 = 0.05/0.05^2 = 20$. At this $\rho$, $\sqrt{0.05/20} = 0.05 = \sigma_1$. ✓

• ADMM linear solve (via random $\mathbf{A}$): $O(MQ)$ (matrix–vector), $\approx 512 \times 1024 \approx 5 \times 10^5$ ops
• DRUNet forward pass (CPU, 4-scale U-Net): $\approx 10^8$ ops per image

Total per-iteration: $\approx 10^8$ flops (dominated by DRUNet). 40 iterations: $\approx 4 \times 10^9$ flops.

FISTA per iteration: one matrix–vector multiply ($O(MQ) \approx 5\times10^5$) + soft-threshold. 200 iterations: $\approx 10^8$ flops (no DRUNet).

PnP costs $\sim 40\times$ more per run, but typically achieves $2$–$3$ dB better MSE for structured scenes. The tradeoff depends on the application's tolerance for computation time. $\square$