高精度計算を用いたポータブルな精度保証付き数値計算法の研究

(1)

高精度計算を用いたポータブルな精度保証付き数値計算法の研究

Studies on Portable Self-validating Methods with Accurate Computations

2007 年 3 月

早稲田大学大学院理工学研究科

情報・ネットワーク専攻情報数理工学研究

尾崎克久

Katsuhisa OZAKI

(2)

1 _序論

精度保証付き数値計算はここ１０年で急激に発展した分野である．シミュレーションを行う際に，物理現象をモデル化・離散化を経て最終的に数値計算による結果が得られるが，シミュレーション結果が実際の現象と異なる場合に，精度保証付き数値計算は数値計算の部分に問題があるのかないのかを明確にすることができ，問題の解決に貢献できる。また通常の数値計算は，得られた近似解に対しての精度が保証されていないが，近年の精度保証付き数値計算は実用的な時間で行われ，本来の意味で数値計算学を確立していると評価され，現在では応用数理学の一分野としての地位を築いており，連立一次方程式，行列の固有値問題，非線形連立方程式，微分方程式，数値積分等の様々な分野で発展し続けている．本論文ではJavaで用いることができるポータブルな連立一次方程式の数値解に対する精度保証法と計算幾何学に現れる行列式の符号を精度保証付きで計算する高速なアルゴリズムについて述べている．

1.1 はじめに

本研究は精度保証付き数値計算に関連するが，精度保証付き数値計算の原理は区間解析と関係が深い．1950年代にSunagaにより区間演算が開発された[25]．区間により定義された2つの数に対する4則演算の結果を，また区間により包含する理論であり，

非常に包括的な論文である．その後Mooreが区間解析を学位論文のテーマとし，出版するにあたり多くの研究者が区間解析に関する研究を行った[26]．

浮動小数点数システムを定める規格としてIEEE 754が登場した[1]．IEEE 754は

W.Kahanを中心とするグループの提案により1977年から制定作業が開始され，1985

年に承認された．IEEE 754が制定された背景としては，浮動小数点数方式が個々に制定され，規格として統一する必要があると考えられたことが原因である．規格が統一され普及されれば，プログラムを別々の計算機で実行した場合にも同じ計算結果を得るという利点がある．現在ではIEEE 754は多くの計算機に標準搭載されている規格であり，浮動小数点演算の結果に対する有向丸めが定義され，精度保証を実現する良い環境が整っている．

(5)

この有向丸めを用いて，浮動小数点数により定義された区間型の4則演算の結果をまた浮動小数点数による区間で求める機械区間演算が定義された．数値計算により真の解を包含できるために，直接解法により解を求めるアルゴリズム中の4則演算をすべて機械区間演算に置き換えて精度保証が行われた．例としては連立一次方程式に対する区間ガウスの消去法がある．ところで区間演算は1回ずつの計算につき有向丸めの変更を行うために非効率な部分があったが，行列乗算の結果の包含には有向丸めを上向き・下向きに一度ずつのみ変更し，また線形計算に最適化されたBLAS等を利用できるために非常に効率良く包含することが可能となった．

精度保証付き数値計算を実現する環境としてはMATLABのツールボックスである

INTLAB [24]が有名である．区間型の演算が定義され，行列計算等の包含を高速に行

うことができ，また三角関数や指数関数等の超越関数も区間型を引数とし，結果を区間により出力する．さらに結果の表示も正しい桁のみを表示するように設計されている．

本論文では連立一次方程式の数値解に対する精度保証付き数値計算法と計算幾何学に現れる行列式の符号に対して高速に精度保証を行う手法について述べる．まず連立一次方程式については

Ax=b, A∈R^n×n, b∈Rⁿ (1.1)

と表記され，この形式の連立一次方程式は科学技術計算においてよく扱われる．通常は計算機に実装されている高速な浮動小数点演算を用いて式(1.1)の解を求めるが，浮動小数点演算には丸め誤差，桁落ち等の計算誤差が常に起こり，得られた数値解（浮動小数点演算により得られた近似解）がどれだけ正しいかは疑問である．式(1.1)を解くアルゴリズムに対して数値解析を行い，アルゴリズムの安定性が論じられ，条件数を用いた式により数値解が持つ誤差の上限を押さえることができたが，その条件数自身を正確に求めることは浮動小数点数を用いた数値計算ではほぼ不可能である．誤差を理論式を用いて表現するのではなく定量的に図るには，よく区間演算という技法を用いるが，

通常の近似計算に要する時間の数10倍以上の手間がかかり，また計算ごとに区間幅の増加を避けられず，誤差半径は次元の増加とともに指数関数的に増えるために満足な結果を得ることが難しかった．ところが近年，OishiとRumpらはIEEE 754が定める浮動小数点演算が満たす性質を用いて，求解に要する時間のわずか数倍の手間により

(6)

連立一次方程式の数値解の持つ誤差の上限に対する定量的評価を与える画期的な精度保証法を開発した[2]．この論文の中ではさらに，連立一次方程式の数値解に必要な計算量と同じ手間をかければ，精度保証が行える高速な手法についても述べられている．

即ち連立一次方程式の数値解に対する精度保証法は数値解を求める時間の2倍で行えることを意味する．その後，連立一次方程式の数値解に対する精度保証法は新たな手法も誕生した．三角行列の事前誤差評価を用いる新しい誤差公式を作成し，有向丸めの変更を利用して精度保証を行う方法が開発された[3]．この手法は高速に精度保証が可能なときは高速に行い，不可能な場合には時間をかけて精度保証を行う2段階方式となっている．またIEEE 754が定める浮動小数点演算における丸めの向きを変更しないポータブルな手法[4]や悪条件な行列を係数に持つ連立一次方程式の精度保証を行うための手法[5]がある（厳密には悪条件用の手法の基礎理論はRumpにより1980年代には開発されていたが，この時点では論文化はされていなかった．また最近Rumpによる悪条件用の近似逆行列の生成法がなぜ上手く機能するかを示す収束証明がOishi，

