logo

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

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

概要

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

ビルドアップ

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

xTKx=[0xT][a11wTwK][0x]>0 \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>0K > 0 だ。便宜上 a11=1a_{11} = 1 と置くと、

A0:=[1wTwK]=[10wI][100K−wwT][1wT0I] 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}

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

A0=L1⋯LmImLmT⋯L1T A_{0} = L_{1} \cdots L_{m} I_{m} L_{m}^{T} \cdots L_{1}^{T}

のように表せる。

一般化

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

A0=[a11wTwK]=[a1100I][11a11wT1a11wK][a1100I] \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*}

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

[11a11wT1a11wK]=[101a11wI][100K−wwTa11][11a11wT0I] \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}

L1:=[a1100I][101a11wI] 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}

そして、

A1:=[100K−wwTa11] A_{1} := \begin{bmatrix} 1 & \mathbb{0} \\ \mathbb{0} & K-{{ \mathbf{w} \mathbf{w}^{T} } \over {a_{11}}} \end{bmatrix}

と置くと、

A0=[a1100I][101a11wI][100K−wwTa11][11a11w0I][a1100I]=L1A1L1T \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*}

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

A0=L1⋯LmImLmT⋯L1T A_{0} = L_{1} \cdots L_{m} I_{m} L_{m}^{T} \cdots L_{1}^{T}

が得られる。L1,L2,⋯LmL_{1} , L_{2} , \cdots L_{m} は下三角行列であり、その積L:=LmLm−1⋯L1L : = L_{m} L_{m-1} \cdots L_{1} も下三角行列であり、A0=A=LLTA_{0} = A = L L^{T} と表される。

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

L:=[l1100⋯0l21l220⋯0l31l32l33⋯0⋮⋮⋮⋱⋮ln1ln2ln3⋯lnn] 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=(aij)=LLTA = (a_{ij}) = LL^{T} を直接計算してみよう。対角成分より後は計算する必要がないので、

aij=∑s=1nlislsjT=∑s=1min⁡(i,j)lisljs a_{ij} = \sum_{s=1}^{n} l_{is} l_{sj}^{T} = \sum_{s=1}^{\min (i,j)} l_{is} l_{js}

i=j=ki=j=k ならば、 akk=∑s=1k−1lks2+lkk2lkk=akk−∑s=1k−1lks2 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=kj=k ならば、i=k+1,⋯ ,ni = k+1, \cdots , n について、 aik=∑s=1k−1lislks+liklkklik=1lkk(aik−∑s=1k−1lislks) 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)

アルゴリズム

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

ステップ1. k=1k = 1

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

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

  • ステップ2-1. lkk=akk−∑s=1k−1lks2 l_{kk} = \sqrt{ a_{kk} - \sum_{s=1}^{k-1} l_{ks}^{2} }

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

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

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

コード

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

5A1D52273

以下のコードは、上記アルゴリズムを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))

実行結果は以下の通り。 5A1D52280

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