Skip to main content
IBM Quantum Platform

クリロフ量子対角化

クリロフ量子対角化(KQD)についてのこのレッスンでは、以下のことにお答えします：

  • 一般的にクリロフ法とは何ですか？
  • なぜクリロフ法が有効なのか、どのような条件下で有効なのか？
  • 量子コンピューターはどのような役割を果たしているのか？

量子計算は文献 [1] に基づいている。

以下のビデオは、古典コンピューティングにおけるクリロフ法の概要、その使用動機、そして量子コンピューティングがそのワークストリームでどのような役割を果たせるかを説明している。 続くテキストでは、より詳細な説明があり、クリロフ法を古典的に、また量子コンピューターを用いて実装している。

1. クリロフ法の導入

クリロフ部分空間法は、 クリロフ部分空間と呼ばれるものを中心に構築されたいくつかの手法のいずれかを指す。 これらの完全なレビューはこのレッスンの範囲外であるが、参考文献 [2-4] はすべて、より実質的な背景を与えることができる。 ここでは、クリロフ部分空間とは何か、それが固有値問題を解く上でどのように、そしてなぜ有用なのか、そして最後に量子コンピューター上でどのように実装できるのかに焦点を当てる。

定義対称的な正の半定値 N×NN\times N 行列 AA が与えられたとき，次数 rr のクリロフ空間 Kr\mathcal{K}^r は，行列 AA の高次のべき乗 r−1≤Nr-1\leq N までを参照ベクトル ∣v⟩\vert v \rangle と掛け合わせることによって得られるベクトルによってスパンされる空間である．

Kr=span{∣v⟩,A∣v⟩,A2∣v⟩,...,Ar−1∣v⟩}\mathcal{K}^r = \text{span}\left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}

上記のベクトルはクリロフ部分空間と呼ばれるものにまたがっているが、直交していると考える理由はない。 多くの場合、 グラム・シュミット直交化と同様の反復直交化処理を用いる。 ここでは、それぞれの新しいベクトルが生成されるときに他のベクトルに対して直交するように作られているため、プロセスが若干異なる。 この文脈では、これをアーノルディ反復と呼ぶ。 最初のベクトル ∣v⟩|v\rangle から始めて、次のベクトル A∣v⟩A|v\rangle を生成し、 ∣v⟩|v\rangle への投影を減算することによって、この2番目のベクトルが最初のベクトルと直交することを保証する。つまり

∣v0⟩=∣v⟩∣∣∣v⟩∣∣∣v1⟩=A∣v⟩−⟨v0∣A∣v⟩∣v0⟩∣∣A∣v⟩−⟨v0∣A∣v⟩∣v0⟩∣∣\begin{aligned} |v_0\rangle &=\frac{|v\rangle}{\left|\left| |v\rangle \right|\right|}\\ |v_1\rangle &=\frac{A|v\rangle-\langle v_0|A|v\rangle |v_0\rangle}{\left|\left|A|v\rangle-\langle v_0|A|v\rangle |v_0\rangle \right|\right|} \end{aligned}

∣v0⟩⊥∣v1⟩|v_0\rangle \perp |v_1\rangle であることが、今や容易に理解できる。なぜなら、

⟨v0∣v1⟩=⟨v0∣A∣v⟩−⟨v0∣A∣v⟩⟨v0∣v0⟩∣∣A∣v⟩−⟨A∣v⟩∣v0⟩∣v0⟩∣∣=0\langle v_0 | v_1\rangle=\frac{\langle v_0 | A|v\rangle-\langle v_0 |A|v\rangle\langle v_0|v_0\rangle}{\left|\left| A|v\rangle-\langle A|v\rangle|v_0\rangle |v_0\rangle \right|\right|}=0

次のベクトルも同様に、前の2つのベクトルと直交するようにする：

∣v2⟩=A∣v1⟩−⟨v0∣A∣v1⟩∣v0⟩−⟨v1∣A∣v1⟩∣v1⟩∣∣A∣v1⟩−⟨v0∣A∣v1⟩∣v0⟩−⟨v1∣A∣v1⟩∣v1⟩∣∣|v_2\rangle=\frac{A |v_1\rangle-\langle v_0|A |v_1\rangle |v_0\rangle-\langle v_1|A |v_1\rangle |v_1\rangle}{\left|\left| A |v_1\rangle-\langle v_0|A |v_1\rangle |v_0\rangle-\langle v_1|A |v_1\rangle |v_1\rangle\right|\right|}

このプロセスをすべてのベクトル（ rr ）について繰り返すと、クリロフ空間の完全な正規直交基底ができる。 mm 直交ベクトルは必然的に全空間にまたがるので、ここでの直交化処理は一度 r>mr>m、ゼロになることに注意。 また、 AA の固有ベクトルとなるベクトルがあれば、それ以降のベクトルはすべてそのベクトルの倍数となるため、この処理はゼロとなる。

1.1 簡単な例：手計算によるクリロフ

ここで、些細に小さな行列に対するクリロフ部分空間の生成の過程を、ステップを踏んで見てみよう。 AA の初期行列から始める：

A=(4−10−14−10−14)A=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}

この小さな例では、手作業でも簡単に固有ベクトルと固有値を求めることができる。 ここでは数値解を示す。

# One might use linalg.eigh here, but later matrices may not be Hermitian. So we use
# linalg.eig in this lesson.

import numpy as np

A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]])
eigenvalues, eigenvectors = np.linalg.eig(A)
print("The eigenvalues are ", eigenvalues)
print("The eigenvectors are ", eigenvectors)

Output:

The eigenvalues are  [2.58578644 4.         5.41421356]
The eigenvectors are  [[ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]
 [ 7.07106781e-01  1.37464400e-16 -7.07106781e-01]
 [ 5.00000000e-01  7.07106781e-01  5.00000000e-01]]

後で比較するために、ここに記録しておく：

a0=2.59,∣0⟩=(1/2−2/21/2)a1=4,∣1⟩=(2/20−2/2)a2=5.41,∣2⟩=(1/22/21/2)\begin{aligned} a_0&=2.59,&|0\rangle&=&\begin{pmatrix}1/2\\-\sqrt{2}/2\\1/2\end{pmatrix}\\ \\ a_1&=4,&|1\rangle&=&\begin{pmatrix}\sqrt{2}/2\\0\\-\sqrt{2}/2\end{pmatrix}\\ \\ a_2&=5.41,&|2\rangle&=&\begin{pmatrix}1/2\\\sqrt{2}/2\\1/2\end{pmatrix} \end{aligned}

クリロフ部分空間 rr の次元を大きくするにつれて、このプロセスがどのように機能するか（あるいは機能しないか）を研究したい。そのために、このプロセスを適用する：

  • ランダムに選ばれたベクトル ∣v⟩|v\rangle （上記のようにすでに正規化されている場合は ∣v0⟩|v_0\rangle ）と呼ぶ）から始まる完全なベクトル空間の部分空間を生成する。
  • その部分空間に完全行列 AA を射影し，その射影行列 A~\tilde{A} の固有値を求める．
  • グラム・シュミット直交化に似たプロセスで、正規直交であることを保証しながら、より多くのベクトルを生成して部分空間のサイズを大きくする。
  • AA をより大きな部分空間に投影し、その結果得られる行列の固有値 A~\tilde{A} を求める。
  • これを固有値が収束するまで（あるいはこのおもちゃのケースでは、元の行列 AA の全ベクトル空間にまたがるベクトルが生成されるまで）繰り返す。

クリロフ法の通常の実装では、クリロフ部分空間を構築するたびに、その部分空間に射影された行列の固有値問題を解く必要はない。 希望する次元の部分空間を構成し、その部分空間に行列を射影し、射影行列を対角化することができる。 各サブスペース次元での投影と対角化は、収束をチェックするためだけに行われる。

r=1r=1 の寸法：

ランダムなベクトルを選ぶ

∣v0⟩=(100)|v_0\rangle=\begin{pmatrix}1\\0\\0\end{pmatrix}

まだ正規化されていない場合は、正規化する。

ここで、行列 AA をこの1つのベクトルの部分空間に射影する：

A~0=⟨v0∣A∣v0⟩=(100)(4−10−14−10−14)(100)=(4)\tilde{A}_0=\langle v_0| A|v_0\rangle=\begin{pmatrix}1&0&0\end{pmatrix}\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}=(4)

これは、クリロフ部分空間への行列の射影であり、この部分空間にはベクトル ∣v0⟩|v_0\rangle が1つだけ含まれている。この行列の固有値はトリビュラルに4である。 これは、 AA の固有値（この場合は1つだけ）の0次推定値と考えることができる。稚拙な推定値ではあるが、桁としては正しい。

r=2r=2 の寸法：

AA、前のベクトルに対して次のベクトルを生成する：

A∣v0⟩=(4−10−14−10−14)(100)=(4−10)A|v_0\rangle=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}=\begin{pmatrix}4\\-1\\0\end{pmatrix}

ここで、直交性を確保するために、このベクトルの前のベクトルへの射影を差し引く。

∣v1⟩=A∣v0⟩−⟨v0∣A∣v0⟩∣v0⟩|v_1\rangle=A|v_0\rangle-\langle v_0 |A|v_0\rangle|v_0\rangle ∣v1⟩=(4−10)−(100)(4−10)(100)=(0−10)|v_1\rangle=\begin{pmatrix}4\\-1\\0\end{pmatrix}-\begin{pmatrix}1& 0& 0\end{pmatrix}\begin{pmatrix}4\\-1\\0\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}=\begin{pmatrix}0\\-1\\0\end{pmatrix}

まだ正規化されていない場合は、正規化する。 この場合、ベクトルはすでに正規化されている

∣v1⟩=(0−10)|v_1\rangle=\begin{pmatrix}0\\-1\\0\end{pmatrix}

ここで、行列Aをこれら2つのベクトルの部分空間に射影する：

A~1=(1000−10)(4−10−14−10−14)(100−100)=(1000−10)(41−1−401)=(4114)\tilde{A}_1= \begin{pmatrix} 1&0&0\\0&-1&0 \end{pmatrix} \begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1&0\\0&-1\\0&0\end{pmatrix}=\begin{pmatrix}1&0&0\\0&-1&0\end{pmatrix}\begin{pmatrix}4&1\\-1&-4\\0&1\end{pmatrix}=\begin{pmatrix}4&1\\1&4\end{pmatrix}

この行列の固有値を求める問題が残っている。 しかし、このマトリックスはフルマトリックスよりもわずかに小さい。 非常に大きな行列を含む問題では、この小さな部分空間で作業することが非常に有利になることがある。

det⁡(A1~−λI)=0\det(\tilde{A_1}-\lambda I)=0 ∣4−λ114−λ∣=(4−λ)2−1=0\begin{vmatrix} 4-\lambda&1\\1&4-\lambda\end{vmatrix} =(4-\lambda)^2-1=0 4−λ=±1→λ=3,54-\lambda=±1→\lambda=3,5

これはまだ良い推定値とは言えないが、0次推定値よりは良い。 プロセスが明確であることを確認するため、もう1度、この反復を実施する。 なぜなら、次の反復で 3x3 行列を対角化することになり、時間も計算量も節約できないからである。

r=3r=3 の寸法：

ここで、前のベクトルに対してAを操作して、部分空間の次のベクトルを生成する：

A∣v1⟩=(4−10−14−10−14)(0−10)=(1−41)A|v_1\rangle=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}0\\-1\\0\end{pmatrix}=\begin{pmatrix}1\\-4\\1\end{pmatrix}

ここで、直交性を確保するために、このベクトルを前のベクトル両方に射影したものを差し引く。

∣v2⟩=A∣v1⟩−⟨v0∣A∣v1⟩∣v0⟩−⟨v1∣A∣v1⟩∣v1⟩∣v2⟩=(1−41)−(100)(1−41)(100)−(0−10)(1−41)(0−10)=(001)\begin{aligned} |v_2\rangle&=A|v_1\rangle-\langle v_0 |A|v_1\rangle|v_0\rangle-\langle v_1 |A|v_1\rangle|v_1\rangle\\ |v_2\rangle&=\begin{pmatrix}1\\-4\\1\end{pmatrix}-\begin{pmatrix}1& 0& 0 \end{pmatrix}\begin{pmatrix}1\\-4\\1\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}-\begin{pmatrix}0&-1& 0\end{pmatrix}\begin{pmatrix}1\\-4\\1\end{pmatrix}\begin{pmatrix}0\\-1\\0\end{pmatrix}=\begin{pmatrix}0\\0\\1\end{pmatrix} \end{aligned}

まだ正規化されていない場合は、正規化する。 この場合、ベクトルはすでに正規化されている

∣v2⟩=(001)|v_2 \rangle=\begin{pmatrix}0\\0\\1\end{pmatrix}

ここで、行列 AA をこれらのベクトルの部分空間に射影する：

A~2=(1000−10001)(4−10−14−10−14)(1000−10001)=(4−101−410−14)(1000−10001)=(410141014)\tilde{A}_2=\begin{pmatrix}1&0&0\\0&-1&0\\0&0&1\end{pmatrix}\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}\begin{pmatrix}1&0&0\\0&-1&0\\0&0&1\end{pmatrix}=\begin{pmatrix}4&-1&0\\1&-4&1\\0&-1&4\end{pmatrix}\begin{pmatrix}1&0&0\\0&-1&0\\0&0&1\end{pmatrix}=\begin{pmatrix}4&1&0\\1&4&1\\0&1&4\end{pmatrix}

ここで固有値を決定する：

