はじめに

量子アニーリングマシンの応用事例として、素因数分解への応用がいくつか既になされている。（例えば、[S. Jiang et al., Scientific Reports 8, 17667 (2018).]） また、2018年後半に新たに公開されたWebインターフェース D-Wave Leap においても乗算器ベースの素因数分解アプリケーションがD-Waveマシン利用のチュートリアルとして紹介されている（T-QARDの日々第12回における紹介動画）。それでは、数論の問題で素因数分解以外に困難な問題でD-Waveマシンが利用できそうな問題はあるだろうか？そのような問題として、 ディオファントス方程式 の求解に取り組んでみた結果を本記事と続く記事によって紹介したい。

ディオファントス方程式とは

整数の係数を持つ変数の整数乗の項からなる方程式であり(例. $x^2+3y^2=10,x^5=1024$)、ディオファントス方程式の名前を用いる時にはその解が有理数または整数に限定する場合が多い。特に、2変数2次のディオファントス方程式 $ax^2+by+c=0$の整数解の存在判定はNP完全問題に属することが知られている。そこで、一連の記事ではこの$ax^2+by+c=0$のタイプの整数解のD-Waveマシンによる求解を目標とする。本記事ではその練習問題として、2変数1次のディオファントス方程式$ax+by+c=0$をD-Waveマシンを用いて解いた。

問題の定式化

二値変数の定義（方程式解の二値表現）

求める整数$x, y$を次のように2進数として表現し、それぞれの桁をD-Waveマシンで求める2値変数とした。（今回は簡単のために$x, y$は正の数のみを探索する。）

$$
x=\sum_i{x_i2^i}, y=\sum_i{y_i2^i}, x_i, y_i \in \{0,1\}
$$

ただし、このままでは無限に変数が増えて計算できないため、適当な範囲で桁数を限定する必要がある。

方程式を解く方針

一般に、方程式$f(x_1, x_2, \dots, x_N)=0$の求解をQUBOソルバーで行うためには、最小化すべき目的関数として、

$$
H(x_1, x_2, \dots, x_N) = (f(x_1, x_2, \dots, x_N))^2
$$

を考えるのが素朴な戦略である。

たとえば、今回のケースでは

$$
H(x, y) = \left( ax+by+c \right)^2
$$

を考えることとなる。イメージとしては以下の図のようになる。

この図は左が元の方程式のグラフ、右がそれを2乗したグラフとなっている。この図を見ると分かるように、正負にかかわらず、0でない数は2乗すると正の値を持つことが分かる。この時、0は2乗しても0のままなので、2乗した後のグラフでは最小値をとる点と0をとる点が一致する。$ax+by+c$も同様に二乗すると$ax+by+c=0$となる点、すなわち方程式の解となる点が最小値になる。そのため、$(ax+by+c)^2$は目的関数として機能することが分かる。

であり、これを2値変数を用いるQUBOの形に直すと、

$$
H(x_0, \dots x_{N-1}, y_0, \dots, y_{N-1}) = \left( a\sum_{i}{x_i2^i}+b\sum_{i}{y_i2^i}+c \right)^2
$$

となる。実際、この式は$x_i, y_j$に関する二次項のみからなっており確かにQUBO型の目的関数となっていることがわかる。

QUBO行列の準備

D-Waveマシンに送信する情報はQUBOの一般形

$$
H(\mathbf{z})=\mathbf{z}^{\mathrm{T}} Q \mathbf{z}
$$

における行列$Q$の要素である。ここで、今回の問題における解の二値表現$x_0, \dots, x_{N-1}, y_0, \dots, y_{N-1}$を繋げたベクトルとして、$\mathbf{z} = (x_0, \dots, x_{N-1}, y_0, \dots, y_{N-1})^{\mathrm{T}}$を用いた。

さて、行列$Q$の要素を計算するためには地道に目的関数を展開すればよく、

