Cholesky factorization in OCaml

Question

I was wondering if someone could help me debug the following OCaml code that is supposed to calculate the upper triangle cholesky factorization of a positive definite matrix.

I know its not very functional and very chunky, and so I apologise in advance. I give some reasons for that below.

Anyway here goes!

let rec calc_S m1 k i = 
if k == i then 
  let float = sum_array m1.(i) in float
else
  begin
    m1.(k).(i) <- m1.(k).(i)**2.;
    calc_S m1 (k+1) i;
  end
;;

let rec calc_S1 m1 k i j len= 
if k == len then 
  let float = sum_array m1.(i) in float
else
  begin
    m1.(k).(j) <- m1.(k).(i)*. m1.(k).(j);
    calc_S1 m1 (k+1) i j len;
  end
;;             

let cholesky m1 =
  let ztol = 1.0e-5 in
  let result = zerof (dimx m1) (dimx m1) in
  let range_xdim = (dimx m1) - 1 in
  let s = ref 0.0 in
  for i=0 to range_xdim do
    begin
      s := calc_S result 0 i;
      let d = m1.(i).(i) -. !s in 
      if abs_float(d) < ztol then
        result.(i).(i) <- 0.0
      else
        if d < 0.0 then
          raise Matrix_not_positive_definite
        else
          result.(i).(i) <- sqrt d;
      for j=(i+1) to range_xdim do
        s:= calc_S1 result 0 i j range_xdim;
        if abs_float (!s) < ztol then
          s:= 0.0;
        result.(i).(j) <- (m1.(i).(j) -. !s) /. result.(i).(i)
      done
    end
  done;
  result;;

Where dimx, dimy are simple functions that return the dimensions of the matrix (2d-array), zerof generates a float matrix of zeroes with the appropriate dimensions, sum_array is a simple function for summing the elements of an array, and the exception is obviously defined previously.

For some reason it calculates the following incorrectly [edit: top-loop was polluted, correct incorrect calculation added]:

 f;;
- : float array array =
[| 1.; 0.; 0.1; 0. |]
[| 0.; 1.; 0.; 0.1 |]
[| 0.; 0.; 1.; 0. |]
[| 0.; 0.; 0.; 1. |]

# cholesky f;;
- : float array array =

 [| 1.; 0.; 0.; 0. |]
 [| 0.; 1.; 0.; 0. |]
 [| 0.; 0.; 1.; 0. |]
 [| 0.; 0.; 0.; 1. |]

which should be according to the python code:

[1.0, 0.0, 0.1, 0.0]
[0, 1.0, 0.0, 0.1]
[0, 0, 0.99498743710662, 0.0]
[0, 0, 0, 0.99498743710662]

But gets the following right:

# r;;
- : float array array = 
[| 0.1; 0. |]
[| 0.; 0.1 |]

# cholesky r;;
- : float array array = 
[| 0.316227766017; 0. |]
[| 0.; 0.316227766017 |]

but with the function bedecked with so many indices, my head is starting to swirl. I am sure this is where the problem is, or the s calculations. Something is definitely wrong ;)

If it helps, I can attach the python code I'm porting it from, since I tried to stay as close to the original, which had no recursion, and so hence no recursion in the attached OCaml version, and thus also some of the awkwardness (in OCaml).

[edit: python code attached]

 def Cholesky(self, ztol=1.0e-5):
    # Computes the upper triangular Cholesky factorization of
    # a positive definite matrix.
    res = matrix([[]])
    res.zero(self.dimx, self.dimx)

    for i in range(self.dimx):
        S = sum([(res.value[k][i])**2 for k in range(i)])
        d = self.value[i][i] - S
        if abs(d) < ztol:
            res.value[i][i] = 0.0
        else:
            if d < 0.0:
                raise ValueError, "Matrix not positive-definite"
            res.value[i][i] = sqrt(d)
        for j in range(i+1, self.dimx):
            S = sum([res.value[k][i] * res.value[k][j] for k in range(self.dimx)])
            if abs(S) < ztol:
                S = 0.0
            res.value[i][j] = (self.value[i][j] - S)/res.value[i][i]
    return res

as per the amazingly fast request :) I'm fairly certain I'm messing up the summing going on inside the statement like this:

S = sum([res.value[k][i] * res.value[k][j] for k in range(self.dimx)])

Answer 1

One definitely wrong thing is you shouldn't modify m1 in calc_S and calc_S1, these functions suppose to return the sums only. Since you didn't provide the full implementation, this is one way to fix it:

let calc_S res i = 
   let sum = ref 0. in
   for k = 0 to i-1 do
      sum := !sum +. (res.[k].[i] ** 2.)
   done;
   !sum

let calc_S1 res i dim = 
   let sum = ref 0. in
   for j = i+1 to dim-1 do
      for k = 0 to dim-1 do
        sum := !sum +. (res.[k].[i] *. res.[k].[j])
      done
   done;
   !sum

Other than this mistake, remaining parts of the code look right to me.