det⁡(A~2−λI)=0\det(\tilde{A}_2-\lambda I)=0 ∣4−λ1014−λ1014−λ∣=(4−λ)((4−λ)2−1)−(4−λ)=0\begin{vmatrix}4-\lambda&1&0\\1&4-\lambda&1\\0&1&4-\lambda\end{vmatrix} = (4-\lambda)((4-\lambda)^2-1)-(4-\lambda)=0\\ 4−λ=0,4−λ=±21/2→λ=4−21/2,4,4+21/2≈2.59,4,5.414-\lambda=0,4-\lambda=±2^{1/2}→\lambda=4-2^{1/2},4,4+2^{1/2}≈2.59,4,5.41

これらの固有値は、まさに元の行列 AA の固有値です。クリロフ部分空間は，元の行列 AA のベクトル空間全体にまたがるように拡張されているので，このようになるに違いない．

この例では、クリロフ法は直接対角化よりも特に簡単には見えないかもしれない。 実際、後の節で述べるように、クリロフ法が有利なのはある行列次元以上だけである。これは、非常に大きな行列の固有値／固有ベクトル問題を解くのに役立つことを目的としている。

非常に大きな行列をクリロフ部分空間に投影した画像。つまり、左側に行列を作るクリロフ・ベクトルの行、ハミルトニアン、右側にクリロフ・ベクトルの列がある。

これは "手で "作業した唯一の例であるが、以下のセクション2には計算例を示す。

用語の明確化

よくある誤解は、与えられた問題に対するクリロフ部分空間は1つしかないというものだ。 しかしもちろん、我々の行列を適用できる初期ベクトルは多数存在するので、可能なクリロフ部分空間も多数存在する。 クリロフ部分空間 」という表現は、特定の例についてすでに定義された特定のクリロフ部分空間を指す場合にのみ使用する。 一般的な問題解決アプローチについては、「a クリロフ部分空間」を参照する。 最後に明確にしておきたいのは、"クリロフ空間 "という言い方が有効だということだ。 初期空間から部分空間への行列の射影という文脈で使われるため、「クリロフ部分空間 」と呼ばれることも多い。 その文脈に沿って、ここではほとんど亜空間と呼ぶことにする。

理解度チェック

なぜ(a)クリロフ部分空間 rr の次元を、注目する行列の次元 NN 以上に拡張することが有用でなく、(b)不可能なのかを説明せよ。

  • (a) ベクトルを生成する際に直交正規化を行っているため、 NN 個のそのようなベクトルの集合は完全基底を成す。つまり、それらの線形結合を用いることで、その空間内の任意のベクトルを生成することができる。

    (b) 直交化のプロセスは、新しいベクトルのすべての先行するベクトルへの射影を差し引くことから成る。 もし、それまでのすべてのベクトルがベクトル空間全体を張っているならば、部分空間全体への射影を引いても、常に零ベクトルが残ることになる。

ある研究者が、小さなおもちゃのような行列にクリロフ法を適用して実演しているとしよう。 行列 AA と初期ベクトル ∣ψ⟩|\psi\rangle の選択に何か問題があるのでしょうか？

A=(213123335)A=\begin{pmatrix}2&1&3\\1&2&3\\3&3&5\end{pmatrix}

および

∣ψ⟩=12(1−10).|\psi\rangle=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\-1\\0\end{pmatrix}.

  • あなたの同僚は誤って初期ベクトルに固有ベクトルを選んでしまった。 初期ベクトルに対して行列を作用させると、同じベクトルを固有値でスケーリングして返すだけである。 これでは、次元が増加する部分空間は生成されない。 同僚に、異なる初期ベクトルを選ぶように助言し、それが固有ベクトルでないことを確認する。

与えられた行列にクリロフ法を適用し、適切な新しい初期ベクトルを選択する。 あなたのクリロフ部分空間の0次および1次における最小固有値の推定値を書き出してください。

