数理解析研究所講究録
第 1815 巻 2012 年 143-151

143

符号情報利用による数値数式 CAD の効率的な実装
岩根秀直穴井宏和(株) 富士通研究所*(株) 富士通研究所/九州大学\dagger
HIDENAO IWANEHIROKAZU ANAI

FUJITSU LABORATORIES LTD

FUJITSU LABORATORIES LTD /KYUSHU UNIVERSITY

屋並仁史(株) 富士通研究所\ddagger
HITOSHI YANAMI
FUJITSU LABORATORIES LTD

1

はじめに

限定記号消去アルゴリズム (quantifier elimination ( $QE$ ) algorithm) とは，与えられた形式理論 (formaltheory) について 「限定記号付きの式 (一階述語論理式)」 を入力とし 「等価で限定記号無しの式」 を出力$QE$ アルゴリズムは最適化問題等の様々な適用分野があり非常に有効であするアルゴリズムのことである．
Algebraic
Decomposition( $CAD$ ) [2] は，与えられた多項式集合に対して，符号が不変でる．Cylindrical
ある領域に分割する手法で，限定記号消去問題の基本的な解法のひとつである．
我々は 2009 年に Dynamic Evaluation ( $DE$ ) と数値数式手法の組み合わせによる効率的な $CAD$ の実装$DE$ は定義多項式の既約性を必要とせず，無を提案した [6]. $DE$ [4,5] は拡大体の表現方法の一つである．

平方性のみを要求する．そのため代数拡大体上での因数分解を回避でき計算の効率化につなげることがで$CAD$ 計算で代数拡大体上での実根の分離が必要となるが，代数的数をそれをきた．数値手法については，含む区間で置き換え，区間係数一変数多項式の実根の分離として扱うようにした．すべての数値演算で区間演算を利用することで正確性を失わずに効率化を実現することができた．
数値手法を利用することで多くの記号計算を回避することができたが，区間が を含む場合には記号計$DE$ が無平方性を要求するために無平方因子分解の処理が回避できない．
算による判定が必要であるし，本稿では $CAD$ が再帰的に射影因子の符号不変な領域を構築することを利用し，数値数式 $CAD$ における記号計算をさらに回避する実装方法について述べる．
$O$

2

$CAD$

アルゴリズム1975 年に，G. E. Collins が，与えられた多項式集合に対して変数空間を各多項式の符号が不変であるセル (cell) とよばれる領域に分割する新しい代数的方法である Cylindrical Algebraic Decomposition( $CAD$ )[2] を提案した．
*iwane@jp.fujitsu.com\dagger anai@jp.fujitsu.com

tyanami@labs.fujitsu.com

144
$CAD$

アルゴリズムは射影段階，底段階，持ち上げ段階の 3 つの段階から成る．射影段階では，入力の多項式集合に射影演算子 PROJ を繰り返し適用することで 1 変数ずつ消去していく．入力の多項式集合をとすると，射影段階により多項式列君 $=F,$ $F_{r-1}=PROJ(F_{r},x_{r})$ ,を既約因子分解して得られる多項式集
PROJを得る．

$F=\{f_{i}(x_{1}, \ldots, x_{r})\}\subseteq \mathbb{Q}[x_{1}, \ldots,x_{r}]$$F_{r-2}=$

PROJ $(F_{r-1}, x_{r-1})$ , . . ., $F_{1}=$$(F_{2}, x_{2})$$\bigcup_{i=1}^{r}F_{i}$を射影因子という．ここで，$P_{k}=\{p\in P|\deg_{x_{k}}(p)>0, \deg_{x}:(p)=0(i>k)\}$ とする．底段階では，の分解を行う．これは射影段階により得られた 1 変数多項式の集合 の実根合

$P\subseteq \mathbb{Q}[x_{1}, \ldots, x_{r}]$$P_{1}$$\mathbb{R}(=\mathbb{R}^{1})$を求めることで$D_{k+1}$を分割する．最後の持ち上げ段階では，$\mathbb{R}^{k}$$\mathbb{R}^{1}$の分解$D_{k}$と$P_{k}$を用いて$\mathbb{R}^{k+1}$の分解を得る $(k=1, \ldots, r-1)$ .

