Exercises

$p(y \mid \gamma) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\left(-\frac{(y-a\gamma)^2}{2\sigma^2}\right)$, $\pi(\gamma) = \frac{1}{\sqrt{2\pi\gamma_0^2}}\exp\!\left(-\frac{\gamma^2}{2\gamma_0^2}\right)$.

$-\log p(\gamma \mid y) \propto \frac{(y-a\gamma)^2}{2\sigma^2} + \frac{\gamma^2}{2\gamma_0^2} = \frac{1}{2}\left(\frac{a^2}{\sigma^2} + \frac{1}{\gamma_0^2}\right)\gamma^2 - \frac{ay}{\sigma^2}\gamma + \text{const}$.

This is $\frac{(\gamma - \hat{\gamma}_{\text{post}})^2}{2\sigma_{\text{post}}^2}$ with $\sigma_{\text{post}}^{-2} = a^2/\sigma^2 + 1/\gamma_0^2$ and $\hat{\gamma}_{\text{post}} = \sigma_{\text{post}}^2 \cdot ay/\sigma^2$.

Substituting: $\hat{\gamma}_{\text{post}} = \frac{a\gamma_0^2}{a^2\gamma_0^2 + \sigma^2}y = \frac{\text{SNR}}{1 + \text{SNR}} \cdot \frac{y}{a}$ where $\text{SNR} = a^2\gamma_0^2/\sigma^2$.

As SNR $\to \infty$: $\hat{\gamma}_{\text{post}} \to y/a$ (ML estimate). As SNR $\to 0$: $\hat{\gamma}_{\text{post}} \to 0$ (prior mean). $\blacksquare$

$-\log p(\boldsymbol{\gamma} \mid \mathbf{y}) \propto \frac{1}{2\sigma^2}\|\mathbf{A}\boldsymbol{\gamma} - \mathbf{y}\|^2 + \frac{1}{2\gamma_0^2}\|\boldsymbol{\gamma}\|^2$. Setting gradient to zero: $\frac{1}{\sigma^2}\mathbf{A}^H(\mathbf{A}\boldsymbol{\gamma} - \mathbf{y}) + \frac{1}{\gamma_0^2}\boldsymbol{\gamma} = 0$, giving $(\mathbf{A}^H\mathbf{A} + \lambda\mathbf{I})\hat{\boldsymbol{\gamma}} = \mathbf{A}^H\mathbf{y}$ with $\lambda = \sigma^2/\gamma_0^2$.

Let $\mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^H$. Then in the SVD basis: $\hat{\boldsymbol{\gamma}}_{\text{MAP}} = \mathbf{V}\,\text{diag}\!\left(\frac{\sigma_k}{\sigma_k^2 + \lambda}\right)\mathbf{U}^H\mathbf{y}$, with filter factors $f_k = \sigma_k^2/(\sigma_k^2 + \lambda)$ — exactly Tikhonov spectral regularization.

As $\gamma_0^2 \to \infty$: $\lambda \to 0$, $f_k \to 1$ — the MAP approaches the pseudoinverse (no regularization). As $\gamma_0^2 \to 0$: $\lambda \to \infty$, $f_k \to 0$ — the MAP approaches $\mathbf{0}$ (prior completely dominates).

$\mathbf{\Gamma}_{\text{post}} = (\sigma^{-2}\mathbf{A}^H\mathbf{A} + \gamma_0^{-2}\mathbf{I})^{-1}$. In the SVD basis of $\mathbf{A}$, $\mathbf{A}^H\mathbf{A} = \mathbf{V}\boldsymbol{\Sigma}^H\boldsymbol{\Sigma}\mathbf{V}^H$ with diagonal entries $\sigma_k^2$. Therefore $[\mathbf{\Gamma}_{\text{post}}]_{kk} = (\sigma_k^2/\sigma^2 + 1/\gamma_0^2)^{-1}$.

For null-space directions: $\sigma_k = 0$, so $[\mathbf{\Gamma}_{\text{post}}]_{kk} = (0 + 1/\gamma_0^2)^{-1} = \gamma_0^2$. The posterior variance equals the prior variance — data provides no information about null-space components.

