数理解析研究所講究録
第 1815 巻 2012 年 133-142

133

安定化理論に基づく ISCZ 法の凸包構成への応用
白柳潔*

関川浩\dagger

東邦大学理学部東海大学理学部

1

KIYOSHI SHIRAYANAGIHIROSHI SEKIGAWA

FACULTY OF SCIENCE, TOHO UNIVERSITYSCHOOL OF SCIENCE, TOKAI UNIVERSITY

はじめに安定化理論 [9] は、 近似計算で実行すると不安定となるアルゴリズムに対し、 それを変形して近似計算で実行しても誤差の影響を抑制し、 安定な出カが得られるようにするための理論である。 [6, 7] では、 その安定化理論の新しい応用として、 なるべく正確計算を軽減させながらも最終的に正確係数の正しい出カを得るための手法 ISCZ 法を提案し、 Buchberger アルゴリズムに適用して得られた実験結果を報告した。 そこでは、有理係数よりも無理係数の方が有効であることはわかったが、 もともとすべてを正確演算で実行した結果と比べて、 必ずしも優位性を主張することができなかった。 [8] では、 その原因を追求し、 Buchbergerアルゴリズムのように、 再帰的な呼び出しが頻発するようなアルゴリズムでは、 ISCZ 法はあまり有効でないことを結論づけた。 さらに、 凸包アルゴリズムのように、 平坦な構造をもったアルゴリズムには有効であることを予想した。本論文では、 実際に ISCZ 法を凸包アルゴリズムに適用し、 無理数を含む場合は格段に有効になることを示す。
まず、 2 節で、 安定化理論と ISCZ 法について復習する。 3 節で、 ISCZ 法を凸包アルゴリズムのーつである Graham のアルゴリズムに適用した計算機実験について導べる。

2

復習

2.1

安定化理論.
.

次のアルゴリズムを対象に安定化理論の復習を簡単に行う。 詳細は [5, 9] を参照されたい。

$\circ$

データは、 すべて多項式環 $R[X_{1}, \ldots , x_{m}]$ の元からなる。 $R$ は実数体の部分体である。

データ間の演算は、 $R[x_{1}, \ldots, x_{m}]$ 内の加減乗である。
データ上の述語は、 不連続点をもつとすればそれは$0$のみである。という意味は、 If $C=0$ ”thenのように、 値が か否かにょって分岐がelse別れることである。 従って、 $C=0$ の代わりに、 $C>0$ や $C\geq 0$ などでもよい。 上記クラスのアルゴリズ

述語の不連続点がムを、 不連続点

$0$$\cdots$$\cdots$$0$の代数的アルゴリズムと呼ぶ。 ほとんどの数式処理のアルゴリズムはこのクラスに入るか、 このクラスのアルゴリズムに変換可能である。
$0$‘lciyoshi.shirayanagi@is.sci.toho-u.ac.jp\dagger sekigawa@tokai-u.jp

134.
.
.

さて、 安定化の 3 つのポイントは、アルゴリズムの構造は変えない。
データ領域において、 ふつうの係数を区間係数に変える。

述語の評価の直前で、 区間係数のゼロ書換えを行なう。

である。 すなわち、 安定化されたアルゴリズムは次のようになる。

区間領域データ領域は区間係数多項式の集合。区間係数は$[A, \alpha]$なる形で、 $A\in R,$$\alpha$は非負の実数。

$[A, \alpha]$

は集合 $\{x\in R||x-A|\leq\alpha\}$ を意味する。

区間演算二項演算 $*\in\{+, -, \cross, /\}$ に対し、$[A, \alpha]*[B, \beta]=[A*B,\gamma_{*}].$

ここに、$\gamma_{*}$

は次を満たす。
$|x-A|\leq\alpha, |y-B|\leq\beta\Rightarrow|x*y-A*B|\leq\gamma_{*}.$

ゼロ書換え不連続点$0$をもつ述語を評価する直前で、 各区間係数$|C|\leq\gamma$

ならば$|C|>\gamma$$[C, \gamma]$$[C, \gamma]$

に対し、を [0,0] に書き換えよ。

