# 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?

• You're sure that's Octave's output? $L$ is not lower triangular. – GNUSupporter 8964民主女神 地下教會 Oct 3 '18 at 19:36
• Thank you @GNUSupporter8964民主女神地下教會. I added the code in the question. I'm probably wrong somewhere but where? – ecjb Oct 3 '18 at 19:45
• The odd output from Octave aside, you do realize that LU factorization is not unique, right? – amd Oct 3 '18 at 19:47
• @ecjb: I wonder if Octave is permutting as it goes. Try [l, u, p] = lu (A) and see if $p$ is a permutation result. – Moo Oct 3 '18 at 19:53
• Thanks @Moo : indeed [l, u, p] = lu(A) proves that octave is permutting. For which sake? I don't know – ecjb Oct 3 '18 at 20:04

## 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:
array([[ 1. ,  0. ,  0. ],
[-0.5,  1. ,  0. ],
[-0.5,  0.5,  1. ]])

U
Out:
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:
array([[-0.5,  0.5,  1. ],
[-0.5,  1. ,  0. ],
[ 1. ,  0. ,  0. ]])

• thanks @RHowe, I agree that Python could be a good alternative but I take a class which uses Octave so I try to do as everybody do. Did you try the code in octave? – ecjb Oct 3 '18 at 19:54
• i just ran it in octave you're right. I'm looking at the documentation – Shogun Oct 3 '18 at 19:56
• updated, its the product of the permutation matrix and the L matrix. its automatically doing it. – Shogun Oct 3 '18 at 20:07
• many thanks @RHowe. Do you have a link to the reference? – ecjb Oct 3 '18 at 20:16
• octave.sourceforge.io/octave/function/lu.html – Shogun Oct 3 '18 at 20:18