The Salty Economist

Things I Should Have Learned in High School
posts - 56, comments - 0, trackbacks - 0

F# Matrix Transpose, Multiply

Well, I have finally decided to write some econometrics routines in F#.  I first did this type of stuff in FORTRAN when I was a grad student, so I figured I'd try to see if I could still figure it all out.

But, first things first.

We need routines to Transpose a matrix and to Multiply two matrices.

My transpose method is this:

--------------------------------------------------------------
//
//  Let's make up some data
//

let n = 6
let k = 3

//  Create data array - n x k

let x0 = Array2D.zeroCreate<float> n k

x0.[0,0] <- 1.0
x0.[1,0] <- 1.0
x0.[2,0] <- 1.0
x0.[3,0] <- 1.0
x0.[4,0] <- 1.0
x0.[5,0] <- 1.0

x0.[0,1] <- 22.0
x0.[1,1] <- 4.0
x0.[2,1] <- 13.0
x0.[3,1] <- 3.0
x0.[4,1] <- 9.0
x0.[5,1] <- -3.0

x0.[0,2] <- 12.0
x0.[1,2] <- 14.0
x0.[2,2] <- 19.0
x0.[3,2] <- 9.0
x0.[4,2] <- 13.0
x0.[5,2] <- 23.0

//---------------------------------------------------
// Transpose a Matrix
// Takes matrix x and returns the transpose as xp
// NEWLY UPDATED VERSION 8/12/2009
//---------------------------------------------------

let rec transposeACol(x : float [,] ) (xp : float [,] ) (i : int) (j : int) =
    match j with
    | -1 -> xp
    | _ -> xp.[j, i] <- x.[i, j]
           //printfn "i: %i j: %i x: %10.5f" i j (if i >=0 && j >=0 then x.[i, j] else -1.0)
           //let userresp = Console.ReadLine()
           let xp = transposeACol x xp i (j-1)
           xp
          
let rec transpose (x : float [,] ) (xp : float [,] ) (i : int) (j : int) =
    match i with
    | -1 -> xp
    | _  -> let xp = transposeACol x xp i j 
            let xp = transpose x xp (i-1) j
            xp

// Create a shell array, k x n
let zero_xp = Array2D.zeroCreate<float> k n

//-----------------------------------------
// Create X Prime
//-----------------------------------------
//
// Handle the fact that matrices are zero based: use n-1, k-1 as starting points
//

let xp = transpose x0 zero_xp (n-1) (k-1)

//
//Print xp
//
// Print it out
for i = 0 to k - 1 do
     printfn "%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f" xp.[i,0] xp.[i,1] xp.[i,2] xp.[i,3] xp.[i,4] xp.[i,5]

--------------------------------------------------------------

The method transpose is a recursive method that starts at the top for both n and k and works down to zero, flipping the elements as it goes.  You send it a blank array and it returns the transpose.  I know I could have done this with a couple of for..do loops, but I wanted to practice writing recursive functions in F#.

Note the Array function Array2D.zeroCreate<float> k n.  This function creates a two dimensional array (floating point) with k rows and n columns.  It is prefilled with zeros.  The array is zero-based.  There is also a method Array.create that allows you to specify the initial value (instead of zero).

 So we run the code, and here's what we get:

      1.00       1.00       1.00       1.00       1.00       1.00
     22.00       4.00      13.00       3.00       9.00      -3.00
     12.00      14.00      19.00       9.00      13.00      23.00

The next thing we need to code is a matrix multiplication routine.

If you aren't familiar with matrix multiplication, here is Wikipedia's page: http://en.wikipedia.org/wiki/Matrix_multiplication

The idea here is that you have a matrix that is dimensioned n x k and you want to multiply it bay a matrix of dimension k x m.  The two k's must be the same.  The resulting matrix is dimensioned n x m.

So, for example, if we have two matrices xp (2 x 4) and x (4 x 2), the resulting matrix (xpx) will have be dimensioned 2 x 2 and the i,jth item will be:

