OCaml中的喬列斯基因式分解

Question

我想知道是否有人可以幫我調試下面的OCaml代碼，它應該計算正定矩陣的上三角形分解。

我知道它不是很實用，很笨重，所以我提前道歉。下面給出一些原因。

無論如何這裏去！

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

但與功能bedecked用許多指數，我的頭開始旋轉。我相信這是問題所在，或者是計算。有些東西肯定是錯的;）

如果有幫助，我可以附上我正在移植它的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