線形方程式の解の精度保証に対する OishiRump 法
2009SE197 中野綾子
指導教員：杉浦洋

1

はじめに行い，計算量を O(n3 ) f lops から O(n2 ) f lops に削減することである [1]．
計算量を削減する方法として，行列 P A の近似 LU 分Ũ + ∼解 P A = L̃Ũ と L̃，Ũ の数値逆行列 L̃+ ∼= L̃−1 ，= Ũ −1+−1を計算し，定理 1 で A = P L̃Ũ とする．このとき本論文では, Shin’Ichi Oishi, Siegfried M.Rump[1] による線形方程式の高速精度保証付き解法を研究し，その有効性を Mathematica による数値実験で検証する. コンピュータでは，高速に数値計算するために，実数を浮動小数点数で近似して計算を行っている．
∥R∥∞ = ∥A+ A − I∥∞ = ∥Ũ + L̃+ P A − I∥∞
本研究では，線形方程式 Ax = b の解を精度保証付き= ∥Ũ + L̃+ (P A − L̃Ũ )∥∞で求める方法を考える．最近，効率的な精度保証付き数+ ∥Ũ + (L̃+ L̃ − I)Ũ ∥∞ + ∥Ũ + Ũ − I∥∞
値計算方法が開発され，元数の大きい方程式に対しても，実用的な精度保証が得られるようになってきた. 本研究となる．
では，その方法を学び，理解を深め，Mathematica 用の第一項は高速な精度保証アルゴリズムの構成を目指す．

2γn ∥ | Ũ + || L̃+ || L̃ || Ũ | e∥∞

精度保証の基本定理［定理 1］方程式 Ax ＝ b の近似解 x̃ と A の近似逆行列
A+ が与えられたとする. 行列 R = A+ A − I が不等式+ δn ∥ | Ũ + || L̃+ | (ne + diag(| Ũ |))∥∞ u,第二項は∥ R ∥∞ < 1γn ∥ | Ũ + || L̃+ || L̃ || Ũ | e∥∞+ nδn ∥ | Ũ + | e∥∞ ||| Ũ | e∥∞ u,が満たされたなら，真の解 x∗ が存在し,∥x∗ − x̃∥∞ ≤∥A+ r∥∞∥A+ ∥∞ ∥r∥∞≤1 − ∥R∥∞
1 − ∥R∥∞(1)が成り立つ．ここで r = Ax̃ − b は残差である．//
式 (1) の右辺を精度保証付きで計算することにより，誤差 ∥x∗ − x̃∥∞ の上界が得られる．重要なのは ∥ R ∥∞ と∥ r ∥∞ の限界の計算である．
低速アルゴリズム∥ R ∥∞ と ∥ r ∥∞ の限界
setround(down);
R = f l(A+ A − I);r = f l(r);setround(up);R = f l(A+ A − I);r = f l(r);rmid = f l(r + r/2);
rrad = f l(rmid − r);

第三項はγn ∥ | Ũ + || Ũ | e∥∞ + δn ∥ne + diag(| Ũ |)∥∞ uでおさえられる. ここで，δn = n/(1 − nu)，γn = uδn ，u = 2−53 ，u = 2−1027 ，e = (1, · · · , 1)T である．これが∥R∥∞ の計算量 O(n2 ) f lops の事前誤差評価である.

Mathematica 用高速精度保証

4

Mathematica の区間演算は非常に遅い. 行列・ベクトルの∞ノルムを事前誤差評価により，区間演算を用いずに求めることを考える.
4.1

ベクトルの 1 ノルム

ベクトル x = (x1 , x2 , …, xn )T ∈ Rn の∞ノルム∥ x ∥1 =n∑| xi |

i=1

ここで，setround(down)，setround(up) は丸めモードのの評価を考える.
切り替えである．しかし，丸めモードの切り替え機能は，(定理１) ベクトル x ∈ Rn の∞ノルムについて，次の評
Mathematica に実装されていないため，数値実験ではこ 価式が成り立つ.
れを等価な区間演算で置き換えた．
1s=M =f l[∥ x ∥1 ].
∥ x ∥1 ≤3 高速精度保証
1 − γn−11 − γn−1
式 (1) の右辺の計算において，計算量の大きい部分は逆行列 A+ の計算と残差 R の計算で，それぞれ計算量はO(n3 ) f lops である．ここで述べる高速精度保証アルゴ
リズムの目的は，∥R∥∞ を直接評価せず事前誤差評価でアルゴリズムは次のようになる.
1. 機械計算 f l[∥ x ∥∞ ] を実行.
2. 評価式 (1) の右辺を区間計算し区間上限を M とする.
☆定理 1 より ∥ x ∥∞ ≤ M である.

表 1 誤差ノルムの厳格な上界n248163264128
2565121024
2048

4.2

上界 (OR 低)
2.13 × 10−14
1.24 × 10−13
5.97 × 10−13
2.79 × 10−12
1.52 × 10−11
1.03 × 10−10
7.68 × 10−10

上界 (OR 高)
2.13 × 10−14
1.24 × 10−13
5.97 × 10−13
2.79 × 10−12
1.52 × 10−11
1.03 × 10−10
7.68 × 10−10
5.92 × 10−9

