スパース行列における連立一次方程式の精度保証付き数値計算

(1)

平成１６年度卒業論文

スパース行列における連立一次方程式の精度保証付き数値計算

平成１７年２月２日

指導教授 : 大石進一 ^教授

早稲田大学理工学部情報学科

1G01P029-9 川崎文裕

(2)

第 1 _章

序論

(6)

1.1 _背景

現在,多様な現象をモデル化し,その現象を支配する基礎方程式を導き,それを解析していくことによって現象を解析していく研究が多方面でなされている.そのため,連続数学の問題を解くということはとても重要な位置を占めている.ところが,連続数学の問題を解くのは難しい場合が多い.そこで,計算機を利用して連続数学を解く研究が始まり,それにともなって数値計算の理論,研究に対して多くの人が携わりその結果,順調に発展してきたといえる.また今日,コンピュータ技術も大幅な発展を遂げ,いまや少し前の大型計算機がパソコンとして個人が使える時代となった.

しかし,そのような時代になっても数値計算には次のような基本問題が残っている.計算機を用いる数値計算では,それぞれの計算を正確に行わなければならない.しかし,実際には四則演算の結果は丸められ,また無限演算はすべて近似して行われる.つまり,計算結果は常にある程度の誤差を伴っているということである.そこで,計算結果から正しい結論を得るためには,この計算結果に含まれている誤差を評価して真の値の範囲を確定することが必要になる.

従来,数値計算を行っている人々の間では,この点についてあまり取り上げられることがなかった.それは,数値計算によって上記のように数学的に正しいことを保証するためには, 計算コストや手間ひまがかかるといった理由のためであろう.

だが,近年では数値計算の誤差解析の精度を厳密に保証する手法が飛躍的に発達し,多くの場合,近似解の精度保証が近似解を得る手間の数倍程度で可能であることがわかった.また,それに加えて計算機の高速化も伴っていくため,近い将来,数値計算において,任意の精度保証を付加することが日常的になるであろうと考えている.実際に去年度,大学に128台のＰＣクラスタが入り,大規模な行列に対しても高速に計算できるようになっている.

(7)

1.2 _{本論文の目的}

精度保証付き数値計算が重要であることは上記の内容などから明らかである.本論では数値計算の中から連立一次方程式を取り上げ,大規模スパース行列の精度保証付き数値計算とその高速化について述べる.

1.3 本論文の構成

本論文の構成は以下の通りである.

第 2 章では,浮動小数点数システム,区間演算,連立一次方程式の精度保証の定理,アルゴリズムを述べる.

第 3 章では,数値計算ツールとしてMATLABについて述べる.

第 4 章では,スパース行列に対する連立一次方程式の数値計算について述べる.

第 5 章では,実際に数値計算を実行し,検証結果について述べる．

第 6 章では,前章の結果を元とした考察を述べる．

第 7 章では,プログラムのソースを載せる.

(8)

第 2 _章

準備

(9)

2.1 _はじめに

本章では,精度保証つき数値計算システムの基礎となる,浮動小数点システム,区間演算, 連立一次方程式の精度保証の定理について述べる.

2.2 浮動小数点数

現代の数値計算では,浮動小数点数が標準として用いられている.浮動小数点数は現在では標準化が進められ,ほとんどのパソコンやワークステーション,ベクトル計算機,並列計算機などで同じ形式の浮動小数点体系が用いられている.この体系に基づいて,高速に浮動小数点演算が実行される回路がコンピュータに実装されているため,数値計算は浮動小数点数体系上で行われることがほとんどである.そこで,このような浮動小数点数の体系の上での数値計算の精度を保証することが重要である.

そのために,まず浮動小数点数について説明する.浮動小数点数についてはIEEE標準754 ( IEEE standard 754 ,以下 IEEE 754 と略記する. ) がパソコンやワークステーションなどをはじめとして,多くのコンピュータで標準的に用いられている.本論文では,以下 IEEE 754 に基づく 2 進数浮動小数点数システムを考えることにする.

2.2.1 浮動小数点数

IEEE 754 では4 つのタイプの数が用意されている.

¶ ³

• 規格化 2 進浮動小数点数

• 零

• 非規格化 2 進浮動小数点数

• NaN ( Not a Number , 非数)

µ ´

(10)

2 進規格化浮動小数点数

2 進規格化浮動小数点数とは a=±

Ã1 2 +d₂

2² + d₃

2³ +· · ·+d_N 2^N

!

×2^e, (di = 0または1). (2.1) と書ける数をいう. e_min を負の整数,e_max を正の整数として,e は e_min ≤e≤e_max となる整数である.

m = 1 2 +d₂

2² +d₃

2³ +· · ·+d_N

2^N. (2.2)

を符合付き仮数(signed mantissa)といい,eを指数(exponent) という.指数eも2進数で表される.通常,単精度,倍精度(8 byte = 64 bit) ,拡張倍精度(10 byte = 80 bit) の浮動小数点システムがあるが,それぞれつぎのような浮動小数点システムである.

N eの範囲

単精度 24 -126≤e≤127

倍精度 53 -1022≤e≤1023 拡張倍精度 64 -16382 ≤e≤ 16383

ただし,拡張倍精度浮動小数点数は標準化されておらず,コンピュータによって異なる.

規格化 2進浮動小数点システムにおいて表される数の絶対値の最大値は xmax=

µ1 2 + 1

2² + 1

2³ +· · ·+ 1 2^N

¶

2^e^max. (2.3)

であり,その最小値は

x_min = 1

22^e^min. (2.4)

である.倍精度数では(2.3),(2.4)はそれぞれ

(2−2⁵²)×2¹⁰²³ = 1.7976931348623157· · · ×10³⁰⁸, (2.5) 2⁻¹⁰²³ = 2.22507358507201· · · ×10⁻³⁰⁸. (2.6)

(11)

である. |x|> x_max のときにオーバーフロー (overflow) が生じたという.

倍精度浮動小数点数においては,仮数部が 53bit である.

