# How to improve the accuracy of lu decomposition?

20 views (last 30 days)

Show older comments

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

##### 0 Comments

### Accepted Answer

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.

##### 0 Comments

### More Answers (2)

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)

##### 5 Comments

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.

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!