[Math] Why do we say SVD can handle singular matrix in least-squares? Comparison of SVD and QR decompositions

least squaresmatricesmatrix decompositionnumerical linear algebrasvd

I don't quite understand why we say that QR decomposition doesn't handle singular matrix, while SVD does when they are used for least square problem?

My example in Matlab seems to support the opposite conclusion.

Here is the information (formula) I have:
Linear function to solve using least square method: $Ax=b$

For QR decomposition: $A = QR$, where $Q$ is orthogonal and $R$ is upper triangular matrix. So we have $Rx = Q^Tb$, and we can solve this using backward substitution.

For SVD method: $A=USV$, where $U$ and $V$ are orthogonal, and $S$ is diagonal. The result is $x=VS^{-1}U^Tb$.

I implemented both methods in Matlab and even though I provided data with high condition number (e+15), QR method still provide accurate result(e-27) while SVD gets much bigger errors (0.8). So I am quite confused why we say SVD is more numerically stable?

I believe I have an incomplete knowledge of SVD(maybe there are more steps to make it better) , could anybody tell me how to let SVD deal with singular matrix? (BTW, do singular matrix, ill conditioning problem, co-linearity, big condition number mean the same thing? )

PS: this is the codes if anybody need to verify. Thank you. 

% we do a linear regression between Y and X
data= [
47.667483331 -122.1070832;
47.667483331001 -122.1070832
];
X = data(:,1);
Y = data(:,2);

X_1 =  [ones(length(X),1),X];

%%
%SVD method
[U,D,V] = svd(X_1,'econ');
beta_svd = V*diag(1./diag(D))*U'*Y;


%% QR method(here one can also use "\" operator, which will get the same result as I tested. I just wrote down backward substitution to educate myself)
[Q,R] = qr(X_1)
%now do backward substitution
[nr nc] = size(R)
beta_qr=[]
Y_1 = Q'*Y
for i = nc:-1:1
    s = Y_1(i)
    for j = m:-1:i+1
        s = s - R(i,j)*beta_qr(j)
    end
    beta_qr(i) = s/R(i,i)
end

svd_error = 0;
qr_error = 0;
for i=1:length(X)
   svd_error = svd_error + (Y(i) - beta_svd(1) - beta_svd(2) * X(i))^2;
   qr_error = qr_error + (Y(i) - beta_qr(1) - beta_qr(2) * X(i))^2;
end

Best Answer

SVD can handle rank-deficiency. In your example, there was a bug. The diagonal matrix D has a near-zero element and you need use pseudoinverse for SVD, i.e. set the 2nd element of 1./diag(D) to 0 other than the original huge value (10^14). You should find SVD and QR have equally good accuracy in this case. For more information, see this document http://www.cs.princeton.edu/courses/archive/fall11/cos323/notes/cos323_f11_lecture09_svd.pdf

Related Question