により得られるセルでは多項式集合の符号が一定なので，標本点 (sample point) と呼ばれる任意の一点の符号を評価することで入力の多項式集合の符号を評価することができる．標本点の各座標は代数的数になり，その定義多項式とその根を唯 1 つ含む分離区間 (isolating interval) を用いて表現し，分離区$CAD$間を求めることを実根の分離と言う．

以下に持ち上げ段階のアルゴリズムを示す．ここでは拡大体上での実根の分離処理をそれぞれの射影因子に対して行い，異なる実根を表す分離区間が重複しないようにする必要があり，分離区間が重複した場合の共通根判定処理を行う．判定処理により異なることが確認できた場合には，標本点 を定める分離区間$S$の区間幅を小さくすることで，分離区間が重複しないようにする．
アルゴリズム 1

LIFT $(P_{k}, S)$
入力$:P_{k}\subseteq \mathbb{Q}[x_{1}, \ldots, x_{k}]$セル $C_{k-1}$ の標本点 $S=(s_{1}, \ldots, s_{k-1})\in \mathbb{R}^{k-1}$

出力: $C_{k-1}$ を持ち上げた結果のセルの標本点

for$f_{i}$

in$L_{i}arrow$for$L_{i},$if$P_{k}$

doREALROOTISO
do$(f_{i}, S)$(分離区間のリスト)$L_{j}$の分離区間に重複がある thenif 同じ代数的数を表している then
一方の分離区間を削除する$L_{i},$$L_{j}$

else
区間の重複がなくなるまで精度を上げる$CAD$ の計算において$P_{k}$

は $P_{k+1}$$D_{k+1}$を構築する場合には疏に対して符号一定な $CADD_{k}$ を用いる．このときに対する射影演算により生成されたもので，$P_{k}$の符号情報を利用することで $P_{k+1}$ に関する情報を得ることができる．本稿ではこれを利用して記号計算を回避することを考える．

2.1

数値数式 $CAD$

本稿は数値数式 $CAD$ [6] の実装に対する効率化であるため，その実装方法について簡単に述べる．
数値数式 $CAD$ では持ち上げ段階における実根の分離処理を数値演算に置き換えることで効率化を実現した．以下に数値数式手法による実根の分離処理を示す．
アルゴリズム 2

REALROOTISO $(f, S)$入力 $:f\in \mathbb{Q}[x_{1}, \ldots, x_{k-1}, x_{k}](k\geq 2)$

145
セル $C_{k-1}$ の標本点 $S=(s_{1}, \ldots, s_{k-1})\in \mathbb{R}^{k-1}$

出力: $f(s_{1}, \ldots, s_{k-1}, x_{k})$ の分離区間$garrow$$Larrow$

SUBSTSAMPLEPOINT $(f, S)$
INTVREALROOTISO $(g)$if error then$farrow$ SYMBOLICSQRFREE $(f, S)$$do$$garrow$

SUBSTSAMPLEPOINT $(f, S)$$Larrow INTvREALRooTIso(g)$if error then$S$の精度をあげる

elsebreakreturn$L$

ここで，SUBSTSAMPLEPOINT は に標本点を表す分離区間を代入し，区間係数の一変数多項式を返す関数である．INTVREALROOTISO は区間係数一変数多項式の実根の分離処理で，Krawczyk 法 [7] を基本とした実装としている．Krawczyk 法は無平方でない場合には正しく実根を分離できないためエラーを返$DE$ が要求する定義多項式の無平方性を満たす必要がある．
す．我々の実装では $DE$ を利用しているため，
そこで Krawczyk 法での実根の分離処理が失敗した場合にのみ記号計算による代数拡大体上での無平方因子分解 SYMBOLICSQRFREE を呼び出すようにすることで不要な演算を回避することを実現した．しかし，INTVREALROOTISO は係数の区間が十分に小さくない場合にも実根の分離ができないためにエラーを返すのでその場合には不要な記号演算を呼び出すことになる．ここで複素根の重根については除去できないが$f$$CAD$ の実装では実根のみを使用するため影響がない．

SUBSTSAMPLEPOINT では分離区間の精度が十分出ない場合には係数が膨張してしまうことがあるため，復帰する多項式の係数は符号が確定しているようにした．っまり，係数が本当に 0 となる場合以外には係数の区間が を含まないようにした．以下に代入処理のアルゴリズムを示す．
$O$

