2021-06-02

グレブナー基底で連立多項式を解く、量子計算機でグレブナー基底を探す

https://arxiv.org/pdf/1903.08270.pdf

pythonのsympyを使うとグレブナー基底が計算できる
sagemath on pythonでもできる

from sympy import *
x,y,z = symbols('x y z')
eqs =[x**2 + y**2 + z**2 - 4, x**2 + 2 * y**2 - 5, x*z-1]
result = groebner(eqs,x,y,z,order='lex')
for eq in list(result):
    print(eq)
    print('\n')

グレブナー基底が以下のように算出される

x + 2*z**3 - 3*z
y**2 - z**2 - 1
2*z**4 - 3*z**2 + 1

これにより、 $x^2 + y^2 + z^2=0,x^2 + 2 y^2-5=0,xz-1=0$ を満足する零点集合が $x + 2z^3 -3z =0,y^2-z^2-1=0,2z^4-3z^2+1=0$ の零点しゅうごうであることがわかる
sagemathではどうなるか

R = PolynomialRing(ZZ,3,"xyz",order="lex") # 3変数の整数係数多項式環を作る。3変数はx,y,zでlex(辞書式)順序
x,y,z = R.gens() # x,y,zという変数をpython/sagemathで使うときに、作成した環 R の生成変数として指定する
eq1 = x^2 + y^2 + z^2 - 4 # イデアルの生成式の1番
eq2 = x^2 + 2 * y^2 - 5 # 同2番
eq3 = x * z -1 # 同3番
I = (eq1,eq2,eq3) * R  # イデアル I は、生成式に多項式環の要素すべてをかけたものなので、このようにする
B = I.groebner_basis() # イデアルIのグレブナー基底を計算する
B

[x + 2*z^3 - 3*z, y^2 - z^2 - 1, 2*z^4 - 3*z^2 + 1]

これを「見て」みる
Rを使って、xyzグリッド点を発生し、連立方程式をだいたい満足する点がかさなることを図示し、また、zの値が0.5のときに、3つの方程式が0に近いかどうかを絶対値で表したときの値が、２つの連立方程式でどのような関係になるかを図示する。値は違うが、０に近い(x,y,z=0.5)座標は共有されることがわかる

f:id:ryamada22:20210602112707p:plain

x <- y <- z <- seq(from=-2,to=2,length=50)
xyz <- as.matrix(expand.grid(x,y,z))
v1 <- xyz[,1]^2 + xyz[,2]^2 + xyz[,3]^2 -4
v2 <- xyz[,1]^2 + 2 *xyz[,2]^2 -5
v3 <- xyz[,1] * xyz[,3] -1

V <- v1**2 + v2**2 + v3**2

library(rgl)

#plot3d(xyz[which(V.<1),])

#plot3d(xyz[which(v1 < 0.0001),])

v1. <- xyz[,1] + 2 * xyz[,3]^2 -3 * xyz[,3]
v2. <- xyz[,2]^2 - xyz[,3]^2 - 1
v3. <- 2 * xyz[,3]^4 - 3 * xyz[,3]^2 + 1

V. <- v1.**2 + v2.**2 + v3.**2

plot(V[which(V<0.5)],V.[which(V<0.5)])

plot3d(xyz[which(V.<0.5),])
spheres3d(xyz[which(V<1),],radius=0.01,color="green")

f:id:ryamada22:20210602112729p:plain

x <- y <- seq(from=-2,to=2,length=50)
xy <- as.matrix(expand.grid(x,y))
z <- 0.5

xyz <- cbind(xy,rep(z,length(xy[,1])))
v1 <- xyz[,1]^2 + xyz[,2]^2 + xyz[,3]^2 -4
v2 <- xyz[,1]^2 + 2 *xyz[,2]^2 -5
v3 <- xyz[,1] * xyz[,3] -1

V <- v1**2 + v2**2 + v3**2

library(rgl)

#plot3d(xyz[which(V.<1),])

#plot3d(xyz[which(v1 < 0.0001),])

v1. <- xyz[,1] + 2 * xyz[,3]^2 -3 * xyz[,3]
v2. <- xyz[,2]^2 - xyz[,3]^2 - 1
v3. <- 2 * xyz[,3]^4 - 3 * xyz[,3]^2 + 1

V. <- v1.**2 + v2.**2 + v3.**2

plot(V[which(V<0.5)],V.[which(V<0.5)])

plot3d(cbind(xyz[,1:2],V))
spheres3d(cbind(xyz[,1:2],V.),radius=0.1,color="green")

{0,1}を取る変数があるときに、それを満足する解を探すには、x(x-1)=0を零点集合の満足条件に加えればよい
- sagemathでやってみる

