Linear least square method for singular matrices

[alyflex]

I have stumbled upon a problem which I have so far been unable to solve.

I we consider a general set of linear equations:
Ax=b,

I know the the system is inconsistent which makes least square method the logical choice.

So the mission is to minimize ||Ax-b||

And the usual way I do this is by setting Ax=p, where p is the projection of b onto Ax.
By isolating:
x=(A^TA)^{-1}A^T \cdot p<br /> =(A^TA)^{-1} \cdot A^T \cdot b

However the product A^T*A is also singular, and thus I am unable to do this.

I pretty sure there is a very simple way to do this, but when i look in my old algebra book I see no solution to the problem.

Anyone know the way?

 

[HallsofIvy]

Yes, that is the standard "least squares" method for solving such a problem with (A^TA)^{-1}A^T being the "generalized inverse". For most such problems A^TA is invertible. If your A is such that A^TA is not invertible then you have a very pathological problem for which there probably is no "simple" way to solve it. What is A?

 

[Wizlem]

Isn't it true that A^TA is only not invertible if there are actually multiple solutions to Ax=b?

 

[Landau]

Isn't it true that A^TA is only not invertible if there are actually multiple solutions to Ax=b?

Correct. The least square solution satisfies the normal equation A^TAx=A^Ty. This solution is unique if and only if A^TA is invertible. See here for a similar discussion.

 

[hotvette]

I can think of a couple of possibilities, both involving miniminum norm solutions. If A has more columns than rows, it is an underdetermined system that has an easy 'least norm' solution. If A has more rows than columns and is rank deficient (my guess in your case), there is a way to minimize the norm of Ax-b and the norm of x at the same time using SVD, but I'm not familiar with the details. Look up LAPACK DGELSS.

* DGELSS computes the minimum norm solution to a real linear least
* squares problem:
*
* Minimize 2-norm(| b - A*x |).
*
* using the singular value decomposition (SVD) of A. A is an M-by-N
* matrix which may be rank-deficient.

 

[alyflex]

Correct. The least square solution satisfies the normal equation A^TAx=A^Ty. This solution is unique if and only if A^TA is invertible. See here for a similar discussion.

This was what I suspected, but thanks for making it clear. Now I just need a method to compute it for non-invertible matrices

I can think of a couple of possibilities, both involving miniminum norm solutions. If A has more columns than rows, it is an underdetermined system that has an easy 'least norm' solution. If A has more rows than columns and is rank deficient (my guess in your case), there is a way to minimize the norm of Ax-b and the norm of x at the same time using SVD, but I'm not familiar with the details. Look up LAPACK DGELSS.

The more I look at this the more afraid I am, that a numerical solution might be the best way. I was really hoping for a clear and illuminating analytical method.

What is A?

The initial system I stumbled upon was this:

\left[ {\begin{array}{ccc}<br /> 1 &amp; 1 &amp; 3 \\<br /> -1 &amp; 3 &amp; 1 \\<br /> 1 &amp; 2 &amp; 4 \\<br /> \end{array} } \right]<br /> \left( {\begin{array}{c}<br /> x \\<br /> y \\<br /> z \\<br /> \end{array} } \right) <br /> = \left( {\begin{array}{c}<br /> -2 \\<br /> 0 \\<br /> 8 \\<br /> \end{array} } \right)

but then I considered more simple system but with the same problem.
One example could be:

y=x
y=x+1

It is clear what the least square solution to this system should be, but I have no idea how to find it methodically

 

[trambolin]

Consider A being column eliminated with R as

<br /> A R R^{-1}x = b<br />

Then,
<br /> \left[ {\begin{array}{ccc}<br /> 1 &amp; 0 &amp; 0 \\<br /> -1 &amp; 4 &amp; 0 \\<br /> 1 &amp; 1 &amp; 0 \\<br /> \end{array} } \right] \left( {\begin{array}{c}<br /> x + y +3z \\<br /> y +z \\<br /> z \\<br /> \end{array} } \right) =\left( {\begin{array}{c}<br /> -2 \\<br /> 0 \\<br /> 8 \\<br /> \end{array} } \right)<br />
now we solve a different linear eq. set, since the last variable does not play a role in the equations we reduce the number of columns.
<br /> \underbrace{\left[ {\begin{array}{cc}<br /> 1 &amp; 0 \\<br /> -1 &amp; 4 \\<br /> 1 &amp; 1 \\<br /> \end{array} } \right]}_{A_r} \left( {\begin{array}{c}<br /> x + y +3z \\<br /> y +z<br /> \end{array} } \right) = \left( {\begin{array}{c}<br /> x + y +3z \\<br /> -x +3y + z \\<br /> x+2y+4z<br /> \end{array} } \right) =\underbrace{\left( {\begin{array}{c}<br /> -2 \\<br /> 0 \\<br /> 8 \\<br /> \end{array} } \right)}_{b}<br />

From the second equation, you can see that it is not solvable. So back your least squares, and call the variables a b i.e.

<br /> \left( {\begin{array}{c}<br /> x + y +3z \\<br /> y +z<br /> \end{array} } \right) = \begin{pmatrix}a\\b\end{pmatrix}<br />
And Matlab gave (from linsolve(Ar,b))

<br /> \begin{pmatrix}a\\b\end{pmatrix} = \begin{pmatrix}3\\1\end{pmatrix}<br />

And from this solution, you can construct admissible least squares solutions. I might have made mistakes so please check the steps but conceptually it should be OK if I didn't mess it up. To your small example, the least squares solution is a = y-x = 0.5 So the whole trick is to embed the underdetermined part inside the x vector and solve the least squares solution. Then you get infinitely many solutions that satisfy the least squares solution.EDIT: By the way, I forgot to write that any full rank decomposition of A = MN^T leads to a similar solution. Similarly with SVD mentioned above.

 

[alyflex]

Thank you very much for the illuminating answer trambolin.

I'm going to work it through in the weekend using your method.

I'll post a longer and more conclusive reply once I have tried the method.

But to all of you, thank you very much for an enlightening discussion on a subject, which was causing me severe headaches yesterday. =)

 

[alyflex]

Well I just went through the calculation myself and as far as I can see everything you wrote was spot on, so thank you very much.

Btw I had no idea that the linsolve function in MATLAB solved a system of linear equations using least squares, I only thought it could find consistent solutions.

 

[mattrb]

alyflex,

Try the truncated SVD of the matrix A (Demmel, Applied Numerical Linear Algebra). Let A=U\Sigma V^T where A,U\in \Re^{m\times n}, \Sigma,V\in \Re^{n \times n}. See en.wikipedia.org/wiki/Singular_value_decomposition for details. Every good linear algebra software should have an algorithm. Let \hat{U}\in \Re^{m\times k},\ \hat{\Sigma} \in \Re^{k \times k},\ \hat{V}\in \Re^{n \times k} where k is the numerical rank of A (ignore smallest singular information).

Your least squares solution is x=A^+b=\hat{V}\hat{\Sigma}^{-1}\hat{U}^Tb.
This should work if your matrix is decently sized (around 1000 unknowns).
Another plus is that of all the residual minimizers, this x has a minimal 2-norm.

Matthew