暂无图片
暂无图片
暂无图片
暂无图片
暂无图片

项目经验分享:保真度算子开发

开源之夏 2021-11-02
627

暑期2021优秀学生的评选已经正式启动了,期待今年闪耀的新星!开源之夏公众号持续欢迎大家分享项目经验,开源心得!



投稿方式:

E-mail :summer@iscas.ac.cn

or 关注公众号,后台回复“投稿”


本期给大家分享的是来自MindSpore社区闫家舜同学的经验分享:保真度算子开发-计算得到量子神经网络演化出的量子态与给定的一个或多个量子态的保真度并能够进行反向传播



摘要
ABSTRACT
本报告为2021暑假开源软件供应链点亮计划项目的的总结报告,项目编号:210610351,项目名称:保真度算子开发。主要完成的任务:第一,在Linux下对mindquantum通过一个更普适的算符测量算子实现了保真度计算功能并和梯度的反向传播;第二,在mindquantum下实现量子相位估计算法的演示。本报告分为三个部分介绍,第一部分结合后端c++代码阐述算子具体实现原理;第二部分在前端python层搭建测试量子线路,通过一个已知结果的末态对算子进行测试和说明;第三部分为独立部分,主要为利用mindquantum现有API实现量子估计算法。本部分作为mindquantum实例教程之一已经并入开源主仓。


保真度算子实现原理

  • 保真度计算

在实际的量子计算中,通常会遇到测量问题,即对算符求某量子态下的期望,例如VQE中需要得到在含参量子态 下的哈密顿量的期望有时问题需要更一般的测量,这时测量左失不等同于右矢 ,且被测量算符为非厄米算符,例如对非厄米哈密顿量进行测量以得到其非对角元 所以,在mindquantum中提供一个算子实现任意算符的基于任意态矢下的测量是有意义的。如果这个算子可以实现,那么计算量子线路末态和目标态的保真度就变成这个测量算子的一个特例,只需要规定测量右矢为初态,被测算符为,测量右矢为演化末态,既可实现保真度计算 或者规定两个测量态矢一个为目标态,一个为线路演化末态 ,被测量算符为单位算符 同样可以实现

总的来说,我们实现了一个更普适的算符测量算子,通过规定测量态矢和算符,保真度计算功能可以当做这个算子的特例得以实现。


  • 梯度计算
一个普适的测量可以表示为

其中 为两个不同的参数化量子线路产生的末态, 为含参量子门,所以测量结果对某一参数的偏导为
对测量右矢的参数求导可以通过迭代依次求得,例如若计算
(a).初始赋值左失,赋值右矢
(b).更新右矢
(c).计算偏微分
(d).更新左失
(f).继续对 求偏导则可重复(b)~(d),以此类推直至右矢所有参数求导完成。
对左失参数的偏微分计算同理,具体实现中,我们将一侧的偏微分计算封装为一个模块 OneSideGrad()
,后端代码为
void OnesideGrad(Projectq<T> &sim_no_grad_side, Projectq<T> &sim_grad_side, 
const VT<BaseGate<T>> &circ, const VT<BaseGate<T>> &herm_circ,
const ParameterResolver<T> &pr, MST<size_t> &p_map, VT<CT<T>> &f_g, bool symmetry = true) {
for (size_t j = 0; j < circ.size(); j++) {
if ((!herm_circ[j].parameterized_) || (herm_circ[j].params_.requires_grad_parameters_.size() == 0)) {
if (herm_circ[j].parameterized_) {
sim_no_grad_side.ApplyGate(herm_circ[j], pr, false);
sim_grad_side.ApplyGate(herm_circ[j], pr, false);
}
else {
sim_no_grad_side.ApplyGate(herm_circ[j]);
sim_grad_side.ApplyGate(herm_circ[j]);
}
}
else {
sim_grad_side.ApplyGate(herm_circ[j], pr, false);
sim_grad_side.run();
Projectq<T> sim_grad_side_tmp = Projectq<T>(1, n_qubits_, sim_grad_side.vec_);
sim_grad_side_tmp.ApplyGate(circ[circ.size() - j - 1], pr, true);
sim_grad_side_tmp.run();
sim_no_grad_side.run();
CT<T> gi = 0;
if (herm_circ[j].ctrl_qubits_.size() == 0) {
gi = ComplexInnerProduct<T, calc_type>(sim_no_grad_side.vec_, sim_grad_side_tmp.vec_, static_cast<Index>(len_));
}
else {
gi = ComplexInnerProductWithControl<T, calc_type>(
sim_no_grad_side.vec_, sim_grad_side_tmp.vec_, static_cast<Index>(len_), GetControlMask(herm_circ[j].ctrl_qubits_));
}
if (symmetry){
for (auto &it : herm_circ[j].params_.requires_grad_parameters_){
f_g[1 + p_map[it]] -= 2 * herm_circ[j].params_.data_.at(it) * gi;
}
}
else{
for (auto &it : herm_circ[j].params_.requires_grad_parameters_){
f_g[1 + p_map[it]] -= herm_circ[j].params_.data_.at(it) * gi;
}
}
sim_no_grad_side.ApplyGate(herm_circ[j], pr, false);
}
}
}

