复现的算法为VQLS算法，用于处理线性方程组Ax=bAx=bAx=b的问题。
这是一种变分量子算法，即用一个经典的优化器来训练一个含参的量子电路，整个算法的思路有些类似机器学习。
代码参考了qiskit提供的VQLS示例

10.17日晚重新看了一遍论文，才发现自己弄错了很多概念，例如：误以为分解AAA是为了作为量子态输入到电路，实际上他是作为门矩阵作用于含参xxx的

流程分析

我大致将该算法分为一下4个流程。

搭建含参电路
计算损失函数
优化器进行梯度下降
求解得标准化后的x

图1 算法整体框架

1. 搭建含参电路

输入量子电路的为试探态，一般选择全∣0⟩|0⟩∣0⟩态，通过含参电路后，则可得：
∣x(α)⟩=V(α)∣0⟩|x(\alpha)⟩ = V(\alpha)|0⟩∣x(α)⟩=V(α)∣0⟩
其中V(α)V(\alpha)V(α)为含参电路的等效矩阵。

将AAA分解为U矩阵和的形式来作为门电路：A=∑cnAnA = \sum c_nA_nA=∑cnAn上式中 cnc_ncn 均为复数，这意味着可以将我们可以将AnA_nAn分别作为量子电路的门电路。

然后让带参数的x(α)x(\alpha)x(α)分别通过AnA_nAn再累加，得到Ax(α)Ax(\alpha)Ax(α)，接着就是构建一个损失函数，用于评估Ax(α)Ax(\alpha)Ax(α)与bbb的接近程度，可以理解为判别两个向量是否在同一方向上。

这一部分对用于图1中的蓝框部分。

2.损失函数

论文中提出两种损失函数，分别对应local和global。
global损失函数如下：
CG^=tr(∣ψ⟩⟨ψ∣(I−∣b⟩⟨b∣))=⟨x∣HG∣x⟩\begin{aligned}\hat{C_G}=&tr(|\psi⟩⟨\psi|(I-|b⟩⟨b|))\\=&⟨x|H_G|x⟩\end{aligned}CG^==tr(∣ψ⟩⟨ψ∣(I−∣b⟩⟨b∣))⟨x∣HG∣x⟩
其中HG=A+(I−∣b⟩⟨b∣)A，∣x⟩=V(α)∣0⟩H_G=A^+(I-|b⟩⟨b|)A，|x⟩ = V(\alpha)|0⟩HG=A+(I−∣b⟩⟨b∣)A，∣x⟩=V(α)∣0⟩
进行归一化则有：
CG=1−∣⟨b∣Ψ⟩∣2\begin{aligned}C_G=1-|⟨b|\Psi⟩|^2\end{aligned}CG=1−∣⟨b∣Ψ⟩∣2
其中∣Ψ⟩=∣ψ⟩⟨ψ∣ψ⟩|\Psi⟩ = \frac{|\psi⟩}{\sqrt{⟨\psi|\psi⟩}}∣Ψ⟩=⟨ψ∣ψ⟩∣ψ⟩。

local损失函数：CL^=⟨x∣HL∣x⟩\begin{aligned}\hat{C_L} = ⟨x|H_L|x⟩\end{aligned}CL^=⟨x∣HL∣x⟩
其中HL=A+U(I−1n∑j=1n(∣0j⟩⟨0j∣⊗Ij))U+AH_L=A^+U(I-\frac{1}{n}\sum_{j=1}^{n}(|0_j⟩⟨0_j|\otimes I_j))U^+AHL=A+U(I−n1∑j=1n(∣0j⟩⟨0j∣⊗Ij))U+A。
对其进行归一化则有：CL=CL^⟨ψ∣ψ⟩\begin{aligned}C_L=\frac{\hat{C_L}}{⟨\psi|\psi⟩}\end{aligned}CL=⟨ψ∣ψ⟩CL^

论文中提出局部损失函数相较于全局损失函数具有更好的梯度下降效果，如图2所示，局部损失函数能够在更少的迭代次数中达到梯度消失。

图2 两种损失函数对比

为了提高运算效率，这一步在量子计算机中进行，对应于图1中的红框部分。

复现过程中，我选择的是global损失函数。

3.优化器

在优化其中，主要实现的功能为：

