不完全コレスキー分解

数値解析において、対称正定値行列の不完全コレスキー分解は、コレスキー分解の疎近似であり、共役勾配法などのアルゴリズムの前処理としてよく使用されます。

正定値行列Aのコレスキー分解はA = LL *で、 L は下三角行列です。不完全なコレスキー分解は、ある意味でLに近い疎な下三角行列Kによって与えられます。対応する前処理はKK *です。

このような行列Kを見つける一般的な方法の 1 つは、 KがAと同じスパース パターンを持つ正確なコレスキー分解を見つけるアルゴリズムを使用することです( Aの対応するエントリも 0 の場合、 Kの任意のエントリは 0 に設定されます)。これにより、行列Aと同じくらいスパースな不完全なコレスキー分解が得られます。

例として次の行列を考えてみましょう。

${\displaystyle \mathbf {A} ={\begin{bmatrix}5&-2&0&-2&-2\\-2&5&-2&0&0\\0&-2&5&-2&0\\-2&0&-2&5&-2\\-2&0&0&-2&5\\\end{bmatrix}}}$

完全な正則コレスキー分解を適用すると、次のようになります。

${\displaystyle \mathbf {L} ={\begin{bmatrix}2.24&0&0&0&0\\-0.89&2.05&0&0&0\\0&-0.98&2.02&0&0\\-0.89&-0.39&-1.18&1.63&0\\-0.89&-0.39&-0.19&-1.95&0.45\\\end{bmatrix}}}$

そして、定義によれば:

