Googleの量子超越性をスパコンで300秒で計算した論文を解説してみる
Googleの量子超越性の実証に用いられた計算をスパコンで5分で達成できたニュースがSNSで話題になっていたので、解説してみます。
背景テンソルネットワーク手法の背景1次元量子状態2次元量子状態テンソルネットワーク計算のパフォーマンス依存性具体例（はじご状テンソル）a_1, a_2, a_4, a_3 の順a_2, a_5, a_8, a_3の順量子超越300秒論文の手法概要手法結果まとめ、感想参考文献
背景
2019年、Googleの開発する53qubit量子コンピュータであるSycamoreが、既存のコンピュータでは計算に途方もない時間がかかる (と思われる)計算を200秒で解いたというニュース *1 が流れ、注目を集めました。 その一方で、IBMが結果に対して「既存のコンピュータでも実現可能だ」という反論 *2 を行ったり、様々な団体が既存のコンピューターでGoogleの行ったタスクを実現した論文を発表するなど、量子デバイスと既存コンピュータのせめぎ合いが行われていました。
しかしながら「量子超越性の指標となっていたタスクがスパコンで300秒で解かれた」という事実のみが独り歩きして話題となってしまっており、具体的にどのような手法で達成したのかに関する情報があまり出てこなかったので、背景知識も含め解説してみることにしました。
テンソルネットワーク
この節ではテンソルネットワークの概要や、手法のメリットについて説明します。
手法の背景
量子計算は、量子状態に対して量子ゲートと呼ばれる演算を作用させることにより様々な演算を行う手法です。以下の図のように量子回路の形式で表されることが多いです。
この計算過程を既存のコンピュータ (量子コンピュータと対比して「古典コンピュータ」と呼ばれます) で再現しようとすると、例えばn個の量子ビットを扱う場合、次のような量子状態を表すベクトルを取り扱う必要があります
は量子状態を表現するベクトルの係数に相当しますが、各々のは の値を取るため、上のの取りうる値は通りにのぼります。この通りの数全てを古典コンピュータで管理しないといけないため、 が比較的小さい値でもすぐに古典コンピュータでの計算が破綻してしまいます。
この状況に対して古典コンピュータも黙って見ているわけもなく、なんとかして量子状態を古典コンピュータで取り扱うための手法が今に至るまで研究されてきました。その研究成果の一つとして、愚直に古典コンピュータで通りの数全てを管理せずとも係数であるを
のようにテンソル積に分割することで見通しがよくなったり計算コストの削減が可能になることが判明し、今日ではテンソルネットワークという手法としてまとめられています。こちらにもテンソルネットワークの説明が記述してあります。
以下では実際にこの記述を用いて、具体的にどのように量子回路を表現するかを見てみましょう。
1次元量子状態
まず、次のように一次元的に量子ビットが並んだ状態を考えてみます。この回路では隣り合う量子ビット間にのみCNOTゲートなどの2量子ビット演算が作用するとします。
この時、テンソルネットワークの記法を用いて量子状態の係数を次のように表現することができます。
上の状態を視覚的にわかりやすくするため、次のように図示します。
ここでの記法のルールとして、
- 四角で表される点が各テンソルを表す。
- 各テンソルから生えている足が、テンソルのindex (やなど)を表す。
- 同じ足が結合している場合、そのindexに対して内積を計算する。
とします。indexに対して内積を計算すると式からそのindexが消えてしまうため、図の中でテンソルの足が結合されずに浮いたままになっているindexが最後に残るindexです。上の例ではにあたります。
の和の範囲を表すの値を変えると、表現できる量子状態の複雑さを制御することができます。例えばが1だと
のように、が行列ではなくただのスカラー値になってしまうため、量子ビット間の相関を表現できません。の大きさを変えていくと異なる量子ビット間との相関を表すことができるようになります。初期状態は量子ビットに相関がない状況が多いため (全て など)、としてしまってよいでしょう。
さて、この状態に、 のみに作用する量子ゲート
をかけてみましょう。次のようになります。
ここで、係数だけに着目すると次のように変換されることがわかります。
これを見ると、 以外のテンソルに対する影響が一切ないことがわかります。このため、 に関わるテンソルを
と置き直すだけで済むことがわかります。量子状態の係数は
となることがわかります。
次に、 に作用する2量子ビット間ゲート
に関しても同様の作業を行うと、量子状態の係数に関して次のように整理できます。
一見すると1, 2量子ビット間にまたがっており、扱いが厄介に見えますが、特異値分解を用いるとこの演算子はそれぞれの量子ビットごとに分解することができ、
とすることができます。改めて量子状態の係数に代入すると、
となり、
と書き直すことにより、量子ゲートを作用させた後の量子状態は次のようになります。
ここで、とが常に一緒にあらわれていることに気づくと、 をまとめてと書き直してしまうことにより、元の記法に戻すことができます。
の和の範囲がになったのは、の自由度も考慮にいれる必要があるためです。異なる量子ビット間の演算を行ったことで相関が発生し複雑な量子状態になったため、の総和の範囲が拡大したと解釈することもできます。
以上をまとめると、テンソルネットワークは量子ゲートを作用させると次のように変化することがわかります。
- 1量子ビットゲートの場合、 の和の範囲である の値は変わりません。
- 2量子ビットゲートの場合、対応するindex (先程の例では)の和の範囲であるの値が2倍になります。
また、2量子ビット間ゲートを作用させている様子を図示すると次のようになります。太い線はの値が2倍になっている様子を指し、1つめと2つめの量子ビット間のindex () の総和の範囲が拡大していることがわかります。
特に、depthが8ごとに2量子ビット間ゲートを作用させるようなランダム量子回路であれば、
と見積もることができます。*5
2次元量子状態
今回は1次元的に配置された量子ビットを扱いましたが、この議論は一般の量子ビットの並べ方に拡張することができます。
たとえば次の図のように2次元的に並ぶ場合は
のようにかけます。1量子ビット、2量子ビット間の量子ゲートの作用に関しては上記と同じです。
テンソルネットワーク計算のパフォーマンス依存性
具体例（はじご状テンソル）
この計算の際、index の和をどの順序で取るかで計算コストが大幅に変わってしまう問題があります。 極端な例としてつぎのようなはしご状のテンソルを考えます。
式を書くと次のようになります。 (以後省略します)
このテンソルを2通りの手法で和をとってみましょう。
a_1, a_2, a_4, a_3 の順
まずは左から右に進むように和を取ってみましょう。次のようにテンソルが変化することがわかります。
この時、計算の過程で生成されるテンソルからは常に足が2本しか生えていないことがわかります。そのため、テンソルのサイズが肥大することがありません。
a_2, a_5, a_8, a_3の順
次に、上から下に進むように和を取ってみます。次のようになります。
この時、計算の過程で生成されるテンソルから足が3本、４本と生えていることがわかります。これは計算の過程でテンソルのサイズが指数的に肥大することを表しています。*6
このように、テンソルネットワーク状態を計算するに当たり、テンソルのサイズを肥大させないようにどの和を優先して取るか
がパフォーマンス向上において非常に重要になってきます。
今回のようにはしご状の模型であればとるべき和の順序はパッと見で判断がつきますが、任意の形状だと最適な和の評価順序は見当がつきません。また、 の2次元格子模型の場合、どうしても縦か横方向に和を取る必要があり、テンソルサイズのでの指数増加は免れません。Googleの量子超越性に用いられた模型の古典コンピュータでのシミュレートが難しい原因はここにあります。
ちなみに、テンソルネットワークのどの和を優先して取れば計算コストを最適にできるかという問題はそれ自体難しいタスクであり、様々なヒューリスティクスを用いて解いた研究 *7 も存在します*8。今回の論文もツールの一つとしてこの研究から生まれているCoTenGra *9 とよばれるソフトウェアを使っています。
量子超越300秒論文の手法
ようやく本題に入りますが、これらの背景を踏まえ、量子超越に用いたタスクを300秒で計算できた論文 Y. Liu, et al. (2021) を紹介します。 結論から言ってしまえばスパコンの数と性能の力で殴り倒したという話なのですが、ただスパコンを用いたところで適切な並列化等を工夫しなければ性能が全く活かせないので、彼らがどのような工夫を用いたかについてまとめます。
概要
本研究ではSunway TaihuLightと呼ばれる中国のスパコン上で大規模なクラスターを使用して計算を行っており、各CPUのアーキテクチャ図が以下のようになります。
詳細は割愛しますが、1つのCPUの中に6つのCore Group, さらにそれぞれのCore Groupのなかに64個のスレッドが動く構成になっています。本研究では107,520個のCPUを用いて量子回路のシミュレートの計算を行ったようです。
このように多数のプロセスを並列して量子回路のシミュレートを行うにあたり、以下の工夫を行う必要があります。
- テンソルネットワーク上のindexの総和を取るタスクを、上記のCPUを無駄なく並列に動作させて計算できるように分割する必要がある。
- 格子状にならんだ量子ビットだけでなく、任意のトポロジーに対しても効率的なindexの総和の取り方を考える必要がある。
手法
この論文では、10x10 の100qubitの格子模型を(1+40+1)depth *10、つまり42回全てのqubitに対して演算を施した回路のシミュレートと、GoogleのSycamoreデバイスで実行された回路のシミュレートの2種類を取り扱っており、それぞれに対して異なるアプローチを取っているため、その2種類に対してなされたテクニックの紹介を行います。
slicing
先程も述べたように、テンソルネットワークの計算において、ネックとなるのはテンソルのindexの総和を計算する際に中間生成物として作られるテンソルのサイズが肥大してしまうことでした。しかし、もし一部のindexの総和を計算しなくて良いのであれば計算コストが大幅に緩和されるケースも存在します。
上で扱ったはしご状のテンソルネットワークを再び見てみましょう。
もしこのindexのうち、例えばの総和をとらなくて良いのであれば、グラフは次の形状になります。
式は次の通りです。
に対して和が取られていないことに注意してください。
この場合、 の順で上から下にindexの総和を計算しても、先程とは違いネットワークがの箇所で切断されているためにテンソルの足は2つ以上に増えないことがわかります。
つまり、 の取りうる値を固定してしまい、各々の値に対してテンソルネットワークのindexの総和を並列に計算し、後で計算したあとの値を全て足せば、元のグラフの形状も復元できますし、テンソルネットワークの計算困難性も緩和でき、さらにスパコンの並列性の恩恵にあずかることもできます。この技法をslicingと呼びます。 要するに空間計算量の困難を時間計算量に押し付けて、それをスパコンの並列計算でゴリ押ししているわけです。
この論文で扱っている10x10 の100qubitの格子模型では次の青色の波線の部分でslicingを行っています。
切断したindexの数 (つまり和を取らないindexの数)は6つで、そのindexの総和の範囲である
がぐらいだとすると、depth が40であることからざっくり の並列計算がこの時点で必要になるので、相当の力技であることがわかります。
ちなみに、この箇所でslicingを行うと、次のようにテンソルのindexの和を取ることができます。四角は全て生成されたテンソルを表します。
このような順序で和を取ると、テンソルの足がつねに5,６本ほどで抑えられていることがわかります。このため、直接格子模型を扱うのに比べて、約で3000万倍程度のメモリ消費量を削減することが可能となります。ちなみにこの各々のテンソルの和の計算も並列化されています。
ヒューリスティクスによる和の順序の最適化
格子状の模型であれば上記の手法で良いのですが、扱う模型が一般の場合に対しては非効率となってしまう可能性があります。本研究ではその際にCoTenGra *12 とよばれるソフトウェアを用いて、適切なslicingや和の順序を模型の形状に応じて計算させています。
その結果が以下の図にあります。
図中の数字の4, 5がそれぞれ格子状模型でslicingを適用した場合とCoTenGraを用いた場合を表し、三角のプロットが格子状模型でのシミュレート結果、丸のプロットがGoogleのSycamore模型でのシミュレート結果になります。
格子状の模型ではCoTenGraを用いても対して性能が変わらないのに比べ、GoogleのSycamore模型では大幅に性能が改善しているのが見て取れます。
その他の工夫
用いているアルゴリズムとしてはほとんど上記で尽きているのですが、以下にその他の工夫を箇条書きで記します。
- テンソルのindexの和計算の際に、メモリアクセスが最適化されるようにDMA (Direct Memory Access)で計算すべき要素を抜き出し、その後にテンソル計算を行う。
- CoTenGraを用いた場合はテンソルの形状に差が出てしまうため、こちらも別途最適化するための手法 *13 を用いてメモリ最適化がかかりやすい形に変換を行っている。
- 単精度浮動小数 (32bit float変数)だと計算速度にボトルネックが生じるため、誤差評価を行い問題がないことを確認した上で 半精度浮動小数 (16bit float)変数を混ぜた計算を行い、FLOPSを約4倍 (1.2Eflopsから4.4Eflopsまで上昇)まで上げている。
結果
これらの数多くの工夫を行い、次の図のようにランダム量子回路で出力されるべき分布であるPorter-Thomas 分布を再現できているという図がこちらになり、この分布を300秒で計算することができたというのが本論文の主張です。
まとめ、感想
本論文の結果のキーとなるポイントは、
- slicingで空間計算量の困難性を時間計算量に押し付けて、スパコンで並列に計算する手法を提案
- 効率的なテンソルの和をとる順序をヒューリスティクスで計算する手法を取り入れた。
- その他メモリアクセスに関連する様々な最適化を取り入れた。
と言えるでしょう。
この手法を読んだ感想ですが、何か全てを変えるような新しい手法を生み出したというよりは、既存の手法を超絶技巧で組み合わせてとんでもない性能を叩き出しているという評価が近いです。 例えばテンソルの和計算に関しても、10x10のqubitで(1+40+1)depthの計算を行う際に、テンソルのサイズが単精度浮動小数を用いても一つのCore Groupのメモリサイズである16GBを超えかけたりと、スパコンのスペックギリギリを渡り歩いていることがわかります。論文に明示的に記していない最適化技法も多々あるでしょう。スパコンのスペックがあと少し足りなかったら今回の結果は実現できなかった可能性もありますし、この論文で「スパコンで300秒で解かれてしまったから量子コンピューターはダメだ」といった感想を持つのは早すぎる気もします。
その一方で、古典コンピュータのアルゴリズムにもupdateがあったら面白いとも感じています。本論文は結局のところテンソル計算を如何に効率的に実行するかに全力を注いでいますが、テンソルネットワークが量子状態を表現するための最適な手法であることは証明されているわけではなく*14、例えば量子状態をボルツマンマシンを用いて表現する手法 *15も存在します。 このような手法も組み合わせることで古典コンピュータではるかに効率の良いシミュレータを作ることもできるかもしれません。
さらなる手法が発展することを期待し、本論文の解説を終わりたいと思います。
参考文献
*1: F. Arute, et al., Quantum supremacy using a programmable superconducting processor. Nature 574, 505–510 (2019)
*3: Y. Liu, et al., Closing the “Quantum Supremacy” Gap, Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (2021)
*4: 作用させるゲートによってはindex ( )の和をまで取らずに適当な値で打ち切ってしまっても大した誤差が生まれない状況もありえます。その場合には例えば を小さくしたテンソルネットワークを用意し、各係数パラメータを最小二乗法で最適化する手法などを取ることができるので、depthが増えても古典コンピュータで容易に扱うことができるのですが、今回扱うようなランダム量子回路の場合はそううまくはいかないはずです。
*5: C. Guo, et al., General-Purpose Quantum Circuit Simulator with Projected Entangled-Pair States and the Quantum Supremacy Frontier, Phys. Rev. Lett. 123, (2019)
*7: J. Gray and S. Kourtis, Hyper-Optimized Tensor Network Contraction, Quantum 5, 410 (2021)
*8:この論文を見てみると、手法の一つとしてコミュニティ検出問題に落とし込んだり、どの和を優先して計算するかというタスクそのものをボルツマン分布を用いてメタヒューリスティクスで選ぶ等を行っているみたいなので、量子アニーリングとの親和性が何となくありそうです。
*10: 最初と最後にあえて1を足しているのは、ランダム量子回路のスキームにおいて最初と最後にアダマールゲートを全てのqubitに作用させるためです。
*11: E. Pednault, J. A Gunnels, G. Nannicini, L. Horesh, T. Magerlein, E. Solomonik, and R. Wisnieff. Breaking the 49-qubit barrier in the simulation of quantum circuits. arXiv preprint arXiv:1710.05867, 15, (2017)
*13: P. Springer and P. Bientinesi, Design of a High-Performance GEMM-like Tensor–Tensor Multiplication, ACM Trans. Math. Softw. 44, 1 (2018)
*15: Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada Phys. Rev. B 96, 205152 (2017)