卒業論文 近似グレブナー基底の計算アルゴリズム 田中 裕子

卒業論文
近似グレブナー基底の計算アルゴリズム
田中 裕子
学籍番号 07251771
2011 年 2 月 14 日
奈良女子大学理学部
情報科学科
近似グレブナー基底の計算アルゴリズム 1
田中 裕子
学籍番号 07251771
内容梗概
本研究では、数式処理システム Asuka を用いた数式処理システムの作成を
行う。
数式処理システムで多変数の連立代数方程式を解く場合、直接解くのは多
くの場合難しい。多変数の連立代数方程式を解くには、等価でかつ解きやす
い形の別の連立代数方程式を求める必要があり、そこに、表れる多項式がグ
レブナー基底である。与えられた方程式のグレブナー基底を計算し、それら
を解くことで元の連立代数方程式の解を求めることができる。しかし、この
計算を厳密に行うには、かなり大変な計算が必要である。また、係数などに
誤差を含む場合には、厳密な代数的アルゴリズムは使えない。
本研究では、近似計算によってグレブナー基底を求めるアルゴリズムを、
数式処理システム Asuka 上で実装を行った。
1
奈良女子大学理学部情報科学科卒業論文, 2011 年 2 月 14 日.
i
目次
第１章 序論
1
第２章 基本事項
2
2.1 数式処理システム Asuka . . . . . . . . . . . . . . . . . . . . .
2
2.2 単項式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2.3 順序 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2.4 グレブナー基底 . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.5 簡約 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.5.1
M-簡約 . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.5.2
S-多項式 . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.6 多項式の精度 . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.6.1
有効浮動小数 . . . . . . . . . . . . . . . . . . . . . . .
6
2.6.2
多項式の計算精度 . . . . . . . . . . . . . . . . . . . . .
6
2.6.3
許容誤差 . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.7 ブッフバーガーのアルゴリズム . . . . . . . . . . . . . . . . .
7
グレブナー基底の計算アルゴリズム . . . . . . . . . . .
8
2.8 近似グレブナー基底 . . . . . . . . . . . . . . . . . . . . . . .
8
近似グレブナー基底の計算アルゴリズム . . . . . . . .
9
2.7.1
2.8.1
第３章 実装
10
3.1 近似グレブナー基底の計算アルゴリズム . . . . . . . . . . . . 10
3.2 プログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2.1
計算精度を求める関数 . . . . . . . . . . . . . . . . . . 11
3.2.2
M-簡約するための関数 . . . . . . . . . . . . . . . . . . 12
3.2.3
S-多項式を求める関数 . . . . . . . . . . . . . . . . . . 12
3.2.4
先頭係数・先頭ベキ乗項を求める関数
3.2.5
ペアをつくる関数 . . . . . . . . . . . . . . . . . . . . . 14
ii
. . . . . . . . . 13
3.2.6
appGB() のプログラム . . . . . . . . . . . . . . . . . . 15
3.2.7
Step1 のプログラム . . . . . . . . . . . . . . . . . . . . 16
3.2.8
Step2 のプログラム . . . . . . . . . . . . . . . . . . . . 17
3.2.9
Step3 のプログラム . . . . . . . . . . . . . . . . . . . . 18
3.2.10 実装例 . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.11 結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
第４章 おわりに
22
謝辞
23
参考文献
24
iii
第１章 序論
本研究では、数式処理システム Asuka を用いた数式処理システムの作成を
行う。
数式処理システム Asuka は、近似代数計算を実行するために作成された
システムである。
連立代数方程式の解法としては、現在のところ、多項式イデアルのグレブ
ナー基底を計算する方法が最も有効であるとされている。しかし、この計算
はしばしば膨大な係数成長を起こすことが知られている。また、丸め誤差や
項消去による桁落ち誤差などで、係数に誤差を含む場合、厳密に計算を行え
なくなる。
本研究では、数式処理システム Asuka を使って、この計算誤差等を制
御しながら、近似計算を行う上で必要な機能を持たせている。
1
第２章 基本事項
2.1
数式処理システム Asuka
Asuka は、近似代数計算を実行するために作成されたシステムである。
厳密な計算を行うことも可能であるが、近似的な計算を行うために、浮動
小数点数や多倍長の浮動小数点数に加えて、区間数、有効浮動小数といった
数値型を用意している。
これらの数値型を使うことによって、近似代数計算による、近似的な（数
値計算）を実行することができる。
2.2
単項式
定義 2.1 (単項式)
α = (α1 , α2 , · · · , αn ) に対して、xα = xα1 1 xα2 2 · · · xαnn ∈ k[x1 , x2 , · · · , xn ] とお
き、これを単項式という。
2.3
順序
定義 2.2 (単項式の順序)
順序集合 Nn に順序 = を一つ定めたとき、α = (α1 , α2 , · · · , αn ), β = (β1 , β2 , · · · , βn ) ∈
Nn とするとき、単項式 xα と xβ の間の順序を
xα = xβ ⇐⇒ α = β
と定義する。Nn の単項式順序 = で決められた単項式の順序を単に単項式順
序という。代表的なものとして、以下のようなものがある。
2
定義 2.3 (辞書式順序 , lex, lexcographic order, >lex )
α = (α1 , α2 , · · · , αn ), β = (β1 , β2 , · · · , βn ) とするとき、
α >lex β
⇐⇒ αi , βi を左から見ていって最初の異なるαi , βi でαi > βi .
定義 2.4 (次数付き辞書式順序 , grlex, graded lexcographic order, >grlex )
α >grlex β
⇐⇒ (|α| > |β|) または (|α| = |β| かつα >lex β).
定義 2.5 (次数付き逆辞書式順序 , grevlex, graded reverse lexcographic order, >grevlex )
α = (α1 , α2 , · · · , αn ), β = (β1 , β2 , · · · , βn ) とするとき、
α >grevlex β ⇐⇒ (|α| > |β|) または (|α| = |β| かつ、αi , βi を右から見て
いって最初の異なるαi , βi でαi < βi ).
{3, 0, 0}, {2, 0, 2}, {1, 2, 1} を大きい順に並べる
lex
grlex
grevlex
(3,0,0) (2,0,2)
(1,2,1)
(2,0,2) (1,2,1)
(2,0,2)
(1,2,1) (3,0,0)
(3,0,0)
例）同じ式の単項式の降べき（大きい順）表示
−5x3 + 7x2 z 2 + 4xy 2 z + 4z 2 (lex)
= 7x2 z 2 + 4xy 2 z − 5x3 + 4z 2 (grlex)
= 4xy 2 z + 7x2 z 2 − 5x3 + 4z 2 (grevlex)
3
2.4
グレブナー基底
多項式 f =
∑
α∈Nn
aα xα について f 6= 0 のとき次のように定義する。
deg(f ) = multideg(f ) = max{α|aα 6= 0}. · · · 多重次数
F を k[x1 , x2 , · · · , xn ] の部分集合とするとき
deg(F ) = {deg(f )|f ∈ F, f 6= 0}
と定義する。
例）f = −5x3 + 7x2 z 2 + 4xy 2 z + 4z 2 について lex で考えると
deg(f ) = (3, 0, 0)
定義 2.6 (グレブナー基底) イデアル I の有限部分集合 G で
deg(I) = h deg(G) i
となるとき G を I のグレブナー基底あるいは標準基底という。
イデアル I のグレブナー基底 G は I の基底である。すなわち I = hGi で
ある。
4
2.5
2.5.1
簡約
M-簡約
多項式 f =
∑
α∈Nn
aα xα について f 6= 0 のとき次のように定義する。
LC(f ) = adeg(f ) · · · f の先頭係数 (leading coef f icient)
LM (f ) = xdeg(f ) · · · f の先頭単項式（leading monomial）
LT (f ) = adeg(f ) xdeg(f ) · · · f の先頭項 (leading term)
とすると、M-簡約の式は以下で表される。
M red(f, g) = LC(g)f − LC(f )
2.5.2
LM (f )
g
LM (g)
S-多項式
LCM を LM (f ) と LM (g) の最小公倍数とするとき、
LCM
LCM
f−
g
LT (f )
LT (g)
= f と g の先頭項を相殺させた式
S(f, g) =
を f と g の S 多項式という。
5
2.6
2.6.1
多項式の精度
有効浮動小数
浮動小数点数は、正確な数値の近似値である。数式処理システム Asuka
では、近似の誤差を扱うために、有効浮動小数点数の型を以下のようにして
いる。
{value, error}
有効浮動小数点数は、値部分 (value) と誤差 (error) から構成されている。
例）
> p := efloat(1.0,0.01);
/* 有効浮動小数を作る */ #e(1.000000000000,0.01000000000000)
> efloat_value(p);
1.000000000000
/* 有効浮動小数の値部分 */
> efloat_error(p); /* 有効浮動小数の誤差部分 */
0.01000000000000
2.6.2
多項式の計算精度
多項式 f の精度は以下の式 acc(f ) で求められる。
ef loat error(f ) acc(f ) = ef loat value(f ) 6
2.6.3
許容誤差
近似計算では、多項式の簡約を何度も繰り返すと、上記の acc(f )、計算誤
差が大きくなる。無視できる程度の小さい誤差であれば、その多項式をその
まま使うことができる。しかし、求める解に影響を与えるくらい大きな誤差
の場合、その多項式を使うことはできない。
そこで、許容誤差を設定する必要ある。計算精度の 10 倍、100 倍などのよ
うに設定する。
例）
> E := 0.001;
/* 多項式の初めの計算精度 */
9.99999999999999999e-4
> Eapp := E*10; /* 許容誤差 （E の 10 倍） */
9.999999999999999998e-3
2.7
ブッフバーガーのアルゴリズム
グレブナー基底を求めるアルゴリズムとしては、ブッフバーガーが発見し
たブッフバーガーアルゴリズムがあり、数式処理の分野での連立代数方程式
の解法として使われている。グレブナー基底の計算はユークリッドの互除法
を多変数多項式に一般化したようなアルゴリズムになっている。
7
グレブナー基底の計算アルゴリズム
2.7.1
Input : F = {f1 , f2 , . . . , fs }
G←F
F を G に入力
repeat
0
0
G ←G
(G = G)
F からペアをつくる (fi , fj ) (1 ≤ i < j < s)
すべての多項式のペア fi , fj ∈ G について
fi , fj の S − 多項式を G で簡約化したものを h とし、
h 6= 0
h=0
ならば G ← G ∪ H として繰り返し
ならば次のペアを処理する
until G = G
0
return G
Output : G = {g1 , g2 , . . . , gt }
この操作は有限回で止まり、得られた G は F で作られるイデアルの基底で
ある。しかし、ブッフバーガー算法で求まるグレブナー基底には一意性がな
く，冗長な元が含まれている（計算途中で生成元が膨れ上がり，最終的にグ
レブナー基底として選ばれる元は，計算途中で求めた生成元の一部である）。
冗長な要素を取り除き、先頭項の係数を 1 になるようにした基底を簡約グレ
ブナー基底という。
2.8
近似グレブナー基底
グレブナー基底の計算は、一般に非常に大変であり、計算途中で出てくる
数も非常に大きくなることがある。このため、単純に浮動小数計算のような
近似計算を行うのではうまく動かない。Asuka では、計算誤差等を制御しな
がら、近似計算を行う上で必要な機能を持たせている。
8
2.8.1
近似グレブナー基底の計算アルゴリズム
Input : F = {f1 , f2 , . . . , fs } εapp = 0.01;（εapp は許容誤差）
G←F
F を G に入力
ε = 0.001; εは多項式の精度
F からペアをつくる (fi , fj ) (1 ≤ i < j < s)
すべての多項式のペア fi , fj ∈ G について
h = Spol(fi , fj ) を計算
h を G で簡約（M − 簡約）する
while hh を G で簡約できる i do
If h が近似的に 0 −→ 次のペアを処理
If h の誤差 > εapp
−→ 零点（次のペアを処理）
else If h の誤差 < ε −→ G ← G ∪ {h}
else −→ Ψ ← Ψ ∪ {h} （ ∗ ∗ ∗の (3) 参照 ）
Ψの中で、誤差が最小の要素を G に追加
このときεapp を追加する式の誤差に変更
od;
return G
Output : G = {g1 , g2 , . . . , gt }
本研究では、近似グレブナー基底を計算するアルゴリズムとして、ブッフバー
ガーのアルゴリズムに
----------------------------- *** ----------------------------(1) ”ゼロ規準”を”近似ゼロ規準”に変更
(2) ”簡約”を”近似的に簡約”に変更
(3) S − 多項式を計算したときに、指定した値より誤差が大きい結果につ
いては、一時的に別のリストとしておき、後で誤差が小さいものから
順に G に追加して、計算を行う
--------------------------------------------------------------という変更を加えたものを考え、Asuka 上に実装した。
9
第３章 実装
3.1
近似グレブナー基底の計算アルゴリズム
手順
appGB(Φ, εapp ==)
Ψ := {};
Step1 :
ε := εinit ;
0
W hile ∀F ∈ ΦがΦ = Φ\F によって簡約できる間
0
Φ
do F · · · F̃ ;
If F 6= F̃ then Φ := Φ − {F };
If acc(F̃ ) > εapp then 零点
else if acc(F̃ ) = ε then Φ := Φ ∪ {F̃ }
else
Ψ := Ψ ∪ {F̃ };
od;
If ](Φ) = 1 then goto Step3;
Step2 :
Φから得られる新しいペア (Fi , Fj ), i < j
Φ
do Spol(Fi , Fj ) · · · F̃ ;
If acc(F̃ ) > εapp then check に行く;
If acc(F̃ ) = ε then
Φ := Φ ∪ {F̃ } で Step1 に行く
else Ψ := Ψ ∪ {F̃ } で check に行く
check :
もし、調べていないペアがあるなら Step2 に行く;
Step3 :
Ψ = {} なら、Φに戻る
F はΨの最も小さな正確性をもった要素とする
F はΨからΦに移動する;
ε := acc(F ) で、Step2 に行く
終了
10
Φ = {F1 , · · · , Fr } ⊂ F[x, y, . . . , z] は与えられた最初の基底とする。
初めにΦ(多項式), εapp（許容誤差）を入力し、手順 appGB が終わったら、グ
レブナー基底を出力する。
3.2
3.2.1
プログラム
計算精度を求める関数
多項式の計算精度を求める関数を以下のように定義した。
proc acc(f);
begin
if zerop(f) then return 0; fi;
if not is_poly(f) then
accF := abs(efloat_error(f)/efloat_value(f));
return accF;
fi;
m := mcoeff(f);
n := length(m)-1;
a := 0;
for i := 0 to n do
b := acc(m[i]);
if a < b then a := b; fi;
od;
return a;
end;
11
3.2.2
M-簡約するための関数
M-簡約の関数を以下のように定義した。
proc M_reduction(f, g);
begin
f := gremainder(f, g); return f;
end;
proc M_red(f, G);
begin
foreach a in G do
f := M_reduction(f, a);
od;
return f;
end;
3.2.3
S-多項式を求める関数
S-多項式の関数を以下のように定義した。
proc Spol(f,g);
begin
lcm := (a,b) -> a*b/gcd(a, b);
F := f; G := g;
Fh := mvar(F);
LMf := lterm(F,Fh);
K := lcm(LMf,LMg);
LCf := lcoeff(F,Fh);
LTf := gquotient(LMf,LCf);
Gh := mvar(G);
LMg := lterm(G,Gh);
LCg := lcoeff(G,Gh);
LTg := gquotient(LMg,LCg);
12
S := gquotient(K,LTf*LCf)*F - gquotient(K,LTg*LCg)*G;
return S;
end;
3.2.4
先頭係数・先頭ベキ乗項を求める関数
proc lc(f);
begin
/* leading coefficient */
while is_poly(f) do
f := minit(f);
od;
return f;
end;
proc mpower(f); /* leading power */
begin
if not is_poly(f) then return 1; fi;
d := mdeg(f);
v := mvar(f);
return v^d;
end;
proc lp(f);
begin
p := 1;
/* leading power product */
while is_poly(f) do
p := p*mpower(f);
f := minit(f);
od;
return p;
end;
13
ペアをつくる関数
3.2.5
proc Pair_make(Pa);
begin
q := Pa; Pair := {};
while
not(q = {}) do
Fi := first(q);
q := rest(q);
foreach Fj in q do
Pair := append(Pair, {{Fi, Fj}});
od;
od;
return Pair;
end;
proc Pair_make1(Ft,P);
begin
F := Ft; Q := {}; q := P;
while not(q = {}) do
fp := first(q);
q := rest(q);
if (fp = F) then ;
else Q := append(Q,{{F,fp}});
fi;
od;
foreach a in P do
Q := append(Q, {{Ft, a}});
od;
return Q;
end;
14
3.2.6
appGB() のプログラム
/* main
*/
proc appGB(Phi, Eapp);
begin
E := Eapp/10;
b := {};
foreach a in Phi do
b := append(b, {efloat(a, E)});
od;
Phi:= b;
R1 := Step1(Phi,Eapp); /* go to Step1 */
while true do
c := Step2(Phi, R1, E, Eapp);
Phi := first(c);
Psi := second(c);
a := Step3(Psi);
/* go to Step2 */
/* go to Step3 */
if a = {} then break; fi;
Phi := append(Phi, {first(a)});
E := second(a);
od;
return
end;
/* from Psi to Phi */
/* E = acc(F)
;
15
go to Step2 */
3.2.7
Step1 のプログラム
/* Step1
proc Step1(Phi,Eap);
begin
P := Phi;
E := Eap/10;
q := P;
*/
Psi := {};
Eapp := Eap;
Pair := {};
Pair := Pair_make(q);
printf("pair =%p\n", Pair);
while not(Pair = {}) do
FF := first(Pair);
f := first(FF);
g := second(FF);
Pair := rest(Pair);
printf("F1 = %p \n",f);
printf("F2 = %p \n",g);
Ft := Mred(f,g);
Ft1 := first(Ft);
accFt := second(Ft);
if not(f = Ft1 or f = -Ft1) then P := setdiff(P,{f});
fi;
if
accFt > Eapp + 0.00000001 then
;
elif
accFt <= E + 0.00000001
then
P := append(P, {Ft1});
else
Psi := append(Psi,{Ft1});
fi;
od;
return P;
end;
16
3.2.8
Step2 のプログラム
/*
Step2
*/
proc Step2(Phi, R1, E, Eapp);
begin
Psi := {};
GP := R1;
P := {};
GPsi := {};
P := append(P,Phi);
P := append(P,GP);
P := kaburi(P);
Pair := {};
Pair := Pair_make(P);
while not(Pair = {}) do
FF := first(Pair);
f := first(FF);
g := second(FF);
Pair := rest(Pair);
printf("F1 = %p \n",f);
printf("F2 = %p \n",g);
Ft := Spol(f,g);
Ft := M_red(Ft, Phi);
accFt := acc(Ft);
if zerop(Ft) then ;
elif
accFt > Eapp+0.00000001 then ;
elif
accFt <= E+0.00000001 then
Q := Pair_make1(Ft, P);
Pair := append(Pair, Pair_make1(Ft,P));
P := append(P, {Ft});
else
Psi := append(Psi, {Ft});
fi;
od;
printf("P = %p \n",P);
17
printf("Psi = %p \n",Psi);
GPsi := append(GPsi,{P});
GPsi := append(GPsi,{Psi});
printf("GPsi = %p \n",GPsi);
return GPsi;
end;
3.2.9
Step3 のプログラム
/* Step3 */
proc Step3(Psi);
begin
if Psi = {} then return {};
fi;
d := first(Psi);
c := acc(d);
Psi := rest(Psi);
foreach a in Psi do
b := acc(a);
if c > b then
d := a;
c := b;
fi;
od;
return {d, c};
end;
18
3.2.10
実装例
以下の、楕円と曲線との交点を求める連立代数方程式について、近似グレ
ブナー基底を求める。なお、項順序は辞書式順序で行った。

