In my last post, I showed how to take a matrix X and decompose it into two triangular matrices L (lower triangular) and U (upper triangular), such that X = LU
Can we get the inverse of X from L and U? Yes, indeedy.
The inverse of X equals the inverse of LU. One of the properties of matrices is that the inverse of the product of two matrices equals the product of the two inverses, where the ordering of the matrices is switched.
That is: if X-1 = (LU)-1 then X-1 = U-1L-1.
It turns out that because U and L are triangular, that calculating their inverses is relatively simple.
For example, if L =
L = [1 0 0]
|0.5 1 0|
[0.75 -0.3 1]
and we know that LL-1= the identity matrix I, then we can calculate the top left value as:
L[1,1] * L-1[1,1] + L[1,2]*L-1[2,1] + L[1,3]*L-1[3,1] = 1
L-1[1,1] = (1 - (L[1,2]*L-1[2,1] + L[1,3]*L-1[3,1])) / L[1,1]
= (1 - (0.0 + 0.0)) / 1 = 1
Moving down the column, we can solve for each of the values, based on the calculation of the previous values.
So for row 2 we have:
L[2,1] * L-1[1,1] + L[2,2]*L-1[2,1] + L[2,3]*L-1[3,1] = 0
L-1[2,1] = (0 - (L[2,1]*L-1[1,1] + L[2,3]*L-1[3,1])) / L[2,2]
= (0 - (0.5*1 + 0.0)) / 1 = -0.5
And for row 3:
L[3,1] * L-1[1,1] + L[3,2]*L-1[2,1] + L[3,3]*L-1[3,1] = 0
L-1[3,1] = (0 - (L[3,1]*L-1[1,1] + L[3,2]*L-1[2,1])) / L[3,3]
= (0 - (0.75*1 + -0.3*-0.5)) / 1 = -.9
For column 2, we have
Row 1:
L[1,1] * L-1[1,2] + L[1,2]*L-1[2,2] + L[1,3]*L-1[3,2] = 0
L-1[1,2] = (0 - (L[1,2]*L-1[2,2] + L[1,3]*L-1[3,2])) / L[1,1]
= (0 - (0.0 + 0.0)) / 1 = 0
Row 2:
L[2,1]*L-1[1,2] + L[2,2]*L-1[2,2] + L[2,3]*L-1[3,2] = 1
L-1[2,2] = (0 - (L[2,1]*L-1[1,2] + L[2,3]*L-1[3,2])) / L[2,2]
= (1 - (0.5*0 + 0.0)) / 1 = 1.0
Row 3:
L[3,1]*L-1[1,2] + L[3,2]*L-1[2,2] + L[3,3]*L-1[3,2] = 0
L-1[3,2] = (0 - (L[3,1]*L-1[1,2] + L[3,2]*L-1[2,2])) / L[3,3]
= (0 - (0.75*0 + -0.3*1.0)) / 1 = 0.3
For column 3, we have
Row 1:
L[1,1] * L-1[1,3] + L[1,2]*L-1[2,3] + L[1,3]*L-1[3,3] = 0
L-1[1,3] = (0 - (L[1,2]*L-1[2,3] + L[1,3]*L-1[3,3])) / L[1,1]
= (0 - (0.0 + 0.0)) / 1 = 0
Row 2:
L[2,1]*L-1[1,3] + L[2,2]*L-1[2,3] + L[2,3]*L-1[3,3] = 0
L-1[2,3] = (0 - (L[2,1]*L-1[1,3] + L[2,3]*L-1[3,3])) / L[2,2]
= (0 - (0.5*0 + 0.0)) / 1 = 0.0
Row 3:
L[3,1]*L-1[1,3] + L[3,2]*L-1[2,3] + L[3,3]*L-1[3,3] = 1
L-1[3,3] = (1 - (L[3,1]*L-1[1,3] + L[3,2]*L-1[2,3])) / L[3,3]
= (1 - (0.75*0 + -0.3*0)) / 1 = 1.0
Thus we have
L-1 = [1 0 0]
|-0.5 1 0|
[-0.9 0.3 1]
for the inverse of L
Jeez, I sure hope I got all those darn little array indices right.
How do we calculate the inverse of U?
Simple. First transpose U into U'. U' is lower triangular. Just apply the above method to calculate the inverse of U' = U'-1. Once the inverse is calculated, just transpose it and the result is U-1.
So now let's write the F# code to do it! I have borrowed a lot of code from my previous post. Here is the new stuff:
-----------------------------------------------------------------------------
let rec doIColumn (L : float [,]) (L0 : float [,]) (X : float [,]) (i : int) (j : int) (k : int) (ret: int) =
match i with
| i when ret <> 0 -> L, ret
| i when i > k -> L, ret
| _ -> let ret0 =
if checkZero L.[i,i] = 0 then
L0.[i,j] <- doCell L L0 X i j (i-1) L.[i,i]
0
else
-2
let L, xRet = doIColumn L L0 X (i+1) j k ret0
L0, xRet
let rec doIColumns (L : float [,]) (L0 : float [,]) (X : float [,]) (j : int) (k : int) =
match j with
| j when j > k -> L0, 0
| _ -> let L0, ret = doIColumn L L0 X 0 j k 0
if ret = 0 then
doIColumns L L0 X (j+1) k
else
L0, ret
let invertL (L : float [,]) =
let k = L.GetUpperBound 0
let L0 = Array2D.zeroCreate<float> (k+1) (k+1)
let I = Array2D.zeroCreate<float> (k+1) (k+1)
//Initialize the Diagonal in L
for j = 0 to k do
I.[j, j] <- 1.0
//Work over the columns one-by-one
let L0, ret = doIColumns L L0 I 0 k
L0, ret
-----------------------------------------------------------------------------
The function invertL takes a lower triangular matrix L as input. We then initialize two square matrices L0 (which will store our inverse) and I which is an identity matrix. Similar to our linear decomposition matrix, we then iterate over the columns from left to right, calling the recursive function doIColumns.
The function doIColumns takes the L matrix, its inverse placeholder L0, the identity matrix I, a column incrementer j, and the maximum dimension k as inputs. This function processes a single column and then moves on to the next column until it reaches its maximum dimension. The function that actually processes the column is doIColumn.
The function doIColumn iterates over each row in a column and calculates its value. It uses the function doCell which is in my prior post. This is the same function that we used to calculate a cell in our linear decomposition function. Anyhow the doCell function accumulates the products of the L and L0 matrix for rows/columns whose index is above the index of the row currently being computed.
So now we need to modify our test program:
-----------------------------------------------------------------------------
let doDecomp =
let ret, L, U = liDecomp X
printfn "ret: %i " ret
let k = X.GetUpperBound 0
for i = 0 to k do
printfn "L Matrix: %10.2f %10.2f %10.2f %10.2f " L.[i,0] L.[i,1] L.[i,2] L.[i,3]
for i = 0 to k do
printfn "U Matrix: %10.2f %10.2f %10.2f %10.2f " U.[i,0] U.[i,1] U.[i,2] U.[i,3]
let zeroX = Array2D.zeroCreate<float> (k+1) (k+1)
let LU = matMult L U zeroX
for i = 0 to k do
printfn "LU Matrix: %10.2f %10.2f %10.2f %10.2f " LU.[i,0] LU.[i,1] LU.[i,2] LU.[i,3]
//--------------------------------
//Now Invert L
//--------------------------------
let LI, ret_LI = invertL L
for i = 0 to k do
printfn "LI Matrix: %10.2f %10.2f %10.2f %10.2f " LI.[i,0] LI.[i,1] LI.[i,2] LI.[i,3]
//-----------------------------------------
//Now Invert U
//First Transpose U - UP is U Transpose
//-----------------------------------------
let Uzero = Array2D.zeroCreate<float> (k+1) (k+1)
let UP = transpose U Uzero k k
for i = 0 to k do
printfn "UP Matrix: %10.2f %10.2f %10.2f " UP.[i,0] UP.[i,1] UP.[i,2]
//-----------------------------------------
//Now Invert UP
//Return UPI as Inverse of UP
//-----------------------------------------
let UPI, ret_UI = invertL UP
for i = 0 to k do
printfn "UPI Matrix: %10.5f %10.5f %10.5f %10.5f " UPI.[i,0] UPI.[i,1] UPI.[i,2] UPI.[i,3]
//-----------------------------------------
//Now Transpose UPI
//UI is now U inverse
//-----------------------------------------
let UI = transpose UPI Uzero k k
for i = 0 to k do
printfn "UI Matrix: %10.5f %10.5f %10.5f %10.5f " UI.[i,0] UI.[i,1] UI.[i,2] UI.[i,3]
//-----------------------------------------
//Multiply UI and LI to calculate XI
//XI is now X inverse
//-----------------------------------------
let XI = Array2D.zeroCreate<float> (k+1) (k+1)
let XI = matMult UI LI XI
for i = 0 to k do
printfn "XI Matrix: %10.5f %10.5f %10.5f %10.5f " XI.[i,0] XI.[i,1] XI.[i,2] XI.[i,3]
//-----------------------------------------
//Now Test
//Multiply XI and X to calculate Identity Matrix
//-----------------------------------------
let XTest = Array2D.zeroCreate<float> (k+1) (k+1)
let XTest = matMult XI X XTest
for i = 0 to k do
printfn "XTest Matrix: %10.5f %10.5f %10.5f %10.5f " XTest.[i,0] XTest.[i,1] XTest.[i,2] XTest.[i,3]
ret, L, U
This program takes the same 4x4 matrix in my prior post and:
(1) Does the linear decomposition of X into L & U
(2) Calculates the inverse of L, returning L0 = L-1
(2) Calculates the transpose of U = UP
(4) Calculates the inverse of UP = UPI
(5) Then takes the transpose of UPI to recover the inverse of U = UI
(6) Multiplies UI and LI to calculate what we really want X-1= XI
(7) Runs a test by multiplying XI and X to make sure it is an identity matrix.
(8) Voilà! Done.
Here is the output:
-------------------------------------------------------------
ret: 0
L Matrix: 1.00 0.00 0.00 0.00
L Matrix: 0.50 1.00 0.00 0.00
L Matrix: 0.75 -0.30 1.00 0.00
L Matrix: 2.25 -2.37 -3.93 1.00
U Matrix: 4.00 7.00 -8.00 6.00
U Matrix: 0.00 7.50 6.00 6.00
U Matrix: 0.00 0.00 -9.20 7.30
U Matrix: 0.00 0.00 0.00 34.42
LU Matrix: 4.00 7.00 -8.00 6.00
LU Matrix: 2.00 11.00 2.00 9.00
LU Matrix: 3.00 3.00 -17.00 10.00
LU Matrix: 9.00 -2.00 4.00 5.00
LI Matrix: 1.00 0.00 0.00 0.00
LI Matrix: -0.50 1.00 0.00 0.00
LI Matrix: -0.90 0.30 1.00 0.00
LI Matrix: -6.97 3.55 3.93 1.00
UP Matrix: 4.00 0.00 0.00
UP Matrix: 7.00 7.50 0.00
UP Matrix: -8.00 6.00 -9.20
UP Matrix: 6.00 6.00 7.30
UPI Matrix: 0.25000 0.00000 0.00000 0.00000
UPI Matrix: -0.23333 0.13333 0.00000 0.00000
UPI Matrix: -0.36957 0.08696 -0.10870 0.00000
UPI Matrix: 0.07547 -0.04168 0.02305 0.02905
UI Matrix: 0.25000 -0.23333 -0.36957 0.07547
UI Matrix: 0.00000 0.13333 0.08696 -0.04168
UI Matrix: 0.00000 0.00000 -0.10870 0.02305
UI Matrix: 0.00000 0.00000 0.00000 0.02905
XI Matrix: 0.17293 -0.07652 -0.07262 0.07547
XI Matrix: 0.14577 0.01158 -0.07704 -0.04168
XI Matrix: -0.06294 0.04915 -0.01800 0.02305
XI Matrix: -0.20261 0.10304 0.11430 0.02905
XTest Matrix: 1.00000 0.00000 0.00000 0.00000
XTest Matrix: 0.00000 1.00000 0.00000 0.00000
XTest Matrix: 0.00000 0.00000 1.00000 0.00000
XTest Matrix: 0.00000 0.00000 0.00000 1.00000
-------------------------------------------------------------