A=(110111011)A=\begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix}

  • 初期ベクトルの選択によって、多くの答えが考えられる。 私たちはそれを選ぶ：

    ∣v0⟩=13(111).|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}.

    ∣v1⟩|v_1\rangle を得るには、 ∣v0⟩|v_0\rangle に AA を一度適用し、次に ∣v1⟩|v_1\rangle を ∣v0⟩|v_0\rangle に対して直交させる。

    A∣v0⟩=(110111011)13(111)=13(232)A|v_0\rangle=\begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}A∣v0⟩−⟨v0∣A∣v0⟩∣v0⟩=13(232)−13(111)13(232)13(111)=13(232)−7313(111)=32(−1/32/3−1/3)A|v_0\rangle - \langle v_0|A|v_0\rangle |v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix} - \frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}-\frac{7}{3}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}=\sqrt{\frac{3}{2}}\begin{pmatrix}-1/3\\2/3\\-1/3\end{pmatrix}

    0次では、クリロフ部分空間への射影は次のようになる

    ⟨v0∣A∣v0⟩=13(111)(110111011)13(111)=73\langle v_0|A|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{7}{3}

    1次では、このクリロフ部分空間への射影は次のようになる

    ⟨V1∣A∣V1⟩=(131313−1623−16)(110111011)(13−16132313−16)\langle V^1|A|V^1\rangle=\begin{pmatrix}\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\-\sqrt{\frac{1}{6}}&\sqrt{\frac{2}{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \begin{pmatrix}\frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\\\frac{1}{\sqrt{3}}& \sqrt{\frac{2}{3}} \\ \frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix}

    これは手作業でもできるが、numpyを使うのが最も簡単である：

    import numpy as np
vstar = np.array([[1/np.sqrt(3),1/np.sqrt(3),1/np.sqrt(3)],[-1/np.sqrt(6),np.sqrt(2/3),-1/np.sqrt(6)]]
)
A = np.array([[1, 1, 0],
              [1, 1, 1],
              [0, 1, 1]])
v = np.array([[1/np.sqrt(3),-1/np.sqrt(6)],[1/np.sqrt(3),np.sqrt(2/3)],[1/np.sqrt(3),-1/np.sqrt(6)]])
proj = vstar@A@v
print(proj)
eigenvalues, eigenvectors = np.linalg.eig(proj)
print("The eigenvalues are ", eigenvalues)
print("The eigenvectors are ", eigenvectors)

    outputs:

    [[ 2.33333333  0.47140452]
 [ 0.47140452 -0.33333333]]
The eigenvalues are  [ 2.41421356 -0.41421356]
The eigenvectors are  [[ 0.98559856 -0.16910198]
 [ 0.16910198  0.98559856]]

    最小固有値の推定値は -0.414 である。

1.2 クリロフ法の種類

「クリロフ部分空間法」は、大規模な線形システムや固有値問題を解くために使われる、いくつかの反復技法のいずれかを指す。 これらに共通するのは、クリロフ部分空間から近似解を構成することである

Kr(A,∣v⟩)=span{∣v⟩,A∣v⟩,A2∣v⟩,...,Ar−1∣v⟩},\mathcal{K}^r(A,|v\rangle ) = \text{span}\{|v\rangle, A|v\rangle, A^2|v\rangle, ..., A^{r-1}|v\rangle\},

ここで、 ∣v⟩|v\rangle は初期推測値である（文献 [5] を参照）。 収束率、メモリ使用量、全体的な計算コストなどの要素をバランスさせながら、この部分空間から最適な近似を選択する方法が異なる。 このレッスンの焦点は、クリロフ部分空間法の文脈で量子コンピューティングを活用することである。 以下の簡単な定義は、あくまでも文脈のためのものであり、これらの方法をさらに調べるための参考文献も含まれている。

共役勾配（CG）法 ：この方法は、対称な正定値線形システムを解くために使用される [6]。 これは、各反復における誤差のAノルムを最小化し、離散化楕円型PDEから生じる系に特に効果的です [7]。 次のセクションでは、このアプローチを使って、なぜクリロフ部分空間が線形システムの解の改善を探るのに効果的な部分空間なのかを説明する。

一般化最小残差（GMRES）法 ：これは、一般的な非対称線形システムを解くために設計されています。 これは、各反復においてクリロフ空間上の残差ノルムを最小化するもので、ロバストであるが、大規模なシステムではメモリを消費する可能性がある [7]。

最小残差(MINRES)法 ：この方法は、対称不定一次方程式を解くのに用いられる。 GMRESと似ていますが、行列の対称性を利用して計算コストを削減します [8]。

その他の注目すべきアプローチには、固有値問題に対するArnoldiの手法と密接な関係にある完全直交化法（FOM）、 両共役勾配法（ BiCG ）、 誘導次元削減法（IDR）などがある。

1.3 クリロフ部分空間法が機能する理由

ここでは、クリロフ部分空間法が、最急降下法というレンズを通して、固有ベクトル近似の反復精密化を通して行列固有値を近似する効率的な方法であることを動機づける。 基底状態の初期推測が与えられたとき、その初期推測に対する連続的な補正のうち、最速の収束をもたらす空間はクリロフ部分空間であることを論じる。 収束動作の厳密な証明には至っていない。

注目の行列 AA が対称かつ正定値であると仮定する。 このため、我々の議論は上記のCG方式に最も適している。 また、 AA がエルミートでなければならないと主張しているわけでもない（ハミルトニアンであればそうである必要がある）。

私たちは通常、次のような形式の問題を解きたいと考える

A∣x⟩=∣b⟩.A|x\rangle=|b\rangle.

∣b⟩=c∣x⟩|b\rangle=c|x\rangle ここで、 cc は固有値問題のようにある定数である。 しかし、今のところ、我々の問題提起はより一般的なものにとどまっている。

まず、近似解であるベクトル ∣x0⟩|x_0\rangle から始めます。 この推測 ∣x0⟩|x_0\rangle と、セクション 1.1 にある ∣v0⟩|v_0\rangle との間には類似点がありますが、ここではこれらを活用していません。 ∣x0⟩|x_0\rangle にはエラーがあると思われます。このエラーを ∣e0⟩|e_0\rangle と呼んでいます：

∣e0⟩:=∣x⟩−∣x0⟩.|e_0\rangle:=|x\rangle−|x_0\rangle.

また、残差 R0R_0 を次のように定義する：

∣R0⟩=∣b⟩−A∣x0⟩.|R_0\rangle=|b\rangle−A|x_0\rangle.

ここでは、クリロフ部分空間 rr の次元と残差を区別するために、大文字の RR。

x とラベル付けされた真の固有ベクトル、x₀ とラベル付けされた推定値、およびこれら2つの間の誤差を表すグラフ。

ここで、次のような補正を加えたい

∣x1⟩=∣x0⟩+∣p0⟩,|x_1\rangle=|x_0\rangle+|p_0\rangle,

これは私たちの近似値を向上させるものである。 ここで ∣p0⟩|p_0\rangle、まだ決定していないベクトルがある。 修正後の誤差を ∣e1⟩|e_1\rangle とする。 結果

∣e1⟩=∣x⟩−∣x1⟩=∣x⟩−(∣x0⟩+∣p0⟩)=∣e0⟩−∣p0⟩.|e_1\rangle=|x\rangle−|x_1\rangle=|x\rangle−(|x_0\rangle+|p_0\rangle)=|e_0\rangle−|p_0\rangle.

真の固有ベクトルと、初期推定値の更新値。 更新された推定値は、真の固有ベクトルにより近い。

私たちは、私たちの誤差が私たちの行列によって変換されたときにどのように振る舞うかに興味がある。 そこで、誤差の AA -normを計算してみよう。 それは

∥∣e0⟩−∣p0⟩∥A2=(⟨e0∣A−⟨p0∣A)(∣e0⟩−∣p0⟩)=⟨e0∣A∣e0⟩−⟨e0∣A∣p0⟩−⟨p0∣A∣e0⟩+⟨p0∣A∣p0⟩=⟨e0∣A∣e0⟩−2⟨e0∣A∣p0⟩+⟨p0∣A∣p0⟩=d−2⟨R0∣p0⟩+⟨p0∣A∣p0⟩,\begin{aligned} ∥|e_0\rangle−|p_0\rangle∥_A^2&=\left(\langle e_0|A−\langle p_0|A\right)\left(|e_0\rangle−|p_0\rangle\right)\\ & = \langle e_0|A|e_0 \rangle − \langle e_0|A|p_0\rangle − \langle p_0|A|e_0\rangle+\langle p_0|A|p_0\rangle\\ & = \langle e_0|A|e_0\rangle−2\langle e_0|A|p_0\rangle+\langle p_0|A|p_0\rangle\\ & = d−2\langle R_0|p_0\rangle +\langle p_0|A|p_0\rangle, \end{aligned}

ここでは、 AA という対称性、および A∣e0⟩=∣R0⟩A |e_0\rangle = |R_0\rangle という関係式を利用している。ここで、 dd は ∣p0⟩|p_0\rangle とは無関係な定数である。セクション 1.2 で述べたように、誤差の AA ノルムは、最小化の対象として選択できる唯一の量ではないが、適切な量の一つである。 この量が、補正ベクトル ∣p0⟩|p_0\rangle の選択によってどのように変化するかを確認したい。そこで、次のように設定して、関数 ff を定義する

f(∣p0⟩)=⟨p0∣A∣p0⟩−2⟨R0∣p0⟩+d.f(|p_0\rangle)=\langle p_0|A|p_0\rangle−2\langle R_0|p_0\rangle+d.

ff は、 AA -normで測定された補正 ∣p0⟩|p_0\rangle の関数としての誤差 ∣e1⟩|e_1\rangle。 したがって、 f(∣p0⟩)f(|p_0\rangle) ができるだけ小さくなるように、 ∣p0⟩|p_0\rangle を選びたい。 この目的のために、 ff の勾配を計算する。 AA の対称性を用いると、次のようになる

∇f(∣p0⟩)=2(A∣p0⟩−∣R0⟩).\nabla f(|p_0\rangle) = 2(A|p_0\rangle−|R_0\rangle).

勾配は上昇が最も急な方向を指しており、つまりその反対方向が、関数が最も減少する方向、すなわち下降が最も急な方向となります。 初期の推測として、 ∣x0⟩|x_0\rangle （ここで、 ∣p0⟩=0|p_0\rangle=0 ）とすると、 ∇f(0)=−2∣R0⟩\nabla f(0) = -2|R_0\rangle が成り立つ。 したがって、関数 ff は、残差 ∣R0⟩|R_0\rangle の方向で最も減少する。したがって、あるスカラー α0\alpha_0 に対して、ベクトル ∣p0⟩=α0∣R0⟩|p_0\rangle=\alpha_0 |R_0\rangle を加えることで、この初期の選択が最も恩恵を受けることになる。

次のステップでは、やはりベクトル ∣p1⟩|p_1\rangle、その値を現在の近似値に加える。 先ほどと同じ論法で、あるスカラー α1\alpha_1 に対して ∣p1⟩=α1∣R1⟩|p_1\rangle = \alpha_1 |R_1\rangle を選ぶ。このように続けると、 kthk^\text{th} の反復ベクトルは次のようになる

∣xk+1⟩=∣x0⟩+α0∣R0⟩+α1∣R1⟩+⋯+αk∣Rk⟩.|x_{k+1}\rangle=|x_0\rangle+\alpha_0 |R_0\rangle+\alpha_1 |R_1\rangle+⋯+\alpha_k |R_k\rangle.

同様に、 ∣R0⟩|R_0\rangle、 ∣R1⟩|R_1\rangle、...と順番に追加していくことで、改善された推定値を選択する空間を構築したい。 kthk^\text{th} 推定されたベクトルは次のようになる

∣xk+1⟩∈∣x0⟩+span{∣R0⟩,∣R1⟩,…,∣Rk⟩}.|x_{k+1}\rangle\in |x_0\rangle+\text{span}\{|R_0\rangle,|R_1\rangle,…,|R_k\rangle \}.

さて、次の関係を使う

∣Rk+1⟩=∣b⟩−A∣xk+1⟩=∣b⟩−A(∣xk⟩+αk∣Rk⟩)=∣Rk⟩−αkA∣Rk⟩,|R_{k+1}\rangle=|b\rangle−A |x_{k+1}\rangle=|b\rangle−A(|x_k\rangle+\alpha_k |R_k\rangle)=|R_k\rangle−\alpha_k A |R_k\rangle,

つまり

span{∣R0⟩,∣R1⟩,…,∣Rk⟩}=span{∣R0⟩,A∣R0⟩,…,Ak∣R0⟩}.\text{span} \{|R_0\rangle,|R_1\rangle,…,|R_k\rangle \}=\text{span} \{|R_0\rangle,A|R_0\rangle,…,A^{k}|R_0\rangle \}.

つまり、正しい解 ∣x⟩|x\rangle を最も効率的に近似する空間は、行列 AA を ∣R0⟩|R_0\rangle に対して順次作用させることによって構築される空間とまったく同じである。クリロフ部分空間とは*、* 最急降下法の各段階の方向ベクトルによって張られる空間のことである。

最後に、我々はこのアプローチのスケーリングについて数値的な主張をしていないこと、また疎行列に対する比較的な利点について議論していないことを再確認する。 これは、クリロフ部分空間法を使う動機付けと、直感的な感覚を加えることだけを意図している。 ここで、これらの方法の数値的な挙動を探ってみることにする。

理解度チェック

上記のワークフローでは、誤差の AA -normを最小化することを提案した。 基底状態とその固有値を求める際に、他にどのような量を最小化することが考えられるか？

  • 誤差の AA -normの代わりに残差ベクトルを使うことも考えられる。 誤差ベクトルそのものを考慮することが有用な場合もあるだろう。

2. クラシカル計算におけるクリロフ法

このセクションでは、固有値問題を解く際にクリロフ部分空間を活用できるように、アーノルディ反復を計算機的に実装する。 まず小規模な例でこれを適用し、次に対象となる行列のサイズが大きくなるにつれて、計算時間がどのようにスケーリングするかを検証する。 ここでの重要な考え方は、クリロフ空間にまたがるベクトルの生成が、必要な総計算時間に大きく寄与するということである。 必要なメモリは、特定のクリロフ法によって異なる。 しかし、メモリの制約により、従来のクリロフ法の使用が制限されることがある。

2.1 単純な小規模な例

クリロフ部分空間を作る過程で、部分空間内のベクトルを正規直交化する必要がある。 私たちの部分空間 vknown （正規化されているとは仮定していない）から確立されたベクトルと、私たちの部分空間 vnext に追加する候補ベクトルを取り、 vnext を vknown に直交させ、正規化する関数を定義しよう。 さらに、クリロフ部分空間内のすべての確立されたベクトルに対して、このプロセスを段階的に実行する関数を定義して、完全に正規直交する集合を保証しよう。

# vknown is some established vector in our subspace. vnext is one we wish to add,
# which must be orthogonal to vknown.


def orthog_pair(vknown, vnext):
    vknown = vknown / np.sqrt(vknown.T @ vknown)
    diffvec = vknown.T @ vnext * vknown
    vnext = vnext - diffvec
    return vnext


# v is the candidate vector to be added to our subspace. s is the existing subspace.


def orthoset(v, s):
    v = v / np.sqrt(v.T @ v)
    temp = v
    for i in range(len(s)):
        temp = orthog_pair(s[i], temp)
    v = temp / np.sqrt(temp.T @ temp)
    return v

ここで、クリロフベクトルの空間が元の行列の全空間に及ぶまで、繰り返しより大きなクリロフ部分空間を構築する関数を定義しよう。 これにより、クリロフ部分空間の次元の関数として、我々のクリロフ部分空間法を用いて得られた固有値が、どの程度正確な値と一致するかを見ることができる。 重要なのは、我々の関数 krylov_full_build 、クリロフ・ベクトル、投影ハミルトニアン、固有値、所要時間を返すことである。

# Necessary imports and definitions to track time in microseconds
import time


def time_mus():
    return int(time.time() * 1000000)


# This function constructs a Krylov subspace that spans the whole space of the original matrix.
#     Input:
#       v0          : initial vector
#       matrix      : original matrix to be diagonalized
#     Output:
#       ks          : Krylov vectors
#       Hs          : projected Hamiltonians
#       eigs        : eigenvalues
#       k_tot_times : time required for the operation


def krylov_full_build(v0, matrix):
    t0 = time_mus()
    b = v0 / np.sqrt(v0 @ v0.T)
    A = matrix
    ks = []
    ks.append(b)
    Hs = []
    eigs = []
    Hs.append(b.T @ A @ b)
    eigs.append(np.array([b.T @ A @ b]))
    k_tot_times = []

    for j in range(len(A) - 1):
        vec = A @ ks[j].T
        ortho = orthoset(vec, ks)
        ks.append(ortho)
        ksarray = np.array(ks)
        Hs.append(ksarray @ A @ ksarray.T)
        eigs.append(np.linalg.eig(Hs[j + 1]).eigenvalues)
        k_tot_times.append(time_mus() - t0)

    # Return the Krylov vectors, the projected Hamiltonians, the eigenvalues,
    # and the total time required.
    return (ks, Hs, eigs, k_tot_times)

まだかなり小さいが、手作業でやるよりは大きい行列でテストしてみよう。

# Define our small test matrix
test_matrix = np.array(
    [
        [4, -1, 0, 1, 0],
        [-1, 4, -1, 2, 1],
        [0, -1, 4, 3, 3],
        [1, 2, 3, 4, 0],
        [0, 1, 3, 0, 4],
    ]
)

# Give the test matrix and an initial guess as arguments in the function defined above.
# Calculate outputs.
test_ks, test_Hs, test_eigs, text_k_tot_times = krylov_full_build(
    np.array([0.5, 0.5, 0, 0.5, 0.5]), test_matrix
)

最後のステップ（クリロフ空間が元の行列の完全なベクトル空間である場合）で、クリロフ法による固有値が厳密な数値対角化の固有値と正確に一致することを確認することで、我々の関数をチェックすることができる：

print(np.linalg.eig(test_matrix).eigenvalues)
print(test_eigs[len(test_matrix) - 1])

Output:

[-1.36956923  8.43756009  2.9040308   5.34436028  4.68361806]
[-1.36956923  8.43756009  2.9040308   4.68361806  5.34436028]

それは成功した。 もちろん、本当に重要なのは、クリロフ部分空間の次元の関数として、我々の近似がどれだけ優れているかということだ。 我々は基底状態やその他の最小固有値を見つけることにしばしば関心を持つので（そして以下に説明する代数的な理由もある）、クリロフ部分空間の次元の関数として最小固有値の推定値を見てみよう。 それは

def errors(matrix, krylov_eigs):
    targ_min = min(np.linalg.eig(matrix).eigenvalues)
    err = []
    for i in range(len(matrix)):
        err.append(min(krylov_eigs[i]) - targ_min)
    return err
import matplotlib.pyplot as plt

krylov_error = errors(test_matrix, test_eigs)

plt.plot(krylov_error)
plt.axhline(y=0, color="red", linestyle="--")  # Add dashed red line at y=0
plt.xlabel("Order of Krylov subspace")  # Add x-axis label
plt.ylabel("Error in minimum eigenvalue")  # Add y-axis label
plt.show()

Output:

Output of the previous code cell

クリロフ部分空間が K2\mathcal{K}^2 まで拡大すると、最小固有値がかなり正確に得られることがわかります。また、 K3\mathcal{K}^3 の時点で完全に一致することがわかります。

2.2 行列の次元による時間スケール

クリロフ法が厳密な数値固有値解法よりも有利であることを、以下の方法で納得させよう：

  • ランダム行列を構成する（スパースではないので、KQDの理想的なアプリケーションではない）
  • NumPy を直接使う方法とクリロフ部分空間を使う方法の2つの方法で固有値を決定する。
  • 我々は、クリロフ推定を受け入れる前に、固有値がどの程度正確でなければならないかのカットオフ値を選択する。
  • この2つの方法で解くのに必要な壁時間を比較する。

注意点以下で詳しく説明するように、クリロフ量子対角化は、行列表現が疎であるか、あるいは少数のパウリ作用素群を使って書くことができる作用素に適用するのが最適である。 ここで使っているランダム行列は、そのような説明には当てはまらない。 これらは、古典的なクリロフ法が役に立つかもしれないスケールを探るのに役立つに過ぎない。 第二に、クリロフ法を使う場合、多くの異なるサイズのクリロフ部分空間を使って固有値を計算する。 基底状態固有値の要求精度を達成する最小次元クリロフ部分空間の所要時間を報告する。 繰り返しになるが、これは厳密固有値ソルバーで実行不可能な問題を解くのとは少し違う。なぜなら、必要な次元を評価するために厳密解を使っているからだ。

まず、ランダムな行列のセットを生成する。

import numpy as np

# Set the random seed
np.random.seed(42)

# how many random matrices will we make
num_matrix = 200

matrices = []
for m in range(1, num_matrix):
    matrices.append(np.random.rand(m, m))

ここで、numpyを使って各行列を直接対角化する。 後で比較するために、対角化に要する時間を計算する。

matrix_numpy_times = []
matrix_numpy_eigs = []
for mm in range(num_matrix - 1):
    t0 = time_mus()
    matrix_numpy_eigs.append(min(np.linalg.eig(matrices[mm]).eigenvalues))
    matrix_numpy_times.append(time_mus() - t0)

plt.plot(matrix_numpy_times)
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Output of the previous code cell

上の画像で、125次元付近で時間が異常に長くなっているのは、行列のランダムな性質によるものか、使用された古典的プロセッサの実装によるものかもしれないが、再現性はない。 コードを再実行すると、異なる異常ピークを持つ異なるプロファイルが得られる。

ここで各行列についてクリロフ部分空間を構築し、段階的に固有値を計算する。 各ステップで、指定した絶対誤差の範囲内で最小固有値が得られたかどうかをチェックする。 指定した誤差内で最初に固有値を与える部分空間が、計算時間を記録する部分空間である。 このセルを実行するには、プロセッサーの速度にもよるが、数分かかるかもしれない。 評価を省略したり、対角化する行列の最大次元を小さくしたりすることは自由である。 事前に計算された結果を見るだけで十分だ。

# Choose the absolute error you can tolerate, and make a list for tracking the Krylov subspace size
# at which that error is achieved.
abserr = 0.05
accept_subspace_size = []

# Lists to store total time spent on the Krylov method, and the subset of that time spent on
# diagonalizing the projected matrix.
matrix_krylov_tot_times = []
matrix_krylov_dim = []

# Step through all our random matrices
for mm in range(0, num_matrix - 1):
    test_ks, test_Hs, test_eigs, test_k_tot_times = krylov_full_build(
        np.ones(len(matrices[mm])), matrices[mm]
    )
    # We have not yet found a Krylov subspace that produces our minimum eigenvalue to
    # within the required error.
    found = 0
    for j in range(0, len(matrices[mm]) - 1):
        # If we still haven't found the desired subspace...
        if found == 0:
            # ...but if this one satisfies the requirement, then record everything
            if (
                abs((min(test_eigs[j]) - matrix_numpy_eigs[mm]) / matrix_numpy_eigs[mm])
                < abserr
            ):
                accept_subspace_size.append(j)
                matrix_krylov_tot_times.append(test_k_tot_times[j])
                matrix_krylov_dim.append(mm)
                found = 1

比較のために、この2つの方法で得られたタイムをプロットしてみよう：

plt.plot(matrix_numpy_times, color="blue")
plt.plot(matrix_krylov_dim, matrix_krylov_tot_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Output of the previous code cell

これらは実際に必要な時間であるが、議論のために、隣接するいくつかの点／マトリックス次元の平均によってこれらの曲線を平滑化してみよう。 これは以下のように行われる：

smooth_numpy_times = []
smooth_krylov_times = []

# Choose the number of adjacent points over which to average forward;
# the same will be used backward.
smooth_steps = 10

# We will do this smoothing for all points/matrix dimensions
for i in range(len(matrix_krylov_tot_times)):
    # Ensure we don't exceed the boundaries of our lists
    start = max(0, i - smooth_steps)
    end = min(len(matrix_krylov_tot_times) - 1, i + smooth_steps)

    # Dummy variables for accumulating an average over adjacent points. This is done for both Krylov
    # and the NumPy calculations.
    smooth_count = 0
    smooth_numpy_sum = 0
    smooth_krylov_sum = 0

    for j in range(start, end):
        smooth_numpy_sum = smooth_numpy_sum + matrix_numpy_times[j]
        smooth_krylov_sum = smooth_krylov_sum + matrix_krylov_tot_times[j]
        smooth_count = smooth_count + 1

    # Appending the averaged adjacent values to our new smooth lists
    smooth_numpy_times.append(smooth_numpy_sum / smooth_count)
    smooth_krylov_times.append(smooth_krylov_sum / smooth_count)
plt.plot(smooth_numpy_times, color="blue")
plt.plot(smooth_krylov_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (smoothed, microsec)")  # Add y-axis label
plt.show()

Output:

Output of the previous code cell

クリロフ部分空間の構築に必要な時間は、当初はnumpyの完全対角化に必要な時間を上回ることに注意してください。 しかし、行列のサイズが大きくなると、クリロフ法が有利になる。 これは許容誤差を小さくしても同じで、マトリックスサイズが大きくなるほど有利になる。 これはバラバラにする価値がある。

数値対角化の時間複雑度は O(n3)O(n^3) （アルゴリズムによって多少の違いはある）。 nn ベクトルの正規直交基底を生成する時間の複雑さも O(n3)O(n^3) である。つまり、クリロフ法の利点は、 some\it{some} の正規直交基底の使用とは関係なく、関心のある固有値を効果的に選び出す特定の正規直交基底の使用にある。 このことは、このレッスンの最初のセクションの証明のスケッチからすでに見たとおりであり、これはクリロフ法の収束保証にとって重要である。

これまでの歩みを振り返ってみよう：

  • 非常に大きな行列の場合、クリロフ部分空間法は、従来の対角化アルゴリズムよりも速く、必要な許容範囲内で近似固有値を得ることができる。
  • このような非常に大きな行列では、クリロフ部分空間の生成がクリロフ部分空間法の最も時間のかかる部分である。
  • したがって、クリロフ部分空間を効率的に生成する方法は非常に価値がある。 ここでようやく量子コンピューターが登場する。

理解度チェック

上記の、対角化時間と行列の次元との関係を示す平滑化されたグラフを参照してください。

(a) このグラフによると、クリロフ法はおよそどの行列次元から計算速度が向上したか

(b) クリロフ法がより高速になる閾値は、計算のどのような側面によって変化しうるか

  • (a) 計算を再実行すると結果が異なる場合がありますが、クリロフ法は次元がおよそ80～85になると処理速度が向上します。

    (b) 答えはいくつもあります。 重要な要素としては、求められる精度や、対角化される行列の疎さが挙げられます。

3. 時間発展によるクリロフ

これまで説明してきたことはすべて、古典的に行うことができる。 では、量子コンピューターをいつ、どのように使うのか？ 非常に大きな行列の場合、クリロフ法は長い計算時間と大量のメモリを必要とする。 ∣v⟩|v\rangle、 HH の行列演算に要する時間は、最悪の場合 O(N2)O(N^2) のような規模になる。 ベクトル上の疎行列の乗算（古典的なクリロフ型ソルバーの典型的なケース）でさえ、 O(N)O(N) のような時間複雑度のスケーリングがある。これは、部分空間内のすべてのベクトルに対して行われます。 通常、部分空間の次元 rr は NN の何分の一にもならず、 log⁡(N)\log(N) のようにスケールすることが多い。したがって、すべてのベクトルを生成することは、最悪の場合、 O(N2log⁡(N))O(N^2 \log(N)) のようにスケールする。 直交化のような他のステップもあるが、このスケーリングは覚えておくべき主要なスケーリングである。

量子コンピューティングは、必要な時間とリソースのスケーリングを決定する問題の属性を変えることを可能にする。 行列の大きさ（ NN ）に全面的に依存するのではなく、ショット数やハミルトニアンを構成する非共約パウリ項の数といったものを見ることになる。 この仕組みを探ってみよう。

3.1 時間発展

量子状態を時間発展させる作用素は e−iHt/ℏe^{-iHt/\hbar} （特に量子コンピューティングでは、表記から ℏ\hbar を削除することが非常に一般的）であることを思い出してほしい。 このような演算子の指数関数を理解し、実現する一つの方法は、そのテイラー級数展開を見ることである。 ある初期ベクトル ∣v⟩|v\rangle に作用するこの操作は、初期状態に適用される HH のべき乗が増加する項の和をもたらすことに注意。 初期推測の状態を時間発展させることで、クリロフ部分空間を作ればよさそうだ！

e−iHt/ℏ→e−iHt≈1−iHt−(H2t2)2+⋯e−iHt∣v⟩≈∣v⟩−iHt∣v⟩−(H2t2)2∣v⟩+⋯\begin{aligned} e^{-iHt/\hbar}→e^{-iHt}&≈1-iHt-\frac{(H^2 t^2)}{2}+⋯\\ e^{-iHt} |v\rangle &≈ |v\rangle-iHt|v\rangle-\frac{(H^2 t^2)}{2}|v\rangle+⋯ \end{aligned}

注意点は、実際の量子コンピューターで時間発展を実現することである。 ハミルトニアンの項の多くは互いに交わらない。 したがって、 e−iZe^{-iZ} のような単純な指数演算子は単純な回路に対応するが、一般的なハミルトニアンはそうではない。 また、これらの項には非整合の項が含まれているため、指数では数のように単純な項の積に分解することはできない。

e−iHt=e−i(H1+H2+⋯+Hn)t≠e−iH1te−iH2t...e−iHnte^{-iHt}=e^{-i(H_1+H_2+⋯+H_n)t}\neq e^{-iH_1 t} e^{-iH_2 t}... e^{-iH_n t}

つまり、これは些細なことではないが、量子コンピューティングではよく研究されているプロセスなのだ。 私たちは量子コンピュータ上でトロッター化と呼ばれるプロセスを用いて時間進化を行うが、それ自体が豊かなテーマである [10]。 しかし、非常に高いレベルでは、時間発展を非常に小さなステップ、例えばサイズ dtdt の mm ステップに分割することによって、項の非可換性の影響を制限することができる。

e−iHt=e−i(H1+H2+⋯+Hn)t=(e−i(H1+H2+⋯+Hn)t/m)m≈(e−iH1dte−iH2dt…e−iHndt)me^{-iHt}=e^{-i(H_1+H_2+⋯+H_n )t} = (e^{-i(H_1+H_2+⋯+H_n )t/m} )^m ≈(e^{-iH_1 dt} e^{-iH_2 dt} …e^{-iH_n dt} )^m

ここで dt=t/mdt = t/m。

Hのべき乗を直接用いて古典的な文脈で生成した、次数rのクリロフ部分空間を「べき乗クリロフ部分空間」と呼ぼう。

KPr(H,∣v⟩)=span{∣v⟩,H∣v⟩,H2∣v⟩…Hr−1∣v⟩}\mathcal{K}_P^r (H,|v\rangle)=\text{span}\{|v\rangle,H|v\rangle,H^2 |v\rangle… H^{r-1} |v\rangle\}

ここで、ユニタリー時間発展演算子 U≡e−iHtU \equiv e^{-iHt} を用いて同様の空間を生成する。これを「ユニタリー・クリロフ空間」 KUr\mathcal{K}_U^r と呼ぶことにする。私たちが古典的に用いているパワークリロフ部分空間 KPr\mathcal{K}_P^r は、 HH がユニタリー演算子ではないため、量子コンピュータ上で直接生成することはできません。 ユニタリー・クリロフ部分空間を用いることで、べき乗クリロフ部分空間と同様の収束保証が得られることが示される。すなわち、初期状態 ∣v⟩|v\rangle、指数関数的に消失しない真の基底状態との重なりがあり、固有値間に十分なギャップがある限り、基底状態の誤差は効率的に収束する。 収束についてのより正確な議論は文献 [1] を参照。

ここで、 UU の累乗は異なる時間ステップになる（ UU の累乗 kthk^\text{th} は時間 k×dtk \times dt だけ進む）。 総時間 kdtk dt に対して時間発展する部分空間の要素を、 ∣ψk⟩|\psi_k\rangle とラベル付けすることができる。

U=e−iHdtUk=e−iH(kdt)KUr=span{∣ψ⟩,U∣ψ⟩,U2∣ψ⟩…Ur−1∣ψ⟩}\begin{aligned} U&=e^{-iHdt}\\ U^k&=e^{-iH(kdt)}\\ \mathcal{K}_U^r&=\text{span}\{|\psi\rangle,U|\psi\rangle,U^2 |\psi\rangle… U^{r-1} |\psi\rangle\} \end{aligned}

ハミルトニアンHをユニタリー・クリロフ部分空間 KUr\mathcal{K}_U^r に射影することができる。言い換えれば、 HH の各行列要素を KUr\mathcal{K}_U^r の基底で計算します。 この投影行列を H~\tilde{H} と呼ぶことにする。

3.2 量子コンピュータでの実装方法

H~\tilde{H} の行列要素は期待値 ⟨ψm∣H∣ψn⟩\langle \psi_m |H| \psi_n\rangle で与えられ、これは量子コンピュータを使って推定できる。 HH、異なる量子ビット上のパウリ演算子の和として書くことができ、すべてのパウリ演算子を同時に測定できるわけではないことに留意してほしい。 パウリ項を通約項のグループに分類し、それらを一度に測定することができる。 しかし、すべての条件をカバーするためには、このようなグループがたくさん必要かもしれない。 そこで、項を分割できる別個の通学グループの数、 NGCPN_\text{GCP} が重要になる。

H=∑α=1NGCPcαPαH=\sum_{\alpha=1}^{N_\text{GCP}} c_\alpha P_\alpha

ここで、 PαP_\alpha は、 Pα∼IZIXII...YZXIXP_\alpha \sim IZIXII...YZXIX の形をとるパウリ文字列、あるいは互いに可換なそのようなパウリ文字列の集合である。 HH をこのような可測演算子の和として表すことができることを踏まえると、 H~\tilde{H} の行列要素に関する以下の式は、 IBM Quantum のプリミティブ推定器を用いて実現することができる。

H~mn=⟨ψm∣H∣ψn⟩=⟨ψeiHtm∣H∣ψe−iHtn⟩=⟨ψeiHmdt∣H∣ψe−iHndt⟩\begin{aligned} \tilde{H}_{mn}&=\langle \psi_m |H| \psi_n\rangle\\ &=\langle \psi e^{iHt_m} |H| \psi e^{-iHt_n}\rangle\\ &=\langle \psi e^{iHmdt} |H|\psi e^{-iHndt}\rangle \end{aligned}

ここで、 ∣ψn⟩=e−iHtn∣ψ⟩\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle はユニタリ・クリロフ空間のベクトルであり、 tn=ndtt_n = n dt は、選択された時間ステップ dtdt の倍数である。 量子コンピュータ上では、各行列要素の計算は、量子状態間の重なりを得ることができる任意のアルゴリズムを用いて行うことができます。 このレッスンでは、アダマール検定について重点的に取り上げます。 KU\mathcal{K}_U の次元が rr であることを踏まえると、この部分空間に射影されたハミルトニアンの次元は r×rr \times r となる。 rr が十分に小さい場合（一般に、固有値の推定値が収束するためには r<<100r<<100 で十分である）、古典的には、射影されたハミルトニアン H~\tilde{H} を容易に対角化することができる。 しかし、クリロフ空間のベクトルが直交していないため、 H~\tilde{H} を直接対角化することはできない。 それらの重なりを測定し、行列を作成する必要があります。 S~\tilde{S}

S~mn=⟨ψm∣ψn⟩\tilde{S}_{mn} = \langle \psi_m \vert \psi_n \rangle

これにより、非直交空間における固有値問題（一般化固有値問題とも呼ばれる）を解くことができる

H~ c⃗=E S~ c⃗\tilde{H} \ \vec{c} = E \ \tilde{S} \ \vec{c}

そして、この一般化された固有値問題の解を見ることによって、 HH の固有値と固有状態の推定値を得ることができる。 例えば、基底状態のエネルギーの推定値は、最小の固有値 EE と、対応する固有ベクトル c⃗\vec{c} から基底状態を取ることによって得られる。 c⃗\vec{c} の係数は、 KU\mathcal{K}_U にまたがるさまざまなベクトルの寄与を決定します。

一般化された固有値問題

なぜ単純に H~\tilde{H} を対角化できないのでしょうか？ S~\tilde{S}、クリロフ基底（非常に特殊な場合を除き、すべて非直交である）の幾何学的形状に関する情報が含まれているので、 H~\tilde{H}、それだけでは完全なハミルトニアンの射影を記述しないので、その固有値は完全なハミルトニアンの固有値とは特に関係がありません。 一般化固有値問題を解くことは、完全なハミルトニアンのクリロフ空間への射影に対応する近似固有値と固有ベクトルを得るために必要である。

修正ハダマード・テストを実行するためには、回路を異なる状態で何度も使用しなければならないことを示す、多くの層を持つ回路図。

図に示すのは、異なる量子状態間の重なりを計算するために用いられる修正ハダマールテストの回路図である。 各行列要素 H~i,j\tilde{H}_{i,j}、状態 ∣ψi⟩\vert \psi_i \rangle、 ∣ψj⟩\vert \psi_j \rangle の間のハダマルド検定が行われる。 図では、行列要素と対応する Prep  ψi\text{Prep} \; \psi_i、 Prep  ψj\text{Prep} \; \psi_j の操作の配色によって、このことが強調されている。 したがって、投影ハミルトニアン H~\tilde{H} のすべての行列要素を計算するために、クリロフ空間ベクトルのすべての可能な組み合わせに対するハダマルド検定のセットが必要となる。ハダマルド検定回路の一番上のワイヤーはアンシラ量子ビットで、XまたはY基底で測定され、その期待値は状態間の重なりの値を決定します。 一番下のワイヤーは、システム・ハミルトニアンのすべての量子ビットを表している。 Prep  ψi\text{Prep} \; \psi_i ∣ψi⟩\vert \psi_i \rangle （ ）、 （ ）はシステム・ハミルトニアンのパウリ分解を表す。量子コンピュータでの実装を以下に詳しく示します。 Prep  ψj\text{Prep} \; \psi_j PP H=∑iPiH = \sum_i P_i

4. 量子コンピュータ上でのクリロフ量子対角化

次に、クリロフ量子対角化を実際の量子コンピュータで実装する。 まずは便利なパッケージをいくつかインポートしてみよう。

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import warnings

from qiskit.quantum_info import SparsePauliOp, Pauli
from qiskit.circuit import Parameter
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

# from qiskit.providers.fake_provider import Fake20QV1
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2 as Estimator, Batch

import itertools as it

warnings.filterwarnings("ignore")

先ほど説明した一般化固有値問題を解く関数を以下に定義する。

def solve_regularized_gen_eig(
    h: np.ndarray,
    s: np.ndarray,
    threshold: float,
    k: int = 1,
    return_dimn: bool = False,
) -> Union[float, List[float]]:
    """
    Method for solving the generalized eigenvalue problem with regularization

    Args:
        h (numpy.ndarray):
            The effective representation of the matrix in our Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of our Krylov subspace
        threshold (float):
            Cut-off value for the eigenvalue of s
        k (int):
            Number of eigenvalues to return
        return_dimn (bool):
            Whether to return the size of the regularized subspace

    Returns:
        lowest k-eigenvalue(s) that are the solution of the regularized generalized eigenvalue problem


    """
    s_vals, s_vecs = sp.linalg.eigh(s)
    s_vecs = s_vecs.T
    good_vecs = np.array([vec for val, vec in zip(s_vals, s_vecs) if val > threshold])
    h_reg = good_vecs.conj() @ h @ good_vecs.T
    s_reg = good_vecs.conj() @ s @ good_vecs.T
    if k == 1:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][0]
    else:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k]

少なくとも最初のベンチマークでは、収束挙動をチェックするために、厳密な古典解を知っておくことは有用である。 以下の関数は、ハミルトニアンと量子ビット数を引数として、ハミルトニアンの基底状態エネルギーを計算します。

def single_particle_gs(H_op, n_qubits):
    """
    Find the ground state of the single particle(excitation) sector
    """
    H_x = []
    for p, coeff in H_op.to_list():
        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))

    H_z = []
    for p, coeff in H_op.to_list():
        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))

    H_c = H_op.coeffs

    print("n_sys_qubits", n_qubits)

    n_exc = 1
    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))
    print("n_exc", n_exc, ", subspace dimension", sub_dimn)

    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)

    sparse_vecs = [
        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)
    ]  # list all of the possible sets of n_exc indices of 1s in n_exc-particle states

    m = 0
    for i, i_set in enumerate(sparse_vecs):
        for j, j_set in enumerate(sparse_vecs):
            m += 1

            if len(i_set.symmetric_difference(j_set)) <= 2:
                for p_x, p_z, coeff in zip(H_x, H_z, H_c):
                    if i_set.symmetric_difference(j_set) == p_x:
                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (
                            (-1) ** len(i_set.intersection(p_z))
                        )
                    else:
                        sgn = 0

                    few_particle_H[i, j] += sgn * coeff

    gs_en = min(np.linalg.eigvalsh(few_particle_H))
    print("single particle ground state energy: ", gs_en)
    return gs_en

