丸めモードを使用しない行列積の精度保証
2010SE040 服部亮輝
指導教員：杉浦洋

1

はじめに行列積は線形代数の最も基本的な演算の一つであり，その精度保証は数値計算における重要な課題である． 行列積の精度保証には，丸めモードの切替えを用いた有用なアルゴリズムがあるが，計算環境によっては，丸めモードの切替えが使用できず，丸めモードは丸め込み (四捨五入) に固定されている． ゆえに，丸め込みモードに固定した環境で働く精度保証アルゴリズムが強く望まれる． 本
研究では，行列 A と B の積 AB を数値計算により包含する手法について考える． ここでの「包含」とは，行列積 AB が含まれる区間の下端と上端となる行列や，その中心と半径となる行列を求め，結果を数値計算を用いて厳密に包み込む「区間」を得ることである． 尾崎らの計算手法 [2] について学び，理解を深め，Mathematica 用のアルゴリズムの構成を目指す．

2

精度保証付き計算について

2.1

丸めモードについて

一般に実数は機械実数 F に丸められ，メモリに格納される．IEEE754 では以下の４つの丸めモードがある．ｃ ∈ Rとする.
・丸め上げ（round upward）ｃ 以上の一番小さい浮動小数点数に丸める．△： R −→ F と表す.
・丸め下げ（round downward）ｃ 以下の一番大きい浮動小数点数に丸める．▽：R −→ F と表す.
・丸め込み（round to nearest）ｃに一番近い F の浮動少数点数に丸める．□：R −→ F と表す.
・丸め捨て（round toward 0）絶対値 | ｃ | が小さい F の浮動小数点数で最も近いものに丸める．

3

丸め込みの誤差の一部の言語では、精度保証に有用な IEEE754 が定める有向丸めの機能がないことが知られている． 現在多くの計算機で，丸め込みモードがデフォルトのモードになっており． このモードの演算のみで機能するアルゴリズムが必要である．

4

簡易法

簡易のため，A, B ∈ Fn×n は正方行列とする．nu < 1のとき，次の不等式が知られている．

|AB − fl(AB)| ≤ γn |A||B|, γn =nu
1 − nuまた，この右辺は|A| · |B| ≤fl(|A| · |B|)
(1 − u)n(2)で評価できる．ゆえに|AB − fl(AB)| ≤γnfl(|A| · |B|)
(1 − u)nで評価できる．さらに 1 ≤ n ≤ u−1 − 3 のとき，()
|a||a|≤ fl，
1 − nu1 − (n + 1)u)11(1 + u)n ≤≤，(1 − u)n
1 − nu|A||B| ≤ (1 + u)n−1 fl(|A||B|)

(3)
(4)
(5)

を用いて，

|AB − fl(AB)| ≤ fl

( γ (|A||B|) )n(1 − (n + 2)u)(6)

R は実数の集合，F は浮動小数の集合とする．以下では， で評価できる．これを簡易法と呼ぶ．
IEEE754 規格の F を浮動小数点に含み， F 上の浮動小5 尾崎らの精密法
数演算が IEEE754 算術規格をみたすことを仮定する. こ定数 λ ∈ F をの仮定は，ほとんどのコンピュータ，ワークステーション⌈ log (n + 1) − log u ⌉メインフレームに用いられている．計算表現のためにあ22λ=る計算式について fl(·) は，丸め込みモードの浮動小数演2算で計算された値を表す．この場合 IEEE754 では．と定める．2 つのベクトル v ∈ Fn と w ∈ Fn をそれぞれ
fl(x op y) = (x op y)(1 + δ),(1)
vi = 2λ・2⌈ log2 max1≤j≤n |aij |⌉ , wj = 2λ・2⌈ log2 max1≤i≤m |bij |⌉|δ| ≤ uと定める．ここで e = (1, 1, …, 1)T ∈ Fn を定義し、が成立する．ここで op ∈ {+, −, ×, /} であり，u = 2−53
は丸め誤差単位である．
A(1) = fl((A + veT ) − veT ), A(2) = fl(A − A(1) )3.1

丸めモード切り替えによる精度保証とその問題点

大石，Rump[1] は IEEE754 の丸めモード切り替えを用いた効果的な行列積の精度保証のアルゴリズムを提案した．しかし，FORTRAN77 や Java，Mathematica など

B (1) = fl((B + ewT ) − ewT ),B (2) = fl(B − B (1) )を実行する．以上の結果

A = A(1) + A(2) ,B = B (1) + B (2)が成立する．ここで行列乗算 AB は

AB = (A(1) + A(2) )(B (1) + B (2) )= A(1) B (1) + A(1) B (2) + A(2) Bと、三つの行列積により変換され、A(1) と B (1) に対してが成立する．

