Exercises

$\mathbf{B} \otimes \mathbf{C} = \begin{bmatrix} 0 & 1 & 0 & 2 \\ 1 & 0 & 2 & 0 \\ 0 & 3 & 0 & 4 \\ 3 & 0 & 4 & 0 \end{bmatrix}.$$

$\mathbf{C} \otimes \mathbf{B} = \begin{bmatrix} 0 & 0 & 1 & 2 \\ 0 & 0 & 3 & 4 \\ 1 & 2 & 0 & 0 \\ 3 & 4 & 0 & 0 \end{bmatrix}.$$

These are different matrices (the blocks are in different positions), confirming non-commutativity.

Eigenvalues of $\mathbf{B}$: $\lambda_1 \approx 5.37$, $\lambda_2 \approx -0.37$. Eigenvalues of $\mathbf{C}$: $\mu_1 = 1$, $\mu_2 = -1$. Products: $\{5.37, -5.37, -0.37, 0.37\}$. Both $\mathbf{B} \otimes \mathbf{C}$ and $\mathbf{C} \otimes \mathbf{B}$ have the same eigenvalues (they are similar). $\square$

$\mathbf{B} \otimes \mathbf{C} = \begin{bmatrix} 1 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 2 & 2 \\ 0 & 0 & 0 & 2 \end{bmatrix}$, $\text{vec}(\mathbf{X}) = [1, 3, 2, 4]^T$. Product: $[4, 3, 12, 8]^T$.

$\mathbf{C}\mathbf{X} = \begin{bmatrix} 4 & 6 \\ 3 & 4 \end{bmatrix}$, $\mathbf{C}\mathbf{X}\mathbf{B}^T = \begin{bmatrix} 4 & 12 \\ 3 & 8 \end{bmatrix}$. $\text{vec}(\mathbf{C}\mathbf{X}\mathbf{B}^T) = [4, 3, 12, 8]^T$. They match. $\square$

$M = 8 \times 16 \times 32 = 4096$. $N = 24^3 = 13824$.

$4096 \times 13824 \times 16$ bytes $= 905$ MB.

$16 \times 24 + 8 \times 24 + 32 \times 24 = 384 + 192 + 768 = 1344$ entries. $1344 \times 16$ bytes $= 21.5$ KB.

$905 \text{ MB} / 21.5 \text{ KB} \approx 43{,}000\times$ reduction. $\square$

$[\mathbf{B} \otimes \mathbf{C}]_{(i,k),(j,l)} = b_{ij} c_{kl}$. The conjugate transpose entry is: $[(\mathbf{B} \otimes \mathbf{C})^H]_{(j,l),(i,k)} = \overline{b_{ij} c_{kl}} = \bar{b}_{ij} \bar{c}_{kl} = [\mathbf{B}^H]_{ji} [\mathbf{C}^H]_{lk} = [\mathbf{B}^H \otimes \mathbf{C}^H]_{(j,l),(i,k)}$. $\square$

$\kappa(\mathbf{A}_{1} \otimes \mathbf{A}_{2}) = 5 \times 12 = 60$.

$\kappa = 5 \times 12 \times 3 = 180$.

Reducing $\kappa(\mathbf{A}_{2})$ from 12 to 6 would reduce overall $\kappa$ from 180 to 90 (saving 90). Reducing $\kappa(\mathbf{A}_{1})$ from 5 to 2.5 would give $\kappa = 90$ (also saving 90). Reducing $\kappa(\mathbf{A}_{3})$ from 3 to 1.5 gives $\kappa = 90$ (saving 90).

In terms of relative improvement, halving any factor halves the total. The key insight is that the multiplicative structure means all factors contribute proportionally --- improving the worst-conditioned factor by a given multiplicative ratio has the same effect as improving any other. $\square$

$\mathbf{A}_{1}^{H}\mathbf{A}_{1} = \begin{bmatrix} 4 & 2 \\ 2 & 2 \end{bmatrix}$. Eigenvalues: $3 \pm \sqrt{5}$, so $\sigma_1 = \sqrt{3+\sqrt{5}} \approx 2.29$ and $\sigma_2 = \sqrt{3-\sqrt{5}} \approx 0.874$.

