solving linear system with decomposition(A,'qr') and qr(A) produce different results
4 次查看(过去 30 天)
显示 更早的评论
load post.mat
x1 = decomposition(CA,'qr')\b_f;
[qq2,rr2,pp2] = qr(CA);
x2= pp2 * (rr2\(qq2'*b_f));
[qq3,rr3,pp3] = qr(CA,"econ","vector");
x3(pp3,:) = rr3\(qq3'*b_f);
any(x1-x2~=0)
any(x3-x2~=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?
0 个评论
采纳的回答
Christine Tobler
2024-7-24
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;
[qq2,rr2,pp2] = qr(CA);
x2= pp2 * (rr2\(qq2'*b_f));
norm(x1 - x2)
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)
2 个评论
Christine Tobler
2024-7-24
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))
min(svd(CA))
cond(CA)
[qq2,rr2,pp2] = qr(CA);
cond(rr2)
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!