アルゴリズム 3

SUBSTSAMPLEPOINT $(f, S)$
入力 $:f\in \mathbb{Q}[x_{1}, \ldots, x_{k-1}, x_{k}](k\geq 2)$
セル $C_{k-1}$ の標本点 $S=(s_{1}, \ldots, s_{k-1})\in \mathbb{R}^{k-1}$

出力: $f$ に標本点の分離区間を代入して得られる区間係数一変数多項式

for $i=0$ to degree of$c_{i}arrow coefficient$ of

$f$

inin

$x_{k}^{i}$$x_{k}$$f$$I_{i}arrow c_{i}(s_{1}, \ldots, s_{k-1})$if$I_{i}$

が$0$を含む thenif SYMBOLICZEROCHECK$(c_{i}, S)$

then$I_{i}arrow[0,0]$else$I_{i}$

returnが$\sum_{i}I_{i}x^{i}$$0$を含まなくなるまで $S$ の精度をあげる

146
標本点の代入した結果が でない場合にも分離区間を代入する際にその精度が十分でない場合に入力の多項式 の係数が大きい場合等に を含んでしまうため，記号演算によりチェックが必要となることがあった．
$O$

$0$

$f$

これまでの話をまとめると，数値数式 $CAD$ の持ち上げ段階において記号計算が発生する箇所は以下の3 箇所である．

1. 係数に標本点を代入した際に発生する零判定 (区間が$O$を含む場合)2. $DE$ が要求する無平方因子分解 (実根の分離に失敗した場合)
3. 異なる射影因子から生成される代数的数の共通根判定 (区間が重複した場合)
次節では $CAD$ の計算結果を利用することによりこれらの記号計算を回避する方法について述べる．

符号情報を利用した数値数式 $CAD$ の実装

3の計算量に大きく影響するため，これまで多くの研究が射影因子を減らすことに注力されてきた．
本稿で提案する手法は使用する射影演算に依存する．ここでは McCallum 射影演算子 [9] を例にして紹介$CAD$する．

アルゴリズム 4

PROJMC$(P_{k}, x_{k})$入力 $:P_{k}\subseteq \mathbb{Q}[x,., xx]$$(k\geq 2)$

主変数$x_{k}$

出力: 射影因子$P_{k-1}\subseteq \mathbb{Q}[x_{1}, \ldots, x_{k-1}]$$Garrow\emptyset$for

$f$

in

$P_{k}$

do$Garrow G\cup$$Garrow G\cup$for$f,$ $g$

in$f$の$x_{k}$$f$の$x_{k}$do$P_{k}$$Garrow G\cup$

return{{{$f$に関するすべての係数}を主変数とみたときの判$B^{1}J$式}$(f\neq g)$と$g$の$x_{k}$を主変数とみたときの終結式}$G$

PROJMC により得られる射影因子には係数，判別式，終結式の情報が含まれている．この符号情報を利用することで多くの記号計算の回避が実現できる．

3.1

標本点の分離区間の代入処理

ここではアルゴリズム 3 における符号情報の利用について述べる．PROJMC では入力の多項式のすべての係数の情報を射影因子として出力する．つまり，持ち上げ段階において多項式 $f$ に標本点 $S$ を代入した結果の符号は一定で，その値は既知である．従ってアルゴリズム 3 において発生する記号演算を完全に回避することができる．

アルゴリズム 5

NEWSUBSTSAMPLEPOINT $(f, S)$
入力 $:f\in \mathbb{Q}[x_{1}, \ldots, x_{k-1}, x_{k}]$$(k\geq 2)$

セル $C_{k-1}$ の標本点 $S=(s_{1}, \ldots, s_{k-1})\in \mathbb{R}^{k-1}$

147
出力: $f$ に標本点の分離区間を代入して得られる区間係数一変数多項式

for $i=0$ to degree ofcoefficient of$c_{i}arrow$$I_{i}arrow$if$x_{k}$$f$$x_{k}^{i}$砺 $(s_{1}, \ldots, s_{k-1})$

が$I_{i}$

inin$f$を含む then$0$

