ブロック三重対角行列

$A_i, B_i, C_i , (i=1,\cdots, n)$ を $m\times m$ 行列、$z$ を実数の定数とするとき、

\begin{align}
M(z) \equiv
\begin{pmatrix}
A_1 & B_1 & & \frac{1}{z}C_1 \\
C_2 & \ddots & \ddots &  \\
& \ddots &  & B_{n-1}  \\
zB_{n} &  & C_{n} & A_{n} 
\end{pmatrix}
\end{align}

の形をした行列を、ブロック三重対角行列と呼びます。このタイプの行列は様々な場面で登場するため、効率的に逆行列や行列式を計算するためのアルゴリズムが色々考えられています。この記事では、ブロック三重対角行列の行列式について成立するMolinariの公式を紹介します。Molinariの公式は、数値計算のアルゴリズムとして見たときには、必ずしも効率が良いわけではありませんが、導出が直感的で解析計算に応用できる場面が色々ありそうです。

Molinariの公式

$n\geq 3$ であるとし、$B_i, C_i$ は正則とします。このときブロック三重対角行列の行列式について、以下のような関係式が成り立つことがMolinariさんによって示されています。

\begin{align}
\det M(z) 
&= \frac{(-1)^{nm}}{(-z)^m} \det\left[ T-zI \right] \det\left[ B_1 \cdots B_{n} \right] \\
&= (-1)^{nm} (-z)^m \det\left[ T^{-1 }-\frac{1}{z} \right] \det\left[ C_1 \cdots C_{n} \right]
\end{align}

ここで $T$ は

\begin{align}
T \equiv 
\begin{pmatrix}
-B_{n}^{-1}A_{n} & -B_{n}^{-1}C_{n} \\ 1 & 0
\end{pmatrix}
\cdots
\begin{pmatrix}
-B_{1}^{-1}A_{1} & -B_{1}^{-1}C_{1} \\ 1 & 0
\end{pmatrix}
\end{align}

で定義されるサイズ $2m$ の行列で、転送行列と呼ばれます。$M(z)$ のサイズは $nm$ だったので、$\det M(z)$ を直接計算するよりも、上の式の右辺を計算する方が、はるかに計算コストが少なくて済みます。

公式の証明

$\Psi = (\psi_1, \psi_2, \dots, \psi_n)^t$ とおき、方程式 $M(z)\Psi = 0$ を考えます。$\psi_i$ はそれぞれ $m$ 成分の縦ベクトルであるとします。この方程式を成分毎に具体的に書き下すと

\begin{align}
\frac{1}{z}C_1 \psi_{n} + A_1 \psi_1 + B_1\psi_2 &= 0 \\
C_2 \psi_1 + A_2\psi_2 + B_2\psi_3 &= 0 \\
&\vdots \\
C_k \psi_{k-1} + A_k\psi_{k} + B_k\psi_{k+1} &= 0 \\
&\vdots \\
%C_{n-1} \psi_{n-2} + A_{n-1}\psi_{n-1} + B_{n-1}\psi_{n-2} &= 0 \\
C_{n}\psi_{n-1} + A_{n}\psi_{n} + zB_{n} \psi_{1} &= 0
\end{align}

のようになります。これを2成分ずつまとめて

\begin{align}
\begin{pmatrix} \psi_2 \\ \psi_{1} \end{pmatrix}
&=
\begin{pmatrix}
-B_1^{-1}A_1 & -B_1^{-1} C_1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix} \psi_1 \\ \frac{1}{z}\psi_{n} \end{pmatrix} \\
%
\begin{pmatrix} \psi_3 \\ \psi_{2} \end{pmatrix}
&=
\begin{pmatrix}
-B_2^{-1}A_2 & -B_2^{-1} C_2 \\
1 & 0
\end{pmatrix}
\begin{pmatrix} \psi_2 \\ \psi_{1} \end{pmatrix} \\
%
&\vdots \\
%
\begin{pmatrix} z\psi_{1} \\ \psi_{n} \end{pmatrix}
&=
\begin{pmatrix}
-B_{n}^{-1}A_{n} & -B_{n}^{-1} C_{n} \\
1 & 0
\end{pmatrix}
\begin{pmatrix} \psi_{n} \\ \psi_{n-1} \end{pmatrix}
\end{align}