Tanabe，Ogita，Rumpらによってなされた[22]）．本論文ではJavaというポータブルでかつ近年は高速化の技術も発展し，数値計算やシミュレーションにも利用される言語を用いて連立一次方程式の数値解に対する精度保証付き数値計算を行う．議論のポイントは，JavaがIEEE 754が定める有向丸めの変更を使用できない点を考慮し，かつ誤差の上限を過大評価しない精度保証法を開発することである．またJava以外に丸めのモードの変更ができない環境があり，例としてはFORTRAN77がある．またライブラリが関数内部で丸めの向きを変更するために，ユーザからの有向丸めの変更命令が無視される環境もある．このような環境では，丸めモードの変更を必要としない最近点への丸めのみを用いた精度保証法が必要である．また将来的に全ての計算機環境において有向丸めが機能するようになったとしても，提案手法はOgita-Rump-Oishiらの精度保証法[4]の誤差半径の過大評価を抑える手法として存在するために，意義を失わない．

次に計算幾何学に現れる行列式の符号を高速かつロバストに求めるアルゴリズムについて述べる．計算幾何学は実社会において幅広く応用を持つ分野であり，凸包・点とポリゴンの内外判定問題・ボロノイ図等はよく議論の対象になる問題である．これ

(7)

らの例を含む応用問題を解くためには，点と直線の位置判定を求めることや点と円の内外判定を用いるなど，基本的な計算の積み重ねが必要である．従来の計算幾何学の問題を解くアルゴリズムは，理論上は高速に実行できるように設計されているが，数値計算による誤差が生じるために正しくない判定結果を利用しながらアルゴリズムが進行するために，正確な結果を返さないことや計算が無限ループに陥ることが問題とされてきた．そこで数値計算に誤差が起こるのを前提に上手く計算が進むアルゴリズムが設計され，数値計算の結果が良いほど求められる結果も比例して良いという方向性を持った優れたアルゴリズムが開発された[38]．

そのために数値計算に求められる期待は高いが，より信頼できる結果を得るために計算の精度を上げると計算時間は膨大に増える．例えば倍精度浮動小数点数の2倍の精度を持つような演算をソフトウェアにより仮想的に実装した場合，通常の浮動小数点演算に比べ30倍から100倍の時間を要すると言われている．そのために多倍長精度の計算を使用する前に，フィルターと呼ばれる浮動小数点演算のみを用いた行列式の符号に対する精度保証法を試すことが効率的である．例えば行列式の符号をLU 分解により求めた後，逆行列を用いて精度保証を行う方式やLU分解自身を区間演算で求める精度保証方式[30], 逆行列を用いない精度保証方式[31]がある．フィルターは問題が悪条件な場合には精度保証に失敗することがあり，その場合には高精度な計算に移行するが，行列式の符号を保証するために必要な精度は問題に依存する．そこでShewchuk は問題を適応的，すなわち問題が簡単な場合には少ない計算時間で，問題が悪条件な場合には難しさに比例した計算時間をかける形式で，非常に優れたアルゴリズムを設計した．行列式の符号が保証されるまで，高精度計算を少しずつ用いていく方式であ

り，またShewchukは，ハードウェアによって最適に実行される浮動小数点演算のみを

用いて高精度に計算を行うアルゴリズムを開発した[14]．その結果，Shewchukは計算幾何学に現れる行列式をロバストに計算する際には，従来の多倍長精度計算ライブラリを用いるよりも浮動小数点演算を用いて高精度に計算を行う方式が高速であることを示した．

数年後にDemmelとHidaは，問題が悪条件な場合に用いられる高精度計算に対して

は，行列式の計算を浮動小数点数の総和を計算する問題に誤差なく変換し，高精度に

(8)

ベクトルの総和を計算するアルゴリズムを適用したほうが，Shewchukのロバストなアルゴリズムより高速であることを示した[13]．

本論文の計算幾何学に現れる行列式の符号を保証する高精度かつ高速なアルゴリズムに関する議論のポイントは，最近開発されたベクトルに対する総和を計算結果の精度を保証付きで求めるアルゴリズムを計算幾何学の計算に特化し，問題が悪条件な場合より高速な形式で提案することである．また，先行研究の優れた結果と合わせて，精度保証が簡単な問題は高速に，難しい問題は先行研究よりも高速に行列式の符号を求めることができるようになった．

1.2 IEEE 754 が定める倍精度浮動小数点数と誤差解析

まずIEEE 754 [1]が定める倍精度浮動小数点数について述べる．倍精度浮動小数点

数は64bitで表現されており，仮数部に52bit，指数部に11bit，符号の表現に1bit使用する．

IEEE 754では，浮動小数点演算において以下の4つの丸めモードを用意している．

• 最近点への丸め(デフォルト)

• +∞方向への丸め

• −∞方向への丸め

• 原点方向への丸め

また2つの浮動小数点数に対する4則演算の際には誤差のない世界で計算してから浮動小数点数に丸めた結果と，浮動小数点演算を行った結果が一致するように規格が定められている．また浮動小数点数に対する平方根も，誤差のない世界で計算してから浮動小数点数に丸めた結果と浮動小数点演算で求めた結果が一致する．ただし指数関数や三角関数などはこのような性質を満たすように規定はされていない．

本論文中では計算はIEEE 754規格に定められた倍精度浮動小数点数を用いる．Fを浮動小数点数の集合とし，fl(·)は括弧中の演算をすべて浮動小数点演算によって実行したときの結果を表す．たとえば，x, y ∈Fⁿに対してfl(x^Ty)は浮動小数点演算によって

(9)

計算した内積の値を意味する．また浮動小数点演算を+∞方向への丸めモードにおいて計算するときはfl₄(·)と表記し，−∞方向への丸めモードにおいて計算するときは fl₅(·)と表記する．有向丸めはベクトルx, y ∈Fⁿに対する内積計算x^TyやA, B ∈F^n×n に対する行列乗算AB 等を下限・上限により包含する際に非常に便利であり，実際には内積・行列乗算の結果は以下のように包含できる．

fl₅(x^Ty)≤x^Ty≤fl₄(x^Ty), fl₅(AB)≤AB ≤fl₄(AB).

これは実数による計算結果は浮動小数点演算による計算を2回行えば包含できることを意味する．上記の計算方式は有向丸めの変更は2回で良く，線形計算ライブラリで

有名なBLASやLAPACKにより最適化された高速な関数をそのまま用いられるので

非常に効率が良い．同じ包含を行う手法として区間演算という手法があるが，区間演算は1回の計算ごとに有向丸めの向きを変更するために，上記の包含方法に比べて効率的ではない．浮動小数点演算において，アンダーフローは起きないと仮定する¹．またu を浮動小数点演算の相対精度(unit roundoff)とする²．

