Cholesky分解を用いた逆行列の計算方法

本稿では， Cholesky分解を用いて正定値対称行列の逆行列を計算する方法と，そのC#による実装例を示します．

#Cholesky分解

正定値対称行列 $A=[a_{ij}] \in \mathbb{R}^{n\times n}$は，適当な下三角行列 $L=[,l_{ij}] \in \mathbb{R}^{n\times n}$（上三角行列でも可）を用いて次のように分解することができる．

A = LL^{\rm T}

これを行列 $A$ のコレスキー分解という．$A^{-1}$ は次のように表される．

A^{-1} = \big(L^{-1}\big)^{\rm T}\big(L^{-1}\big)

三角行列の逆行列は容易に計算することが可能なため，行列 $A$ を下三角行列 $L$ に分解することができれば，$A^{-1}$ を求めることができる．よって，以降では行列 $A$ から下三角行列 $L$ を求めるための方法を示す．まず，視覚的に分かり易くするため，$A$, $L$ の各成分を行列形式で示す．

A = 
\begin{bmatrix}
a_{11} & a_{21} & \cdots & a_{n1} \\
a_{21} & a_{22} & \cdots & a_{n2} \\
\vdots & \vdots & \ddots & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn} \\
\end{bmatrix},\quad
L = 
\begin{bmatrix}
\, l_{11} &   0    & \cdots &    0   \\
\, l_{21} & l_{22} & \cdots & \vdots \\
\vdots & \vdots & \ddots &    0   \\
\, l_{n1} & l_{n2} & \cdots & l_{nn} \\
\end{bmatrix}

$LL^{\rm T}$ の各成分は次のようになる．

LL^{\rm T} = 
\begin{bmatrix}
\, l_{11} &   0    & \cdots &    0   \\
\, l_{21} & l_{22} & \cdots & \vdots \\
\vdots & \vdots & \ddots &    0   \\
\, l_{n1} & l_{n2} & \cdots & l_{nn} \\
\end{bmatrix}
\begin{bmatrix}
\, l_{11} & l_{21} & \cdots & l_{n1} \\
\,   0    & l_{22} & \cdots & l_{n2} \\
   \vdots & \vdots & \ddots & \vdots \\
\,   0    & \cdots &    0   & l_{nn} \\
\end{bmatrix} \\
 =
\begin{bmatrix}
\, l_{11}^2 & l_{11}l_{21}    & \cdots & l_{11}l_{n1} \\
\, l_{11}l_{21} & l_{21}^2+l_{22}^2 & \cdots & l_{21}l_{n1}+l_{22}l_{n2} \\
\vdots & \vdots & \ddots & \vdots \\
\, l_{11}l_{n1} & l_{21}l_{n1}+l_{22}l_{n2} & \cdots & l_{n1}^2+l_{n2}^2+\cdots+l_{nn}^2 \\
\end{bmatrix}

$LL^{\rm T}$ も当然，対称行列となる．$A=LL^{\rm T}$ の両辺を比較することで，次式が得られる．

\begin{array}{cccc}
a_{11} = l_{11}^2 &     &  &  \\
a_{21} = l_{11}l_{21} & a_{22} = l_{21}^2+l_{22}^2 &  &  \\
\vdots & \vdots & \ddots &  \\
a_{n1} = l_{11}l_{n1} & a_{n2} = l_{21}l_{n1}+l_{22}l_{n2} & \cdots & a_{nn} = l_{n1}^2+l_{n2}^2+\cdots+l_{nn}^2 \\
\end{array}

$L$ の対角成分を求める際には，平方根の計算が含まれるため，$\pm$の任意性がある．すなわち，$L$ の選び方は一意に定まらない．しかし，$L$ の対角成分を全て正と見なして計算を行えば，$L$ の選び方を一意に定めることができ，逆行列も問題なく求めることができる．$L$ の各成分は，上式を1列目から順に，左から計算することで求められる．

\begin{array}{cccc}
l_{11} = \sqrt{a_{11}} &     &  &  \\[1.5ex]
l_{21} = \dfrac{a_{21}}{l_{11}} &  l_{22} = \sqrt{a_{22} - l_{21}^2} &  &  \\
\vdots & \vdots & \ddots &  \\
l_{n1} = \dfrac{a_{n1}}{l_{11}} & l_{n2} = \dfrac{a_{n2} - l_{21}l_{n1}}{l_{22}} & \cdots & l_{nn} = \sqrt{a_{nn} - \big(l_{n1}^2+l_{n2}^2+\cdots+l_{n(n-1)}^2 \big)} \\
\end{array}

