Applications of LU decomposition

1096 days ago by ncr006

Solving systems of equations via the LU decomposition

To be able to do LU decomposition in basic Sage (not numpy or scipy, say), you have to make matrices over exact rings.  AA is the set of exact algebraic reals (i.e., reals that are roots to polynomial equations with integer or rational coefficients).

A = matrix(AA, 3,3, [2,7,-3, -1,-1,5,1,2,-4 ]) b = matrix(AA,3,1,[1,2,3]) 
       
P,L,U = A.LU() 
       

Now we solve for $A\mathbf{x}=\mathbf{b}$ using this decomposition.  We start by solving $L\mathbf{z} = \mathbf{b}$.  I.e., we find a vector $\mathbf{z}=\begin{pmatrix} z_1\\z_2\\z_3\end{pmatrix}$ so that $L\mathbf{z} = \mathbf{b}$.

       
[   1    0    0]
[-1/2    1    0]
[ 1/2 -3/5    1]
[   1    0    0]
[-1/2    1    0]
[ 1/2 -3/5    1]
b1 = b[0][0] b2 = b[1][0] b3 = b[2][0] 
       

Since $L$ is lower triangular we can easily find the $\mathbf{z}$ we are looking for

z1 = b1 z2 = b2 + (1/2)*z1 z3 = b3 + (3/5)*z2- (1/2)*z1 print z1,z2,z3 
       
1 5/2 4
1 5/2 4

Now we find $\mathbf{x}=\begin{pmatrix} x_1\\x_2\\x_3\end{pmatrix}$ so that $U\mathbf{x}=\mathbf{z}$, using the $\mathbf{z}$ we just found using $L$.

       
[   2    7   -3]
[   0  5/2  7/2]
[   0    0 -2/5]
[   2    7   -3]
[   0  5/2  7/2]
[   0    0 -2/5]

Since $U$ is upper triangular you can find $\mathbf{x}$ by back substitution.

x3 = (-5/2)*z3 x2 = (2/5)*(z2 - (7/2)*x3) x1 = (1/2)*(z1 + 3*x3 - 7*x2) print x1, x2, x3 
       
-67 15 -10
-67 15 -10

Let's check our solution.

xx = vector(AA,[x1,x2,x3]) print A*xx 
       
(1, 2, 3)
(1, 2, 3)

Finding inverses via the LU decomposition

A = matrix(AA, 3,3, [2,7,-3, -1,-1,5,1,2,-4 ]) 
       

A way to think about finding the matrix inverse of $A$ is that you want the matrix $B$ so that $A\mathbf{c}_1= \begin{pmatrix} 1\\0\\\vdots\\ 0\end{pmatrix}$ where $c_1$ is the first column of the of $B$.  In general, $A\mathbf{c}_i=\mathbf{e}_i$ where $c_i$ is the $i$th column of $B$ and $e_i$ is the vector with zeros everywhere except in the $i$th row.  If in turn you can solve for $c_1,c_2,\dots$ you can construct the matrix $B$ which will be the inverse of $A$.  Let's see an example to make this believable.

B = A.inverse() show(matrix(AA,3,1,A*B.columns()[0])) show(matrix(AA,3,1,A*B.columns()[1])) show(matrix(AA,3,1,A*B.columns()[2])) 
       

                                
                            

                                

Let's find the LU decomposition.

P, L,U = A.LU() 
       

Let's construct the vector $\mathbf{e}_1$.

e = matrix(AA,3,1,[1,0,0]) e1 = e[0][0] e2 = e[1][0] e3 = e[2][0] 
       

Let's solve for $L\mathbf{z}=\mathbf{e}_1$ as we did above.

z1 = e1 z2 = e2 + (1/2)*z1 z3 = e3 + (3/5)*z2- (1/2)*z1 print z1,z2,z3 
       
1 1/2 -1/5
1 1/2 -1/5

Let's solve for $\mathbf{c}_1$ as above.

c3 = (-5/2)*z3 c2 = (2/5)*(z2 - (7/2)*c3) c1 = (1/2)*(z1 + 3*c3 - 7*c2) print c1, c2, c3 
       
3 -1/2 1/2
3 -1/2 1/2

Notice that the three numbers I just got make up the first column of the inverse of $A$.  Proceeding in this way, you'd get the second column and the third column of the inverse.

A.inverse() 
       
[   3  -11  -16]
[-1/2  5/2  7/2]
[ 1/2 -3/2 -5/2]
[   3  -11  -16]
[-1/2  5/2  7/2]
[ 1/2 -3/2 -5/2]