2⁻⁵³= 1.1102230246· · · ×10⁻¹⁶. (2.7) より,倍精度浮動小数点数は 10進数で約 16桁の精度がある.

零

零は規格化されて

+

µ0 2+ 0

2² + 0

2³ +· · ·+ 0 2^N

¶

2^e^min. (2.8)

と表される.

非規格化 2 進浮動小数点数

IEEE 754 では浮動小数点数は指数部がe_min となったとき,仮数部の最初の桁が 1よ

り小さい数を表すために,デフォルトで最初の桁を1とすることをやめ,ここが 0となる数を置くことを許す規格となっている.これを非規格化数 (denormalized number) という.非規格化数の範囲に数が入ることを,漸近アンダーフロー(gradual underflow) という.

このような非規格化浮動小数点数の正で最も小さな数は

2⁻¹⁰⁷⁴ = 4.94065645841246544· · · ×10⁻³²⁴. (2.9) である.これ以下の数になると,アンダーフロー (underflow) が生じたという.

NaN

このほかに, IEEE 754 ではつぎのような特別な数が用意されている.

(a) NaN(Not a Number)は √

−5,∞

∞,+∞+ (−∞) など不当な演算の結果として得

られる.

(b) ±∞ はオーバーフローの結果や零で割った結果として得られる.

(c) ±0はアンダーフローか ±∞での割り算の結果として得られる.

(12)

2.2.2 浮動小数点数と丸め

IEEE 754 では,浮動小数点数の集合F 上での演算は丸目を用いて定義されている.

また,つぎの 4つの丸めのモードが指定できるように,コンピュータが設計されていることを要請している.これは現在のパソコン,ワークステーションを始めとして,ほとんどすべてのコンピュータで実現されている. c を実数(c∈R) とする.

¶ ³

上向きの丸め(round upward)

c 以上の浮動小数点数の中で最も小さい数に丸める.これを 4:R→F と表す.

下向きの丸め(round downward)

c 以下の浮動小数点数の中で最も大きい数に丸める.これを 5:R→F と表す.

最近点への丸め(round to nearest)

c に最も近い浮動小数点数に丸める.これを □:R→F と表す.アルゴリズムでの表記は NEARとする.もし,このような点が2点ある場合には,仮数部の最後のビットが偶数である浮動小数点数に丸める.これを偶数丸め方式 (round to even) という.

切り捨て(round toward 0)

絶対値が c 以下の浮動小数点数の中で, cに最も近いものに丸める.

µ ´

2.2.3 浮動小数点数の性質

丸めの演算を写像として°:R→F と書く.すなわち, ° は 4,5,□ のいずれかと考える. IEEE 754 では,丸めの演算はつぎの条件を満たすが要請されている.

°x=x, (任意のx∈Fについて),

x <=y→ °x <=°y, (任意のx, y ∈Rについて). (2.10)

(13)

また, x∈F のとき,符号を変えることにより, −x や |x| が得られるので,これらは正確に計算される.

IEEE 754 では,つぎの性質が成立する.

□(−x) = −□x, (任意のx∈Rについて), 4(−x) = − 5x, (任意のx∈Rについて),

5(−x) = − 4x, (任意のx∈Rについて). (2.11)

IEEE 754では,浮動小数点数演算(F 上での四則演算)は丸めとの関係により,つぎのよ

うに定義されている. · ∈ {+,−,×, /} , ° ∈ {4,5,□} のとき

x°· y=°(x·y), (任意のx, y ∈Rについて). (2.12) この式は,左辺の浮動小数点数の四則演算の結果 x°· y は,右辺の数学的に正しい (実数としての) 四則演算の結果 x·y を指定された丸めを行って得られた数 °(x·y) に一致するように計算することを表している.

また,平方根も

(√

x)_{f p}=°(√

x), (任意のx∈Fについて). (2.13) と,浮動小数点数演算によって計算された平方根 (√

x)_{f p} は,正確な実数演算で計算された平方根 √

xを指定された丸めの方向へ丸めた数となる.注意すべきことは,指数関数や三角関数などはこのような規格を満たしていない.つまり,初等関数の値を精度保証付きで求めるためには工夫が必要である.

2.3 _区間解析

以前に述べたような精度保証付き数値計算の原理は，実数値で与えられる真の解の上限と下限を浮動小数点演算により計算しようとするものであった．このような考え方は区間解析と呼ばれる数学的理論と密接な関係がある．この区間解析についてまず簡単に述べる．

(14)

区間解析では区間を数の拡張と考える．そして，その四則演算が定義される．これを区間演算という．区間演算を浮動小数点システムの上で展開するものを機械区間演算という．

ここで注意すべきことは，区間演算および機械区間演算はそれを基本演算として直接的に考えると実用上大きな問題を生じることが多いということである．このような問題点を解決するために，前述の精度保証付き数値計算の原理が考えられた．

以下，区間演算について述べていく．

2.4 区間演算

実数上での演算を仮定する区間演算の概念を説明するここでの議論において四則演算は実数の四則演算とする．

区間解析においては，区間とは

[x, x] = {x∈R| x≤x≤x} (2.14)

と表される閉区間であるとする．ただし，x≤x ∈Rで，それぞれの区間の下端，上端とする．以下，記号の節約のため，

[x] = [x, x] (2.15)

と区間を表すことにし，区間[x]と単に書けば，それは[x] = [x, x]を表すものとする．x=x となる区間[x]は点区間という．点区間は実数であるので，点区間[x]の表す実数をxと書くことにする．

区間[x]について以下を定義する．

d([x]) = x−x, r([x]) = x−x

2 , m([x]) = x+x

2 (2.16)

d([x]), r([x]), m([x])をそれぞれ区間[x]の直径，半径，中心という．

区間[x]の最小絶対値と最大絶対値をそれぞれ次で定義する．

h[x]i= min{|x| |x∈[x]}

|[x]|= max{|x| | x∈[x]}= max{|x|,|x|} (2.17)

(15)

