行列指数関数とベクトルとの積を高速に計算する方法について。

導入

問題設定

以下の行列指数関数とベクトルとの積を計算する。

exp ( 𝐀t ) 𝐛

記号の定義

t の存在理由

多くの応用で、複数の t の値について計算する。

そのため、あえて 𝐀 と t に分け、一つの記号にまとめていない。

nnz の意味

nnz は "the Number of None Zero" の略で、行列の非ゼロ要素の個数を表す。

そのため、nnz(𝐀)⁄n2 は行列 𝐀 が n 次の正方行列なら、𝐀 の疎性を表す。

用語の定義

行列指数関数

行列指数関数は以下の Taylor 展開で定義される。

exp ( 𝐌 ) = ∑ k=0 ∞ 1 k! 𝐌 k

これは以下の通常指数関数の Taylor 展開に由来する。

exp ( x ) = ∑ k=0 ∞ 1 k! x k

両者は引数の型をのぞいて全く同じ構造である。

近似精度

結果は Taylor 展開の m-1 次まで求める。

exp ( 𝐀t ) 𝐛 ≈ ∑ k=0 m-1 1 k! t k 𝐀 k 𝐛

ただし、m≪n とする。

各種アルゴリズム

二つの手法があり、それぞれ長所と短所がある。

推奨手法

少数の t について計算するなら手法 A が、そうでなければ手法 B が推奨される。なぜなら、手法 A は手法 B にある重い前処理 (後述の Arnordi 法) がない反面、t を指定して計算する本処理が重くなっている。

時間計算量

それぞれの手法と処理部分の時間計算量は以下の通り。

本処理 前処理
手法 A O ( m · nnz(𝐀) ) O ( 0 )
手法 B O ( m·n ) O ( m · nnz(𝐀) + m2n )

ただし、行列 𝐀 に密行列のためのデータ構造を利用している場合、nnz(𝐀) の箇所は n2 に増える。しかし、これを嫌って疎行列のためのデータ構造を利用した場合、メモリ効率が悪くなり計算量に表れない比例係数が悪化しやすいため、推奨されない。

手法 A

計算順序の単純な工夫

計算には Taylor 展開が含まれるため、有限級数近似を行う。

ここで、計算順序を少し工夫するだけで、計算量を改善できる。

そのまま計算した場合

行列指数関数の Taylor 展開の計算後に、ベクトルとの積を行う。

計算量は密行列でも疎行列でも O(m·n3) になる。なお、n3 の箇所は n 次正方行列どうしの行列積の計算量に由来する。つまり、Taylor 展開の各回の計算では、保存した前回の結果との行列積が主な計算となる。なお、疎行列でも密行列でも計算量が同じなのは、行列積のたびに行列が密化していくためである。

ちなみに、行列指数関数の計算には Pade 近似など他の方法もある (詳細は割愛)。これらはより少ない級数で近似できるなど効率的ではあるが、影響は高々定数倍である。

工夫して計算した場合

行列指数関数の Taylor 展開にベクトルとの積を混ぜる。

計算量は疎行列で O(m·nnz(𝐀))、密行列で O(m·n2) になる。これは n 次正方行列どうしの行列積の計算量より小さい。つまり、Taylor 展開の各回の計算では、保存した前回の結果との積が必要になるが、これは行列とベクトルとの積であり、行列積が不要になる。

手法 B

Krylov 部分空間での近似

部分空間を定義し、そこに対象を圧縮してから計算する。

Krylov 部分空間

上述の近似式は m 個のベクトル 𝐀k𝐛 の線形結合になっている。

これらのベクトルで張られる部分空間は Krylov 部分空間と呼ばれる。

𝒦 m = span { 𝐛, 𝐀𝐛, 𝐀2𝐛, ⋯, 𝐀m-1𝐛 }

この空間の次元数は m で、これは元の問題の n より小さい。

そのため、計算前に座標変換して計算後に戻せば、計算が楽になる。

座標変換と計算

座標変換を挟んだ近似計算は次式で表せる。

exp ( 𝐀t ) 𝐛 ≈ 𝐕m exp ( 𝐇m t ) 𝐞1

ここで、新しく以下の概念を導入している。

  • 𝒦m の正規直交基底: {𝐯1,𝐯2,⋯,𝐯m}
  • 新座標系への m×n 変換行列: 𝐕m={𝐯1,𝐯2,⋯,𝐯m}
  • 新座標系での m×m 指数部行列: 𝐇m=𝐕m⊤𝐀𝐕m
  • 最初の要素が 1 で他が 0 の m 次元ベクトル: 𝐞1

ただし、𝐯1=𝐛/|𝐛| とする。

なお、𝐕m と 𝐇m は Arnordi 法で計算できる。

時間計算量

まず、前処理は Arnordi 法の計算量になる。

次に、本処理は以下の計算なので O(mn) になる。

𝐕m exp ( 𝐇m t ) 𝐞1

