量子アニーリングでPell方程式の解を求める

こんにちは、@snuffkin と言います。
これは量子コンピュータ Advent Calendar 2018の22日目のブログです。

はじめに

昔から研究されている数学の問題にPell方程式と呼ばれるものがあります。これは、

平方数でない(= n^2 の形で書けない)自然数 D に対して、\\
x^2 - Dy^2 = 1 \\
を満たす整数解 (x, y) を求めよ

という問題です。
Pell方程式の性質をいくつか挙げてみます。

• $x = 0$ とすると、$1 = -Dy^2 \le 0$ となり矛盾が生じるため、$x \neq 0$ です。
• $y = 0$ とすると、$x^2 = 1$ なので、$(x, y) = (\pm1, 0)$ がどんな $D$ に対しても解となります。これを自明な解と呼びます。
• $x$ と $y$ は2次の項しかないため、$(\alpha, \beta)$ が解の場合は、$(-\alpha, \beta), (\alpha, -\beta), (-\alpha, -\beta)$ も解になります。

これらの性質により、自然数解($x, y \ge 1$ となる解)が分かれば他の整数解が分かります。

実はQuantum Algorithm Zoo全訳の「Algorithm: ペル方程式」の部分に記載されているように、Pell方程式の解を多項式時間で求める一般的な方法は見つかっていません(古典コンピュータで高速に求めるアルゴリズムは見つかっていない)。

そこで、最近話題の量子アニーリングでPell方程式の解を求めてみましょう(といっても、実機ではなくシミュレータを利用します)。
このブログでは、$D=2$ とした、

x^2 - 2y^2 = 1

を解きます。手計算で求められる自然数解の1つとしては $(x, y) = (3, 2)$ があります($3^2 - 2\cdot2^2 = 1$)。

Pell方程式の解を求めてみましょう

Pell方程式をQUBO(Quadratic Unconstrained Binary Optimization)と呼ばれる表現に変換し、QUBOの最適化問題を量子アニーリングで解く、という方針で進めます。
QUBOについては、こちらの記事が分かりやすいです。

QUBOを使って日本地図の彩色問題やってみた

最終的には$D=2$ の場合を解きますが、一般化を意識して議論を進めます。

ステップ1: 最適化問題に変換する

量子アニーリングで解けるのは最適化問題(最小値を取る箇所を求める問題)であるため、「方程式の解を求める問題」を「最小値を取る箇所を求める問題」に変換します。
Pell方程式を「$○○=0$」の形に変形すると、

x^2 - Dy^2 -1 = 0

となります。両辺を二乗すると、

(x^2 - Dy^2 -1)^2 = 0

となります。二乗したため、どんな $x, y$ に対しても「$左辺 \ge 0$」です。そのため、

(x^2 - Dy^2 -1)^2 が最小値を取る (x, y) を求めよ

という最適化問題を解けば良い事になります。

ステップ2: x, yをビットで表現する

QUBOで表現する変数は0もしくは1を値として取ることができます。要するにビットを表現することができます。
$x, y$ が2ビット以下で表せる自然数の場合、$x_i, y_i \in \{0, 1\}$を使って、

x = x_0 + 2x_1 \\
y = y_0 + 2y_1 \\

と書けるため、この表現を使いましょう(ここでは、2ビットとしましたが、解の大きさによってはもっと多くのビットが必要です)。
これをPell方程式に代入すると、

