MINRES法を実装してみた。

MINRES法とは？

MINRES法（最小残差法）とは、krylov部分空間法と言われる手法の一つです。
以下のように、正方対称行列の正規方程式を扱う計算です。

Ax = b,\;\; A \in \mathbb{R}^{n \times n} :: \text{Symmetric}

実装

JuliaでのシンプルなMINRES法の実装は以下の通りです。
この実装はちゃんと収束しますが、かなりNaiveです。
また、実際には、k_maxまでで収束しない場合に解の改善を行うようにすべきですが、この実装は短いコードで作りたかったので削減しています。

参考
自ページ

using LinearAlgebra

function simple_minres(
      A::AbstractMatrix
    , b::AbstractMatrix
    ; k_max = 10
    , tol   = 1e-7
    , n     = size(A, 2)
)
    pprev, prev, cur, next = 1, 2, 3, 4

    #    k-2           k-1           k             k+1
    #    -1            0             1             2           (k=1)
    x = [zeros(n, 1) , rand(n, 1)  , zeros(n, 1) , zeros(n, 1) ]
    v = [zeros(n, 1) , zeros(n, 1) , b-A*x[prev] , zeros(n, 1) ]
    β = [0.0         , 0           , norm(v[cur]), 0           ]
    γ = [0.0         , 1           , 1           , 0           ]
    σ = [0.0         , 0           , 0           , 0           ]
    ω = [zeros(n, 1) , zeros(n, 1) , zeros(n, 1) , zeros(n, 1) ]
    α = [0.0         , 0           , 0           , 0           ]

    rknorm = η = β[cur]
    
    for k = 1:k_max
        v[cur]  = (1 / β[cur])*v[cur]
        Av      = A*v[cur]
        α[cur]  = ((v[cur]')*Av)[]
        v[next] = Av - α[cur]*v[cur] - β[cur]*v[prev]
        β[next] = norm(v[next])
        δ       = γ[cur]*α[cur] - γ[prev]*σ[cur]*β[cur]
        ρ1      = √(δ^2 + β[next]^2)
        ρ2      = σ[cur]*α[cur] + γ[prev]*γ[cur]*β[cur]
        ρ3      = σ[prev]*β[cur]
        γ[next] = δ / ρ1
        σ[next] = β[next] / ρ1
        ω[cur]  = (v[cur] - ρ2*ω[prev] - ρ3*ω[pprev]) / ρ1
        x[cur]  = x[prev] + γ[next]*η*ω[cur]
        rknorm  = abs(σ[next])*rknorm
        η       = - σ[next]*η

        if (rknorm < tol) println("n_ittr = $k"); break; end

        pprev, prev, cur, next = prev, cur, next, pprev
    end
    return x[prev]
end
;

以下のようにすると実際に収束することが確認できる。

A = Symmetric(rand(10, 10)); b = rand(10, 1); 
x = simple_minres(A, b, k_max=20)

norm(A*x - b)