# How to improve the accuracy of lu decomposition?

14 views (last 30 days)
Benson Gou on 6 Sep 2019
Commented: John D'Errico on 7 Sep 2019
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))
1. inv(L)*[A*x_sol - b]
2. 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)

Benson Gou on 6 Sep 2019
A = [
Columns 1 through 4
0.381038530261726 -21.045732776083682i 0.000000000000000 + 0.000000000000000i 0.000000000000000 + 0.000000000000000i 0.840826939016329 -28.151837797820306i
0.000000000000000 + 0.000000000000000i 0.082816227356801 - 2.588225042340450i 0.000000000000000 + 0.000000000000000i -0.082816227356801 + 2.588225042340450i
-0.076328861099904 + 2.566244861142124i 0.000000000000000 + 0.000000000000000i 0.076328861099904 - 2.566244861142124i 0.000000000000000 + 0.000000000000000i
0.000000000000000 + 0.000000000000000i -0.082816227356801 + 2.588225042340450i 0.000000000000000 + 0.000000000000000i -0.686700543451740 +23.086306110927982i
0.000000000000000 + 0.000000000000000i 0.000000000000000 + 0.000000000000000i 0.000000000000000 + 0.000000000000000i -0.071310168207788 + 2.477306644551874i
Column 5
0.000000000000000 + 0.000000000000000i
0.000000000000000 + 0.000000000000000i
0.000000000000000 + 0.000000000000000i
-0.071310168207788 + 2.477306644551874i
0.071310168207788 - 2.477306644551874i]
x_sol = [ 1.005058882165182 + 0.209198956452632i
1.001631904221384 + 0.166645397313698i
1.010699880863042 + 0.157512414823176i
1.005697940832339 + 0.207098748923282i
1.010794466085600 + 0.135227908845486i]
Thanks a lot. You have a good weekend!
Benson
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
Bruno is correct. The correct way to give us your data is by putting it into a .mat file. That way we can get the exact numbers, not a 15 digit approximation as written on the command window. And you needed to give us b, not x_sol.

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.