大規模な三角 Toeplitz 線形方 程式の高速解法とその応用 ○ 安村 修一（法政大学 4 年） 李 磊（法政大学） 日本応用数理学会「行列・固有値の解法とその応用」研究部会 第６回研究会

目次 研究背景 Toeplitz 行列 ・定義 ・連立一次方程式の解法 三角 Toeplitz 行列 ・定義 ・連立一次方程式の解法 性能評価 ・計算量 ・アルゴリズム ・数値実験 応用

研究背景 Toeplitz 行列は様々な分野, 例えば信号処理分野 の線形システムや自己相関行列などで現れる 工学上有用な行列である. Toeplitz 行列は特殊な構造をしているため, 一般 的な行列より計算しやすい性質がある 連立一次方程式や逆行列を高速に解くこと が出来れば応用解析の効率があがる

Toeplitz 行列 n×n の 行列 A が以下の形式であれば Toeplitz 行列 と呼ぶ A の i 行 j 列の要素を A i, j と表すとき以下が成り立つ

Toeplitz 連立一次方程式 A : n×n Toeplitz 行列 b, x : n 次ベクトル Ax = b 直接法 Levinson-Durbin アルゴリズム [2] : O(n 2 ) 最近の研究 （ superfast, asymptotically algorithms ） 高速フーリエ変換を利用した方法 : O(nlog 2 3 n) 三重対角化を利用した方法 [3] : O(nlog 2 2 n) 反復法 前処理付き共役勾配法 (PCG 法 ) 一般化共役残差法（ GCR 法）

三角 Toeplitz 行列 行列 A が t  k =0 （ k =1,2, ･･･ n  1 ）のとき 三角 Toeplitz 行列 と呼ぶ 行列 T は第 1 列の要素のみで表すことができる O

三角 Toeplitz 連立一次方程式 直接法 前進代入 ： O(n 2 ) 多項式の除算から求める方法 ： O(nlog 2 2 n) 今回紹介する手法 [1] 分割統治法＋高速フーリエ変換 : O(nlog 2 n) 1987 年 Z. You, L. Li によって提案された T : n×n 三角 Toeplitz 行列 b, x : n 次ベクトル Tx = b

計算量（補題） T 1 (n) : 三角 Toeplitz 連立一次方程式を解くために必要な計算量 T 2 (n) : 三角 Toeplitz 行列の逆行列を求めるために必要な計算量 T 2 (n) ≦ T 1 (n) ≦ T 2 (n) + O(nlog 2 n) 三角 Toeplitz 行列 T の逆行列 T  1 は三角 Toeplitz 行列 T  1 の全要素は第 1 列によって決定 T  1 の第 1 列をベクトル x’ としたとき Tx’ = (1, 0, ･･･, 0) T を満たす方程式を解けば逆行列が求められる 三角 Toeplitz 連立一次方程式 Tx = b を解く 高速フーリエ変換を使って巡回畳み込みを計算 : O(nlog 2 n) x = T  1 b ： 逆行列とベクトルの積 T 2 (n) ≦ T 1 (n) T 1 (n) ≦ T 2 (n)+O(nlog 2 n) （1）（1）

アルゴリズム（解） n × n の三角 Toeplitz 連立一次方程式： T, b を以下のように 2 k  1 次に分解 三角 Toeplitz 行列の逆行列 三角 Toeplitz 連立一次方程式の解 （2）（2） T 1 : 三角 Toeplitz 行列 T 2 : Toeplitz 行列

アルゴリズム（逆行列） 逆行列を求めるとき b = (1, 0, ･･･, 0) T b 1 = (1, 0, ･･･, 0) T b 2 = (0, ･･･, 0) T (2) 式に代入して整理すると : 逆行列の 1 列目 （3）（3） Toeplitz 行列とベクトルの積（畳み込み）の回数 解を求めるとき : 3 回 逆行列を求めるとき : 2 回

アルゴリズム procedure ttoeplitz(T[1..n], b[1..n], x[1..n]) begin T 1 [1..n/2] ← T[1..n/2]; T 2 [1..n/2] ← T[n/2+1..n]; if n > 2 then ttoeplitz(T 1, NULL, T 1  1 ); else T 1  1 [1] = 1 / T 1 [1]; if b == NULL then begin x[1..n/2] = T 1  1 ; x[n/2+1..n] =  T 1  1 * T 2 * T 1  1 ; else b 1 [1..n/2] ← b[1..n/2]; b 2 [1..n/2] ← b[n/2+1..n]; x 1 [1..n/2] = T 1  1 * b 1 ; x[1..n/2] = x 1 ; x[n/2+1..n] =  T 1  1 * (T 2 * x 1  b 2 ); end (3) 式 逆行列の計算 (2) 式 解の計算 再帰呼び出し