计算损失函数
进行进行梯度下降

由步骤2分析可知，损失函数需要计算的有两个式子：
C^G=1−∣⟨b∣ψ⟩∣⟨ψ∣ψ⟩⟨ψ∣ψ⟩=⟨0∣V(α)+A+AV(α)∣0⟩=∑m∑ncm∗cn⟨0∣V(α)+Am+AnV(α)∣0⟩∣⟨b∣ψ⟩∣2=∑m∑ncm∗cn⟨0∣U+AnV(α)∣0⟩⟨0∣V(α)+AmU∣0⟩=∑m∑ncmcn⟨0∣U+AnV(α)∣0⟩⟨0∣U+AmV(α)∣0⟩\begin{aligned} \hat{C}_G&=1-\frac{|⟨b|\psi⟩|}{⟨\psi|\psi⟩}\\ ⟨\psi|\psi⟩&=⟨0|V(\alpha)^+A^+AV(\alpha)|0⟩\\ &=\sum_{m}\sum_{n}c^*_mc_n⟨0|V(\alpha)^+A^+_mA_nV(\alpha)|0⟩\\ |⟨b|\psi⟩|^2&=\sum_{m}\sum_{n}c_m^*c_n⟨0|U^+A_nV(\alpha)|0⟩⟨0|V(\alpha)^+A_mU|0⟩\\ &=\sum_{m}\sum_{n}c_mc_n⟨0|U^+A_nV(\alpha)|0⟩⟨0|U^+A_mV(\alpha)|0⟩ \end{aligned}C^G⟨ψ∣ψ⟩∣⟨b∣ψ⟩∣2=1−⟨ψ∣ψ⟩∣⟨b∣ψ⟩∣=⟨0∣V(α)+A+AV(α)∣0⟩=m∑n∑cm∗cn⟨0∣V(α)+Am+AnV(α)∣0⟩=m∑n∑cm∗cn⟨0∣U+AnV(α)∣0⟩⟨0∣V(α)+AmU∣0⟩=m∑n∑cmcn⟨0∣U+AnV(α)∣0⟩⟨0∣U+AmV(α)∣0⟩
其中U∣0⟩=∣b⟩U|0⟩=|b⟩U∣0⟩=∣b⟩

Hadamard Test

该电路主要作用是对任给的幺正算符UUU和量子态∣Ψ⟩|\Psi⟩∣Ψ⟩，可以给出该幺正算符在量子态上的投影期望 ⟨Ψ∣U∣Ψ⟩⟨\Psi|U|\Psi⟩⟨Ψ∣U∣Ψ⟩，如图3所示。

图3 Hadamard Test电路

由推导可知：
P(0)=12(1+Re⟨ψ∣U∣ψ⟩)P(1)=12(1−Re⟨ψ∣U∣ψ⟩)\begin{aligned} P(0)=\frac{1}{2}(1+Re\displaystyle \left\langle \psi|U|\psi \right\rangle)\\ P(1)=\frac{1}{2}(1-Re\displaystyle \left\langle \psi|U|\psi \right\rangle)\\ \end{aligned}P(0)=21(1+Re⟨ψ∣U∣ψ⟩)P(1)=21(1−Re⟨ψ∣U∣ψ⟩)
则有：
Re⟨ψ∣U∣ψ⟩=P(0)−P(1)\begin{aligned} Re\displaystyle \left\langle \psi|U|\psi \right\rangle=P(0)-P(1) \end{aligned}Re⟨ψ∣U∣ψ⟩=P(0)−P(1)

我们需要做的是，将x(α)x(\alpha)x(α)作为图3中的∣ψ⟩|\psi⟩∣ψ⟩输入到电路，选取AnA_nAn作为UUU，通过测量第一个qubit的0、1的概率再分别乘以系数累加则可得⟨ψ∣ψ⟩⟨\psi|\psi⟩⟨ψ∣ψ⟩

Hadamard Overlap Test

Hadamard Overlap Test电路是用于计算：
⟨0∣U+AnV(α)∣0⟩⟨0∣V(α)+AmU∣0⟩⟨0|U^+A_nV(\alpha)|0⟩⟨0|V(\alpha)^+A_mU|0⟩⟨0∣U+AnV(α)∣0⟩⟨0∣V(α)+AmU∣0⟩

