テンソルネットワークで量子フーリエ変換してみた

この記事は量子コンピューター Advent Calendar 2022の25日目の記事です。
 
はじめまして、株式会社Jij の西村と申します。qiitaではj-i-k-o というIDです。
量子フーリエ変換がテンソルネットワークで効率よく計算できる論文[1]に興味を持ったので、論文の概略を追う + テンソルネットワークシミュレータで量子フーリエ変換を動かしてみました。
この記事では
  • 量子フーリエ変換においてエンタングルメントが小さく抑えることのできる論文[1]の概略の紹介
  • テンソルネットワークシミュレータでの実験
を行おうと思います。
 

量子フーリエ変換

量子コンピュータを用いるアルゴリズムにおいて、サブルーチンとしてよく用いられるアルゴリズムとして、主にGroverの探索アルゴリズムと量子フーリエ変換があります。特に量子フーリエ変換は、ショアのアルゴリズム [2,3]で素因数分解を行う際に用いる位相推定のために用いられることからも有名です。
詳細な解説はこちらのQuantum Native Dojoの資料に譲りますが、量子フーリエ変換は要素のの離散フーリエ変換に相当する
に対応する状態へ量子状態を変化させるアルゴリズムとなります。つまり、
の変換を行うアルゴリズムです。
ここで、は整数の二進数表示に対応する量子状態を省略して表記しています。例えば4qubitの場合 (の場合)、 = になります。
特にこの逆変換である
を見ると、波動関数の係数であるの中に存在する位相 に相当する情報を量子ビット へ移すことができるため、量子計算の途中で発生した位相の情報を観測によって読み取る際に多く用いられます。
 
数式より、量子フーリエ変換を行う量子回路を作成するには
のような演算を行う量子回路を作成すれば良いことがわかります。
これを実現する量子回路は例えば4qubitだと次のような、アダマールゲートと制御位相ゲートを用いたものになります。[4]
 
ここで、を2進数表記である ()で書き直して各量子ゲート毎の量子状態の変化を追うと次のようになります。
ここで、 は2進小数を表します。
これを量子ビットが一般のの場合に拡張すると、最終結果は
 
となり、最後の量子ビットをSWAPゲートを用いて逆順に並び替えると量子フーリエ変換と同じ回路が出来上がります。

量子回路のOperator Schmidt Decomposition

量子回路を取り扱うにあたり、テンソルネットワーク形式で効率よく計算できるかどうかの指標となる量として、Operator Schmidt Decomposition を用いて得られる量の分布があります。
テンソルネットワーク形式の計算、文献についてはこちらに筆者が書いた記事があるのでそちらを参照ください。
 
今、ある量子回路が下の図のようにあるとして、
これを下の図のように量子ビットが分かれている2つの演算子に分割することを考えます。
これを式で表すと次のようになります。この分解をOperator Schmidt Decomposition と呼びます。
ここで、はそれぞれ規格化のため、
を満たすとします。
この中での分布がとても重要な量となります。というのも、もしが全てのに対して同じ程度の大きさで、かつがとても大きくなる場合、上の図の真ん中の垂直な二重線に相当する箇所の結合が無視できず、テンソルネットワーク形式で扱った際に計算コストの増大につながってしまうからです。
その一方。もしに対して素早く (例えば指数的に)減少するのであれば、 のような量を用意して、
としてもが大きい箇所でがほぼ0に近い値を取ることから計算誤差が無視できるので、テンソルネットワーク形式で計算誤差を抑えたまま効率的な計算が可能となります。[5]

The Quantum Fourier Transform Has Small Entanglement

量子フーリエ変換、Operator Schimidt Decompositionについて整理したところで、メインの論文[1] の主張に触れようと思います。
量子フーリエ変換は
  1. 上記で示した以下の量子回路
  1. SWAPゲートを用いて量子ビットを逆順に並び替える操作
 
