MikeyW said:
I think I'm understanding it slightly. My problem is mostly concerned with numerical solutions to under-determined problems, and I think I do have the issue that if A is rank-deficient, then rounding errors would appear to make A'*A "singular to working precision" as MATLAB phrases it.
If you are solving an underdetermined linear problem Ax=b, then by definition your A matrix is rank deficient. If rank(A)=r then you only have r nonzero singular values. Now, as indicated by D H, when you use a numerical method to determine the singular values, you may not find any of them that are exactly zero due to the finite precision arrithmatic that is used in the solution - instead you will find r "large" ones, and the rest will be "small."
I'm assuming you are doing rank-deficient least-squares? If so then the SVD is one useful tool. I think of it this way, if U^\dagger A V = \Sigma is the SVD of a matrix A with rank(A)=r, then you can write
<br />
A = \sum_{i=1}^{r} \sigma_i u_i v_i^\dagger,<br />
where \sigma_i is the i^{th} singular value, u_i is the i^{th} left singular vector, andv_i is the i^{th} right singular vector. Since your rank deficient problem, you can add any vector in the nullspace of A to any least-squares solution x_{LS} and maintain the same residual - that is, the rank-deficient problem has an infinite number of solutions. Sometimes you want the solution that has the smallest 2-norm. In these cases the SVD can used to give it to you and the resulting solution is:
<br />
x_{LS} = \sum_{i=1}^{r} \sigma_i^{-1} \left( u_i^\dagger b\right) v_i.<br />
This is related to the so-called pseudo-inverse (
http://en.wikipedia.org/wiki/Moore–Penrose_pseudoinverse). Note that other "complete orthogonal factorizations" can also give you the minimum 2-norm solution. I wrote my own QR based code that does this (see Golub and Van-Loan).
Since you have access to matlab, here is some simple code that solves an underdetermined system via least-squares:
A = randn(15,20);
b = ones(15,1);
x_mldivide = A\b;
x_pinv = pinv(A)*b;
you will find that x_mldivide has 5 zeros, as it provides the solution with the fewest number of non-zero elements (often this is called a "basic" solution), and that x_pinv is actually the minimum norm solution that the svd would give you. A more explicit computation is,
x_minnorm = 0;
[U,S,V] = svd(A);
for ii=1:rank(A)
x_minnorm = x_minnorm + (U(:,ii)'*b)*V(:,ii)/S(ii,ii);
end
where I have used the formula above. You should find that x_minnorm=x_pinv. Note that the command rank(A) is likely using an SVD to estimate rank - I think it looks to find the number of singular values that are above some numerical error tolerance.
enjoy!
jason