二つの区間[x],[y]が与えられたときその二つの区間の四則演算を次で定義する．

[x]◦[y] = {x◦y | x∈[x], y ∈[y]} (2.18) ただし，◦ ∈ {+,−,×, /}．これを区間演算という．この定義では区間演算は無限回の実数 演算をしないと実行できないような印象を与えるが，実は次が成立する．

¶ ³

[x] + [y] = [x+y, x+y]

[x]−[y] = [x−y, x−y]

[x]×[y] = [min{xy, xy, xy, xy},max{xy, xy, xy, xy}]

[x]/[y] = [x]×

"

1 y,1

y

#

, (0∈/[y]) (2.19)

µ ´

これから区間演算が有限回の実数演算で実行できることがわかる．かけ算(表 2.1)と割算(表 2.2)においては，場合分けにより，より少ない手間で計算する次のような計算規則も簡単に導ける．

Table 2.1: 区間[x],[y]の乗算[x][y]

y≥0 y30 y≤0

x≥0 [xy, xy] [xy, xy] [xy, xy]

x30 [xy, xy] [min{xy, xy},max{xy, xy}] [xy, xy]

x≤0 [xy, xy] [xy, xy] [xy, xy]

区間演算については，包含関係における単調性

[x]⊆[x⁰], [y]⊆[y⁰]⇒[x]◦[y]⊆[x⁰]◦[y⁰], ◦ ∈ {+,−,×, /} (2.20) が成立する．

(16)

Table 2.2: 区間[x],[y]の除算[x]/[y]

y≥0 y≤0 x≥0 [x/y, x/y] [x/y, x/y]

x30 [x/y, x/y] [x/y, x/y]

x≤0 [x/y, x/y] [x/y, x/y]

また，加法と乗法に関し交換則と結合則が成立する．

[x]◦[y] = [y]◦[x], ◦ ∈ {+,×}

[x]◦([y]◦[z]) = ([x]◦[y])◦[z], ◦ ∈ {+,×} (2.21) しかし，加法と乗法の逆元は存在しない．すなわち，−[x] = [−x,−x]であるが

0 = [0]⊆[x]−[x] = [x−x, x−x]

1 = [1]⊆[x]/[x] (2.22)

となることがわかる．上式で等号は[x]が点区間のときのみ成立する．

また，分配則も区間演算に対しては成立しない．そのかわり次の劣分配則が成立する．

[x]×([y] + [z])⊆[x]×[y] + [x]×[z] (2.23) この式で等号は例えば区間[y]と[z]が同じ符合をもつときに成立する．

関数f :D ⊂R→Rを領域Dの任意の閉区間の上で連続な初等関数とする．関数fを次のようにしてIR上の関数に拡張することができる．

f([x]) ={f(x) | x∈[x]} (2.24)

IRを実数上の有界閉区間の集合とする．関数f : D ⊂ R → Rが初等関数の合成や四則演算に基づく演算規則で定義されているとき，f をIR上の関数に拡張したい．しかし，

f([x])を厳密に計算することは非常な計算時間を要することがある．そこで，fの演算規則

(17)

を単純に区間演算で置き換えた演算規則によって定義される関数をfの区間拡張といいf_□ で表す．

f([x])⊆f_□([x]) (2.25)

が成立する．区間拡張は一意的に定まるわけではなく，f を定義する数式を同等な他の表現形式に書き直したとき，それらの表現に基づく区間拡張が一般に異なる評価を与えることに注意する．

上式で等号が成立するとき，厳密な包み込みができたという．厳密な包み込みは一般に難しいが，fの数学的な表現式に区間値パラメータが現れず，また，[x]が一度しかその表現に現れなければ，その表現に基づく区間拡張が厳密な包み込みを与えることが知られている．一般に式の表現に[x]が少なく現れれば現れるほど，区間包囲は厳密な包み込みに近くなることが多いことが経験されている．区間[x]の幅d([x])が小さいとき，

q(f([x]), f_□([x]))≤αd([x]) (2.26) となることが知られている．ただし，αは正の[x]によらない定数である．

なお，一般に関数の値域f([a, b])をf([a, b]) ⊂ [c, d]と区間[c, d]で評価するとき，区間

[c, d]のことをf([a, b])の区間包囲という．区間拡張は区間包囲の一種である．

2.5 機械区間演算

前節で説明した区間演算においては，区間の両端の数は実数であるとし，厳密な実数演算に基づいて区間演算が定義された．しかし，浮動小数点数システム上で数値計算するときには，厳密な区間演算は実行できない．そこで，数学的な正解が含まれる区間が演算の結果得られることを念頭に区間演算で与えられた浮動小数点数システムR上に展開する方法を説明する．

IR={[x]∈IR | x, x∈R} (2.27)

(18)

とする．IRの要素を機械区間という．機械区間は区間の定義で上端と下端の数が浮動小数点数で与えられるものをいう．区間[a, b] ∈ IRを区間[c, d] ∈ IRに丸める区間の丸め

♦:IR→IRを

[x]⊂ ♦[x] (2.28)

および，次を満たすように定義する．

♦[x] = [x] (すべての[x]∈IRについて)

[x]⊂[y]⇒ ♦[x]⊂ ♦[y] (すべての[x],[y]∈IRについて)

♦(−[x]) = −♦[x] (すべての[x]∈IRについて) (2.29)

IEEE754上の浮動小数点数システムにおいては

♦[x] = [5x,4x] (2.30)

と定義すれば上記の条件は満たされる．機械区間演算の四則演算を

[x]♦[y] = ♦([x]◦[y]) (2.31)

と定義する．ただし，◦ ∈ {+,−,×, /}である．

具体的には機械区間演算は次のように定義される．

¶ ³

[x] + [y] = [5(x+y),4(x+y)]

[x]−[y] = [5(x−y),4(x−y)]

