Download presentation
Presentation is loading. Please wait.
1
大規模な三角 Toeplitz 線形方 程式の高速解法とその応用 ○ 安村 修一（法政大学 4 年） 李 磊（法政大学） 日本応用数理学会「行列・固有値の解法とその応用」研究部会 第６回研究会
2
目次 研究背景 Toeplitz 行列 ・定義 ・連立一次方程式の解法 三角 Toeplitz 行列 ・定義 ・連立一次方程式の解法 性能評価 ・計算量 ・アルゴリズム ・数値実験 応用
3
研究背景 Toeplitz 行列は様々な分野, 例えば信号処理分野 の線形システムや自己相関行列などで現れる 工学上有用な行列である. Toeplitz 行列は特殊な構造をしているため, 一般 的な行列より計算しやすい性質がある 連立一次方程式や逆行列を高速に解くこと が出来れば応用解析の効率があがる
4
Toeplitz 行列 n×n の 行列 A が以下の形式であれば Toeplitz 行列 と呼ぶ A の i 行 j 列の要素を A i, j と表すとき以下が成り立つ
5
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 : 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 法）](http://images.slidesplayer.net/39/11031470/slides/slide_5.jpg "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 法）")
6
三角 Toeplitz 行列 行列 A が t k =0 （ k =1,2, ･･･ n 1 ）のとき 三角 Toeplitz 行列 と呼ぶ 行列 T は第 1 列の要素のみで表すことができる O
7
三角 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
![三角 Toeplitz 連立一次方程式 直接法 前進代入 ： O(n 2 ) 多項式の除算から求める方法 ： O(nlog 2 2 n) 今回紹介する手法 [1] 分割統治法＋高速フーリエ変換 : O(nlog 2 n) 1987 年 Z.](http://images.slidesplayer.net/39/11031470/slides/slide_7.jpg "You, L. Li によって提案された T : n×n 三角 Toeplitz 行列 b, x : n 次ベクトル Tx = b.")
8
計算量（補題） 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）
9
アルゴリズム（解） n × n の三角 Toeplitz 連立一次方程式： T, b を以下のように 2 k 1 次に分解 三角 Toeplitz 行列の逆行列 三角 Toeplitz 連立一次方程式の解 （2）（2） T 1 : 三角 Toeplitz 行列 T 2 : Toeplitz 行列
10
アルゴリズム（逆行列） 逆行列を求めるとき b = (1, 0, ･･･, 0) T b 1 = (1, 0, ･･･, 0) T b 2 = (0, ･･･, 0) T (2) 式に代入して整理すると : 逆行列の 1 列目 （3）（3） Toeplitz 行列とベクトルの積（畳み込み）の回数 解を求めるとき : 3 回 逆行列を求めるとき : 2 回
11
アルゴリズム 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) 式 解の計算 再帰呼び出し
![アルゴリズム 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) 式 解の計算 再帰呼び出し](http://images.slidesplayer.net/39/11031470/slides/slide_11.jpg "アルゴリズム 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) 式 解の計算 再帰呼び出し")
12
計算量 連立一次方程式・逆行列を求めるために必要な計算量は O(nlog 2 n) を超えない 補題 T 2 (n) ≦ T 1 (n) ≦ T 2 (n) + O(nlog 2 n) （2）式（2）式
13
性能評価 OS: Microsoft Windows XP Professional CPU: Intel Core2 Duo 2.60GHz RAM: 3.0GB Compiler: Microsoft Visual C++ 2005 （ 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 を設定し連立一次方程式を解く 比較手法として前進代入を使い, （絶対値）平均誤差, 最大誤差, 計算時間 [ 秒 ] を測定する
![性能評価 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 を設定し連立一次方程式を解く 比較手法として前進代入を使い, （絶対値）平均誤差, 最大誤差, 計算時間 [ 秒 ] を測定する](http://images.slidesplayer.net/39/11031470/slides/slide_13.jpg "性能評価 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 を設定し連立一次方程式を解く 比較手法として前進代入を使い, （絶対値）平均誤差, 最大誤差, 計算時間 [ 秒 ] を測定する")
14
性能評価（結果 1 ） サイズ 従来手法 O(n 2 ) 提案手法 O(nlog 2 n) 平均誤差最大誤差計算時間 [ 秒 ] 平均誤差最大誤差計算時間 [ 秒 ] 1281.839E-164.441E-160.0002.478E-166.724E-160.000 2562.030E-164.441E-160.0002.294E-166.053E-160.000 5122.125E-164.441E-160.0002.435E-166.755E-160.000 1,0242.173E-164.441E-160.0002.925E-169.441E-160.000 2,0482.197E-164.441E-160.0152.416E-169.287E-160.000 4,0962.209E-164.441E-160.0472.756E-161.004E-150.000 8,1922.214E-164.441E-160.1563.083E-161.111E-150.000 16,3842.217E-164.441E-160.5473.280E-161.119E-150.031 32,7682.219E-164.441E-162.1092.843E-161.337E-150.047 65,5362.220E-164.441E-168.2973.012E-161.425E-150.094 131,0722.220E-164.441E-1636.5943.368E-161.362E-150.187 262,1442.220E-164.441E-16207.8283.515E-161.444E-150.453 524,2882.220E-164.441E-16955.1253.856E-161.581E-150.969 t k = exp( k), x = 1.0 のときの結果
![性能評価（結果 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 のときの結果](http://images.slidesplayer.net/39/11031470/slides/slide_14.jpg "性能評価（結果 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 のときの結果")
15
性能評価（結果 1 ） t k = exp( k), x = 1.0
16
サイズ 従来手法 O(n 2 ) 提案手法 O(nlog 2 n) 平均誤差最大誤差計算時間平均誤差最大誤差計算時間 1286.465E-152.531E-140.0007.409E-155.943E-140.000 2562.054E-141.097E-130.0001.564E-143.363E-130.000 5125.004E-142.944E-130.0003.418E-147.799E-130.000 1,0241.293E-137.607E-130.0006.099E-142.592E-120.000 2,0483.724E-132.464E-120.0161.161E-133.932E-120.000 4,0961.055E-127.039E-120.0312.731E-133.332E-110.000 8,1922.976E-122.287E-110.1254.754E-132.337E-100.016 16,3848.562E-126.722E-110.5161.430E-121.338E-090.015 32,7682.428E-111.962E-102.3122.848E-125.029E-090.047 65,5366.939E-115.142E-108.4695.928E-123.444E-080.094 131,0721.961E-101.526E-0938.2669.913E-127.870E-080.187 262,1445.571E-104.839E-09216.7502.562E-119.894E-070.437 524,2881.577E-091.313E-081036.8594.288E-112.163E-060.953 性能評価（結果 2 ） t k = cos(k), x = 1.0 のときの結果
17
性能評価（結果 2 ） t k = cos(k), x = 1.0
18
考察 計算時間 ・提案手法は従来手法 O ( n 2 ) に比べて, 大規模な連立一次 方 程式でも高速に解を求めることができる. ・比較的規模の小さな問題に対しては、計算にかかる時 間は 従来手法とほとんど変わらない 誤差 ・従来手法とほとんど変わらないが, 問題によっては誤 差に ばらつきが出てしまう. これは畳み込みの丸め誤差の影響であると考えられる
19
応用 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 法の反復式より
![応用 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 法の反復式より](http://images.slidesplayer.net/39/11031470/slides/slide_19.jpg "応用 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 法の反復式より")
20
応用（数値実験） 以下の式から生成された Toeplitz 行列とベクトルを係数とし Ax = b を解く 三角 Toeplitz 行列の高速解法を SOR 法に適用したもの, 通常の SOR 法 BiCG 法の計算時間と反復回数を測定する サイズ SOR 法（ ω=1.1 ）＋高速解法 SOR 法（ ω=1.2 ） BiCG 法 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ] 10241710.0941440.6401190.109 20481710.2041442.5461220.234 40961710.45314410.1561240.500 81921710.95314440.2031240.969 163841712.530144162.2031252.500
![応用（数値実験） 以下の式から生成された Toeplitz 行列とベクトルを係数とし Ax = b を解く 三角 Toeplitz 行列の高速解法を SOR 法に適用したもの, 通常の SOR 法 BiCG 法の計算時間と反復回数を測定する サイズ SOR 法（ ω=1.1 ）＋高速解法 SOR 法（ ω=1.2 ） BiCG 法 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ]](http://images.slidesplayer.net/39/11031470/slides/slide_20.jpg "応用（数値実験） 以下の式から生成された Toeplitz 行列とベクトルを係数とし Ax = b を解く 三角 Toeplitz 行列の高速解法を SOR 法に適用したもの, 通常の SOR 法 BiCG 法の計算時間と反復回数を測定する サイズ SOR 法（ ω=1.1 ）＋高速解法 SOR 法（ ω=1.2 ） BiCG 法 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ] 反復回数計算時間 [ 秒 ]")
21
応用（考察） 高速解法を適用した SOR 法は他の 2 つに比べて反復回数 が多くなってしまったが, 計算時間は BiCG 法と同じくら いで通常の SOR 法よりはるかに早い結果になった SOR 法を単独で用いるのではなく PCG 法, BiCG 法, GCR 法などの前処理として今回の高速解法を適用した SOR 法 を用いれば更なる高速化が期待できる
22
まとめ 参考文献 [1] の手法を元に三角 Toeplitz 連立一次方程式を 分割統治法と高速フーリエ変換を使って O(nlog 2 n) で計算 できる高速解法の性能評価を行った 数値実験より提案手法が大規模な三角 Toeplitz 連立一次方 程式に対して高速, かつ, 精度良く計算できることを示し た 提案手法を Gauss-Seidel 法（ SOR 法）に適用することに よって Toeplitz 連立一次方程式を高速に計算できることを 示した 他の手法の前処理として提案手法を使うことによって更 なる高速化が期待できる
![まとめ 参考文献 [1] の手法を元に三角 Toeplitz 連立一次方程式を 分割統治法と高速フーリエ変換を使って O(nlog 2 n) で計算 できる高速解法の性能評価を行った 数値実験より提案手法が大規模な三角 Toeplitz 連立一次方 程式に対して高速, かつ, 精度良く計算できることを示し た 提案手法を Gauss-Seidel 法（ SOR 法）に適用することに よって Toeplitz 連立一次方程式を高速に計算できることを 示した 他の手法の前処理として提案手法を使うことによって更 なる高速化が期待できる](http://images.slidesplayer.net/39/11031470/slides/slide_22.jpg "まとめ 参考文献 [1] の手法を元に三角 Toeplitz 連立一次方程式を 分割統治法と高速フーリエ変換を使って O(nlog 2 n) で計算 できる高速解法の性能評価を行った 数値実験より提案手法が大規模な三角 Toeplitz 連立一次方 程式に対して高速, かつ, 精度良く計算できることを示し た 提案手法を Gauss-Seidel 法（ SOR 法）に適用することに よって Toeplitz 連立一次方程式を高速に計算できることを 示した 他の手法の前処理として提案手法を使うことによって更 なる高速化が期待できる")
23
参考論文 [1] Z. You, L. Li, The Time Complexity of Triangular Toeplitz Systems, Mathematica Numerica Sinica, 9:3, pp282-285, 1987 [2] Levinson, N, The Wiener RMS error criterion in filter design and prediction. J. Math. Phys., v. 25, pp. 261-278. 1947 [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
![参考論文 [1] Z. You, L.](http://images.slidesplayer.net/39/11031470/slides/slide_23.jpg "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,")
Similar presentations
