3
私は非可逆残差行列からの崩壊の影響を受けないブロック共役勾配アルゴリズムを実装しようとしています。しかし、私は無意味な結果を得ています（各反復において、Rcurrentのランクは、小さくならず、増加しなくてはなりません）。これは、Hao JiおよびYaohang Liによる論文「A breakdown-free block conjugate gradient法」に示されています。このアルゴリズムを正しく実装していますか？
これはジュリアの私の実装です：
function orth(M::Matrix)
matrixRank = rank(M)
Ufactor = svdfact(M)[:U]
return Ufactor[:,1:matrixRank]
end
function BFBCG(A::Matrix, Xcurrent::Matrix, M::Matrix, tol::Number, maxit::Number, Rcurrent::Matrix)
# initialization
#Rcurrent = B - A*Xcurrent;
Zcurrent = M*Rcurrent;
Pcurrent = orth(Zcurrent);
Xnext::Matrix = ones(size(Xcurrent))
# iterative method
for i = 0:maxit
Qcurrent = A*Pcurrent
acurrent = (Pcurrent' * Qcurrent)\(Pcurrent'*Rcurrent)
Xnext = Xcurrent+Pcurrent*acurrent
Rnext = Rcurrent-Qcurrent*acurrent
# if Residual norm of columns in Rcurrent < tol, stop
Znext = M*Rnext
bcurrent = -(Pcurrent' * Qcurrent)\ (Qcurrent'*Znext)
Pnext = orth(Znext+Pcurrent*bcurrent)
Xcurrent = Xnext
Zcurrent = Znext
Rcurrent = Rnext
Pcurrent = Pnext
@printf("\nRANK:\t%d",rank(Rcurrent))
@printf("\nNORM column1:\t%1.8f",vecnorm(Rcurrent[:,1]))
@printf("\nNORM column2:\t%1.8f\n=============",vecnorm(Rcurrent[:,2]))
end
return Xnext
end
それらの入力のための論文の結果：
A = [15 5 4 3 2 1; 5 35 9 8 7 6; 4 9 46 12 11 10; 3 8 12 50 14 13; 2 7 11 14 19 15; 1 6 10 13 15 45]
M = eye(6)
guess = rand(6,2)
R0 = [1 0.537266261211281;2 0.043775211060964;3 0.964458562037146;4 0.622317517840541;5 0.552735938776748;6 0.023323943544997]
X = BFBCG(A,guess,M,tol,9,R0)
はここ
は、アルゴリズムであります第3の反復でゼロに達するランク。
広すぎます。コンテキストのない擬似コードを表示することはそれほど役に立ちません。また、前提条件という言葉と行の検索がないと、反復の中に非減少がある理由を説明するかもしれません（しかし、それはちょうど推測です）。 – sascha
私は、アルゴリズムの目的を明確にする/コンテキストを追加するために含めることを試みた。前提条件については、論文中の恒等行列として残されている。 –
コードを少し最適化しようと思っているなら、ここで一時配列の数を減らすために使用できる事前割り振りと埋め込み演算のトンがあります。 'A_mul_B！'と 'At_mul_B！ 'を見てください。 –