$\mathbf{A}_{2}$ is diagonal: $\sigma_1 = 3$, $\sigma_2 = 1$.

All pairwise products: $\{2.29 \times 3, 2.29 \times 1, 0.874 \times 3, 0.874 \times 1\}$ $= \{6.87, 2.29, 2.62, 0.874\}$. Sorted: $6.87, 2.62, 2.29, 0.874$. $\kappa = 6.87/0.874 \approx 7.86 = (2.29/0.874) \times (3/1) = 2.62 \times 3 \approx 7.86$. $\square$

$[\mathbf{A}]_{n,q} = \frac{1}{\sqrt{N}} e^{-j\pi n u_q}$ for $n = 0,\ldots,15$ and $u_q = -1 + (2q-1)/32$.

$|\langle \mathbf{A}_{p}, \mathbf{A}_{p+1}\rangle| = \frac{1}{N}\left|\sum_{n=0}^{15} e^{j\pi n/16}\right| = \frac{1}{16}\left|\frac{\sin(16\pi/32)}{\sin(\pi/32)}\right| = \frac{|\sin(\pi/2)|}{16\sin(\pi/32)} = \frac{1}{16 \times 0.098} \approx 0.637$.

$\mu_{\text{Welch}} = \sqrt{(32-16)/(16 \times 31)} = \sqrt{16/496} \approx 0.180$.

The physical matrix coherence (0.64) greatly exceeds the Welch bound (0.18). $\square$

$k \approx \frac{1}{2}\sqrt{10000}\ln(2 \times 10^4) = 50 \times 9.9 \approx 495$ iterations.

$k \approx \frac{1}{2}\sqrt{25}\ln(2 \times 10^4) = 2.5 \times 9.9 \approx 25$ iterations.

Unpreconditioned total cost: $495 \times C$. Preconditioned total cost: $25 \times 3C = 75C$. Speedup: $495/75 = 6.6\times$. Preconditioning is highly worthwhile. $\square$

$\boldsymbol{\mathcal{Y}} \in \mathbb{R}^{3 \times 2 \times 5}$.

The mode-2 unfolding $\mathbf{X}_{(2)} \in \mathbb{R}^{4 \times 15}$ stacks the fibers along mode 2. $\mathbf{Y}_{(2)} = \mathbf{M}\mathbf{X}_{(2)} \in \mathbb{R}^{2 \times 15}$. Refolding $\mathbf{Y}_{(2)}$ gives $\boldsymbol{\mathcal{Y}}$.

$2 \times 4 \times 15 = 120$ multiplications. Compare to the full Kronecker product which would involve a matrix of size $(m_1 \times 2 \times m_3) \times (3 \times 4 \times 5)$. $\square$

$\mathbf{G} = (\mathbf{A}_{\text{ang}}^H\mathbf{A}_{\text{ang}}) \otimes (\mathbf{A}_{f}^{H}\mathbf{A}_{f})$. Each factor Gram matrix is Hermitian positive semi-definite.

For a DFT-like matrix with columns $[\mathbf{A}]_{n,q} = e^{-j2\pi nq/N}/\sqrt{M}$, $[\mathbf{G}]_{qq} = \sum_n |[\mathbf{A}]_{nq}|^2 = M/M = 1$ (normalized columns). Since both factors have unit-norm columns, $[\mathbf{G}]_{(q_1,q_2),(q_1,q_2)} = 1$ for all voxels.

Angular PSF: width $\approx 2/N_{\text{virt}} = 2/16 = 0.125$ in $u = \sin\theta$. Range PSF: width $\approx 1/N_f = 1/8$ in normalized delay. The 2D PSF is the outer product of these 1D PSFs. $\square$

$\delta_{s_1 s_2} \leq 0.1 + 0.15 + 0.1 \times 0.15 = 0.265$.

$0.265 < 0.414$, so the RIP condition for LASSO recovery is satisfied.

