MATLAB: How to improve the accuracy of lu decomposition

lu decomposition

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

Best Answer

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.
Related Question