次に，本論文で利用する有向丸めを用いずに誤差の上限を求めるための関係式を以下に紹介する[4]．a, b∈F に対して

|a+b| ≤ (1 +u)fl(|a+b|) (1.2) (1 +u)ⁿ|a| ≤ fl

µ |a|

1−(n+ 1)u

¶

(1.3) が成り立つことが知られている．また，x, y ∈Fⁿ に対して

kxk∞ = fl(kxk∞) (1.4)

|x^T||y| ≤ fl

µ |x^T||y|

1−(n+ 1)u

¶

(1.5)

|x^Ty| ≤ |fl(x^Ty)|+ fl(eγn,2n+2|x^T||y|) (1.6) が成立することが知られている．ただし，m, n∈N に対して eγm,nを

e

γm,n := fl

µ mu 1−nu

¶

(1.7)

1アンダーフローを考慮しても議論は大きくは変わらない．

2IEEE 754倍精度では，u= 2⁻⁵³ となる．本論文では一番大きい次数の項のみを表記する．

(10)

と定義し，行列とベクトルに対する絶対値は，すべての成分に絶対値をつけたベクトル

·行列を意味する．また，ベクトルと行列の最大値ノルムについては，x= (x₁, . . . , x_n)^T 及びA= (a_ij)∈F^m×n に対して

kxk_∞ := max

1≤i≤n|x_i|, kAk_∞ := max

1≤i≤m

Xn

j=1

|a_ij| (1.8)

である．

1.3 連立一次方程式の数値解に対する精度保証法

精度保証付き数値計算の分野では近年，OishiとRumpによって連立一次方程式の数値解に対する高速な精度保証法が提案された[2]．Oishi-Rumpの方法では，式(1.1)の近似解をex，Aの近似逆行列をR，I を n×n 単位行列とすると

kRA−Ik_∞<1 (1.9)

を証明できた場合に

°°A⁻¹°

°∞≤ kRk_∞

kRA−Ik_∞ (1.10)

が成立し

kex−A⁻¹bk_∞≤ kR(Aex−b)k_∞

1− kRA−Ik_∞ (1.11)

という誤差評価式を用いて精度保証を行っている．ここで(1.11)の右辺の計算を前節

で述べたIEEE 754有向丸めを適宜変更しながら計算することで誤差の上限を求める

ことができる．実際にkRA−Ik_∞の上限は

kRA−Ik_∞≤ kmax(|fl5(RA−I)|,|fl4(RA−I)|)k_∞ (1.12) と求められる．浮動小数点演算の回数をflops(floating point operations)と表記すると，

連立一次方程式を精度保証付き数値計算で行う場合は近似逆行列を求めるのに2n³flops 必要であり，また2n³flops必要な行列乗算RAを+∞方向への丸めモードと−∞方向への丸めモードにより計算するために，全体で6n³flopsの計算量が必要となる．係数行列が密行列である場合によく近似解法として用いられるLU分解法に必要な計算量

(11)

は2/3n³flops であることから，精度保証付き数値計算は近似計算の9倍の手間をかければ実行可能なことを意味する．また係数行列に何も制約をつけない場合において最も高速な連立一次方程式の数値解に対する精度保証法はOishiとRumpにより開発された[2]．P ∈F^n×nを置換行列とした軸交換付きLU分解

P A≈LˆUˆ

を利用し，L, Uの近似逆行列をそれぞれXL, XUとしたときに式(1.11)は

°°ex−A⁻¹b°

°∞≤ kX_U ·X_L·P(Aex−b)k_∞

1− kX_U ·X_L·P ·A−Ik_∞ (1.13) と変形できる．ここでL, Uの近似行列をL,ˆ Uˆ としγ_nを

γ_n:= nu 1−nu と定義すると

kX_U·X_L·P ·A−Ik_∞≤γ_2n

°°

°|X_U||X_L||P||L||ˆ Uˆ|

°°

°∞+γ_n

°°

°|X_U||Uˆ|

°°

°∞ (1.14) であることが証明される．上式は三角行列をW，W に対する近似逆行列XW を消去法によって計算された場合には以下の不等式

|X_WL−I| ≤γ_n|X_W||W|

が成立し，この関係式を用いて得られる（三角行列における逆行列計算の事前誤差評価と言われている）．また行列X, Y ∈F^n×nに対しては

k|X||Y|k_∞ =k|X|(|Y|e)k_∞, e= (1,1, . . . ,1)^T ∈Fⁿ

という関係を利用することにより，式(1.14)における右辺の評価をそれぞれ γ_2n

°°

°|X_U||X_L||P||L||ˆ Uˆ|

°°

°∞ =γ_2n

°°

°|X_U|(|X_L|(|P|(|L|(|ˆ Uˆ|e))))

°°

°∞, γ_n

°°

°|X_U||Uˆ|

°°

°∞=γ_n

°°

°|X_U|(|Uˆ|e)

°°

°∞

と評価することにより次元の2乗に比例する計算量で実行できる．この精度保証法は三角行列L, Uに対する逆行列の計算に要するコストが大部分を占め，合計の計算量は

(12)

2/3n³flopsとなるために，数値解を求めるための近似計算と同じコストをもう一度かけると誤差の上限が求められる高速な精度保証法である．

また荻田と大石により，上記の最も高速な精度保証法よりも扱える問題の幅が広い精度保証法も開発された[3]．

°°ex−A⁻¹b°

°∞≤ kX_UX_LP(Aex−b)k_∞

1−(kX_LP A−Uˆk_∞+γ_nkUkˆ _∞)kX_Uk_∞ (1.15) という評価を行う．ここで

°°

°X_LP A−Uˆ

°°

°∞の上限は+∞方向への丸めモードと−∞

方向への丸めモードを用いて

°°

°X_LP A−Uˆ

°°

°∞ ≤

°°

°max(fl₅(|X_LP A−Uˆ|),fl₄(|X_LP A−Uˆ|))

°°

°∞ (1.16) と計算され，精度保証に必要な計算量は8/3n³flopsとなる．この方式は2段階方式で提案されており，はじめに式(1.13)を用いた精度保証を実行し，精度保証が失敗した

