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)
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)
    beta_qr(i) = s/R(i,i)

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;
| cite | improve this question | | | | |
  • $\begingroup$ you don't understand why with $USV^*$ the SVD of $A$, i.e. $U SS^* U^*$ the SVD of $A A^T$ : $\min_x \|Ax-b\|^2$ with $x = A^T y$ becomes $\min_x \|AA^Ty -b\|^2 = \min_y \|U S S^*U^* y -b\|^2 = \min_y \|S S^*U^* y -U^* b\|^2 = \min_z \|S S^* z -U^* b\|^2 = \min_z \sum_{k=1}^K |\, |S_{k,k}|^2 z_k- (U^*b)_k|^2$ which is easy to solve even when $S$ is singular $\endgroup$ – reuns Apr 6 '16 at 4:02
  • $\begingroup$ A singular matrix means it has 0 for an eigenvalue, its determinant is zero, and that the matrix is non-invertible. The condition number of a matrix measures how the effect of a small change in the elements of a matrix changes the output for an inverse. $\endgroup$ – Xoque55 Apr 6 '16 at 4:09
  • $\begingroup$ @Xoque55 Are cond number and singularity somehow related? Say if a matrix have a big cond number then it is near singularity? $\endgroup$ – daydayup Apr 6 '16 at 13:13
  • $\begingroup$ @user1952009 hmmmm.... I used the formula $x = VS^{-1}U^Tb$, which means $S$ can not be singular. How should I solve it if $S$ is non-invertible? Thank you. $\endgroup$ – daydayup Apr 6 '16 at 13:21
  • $\begingroup$ this was cross-posted on Stack Overflow: stackoverflow.com/q/36452701/97160 $\endgroup$ – Amro Apr 11 '16 at 18:51

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

| cite | improve this answer | | | | |

While in theory SVD can solve Linear System which is singular we live in a practical world where a computer solves the equation.

Solving linear equation deals with the conditional number of the equation and if it high numerical issues arise (stability as well).

So in practice you need to use the SVD but also not to try inverse small singular values (Small is relative to the largest one) bu clamp them before inverting the Singular Values matrix.

This is basically what being done when one use the tol parameter of MATLAB's pinv() function (Which solves Linear System using the SVD).

| cite | improve this answer | | | | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.