I was wondering if anyone could help me debug the following OCaml code, which should calculate the cholez factorization of the upper triangle of a positive definite matrix.
I know that it is not very functional and very stocky, so I apologize in advance. I have given several reasons for this below.
Anyway, here it goes!
let rec calc_S m1 ki = 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 kij 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) ij 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 ij 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 floating-point matrix of zeros with the corresponding sizes, sum_array is a simple function for summing the elements of the array, and the exception is obviously defined earlier.
For some reason, it incorrectly calculates the following: [edit: top-loop was dirty, fixed 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. |]
which should match 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 has the following right:
# r;; - : float array array = [| 0.1; 0. |] [| 0.; 0.1 |]
but with a function adorned with a multitude of indexes, my head begins to swirl. I am sure this is a problem, or s calculations. Something is definitely wrong;)
If this helps, I can attach the python code from which I am porting it, since I tried to stay as close as possible to the original, which did not have recursion, and therefore there was no recursion in the attached version of OCaml, and therefore some of the awkwardness ( at 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
according to an amazingly quick request :) I'm pretty sure I was messing around with the summation inside the statement like this:
S = sum([res.value[k][i] * res.value[k][j] for k in range(self.dimx)])