Мне было интересно, может ли кто-нибудь помочь мне отладить следующий код 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 - это простые функции, которые возвращают размеры матрицы (2d-массив), 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 |]
но с этой функцией, имеющей так много индексов, моя голова начинает кружиться. Я уверен, что это где проблема, или с расчетами. Что-то определенно не так;)
Если это поможет, я могу прикрепить код 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)])