Inverse Problems and Regularization

Definition Let $A \in \mathbb{F}^{m \times n}$ be a known matrix that models a measurement. A linear inverse problem asks for an estimate of an unknown input $x_{\text{true}} \in \mathbb{F}^{n}$ from observed data

where the noise $e$ is unknown too. Take a full SVD $A = U\Sigma V^{*}$ with $U$ and $V$ unitary, columns $u_1, \dots, u_m$ of $U$ and $v_1, \dots, v_n$ of $V$, rank $r \ge 1$ and positive singular values $\sigma_1 \ge \dots \ge \sigma_r > 0$. Then $Av_j = \sigma_j u_j$ for $j \le r$ and $Av_j = 0$ for $j > r$. Write $\beta_j = u_j^{*}b$ for the coordinates of the data. The baseline reconstruction is the pseudoinverse solution

the least-squares solution of smallest norm. The least-squares solutions are exactly the vectors $x^{+} + z$ with $z \in \ker A$. Noise amplification, the weakness of this baseline, means that a small error in the data can become a much larger error in the reconstruction. The noise passes through linearly, as

so the component of the noise along $u_j$ is multiplied by $1/\sigma_j$ for $j \le r$, and the components along $u_j$ with $j > r$ are dropped. Because the $v_j$ are orthonormal and $\sum_{j} \lvert u_j^{*}e \rvert^2 = \lVert e \rVert^2$, this gives

with equality when $e$ is a multiple of $u_r$. For square invertible $A$ and $x_{\text{true}} \ne 0$, the reconstruction $x^{+} = A^{-1}b$ has error $x^{+} - x_{\text{true}} = A^{-1}e$, and the same bound together with $\lVert Ax_{\text{true}} \rVert \le \sigma_1 \lVert x_{\text{true}} \rVert$ gives the relative bound

Definition In the setting above, choose a threshold $\tau > 0$. The truncated SVD reconstruction keeps only the singular directions with $\sigma_j > \tau$,

If exactly $k$ singular values exceed $\tau$, then $x_{\tau} = A_k^{+}b$ for the rank $k$ truncation $A_k = \sum_{j=1}^{k} \sigma_j u_j v_j^{*}$, so $x_{\tau}$ is the minimum-norm least-squares solution for the truncated matrix. Since $v_j^{*}A^{*}b = \sigma_j \beta_j$, the reconstruction equals $\sum_{\sigma_j > \tau} \sigma_j^{-2}\, v_j v_j^{*} A^{*}b$, and grouping equal singular values turns the sum of the $v_j v_j^{*}$ into orthogonal projectors onto eigenspaces of $A^{*}A$. So $x_{\tau}$ does not depend on which SVD is used, even when singular values repeat. Every retained data coordinate $\beta_j$ is multiplied by $1/\sigma_j < 1/\tau$, and discarded directions contribute nothing. The threshold is a modeling decision tied to the noise level, the units or the resolution needed. Linear algebra alone does not supply it.

Definition Let $A \in \mathbb{F}^{m \times n}$, $b \in \mathbb{F}^{m}$ and $\lambda > 0$. Here $\lambda$ is a penalty parameter, not an eigenvalue. Tikhonov regularization chooses $x \in \mathbb{F}^{n}$ to minimize

The theorem says that $J_{\lambda}$ has exactly one minimizer $x_{\lambda}$, that $x_{\lambda}$ is the unique solution of the regularized normal equations

and that in a full SVD of $A$

For the proof, put $c = V^{*}x$. Unitary matrices preserve norms, so

Each term involves one coefficient. In a term with $j \le r$, put $\sigma = \sigma_j$, $\beta = \beta_j$, $a = \sigma^2 + \lambda > 0$ and $c_0 = \sigma\beta / a$, and complete the square,

Expanding $a \lvert c - c_0 \rvert^2$ gives back the same cross term, and the constants agree because $\lvert \beta \rvert^2 - \sigma^2 \lvert \beta \rvert^2 / a = \lambda \lvert \beta \rvert^2 / a$. So the term is smallest exactly at $c_j = c_0$. A term $\lambda \lvert c_j \rvert^2$ with $j > r$ is smallest exactly at $c_j = 0$, because $\lambda > 0$, and the last sum does not involve $x$. This proves existence, uniqueness and the formula without any calculus. For the linear system, every $z \ne 0$ gives

