16 views (last 30 days)

Dear All,

I did the following calculation and found lu decomposition is Not accurate as I expected. For a set of linear equations (Ax=b), I know the accurate solution x_sol of it. I compared the following two cases and found lu decomposition is very Bad in accuracy: ([L,U] = lu(A))

- inv(L)*[A*x_sol - b]
- U*x_sol - inv(L)*b

I found the first case is very accurate (10e-14), but the second case has a difference as follows:

[U*x_sol - invL*b]

=

0.000000000000000 - 0.000000000000007i

-0.000000000000000 + 0.000000000000001i

-0.000037006789095 - 0.000015944272153i

-0.000970975504200 + 0.004054164127652i

0.000106135310812 + 0.000332352007079i

I hope somebody could help me to fix this problem. You have a good weekend.

Benson

Bruno Luong
on 7 Sep 2019

Edited: Bruno Luong
on 7 Sep 2019

If you want to check accuracy of LU decomposition you should check unitless quantity

error = norm(L*U-A) / norm(A)

norm(.) is the 2-norm or spectral norm.

It should not involve another vector and matrix inversion, which related to matrix conditioning.

If you want to involve a vector, this quantity error is theoretically larger than

q = norm(L*U*x-A*x) / norm(A*x)

for all x (q is the square root of the Rayleigh quotient). So the second test must meet (q is small) for all arbitrary x. This is necessary condition to show the decomposition is accurate (but not sufficient).

NOTE: there are actually 2 ways to compute L*U*x

L*(U*x)

(L*U)*x % is equal to non-pathesis expression L*U*x with MATLAB parser

They might return different results with finite floating point arithmetic. Prefer the first one for better cpu effort.

Fabio Freschi
on 6 Sep 2019

It's a weird way to check the accuracy. What about checking the norm of the residual

% data

N = 5; A = rand(N)+eye(N)+1j*rand(N); b = ones(N,1);

% factor

[L,U] = lu(A);

x0 = U\(L\b)

norm(A*x0-b)

% inverse

invA = inv(A);

x1 = invA*b;

norm(A*x1-b)

Bruno Luong
on 7 Sep 2019

Sorry: copy/past values on screen is NOT a right/correct to provide your data. We lost the precision, and a matrix can suddenly become singular/regular with such "sloppy" (to quote John) method.

I won't even bother to look at before you attach your data in MAT form (beside the fact that we need to edit your screen-capture to enter the matrix, which can lead to again error).

I hope you did correctly your test using the real accurate data and not copy/past with your data A, b, L, U etc.... before drawing a conclusion.

I also expect as John: you must make a mistake somewhere, and posting what you observe instead of telling us exactly what you did.

In computer, bugs and mistakes happen all the time, so please do provide the data/code you run for us to check and instead of drawing a conclusion with wording. This is waste of time.

John D'Errico
on 7 Sep 2019

John D'Errico
on 6 Sep 2019

Edited: John D'Errico
on 6 Sep 2019

Sorry, but it is often the case that people are sloppy in their work. That is my expectation here.

You gave us A and x_sol, but not b. But since A is non-singular, and it is not too poorly conditioned, I can get b as:

b = A*x_sol

b =

11.4615398028092 - 49.2105988921614i

-0.105039111494943 + 0.00717362381530284i

-0.132209750840726 - 0.0184213587768616i

-6.39310652188872 + 28.1486882505365i

-0.177682675601164 - 0.0177507775692129i

So I'll assume that is what you stared with.

Now , compute the LU, as

[L,U] = lu(A);

Finally, while I have no idea why you wanted to do a comparison this way (it seems to be a bit silly to me.) norm(A*x_sol - b) is the common way to do a comaprison of the accuracy. Why you are using an LU anyway is beyond me. Just use backslash.

inv(L)*[A*x_sol - b]

ans =

0

0

0

0

0

U*x_sol - inv(L)*b

ans =

0 + 0i

0 + 0i

0 - 8.88178419700125e-16i

-8.88178419700125e-16 + 0i

0 + 8.88178419700125e-16i

I fail to see the problem. There is a difference in the least significant bits. You should never trust the least significant bits of a vector when it is computed in different ways, even if the mathematics tell you the result should be the same.

Remember that floating point arithmetic is only an APPROXIMATION of real mathematics. It is not exact.

But if you got something different, then you made a mistake in what you are telling us that you did. For example, I see this in your question:

[U*x_sol - invL*b]

=

0.000000000000000 - 0.000000000000007i

-0.000000000000000 + 0.000000000000001i

-0.000037006789095 - 0.000015944272153i

-0.000970975504200 + 0.004054164127652i

0.000106135310812 + 0.000332352007079i

I'd bet a decent sum of money that the invL you used there is NOT the same matrix as inv(L), perhaps computed from a different problem. It happens. It happens a LOT, especailly with novice users, who tend to be sloppy in their code. Are you one of them? Maybe. After all, you used invL there instead of inv(L).

Again, that is a mistake that people make all the time, then they get all upset, and think there is a problem in MATLAB, when it was in fact the user who screwed things up.

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.