I am trying to solve an ill-posed linear system of equations. The particular system has 160 equations and 400 variables. Moreover, the condition number of the left hand side matrix is of order $10^{16}$.
I came across two methods to solve this problem: constrained optimization and regularization methods. My question is what are the pros and cons of each method; which one should I prefer? Any website or paper explaining the same will be much appreciated.
We have the following linear system in $\mathrm x \in \mathbb R^n$
$$\rm A x = b$$
where $\mathrm A \in \mathbb R^{m \times n}$ is fat (i.e., $n > m$) and $\mathrm b \in \mathbb R^m$.
If the linear system is consistent, we look for the least-norm solution via the following (convex) quadratic program
$$\begin{array}{ll} \text{minimize} & \| \mathrm x \|_2^2\\ \text{subject to} & \mathrm A \mathrm x = \mathrm b\end{array}$$
Let the Lagrangian be
$$\mathcal L (\mathrm x, \lambda) := \frac 12 \mathrm x^{\top} \mathrm x + \lambda^{\top} (\mathrm A \mathrm x - \mathrm b)$$
Taking the partial derivatives of $\mathcal L$ and finding where they vanish, we obtain the linear system
$$\begin{bmatrix} \mathrm I_n & \mathrm A^\top\\ \mathrm A & \mathrm O_m\end{bmatrix} \begin{bmatrix} \mathrm x\\ \lambda \end{bmatrix} = \begin{bmatrix} 0_n\\ \mathrm b\end{bmatrix}$$
If $\mathrm A$ has full row rank, then $\rm A A^\top$ is invertible and we can conclude that the least-norm solution is
$$\mathrm x_{\text{LN}} := \color{blue}{\mathrm A^{\top} \left( \mathrm A \mathrm A^{\top} \right)^{-1} \mathrm b}$$
If the linear system is inconsistent, we can look for the least-squares solution via the following unconstrained (convex) quadratic program
$$\begin{array}{ll} \text{minimize} & \| \mathrm A \mathrm x - \mathrm b \|_2^2\end{array}$$
Taking the gradient of the objective function and finding where it vanishes, we obtain the normal equations $\mathrm A^{\top} \mathrm A \,\mathrm x = \mathrm A^{\top} \mathrm b$. However, since $\rm A$ is fat, its rank is at most $m$ and, thus,
$$\mbox{rank} (\mathrm A^{\top} \mathrm A) \leq m < n$$
Hence, $\mathrm A^{\top} \mathrm A$ is never invertible and, thus, the normal equations have infinitely many solutions. Thus, let us add a regularization term to the objective function, i.e.,
$$\begin{array}{ll} \text{minimize} & \| \mathrm A \mathrm x - \mathrm b \|_2^2 + \gamma \| \mathrm x \|_2^2\end{array}$$
where $\gamma \geq 0$. The new normal equations are
$$\left( \mathrm A^{\top} \mathrm A + \gamma \mathrm I_n \right) \mathrm x = \mathrm A^{\top} \mathrm b$$
If $\color{blue}{\gamma > 0}$, then $\mathrm A^{\top} \mathrm A + \gamma \mathrm I_n$ is positive definite and, thus, invertible. Hence, the regularized least-squares solution is
$$\mathrm x_{\text{LS}} := \color{blue}{\left( \mathrm A^{\top} \mathrm A + \gamma \mathrm I_n \right)^{-1} \mathrm A^{\top} \mathrm b}$$
Using a matrix inversion lemma, the regularized least-squares solution can be rewritten as follows
$$\mathrm x_{\text{LS}} := \color{blue}{\mathrm A^{\top} \left( \mathrm A \mathrm A^{\top} + \gamma \mathrm I_m \right)^{-1} \mathrm b}$$
which resembles the least-norm solution. However, we now invert $\mathrm A \mathrm A^{\top} + \gamma \mathrm I_m$, which is positive definite whenever $\gamma > 0$, rather than $\mathrm A \mathrm A^{\top}$ (which may be ill-conditioned or even non-invertible).