ならばそのままとせよ。

区間演算の具体的な定義については、文献 [1] を参照されたい。
今、 入力 $f\in R[x_{1}, \ldots, x_{m}]$ を$f= \sum_{i_{1},\ldots,i_{m}}a_{i_{1}\ldots i_{m}}x_{1}^{i_{1}}\cdots x_{m}^{i_{m}}$と表したとき、$f$に対する近似列 Int$\{$Intで定義する。 ここに、 すべての$(f)_{j}\}_{j}$を...

$x_{m}^{i_{m}}$$(f)_{j}= \sum_{i_{1},\ldots,i_{m}}[(a_{i_{1}\ldots i_{m}})_{j}, (\alpha_{i_{1}\ldots i_{m}})_{j}]x_{1}^{i_{1}}$$i_{i},$$\ldots,$$i_{m}$

について、

for

$|a_{i_{1}\ldots i_{m}}-(a_{i_{1}\ldots i_{m}})_{j}|\leq(\alpha_{i_{1}\ldots i_{m}})_{j}$$(\alpha_{i_{1}\ldots i_{m}})_{j}arrow 0$

as$\forall j,$$jarrow\infty.$

このとき、 単に

Int $(f)_{j}arrow f$と書く。
さて、$\mathcal{A}$を安定化したアルゴリズムを Stab定理 1(安定化理論の基本定理)$\mathcal{A}$$(\mathcal{A})$司に対しは不連続点 の代数的アルゴリズムで、 入力 $f\in R[x_{1},$$i\geq n$が存在して、あるに対し、Intに対する任意の近似列$0$正常終了するとせよ。 このとき、に対し正常終了し、は Int
ならば、 Stab$f$$(\mathcal{A})$と書くと、 次が安定化理論の基本定理である。
$\ldots,$$\{$$(f)_{j}\}_{j}$$x$$n$$(f)_{j}$Stab$(\mathcal{A})(Int(f)_{j})arrow \mathcal{A}(f)$.

簡明を期すため、 入力は一つだけの多項式にしているが、 入力はもちろん、 多項式の有限集合でもよい。

主題のグレブナ基底の場合も、 入力は多項式の集合である。本定理の Buchberger アルゴリズムに特化した証明は [4] に、より一般的な場合に対する厳密な証明は [9] にある。

1352.2ISCZ 法

2.2.1

シンボル付き区間

区間と形式的なシンボルを組み合わせた係数 (シンボル付き区間) を導入する。 区間は、 従来と同様の意味の区間である。 シンボルは、 アルゴリズム実行中に現れる係数の(記録) を取るのに使われる。 例
えば、 入力係数が 1/3 と 1/9 だったとしよう。 これらの精度 3 の (円形) 区間はそれぞれ [.333, .0005] と$\log$[.111, .0005] である。 さて、 1/3 に対するシンボルをせると、 それぞれ、 [.333, .0005], $s]$ と [.111, .0005],$[$$[$$[$$s$$t]$、$1/9$に対するシンボルを$t$として、 区間と組み合わとなる。次に、 これらの間の演算、 例えば、 加算を、[.333, .0005], $s]+[[.111, .0005], t]=[[.333, .0005]+[.111, .0005], \dotplus(s, t)]$と定義する。 $[.333, .0005]+$ [.111, .0005] に対しては通常の区間演算を使う。 シンボル部分$\dotplus(s, t)$は再び形式的なシンボルで、 加算を実施したことを記録できれば何でもよい。

アルゴリズムの実行中あるいは実行後に、 必要に応じて、 シンボルを正確係数 (実数値) に復元することができる。先の簡単な例で言えば、 もしシンボルがであったとすれば、 に 1/3 を、 に 1/9 を代入$\dotplus(s, t)$し、$\dotplus$には加算の意味を与えて、 $1/3+1/9=4/9$ と復元する。

シンボル付き区間のことを interval

2.2.2

$t$

$s$

$-symbo1_{\backslash }$あるいは単に $IS$ と呼ぶ。

手続き

従来の ISZ 法 [10, 11, 6] は、 述語を評価する直前で、 任意の $IS[[E, \epsilon], s]$ に対し、であれ、$[[E, \epsilon], s]$を$0$に置き換える手法であった。 それに対して、 ISCZ 法は、を実数値に復元して、 実際にそれが$\mathcal{A}$を不連続点

$0$

$0$

$|E|\leq\epsilon$$|E|\leq\epsilon$

ならば、 が何$s$ならばその都度、$s$であるかどうかを確かめるものである。 詳細は以下の通りである。の代数的アルゴリズムとする。

[ISCZ 法]-to- $IS$ 各入力係数$R$$Symbol_{a}$$IS$は$a$$a$をペア $[[\tilde{a}, \alpha], Symbol_{a}]$ に変換する。 ここに、$[\tilde{a}, \alpha]$

は$a$の指定された精度の区間、を表すシンボル (以下、 入カシンボルと呼ぶ) である。

演算 $IS$ 間の演算を次のように実行する :$[[A, \alpha], s]+[[B, \beta], t]=[[A, \alpha]+[B, \beta], \dotplus(s, t)]$$[[A, \alpha], s]-[[B, \beta], t]=[[A, \alpha]-[B, \beta], -(s, t)]$$[[A, \alpha], s]\cross[[B, \beta], t]=[[A, \alpha]\cross[B, \beta], \cross(s, t)]$すなわち、 区間部分については区間演算を用い、 シンボル部分については加算、 減算、 乗算の形式的なシンボル十，$-,$$\cross$を使って、 どういう演算が行なわれたかを記録する。

正しいゼロ書換え 述語を評価する直前で、 任意の $IS[[E, \epsilon], s]$ に対し、ならば、 をそれに対応する実数値 $r(s)$ に復元する。 もし、 $r(s)=0$ ならば、 次のステップに進む。 そうでなければ、 精度を上$|E|\leq\epsilon$$s$

げて $R$-to- $IS$ に戻る。
$IS$-to- $R$ 出力のシンボル部分の中の各入カシンボルにそれぞれ対応する入力係数を代入し、 演算シンボルに演算の意味を与えて実数値に復元する。

136
この手法を ISCZ 法 ( $IS$ method with correct zero rewriting) と呼ぶ。に入力したときの実行過程は、 もしアルゴリズム中のすべてのを Stab
さて、 ある精度 で Int
ゼロ書換えが正しいならば、 真の出力 $f$ を に入力したときの実行過程と完全に一致する。 ISCZ 法では、$(f)_{j}$$i$$(\mathcal{A})$$\mathcal{A}$

各ゼロ書換えにおいて、それが正しいかどうかを確認する。さらに、定理 1 により、すべてのゼロ書換えが正しくなる精度が存在する。従って、 ISCZ 法は有限ステツプで終了し、 その出力の各 $IS$ 係数のシンボルは正しい正確係数を与える。
これを定理にまとめる。
に対する ISCZ 法が入力 で正常終了するとせよ。 このとき、
定理 2(ISCZ 法の停止性と正当性)の出力と同じ結果を与える。
は、 常に有限ステツプで終了し、 正しい結果、 すなわち、$\mathcal{A}$$\mathcal{A}$$I$$\mathcal{A}(I)$

ISCZ 法の利点は、 ISZ 法と違い、 出力の正当性を確認する必要がないことである。任意の$IS[[E, \epsilon], s]$にでない限り、 正確計算をスキツプすることができ、 浮動小数点計算だけで済む。 換言すれば、対し、ゼロでない係数についての正確計算を省略することができる。従って、 本手法は、 $|E|>\epsilon$ の場合が$|E|\leq\epsilon$$|E|\leq\epsilon$の場合よりも多ければ多いほど有効であるということができる。

3

凸包アルゴリズムへの応用[8] では、 ISCZ 法が Buchberger のアルゴリズムに対してあまり有効でないのは、 簡約など再帰的な呼び出しが頻繁に行なわれるため、シンボル係数のシンボルの膨張が激しいからであることを究明した。シンボルの膨張を防ぐために導入されたシンボルリストも功を奏さなかった。 (シンボルリストの詳細については、 [8] を参照されたい)。 そこで、 本論文では、 再帰的な呼び出しのない、 構造が平坦な凸包アルゴリズムに適用する。凸包アルゴリズムは複数存在するが、 ここでは簡明な Graham のアルゴリズム [2, 3] を選定した。 そのアルゴリズムを記述する。

Graham のアルゴリズム
Input: 有限集合
Output: の凸包の頂点集合$S=\{P_{1}, P_{2}, \ldots, P_{m}\}\subset \mathbb{R}^{2}$$S$

1. [基準点]$S$の点を一つ選び、 $O$ とする。の点を、 $0$ からの偏角が小さい順にソートする。 ただし、 偏角が同じ点が複数ある場合は、とし、 $Q_{0}=O$ とおく。
$O$ から最も遠い点のみをとり、 その他は捨てる。 その結果を

2. [ソート]$S$$Q_{1},$$Q_{2},$$\ldots,$$Q_{n}$

3. [走査] 内角 $Q_{i}Q_{i-1}Q_{i-2}$ が $180 °$ 未満かどうか調べながら凸包を構成する。 厳密に書くと次のようになる :$j:=0$

forfrom to$j:=j+1$
$i$

$0$

$n$

do

$R_{j}:=Q_{i}$while $j\geq 4$ and $\angle R_{j}R_{j-1}R_{j-2}\geq 180^{o}$ do$R_{j-1}:=R_{j}$$j:=j-1$
return$[R_{1}, R_{2}, \ldots, R_{j}]$137
ここで、 一般に 3 点 $P_{1}=(x_{1}, y_{1}),$ $P_{2}=(x_{2}, y_{2}),$ $P_{3}=(x_{3}, y_{3})$ に対し、$\angle P_{3}P_{2}P_{1}\geq 180^{o}\Leftrightarrow$が成立するので、 Graham のアルゴリズムは不連続点$0$$x_{1}$$y_{1}$$x_{2}$$y_{2}$$x_{3}$$y_{3}$

111$\leq 0$の代数的アルゴリズムである。従って、 安定化理論の適用条件を満たす。
我々は、 このアルゴリズムを適用するための入力例を 4 つの群に分けた。 第一群は 1000 個の整数点、 第二群は 1000 個の無理数点、第三群は 5000 個の整数点、 第四群は 5000 個の無理数点からなる。各群に対し、4 つの例を設けた。
第一群は、 次の 4 例である。
例 $10\leq X,$ $Y\leq 100$ を満たすランダムに生成された整数の組 $(X, Y)1000$ 点。

例 2 次を満たすランダムに生成された整数の組 $(X, Y)$ 1000 点。
$-100\leq X, Y\leq 100, X^{2}+Y^{2}\leq 100^{2}.$例 3 次を満たすランダムに生成された整数の組 $(X, Y)$ 1000 点。
$0\leq X, Y\leq 100, X^{2}+Y^{2}\leq 100^{2}, Y^{2}\leq 3X^{2}.$すなわち、 中心が原点、 半径 $100$、角度 $60 °$ の扇形の内部か境界線にある 1000 点である。

例 4 原点と、 次を満たすランダムに生成された整数の組 $(X, Y)999$ 点。
$1 \leq X, Y\leq 3000, X^{2}+Y^{2}\leq 3000^{2}, \frac{9}{10}\leq\frac{Y}{X}\leq1.$

すなわち、 中心が原点、 半径 3000 の長細い扇形の内部か境界線にある 1000 点である。

第二群は、 次の 4 例である。
例 5 $0\leq X,$ $Y\leq 100$ を満たすランダムに生成された整数の組 $(X, Y)1000$ 点に対し、 $(\sqrt{X}, \sqrt{Y})$ を作ったもの。

例 6 次を満たすランダムに生成された整数の組 $(X, Y)$ 1000 点に対し、 (sign $(X)\sqrt{|X|}$, sign $(Y)\sqrt{|Y|}$ ) を作ったもの。 ただし、 sign(X) は $X$ の符号である。
$-100\leq X, Y\leq 100, X+Y\leq 100.$

例 7 次を満たすランダムに生成された整数の組 $(X, Y)$ 1000 点に対し、$(\sqrt{X}, \sqrt{Y})$を作ったもの。

$0\leq X, Y\leq 100, X+Y\leq 100, Y\leq 3X.$

すなわち、 中心が原点、 半径 $10$、角度$60^{O}$の扇形の内部か境界線にある 1000 点である。

例 8 原点と、 次を満たすランダムに生成された整数の組 $(X, Y)999$ 点に対し、$(\sqrt{X}, \sqrt{Y})$$1 \leq X, Y\leq 3000, X+Y\leq 3000, \frac{9}{10}\leq XY$$-\leq 1.$

すなわち、 中心が原点、 半径 $10\sqrt{30}$ の長細い扇形の内部か境界線にある 1000 点である。

を作ったもの。

138
第三群は、 次の 4 例である。
例 9 $0\leq X,$ $Y\leq 200$ を満たすランダムに生成された整数の組 $(X, Y)$ 5000 点。

例 10 次を満たすランダムに生成された整数の組 $(X, Y)$ 5000 点。
$-200\leq X, Y\leq 200, X^{2}+Y^{2}\leq 200^{2}.$例 11 次を満たすランダムに生成された整数の組 $(X, Y)$ 5000 点。
$0\leq X, Y\leq 400, X^{2}+Y^{2}\leq 400^{2}, Y^{2}\leq 3X^{2}.$すなわち、 中心が原点、 半径 $400$、角度 $60 °$ の扇形の内部か境界線にある 5000 点である。

例 12 原点と、 次を満たすランダムに生成された整数の組 $(X, Y)$ 4999 点。
$1 \leq X, Y\leq 6000, X^{2}+Y^{2}\leq 6000^{2}, \frac{9}{10}\leq\frac{Y}{X}\leq1.$

すなわち、 中心が原点、 半径 6000 の長細い扇形の内部か境界線にある 5000 点である。

第四群は、 次の 4 例である。
例 13 $0\leq X,$ $Y\leq 200$ を満たすランダムに生成された整数の組 $(X, Y)$ 5000 点に対し、$(\sqrt{X}, \sqrt{Y})$を作ったもの。

例 14 次を満たすランダムに生成された整数の組 $(X, Y)$ 5000 点に対し、 (sign $(X)\sqrt{|X|}$ , sign $(Y)\sqrt{|Y|}$) を作ったもの。 ただし、 sign(X) は $X$ の符号である。
$-200\leq X, Y\leq 200, X+Y\leq 200.$

例 15 次を満たすランダムに生成された整数の組 $(X, Y)$ 5000 点に対し、$(\sqrt{X}, \sqrt{Y})$を作ったもの。

$0\leq X, Y\leq 400, X+Y\leq 400, Y\leq 3X.$

例 16 原点と、次を満たすランダムに生成された整数の組 $(X, Y)4999$ 点に対し、$(\sqrt{X}, \sqrt{Y})$を作ったもの。

$1 \leq X, Y\leq 6000, X+Y\leq 6000, \frac{9}{10}\leq\frac{Y}{X}\leq 1.$

すなわち、 中心が原点、 半径 $20\sqrt{15}$ の長細い扇形の内部か境界線にある 5000 点である。

例 4 と例 12 の入力を図 1 に示す。
計算機は、 dual core AMD(R) Opteron(R) processor $(2.85GHz),$ $8GB$ RAM, Linux(R) $OS$ である。
我々は、 次の 2 種類の方法を試みた。

1. Maple 12 で実装した ISCZ 法。 初期精度は 1 で、 精度の増加幅も 1。

2. Maple 12 で実装した ISCZ 法。初期精度は 10 で、 精度は倍々で増加。

139

図 1: 例 4(左) と例 12(右).

表 1: 方法 1–1000 点 (整数)

表 2: 方法 1–1000 点 (無理数)

表 3: 方法 1–5000 点 (整数)$\ovalbox{\tt\small REJECT}_{R_{-}CHJ}ffi\mathfrak{l}\ovalbox{\tt\smallREJECT} j\iota B$セ
きえ $ffl$スキッ$\mathfrak{o}$表 4: 方法 1–5000 点 (無理数)$\ovalbox{\tt\small REJECT}_{R_{-}CH}fflJ\ovalbox{\tt\smallREJECT}\# y\iota*\mathfrak{o}$
きえスキッ$\# j$

140

表 5: 方法 2–1000 点 (整数)

表 6: 方法 2–1000 点 (無理数)

表 7: 方法 2–5000 点 (整数)

表 8: 方法 2–5000 点 (無理数)141
方法 1 と方法 2 による実験結果をそれぞれ、 表 1-4 と表 5-8 に示す。各表において、 時間は、 初期精度から成功精度に至るまでにかかった cpu 時間 (単位は秒) の合計を表す。 成功精度とは、 ISCZ 法が終了した (成功に至った) ときの入力の浮動小数点精度である。 ゼロ書き換えは、 成功精度において、 ゼロ書換えによって実際にゼロに書き換えられた区間の個数を表す。 スキップ数は、 成功精度において、 “ $|E|>\epsilon$ ” であるゆえに正確計算がスキップされたゼロでない係数の個数を表す。 R-$CH$ は、 比較のために Maple12 で自前で Graham のアルゴリズムを実装したもので、 全過程を通して正確演算で計算したときの cpu 時間である。

これらの表により、 次のことが結論づけられる。

ISCZ 法は、 Graham のアルゴリズムに対し、 入力が有理数点の場合は、 正確計算の直接法よりも 1.5 倍から 2 倍程度遅くなるが、 無理数点の場合は、 5 倍から 8 倍程度高速となる。 従って、正確計算回避の効果が十分に現れているといえる。
さらに、 凸包構成の場合、 得られた結果が真に正しい凸包であるかを、 正確計算で求めることなしに検証することは一般に難しいが、 ISCZ 法を使えば、 後の検証の必要もなく結果が真に正しいことを保証してくれる。

4

おわりにGraham のアルゴリズムに対する ISCZ 法は、 入力が無理数点の場合に大変有効であることが実証された。今後は、 他のアルゴリズムにも適用し、 ISCZ 法が有効となるアルゴリズムを特定していく予定である。

謝辞
本研究は科研費 21500026 の助成を受けたものである。

[1] G. Alefeld and J. Herzberger. Introduction to Interval Computations, Computer Science and Applied
Mathematics. Academic Press, 1983.
[2] R. L. Graham. An efficient algorithm for determining the convex hull of a finite planar set. Information Processing Letters, $1(4):132-133$ , 1972.
[3] F. P. Preparata and M. I. Shamos. Computational Geometry: An Introduction. Springer-Verlag,1985.
[4] K. Shirayanagi. Floating point $Gr6bner$ bases. Mathematics and Computers in Simulation, $42(4-$$6):509-528$ , 1996.
[5] 白柳潔．アルゴリズムの安定化理論．数式処理，$5(2):2-21$ ,1997.

[6] 白柳潔，関川浩．安定化理論に基づく log method について．京都大学数理解析研究所講究録掲載予定[7] K. Shirayanagi and H. Sekigawa. Reducing Exact Computations to Obtain Exact Results Based onStabilization Techniques. In Proc. International Workshop on Symbolic-Numenc Computation 2009
(SNC2009), pages 191-197, 2009.

142[S] 白柳潔，関川浩．安定化理論に基づく ISCZ 法の有効性について．京都大学数理解析研究所講究録掲載予定
theory of stabilizing algebraic algorithms. Technical Report[9] K. Shirayanagi and M. Sweedler.
95-28, Mathematical Sciences Institute, Comell University, 1995. 92 pages.
. ac.jp/sci/is/shirayanagi/msitr95-28. pdf)
(http: $//www$ . lab.
$A$$toh\triangleright u$[10] K. Shirayanagi and M. Sweedler. Automatic algorithm stabilization. In ISSAC’96 poster session
abstracts, pages 75-78, 1996.
[11] K. Shirayanagi and M. Sweedler. Remarks on automatic algorithm stabilization. J. Symbolic Computation, $26(6):761-766$ , 1998.