\begin{align}
(x^2 - Dy^2 -1)^2 =& ((x_0 + 2x_1)^2 - D(y_0 + 2y_1)^2 -1)^2 \\
=& D^{2} y_{0}^{4} + 8 D^{2} y_{0}^{3} y_{1} + 24 D^{2} y_{0}^{2} y_{1}^{2} + 32 D^{2} y_{0} y_{1}^{3} + 16 D^{2} y_{1}^{4} - 2 D x_{0}^{2} y_{0}^{2} - 8 D x_{0}^{2} y_{0} y_{1} - 8 D x_{0}^{2} y_{1}^{2} \\
& - 8 D x_{0} x_{1} y_{0}^{2} - 32 D x_{0} x_{1} y_{0} y_{1} - 32 D x_{0} x_{1} y_{1}^{2} - 8 D x_{1}^{2} y_{0}^{2} - 32 D x_{1}^{2} y_{0} y_{1} - 32 D x_{1}^{2} y_{1}^{2} + 2 D y_{0}^{2} + 8 D y_{0} y_{1} \\
& + 8 D y_{1}^{2} + x_{0}^{4} + 8 x_{0}^{3} x_{1} + 24 x_{0}^{2} x_{1}^{2} - 2 x_{0}^{2} + 32 x_{0} x_{1}^{3} - 8 x_{0} x_{1} + 16 x_{1}^{4} - 8 x_{1}^{2} + 1
\end{align}

と書くことができます。(さて、計算が激しくなってきましたね)

ステップ3: 二乗以上の変数の次数を減らす

量子アニーリングでQUBOを解くには、2次以下の多項式にする必要があります。そのため、多項式の次数を減らす必要があります。

$x_i, y_i \in \{0, 1\}$であるため、必ず$x_i^2 = x_i$ となります。したがって、$x_i^4 = x_i^3 = x_i^2 = x_i$ となります。同じ変数のべき乗は無視しても構わないです。これは、$y_i$ も同様です。
これをステップ2で求めた式に適用すると、次のようになります。

\begin{align}
(x^2 - Dy^2 -1)^2 =& 
64 D^{2} y_{0} y_{1} + D^{2} y_{0} + 16 D^{2} y_{1} - 32 D x_{0} x_{1} y_{0} y_{1} - 8 D x_{0} x_{1} y_{0} - 32 D x_{0} x_{1} y_{1} - 8 D x_{0} y_{0} y_{1} - 2 D x_{0} y_{0} \\
& - 8 D x_{0} y_{1} - 32 D x_{1} y_{0} y_{1} - 8 D x_{1} y_{0} - 32 D x_{1} y_{1} + 8 D y_{0} y_{1} + 2 D y_{0} + 8 D y_{1} + 56 x_{0} x_{1} \\
& - x_{0} + 8 x_{1} + 1
\end{align}

ステップ4: 3個以上の変数の積を減らす

同じ変数のべき乗はなくなりましたが、$- 32 D x_{0} x_{1} y_{0} y_{1}$ や $- 8 D x_{0} x_{1} y_{0}$ のように3個以上の変数の積が残っています。
そのため、新たな変数に置き換えことで、2個以下の変数の積になるようにします。
ここでは、以下で紹介されている方法を利用します。

N次式をQUBOに変換してアニーリングマシンで扱えるようにする

Pell方程式で必要なのは、$N=4$ の場合です。
新たな変数 $z_0, z_1$ を導入します。
$H_c$ として $H_c(k, m, n) := km -2(k+m)n +3n$ を選びます。$H_c$ を使って $x_{0} x_{1} y_{0} y_{1}$ を表すと、次のようになります(上記のブログの変数 $\lambda$ はここでは $L$ にしています)。

\begin{align}
x_{0} x_{1} y_{0} y_{1} &= z_1 y_1 + L H_c(x_0, x_1, z_0) + L H_c(z_0, y_0, z_1) \\
&= L x_{0} x_{1} - 2 L x_{0} z_{0} - 2 L x_{1} z_{0} + L y_{0} z_{0} - 2 L y_{0} z_{1} - 2 L z_{0} z_{1} + 3 L z_{0} + 3 L z_{1} + y_{1} z_{1}
\end{align}

ただし、$L > 0$ です。
これで「4個の変数の積」を「2個以下の変数の積」で表すことができます。
この式に、$y_1 = 1$ を代入すれば、$x_{0} x_{1} y_{0}$ といった「3個の変数の積」を「2個以下の変数の積」で表すことができます。