if 射影因子の符号から$c_{i}(s_{1}, \ldots, s_{k-1})$が$0$である then$I_{i}arrow[0,0]$else
が$I_{i}$return3.2$0$を含まなくなるまで $S$ の精度をあげる$\sum_{i}I_{i}x^{i}$

代数拡大体上での無平方因子分解

数値数式 $CAD$ における実根の分離処理の部分では符号情報を利用することで不要な無平方因子分解を回避することができる．
アルゴリズム 6

NEWREALROOTISO $(f, S)$
入力 $:f\in \mathbb{Q}[x_{1}, \ldots, x_{k-1}, x_{k}]$セル$C_{k-1}$の標本点$(k\geq 2)$$S=(s_{1}, \ldots, s_{k-1})\in \mathbb{R}^{k-1}$

出力: $f(s_{1}, \ldots, s_{k-1}, x_{k})$ の分離区間$garrow SuBSTSAMPLEPoINT(f, S)$$Larrow INTvREALRooTIso(g)$if error then$c_{i}arrow le$

if

$f$

ding coefficient of $f$ inの判別式が$0$$x_{k}$

then$farrow$

SYMBOLICSQRFREE $(f, S)$$garrow$

SUBSTSAMPLEPOINT $(f, S)$$do$$Larrow INTvREALRooTIso(g)$if error then$S$の精度をあげる

else
break
return$L$に関する主係数が $S$ でになる場合に，判別式の値は必ず になるために不要な無平方因子分解が発生することがある．それでも，従来の方法では発生していた不要な無平方因子分解の多くを回避できるため，計算効率の改善が実現できる．
ここでは $f$ の$x_{k}$$0$$0$

148
次数 2 の多項式の無平方因子分解

3.2.1

SYMBOLICSQRFREE は代数拡大体上での無平方因子分解を行う．次数 2 の多項式が平方因子をもつと分かっている場合には以下のように置き換えることができる．
$c_{2}x^{2}+c_{1}x+c_{0}=w(2c_{2}x+c_{1})^{2}$

ここで， は$w$$CAD$(1)と符号が同じ値となる定数である．

計算では符号情報のみに興味があるため定数倍の違いは結果に影響しないので，$|w|=1$ として良い．また3.3$c_{2}$.

$c_{2}$の符号は射影因子の符号からすぐに得ることができる．

共通根判定処理

ここでは異なる多項式から得られた分離区間が重複した場合の記号計算回避について述べる．
$f,$$g\in \mathbb{Q}[x_{1}, \ldots, x_{k}]$に標本点$S\in \mathbb{R}^{k-1}$を代入して得られる一変数多項式の分離区間が重複する場合には，記号計算で共通根であるか判定する必要がある．ここでは終結式の符号情報が得られているため，終結式の値が 0 でない場合には記号計算を行うことなく標本点の精度を上げることで実根を分離することが可能となる．

4計算機実験

ここで提案した手法は SyNRAC [6] 上に実装した．表の $SyN2009$ は 2009 年に提案した SyNRAC の実言語で実装したものである．
$SyN2011$ は本稿で提案した内容を装で Maple ユーザ言語で実装している．
他の $CAD$ を基本とする $QE$ ツールである QEPCAD $B1.58[1]$ , Mathematica 8.0 [10] との比較も行った．
良く知られた $QE$ のベンチマーク問題 (adaml から adam3) と，多目的最適化問題を $QE$ を用いて解く$C$問題 (mooea, wilson, lampinen) , 社内で発生した最適化問題 (Optl, opt2) を比較のために使用した．
表 1 は計算時間 (CPU 時間 (秒)) を表している．使用した環境は Intel(R) Core $(TM)i3$ CPU $U330$

1.20 GHz で 2.92 GByte メモリである．2009 年に提案した $SyN$ RAC の実装に比べても記号計算の回避による効率化の実現が確認できる．またいくつかの間題では Mathematica よりも良い結果を得ることができた．

表 2 は $CAD$ 計算における記号計算の呼出回数を表している．符号情報の利用により多くの記号計算を削減できていることが確認できる．
表 3 は SyNRAC 上で実装した，次数 2 の場合の無平方因子分解処理を通常の拡大体上での演算と，式(1) による演算による実行時間を表している．扱った問題では多くの呼出が発生し効率化が実現できていることを確認できる．

