LU Decomposition: difference between between hand calculation and solver?

Question

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?

Answer 1

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. ]])