基于此,一个广义的测量和对应梯度的返回可以以下代码实现 
 VT<VT<CT<T>>> NonHermitianMeasureWithGrad(const VT<Hamiltonian<T>> &hams, const VT<Hamiltonian<T>> &herm_hams,
const VT<BaseGate<T>> &left_circ, const VT<BaseGate<T>> &herm_left_circ,
const VT<BaseGate<T>> &right_circ, const VT<BaseGate<T>> &herm_right_circ,
const ParameterResolver<T> &pr, const VT<std::string> &params_order,
size_t mea_threads) {
auto n_hams = hams.size();
auto n_params = pr.data_.size();
bool symmetry = false;
MST<size_t> p_map;
for (size_t i = 0; i < params_order.size(); i++) {
p_map[params_order[i]] = i;
}
VT<VT<CT<T>>> output(n_hams);
Projectq<T> sim = Projectq<T>(1, n_qubits_, vec_);
sim.ApplyCircuit(right_circ, pr);
Projectq<T> sim2 = Projectq<T>(1, n_qubits_, vec_);
sim2.ApplyCircuit(left_circ, pr);
#pragma omp parallel for schedule(static) num_threads(mea_threads)
for (size_t i = 0; i < n_hams; i++) {
VT<CT<T>> f_g(n_params + 1, 0);
{
Projectq<T> sim_left = Projectq<T>(1, n_qubits_, sim2.vec_);
sim_left.ApplyHamiltonian(herm_hams[i]);
f_g[0] = ComplexInnerProduct<T, calc_type>(sim.vec_, sim_left.vec_, static_cast<Index>(len_));
Projectq<T> sim_right = Projectq<T>(1, n_qubits_, sim.vec_);
Projectq::OnesideGrad(sim_left, sim_right, right_circ, herm_right_circ, pr, p_map, f_g, symmetry);
}
{
Projectq<T> sim_right = Projectq<T>(1, n_qubits_, sim.vec_);
sim_right.ApplyHamiltonian(hams[i]);
f_g[0] = ComplexInnerProduct<T, calc_type>(sim2.vec_, sim_right.vec_, static_cast<Index>(len_));
Projectq<T> sim_left = Projectq<T>(1, n_qubits_, sim2.vec_);
Projectq::OnesideGrad(sim_right, sim_left, left_circ, herm_left_circ, pr, p_map, f_g, symmetry);
}
output[i] = f_g;
}
return output;
}


  • 保真度算子实例演示
    保护

  • 实例说明
在本次说明中,我们将对Simulator类中的成员函数nonhermitian_measure_with_grad()进行说明,并通过一个简单的两比特激发转换实例进行演示。
nonhermitian_measure_with_grad()是用于计算非厄米哈密顿量 H 在两自定义态矢下的期望值
实例需要传递的实参有四个:两个不同的线路circ1circ2、厄米或非厄米哈密顿量hams、和线路参量parameter_resolvers
.

做一个简单的演示。我们利用旋转Y门和SWAP门组成线路 circ1
,当旋转Y门
的旋转角为 时就变成了一个单比特X门,作用到 qubit[0]
上使其处于激发态 ,再用两比特SWAP门将激发转换到 quibt[1]
上,所以线路 circ1
得到的末态为此时如果 circ2
我们不加任何门操作,则测量左失为初态 为了计算保真度,将被测量的哈密顿量定义为, 则测量得到的期望值变为末态和目标态
之间的保真度.


  • 算子测试
