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

大規模な三角Toeplitz線形方
程式の高速解法とその応用
○安村 修一（法政大学 4年）
李 磊（法政大学）
日本応用数理学会「行列・固有値の解法とその応用」研究部会 第６回研究会
目次


研究背景
Toeplitz行列
・定義
・連立一次方程式の解法

三角Toeplitz行列
・定義
・連立一次方程式の解法

性能評価
・計算量
・アルゴリズム
・数値実験

応用
研究背景
Toeplitz行列は様々な分野, 例えば信号処理分野
の線形システムや自己相関行列などで現れる工
学上有用な行列である.
Toeplitz行列は特殊な構造をしているため, 一般的
な行列より計算しやすい性質がある
連立一次方程式や逆行列を高速に解くことが
出来れば応用解析の効率があがる
Toeplitz行列
n×n の 行列A が以下の形式であれば Toeplitz行列 と呼ぶ
 t0

 t1
t
A   2

 

 t n1
t 1
t0
t1

t 2
t 1



 




t1
t2
 t (n1) 

 

 
, t0  0

t 1 t 2
t0
t 1 

t1
t0 
A のi 行 j 列の要素を Ai, j と表すとき以下が成り立つ
Ai, j  t i  j
Toeplitz連立一次方程式
Ax = b
A : n×n Toeplitz行列
b, x : n次ベクトル
直接法
• Levinson-Durbin アルゴリズム[2] : O(n2)
• 最近の研究 （superfast, asymptotically algorithms）
高速フーリエ変換を利用した方法 : O(nlog23n)
三重対角化を利用した方法[3] : O(nlog22n)
反復法
• 前処理付き共役勾配法 (PCG法)
• 一般化共役残差法（GCR法）
三角Toeplitz行列
行列A が tk=0（k =1,2,･･･n1）のとき 三角Toeplitz行列 と呼ぶ
 t0
 t
 1
t2

T
 
 

 t n1
t0
O
t1 
  
 t1 t 0
  t 2 t1



　, t0  0



t0 
行列T は第1列の要素のみで表すことができる
T  (t 0 , t1 , , t n1 )T
三角Toeplitz連立一次方程式
Tx = b
T : n×n 三角Toeplitz行列
b, x : n次ベクトル
直接法
• 前進代入 ： O(n2)
• 多項式の除算から求める方法 ： O(nlog22n)
今回紹介する手法 [1]
分割統治法＋高速フーリエ変換 : O(nlog2n)
1987年 Z. You, L. Li によって提案された
計算量（補題）
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)
T1(n) ≦
T2(n)+O(nlog2n)
T1の第1列をベクトルx’としたとき
Tx’ = (1, 0, ･･･, 0)T
を満たす方程式を解けば逆行列が求められる
三角Toeplitz連立一次方程式Tx = bを解く
x = T1b ： 逆行列とベクトルの積
高速フーリエ変換を使って巡回畳み込みを計算 : O(nlog2n)
アルゴリズム（解）
n×nの三角Toeplitz連立一次方程式： Tx  b, n  2k
T, bを以下のように2k1次に分解
T

b 
T  1
b   1 

T2 T1 
 b2 
T1: 三角Toeplitz行列
T2: Toeplitz行列
三角Toeplitz行列の逆行列
T
1
1
1
1






T

T

T
1
　  1
T  1

 
1
1
1 
1  
1 

T1 
T1  
T1 
  T1 T2T1

三角Toeplitz連立一次方程式の解
1

 b1 
T
1

x  T b  
1
1
1  b 
  T1 T2T1 T1  2 
x1


 
　, x1  T11b1
1
  T1 (T2 x1  b2 ) 
1
（2）
アルゴリズム（逆行列）
逆行列を求めるとき b = (1, 0, ･･･, 0)T
b1= (1, 0, ･･･, 0)T
b2= (0, ･･･, 0)T
(2)式に代入して整理すると
x1


x  
　1

  T1 T2 x1 