3
2
 f 1 = 1.02x + 0.92y − 3.11x − 3.12y − 10.92

f 2 = 1.06x2 + 1.20y 2 − 8.7
/*
ターミナルで入力
*/
>appGB({1.02*x^3+0.92*y^2-3.11*x-3.12*y-10.92,
1.06*x^2+1.20*y^2-8.7},0.01);
3.2.11
結果
実装すると、以下の 6 つのグレブナー基底が出力される。値部分のみ示し、
誤差部分は省略してある。
G1 = 1.020000000000*x^3 - 3.110000000000*x + 0.9200000000000*y^2
- 3.120000000000*y -10.92000000000
G2 = 1.060000000000*x^2 + 1.200000000000*y^2 - 8.700000000000
G3 = (-1.132075471698*y^2 + 5.158527561968)*x
+ 0.9019607843137*y^2
- 3.058823529412*y - 10.70588235294
G4 = (+
+
+
3.602716981132*y - 7.768747169811)*x^2
(- 3.219075258099*y^2 - 2.870400000000*y
8.478748451406)*x + 1.202648629405*y^4
4.078547525810*y^3 - 19.75502456390*y^2
18.58471484514*y + 65.04650195799
G5 = (3.890823529412*y + 8.390007996924)*x
- 1.630188679245*y^4 + 18.33306164683*y^2
+ 3.099943713956*y - 43.00522474810
19
G6 =
/*
+
+
563.9368745249*y^7
8859.163330330*y^5
49133.84856210*y^3
93296.65630034*y -
+ 1216.049726043*y^6
- 16958.77794822*y^4
+ 70719.95131383*y^2
85189.32431881
以下 G3 ~ G6 の計算誤差
全て Eapp(0.01) 以下
> acc(G3);
0.001086956521739
> acc(G4);
0.002142050840471
> acc(G5);
0.001531233203620
> acc(G6);
0.004194422248356
20
*/
よって以下