のように書き換えてみます。これらの式を見ると、$(\psi_{k} ,, \psi_{k-1})^t$ が一つ前のベクトル $(\psi_{k-1} ,, \psi_{k-2})^t$ によって定まるという関係が延々連なっていることがわかります。よってこれらをすべてつなぎ合わせることで、

\begin{align}
z\begin{pmatrix} \psi_{1} \\ \frac{1}{z}\psi_{n} \end{pmatrix}
=
\begin{pmatrix} z\psi_{1} \\ \psi_{n} \end{pmatrix}
=
T
\begin{pmatrix} \psi_1 \\ \frac{1}{z}\psi_{n} \end{pmatrix}
\end{align}

すなわち、

\begin{align}
(T -zI) \begin{pmatrix} \psi_1 \\ \frac{1}{z}\psi_{n} \end{pmatrix} = 0
\end{align}

という方程式が得られます。これが非自明な解をもつのは $\det(T -zI) = 0$ の場合に限られます。一方、もともとの方程式が非自明な解をもつのは $\det M(z) = 0$ の場合に限られます。$\det(T -zI)$ は $z$ の $2m$ 次多項式、$\det M(z)$ は $m$ 次多項式だから、

z^m \det M(z) \propto \det(T -zI)

でなければなりません。これでほぼほぼ証明が終わりました。

最後に比例係数を求めます。そのためには、特殊な $z$ に限って考えれば十分なので、$z \gg 1$ の場合を考えてみます。すると、$M(z)$ の右上の隅のブロックが無視できて、

\begin{align}
\det M(z) 
\sim
\det \begin{pmatrix}
A_1 & B_1 & & 0 \\
C_2 & \ddots & \ddots &  \\
& \ddots &  & B_{n-1}  \\
zB_{n} &  & C_{n} & A_{n} 
\end{pmatrix}
=
(-1)^{nm} (-z)^m \det\left( B_0 \cdots B_{n-1}\right)
\end{align}

と振る舞うことがわかります。従って、比例係数は $(-1)^{nm} (-1)^m \det\left( B_0 \cdots B_{n-1}\right)$ と決まります。公式の2段目も同様に示すことができます。

Julia で計算してみる

数値計算でチェックしてみます。
以下はランダムな要素をもつブロック三重対角行列を作る関数です。

function create_block_tridiagonal_matrix(n, m)
    M = zeros(Float64, n*m, n*m)
    
    # A1, A2, ... An
    @inbounds for j in 1:n
        @views M[(j-1)*m+1:j*m,(j-1)*m+1:j*m] = rand(m, m)
    end

    # B1, B2, ... Bn
    @inbounds for j in 2:n
        @views M[(j-2)*m+1:(j-1)*m, (j-1)*m+1:j*m] = rand(m, m)
    end
    @views M[(n-1)*m+1:n*m, 1:m] = rand(m, m)

    # C1, C2, ... Cn
    @views M[1:m,(n-1)*m+1:n*m] = rand(m, m)
    @inbounds for j in 2:n
        @views M[(j-1)*m+1:j*m,(j-2)*m+1:(j-1)*m] = rand(m, m)
    end
    return M
end

以下は部分行列$A_n, B_n, C_n$を取り出す関数です。上とやっていることが重複しているのですが、あまり深く考えないことにします。

function get_An(n, m, M)
    @views M[(n-1)*m+1:n*m,(n-1)*m+1:n*m]
end

function get_Bn(n, m, M)
    if n*m == size(M)[1]
        @views M[(n-1)*m+1:n*m, 1:m]
    else
        @views M[(n-1)*m+1:n*m, n*m+1:(n+1)*m]
    end
end

function get_Cn(n, m, M)
    if n == 1
        nm = size(M)[1]
        @views M[1:m,nm-m+1:nm]
    else
        @views M[(n-1)*m+1:n*m,(n-2)*m+1:(n-1)*m]
    end
end

これらを用いると、転送行列は以下のようにして計算することができます。

function transfer_matrix(n, m, M)
    T = Matrix{Float64}(I, 2m, 2m)
    for j in 1:n
        Tj = zeros(2m, 2m)
        @views Tj[1:m, 1:m] = - inv(get_Bn(j, m, M)) * get_An(j, m, M)
        @views Tj[1:m, m+1:2m] = - inv(get_Bn(j, m, M)) * get_Cn(j, m, M)
        @views Tj[m+1:2m, 1:m] = Matrix{Float64}(I, m, m)
        T = Tj * T
    end
    return T
