Jij Tech Blog

Jij inc.の開発日記です

断熱量子線形回帰

序章: 線形回帰

 N>dとして、 N組のデータ (Y_i, x_{i2}, x_{i3}, \dots , x_{id}) \ (i=1, 2, \dots, N)に対する、線形回帰モデルは

 \displaystyle{
Y_i = \omega_1 + \omega_2 x_{i2} + \omega_3 x_{i3} + \cdots + \omega_d x_{id} + \epsilon_i \ (i = 1, 2, \dots, N)
}

となります。ここで

 \displaystyle{
{\bf Y} = \left( \begin{array}{c}
Y_1 \\
Y_2 \\
\vdots \\
Y_N
\end{array} \right), \ 
{\boldsymbol \omega} = \left( \begin{array}{c}
\omega_1 \\
\omega_2 \\
\vdots \\
\omega_d
\end{array} \right), \ 
{\boldsymbol \epsilon} = \left( \begin{array}{c}
\epsilon_1 \\
\epsilon_2 \\
\vdots \\
\epsilon_N
\end{array} \right)
}

と、以下の N \times dの行列を定義します。

 \displaystyle{
X = ({\bf 1}, {\bf x}_2, \dots, {\bf x}_d) 
= \left( \begin{array}{cccc}
1 & x_{12} & \cdots & x_{1d} \\
1 & x_{22} & \cdots & x_{2d} \\
\vdots & \vdots & \ddots & \vdots \\
1 & x_{N2} & \cdots & x_{Nd}
\end{array} \right)
}

これらを用いて線形回帰モデルの式を表現すると

 \displaystyle{
{\bf Y} = X {\boldsymbol \omega} + {\boldsymbol \epsilon}
}

となります。ここから誤差の2乗和

 \displaystyle{
\| {\boldsymbol \epsilon} \|^2 
= \| {\bf Y}- X {\boldsymbol \omega} \|^2 
= {\bf Y}^T {\bf Y} - 2 {\boldsymbol \omega}^T X^T {\bf Y} + {\boldsymbol \omega}^T X^T X {\boldsymbol \omega}
}

の最小値を与える {\boldsymbol \omega}を推定することを最小二乗法と呼びます。最適化において {\bf Y}^T {\bf Y}は定数項なので、これを無視した

 \displaystyle{
\| {\boldsymbol \epsilon} \|^2 
= - 2 {\boldsymbol \omega}^T X^T {\bf Y} + {\boldsymbol \omega}^T X^T X {\boldsymbol \omega} \tag{1}
}

を最小にすることを目指します。

QUBO化

精度ベクトル

(1)式をQUBO化します。そのために精度ベクトル {\bf P} = (p_1, p_2, \dots, p_K)^Tを導入します。ここで {\bf P}の成分は2の冪乗で表現し、プラス・マイナス両方取りうるとします。ただし、 {\bf P}の成分は大きさ順に並べられていなければなりません。

例として次のような精度ベクトルが考えられます:  {\bf P} = (-2, -1, -1/2, 1/2, 1, 2)

バイナリ変数

さらに d \times Kのバイナリ変数の行列 \hat{\omega}を導入し、これと精度ベクトルを用いて \omegaを表現します。

 \displaystyle{
\omega_i = \sum_{k=1}^K p_k \hat{\omega}_{ik}
}

例として {\bf P} = (-1, -1/2, 1/2, 1)を用意すると、 \omega = (-3/2, -1, -1/2, 0, 1/2, 1, 3/2)を表現することができます。

 \mathcal{P} = {\bf I}_d \otimes {\bf P}^Tを用いることで、上式を行列形式に書き換えましょう。

 \displaystyle{
\omega = \mathcal{P} \hat{\omega}
}

QUBO化

以上より(1)式は

 \displaystyle{
\min \| {\boldsymbol \epsilon} \|^2 
= - 2 {\hat \omega}^T {\mathcal P}^T X^T {\bf Y} + {\hat \omega}^T {\mathcal P}^T X^T X {\mathcal P} {\hat \omega}
}

となります。これをD-Waveで解くことを目指します。この式は

 \displaystyle{
\min_z (z^T A z + z^T b)
}

と同じ形をしており、これはNP困難であることが知られています。

数値実験

比較対象

D-Wave 2000Qマシンを用いた量子アニーリングがどれほど有用であるかを検証するために、古典的なアプローチとして以下のようなものを用います。

  • 3.6 GHz 8-Core Intel i9 プロセッサ
  • 64GB 2.666 MHz DDR4 メモリ
  • 用いた数値計算ライブラリ: Scikit-learn

数値実験結果

回帰エラー

回帰エラーを元に2つの比較結果を見てみましょう。

Experimental Runs Where Scikit-learn Error D-Wave Error
D-Wave fit the data (80% runs) 5.0261 5.0362
D-Wave fit the data (20% runs) 4.8188 15.0657
Overall 4.9846 7.0421

これを見るとD-Waveでの実行の方がエラーが大きくなっています。これは量子アニーリング実行中に量子ビット間の接続が切れると、誤った結果が出ることが知られているD-Waveマシンのハードウェアの問題だと考えられます。

データ数Nに対するスケーラビリティ
Scalability of N comparison of Scikit-learn and D-Wave
Nに対するスケーラビリティの比較

上図は横軸をデータ点の数、縦軸を計算時間にしてグラフを書いたものです。左図はデータ点数の小さいときの比較図、そそして右図はデータ点数が大きい時の比較図です。データ点数が1600万になると、D-Waveでの計算時間がScikit-learnでのそれよりも2.8倍高速であることがわかりました。

特徴量数dに対するスケーラビリティ
Scalability of d comparison of Scikit-learn and D-Wave
dに対するスケーラビリティの比較

こちらの図は横軸を特徴量の数、縦軸を計算時間にしてグラフを書いたものです。特徴量数が32の時点で、D-Wave実行がScikit-learnでの計算時間よりも2.8倍高速であることがわかりました。

結言

古典コンピュータ上でScikit-learnによる線形回帰と、D-Wave上で同じ問題を解いた場合の比較実験を行いました。その結果、D-Waveのハード性能の問題から、精度には古典コンピュータに軍配が上がりました。一方で、データ点数が大量になるとD-Waveでの量子アニーリングの実行の方が高速であることがわかりました。

文責

中村翔、株式会社Jij
Sho K. NAKAMURA, Jij inc.
憂いの篩 -Pensieve-

※ 疑問点や間違いなどございましたら、お気軽にブログ上でご質問ください。