What is this regularization technique?

Question

I'm currently working on understanding some very old code. They are trying to solve an underdetermined system and the comments say that they want the minimum norm solution.

The system they solve is: $A x = b$, with $A \in \mathbb{R}^{5 \times n} (n > 5)$

Unless I'm mistaken, you can obtain a minimum norm solution via an LQ factorization, for example using the LAPACK function dgels.

The code instead constructs the matrix

$B = \begin{bmatrix} 2 & A^T \\ A & 0 \end{bmatrix}$

and then solves:

$B \begin{bmatrix} x \\ \tau \end{bmatrix} = \begin{bmatrix} 0 \\ b \end{bmatrix}$

This appears to result in the same minimum norm solution as with an LQ. What is this technique?

Answer 1

All they are doing is taking the problem of minimizing $x^T x$ subject to $A x = b$ and then forming a Lagrangian, like so: $$L(x, \lambda) = x^T x + \lambda^T (A x - b)$$

From here, you want to find optimal values $x^*$ and $\lambda^*$ that make the gradient of the Lagrangian the zero vector. This leads to the following condition that must hold for the optimal values: $$\begin{align*} 2 x^* + A^T \lambda^* &= 0 \\ A x^* &= b \end{align*}$$

which gives you the system of equations you see.

Answer 2

I have never seen it used for underdetermined systems, but for overdetermined systems (tall thin $A$) the version obtained replacing $A$ with $A^T$ in your block matrix a common reformulation, known as "augmented system"; Golub and Van Loan (4th ed, sec. 5.3.8) for instance suggest using it for iterative improvement.

It's not a great idea to use it as a solver instead of a QR / LQ factorization: the main drawback is that $B$ is very large (the first block alone is $n \times n$), so you're replacing a problem that can be solved in $O(n)$ with one that can be solved in $O(n^3)$ (where $n$ is the largest dimension).

Answer 3

Spektr's answer is already great but note that you can also argue about this in the following manner $2x + A^T \tau = 0$ and $Ax=b$ implies that $x = -\frac{1}{2} A^T\tau$ and plugging it in the second equation you get $AA^T \tau = -2b$ and thus $x = A^T(AA^T)^{-1}b$ which is a right inverse for $A$ provided of course that $A$ has linearly independent rows (this is the necessary and sufficient condition for $(AA^T)^{-1}$ to exist). If $A$ doesn't have linearly independent rows then it can happen that $b$ is not in the range of $A$, in which case this has no solution. If $b$ is in the range of $A$ then this will still work but it will be equivalent to using a generalized equation-solving inverse $(AA^T)^{-}$ of $AA^T$. The matrix $A^T(AA^T)^{-}$ would then necessarily be a $\{1,4\}$ inverse, i.e. an equation solving inverse minimizing the norm of the solution in the Euclidean norm. This means that the above solves the problem $\min_{Ax=b} \|x\|^2$. If you are interested in iterative algorithms doing this, then look up the CGNE/Craig/CGLS algorithm.

Note that the above will fail if $b$ is not in the range of $A$. To rectify this one often considers $\min_{x} \|Ax-b\|^2$. If $b$ is not in the range of $A$ this will give you a $x$ such that $Ax$ is the projection of $b$ on the span of the columns of $A$. The minimizers solve the normal equations $A^TAx = A^Tb$, or equivalently: $$\begin{bmatrix} -I & A\\ A^T & 0 \end{bmatrix} \begin{bmatrix} y \\ x \end{bmatrix} = \begin{bmatrix} 0 \\ A^Tb\end{bmatrix}.$$ You can of course try to filter again the solution with minimal norm $\min_{x \in S} \|x\|^2, \, S = \arg\min_x \|Ax-b\|^2$. If you initialize the CGNR algorithm with an intial guess of zero you get this solution, it is the Moore-Penrose pseudoinverse solution.