基于 MindQuantum 0.9 实现 expm 门

基于 MindQuantum 0.9 实现 expm 门

在基于变分算法实现量子优化控制时，我们往往需要执行 expm 门。即，设计一个量子门实现对于一个给定的力学量 $H$ 执行操作 $U(t)=\text{expm}(-iHt)$，其中 $H$ 具有厄米性。但 mindquantum 0.9 （以下简称 mq09）还没有 API 来直接实现此功能。本文中我们将介绍如何基于 mq09 的自定义单参数门来实现该功能。

由算子函数的定义，设
$$A=\sum _{a}a|a⟩⟨a|$$
是正规矩阵 A（厄米矩阵属于正规矩阵）的一个谱分解。则有
$$f(A)\equiv \sum _{a}f(a)|a⟩⟨a|$$

如正算子的平方根、正定算子的的对数或正规算子的指数 expm 等都可以如此定义。

以 $\text{expm}(-i{\sigma }_{x}t/2)$ 为例（也就是 mq09 中的 RX 门）。根据上述定义，则可写为
$$\text{expm}(-i{\sigma }_{x}t/2)={\text{e}}^{-i{\sigma }_{x}t/2}\equiv \sum _i{\text{e}}^{-i{E}_{i}t/2}|E_i⟩⟨E_i|$$

其中 $E_i$ 和 $|E_i⟩$ 分别为 ${\sigma }_x$ 第 $i$ 个本征值和对应的本征态。

考虑到，$E_0=1$，$E_1=-1$，$|E_0⟩=[0.5,0.5{]}^T$，$|E_1⟩=[0.5,-0.5{]}^T$

$$\text{expm}(-i{\sigma }_{x}t)={\text{e}}^{-it/2}\left(\begin{array}{cc}0.5& 0.5\\ 0.5& 0.5\end{array}\right)-{\text{e}}^{it/2}\left(\begin{array}{cc}0.5& -0.5\\ -0.5& 0.5\end{array}\right)$$

我们也可以使用 numpy 来计算各本征值，以及本征向量组成的矩阵。

import numpy as np

sx = np.array([[0, 1], [1, 0j]])
alpha = 1
eig_vs, eig_vecs = np.linalg.eig(sx)
expm = 0
for i in range(len(eig_vs)):
    eig_mat = eig_vecs[:,i].reshape(-1, 1).conj() @ eig_vecs[:,i].reshape(1, -1)
    print(f'eig_v{i}\n', eig_vs[i])
    print(f'eig_mat{i}\n', eig_mat, '\n')
    expm = expm + np.exp(-1j * eig_vs[i] * alpha / 2) * eig_mat
print('expm:\n', expm)

eig_v0
 (0.9999999999999996+0j)
eig_mat0
 [[0.5+0.j 0.5+0.j]
 [0.5+0.j 0.5+0.j]] 

eig_v1
 (-0.9999999999999999+0j)
eig_mat1
 [[ 0.5+0.j -0.5+0.j]
 [-0.5+0.j  0.5+0.j]] 

expm:
 [[8.77582562e-01+1.38777878e-16j 5.55111512e-17-4.79425539e-01j]
 [5.55111512e-17-4.79425539e-01j 8.77582562e-01-2.77555756e-17j]]

我们检查该计算结果是否与 mq09 的结果一致。

from mindquantum import *
np.allclose(expm, RX(1).on(0).matrix())

True

接下来，我们根据以上方法，基于 mq09 的自定义单参数门，实现 expm 门的构造。

sx_p = np.array([[1, 1], [1, 1]]) / 2
sx_n = np.array([[1, -1], [-1, 1]]) / 2

def matrix(alpha):
    return np.exp(-1j*alpha/2) * sx_p + np.exp(1j*alpha/2) * sx_n

def diff_matrix(alpha):
    return -1j * sx @ (np.exp(-1j*alpha/2) * sx_p + np.exp(1j*alpha/2) * sx_n)
    
Expm = gene_univ_parameterized_gate('Expm', matrix, diff_matrix)

接下来就可以使用该量子门构建量子线路，以及执行变分量子算法了。

circ = Circuit()
circ += Expm('a').on(0)
print(circ)

q0: ──Expm(a)──

import mindspore as ms
from mindspore.nn import Adam, TrainOneStepCell
from mindspore.common.initializer import initializer
from mindspore.common.parameter import Parameter

sim = Simulator('mqvector', 1)
ham = Hamiltonian(QubitOperator('Z0'))
grad_ops = sim.get_expectation_with_grad(ham, circ)

ms.set_context(mode=ms.PYNATIVE_MODE, device_target="CPU")

qnet = MQAnsatzOnlyLayer(grad_ops)

opti = Adam(qnet.trainable_params(), learning_rate=0.5) 
net = TrainOneStepCell(qnet, opti)

for i in range(200):
    res = net()
    if i % 20 == 0:
        print(i, ': ', res)
        
print(qnet.weight.asnumpy())

0 :  [0.99999225]
20 :  [-0.89419425]
40 :  [-0.98503464]
60 :  [-0.99691594]
80 :  [-0.9999923]
100 :  [-0.99994886]
120 :  [-0.99999577]
140 :  [-0.9999998]
160 :  [-1.]
180 :  [-1.]
[3.1415365]

以上得到的结果与 RX 门的结果完全一致。其他 expm 门也可以照此方法构造。

注意事项：由于编译问题，matirx() 和 diff_matrix() 中不能使用 for 循环。即 不能 构造为下面的代码。

def matrix(alpha):
    eig_vs, eig_vecs = np.linalg.eig(sx)
    expm = 0
    for i in range(len(eig_vs)):
        eig_mat = eig_vecs[:,i].reshape(-1, 1) @ eig_vecs[:,i].reshape(1, -1)
        expm = expm + np.exp(-1j * alpha * eig_vs[i] / 2) * eig_mat
    return expm

def diff_matrix(alpha):
    eig_vs, eig_vecs = np.linalg.eig(sx)
    expm = 0
    for i in range(len(eig_vs)):
        eig_mat = eig_vecs[:,i].reshape(-1, 1) @ eig_vecs[:,i].reshape(1, -1)
        expm = expm + np.exp(-1j * alpha * eig_vs[i] / 2) * eig_mat
    return -1j * sx @ expm