場合は式(1.16)を用いて精度保証を行う形式となっている．式(1.13)の方式に必要な

計算量はX_LとX_Uを求めた後はO(n²)flopsのみであるために，失敗した場合の時間のロスも無視できるほどに小さいために有効な2段階方式になっている．

またOgita,Rump,Oishiらにより，有向丸めを用いず，最近点への丸めのみを用いた

浮動小数点演算により誤差の上限を求める精度保証法が開発された．実際にkRA−Ik_∞ の上限は

kRA−Ik_∞≤fl

µα1+ fl(eγ3n+2(α2+ 2)) 1−2u

¶

と計算される．ここでα₁,α₂は

α₁ := fl(kRA−Ik_∞), α₂ := fl(k|R|(|A|e)k_∞)

と計算される．また精度保証式の分子であるkR(Ax−b)k_∞の上限は r_mid := fl(Aex−b), r_rad := fl(eγ_2n+4(|A||x|+|b|))

としたときに，以下の関係式により求められる．

kR(Ax−b)k_∞ ≤fl

µRr_mid+q 1−2u

¶

(13)

ここで上式にあるq，tは q:= fl

µ|R|(t+r_rad) 1−(n+ 3)u

¶

, t:= fl(eγ_n+1|r_mid|)

と定義する．精度保証の計算においては逆行列，行列乗算に時間を要するが，有向丸めを用いる方式が行列乗算を2回行うのに対し，有向丸めを用いない形式方式は行列乗算が1度のみであるので高速である．ただし，有向丸めを用いた手法よりも扱える問題の幅が狭く，誤差の上限を過大評価する傾向にある．

計算量と計算時間は一概には比例しないことが良く知られている．例えば逆行列を求める際に必要な計算量は一般に2n³flopsと言われるが，分割倒置法を繰り返し適用するとおよそn^log²⁷に比例する計算量まで削減できる．ただし再帰アルゴリズムの複雑さやオーバヘッド等が原因で次元が低いときには高速にならないことが知られている．

即ち，flopsで示される量が小さい場合でも，多くの計算時間を要してしまうこともある．ただし，これまでに紹介した連立一次方程式の精度保証法はすべてBLAS等で最適化された高速な関数をそのまま使用するために計算時間は計算量にほぼ比例し，高速に実装が可能な手法である．ここで実際に各手法の計算時間を計測した時間を示す．

計算機環境として，CPUはPentium M 1.7GHz，ソフトウェア環境としてはMATLAB 7.1を用いて以下の手法に要した計算時間を比較する．

• 手法A 式(1.13)による精度保証法[2]．

• 手法B 式(1.15)による精度保証法[3]．

• 手法C 式(1.12)による精度保証法[2]．

表1は上記の手法に要した時間を表し，その結果からほぼ計算量に比例した計算時間がかかることがわかる．

係数行列が悪条件である連立一次方程式に対する精度保証法は太田·荻田·大石·Rump らにより開発された[5]．係数行列が悪条件な場合には，kRA−Ik_∞ < 1を満たすようなAの近似逆行列Rを得ることが難しい．この問題に対し，従来の多倍長精度計算を用いてRを高精度に求めることは可能であるが，膨大な計算時間を要する．また，

kRA−Ik_∞ < 1を満たすために用意する精度は条件数等に依存するために，事前に

(14)

表 1: 各手法に要する計算時間の比較(秒).

n A B C

1000 0.74 2.47 4.53 2000 5.08 18.84 36.39

はわからないことも問題となる．そのために，Rumpは浮動小数点演算のみを用いて kRA−Ik_∞を満たすRを適応的に求めるアルゴリズムを開発した．Rを複数の浮動小数点数の和

R=R1+R2+R3+. . . で表現し，このようなP

R_iを求めるためにkRA−Ik_∞<1を満たすまで以下の計算

• C =RA

• T = inv(C)

• R1,2,3,...,k =T R

を繰り返して計算する．保持する行列の個数kは，行列の条件数の大きさに依存し，適応的に増えていく．ここでC = RAに対しては高精度な計算を用い(Cの計算結果は 1つの浮動小数点数による行列)，R =T Rの計算では高精度な計算を用い，結果を複数のRの和として保存する．太田らはここで用いる高精度な計算を近年開発された浮動小数点数のみを用いた高速かつ高精度な内積の計算法[8]により実装し精度保証を行った．

精度保証を行う際に使用するメモリを削減する方式も提案されている[27]．通常密行列を係数行列とする場合には係数行列，ＬＵ分解による上三角行列・下三角行列，近似逆行列を保持する必要があり，また途中の計算に対してもワークメモリが必要なために，精度保証は係数行列の格納に必要なメモリの数倍必要であることから，大規模な行列を扱うときに問題となる．この手法では係数行列の格納に必要なメモリの約２倍のみのメモリを用いて節約するが，最適化が優れているLebel 3のBLASを可能な限り用いて精度保証を行っており，ブロック化のサイズによる計算時間の比較が論文

(15)

に書かれている．メモリを節約しない場合に比べて速度は低下するが，約3割程度の計算時間の増加により抑えられている．

また連立一次方程式の数値解の1つ1つの要素に対し誤差上限を与える精度保証法，

すなわち成分毎の精度保証法も開発された．成分毎に精度保証を行うための有効な関係式は1980年代にYamamotoにより開発された[6]．式(1.1)の近似解をx，Ae の近似逆行列をR，I を n×n 単位行列とすると

|ex−A⁻¹b| ≤ |R(Aex−b)|+ kR(Aex−b)k_∞

1− kRA−Ik_∞κ (1.17) が成立する．ここでκ≥

Ã _n X

j=1

|G_1j|, . . . , Xn

j=1

|G_nj|

!

, G=RA−Iとする．ただし，ベクトル v = (v₁, . . . , v_n)^T に対して |v| を |v|:= (|v₁|, . . . ,|v_n|)^T と定義し，ベクトルにおける不等式は各成分ごとに不等式が成立することを意味する．また，行列に対しても同様の表記を用いる．この手法に必要な計算量は6n³flopsであり，連立一次方程式の数値解の持つ誤差の上限をノルムにより評価する式(1.11)と同様の計算量で実行できることがわかる．

その後Ogita，Oishi，Ushiroらにより誤差上限の過大評価を抑える精度保証法が開