なぜなら、𝐕m と exp(𝐇mt)𝐞1 の積は n×m 行列と m 次元ベクトルの積のため、その計算量は O(mn) である。そして、exp(𝐇mt)𝐞1 の計算に手法 A を適用すると、その計算量は O(m2) となるが、これは前者に比べて小さく無視できる。

Arnordi 法

Arnordi 法は Krylov 部分空間の正規直交基底を求めるアルゴリズムである。

計算手順

以下の疑似コードの手順で計算する。 



	v1
	:=
	𝐛/
	
		|
		𝐛
		|
	


	for
	i:=1
	to
	m-1 {

	
		𝒘
		:=
		𝐀𝐯i
	
	
		for
		j:=1
		to
		i {
	
		
			hj,i
			:=
			
				⟨
				𝒘
				,
				𝐯j
				⟩
			
		
		
			𝒘
			:=
			𝒘
			-
			hj,i
			𝐯j
		
	}
	
		
			𝐯
			i+1
		
		:=
		𝒘
		/
		
			|
			𝒘
			|
		
	
	
		
			h
			i+1,j
		
		:=
		
			|
			𝒘
			|
		
	
}

Gram-Shmidt 法との関係

正規直交化の基本的な考え方は、Gram-Schmidt 法と同じである。なお、やや非効率だが、目的の計算は Gram-Schmidt 法でもできる。最初に 𝐕m を求め、そこから 𝐇m を計算すればよい。Arnordi 法では 𝐇m を同時に求められる。

基本方針

後述する『部分空間の別表現』と『基底変換行列』から、以下の性質を利用する。

𝒦m+1 = span { 𝐯1, 𝐀𝐯1, 𝐀𝐯2, ⋯, 𝐀𝐯m }
𝐀 𝐕m = 𝐕m 𝐇m

これらにより、対象の部分空間を 𝐯1 と 𝐀𝐕m の各列で表現すると、その正規直交化により得られる正規直交基底が 𝐕m、基底変換行列が 𝐇m になる事が分かる。

なお、𝐇m が 𝒦m+1 への変換行列となるには次元数が一つ足りない。しかし、最初の列は 𝐯1 から 𝐯1 への変換に対応するため固定値となる。また、最後の行も 𝒦m までの計算には必要ない。そのため、基底変換行列 𝐇m は m+1 次元の上三角行列からそれらを削った Hessenberg 標準系と呼ばれる形式で表せる。

部分空間の別表現

以下の各行はどれも同じ部分空間の別表現である。

𝒦m = span { 𝐯1, 𝐯2, 𝐯3, ⋯, 𝐯m } = span { 𝐛, 𝐀𝐛, 𝐀2𝐛, ⋯, 𝐀 m-1 𝐛 } = span { 𝐯1, 𝐀𝐯1, 𝐀𝐯2, ⋯, 𝐀 𝐯 m-1 }

一行目は 𝐯i の定義から明らかである。

二行目は 𝒦m の定義そのものである。

三行目を補題 P(m) とすると、これは数学的帰納法で証明できる。

数学的帰納法の導入部: P(2) は span{𝐛,𝐀𝐛}=span{𝐯1,𝐀𝐯1} から自明である。

数学的帰納法の連鎖部: P(m)⇒P(m+1) は以下の手順により示せる。

まず、P(m) の仮定により、i≤m で以下が成り立つ。

𝐀i-1 𝐛 ∈ span { 𝐯1, 𝐯2, 𝐯3, ⋯, 𝐯m } 𝐀𝐯i-1 ∈ span { 𝐯1, 𝐯2, 𝐯3, ⋯, 𝐯m }

また、𝐀m𝐛=𝐀·𝐀m-1𝐛 と分解できるため、次が成り立つ。

𝐀m𝐛 ∈ span { 𝐀𝐯1, 𝐀𝐯2, 𝐀𝐯3, ⋯, 𝐀𝐯m-1, 𝐀𝐯m }

これをなるべく 𝐯1 から 𝐯m の線形結合で表すと、次のように変形できる。

𝐀m𝐛 ∈ span { 𝐯1, 𝐯2, 𝐯3, ⋯, 𝐯m, 𝐀𝐯m }

以上より、次が成り立ち、P(m+1) が証明できる。

𝒦m+1 = span { 𝐯1, 𝐯2, 𝐯3, ⋯, 𝐯m, 𝐀𝐯m } = span { 𝐛, 𝐀𝐛, 𝐀2𝐛, ⋯, 𝐀m-1 𝐛, 𝐀m 𝐛 } = span { 𝐯1, 𝐀𝐯1, 𝐀𝐯2, ⋯, 𝐀 𝐯 m-1 , 𝐀𝐯m }

基底変換行列

𝐀𝐕m の各ベクトルは、正規直交基底 𝐕m と基底変換行列 𝐇m で表現できる。

𝐀 𝐕m = 𝐕m 𝐇m

なぜなら、両辺の左側から 𝐕m⊤ を掛けると、右辺の 𝐕m⊤𝐕m は 𝐕m が正規直交基底のため単位行列となり式から消去でき、この表現は 𝐇m の定義と一致する。