$$
\begin{equation*}
\begin{split}
H
&= a^2\left(\sum_{i}{x_i2^i}\right)^2+b^2\left(\sum_{i}{y_i2^i}\right)^2+c^2 \\
& + 2ab\left(\sum_{i}{x_i2^i}\right)\left(\sum_{i}{y_i2^i}\right)+2ac\left(\sum_{i}{x_i2^i}\right)+2bc\left(\sum_{i}{y_i2^i}\right) \\
&= a^2\left(\sum_{i}{x_i2^{2i}}+2\sum_{i<j}{x_ix_j2^{i+j}}\right)+b^2\left(\sum_{i}{y_i2^{2i}}+2\sum_{i<j}{y_iy_j2^{i+j}}\right)+c^2 \\
& + 2ab\left(\sum_{i,j}{x_iy_j2^{i+j}}\right)+2ac\left(\sum_{i}{x_i2^i}\right)+2bc\left(\sum_{i}{y_i2^i}\right) \\
&= \sum_{i}\left(a^22^{2i}+2ac2^i\right)x_i^2+\sum_{i}\left(b^22^{2i}+2bc2^i\right)y_i^2+\sum_{i<j}{(2a^22^{i+j})x_ix_j} \\
& + \sum_{i<j}{(2b^22^{i+j})y_iy_j}+\sum_{i,j}{(2ab2^{i+j})x_iy_j}
\end{split}
\end{equation*}
$$

と変形できる。続いて、例として桁数を$N$とした時、

$$
z_i=\left\{\begin{array}{11}x_i & (i\leq N-1) \\ y_{i-N} & (i\geq N)\end{array} \right.
$$

であるから、目的関数を$z_i$で表すと、

$$
\begin{equation*}
\begin{split}
Obj
&= \sum_{i=0}^{N-1}\left(a^22^{2i}+2ac2^i\right)z_i^2+\sum_{i=0}^{N-1}\left(b^22^{2i}+2bc2^i\right)z_{i+N}^2+\sum_{i=0}^{N-1}\sum_{j>i}^{N-1}{(2a^22^{i+j})z_iz_j} \\
& + \sum_{i=0}^{N-1}\sum_{j>i}^{N-1}{(2b^22^{i+j})z_{i+N}z_{j+N}}+\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}{(2ab2^{i+j})z_iz_{j+N}}
\end{split}
\end{equation*}
$$

したがって、Qの対角成分$Q_{ii}$と非対角成分$Q_{ij}$は次のようになる。

$$
\begin{align}Q_{ii}
&= a^22^{2i}+2ac2^i\ \ \ (0\leq i \leq N-1)\\Q_{(i+N)(i+N)}
&= b^22^{2i}+2bc2^i\ \ \ (0\leq i \leq N-1)\\Q_{ij}
&= 2a^22^{i+j}\ \ \ (0\leq i<j\leq N-1)\\Q_{(i+N)(j+N)}
&= 2b^22^{i+j}\ \ \ (0\leq i<j\leq N-1)\\Q_{i(j+N)}
&= 2ab2^{i+j}\ \ \ (0\leq i \leq N-1,0\leq j\leq N-1)
\end{align}
$$

問題の実行例

jupyter notebookを用いたD-Waveマシンの使用

上の定式化に従って、jupyter notebookを介したD-Waveマシンを使って実際に解いてみた。

使用したプログラム

今回、大きく分けてD-Waveマシンへ入力する行列を生成するプログラムとD-Waveマシンから返ってきた出力を解析するプログラムの2つを使用した。

行列生成プログラム

#方程式の入力
a = 3
b = 8
c = -73

#x,yの桁数の設定
xmax=(0-c)/a #xの最大値
ymax=(0-c)/b #yの最大値

for i in range(100):
        if(2**(i-1)<xmax<=2**i):
            length = i
            break
for i in range(100):
        if(2**(i-1)<ymax<=2**i and i > length):
            length = i
            break

#行列のサイズの決定            
length2 = length * 2
Q1 = [[0] * length2 for i in range(length2)]

#行列Q1の要素の計算
for i in range(0,length):
    Q1[i][i] += (2**(2*i))*(a*a)+(2**(i+1))*a*c
    Q1[length+i][length+i] += (2**(2*i))*(b*b)+(2**(i+1))*b*c
    for j in range(i+1,length):
        Q1[i][j] += (2**(i+j+1))*(a*a)
        Q1[length+i][length+j] += (2**(i+j+1))*(b*b)
    for j in range(0,length):
        Q1[i][length+j] += (2**(i+j+1))*a*b