の2段階によってなされますが、1番目と2番目の手順を合わせた量子フーリエ変換全体の回路についてOperator Schimidt Decompositionを行った結果、に対して一様になる、つまりテンソルネットワーク形式で効率よく扱えるわけでもないという結果が [6,7] にて示されています。しかし本当に量子フーリエ変換として重要な箇所はこのうち1番目の操作だけであり、この操作だけに限ればに対して指数的に減少する、つまりテンソルネットワーク形式で効率よく扱える可能性があることが数値計算で示唆されていました [8]。本論文 [1] ではこの結果に着目して、この1番目の手順についてみると、
となる結果を得ています。これより、に対して指数的に減少しているため、テンソルネットワーク形式で量子フーリエ変換が効率よく扱えるという結論を出しています。この証明については論文 [1] で詳しく語られていますが、かいつまんでまとめると次のようになります。
  • 量子フーリエ変換の2段階の手順のうち1番目の手順をとして、が次の要素に分割できることを証明
    • qubitを用いた量子フーリエ変換
    • 制御位相ゲートの集まりである
    • qubitを用いた量子フーリエ変換
  • 以上を用いて、 と分割できる。任意のに対してこの分割が成り立つ。[9]
    • 文献 [1]のFig. 2より引用
      文献 [1]のFig. 2より引用
  • このなかで一番Operator Schimidt Decompositionに効いてくるのがのため、 だけに着目すれば良い
  • をスペクトル分解し、Operator Schmidt Decompositionの以下の数式へ当てはめることで の分布を計算
  • の特異値を計算することが の固有値を計算することに対応するため、その固有値を計算する過程で次のような を固有値とする固有方程式が現れる。
  • この固有方程式がperiodic discrete prolate spheroidal sequences と呼ばれるクラスに属し、特にqubit数であるに対し の極限をとることで、この方程式がdiscrete prolate spheroidal sequences と呼ばれるクラスに属するため、このクラスに対する既存の結果を使うことでの上限がに対し指数的に減少する結果が得られ、これを用いることで最終的に求めたかった以下の結果を得ることができる。
    • この結果により、Operator Schimidt Decompositionにおける特異値が指数的に減少することが分かります。
      この後に、実際に量子フーリエ変換に相当するユニタリ演算子を、ボンド次元が有限のMatrix Product Operator形式で実際に構築することも行っており、むしろこっちのほうが論文のメインの内容となるのですが、これについては本記事では割愛させていただきます。
      💡
      2022/12/25 追記: この記事を書いたあとに気づいたのですが、既に同じ論文を解説した記事が量子コンピューター Advent Calendar 2022の一日目にありました。 こちらの記事では量子インスパイヤド離散フーリエ変換としての本アルゴリズムの側面にも触れており、面白いので合わせてご覧ください。
       

テンソルネットワークシミュレータを用いた計算

ここまで論文[1] の主張、証明を見てきましたが、実際にテンソルネットワークシミュレータを用いて量子フーリエ変換を行って、量子ビット数に対する実行速度がどの程度になるかを見ていきたいと思います。
この計算を行うにあたって、現在開発中のテンソルネットワークシミュレータを用いてみようと思います。以前に書いた記事であるこちら とほとんど同じ、time-dependent DMRGを用いた計算を内部で行っています。
論文 [1] では量子フーリエ変換に相当する演算子をMPO形式に書き直して数値計算を行っていますが、本テンソルネットワークシミュレータを使い、演算子をMPO形式を使わず直接計算を行った場合に、パフォーマンスがどのようになるのかを見るのが目的となります。
このテンソルネットワークシミュレータは株式会社Jijのインターンの方によって作成されています。
こちらのinstallガイドに従えばインストールできるので、みなさんの環境でお試しいただくことが可能です。
 
まだドキュメントをきちんと用意していない段階なので本記事に書く形式となりますが、使い方は以下の通りとなります。pythonで動きます。
 
pythonで試しに全てのqubitが重ね合わさった状態を、システムサイズ100qubit程度で量子フーリエ変換してみます。コードは以下の通りです。もし量子フーリエ変換が意図通り行われているならば が出力されるはずです。[10]
 
