LU decomposition in Python

In linear algebra, we define LU (Lower-Upper) decomposition as the product of lower and upper triangular matrices. In this tutorial, we will learn LU decomposition in Python. Computers use LU decomposition method to solve linear equations.

How to solve LU decomposition?

Let us, first see some algebra. Mainly two methods are used to solve linear equations: Gaussian elimination and Doolittle method/ LU decomposition method. As defined, LU is a product of upper and lower triangular matrices. At times, permutation matrix is included as well. Hence, the equation looks something like this:

A = PLU, where A is a square matrix, L and U are it’s upper and lower triangular matrices respectively and P is the permutation matrix.

When linear equations are given, we write in the form of Ax = B. Using LU decomposition, we know that PA = LU.

Permutation matrix : There should be single 1 in each row and column. Rest of the elements are 0. This matrix is needed to solve some singularity issues.

Upper triangular matrix : All the elements below the main diagonal should be 0. To calculate upper triangle, we use the formula:

uij = aij − ∑k=1i=1  (ukj lik)

Lower triangular matrix : Similar to upper triangular matrix, in lower triangular matrix, all the elements above the main diagonal should be 0. The formula for this matrix is also similar. Each element is just divided by corresponding diagonal element of U.

uij = 1/ujj [aij − Σk=1i=1  (ukj lik)]

A matrix is a 2D structure consisting of rows and columns. Python does not have a built-in function called matrix. Nonetheless, we can create lists or arrays instead of matrix. This can be done by using array() method.

LU decomposition in Python with SciPy Library

Scipy library-Scientific library for Python

Scipy is an open source library in Python used for mathematical calculations, scientific computing and engineering. It contains all the features of numpy including some additional features. Hence, it is faster and more preferred than numpy. Since it is built on numpy, both of them can work hand in hand. We would need this library to prove LU decomposition.

The Scipy library holds many packages available to help in scientific computing. One such built-in package is linalg. Linalg enables solving linear algebra routines very quickly. One such linear algebra function is solving LU. It can easily be computed using lu() method. This method automatically computes P, L and U.

We have already seen how to create matrix above. To multiply two matrices, we use dot() method. Matmul() can also be used. It is to be noted that “*” only multiplies corresponding elements of the matrices and cannot be used for matrix multiplication.

Let us see the formula of LU decomposition in Python part by part. A = PLA

import scipy.linalg
A = scipy.array([[1, 2, 3],
             [4, 5, 6],
             [10, 11, 9]])

P, L, U = scipy.linalg.lu(A)
print(P)
print(L)
print(U)
[[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]

[[1.         0.         0.        ]
 [0.1        1.         0.        ]
 [0.4        0.66666667 1.        ]]

[[10.  11.   9. ]
 [ 0.   0.9  2.1]
 [ 0.   0.   1. ]]

The above outputs P, L and U in order. All of them fits their definition perfectly. Now, let us multiply them and check if we get the original matrix A back. Here is the whole code:

A = scipy.array([[1, 2, 3],
             [4, 5, 6],
             [10, 11, 9]])

P, L, U = scipy.linalg.lu(A)
mult = P.dot((L.dot(U)))
print(mult)

Output:

[[ 1.  2.  3.]
 [ 4.  5.  6.]
 [10. 11.  9.]]

Thus, we successfully get A back by multiply P, L and U. This proves the LU decomposition.

Also read:

Leave a Reply

Your email address will not be published. Required fields are marked *