end

Molinariの公式には、$\det B_1 \cdots B_n$ を計算する箇所があるので、それを行うコードも用意します。

function detB1B2Bn(n, m, M)
    ΠBj = Matrix{Float64}(I, m, m)
    for j in 1:n
        ΠBj *= get_Bn(j, m, M)
    end
    return det(ΠBj)
end

以上で、Molinariの公式に従って行列式を計算する準備ができました。

function MolinariFormula(n, m, M)
    T = transfer_matrix(n, m, M)
    (-1)^(n*m) / (-1)^m * det(T-I) * detB1B2Bn(n, m, M)
end

実際に計算してみると

M = create_block_tridiagonal_matrix(4, 2)
det(M) ≈ MorinariFormula(4, 2, M)

などが true であることが確認できます。計算速度はどんなものかというと

M = create_block_tridiagonal_matrix(8, 4)
@time det(M)                   # 0.000463 seconds (4 allocations: 8.500 KiB)
@time MolinariFormula(8, 4, M) # 0.000095 seconds (194 allocations: 64.562 KiB)

M = create_block_tridiagonal_matrix(16, 8)
@time det(M)                    # 0.004803 seconds (5 allocations: 129.266 KiB)
@time MolinariFormula(16, 8, M) # 0.000197 seconds (378 allocations: 297.516 KiB)

となりました。(Julia 1.4.0 で実行。)行列サイズが大きければ大きいほど、Molinariの公式を用いたほうが効率が良くなります。Molinariの公式の計算量については参考文献を御覧ください。サンプルコードはGistで公開しています。

おまけ：フェルミオン場の理論

Molinariの公式に登場した転送行列は物理でもおなじみの概念ですが、この公式自身フェルミオン場の理論でよく知られているある性質の離散バージョンであることがわかります。

フェルミオン場の理論として、作用が

\begin{align}
S = S_\text{B} + \int d\tau \int d^dx \psi^*(x,\tau) O \psi(x,\tau)
\end{align}

で与えられる場合を考えます。$S_\text{B}$はボソン場の作用、$\psi$はフェルミオン場、$d$は空間次元、$O$は微分、質量、化学ポテンシャル、ボソン場との結合項などを含むもので、一般に以下のように書けます

\begin{align}
O = \frac{\partial}{\partial \tau} + H
\end{align}

$H$の具体例は

\begin{align}
H = -\nabla^2 + \omega + \lambda \phi(x,\tau)
\end{align}

などです。$A$を物理量とすると、その期待値は経路積分を用いて

\begin{align}
\langle A \rangle = \frac{1}{Z} \int \mathcal{D}\phi e^{-S_\text{B}}\det\left( \frac{\partial}{\partial \tau} + H \right) A
\end{align}

と表されます。ただし$Z$は大分配関数で

\begin{align}
Z = \int \mathcal{D}\phi e^{-S_\text{B}}\det\left( \frac{\partial}{\partial \tau} + H \right)
\end{align}

で与えられます。この式中の行列式はフェルミオン場を積分した結果として出てくるもので、フェルミオン行列式と呼ばれます。このフェルミオン行列式について、一般に次が成り立ちます。

\begin{align}
\det_{x,\tau} \left[ \frac{\partial}{\partial \tau} + H \right]
\propto
\det_{x} \left[ I + \mathrm{T} \exp\left( -\int_{0}^{\beta} d\tau H(\tau) \right) \right]
\end{align}

ここで$\mathrm{T}$は虚時間順序積を表します。左辺の $\det_{x,\tau}$ は中身の関数の引数が時空全体であること、左辺の $\det_{x}$ は空間だけであることを表します。これはちょうど、Molinariの公式の連続理論バージョンになっていると言えます。

参考文献

Molinariの原論文: L. G. Molinari, Linear Algebra and its Applications 429, 2221 (2008), arXiv:0712.0681.
Molinariの方法と他の計算量を比較した論文: Filipiak, Markiewicz and Sawikowska, Electronic Journal of Linear Algebra 25, 102-118 (2012).
場の理論におけるフェルミオン行列式の計算について: R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).