АЛГОРИТМ ХОЛЕЦКОГО

Алгоритм Холецкого. Алгоритм нижнего треугольного разложения Холецкого P=LL' состоит в том, что последовательно обрабатываются диагональный элемент и убывающий по высоте столбец нижнего треугольного блока матрицы P, результат L размещается по месту матрицы P.


A=rands(4), P={A'*A}, 

L=ch(P), mesh(L), tr(L), P2={L'*L}, N={P-P2}, N=norm(N), N=?,  

function: ch(P),
var i j k s,
L={P}, L[0][0]=sqrt(P[0][0]), s=rowm(L), 
for i=1:s, L[0][i]=L[0][i]/L[0][0], L[i][0]=0, end, 
for k=1:s, 
for j=0:k-1, L[k][k]=L[k][k]-L[j][k]*L[j][k], end,
L[k][k]=sqrt(L[k][k]),
for i=k+1:s, L[i][k]=0,
for j=0:k-1, L[k][i]=L[k][i]-L[j][i]*L[j][k], end,
L[k][i]=L[k][i]/L[k][k], 
end, end,
return L,
end,

Разместим дополнительной нижней строкой матрицы P вектор правой части r, алгоритм Холецкого подготовит новый вектор правой части L-1r. Убедимся в этом, разместив этот же вектор, вычисленный по формуле, еще ниже, самой нижней строкой.


A=rands(4), P={A'*A}, r=[1 1 1 1], 

Pr={P}, Pr={[P r]}, tr(Pr), L2=ch2(Pr), 
% СРАВНЕНИЕ РЕШЕНИЯ, 
L=chol(P), r={L\r}, 
tr(L2), L2={[L2 r]}, tr(L2), mesh(L2), 

function: ch2(P),
var i j k s,
L={P}, L[0][0]=sqrt(P[0][0]), s=rowm(L), 
for i=1:s, L[0][i]=L[0][i]/L[0][0], if i<>s, L[i][0]=0, end, end, 
for k=1:s-1, 
for j=0:k-1, L[k][k]=L[k][k]-L[j][k]*L[j][k], end,
L[k][k]=sqrt(L[k][k]),
for i=k+1:s, if i<>s, L[i][k]=0, end,
for j=0:k-1, L[k][i]=L[k][i]-L[j][i]*L[j][k], end,
L[k][i]=L[k][i]/L[k][k], 
end, end,
return L,
end,

Таким образом, видно, что алгоритм Холецкого попутно едва ли не решает систему нормальных уравнений, сведенную к L'x=L-1r.



Rambler's Top100