图4 Hadamard Overlap Test电路

推导图4所示量子电路，可得：
P(0)−P(1)=Re(⟨0∣U+AlV(α)∣0⟩⟨0∣V(α)+Al+U∣0⟩)P(0)-P(1)=Re(⟨0|U^+A_lV(\alpha)|0⟩⟨0|V(\alpha)^+A_l^+U|0⟩)P(0)−P(1)=Re(⟨0∣U+AlV(α)∣0⟩⟨0∣V(α)+Al+U∣0⟩)
同样，通过测量第一个qubit的0、1的概率再分别乘以系数累加则可得∣⟨b∣ψ⟩∣2|⟨b|\psi⟩|^2∣⟨b∣ψ⟩∣2

梯度下降

对于global损失函数，我们对其求梯度：
∂CG^∂αi=∑m,nm∗n(∂⟨0∣V(α)+Am+AnV(α)∣0⟩∂αi−∂⟨0∣U+AnV(α)∣0⟩⟨0∣V(α)+AmU∣0⟩∣0⟩∂αi)\begin{aligned} \frac{\partial \hat{C_G}}{\partial \alpha_i}=\sum_{m,n}m^*n\displaystyle \left( \frac{\partial ⟨0|V(\alpha)^+A^+_mA_nV(\alpha)|0⟩}{\partial \alpha_i} - \frac{\partial ⟨0|U^+A_nV(\alpha)|0⟩⟨0|V(\alpha)^+A_mU|0⟩|0⟩}{\partial \alpha_i}\right) \end{aligned}∂αi∂CG^=m,n∑m∗n(∂αi∂⟨0∣V(α)+Am+AnV(α)∣0⟩−∂αi∂⟨0∣U+AnV(α)∣0⟩⟨0∣V(α)+AmU∣0⟩∣0⟩)

最后经过迭代，得到最终的参数。这一部分对应图1中的绿框部分。

4.求解得标准化后的x

经过步骤3我们得到可以得到最终的量子电路参数，所有qubit位输入∣0⟩\lvert 0⟩∣0⟩，我们得到的是标准化后的∣x⟩=x∣∣x∣∣2|x⟩=\frac{x}{||x||_2}∣x⟩=∣∣x∣∣2x

实践部分

量子电路搭建

含参电路

首先需要搭建的是含参电路，以3qubits电路，参数初始化全初始化为1为例，如图5所示。

图5 含参电路

代码如下，加入ry门以及cz门，ry门的参数作为可变参数，cz门使量子之间产生纠缠。

def apply_fixed_ansatz(circ, qubits, parameters):# 必须写为张量的形式parameters = paddle.to_tensor(parameters)parameters = paddle.reshape(parameters, [3, -1])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])

Hadamard Test电路

再含参电路的基础上，搭建不同的U矩阵，计算⟨x(α)∣Am+An∣x(α)⟩⟨x(\alpha)|A_m^+A_n|x(\alpha)⟩⟨x(α)∣Am+An∣x(α)⟩
复现过程中取A=0.45Z3+0.55IA=0.45Z_3+0.55IA=0.45Z3+0.55I
故一共需要构建四种电路，分别对应Z3、IZ_3、IZ3、I的四种组合。
生成电路如图6所示（以Am=I,An=Z3A_m=I,A_n=Z_3Am=I,An=Z3为例）：

图6 Hadamard Test

代码如下所示：

def had_test(circ, gate_type, qubits, auxiliary_index, parameters):circ.h(auxiliary_index)apply_fixed_ansatz(circ, qubits, parameters)for ie in range (0, len(gate_type[0])):if (gate_type[0][ie] == 1):circ.cz([auxiliary_index, qubits[ie]])for ie in range (0, len(gate_type[1])):if (gate_type[1][ie] == 1):circ.cz([auxiliary_index, qubits[ie]])circ.h(auxiliary_index)

Hadamard Overlap Test电路