[x]∗[y] = [5(min{xy, xy, xy, xy),4(max{xy, xy, xy, xy)]

[x]/[y] = [5(min{x/y, x/y, x/y, x/y),4(max{x/y, x/y, x/y, x/y)] (2.32)

µ ´

解を有限回の四則演算の組合せで求める方法(直接解法という)において，浮動小数点数を機械区間に置き換え，四則演算を機械区間演算に置き換えて計算をすれば，真の解を含む区間が得られる．こうして,数値計算により数学的に正しい結論を導き出すために，直接

(19)

解法が知られている問題については, このような置き換えを行なうことが盛んに行なわれた．ガウスの消去法において，途中の計算における数を機械区間に，四則演算を機械区間演算に置き換えたアルゴリズムを区間ガウス消去法という．区間ガウス消去法はこのようなアルゴリズムの代表的なものである．

2.6 中心と半径による区間

区間演算におけるかけ算は,場合分けも多く複雑である.これをシンプルに計算するために区間の中心と半径による表示を考える.まず,点をかける区間の計算が簡単にできることを示す.αを実数とする.区間を[β−δ, β+δ]と表す.βは区間の中心でδ >= 0は半径である.

このとき¶ ³

α[β−δ, β+δ] = [αβ − |α|δ, αβ+|α|δ] (2.33)

µ ´

が成立する.ただし,四則演算は厳密にできるとする.実際

min{αβ−αδ, αβ+αδ}=αβ− |αβ|=αβ − |α|δ

max{αβ−αδ, αβ+αδ}=αβ+|αβ|=αβ+|α|δ (2.34) より式(2.33)が示された.

(20)

2.7 _{連立一次方程式}

以下は,連立一次方程式の近似解の精度保証を行うための原理となる定理である.

¶ ³

方程式Ax=bの近似解x˜とAの逆行列の近似行列Rが求められたとき,行列G=RA−I が不等式

||G||<1 (2.35)

を満たすときは,逆行列A⁻¹が存在し

||A⁻¹||<= ||R||

1− ||RA−I|| (2.36)

およびx˜を真の解として

||x^∗−x||˜ <= ||R(Ax˜−b)||

1− ||RA−I|| (2.37)

が成り立つ.

µ ´

2.7.1 精度保証アルゴリズム

A˜x−bに対する正確な区間を求める方法を示す.

setround(down);

res=A˜x−b; %区間の下端 setroun(up);

res=A˜x−b; %区間の上端

よって

res≤A˜x−b ≤res (2.38)

(21)

次に,RA-IとAx˜−bの正確は区間を区間の中心res_midと半径res_radから求める.

setround(down) %Rとresの下方向丸め G=RA−I

res=A˜x−b

setround(up) %Gとresの上方向丸め

G=RA−I res=A˜x−b

res_mid = ((res+res)/2) %中心と半径のかたちに変換 res_rad = (res_mid−res)

よって,

G=RA−I ∈ {G, G}

res_mid−res_radı <=A˜x−b <=res_mid+res_rad (2.39) 次は,||x^∗−x||˜ _∞に対する正確な区間を求める.

G= max{|G|,|G|}

setround(down) %下方向丸め計算

4=R∗resmid

setround(up) %上方向丸め計算

4=R∗res_mid

4= max(|4|,|4|) +|R| ∗resrad

G_norm= max

i=1,2,···,n

Xn j=1

G_ij

(22)

R_norm = max

i=1,2,···,n

Xn j=1

|R|_ij 4_norm = max

i=1,2,···,n{4_i} D=−(Gnorm−1) if (D >0) err_norm = 4_norm

D else

printf(”vrif icationf ailed”) end

R∗res_mid− |R| ∗res_rad <=R(A˜x−b)<=R∗res_mid+|R| ∗res_rad

となることを示している。

|R(Ax˜−b)|<= max(|4|,|4|) +|R| ∗resrad <=4. (2.40)

2.8 むすび

本章では，浮動小数点システム，区間演算,機械区間演算,連立一次方程式の近似解の精度保証アルゴリズムについて述べた.

(23)

第 3 _章

数値計算ツール

MATLAB

(24)

3.1 _はじめに

本章では，大規模スパース行列の計算実験に必要な数値計算ツールを紹介する.

3.2 MATLAB

MATLABは、テクニカルコンピューティングのための高性能プログラミング言語であ

る.MATLABは、問題やそれに対する解をよく知られた数学の記法で表現し,計算,可視化, プログラミングなどを利用しやすい環境で統合する. if文やfor文といた制御構造を持ち, コマンド入力することも可能である.なお,基本となるエンジンはＣ言語で作成されているため,インタプリタとしての実行速度も決して遅くなく,快適に仕事をこなすことができる.

典型的な利用形態としては、次のようなものがある.

• 数学と計算

• アルゴリズムの開発

• モデリング、シミュレーション、プロトタイピング

• データ解析、データ補間、データ可視化

• 科学、工学でのグラフィックス

• グラフィカルユーザインタフェース構築などのアプリケーションの開発

3.2.1 MATLAB の特徴

• 行列操作が簡単

MATLABは行列操作が容易に行える.配列宣言は不要で,行列の拡大や一部取り出しが簡

単に行える.また,+,-,*などの演算子はそのまま行列の四則演算に使用することができる.

• グラフィック機能が充実

(25)

さまざまな２次元,３次元グラフ出力関数を持っており,MATLABはアスキーセーブされたファイルからデータを読み込むことができるので,実験等で得られたデータを表示するためにも使用できる.さらに,アニメーションもサポートされている.

• ライブラリが充実

制御系の解析や設計に必要な関数をまとめたtoolboxと呼ばれるライブラリが充実している. toolboxとは,MATLAB関数を広く集めたもので,特定分野の問題を解くのにMATLAB の環境を拡張したものである.

3.3 _むすび

本章では,数値計算ツールMATLABに関して述べた.本論の実験ではMATLABを利用し精度保証実験を行った.

(26)

第 4 _章

スパース行列に対する連立一次方程式の

精度保証付き数値計算

(27)

4.1 _はじめに

本章では,連立一次方程式の精度保証とそのMATLABでの扱いについて述べる.

4.2 連立１次方程式

連立一次方程式

Ax=b (4.1)

の解の精度保証付き数値計算法を考える.ただし,Aはn×n行列,x,bはn次元ベクトルである.行列Aと右辺ベクトルbが与えられたとして,ベクトルxを求める.

4.3 matlab

matlabでAx=bを解く.

4.3.1 連立一次方程式の直接解法

連立一次方程式の直接解法というと,多くの場合,ガウスの消去法に基づいたLU 分解法を指す.この方法は係数行列が特異でなければ適用可能であるという意味で凡用性が高い.

直接解法は,必要とする計算量と計算容量の点から,比較的小規模な問題に適用されることが多い.

MATLABにおける直接解法は,実は大変使いやすく設計されており,

>> x =A\b;

とするだけで,簡単に解を求めることができる.

(28)

しかし,Aが大規模スパース行列であるような連立１次方程式を解くとすると,逆行列を直接求めることによって解を求めることは,あまり望ましいとはいえない.なぜなら,係数行列Aがスパースであっても,その逆行列がスパースであるとは限らず,逆行列が密行列となる結果,莫大なメモリを使ってしまうので,コストパフォーマンスが良くない.

更に,直接解法の場合,計算が完了するまで解の予想が困難であるため,解に関する何らかの情報を得られるまでに必要な計算量は常に一定（問題サイズの3乗オーダ）であり,得られる計算解の精度は高い傾向にあるが,ユーザが常にそこまでの精度を求めているかどうかは疑問で,せいぜい数桁正しい解が得られればよいという場合でも,やはり同じ計算量が必要である.

一般的に直接解法は大規模スパース系には向いていないと考えられている.

4.3.2 スパース行列

スパース行列とは,大部分がゼロ要素である行列の特別なクラスのことである.この行列の特性によって,MATLABではつぎのようなことができる.

¶ ³

1. 行列の非ゼロ要素のみをインデックスと共に格納 2. ゼロ要素への演算を省略して計算時間を短縮

µ ´

4.3.3 スパース行列ストレージ

フル行列の場合,MATLABではすべての行列要素を内部に格納する.ゼロ要素部分も,その他の行列要素部分と同じストレージ容量を必要とする.ただし,スパース行列の場合は,上記のとおり,非ゼロ要素とそのインデックスのみを格納する.この時,ゼロ要素が含まれる割合の高い大きな行列に対しては,データストレージに必要なメモリ総量をかなり削減できる.

(29)

4.3.4 連立一次方程式の反復解法

反復解法とは,ある決められた一連の手続きを何度も繰り返すことによって解の精度を上げていく方法である.反復解法の長所の一つは,係数行列のスパース性を保持できることである. さらに直接解法に対し,反復解法では繰り返しの過程において解の制度をある程度は推測できるため,それに応じて適当なところで計算を打ち切ることができる.

しかしながら,数値計算における反復解法の最大の弱点は,数学的に解の収束が保証されていたとしても,係数行列の性質が悪いと丸め誤差によって必ずしも計算解が真の解に収束するとは限らないことである. 対称正定値行列などの一部の特殊な行列を除いては,決定的なアルゴリズムをつくることは困難とされている.したがって,問題に応じてアルゴリズムを変える必要がある,というのが現状である.

4.4 MATLAB _{で使用した反復解法}

連立一次方程式の反復解法は,ガウス・ザイデル法やSOR法（逐次過剰緩和法）を代表とする定常反復法系統と,CG法（共役勾配法）を代表とする非定常反復法系統の大きく二つに分類できる.MATLABには現時点（virsion7.0）では,連立方程式のために九つの反復解法が用意されているが,それらはすべて非定常反復法である.

今回取り上げる反復解法は,PCG法,Bi-CGSTAB法,GMRES法である.

4.4.1 PCG 法

PCG法は係数行列Aが対称正定値の場合に最も効果的な解法である.ただし,与えられた係数行列が正定値かどうかを判定するのは困難であるため,あらかじめ係数行列が正定値になることが分かっている問題に有効である.

MATLABでは関数pcgが用意されており,反復法の収束においては,同値で条件の良い

方程式に変換してから解くと,速く収束することが多い.この操作を前処理という.前処理行列としてAを近似する行列MあるいはM=M(1)M(2)のように下三角行列M(1)と上三角

(30)

行列M(2)を用意し,Ax =bの代わりにM⁻¹Ax =M⁻¹bを解く.このように条件の良い方程式に変換してから解く前処理付き共役勾配法をPCG法という.

4.4.2 Bi-CGSTAB 法

Bi-CGSTAB法はランチョス原理に基づく反復解法で,係数行列Aが非対称行列の場合に

よく使われる解法である.Bi-CG法系統の解法は,後述のGMRES法系統の解法と違い,リスタートなどのパラメータ設定が基本的に不要であるため使い勝手が良く,ユーザに好まれる傾向がある.

MATLABでは,関数bicgstabが用意されており,基本的な使い方は関数pcgと同様であ

る.Bi-CGSTABの前処理法はいくつかあるが,今回はILU分解（不完全LU分解）を用いている.

本論での実験では,主にbicgstab法を用いた.

4.4.3 GMRES 法

GMRES法はアーノルディ原理に基づく反復解法で,Bi-CGSTAB法と同様に,係数行列A

が非対称行列の場合によく使われる解法である.オリジナルのGMRES法は計算量及び計算容量の点で実用的でなく,適当な正整数kに対して,k回ごとにリスタートすることによって必要な記憶容量を減らしたGMRES(k)法が実際に用いられる.

MATLABでは,関数gmresが用意されており,基本的な使い方は関数bicgstabと同じで

ある.ただし,リスタートのためのパラメータ設定が必要である.

リスタートを行うGMRES(k)法において,kの値は本来,大きいほどオリジナルのGMRES 法に近づき収束も速くなる傾向にあるが,その分必要なメモリ量も増加する.これは,解きたい問題に応じて対応すれば良い.

(31)

4.5 _前処理

ここでは,反復解法の前処理として,不完全LU分解を用いることにする.前処理は,係数行列Aの条件を良くして問題を解きやすくするための技法である.

4.6 外部データを MATLAB で使う方法

大規模スパース行列のうち,MATLABであらかじめ用意されているものではなく,他のアプリケーションなどを使って生成した行列のデータをファイルに保存しておくことで,MAT- LAB上でその行列をスパース型行列として利用することができる. 外部データの取り込みは,パスの設定のすんでいるフォルダ(例えばwork)に当該データを保存しておいて,

>>load filename.mtx

とし,MATLAB上に読み込む.さらに,

A = spconvert(filename);

とすれば,スパース型の行列Aが生成される.

4.7 _精度保証

MATLAB上でメモリに積まれているのは,係数行列Aと近似解x˜とベクトルbである.

この3つを利用し,真の解を見積もる事を考える. 具体的には,真の解x^∗と数値計算による近似解x˜との誤差を見積もる. この為にいくつか,数学的な準備を行う. まず,Ax = bの逆行列が存在しないことには誤差解析もできない.

逆行列が存在するか否かは,左近似逆行列Rとすると,

||RA−I||_∞ <1 (4.2)

(32)

であれば逆行列が存在(右逆行列ならAR)し,

||A⁻¹||∞ <= ||R||∞

1− ||RA−I||_∞ (4.3)

ここで,||˜x−x^∗||_∞については,

||˜x−x^∗||_∞<=||A⁻¹A˜x−A⁻¹b||_∞ (4.4)

=||A⁻¹(A˜x−b)||∞ (4.5)

<=||A⁻¹||_∞||A˜x−b||_∞ (4.6) である.以上により,

||˜x−x^∗||_∞<= ||R(A˜x−b)||_∞ 1− ||RA−I||∞

(4.7)

によって真の解と数値計算によって求められた解の誤差を見積もることができる. よって, 左または右逆行列を求めれば,それによって精度保証を行うことができる. したがって,Rが必要となるが,係数行列が大規模行列の場合,Rを求めると,前述のように, 莫大なメモリを消費するため得策ではない.ここで,||R||∞とは

||R||_∞ = max

||x||∞=1||Ax||_∞= max

i=1,2,···,n

Xn j=1

|R_ij|

であるから,Rの行ベクトル単位で分かれば（時間はかかるが）計算することができ,メモリの消費も抑えることができる.Rのi行の行ベクトルy_iとし,単位行列Iのi行目の行ベクトルをe_iとすると

y_iA=e_i (4.8)

A^ty_i^t=e^t_i (4.9)

つまり,y_i^tさえ求めればよい.

(33)

4.7.1 y

^t_i

の求め方

y_i^tの求め方は,matlab上に用意されている,反復解法を使用することにする.

matlabでは,多くの反復解法が用意されており,下記プログラム中funkに反復解法を指定

することにより,指定の反復解法を使う.ただし,gmresを指定した場合はリスタートの設定をする必要がある.

4.7.2 近似逆行列 R について

上記のように,Aの逆行列の行ベクトルをそのまま求める方法もあるが,Rがよい近似逆行列でなくとも,とりあえず,||RA−I||_∞<1であれば,逆行列の存在を示す事ができるので, あまり解の精度を要求しない場合には,||RA−I||_∞ <1を満たすようなRを場合に応じて選んで計算すればよい.

4.8 むすび

本章では,連立一次方程式のMATLABによる精度保証付き数値計算について述べた.

(34)

第 5 _章

結果

(35)

5.1 _はじめに

本章では，大規模スパース行列系に向いた反復解法の実験結果について述べる．

5.2 実行環境

1. CPU : Intel(R) Xeon(TM) CPU 3.06GHz 2. Memory : 1.0 GB

3. OS : Windows XP 4. MATLAB 7.0

5.3 使用するスパース行列

今回は,テスト行列として,Matrix Marketにある行列を使用することにした. 実験したテスト行列の一覧を下に示す.いずれもSPARSKIT Collectionにあるものである.

5.3.1 表に関する準備

下記,表中のtoliluについて説明する.

• tolilu

前処理(前章参照)としての不完全LU分解においてf illinの基準となるしきい値を

表す.この値を小さくすればするほど,完全なLU分解に近づくが,f ill_inがその分増加して,必要なメモリが莫大になる.

データはすべて反復解法bicgstabで実験しているので表には加えないものとする. bicgstab を使用した目的としてはリスタートのパラメータ設定が基本的に不要である事、非対称行列に対しても利用できる事、データが他の方法に比べ比較的安定かつ高速にとれた点がある.

(36)

5.3.2 実験に用いた行列

T able5.1からT able5.4に示す.

Table 5.1 実験に使用したDRIVCAVセットに用意されている行列について

行列名次元数 tolilu 近似解の反復数近似解算出時間[s] 非ゼロ要素

E20R0000 4241 1.00E-03 31 4.984 131412

E20R1000 4241 1.00E-06 2 88.578 131430

E30R0000 9661 1.00E-06 2 65.281 305794

E30R1000 9661 1.00E-06 2.5 261.125 306002

E40R0000 17281 1.00E-06 2.5 521.921 553216

E40R0100 17281 1.00E-06 2.5 580.547 553562

Table 5.2実験に使用したDRIVCAV OLDのセットに用意されている行列について

CAVITY16 4562 1.00E-06 39.5 5.625 137887

CAVITY17 4562 1.00E-03 53 7.047 131735

CAVITY18 4562 1.00E-06 49 6.937 138040

CAVITY19 4562 1.00E-03 52.5 6.984 131735

CAVITY20 4562 1.00E-03 83.5 9.562 138040

CAVITY21 4562 1.00E-06 83 9.906 131735

CAVITY23 4562 1.00E-06 150.5 15.188 131735

CAVITY24 4562 1.00E-06 155.5 16.515 138040

CAVITY25 4562 1.00E-06 260 25.093 131735

CAVITY26 4562 1.00E-06 2.5 72.016 138040

(37)

Table 5.3実験に使用したTOKAMAKのセットに用意されている行列について

UTM3060 3060 1.00E-06 704 20.547 42211

UTM5940 5940 1.00E-09 7.5 74.281 83842

Table 5.4実験に使用したFIDAPのセットに用意されている行列について

FIDAP014 3251 1.00E-09 5.5 2.297 65747

FIDAP015 6867 1.00E-09 1 69.172 96421

FIDAP018 5773 1.00E-09 4.5 1.594 69231

FIDAP019 12005 1.00E-12 1 3.796 259561

FIDAP031 3909 1.00E-06 7.5 0.782 91165

FIDAP035 19716 1.00E-09 6.5 3.828 217972

FIDAP037 3565 0.001 3 0.36 67591

FIDAPM08 3876 1.00E-06 18.5 1.578 86793

FIDAPM09 4683 1.00E-06 2 1.64 93733

FIDAPM11 22294 1.00E-06 38.5 49.265 617874

FIDAPM29 13668 1.00E-06 2.5 35.047 183394

FIDAPM37 9152 1.00E-09 17.5 17.219 765944

(38)

5.3.3 手順

上に述べた行列を使用するときの手順を示す.Ax = bにおける右辺ベクトルbはすべての要素が1のベクトルと係数行列Aとの積とする. そうすることで, 求める近似解x˜を ones(n,1)ベクトルとして指定する. また,前処理としてはLU分解を用い,fill-inの基準となるしきい値を10⁻³から始め,状況に応じて10⁻⁶,10⁻⁹,10⁻¹²としていく.反復解法はbi-cgstab 法を利用.最大反復回数をloopで指定する.これらの条件で解が収束しない場合,そこで計算を中止する.

以下,具体的なコマンドを示す.

>>load filename.mtx

>>A = spconvert(filename);

>>n = length(A) ; b = A*ones(n,1);

>>clear filename

これにより得られた数値解に精度保証をかける.

>>err=vlin_test(A,b,’bicgstab’,1e-12,1e-3,loop)

5.3.4 誤差

表において,誤差(error)となっているのは

||R(A˜x−b)||_∞

1− ||RA−I||_∞ (5.1)

である.

5.3.5 結果

T able5.5からT able5.8に示す.

(39)

Table 5.5 DRIVCAVセット行列の精度保証結果

行列名 nnz_L nnz_U 精度保証の反復数精度保証の時間[s] 誤差

E20R0000 646885 902660 32100 2440 1.74E-08

E20R1000 853176 986803 2643 685.266 8.97E-09

E30R0000 3068851 3209257 5336 5040 5.98E-09

E30R1000 2915776 3416928 5680 5090 5.20E-09

E40R0000 7310358 7750614 9605 22200 8.21E-10

E40R0100 6909118 8143949 10907 22700 2.94E-08

nnz=非ゼロ要素数

Table 5.6 DRIVCAV OLDセット行列の精度保証結果

CAVITY16 954590 1046092 2481 747.063 1.53E-08

CAVITY17 687110 926318 44100 3210 7.98E-09

CAVITY18 948634 1046760 2481 747 8.96E-08

CAVITY19 654852 942218 63588 4430 6.91E-09

CAVITY20 941387 1061353 2480 752 5.91E-08

CAVITY21 646484 921030 112000 7430 9.98E-08

CAVITY23 877799 1004653 3114 778.484 1.25E-07

CAVITY24 926996 1083525 2530 751.563 2.57E-07

CAVITY25 875442 991848 3461 846.578 2.22E-07

CAVITY26 918345 1068615 2580 756.907 4.04E-10

(40)

Table 5.7 TOKAMAKセット行列の精度保証結果

UTM3060 373282 416872 2280 255.344 4.79E-08

UTM5940 1512912 1574201 4332 1410 1.46E-06

Table 5.8 FIDAPセット行列の精度保証結果

行列名 nnzL nnzU 精度保証の総反復数精度保証時間[s] 誤差

FIDAP015 193593 367392 6862.00 747.844 3.37E-04

FIDAP018 209548 292409 19760.00 946.516 1.30E-03

FIDAP019 523293 486137 6013.50 1610.700 7.36E-05

FIDAP031 210115 248711 1985.50 221.500 1.02E-09

FIDAP035 607184 662599 33091.00 5200.500 6.83E-04

FIDAP037 37396 38695 2675.00 97.531 3.25E-13

FIDAPM08 397410 525357 2161.50 365.281 2.02E-10

FIDAPM09 232706 431041 2362.00 362.375 1.69E-11

FIDAPM11 3974801 5081888 39035.00 11341.000 1.50E-08

FIDAPM29 697698 868566 7993.00 2567.000 1.76E-11

FIDAPM37 1860738 1907389 4749.50 3954.800 4.95E-08

(41)

5.4 _むすび

本章では,MATLABによる連立一次方程式(大規模スパース行列)の精度保証の実験結果をまとめた.

次章では実験の考察を述べる.

(42)

第 6 _章

統括

(43)

6.1 _はじめに

本章では，前章の結果について述べる．

6.2 実行結果について

6.2.1 実行時間

今回の実験では条件を揃えるために主に反復解法bicgstab法を用いた.しかし,各行列によって適した他の反復解法があるはずである.それは,DRICAV OLDセットの行列を見たら分かるように,次元数,非ゼロ要素数が似通っているにも関わらず,実行時間に違いが表れていることから明白である.そういった意味では,行列毎に相性のよい反復解法を試すことにより計算の高速化を図ることができる.それらは今後の展望として捉えておきたい.

6.2.2 次元数

前章の結果の通り,最大で22294×22294の行列(FIDAPM11)を係数行列Aとする精度保証を実現した. これには今回の実験において,メモリの節約を優先したことと,実行環境であるPCのメモリが1GBと(一台のPCとしては)大きかったことが関わっているだろう.

6.2.3 精度保証

全体的に誤差レベルとして10⁻³から10⁻¹³と,高い精度で近似解の誤差を求めることができた. 精度として差が表れるのは,行列と反復解法の相性や,不完全LU分解によるしきい値を10⁻³としてはじめたことで分解が完全LU分解と離れてしまい,行列が特異に近くなってしまったのではないかと考える.

(44)

6.2.4 展望

本実験で最大となった行列は2万次元を超える程度のものであったが,さらに大きな次元の行列の精度保証を実現したい.それには第一に実行環境の改善があげられる.序論でも述べたように,大学に導入された128台のPCクラスタを利用すれば,数万,数十万行といった行列の精度保証も可能となるであろう.

6.3 むすび

本論文では,スパース行列における連立一次方程式の精度保証付き数値計算について述べた.

(45)

第 7 _章

ソース

(46)

7.1 vlin test.m

function err = vlin_test(A,b,func,tol1,tol2,loop) system_dependent(’setround’,0);

if isequal(func,’gmres’)

funcin = ’(P*A,P*b,10,tol1,loop,L,U)’;

else

funcin = ’(P*A,P*b,tol1,loop,L,U)’;

restart = 0;

end

nnz_A = nnz(A) tic;for ii=2:5

switch ii case 1

tolilu = ’0’

case 2

tolilu = 1e-3 case 3

tolilu = 1e-12

end% incomplete LU factorization [L,U,P] = luinc(A,tolilu);

% solve Ax=b

[x,flag,relres,iter_sol] = eval([func funcin]);

if flag == 0

nnz_L = nnz(L) nnz_U = nnz(U) break

elseii, flag endclear L U end

t_sol = toc iter_sol

%return

% verification

disp(’** verification **’)

tic;[err,iter_ver] = vlinsparse(A,b,x,func,tol2,loop);

t_ver = toc iter_ver

(47)

7.2 vlinsparse.m

function [err,iter_ver] = vlinsparse(A,b,x,func,tol2,loop) [m,n] = size(A);

if m ~= n

disp(’rectangular matrix’) return

end

if isequal(func,’gmres’)

funcin = ’(B,P*e,10,tol2,loop,L,U)’;

else

funcin = ’(B,P*e,tol2,loop,L,U)’;

end

% incomplete LU factorization [L,U,P] = luinc(A’,1e-3);

B = P*A’;

system_dependent(’setround’,-Inf);

rl = A*x - b;

system_dependent(’setround’,Inf);

ru = A*x - b;

rmid = 0.5*(rl + ru);

rrad = rmid - rl;

G_norm = 0;

Rr_norm = 0;

e = zeros(n,1);

iter_ver = 0;

ii = 2;

for i=1:n

system_dependent(’setround’,0);

e(i) = 1;

% solve A’y=e

[v,flag,relres,iter] = eval([func funcin]);

iter_ver = iter_ver + iter;

if flag ~= 0 flag

ii = ii + 1;

switch ii case 2

tolilu = 1e-12 otherwise

err = Inf;

error(’impossible!’)

(48)

end

clear L U P B

[L,U,P] = luinc(A’,tolilu);

nnz_L = nnz(L) nnz_U = nnz(U) B = P*A’;

i = i - 1;

else

system_dependent(’setround’,-Inf);

tg_l = v’*A - e’;

yl = v’*rmid;

system_dependent(’setround’,Inf);

tg_u = v’*A - e’;

tg = max(abs(tg_l), abs(tg_u));

norm_tg = sum(tg);

yu = v’*rmid;

yu = max(abs(yl),abs(yu));

norm_rr = yu + abs(v’)*rrad;

if norm_tg >=1 i, norm_tg err = Inf;

error(’norm_tg is larger than 1.’) end

G_norm = max(G_norm, norm_tg);

Rr_norm = max(Rr_norm, norm_rr);

e(i) = 0;

end end

if G_norm < 1

d = -(G_norm - 1);

err = Rr_norm/d;

else

% disp(’verification failed’);

err = Inf;

end

system_dependent(’setround’,0);

(49)

謝辞

(50)

この卒業論文の完成にあたりたくさんの方々から多くの助言や指導を頂きました.

まず,本研究を進めるに当たり,終始丁寧な御指導及び御激励を賜り,その他多くの面でも色々と御面倒を見て下さり御助言を与えて下さいました大石進一教授に深く感謝いたします.

お忙しい中でも私の質問に熱心に答えて下さいました,理工学研究科客員講師,丸山晃佐氏,荻田武史氏に深く感謝いたします.

研究室内や合宿など,日常的な場面でも大変お世話になりました,大石研究室博士課程1 年,尾崎克久氏,並びに大石研究室修士課程2年,大山博毅氏,森山敦史氏,大石研究室修士課程1年,坂内太郎氏,島本誠氏,関本竜平氏,徳永克久氏,福地健氏,細田幸裕氏,横山大輔氏に深く感謝いたします.

最後に,研究内容だけでなく日常でも意見の交換や協力などをして下さいました,大石研究室4年, 有滝貴広氏,岩田揚平氏,大成高顕氏,栗原崇氏,佐藤友規氏,平野晃氏,ト詣各氏,水島裕氏,吉岡毅氏,吉田昂央氏に感謝いたします.

(51)

参考文献

(52)

1. MATLAB HOME PAGE http://www.cybernet.co.jp/matlab/

2. MATRIX＿MARKET

http://math.nist.gov/MatrixMarket/

3. 大石進一:“応用数学セミナー数値計算”,裳華房, (1999), 4. 大石進一:“Linux数値計算ツール”,コロナ社,(2000), 5. 大石進一:“精度保証付き数値計”,コロナ社,(2000), 6. 大石進一:“MATLABによる数値計算”,培風館,(2000),

7. 山下栄吉:“マイクロ波シュミレータの基礎”,電子情報通信学会,(2004),

8. S.Oishi,S.M.Rump:“Fast verification of solutions of matrix equations”,Numer,Math 90 (2002),pp755–773,

スパース行列における 連立一次方程式の 精度保証付き数値計算

平成１６年度 卒業論文