下面将演示如何实现上述过程,首先搭建线路 circ1,并定义线路参数pr
from mindquantum import Circuit
from mindquantum.gate import RY, SWAPGate
from mindquantum import ParameterResolver
import numpy as np
circ1 = Circuit()
circ1 += RY('p1').on(0)
circ1 += SWAPGate().on([0, 1])
theta = np.pi
pr = ParameterResolver({'p1': theta})

搭建一个空线路circ2,不添加任何门操作
circ2 = Circuit()
将目标态矢编码成稀疏矩阵,并转换成 Hamiltonian 格式
from mindquantum import *
from scipy.sparse import csr_matrix
zero_state = np.array([1, 0, 0, 0], dtype=complex)
target_state = np.array([0, 0, 1, 0], dtype=complex)
target_state_ham = zero_state.reshape(4, 1).dot(target_state.reshape(1, 4))
hams = Hamiltonian(csr_matrix(target_state_ham))
生成量子模拟器,将右矢线路 circ1,左失线路circ2
、哈密顿量hams、线路参数 pr一起传入函数 nonhermitian_measure_with_grad()并接受返回值,保真度 f
、梯度g,最后对其打印显示。

from mindquantum.simulator import Simulator
s = Simulator('projectq', 2)
f, g = s.nonhermitian_measure_with_grad(circ2, circ1, [hams], [pr])
print('The fidelity is',f)
print('The gradient is', np.round(g, 8))
The fidelity is [[1.]]
The gradient is [[[0.]]]

可见线路 circ1得到的末态与目标态保真度为1,返回梯度为0。由于我们预先已经将目标态矢转换成非厄米哈密顿量,又将circ1设计成能产生该态的线路,所以保真度为1符合我们的预期。又因为此时保真度处于最大值,其对线路参数的梯度应为0,也符合预期。

下面为了不失一般性,我们对保真度和梯度值进行进一步验证,我们将第一个旋转Y门 的旋转角度加一个小的偏移,这时我们仍以作为目标态,此时得到的结果为
delta = 0.1 * theta
theta = theta + delta
pr = ParameterResolver({'p1': theta})
f, g = s.nonhermitian_measure_with_grad(circ2, circ1, [hams], [pr])
print('The fidelity is',f)
print('The gradient is', np.round(g, 8))
The fidelity is [[0.98768834]]
The gradient is [[[-0.07821723]]]

由于我们引入了偏移角,所以保真度不再为1,梯度不再为零。下面我们对两个值进行验证:首先对量子门 RY
门、SWAP
门进行定义
RY = np.array([[np.cos(theta / 2), -np.sin(theta / 2)], [np.sin(theta / 2), np.cos(theta / 2)]],dtype=np.complex128)
RY_twoqubits = np.kron(np.array([[1,0], [0,1]]), RY)
SWAP = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]], dtype=np.complex128)
并将其作用到初态上得到 right_state,定义目标态 target_state,计算内积得到保真度
right_state = np.dot(np.dot(SWAP, RY_twoqubits), np.array([1,0,0,0]).reshape(4,1))
target_state = np.array([0,0,1,0])
print('The fidelity between the target state and the final state is', np.round(np.dot(target_state, right_state), 8))
The fidelity between the target state and the final state is [0.98768834+0.j]
同样地,定义 门的对参求导的矩阵形式
RY_diff = 0.5 * np.array([[-np.sin(theta / 2), -np.cos(theta / 2)], [np.cos(theta / 2), -np.sin(theta / 2)]], dtype=np.complex128)
RY_diff_twoqubits = np.kron(np.array([[1,0], [0,1]]), RY_diff)
演化得到右矢,计算梯度
right_state_grad = np.array([1,0,0,0]).reshape(4,1)
right_state_grad = np.dot(np.dot(SWAP, RY_diff_twoqubits), right_state_grad)
print('The gradient is', np.round(np.dot(target_state, right_state_grad), 8))


The gradient is [-0.07821723+0.j]

结果相符,验证完毕。


量子相位估计算法演示

  • 概述
量子相位估计算法(Quantum hase Estimation Algorithm,简称QPE),是很多量子算法的关键。假设一个幺正算符 U,这个幺正算符作用在其本征态上会出现一个相位,现在我们假设U算符的本征值未知,也就是 未知,但是U 算符和本征态 已知,相位估计算法的作用就是对这个相位 进行估计。 

  • 算法解析
