安定化手法の発展形 ISCZ 法における
シンボルリストの新しい評価法
伊井誠和
$*$東邦大学大学院理学研究科
TOMOKAZU II
GRADUATE SCHOOL
OFSCIENCE, TOHO UNIVERSITY
白柳潔
$\dagger$東邦大学理学部
KIYOSHI
SHIRAYANAGI
FACULTY OF SCIENCE, TOHO UNIVERSITY
1
はじめに
安定化理論[1, 2] は、近似計算で実行すると不安定となるアルゴリズムに対し、 それを変形して近似計算 で実行しても誤差の影響を抑制し、安定な出力が得られるようにするための理論である。 本論では、 その安 定化手法の発展形であるISCZ
法[3,
4, 5] で用いるシンボルリストの評価方法について提案する。2節で、 安定化理論とISCZ法について復習する。 3 節で、 提案するシンボルリストの評価方法を、 4 節で、 この評 価方法をISCZ法のユークリッドの互除法アルゴリズムに適用した実験結果を、 計算時間とメモリ使用量に ついて分析し、 本手法の有効性について検討する。2
復習
2.1
安定化理論
次のアルゴリズムを対象に安定化理論の復習を簡単に行う。 詳細はを参照されたい。 $\bullet$ データは、すべて多項式環$R[X_{1}, \cdots, x_{m}]$ の元からなる。$R$ は実数体の部分体である。 $\bullet$ データ間の演算は、$R[X_{1}, \cdots, x_{m}]$ 内の加減乗である。 $\bullet$ データ上の述語は、 不連続点をもつとすればそれは$0$のみである。 $\[email protected] $\dagger$ [email protected]述語の不連続点が$0$ という意味は、
If
$C=0$ ”then.
.
. else.
. .
のように、値が$0$ か否かによって分岐が 別れることである。従って、$C=0$の代わりに$C>0$ や$C\leq 0$などでもよい。上記クラスのアルゴリズム を、不連続点$0$の代数的アルゴリズムと呼ぶ。 ほとんどの数式処理のアルゴリズムはこのクラスに入るか、 このクラスのアルゴリズムに変換可能である。 さて、 安定化の3つのポイントは、 $\bullet$ アルゴリズムの構造は変えない。 $\bullet$ データ領域において、ふつうの係数を区間係数に変える。 $\bullet$ 述語の評価の直前で、 区間係数のゼロ書換えを行なう。 である。すなわち、安定化されたアルゴリズムは次のようになる。区間領域データ領域は区間係数多項式の集合。
区間係数は $[A, B]$ なる形で、$A,$$B\in \mathbb{R},$ $[A, B]$ は集合$\{r\in \mathbb{R}|A\leq r\leq B\}$ を意味する。
区間演算二項演算$\star\in\{+, -, \cross, \div\}$ に対し、
$[A, B] \star[C, D]=[\min(A\star C, A\star D, B\star C, B\star D), \max(A\star C, A\star D, B\star C, B\star D)]$
ゼロ書換え不連続点$0$ をもつ述語を評価する直前で、各区間係数$[A, B]$ に対し、
$A\leq 0\leq B$ ならば$[A, B]$ を $[0, 0]$ に書き換えよ。
そうでないならばそのままとせよ。 今、入力 $f\in R[x_{1}, \cdots, x_{m}]$ を
$f= \sum_{i_{1},\cdots,i_{m}}r_{i1\cdots i.m}x_{1}^{i_{1}}\cdots x_{m}^{i_{m}}$
と表したとき、$f$ に対する近似列
Int(f) 乃を
Int
$(f)_{j}= \sum_{i_{1},\cdots,i_{m}}[(a_{i1\cdots im})_{j}, (b_{i1\cdots im})_{j}]x_{1}^{i_{1}}$ $x_{m}^{i_{m}}$で定義する。 ここに、すべての$i_{1},$
$\cdots,$$i_{m}$ について、
$(a_{i1\cdots im})_{j}\leq r_{11\cdots im}\leq(b_{i1\cdots im})_{j}for\forall j$
$(b_{i1\cdots im})_{j}-(a_{i1\cdots im})_{j}arrow 0$
as
$jarrow\infty$このとき、 単に
Int$(f)_{j}arrow f$ と書く。
さて、$A$を安定化したアルゴリズムをStab($\mathcal{A}$
) と書くと、 次が安定化理論の基本定理[2] である。
定理
1(
安定化理論の基本定理)
$A$は不連続点$0$の代数的アルゴリズムで、入力$f\in R[x_{1}, \cdots, x_{m}]$ に対し正常終了するとせよ。 このとき、$f$ に対する任意の近似列 $\{$Int$(f)_{j}\}_{j}$ に対し、 ある$n$ が存在して、$j\geq n$
ならば、Stab($\mathcal{A}$) $|$は$\{$Int$(f)_{j}\}_{j}$ に対し正常終了し、
Stab$(\mathcal{A})(Int(f)_{j})arrow \mathcal{A}(f)$
2.2
ISCZ
法
2.2.1
シンボル付き区間 区間と形式的なシンボルを組み合わせた係数 (シンボル付き区間) を導入する。区間は、従来と同様の意 味の区間である。 シンボルは、 アルゴリズム実行中に現れる係数の$\log$ (記録) を取るのに使われる。例え ば、入力係数が 1/3 と 1/7 だったとしよう。これらの精度 3 の区間はそれぞれ[.333, 334] と[.142, 143] で ある。 さて、1/3 に対するシンボルを$s$、 $1/7$に対するシンボルを $t$ として、区間と組み合わせると、それ ぞれ、$[[0.333, 0.334], s]$ と $[[0.142, 0.143], t]$ となる。次に、これらの間の演算、$[[0.333, 0.334], s]+[[O.142, 0.143], t]=[[O.333, 0.334]+[0.142, 0.143], s\dotplus t]$
と定義する。$[$0.333,$0.334]+[0.142$,
0.143
$]$ に対しては通常の区間演算を使う。シンボル部分$s\dotplus t$は再び形 式的なシンボルで、加算を実施したことを記録できれば何でもよい。 効率的なシンボル付けについては 2.3 節で議論する。 アルゴリズム終了後、最終的なシンボルを正確係数に復元する。先の簡単な例で言えば、 もし最終的なシン ボルが$s\dotplus t$であったとすれば、$s$ に1/3を、$t$に1/7を代入し、$\dotplus$ には加算の意味を与えて、 $1/3+1/7=10/21$ と復元する。 シンボル付き区間のことを interval$-symbo1_{Y}$ あるいは単にIS
と呼ぶ。 2.2.2 手続き $A$ を不連続点$0$の代数的アルゴリズムとする。ISCZ 法の手続きは次の通りである。$R-to$
-IS
各入力係数$r$ を $[[A, B], Symbol_{r}]$ に変換する。ここに、$[A, B]$ は$r$の予め定められた精度の区間、$Symbol_{r}$ は$r$ を表すシンボル(以下、入カシンボルと呼ぶ) である。
IS演算
IS
間の演算を次のように実行する:$[[A, B], s]+[[C, D], t]=[[A, B]+[C, D], s\dotplus t]$
$[[A, B], s]-[[C, D], t]=[[A, B]-[C, D], s-t]$
$[[A, B], s]\cross[[C, D], t]=[[A, B]\cross[C, D], s\cross t]$
すなわち、区間部分については区間演算を用い、 シンボル部分については加算、 減算、乗算の形式的
なシンボル十,
$\cross$ を使って、 どういう演算が行なわれたかを記録する。正しいゼロ書換え任意の
IS
$[[A, B], s]$ に対し、$A\leq 0\leq B$ならば、$s$ をそれに対応する実数$r(s)$ に復元する。 もし、 $r(s)=0$ならば、 次のステップに進む。 そうでなければ、 精度を上げて$R-to$
-IS
に戻る。$IS-to-R$ 出力のシンボル部分の中の各入カシンボルにそれぞれ対応する入力係数を代入し、 演算シンポル
に演算の意味を与えて実数値に復元する。
この手法を
ISCZ
法 (IS method withcorrectzero
rewriting) と呼ぶ。さて、 ある精度$j$でInt$(f)_{j}$ を Stab($\mathcal{A}$
) に入力したときの実行過程は、もしアルゴリズム中のすべての ゼロ書換えが正しいならば、真の出力$f$ を$A$に入力したときの実行過程と完全に一致する。ISCZ 法では、 各ゼロ書換えにおいて、それが正しいかどうかを確認する。 さらに、定理 1 により、すべてのゼロ書換えが 正しくなる精度が存在する。従って、ISCZ法は有限ステップで終了し、 その出力の各
IS
係数のシンボル は正しい正確係数を与える。 これを定理にまとめる。定理 2(ISCZ 法の停止性と正当性) $A$ が入力$I$ で正常終了するとせよ。 このとき、$A$ に対する
ISCZ
法は、常に有限ステップで終了し、 正しい結果、 すなわち、$\mathcal{A}(I)$ の出力と同じ結果を与える。
ISCZ
法の利点は、ISZ
法と違い、出力の正当性を確認する必要がないことである。任意のIS
$[[A, B], s]$ に対し、$A\leq 0\leq B$でない限り、正確計算をスキップすることができ、 浮動小数点計算だけで済む。 換言す れば、ゼロでない係数についての正確計算を省略することができる。従って、本手法は、$A\leq 0\leq B$ でな い場合が$A\leq 0\leq B$である場合よりも多ければ多いほど有効であるということができる。
2.3
シンボルリスト
ISCZ
法では、IS
のシンボル部分が膨張してしまうことがある。それを防ぐため、IS
に用いるシンボルは 全て整数とし、計算の経過を行番号に対応させて保存する行列 (シンボルリスト) を用いることとする。シ ンボルリストについては[4] において提案されている。 シンボルリストには$n$行3列の行列を用い、$0$列目1列目には数値もしくはシンボル、2列目には演算 を保存する。 もし、 シンボルが表すものが実数であれば、$0$列目に実数、1 列目と 2 列目に$-1$ を保存し、 演算ではないことを示す。 シンボル番号は、IS
演算を行うごとに新しい番号をつけていく。 例えば $\frac{1}{3}+\frac{1}{7}\cross\frac{2}{9}$ $=$ $[[0.333, 0.334], 0]+[[0.142$, 0.143],1
$]$ $\cross[[0.222$,0.223],2
$]$ $=$ $[[0.333, 0.334], 0]+[[0.031524$,0.031889],3
$]$ $=$ $[[0.364524$,
0.365889],4
$]$ の計算経過をシンボルリストに保存する場合、$0 1 2$
$SL=40231(^{\frac{1}{31}} \frac{}{\frac{72}{091}} -1-1-132 -1-1-1+\cross)$
となる。
3
シンボルリストの評価方法
3.1
行列法
2.3 節で保存したシンボルリストを用いて、シンボル番号4を評価する手順によって説明する。$0 1 2$
$SL=43021(^{\frac{1}{\frac{}{019}\frac{31}{27}}} -1-1-132 -1-1-1+\cross)$流れとしては、 シンボルリストから評価行列を作り、 その評価行列を用いて計算することになる。 各行の$0$列目にはその行にいくつのシンボルリストから取つた3つ組が存在するかを保存する。例えば $0$行目には評価を行うシンボル番号4に対応するシンボルリストの行のみを入れるため、 3つ組の数は1と なり、$0$行$0$列目に 1 が保存される。
$0 1 2 3$
$eva=0(1 0 3 +)$
eva
行列の$n$行$3m+3$列目 $(m\geq 0, m\in \mathbb{Z})$ を読み、 演算記号ならばシンボルリストから $3m$列目 3m$+$1列目のシンボルに対応するシンボルリストの行を$n+1$行目に左詰めで入れる。演算記号ではなく$-1$ ならば
何も行わずスキップする。$n$行目の全ての列が終わったら $n+1$行目以降にも同様の操作を行い、$n$行目に
演算シンボルがなくなるまで繰り返す。
$012 3 456$
$eva=01(\begin{array}{lllllll}1 0 3 + 2 \frac{1}{3} -1 -1 1 2 \cross\end{array})$
$012 3 4 5 6$
$eva=201(\begin{array}{lllllll}1 0 3 + 2 \frac{1}{3} -1 -1 1 2 2 \frac{1}{7} -1 -1 \frac{2}{9} -1 -1\cross\end{array})$
評価行列が完成してから演算を行う。$n$行の評価行列での演算は$n$行目,n–l行目,$\cdots$,0 行目の順に行い、
演算シンボルがあれば 1 行下に保存されている実数を用いて計算を行うこととする。計算結果は演算シン
ボルが属する3つ組に上書きする。
$012 3 4 5 6$
$eva=201(\begin{array}{lllllll}1 0 3 + 2 \frac{1}{3} -1 -l 1 2 2 \frac{1}{7} -1 -1 \frac{2}{9} -1 -1\cross\end{array})$
$\downarrow$
$012 3 4 5 6$
$eva=201(_{2}^{1}2 \frac{1}{\frac{0}{}3,71} -1-13 -1-1+ \frac{2}{\frac{632}{9}}-1-1-- -1-1)$
また、同じシンボルを複数回評価するのは無駄であるため、 一度評価されたシンボルはシンボルリストに
実数として上書きする。
$01 2012$
$SL=43201 ( \frac{1}{021\frac {}{}9731} -1-1-123 -1-1-1+\cross)arrow SL=43201(\begin{array}{lll}\frac{1}{3} -1 -1\frac{1}{7} -1 -1\frac{2}{9} -1 -1\frac{2}{63} -1 -10 3 +\end{array})$
3.2
式列法
式列法の基本形は次のようになる。
$012 3 456 789$
$eva= (a” s =,, s” t \star s" u)$
ここで、$s,$$t,$ $u$ はシンボル番号を、 $a$” は次の数が評価されるシンボルであることを、 $s$” は次の数がシ ンボルであることを、$\star$は演算記号を表す。この行列を$0$列目から走査し、以下の条件に当てはまったら操 作を加える。 $\bullet$ $n$タリ目が $s$ かつ $n+1$列目が表すシンボルが演算のとき $n$列目を起点として基本形に展開する。$n+2$列目以降を右に8つシフトする。
$n n+1 n+2 n+3 n+4 n+5$
.
.
.
$eva=$ $(\cdots$ “$s$” $s$ $\star$ $s$ $t$ $\cdots)$
$\downarrow$
$n n+1 n+2 n+3 n+4 n+5 n+6 n+7 n+S n+9 n+10$
. . .
$eva=$ $(\cdots$ $a$” $s$ $=$ $s$ $u$ $\star$
“s”
$v$ $\star$ $\cdots)$$\bullet$ $n$夕1川が $s$” かつ $n+1$列目が表すシンボルが実数のとき シンボルを表す $s$ ” を実数を表す $r$” に変え、 後ろに実数を入れる。
$n n+1 n+2 n+3 n+4 n+5$
.
.
.
$eva=$ $(\cdots$ “$s$ ”’ $s$ $\star$ $s$” $t$ $\cdots)$ $\downarrow$$n n+1 n+2 n+3 n+4 n+5$
. .
.
$eva=$ $(\cdots$ $r$” $p$ $\star$ $s$” $t$ $\cdots)$
$\bullet$ $n+4$ 列目と $n+7$列目が共に $r$” のとき
演算を行い、$n+10$列目以降を左に8つシフトする。
$n n+1 n+2 n+3 n+4 n+5 n+6 n+7 n+8 n+9 n+10$
. . .
$eva=$
$(\cdots " a" s =,, r p \star r q \star \cdots)$
$\downarrow$
$n n+1 n+2$
.
.
.
$eva=$ $(\cdots$ “$r$” $r$ $\star$ $\cdots)$
また、
演算後は行列法と同様にシンボルリストにも演算結果を上書きする。
この操作を繰り返し、$0$列目が $r$” になった時の 1 列目が評価結果となる。
4
ユークリッドの互除法への適用
表 1: 実験結果(安定化計算有理係数)
4.1
実験方法
$\bullet$ 1変数多項式のユークリッドの互除法を行う。 $\bullet$ 以下の4つの方法で計算開始から出力までの時間メモリ使用量を計測する。ISCZ法については評 価行列の最大要素数も計測する。1.
ISCZ法のユークリッドの互除法アルゴリズムに、 行列法の評価方法を組み込んだもの2.
ISCZ法のユークリッドの互除法アルゴリズムに、 式列法の評価方法を組み込んだもの 3. 自作の正確演算でのユークリッドの互除法アルゴリズム4.
Maple に組み込まれているGCD
関数(正確演算) $\bullet$ 実験環境は以下の通り。 使用ソフト:数式処理ソフト Maple14 使用コンピュータ OS:Windows7
ProfessionalCPU:Intel(R) Pentium(R)
CPU
G840
@2.
$80GHz$メモリ$:4.00GB$
4.2
実験結果
実験結果表中の略号は以下の通り。 入出力次数:(入力した多項式$A$の次数)-(入力した多項式$B$の次数)-(出力された多項式の次数) SP:開始精度から1桁ずつ増やした時の出力精度(開始精度は有理係数では 3 桁、 無理係数では 1 桁) ZR:ゼロ書き換え回数 SL:シンボルリストの長さ Time:ミリ秒 Memory:k$=kibM=Mib$ Elements:評価行列の要素数 表1と表2は入力多項式の係数が有理数だけのものである。 有理係数での実験では行列法と式列法につ いてはどちらが優位とは決められなかった。 また、 自作の正確演算アルゴリズム、Mapleに組み込まれているGCD
関数と比べると時間メモリ使用 量共に ISCZ 法を用いない方が優位であった。表 2: 実験結果
(
厳密計算有理係数)
表 3: 実験結果(安定化計算無理係数)