∥ A ∥∞ = max

1≤i≤mn∑| aij |

j=1

は，上で述べたベクトルの 1-ノルムの評価技法を用いて次のように評価できる．
1. 機械計算 f l[∥ A ∥∞ ] を実行.
2. 評価式 (7) の最右辺を区間計算し区間上限を M とする.
☆ ∥ A ∥∞ ≤ M が成立する.
線形方程式の残差ノルム

線形方程式 Ax = b，A = (aij ) ∈ Rn×n ，b = (bi )，x =(xi ) ∈ Rn ) の残差 r = Ax − b = (ri ) ∈ Rn の∞ノルム∥ r ∥∞ の評価も以下のように行える.
1.| r˜i | ，ρ˜i を機械計算する.
2.(9) を区間計算し, 区間 [Mi ] の上限を Mi とする.
3.M = max Mi を計算する.
1≤i≤n☆ ∥ r ∥∞ ≤ M である.

4.4

∥ R ∥∞ の上界∥ R ∥∞ の上界も次のように計算できる．
1.R̃ = f l[A+ A − I] を計算
2.R̃1 = f l[| A+ || A | +I] を計算
3.M1 = f l[R̃ ∥∞ ] を計算
4.M2 = f l[∥ R̃1 ∥∞ ] を計算
γn+1
5.M = 1−γ1n−1 (M1 + 1−γM2 ) を区間計算.
n+1

5n2
4
8
163264128
25651210242048

時間 (OR 低)
01.50 × 10−2
1.60 × 10−2
1.09 × 10−1
3.43 × 10−1
2.64
2.07 × 101

時間 (OR 高)
1.60 × 10−2
1.50 × 10−2
4.70 × 10−2
1.09 × 10−1
2.03 × 10−1
7.64 × 10−1
3.25
3.50 × 101

時間 (M 高速)2.49 × 10−4
2.19 × 10−4
2.96 × 10−4
2.97 × 10−4
4.83 × 10−4
1.12 × 10−3
3.84 × 10−3
1.94 × 10−2
1.03 × 10−1
7.80 × 10−1
4.66

行列の∞ノルム

行列 A = (aij ) ∈ Rm×n の∞ノルム

4.3

表 2 計算時間

上界 (M 高速)1.88 × 10−14
1.40 × 10−13
9.04 × 10−13
6.45 × 10−12
1.96 × 10−10
3.53 × 10−97.93 × 10−9
6.67 × 10−8
3.01 × 10−7
2.19 × 10−5
1.22 × 10−4

OR 低と OR 高を比べると誤差ノルム上界は等しい．実行時間では，OR 高の方が短く，比は n が大きくなるほ
ど大きい．n=128 のときは 1/6 となる．今回開発した
Mathematica 用のアルゴリズム (M 高速) は，誤差ノルム上界が OR 低，OR 高の約 10 倍大きめにでる．しかし，非常に高速で例えば，n=256 のとき，OR 高に対する実行時間の比が 1/1800 となる．誤差ノルムの上界は，値が小さいほど正確で，大きいほど不正確である．また，時間が短いほど良い．

6終わりに

本研究では，Shin’Ichi Oishi，Siegfried M.Rump[1] により，線形方程式の解を精度保証付きで求める方法を学び，Mathematica 上でプログラムを実装し数値実験を行った．
三輪 [2] は正規数領域で計算が行われると仮定し研究を行った．本研究では，正規数領域と副正規数領域で計算を行われる場合に拡張し，研究を行った．しかし，Mathematicaの数システムを調査し，副正規数が使われていないことが数値実験により分かり，副正規数領域の実験は行わなかった．精度保証付き計算では，Oishi，Rump の低速アルゴリズム，高速アルゴリズムを実装し，数値実験を行った．高速アルゴリズムは，方程式の次元が 128 のとき，約 6 倍高速であったが，きわめて遅い．その原因は，Mathematicaの区間演算の低速性にある．そこで，事前誤差解析に基づき，Mathematica 専用の高速アルゴリズムを開発した．
このアルゴリズムは，誤差ノルムの上界が約 10 倍大きめにでるが，非常に高速で，方程式の次元が 256 のとき，Oishi，Rump の高速版に比べて約 1800 倍高速である．

数値実験

三種類のアルゴリズムを Mathematica 上に実現し，数 参考文献
値実験を行った．n 次 Frank 行列を係数行列とする線形 [1] Shin’Ichi Oishi, Siegfried M.Rump:Fast verification
方程式の解の精度保証を行った．Oishi，Rump の低速版
of solutions of matrix equations, Numer. Math. Vol.
(OR 低)，高速版 (OR 高) と今回開発した Mathematica
90, pp. 755-773(2002).
用のアルゴリズム (M 高速) の近似解の誤差ノルムの厳格[2] 三輪嵩：「線形方程式の精度保証付き解法」．南山大な上界を計算した．その結果を表 1，2 に示す．表 1 は，学数理情報学部数理科学科卒業論文（2010)．
計算された上界の値である．n は方程式の次元である．表2 は，計算時間を示す．