まず、

$$\begin{align*}x=\sum_i{x_i2^i}&&y=\sum_i{y_i2^i} \\ x_i=\{0,1\}&&y_i=\{0,1\}\end{align*}$$

と置いている都合上x,yは正の値しか取れない。それを考慮し、このプログラムではa>0,b>0,c<0に係数を限定することで解の個数を有限にしてD-Waveマシンで解きやすいプログラムにした。

xが最大値を取るのはy=0の時であり、そのとき$ax+c=0$となるためxの最大値$x_{max}=-c/a$。同様にyの最大値$y_{max}=-c/b$となる。このように、x,yが最大値で抑えられることを利用してxやyは2進数で何桁必要かを求める。$2^{i-1} \leq x_{max}<2^i$ならば全てのxはi桁あれば表現でき、yにも同様のことが言える。プログラム中盤のfor文で今の計算を行なっており、x,yのうち必要な桁数の多い方に合わせてxとyに用意する二値変数の数(length)を決めている。

先ほど決めたlength (桁数$N$)の2倍だけ変数を準備すればいいため、length*2×length*2の正方行列Q1をD-Waveマシンに送る行列として用意した。

その後、問題の定式化で説明した通りに行列の要素を生成し、行列をD-Waveマシンに送り、返ってきた出力を次のプログラムで処理した。

出力解析プログラム

#得られた解の数に応じた繰り返し
for i in range(answers):
    x=0
    y=0
#2進数から10進数へ
    for j in range(length):
        if(decode_answer[i][j]==1):
            x += decode_answer[i][j]*(2**i)
        if(decode_answer[i][length+j]==1):
            y += decode_answer[i][length+j]*(2**i)
#出力
    print 'x =',
    print x,
    print ' y =',
    print y,
#ダブルチェック
    if(a*x+b*y+c==0):
        print 'true'
    else:
        print 'false'

D-Waveマシンは何回も同じ問題を解き、その解答の中から良い解をいくつか選び出力しているため、出力される解は複数存在する。具体的には次のように行列の形で出力されている。

上図のようにiで何番目に得られた解を調べるかを指定し、jで何番目の変数を見るか指定することができる。今回のプログラムでは得られた解の個数だけプログラムを繰り返して、全ての解に対して解析を行なった。

次に、解として得られた解は2進数の形で得られるため、下式にしたがってx,yを10進数に戻している。

$$
\begin{align*}
x=\sum_{i=0}^{length-1}{x_i2^i}=\sum_{i=0}^{length-1}{z_i2^i}&&y=\sum_{i=0}^{length-1}{y_i2^i}=\sum_{i=0}^{length-1}{z_{length+i}2^i}
\end{align*}
$$

D-Waveマシンの特性上、ここで得られた解は近似解であるので必ずしも正しい解とは限らないので、元の方程式に代入して満たすかどうかをダブルチェックしている。

プログラムの実行例

実際に7x+3y=88を解かせてみると、(x,y)=(10,6),(0,29)の2通りの解が得られた。このうち、(10,6)の方は正答であるが、(0,29)は7x+3y=88を満たさないことが分かる。このように、D-Waveマシンの結果は絶対に合っているわけではないので、先ほど述べたダブルチェックが必要となる。

問題点としてはもう一つ、7x+3y=88の解は(x,y)=(10,6)以外にも(x,y)=(1,27),(4,20),(7,13)が存在しているが、それは出力されていない。これらの解が求められない理由は2つある。

• 試行回数
  D-Waveマシンは近似解を出力するため、同じ問題を解かせても毎回異なる解が出力される。一度の結果が全てというわけではなく、何回か繰り返すことで他の解が得られることがあるので解のバリエーションを得たい場合複数回D-Waveマシンに解かせる必要がある。ただし、次の理由によって試行回数を重ねても出ない解が存在する可能性がある。
• 出力結果の偏り
  Waveマシンの解には出やすいものと出にくいものが存在することがわかっている。

最後に

今回は定式化ととりあえずの試行の様子を述べた。次回は二変数二次のディオファントス方程式の求解に取り組み、また、今回と合わせて得られる近似解の振る舞いについて議論したい。

本記事の担当者

金垣 葵（編集：観山正道）