OCaml 中的 Cholesky 分解

Question

我想知道是否有人可以帮助我调试以下 OCaml 代码，该代码应该计算正定矩阵的上三角 Cholesky 分解。

我知道它不是很实用而且很厚实，所以我提前道歉。我在下面给出了一些原因。

不管怎么样！

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;;

其中 dimx、dimy 是返回矩阵维度（二维数组）的简单函数，zerof 生成具有适当维度的零浮点矩阵，sum_array 是用于对数组元素求和的简单函数，而异常显然是之前定义的。

由于某种原因，它错误地计算了以下内容[编辑：顶部循环被污染，添加了正确的错误计算]：

 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. |]

应该根据python代码：

[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]

但得到以下权利：

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

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

但是这个函数有这么多索引，我的头开始旋转。我确定这是问题所在，或者是 s 计算。肯定有问题;)

如果有帮助，我可以附加要从中移植它的 python 代码，因为我试图保持与原始代码尽可能接近，它没有递归，因此在附加的 OCaml 版本中没有递归，因此也一些尴尬（在 OCaml 中）。

[编辑：附上python代码]

 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

根据惊人的快速请求 :) 我相当肯定我搞砸了这样的语句中的求和：

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

Answer 1

一个绝对错误的事情是你不应该在calc_S 和calc_S1 中修改m1，这些函数假设只返回总和。由于您没有提供完整的实现，这是修复它的一种方法：

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

除了这个错误，其余部分代码对我来说都是正确的。

Answer 2

你在这里错过了一个重要的点。只能为 对称（或厄米特）正定矩阵计算 Cholesky 分解。

你失败的例子中的矩阵不是对称的，所以你不能在这里期待任何东西。您的工作示例中的矩阵 是对称的，并且可以正常工作。

如果您需要对非对称矩阵进行分解，请考虑 QR decomposition