For $\sigma_k \gg \sigma/\gamma_0$: $[\mathbf{\Gamma}_{\text{post}}]_{kk} \approx ({\sigma_k^2/\sigma^2})^{-1} = \sigma^2/\sigma_k^2$. Uncertainty in mode $k$ is determined by the noise level $\sigma^2$ relative to the singular value $\sigma_k^2$. $\blacksquare$

$-\log p(\gamma|y) \propto (y-\gamma)^2/(2\sigma^2) + \lambda|\gamma|$. Subgradient optimality: $(\gamma-y)/\sigma^2 + \lambda\,\partial|\gamma| \ni 0$. For $\gamma > 0$: $\gamma = y - \lambda\sigma^2 > 0$ iff $y > \lambda\sigma^2$. For $\gamma < 0$: $\gamma = y + \lambda\sigma^2 < 0$ iff $y < -\lambda\sigma^2$. For $|y| \leq \lambda\sigma^2$: $\gamma = 0$ satisfies the inclusion. Therefore $\hat{\gamma}_{\text{MAP}} = \mathcal{S}_{\lambda\sigma^2}(y)$.

The unnormalized posterior is $\exp(-(y-\gamma)^2/(2\sigma^2) - \lambda|\gamma|)$. Splitting on $\gamma > 0$ and $\gamma < 0$: $\hat{\gamma}_{\text{MMSE}} = [I_+ - I_-]/Z$ where $I_\pm = \int_0^\infty \gamma \exp(-(y \mp \gamma)^2/(2\sigma^2) - \lambda\gamma)\mathrm{d}\gamma$. For $\sigma=1, \lambda=1$: $\hat{\gamma}(y=0.1) \approx 0.12$, $\hat{\gamma}(y=0.5) \approx 0.21$, $\hat{\gamma}(y=1.0) \approx 0.47$, $\hat{\gamma}(y=3.0) \approx 2.82$.

MAP and MMSE agree for large $|y|$ (both approximately $y - \lambda\sigma^2$ for $y \gg \lambda\sigma^2$). They differ near the threshold: MAP gives exact zeros for $|y| < \lambda\sigma^2$, MMSE gives small nonzero values. For support identification, MAP is preferred (exact sparsity). For prediction/reconstruction accuracy (MSE), MMSE is optimal by construction.

$\boldsymbol{\gamma}$ and $\mathbf{w}$ are independent Gaussians, so $\mathbf{y} = \mathbf{A}\boldsymbol{\gamma} + \mathbf{w}$ has covariance $\alpha^{-1}\mathbf{A}\mathbf{A}^H + \sigma^2\mathbf{I}$. The marginal is $\mathbf{y} \sim \mathcal{N}(\mathbf{0}, \mathbf{C}_\alpha)$ with $\mathbf{C}_\alpha = \sigma^2\mathbf{I} + \alpha^{-1}\mathbf{A}\mathbf{A}^H$.

$\log\mathcal{Z}(\mathbf{y} \mid \alpha) = -\frac{1}{2}[\mathbf{y}^H\mathbf{C}_\alpha^{-1}\mathbf{y} + \log\det\mathbf{C}_\alpha + m\log(2\pi)]$. Using the determinant lemma: $\log\det\mathbf{C}_\alpha = m\log\sigma^2 + \sum_k\log(1 + \sigma_k^2/(\alpha\sigma^2))$.

Differentiating $\log\mathcal{Z}$ w.r.t. $\alpha$ and setting to zero yields the MacKay update: $\hat{\alpha}^{-1} = \|\boldsymbol{\mu}\|^2/(n - \hat{\alpha}\operatorname{tr}(\mathbf{\Sigma}))$ where $n - \hat{\alpha}\operatorname{tr}(\mathbf{\Sigma})$ is the effective number of "well-determined" parameters (the "gamma number" in the original evidence approximation). $\blacksquare$

From Theorem (Gaussian Prior-Posterior): $\mathbf{\Sigma} = (\sigma^{-2}\mathbf{A}^H\mathbf{A} + \mathbf{D}_\alpha)^{-1}$ where $\mathbf{D}_\alpha = \text{diag}(\boldsymbol{\alpha})$, and $\boldsymbol{\mu} = \sigma^{-2}\mathbf{\Sigma}\mathbf{A}^H\mathbf{y}$.

