我是否正确实施此算法？

我试图实现一个块不共轭梯度算法，不受从不可逆残差矩阵的分解;但是我得到了无意义的结果（在每次迭代中，排名Rcurrent应该变小，而不是增加）。它被提供在Hao Ji和Yaohang Li的论文“A breakdown-free block conjugate gradient method”中。我是否正确实施此算法？

这里的算法：

这是我在茱莉亚实现：

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)

是等级在第三次迭代中达到零。

来源

2017-09-14 K. Gkinis

太宽了。在没有上下文的情况下显示一些伪代码并不有用。另外：单词preconditioner和没有线搜索可能会解释为什么在某些迭代中非减少（但这只是猜测）。 – sascha

我试图包括算法的目标，以澄清/添加上下文。至于预处理器，它作为文件中的单位矩阵。 –

如果您希望优化代码，那么可以使用大量的预分配和就地操作来减少临时数组的数量。查找'A_mul_B！'和'At_mul_B！'。 –

该算法起作用，并且在第三次迭代中排名变为零。问题在于数字不准确，这将使任何矩阵完全排名。为了获得更好的结果，请使用rank(Rcurrent, tol)而不是rank(Rcurrent)这是一个考虑到容差的版本。之后，至少在我的机器上，等级降到零。

julia> X = BFBCG(A,guess,M,tol,9,R0) 

RANK: 2 
NORM column1: 1.78951939 
NORM column2: 0.41155080 
============= 
RANK: 2 
NORM column1: 0.97949620 
NORM column2: 0.16170799 
============= 
RANK: 0 
NORM column1: 0.00000000 
NORM column2: 0.00000000 
============= 
RANK: 0 
NORM column1: 0.00000000 
NORM column2: 0.00000000 
=============

来源

2017-09-14 17:12:50

有趣。你的评论似乎是正确的，我会进一步调查。谢谢您的回答。（我迟早会接受） –

我是否正确实施此算法？

回答

相关问题