Error in lu-decomposition function?

I need to do a lu-decomposition, receiving a lower triangular matrix with unit diagonal. The MATLAB-Function [L,R,P] = lu(A) should do exactly that.

A=[ 28   -35    21
   -12    19   -15
    16   -30    48];
[L,R,P]=lu(A);

Testing the result:

P*A==L*R

ans =

     1     1     1
     1     1     1
     1     1     0

What am I overlooking, why do the results not match?

Answer 1

[John D'Errico]

No. This is NOT an error in the lu function, but an error in how you think about floating point arithmetic.

Instead of testing for EXACT equality, why not subtract the two? What is the difference?

Do you want to make a bet that the differences are tiny numbers? In fact, it looks like only one of the elements was not exactly the same. I'll save you some time.

P*A-L*R
ans =
          0            0            0
          0            0            0
          0            0  -1.7764e-15

So welcome to the world of floating point arithmetic, where you should virtually NEVER test for exact equality.

[Benjamin Lender]

Errors are small (here 10^-14), but grow with larger matrizes with a wider range of values. LU-decomposition should be stable. Ultimately I need to do this for a 1073x1073 matrix and the results can be guessed from this 3x3. I have a "brut-force" lu-decomposition that get's there but is very slow. Any ideas?

[John D'Errico]

I'm sorry, but if you think that your brute force LU is doing much better, unless it is written in higher precision than double precision, it will not do better in general.

Even then you will ALWAYS have floating point errors when working with floating point arithmetic. It is just the size of those errors. You need to learn about floating point arithmetic if you are working with it.

Note that even though you started with an integer matrix, the LU factors are NOT integers. And almost all fractions are not representable exactly in terms of a floating point number.

For example, consider the L factor as generated.

L
L =
          1            0            0
    0.57143            1            0
   -0.42857         -0.4            1

Those numbers are not exact. That 0.57143 is only an approximation to the fraction 4/7. In fact, the actual number as stored as a double is:

0.5714285714285713968507707249955274164676666259765625000

See that it is not truly 4/7, but only an approximation.

http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

[Benjamin Lender]

Thank you very much for your hint, I am confused because I am aware of what you are floating point errors.

I am using this code in Finite Elements, and for an artificial test case both inv() and the brute force LU yield errors of below 10^-15. The MATLAB lu-function currently yields 80% error.

If this was just floating point error, I should not be able to get such low an error with my brute force LU.

[John D'Errico]

As far as the errors growing with size, again, this is completely expected. The errors involved in ANY linear algebra solve like an LU factor (not just that in MATLAB) will grow with the condition number of the matrix. And generally the condition number of some random matrix will tend to grow with the size of that matrix.

[John D'Errico]

It IS just floating point error. The error that you found when you used your code was still non-zero, just arbitrarily different. I'm not sure what you mean by that 80% number. And of course, I have no way to know what your own code does. For example, perhaps you used higher precision to accumulate things. This takes more time of course.

You might also want to read section 5.1 of the MATLAB FAQ.

http://www.mit.edu/~pwb/cssm/matlab-faq_6.html#eps

[Walter Roberson]

Your "brute-force" decomposition is going to have the same problem unless you convert all of your entries from floating point into rational values and use a symbolic solver that calculates in rational form.

For example the exact L matrix is [[1, 0, 0]; [4/7, 1, 0]; [-3/7, -2/5, 1]] and you know that there is no finite binary floating point representation for 1/7 so of course errors are going to accumulate.

[Benjamin Lender]

That is why I wonder. I am getting a distinct error with the lu-function, and I should be getting the same kind of result.

These are my results that make me suspect some other error. For each pair the smooth comes from the "brute-force" lu, the other from the lu-function, that is the only difference in calculation. If this doesn't mean anything, please let me know. I am learning to use MATLAB.

I hope the image upload works as intended.

The global errors are closer than I first thought they were.

[Benjamin Lender]

[Benjamin Lender]

[Benjamin Lender]

[Benjamin Lender]

I know code errors won't be identified by these photos, I want to find out whether this difference is considered significant.

Only difference is using "brute force" (slow) lu-decomposition and using the lu-function from MATLAB. To my understanding of floating point arithmetic, both should "suffer" from those errors equally.

The code for the straight approach comes from the Book "Höhere Mathematik in Rezepten" by Christian Karpfinger

[n,n] = size(A);
for j = 1:n-1
    I = j+1:n;
    A(I,j) = A(I,j)/A(j,j);
    A(I,I) = A(I,I)-A(I,j)*A(j,I);
end

Ri = triu(A);
Le = eye(n,n) + tril(A,-1);

I am asking to learn

[Christine Tobler]

The problem with this brute-force lu decomposition is that it doesn't do permutation, which can lead to completely wrong results. For example, try this matrix:

A = [0 0 1; 1 0 0; 0 1 0];

The brute-force LU will return matrices L and U with Inf and NaN values, while the lu function's output is still correct. In less extreme examples, not using permutation will simply lead to high levels of numerical error.

For the plots you show, I assume from your answer that the first is using your brute-force algorithm and the second the output from lu? In that case, I would assume there is something wrong with the way you are using L, U and P to solve the linear system. Can you post your formula for this?

If you only need to solve one linear system, you can also just use

x = A \ b;

which will solve the linear system A*x = b directly without any worries about the decomposition used.

[John D'Errico]

Yes. the pivoted LU will be more accurate, although still not perfect. Nothing can achieve that, unlesss you use exact arithmetic.