上式をより一般化した形式にまとめることを考える．
$L$ の非対角成分を用いて次のような行ベクトルを定める．

\boldsymbol{l}_{1} = \boldsymbol{[} \ 0 \quad 0 \quad 0 \quad \cdots \quad 0 \ \boldsymbol{]} \in \mathbb{R}^n \\
\boldsymbol{l}_{2} = \boldsymbol{[} \  l_{21} \ 0 \quad 0 \quad \cdots \quad 0 \ \boldsymbol{]} \in \mathbb{R}^n \\
\boldsymbol{l}_{3} = \boldsymbol{[} \  l_{31} \ l_{32} \ 0 \quad \cdots \quad 0 \ \boldsymbol{]} \in \mathbb{R}^n \\
  \vdots \nonumber \\
\boldsymbol{l}_{n} = \boldsymbol{[} \  l_{n1} \ l_{n2} \ \cdots \ l_{n(n-1)} \ 0 \ \boldsymbol{]} \in \mathbb{R}^n

このとき，$A$ と $L$ の各成分の関係は次のように書くことができる．

$i>j$ のとき

l_{ij} = \dfrac{a_{ij} - \boldsymbol{l}_i \boldsymbol{l}_j^{\rm T}}{l_{jj}}

$i=j$ のとき

l_{ii} = \sqrt{a_{ii} - \boldsymbol{l}_i \boldsymbol{l}_i^{\rm T}}

上記 $l_{ij}$ の計算式は，他の文献で示されているものと比べて簡潔に書けているはず．

#C#によるCholesky分解の実装例

Cholesky分解を用いて逆行列を計算する関数です．行列の積や転置を求める関数も下部にあります．

C#

        float[,] inverse(float[,] A)
        {
            float[,] Ainv = new float[A.GetLength(0), A.GetLength(1)];
            float[,] L = new float[A.GetLength(0), A.GetLength(1)];
            float[,] Linv = new float[A.GetLength(0), A.GetLength(1)];
            
            //L生成
            for (int j = 0; j < A.GetLength(1); j++)
            {

                L[j, j] = (float)Math.Sqrt(A[j, j] - calcllT(j, j, L));
                
                for (int i = j + 1; i < A.GetLength(0); i++)
                {
                    L[i, j] = (A[i, j] - calcllT(i, j, L)) / L[j, j];
                }

            }
           
            //掃き出し法でLinv計算
            for (int j = 0; j < A.GetLength(1); j++)
            {
                for (int i = j + 1; i < A.GetLength(0); i++)
                {                                        
                    L[i, j] = L[i, j] / L[i, i];                    
                }
                                
            }
            for (int i = 0; i < A.GetLength(0); i++)
            {
                Linv[i, i] = 1 / L[i, i];
                L[i, i] = 1;
            }
            for (int j = 0; j < A.GetLength(1); j++)
            {
                for (int k = 0; k < j + 1; k++)
                {
                    for (int i = j + 1; i < A.GetLength(0); i++)
                    {
                        Linv[i, k] = Linv[i, k] - L[i, j] * Linv[j, k];
                    }
                }
            }
                       
            Ainv = MatrixTimesMatrix(Transpose(Linv), Linv);

            return Ainv;
        }

        
        float calcllT(int i, int j, float[,] L)
        {
            float llT = 0;

            for (int k = 0; k < i && k < j; k++)
            {
                llT = llT + L[i, k] * L[j, k];
            }

            return llT;

        }


        float[,] MatrixTimesMatrix(float[,] A, float[,] B)
        {
            float[,] product = new float[A.GetLength(0), B.GetLength(1)];

            for (int i = 0; i < A.GetLength(0); i++)
            {
                for (int j = 0; j < B.GetLength(1); j++) 
                {
                    for (int k = 0; k < A.GetLength(1); k++) 
                    {
                        product[i, j] += A[i, k] * B[k, j];
                    }
                }
            }

            return product;
        }

        float[,] Transpose(float[,] A)
        {
            float[,] AT = new float[A.GetLength(1), A.GetLength(0)];

            for (int i = 0; i < A.GetLength(1); i++)
            {
                for (int j = 0; j < A.GetLength(0); j++)
                {
                    AT[i, j] = A[j, i];
                }
            }

            return AT;
        }

引数の行列が正定値対称行列でない場合，ゼロ除算が起こってしまうので注意．