$\mathbb{E}[\log\pi(\boldsymbol{\gamma} \mid \boldsymbol{\alpha})] = \sum_i \frac{1}{2}\log\alpha_i - \frac{\alpha_i}{2}\mathbb{E}[\gamma_i^2] = \sum_i \frac{1}{2}\log\alpha_i - \frac{\alpha_i}{2}(\mu_i^2 + [\mathbf{\Sigma}]_{ii})$. Differentiating w.r.t. $\alpha_i$ and setting to zero: $\frac{1}{2\alpha_i} = \frac{\mu_i^2 + [\mathbf{\Sigma}]_{ii}}{2}$, giving $\alpha_i^{-1} = \mu_i^2 + [\mathbf{\Sigma}]_{ii}$. Rearranging using $[\mathbf{\Sigma}]_{ii} = \alpha_i^{-1}(1 - \alpha_i[\mathbf{\Sigma}]_{ii})$... yields the stated formula.

$\alpha_i^{\text{new}} = (1 - \alpha_i[\mathbf{\Sigma}]_{ii})/\mu_i^2$. Note $\alpha_i[\mathbf{\Sigma}]_{ii} \in (0,1)$ (fraction of posterior variance from prior). If $\mu_i^2 < (1 - \alpha_i[\mathbf{\Sigma}]_{ii})$ then $\alpha_i^{\text{new}} > 1/\mu_i^2 > \alpha_i$ (grows). If additionally $\mu_i^2 \leq \alpha_i[\mathbf{\Sigma}]_{ii}$, then $\alpha_i^{\text{new}} \leq 0$, which is infeasible — the EM update drives $\alpha_i \to \infty$ (prune). $\blacksquare$

With $\gamma \mid \lambda \sim \mathcal{N}(0, \lambda^2\tau^2)$ and $y \mid \gamma \sim \mathcal{N}(\gamma, 1)$, the conjugate formula gives: posterior variance = $(1/(\lambda^2\tau^2) + 1)^{-1} = \lambda^2\tau^2/(1 + \lambda^2\tau^2) = 1 - \kappa$, posterior mean $= (1 - \kappa) \cdot 1 \cdot y = (1-\kappa)y$.

$\lambda \sim \text{Half-Cauchy}(0,1)$ has density $p(\lambda) = 2/(\pi(1+\lambda^2))$ for $\lambda > 0$. Change variables: $\kappa = 1/(1+\lambda^2)$, $\lambda = \sqrt{(1-\kappa)/\kappa}$, $|d\lambda/d\kappa| = 1/(2\kappa^{3/2}(1-\kappa)^{1/2})$. $p(\kappa) \propto \kappa^{-1/2}(1-\kappa)^{-1/2}$ — poles at $\kappa = 0$ and $\kappa = 1$.

$\mathbb{E}[\gamma \mid y] = \int_0^\infty (1-\kappa(\lambda))y \cdot p(\lambda \mid y) \mathrm{d}\lambda$. For $y=0.1$: horseshoe $\approx 0.04$ (aggressive shrinkage), Laplace MAP $= 0$ (threshold $\approx 1$). For $y=4$: horseshoe $\approx 3.98$ (minimal shrinkage), Laplace MAP $= 3$ (still shrunk). $\blacksquare$

The operator $-\partial_{xx}$ on $L^2([0,1])$ with periodic BC has eigenfunctions $\phi_k(x) = e^{2\pi ikx}$ with eigenvalues $(2\pi k)^2$, $k \in \mathbb{Z}$. Therefore $(\kappa^2 I - \partial_{xx})\phi_k = (\kappa^2 + (2\pi k)^2)\phi_k$, so $\lambda_k = (\kappa^2 + (2\pi k)^2)^{-s}$.

$\operatorname{tr}(\mathcal{C}_0) = \sum_{k=-\infty}^\infty (\kappa^2 + (2\pi k)^2)^{-s}$. For large $|k|$, $\lambda_k \sim |k|^{-2s}$. The series converges iff $\sum_k |k|^{-2s} < \infty$, which requires $2s > 1$. Hence trace class iff $s > 1/2$.

$\mathcal{H} = \operatorname{Range}(\mathcal{C}_0^{1/2})$ consists of functions $h = \sum_k h_k \phi_k$ with $\|h\|_{\mathcal{H}}^2 = \sum_k |h_k|^2/\lambda_k = \sum_k |h_k|^2(\kappa^2 + (2\pi k)^2)^s < \infty$. This is the Sobolev space $H^s([0,1])$. The Cameron-Martin norm is the $H^s$ seminorm weighted by $(\kappa^2 + |\cdot|^2)^s$. $\blacksquare$