4.1 ステップ1: 問題を量子回路と演算子にマッピングする

ここでハミルトニアンを定義する。 これは上の関数とは異なり、上の関数はハミルトニアンを引数にとり、基底状態のみを返す。 ここで定義するハミルトニアンは、すべてのエネルギー固有状態のエネルギー準位を決定するものであり、このハミルトニアンはパウリ作用素を用いて構成することができ、量子コンピューター上で実装することができる。

我々は、「ハイゼンベルグ鎖」と呼ばれる、空間上でどのような向きにもなりうるスピンの鎖に対応するハミルトニアンを選ぶ。 我々は、 ithi^\text{th} のスピンは、最も近い隣人（ (i−1)th(i-1)^\text{th} と (i+1)th(i+1)^\text{th} のスピン）の影響を受けるが、より遠い隣人の影響を受けることはないと仮定している。 また、スピンが異なる軸を向いている場合、スピン間の相互作用が異なる可能性も考慮する。 この非対称性は、例えばスピンが埋め込まれた結晶格子の構造によって生じることがある。

# Define problem Hamiltonian.
n_qubits = 10
# coupling strength for XX, YY, and ZZ interactions
JX = 1
JY = 3
JZ = 2

# Define the Hamiltonian:
H_int = [["I"] * n_qubits for _ in range(3 * (n_qubits - 1))]
for i in range(n_qubits - 1):
    H_int[i][i] = "Z"
    H_int[i][i + 1] = "Z"