R = PolynomialRing(ZZ,6,"x1,x2,x3,a2,b1,b2",order="lex")
x1,x2,x3,a2,b1,b2 = R.gens()
eqs = [x1+x2+x3-b1,x1+a2*x2-b2]
eqs = eqs + [x*(x-1) for x in [x1,x2,x3]] # 3変数 x1,x2,x3につき xi(xi-1)=0の条件を付与
eqst = tuple(eqs) # リストをタプルにする
I = eqst * R 
B = I.groebner_basis()
for eq in list(B):
    print(str(eq)+",")

- グレブナー基底は「変数順序での掃き出し」になっているので、後半には、xi変数は入っていない(最後の７つの式)。これが、 $x1,x2,x3 \in \{0,1\}$ を満足する条件の下での、a2,b1,b2の満足条件となる

x1 + x2 + x3 - b1,
x2^2 - x2,
x2*x3 - x2*b1 + x2*b2 - x3 + b1 - b2,
x2*a2 - x2 - x3 + b1 - b2,
x2*b1^2 - 3*x2*b1 + 2*x2 - 4*x3*b1 + 4*x3*b2 + 4*x3 + 2*a2*b1^2 - 4*a2*b1*b2 - 2*a2*b1 + 4*a2*b2 + b1^3 - 4*b1^2*b2 - b1^2 + 4*b1*b2^2 + 4*b1*b2 - 4*b2^2,
x2*b2^2 + x2*b2 - 2*x2 - 4*x3*b1 + 6*x3*b2 + 2*x3 + 3*a2*b1^2 - 6*a2*b1*b2 - 3*a2*b1 + 6*a2*b2 + 2*b1^3 - 6*b1^2*b2 - 4*b1^2 + 6*b1*b2^2 + 6*b1*b2 + 4*b1 - 7*b2^2 - b2,
2*x2*b2 - 2*x2 + 2*x3*b1 - 4*x3 - b1^2 + 3*b1 - 2*b2,
x3^2 - x3,
x3*a2 + 2*x3*b1 - 2*x3*b2 - 2*x3 - a2*b1 + a2*b2 - b1^2 + 2*b1*b2 + b1 - b2^2 - b2,
x3*b1^2 - x3*b1 - 2*x3*b2 - a2*b1^2 + 2*a2*b1*b2 + a2*b1 - 2*a2*b2 - b1^3 + 2*b1^2*b2 + 2*b1^2 - 2*b1*b2^2 - 2*b1*b2 - b1 + 2*b2^2,
2*x3*b1*b2 - 4*x3*b2 - b1^2*b2 + 3*b1*b2 - 2*b2,
6*x3*b1 - 6*x3*b2 - 6*x3 - 3*a2*b1^2 + 6*a2*b1*b2 + 3*a2*b1 - 6*a2*b2 - 2*b1^3 + 6*b1^2*b2 + 3*b1^2 - 6*b1*b2^2 - 6*b1*b2 - b1 + 6*b2^2,
x3*b2^2 - x3*b2 - a2*b2^2 + a2*b2 - b1*b2^2 + b1*b2 + b2^3 - b2,
a2^2*b1^2 - 2*a2^2*b1*b2 - a2^2*b1 + 2*a2^2*b2 + a2*b1^2 - 2*a2*b1*b2 - a2*b1 + 4*a2*b2^2 - 2*a2*b2 - 3*b1^2*b2^2 + 3*b1^2*b2 + 2*b1*b2^3 + 9*b1*b2^2 - 11*b1*b2 - 6*b2^3 - 2*b2^2 + 8*b2,
a2^2*b2^2 - a2^2*b2 - 2*a2*b2^3 + 3*a2*b2^2 - a2*b2 + b2^4 - 2*b2^3 + b2^2,
a2*b1^3 - 3*a2*b1^2 + 2*a2*b1 + b1^3 - 3*b1^2*b2 - 3*b1^2 + 9*b1*b2 + 2*b1 - 6*b2,
a2*b1^2*b2 - 3*a2*b1*b2 + 2*a2*b2 - b1^2*b2^2 + b1^2*b2 + 3*b1*b2^2 - 3*b1*b2 - 2*b2^2 + 2*b2,
2*a2*b1*b2^2 - 2*a2*b1*b2 - 4*a2*b2^2 + 4*a2*b2 + b1^2*b2^2 - b1^2*b2 - 2*b1*b2^3 - b1*b2^2 + 3*b1*b2 + 4*b2^3 - 2*b2^2 - 2*b2,
b1^4 - 6*b1^3 + 11*b1^2 - 6*b1,
b1^3*b2 - 6*b1^2*b2 + 11*b1*b2 - 6*b2,

$y1+2y2+3y3+3y4-z=0,y1+y+2y3+y4-3=0$ , $yi \in \{0,1\}$ を満足するような整数zを求め、その最小値を知りたい…

