LU decomposition

Simple matrix decomposition

Source

Intro

LU (lower–upper) decomposition (factorization) outputs (factors original matrix into) lower and upper triangular matrix. These matrices could be used to efficiently solve system of non-sparse linear systems or find inverse of a matrix. It might sound a bit difficult, but we’ll have an example later. Goal of this post is to show what LU decomposition is about, how it could be calculated and used solving practical tasks.

Gaussian elimination

In LU decomposition we want to decompose original into upper and lower triangular matrices, so that:

A = LU, where:

• A is original matrix we want to decompose
• L is lower triangular matrix (we assume it has 1-s in diagonal)
• U is upper triangular matrix

In simplest form LU decomposition could be calculated using Gaussian elimination. Gaussian elimination allows three types of operations:

• swapping two rows
• adding multiple of one row to other
• multiplying a row with a nonzero number

Here we have a simple example of how Gaussian elimination gives upper triangular matrix. We start with simple matrix:

[[ 1,  1,  0],
 [ 2,  1, -1],
 [ 3, -1, -1]]

We want to have all elements below main diagonal to be zero, lets do following steps:

• subtract 2 times row 1 from row 2:

[[ 1,  1,  0],
 [ 0, -1, -1],
 [ 3, -1, -1]]

• subtract 3 times row 1 from row 3 (1st column has now 0s except first element):

[[ 1,  1,  0],
 [ 0, -1, -1],
 [ 0, -4, -1]]

• subtract 4 times row 2 from row 3 and we are done (now matrix is in row echelon form):

[[ 1,  1,  0],
 [ 0, -1, -1],
 [ 0,  0,  3]]

We have now upper triangular matrix. Getting lower triangular matrix is not very difficult, we have to organize coefficients from the previous steps into identity matrix:

• start with identity matrix:

[[ 1,  0,  0],
 [ 0,  1,  0],
 [ 0,  0,  1]]

• from backwards with previous Gaussian elimination steps to go from upper triangular matrix to original matrix, we should first add 4 times row 2 to row 3:

[[ 1,  0,  0],
 [ 0,  1,  0],
 [ 0,  4,  1]]

• then we should add 3 times row 1 to row 3:

[[ 1,  0,  0],
 [ 0,  1,  0],
 [ 3,  4,  1]]

• finally we should add 2 times row 1 to row 2 and we are done:

[[ 1,  0,  0],
 [ 2,  1,  0],
 [ 3,  4,  1]]

And now we have lower and upper triangular matrices. In python Gaussian elimination could be implemented like that:

Example usage for this function:

A = np.array([[60, 91, 26], [60, 3, 75], [45, 90, 31]], dtype='float')
L, U = lu(A)
L
#[[ 1.  ,  0.        ,  0. ],
# [ 1.  ,  1.        ,  0. ],
# [ 0.75, -0.24715909,  1. ]]
U
#[[ 60. ,  91.  ,  26.        ],
# [  0. , -88.  ,  49.        ],
# [  0. ,   0.  ,  23.61079545]]
#we should get original matrix L@U
np.allclose(A, L@U)
#True

You should compare the results with scipy.linalg.lu and see that the results are same. This is very simple LU decomposition and it has a problem. For example following matrix doesn’t return usable result.

A = np.array([[0, 91, 26], [60, 3, 75], [45, 90, 31]], dtype='float')
L, U = lu(A)
L
#[[ 1.,  0.,  0.],
# [inf,  1.,  0.],
# [inf, nan,  1.]]
U
#[[  0.,  91.,  26.],
# [ nan, -inf, -inf],
# [ nan,  nan,  nan]]

For example if we have 0 on diagonal element, our simple function doesn’t work. For that example to work we need a permutation matrix.

PLU

From PLU decomposition we create permutation matrix with upper and lower diagonal matrix. Permutation stands for swapping rows. PLU formula is:

A = PLU, where:

• P stands for permutation matrix

In our case we could swap rows so that PLU decomposition doesn’t break if there is 0 on diagonal. In python PLU implementation might look like this:

We just reused lu function and added row permutation part. We could check if now we get something useful if there is 0 on the diagonal:

A = np.array([[0, 91, 26], [60, 3, 75], [45, 90, 31]], dtype='float')
P, L, U = plu(A)
P
#[[0., 1., 0.],
# [1., 0., 0.],
# [0., 0., 1.]]
L
#[[1.        , 0.        , 0.  ],
# [0.        , 1.        , 0.  ],
# [0.75      , 0.96428571, 1.  ]]
U
#[[ 60.      ,   3.        ,  75.        ],
# [  0.      ,  91.        ,  26.        ],
# [  0.      ,   0.        , -50.32142857]]

You can verify that you should get similar result using scipy.linalg.lu function. We have looked LU decomposition from the Gaussian elimination perspective. As in many cases in linear algebra, there some other way.

Doolittle algorithm

We could get LU decomposition by using lower and upper triangular matrices specifics. From the A = LU we know:

LU decomposition matrix elements. Source