5

まとめ
数値数式 $CAD$ 計算における符号情報の利用による記号計算の回避方法について述べた．計算結果を利用することによる記号計算の回避であるため負荷はほとんど発生せず，これにより多くの不要な記号計算を回避でき，大きな効率化を実現できた．
本稿では McCallum 射影演算子に対する手法について述べた力 , 部分終結式主係数 (principal subresultant$\grave{}\grave{}$coefficient) を用いた Collins 射影演算子等を利用するとより多くの情報が得らる．例えば，重根の数が正149

表 1: Computing time $(\sec)$

表 2: 関数の呼出回数確に分かるため数値計算のみで $CAD$ 構築を実現できる可能性がある．ただし，部分終結式主係数を用いた射影演算子は McCallum が提案したものにくらべ，一般に大きな射影因子集合を構築するために全体としてのセルの数が増加してしまう傾向にある．扱う問題や数値手法の効率化にょっては部分終結式主係数を用いたもののほうが良い結果をもたらす可能性もあるが，我々の実験では McCallum 射影演算子を用いたものが良い結果が確認できている．

参考文献[1] C. W. Brown. QEPCAD : program for computing with semi-algebraic sets using CADs. SIGSAM
BULLETIN, 37:97-108, 2003.
$B$

$A$

[2] G. E. Collins. Quantifier elimination and cylindncal algebraic decomposition, pages 8-23. Texts andMonographs in Symbolic Computation. Springer, 1998.
[3] K. Deb. Multi-Objective optimization using Evolutionary Algorethms. Wiley-Interscience Series inSystems and optimization. John Wiley & Sons, Chichester, 2001.

[4] J. D. Dora, C. Dicrescenzo, and D. Duval. About a new method for computing in algebraic numberfields. In Research Contributions from the European Conference on Computer Algebra-Volume 2,EUROCAL 85, pages 289-290, London, $UK$ , 1985. Springer-Verlag.
[5] D. Duval. Algebraic numbers: An example of dynamic evaluation. Journal of Symbolic Computation,
18(5):429-445, 1994.

150

表 3: 次数 2 の無平方因子分解処理[6] H. Iwane, H. Yanami, H. Anai, and K. Yokoyama. An effective implementation of a symbolicnumeric cylindrical algebraic decomposition for quantifier elimination. In Proceedings of the 2009
Intemational Workshop on Symbolic-Numenc Computation, volume 1, pages $5k64$ , 2009.
[7] R. Krawczyk. Newton-algorithmen zur bestimmung von nullstellen mit fehlerschranken. Computing,pages 187-201, 1969.
[8] J. Lampinen. Multiobjective Nonlinear Pareto-optimization. Laboratory of Information Processing,2000.
$F.$[9] S. McCallum. An improved projection operator for cyhndrical algebraic decomposition. In B.
Caviness and J. R. Johnson, editors, Quantifier elimination and cylindncal algebraic decomposition,
Texts and Monographs in Symbolic Computation, pages 242-268. Springer, 1998.

[10] A. W. Strzebo\’{n}ski. Cylindrical algebraic decomposition using validated numerics. Joumal of Symbolic Computation, $41(9):1021-1038$ , 2006.
[11] B. Wilson, D. Cappelleri, T. W. Simpson, and M. Frecker. Efficient Pareto frontier exploration using
surrogate approximations. optimization and Engineemng, 2:31-50, 2001.

A

例

実験で使用した問題を以下に示す．

adaml [10] Stability of Dormand-Prince fifth-order embedded seven-stage method (Example 4.4 fromHong, et al. (1997)$)$$\forall x\forall y((x<0\wedge x^{2}+y^{2}<\frac{99438}{100000})\Rightarrow R(x+iy)R(x-iy)<1)$

,

where
$R(z)=1+z+ \frac{z^{2}}{2}+\frac{z^{3}}{6}+\frac{z^{4}}{24}+\frac{z^{5}}{120}+\frac{z^{6}}{600}.$