${\displaystyle \mathbf {A} =\mathbf {L} \mathbf {L'} }$

しかし、コレスキー分解を適用すると、元の行列のゼロ要素の一部が、分解された行列では非ゼロ要素になることがわかります。この例では、要素 (4,2)、(5,2)、(5,3) がこれに該当します。これらの要素は「フィルイン」と呼ばれます。

これは、それ自体は問題ではありませんが、スパース行列を扱うときには非常に問題になります。これは、フィルイン生成がほとんど予測不可能であり、行列のスパース性が低下し、スパース行列アルゴリズムの効率に影響するためです。

したがって、行列計算におけるコレスキー分解の重要性を考えると、通常の方法を再利用してフィルイン生成を排除することが非常に重要です。不完全コレスキー分解は、通常の方法で生成されたものと同様の行列 L を生成しながら、元の行列のゼロ要素を保存することで、まさにそれを実現します。

当然：

${\displaystyle \mathbf {A} \neq \mathbf {L_{ichol}} \mathbf {L_{ichol}'} }$

不完全コレスキー分解によって生成された行列 L にその転置行列を乗算すると、元の行列は生成されませんが、類似の行列が生成されます。

からまでの場合: ${\displaystyle i}$${\displaystyle 1}$${\displaystyle N}$

${\displaystyle L_{ii}=\left({a_{ii}-\sum \limits _{k=1}^{i-1}{L_{ik}^{2}}}\right)^{1 \over 2}}$

からまでの場合: ${\displaystyle j}$${\displaystyle i+1}$${\displaystyle N}$

${\displaystyle L_{ji}={1 \over {L_{ii}}}\left({a_{ji}-\sum \limits _{k=1}^{i-1}{L_{ik}L_{jk}}}\right)}$

GNU Octave言語での不完全コレスキー分解の実装。分解は下三角行列として保存され、上三角の要素はゼロに設定されます。

関数 a = ichol ( a ) n = size ( a , 1 );  
	  

	 k  =  1の場合: n 
		a ( k , k )  =  sqrt ( a ( k , k )); 
		i = ( k + 1 )の場合 : n if ( a ( i , k ) != 0 ) a ( i , k ) = a ( i , k ) / a ( k , k ); endif endfor j = ( k + 1 )の場合: n i = j の場合: n if ( a ( i , j ) != 0 ) a ( i , j ) = a ( i , j ) - a ( i , k ) * a ( j , k ); endif endfor endfor endfor  
		       
		                      
		    
		
		   
		       
		           
		                  
		        
		    
		
	

     i  =  1 の場合: n 
        j = i + 1 の場合 : n a ( i , j ) = 0 ; endfor endfor endfunction  
              
        
                

この記事の冒頭で示した行列をもう一度考えてみましょう。この行列は対称であり、この方法では下三角要素のみを使用するため、次のように表すことができます。

${\displaystyle \mathbf {A_{tri}} ={\begin{bmatrix}5&0&0&0&0\\-2&5&0&0&0\\0&-2&5&0&0\\-2&0&-2&5&0\\-2&0&0&-2&5\end{bmatrix}}}$

より具体的には、座標リストとしてのスパース形式で、最初に行をスイープします。

値 5 -2 -2 -2 5 -2 5 -2 5 -2 5 -2 5
行 1 2 4 5 2 3 3 4 4 5 5
列 1 1 1 1 2 2 3 3 4 4 5

次に、(1,1) の平方根をとり、その結果で他の (i,1) 要素を割ります。

値 2.24 -0.89 -0.89 -0.89 | 5 -2 5 -2 5 -2 5
行 1 2 4 5 | 2 3 3 4 4 5 5
列 1 1 1 1 | 2 2 3 3 4 4 5

その後、列が 1 より大きい他のすべての要素について、(i,1) と (j,1) が存在する場合は (i,j)=(i,j)-(i,1)*(j,1) を計算します。たとえば、(5,4) = (5,4)-(5,1)*(4,1) = -2 -(-0.89*-0.89) = -2.8 です。

値 2.24 -0.89 -0.89 -0.89 | 4.2 -2 5 -2 4.2 -2.8 4.2
行 1 2 4 5 | 2 3 3 4 4 5 5
列 1 1 1 1 | 2 2 3 3 4 4 5
                                  ↑ ↑ ↑ ↑

要素 (2,2)、(4,4)、(5,4)、(5,5) (矢印でマーク) は、このルールに従っているため、再計算されています。一方、要素 (3,2)、(3,3)、(4,3) は、要素 (2,1) と (4,1) が存在しても、要素 (3,1) が存在しないため、再計算されません。

ここで、(i,2) について同じ手順を繰り返します。(2,2) の平方根を取り、その結果で他の (i,2) 要素を割ります。

値 2.24 -0.89 -0.89 -0.89 | 2.05 -0.98 | 5 -2 4.2 -2.8 4.2
行 1 2 4 5 | 2 3 | 3 4 4 5 5
列 1 1 1 1 | 2 2 | 3 3 4 4 5

再び、列が2より大きい要素については、(i,2)と(j,2)が存在する場合、(i,j)=(i,j)-(i,2)*(j,2)を計算します。

値 2.24 -0.89 -0.89 -0.89 | 2.05 -0.98 | 4.05 -2 4.2 -2.8 4.2
行 1 2 4 5 | 2 3 | 3 4 4 5 5
列 1 1 1 1 | 2 2 | 3 3 4 4 5
                                               ↑

(i,3)についても繰り返します。(3,3)の平方根を取り、もう一方の(i,3)を割ります。

値 2.24 -0.89 -0.89 -0.89 2.05 -0.98 | 2.01 -0.99 | 4.2 -2.8 4.2
行 1 2 4 5 2 3 | 3 4 | 4 5 5
列 1 1 1 1 2 2 | 3 3 | 4 4 5

列が3より大きい要素については、(i,3)と(j,3)が存在する場合、(i,j)=(i,j)-(i,3)*(j,3)を計算します。

値 2.24 -0.89 -0.89 -0.89 2.05 -0.98 | 2.01 -0.99 | 3.21 -2.8 4.2
行 1 2 4 5 2 3 | 3 4 | 4 5 5
列 1 1 1 1 2 2 | 3 3 | 4 4 5
                                                          ↑

(i,4)についても繰り返します。(4,4)の平方根を取り、もう一方の(i,4)を割ります。

値 2.24 -0.89 -0.89 -0.89 2.05 -0.98 2.01 -0.99 | 1.79 -1.56 | 4.2
行 1 2 4 5 2 3 3 4 | 4 5 | 5
列 1 1 1 1 2 2 3 3 | 4 4 | 5

列が4より大きい要素については、(i,4)と(j,4)が存在する場合、(i,j)=(i,j)-(i,4)*(j,4)を計算します。

値 2.24 -0.89 -0.89 -0.89 2.05 -0.98 2.01 -0.99 | 1.79 -1.56 | 1.76
行 1 2 4 5 2 3 3 4 | 4 5 | 5
列 1 1 1 1 2 2 3 3 | 4 4 | 5
                                                                     ↑

最後に(5,5)の平方根をとります。

値 2.24 -0.89 -0.89 -0.89 2.05 -0.98 2.01 -0.99 1.79 -1.56 | 1.33
行 1 2 4 5 2 3 3 4 4 5 | 5
列 1 1 1 1 2 2 3 3 4 4 | 5

マトリックスを完全な形に拡張します。

${\displaystyle \mathbf {L} ={\begin{bmatrix}2.24&0&0&0&0\\-0.89&2.05&0&0&0\\0&-0.98&2.01&0&0\\-0.89&0&-0.99&1.79&0\\-0.89&0&0&-1.56&1.33\end{bmatrix}}}$

この場合、元の行列と比較して、フィルインは生成されていないことに注意してください。要素 (4,2)、(5,2)、(5,3) は依然としてゼロです。

しかし、L をその転置に乗算すると次のようになります。

${\displaystyle \mathbf {LL'} ={\begin{bmatrix}5&-2&0&-2&-2\\-2&5&-2&0.8&0.8\\0&-2&5&-2&0\\-2&0.8&-2&5&-2\\-2&0.8&0&-2&5\end{bmatrix}}}$

分解では、埋め込んだ要素を排除するためにすべての要素が考慮されなかったため、元の行列とは少し異なる行列が得られます。

MATLAB で実装された不完全コレスキー分解のスパース バージョン (上記と同じ手順) を以下に示します。当然、MATLAB にはスパース行列を処理する独自の手段がありますが、以下のコードは教育目的で明示的に作成されています。このアルゴリズムは、行列を連続した 1D 配列として扱い、ゼロ要素を自動的に回避するため、効率的です。

function  A = Sp_ichol ( A ) 
n = size ( A , 1 ); ncols = A ( n ) .col ; c_end = 0 ; for col = 1 : ncols is_next_col = 0 ; c_start = c_end + 1 ; for i = c_start : n if A ( i ) .col == col % 現在の列 (col) 内: if A ( i ) .col == A ( i ) .row A ( i ) .val = sqrt ( A ( i ) .val ); % 現在の列の対角要素の平方根を取るdiv = A ( i ) .val ; else A ( i ) .val = A ( i ) .val / div ; % 現在の列の他の要素を対角要素の平方根で割るend end if A ( i ). col > col % 次の列 (col+1 ... ncols) で: if is_next_col == 0 c_end = i - 1 ; is_next_col = 1 ; end v1 = 0 ; v2 = 0 ; for j = c_start : c_end if A ( j ). col == col if A ( j ). row == A ( i ). row % 行インデックスが現在の要素 A(i) の行インデックスと同じである現在の列 (col) 要素 A(j) を検索します。v1 = A ( j ). val ;	
	
    
     
        
        
         
              
                  
                     
                    
                
                     
                
            
              
                 
                    
                    
                
                
                
                 
                     
                          
                            
                        end 
if A ( j ). row == A ( i ). col % 行インデックスが現在の要素 A(i) の列インデックスと同じである現在の列 (col) の要素 A(j) を検索しますv2 = A ( j ). val ; end if v1 ~= 0 && v2 ~= 0 % これらの要素が現在の列 (col) に存在する場合は、現在の要素 A(i) を再計算します: A ( i ). val = A ( i ). val - v1 * v2 ; break ; end end end end end end end                          
                            
                        
                            
                            
                            
                        
                    
                
            
        
    

• CFD Online wiki の不完全コレスキー分解
• ゴルブ、ジーン H. ;ヴァン・ローン、チャールズ F. (1996)、マトリックス計算(第 3 版)、ジョンズ・ホプキンス、ISBN 978-0-8018-5414-9セクション10.3.2を参照してください。