\begin{align}
x_{0} x_{1} y_{0} &= L x_{0} x_{1} - 2 L x_{0} z_{0} - 2 L x_{1} z_{0} + L y_{0} z_{0} - 2 L y_{0} z_{1} - 2 L z_{0} z_{1} + 3 L z_{0} + 3 L z_{1} + z_{1} \\
x_{0} x_{1} y_{1} &= L x_{0} x_{1} - 2 L x_{0} z_{0} - 2 L x_{1} z_{0} - 2 L z_{0} z_{1} + 4 L z_{0} + L z_{1} + y_{1} z_{1} \\
x_{0} y_{0} y_{1} &= - 2 L x_{0} z_{0} + L x_{0} + L y_{0} z_{0} - 2 L y_{0} z_{1} - 2 L z_{0} z_{1} + L z_{0} + 3 L z_{1} + y_{1} z_{1} \\
x_{1} y_{0} y_{1} &= - 2 L x_{1} z_{0} + L x_{1} + L y_{0} z_{0} - 2 L y_{0} z_{1} - 2 L z_{0} z_{1} + L z_{0} + 3 L z_{1} + y_{1} z_{1}
\end{align}

これをステップ3で求めた式に適用すると、次のようになります。

\begin{align}
(x^2 - Dy^2 -1)^2 =& 
64 D^{2} y_{0} y_{1} + D^{2} y_{0} + 16 D^{2} y_{1} - 72 D L x_{0} x_{1} + 160 D L x_{0} z_{0} - 8 D L x_{0} + 208 D L x_{1} z_{0} - 32 D L x_{1} \\
& - 80 D L y_{0} z_{0} + 160 D L y_{0} z_{1} + 224 D L z_{0} z_{1} - 288 D L z_{0} - 272 D L z_{1} - 2 D x_{0} y_{0} - 8 D x_{0} y_{1} - 8 D x_{1} y_{0} \\
& - 32 D x_{1} y_{1} + 8 D y_{0} y_{1} + 2 D y_{0} - 104 D y_{1} z_{1} + 8 D y_{1} - 8 D z_{1} + 56 x_{0} x_{1} - x_{0} \\
& + 8 x_{1} + 1
\end{align}

これで晴れて2次式になりましたね。(ちなみに、手計算が大変なので実際にはsympyを使っています)

ステップ5: LとDに具体的な値を入れる

そろそろ目的を見失いそうですが、今回やろうとしていた事を思い出すと、$D=2$ でのPell方程式を解こうとしていたのでした。そこで、$D=2$ としましょう。また、$L=1$ とします。
これをステップ4で求めた式に適用すると、次のようになります。

\begin{align}
(x^2 - Dy^2 -1)^2 =& 
- 88 x_{0} x_{1} - 4 x_{0} y_{0} - 16 x_{0} y_{1} + 320 x_{0} z_{0} - 17 x_{0} - 16 x_{1} y_{0} - 64 x_{1} y_{1} + 416 x_{1} z_{0} \\
& - 56 x_{1} + 272 y_{0} y_{1} - 160 y_{0} z_{0} + 320 y_{0} z_{1} + 8 y_{0} - 208 y_{1} z_{1} + 80 y_{1} + 448 z_{0} z_{1} \\
& - 576 z_{0} - 560 z_{1} + 1
\end{align}

これでようやく2次式になり、QUBOの最適化問題を解くことができます。

ステップ6: 量子アニーリングで解を求める

このブログでは、量子アニーリングのシミュレータとして、D-Wave社が開発しているD-Wave's Ocean Softwareを利用します。

D-Wave's Ocean Software

これはPythonで動作するライブラリで、次のコマンドでインストールできます。

pip install dwave-ocean-sdk

アカウントを取得すればD-Wave社の実機に接続して実行することもできますが、ここでは自分のPCでシミュレータとして実行します。
変数 qubo にはsympyで作成したQUBOの多項式が設定されているとします。

