📂行列代数

正定値行列のコレスキー分解

概要

可逆行列から対角行列への条件が強化され、LDU分解ができたように、正定値行列への条件強化で、より効率的な行列分解であるコレスキー分解^{cholesky decomposition}ができる。

ビルドアップ

$m \times m$ 正定値行列 $A : = \begin{bmatrix} a_{11} & \mathbf{w}^{T} \\ \mathbf{w} & K \end{bmatrix} > 0$ を考えると、$a_{11}$ は正の数であり、$\mathbf{w} \in \mathbb{R}^{m-1}$ 且つ $K \in \mathbb{R}^{(m-1) \times (m-1)}$ だ。もし $a_{11} \le 0$ ならば、$\mathbf{e}_{1} = (1, 0, \cdots , 0)$ に対して、$\mathbf{e}_{1}^{T} A \mathbf{e}_{1} = a_{11} \le 0$ なので、$A$ は正定値になれない。一方で、任意のベクトル$\mathbb{0 } \ne \mathbf{x} \in \mathbb{R}^{m-1}$ について

$$ \mathbf{x}^{T} K \mathbf{x} = \begin{bmatrix} 0 & \mathbf{x}^{T} \end{bmatrix} \begin{bmatrix} a_{11} & \mathbf{w}^{T} \\ \mathbf{w} & K \end{bmatrix} \begin{bmatrix} 0 \\ \mathbf{x} \end{bmatrix} > 0 $$

となるので、$K > 0$ だ。便宜上 $a_{11} = 1$ と置くと、

$$ A_{0} : = \begin{bmatrix} 1 & \mathbf{w}^{T} \\ \mathbf{w} & K \end{bmatrix} = \begin{bmatrix} 1 & \mathbb{0} \\ \mathbf{w} & I \end{bmatrix} \begin{bmatrix} 1 & \mathbb{0} \\ \mathbb{0} & K - \mathbf{w} \mathbf{w}^{T} \end{bmatrix} \begin{bmatrix} 1 & \mathbf{w}^{T} \\ \mathbb{0} & I \end{bmatrix} $$

と表せる。ここで、$L_{1} := \begin{bmatrix} 1 & \mathbf{w}^{T} \\ \mathbf{w} & K \end{bmatrix}$ は上三角行列、$A_{1} := \begin{bmatrix} 1 & \mathbb{0} \\ \mathbb{0} & K - \mathbf{w} \mathbf{w}^{T} \end{bmatrix}$ は正定値なので、$K - \mathbf{w} \mathbf{w}^T > 0$ だ。したがって、このプロセス $A_{n-1} = L_{n} A_{n} L_{n}^{T}$ を $n=m$ まで繰り返すと、

$$ A_{0} = L_{1} \cdots L_{m} I_{m} L_{m}^{T} \cdots L_{1}^{T} $$

のように表せる。

一般化

今、$a_{11}>0$ の場合について一般化してみよう。$A := A_{0}$ とすると、

$$ \begin{align*} A_{0} &= \begin{bmatrix} a_{11} & \mathbf{w}^{T} \\ \mathbf{w} & K \end{bmatrix} \\ =& \begin{bmatrix} \sqrt{a_{11}} & \mathbb{0} \\ \mathbb{0} & I \end{bmatrix} \begin{bmatrix} 1 & {{1} \over {\sqrt{a_{11}}}} \mathbf{w}^{T} \\ {{1} \over {\sqrt{a_{11}}}} \mathbf{w} & K \end{bmatrix} \begin{bmatrix} \sqrt{a_{11}} & \mathbb{0} \\ \mathbb{0} & I \end{bmatrix} \end{align*} $$

$\begin{bmatrix} 1 & {{1} \over {\sqrt{a_{11}}}} \mathbf{w}^{T} \\ {{1} \over {\sqrt{a_{11}}}} \mathbf{w} & K \end{bmatrix}$は$a_{11}=1$ の場合で扱った形なので、

$$ \begin{bmatrix} 1 & {{1} \over {\sqrt{a_{11}}}} \mathbf{w}^{T} \\ {{1} \over {\sqrt{a_{11}}}} \mathbf{w} & K \end{bmatrix} = \begin{bmatrix} 1 & \mathbb{0} \\ {{1} \over {\sqrt{a_{11}}}} \mathbf{w} & I \end{bmatrix} \begin{bmatrix} 1 & \mathbb{0} \\ \mathbb{0} & K-{{ \mathbf{w} \mathbf{w}^{T} } \over {a_{11}}} \end{bmatrix} \begin{bmatrix} 1 & {{1} \over {\sqrt{a_{11}}}} \mathbf{w}^{T} \\ \mathbb{0} & I \end{bmatrix} $$

$$ L_{1} := \begin{bmatrix} \sqrt{a_{11}} & \mathbb{0} \\ \mathbb{0} & I \end{bmatrix} \begin{bmatrix} 1 & \mathbb{0} \\ {{1} \over {\sqrt{a_{11}}}} \mathbf{w} & I \end{bmatrix} $$

そして、

$$ A_{1} := \begin{bmatrix} 1 & \mathbb{0} \\ \mathbb{0} & K-{{ \mathbf{w} \mathbf{w}^{T} } \over {a_{11}}} \end{bmatrix} $$

