本教程復現論文Variational Quantum Linear Solver中的圖四，圖四使用了文中提出的VQLS演算法求解文中II.B.1中給出的問題Ising-inspired QLSP，給出了引數\(\kappa\)與線路運行次數的關系，
VQLS演算法用于求解線性方程的解，即對方程\(Ax=b\)，已知\(A\)和\(b\)，得出方程的解\(x\)，如上圖所示，在VQLS演算法中，作者利用量子線路來代替\(A\)，使用含參量子線路\(V(\alpha)\)來制備\(x\)，即\(|x\rangle=V(\alpha)|0\rangle\)，使用量子線路\(U\)來制備\(b\)，即\(U|0\rangle=|b\rangle\)，通過對\(V(\alpha)\)中的引數\(\alpha\)進行優化，使得\(AV|0\rangle=A|x\rangle\)逼近\(U|0\rangle=|b\rangle\)，從而得到方程的解\(x\)，
損失函式的匯出
文中使用了量子網路得到損失函式的值，并根據損失函式的值對量子網路的引數\(\alpha\)進行優化，下面對文中的損失函式進行簡要說明，
損失函式值的形式
文中通過Ansatz線路生成\(|x(\alpha)\rangle\)，將\(|\psi\rangle=A|x\rangle\)逼近\(|b\rangle\)來得到引數\(\alpha\)，
因此損失函式定義為\(C_G=1-|\langle b|\psi\rangle|^2=1-\langle\psi|b\rangle\langle b|\psi\rangle=1-\langle x|A^+U|\mathbf0\rangle\langle\mathbf0|U^+A|\psi\rangle\).
但是\(C_G\)在優化時容易遇到barren plateaus，
因此引入了新的損失函式
\(C_L=1-\langle x|A^+U(\frac{1}{n}\sum_{j=1}^n|0_j\rangle\langle0_j|\bigotimes\Bbb1_{\overline j})U^+A|x\rangle\)
代入\(|0_j\rangle\langle0_j|=\frac{Z_j+I_j}{2}\)
即
\[C_L=\frac{1}{2}-\frac{1}{2n}\sum_{j=1}^n\langle x|A^+UZ_jU^+A|x\rangle \]其中\(A\)可以分解成一系列幺正矩陣之和\(A=\sum_{l=1}^Lc_lA_l\)，其中\(A_l\)對應了簡單的門電路，方便使用Hadamard Test求解，于是\(C_L\)可以分解為
\[C_L=\frac{1}{2}-\frac{1}{2n}\sum_{ll'j}c_lc_{l'}\langle x|A^+_{l'}UZ_jU^+A_l|x\rangle \]Hadamard Test線路計算損失函式
文中提供了Hadamard Test線路和Hadamard Overlap Test線路計算上面提到的損失函式\(C_L\)，這里使用Hadamard Test線路進行計算，文中的Fig. 9(c)為對應的線路圖，
對于中間的受控\(Z_j\)，由于需要作用到每一個位元上，那么線路將多生成\(n_{qubits}\)個，我們可以使用受控PS門來代替受控\(Z_j\)，將PS門的引數以encoding的形式傳入到量子網路中，首先我們在\(Z_j\)的位置對所有位元都施加一個受控PS門，當想要對第\(j\)個位元實作\(Z_j\)作用時，我們將作用在第\(j\)個位元的\(PS\)門中傳入引數\(\pi\)，這樣受控\(PS_j\)門會變成\(Z_j\)，對于其他所有位元的\(PS_{\overline j}\)門傳入引數0，那么其他位元的的\(PS\)門會變成單位矩陣\(I_\overline j\)，即可實作想要的效果，
推導
下面簡要推導從Hadamard Test得到\(C_L\)的程序，
系統初態為\(|\mathbf0\rangle=|0\rangle_a\bigotimes|0\rangle^{\bigotimes n}\)，經過H門作用，系統變為\(\frac{1}{\sqrt2}|0\rangle_a\bigotimes|0\rangle^{\bigotimes n} + \frac{1}{\sqrt2}|1\rangle_a\bigotimes|0\rangle^{\bigotimes n}\)，
經過各種門電路作用，系統的態最終變成了
\[|\phi\rangle=\frac{1}{\sqrt2}H|0\rangle_a\bigotimes V|0\rangle^{\bigotimes n} -j \frac{1}{\sqrt2}H|1\rangle_a\bigotimes A_{l'}^+UZ_jU^+A_lV|0\rangle^{\bigotimes n}=\frac{1}{2}(|0\rangle_a+|1\rangle_a)\bigotimes V|0\rangle^{\bigotimes n} -j \frac{1}{2}(|0\rangle_a-|1\rangle_a)\bigotimes A_{l'}^+UZ_jU^+A_lV|0\rangle^{\bigotimes n} \]線路最后對ancilla位元進行了測量，得到0的概率為
\[P(|0\rangle_a)=\langle\phi|(|0\rangle_a\langle0|_a\bigotimes I)|\phi\rangle =\frac{1}{4}\langle0|^{\bigotimes n}(V^++jV^+A_l^+UZ_j^+U^+A_{l'})(V-jA_{l'}^+UZ_jU^+A_lV)|0\rangle^{\bigotimes n} =\frac{1}{2}-j\langle x|A^+_{l'}UZ_jU^+A_l|x\rangle+j\langle x|A^+_lUZ^+jU^+A_{l'}|x\rangle =\frac{1}{2}+\frac{1}{2}Im[\langle x|A^+_{l'}UZ_jU^+A_l|x\rangle]\]同理可得，得到1的概率為
\[P(|1\rangle_a) =\frac{1}{2}-\frac{1}{2}Im[\langle x|A^+_{l'}UZ_jU^+A_l|x\rangle]\]因此\(C_L\)中\(\langle x|A^+_{l'}UZ_jU^+A_l|x\rangle\)的虛部可由\(P(|0\rangle_a)-P(|1\rangle_a)\)得到，如要求實部將線路中\(S^+\)去掉，同樣使用\(P(|0\rangle_a)-P(|1\rangle_a)\)得到，
\[P(|0\rangle_a)-P(|1\rangle_a) =\langle\phi|(|0\rangle_a\langle0|_a\bigotimes I)|\phi\rangle-\langle\phi|(|1\rangle_a\langle1|_a\bigotimes I)|\phi\rangle =\langle\phi|(Z_a\bigotimes I)|\phi\rangle\]在Mindquantum中要求\(P(|0\rangle_a)-P(|1\rangle_a)\)的值，只需要使用量子線路生成\(|\phi\rangle\)后，再求哈密頓量\(Z_a\)在\(|\phi\rangle\)上的期望即可，
VQLS演算法實體
本教程使用MindQuantum首先復現https://pennylane.ai/qml/demos/tutorial_vqls.html 中的一個簡單算例，詳細說明VQLS的運行程序，
網路結構
為了得到
\[C_L=\frac{1}{2}-\frac{1}{2n}\sum_{ll'j}c_lc_{l'}\langle x|A^+_{l'}UZ_jU^+A_l|x\rangle \]其中的
\[Im[\langle x|A^+_{l'}UZ_jU^+A_l|x\rangle]=\langle\phi(ll'j)||0\rangle_a\langle 0|_a|\phi(ll'j)\rangle \]對于\(l,l'\)不同的取值，每個都要生成兩個網路，用于計算實部和虛部，\(j\)可以通過對PS門傳參控制，
本教程中，先構造網路生成\(|\phi(ll'j)\rangle\),通過對其求期望的方式構造梯度算子，從而得到了各個量子網路的輸出，再使用\(c_lc_{l'}\)加權疊加，計算得到\(C_L\)，并使用Mindspore中的自動微分功能，對\(C_L\)自動求導并優化，
構造\(A, U, V\)線路
from mindquantum.core.circuit import Circuit, UN, controlled, dagger
from mindquantum.core.gates import PhaseShift, I, RX, RY, RZ, H, X, Z
import numpy as np
# 引數定義
n_qubits = 3
# 構造A_l
A_list = [[1, Circuit()], [0.2, Circuit(X.on(0))+Z.on(1)], [0.2, Circuit(X.on(0))]]
print('A線路:')
for n, A in enumerate(A_list):
print(f'A_{n}:\n')
print(A[1],'\n')
# 構造U
U = UN(H, n_qubits)
print('U線路')
print(U, '\n')
# 構造V
V = Circuit(RY(f'omega_{i}').on(i) for i in range(n_qubits))
print('V線路')
print(V, '\n')
/home/ma-user/modelarts-dev/modelarts-sdk/requests/__init__.py:104: RequestsDependencyWarning: urllib3 (1.26.13) or chardet (5.0.0)/charset_normalizer (2.0.12) doesn't match a supported version!
RequestsDependencyWarning)
A線路:
A_0:
A_1:
q0: ──X──
q1: ──Z──
A_2:
q0: ──X──
U線路
q0: ──H──
q1: ──H──
q2: ──H──
V線路
q0: ──RY(omega_0)──
q1: ──RY(omega_1)──
q2: ──RY(omega_2)──
構造Hadmard Test線路并生成梯度算子
# 構造Hadamard Test線路
def generate_hdt(n_qubits, V, U, A_l, A_lp, is_imag=True):
ancilla_ind = n_qubits
A_l = Circuit(A_l)
A_lp = Circuit(A_lp)
ancilla = Circuit(H.on(ancilla_ind))
if is_imag:
ancilla += PhaseShift(-np.pi/2).on(ancilla_ind)
ancilla += controlled(A_l)(ancilla_ind)
hdt = Circuit()
hdt += ancilla
hdt += dagger(U)
for j in range(0, n_qubits):
hdt += controlled(Circuit(PhaseShift(f'Z_{j}').on(j)))(ancilla_ind)
hdt += U
hdt += controlled(dagger(A_lp))(ancilla_ind)
hdt += H.on(ancilla_ind)
hdt = hdt.no_grad()
hdt.as_encoder()
return V + hdt
hdt_imag = generate_hdt(n_qubits, V, U, A_list[0][1], A_list[0][1])
hdt_real = generate_hdt(n_qubits, V, U, A_list[0][1], A_list[0][1], False)
print(hdt_imag)
hdt_imag.summary()
print(hdt_real)
hdt_real.summary()
# 構造梯度算子以得到測量概率
from mindquantum.core.operators import Hamiltonian, QubitOperator
from mindquantum.simulator import Simulator
# 構造哈密頓量
# 當ancilla測量結果為1，波函式系數乘-1；當測量結果為0，系數不變.
Qops = QubitOperator('Z%d'%n_qubits)
hams = Hamiltonian(Qops)
sim = Simulator('projectq', n_qubits+1)
# 得到<0|hams hdt_imag|0>=P(0)-P(1)
grad_ops_imag = sim.get_expectation_with_grad(hams, hdt_imag)
grad_ops_real = sim.get_expectation_with_grad(hams, hdt_real)
# 運行多次均為實數，證明可以得到一個概率
weight = np.pi*np.array(np.random.rand(n_qubits))
print(grad_ops_imag(np.pi*np.array([[1.,0, 0]]), weight))
print(grad_ops_real(np.pi*np.array([[1.,0, 0]]), weight))
q0: ──RY(omega_0)───────H────────PS(Z_0)───────H─────────────────────
│
q1: ──RY(omega_1)───────H───────────┼───────PS(Z_1)───────H──────────
│ │
q2: ──RY(omega_2)───────H───────────┼──────────┼───────PS(Z_2)────H──
│ │ │
q3: ───────H─────────PS(-π/2)───────●──────────●──────────●───────H──
==========================Circuit Summary==========================
|Total number of gates : 15. |
|Parameter gates : 6. |
|with 6 parameters are : |
|omega_0, omega_1, omega_2, Z_0, Z_1, Z_2 .|
|Number qubit of circuit: 4 |
===================================================================
q0: ──RY(omega_0)────H────PS(Z_0)───────H─────────────────────
│
q1: ──RY(omega_1)────H───────┼───────PS(Z_1)───────H──────────
│ │
q2: ──RY(omega_2)────H───────┼──────────┼───────PS(Z_2)────H──
│ │ │
q3: ───────H─────────────────●──────────●──────────●───────H──
==========================Circuit Summary==========================
|Total number of gates : 14. |
|Parameter gates : 6. |
|with 6 parameters are : |
|omega_0, omega_1, omega_2, Z_0, Z_1, Z_2 .|
|Number qubit of circuit: 4 |
===================================================================
(array([[0.+0.j]]), array([[[0.+0.j, 0.+0.j, 0.+0.j]]]), array([[[-1.23010388e-17+0.j, 4.23704588e-33+0.j, -1.54074396e-33+0.j]]]))
(array([[0.56641042+0.j]]), array([[[0.+0.j, 0.+0.j, 0.+0.j]]]), array([[[ 8.24123313e-01+0.j, 2.77555756e-17+0.j, -1.77809156e-17+0.j]]]))
# 對每對$A_l$, $A_{l'}^+$生成量子線路，得到所有線路對應的梯度算子
def generate_grad_ops(n_qubits, V, U, A):
circuits_list = []
grad_ops_list = []
coeff_list = []
# 構造哈密頓量
Qops = QubitOperator('Z%d'%n_qubits)
hams = Hamiltonian(Qops)
# 生成模擬器
sim = Simulator('projectq', n_qubits+1)
for c_l, A_l in A:
for c_lp, A_lp in A:
hdt_imag = generate_hdt(n_qubits, V, U, A_l, A_lp)
hdt_real = generate_hdt(n_qubits, V, U, A_l, A_lp, False)
# 得到<0|hams hdt_imag|0>=P(0)-P(1)
grad_ops_imag = sim.get_expectation_with_grad(hams, hdt_imag)
grad_ops_real = sim.get_expectation_with_grad(hams, hdt_real)
grad_ops_list.extend([grad_ops_real, grad_ops_imag])
coeff_list.append(c_l*np.conj(c_lp))
circuits_list.append(hdt_real)
return coeff_list, grad_ops_list, circuits_list
coeff_list, grad_ops_list, circuits_list = generate_grad_ops(n_qubits, V, U, A_list)
# 輸出所有網路檢查
for ind, circuit in enumerate(circuits_list):
print(f'第{ind+1}個網路')
print(circuit, '\n')
第1個網路
q0: ──RY(omega_0)────H────PS(Z_0)───────H─────────────────────
│
q1: ──RY(omega_1)────H───────┼───────PS(Z_1)───────H──────────
│ │
q2: ──RY(omega_2)────H───────┼──────────┼───────PS(Z_2)────H──
│ │ │
q3: ───────H─────────────────●──────────●──────────●───────H──
第2個網路
q0: ──RY(omega_0)────H────PS(Z_0)───────H────────────────────────────X───────
│ │
q1: ──RY(omega_1)────H───────┼───────PS(Z_1)───────H────────────Z────┼───────
│ │ │ │
q2: ──RY(omega_2)────H───────┼──────────┼───────PS(Z_2)────H────┼────┼───────
│ │ │ │ │
q3: ───────H─────────────────●──────────●──────────●────────────●────●────H──
第3個網路
q0: ──RY(omega_0)────H────PS(Z_0)───────H───────────────────────X───────
│ │
q1: ──RY(omega_1)────H───────┼───────PS(Z_1)───────H────────────┼───────
│ │ │
q2: ──RY(omega_2)────H───────┼──────────┼───────PS(Z_2)────H────┼───────
│ │ │ │
q3: ───────H─────────────────●──────────●──────────●────────────●────H──
第4個網路
q0: ──RY(omega_0)────X────H─────────PS(Z_0)───────H─────────────────────
│ │
q1: ──RY(omega_1)────┼────Z────H───────┼───────PS(Z_1)───────H──────────
│ │ │ │
q2: ──RY(omega_2)────┼────┼────H───────┼──────────┼───────PS(Z_2)────H──
│ │ │ │ │
q3: ───────H─────────●────●────────────●──────────●──────────●───────H──
第5個網路
q0: ──RY(omega_0)────X────H─────────PS(Z_0)───────H────────────────────────────X───────
│ │ │
q1: ──RY(omega_1)────┼────Z────H───────┼───────PS(Z_1)───────H────────────Z────┼───────
│ │ │ │ │ │
q2: ──RY(omega_2)────┼────┼────H───────┼──────────┼───────PS(Z_2)────H────┼────┼───────
│ │ │ │ │ │ │
q3: ───────H─────────●────●────────────●──────────●──────────●────────────●────●────H──
第6個網路
q0: ──RY(omega_0)────X────H─────────PS(Z_0)───────H───────────────────────X───────
│ │ │
q1: ──RY(omega_1)────┼────Z────H───────┼───────PS(Z_1)───────H────────────┼───────
│ │ │ │ │
q2: ──RY(omega_2)────┼────┼────H───────┼──────────┼───────PS(Z_2)────H────┼───────
│ │ │ │ │ │
q3: ───────H─────────●────●────────────●──────────●──────────●────────────●────H──
第7個網路
q0: ──RY(omega_0)────X────H────PS(Z_0)───────H─────────────────────
│ │
q1: ──RY(omega_1)────┼────H───────┼───────PS(Z_1)───────H──────────
│ │ │
q2: ──RY(omega_2)────┼────H───────┼──────────┼───────PS(Z_2)────H──
│ │ │ │
q3: ───────H─────────●────────────●──────────●──────────●───────H──
第8個網路
q0: ──RY(omega_0)────X────H────PS(Z_0)───────H────────────────────────────X───────
│ │ │
q1: ──RY(omega_1)────┼────H───────┼───────PS(Z_1)───────H────────────Z────┼───────
│ │ │ │ │
q2: ──RY(omega_2)────┼────H───────┼──────────┼───────PS(Z_2)────H────┼────┼───────
│ │ │ │ │ │
q3: ───────H─────────●────────────●──────────●──────────●────────────●────●────H──
第9個網路
q0: ──RY(omega_0)────X────H────PS(Z_0)───────H───────────────────────X───────
│ │ │
q1: ──RY(omega_1)────┼────H───────┼───────PS(Z_1)───────H────────────┼───────
│ │ │ │
q2: ──RY(omega_2)────┼────H───────┼──────────┼───────PS(Z_2)────H────┼───────
│ │ │ │ │
q3: ───────H─────────●────────────●──────────●──────────●────────────●────H──
構建量子網路和損失函式\(C_L\)
由于所有梯度算子共享同一套優化引數omega_1, omega_2, omega_3，我們需要把這些梯度算子封裝到一層神經網路中，量子網路的輸出為一個list，list中每一項都為一個Hadamard Test線路的輸出，
圖中所示的Parameters，就是需要優化的引數omega_1, omega_2, omega_3，所有量子網路共享這些引數，以便參與到對\(C_L\)的求導優化，這是通過下面代碼構造的MQMultiGradsLayer類來完成的，
# 根據梯度算子構造神經網路，以便參與到優化
# 由于所有梯度算子共享同一套優化引數omega_1, omega_2, omega_3，我們需要把這些梯度算子封裝到一層神經網路中
import mindspore as ms
ms.context.set_context(mode=ms.context.PYNATIVE_MODE, device_target="CPU")
from mindspore import nn
from mindquantum.framework import MQOps
from mindspore.common.initializer import initializer
from mindspore import ops, Tensor, Parameter
class MQMultiGradsLayer(nn.Cell):
def __init__(self, grads_ops_list, weight='normal'):
"""Initialize a MQLayer object."""
super().__init__()
self.evo_list = []
weight_name_set = set()
# 相同名稱的ansatz為同一個權重
ansatz_params_list = []
for grads in grads_ops_list:
self.evo_list.append(MQOps(grads))
# print(grads.circ_right)
weight_name_set = weight_name_set | set(grads.ansatz_params_name)
ansatz_params_list.append(grads.ansatz_params_name)
weight_name_list = list(weight_name_set)
weight_name_list = np.sort(weight_name_list).tolist()
self.weight_ind_list = []
for params_name in ansatz_params_list:
weight_ind = [weight_name_list.index(p_name) for p_name in params_name]
self.weight_ind_list.append(Tensor(weight_ind, dtype=ms.int32))
weight_size = len(weight_name_list)
if isinstance(weight, ms.Tensor):
if weight.ndim != 1 or weight.shape[0] != weight_size:
raise ValueError(f"Weight init shape error, required ({weight_size}, ), but get {weight.shape}.")
self.weight = Parameter(initializer(weight, weight_size, dtype=ms.float32), name='ansatz_weight')
def construct(self, x):
result = []
for evo, weight_ind in zip(self.evo_list, self.weight_ind_list):
# print(evo.expectation_with_grad.circ_right)
weight = ops.gather(self.weight, weight_ind, axis=0) # 挑選網路中需要的的引數
result.append(evo(x, weight))
return result
# 生成網路并嘗試運行，此時我們會得到2x3x3個網路的輸出
QuantumNets = MQMultiGradsLayer(grad_ops_list)
result = QuantumNets(Tensor(np.pi*np.array([[1.,0, 0], [0, 0, 1.]])))
print('共%d個網路' % len(result))
print(result)
共18個網路
[Tensor(shape=[2, 1], dtype=Float32, value=
https://www.cnblogs.com/liniganma/archive/2023/04/16/[[-8.28392163e-04],
[ 1.70038268e-02]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-5.51035729e-17],
[-2.74975515e-17]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 9.99993265e-01],
[-1.40857428e-05]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-1.47796430e-17],
[-2.82190218e-17]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 1.00000000e+00],
[-1.40858374e-05]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-1.47796430e-17],
[-2.82190218e-17]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 9.99993265e-01],
[-1.40857428e-05]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-1.50369954e-17],
[-2.83241233e-17]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-8.28392163e-04],
[ 1.70038268e-02]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 5.50419082e-17],
[ 1.39000767e-16]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-8.28386634e-04],
[ 1.70037132e-02]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 5.50419082e-17],
[ 1.39000767e-16]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 1.00000000e+00],
[-1.40858374e-05]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-1.39317067e-17],
[-2.78645980e-17]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-8.28386634e-04],
[ 1.70037132e-02]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 5.54116745e-17],
[ 1.38983191e-16]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[-8.28392163e-04],
[ 1.70038268e-02]]), Tensor(shape=[2, 1], dtype=Float32, value=
[[ 5.54116745e-17],
[ 1.38983191e-16]])]
# 構建C_L
class NetworkCL(nn.Cell):
def __init__(self, n_qubits, grad_ops_list, coeff_list):
super().__init__()
self.vqcs = MQMultiGradsLayer(grad_ops_list)
self.coeff_list = Tensor(coeff_list)
self.n_qubits = n_qubits
def construct(self):
x = np.pi*np.eye(n_qubits)
x = np.row_stack((x, [0]*n_qubits))
x = Tensor(x)
qnets = self.vqcs(x)
# 此時qnets按照u_{ll',j}的實部和虛部依次排布
real = Tensor(0)
imag = Tensor(0)
for n in range(0, len(self.coeff_list)):
real += qnets[2*n]*self.coeff_list[n]
imag += qnets[2*n+1]*self.coeff_list[n]
# 歸一化
norm = ops.pow(ops.square(sum(real[0:n_qubits]))+ops.square(sum(imag[0:n_qubits])), 0.5)
norm /= ops.pow(ops.square(real[-1])+ops.square(imag[-1]), 0.5)
return 0.5 - 0.5/n_qubits*norm
CL = NetworkCL(n_qubits, grad_ops_list, coeff_list)
print('初始C_L:', CL())
print('優化引數為:', CL.trainable_params())
初始C_L: [0.38568467]
優化引數為: [Parameter (name=vqcs.weight, shape=(3,), dtype=Float32, requires_grad=True)]
優化\(C_L\)
%matplotlib inline
import matplotlib.pyplot as plt
import mindspore
from mindspore import ops, Tensor, Parameter
from mindspore.nn import TrainOneStepCell
# 優化
learning_rate = 2
optimizer = nn.SGD(CL.trainable_params(), learning_rate=learning_rate)
net = TrainOneStepCell(CL, optimizer)
cost_history = []
for it in range(30):
res = net()
cost_history.append(float(res))
plt.plot(np.arange(30), cost_history)
plt.xlabel('Optimization Steps')
plt.ylabel('Cost Function')
結果驗證
# 得到x并驗證結果
from mindquantum.core.gates import I
I_on_all = sum(Circuit(I.on(n)) for n in range(n_qubits))
sim = Simulator('projectq', n_qubits)
sim.apply_circuit(V, CL.trainable_params()[0].asnumpy())
x = sim.get_qs()
sim.reset()
b = sim.apply_circuit(U)
b = sim.get_qs()
A = sum([c*(A_l+I_on_all).matrix(big_end=True) for c, A_l in A_list])
x_classical = np.dot(np.linalg.inv(A.astype(np.float64)), b)
x_classical /= np.linalg.norm(x_classical)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 4))
ax1.bar(np.arange(0, 2 ** n_qubits), np.abs(x_classical)**2, color="blue")
ax1.set_xlim(-0.5, 2 ** n_qubits - 0.5)
ax1.set_xlabel("Vector space basis")
ax1.set_title("Classical probabilities")
ax2.bar(np.arange(0, 2 ** n_qubits), np.abs(x)**2, color="green")
ax2.set_xlim(-0.5, 2 ** n_qubits - 0.5)
ax2.set_xlabel("Hilbert space basis")
ax2.set_title("Quantum probabilities")
plt.show()
/home/ma-user/anaconda3/envs/MindQuantum/lib/python3.7/site-packages/ipykernel/__main__.py:13: ComplexWarning: Casting complex values to real discards the imaginary part
原文復現
構造\(A\)
原文中要求\(J=0.1\)，\(\xi\)和\(\eta\)的選擇應滿足讓\(A\)的最大特征值為1，最小特征值為\(1/\kappa\)，其中\(\kappa\)是給定的引數，
假設\(A\)中的前兩個求和項的\(\sum_{j=1}^nX_j+J\sum_{j=1}^{n-1}Z_jZ_{j+1}\)最大最小特征值為\(E_{max}, E_{min}\)，\(\Bbb1\)為單位陣，其特征值為1，因此\(A\)的最大、最小特征值分別為\(\frac{E_{max}+\eta}{\xi}=1\)，\(\frac{E_{min}+\eta}{\xi}=\frac{1}{\kappa}\)，
我們首先計算\(E_{max}\)和\(E_{min}\)，然后代入上式計算\(\xi\)和\(\eta\)，
n_qubits = 10
def generate_A_list(kappa, n_qubits):
J = 0.1
I_on_all = sum(Circuit(I.on(n)) for n in range(n_qubits))
A_list = [[1, Circuit(X.on(j))] for j in range(n_qubits)]
A_list.extend([[J, Circuit(Z.on(j)) + Z.on(j+1)] for j in range(n_qubits-1)])
A = sum([c*(A_l+I_on_all).matrix(big_end=True) for c, A_l in A_list])
eigvals = np.linalg.eigvals(A)
xi = (max(eigvals)-min(eigvals))/(1-1/kappa)
eta = (max(eigvals)-kappa*min(eigvals))/(kappa-1)
A_list.append([eta, Circuit()])
A_list = [[float(c/xi), A_l] for c, A_l in A_list]
return A_list
# 求$\kappa=20$時的$A$
kappa = 20
A_list = generate_A_list(kappa, n_qubits)
I_on_all = sum(Circuit(I.on(n)) for n in range(n_qubits))
A = sum([c*(A_l+I_on_all).matrix(big_end=True) for c, A_l in A_list])
eigvals = np.linalg.eigvals(A)
print('最大特征值為%f, 最小特征值為%f' % (max(eigvals), min(eigvals)))
/home/ma-user/anaconda3/envs/MindQuantum/lib/python3.7/site-packages/ipykernel/__main__.py:13: ComplexWarning: Casting complex values to real discards the imaginary part
最大特征值為1.000000, 最小特征值為0.050000
/home/ma-user/anaconda3/envs/MindQuantum/lib/python3.7/site-packages/ipykernel/__main__.py:22: ComplexWarning: Casting complex values to real discards the imaginary part
構造\(V\)
文中使用的是Hardware-Efficient Ansatz，
from mindquantum.core.circuit import shift, add_prefix, Circuit
from mindquantum.core.gates import RY, X
n_layers = 4
n_qubits = 10
def construct_ansatz(n_layers, n_qubits):
ansatz = Circuit()
ansatz += Circuit(RY(f'0_ry{i}').on(i) for i in range(n_qubits))
# 構造每一層的線路模板
template = Circuit([X.on(1, 0)])
single_layer = sum(shift(template, i) for i in range(0, n_qubits, 2))
single_layer += Circuit(RY(f'1_ry{i}').on(i) for i in range(n_qubits))
single_layer += sum(shift(template, i) for i in range(1, n_qubits-1, 2))
single_layer += Circuit(RY(f'2_ry{i}').on(i) for i in range(1, n_qubits-1))
# 總共有n_layers層
ansatz += sum(add_prefix(single_layer, f'layer{i}') for i in range(0, n_layers))
return ansatz
eff_ansatz = construct_ansatz(n_layers, n_qubits)
print(eff_ansatz)
q0: ──RY(0_ry0)────●────RY(layer0_1_ry0)─────────────────────────────●────RY(layer1_1_ry0)─────────────────────────────●────RY(layer2_1_ry0)─────────────────────────────●────RY(layer3_1_ry0)───────────────────────────
│ │ │ │
q1: ──RY(0_ry1)────X────RY(layer0_1_ry1)────●────RY(layer0_2_ry1)────X────RY(layer1_1_ry1)────●────RY(layer1_2_ry1)────X────RY(layer2_1_ry1)────●────RY(layer2_2_ry1)────X────RY(layer3_1_ry1)────●────RY(layer3_2_ry1)──
│ │ │ │
q2: ──RY(0_ry2)────●────RY(layer0_1_ry2)────X────RY(layer0_2_ry2)────●────RY(layer1_1_ry2)────X────RY(layer1_2_ry2)────●────RY(layer2_1_ry2)────X────RY(layer2_2_ry2)────●────RY(layer3_1_ry2)────X────RY(layer3_2_ry2)──
│ │ │ │
q3: ──RY(0_ry3)────X────RY(layer0_1_ry3)────●────RY(layer0_2_ry3)────X────RY(layer1_1_ry3)────●────RY(layer1_2_ry3)────X────RY(layer2_1_ry3)────●────RY(layer2_2_ry3)────X────RY(layer3_1_ry3)────●────RY(layer3_2_ry3)──
│ │ │ │
q4: ──RY(0_ry4)────●────RY(layer0_1_ry4)────X────RY(layer0_2_ry4)────●────RY(layer1_1_ry4)────X────RY(layer1_2_ry4)────●────RY(layer2_1_ry4)────X────RY(layer2_2_ry4)────●────RY(layer3_1_ry4)────X────RY(layer3_2_ry4)──
│ │ │ │
q5: ──RY(0_ry5)────X────RY(layer0_1_ry5)────●────RY(layer0_2_ry5)────X────RY(layer1_1_ry5)────●────RY(layer1_2_ry5)────X────RY(layer2_1_ry5)────●────RY(layer2_2_ry5)────X────RY(layer3_1_ry5)────●────RY(layer3_2_ry5)──
│ │ │ │
q6: ──RY(0_ry6)────●────RY(layer0_1_ry6)────X────RY(layer0_2_ry6)────●────RY(layer1_1_ry6)────X────RY(layer1_2_ry6)────●────RY(layer2_1_ry6)────X────RY(layer2_2_ry6)────●────RY(layer3_1_ry6)────X────RY(layer3_2_ry6)──
│ │ │ │
q7: ──RY(0_ry7)────X────RY(layer0_1_ry7)────●────RY(layer0_2_ry7)────X────RY(layer1_1_ry7)────●────RY(layer1_2_ry7)────X────RY(layer2_1_ry7)────●────RY(layer2_2_ry7)────X────RY(layer3_1_ry7)────●────RY(layer3_2_ry7)──
│ │ │ │
q8: ──RY(0_ry8)────●────RY(layer0_1_ry8)────X────RY(layer0_2_ry8)────●────RY(layer1_1_ry8)────X────RY(layer1_2_ry8)────●────RY(layer2_1_ry8)────X────RY(layer2_2_ry8)────●────RY(layer3_1_ry8)────X────RY(layer3_2_ry8)──
│ │ │ │
q9: ──RY(0_ry9)────X────RY(layer0_1_ry9)─────────────────────────────X────RY(layer1_1_ry9)─────────────────────────────X────RY(layer2_1_ry9)─────────────────────────────X────RY(layer3_1_ry9)───────────────────────────
優化終止條件
根據文中式(9)的定義，可以將\(\gamma=\frac{\epsilon^2}{n_{qubits}\kappa^2}\)作為\(C_L\)的下界，當\(C_L\le\gamma\)時終止優化程式，得到問題的解，
# 構造U
U = UN(H, n_qubits)
# 構造V
eff_ansatz = construct_ansatz(n_layers, n_qubits)
# 構造所有線路
coeff_list, grad_ops_list, circuits_list = generate_grad_ops(n_qubits, eff_ansatz, U, A_list)
# 構造損失函式
CL = NetworkCL(n_qubits, grad_ops_list, coeff_list)
# 優化
learning_rate = 5
epsilon = 10
optimizer = nn.SGD(CL.trainable_params(), learning_rate=learning_rate)
net = TrainOneStepCell(CL, optimizer)
cost_history = []
time_cost = 0
gamma = 1/n_qubits*(epsilon/kappa)**2
res = gamma + 1
while gamma < res:
res = net()
cost_history.append(float(res))
print('Loss: ', float(res))
time_cost += 1
plt.plot(np.arange(len(cost_history)), cost_history)
plt.xlabel('Optimization Steps')
plt.ylabel('Cost Function')
print('Time Cost:', time_cost)
Loss: 0.42248088121414185
Loss: 0.21629470586776733
Loss: 0.09381318092346191
Loss: 0.04756203293800354
Loss: 0.03003239631652832
Loss: 0.021059304475784302
Time Cost: 6
Time-to-solution的計算
圖四描述了\(\kappa\)取值與Time-to-solution的關系，文中的Time-to-solution為計算損失函式的次數(the number of exact cost function evaluations during the optimization)，
def time_to_solution(A_list, n_qubits, epsilon):
n_layers = 4
U = UN(H, n_qubits)
eff_ansatz = construct_ansatz(n_layers, n_qubits)
coeff_list, grad_ops_list, circuits_list = generate_grad_ops(n_qubits, eff_ansatz, U, A_list)
CL = NetworkCL(n_qubits, grad_ops_list, coeff_list)
learning_rate = 5
optimizer = nn.SGD(CL.trainable_params(), learning_rate=learning_rate)
net = TrainOneStepCell(CL, optimizer)
time_to_solution = 0
gamma = 1/n_qubits*(epsilon/kappa)**2
res = gamma + 1
while gamma < res:
res = net()
time_to_solution += 1
return time_to_solution
Fig.4復現
\(\kappa\)取不同值時，我們以\(\kappa=20\)為基準，對\(A\)中每項的系數進行放縮，以滿足其最大特征值和最小特征值為\(1,1/\kappa\)，由于\(n_{qubits}>10\)時運行時間過長，這里僅復現\(n_{qubits}=10\)的情況，這里取\(\epsilon\)為50和10，以讓程式更快執行完畢，
epsilon_list = (50, 10)
all_time_cost = []
n_qubits = 10
kappa_list = np.arange(20, 180, 20)
A_list = generate_A_list(kappa_list[0], n_qubits)
for epsilon in epsilon_list:
time_cost = []
print('epsilon:', epsilon)
for kappa in kappa_list:
# 放縮A_list
xi = (1-1/kappa_list[0])/(1-1/kappa)
eta = (1-kappa*1/kappa_list[0])/(kappa-1)
A_resized = [[float(c/xi), A_l] for c, A_l in A_list]
A_resized[-1][0] += eta/xi
time_cost.append(time_to_solution(A_resized, n_qubits, epsilon))
print('kappa:', kappa, 'Time-to-solution:', time_cost[-1])
# A = sum([c*(A_l+I_on_all).matrix(big_end=True) for c, A_l in A_resized])
# eigvals = np.linalg.eigvals(A)
# print('最大特征值為%f, 最小特征值為%f' % (max(eigvals), min(eigvals)))
all_time_cost.append(time_cost)
all_time_cost = np.asarray(all_time_cost)
/home/ma-user/anaconda3/envs/MindQuantum/lib/python3.7/site-packages/ipykernel/__main__.py:13: ComplexWarning: Casting complex values to real discards the imaginary part
epsilon: 50
kappa: 20 Time-to-solution: 1
kappa: 40 Time-to-solution: 3
kappa: 60 Time-to-solution: 4
kappa: 80 Time-to-solution: 5
kappa: 100 Time-to-solution: 6
kappa: 120 Time-to-solution: 7
kappa: 140 Time-to-solution: 8
kappa: 160 Time-to-solution: 9
epsilon: 10
kappa: 20 Time-to-solution: 6
kappa: 40 Time-to-solution: 12
kappa: 60 Time-to-solution: 19
kappa: 80 Time-to-solution: 26
kappa: 100 Time-to-solution: 36
kappa: 120 Time-to-solution: 46
kappa: 140 Time-to-solution: 54
kappa: 160 Time-to-solution: 67
%matplotlib inline
plt.plot(kappa_list, all_time_cost[0, :], 'o--')
plt.plot(kappa_list, all_time_cost[1, :], 'o-')
plt.legend([r'$\epsilon=$'+str(eps) for eps in epsilon_list])
plt.xlabel(r'$\kappa$')
plt.ylabel('Time-to-solution')
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/550260.html
標籤：其他