3
2
 f 1 = 1.02x + 0.92y − 3.11x − 3.12y − 10.92

f 2 = 1.06x2 + 1.20y 2 − 8.7
の連立代数方程式のグレブナー基底を求めると、





































G1 = 1.02x3 + 0.92y 2 − 3.11x − 3.12y − 10.92
G2 = 1.06x2 + 1.20y 2 − 8.7
G3 = (−1.132075471698y 2 + 5.158527561968)x
+0.9019607843137y 2 − 3.058823529412y − 10.70588235294
G4 = (−3.602716981132y − 7.768747169811)x2
+(−3.219075258099y 2 − 2.870400000000y + 8.478748451406)x
+1.202648629405y 4 − 4.078547525810y 3
−19.75502456390y 2 + 18.58471484514y
+65.04650195799
















G5 = (3.890823529412y + 8.390007996924)x − 1.630188679245y 4




+18.33306164683y 2 + 3.099943713956y − 43.00522474810









G6 = −563.9368745249y 7 + 1216.049726043y 6 + 8859.163330330y 5




−16958.77794822y 4 − 49133.84856210y 3 + 70719.95131383y 2



+93296.65630034y − 85189.32431881
上記の６つの多項式（基底）が出力される。
21
第４章 おわりに
本研究では、近似グレブナー基底を求める計算アルゴリズムについて、数
式処理システム Asuka を使って実装を行った。このプログラムを使って計
算を行ったが、計算誤差によって正しく計算できない場合が存在した。
本質的にグレブナー基底が不安定になる場合に、どのような処理を行うか
について、より詳しく調べる必要がある。
22
謝辞
本論文の執筆及び研究にあたり、様々な御指導を賜わった本学科加古 富志
雄 教授 に厚く御礼申し上げます。 また、助言、 ご協力をして頂いた加古研
究室の皆さんに深く感謝いたします。
23
参考文献
[1] 加古富志雄（著） 「数式処理システム Asuka」 講義資料
[2] 佐々木建昭（著）（筑波大学教授）「近似グレブナー基底とその応用」
講義資料
[3] D. コックス, J. リトル, D. オシー（著）, 落合啓之, 示野信一, 西山亨, 室
政和, 山本敦子（訳）「グレブナ基底と代数多様体入門（上）」, 第２版,
シュプリンガー・フェアラーク東京株式会社, 2000. 参考図書
24