for i in range(n_qubits - 1):
    H_int[n_qubits - 1 + i][i] = "X"
    H_int[n_qubits - 1 + i][i + 1] = "X"
for i in range(n_qubits - 1):
    H_int[2 * (n_qubits - 1) + i][i] = "Y"
    H_int[2 * (n_qubits - 1) + i][i + 1] = "Y"
H_int = ["".join(term) for term in H_int]
H_tot = [
    (term, JZ)
    if term.count("Z") == 2
    else (term, JY)
    if term.count("Y") == 2
    else (term, JX)
    for term in H_int
]

# Get operator
H_op = SparsePauliOp.from_list(H_tot)
print(H_tot)

Output:

[('ZZIIIIIIII', 2), ('IZZIIIIIII', 2), ('IIZZIIIIII', 2), ('IIIZZIIIII', 2), ('IIIIZZIIII', 2), ('IIIIIZZIII', 2), ('IIIIIIZZII', 2), ('IIIIIIIZZI', 2), ('IIIIIIIIZZ', 2), ('XXIIIIIIII', 1), ('IXXIIIIIII', 1), ('IIXXIIIIII', 1), ('IIIXXIIIII', 1), ('IIIIXXIIII', 1), ('IIIIIXXIII', 1), ('IIIIIIXXII', 1), ('IIIIIIIXXI', 1), ('IIIIIIIIXX', 1), ('YYIIIIIIII', 3), ('IYYIIIIIII', 3), ('IIYYIIIIII', 3), ('IIIYYIIIII', 3), ('IIIIYYIIII', 3), ('IIIIIYYIII', 3), ('IIIIIIYYII', 3), ('IIIIIIIYYI', 3), ('IIIIIIIIYY', 3)]

以下のコードでは、ハミルトニアンを単一粒子状態に制限し、スペクトルノルムを用いて時間ステップ dtdt の適切なサイズを設定しています。ハミルトニアンノルムの上限に基づいて）時間ステップ dt の値を発見的に選択します。 文献 [9] は、十分に小さいタイムス テップは π/∣∣H∣∣\pi/\vert \vert H \vert \vert、ある点までは、この値を過大評価するよりも過小評価する方が望ましいことを示した。過大評価すると、高エネルギー状態からの寄与がクリロフ空間の最適状態さえも破壊してしまうからである。 一方、 dtdt を小さくしすぎると、クリロフ基底ベクトルのタイムステップごとの違いが小さくなるため、クリロフ部分空間のコンディショニングが悪くなる。

# Get Hamiltonian restricted to single-particle states
single_particle_H = np.zeros((n_qubits, n_qubits))
for i in range(n_qubits):
    for j in range(i + 1):
        for p, coeff in H_op.to_list():
            p_x = Pauli(p).x
            p_z = Pauli(p).z
            if all(p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)):
                sgn = ((-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))) * (
                    (-1) ** p_z[i]
                )
            else:
                sgn = 0
            single_particle_H[i, j] += sgn * coeff
for i in range(n_qubits):
    for j in range(i + 1, n_qubits):
        single_particle_H[i, j] = np.conj(single_particle_H[j, i])

# Set dt according to spectral norm
dt = np.pi / np.linalg.norm(single_particle_H, ord=2)
dt

Output:

np.float64(0.17453292519943295)

時間発展に使用するトロッターステップの数を指定します。 また、最大クリロフ次元を4と指定します。 このクリロフ次元は、現実的な応用には十分ではない。 しかし、この例にはこれで十分である。 さらに、より小さな次元での収束性を検証する。 後のレッスンでは、ハミルトニアンをより大きな部分空間にスケーリングおよび射影する方法を探求します。

# Set parameters for quantum Krylov algorithm
krylov_dim = 4  # size of krylov subspace
num_trotter_steps = 4
dt_circ = dt / num_trotter_steps

状態準備

基底状態と重なる参照状態（ ∣ψ⟩\vert \psi \rangle ）を選ぶ。 このハミルトニアンでは、中央の量子ビット（ ∣00..010...00⟩\vert 00..010...00 \rangle ）が励起された状態を参照状態として用いる。

qc_state_prep = QuantumCircuit(n_qubits)
qc_state_prep.x(int(n_qubits / 2) + 1)
qc_state_prep.draw("mpl", scale=0.5)

Output:

Output of the previous code cell

時間発展

我々は、与えられたハミルトニアンによって生成される時間発展演算子を実現することができる。 U=e−iHtU=e^{-iHt}、 リー・トロッター近似を経由する。 簡単のため、時間発展回路に内蔵されている PauliEvolutionGate 。 一般的な構文はこうだ。

t = Parameter("t")

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)
)

qr = QuantumRegister(n_qubits)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

Output:

<qiskit.circuit.instructionset.InstructionSet at 0x7ccaa4664250>

以下、ハダマードテストでこのバージョンを使用するが、時間 dtdt に対して前進する。

ハダマール検定

ここで、ハダマール検定を用いて、 H~\tilde{H} およびグラム行列 S~\tilde{S} の両方の行列要素を計算したいことを思い出してください。 この文脈において、これがどのように機能するのかを振り返ってみましょう。まずは、 H~\tilde{H} の構成に焦点を当てます。全体的なプロセスを以下の図に示します。 色分けされた状態準備ブロックの層（ Prep∣ψi⟩\text{Prep}|\psi_i\rangle ）は、この処理が、我々の部分空間における ∣ψi⟩|\psi_i\rangle および ∣ψj⟩|\psi_j\rangle のすべての組み合わせに対して行われることを示唆している。

ハダマルド・テストを行うためには、多くの異なる状態について回路を評価する必要があることを示す、多くの層を持つ量子回路図のイメージ。

各ステップにおけるシステムの状態は以下の通りである：

Step 0:∣Ψ⟩=∣0⟩∣0⟩NStep 1:∣Ψ⟩=12(∣0⟩+∣1⟩)∣0⟩NStep 2:∣Ψ⟩=12(∣0⟩∣0⟩N+∣1⟩∣ψi⟩)Step 3:∣Ψ⟩=12(∣0⟩∣0⟩N+∣1⟩P∣ψi⟩)Step 4:∣Ψ⟩=12(∣0⟩∣ψj⟩+∣1⟩P∣ψi⟩)\begin{aligned} \text{Step 0:}\qquad|\Psi\rangle & = |0\rangle|0\rangle^N \\ \text{Step 1:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle + |1\rangle \Big)|0\rangle^N \\ \text{Step 2:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle|0\rangle^N+|1\rangle |\psi_i\rangle\Big)\\ \text{Step 3:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle |0\rangle^N+|1\rangle P |\psi_i\rangle\Big) \\ \text{Step 4:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\Big(|0\rangle |\psi_j\rangle+|1\rangle P|\psi_i\rangle\Big) \end{aligned}

ここで、 PP はハミルトニアンの分解におけるパウリ項である（ユニタリーでなくなるため、複数のパウリ項の線形結合にはなりえないことに注意--後で示す別の構成を使ってグループ化は可能）。 Prep  ψi\text{Prep} \; \psi_i、 Prep  ψj\text{Prep} \; \psi_j は、ユニタリー・クリロフ空間のベクトル ∣ψi⟩|\psi_i\rangle、 ∣ψj⟩|\psi_j\rangle を準備する制御された演算であり、 ∣ψk⟩=e−iHkdt∣ψ⟩=e−iHkdtUψ∣0⟩N|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N。この回路に XX と YY の測定を適用すると、必要な行列要素の実部と虚部がそれぞれ計算される。

上記のステップ4から始めて、0番目の量子ビットにハダマードゲート HH。

∣Ψ⟩⟶12∣0⟩(∣ψj⟩+P∣ψi⟩)+12∣1⟩(∣ψj⟩−P∣ψi⟩)\begin{equation*} |\Psi\rangle \longrightarrow\quad\frac{1}{2}|0\rangle\Big( |\psi_j\rangle + P|\psi_i\rangle\Big) + \frac{1}{2}|1\rangle\Big(|\psi_j\rangle - P|\psi_i\rangle\Big) \end{equation*}

次に、 XX または YY を測定する。

⇒⟨X⟩=14(∥∣ψj⟩+P∣ψi⟩∥2−∥∣ψj⟩−P∣ψi⟩∥2)=Re[⟨ψj∣P∣ψi⟩].\begin{equation*} \begin{split} \Rightarrow\quad\langle X\rangle &= \frac{1}{4}\Bigg(\Big\|| \psi_j\rangle + P|\psi_i\rangle \Big\|^2-\Big\||\psi_j\rangle - P|\psi_i\rangle\Big\|^2\Bigg) \\ &= \text{Re}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{split} \end{equation*}

恒等式から ∣a+b∥2=⟨a+b∣a+b⟩=∥a∥2+∥b∥2+2Re⟨a∣b⟩|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle。同様に、 YY を測定すると、次のようになる

⟨Y⟩=Im[⟨ψj∣P∣ψi⟩].\begin{equation*} \langle Y\rangle = \text{Im}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{equation*}

先に設定した時間進化にこれらのステップを加えると、次のようになる。

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)

