我是否正確實施此算法？

我試圖實現一個塊不共軛梯度算法，不受從不可逆殘差矩陣的分解;但是我得到了無意義的結果（在每次迭代中，排名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

有趣。你的評論似乎是正確的，我會進一步調查。謝謝您的回答。（我遲早會接受） –

我是否正確實施此算法？

回答

相關問題