The Mad Economist

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