With equal constants $\delta$: $2\delta + \delta^2 < \sqrt{2} - 1$. $\delta^2 + 2\delta - 0.414 < 0$. Solving: $\delta < \frac{-2 + \sqrt{4 + 1.656}}{2} = \frac{-2 + 2.378}{2} \approx 0.189$. Each factor can have RIP constant up to 0.189. $\square$

$\lambda = 3 \times 10^8 / (28 \times 10^9) = 10.71$ mm. $D = 7 \times 10.71/2 = 37.5$ mm. $R_F = 2 \times (0.0375)^2 / 0.01071 = 0.263$ m.

$R = 5$ m $\gg R_F = 0.26$ m. The far-field assumption is valid for this small array. For a larger array ($N_t = 64$, $D = 340$ mm), $R_F = 21.6$ m, and the near-field correction becomes essential.

The exact distance from element $n$ to a target at $(r, \theta)$: $d_n = \sqrt{r^2 - 2rn d\sin\theta + n^2 d^2} \approx r - nd\sin\theta + \frac{n^2 d^2\cos^2\theta}{2r}$. The quadratic term $n^2 d^2\cos^2\theta/(2r)$ breaks the linear phase assumption and hence the Kronecker factorization.

Include the quadratic phase as a diagonal correction: $\mathbf{A}_{\text{exact}} = \mathbf{A}_{\text{Kron}} \cdot \text{diag}(e^{j\phi_{\text{quad}}})$. This preserves the computational structure (Kronecker + diagonal) while correcting the near-field phase. $\square$

$\mathbf{P}_{\text{diag}} = \text{diag}(\mathbf{A}_{2}^{H}\mathbf{A}_{2}) \otimes \text{diag}(\mathbf{A}_{1}^{H}\mathbf{A}_{1})$. For DFT-like matrices with equal column norms, $\mathbf{P}_{\text{diag}} = c\mathbf{I}$ and the diagonal preconditioner does nothing.

Full singular values: $\{12, 4, 6, 2, 3, 1\}$ (all products). $\kappa = 12/1 = 12$. For unequal column norms (non-DFT case), diagonal preconditioning corrects the norm variation but not the off-diagonal coherence, so $\kappa_{\text{eff}} \leq \kappa$ but improvement is limited.

$\mathbf{P}_{\text{Kron}} = (\mathbf{A}_{2}^{H}\mathbf{A}_{2}) \otimes (\mathbf{A}_{1}^{H}\mathbf{A}_{1})$. $\mathbf{P}_{\text{Kron}}^{-1}\mathbf{A}^{H}\mathbf{A} = \mathbf{I}$ (perfect preconditioning when factors are invertible). The Kronecker preconditioner is strictly better than diagonal. $\square$

Using the factored eigendecomposition, the eigenvalues of $(\mathbf{A}^{H}\mathbf{A} + \tau^{-1}\mathbf{I})^{-1}\mathbf{A}^{H}\mathbf{A}$ are $\frac{\lambda_i\mu_j\nu_l}{\lambda_i\mu_j\nu_l + \tau^{-1}}$.

The trace is $\sum_{i,j,l} \frac{\lambda_i\mu_j\nu_l}{\lambda_i\mu_j\nu_l + \tau^{-1}}$.

This triple sum has $n_1 n_2 n_3 = N$ terms. But if we precompute the per-factor eigenvalues ($n_1 + n_2 + n_3$ values), the sum is $O(N)$ in the worst case.

The Hutchinson estimator uses $R$ random probe vectors: $\text{tr}(\mathbf{M}) \approx \frac{1}{R}\sum_{r=1}^R \mathbf{z}_r^H\mathbf{M}\mathbf{z}_r$. Each probe requires one matvec ($O(N^{4/3})$ with Kronecker) plus one LMMSE solve ($O(N^{4/3})$). Total: $O(RN^{4/3})$.