## Create the time-evo op dagger circuit
evol_gate_d = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)
evol_gate_d = evol_gate_d.inverse()

# Put pieces together
qc_reg = QuantumRegister(n_qubits)
qc_temp = QuantumCircuit(qc_reg)
qc_temp.compose(qc_state_prep, inplace=True)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate, qargs=qc_reg)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate_d, qargs=qc_reg)
qc_temp.compose(qc_state_prep.inverse(), inplace=True)

# Create controlled version of the circuit
controlled_U = qc_temp.to_gate().control(1)

# Create hadamard test circuit for real part
qr = QuantumRegister(n_qubits + 1)
qc_real = QuantumCircuit(qr)
qc_real.h(0)
qc_real.append(controlled_U, list(range(n_qubits + 1)))
qc_real.h(0)

print("Circuit for calculating the real part of the overlap in S via Hadamard test")
qc_real.draw("mpl", fold=-1, scale=0.5)

Output:

Circuit for calculating the real part of the overlap in S via Hadamard test
Output of the previous code cell

トロッター・サーキットの奥深さについてはすでに警告した。 このような条件でハダマルドテストを実行すると、特にネイティブゲートに分解すると、さらに深い回路が得られる。 もしデバイスのトポロジーを考慮すれば、これはさらに増えるだろう。 従って、量子コンピュータを使う前に、回路の2量子ビットの深さをチェックするのは良い考えだ。

print(
    "Number of layers of 2Q operations",
    qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),
)

Output:

Number of layers of 2Q operations 14401

これほど複雑な回路では、現代の量子コンピュータでは実用的な結果を得ることができません。 H~\tilde{H} と S~\tilde{S} を構築するには、もっと良い方法が必要です。 これが、後述する効率的なアダマール検定が導入された理由である。

4. 2 ステップ2. 対象ハードウェア向けに回路と演算子を最適化する

効率的なハダマール検定

いくつかの近似を導入し、モデルのハミルトニアンに関するいくつかの仮定に依存することで、得られたハダマルド・テストの深層回路を最適化することができる。 例えば、次のようなハダマード・テストの回路を考えてみよう：

効率的なハダマルドテストを行うためには、多くの異なるユニタリー演算子に対して回路を評価する必要があることを示す。

ハミルトニアン HH の下での ∣0⟩N|0\rangle^N の固有値 E0E_0 を古典的に計算できると仮定する。 これはハミルトニアンがU(1)対称性を保存しているときに満たされる。 これは強い仮定のように思えるかもしれないが、ハミルトニアンの作用に影響されない真空状態（この場合、 ∣0⟩N|0\rangle^N 状態に対応する）が存在すると仮定しても安全な場合がたくさんある。 これは例えば、安定な分子（電子の数が保存されている）を記述する化学のハミルトニアンに当てはまります。 ゲート Prep  ψ0\text{Prep} \; \psi_0、所望の参照状態 ∣ψ0⟩=Prep  ψ0∣0⟩=e−iH0dtUψ0∣0⟩\ket{\psi_0} = \text{Prep} \; \psi_0 \ket{0} = e^{-i H 0 dt} U_{\psi_0} \ket{0} を準備することを考えると、例えば、化学のためのHF状態 Prep  ψ0\text{Prep} \; \psi_0 を準備することは、単一量子ビットNOTの積になるので、制御された- Prep  ψ0\text{Prep} \; \psi_0 は単なるCNOTの積である。 すると、上記の回路は、測定前の以下の状態を実現する：

Step 0:∣Ψ⟩=∣0⟩∣0⟩NStep 1:∣Ψ⟩=12(∣0⟩∣0⟩N+∣1⟩∣0⟩N)Step 2:∣Ψ⟩=12(∣0⟩∣0⟩N+∣1⟩∣ψ0⟩)Step 3:∣Ψ⟩=12(eiϕ∣0⟩∣0⟩N+∣1⟩U∣ψ0⟩)Step 4:∣Ψ⟩=12(eiϕ∣0⟩∣ψ0⟩+∣1⟩U∣ψ0⟩)=12(∣+⟩(eiϕ∣ψ0⟩+U∣ψ0⟩)+∣−⟩(eiϕ∣ψ0⟩−U∣ψ0⟩))=12(∣+i⟩(eiϕ∣ψ0⟩−iU∣ψ0⟩)+∣−i⟩(eiϕ∣ψ0⟩+iU∣ψ0⟩))\begin{aligned} \text{Step 0:}\qquad|\Psi\rangle & = \ket{0} \ket{0}^{N}\\ \text{Step 1:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(\ket{0}\ket{0}^N+ \ket{1} \ket{0}^N\right)\\ \text{Step 2:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(|0\rangle|0\rangle^N+|1\rangle|\psi_0\rangle\right)\\ \text{Step 3:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0}\ket{0}^N+\ket{1} U\ket{\psi_0}\right)\\ \text{Step 4:}\qquad|\Psi\rangle & = \frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0} \ket{\psi_0}+\ket{1} U\ket{\psi_0}\right)\\ & = \frac{1}{2}\left(\ket{+}\left(e^{i\phi}\ket{\psi_0}+U\ket{\psi_0}\right)+\ket{-}\left(e^{i\phi}\ket{\psi_0}-U\ket{\psi_0}\right)\right)\\ & = \frac{1}{2}\left(\ket{+i}\left(e^{i\phi}\ket{\psi_0}-iU\ket{\psi_0}\right)+\ket{-i}\left(e^{i\phi}\ket{\psi_0}+iU\ket{\psi_0}\right)\right) \end{aligned}

ここでは、ステップ2からステップ3まで、古典的なシミュラブル位相シフト U∣0⟩N=eiϕ∣0⟩N U\ket{0}^N = e^{i\phi}\ket{0}^N。 したがって期待値は

⟨X⊗P⟩=14((e−iϕ⟨ψ0∣+⟨ψ0∣U†)P(eiϕ∣ψ0⟩+U∣ψ0⟩)−(e−iϕ⟨ψ0∣−⟨ψ0∣U†)P(eiϕ∣ψ0⟩−U∣ψ0⟩))=Re[e−iϕ⟨ψ0∣PU∣ψ0⟩],\begin{aligned} \langle X\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi_0}+\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi}\ket{\psi_0}+U\ket{\psi_0}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi_0}-\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi}\ket{\psi_0}-U\ket{\psi_0}\right) \Big)\\ &=\text{Re}\left[e^{-i\phi}\bra{\psi_0}PU\ket{\psi_0}\right], \end{aligned} ⟨Y⊗P⟩=14((e−iϕ⟨ψ0∣+i⟨ψ0∣U†)P(eiϕ0∣ψ0⟩−iU∣ψ0⟩)−(e−iϕ⟨ψ0∣−i⟨ψ0∣U†)P(eiϕ∣ψ0⟩+iU∣ψ0⟩))=Im[e−iϕ⟨ψ0∣PU∣ψ0⟩]. \begin{aligned} \langle Y\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi_0}+i\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi_0}\ket{\psi_0}-iU\ket{\psi_0}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi_0}-i\bra{\psi_0}U^\dagger\right)P\left(e^{i\phi}\ket{\psi_0}+iU\ket{\psi_0}\right) \Big)\\ &=\text{Im}\left[e^{-i\phi}\bra{\psi_0}PU\ket{\psi_0}\right]. \end{aligned}

これらの仮定を用いることで、より少ない制御演算で、関心のある演算子の期待値を書くことができた。 実際には、制御された状態準備（ Prep  ψ0\text{Prep} \; \psi_0 ）だけを実装すればよく、制御された時間発展は必要ない。 上記のように計算を組み直すことで、結果として得られる回路の深さを大幅に減らすことができる。

ボーナスとして、パウリ演算子は回路の途中で制御されたゲートとしてではなく、回路の最後で測定されるようになったので、上で示した分解 H=∑α=1NGCPcαPαH=\sum_{\alpha = 1}^{N_\text{GCP}}c_\alpha P_\alpha のように、他の通学パウリ演算子と一緒に測定することができる。

トロッター分解を用いて時間発展演算子を分解する

時間発展演算子を厳密に実装する代わりに、トロッター分解を用いてその近似を実装することができます。 ある順序のトロッター分解を数回繰り返すことで、近似によって生じる誤差をさらに低減することができます。 以下では、検討しているハミルトニアンの相互作用グラフ（最近接相互作用のみ）に対して、最も効率的な方法でトロッター法の実装を直接構築する。 実際には、結合強度 JxJ_x、 JyJ_y、 JzJ_z およびパラメータ化された角度 tt を持つパウリ回転 RxxR_{xx}、 RyyR_{yy}、 RzzR_{zz} を挿入します。これらは、 e−i(JxXX+JyYY+JzZZ)te^{-i (J_x XX + J_y YY + J_z ZZ) t} の近似的な実装に対応しています。パウリ回転の定義と、実装しようとしている時間発展との違いを考慮すると、 dtdt の時間発展を実現するには、パラメータ 2∗dt2*dt を使用する必要があります。 さらに、トロッターステップの繰り返し回数が奇数の場合、演算の順序を逆転させます。これは機能的には同等ですが、隣接する演算を単一の SU(2)SU(2) ユニタリー演算で合成できるようになります。 これにより、汎用機能 PauliEvolutionGate() を使用した場合よりもはるかに浅い回路が得られます。

t = Parameter("t")

# Create instruction for rotation about XX+YY-ZZ:
Rxyz_circ = QuantumCircuit(2)
Rxyz_circ.rxx(2 * JX * t, 0, 1)
Rxyz_circ.ryy(2 * JY * t, 0, 1)
Rxyz_circ.rzz(2 * JZ * t, 0, 1)
Rxyz_instr = Rxyz_circ.to_instruction(label="R J_x XX + J_y YY + J_z ZZ")

interaction_list = [
    [[i, i + 1] for i in range(0, n_qubits - 1, 2)],
    [[i, i + 1] for i in range(1, n_qubits - 1, 2)],
]  # linear chain

qr = QuantumRegister(n_qubits)
trotter_step_circ = QuantumCircuit(qr)
for i, color in enumerate(interaction_list):
    for interaction in color:
        trotter_step_circ.append(Rxyz_instr, interaction)
    if i < len(interaction_list) - 1:
        trotter_step_circ.barrier()
reverse_trotter_step_circ = trotter_step_circ.reverse_ops()

qc_evol = QuantumCircuit(qr)
for step in range(num_trotter_steps):
    if step % 2 == 0:
        qc_evol = qc_evol.compose(trotter_step_circ)
    else:
        qc_evol = qc_evol.compose(reverse_trotter_step_circ)

qc_evol.decompose().draw("mpl", fold=-1, scale=0.5)

Output:

Output of the previous code cell

この効率的なハダマルド・テストのために、再び初期状態を準備する。

control = 0
excitation = int(n_qubits / 2) + 1
controlled_state_prep = QuantumCircuit(n_qubits + 1)
controlled_state_prep.cx(control, excitation)
controlled_state_prep.draw("mpl", fold=-1, scale=0.5)

Output:

Output of the previous code cell

S~\tilde{S} および H~\tilde{H} の行列要素をHadamardテストを用いて計算するためのテンプレート回路

ハダマード・テストで使用される回路の唯一の違いは、時間発展演算子の位相と測定される観測値である。 したがって、時間発展演算子に依存するゲートのプレースホルダーを持つ、ハダマード・テストの一般的な回路を表すテンプレート回路を用意することができる。

# Parameters for the template circuits
parameters = []
for idx in range(1, krylov_dim):
    parameters.append(dt_circ * (idx))
# Create modified hadamard test circuit
qr = QuantumRegister(n_qubits + 1)
qc = QuantumCircuit(qr)
qc.h(0)
qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)
qc.barrier()
qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)
qc.barrier()
qc.x(0)
qc.compose(controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True)
qc.x(0)

qc.decompose().draw("mpl", fold=-1)

Output:

Output of the previous code cell 
print(
    "The optimized circuit has 2Q gates depth: ",
    qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),
)

Output:

The optimized circuit has 2Q gates depth:  50

この深さは、オリジナルのハダマード・テストに比べて大幅に減少している。 この深さは、最新の量子コンピューターでも管理可能なものだが、それでもかなり高い。 有用な結果を得るためには、最先端のエラー緩和を用いる必要がある。

量子クリロフ計算を実行するバックエンドを選択し、その量子コンピューター上で実行できるように回路をトランスパイルする。

# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)

# Or you may choose a specify backend and channel if necessary for your workflow.
# service = QiskitRuntimeService(channel="ibm_quantum_platform")
# backend = service.backend("ibm_fez")

ここで回路と演算子をトランスパイルする。

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

target = backend.target
basis_gates = list(target.operation_names)
pm = generate_preset_pass_manager(
    optimization_level=3, backend=backend, basis_gates=basis_gates
)

qc_trans = pm.run(qc)
print(qc_trans.depth(lambda x: x[0].num_qubits == 2))
print(qc_trans.count_ops())
qc_trans.draw("mpl", fold=-1, idle_wires=False, scale=0.5)

Output:

36
OrderedDict([('rz', 410), ('sx', 361), ('cz', 156), ('x', 18), ('barrier', 6)])
Output of the previous code cell

最適化の後、我々のトランスパイルされた2量子ビットの深さはさらに減少した。