発された[7]．式(1.1)の近似解をx，Ae の近似逆行列をR，r :=Aex−b，ezをAez =r の近似解とすると，

|ex−A⁻¹b| ≤ |ez|+kA⁻¹k_∞kAez−rk_∞e (1.18) が成立する．ここでe= (1, . . . ,1)^T ∈Rⁿとし，kA⁻¹k_∞の上限は式(1.10)から求める．

式(1.18)における第2項のkA⁻¹k_∞kAez−rk_∞eの値が大きい場合は，絶対値の小さな成分に対しても大きな成分に対しても同じ誤差限界が加えられるために，成分毎に良い誤差限界を与えているとはいえないが，kAez−rk_∞ を小さくできるために，第2項の過大評価を抑えて成分毎に良い誤差限界を与えるように設計されている．計算量は 6n³flopsとなり，式(1.17)に要する計算量とほぼ同じである．

このように連立一次方程式の数値解に対する精度保証法は様々な方式で開発されてきたが，基本的に計算コストと扱える問題の幅にはトレードオフがある．高速な手法は次元や条件数が高い場合には精度保証に失敗しやすくなり，コストをかける手法ほど次元

(16)

10² 10⁴ 10⁶ 10⁸ 10¹⁰ 10¹² 10¹⁴

MethodA MethodB MethodC

Fig. 1: kRA−Ik_∞の上限の比較．

や条件数が高くても精度保証に成功する傾向にある．ここで各手法におけるkRA−Ik_∞ の上限を比較する．

Fig.1は以下の手法

• 手法A 式(1.13)による精度保証法[2]．

• 手法B 式(1.15)による精度保証法[3]．

• 手法C 式(1.12)による精度保証法[2]．

を実行した結果であり，係数行列はMATLABのgalleryにあるrandsvdを用いて，1000 次元の行列を生成した(特異値の分布モードは3とした)．表1からコストが少ない手法ほど，kRA−Ik_∞の上限が1を超えやすく，コストをかける手法ほど，高い条件数の問題まで精度保証が可能となることがわかる．

また係数行列が特殊な構造を持つときには，その特殊性を考慮した高速な精度保証法が提案されている．係数行列が対角優位であれば不動点定理に基づいた手法が提案されている[34]．係数行列Aを対角行列D，対角成分が0である上三角行列A_U，対角

(17)

成分0である下三角行列A_Lを用いてA=D+A_U +A_Lとしたときに，真の解xと近似解xeの間の誤差上限は以下の式により抑えられる．

kex−xk_∞ ≤ 1 1−q

°°ex+D⁻¹(AL+AU)ex−D⁻¹b°

°∞.

ここでq =kD⁻¹(A_L+A_R)k_∞であり，O(n²)flopsで誤差の上限が計算できる(対角優位という仮定よりq <1は満たされ，精度保証は成功する)．

係数行列が対称の場合にLDL^T 分解を用いる精度保証法が開発された[7]．逆行列の最大値ノルムの上限値を，固有値に関する定理から求めて精度保証を行う手法である．

近似計算に要する時間よりも高速に精度保証が行える場合もあることが示されている．

また係数行列が単調行列な場合に用いる精度保証法も開発されている[11]．単調行列の逆行列A⁻¹の最大値ノルムの上限は，s=Aey−e, ey∈Fⁿとしたときに

°°A⁻¹°

°∞ ≤ keyk_∞ 1− ksk_∞ と抑えられる．ここでA⁻¹の最大値ノルムを以下の不等式

kex−xk_∞≤°

°A⁻¹°

°∞kAex−bk_∞

に代入することで，精度保証が可能となる．精度保証が成功するにはksk_∞<1でなければならないためにeyはAy =eの近似解として与えられる．

また従来の精度保証では近似解・近似逆行列の計算に直接解法を用いていたが，反復解法を用いた精度保証法も開発された[28]．近似逆行列を一度に全て求めずに，n回の連立一次方程式の求解に帰着させて近似逆行列を行単位に求める方式でその処理は並列計算にも向いている．

係数行列が対称正定値な行列であれば近似解を求める際に用いるコレスキー分解の

後にO(n²)flopsで精度保証が行える手法が開発された[10]．数値解を得るために必要

な計算量と比べて精度保証に必要な計算量はほぼ無視できるために，あたかも数値計算に要する時間のみで精度保証ができるという優れた手法である．

このように連立一次方程式の精度保証法は係数行列に何も制約をつけない場合・係数行列が特殊な構造を持つ場合・用いるメモリを削減して実行する場合など多角的に議論がされている．

(18)

1.4 浮動小数点数のみを用いた高速な高精度計算

本節では，本論文中に用いる浮動小数点演算を用いた高精度計算について紹介する．

浮動小数点演算により良い精度を持つ結果を得ることの難しさを表す指標として条件数があり，x, y ∈Fⁿとしたときに，ベクトルの総和に対する条件数は

cond(X x_i) =

P|x_i|

|P x_i| と定義され，またベクトルの内積は

cond(x^Ty) = 2|x^T||y|

|x^Ty|

と定義される．例えば，IEEE 754による倍精度浮動小数点数を用いた場合，条件数が 10¹⁶を超える問題では，浮動小数点演算による結果はほとんど正しい桁をもっていないことが問題となる．

計算機によるベクトルの総和や内積といった基本的な計算を高速かつ高精度に行うアルゴリズムについては昔から多くの議論があったが，近年，浮動小数点演算による高速かつ高精度なベクトルの総和及び内積計算をOgita-Rump-Oishi[8]，Rump-Ogita-

Oishi[9]らが開発した．この手法はソフトウェアにより仮想的に浮動小数点数の仮数

部・指数部を拡張して計算する方式ではなく，ハードウェアが通常サポートする浮動小数点演算のみを用いて高精度な計算結果を得る手法である．特徴として

• 仮数部や指数部への直接的な参照を必要としない．

• 計算機環境に依存しない．

• メインループに分岐がないためにコンパイラが最適化しやすい．

• 入力されるデータに対してソートの必要がない．

• 誤差評価に対する証明を行っている．

といった非常に優れた性質を持つ．結果として現在主流として使用されているライブラリよりも高速に実行されることが知られている．ここで，ある自然数Kに対して，

(19)