Matérn-3/2 eigenvalues: $\lambda_k = (\kappa^2 + (2\pi k/n)^2)^{-3/2}$ (discrete Fourier modes). Prior sample: draw $\xi_k \sim \mathcal{N}(0,1)$ for $k = 0, \ldots, 49$, compute $\hat{\gamma}_k = \sqrt{\lambda_k}\xi_k$, pad with zeros for $k \geq 50$, IFFT to get $\boldsymbol{\gamma}$.

At each step: propose $\boldsymbol{\gamma}' = \sqrt{1-\beta^2}\boldsymbol{\gamma} + \beta\boldsymbol{\xi}_{\text{prior}}$. Compute $\Delta\Phi = \Phi(\boldsymbol{\gamma}') - \Phi(\boldsymbol{\gamma})$. Accept if $u < \exp(-\Delta\Phi)$, $u \sim U[0,1]$. Near-optimal $\beta$ gives acceptance rate $\approx 0.234$.

From samples $\{\boldsymbol{\gamma}^{(t)}\}_{t=1001}^{5000}$: $\hat{\boldsymbol{\mu}} = \frac{1}{4000}\sum_t \boldsymbol{\gamma}^{(t)}$, $\hat{\sigma}_i^2 = \frac{1}{4000}\sum_t(\gamma_i^{(t)} - \hat{\mu}_i)^2$. Credible band: $[\hat{\mu}_i - 1.96\hat{\sigma}_i,\, \hat{\mu}_i + 1.96\hat{\sigma}_i]$. With correct $\beta$ and sufficient samples, $\approx 95\%$ of true values should be covered. $\blacksquare$

Proposal: $\boldsymbol{\gamma}' = \boldsymbol{\gamma} + \delta\boldsymbol{\varepsilon}$, $\boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. Acceptance: $\min(1, \exp(-\|\boldsymbol{\gamma}'\|^2/2 + \|\boldsymbol{\gamma}\|^2/2))$. With $\delta = 2.38/\sqrt{n}$: acceptance $\approx 0.234$, ESS $\approx T/n$ — scales as $n^{-1}$.

Introduce momentum $\mathbf{p} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. Leapfrog with $L = 10$ steps, step size $\varepsilon$ tuned for 65% acceptance. ESS $\approx T/n^{1/4}$ — confirmed numerically. Each HMC proposal costs $L$ gradient evaluations.

For $n = 500$: random-walk MH has ESS/eval $\approx 1/500 = 0.002$ while HMC has ESS/eval $\approx 1/4.7 \approx 0.2$. HMC is $\approx 100\times$ more efficient. The tuning difficulty: MH requires only $\delta$; HMC requires $\varepsilon$ and $L$ (mitigated by NUTS which auto-tunes $L$). $\blacksquare$

$\mathbf{\Gamma}_{\text{post}} = (\sigma^{-2}\mathbf{A}^H\mathbf{A} + \gamma_0^{-2}\mathbf{I})^{-1}$. For DFT rows, $\mathbf{A}^H\mathbf{A}$ has a circulant structure exploitable via FFT.

$\min_{S \subseteq [2N],\, |S|=m} \operatorname{tr}\!\left((\sigma^{-2}\sum_{k\in S}\mathbf{a}_k\mathbf{a}_k^H + \gamma_0^{-2}\mathbf{I})^{-1}\right)$. This is NP-hard in general; the greedy algorithm has a $(1-1/e)$ guarantee for the related D-optimal criterion.

Greedy: initialize $\mathbf{B}_0 = \gamma_0^{-2}\mathbf{I}$; at step $j$, add frequency $k^* = \arg\min_k \operatorname{tr}((\mathbf{B}_{j-1} + \sigma^{-2}\mathbf{a}_{k^*}\mathbf{a}_{k^*}^H)^{-1})$ via rank-1 update formula. For $N=64$, $m=16$, greedy achieves $\operatorname{tr}(\mathbf{\Gamma}_{\text{post}}) \approx 0.7\times$ that of random selection — confirming the practical value of optimal design. $\blacksquare$