4.3 ステップ3。 IBM Quantum プリミティブを使用して実行する

次に、Estimatorで実行するためのPUBを作成する。

# Define observables to measure for S
observable_S_real = "I" * (n_qubits) + "X"
observable_S_imag = "I" * (n_qubits) + "Y"

observable_op_real = SparsePauliOp(
    observable_S_real
)  # define a sparse pauli operator for the observable
observable_op_imag = SparsePauliOp(observable_S_imag)

layout = qc_trans.layout  # get layout of transpiled circuit
observable_op_real = observable_op_real.apply_layout(
    layout
)  # apply physical layout to the observable
observable_op_imag = observable_op_imag.apply_layout(layout)
observable_S_real = (
    observable_op_real.paulis.to_labels()
)  # get the label of the physical observable
observable_S_imag = observable_op_imag.paulis.to_labels()

observables_S = [[observable_S_real], [observable_S_imag]]


# Define observables to measure for H
# Hamiltonian terms to measure
observable_list = []
for pauli, coeff in zip(H_op.paulis, H_op.coeffs):
    # print(pauli)
    observable_H_real = pauli[::-1].to_label() + "X"
    observable_H_imag = pauli[::-1].to_label() + "Y"
    observable_list.append([observable_H_real])
    observable_list.append([observable_H_imag])

layout = qc_trans.layout

observable_trans_list = []
for observable in observable_list:
    observable_op = SparsePauliOp(observable)
    observable_op = observable_op.apply_layout(layout)
    observable_trans_list.append([observable_op.paulis.to_labels()])

observables_H = observable_trans_list


# Define a sweep over parameter values
params = np.vstack(parameters).T


# Estimate the expectation value for all combinations of
# observables and parameter values, where the pub result will have
# shape (# observables, # parameter values).
pub = (qc_trans, observables_S + observables_H, params)

t=0t=0 の回路は古典的に計算可能である。 量子コンピューターを使った t≠0t\neq 0。

from qiskit.quantum_info import StabilizerState, Pauli


qc_cliff = qc.assign_parameters({t: 0})


# Get expectation values from experiment
S_expval_real = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "X")
)
S_expval_imag = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "Y")
)

# Get expectation values
S_expval = S_expval_real + 1j * S_expval_imag

H_expval = 0
for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Get expectation values from experiment
    expval_real = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "X")
    )
    expval_imag = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "Y")
    )
    expval = expval_real + 1j * expval_imag

    # Fill-in matrix elements
    H_expval += coeff * expval


print(H_expval)

Output:

(10+0j)

効率的なハダマード・テストを用いることで、ゲートの深さを1桁減らすことができたが、それでも最先端のエラー緩和を必要とするには十分な深さである。 以下では、使用されるミティゲーションの属性を指定する。 使用される手法はどれも重要だが、特に確率的エラー増幅（PEA ）について言及する価値がある。 この強力なテクニックには、大量の量子オーバーヘッドが伴う。 ここで行われた計算は、実際の量子コンピューターで実行すると20分以上かかる。 精度とオーバーヘッドを増減させるために、以下のパラメータを調整するとよいだろう。 以下のデフォルト設定では、忠実度の高い結果が得られる。

# Experiment options
num_randomizations = 300
num_randomizations_learning = 20
max_batch_circuits = 20
shots_per_randomization = 100
learning_pair_depths = [0, 4, 24]
noise_factors = [1, 1.3, 1.6]

# Base option formatting
options = {
    # Builtin resilience settings for ZNE
    "resilience": {
        "measure_mitigation": True,
        "zne_mitigation": True,
        "zne": {"noise_factors": noise_factors},
        # TREX noise learning configuration
        "measure_noise_learning": {
            "num_randomizations": num_randomizations_learning,
            "shots_per_randomization": shots_per_randomization,
        },
        # PEA noise model configuration
        "layer_noise_learning": {
            "max_layers_to_learn": 10,
            "layer_pair_depths": learning_pair_depths,
            "shots_per_randomization": shots_per_randomization,
            "num_randomizations": num_randomizations_learning,
        },
    },
    # Randomization configuration
    "twirling": {
        "num_randomizations": num_randomizations,
        "shots_per_randomization": shots_per_randomization,
        "strategy": "all",
    },
    # Experimental settings for PEA method
    "experimental": {
        # # Just in case, disable any further qiskit transpilation not related to twirling / DD
        # "skip_transpilation": True,
        # Execution configuration
        "execution": {
            "max_pubs_per_batch_job": max_batch_circuits,
            "fast_parametric_update": True,
        },
        # Error Mitigation configuration
        "resilience": {
            # ZNE Configuration
            "zne": {
                "amplifier": "pea",
                "return_all_extrapolated": True,
                "return_unextrapolated": True,
                "extrapolated_noise_factors": [0] + noise_factors,
            }
        },
    },
}

最後に、 S~\tilde{S} と H~\tilde{H} の回路をEstimatorで実行する。

# This job required 17 minutes of QPU time to run on a Heron r2 processor. This is only an estimate.
# Your execution time may vary.

with Batch(backend=backend) as batch:
    # Estimator
    estimator = Estimator(mode=batch, options=options)

    job = estimator.run([pub], precision=1)

4.4 ステップ4. 結果の後処理と分析

私たちが量子コンピュータから得たものは、 S~\tilde{S} の個々の行列要素と、 H~\tilde{H} の行列要素を構成する通約パウリ群である。これらの項を組み合わせて行列を復元し、一般化固有値問題を解くことができるようにしなければならない。

# Store the outputs as 'results'.
results = job.result()[0]

有効ハミルトニアンとオーバーラップ行列を計算する

まず、制御されていない時間発展中に ∣0⟩\vert 0 \rangle 状態によって蓄積された位相を計算する

prefactors = [
    np.exp(-1j * sum([c for p, c in H_op.to_list() if "Z" in p]) * i * dt)
    for i in range(1, krylov_dim)
]

回路の実行結果が得られたら、データを後処理して、以下の行列要素を計算することができる。 SS

# Assemble S, the overlap matrix of dimension D:
S_first_row = np.zeros(krylov_dim, dtype=complex)
S_first_row[0] = 1 + 0j

# Add in ancilla-only measurements:
for i in range(krylov_dim - 1):
    # Get expectation values from experiment
    expval_real = results.data.evs[0][0][i]  # automatic extrapolated evs if ZNE is used
    expval_imag = results.data.evs[1][0][i]  # automatic extrapolated evs if ZNE is used

    # Get expectation values
    expval = expval_real + 1j * expval_imag
    S_first_row[i + 1] += prefactors[i] * expval

S_first_row_list = S_first_row.tolist()  # for saving purposes


S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        S_circ[j, i] = S_first_row[i - j]
    else:
        S_circ[j, i] = np.conj(S_first_row[j - i])
from sympy import Matrix

Matrix(S_circ)

Output:

[1.00.149322296177984−0.283023058106896i0.185815978760175−0.0910521940394691i0.0940509850777074−0.094154537369141i0.149322296177984+0.283023058106896i1.00.149322296177984−0.283023058106896i0.185815978760175−0.0910521940394691i0.185815978760175+0.0910521940394691i0.149322296177984+0.283023058106896i1.00.149322296177984−0.283023058106896i0.0940509850777074+0.094154537369141i0.185815978760175+0.0910521940394691i0.149322296177984+0.283023058106896i1.0]\displaystyle \left[\begin{matrix}1.0 & 0.149322296177984 - 0.283023058106896 i & 0.185815978760175 - 0.0910521940394691 i & 0.0940509850777074 - 0.094154537369141 i\\0.149322296177984 + 0.283023058106896 i & 1.0 & 0.149322296177984 - 0.283023058106896 i & 0.185815978760175 - 0.0910521940394691 i\\0.185815978760175 + 0.0910521940394691 i & 0.149322296177984 + 0.283023058106896 i & 1.0 & 0.149322296177984 - 0.283023058106896 i\\0.0940509850777074 + 0.094154537369141 i & 0.185815978760175 + 0.0910521940394691 i & 0.149322296177984 + 0.283023058106896 i & 1.0\end{matrix}\right]

そして H~\tilde{H}

import itertools

# Assemble S, the overlap matrix of dimension D:
H_first_row = np.zeros(krylov_dim, dtype=complex)
H_first_row[0] = H_expval

for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Add in ancilla-only measurements:
    for i in range(krylov_dim - 1):
        # Get expectation values from experiment
        expval_real = results.data.evs[2 + 2 * obs_idx][0][
            i
        ]  # automatic extrapolated evs if ZNE is used
        expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][
            i
        ]  # automatic extrapolated evs if ZNE is used

        # Get expectation values
        expval = expval_real + 1j * expval_imag
        H_first_row[i + 1] += prefactors[i] * coeff * expval

H_first_row_list = H_first_row.tolist()

H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in itertools.product(range(krylov_dim), repeat=2):
    if i >= j:
        H_eff_circ[j, i] = H_first_row[i - j]
    else:
        H_eff_circ[j, i] = np.conj(H_first_row[j - i])
from sympy import Matrix

Matrix(H_eff_circ)

Output:

[10.0−3.02044405310714−2.80721615865252i0.496487054782717+0.188101957039621i1.0770511571923+0.104340737159455i−3.02044405310714+2.80721615865252i10.0−3.02044405310714−2.80721615865252i0.496487054782717+0.188101957039621i0.496487054782717−0.188101957039621i−3.02044405310714+2.80721615865252i10.0−3.02044405310714−2.80721615865252i1.0770511571923−0.104340737159455i0.496487054782717−0.188101957039621i−3.02044405310714+2.80721615865252i10.0]\displaystyle \left[\begin{matrix}10.0 & -3.02044405310714 - 2.80721615865252 i & 0.496487054782717 + 0.188101957039621 i & 1.0770511571923 + 0.104340737159455 i\\-3.02044405310714 + 2.80721615865252 i & 10.0 & -3.02044405310714 - 2.80721615865252 i & 0.496487054782717 + 0.188101957039621 i\\0.496487054782717 - 0.188101957039621 i & -3.02044405310714 + 2.80721615865252 i & 10.0 & -3.02044405310714 - 2.80721615865252 i\\1.0770511571923 - 0.104340737159455 i & 0.496487054782717 - 0.188101957039621 i & -3.02044405310714 + 2.80721615865252 i & 10.0\end{matrix}\right]

最後に、 H~\tilde{H} の一般化固有値問題を解くことができる：

H~c⃗=cSc⃗\tilde{H} \vec{c} = c S \vec{c}

を計算し、基底状態のエネルギーの推定値を得る。 cminc_{min}

gnd_en_circ_est_list = []
for d in range(1, krylov_dim + 1):
    # Solve generalized eigenvalue problem
    gnd_en_circ_est = solve_regularized_gen_eig(
        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=1e-1
    )
    gnd_en_circ_est_list.append(gnd_en_circ_est)
    print("The estimated ground state energy is: ", gnd_en_circ_est)

Output:

The estimated ground state energy is:  10.0
The estimated ground state energy is:  5.933953916292923
The estimated ground state energy is:  4.4101773995740645
The estimated ground state energy is:  3.921288588521255

一粒子セクターの場合、ハミルトニアンのこのセクターの基底状態を古典的に効率よく計算することができる

gs_en = single_particle_gs(H_op, n_qubits)

Output:

n_sys_qubits 10
n_exc 1 , subspace dimension 11
single particle ground state energy:  2.391547869638771

len(H_op)

Output:

27

plt.plot(
    range(1, krylov_dim + 1),
    gnd_en_circ_est_list,
    color="blue",
    linestyle="-.",
    label="KQD estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.title("Estimating Ground state energy with Krylov Quantum Diagonalization")
plt.show()

Output:

Output of the previous code cell

5. 考察と発展

要約すると、まず参照状態から始め、それをさまざまな期間進化させてユニタリー・クリロフ部分空間を生成する。 ハミルトニアンをその部分空間に投影する。 また、部分空間ベクトルの重なりも推定する。 最後に、低次元の一般化固有値問題を古典的に解く。

QKDのフローチャートの概要：参照状態から開始し、状態を近似クリロフベクトルに進化させ、クリロフ部分空間に投影し、投影された部分空間を古典的に対角化し、基底状態の特性を決定する。

古典的なクリロフ技術と量子力学的なクリロフ技術の計算コストを決定するものを比較してみよう。 すべてのステップにおいて、古典的アプローチと量子的アプローチの間に完全な類似性があるわけではない。 この表は、さまざまなステップのスケーリングを集めたものである。

クリロフ法に対する古典的および量子的アプローチにおける様々なプロセスのスケーリングを説明する表。 量子のステップの中には、アナログのないものもある。 スケーリングは本文に記載されているものと同じ。

ハミルトニアンには、一般的に（互いに可換でないため）同時に測定できない項が含まれていることを思い出してください。 ハミルトニアンの項を、すべて同時に測定可能な、互いに可換なパウリ演算子のグループに分類します。また、互いに可換でないすべての項を説明するためには、このようなグループを多数必要とする場合があります。 量子コンピュータ上で H~\tilde{H} を構築するには、ハミルトニアン内の可換なパウリ文字列の各グループごとに個別の測定を行う必要があり、そのそれぞれに多数のショットが必要となる。 r2r^2 個の異なる行列要素について、これを行う必要があり、これらは r2r^2 通りの異なる時間発展因子の組み合わせに対応する。 これを削減する方法がいくつかある場合もありますが、この大まかな処理では、これに必要な時間は Nshots×NGCP×r2N_\text{shots}\times N_\text{GCP} \times r^2 のオーダーで増加します。 SS の各要素は推定する必要がありますが、その計算量は O(Nshots×r2)O(N_\text{shots}\times r^2) のオーダーで増加します。最後に、投影空間における一般化固有値問題を解くには、従来の手法では O(r3)O(r^3) の計算量が必要です。

量子クリロフ対角化は、ハミルトニアン中のパウリ群の数が比較的少ない場合に有効であることがわかる。 これらのスケーリング依存性は、クリロフ法が有用なアプリケーションとそうでないアプリケーションを示唆している。 ハミルトニアンの中には、量子ビットに写像したときに複雑さを増すものがあり、多数の非共通パウリ文字列を含み、少数の共通群に簡単に分割できないものがある。 例えば、量子化学の問題ではよくあることだ。 この複雑さは、近い将来の量子コンピューターに2つの主要な課題をもたらす：

  • H~\tilde{H} の各要素の推定は、項数が多いため計算コストが高くなる。
  • 必要なトロッター回路は法外に深くなる。

上記の2点は、量子コンピューターがフォールトトレランス（耐障害性）に到達すれば問題は少なくなるが、当面は考慮しなければならない。 量子化学の系よりも "単純な "写像を持つ系であっても、ハミルトニアンに非共約項が多すぎる場合には、同じ障害が発生する可能性がある。 クリロフ法は、ハミルトニアンが比較的少数の通約パウリ群に分割でき、 HH、トロッター回路の実装が容易な場合に最も有用である。 例えば、物理学で注目されている多くの格子モデルでは、これらの条件が両方とも満たされている。 KQDは、基底状態についてほとんど知られていない場合に特に有効である。 これは、その固有の収束保証と、地上状態に関する知識が不十分であるために他の方法が使えないようなシナリオで適用できることに起因する。

KQDは強力なツールであるが、プロトコルの時間のかかる側面、特に投影ハミルトニアンの各要素の推定とクリロフ状態のオーバーラップには改善の余地がある。 別のアプローチとしては、クリロフ法をサンプリングベースの手法と組み合わせて活用する方法がある。

6. 付録

付録I：実時間発展からのクリロフ部分空間

ユニタリー・クリロフ空間は次のように定義される

KU(H,∣ψ⟩)=span{∣ψ⟩,e−iH dt∣ψ⟩,…,e−irH dt∣ψ⟩}\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ |\psi\rangle, e^{-iH\,dt} |\psi\rangle, \dots, e^{-irH\,dt} |\psi\rangle \right\}

