1 / 23

大規模な三角 Toeplitz 線形方程式の高速解法とその応用

大規模な三角 Toeplitz 線形方程式の高速解法とその応用. ○ 安村 修一（法政大学 4 年） 　 李 磊（法政大学）. 日本応用数理学会「行列・固有値の解法とその応用」研究部会　第６回研究会 . 目次. 研究背景 Toeplitz 行列 ・定義 ・連立一次方程式の解法 三角 Toeplitz 行列 ・定義 ・連立一次方程式の解法 性能評価 ・計算量 ・アルゴリズム ・数値実験 応用. 研究背景. Toeplitz 行列 は様々な分野 , 例えば信号処理分野の線形システムや自己相関行列などで現れる工学上有用な行列である .

hana
Download Presentation

大規模な三角 Toeplitz 線形方程式の高速解法とその応用

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

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

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

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

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

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

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

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

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

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

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

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

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

  13. 性能評価 三角Toeplitz連立一次方程式 Tx = b （1） tk = exp(-k), k = 0,1, ･･･, n-1 のとき x = 1.0 となるように bを設定し連立一次方程式を解く （2） tk = 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++ 2005 （C言語）

  14. 性能評価（結果1） tk = exp(-k), x = 1.0 のときの結果

  15. 性能評価（結果1） tk = exp(-k), x = 1.0

  16. 性能評価（結果2） tk = cos(k) , x = 1.0 のときの結果

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

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

  19. 応用 Toeplitz連立一次方程式 Ax = b Gauss-Seidel法に適用 SOR法に適用 A = L+U　より Gauss-Seidel法の反復式より よって反復式は 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 xk+1 = L-1 * (b – U * xk); end

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

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

  22. まとめ • 参考文献[1]の手法を元に三角Toeplitz連立一次方程式を分割統治法と高速フーリエ変換を使ってO(nlog2n)で計算できる高速解法の性能評価を行った • 数値実験より提案手法が大規模な三角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

More Related