と置くと、

$$ \begin{align*} A_{0} =& \begin{bmatrix} \sqrt{a_{11}} & \mathbb{0} \\ \mathbb{0} & I \end{bmatrix} \begin{bmatrix} 1 & \mathbb{0} \\ {{1} \over {\sqrt{a_{11}}}} \mathbf{w} & I \end{bmatrix} \begin{bmatrix} 1 & \mathbb{0} \\ \mathbb{0} & K-{{ \mathbf{w} \mathbf{w}^{T} } \over {a_{11}}} \end{bmatrix} \begin{bmatrix} 1 & {{1} \over {\sqrt{a_{11}}}} \mathbf{w} \\ \mathbb{0} & I \end{bmatrix} \begin{bmatrix} \sqrt{a_{11}} & \mathbb{0} \\ \mathbb{0} & I \end{bmatrix} \\ =& L_{1} A_{1} L_{1}^{T} \end{align*} $$

同様にプロセス $A_{n-1} = L_{n} A_{n} L_{n}^{T}$ を $n=m$ まで繰り返すと、

$$ A_{0} = L_{1} \cdots L_{m} I_{m} L_{m}^{T} \cdots L_{1}^{T} $$

が得られる。$L_{1} , L_{2} , \cdots L_{m}$ は下三角行列であり、その積$L : = L_{m} L_{m-1} \cdots L_{1}$ も下三角行列であり、$A_{0} = A = L L^{T}$ と表される。

こうして、下三角行列 $L$ だけで $A = L L^{T}$ を表すことが、コレスキー分解と言われる。正定値という条件はかなり強力だが、その分解は驚くほど単純で、非常に少ない計算で求めることができる。これから具体的に、

$$ L : = \begin{bmatrix} l_{11} & 0 & 0 & \cdots & 0 \\ l_{21} & l_{22} & 0 & \cdots & 0 \\ l_{31} & l_{32} & l_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \cdots & l_{nn} \end{bmatrix} $$

を見つけるために、$A = (a_{ij}) = LL^{T}$ を直接計算してみよう。対角成分より後は計算する必要がないので、

$$ a_{ij} = \sum_{s=1}^{n} l_{is} l_{sj}^{T} = \sum_{s=1}^{\min (i,j)} l_{is} l_{js} $$

$i=j=k$ ならば、 $$ a_{kk} = \sum_{s=1}^{k-1} l_{ks}^{2} + l_{kk}^2 \\ l_{kk} = \sqrt{ a_{kk} - \sum_{s=1}^{k-1} l_{ks}^{2} } $$

$j=k$ ならば、$i = k+1, \cdots , n$ について、 $$ a_{ik} = \sum_{s=1}^{k-1} l_{is} l_{ks} + l_{ik} l_{kk} \\ l_{ik} = {{1} \over {l_{kk}}} \left( a_{ik} - \sum_{s=1}^{k-1} l_{is} l_{ks} \right) $$

アルゴリズム

$(a_{ij}) \in \mathbb{R}^{n \times n}$ が正定値行列だとする。

ステップ1. $k = 1$

$d_{11} = \sqrt{a_{11}}$ を代入して、$\displaystyle l_{i1} = {{1} \over {d_{11}}} a_{i1}$ を計算する。

---

ステップ2. $k = 2, 3, \cdots , n-1$ に対して、以下を計算する。

• ステップ2-1. $$ l_{kk} = \sqrt{ a_{kk} - \sum_{s=1}^{k-1} l_{ks}^{2} } $$

• ステップ2-2. $i = k+1 , k+2 , \cdots , n-1$ に対して、以下を計算する。 $$ l_{ik} = {{1} \over {l_{kk}}} \left( a_{ik} - \sum_{s=1}^{k-1} l_{is} l_{ks} \right) $$

---

ステップ3. $k = n$ に対して、以下を計算する。 $$ l_{nn} = \sqrt{ a_{nn} - \sum_{s=1}^{n-1} l_{ns}^{2} } $$

---

• LU分解やLDU分解と比べて、計算量が目に見えて減ったのを確認できる。さらに、行なっている計算も複雑ではなく、内積やノルムを求めるのに似ていることに注目しよう。

コード

さあ、実際にコレスキー分解を試してみよう。

以下のコードは、上記アルゴリズムをRで実装したものである。

B <- matrix(c(2,1,0,1,2,1,0,1,2), ncol= 3); B
Cholesky<-function(A){ 
    n=length(A[,1])
    L<-diag(0,n)
    k=1
    L[1,1] = sqrt(A[1,1])
    L[,1] = A[,1]/L[1,1]
    for(k in 2:(n-1))  {
        L[k,k] = sqrt(A[k,k] - sum(L[k,1:(k-1)])^2 )
        for(i in (k+1):n)    {
            L[i,k] = ( A[i,k] - sum( L[i,1:(k-1)] * L[k,1:(k-1)] ) )/L[k,k]
        }
    }
    L[n,n] = sqrt( A[n,n] - sum(L[n,1:(n-1)]^2) )
    return(list(L=L))
}

Cholesky(B)
t(chol(B))

実行結果は以下の通り。

Rに内蔵されているコレスキー分解の結果と比較して、正確に一致していることが確認できる。