Talk:LU decomposition

Python example contains error

The Python example has a divide by zero error for the matrix

b = [[1, 1, 1, 1], [1, 1, -1, -1], [1, -1, 0, 0], [0, 0, 1, -1]]

although a LUP decomposition exists: [1]

The permutation matrix has to be updated at each step, but that will make the code a lot more complicated.

Example 2 pivot matrix seems to be wrong

The pivot matrix in example 2 should also swap the last two rows of the current resulting pivoted matrix: A'(3, 3) = 2 while there is a 7 right beneath it. Here is the result of the multiplication of the two matrices on wolfram alpha

The pivot matrix I propose is {{1,0,0,0},{0,0,1,0},{0,0,0,1},{0,1,0,0}}.

I did not change the article because it seemed very strange that I would be the first to see this and I therefore wonder if I'm not wrong.

Re: Example 2 pivot matrix seems to be wrong

I noticed this too. In the implementations it appears that the max pivot values are taken from the original A matrix without swapping. When the max pivot is being found for row 3 of A the choices presented are row 3 [3, 7, 18, 1] and row 4 [2, 5, 7, 1], since row 3 contains the largest pivot value the third row of the current permutation matrix is swapped with itself, which contains [0, 1, 0, 0] from the initial swap with row 2. It's strange, but a lot of the algorithm implementations seem to follow that trend and multiplying by the permutation matrix instead of swapping inline with the permutation matrix.

I think the permutation matrix you proposed makes more sense as well.

Matlab code is wrong

Please lock at the following example A=[1 2 -1 0;2 4 -2 -1; -3 -5 6 1; -1 2 8 -2] the code returns nonsense for U, namely

U =

        0         0         0   -1.0000
   0.5000    1.0000   -0.5000         0
  -0.1429         0    1.0000    0.2857
   0.1563         0         0    1.0000

I presume the reason is the permutation matrix. I will try to correct it.

Different output as per common literature

If one sees the document at: https://ask.fxplus.ac.uk/tools/HELM/pages/workbooks_1_50_jan2008/Workbook30/30_3_lu_decmp.pdf

For A= [{1,2,4;3,8,14;2,6,13}] L= [{1,0,0;3,1,0;2,1,1}] and U= [{1,2,4;0,2,2;0,0,3}].

But on applying the LU VBA code the output is different. We get, L=[{1,0,0;0.333333333333333,1,0;0.666666666666667,-1,1}] and H=[{3,8,14;0,-0.666666666666667,-0.666666666666666;0,0,3}].

Why is this so?

First, the LU decomposition is not unique. Notably, the algorithm usually force the main diagonal of L **or** U to be all 1. But there is an additional twist: the "A=LU" decomposition may not exist. This happens if one of the pivots is zero in Gauss algorithm. Then you have to swap rows to continue. And like with Gauss algorithm, it's possible to always apply row pivoting, to get at each step the maximum pivot in absolute value, to enhance the numerical precision. The decomposition is then PA=LU, with P a permutation matrix, thus LU is the same as A with possibly some rows swapped. It's what the VBA code does, as well as other implementations. However, even if the result is on the test cases, there is a bug in the code: the maximum pivot is selected, but not in absolute value. Which means the zero pivot would be selected if all other choices were < 0. I'll change this. Eoraptor (talk) 13:46, 15 August 2020 (UTC)

Addition: the VBA code, while not incorrect per se with this modification, is still awfully inefficient (one matrix multiply for each pivot!). All the code needs to be changed. Eoraptor (talk) 13:51, 15 August 2020 (UTC)

Problems in solutions

After correcting the VBA code, I noticed the Ada, C++ and Lisp codes have the same problem, and probably others as well: when selecting the maximum pivot, it's necessary to get the absolute value, otherwise the zero pivot will be selected if it's the only nonnegative element in the possible choices. This is a **bug**.

There are other problems: the Fortran code did the permutation in a blatantly inefficient way (permute all the lower rows at each step), and the author didn't notice that if you permute the whole rows, you don't have to permute anything afterwards to get the L and U components in A: they are just the lower and upper part of A. It's a matrix trick that's easy to prove.

Overall, I fear several implementations here are badly written. Sadly, it's not the first time I notice this on RC: numerical analysis tasks are not correctly handled. Eoraptor (talk) 14:21, 15 August 2020 (UTC)