あたかも使用している浮動小数点数のK倍の精度で計算した結果を通常の精度に丸めたような結果が得られるDotKと呼ばれるアルゴリズム[8]とベクトルの総和に対する計算結果がfaithful³であることを保証するAccSumと呼ばれるアルゴリズム[9]がある．解きたい問題の条件数がはじめからわかる問題であれば，使用する精度を決めて高精度計算を行うことも可能だが，わからない場合には精度を決められない問題がある．AccSumは悪条件な問題に対しても精度が保証されるまで自動的に反復を繰り返す手法であるために，この問題点に対し何も心配はいらず，速度面に関しても優れていることが示されている．これらの高精度な計算法はこれより本論文で述べる手法において核としての役割を担う．

1.5 計算幾何学と行列式

計算幾何学は実社会において幅広く応用される分野である．点とポリゴンの内外判定・凸包・ボロノイ図等は特に話題になるテーマである．従来の計算幾何のアルゴリズムは理論的には確実に問題の解を求めることができるが，計算機による誤差を前提としていないためにアルゴリズムが暴走することもあった．この問題点を解決するために問題を整数演算のみに変換する整数帰着法や，誤差を前提にしてもアルゴリズムが正しく機能する位相優先法などが開発されてきた[38]．応用問題を解く際には，点と直線の位置関係・点と平面の位置関係・点と円の内外判定・点と球面の内外判定問題等はよく基本として用いられるが，それぞれ低次元の行列式に帰着されることが知られている．

例えば2次元平面上における点と直線の位置関係については，2点(ax, ay),(bx, by)から生成される直線と点(c_x, c_y)の位置関係は以下の行列式

¯¯

¯

ax ay 1 b_x b_y 1 c_x c_y 1

¯¯

¯

の符号により判定される．また3次元空間上における点と平面の位置関係では 3点

3誤差を含まない正確な結果を隣り合うどちらかの浮動小数点に丸めたもの．

(20)

(a_x, a_y, a_z),(b_x, b_y, b_z),(c_x, c_y, c_z)による平面と点(d_x, d_y, d_z)の位置関係は以下の行列式

¯¯

¯

ax ay az 1 b_x b_y b_z 1 c_x c_y c_z 1 dx dy dz 1

¯¯

¯

の符号により判定される．また 2 次元平面上における点と円の内外判定では 3 点 (a_x, a_y),(b_x, b_y),(c_x, c_y)により定義される円と点(d_x, d_y) の内外判定は以下の行列式

¯¯

¯

a_x a_y a²_x+a²_y 1 bx by b²_x+b²_y 1 c_x c_y c²_x+c²_y 1 d_x d_y d²_x+d²_y 1

¯¯

¯

の符号により判定される．また3次元空間上における点と球面の内外判定では 4点 (a_x, a_y, a_z),(b_x, b_y, b_z),(c_x, c_y, c_z),(d_x, d_y, d_z)により定義される球面と点(e_x, e_y, e_z)の内外判定は以下の行列式

¯¯

¯

a_x a_y a_z a²_x+a²_y+a²_z 1 b_x b_y b_z b²_x+b²_y+b²_z 1 cx cy cz c²_x+c²_y+c²_z 1 d_x d_y d_z d²_x+d²_y+d²_z 1 e_x e_y e_z e²_x+e²_y +e²_z 1

¯¯

¯

の符号により判定される．これらの計算幾何学の判定問題に対し，符号の精度保証付き数値計算が考えられてきた．通常，行列式を求める数値計算法としてはLU分解を利用する方法が有力であり，LU分解を利用した精度保証法が開発された．例えば[30]

ではAの部分軸交換LU分解を

P A≈LˆUˆ

とする( ˆL，Uˆ はL，Uの近似行列)．またL，U の近似行列をXL,XU としたときの近似逆行列を

R :=XUXLP としたときに

kI−RAk_∞<1

(21)

が示せた場合は，行列式の符号は

sign(det(A)) = det(P)Πsign(ˆu_i,i) と保証されることが示されている．

また上記の手法が提案されている論文には同時に区間演算を使用する精度保証法も紹介されている．部分軸交換付きLU分解法を機械区間演算により求め，U のすべての対角成分の区間に0が含まれていなければ，置換行列の行列式とUの対角成分のすべての積を考慮すれば，行列式の符号は保証される．このような区間演算を使用する手法は与えられる行列の次元が高いときには精度保証に失敗しやすいが，計算幾何学の問題は低次元の行列を扱うことも多いので，応用可能である．

また区間型を用いた別の精度保証法も提案されている[31]．部分軸交換LU分解を P A ≈ LˆUˆ とし，LとU の近似逆行列をX_L，X_U とすると，Aに対する行列式dは

INTLABを用いて以下のステップにより精度保証付き(区間)で求まる．

B =X_L∗intval(P ∗A)∗X_U; M = mid(diag(B));

G= midrad(M,abss(X

(B−diag(M),2)));

d= prod(G)/prod(intval(diag(X_U)))/det(P);

上記の手法は区間ガウスの消去法よりもはるかに高い次元まで精度保証が可能である．

また，LやU の近似逆行列を求めずに精度保証を行う手法も提案され，以下の不等式が成立した場合に行列式の符号が保証される[31]．

|det(U)|> n²e^∗D⁺.

ここでe^∗ =γnmaxi,j(|L| · |U|)，D⁺ = Πⁿ_k=1(kA(:, k)k2+ne^∗)と定義される．

これまで紹介した行列式の符号に対する精度保証付き数値計算法は係数行列が悪条件な場合には精度保証が失敗することがある．ここで係数行列が悪条件であるということは，点と平面の位置関係の判定問題では点が平面上に非常に近い位置に存在する

(22)

ことを意味する．この場合には浮動小数点演算により正しい判定が行われないことがある．そこで，このような悪条件な場合にも精度保証を行う手法が開発されてきた．

線分と線分の交差判定問題に現れる計算を除算のない形式で与え，整数演算を用いて計算に誤差を起こさない方式が提案された[33]．計算に割り算がない場合には，整数演算に必要な精度を用意すれば無誤差で計算ができるために必ず判定は正しくでき

る(幾何的な問題点にも触れ，全体的に矛盾のない手法を提案している)．この手法は

点とポリゴンの内外判定に対する偶奇テスト法に応用され，正しい結果が得られるように提案されている．

