Variational Quantum Linear Solver 的MindQuantum復現

本教程復現論文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=
[[-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')