x　: 逆行列の1列目
（3）
Toeplitz行列とベクトルの積（畳み込み）の回数
解を求めるとき:
3回
逆行列を求めるとき: 2回
アルゴリズム
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)式 解の計算
計算量
T1 (2k )  T2 (2k 1 )  3  O(2k 1 (k  1))
 T2 (2
k 1
)  O(2
k 1
(k  1))
 T1 (2k 1 )  O(2k 1 (k  1))
 T1 (2k 2 )  O(2k 1 (k  1))  O(2k 1 (k  1))
 T1 (2k 2 )  O(2k 1 (k  1))
 
（2）式

 b1 
T11

x  T 1b  
1
1
1  b 
  T1 T2T1 T1  2 
x1


 
　, x1  T11b1
1
  T1 (T2 x1  b2 ) 
補題
T2(n) ≦ T1(n) ≦ T2(n) + O(nlog2n)
 T1 (1)  O(2 1)  O(22  2)    O(2k 1 (k  1))
 1  O(2  22  2    2k 1 (k  1))
 O(k  2k  2  2k  2)
 O(k  2k )
 O(n log 2 n)
T1 (n)  O(n log 2 n)
T2 (n)  O(n log 2 n)
連立一次方程式・逆行列を求めるために必要な計算量はO(nlog2n)を超えない
性能評価
三角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言語）
性能評価（結果1）
tk = exp(k), x = 1.0 のときの結果
サイズ
従来手法 O(n2)
平均誤差
最大誤差
提案手法 O(nlog2n)
計算時間 [秒]
平均誤差
最大誤差
計算時間[秒]
128
1.839E-16
4.441E-16
0.000
2.478E-16
6.724E-16
0.000
256
2.030E-16
4.441E-16
0.000
2.294E-16
6.053E-16
0.000
512
2.125E-16
4.441E-16
0.000
2.435E-16
6.755E-16
0.000
1,024
2.173E-16
4.441E-16
0.000
2.925E-16
9.441E-16
0.000
2,048
2.197E-16
4.441E-16
0.015
2.416E-16
9.287E-16
0.000
4,096
2.209E-16
4.441E-16
0.047
2.756E-16
1.004E-15
0.000
8,192
2.214E-16
4.441E-16
0.156
3.083E-16
1.111E-15
0.000
16,384
2.217E-16
4.441E-16
0.547
3.280E-16
1.119E-15
0.031
32,768
2.219E-16
4.441E-16
2.109
2.843E-16
1.337E-15
0.047
65,536
2.220E-16
4.441E-16
8.297
3.012E-16
1.425E-15
0.094
131,072
2.220E-16
4.441E-16
36.594
3.368E-16
1.362E-15
0.187
262,144
2.220E-16
4.441E-16
207.828
3.515E-16
1.444E-15
0.453
524,288
2.220E-16
4.441E-16
955.125
3.856E-16
1.581E-15
0.969
性能評価（結果1）
1000
1.0E-13
計算時間（従来）
計算時間（提案）
平均誤差（従来）
最大誤差（従来）
平均誤差（提案）
最大誤差（提案）
900
800
1.0E-14
600
誤差
計算時間 [秒]
700
500
400
1.0E-15
300
200
100
0
1.0E-16
7
8
9
10
11
12
13
14
15
16
行列サイズ [2^x]
tk = exp(k), x = 1.0
17
18
19
性能評価（結果2）
tk = cos(k) , x = 1.0 のときの結果
サイズ
従来手法 O(n2)
平均誤差
最大誤差
提案手法 O(nlog2n)
計算時間
平均誤差
最大誤差
計算時間
128
6.465E-15
2.531E-14
0.000
7.409E-15
5.943E-14
0.000
256
2.054E-14
1.097E-13
0.000
1.564E-14
3.363E-13
0.000
512
5.004E-14
2.944E-13
0.000
3.418E-14
7.799E-13
0.000
1,024
1.293E-13
7.607E-13
0.000
6.099E-14
2.592E-12
0.000
2,048
3.724E-13
2.464E-12
0.016
1.161E-13
3.932E-12
0.000
4,096
1.055E-12
7.039E-12
0.031
2.731E-13
3.332E-11
0.000
8,192
2.976E-12
2.287E-11
0.125
4.754E-13
2.337E-10
0.016
16,384
8.562E-12
6.722E-11
0.516
1.430E-12
1.338E-09
0.015
32,768
2.428E-11
1.962E-10
2.312
2.848E-12
5.029E-09
0.047
65,536
6.939E-11
5.142E-10
8.469
5.928E-12
3.444E-08
0.094
131,072
1.961E-10
1.526E-09
38.266
9.913E-12
7.870E-08
0.187
262,144
5.571E-10
4.839E-09
216.750
2.562E-11
9.894E-07
0.437
524,288
1.577E-09
1.313E-08
1036.859
4.288E-11
2.163E-06
0.953
性能評価（結果2）
1200
1.0E+00
計算時間（従来）
計算時間（提案）
平均誤差（従来）
最大誤差（従来）
平均誤差（提案）
最大誤差（提案）
1000
1.0E-02
1.0E-03
1.0E-04
1.0E-05
1.0E-06
1.0E-07
600
1.0E-08
1.0E-09
400
1.0E-10
1.0E-11
1.0E-12
200
1.0E-13
1.0E-14
0
1.0E-15
7
8
9
10
11
12
13
14
15
行列サイズ [2^x]
tk = cos(k), x = 1.0
16
17
18
19
誤差
計算時間 [秒]
800
1.0E-01
考察