adam2 [10] Stability of a six-point upwind-based second-order accurate scheme for approximating a
two-dimensional advection equation (Example 5.4 from Hong, et al. (1997))$2-1 \forall\alpha\forall\beta\forall C_{2}((\alpha\geq 0\wedge\beta\geq 0\wedge4(\alpha^{2}+\beta^{2})<1)\Rightarrow(B\leq 0\vee D\leq 0))$

2–2,$\forall\alpha\forall\beta\forall C_{2}((0\leq\alpha\leq 1\wedge 0\leq\beta\leq 1)\Rightarrow A\leq 0\wedge C\leq 0\wedge(B\leq 0\vee D\leq 0))$,

151
where
$A = C_{2}^{4}(\alpha-\beta+1)(\alpha-\beta-1)(\alpha-\beta)^{2},$$B = 2C_{2}^{4}\beta(3\alpha^{2}\beta-2\alpha^{2}-2\alpha\beta^{2}+\alpha+\beta^{3}-\beta)$$+4C_{2}^{3}\alpha\beta(\alpha^{2}-\alpha+\beta^{2}-\beta)+2C_{2}^{2}\alpha(\alpha^{3}-2\alpha^{2}\beta+3\alpha\beta^{2}-\alpha-2\beta^{2}+\beta)$,$C = C_{2}^{4}\beta^{2}(\beta^{2}-1)+4C_{2}^{3}\alpha\beta^{2}(\beta-1)+\alpha^{2}(\alpha^{2}-1)$$+2C_{2}^{2}\alpha\beta(3\alpha\beta-2\alpha-2\beta+1)+4C_{2}\alpha^{2}\beta(\alpha-1)$,$D = C_{2}^{2}R+2C_{2}S+T,$$R = 8\alpha^{2}\beta^{2}-12\alpha^{2}\beta+5\alpha^{2}-8\alpha\beta^{3}+8\alpha\beta^{2}+2\alpha\beta-4\alpha+4\beta^{4}-4\beta^{3}-3\beta^{2}+4\beta,$$S = 4\alpha^{3}\beta-2\alpha^{3}-4\alpha^{2}\beta^{2}-2\alpha^{2}\beta+\alpha^{2}+4\alpha\beta^{3}-2\alpha\beta^{2}+2\alpha\beta-2\beta^{3}+\beta^{2},$$T = 4\alpha^{4}-8\alpha^{3}\beta-4\alpha^{3}+8\alpha^{2}\beta^{2}+8\alpha^{2}\beta-3\alpha^{2}-12\alpha\beta^{2}+2\alpha\beta+4\alpha+5\beta^{2}-4\beta.$

adam3 [10] Robust multi-objective feedback design (Example 4.2 from Dorato, el at. (1997))
Find the set of satisfying:$n$$\exists q_{1}\exists q_{2}\forall w(q_{1}>1\wedge q_{2}>0\wedge n>0\wedge$$(n-q_{1}^{2})w^{4}+(n((q_{1}+1)^{2}-2q_{2})-(q_{1}^{2}+q_{2}^{2}))w^{2}+(n-1)q_{2}^{2}\geq 0\wedge$$(n-q_{1}^{2})w^{4}+(n((q_{1}-1)^{2}-2q_{2})-(q_{1}^{2}+q_{2}^{2}))w^{2}+(n-1)q_{2}^{2}\geq 0)$.

mooea: extended problem of example 1 in [3, p. 11]$\exists x_{1}\exists x_{2}\exists x_{3} (f=x_{1}^{2}+x_{2}^{2}+x_{3}\wedge g=(x_{1}-1)^{2}+x_{2}^{2}+x_{3}\wedge$$-2\leq x_{1}\leq 2\wedge-2\leq x_{2}\leq 2\wedge-1\leq 10x_{3}\leq 1)$.

wilson [11]$\exists x_{1}\exists x_{2} (f=(x_{1}-2)^{2}+(x_{2}-1)^{2}\wedge g=x_{1}^{2}+(x_{2}-6)^{2}\wedge$$2/5\leq x_{1}\leq 8/5\wedge 2\leq x_{2}\leq 5)$.

lampinen [8, p. 6]$\exists x_{1}\exists x_{2} (f=x_{1}^{2}+x_{2}\wedge g=x_{1}+x_{2}^{2}\wedge-10\leq x_{1}\leq 10\wedge-10\leq x_{2}\leq 10)$.