The exact trace is $O(N)$, which is asymptotically faster than Hutchinson's $O(RN^{4/3})$. The exact computation is preferable whenever the Kronecker eigendecomposition is available (precomputed once per geometry). Hutchinson is useful when the matrix does not admit exact Kronecker factorization. $\square$

$\kappa_{k} = 2\pi(f_0 + f_k)/c$ for each subcarrier. The spatial frequency sampled is $\mathbf{k}_{i,j,k} = \kappa_{k}(\hat{\mathbf{r}}_j - \hat{\mathbf{s}}_i)$. With 3 Tx and 4 Rx, there are $3 \times 4 = 12$ Tx-Rx pairs, each contributing 32 frequency points: 384 total k-space samples.

The radial (range) extent in k-space is $\Delta k_r = 2\piW/c = 2\pi \times 200\times 10^6/3\times 10^8 \approx 4.19$ rad/m. The angular extent depends on the array placement geometry.

If the arrays are clustered on one side of the target (limited angular diversity), the angular k-space coverage is narrow, producing a large condition number for $\mathbf{A}_{\text{ang}}$. The frequency factor is a well-conditioned partial DFT. The overall conditioning is dominated by the angular factor. $\square$

Define $[\mathcal{R}(\mathbf{A})]_{(i-1)n_1+j, (k-1)n_2+l} = [\mathbf{A}]_{(i-1)m_2+k, (j-1)n_2+l}$. This rearranges the $m_1 m_2 \times n_1 n_2$ matrix $\mathbf{A}$ into an $m_1 n_1 \times m_2 n_2$ matrix whose rank-1 SVD gives the NKP.

Let $\mathcal{R}(\mathbf{A}) = \sum_r \sigma_r \mathbf{u}_r \mathbf{v}_r^H$. The best rank-1 approximation gives $\mathbf{B}^* = \sigma_1 \text{reshape}(\mathbf{u}_1, m_1 \times n_1)$, $\mathbf{C}^* = \text{reshape}(\mathbf{v}_1, m_2 \times n_2)$.

$\mathbf{A} \approx \sum_{r=1}^R \mathbf{B}_r \otimes \mathbf{C}_r$. Each matvec requires $R$ Kronecker matvecs: cost $O(R \cdot N^{4/3})$. For $R = 3$: $3\times$ overhead compared to exact Kronecker, but still $\sim 170\times$ faster than the full matvec. $\square$

Let $\mathbf{G} = \mathbf{A}^{H}\mathbf{A}$ and $\hat{\mathbf{G}} = \hat{\mathbf{P}}$. Then $\hat{\mathbf{P}}^{-1}\mathbf{G} = \mathbf{I} + \hat{\mathbf{P}}^{-1}(\mathbf{G} - \hat{\mathbf{G}})$. By Weyl's inequality: $\kappa(\hat{\mathbf{P}}^{-1}\mathbf{G}) \leq \frac{1 + \|\hat{\mathbf{P}}^{-1}\|\|\mathbf{G} - \hat{\mathbf{G}}\|}{1 - \|\hat{\mathbf{P}}^{-1}\|\|\mathbf{G} - \hat{\mathbf{G}}\|}$ which depends on $\epsilon$ and the regularization $\alpha_k$.

The NKP preconditioner improves over unpreconditioned CG when $\kappa(\hat{\mathbf{P}}^{-1}\mathbf{G}) < \kappa(\mathbf{G})$. This holds for $\epsilon \lesssim 1/\sqrt{\kappa(\mathbf{G})}$, i.e., when the NKP approximation error is smaller than the spectral spread of the Gram matrix. $\square$

Classical OAMP: $T_{\text{class}} \sim \sqrt{\kappa}$ iterations. Learned OAMP: $T_{\text{learn}} \sim 5$--$10$ layers (fixed). Per-image cost ratio: $T_{\text{class}}/T_{\text{learn}} \sim \sqrt{\kappa}/5$. For $\kappa = 100$: learned is $2\times$ faster. For $\kappa = 10{,}000$: learned is $20\times$ faster. The advantage grows with the condition number, which is precisely the regime where learned methods shine. $\square$