後で決めるタイムステップ dtdt。 仮に rr が偶数であると仮定する。 d=r/2d=r/2 を定義する。ハミルトニアンを上記のクリロフ空間に射影すると，クリロフ空間と区別がつかないことに注意してください．

KU(H,∣ψ⟩)=span{ei d H dt∣ψ⟩,ei(d−1)H dt∣ψ⟩,…,e−i(d−1)H dt∣ψ⟩,e−i d H dt∣ψ⟩},\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ e^{i\,d\,H\,dt}|\psi\rangle, e^{i(d-1)H\,dt} |\psi\rangle, \dots, e^{-i(d-1)H\,dt} |\psi\rangle, e^{-i\,d\,H\,dt} |\psi\rangle \right\},

つまり、すべての時間進化が dd タイムステップだけ後ろにシフトしている。 見分けがつかない理由は、行列の要素にある

H~j,k=⟨ψ∣ei j H dtHe−i k H dt∣ψ⟩=⟨ψ∣Hei(j−k)H dt∣ψ⟩\tilde{H}_{j,k} = \langle\psi|e^{i\,j\,H\,dt}He^{-i\,k\,H\,dt}|\psi\rangle=\langle\psi|He^{i(j-k)H\,dt}|\psi\rangle

は、進化時間の全体的なシフトに対して不変であり、時間進化はハミルトニアンと通約するからである。 奇数 rr については、 r−1r-1 の分析を使うことができる。

このクリロフ空間のどこかに、低エネルギー状態が必ず存在することを示したい。 これは [3] の定理 3.1 ：

請求項1： ハミルトニアンのスペクトル範囲（つまり、基底状態エネルギーと最大エネルギーの間）のエネルギー EE に対して、次のような関数 ff が存在する。

  1. f(E0)=1f(E_0)=1
  2. ∣f(E)∣≤2(1+δ)−d|f(E)|\le2\left(1 + \delta\right)^{-d} E0E_0 から 離れたところにある のすべての値に対して、つまり指数関数的に抑制される。 ≥δ\ge\delta EE
  3. f(E)f(E) に対する eijE dte^{ijE\,dt} の線形結合である。 j=−d,−d+1,...,d−1,dj=-d,-d+1,...,d-1,d

以下に証明を示すが、完全で厳密な議論を理解したい人以外は読み飛ばしても構わない。 今のところは、上記の主張の意味に焦点を当てる。 上記の性質3により、上記のシフト・クリロフ空間は状態 f(H)∣ψ⟩f(H)|\psi\rangle を含むことがわかる。これが我々の低エネルギー状態である。 その理由を知るために、 ∣ψ⟩|\psi\rangle をエネルギー固有ベーシスで書いてみよう：

∣ψ⟩=∑k=0Nγk∣Ek⟩,|\psi\rangle = \sum_{k=0}^{N}\gamma_k|E_k\rangle,

ここで ∣Ek⟩|E_k\rangle はk番目のエネルギー固有状態であり、 γk\gamma_k は初期状態 ∣ψ⟩|\psi\rangle におけるその振幅である。これによって表現すると、 f(H)∣ψ⟩f(H)|\psi\rangle は次式で与えられる

f(H)∣ψ⟩=∑k=0Nγkf(Ek)∣Ek⟩,f(H)|\psi\rangle = \sum_{k=0}^{N}\gamma_kf(E_k)|E_k\rangle,

が固有状態 ∣Ek⟩|E_k\rangle に作用するとき、 HH を EkE_k で置き換えることができるという事実を利用している。したがって、この状態のエネルギー誤差は

energy error=⟨ψ∣f(H)(H−E0)f(H)∣ψ⟩⟨ψ∣f(H)2∣ψ⟩\text{energy error} = \frac{\langle\psi|f(H)(H-E_0)f(H)|\psi\rangle}{\langle\psi|f(H)^2|\psi\rangle} =∑k=0N∣γk∣2f(Ek)2(Ek−E0)∑k=0N∣γk∣2f(Ek)2.= \frac{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.

これを理解しやすい上界にするために、まず分子の和を Ek−E0≤δE_k-E_0\le\delta の項と Ek−E0>δE_k-E_0>\delta の項に分ける：

energy error=∑Ek≤E0+δ∣γk∣2f(Ek)2(Ek−E0)∑k=0N∣γk∣2f(Ek)2+∑Ek>E0+δ∣γk∣2f(Ek)2(Ek−E0)∑k=0N∣γk∣2f(Ek)2.\text{energy error} = \frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} + \frac{\sum_{E_k> E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.

第一項を δ\delta で上界することができる、

∑Ek≤E0+δ∣γk∣2f(Ek)2(Ek−E0)∑k=0N∣γk∣2f(Ek)2<δ∑Ek≤E0+δ∣γk∣2f(Ek)2∑k=0N∣γk∣2f(Ek)2≤δ,\frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} < \frac{\delta\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} \le \delta,

ここで、最初のステップは、和のすべての EkE_k について Ek−E0≤δE_k-E_0\le\delta、2番目のステップは、分子の和が分母の和の部分集合であるため、続く。 第2項については、まず分母を ∣γ0∣2|\gamma_0|^2。 f(E0)2=1f(E_0)^2=1 ：すべてを足し合わせると、次のようになる

energy error≤δ+1∣γ0∣2∑Ek>E0+δ∣γk∣2f(Ek)2(Ek−E0).\text{energy error} \le \delta + \frac{1}{|\gamma_0|^2}\sum_{E_k>E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0).

残されたものを単純化するために、これらのすべての EkE_k、 ff の定義によって、 f(Ek)2≤4(1+δ)−2df(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}。さらに、 Ek−E0<2∥H∥E_k-E_0<2\|H\| を上界し、 ∑Ek>E0+δ∣γk∣2<1\sum_{E_k>E_0+\delta}|\gamma_k|^2<1 を上界すると、次のようになる

energy error≤δ+8∣γ0∣2∥H∥(1+δ)−2d.\text{energy error} \le \delta + \frac{8}{|\gamma_0|^2}\|H\|\left(1 + \delta\right)^{-2d}.

これは、どのような δ>0\delta>0 についても成り立つので、 δ\delta を目標誤差に等しく設定すれば、上の誤差境界は、クリロフ次元 2d=r2d=r に従って指数関数的に収束する。また、 δ<E1−E0\delta<E_1-E_0 の場合、 δ\delta の項は、実際には、上記の境界から完全になくなることに注意してください。

議論を完結させるために、上記はクリロフ空間における最低エネルギー状態のエネルギー誤差ではなく、特定の状態 f(H)∣ψ⟩f(H)|\psi\rangle のエネルギー誤差に過ぎないことにまず注意する。 しかし、(Rayleigh-Ritz) 変分原理により、クリロフ空間内の最低エネルギー状態のエネルギー誤差は、クリロフ空間内の任意の状態のエネルギー誤差によって上界されるため、上記は最低エネルギー状態のエネルギー誤差、つまりクリロフ量子対角化アルゴリズムの出力に対する上界でもある。

上記と同様の分析が、ノイズとノートで説明した閾値処理手順を追加して実施できる。 この分析については [2] と [4] を参照のこと。

付録II：請求項１の証明

以下はほとんど [3] の定理 3.1: 0<a<b0 < a < b とし、 Πd∗\Pi^*_d を最大次数 dd の残差多項式（0での値が1である多項式）の空間とする。に対する解は

β(a,b,d)=min⁡p∈Πd∗max⁡x∈[a,b]∣p(x)∣\beta(a, b, d) = \min_{p \in \Pi^*_d} \max_{x \in [a, b]} |p(x)| \quad

次と同一である

p∗(x)=Td(b+a−2xb−a)Td(b+ab−a),p^*(x) = \frac{T_d\left(\frac{b + a - 2x}{b - a}\right)}{T_d\left(\frac{b + a}{b - a}\right)}, \quad

となり、対応する最小値は

β(a,b,d)=Td−1(b+ab−a).\beta(a, b, d) = T_d^{-1}\left(\frac{b + a}{b - a}\right).

我々はこれを複素指数で自然に表現できる関数に変換したい。なぜなら、それが量子クリロフ空間を生成する実時間進化だからである。 そのためには、ハミルトニアンのスペクトル範囲内のエネルギーを、 [0,1][0,1] ：defineの範囲内の数に変換する以下の変換を導入するのが便利である

g(E)=1−cos⁡((E−E0)dt)2,g(E) = \frac{1-\cos\big((E-E_0)dt\big)}{2},

ここで dtdt は −π<E0dt<Emaxdt<π-\pi < E_0dt < E_\text{max}dt < \pi のようなタイムス テップである。 EE が E0E_0 から遠ざかるにつれて、 g(E0)=0g(E_0)=0 と g(E)g(E) が成長することに注意。

ここで、パラメータa, b, dを a=g(E0+δ)a = g(E_0 + \delta), b=1b = 1, d = int( r/2 ) に設定した多項式 p∗(x)p^*(x) を使い、関数を定義する：

f(E)=p∗(g(E))=Td(1+2cos⁡((E−E0)dt)−cos⁡(δ dt)1+cos⁡(δ dt))Td(1+21−cos⁡(δ dt)1+cos⁡(δ dt))f(E) = p^* \left( g(E) \right) = \frac{T_d\left(1 + 2\frac{\cos\big((E-E_0)dt\big) - \cos\big(\delta\,dt\big)}{1 +\cos\big(\delta\,dt\big)}\right)}{T_d\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)}

ここで、 E0E_0 は基底状態のエネルギーである。 f(E)f(E) が次数 dd の三角多項式、つまり j=−d,−d+1,...,d−1,dj=-d,-d+1,...,d-1,d に対する eijE dte^{ijE\,dt} の線形結合であることは、 cos⁡(x)=eix+e−ix2\cos(x)=\frac{e^{ix}+e^{-ix}}{2} を挿入することでわかる。さらに、上記の p∗(x)p^*(x) の定義から、 f(E0)=p(0)=1f(E_0)=p(0)=1 と、 ∣E−E0∣>δ\vert E-E_0 \vert > \delta のようなスペクトル範囲内の任意の EE に対して、次のようになる

∣f(E)∣≤β(a,b,d)=Td−1(1+21−cos⁡(δ dt)1+cos⁡(δ dt))|f(E)| \le \beta(a, b, d) = T_d^{-1}\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right) ≤2(1+δ)−d=2(1+δ)−⌊k/2⌋.\leq 2\left(1 + \delta\right)^{-d} = 2\left(1 + \delta\right)^{-\lfloor k/2\rfloor}.

参照資料：

[1 ]https://arxiv.org/abs/2407.14431

[2 ]https://arxiv.org/abs/1811.09025

[3 ]https://people.math.ethz.ch/~mhg/pub/biksm.pdf

[4 ]https://academic.oup.com/book/36426

[5] https://en.wikipedia.org/wiki/Krylov

[6] クリロフ部分空間法：原理と分析, イェルク・リーゼン, ズデネク・ストラコス https://academic.oup.com/book/36426

[7] スパース線形システムの反復法」 ユセフ・サード著

[8] "MINRES-QLP: A Krylov Subspace Method for Indefinite or Singular Symmetric Systems" by Sou-Cheng Choi, Christopher Paige, and Michael Saunders. マイケル・サンダース( https://epubs.siam.org/doi/10.1137/100787921 )

[9] イーサン・エッパリー、リン・リン、中務裕司。 「量子部分空間対角化の理論」。 SIAM Journal on Matrix Analysis and Applications 43, 1263-1290 (2022).

[10 ]https://link.aps.org/doi/10.1103/PRXQuantum.4.030319

このページは役に立ちましたか？
バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。