R = PolynomialRing(QQ,5,"y1,y2,y3,y4,z",order="lex")
y1,y2,y3,y4,z = R.gens()
eqs = [y1+2*y2+3*y3+3*y4-z,y1+y2+2*y3+y4-3]
eqs = eqs + [x*(x-1) for x in [y1,y2,y3,y4]]
eqst = tuple(eqs)
I = eqst * R 
B = I.groebner_basis()
lasteq = B[-1]
lasteq
Z = var("z")
solve(Z^3-15*Z^2+74*Z-120,Z)

[y1 + y3 - 1/2*z^2 + 11/2*z - 16, y2 + y3 + z^2 - 10*z + 23, y3^2 - y3, y3*z - 6*y3 - z + 6, y4 - 1/2*z^2 + 9/2*z - 10, z^3 - 15*z^2 + 74*z - 120]

z^3 - 15*z^2 + 74*z - 120

[z == 4, z == 5, z == 6]

- なので、満足するzのうち、最小なのは4
トーリックイデアルのグレブナー基底を探す
- 行列 $A= \begin{pmatrix} 4 & 5 & 1 & 0\\ 2 & 3 & 0 & 1 \end{pmatrix}$ が与えられたとき、これがトーリックイデアルを指定する。Aは配置と呼ばれる
- $A \begin{pmatrix}x \\ y \\ z \\ w \end{pmatrix} = \begin{pmatrix} 4x+5y + z \\ 2x + 3y + w \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$ を考える
- このような $(x,y,z,w)$ は $(x,y,z,w) = a(1,0,-4,-2) + b(0,1,-5,-3)=(a,b,-4a-5b,-2a-3b)$ と、a,bを使って表せる。a,bを整数に取れば、特定の格子点であって、式を満足する点ということになる
- ここで、 $(a,b,-4a-5b,-2a-3b)$ が、例えば、 $(-2,1,3,1)$ だったときに、この長さ４のベクトルの正の要素と負の要素を別々に扱って、 $y z ^3 w - x^2$ というような２個の単項式の差の形をとる多項式に対応付けることにする
- このような多項式(が複数あるときに、そ）の零点集合を考えることがトーリックイデアル
- 整数ベクトルはその要素の符号を気にすることで、２個の単項式の差の式を取れることが、トーリックイデアルの特徴である
- 今、このように制約のある多項式集合(イデアル)を、いい感じの生成元多項式のセットに結び付けたい
- グレブナー基底はそういういい感じの生成元多項式のセット
- グレブナー基底は、単項式に順序をつけることで、取られ方が変わる
量子計算機では、この探索を、整数格子点最適化問題にして、探索する(らしい)
- 格子点を量子ビット表現して、あるコスト関数を最小にするような格子点を探す問題として求解する
- コスト関数としては、変数(x,y,z,w)に重みをつけることで、単項式の辞書式順序がコストになるようにする。また、量子計算機が計算しやすいように、辞書式順序がコストになるような値を、(x,y,z,w)の重みベクトルに行列を掛けて変換し、その変換してできたベクトルの二乗ノルムをコスト関数とする
- この話が、本記事の冒頭で引用したペイパーの節６の話
- ちなみに、この配置Aにまつわるグレブナー基底をパイソン計算するコードは以下。w,z,y,xの順のグレブナー基底を多項式セットとしてw,z,y,x順でのグレブナー基底を求めると、グレブナー基底そのものが得られるし、その多項式セットに、元のイデアルに含まれる多項式を付け加えても同じこと。x,y,z,wの順のグレブナー基底を求めると、もとの設定が単純だったこともあり、２つの多項式が基底として得られている

from sympy import *
x,y,z,w = symbols('x y z w')
eqs =[x*z*w - y,y**2*z**2-x**3,x**4*w-y**3*z,z**4*w**2 - x,y*z**3*w -x**2]
eqs1 = eqs + [x - z**4 * w**2,y - z**5 * w**3,x * y - z**9 * w**5,y - x*z*w]
result = groebner(eqs,w,z,y,x,order='lex')
result1 = groebner(eqs1,w,z,y,x,order='lex')
for eq in list(result):
    print(eq)
    print('\n')

print("----")
for eq in list(result1):
    print(eq)
    print('\n')
print("=====")
result = groebner(eqs,x,y,z,w,order='lex')
result1 = groebner(eqs1,x,y,z,w,order='lex')
for eq in list(result):
    print(eq)
    print('\n')
print("----")
for eq in list(result1):
    print(eq)
    print('\n')

w**2*z**4 - x


w*y*z**3 - x**2


w*x*z - y


w*x**4 - y**3*z


-x**3 + y**2*z**2


----
w**2*z**4 - x


w*y*z**3 - x**2


w*x*z - y


w*x**4 - y**3*z


-x**3 + y**2*z**2


-w**2*z**4 + x


-w**3*z**5 + y


----
-w**2*z**4 + x


-w**3*z**5 + y