はじめに
こんにちは。量子アルゴリズム開発のスタートアップであるQunaSysで量子情報エンジニア兼エバンジェリストをしている井辺と申します。
今回は量子コンピュータAdvent Calendar 2019 の10日目として、本年9~10月に開催された量子アルゴリズムのハッカソンであるIBM Quantum Challenge に社内外のメンバーで構成される3名のチームで参加し、結果準優勝したのでその記録を書きます。
ハッカソンで上位に入るために工夫した点や、Groverのアルゴリズムに対する理解が深まったのでそれらの点について共有できれば幸いです。
なお、最終的なコードはこちらで公開されています。
$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$
ハッカソン課題：自治体のコンビニ出店プランを提示せよ¶
都内某市に11の区域があり、そこに4社（A~D社）のコンビニが出店する配置を考えなさいという問題です。
下図中の11個のノード（点）が各区域に対応しており、色付きのA,B,C,Dで表示された区域は対応するコンビニがすでに出店しているため、残りの7つの区域（灰色の0,...,6で示されたノード）におけるコンビニの出店パターンを考えます。
コンビニの配置の方法には条件があり、隣接する区域（図中のエッジ（線）で繋がった区域）には異なるコンビニが出店する必要があります。たとえば2で示された地域にはBまたはDのコンビニのみ出店可能です。
なお、これは数学的にはグラフ彩色という古くから研究されている問題の一種のようです。
課題の詳細な説明は公式のリポジトリをご覧ください。
本ハッカソンの課題は、Groverのアルゴリズムをイテレーション回数（後述）5回で用いて、条件を満たす出店パターンを全て列挙することです。
ちなみに条件を満たす全出店パターンは以下の9通りとなります。
(図は公式の解説・講評より)
ハッカソンの順位は以下で定義される「コスト」の小ささで決まります。
{\rm cost} = 1 {\rm qubit}ゲート数 + 10 \times {\rm CNOT}ゲート数
このコストはGroverのアルゴリズムを動かすために用いた量子ゲートの多さに対応しています。量子ゲートが多いほど計算時間が長くなり、かつ量子コンピュータ実機での実行が困難になる（i.e., エラーの影響が大きくなる）傾向にあります。本ハッカソンでは実機ではなくシミュレータ（量子回路の計算を再現するための古典コンピュータ）で計算を実行しますが、将来実機で実行する上でコストは非常に重要な指標です。
なお、CNOTゲート数にのみ10が掛かっているのは、実機での実行が1qubitゲートに比べそれだけ困難であることを反映しています（難しさが10倍かどうかは議論の余地あり）。
また、本ハッカソンではIBM謹製のシミュレータであるibmq_qasm_simulator(32qubit)を用いるルールのため、最大32qubitという縛りのもとでアルゴリズムを組むことに注意が必要です。
Groverのアルゴリズム
Groverのアルゴリズムは、数多くの選択肢から正答を選び出す作業を高速化できるアルゴリズムです。
N個のダンボール箱の中に1つだけリンゴが入っているとします。
愚直にリンゴを探そうとすると、最大でN-1回箱を開けるという操作が必要になります。すなわち、箱を開ける回数は$N$に比例します。
一方でGroverのアルゴリズムでは、$\sqrt{N}$に比例する回数で正答が得られます（答えが$k$個ある場合は$\sqrt{N/k}$に比例）。
ダンボール箱とリンゴは比喩ですが、本ハッカソンの問題ではコンビニの$N=4^7=16384$種類の全配置パターンが各ダンボール箱に対応しており、条件を満たす配置パターンがリンゴの入っている箱に対応しています。
Groverのアルゴリズムでは量子的な重ね合わせ状態を用い、かつ正答が得られる確率を巧妙に増幅していくことで、計算回数（i.e., ダンボール箱を開ける回数）を$N\to\sqrt{N}$のオーダーで削減することができます。
本ハッカソンではこの計算回数（Grover iterationと呼ばれる）は5回に定められています。今の場合は5回程度でも十分正答が得られる事が後ほどわかります。
Groverのアルゴリズムに関するより詳しい説明はQuantum Challenge Week2やQuantum Native Dojo 8.2をご覧ください。
Groverのアルゴリズムを用いたコンビニ出店問題解決の手順
本ハッカソンの我々のチームの解法について、以下で構成されるGroverのアルゴリズムの手順に従って説明していきます。
- 入力状態生成
- オラクル
- ディフュージョン
- 測定
入力状態生成
Groverのアルゴリズムでは、まず考えられる全配置パターンを量子状態に1:1で対応させ、そららの重ね合わせ状態を作る必要があります。
つまり、入力状態として、$N=4^7=16384$個の項をもつ
\ket{{\rm Input}} = \frac{1}{\sqrt{N}}\ket{{\rm AAAAAAA}} + \frac{1}{\sqrt{N}}\ket{{\rm AAAAAAB}} + \cdots + \frac{1}{\sqrt{N}}\ket{{\rm DDDDDDD}}
のような状態を作成します。ここで $\ket{{\rm AAAAAAA}}$は未出店の全7区域にコンビニAが出店する配置に対応（ほかも同様）しています。
（これに後述のGrover iterationを適用することで、正解の配置に対応する9個の状態（例：$\ket{{\rm BADCDAB}}$）のみ係数が増幅していきます。）
各ノードの$\ket{{\rm A}}, \ket{{\rm B}}$ といった情報をqubitで表すためには、qubitを2つ使って
\ket{{\rm A}} = \ket{00} \\
\ket{{\rm B}} = \ket{01} \\
\ket{{\rm C}} = \ket{10} \\
\ket{{\rm D}} = \ket{11}
と対応付ければOKです。したがって、$2*7=14$qubitでコンビニの配置の全パターンを網羅した入力状態を作成できます。
一方、冒頭の図を見ると上記の表現は少し冗長であることに気づきます。
例えばノード2はコンビニBorDの2択なので、$\ket{{\rm B}} = \ket{0}, \ket{{\rm D}} = \ket{1}$と対応付ければ qubitは1つで十分です。
さらに、ノード0をDと仮定するとノード6で矛盾が生じることがすぐに分かるので、ノード0はBorCの2択です。同様に、ノード1もAorDに制限されます。
以上より、ノード0, 1, 2には1qubitずつを割り振れば十分であるため、入力状態は$14-3=11$ qubitで表現できる事がわかります。
これによって考える配置パターン数が$2^{11} = 2048$個に減り、探索する空間が狭まりました。
$\ket{0}$に初期化された11個のqubit全てにHadamardゲートをかけ、重ね合わせ状態を生成します。
\begin{eqnarray}
\ket{{\rm Input}} &=& H^{\otimes 11} \ket{00000000000} \\
&=& \frac{1}{\sqrt{2^{11}}}(\ket{0} + \ket{1}) \otimes \cdots \otimes (\ket{0} + \ket{1})\\
&=& \frac{1}{\sqrt{2^{11}}}\sum_{x:全配置パターン} \ket{x}
\end{eqnarray}
これで入力状態が完成です。1
オラクル
次にオラクルを構成します。
オラクル$U$とは「求める量子状態に対してのみ$-1$をかける量子回路」のことです。
すなわち、求める量子状態を$\ket{w}$とすると、
\begin{eqnarray}
U\ket{x} =\left\{ \begin{array}{ll}
-\ket{x} & (x=w) \\
\ket{x} & (x\neq w) \\
\end{array} \right.
\end{eqnarray}
と作用します。
このように量子回路（ユニタリ行列）としてオラクルを構成しておくことで、先程作成した2048個の重ね合わせ状態のそれぞれが条件を満たしているかどうかを一度に判定できます。
私は当初「オラクルが分かっているなら答え分かってるのと同じでは」と思っていましたが、実はそうではありません。
オラクルは単に答えが合っているかをチェックする「判定機」のようなもので、これを作るために答えの各パターンを全て知っていないといけないというわけではありません。
まず小さい部分ごとの条件（コンビニ問題の場合、例えばノード0と1に異なるコンビニが出店しているか）が満たされているかをチェックする「小さいオラクル」を作ります。それらを統合することで、条件が全て満たされている場合のみ「正解」という判断を下す回路ができます。これがGroverのアルゴリズムを動かすためのオラクルです。
ちょっと何を言っているのか伝わりにくいと思うので、具体的に我々のオラクルの構成法を例に取り説明します。
本ハッカソン問題はオラクルの設計に最も自由度があり、各チームでコストに差が付きやすい重要な箇所だと思います。
我々のチームのオラクル設計戦略
まず、オラクル構成法として愚直には以下の手順が考えられます。
- 各エッジで繋がっているノードどうしが異なるかを判定
- 判定結果を各エッジごとに1qubitに格納し（条件が満たされていれば$\ket{1}$を格納）
- そらら全ての
ANDを取り、結果を更に別のqubit (target)に格納
しかし、これではqubit数の上限である32qubitを超えてしまいます。
冒頭のグラフ上で考慮すべきエッジは（自明なC,D間を除いて）21本あり、ノード用の11qubitとtarget用の1qubitで既に33qubitとなってしまうためです。
他にも途中の計算で用いる補助qubit（アンシラと呼ばれる）に用いるqubitも必要なため、オラクルに用いるqubit数を何とかして削減する必要があります。
我々の実装では冒頭のグラフのエッジを以下の7つのグループ（それぞれを「パート」と呼ぶ）に分解し、それぞれに別々の条件判定回路を作成、集計を行うというアプローチを取りました。こうすることで、前述の愚直なオラクル実装法で21qubit(i.e., エッジの数）だけ必要だった部分を8qubit(注：パート4のみ判定用qubitが2個必要)に削減できます。
この実装法は問題に特化したやり方で汎用性は低いですが、このようなパート分けを自動で行うアルゴリズムを考えるのも面白いかもしれません。
例えばパート2（ノード0, 1, 3で構成される三角形）の判定回路はこのように実装しました。
def part2(circuit, q0, q1, q2, q3, t0):
"""
for nodes 0, 1, 3
:param circuit: QuantumCircuit
:param q0: node 0 (|0> = B, |1> = C)
:param q1: node 1 (|0> = A, |1> = D)
:param q2: node 3
:param q3: node 3
:param t0: target
:return:
"""
circuit.cx(q1, q0)
circuit.cx(q2, q3)
circuit.ccx(q0, q3, q2)
circuit.cx(q2, q3)
# target であるt0に判定結果を格納
circuit.cx(q1, t0)
circuit.cx(q3, t0)
# 逆操作（uncomputation）を行いtarget以外のqubitを現状復帰
circuit.cx(q2, q3)
circuit.ccx(q0, q3, q2)
circuit.cx(q2, q3)
circuit.cx(q1, q0)
回路図：
各ノード0, 1, 3に異なるコンビニが割り当てられている場合にのみ、target qubit（図中のq0_4）に$\ket{1}$が格納されることが確かめられるはずです。
重要なことは、各パートの実装の内部で逆操作によってノード用のqubitを原状復帰（uncomputation）していることです。
各ノードは一般に複数のパートに共有されているため、パートの条件判定でノードの状態が変化してしまうと正しい結果が得られなくなるためです。
オラクルの実装で工夫したのは、上記の実装で通常のToffoliゲート circuit.ccx()の代わりによりコストが低い「相対位相つきToffoliゲート」circuit.rccx()を最終的に用いた点です（公式のドキュメント）。Toffoliと比べCNOTの数が半分で済みます。コストが低い代償に、量子状態の位相が異なる部分が出てくるという欠点がありますが、最終的に逆操作をしてqubitの状態をもとに戻すので、ここでは通常のToffoliゲートと同様に用いる事ができます。
全体のオラクルは各パートの8個のtarget qubitをcontrolとするmct（multi-control Toffoli）ゲートを用いることで、それらが全て$\ket{1}$の場合（i.e., 全パートについて条件判定をクリアした状態）にのみ位相反転を行うようにします。
ここでもuncomputationを行い、各パーツのtarget用qubitを含む全てのqubitの状態をリセットしておきます。
ところで、uncomputationを各パート中と全体のオラクルの実装の両方の箇所で行っていることに気づきます。
詳細な説明は省きますが、我々の実装ではパートの判定の順番を工夫することで、余分なuncomputationを極限まで減らし、コストを削減しました。
ハッカソンの終盤ではこのような冗長な操作がないかのチェック作業と修正に主に時間を割きました。
ディフュージョン
さて、オラクルを作用させただけでは、欲しい状態の係数の符号が反転するだけです。
最終的に正しい解に収束する確率を高めるためには、オラクルに加えここで述べるディフュージョンという操作が必要です。
ディフュージョンでは入力の重ね合わせ状態に関する反転を行います。
上記のオラクルとこのディフュージョンを交互に繰り返す（Grover iteration）ことで、正しい解に対応する状態の振幅が大きくなっていきます。
図的な説明は Quantum Native Dojo 8.2 にあるのでご参照下さい。
実装で工夫した点は、オラクルの節でも出てきたmct（multi-control Toffoli）ゲートで、最も低コストなモードを指定することです(mctゲートの実装はこちら)。
mctゲートには4種類の実装があり、modeという引数で指定できます。この内'basic'は最も多くの補助qubit（アンシラ）を用いますが、同時に最も少ないコストで動作します。
多数のアンシラが必要なため、使えそうなqubitを探して用いました。
本ハッカソン問題を解き始めた当初はより使用qubit数の少ない（i.e., コストが大きい）モードを指定したり、rcccx(位相つきcccxゲート)を用いる等していましたが、色々試すうちにmctゲートの'basic'モードが最もコストが低いことが分かり、これを用いました。
測定
上記のGrover iterationを経て、得られた状態を基底測定します。
実はここで若干苦戦しました。
問題文ではおそらく各ノードに2つずつのqubitを割り当てることを前提として
古典レジスタc[2n], c[2n+1]にn番目の区域に出店するコンビニをマッピングしてください
という指示が出ています。
一方で我々の実装ではqubit数を削減するために一部のnodeを1qubitで表しており、測定によって11個の古典レジスタに結果が格納されます。これを問題文で求められている14個の古典レジスタにmapするのに少し手間取りました。
qubitの測定以外で古典レジスタを操作する方法がQiskit上で見つけられなかったためです。
結局量子回路上で追加のゲートを用いて対処しましたが、古典レジスタを直接操作できる機能も欲しいと思いました。
qubitの割当一覧
最終的に我々の解法では27個のqubitを用いました。各qubitの用途一覧は以下の通りです。
| qubit | node |
|---|---|
| 0 | 0 (one qubit to represent B/C) |
| 1 | 1 (one qubit to represent A/D) |
| 2 | 2 (one qubit to represent B/D) |
| 3 | 3 (the right digit of node 3 e.g., $\ket{1}$ if $\ket{{\rm node3}}=\ket{B}=\ket{01}$ and so forth) |
| 4 | 3 (the left digit of node 3 e.g., $\ket{0}$ if $\ket{{\rm node3}}=\ket{B}=\ket{01}$ and so forth) |
| 5 | 4 |
| 6 | 4 |
| 7 | 5 |
| 8 | 5 |
| 9 | 6 |
| 10 | 6 |
| 11 | target for part0 and diffusion MCT ancilla 0 |
| 12 | target for part1 and diffusion MCT ancilla 1 |
| 13 | target for part2 and diffusion MCT ancilla 2 |
| 14 | target for part3 and diffusion MCT ancilla 3 |
| 15 | target for part4 and diffusion MCT ancilla 4 |
| 16 | target for part4 and diffusion MCT ancilla 5 |
| 17 | target for part5 and diffusion MCT ancilla 6 |
| 18 | target for part6 and diffusion MCT ancilla 7 |
| 19 | part4, ancillary |
| 20 | part5, ancillary |
| 21 | part6, ancillary |
| 22 | phase MCT ancilla0 |
| 23 | phase MCT ancilla1 |
| 24 | phase MCT ancilla2 |
| 25 | phase MCT ancilla3 |
| 26 | phase MCT ancilla4 |
計算結果とコスト
計算結果は下図のとおりです。
ここでGrover iteration数とshot数(量子回路の実行・測定という一連の流れの繰り返し回数)は指示どおりそれぞれ5回、8000回としています。また、ここでは見やすさのため測定によって得られた回数の上位30パターンのみ図示しています。
正解の状態に対応する9つの測定結果の頻度が、他と比べ有意に高いことが見て取れます。
最終的なコストは17053となりました。
締め切り前ラスト数時間はランキングがかなり変動していましたが、締め切り時点で我々のチームは5位に着けました。
その後judgeによる審査を経て、1週間後に最終的な結果が発表されます。
ハッカソン最終結果
結果は2位でした。
(cf. Leader Board のQunaVillageチーム)
締め切り時点では我々のチームは5位でしたが、上位4チームの内3チームが失格の扱いとなり、我々が繰り上がっていました。
後でそれらの内一つのチームの方に聞いたところ、5回のiterationの内1,2回目と3,4,5回目とで異なる実装にする、といった手法でゲート数を削減していたそうです。最終的な答えは合っており、なかなか際どいラインではありますが、今回は失格と見なされたようです。
問題文には
week2,3で使ったグローバーのアルゴリズムをiteration回数＝5で用いて、
とあるのみなので「全く同じ回路を繰り返す必要は無い」とも取れなくはないですが、これを許してしまうと極端ですが例えば「1,2,3回目は通常のGrover iteration, 4,5回目のiterationではidentityをかける」といった解法もまかり通ってしまうので、失格は妥当かもしれないと思いました。なお当然ながら他の失格のパターンもあったのかもしれません。
一方で1位に輝いたWhit3zチームは、同じqubitに連続して掛かっている1qubitゲートをまとめる等、我々よりも細部まで行き届いた最適化を行っていたようです。
（コードはこちら）
ハッカソン全体のまとめと感想
本ハッカソンでは、条件を満たすコンビニの出店パターンをGroverのアルゴリズムを用いて探索し、更に計算に用いる量子ゲート数（コスト）を可能な限り削減するという問題に挑戦しました。
普段は専らVQEを始めとするNISQアルゴリズムの実装を行っているので、Groverのアルゴリズムのようなlong-termのアルゴリズムで問題を解くというのは、新鮮で得難い経験だったと思います。
Groverのアルゴリズムに関して、具体的には以下の点で特に勉強になりました。
- 正答が分かっていなくてもオラクルを作れることが分かったこと（実はこのハッカソンに参加するまでは「オラクルが分かっているなら答えが分かってるのと同じでは？」と思っていました）
- Qubitの再利用。uncomputationを行い「ゴミ」を残さず可逆計算のみで欲しい結果を得るための手順
- Grover iterationを繰り返していくことで、欲しい解に対応する状態の確率（振幅）が徐々に上がっていく様子を実感できたこと
また、本ハッカソンで初めてQiskitに触れましたが、QiskitはPythonでシンプルに記述できるうえ、回路のプロット機能やチュートリアルも充実しているため、入門者にも非常にとっつきやすい言語だと感じました。
クラウドで公開されている実機やシミュレータで、コードをそのまま実行できる点も素晴らしいです。
追加で欲しい機能を一つ挙げるとするならば、「測定」で触れた通り、古典レジスタを直接操作する機能があればより便利かもしれません。
なお本ハッカソン準優勝で招待されたQiskit Camp Asiaでは、上記のコンビニ問題とほぼ同様の考え方でGroverのアルゴリズムを応用し、優秀な学生の方々やIBMのコーチの方とチームを組み、有名な数字パズルである数独を解きました（Github repoはこちら）。
本記事の続編として後日まとめられればと思っています。