Ok, so here is my F# code to perform a linear decomposition of a matrix:
First, let's make up some data:
open System
open System.Collections.Generic
open System.IO
open Xunit
open MathMod
open Microsoft.FSharp.Math
let X = Array2D.zeroCreate<float> 4 4
X.[0,0] <- 4.0
X.[1,0] <- 2.0
X.[2,0] <- 3.0
X.[3,0] <- 9.0
X.[0,1] <- 7.0
X.[1,1] <- 11.0
X.[2,1] <- 3.0
X.[3,1] <- -2.0
X.[0,2] <- -8.0
X.[1,2] <- 2.0
X.[2,2] <- -17.0
X.[3,2] <- 4.0
X.[0,3] <- 6.0
X.[1,3] <- 9.0
X.[2,3] <- 10.0
X.[3,3] <- 5.0
X is just a 4 x 4 array.
We next create a function that takes the X array as an input and does the decomposition, returning both L and U as the Lower and Upper matrices.
let deComp (X : float [,]) (k : int) =
//-------------------------------------------------------
//Instantiate the Upper and Lower Diagonal Arrays
//-------------------------------------------------------
let L = Array2D.zeroCreate<float> (k+1) (k+1)
let U = Array2D.zeroCreate<float> (k+1) (k+1)
//Initialize the Diagonal in L
for j = 0 to k do
L.[j, j] <- 1.0
//Work over the columns one-by-one
let L, U, ret = doColumns L U X 0 k
ret, L, U
let liDecomp (X : float [,]) =
let k, ret = isSquare X
if ret = -1 then
deComp X k
else
0, X, X
The function liDecomp is the main function. The function first checks to make sure that X is in fact square (isSquare is a function that does the checking). If the X array is square, then deComp performs the decomposition, returning L, U and a return code (ret) to indicate success or failure. These value are returned in a tuple.
The function isSquare just uses the method .GetUpperbound for both dimensions of X and makes sure they are the same:
//------------------------------------------
//Check to make sure array is square
//------------------------------------------
let isSquare (X : float [,]) =
let k0 = X.GetUpperBound 0
let k1 = X.GetUpperBound 1
let ret =
if k0 <> k1 then
0
else
-1
k0, ret
The function returns the dimension of the matrix, as well 0 or -1 depending if the matrix is square or not. The values are returned in a tuple.
If the matrix is square, the liDecomp function calls the deComp function.
Let's now look at it:
The deComp function first instantiates the L & U arrays with dimensions k+1, where k+1 is the number of elements for a matrix with an upper bound of k. We then initialize the main diagonal of L to be all ones.
We then call a function doColumns to iterate over columns of the X matrix.
Here is the doColumns function:
let rec doColumns (L : float [,]) (U : float [,]) (X : float [,]) (j : int) (k : int) =
match j with
| j when j > k -> L, U, 0
| _ -> let L, U, ret = doColumn L U X 0 j k 0
if ret = 0 then
doColumns L U X (j+1) k
else
L, U, ret
The doColumns function is a recursive function that calls itself, starting with the first column j (initialized equal to 0), and ending when j reaches the total number of columns k. The column number (j) is incremented in each successive call. For each value of j, the function calls another function called doColumn. When j is greater than k, the recursion stops and we return L, U, and a return code (indicating success or failure) in a tuple. The return code is required to guard against a divide by zero occurrence in the doColumn function
The function doColumn is:
let rec doColumn (L : float [,]) (U : float [,]) (X : float [,]) (i : int) (j : int) (k : int) (ret : int) =
match i with
| i when ret <> 0 -> L, U, ret
| i when i > k -> L, U, ret
| _ -> let ret0 =
if i > j then
if checkZero U.[j,j] = 0 then
L.[i,j] <- doCell L U X i j (j-1) U.[j,j]
0
else
-2
else
if checkZero L.[j,j] = 0 then
U.[i,j] <- doCell L U X i j (i-1) L.[j,j]
0
else
-2
let L, U, xRet = doColumn L U X (i+1) j k ret0
L, U, xRet
The function takes as input parameters the matrices L, U and X, the row index i (which starts at zero and increments by one as the function is called recursively, the column index j, and the maximum dimension of the array k. When i hits k, the recursion stops. The return code ret is also passed in with an initial value of zero. If we hit a divide by zero occurrence, the value of ret will be set at -2 and the recursion will quit.
The match statement has three parts:
1. When ret <> 0, we have a divide by zero, so we just quit and return whatever is in the matrices Land U and the return code (should be -2)
2. When i > k, then return the current state of the matrices Land U and the current return code
3. Otherwise, compute the value for the cell. Before we calculate the value for a cell, we check to see if the denominator is approximately zero.
We use the function checkZero:
let checkZero (v: float) =
if Math.Abs(v) < 0.000001 then
-1
else
0
to accomplish this.
We need to first check to see if we are working above or below the main diagonal. If the row i is equal to or less than the column j, we must be working on the values in the U (upper diagonal) array; for i greater than j we must be working on the values in the L (lower diagonal) array. For a nonzero denominator, we call the function doCell:
let doCell (L : float [,]) (U : float [,]) (X : float [,]) (i : int) (j : int) (iMx : int) (denom : float) =
let xSum = USum L U 0 i j iMx 0.0
let aVal = (X.[i,j] - xSum) / denom
aVal
This function does all the math work and returns the value for the cell in the variable aVal. The input parameters here are the matrices L, U and X, the row and column indices i & j, as well as the stopping point iMx. For the L matrix, the stopping point will be j-1; for the U matrix, the stopping point will be i-1. The value of aVal:
let aVal = (X.[i,j] - xSum) / denom
is the general equation for calculating the value of a cell (see my previous post for the logic behind this). I use the temporary variable xSum to hold the sum of the products used in the calculation. I could have substituted out this variable, but I thought the function would be easier to read by using two lines rather than one. I could just have easily written:
let aVal = (X.[i,j] - (USum L U 0 i j iMx 0.0)) / denom
instead.
The USum function sums up the cross products between the two arrays, L and U.
let rec USum (L : float [,]) (U : float [,]) (ii : int) (i : int) (j : int) (iMx : int) (v: float) =
// start at 0 and go to iMx
// ii is the current index
match ii with
| ii when ii > iMx -> 0.0
| ii when ii = iMx -> v + L.[i,ii] * U.[ii,j]
| _ -> let tmp = v + L.[i,ii] * U.[ii,j]
USum L U (ii+1) i j iMx tmp
The USum function takes as inputs the matrices L and U, a row/column index ii, the row and column indices i & j, the stopping point iMx, and an accumulator variable v, initialized with a value of zero. The index ii runs from zero the maximum value and accumulates the products of the elements of the row i and the column j:
L.[i,ii] * U.[ii,j]
Note that the index ii is used as the column index for the L array and as the row index for the U array. We also use a placeholder variable tmp to store the result before passing it on to the next recursive call. Again, we could have written this as:
USum L U (ii+1) i j iMx (v + L.[i,ii] * U.[ii,j])
but it would not be quite as readable.
At last, we are finished with the gory details. One last thing I want to do is write a function to test it out:
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]
ret, L, U
This function does the decomposition of X into the matrices L and U and then multiples them to make sure their product does indeed equal out original matrix. I sure hope I haven't messed up any of the row & column indices, they were quite a pain to keep track of. So 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
So, what's next? Calculate the inverse of the matrix from the L & U matrices, of course. That's the topic of my next post.
-----------------------------------------------------------------------------
Below is the complete sorce code
-----------------------------------------------------------------------------
#light
open System
open System.Collections.Generic
open System.IO
open Xunit
open MathMod
open Microsoft.FSharp.Math
let X = Array2D.zeroCreate<float> 4 4
X.[0,0] <- 4.0
X.[1,0] <- 2.0
X.[2,0] <- 3.0
X.[3,0] <- 9.0
X.[0,1] <- 7.0
X.[1,1] <- 11.0
X.[2,1] <- 3.0
X.[3,1] <- -2.0
X.[0,2] <- -8.0
X.[1,2] <- 2.0
X.[2,2] <- -17.0
X.[3,2] <- 4.0
X.[0,3] <- 6.0
X.[1,3] <- 9.0
X.[2,3] <- 10.0
X.[3,3] <- 5.0
//------------------------------------------
//Check to make sure array is square
//------------------------------------------
let isSquare (X : float [,]) =
let k0 = X.GetUpperBound 0
let k1 = X.GetUpperBound 1
let ret =
if k0 <> k1 then
0
else
-1
k0, ret
//------------------------------------------
//Accumulate cross products
//------------------------------------------
let rec USum (L : float [,]) (U : float [,]) (ii : int) (i : int) (j : int) (iMx : int) (v: float) =
// start at 0 and go to iMx
// ii is the current index
match ii with
| ii when ii > iMx -> 0.0
| ii when ii = iMx -> v + L.[i,ii] * U.[ii,j]
| _ -> let tmp = v + L.[i,ii] * U.[ii,j]
USum L U (ii+1) i j iMx tmp
//------------------------------------------
//Compute a Cell Value
//------------------------------------------
let doCell (L : float [,]) (U : float [,]) (X : float [,]) (i : int) (j : int) (iMx : int) (denom : float) =
//let iMx = j-1
let xSum = USum L U 0 i j iMx 0.0
//let aVal = (X.[i,j] - xSum) / U.[j,j]
let aVal = (X.[i,j] - xSum) / denom
aVal
//------------------------------------------
//Check for Zero
//The value 0.000001 is arbitrary
//------------------------------------------
let checkZero (v: float) =
if Math.Abs(v) < 0.000001 then
-1
else
0
//-----------------------------------------------------------
//Calculate the Values for a Column
//Work on the Lower or Upper Matrix depending on the row i
//J is the current column
//k is the maximum array dimension
//-----------------------------------------------------------
let rec doColumn (L : float [,]) (U : float [,]) (X : float [,]) (i : int) (j : int) (k : int) (ret: int) =
match i with
| i when ret <> 0 -> L, U, ret
| i when i > k -> L, U, ret
| _ -> let ret0 =
if i > j then
if checkZero U.[j,j] = 0 then
L.[i,j] <- doCell L U X i j (j-1) U.[j,j]
0
else
-2
else
if checkZero L.[j,j] = 0 then
U.[i,j] <- doCell L U X i j (i-1) L.[j,j]
0
else
-2
let L, U, xRet = doColumn L U X (i+1) j k ret0
L, U, xRet
//------------------------------------------
//Calculate the Array values, Column by Column
//------------------------------------------
let rec doColumns (L : float [,]) (U : float [,]) (X : float [,]) (j : int) (k : int) =
match j with
| j when j > k -> L, U, 0
| _ -> let L, U, ret = doColumn L U X 0 j k 0
if ret = 0 then
doColumns L U X (j+1) k
else
L, U, ret
let deComp (X : float [,]) (k : int) =
//-------------------------------------------------------
//Instantiate the Upper and Lower Diagonal Arrays
//-------------------------------------------------------
let L = Array2D.zeroCreate<float> (k+1) (k+1)
let U = Array2D.zeroCreate<float> (k+1) (k+1)
//Initialize the Diagonal in L
for j = 0 to k do
L.[j, j] <- 1.0
//Work over the columns one-by-one
let L, U, ret = doColumns L U X 0 k
//let ret = -1
ret, L, U
let liDecomp (X : float [,]) =
let k, ret = isSquare X
if ret = -1 then
deComp X k
else
-1, X, X
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" L.[i,0] L.[i,1] L.[i,2]
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 " U.[i,0] U.[i,1] U.[i,2]
printfn "U Matrix: %10.2f %10.2f %10.2f %10.2f " U.[i,0] U.[i,1] U.[i,2] U.[i,3]
let userresp = Console.ReadLine()
ret, L, U