solving linear system with decomposition(A,'qr') and qr(A) produce different results

load post.mat
x1 = decomposition(CA,'qr')\b_f;

Warning: Rank deficient, rank = 44, tol = 1.470347e-01.

[qq2,rr2,pp2] = qr(CA);
x2= pp2 * (rr2\(qq2'*b_f));

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.100607e-25.

[qq3,rr3,pp3] = qr(CA,"econ","vector");
x3(pp3,:) = rr3\(qq3'*b_f);

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.100607e-25.

any(x1-x2~=0)

ans = 1x3 logical array
   1   1   1

any(x3-x2~=0)

ans = 1x3 logical array
   0   0   0

I know that the CA matrix is ill-conditioned. But that doesn't explain the difference in solution, right?

Also, solving the system using decomposition(CA,'lu') and lu(CA) produce the same results. So why not the 'qr' pair?

Answer 1

[Christine Tobler]

There are two reasons that the results don't match:

1) When the matrix is not full-rank, the QR-based solver in decomposition only uses the upper left triangle of the R matrix returned from QR. I have replicated this in output x3 in the code below.

2) The matrix Q in decomposition is stored in a more compact format (Householder reflectors) than the one returned by the QR function. This results in some differences on the level of round-off error.

load post.mat
x1 = decomposition(CA,'qr')\b_f;

Warning: Rank deficient, rank = 44, tol = 1.470347e-01.

[qq2,rr2,pp2] = qr(CA);
x2= pp2 * (rr2\(qq2'*b_f));

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.100607e-25.

norm(x1 - x2)

ans = 4.3604e+09

rk = 44;
m = size(CA, 1);
nrhs = size(b_f, 2);
[qq2,rr2,pp2] = qr(CA);
x3 = pp2 * [(rr2(1:rk, 1:rk)\(qq2(:, 1:rk)'*b_f)); zeros(m-rk, nrhs)];
norm(x1 - x3)

ans = 1.2555e-13

[Qiang Li]

Hi Christine,

I used the 'RankTolerance' option of decomposition, and compared the results. Is the 0.0024 difference coming from your second point? Thanks.

load post.mat
min(svd(CA))

ans = 3.2181e-12

x1 = decomposition(CA,'qr','RankTolerance',1e-12)\b_f;
rk = 60;
m = size(CA, 1);
nrhs = size(b_f, 2);
[qq2,rr2,pp2] = qr(CA);
x3 = pp2 * [(rr2(1:rk, 1:rk)\(qq2(:, 1:rk)'*b_f)); zeros(m-rk, nrhs)];

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.100607e-25.

norm(x1 - x3)

ans = 0.0024

[Christine Tobler]

Yes, it's initially coming from point (2). This is amplified by the fact that matrix CA is ill-conditioned (max / min singular value is about 3e24), which means that the full matrix R is also ill-conditioned, and so the application of R\(Q'*b) takes a small difference in Q'*b and amplifies it a lot.

While RankTolerance=1e-12 doesn't seem particularly tiny, this is an absolute tolerance (not scaled by the maximum singular value of A), meaning it corresponds to about 1e-25 in relative terms.

load post.mat
max(svd(CA))

ans = 1.1161e+13

min(svd(CA))

ans = 3.2181e-12

cond(CA)

ans = 3.4681e+24

[qq2,rr2,pp2] = qr(CA);
cond(rr2)

ans = 3.4532e+24