Lanczos法とBlock-Lanczos法の実装

Lanczos法（Lanczos algorithm）

Lanczos法とは，クリロフ部分空間法の一種です。
対称行列を以下のように三重対角化します。

$$
A = P D P^\mathrm{T}, AP = PD
$$

なお、固有値の計算などで一般的に３重対角化に用いられる手法はLanczos法ではなく、Householder法です。したがって、Lanczos法はHouseholderが不安定な場合の代替手段などに用いられているらしいです。

Lanczos法のアルゴリズム

$P = [p_1, p_2, \dots, p_n]$ とすると、以下のように示されます。

A[p_1, p_2, \dots, p_n] = 
[p_1, p_2, \dots, p_n]
\begin{bmatrix}
  \alpha_1 & \beta_1  &         & & & O \\
  \beta_1  & \alpha_2 & \beta_2 & & & \\
           & \beta_2  & \alpha_3& & & \\
           &          &         & \dots & & \\
           &          &         & & \alpha_{n-1} & \beta_{n-1} \\
  O        &          &         & & \beta_{n-1}  & \alpha_{n}  
\end{bmatrix}

これらの両辺を比較すると、以下の式が導かれます。

\begin{matrix}
  A p_1 = \alpha_1 p_1 + \beta_1 p_2 \\
  A p_2 = \beta_1 p_1 + \alpha_2 p_2 + \beta_2 p_3 \\
  A p_3 = \beta_2 p_2 + \alpha_3 p_3 + \beta_3 p_4 \\
  \dots \\
  A p_{n-1} = \beta_{n-2} p_{n-2} + \alpha_{n-1} p_{n-1} + \beta_{n-1} p_n \\
  A p_n = \beta_{n-1} p_{n-1} + \alpha_n p_n
\end{matrix}

これらの式に、 $p_1, p_2, p_3, \dots, p_{n-1}, p_n$ を左から転置して作用すると、以下の式が得られる。

$$
p_k^* A p_k = \alpha_k
$$

ここで、$u_1$ を適当な大きさ $1$ のベクトルとして、初期値として計算していくと、以下の式によって順次 $p_k$ が求められる。

$$
p_{k+1} = A p_k - \beta_{k-1} p_{k-1} - \alpha_k p_k
$$

直交条件 $p_k^* p_k = 1$ を満たすように、 $\beta_k = | v_{k+1} |$ とする。

Lanczos法の実装

Lanczos法の実装は以下の通りです。※ωとwの違いに気を付けてください。

using LinearAlgebra

function simple_lanczos(
      A::AbstractMatrix
    ; n = size(A, 2)
    , m = n
)
    # 確保
    P = zeros(n, n)
    α = zeros(n)
    β = zeros(n)

    # 初期化
    P[:, 1] .= (x -> x / norm(x))(rand(n, 1))
    w    = A * P[:, 1]
    α[1] = w' * P[:, 1]
    ω    = w - α[1] * P[:, 1]

    # 反復
    for k = 2:m
        β[k] = norm(ω)
        if β[k] != 0   P[:,k] = ω / β[k]  end
        w    = A * P[:, k] # tmp vector
        α[k] = w' * P[:, k]
        ω    = w - α[k] * P[:, k] - β[k] * P[:, k-1]
    end

    # 行列化
    D = Tridiagonal(β[begin+1:end], α, β[begin+1:end])

    return P, D
end
;

以下のコードを実行すると、分解がちゃんとできていることが分かります。

A = Symmetric(rand(10, 10))
P, D = simple_lanczos(A)

norm(A*P - P*D)

Block-Lanczos法のアルゴリズム

$V_i \in \mathbb{R}^{n \times p}, V_j^\mathrm{T} V_j = I_p$ として、$B_1, M_1 \in \mathbb{R}^{p \times p}$ のような行列に対して、以下のようにする。

A
[V_1, V_2, \dots, V_l]
= 
[V_1, V_2, \dots, V_l]
\begin{bmatrix}
  M_1 & B_1^\mathrm{T} &                    & O \\
  B_1 & M_2            & B_2^\mathrm{T}     &   \\
      & B_2            & \dots              &   \\
      &                & \dots              & B_{j-1}^\mathrm{T} \\
  O   &                & B_{j-1}^\mathrm{T} & M_j
\end{bmatrix}

Block-Lanczos法の実装

using LinearAlgebra

orth(n, p) = (rand(n, p) |> qr |> x -> x.Q |> Matrix) # n \times p の直行行列を生成
qr2(R) = begin
    F = qr(R)
    Matrix(F.Q), F.R
end

# 
function make_ret_vals(M, V, B, n, p, q)
    T = zeros(n, n)
    VV= zeros(n, n)
    for k = 1:q
        kk = (k-1)*p+1:k*p; 
        jj = k*p+1:(k+1)*p
        T[kk, kk] = M[k]
        VV[:, kk] = V[k]
        if k == q continue end
        T[jj, kk] = B[k]
        T[kk, jj] = B[k]'
    end
    return VV, T
end

function simple_block_lanczos(
      A::AbstractMatrix
    ; n = size(A, 2)
    , q = 4
    , p = n ÷ q
)
    # 確保
    V = fill(zeros(n, p), q)
    M = fill(zeros(p, p), q)
    B = fill(zeros(p, p), q-1)

    # 初期化
    V[1] = orth(n, p)
    U    = A * V[1]
    M[1] = U' * V[1]
    R    = U - V[1] * M[1]

    # 反復
    for k = 2:q
        V[k], B[k-1] = qr2(R)
        U    = A * V[k] - V[k-1] * B[k-1]'
        M[k] = U' * V[k]
        R    = U - V[k] * M[k]
    end
    
    # 行列化
    P, T = make_ret_vals(M, V, B, n, p, q)

    return P, T
end
;

以下のコードを実行すれば、分解がしっかりされていることが分かります。

A = Symmetric(rand(20, 20)); 
V, T = simple_block_lanczos(A)

norm(A * V - V * T)

参考文献

自ページ