結果は次のようになるため、たしかに が計算できていることが分かります。手元の環境 (Intel Core i7-8565U @ 8x 4.6GHz) ではおよそ20秒ぐらいで計算が完了しました。
このアルゴリズムでは特異値分解のcutoffのみを設定しているため、特異値が一様な値を取るような状況が発生するとtruncationがうまく行かず、計算時間が指数的に増大する要因となります。しかし、計算の途中でいきなり計算速度が低下するような振る舞いもないため、エンタングルメントが小さい状態を維持てきていると考えられます。
また、回路を少し変えてみて、SWAPゲートの一部をコメントアウトして次のような構成にしてみた場合、
実行時間を測ってみると途中で計算速度が大幅に減少してしまい、測定不能レベルで遅くなってしまうことが分かります。これにより、量子フーリエ変換とは違う適当な演算をおこなうと、エンタングルメントが増えていってしまい、cutoffで指定した誤差1e-5を達成するための特異値分解のtruncationもうまくいかないため計算時間が大きくなってしまうことが示唆されます。
 
初期状態として(直積状態ではない) HゲートとCNOTゲートを作用させた状態から初め、量子フーリエ変換を完了するまでにかかる時間を以下のコードで簡単にベンチマークを取ってみました。
 
グラフは次のようになり、計算時間の依存性はなく 依存性があることが分かります。

まとめ

今回の記事にて、
  • 量子フーリエ変換においてエンタングルメントが小さく抑えることのできる論文[1]の概略の紹介
  • テンソルネットワークシミュレータでの実験
を行いました。
量子フーリエ変換が古典コンピュータでシミュレートできてしまうといった内容だけみると、「Shorの素因数分解アルゴリズムも同じようにできるの?」といった疑問が生まれますが、Shorの素因数分解アルゴリズムはorder-findingアルゴリズムのために
のようなユニタリ行列の位相推定 (解説はこちら)を行う必要があり、おそらくこのユニタリ行列の計算でエンタングルメントが増大してしまうのではという可能性があります。実際に検証してみると面白そうです。
また、本テンソルネットワークシミュレータは注釈 [10] にもあるように論文[1]での状況設定と違う形ではありますが、それでも計算時間が多項式で抑えられることが確認できました。論文[1]のようにMatrix Product Operatorをちゃんと作って実験すればもっとパフォーマンスが上がるかもしれません。とはいえ真面目にMatrix Product Operatorを作らずにtime-dependent DMRGを用いて計算を行う本シミュレータでもある程度のパフォーマンスを出せたのは意外な結論でした。
 

[1] J. Chen, E. M. Stoudenmire, and S. R. White, (2022).
[2] P. W. Shor, Proceedings 35th Annual Symposium on Foundations of Computer Science (1994).
[3] P. W. Shor, SIAM Journal on Computing 26, 1484 (1997).
[4] この回路は論文[1]から取ってきたものですが、実はNielsen & Chuang やQuantum Native Dojoにある回路と若干異なります。とはいえ記事中で示している通り結局最終的には同じ式になるので問題ありません。実際、制御位相ゲートにてcontrol qubitとtarget qubitを入れ替えても同じゲートになる性質を利用するとNielsen & Chuang やQuantum Native Dojoにある回路と同じ回路にすることができます。
[5] ここで、テンソルネットワーク形式における特異値分解のプロセスで値の小さい特異値を無視するtruncationと呼ばれるプロセス、つまりに置き換えてしまうプロセスを経ないと全ての特異値を保持してしまうため、この方法に意味がないことに注意です。
[6] M. A. Nielsen et al., Physical Review A 67, (2003).
[7] J. Tyson, Journal of Physics A: Mathematical and General 36, 6813 (2003).
[8] K. Woolfe, C. Hill, and L. Hollenberg, Quantum Information and Computation 17 (2014)
[9] 一般の証明についての詳細は論文[1]のAppendix Aにまとめられていますが、アダマールゲートや制御位相ゲートを他の量子ゲートとぶつからないように順序を入れ替えることで直感的に示すことができます。
[10] 本シミュレータでは制御位相ゲートをかける際に量子ビットの移動のためにSWAPゲートを併用していますが、論文[1] では制御位相ゲートをそのままMatrix Product Operatorに変換してから量子状態に作用させる手法を提案しているため、両者の計算量に違いが出てくる可能性があります。おそらく本記事で紹介した内容だとSWAPゲートを作用させている間に計算時間が増大し、SWAPゲートを処理し終わった段階でエンタングルメントが減少し計算が速くなる振る舞いがみられると予想されます。