はじめから多倍長計算を用いると，通常の浮動小数点演算で間違わない判定ができる問題にまで計算時間をかけるために，非効率な部分があり，問題を適応的に解く手法が開発された[14]．手法を複数用意し(判定問題により用意する手法の数は異なる)，

計算量の少ない精度保証法から順に試していく方式である．簡単に行列式の符号が保証される場合には少ない計算時間で，問題が悪条件な場合には問題の難しさに依存した計算時間をかけて精度保証が実行されるために優れている．そのためにShewchukは独自の浮動小数点数を用いた高精度な計算方法を定義し高速に各手法を作成した．その結果，多倍長精度を用いた計算法よりも高速に実行されることを示している．

またDemmelとHidaは点と平面の位置関係に対して問題が悪条件な場合にShewchuk

の手法より高速に行列式の符号を保証する手法を開発した[13]．与えられる行列式を浮動小数点数の総和に誤差なく変換し，高精度にベクトルの総和を計算するアルゴリズムを適用する方法である．この手法は入力データにソートが必要であり，また拡張倍精度浮動小数点数を利用する．一般にソートは時間のかかる作業であるが，浮動小数点数の指数部によるソートが高速に機能することと，拡張倍精度を用いた計算は近年の計算機環境において高速に実行できるために，優れた手法であることが知られている．

計算幾何学において，入力に誤差がある場合の判定法も開発された[32]．例えば線分と線分の交差判定において，線分を構成する2点の座標それぞれに誤差がある場合にも正しく判定を行うという手法である．ただし，点の座標に誤差半径が与えられた場合には点は誤差半径内を自由に動くために，点の位置により線分が交差する・しな

(23)

いと変わることがある．そのために

• 入力に誤差がある場合にも必ず線分と線分は交差する．

• 入力に誤差がある場合にも必ず線分と線分は交差しない．

• 入力に誤差がある場合には判定不能．

のいずれかを出力する形式で提案されている．また入力誤差がある場合の凸胞やボロノイ図の扱いについても応用されている．

1.6 本論文の構成

この論文では今までに作成してきた論文をまとめる．第1章では第2章・第3章・第 4章での準備を行う．第2章ではJavaで実装可能なポータブルな連立一次方程式の数値解に対する精度保証法について述べる．誤差上限の過大評価を抑制するために，適宜高精度計算を用いた．第3章では連立一次方程式の数値解に対して，従来手法では近似解の絶対値の大きさに差がある場合，絶対値の小さな数値解の成分に対しては良い相対精度を保証することができない点について，ポータビリティを保ちながら改良した手法について述べる．第4章では計算幾何学に現れる「点と平面の位置関係」を判定する高速なアルゴリズムについて述べる．第5章では結論を述べる．

(24)

2 Java による連立一次方程式の数値解に対する高精度な精度保証法

概要

Javaは近年の最適化技術などによって数値計算や数値シミュレーションを行う言語としても注目されてきている．従来の連立一次方程式の数値解の精度保証法は浮動小数点演算の丸めモードの変更を必要とするが，JavaではIEEE 754規格で定められている丸めモードの変更がサポートされていない．そこで本報告では丸めモードの変更を必要としない連立一次方程式の数値解の高速かつ高精度な精度保証法を提案する．最後に，数値実験によって提案手法の有効性を示す．

2.1 はじめに

本章では，連立一次方程式(1.1)における数値解の精度保証をJavaにおいて実行する方法について考える．ここで，数値解とは計算機によって得られる近似解のことであり，精度保証とは厳密解に対する近似解の定量的誤差評価のことを指す．

JavaはOS等に依存しないポータブルな言語であり，一度実装すれば異なるプラットフォームにおいて同一の結果を得ることができる．従来，Javaは数値計算の分野においては処理速度の点で問題視されていたが，近年のJust-In-Time (JIT)コンパイラ

やHot Spot VM等の最適化技術により，その点もかなり改善することが可能となって

きた[36]．したがって，Javaは数値計算や数値シミュレーションを行う言語としても注

目されてきている[16, 35, 37]．しかしながら，一方ではJavaにおける数値計算の様々な問題点も指摘されている[18]．

連立一次方程式の分野において代表的な手法であるOishi-Rumpの方法[2]では「+∞

方向への丸め」及び「−∞方向への丸め」を用いている．ところが，現在のJavaにおける浮動小数点演算の規格では，丸めモードの変更が定められていない．そのため，Java で丸めモードの変更を実装するには，C言語などで丸めモードの変更を行うライブラリを生成し，Java Native Interface (JNI)を用いなければならない．しかしながら，他

(25)

言語で生成したライブラリを用いると，Javaの大きな利点の1つであるソースコードのポータビリティ(アーキテクチャ，OS，コンパイラなどに依存しない)が失われてしまう．Javaのポータビリティを保持するためにはデフォルトの丸めモードである「最近点への丸め」のみで精度保証を行う必要があるため，従来の丸めモードの変更を必要とする精度保証法を用いることができない．

これに対し，浮動小数点演算の事前誤差評価を用いて丸めモードを変更せずに最近点への丸めのみで高速に連立一次方程式の数値解を精度保証する手法をOgitaらが提案している[4]．本論文ではその手法を拡張し，高精度内積計算アルゴリズム[8]を用いた新しい手法を提案する．これによりJavaのポータビリティが損なわれない上に，精度保証が失敗しない限り，従来の手法よりも数値解の持つ精度をできるだけ過大評価することなく保証することが可能となる．本章の最後では，数値実験によって提案手法の有効性を示す．

2.2 準備

本節では，本章で用いるJavaの浮動小数点システム[17]について簡潔に説明をする．

Javaでは浮動小数点数にIEEE 754の単精度及び倍精度の表現形式を用いている．デフォルトの丸めモードはIEEE 754の最近点への丸めモードであり，Java単独では丸めモードの切り替えが出来ない．またIEEE 754は浮動小数点演算の計算過程において，

単精度もしくは倍精度を超えた拡張精度での内部演算を認めている．拡張精度に関してはCPUに依存するため，浮動小数点演算の結果が異なることがある．現在のJava における浮動小数点演算の規格では，この拡張精度で内部演算を行うwidefpモードがデフォルトとなっている．ただしJavaのstrictfpモードを使用すると，拡張精度を用いずにすべてIEEE 754の単精度及び倍精度で演算が行われる．そのためstrictfpモードを用いて得られる演算結果はCPUに依存しないためポータビリティを持つ．