And from there we can write down all the element calculations:

LU decomposition calculations. Source

Now because of the triangular matrices we could “hack” our way to solution. Some calculations are very simple and we could move from getting first elements to calculating next elements. We know that l₁₁=1 (all elements in lower diagonal matrix are 1s). From that we could get:

• u₁₁=a₁₁, u₁₂=a₁₂ and u₁₃=a₁₃, we have first row of U
• as we now know u₁₁, we can find l₂₁=a₂₁/u₁₁ (which is same as a₂₁/a₁₁) and l₃₁=a₃₁/u₁₁ (which is same as a₃₁/a₁₁)
• now we can move on to second row of U. we know that l₂₁u₁₂+l₂₂u₂₂=a₂₂ and that l₂₂=1, l₂₁=a₂₁/a₁₁ and u₁₂=a₁₂ then u₂₂=a₂₂-a₂₁a₁₂/a₁₁, we can continue this approach and find all the elements.

This approach is known as Doolittle algorithm. In python example implementation might look like so (you might want to have a deeper look at it because this code is compact and uses dot-product instead of looping elements separately):

I’ll leave it to users to verify that this algorithm produces valid LU decomposition. Note that this implementation assumes that matrix is square and we don’t do the permutation part. Now let’s see some examples where LU decomposition could have practical value.

Solving linear systems

LU decomposition is used for solving equation of linear systems. We have:

Ax = b, where

• A and b are known
• x is unknown, we want to find it

With LU decomposition we could do:

LUx = b and solve it in two parts:

• Ux = y using backward substitution
• Ly = b using forward substitution

If we go from back to the beginning we could see that again using triangular matrix we can “hack” the solution. We have Ly = b and L, b are known:

[[L11, 0,     0],  [[y1],    [[b1],
 [L21, L22,   0],   [y2],  =  [b2],
 [L31, L32, L33]]   [y3]]     [b3]]

We can find values of y:

• L11y1=b1, so y1=b1/L11,
• L21y1+L22y2=b2, so y2=(b2-L21y1)/L22
• etc.

This is forward substitution. Sample python implementation:

Now for last part, doing backward substitution. We have Ux = y, we know U and y and need to find x:

[[U11, U12,  U13],  [[x1],    [[y1],
 [0,   U22,  U23],   [x2],  =  [y2],
 [0,     0,  U33]]   [x3]]     [y3]]

For finding values of x, we can start from the beginning of the matrix:

• U33x3=y3, so x3=y3/U33
• U22x2+U23x3=y2, so x2=(y2-U23x3)/U22
• etc.

This is called backward substitution. In python might be implemented like this:

We could test this code with the following example:

def plu_solve(A, b):
    P, L, U = plu(A)
    y = forward_substitution(L, P@ b)
    return back_substitution(U, y)
A = np.array([[1, 5, 5], [6, 9, 22], [32, 5., 5]])
b = np.array([1, 2, 7])
plu_solve(A, b)
#[ 0.19354839,  0.20843672, -0.0471464 ]
np.allclose(plu_solve(A, b), scipy.linalg.solve(A, b))
#True

So we could use LU decomposition to solve system of linear equations. Next we’ll see how to use it finding matrix inverse.

Matrix inverse

LU decomposition helps to calculate matrix inverse. We know that:

AA¹ = I, where

• A stands for original matrix
• A¹ stands for inverse of original matrix
• I stands for identity matrix

Using LU we could write it as:

LUA¹ = I

Which could be written as (source):

[[1,     0, 0],  [[U11, U12, U13],
 [L21,   1, 0],   [  0, U22, U23],  [v1, v2, v3] = [e1, e2, e3]
 [L31, L32, 1]]   [  0,   0, U33]]

Elements v1, v2, v3 are columns of A¹ and e1, e2, e3 are columns of I (in our example e1=[1, 0, 0]). So we know that L⋅U⋅vi=e. Now we can use forward and backward substitution to find all values of v.

We can check if it gives correct values:

A = np.array([[2, -1, 0],[-1, 2, -1.], [0, -1, 2.]])
plu_inverse(A)
#[[0.75, 0.5 , 0.25],
# [0.5 , 1.  , 0.5 ],
# [0.25, 0.5 , 0.75]]
np.allclose(plu_inverse(A), np.linalg.inv(A))
#True

LU decomposition has many practical usages. It is relatively fast algorithm which makes it useful (source).

Conclusion

LU factorization is a fast decomposition which is used for solving linear equations and finding inverses of matrices. Calculations behind this algorithm is not very complicated. In practical cases you want to use scipy.linalg.lu because it is faster and supports non-square matrices. Examples here have only educational purpose.

References

Backward substitution. AlgoWiki

Easy way to calculate inverse of an LU decomposition., math.stackexchange.com

Forward substitution, AlgoWiki

Gaussian elimination, Wikipedia

How can LU factorization be used in non-square matrix?, math.stackexchange.com

LU decomposition, Wikipedia

LU factorization, johnfoster.pge.utexas.edu

Row echelon form, Wikipedia