from dimod import *
from dwave_qbsolv import QBSolv
var_list = [x0, x1, y0, y1, z0, z1]

linear = dict()
for var1 in var_list:
    linear[var1] = qubo.coeff_monomial(var1)

quadratic = dict()
for index, var1 in enumerate(var_list):
    for var2 in var_list[index + 1:]:
        quadratic[(var1, var2)] = qubo.coeff_monomial(var1*var2)

bqm = BinaryQuadraticModel(linear, quadratic, qubo.coeff_monomial(1), dimod.BINARY)
result = QBSolv().sample(bqm)
for s, e in result.data(['sample','energy']):
    print(s,'\t E = ', e)

実行結果はこちら

{x0: 1, x1: 1, y0: 0, y1: 1, z0: 0, z1: 1} 	 E =  -928.0

この結果から $x, y$ を求めると、

x = x_0 + 2x_1 = 1 + 2\cdot1 = 3\\
y = y_0 + 2y_1 = 0 + 2\cdot1 = 2

となり、期待した結果である $(x, y) = (3, 2)$ が得られました！
(想定したビットの範囲に解が存在しない場合もあるので、検算は必要です)

$D=2$ を解くにはオーバースペックですが、大きな $D$ に対しては量子アニーリングはもっと優位になる？？

と言ったものの、この手法には「QUBOの係数の数字が大きくなる」という難点があります。
$x$ や $y$ が $n$ ビットの数だとすると、(4乗しているので)$4n$ ビット程度の数値を係数として扱える必要があります。64ビットの係数を扱えるとしても、$n=16$ が限界です。$x$ や $y$ が 16ビットだとすると、65535までしか表せません(実際には $D$ や$ L$ をかけているので、限界値はもっと低いでしょう)。
Wikipediaの「ペル方程式」を見ると、$D=53$ の場合は求められないことになります。

Wikipedia: ペル方程式

根本的な問題が。。。

試しに、$D=3$ のケースも実行してみます。Pell方程式は $x^2 - 3y^2 = 1$ なので、解 $(x, y) = (2, 1)$ が得られることを期待しましょう。
ステップ5のところで $D=3$ とすると、次のようになります。

\begin{align}
(x^2 - Dy^2 -1)^2 =& 
- 160 x_{0} x_{1} - 6 x_{0} y_{0} - 24 x_{0} y_{1} + 480 x_{0} z_{0} - 25 x_{0} - 24 x_{1} y_{0} - 96 x_{1} y_{1} + 624 x_{1} z_{0} \\
& - 88 x_{1} + 600 y_{0} y_{1} - 240 y_{0} z_{0} + 480 y_{0} z_{1} + 15 y_{0} - 312 y_{1} z_{1} + 168 y_{1} + 672 z_{0} z_{1} \\
& - 864 z_{0} - 840 z_{1} + 1
\end{align}

ステップ6のプログラムに適用した結果はこちら。

{x0: 1, x1: 1, y0: 0, y1: 1, z0: 0, z1: 1} 	 E =  -1376.0

$x, y$ を求めると、

x = x_0 + 2x_1 = 1 + 2\cdot1 = 3\\
y = y_0 + 2y_1 = 0 + 2\cdot1 = 2

なので、 $(x, y) = (2, 1)$ にならない。。。
う～ん、これでは $D=2$ が運良く当たっただけの可能性もあるし、信頼性に欠けますね。私の計算(プログラム)が誤っている可能性もありますが、解決しませんでした。
以下にプログラムを公開しているので、分かる方、せひ教えてください！

pell.ipynb

今後の期待

実際、量子アニーリングを行おうとすると、次数の大きい方程式が頻出すると思うんですよね。
次数を減らして2次以下にする必要がありますが、Pell方程式程度でもこんなに大変な作業です。
こういうのを簡単に行うライブラリが求められますね😁
(自分で作るフラグではありません)