(26)

2.3 高精度演算アルゴリズム

本節では，提案する精度保証法において用いる内積や行列・ベクトル積の高精度演算について述べる．

まず，浮動小数点演算の加算及び乗算に対する誤差の計算方法について説明する．これらは，内積計算や行列・ベクトル積を高精度に計算するアルゴリズムの中で用いる．

Knuthは2つの浮動小数点数の和a+bを近似値xと誤差yの和x+y (x, y ∈F)に変換するアルゴリズムを提案している[19]．このとき，a+b =x+y が厳密に成り立つ．

そのアルゴリズムを以下に記す．本報告中のアルゴリズムはコード記述の簡略化のた

めにMatlabに似たコードで記述する．

アルゴリズム 1 (Knuth [19]) a, b∈Fに対してa+bの近似値xとその誤差y (x, y ∈ F)を求めるアルゴリズム．

function [x, y] = TwoSum(a, b) x= fl(a+b)

z = fl(x−a)

y= fl((a−(x−z)) + (b−z))

2

次にDekkerによる仮数部がpビットの浮動小数点数aに対しs =dp/2eとし上位s ビットに対応するxと下位p−sビットに対応するy (x, y ∈F)に誤差なく分解するアルゴリズム[15]を以下に記す⁴．

アルゴリズム 2 (Dekker [15]) a ∈Fに対して，仮数部がpビットである浮動小数点を上位s =dp/2eビットに対応するxと下位p−sビットに対応するy (x, y ∈F)に分解するアルゴリズム．

function [x, y] =Split(a)

c= fl(factor·a) % factor= 2^s+ 1 x= fl(c−(c−a))

y = fl(a−x)

2

4たとえば，IEEE 754倍精度の場合，p= 53,s= 27となる．

(27)

上記アルゴリズムを利用して，DekkerとVeltkampは2つの浮動小数点数の積a·b をその近似値xと誤差y (x, y ∈ F)の和に変換するアルゴリズムを提案している[15]．

このとき，a·b =x+y が厳密に成り立つ．そのアルゴリズムを以下に記す．

アルゴリズム 3 (Dekker [15]) a, b∈Fに対しa·bの近似値xと誤差y(x, y ∈F)を求めるアルゴリズム．

function [x, y] = TwoProduct(a, b) x= fl(a·b)

[a₁, a₂] = Split(a) [b₁, b₂] =Split(b)

y= fl(a₂·b₂−(((x−a₁·b₁)−a₂·b₁)−a₁·b₂))

2

ここで上記アルゴリズムを実行する際には，IEEE 754で定められた浮動小数点数の規格で厳密に計算される必要がある．CPUによる拡張精度で計算を行うとアルゴリズムが正常に機能しないことがあることに注意をしておく．またコンパイラの最適化により，アルゴリズム中の括弧をはずして計算されると正常に機能しないことにも注意をしておく．

上記アルゴリズムを用いてOgitaらはハードウェアでサポートされている通常の浮動小数点演算，特に，加減算及び乗算のみを用いて内積計算を高速かつ高精度に行う手法を提案している[8]．これはTwoSumやTwoProductにより得られる誤差を計算値に還元していく手法であり，実行レベルでの最適化に優れていることから高速に行えることが示されている．また高精度計算アルゴリズムによって計算された値に関する誤差解析が行われている．その誤差解析の結果から，使用している精度のk倍精度で計算を行い，元の精度に丸めた結果が得られることが知られている．その中で，提案手法で用いるアルゴリズムDot2Errを以下に記す．Dot2Errは，使用している精度に対してその倍の精度(倍精度であれば4倍精度)で内積を計算し，その計算値と誤差の上限を返す関数である．

アルゴリズム 4 (Ogita-Rump-Oishi [8]) x, y ∈Fⁿに対し，内積x^Tyを2倍の精度

(28)

で計算したときの計算値res ∈Fと誤差の上限err ∈Fを求めるアルゴリズム．

function [res,err] =Dot2Err(x, y)

if 2nu≥1, error(’inclusion failed’), end [p, s] = TwoProduct(x₁, y₁)

e=|s|

for i= 2 : n

[h, r] =TwoProduct(x_i, y_i) [p, q] = TwoSum(p, h) t = fl(q+r)

s= fl(s+t) e= fl(e+|t|) end

res= fl(p+s)

δ = fl((nu)/(1−2nu)) α = fl(u|res|+δe) err= fl(α/(1−2u))

2

表2は，各成分が [−1,1]の一様乱数であるベクトル x, y ∈Fⁿ に対して nを10から 10000まで変化させて Dot2Err を用いたときの，相対誤差 err/|res| の値である．数値実験ではJavaのstrictfpモードで倍精度を使用した．

この結果から，内積計算における誤差の上限を相対的に非常に小さく計算できることがわかる．

次に，アルゴリズム 4 (Dot2Err)を行列とベクトルの積に適用できるように拡張し，

A∈F^n×n，x∈Fⁿに対して

y−r≤Ax≤y+r (2.1)

(29)

表 2: 乱数ベクトルに対するDot2Errの相対誤差 n err/|res|

10 1.12e-16 100 1.12e-16 1000 1.12e-16 10000 1.12e-16

を満たすような中心y, 半径r(y, r∈Fⁿ)を求めるアルゴリズムを以下に記す．

アルゴリズム 5 A∈F^m×n, x∈Fⁿ に対しy−r ≤Ax≤y+r を満たす Ax の近似ベクトル y と誤差ベクトルr (y, r ∈Fⁿ)を求めるアルゴリズム．

function [y, r] =MVErr(A, x) [m, n] =size(A) % m行n列 for i= 1 :m

w=A(i,:)

[y_i, r_i] = Dot2Err(w^T, x) end

2

2.4 提案手法による精度保証

本節では，Javaにおいて連立一次方程式の数値解の精度保証をするために，丸めモードの制御なしで行う高速かつ高精度な精度保証法を提案する．

2.4.1 従来手法

Ogitaらは，丸めモードの変更なしで連立一次方程式 Ax=bの数値解の精度保証を

する手法を提案している[4]．この手法は浮動小数点演算の事前誤差評価を用いた手法

高精度計算を用いたポータブルな 精度保証付き数値計算法の研究