xpx[i, j] = (xp[i,0] * x[0,j]) + (xp[i,1] * x[1,j]) + (xp[i,2] * x[2,j]) + (xp[i,3] * x[3,j])

So here's my code:

-----------------------------------------------------------------------------

let accum (x1 : float) (x2 : float) (ac : float) =
    ac + (x1 * x2)

let rec xply (xp : float [,] ) (x : float [,] ) (j0 : int) (j1 : int) (i : int) (cell_accum : float) =
    match i with
    | 0 -> accum xp.[j0, i] x.[i, j1]  cell_accum
    | _  -> xply xp x j0 j1 (i-1) (accum xp.[j0, i] x.[i, j1]  cell_accum)
       

//------------------------------------------------
// Multiply two Matrices:  [n,k] x [k,m] = [k, k]
//------------------------------------------------
// j0 is the strating iterator for the row
// j1 is the starting iterator for the column
// k is the size of k -> stays fixed
// m is the column count of the second array
// It is used to reset j1 as we work down the rows
let rec xMultiply (xp : float [,] ) (x : float [,] ) (xpx : float [,] )(j0 : int) (j1 : int) (k : int) (m : int) =

    match j0 with
    | -1 -> xpx
    | _ -> match j1 with
        | -1 -> xMultiply xp x xpx (j0-1) m k m
        | _  -> let tmp = xply xp x j0 j1 k 0.0
                xpx.[j0, j1] <- tmp
                printfn "xpx j0: %i j1 %i = %10.5f" j0 j1 xpx.[j0, j1]
                xMultiply xp x xpx j0 (j1-1) k m


//------------------------------------------------
// Multiply two Matrics:  [n,k] x [k,m] = [n, m]
//------------------------------------------------
let matMult (xp : float [,] ) (x : float [,] ) (xpx : float [,] ) =
    let n = xp.GetUpperBound 0
    let k = xp.GetUpperBound 1
    let m = x.GetUpperBound 1
    xMultiply xp x xpx n m k m


//-----------------------------------------
// Create X Prime X
//-----------------------------------------

let zero_xpx = Array2D.zeroCreate<float> k k
//let xpx = matMult xp x0 zero_xpx


-----------------------------------------------------------------------------

The routine goes as follows:

The function accum, just takes two values, multiplies them, and adds then to an accumulator value (ac) that is passed in.  Obviously the accumulator value would be initialized at 0.0

The function xply iterates over a row (row j0)  in the left matrix and a column (j1) in the right matrix. calling accum to accumulate the result.

Note that in the xply function, this line

    | _  -> xply xp x j0 j1 (i-1) (accum xp.[j0, i] x.[i, j1]  cell_accum)

could be written in two lines as:

    | _  -> let tmp = accum xp.[j0, i] x.[i, j1]  cell_accum
            xply xp x j0 j1 (i-1) tmp

but the former is somewhat more 'elegant.'  The return value from xply is the value for one cell (j0, j1) the resultant matrix from the multiplication operation.

The function xMultiply is a recursive function that iterates over all the cells in the matrix, calling xply each time to fill in a cell.  It has two matching levels, one outer level that iterates over the rows; one inner level that iterates over the rows in a given column.

The top level function  matMult just takes the two input matrices (left & right) and a empty output matrix as inputs.  Note the method GetUpperBound i.  This method returns the upper bound for a matrix for the dimension you pass in.  In a 2D array, The 0th dimension is for the rows, the 1st dimension is for the columns.  This is a very handy function because it means the user does not have to pass in the dimensions of the matrices to be multiplied.

The out of the matrix multiplication is:

  6.00      48.00      90.00
 48.00     768.00     642.00
 90.00     642.00    1480.00

 

 

Print | posted on Wednesday, July 22, 2009 2:17 PM |

Powered by:
Powered By Subtext Powered By ASP.NET