計算時間
・提案手法は従来手法O(n2)に比べて, 大規模な連立一次方
程式でも高速に解を求めることができる.
・比較的規模の小さな問題に対しては、計算にかかる時間は
従来手法とほとんど変わらない

誤差
・従来手法とほとんど変わらないが, 問題によっては誤差に
ばらつきが出てしまう.
これは畳み込みの丸め誤差の影響であると考えられる
応用
Toeplitz連立一次方程式 Ax = b
Gauss-Seidel法に適用
A = L+U より
(L  U ) x  b
Lx  b  Ux
x  L1 (b  Ux)
よって反復式は
SOR法に適用
Gauss-Seidel法の反復式より
x  L1 (b  Uxk )
xk 1  xk  ( x  xk )
0   2
xk 1  L1 (b  Uxk )
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
応用（数値実験）
以下の式から生成されたToeplitz行列とベクトルを係数とし Ax = b を解く
1 
ak 
f ( )e ik d , f ( )   4  1
2 
b  (1, ,1)T

三角Toeplitz行列の高速解法をSOR法に適用したもの, 通常のSOR法
BiCG法の計算時間と反復回数を測定する
サイズ
SOR法（ω=1.1）＋高速解法
反復回数
計算時間[秒]
SOR法（ω=1.2）
反復回数
BiCG法
計算時間[秒]
反復回数
計算時間[秒]
1024
171
0.094
144
0.640
119
0.109
2048
171
0.204
144
2.546
122
0.234
4096
171
0.453
144
10.156
124
0.500
8192
171
0.953
144
40.203
124
0.969
16384
171
2.530
144
162.203
125
2.500
応用（考察）

高速解法を適用したSOR法は他の2つに比べて反復回数が
多くなってしまったが, 計算時間はBiCG法と同じくらいで通常
のSOR法よりはるかに早い結果になった

SOR法を単独で用いるのではなくPCG法, BiCG法, GCR法な
どの前処理として今回の高速解法を適用したSOR法を用い
れば更なる高速化が期待できる
まとめ

参考文献[1]の手法を元に三角Toeplitz連立一次方程式を分
割統治法と高速フーリエ変換を使ってO(nlog2n)で計算でき
る高速解法の性能評価を行った

数値実験より提案手法が大規模な三角Toeplitz連立一次方
程式に対して高速, かつ, 精度良く計算できることを示した

提案手法をGauss-Seidel法（SOR法）に適用することによって
Toeplitz連立一次方程式を高速に計算できることを示した

他の手法の前処理として提案手法を使うことによって更なる
高速化が期待できる
参考論文
[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