4

I wrote a code to perform some basic Matrix calculation using F#. I would like to know if there are some possible improvements on this code in order to decrease the calculation time. Indeed the operations performed are quite basic (multiplication of 2 matrices and transpose mainly) however sizes of the matrix is high (around 10000 * 100000) leading to a huge calculation duration (few hours).

My questions/remarks are the following :

  1. Is there any way to improve the following code? There are many "for loops" which may lead to a serious slow down of the algorithm, but I don't know how to avoid these "for loops".
  2. I created some initial Matrix with initial values at 0 and, in a second time, filled their elements with results. Maybe it is possible to avoid the first step of initialisation.

Here is the algorithm :

// I use the #time function to calculate the calculation duration of the algorithm
#time

#r "Microsoft.Office.Interop.Excel"
#r "FSharp.PowerPack.dll"

open System
open System.IO

open Microsoft.FSharp.Math
open System.Collections.Generic

// Algorithm
let matrixCalculation (matA : matrix) (matB : matrix) (matC : matrix) =  

    // First step : Renamed the matrix A and B size to initialize the matrix "matrixCalcul" 
    let nbrOfElementsA = matA.NumRows
    let nbrOfElementsB = matB.NumRows
    let nbrOfCaracteristicsA = matA.NumCols
    let nbrOfCaracteristicsB = matB.NumCols

    // Second step : MatB has to be transposed 
    let tmatB = matB.Transpose

    // Initialisation of the final output named matrixCalcul. A weighted vector is also initialised 
    let mutable matrixCalcul = Matrix.create (nbrOfElementsA + 1) (nbrOfElementsB + 1) 0.            
    let mutable weightedVector = Matrix.create nbrOfCaracteristicsA 1 0.                   

    // The first column of matA and matB represents IDs, and are "copy/past" in matrixCalcul's first colum and first row respectively
    matrixCalcul.[1.. ,0..0] <- matA.[0..,0..0]
    matrixCalcul.[0..0,1 ..] <- matB.[0..,0..0].Transpose

    // Then the core of the matrix named "matrixCalcul" can be calculated
    for j = 0 to (nbrOfElementsB - 1) do
        weightedVector <- matC * tmatB.[1..(nbrOfCaracteristicsB - 1),0..(nbrOfElementsB-1)].Columns(j,1)                       
        for i = 0 to (nbrOfElementsA - 1) do
            let mutable acc  = matA.[0..(nbrOfElementsA - 1),1..(nbrOfCaracteristicsA-1)].Rows(i,1) * weightedVector                
            matrixCalcul.[i+1,j+1] <- (acc.[0,0])
    matrixCalcul


// Two matrix generators (one for matA and matB and another one for matC)

let matrixTestGeneratorAandB nbrOfElements nbrOfCaracteristics = 
    let matrixTestGeneratedAandB = Matrix.create nbrOfElements nbrOfCaracteristics 0.
                                   |> Matrix.mapi (fun i j value -> if j = 0 then float(i + 1) elif j % 2 = 0 then 1. else 0.)
    matrixTestGeneratedAandB

let matrixTestGeneratorC nbrOfElements nbrOfCaracteristics = 
    let matrixTestGeneratedC = Matrix.create nbrOfElements nbrOfCaracteristics 0.
                               |> Matrix.mapi (fun i j value -> if j = 0 then 0. elif j % 2 = 0 then 1. else 0.)
    matrixTestGeneratedC


// Generation of matrixA, matrixB and matrixC

let matrixA = matrixTestGeneratorAandB 100 179

let matrixB = matrixTestGeneratorAandB 100 639

let matrixC = matrixTestGeneratorC 178 638

// Calculation 
matrixCalculation matrixA matrixB matrixC

Basically the calculation duration amounts to circa 2 seconds, but if you change the number of matrixA and matrixB up to 10000, it can take hour. Just for information, in my algorithm, matrixC's size will remain constant, only matrix A and B can have an increasing number of rows.

If you have any idea of improvement, I take it.

5
  • OK no problem I am on it
    – fabco63
    Dec 19, 2011 at 18:30
  • Could you provide a test script for benchmarking and explain more clearly what your code is actually doing?
    – pad
    Dec 19, 2011 at 18:31
  • My gut reaction would be to remove the matrix creation from the inner most loop.
    – gradbot
    Dec 19, 2011 at 18:41
  • 1
    I'm thinking this question is better asked on Code Review. Dec 19, 2011 at 19:02
  • Have you tried using some matrix library instead?
    – svick
    Dec 19, 2011 at 19:29

1 Answer 1

9

From your code, it's quite difficult to understand what you're trying to achieve. I think you mean to compute a matrix d[0..m, 0..n] as follows:

  +---------+-------------------------+
  | 0.0     | b00 b10 ......  b(n-1)0 |
  +---------+-------------------------+
  | a00     | d11 d12 ......  d1n     |
  | a10     | d21 d22 ......  d2n     |
  | ...     | ... ... ......  ...     |
  | ...     | ... ... ......  ...     |
  | ...     | ... ... ......  ...     |
  | a(m-1)0 | dm1 dm2 ......  dmn     |
  +---------+-------------------------+

where the core part (the inner matrix d[1..m, 1..n]) is a multiplication of three matrices matA1 (matA after trimmed the first columns), matC and matB1 (matB after trimmed the first column and transposed).

To understand matrix operation, a good way to go is reasoning on matrix size. Let ra, ca, rb, cb, rc and cc denote the number of rows and columns in matA, matB and matC respectively. The multiplication is among three matrices of size ra x (ca-1), rc x ccand (cb-1) x rb; this only makes sense if rc = ca-1 and cc = cb-1. We have the resulting matrix d of size (ra+1) x (rb+1).

Here is my attempt without using any for loop:

let calculate (matA : matrix) (matB : matrix) (matC : matrix) = 
    let ra = matA.NumRows
    let ca = matA.NumCols
    let rb = matB.NumRows
    let cb = matB.NumCols
    let matrixCalcul = Matrix.zero (ra+1) (rb+1)

    matrixCalcul.[1.., 0..0] <- matA.[0.., 0..0]
    matrixCalcul.[0..0, 1..] <- matB.[0.., 0..0].Transpose

    matrixCalcul.[1.., 1..] <- (matA.Columns(1, ca-1) * matC) * matB.Columns(1, cb-1).Transpose
    matrixCalcul

I have tested with matA, matB and matC of size 200x279, 200x1279 and 278x1238 respectively. Two versions yield the same result and my function is 40x faster than the original one. There are many reasons for that, but in general a vectorized version has much better performance when it comes to matrix calculation.

2
  • 2
    Deserves +10 just for the effort!
    – Benjol
    Dec 20, 2011 at 6:32
  • Thank you pad, it's perfectly working. Moreover my question was not very easy to decript... Thank you again
    – fabco63
    Dec 20, 2011 at 10:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.