$\mathbf{A} \in \mathbb{R}^{200 \times 4096}$ with i.i.d. $\mathcal{N}(0, 1/200)$ entries. Place 3 targets at random positions with $|\gamma_i| \sim \text{Uniform}(0.5, 2)$, angle uniform on $[0, 2\pi)$.

Use Woodbury: $(\sigma^{-2}\mathbf{A}^H\mathbf{A} + \mathbf{D}_\alpha)^{-1} = \mathbf{D}_\alpha^{-1} - \mathbf{D}_\alpha^{-1}\mathbf{A}^H(\sigma^2\mathbf{I}_m + \mathbf{A}\mathbf{D}_\alpha^{-1}\mathbf{A}^H)^{-1}\mathbf{A}\mathbf{D}_\alpha^{-1}$. Cost per EM iteration: $O(m^2n + m^3)$ — much cheaper than $O(n^3)$ for $m \ll n$.

Typical results (SNR 10 dB): SBL NMSE $\approx -18$ dB, detection rate $\approx 97\%$ at PFA $\leq 1\%$. LASSO NMSE $\approx -14$ dB, detection rate $\approx 88\%$ at same PFA. SBL outperforms due to automatic regularization and sparse-promoting hyperparameter updates. $\blacksquare$

For a uniformly sampled circular aperture with $m$ Tx-Rx pairs evenly spaced in angle, the PSF (point spread function) $[\mathbf{A}^H\mathbf{A}]_{ij}$ is a function only of $|x_i - x_j|$ (in the 2D sense) — approximately shift-invariant. The posterior variance is approximately constant: $[\mathbf{\Gamma}_{\text{post}}]_{ii} \approx \sigma_{\text{post}}^2$ for all $i$ near the center.

Linear aperture: k-space coverage is limited to a sector, leading to non-uniform resolution. Pixels at scene edges require large-angle components of $\mathbf{A}$, which are undersupported. $[\mathbf{\Gamma}_{\text{post}}]_{ii}$ grows toward the edges, revealing reduced resolution and increased uncertainty. This motivates MIMO configurations ([?ch33:s01]) that provide more uniform k-space coverage. $\blacksquare$

For the correctly specified model, by Bayes optimality the credible intervals are calibrated: empirical coverage should be $\approx 95\%$ (with $\pm 1\%$ Monte Carlo error for 500 trials and 20 pixels $\approx 10000$ coverage events).

With $\boldsymbol{\gamma} \sim \mathcal{N}(\mathbf{0}, 4\mathbf{I})$ but prior $\mathcal{N}(\mathbf{0}, \mathbf{I})$: the true prior variance is $4\times$ larger than assumed. The posterior covariance $\mathbf{\Gamma}_{\text{post}}$ is too small (overconfident), so credible intervals are too narrow. Empirical coverage $\approx 72\%$ instead of $95\%$. $\blacksquare$

$\mathcal{R}(\hat{\boldsymbol{\gamma}}) = \mathbb{E}[\|\boldsymbol{\gamma} - \hat{\boldsymbol{\gamma}}(\mathbf{y})\|^2] = \mathbb{E}_{\mathbf{y}}\bigl[\mathbb{E}[\|\boldsymbol{\gamma} - \hat{\boldsymbol{\gamma}}(\mathbf{y})\|^2 \mid \mathbf{y}]\bigr]$. Minimizing over all estimators $\hat{\boldsymbol{\gamma}}(\cdot)$ is equivalent to minimizing the inner expectation pointwise for each $\mathbf{y}$.

Fix $\mathbf{y}$ and minimize $f(\mathbf{c}) = \mathbb{E}[\|\boldsymbol{\gamma} - \mathbf{c}\|^2 \mid \mathbf{y}]$ over $\mathbf{c} \in \mathbb{R}^n$. $\nabla_{\mathbf{c}} f = -2\mathbb{E}[\boldsymbol{\gamma} - \mathbf{c} \mid \mathbf{y}] = 0$ gives $\mathbf{c} = \mathbb{E}[\boldsymbol{\gamma} \mid \mathbf{y}] = \hat{\boldsymbol{\gamma}}_{\text{MMSE}}$. The Hessian $\nabla^2 f = 2\mathbf{I} \succ 0$, confirming this is a global minimum. $\blacksquare$