我是否正确地实现了这个算法？

Question

我正在尝试实现一个块共轭梯度算法，它不受不可逆残差矩阵的分解；但是我得到了荒谬的结果（在每次迭代中，Rcurrent 的排名应该越来越小，而不是增加）。它在 Hao Ji 和 Yaohang Li 的论文 "A breakdown-free block conjugate gradient method" 中提出。

算法如下：

这是我在 Julia 中的实现：

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)

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

Answer 1

算法有效，并且在第三次迭代中排名变为零。问题是数值不准确，这会使任何矩阵完全排名。为了获得更好的结果，请使用 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
=============

我是否正确地实现了这个算法？

Do I implement this algorithm correctly?

algorithm

linear-algebra

julia