Hadamard Overlap Test电路是对Hadamard Test的推广，用于计算∣⟨b∣ψ⟩∣2=∑m∑ncmcn⟨0∣U+AnV(α)∣0⟩⟨0∣U+AmV(α)∣0⟩|⟨b|\psi⟩|^2=\sum_{m}\sum_{n}c_mc_n⟨0|U^+A_nV(\alpha)|0⟩⟨0|U^+A_mV(\alpha)|0⟩∣⟨b∣ψ⟩∣2=m∑n∑cmcn⟨0∣U+AnV(α)∣0⟩⟨0∣U+AmV(α)∣0⟩

qiskit提供的VQLS示例中，取A为8∗8A为8*8A为8∗8矩阵，b=[0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.125]Tb=[\sqrt{0.125},\sqrt{0.125},\sqrt{0.125}, \sqrt{0.125}, \sqrt{0.125}, \sqrt{0.125}, \sqrt{0.125}, \sqrt{0.125}]^Tb=[0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.125]T
故取的U=H1H2H3U=H_1H_2H_3U=H1H2H3，如图7所示，

图7 U电路

最终生成的Hadamard Overlap Test电路如图8（由于jupyter生成的图片显示不全，选择在pycharm生成，建议paddle优化一下可视化效果）：

图8 Hadamard Overlap Test电路

同样也是对第一个qubit进行测量，然后乘系数累加则可得到∣⟨b∣ψ⟩∣2|\displaystyle \left\langle b|\psi \right\rangle|^2∣⟨b∣ψ⟩∣2
代码如下所示：

def special_had_test(circ, gate_type, qubits, auxiliary_index, parameters):circ.h(auxiliary_index)control_fixed_ansatz(circ, qubits, parameters, auxiliary_index)for ty in range (0, len(gate_type)):if (gate_type[ty] == 1):circ.cz([auxiliary_index, qubits[ty]])control_b(circ, auxiliary_index, qubits)circ.h(auxiliary_index)

优化器构建

这一部分同时参考paddle中VQE示例以及qiskit中示例，由于对双方的api都不熟悉，这一部分代码编写花费了很长时间，还存在着参数不能优化的bug以及我也无法验证loss计算过程是否存在问题。

loss计算

最开始使用paddle的优化器，由于api不熟悉，于是更换了scipy.optimize的minimize作为优化器
代码如下：

# 定义损失函数和前向传播机制
def calculate_cost_function(parameters):overall_sum_1 = 0# 参数1for i in range(0, len(gate_set)):for j in range(0, len(gate_set)):# 用hadmard-test电路计算loss的一个值cir_loss1 = UAnsatz(4)multiply = coefficient_set[i]*coefficient_set[j]had_test(cir_loss1, [gate_set[i], gate_set[j]], [1, 2, 3], 0, parameters)cir_loss1.run_state_vector()result = cir_loss1.measure(shots = 0, which_qubits = [0])# print(result['1'])overall_sum_1+=multiply*(1-(2*(result['1'])))print(overall_sum_1)overall_sum_2 = 0for i in range(0, len(gate_set)):for j in range(0, len(gate_set)):multiply = coefficient_set[i]*coefficient_set[j]mult = 1for extra in range(0, 2):cir_loss2 =  UAnsatz(5)if (extra == 0):special_had_test(cir_loss2, gate_set[i], [1, 2, 3], 0, parameters)if (extra == 1):special_had_test(cir_loss2, gate_set[j], [1, 2, 3], 0, parameters)cir_loss2.run_state_vector()result = cir_loss2.measure(shots = 0, which_qubits = [0])# print(result['1'])mult = mult*(1-(2*(result['1'])))overall_sum_2+=multiply*mult# print(overall_sum_2)loss = 1-float(overall_sum_2/overall_sum_1)# print(loss)return loss

如同前面分析

先对⟨ψ∣ψ⟩、∣⟨b∣ψ⟩∣2⟨\psi|\psi⟩、|⟨b|\psi⟩|^2⟨ψ∣ψ⟩、∣⟨b∣ψ⟩∣2分开求解：搭建两个电路，输入全0态
测量得到两个电路第一个qubit为1得概率P(1)P(1)P(1)，分别得到
⟨ψ∣ψ⟩=∑m∑ncm∗cn(1−2P(1))∣⟨b∣ψ⟩∣2=∑m∑ncmcn(1−2P(1))\begin{aligned} ⟨\psi|\psi⟩&=\sum_{m}\sum_{n}c^*_mc_n(1-2P(1))\\ |⟨b|\psi⟩|^2&=\sum_{m}\sum_{n}c_mc_n(1-2P(1)) \end{aligned}⟨ψ∣ψ⟩∣⟨b∣ψ⟩∣2=m∑n∑cm∗cn(1−2P(1))=m∑n∑cmcn(1−2P(1))
带入loss=1−∣⟨b∣ψ⟩∣2⟨ψ∣ψ⟩loss=1-\frac{|⟨b|\psi⟩|^2}{⟨\psi|\psi⟩}loss=1−⟨ψ∣ψ⟩∣⟨b∣ψ⟩∣2

