LU Decomposition¶

Every sqaure matrix $A$ can be written as a product of $PLU$ where $P$ is a permutation matrix and $L$ is a lower triangluar matrix with $1$'s on its diagonal and $U$ is an upper triangular matrix.

In [1]:

%display latex
A = matrix(QQ,[[-3,-1,-2],[-1,5,-2],[1,-1,0]]); A

Out[1]:

$\displaystyle \left(\begin{array}{rrr} -3 & -1 & -2 \\ -1 & 5 & -2 \\ 1 & -1 & 0 \end{array}\right)$

In [2]:

P,L,U = A.LU() #First we show that sage can compute the LU decomposition.

In [3]:

P, L, U

Out[3]:

$\displaystyle \left(\left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right), \left(\begin{array}{rrr} 1 & 0 & 0 \\ \frac{1}{3} & 1 & 0 \\ -\frac{1}{3} & -\frac{1}{4} & 1 \end{array}\right), \left(\begin{array}{rrr} -3 & -1 & -2 \\ 0 & \frac{16}{3} & -\frac{4}{3} \\ 0 & 0 & -1 \end{array}\right)\right)$

In [4]:

P*L*U #this checks that PLU = A

Out[4]:

$\displaystyle \left(\begin{array}{rrr} -3 & -1 & -2 \\ -1 & 5 & -2 \\ 1 & -1 & 0 \end{array}\right)$

So how are we going to find this decomposition? And what is it good for?

First we tag the identity to the right of $A$ just like the first step for finding inverse using row reduction.

In [5]:

I3 = identity_matrix(3)
M = block_matrix([A,I3], ncols =2)
M

Out[5]:

$\displaystyle \left(\begin{array}{rrr|rrr} -3 & -1 & -2 & 1 & 0 & 0 \\ -1 & 5 & -2 & 0 & 1 & 0 \\ 1 & -1 & 0 & 0 & 0 & 1 \end{array}\right)$

Next we turn the left-half to an upper triangluar matrix using only the 3rd type of ERO.

In [6]:

M.add_multiple_of_row(1,0,-1/3); M

Out[6]:

$\displaystyle \left(\begin{array}{rrr|rrr} -3 & -1 & -2 & 1 & 0 & 0 \\ 0 & \frac{16}{3} & -\frac{4}{3} & -\frac{1}{3} & 1 & 0 \\ 1 & -1 & 0 & 0 & 0 & 1 \end{array}\right)$

In [7]:

M.add_multiple_of_row(2,0,1/3); M

Out[7]:

$\displaystyle \left(\begin{array}{rrr|rrr} -3 & -1 & -2 & 1 & 0 & 0 \\ 0 & \frac{16}{3} & -\frac{4}{3} & -\frac{1}{3} & 1 & 0 \\ 0 & -\frac{4}{3} & -\frac{2}{3} & \frac{1}{3} & 0 & 1 \end{array}\right)$

In [8]:

R=M.with_added_multiple_of_row(2,1,1/4);
R

Out[8]:

$\displaystyle \left(\begin{array}{rrr|rrr} -3 & -1 & -2 & 1 & 0 & 0 \\ 0 & \frac{16}{3} & -\frac{4}{3} & -\frac{1}{3} & 1 & 0 \\ 0 & 0 & -1 & \frac{1}{4} & \frac{1}{4} & 1 \end{array}\right)$

Now the left-half is an upper triangular matrix. And since we used only the 3rd type of ERO, the right-side is now a lower triangular matrix with 1 on its diagonal.

In [9]:

U = R.submatrix(0,0,ncols=3); U #U is the left-half of R.

Out[9]:

$\displaystyle \left(\begin{array}{rrr} -3 & -1 & -2 \\ 0 & \frac{16}{3} & -\frac{4}{3} \\ 0 & 0 & -1 \end{array}\right)$

The right-half $R$ is not the $L$ we want to find but the inverse of that $L$.

In [10]:

Linv=R.submatrix(0,3); Linv #this is L^{-1}

Out[10]:

$\displaystyle \left(\begin{array}{rrr} 1 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 \\ \frac{1}{4} & \frac{1}{4} & 1 \end{array}\right)$

So $L$ is the inverse of the matrix above.

In [11]:

L = Linv.inverse(); L

Out[11]:

$\displaystyle \left(\begin{array}{rrr} 1 & 0 & 0 \\ \frac{1}{3} & 1 & 0 \\ -\frac{1}{3} & -\frac{1}{4} & 1 \end{array}\right)$

In [12]:

L*U==A #check that $A = LU$.

Out[12]:

$\displaystyle \mathrm{True}$

In [13]:

L*U

Out[13]:

$\displaystyle \left(\begin{array}{rrr} -3 & -1 & -2 \\ -1 & 5 & -2 \\ 1 & -1 & 0 \end{array}\right)$