so $A^{*}A + \lambda I$ is positive definite and invertible even when $A$ is rank deficient. In the basis of the $v_j$ it is diagonal with entries $\sigma_j^2 + \lambda$, where $\sigma_j = 0$ for $j > r$, and $V^{*}A^{*}b$ has entries $\sigma_j \beta_j$ for $j \le r$ and $0$ after that, so its solution has the coordinates found above. At $\lambda = 0$ the coordinates $c_j$ with $j > r$ drop out of the function, so the minimizer is not unique when $r < n$.

Definition All three reconstructions have the form

with numbers $\varphi_j$ between $0$ and $1$, called filter factors. The pseudoinverse uses $\varphi_j = 1$. Truncation uses $\varphi_j = 1$ when $\sigma_j > \tau$ and $\varphi_j = 0$ otherwise. Tikhonov regularization uses

which is near $1$ when $\sigma_j^2$ is much larger than $\lambda$ and near $0$ when it is much smaller. The number $\varphi_j / \sigma_j$ is the gain applied to the data coordinate $\beta_j$. A data error $e$ changes the reconstruction by $\sum_j (\varphi_j/\sigma_j)(u_j^{*}e)\, v_j$, whose norm is at most $\max_j (\varphi_j/\sigma_j) \lVert e \rVert$. The largest pseudoinverse gain $1/\sigma_r$ grows without limit as $\sigma_r$ shrinks toward $0$, the truncated gains stay below $1/\tau$, and every Tikhonov gain satisfies

because $\sigma^2 + \lambda - 2\sigma\sqrt{\lambda} = (\sigma - \sqrt{\lambda})^2 \ge 0$, with equality at $\sigma = \sqrt{\lambda}$. As $\lambda \to 0^{+}$ every Tikhonov factor tends to $1$ and $x_{\lambda} \to x^{+}$. As $\lambda \to \infty$ every factor tends to $0$ and $x_{\lambda} \to 0$.

Definition For a reconstruction with filter factors $\varphi_j$ between $0$ and $1$, orthonormality of the $u_j$ and of the $v_j$ gives

Raising a filter factor $\varphi_j$ with $\beta_j \ne 0$ lowers the residual and raises the solution norm. This is the residual and solution norm tradeoff. For every input $x$ whatsoever, the squared residual is at least the sum of $\lvert \beta_j \rvert^2$ over $j > r$, which does not depend on $x$, so the pseudoinverse, with all $\varphi_j = 1$, has the smallest residual possible. In Tikhonov regularization every $\varphi_j$ strictly decreases as $\lambda$ grows, so $\lVert x_{\lambda} \rVert$ decreases and $\lVert b - Ax_{\lambda} \rVert$ increases, strictly when $\beta_j \ne 0$ for some $j \le r$. No input does better than $x_{\lambda}$ on both counts. If $x \ne x_{\lambda}$ and $\lVert x \rVert \le \lVert x_{\lambda} \rVert$, then $J_{\lambda}(x) > J_{\lambda}(x_{\lambda})$ gives

Definition For $A \in \mathbb{F}^{m \times n}$, the unrecoverable component of an input $x$ is its orthogonal projection onto $\ker A$, which is $(I - A^{+}A)x$. Because $A(x + z) = Ax$ for every $z \in \ker A$, inputs that differ only in this component produce identical noiseless data, so when $\ker A \ne \{0\}$ no reconstruction computed from $A$ and $b$ alone can determine it. All three reconstructions lie in $(\ker A)^{\perp} = \operatorname{span}\{v_1, \dots, v_r\}$. They set the unrecoverable component to zero by a convention or a penalty, not by a measurement. For the pseudoinverse the error splits into amplified noise and the lost kernel part,

A different loss lives in the data space. The component $(I - AA^{+})b$ of the data in $(\operatorname{im} A)^{\perp} = \ker A^{*}$ is orthogonal to every output $Ax$, so it stays in every residual and $\lVert b - Ax \rVert \ge \lVert (I - AA^{+})b \rVert$ for every $x$. The kernel $\ker A$ has dimension $n - r$ and lies in the input space, while $\ker A^{*}$ has dimension $m - r$ and lies in the output space. The first consists of input directions that the measurement never records, and the second of data directions that no input can explain.