参数优化器配置

out = minimize(calculate_cost_function, x0=[float(random.randint(0,3000))/1000 for i in range(0, 9)], method="COBYLA", options={'maxiter':100})
print(out)

运行结果如下：

图9 优化器运行结果

得到标准化后的x

图10 VQLS和numpy求解x对比（标准化后）

除了相差π\piπ的全局相位，我复现求解的xxx与真实标准化的xxx存在着排列顺序上以及的区别(我将示例中求解出来的电路参数带入也有相同的问题)

小结

复现过程中，发现的问题如下：

对机器学习部分基础知识掌握不到位，对paddle的api调用不熟练
量子计算理论部分还相当薄弱，对于测量、纠缠以及量子门对应的作用（尤其是多比特量子门）等部分不熟悉。

VQLS：变分量子算法解线性方程组相关推荐

matlab sor解线性方程组,SOR算法解线性方程组的matlab程序
编写用SOR 方法求解线性方程组Ax=B 的标准程序,并求下列方程组的解,并比较松弛因子取 1.0.1.25.1.5时所需迭代的次数. 可取初始向量 (0) (1,1,1) T =x ,迭代终止条件 ...
java实现高斯赛德尔算法解线性方程组
packagelinear_equation; importjava.util.Scanner; /*使用高斯赛德尔迭代法求解线性方程组*/ publicclass Gauss_Seidel_Iter ...
量子计算机算法详解,量子计算机量子算法以及物理实现.pdf
CN43-1258/TP 计算机工程与科学 2012年第34卷第8期 V01．34,No．8,2012 ISSN1007-130X ENGINEERING&SCIENCE COMPUTER 文 ...
通识：量子算法发展的四大阶段
尽管对于绝大多数人来说,"量子计算"是一个属于科幻世界的遥远词汇,但实际上它已经走过了近四十年的历史.在这四十年里,量子计算及其算法从无到有,从理论概念逐渐变成了现实.我们从量子算 ...
量子计算机个人化时间,科学家发现量子算法可以停止时间
停止时间,这对于人类来说几乎是不可能完成的任务,但对于计算机来说却是有可能的. 据外媒thenextweb报道,来自麻省理工学院和马里兰大学的两项研究表明,量子算法可以解决非线性微分方程. 比如麻省理 ...
量子计算机模拟,新量子算法将量子模拟器变成量子计算机，可以进行量子计算...
一项研究发现了一种新的量子计算通用模型,即变分模型,填补了量子模拟器与传统量子计算机之间的空白. 研究于3月10日发表在<物理评论A>上,标题为"Universal variat ...
matlab lu分解求线性方程组_计算方法（二）直接三角分解法解线性方程组
封面是WH2里春希在编辑部的上司麻理前辈,有一说一,这条线的第一次H有点恶趣味,不是很喜欢. 一:概述矩阵分解我学过的挺多种,比如极分解,谱分解,满秩分解,正交三角分解还有这里的直接三角分解大部分我 ...
超松弛迭代法解线性方程组c语言,超松弛迭代法解线性方程组.doc
PAGE PAGE 2 姓名:___________________________ 设计题目:超松弛迭代法解线性方程组专业: 摘要本文是在matlab环境下熟悉的运用计算机编程语言并结合超松弛变 ...
选主元的高斯-约旦（Gauss-Jordan）消元法解线性方程组/求逆矩阵
选主元的高斯-约当(Gauss-Jordan)消元法在很多地方都会用到,例如求一个矩阵的逆矩阵.解线性方程组(插一句:LM算法求解的一个步骤),等等.它的速度不是最快的,但是它非常稳定(来自网上的定义 ...

VQLS：变分量子算法解线性方程组