計算量 連立一次方程式・逆行列を求めるために必要な計算量は O(nlog 2 n) を超えない 補題 T 2 (n) ≦ T 1 (n) ≦ T 2 (n) + O(nlog 2 n) （2）式（2）式

性能評価 OS: Microsoft Windows XP Professional CPU: Intel Core2 Duo 2.60GHz RAM: 3.0GB Compiler: Microsoft Visual C （ C 言語） 計算機環境 （ 1 ） t k = exp(  k), k = 0,1, ･･･, n-1 のとき x = 1.0 となるように b を設定し連立一次方程式を解く 三角 Toeplitz 連立一次方程式 Tx = b （ 2 ） t k = cos(k), 0 ≦ k ≦ 2π のとき x = 1.0 となるように b を設定し連立一次方程式を解く 比較手法として前進代入を使い, （絶対値）平均誤差, 最大誤差, 計算時間 [ 秒 ] を測定する

性能評価（結果 1 ） サイズ 従来手法 O(n 2 ) 提案手法 O(nlog 2 n) 平均誤差最大誤差計算時間 [ 秒 ] 平均誤差最大誤差計算時間 [ 秒 ] E E E E E E E E E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E t k = exp(  k), x = 1.0 のときの結果

性能評価（結果 1 ） t k = exp(  k), x = 1.0

サイズ 従来手法 O(n 2 ) 提案手法 O(nlog 2 n) 平均誤差最大誤差計算時間平均誤差最大誤差計算時間 E E E E E E E E E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E , E E E E 性能評価（結果 2 ） t k = cos(k), x = 1.0 のときの結果

性能評価（結果 2 ） t k = cos(k), x = 1.0

考察 計算時間 ・提案手法は従来手法 O ( n 2 ) に比べて, 大規模な連立一次 方 程式でも高速に解を求めることができる. ・比較的規模の小さな問題に対しては、計算にかかる時 間は 従来手法とほとんど変わらない 誤差 ・従来手法とほとんど変わらないが, 問題によっては誤 差に ばらつきが出てしまう. これは畳み込みの丸め誤差の影響であると考えられる

応用 Gauss-Seidel 法に適用 Toeplitz 連立一次方程式 Ax = b A = L+U より よって反復式は procedure gauss_seidel(L[1..n], U[1..n], b[1..n], x[1..n]) begin ttoeplitz(L, NULL, L -1 ); for k ← 0 step 1 until convergence x k+1 = L -1 * (b – U * x k ); end SOR 法に適用 Gauss-Seidel 法の反復式より

応用（数値実験） 以下の式から生成された Toeplitz 行列とベクトルを係数とし Ax = b を解く 三角 Toeplitz 行列の高速解法を SOR 法に適用したもの, 通常の SOR 法 BiCG 法の計算時間と反復回数を測定する サイズ SOR 法（ ω=1.1 ）＋高速解法 SOR 法（ ω=1.2 ） BiCG 法 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ]

応用（考察） 高速解法を適用した SOR 法は他の 2 つに比べて反復回数 が多くなってしまったが, 計算時間は BiCG 法と同じくら いで通常の SOR 法よりはるかに早い結果になった SOR 法を単独で用いるのではなく PCG 法, BiCG 法, GCR 法などの前処理として今回の高速解法を適用した SOR 法 を用いれば更なる高速化が期待できる

まとめ 参考文献 [1] の手法を元に三角 Toeplitz 連立一次方程式を 分割統治法と高速フーリエ変換を使って O(nlog 2 n) で計算 できる高速解法の性能評価を行った 数値実験より提案手法が大規模な三角 Toeplitz 連立一次方 程式に対して高速, かつ, 精度良く計算できることを示し た 提案手法を Gauss-Seidel 法（ SOR 法）に適用することに よって Toeplitz 連立一次方程式を高速に計算できることを 示した 他の手法の前処理として提案手法を使うことによって更 なる高速化が期待できる

参考論文 [1] Z. You, L. Li, The Time Complexity of Triangular Toeplitz Systems, Mathematica Numerica Sinica, 9:3, pp , 1987 [2] Levinson, N, The Wiener RMS error criterion in filter design and prediction. J. Math. Phys., v. 25, pp [3] G. Codevico, G. Heinig, and M. Van Barel, A superfast solver for real symmetric Toeplitz systems using real trigonometric transformations, Department of Computer Science, K.U.Leuven, Report TW 386, Leuven, Belgium, February, 2004