相位估计算法的实现需要两个寄存器(register),第一个寄存器包含t个初始在的
的量子比特,比特数和最后相位估计的结果的精度和算法的成功概率相关;第二个寄存器初始化在幺正算符 U 的而本征态 上。相位估计算法主要分为三步:
(a). 对第一个寄存器的所有量子比特进行 Hadamard
门操作,对第二寄存器连续进行 控制U
门操作,其中U门的幂次依次为
,控制比特依次为这时第一寄存器中的态就会变为

其中k为直积态的十进制表示,比如 k=0 表示第一寄存器中t个比特全部在基态 ,  k=2 表示,以此类推。
(b). 对第一寄存器的进行量子傅里叶变换的逆变换(Inverse Quantum Fourier Transform),在线路中表示成 , 对 进行逆量子傅里叶变换可得
其中
为本征基矢 对应的概率幅 。由上式可得,当为整数,且满足 时,概率幅取最大值1,此时第一寄存器的末态可以精确反映
不是整数时,的估计,且t越大,估计精度越高。
(c). 对第一寄存器的量子比特进行测量,得到第一寄存器的末态;从中找到最大的概率幅 ,其对应的本征基矢中的在除以 即为相位的估计值。

  • QPE代码实现

下面用一个实例来演示如何在mindquantum实现相位估计算法,选择 T门作为进行估计的幺正算符,由定义
可知需要估计的相位角为
首先导入相关依
from mindquantum import Circuit
from mindquantum import Simulator
from mindquantum import UN, PhaseShift, qft, H, X, BARRIER
import numpy as np




Mindquantum.UN
可以指定量子门,目标比特和控制比特,从而在线路中搭建门操作。因为我们已知,当 时可令为整数,所以第一寄存器只需要3个比特即可准确估计;又已知 T
门的本征态为
,所以第二寄存器选择一个比特,即:我们需要搭建4比特线路,前比特用于估计,属于第一寄存器;属于第二寄存器用于传入 T 算符的本征态。
利用 UN
进行 Hadamard
门操作, 用 X
门对 进行翻转,得到 T
门的本征态
n = 3c = Circuit()c += UN(H, n)c += X.on(n)c
q0: ──H──
q1: ──H──
q2: ──H──
q3: ──X──</pre>

为目标比特,添加 控制PhaseShift
for i in range(n):
c += PhaseShift({'phi': 2**i}).on(n, n-i-1)
对第一寄存器比特进行逆傅里叶变换
c += BARRIER
c += qft(range(n)).hermitian

选择后端、传入总比特数创建模拟器,将 值传入并进行演化,得到末态
sim = Simulator('projectq', c.n_qubits)phi = 0.125sim.apply_circuit(c,{'phi': 2*np.pi*phi})qs = sim.get_qs()

找出概率幅最大值的位置
index = np.argmax(np.abs(qs))print(index)
12

注意此时的 index
对应的
并不是真正的估计值,被除之后也不是,因为测量结果中包括第二寄存器中的辅助比特,需要将index
转成二进制后将辅助位剔除
bit_string = bin(index)[2:].zfill(c.n_qubits)[1:]
print(bit_string)
100
在将二进制转回十进制,得到我们最终的估计值
theta_exp = int(bit_string[::-1], 2) / 2**n


0.125

可见得到的估计相位和 近似相等。


总结


这次暑假点亮计划对我来说是一次很特殊的学习和体验,作为物理系学生,在日常的科研生活中很少有机会接触到开发课题,所以十分缺乏对开发工具和开发流程的基本认识。随着课题的逐渐推进,对基于c++加python的开发流程逐渐熟悉,同时对不同语言的下的计算逻辑也形成了一个清晰的认识。在项目开发期间,个人都处在一个充分学习的状态,非常奇妙。

另外感谢导师徐旭升老师,他在过程中非常耐心,乐观的性格并没有让在线会议显得非常局促,交流过程自然且有效,专业的能力给予我很多帮助。
总而言之,感谢本次开源之夏活动,这个暑假让我受益匪浅。


/

本文为浙江大学闫家舜同学原创文章,欢迎大家投稿分享


END


文章转载自开源之夏,如果涉嫌侵权,请发送邮件至:contact@modb.pro进行举报,并提供相关证据,一经查实,墨天轮将立刻删除相关内容。

评论