# Regularization vs constrained optimization of an ill posed tomography problem

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.

• Have you read this post in Mathematics.SO? – nicoguaro Jun 17 '17 at 15:46
• For an ill-posed, underdetermined problem, the condition number is of course infinite. Your statement of $10^{16}$ is only a numerical representation of this infinity. – Wolfgang Bangerth Jun 17 '17 at 16:30
• I think you will want to read the book by Engl, Hanke, Neubauer: "Regularization of Inverse Problems". – Wolfgang Bangerth Jun 17 '17 at 16:31
• Do you know anything else about the desired solution? Is it nonnegative? Smooth? Sparse? Are there constraints that the solution should satisfy?you need a reason to pick out one solution as best. – Brian Borchers Jun 20 '17 at 2:16
• @nicoguaro: Thanks for sharing this post. It helps in understanding the solution from different approaches, but I am confused about which one would be better. – Sahil Gupta Jun 21 '17 at 10:41

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$$.

### Least-norm

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}$$

### Least-squares

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).

• Thanks for the answer, but this does not answer the original question: I wanted to understand the differences between constrained optimization and regularization. Your approach to this problem is similar to optimization, where you define a Lagrangian and minimize the function. Theoretically, the final solution by you is correct, provided the data noise was absent. In my case, as the condition number of A is very high, it blows-up the error. This is where the regularization of problem helps in obtaining a better solution. – Sahil Gupta Jun 21 '17 at 12:18
• A further note, I have tried constrained optimization and regularization of the problem. Tikhonov Regularization fails for complicated profiles (I have not yet tried other regularization techniques) and constrained optimization takes too long to solve the problem. Before moving further, I wanted to learn more about the differences in these approaches. – Sahil Gupta Jun 21 '17 at 12:21
• @Sahil Gupta write in the preceding comment: "constrained optimization takes too long to solve the problem". What constraints did you use? What is the formulation of your optimization problem? What method did you use to solve it? – Mark L. Stone Jun 21 '17 at 15:12
• @MarkL.Stone I used the fmincon function available in MATLAB link. The additional constraints along with the linear system of equations were lower and upper bounds of 0 and 1. – Sahil Gupta Jun 23 '17 at 5:55
• Thanks for the updated answer; this sums up pretty much the complete regularization approach. Tikhnov Regularization is similar to your least-squares approach; the only difference being the regularization term. I have a doubt regarding this approach: how does one find a suitable value for gamma? The usual approach is the L-Curve corner. But in my case, as I increase the complexity of the image, regularization methods fail to obtain a suitable solution for the corner value. Although the constrained optimization captures a reasonable solution. What could be the reason for this? – Sahil Gupta Jun 23 '17 at 6:03