LU Decomposition: difference between between hand calculation and solver?

I have a $3 \times 3$ matrix $A$ and have to perform the $LU$ Factorization (1)

$$ A = \begin{bmatrix} 1 & 2 & 3 \\ 1 & -1 & 3 \\ -2 & -10 & -2 \end{bmatrix}$$

Using row reduction, I would first substract the second row by (1/1) the first row which gives (2):

$$ \begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & 0 \\ -2 & -10 & -2 \end{bmatrix}$$

I would then substract the third row by (-2/1) the first row which gives (3)

$$ \begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & 0 \\ 0 & -6 & 4 \end{bmatrix}$$

Finally I would substract the third row by (-6/2) the second row which would give (4):

$$ \begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & 0 \\ 0 & 0 & 4 \end{bmatrix} = U$$

with the factor found in (2) (3) and (4) I can get the matrix $U$

$$L = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ -2 & 2 & 1 \end{bmatrix}$$

I can verify that $L \times U = A$

However when I compute [L, U] = lu(A) with Octave, I get $$ U = \begin{bmatrix} -2 & -10 & -2 \\ 0 & -6 & 2 \\ 0 & 0 & 1 \end{bmatrix}$$

and

$$ L= \begin{bmatrix} -0.5 & 0.5 & 1 \\ -0.5 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}$$

Here is the octave / matlab code:

A = [1 2 3; 1 -1 3; -2 -10 -2];
U_hand = [1 2 3; 0 -3 0; 0 0 4];
L_hand = [1 0 0; 1 1 0; -2 2 1];
L_hand * U_hand
[L, U] = lu(A)

How can how explain the differences? I'm probably wrong somewhere but where?

1 Answer

You shouldn't have got that for your LU decomp. I used python which uses the same LAPACK

import scipy.linalg
import
A = scipy.array([[1 ,2,3],[1, -1, 3 ] ,[-2,-10,-2]])
P,L,U = scipy.linalg.lu(A)
L
Out[3]:
array([[ 1. , 0. , 0. ], [-0.5, 1. , 0. ], [-0.5, 0.5, 1. ]])
U
Out[4]:
array([[ -2., -10., -2.], [ 0., -6., 2.], [ 0., 0., 1.]])

Upon further inspection. The reference says

It's the product of the permutation matrix and the L matrix

With two output arguments, returns the permuted forms of the upper and lower triangular matrices, such that A = L * U. With one output argument y, then the matrix returned by the LAPACK routines is returned. If the input matrix is sparse then the matrix L is embedded into U to give a return value similar to the full case. For both full and sparse matrices, lu loses the permutation information.

Which leads to testing it in Python.

Atest = np.dot(P,L)
Out[2]:
array([[-0.5, 0.5, 1. ], [-0.5, 1. , 0. ], [ 1. , 0. , 0. ]])