AB ∈ ⟨M + T1 + T2 , R(1) + R(2) ⟩ ⊆ ⟨M, |T1 | + |T2 | + R(1) + R(2) ⟩を得る．さらに

fl(|T1 | + |T2 | + R(1) + R(2) )1 − 3u( |T | + |T | + R(1) + R(2) )12≤ fl
1 − 4u:= R|T1 | + |T2 | + R(1) + R(2) ≤fl(A(1) B (1) ) = A(1) B (1)が成立する．

fl(|A||B|)fl(γ n (|A||B|))≤ (1 + u)(1 − u)n(1 − u)nfl(γ n (|A||B|))fl(γ n (|A||B|))
1)·
=
≤(
1−u
(1 − u)n(1 − u)n+1(
fl(γ n (|A||B|))
γ n (|A||B|) )
≤
≤ fl
(7)
(1 − (n + 1)u)
(1 − (n + 2)u)γn |A||B| ≤ γ nを得る．よって，AA(1) B (2) ∈⟨

(1)

B(2)

(2)とA

B の包含を( γ (|A(1) ||B (2) |) )⟩fl(A(1) B (2) ), fl n(1 − (n + 2)u)=: ⟨M (1) , R(1) ⟩⟨( γ (|A(2) ||B|) )⟩nA(2) B ∈ fl(A(2) B), fl(1 − (n + 2)u)=: ⟨M (2) , R(2) ⟩と得る．ここで M (0) をfl(A(1) B (1) ) とすれば AB は下記のように包含される．

AB ∈ M (0) + ⟨M (1) , R(1) ⟩ + ⟨M (2) , R(2)⟩= ⟨M (0) + M (1) + M (2) , R(1) + R(2) ⟩上式の中心と半径を最近点への丸めで計算を行い，最終的にはすべての成分が浮動小数点数である中心と半径の行列をそれぞれ求めなければならない．
ここで，中心の計算を高精度に行うために，浮動小数点の和に関するアルゴリズムを紹介する． Knuth は a, b ∈ Fに対してと計算することにより，最終的な包含AB ∈ ⟨M, R⟩, M, R ∈ Fm×pを得る．

6数値実験

前章のアルゴリズムに従って，数値実験を行った．計算機を FUJITSU のノートパソコン FMV-S8390，CPU はintel Core2 Duo P8700 2.53GHz メモリ 2.00GB である．
OS は Windows7 で Mathematica ver.8.0.1.0 上でプログラムを作成した．今回は，前章で述べたアルゴリズム (尾崎) を Mathematica に実現し、簡易法の製作した精度保証付きアルゴリズムとで Mathematica 上での数値計算の精度と計算速度の比較を行った．
時間について精密法との結果を比較すると n = 256 の時，精密法の結果は簡易法に比べ，約 3000 の時間がかかった． また n の値が上がるにつれ，計算速度も比例して大きくなることが分かる.．
精度において，n の値が小さいときは大きな精度の差は見られないが n = 256 において精密法の計算結果は簡易法の計算結果よりも約 3 桁数値が正確に求められていることが分かる．また n の値が大きくなるにつれ，精密法の数値計算は簡易法の数値計算よりもより精度の高い計算結果を得ることが可能となる．

7終わりに

本研究では，尾崎ら [2] により，成分がすべて浮動小数点数で表現される行列 A と B の積 AB を包含する方法
a + b = x + y, x = fl(a + b)を学び，Mathematica 上でのアルゴリズムの構築を実現となるような x, y ∈ F を求めるアルゴリズムを開発した． した．計算速度において簡易法の方が尾崎よりも高速に数値計算の実行ができた．しかし，精度保証において nすなわち，|a| ≥ |b| となるように並べかえてからの値が大きくなるにつれ，尾崎の数値計算は簡易法よりx = fl(a + b),精密な精度保証ができると分かった．

y = fl(a − x) + b.
[]
このアルゴリズム (尾崎) を x, y = TS(a, b) とし，行列の和に拡張して次のように適用する．
[][]
H1 , H2 = TS(M (0) , M (1) ), H3 , T1 = TS(H2 , M (2) ),[]M, T2 = TS(H3 , H1 ).
この結果

M

(0)+M

(1)+M

(2)

= M + T1 + T2

参考文献
[1] Oishi,S.and Rump,S.M.：「Fast venefication of solutions of matrix equations,Numerische Mathematik」，
90[4]．(2002)．
[2] 尾崎克久，荻田武史，大石進一：「有向丸めの変更を使用しないタイトな行列積の包含方法」．応用数
理,Vol21,No3,PP186〜196(2011)．
[3] 杉浦洋：「数値計算の基礎と応用-数値解析学への入門」．サイエンス社 (1997)．