はじめに
量子コンピューターアドベントカレンダー19日目の記事になります。
今回はタイトルの通り、VQLSというアルゴリズムについて調査したので、それを少し詳しめに解説します。去年の「時を巻き戻す量子アルゴリズム」同様、割と本気で書いた結果、かなり長くなっていますのでゆっくり読み進めてください。難易度は中級レベルです。なお、コード実装にはQiskitを利用します。
原論文はこちらです。
「Variational Quantum Linear Solver」（arXiv:1909.05820 [quant-ph]）
$$
\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}}
$$
量子コンピューターで連立一次方程式を解く方法
実は、量子コンピューターで「高速に解くことができる」連立一次方程式を解く方法として、Harrow-Hassidim-Lloyd (HHL) アルゴリズムという手法が2009年に提唱されています。アルゴリズム名は提唱者の３名の頭文字が使われています。
このアルゴリズム自体は、量子位相推定アルゴリズムを使って解く手法となっています。ただし、HHLアルゴリズムを効果的に使うためにはいくつか前提があります。
・位相推定を使うため、十分なエラー耐性があるデバイス上で動くことを前提としている
・「高速に解く」ためには、問題となる連立一次方程式の係数の多くが0であること（スパース）であることが望ましい
つまり、まだこのアルゴリズムを上手く使うためのハードルが高く、効果を出すには条件が厳しいという話になります。
HHLとは別のアプローチで実現するアルゴリズムが2019年に提唱されました。それが今回の主役のVariational Quantum Linear Solver（VQLS）になります。
Variational Quantum Linear Solver（VQLS）とは
Variationalという単語がついているのでわかる人にはわかると思いますが、変分法を使うアルゴリズムになります。
変分法を用いることで、エラー耐性の低いNISQでも動くことを前提としたアルゴリズムであることが、HHLアルゴリズムとの大きな違いと言えます。
変分法とは（ご存知の方は読み飛ばして大丈夫です）
すごく簡単に説明すると、問題設定に合わせてパラメーターとコスト関数を用意し、コスト関数を最小化するようにパラメーターの微調整を繰り返す手法です。機械学習などではおなじみの方法です。
量子計算になった場合でも、量子回路の一部のゲートをパラメーターで調整できるようにして、観測と調整を繰り返して、イイ感じのパラメーターを見つけ出します。観測と調整は古典コンピューターを使うのである意味ではハイブリッド的に動かすことになります。
VQLSの実装方針
問題設定を以下のように考えます。
次のような$n$次連立一次方程式を考えます。各$a_{ij}$および$b_{k}$は全て既知であり何らかの値があるとして、$n$個ある$x_{l}$を求めたいとします。
\left\{
\begin{array}{ll}
a_{11}x_{1} + \cdots + a_{1n}x_{n} = b_{1} \\
\vdots \\
a_{n1}x_{1} + \cdots + a_{nn}x_{n} = b_{n}
\end{array}
\right.
この連立方程式は、線形代数の言葉で行列とベクトルを使って、次のように書き直すことができます。
\begin{pmatrix}
a_{11} & \cdots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \cdots & a_{nn}
\end{pmatrix}
\begin{pmatrix}
x_{1} \\
\vdots \\
x_{n}
\end{pmatrix}
=
\begin{pmatrix}
b_{1} \\
\vdots \\
b_{n}
\end{pmatrix}
ここで、以下の式になるように$A$、$b$、$x$を行列とベクトルとして置き換えます。
A =
\begin{pmatrix}
a_{11} & \cdots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \cdots & a_{nn}
\end{pmatrix}
、x =
\begin{pmatrix}
x_{1} \\
\vdots \\
x_{n}
\end{pmatrix}
、
b=
\begin{pmatrix}
b_{1} \\
\vdots \\
b_{n}
\end{pmatrix}
上記のように見なすと、以下のようにスッキリした表現にできます。
$$ Ax = b $$
まずは問題の設定はここまでです。この時点では量子かどうかは考えていません。
次に処理の大まかな流れは以下の図（原論文より）で表現されています。ステップを１つずつ紹介していきます。
ステップ①
以下の３つのことをやります。
1.行列$A$をユニタリ行列の線形和に分解する
2.ベクトル$b$を規格化し、量子状態$\ket{\psi}$とみなして、$\ket{0}$からその$\ket{b}$を生成するユニタリ行列$U$を用意する
3.コスト関数を定義しておく
まず、問題の前提として2つ気にしておくことがあります。
1つ目は問題設定に出てくる行列$A$は任意の行列となりますので、ユニタリ行列でない場合があります。むしろユニタリ行列ではない方が一般的と言えます。
それに関連して2つ目には、量子計算で行列を扱う場合、扱えるのはユニタリ行列のみになります。
またベクトル$b$と$x$も一般の行列であり量子計算上は、規格化されているベクトルでなければなりません。これは量子論の要請です。なので、それに合わせた事前処理が必要になりますので、次のようなことをしておく必要があります。
ステップ①-1 行列Aのユニタリ行列の線形和に分解
行列$𝐴$をユニタリ行列の線形和に分解します。つまり、何らかの$n$次元ユニタリ行列$A_{k}$と複素数$c_{k}$を用いて以下のように表現をし直します。
A = \sum_{k}^{} c_{k}A_{k}
ステップ①-2 ベクトルの規格化
$ \ket{𝑏}=𝑈 \ket{0}$となる$U$を定義しておく必要があります。定義というと畏まり過ぎているので、どちらかと言えば$U$を探すという方が適切かもしれません。
後述実験では、簡単なものを用意して使っているので探すようなことはしていません。
ステップ①-3 コスト関数の定義
ここから少し長い議論と計算を行うので、急いで読みたい方は「コスト関数計算まとめ」まで飛んでください。
コスト関数のベースとなる考え方は、「正解から遠いほどコストが大きく、正解に近づくとコストが小さくなる」とすることにあります。この考え方に基づいて今回はコスト関数を比較的シンプルに構成することができます。
まず、求めたいベクトル$x$の規格化されたベクトルを$\ket{\psi} = \frac{\ket{x}}{|x|}$としましょう。この$\ket{\psi}$となるような何からのユニタリ行列$V$があるとします。つまり、
$$ \ket{\psi} = V\ket{0} $$
となるとしましょう。この$V$を回路で表現できると、解を求められたことに相当します。実際にはこの$V$を再現できるパラメーターを探すことになります。
では、いったん何らかのパラメーター$k$を用いて$V(k)$としているとして、$\ket{\psi(k)} = V(k)\ket{0}$としておきます。$\ket{\psi}$は$\ket{x}$を再現しようとしているベクトルとなります。ここに$A$を左からかけると、それは$\ket{b}$を再現しようとするベクトルになっているはずです。つまり、$AV(k)\ket{0}$がどれだけ規格化する前のベクトル$\ket{b}$に似ているかどうかが、コスト関数を考える上での焦点になります。
そのために、$ H_{p} = I - \ket{\psi}\bra{\psi} $という射影演算子を考えます。これを任意の状態$\ket{\psi(k)}$での期待値を計算してみます。
$$ C_{p} = \bra{\psi(k)}H_{p}\ket{\psi(k)} = \braket{\psi(k)}{\psi(k)} - \braket{\psi(k)}{b}\braket{b}{\psi(k)} $$
これを見ると、$\ket{\psi(k)}$が$\ket{b}$に近いほど第二項目が大きくなり全体の値が小さくなることがわかります。つまりこれをコスト関数として採用できそうです。
ただし、1点注意しないといけないのは、行列$A$はユニタリ行列とは限らないので、それを作用させたベクトルの長さが変化する可能性が出てきます。（なので、$\braket{\psi(k)}{\psi(k)}$を1として計算しなかったという意図があります。）この可能性が意味するところは、ベクトルの長さが変化すると、たとえ同じ向きのベクトルであってもコストが変化します。これではあくまでベクトルの向きに着目してコスト関数を定義したいという話になりますので、もうひと手間加えます。
シンプルに、$\braket{\psi(k)}{\psi(k)}$で全体を割ってみることにします。（ただし、$\ket{\psi(k)}\neq0$であることが必要です。冷静に考えればこれは妥当な前提です。もし零ベクトルであったとすると、そもそも$\ket{b}$も零ベクトルになるので、それは自明な解になるわけですから、わざわざ苦労して量子回路を組む意味が無くなるケースなので・・・）
こうすると、
$$ \hat{C_{p}} = 1- \frac{\braket{\psi(k)}{b}\braket{b}{\psi(k)}}{\braket{\psi(k)}{\psi(k)}} = 1- \frac{|\braket{\psi(k)}{b}|^2}{\braket{\psi(k)}{\psi(k)}} $$
になるので、長さが変わってしまうことを心配する必要は無くなります。あくまでも$\ket{b}$と$\ket{\psi(k)}$の向きだけでコスト関数（正解に近いベクトルであるかどうか）を定義することができました。
さて、コスト関数は定義できたのですが、1つ課題が出てきました。このコスト関数を具体的に計算するにはどうした良いのでしょうか？式の構造はわかりましたが、このままではよくわかりませんし、使えません。まずはいったん落ち着いて、分母と分子に分けてそれぞれがどんな風になっているかを見ていきます。
コスト関数の分母をどうやって求めるか
まずは分母$\braket{\psi(k)}{\psi(k)}$から。これは、①-1で準備でした行列$A$と$V(k)$で書き直します。
$$ \braket{\psi(k)}{\psi(k)} = \bra{\psi(k)} A^{\dagger}A\ket{\psi(k)} = \bra{0}V(k)^{\dagger}A^{\dagger} A V(k) \ket{0}
= \sum_{i}^{}\sum_{j}^{} c_{j}^{\ast}c_{i}
\bra{0} V(k)^{\dagger}A_{j}^{\dagger} A_{i} V(k) \ket{0} = \sum_{i}^{}\sum_{j}^{} \bra{\psi(k)} c_{j}^{\ast}c_{i} A_{j}^{\dagger} A_{i} \ket{\psi(k)} $$
と書き下せます。「量子状態$\ket{0}$に対するユニタリ行列$V(k)^{\dagger} A_{j}^{\dagger} A_{i}V(k)$の期待値」または「量子状態$\ket{\psi(k)}$に対するユニタリ行列$A_{j}^{\dagger}A_{i}$の期待値」とみなせて、それらに複素係数をかけた和から求めることができようになります。ここまでくると「XXの期待値」の部分さえ計算できれば、古典コンピューターで算出可能です。そうなると「XXの期待値」はどうやって求めるんだ？という疑問が湧いてきます。
これの答えとしては、Hadamrdテストと呼ばれる回路から計算できます。
Hadamardテストとは
ここでは少し脱線してHadamardテストとは何かを見ていきます。この後の議論でも出てくるのでここで一度整理しておきます。
結論から申し上げると、「状態$\ket{\psi}$に関するユニタリ行列$U$の期待値」を求めることができます。まさに、今回のにピッタリですね！
実際には以下のような回路で構成されています。
ここではHadamardテストを見ていきましょう！順にゲートを通していきます。
まずは量子ビット0にHadamardゲートを作用させます。
\ket{0} \otimes \ket{\psi} \xrightarrow{H\otimes I} H\ket{0} \otimes \ket{\psi} = \frac{1}{\sqrt{2}}(\ket{0}+\ket{1}) \otimes \ket{\psi} = \frac{1}{\sqrt{2}} (\ket{0}\otimes\ket{\psi} + \ket{1}\otimes\ket{\psi})
次にcontrolled-Uゲートを作用させます。
\xrightarrow{controlled-U} \frac{1}{\sqrt{2}} (\ket{0}\otimes\ket{\psi} + \ket{1}\otimes U\ket{\psi})
もう一度、量子ビット0にHadamardゲートを作用させて、量子ビット0の状態が$\ket{0}$か$\ket{1}$の塊で式を整理します。
\xrightarrow{H\otimes I} \frac{1}{\sqrt{2}} (H\ket{0}\otimes\ket{\psi} + H\ket{1}\otimes U\ket{\psi}) = \frac{1}{\sqrt{2}} ( \frac{1}{\sqrt{2}}(\ket{0}+\ket{1}) \otimes \ket{\psi} + \frac{1}{\sqrt{2}}(\ket{0}-\ket{1}) \otimes U \ket{\psi} ) = \ket{0} \otimes \frac{I+U}{2} \ket{\psi} + \ket{1} \otimes \frac{I-U}{2} \ket{\psi}
割とスッキリする所まできました。
ここで、量子ビット0が状態$\ket{0}$が測定される確率と、状態$\ket{1}$が測定される確率を求めると、「状態$\ket{\psi}$に関するユニタリ行列$U$の期待値」が現れくるのです！
確率はどうやって求めるかというと、量子ビットのそれぞれの状態に付随している量子ビット1の方の状態の内積を計算すれば良いのです。なので、そういった式変形をしていました。
それぞれの状態の確率を$P(x)$と定義しましょう。
P(0) = \bra{\psi}\frac{(U^{\dagger}+I)(U+I)}{4}\ket{\psi} = \bra{\psi}\frac{(2I
+ U^{\dagger} + U)}{4}\ket{\psi} =\frac{1}{4}(2\braket{\psi}{\psi} + \bra{\psi} U \ket{\psi} + \bra{\psi} U^{\dagger} \ket{\psi}) = \frac{1}{2}(1+Re(\bra{\psi}U\ket{\psi}))
同様にP(1)を計算すると、
P(1) = \frac{1}{2}(1-Re(\bra{\psi}U\ket{\psi}))
これらを$Re(\bra{\psi}U\ket{\psi})$について整理すると、
Re(\bra{\psi}U\ket{\psi}) = P(0)-P(1)
となります。この式の意味するところは、状態$\ket{0}$が観測される確率と状態$\ket{1}$が観測される確率を見て、それらを使って計算できることを示しています。
イメージを図にすると、以下の矢印の部分を求めたかった期待値と捉えることができます。
「なるほど、これで期待値の実部は得られそうだけど、虚部はどうするの？」という疑問があるかもしれません。これの回答としては、今回の問題設定では実数のみを扱うことにするため、ユニタリ行列も状態ベクルも実数のみで構成できるので、無視して良くできます。なにで今回は求めないことにします。一般のユニタリ行列と状態ベクトルで考える場合は、虚部を求める場合もあるでしょう。その場合は、もう一手間加える必要があります。
その答えはこちらの記事に書いてあるので、是非読んでてみください。
https://qiita.com/SamN/items/917e74bcb794b4c6e776
話を戻すと、Hadamardテストを使って、求めたい期待値を得ることができそうです。Hadamradテストの回路を参考にして、今回はこんな回路を組むと「量子状態$\ket{\psi(k)}$に対するユニタリ行列$A_{j}^{\dagger}A_{i}$の期待値」が求められそうです。これで分母の計算はこれらを１つずつ求めて係数かけて足し上げるだけで行けそうです。
コスト関数の分子を求めてみよう
一方で、分子$|\braket{b}{\psi}|^2$の計算方法を考えます。方針はこれまでと全く同じく$A_{i}$で展開していきます。
|\braket{b}{\psi}|^2 = |\bra{b}AV(k)\ket{0}|^2 = |\bra{0}U^{\dagger}AV(K)\ket{0}|^2 = \bra{0}V^{\dagger}(k)A^{\dagger}U\ket{0}\bra{0}U^{\dagger}AV(K)\ket{0} = \sum_{i}^{}\sum_{j}^{} c_{i}c^{*}_{j} \bra{0}V^{\dagger}(k)A^{\dagger}_{j}U\ket{0}\bra{0}U^{\dagger}A_{i}V(k)\ket{0}
ここで、$\bra{0}U^{\dagger}A_{i}V(k)\ket{0}$を求める方法を考えますが、同じようにHadamardテストを構成することで求めることができます。
これらを全ての$A_{i}$について、１つずつ求めて、かけては足しを繰り返せば良いのです
コスト関数計算方法まとめ
長い長い計算の果てに、コスト関数の計算方法を議論してきました。ここでまとめておくと以下のような関係性になっています。一見すると複雑な式になっていますが、１つずつHadamardテストでパーツを求めて足し上げていくことになります。
後でこれを実装していきます。
さて、これでやっとコスト関数を求める方法も見つかりました。
これで準備段階は終わりました。
ステップ②
次のステップでは、量子状態$\ket{x}=V\ket{0}$を生成する回路を構成します。ここでは、3量子ビットのケースに限定して話をします。
ここでは、以下のような回路を考えてみることにします。
def apply_fixed_ansatz(qubits, parameters):
for iz in range (0, len(qubits)):
circ.ry(parameters[0][iz], qubits[iz])
circ.cz(qubits[0], qubits[1])
circ.cz(qubits[2], qubits[0])
for iz in range (0, len(qubits)):
circ.ry(parameters[1][iz], qubits[iz])
circ.cz(qubits[1], qubits[2])
circ.cz(qubits[2], qubits[0])
for iz in range (0, len(qubits)):
circ.ry(parameters[2][iz], qubits[iz])
circ = QuantumCircuit(3)
apply_fixed_ansatz([0, 1, 2], [[1, 1, 1], [1, 1, 1], [1, 1, 1]]) # 全て変分パラメーターが1としてセットした例
これを表示すると、こんな回路になります。
ポイントはこれのパラメーターは9つで構成されていて、実行の際にはこの回路の構成を一切変えずにRyゲートのパラメーターのみを調整するようにしていきます。そういう用途もあってか、Fixed Hardware Ansatzという名称と呼んでいます。（一般的な言葉ではないと思いますが）
ステップ③
コスト関数が閾値γよりも大きい場合は、𝛼の値を変分法により再計算し、新たに量子回路に入れて計算を行い、これを閾値よりも小さくなるまで繰り返します。
ステップ④
コスト関数の値が閾値以下になった𝛼を$α_{opt}$として、再度状態$\ket{x}$を作る回路に投入することで、求めたい連立一次方程式の解$x$の規格化された量子状態$\ket{x}$を得ることができます。
VQLSの実装と実験
実験の設定
今回は簡易的に動かすために、以下のような設定を前提とします。
・3量子ビット（最大8変数）とする
・$U=H\otimes H\otimes H$とする ⇒ $\ket{b} = H\otimes H\otimes H \ket{0} = (\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{t}$
・係数は全て実数とする ⇒ $行列Aの要素は全て実数とする$
・行列$A$は、$I$と$Z$で分解できるものに限定する
・初期変分パラメーターαを、$\frac{k}{1000} (k\in {0,3000})$ としてランダムに選ぶ
・変分パラメーターの最適化にはCOBYLA（ constrained optimization by linear approximation method ）法を使う
実装してみよう！
それでは実装を始めていきます！
コスト関数の分子のHadamardテストの実装
では、まずはHadamardテストを実装してみましょう。$V(k)$は上述したFixed hardware Ansatzの構造を使います。
コードはこうなります。関数として定義する部分と、可視化する処理として書いておきます。
def had_test(gate_type, qubits, ancilla_index, parameters):
circ.h(ancilla_index)
apply_fixed_ansatz(qubits, parameters)
for ie in range (0, len(gate_type[0])):
if (gate_type[0][ie] == 1):
circ.cz(ancilla_index, qubits[ie])
for ie in range (0, len(gate_type[1])):
if (gate_type[1][ie] == 1):
circ.cz(ancilla_index, qubits[ie])
circ.h(ancilla_index)
circ = QuantumCircuit(4)
had_test([[0, 0, 0], [0, 0, 1]], [1, 2, 3], 0, [[1, 1, 1], [1, 1, 1], [1, 1, 1]])
circ.draw()
had_test([[0, 0, 0], [0, 0, 1]], [1, 2, 3], 0, [[1, 1, 1], [1, 1, 1], [1, 1, 1]])の引数について説明します。
最初の[[0,0,0],[0,0,1]]はどのゲートをどこに入れるかをフラグ化しています。前がIゲートで、後ろがZゲートです。3つの要素は、どの量子ビットに入れるかです。つまり[0,0,1]は量子ビット3にZゲートを作用させることを意味しています。
２つ目の引数[1,2,3]は実際の連立方程式の式に相当する部分に対してどの量子ビットを使うかという宣言です。量子ビット１～３を連立方程式の式として使われます。
３つの引数0は、どれを補助ビットとするかを宣言しています。
４つの引数[[1,1,1],[1,1,1],[1,1,1]]はFixed Hardware Ansatzに入る変分パラメーターの設定です。ここでは可視化のために変分パラメーターを1にセットしていますが、本来は不定の変分パラメーターが入ることになります。
可視化するとこんな感じです。
コスト関数の分母部分のHadamardテストの実装
同じくコードはこのようになります。
def special_had_test(gate_type, qubits, ancilla_index, parameters, reg):
circ.h(ancilla_index)
control_fixed_ansatz(qubits, parameters, ancilla_index, reg)
for ty in range (0, len(gate_type)):
if (gate_type[ty] == 1):
circ.cz(ancilla_index, qubits[ty])
control_b(ancilla_index, qubits)
circ.h(ancilla_index)
q_reg = QuantumRegister(5)
circ = QuantumCircuit(q_reg)
special_had_test([0, 0, 1], [1, 2, 3], 0, [[1, 1, 1], [1, 1, 1], [1, 1, 1]], q_reg)
circ.draw(output='mpl')
引数についてちょっとだけ補足すると、最初の引数[0,0,1]はZゲートをどの量子ビットに配置するかを指定しています。
可視化するとこんな感じです。
分子部分の計算のときに回路は似ているのですが、$V(k)$の部分が制御$V(k)$になっている点に注意が必要です。分子部分のときはCZゲートが出てきた部分はCCZに、Ryゲートは制御Ryゲートに置き換える必要があります。その点が回路を複雑に見せています。
コスト関数の実装
上記の関数を使って、以下のように１つの関数にまとめます。
def calculate_cost_function(parameters):
global opt
overall_sum_1 = 0
parameters = [parameters[0:3], parameters[3:6], parameters[6:9]]
for i in range(0, len(gate_set)):
for j in range(0, len(gate_set)):
global circ
qctl = QuantumRegister(5)
qc = ClassicalRegister(5)
circ = QuantumCircuit(qctl, qc)
backend = Aer.get_backend('statevector_simulator')
multiply = coefficient_set[i]*coefficient_set[j]
had_test([gate_set[i], gate_set[j]], [1, 2, 3], 0, parameters)
job = execute(circ, backend)
result = job.result()
outputstate = np.real(result.get_statevector(circ, decimals=100))
o = outputstate
m_sum = 0
for l in range (0, len(o)):
if (l%2 == 1):
n = o[l]**2
m_sum+=n
overall_sum_1+=multiply*(1-(2*m_sum))
overall_sum_2 = 0
for i in range(0, len(gate_set)):
for j in range(0, len(gate_set)):
multiply = coefficient_set[i]*coefficient_set[j]
mult = 1
for extra in range(0, 2):
qctl = QuantumRegister(5)
qc = ClassicalRegister(5)
circ = QuantumCircuit(qctl, qc)
backend = Aer.get_backend('statevector_simulator')
if (extra == 0):
special_had_test(gate_set[i], [1, 2, 3], 0, parameters, qctl)
if (extra == 1):
special_had_test(gate_set[j], [1, 2, 3], 0, parameters, qctl)
job = execute(circ, backend)
result = job.result()
outputstate = np.real(result.get_statevector(circ, decimals=100))
o = outputstate
m_sum = 0
for l in range (0, len(o)):
if (l%2 == 1):
n = o[l]**2
m_sum+=n
mult = mult*(1-(2*m_sum))
overall_sum_2+=multiply*mult
print(1-float(overall_sum_2/overall_sum_1))
return 1-float(overall_sum_2/overall_sum_1)
長くて読むのが大変そうに見えますが、これはこれまでの流れにそった処理を１つずつやっているに過ぎません。繰り返しになるので、ここでは説明しませんが、もう一度流れを押さえた上で、コードを読んでみてください。このコスト関数をCOBLAに突っ込んでパラメーターを最適化します。
実問題を解いてみよう！
では２つほど問題を設定してコードを動かしみましょう。
問題1
$A=0.45Z_{3}+0.55I$
こちらの問題を解く処理をコード化したのが以下です。
coefficient_set = [0.55, 0.45]
gate_set = [[0, 0, 0], [0, 0, 1]]
out = minimize(calculate_cost_function, x0=[float(random.randint(0,3000))/1000 for i in range(0, 9)], method="COBYLA", options={'maxiter':200})
print(out)
out_f = [out['x'][0:3], out['x'][3:6], out['x'][6:9]]
circ = QuantumCircuit(3, 3)
apply_fixed_ansatz([0, 1, 2], out_f)
backend = Aer.get_backend('statevector_simulator')
job = execute(circ, backend)
result = job.result()
o = result.get_statevector(circ, decimals=10)
a1 = coefficient_set[1]*np.array([[1,0,0,0,0,0,0,0], [0,1,0,0,0,0,0,0], [0,0,1,0,0,0,0,0], [0,0,0,1,0,0,0,0], [0,0,0,0,-1,0,0,0], [0,0,0,0,0,-1,0,0], [0,0,0,0,0,0,-1,0], [0,0,0,0,0,0,0,-1]])
a2 = coefficient_set[0]*np.array([[1,0,0,0,0,0,0,0], [0,1,0,0,0,0,0,0], [0,0,1,0,0,0,0,0], [0,0,0,1,0,0,0,0], [0,0,0,0,1,0,0,0], [0,0,0,0,0,1,0,0], [0,0,0,0,0,0,1,0], [0,0,0,0,0,0,0,1]])
a3 = np.add(a1, a2)
b = np.array([float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8))])
print((b.dot(a3.dot(o)/(np.linalg.norm(a3.dot(o)))))**2)
では動かしてみると、１回の試行ごとのコスト関数の値は以下のように推移していきます。
0.7260295360394795
0.7752900217571115
0.904521042497052
0.6724723336557688
0.6868949431101335
・
・
・
0.04575109909207409
0.04639360288121852
0.04538357433638929
fun: 0.04538357433638929
maxcv: 0.0
message: 'Maximum number of function evaluations has been exceeded.'
nfev: 200
status: 2
success: False
x: array([1.88514537, 0.00847878, 3.04330519, 1.75327891, 0.08770506,
2.43756353, 1.63178745, 2.09752409, 3.62520129])
(0.9546164256514695-0j)
最初のコストが0.7以上に対して、最終的なコスト関数の値が0.04まで下げることができました。
２つのベクトルの類似度は、一番最後に書かれている0.95XXXの部分ですが、これが1に近いほどベクトルが似ていることが客観的にわかります。かなり近い値になったことがわかります。
問題2
$A=0.55I+0.225Z_{2}+0.225Z_{3}$
同様にコード化してみます。大きな違いはcoefficient_setの部分になります。
coefficient_set = [0.55, 0.225, 0.225]
gate_set = [[0, 0, 0], [0, 1, 0], [0, 0, 1]]
out = minimize(calculate_cost_function, x0=[float(random.randint(0,3000))/1000 for i in range(0, 9)], method="COBYLA", options={'maxiter':200})
print(out)
out_f = [out['x'][0:3], out['x'][3:6], out['x'][6:9]]
circ = QuantumCircuit(3, 3)
apply_fixed_ansatz([0, 1, 2], out_f)
backend = Aer.get_backend('statevector_simulator')
job = execute(circ, backend)
result = job.result()
o = result.get_statevector(circ, decimals=10)
a1 = coefficient_set[2]*np.array([[1,0,0,0,0,0,0,0], [0,1,0,0,0,0,0,0], [0,0,1,0,0,0,0,0], [0,0,0,1,0,0,0,0], [0,0,0,0,-1,0,0,0], [0,0,0,0,0,-1,0,0], [0,0,0,0,0,0,-1,0], [0,0,0,0,0,0,0,-1]])
a0 = coefficient_set[1]*np.array([[1,0,0,0,0,0,0,0], [0,1,0,0,0,0,0,0], [0,0,-1,0,0,0,0,0], [0,0,0,-1,0,0,0,0], [0,0,0,0,1,0,0,0], [0,0,0,0,0,1,0,0], [0,0,0,0,0,0,-1,0], [0,0,0,0,0,0,0,-1]])
a2 = coefficient_set[0]*np.array([[1,0,0,0,0,0,0,0], [0,1,0,0,0,0,0,0], [0,0,1,0,0,0,0,0], [0,0,0,1,0,0,0,0], [0,0,0,0,1,0,0,0], [0,0,0,0,0,1,0,0], [0,0,0,0,0,0,1,0], [0,0,0,0,0,0,0,1]])
a3 = np.add(np.add(a2, a0), a1)
b = np.array([float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8)),float(1/np.sqrt(8))])
print((b.dot(a3.dot(o)/(np.linalg.norm(a3.dot(o)))))**2)
こんな結果になります。
0.9799856290093593
0.9191792111360295
0.973364494052271
0.8305003425889084
0.7516802955839005
・
・
・
0.0043649061738640915
0.0043229904460660995
0.004229819596867079
0.004181852882589876
fun: 0.004181852882589876
maxcv: 0.0
message: 'Maximum number of function evaluations has been exceeded.'
nfev: 200
status: 2
success: False
x: array([2.85165255, 0.21129123, 3.1363429 , 2.89342992, 0.37618081,
2.28115687, 1.72738044, 2.64220854, 3.60407879])
(0.9958181471157969-0j)
こちらも最終コストが0.004まで落ちています。類似度も0.99とほぼ解と同じくらいのベクトルが再現できています。
やってみた疑問点
実際に利用することを想定した場合、行列Aをユニタリ行列での線形結合を古典コンピューターないしは人手で探してから動かさないといけないのでは？
何らかの古典処理を前提とするため、その分の時間がかかることを想定しておく必要がありそうです。変数が増えるとこの処理が大きな問題になる可能性があります。
一般のn変数での連立方程式の場合のFixed Hardware Ansatzはどうやって構成するのだろうか？
今回は3変数に限定しましたが、実際に$n$とした場合、Fixed Hardware Ansatzをどういった思想で構成する必要があるのかについては、特に言及はありませんでした。要検討です。
HHLアルゴリズムと異なり、スパース（行列の成分が0）でなくても安定するか？
特に言及はありませんが、アルゴリズム的には、安定するとると想像しています。
このアルゴリズムは計算量的に古典よりも効率的になりうるか？
こちらも言及されていませんでした。。。行列$A$の線形結合を見つける処理も踏まえる必要はありますし、線形結合で使われるユニタリ行列$A_{i}$の数×Hadamardテストのサンプリング回数にも依存するでしょう。
連立一次方程式には、Gaussの消去法（$O(n^3)$）やLU分解（$O(n^2)$）などの多項式時間で解ける古典アルゴリズムはすでにあるので、それらの方が優位性がある可能性は否定できないと思います。
最後に
こちらの記事の内容はQiskit Textbook 4.2.1に紹介されていたのをきっかけに解説しています。内容的にも読んでいて行間を埋めないと苦しい箇所が多かったので、それらを整理して、補足的な位置づけになるような意味で書き溜めておきました。
ちなみに、こちらの日本語訳を筆者が担当していますｗ
翻訳はほぼ直訳なので、ここで紹介した行間はあまり語られていないです。
手前味噌ながら、解説動画も公開しています。内容はこの記事とほぼ同等ですが、１時間と長いので覚悟して視聴してください！かなり丁寧には解説しているつもりです。
https://www.youtube.com/watch?v=j3gYU4-7VFc&t=59s
Qiskit Textbookは近いうちに日本語版が正式に公開されるみたいなので、是非読んでみてください！
また、このアルゴリズムの詳細が気になる方は、最初に紹介した原論文から読んでみられるとよろしいかと。