29 views (last 30 days)

Show older comments

Qiang Li on 24 Jul 2024 at 11:35

Commented: Christine Tobler on 24 Jul 2024 at 16:41

Accepted Answer: Christine Tobler

- post.mat

Open in MATLAB Online

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)

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 Comments Show -2 older commentsHide -2 older comments

Show -2 older commentsHide -2 older comments

Sign in to comment.

Sign in to answer this question.

### Accepted Answer

Christine Tobler on 24 Jul 2024 at 12:23

Open in MATLAB Online

- post.mat

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

##### 2 Comments Show NoneHide None

Show NoneHide None

Qiang Li on 24 Jul 2024 at 14:05

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2139861-solving-linear-system-with-decomposition-a-qr-and-qr-a-produce-different-results#comment_3219616

Open in MATLAB Online

- post.mat

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)];

norm(x1 - x3)

ans = 0.0024

Christine Tobler on 24 Jul 2024 at 16:41

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2139861-solving-linear-system-with-decomposition-a-qr-and-qr-a-produce-different-results#comment_3219751

Open in MATLAB Online

- post.mat

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

Sign in to comment.

### More Answers (0)

Sign in to answer this question.

### See Also

### Categories

Signal ProcessingDSP System ToolboxStatistics and Linear AlgebraLinear Algebra

Find more on **Linear Algebra** in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français

- United Kingdom(English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)

Contact your local office