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 =(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} 1 & 1 & 3 \\ -1 & 3 & 1 \\ 1 & 2 & 4 \\ \end{array} } \right] \left( {\begin{array}{c} x \\ y \\ z \\ \end{array} } \right) = \left( {\begin{array}{c} -2 \\ 0 \\ 8 \\ \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

$$A R R^{-1}x = b$$

Then,
$$\left[ {\begin{array}{ccc} 1 & 0 & 0 \\ -1 & 4 & 0 \\ 1 & 1 & 0 \\ \end{array} } \right] \left( {\begin{array}{c} x + y +3z \\ y +z \\ z \\ \end{array} } \right) =\left( {\begin{array}{c} -2 \\ 0 \\ 8 \\ \end{array} } \right)$$
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.
$$\underbrace{\left[ {\begin{array}{cc} 1 & 0 \\ -1 & 4 \\ 1 & 1 \\ \end{array} } \right]}_{A_r} \left( {\begin{array}{c} x + y +3z \\ y +z \end{array} } \right) = \left( {\begin{array}{c} x + y +3z \\ -x +3y + z \\ x+2y+4z \end{array} } \right) =\underbrace{\left( {\begin{array}{c} -2 \\ 0 \\ 8 \\ \end{array} } \right)}_{b}$$

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.

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

$$\begin{pmatrix}a\\b\end{pmatrix} = \begin{pmatrix}3\\1\end{pmatrix}$$

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