基于 MindQuantum 实现振幅编码

MindQuantum 当前版本 0.8.0 已拥有振幅编码 API：amplitude_encoder，不过其只能处理不含相位信息的输入数据。比如 [0.5, -0.5, 0.5, 0.5] 这样的输入数据。实践中，我们有时需要实现 $\mathbf{x}=(|x_0|{\mathrm{e}}^{i{\omega }_0},\cdots ,|x_{N-1}|{\mathrm{e}}^{i{\omega }_{N-1}})\to {\sum }_p|x_p|{\mathrm{e}}^{i{\omega }_p}|p⟩$。也即输入的数据是带有相位的。

本文采用“从上到下”振幅编码策略。
参考文献：Configurable sublinear circuits for quantum state preparation, Quantum Information Processing (2023) 22:123. https://doi.org/10.1007/s11128-023-03869-7

先来看看 MindQuantum 的方案。编码线路及参数为：

from mindquantum import *

vec_in = [0.5, -0.5, 0.5, 0.5] # 输入数据矢量
encoder, parameterResolver = amplitude_encoder(vec_in, 2) # 会自动对输入矢量进行归一化、补零或截取操作
sim = Simulator('mqvector', 2)
sim.apply_circuit(encoder, parameterResolver)

print('\n编码后的量子态为：\n')
print( sim.get_qs(True))
print('\n旋转角度参数为：\n\n', parameterResolver)
print('\n编码线路为：\n')
encoder.svg()

编码后的量子态为：

1/2¦00⟩
-1/2¦01⟩
1/2¦10⟩
1/2¦11⟩

旋转角度参数为：

 {'alpha0': 1.5708, 'alpha1': -1.5708, 'alpha2': 1.5708}, const: 0

编码线路为：

接下来，我们来完成可实现带有相位信息的振幅编码方案。我们先介绍一下状态树和角度树的概念。

状态树和角度树分别用来统计振幅相位信息和分配旋转角度。树的每一个节点都代表了一个受控门。树的高度代表了所用比特数。从每个点分出来的两条边表明了每一个受控门操作将希尔伯特空间分为两个子空间。因此，经过 $n$ 层操作后，会有 $2^n$ 个不同概率振幅的子空间。

为解释流程，需要根据目标矢量 $\mathbf{x}$ 定义四个参数：

$${\Omega }_{i,k}=\sum _{l=0}^{2^k-1}{\omega }_{i\cdot 2^k+l}/2^{k-1}相位统计\text{\hspace{1em}(1)}$$
$${\eta }_{i,k}=\sqrt{\sum _{l=0}^{2^k-1}|x_{i\cdot 2^k+l}{|}^2}振幅统计\text{\hspace{1em}(2)}$$
$${\lambda }_{j,v}={\Omega }_{2j,v-1}-{\Omega }_{j,v}z轴旋转角度，编码相位\text{\hspace{1em}(3)}$$
$${\beta }_{j,v}={\eta }_{2j,v-1}/{\eta }_{j,v}{\alpha }_{j,v}=2\arcsin ({\beta }_{j,v})y轴旋转角度，编码振幅\text{\hspace{1em}(4)}$$

其中，$j=0,1,2,...,2^{n-v}-1$，$v=1,2,...,n$，并且 $n={\log }_2(N)$。
这些参数用来构造树表示。指标 $k,v$ 是层指标，而 $i,j$ 是给定层下的节点指标。状态树的节点值标志着待编码量子态的振幅和用于构建编码线路的旋转角度。

状态树：统计振幅和相位。节点值：${\eta }_{点,层}{\Omega }_{点,层}$。 ${\eta }_{i,k}$ 振幅， ${\Omega }_{i,k}$ 相位。计算方向为自下而上。

对于输入矢量：$\mathbf{x}=(|x_0|{\mathrm{e}}^{i{\omega }_0},\cdots ,|x_{N-1}|{\mathrm{e}}^{i{\omega }_{N-1}})$。

状态树的最底层对应于 $k=0$，此时，参数点就作用到输入矢量 $\mathbf{x}$ 上了，${\eta }_{i,0}=|x_i|$ （${\eta }_{0,0}=|x_0|$，${\eta }_{1,0}=|x_1|$）和 ${\Omega }_{i,0}=2{\omega }_i$ （${\Omega }_{0,0}=2{\omega }_0$, ${\Omega }_{1,0}=2{\omega }_1$）。向上的每一层的节点中的 $\eta$ 是该节点下辖所有 $|x_i|$ 的振幅和(模方和的开方)（统计下辖节点的振幅）。如 ${\eta }_{0,1}=\sqrt{|x_0{|}^2+|x_1{|}^2}$ ， ${\eta }_{1,1}=\sqrt{|x_2{|}^2+|x_3{|}^2}$ 。

而每个节点中相位 $\Omega$ 则统计该节点下辖所有底层项相位信息 ${\Omega }_{点,层}=下辖所有{\omega }_i的和/2^{层-1}$

角度树：分配振幅和相位。公式 (3) 和 (4) 决定了角度树的旋转角度值。节点值：${\alpha }_{点,层}{\lambda }_{点,层}$。$\alpha$ 绕 $y$ 轴旋转角度，$\lambda$ 绕 $z$ 轴旋转角度。计算方向为自上而下。

$${\alpha }_{点,层}=2\arcsin (\frac{右下位置的{\eta }_{{点}^{\prime },层-1}}{状态树同位置的{\eta }_{点,层}})$$

$${\lambda }_{点,层}=右下位置的{\Omega }_{{点}^{\prime },层-1}-状态树同位置节点的{\Omega }_{点,层}$$

矢量 $|x⟩$ 的相位参量通过绕 $z$ 轴的旋转 ${\lambda }_{j,v}$ 来实现，而大小则通过绕 $y$ 轴旋转角度 ${\alpha }_{j,v}=2\arcsin ({\beta }_{j,v})$ 来实现。

然后介绍“自上而下法”的工作流程。

该方案的量子线路图在接近文尾处有显示 (以两比特情况为例)：

从角度树的 root 开始，先生成
$$|{\psi }_n⟩={\text{e}}^{-i\frac{{\lambda }_{1,n}}{2}}\sqrt{1-|{\beta }_{1,n}{|}^2}|0⟩+{\text{e}}^{i\frac{{\lambda }_{1,n}}{2}}{\beta }_{1,n}|1⟩$$
再生成：
$$|{\psi }_v⟩=\sum _{j=1}^{2^{n-v}}|j-1⟩⟨j-1|{\psi }_{v+1}⟩({\text{e}}^{-i\frac{{\lambda }_{j,v}}{2}}\sqrt{1-|{\beta }_{j,v}{|}^2}|0⟩+{\text{e}}^{i\frac{{\lambda }_{j,v}}{2}}{\beta }_{j,v}|0⟩)$$
如此逐层下分，$v=(n-1),...,1$，就可以得到目标态
$$|{\psi }_1⟩=|x_0|{\text{e}}^{i{\omega }_0}|0⟩+..+|x_{N-1}|{\text{e}}^{i{\omega }_{N-1}}|N-1⟩$$

接下来，我们实基于上述参考文献中的方案现振幅编码。我们具有相位信息的待编码数据为 [0.5, -1+1j, 0.5, 0.5]。 先对输入矢量进行归一化。

vec_in = [0.5, -1+1j, 0.5, 0.5]
vec = normalize(vec_in)

print('归一化后的矢量为：\n\n', vec)

归一化后的矢量为：

 [ 0.30151134+0.j         -0.60302269+0.60302269j  0.30151134+0.j
  0.30151134+0.j        ]

计算树的层数 $n$（等于编码线路所用比特数 n_qubits）。

import numpy as np

n = int(np.log2(len(vec)))

print('树的层数为：\n\n', n)

树的层数为：

 2

先计算状态树中的 ${\eta }_{点,层}$ 部分。注意，在参考文献中，${\eta }_{点,层}$、${\Omega }_{点,层}$、${\alpha }_{点,层}$、${\lambda }_{点,层}$ 中点的取值范围为 [1,2,3…]，而本文中为 [0,1,2,…]。

# 先将 vec 中的各元素 vec_i 标记为 eta_i0（eta_点层）
for i in range(len(vec)):
    exec(f'eta_{i}0=abs(vec[i])')
    
print('0 层诸 eta_i0 的值为:\n')
for i in range(len(vec)):
    print(f'eta_{i}0 =',eval(f'eta_{i}0'))

# 计算其他层的 eta_iv （eta_点层） 值
print('\n其它层的 eta_点层 的值为：\n')

for v in range(1,n+1): # v 层，从小到大，标志方向为从下到上
    for i in range(int(2**(n-v))):
        exec(f'eta_{i}{v}=np.sqrt((eta_{int(2*i)}{v-1})**2 + (eta_{int(2*i+1)}{v-1})**2)')
        print(f'eta_{i}{v} =', eval(f'eta_{i}{v}'))

0 层诸 eta_i0 的值为:

eta_00 = 0.30151134457776363
eta_10 = 0.8528028654224418
eta_20 = 0.30151134457776363
eta_30 = 0.30151134457776363

其它层的 eta_点层 的值为：

eta_01 = 0.9045340337332909
eta_11 = 0.4264014327112209
eta_02 = 1.0

计算状态树的 ${\Omega }_{点,层}$ 部分。

# 先根据 vec 中的各元素 vec_i 计算得到 omega_i(每个数据自带的相位)，该数据用于后面计算 omega_点层
    
for i in range(len(vec)):
    exec(f'omega_{i}=np.angle(vec[i])')
    
print('0 层诸 omega_i 的值为:\n')
for i in range(len(vec)):
    print(f'omega_{i} =',eval(f'omega_{i}'))
    
# 再根据 vec 中的各元素 vec_i 计算得到 omega_i0(omega_点层)，该数据用于后面计算绕 z 的旋转角度 lambda_点层
for i in range(len(vec)):
    exec(f'omega_{i}0=2*np.angle(vec[i])')
    
print('\n0 层诸 omega_i0 的值为:\n')
for i in range(len(vec)):
    print(f'omega_{i}0 =',eval(f'omega_{i}0'))
    
# 计算其他层的 omega_iv （omega_点层） 值
print('\n其它层的 omega_点层 的值为：\n')
n = 2
for v in range(1,n+1): # v 层，从小到大，标志方向为从下到上
    for i in range(int(2**(n-v))):
        exec(f'omega_{i}{v}=0')
        for j in range(int(2**v)):
            exec(f'omega_{i}{v}+=omega_{int(i*2**v+j)}/2**{v-1}')
        print(f'omega_{i}{v} =', eval(f'omega_{i}{v}'))

0 层诸 omega_i 的值为:

omega_0 = 0.0
omega_1 = 2.356194490192345
omega_2 = 0.0
omega_3 = 0.0

0 层诸 omega_i0 的值为:

omega_00 = 0.0
omega_10 = 4.71238898038469
omega_20 = 0.0
omega_30 = 0.0

其它层的 omega_点层 的值为：

omega_01 = 2.356194490192345
omega_11 = 0.0
omega_02 = 1.1780972450961724

计算角度树中 ${\alpha }_{点,层}$ 部分。

alphas = {} # 用列表盛装，方便传入 mindquantum 接口

print('角度树诸节点的数值为：\n')
for v in range(n, 0, -1): # v 层，从大到小，标志方向从上到下
    for i in range(n-v + 1): # i 点
        exec(f'alpha_{i}{v}=2*np.arcsin(eta_{int(2*i+1)}{v-1}/eta_{i}{v})')
        alphas[f'alpha_{i}{v}'] = eval(f'alpha_{i}{v}')
        print(f'alpha_{i}{v} =', eval(f'alpha_{i}{v}'))
        
print('\n汇总：\n')
print(alphas)

角度树诸节点的数值为：

alpha_02 = 0.881021326009397
alpha_01 = 2.4619188346815495
alpha_11 = 1.5707963267948963

汇总：

{'alpha_02': 0.881021326009397, 'alpha_01': 2.4619188346815495, 'alpha_11': 1.5707963267948963}

根据状态树中 ${\Omega }_{点,层}$ 计算角度树中绕 $z$ 轴的旋转角度 ${\lambda }_{点,层}$。

lambdas = {} # 用列表盛装，方便传入 mindquantum 接口

print('角度树诸节点的 lambda_点层 数值为：\n')
for v in range(n, 0, -1): # v 层，从大到小，标志方向从上到下
    for i in range(n-v + 1): # i 点
        exec(f'lambda_{i}{v}=omega_{2*i+1}{v-1} - omega_{i}{v}')
        lambdas[f'lambda_{i}{v}'] = eval(f'lambda_{i}{v}')
        print(f'lambda_{i}{v} =', eval(f'lambda_{i}{v}'))
        
print('\n汇总：\n')
print(lambdas)

角度树诸节点的 lambda_点层 数值为：

lambda_02 = -1.1780972450961724
lambda_01 = 2.356194490192345
lambda_11 = 0.0

汇总：

{'lambda_02': -1.1780972450961724, 'lambda_01': 2.356194490192345, 'lambda_11': 0.0}

搭建振幅编码线路。

def amp_encoder(n): # 角度树层数 n，其和所用比特数是相同的
    _circ = Circuit()
    _circ += RY(f'alpha_0{n}').on(n-1)
    _circ += RZ(f'lambda_0{n}').on(n-1)
    _circ += BarrierGate()
    for i in range(1,n): # 用于确定控制位的字符串位数
        for j in range(int(2**i)):
            string = bin(j)[2:].zfill(i) # 根据字符串的比特来确定是实控还是虚控
            for tem_qubit, bit in enumerate(string):
                qubit = int(n - tem_qubit)
                if bit == '0': # 如果比特为 0，就是虚空，要加 X 门
                    _circ += X.on(qubit-1)
                    
            _circ += RY(f'alpha_{j}{int(n-i)}').on(int(n-1-i), list(range(n-i,n))) # 增加受控 RY 门
            _circ += RZ(f'lambda_{j}{int(n-i)}').on(int(n-1-i), list(range(n-i,n))) # 增加受控 RZ 门
            string = bin(j)[2:].zfill(i)
            
            # 虽然显着重复，但这一句必须要有，因为上个循环汇总，string 已经被 enumerate 完了。
            for tem_qubit, bit in enumerate(string): # 补上虚控位的 X 门
                qubit = int(n - tem_qubit)
                if bit == '0':
                    _circ += X.on(qubit-1)
            _circ += BarrierGate()
    return _circ

circ = amp_encoder(2)
circ.svg()

将计算得到的角度输入到编码线路中，并得到计算结果。

sim = Simulator('mqvector',2)
params = {**alphas, **lambdas}
sim.apply_circuit(circ, pr=params)
print(sim.get_qs(True))

(0.2506975207808922-0.16751072796512168j)¦00⟩
(-0.16637358563154098+0.8364164974920278j)¦01⟩
(0.2506975207808923-0.16751072796512173j)¦10⟩
(0.25069752078089225-0.1675107279651217j)¦11⟩

由于有整体相位，不容易检查结果，我们下面就去掉整体相位。

state0 = sim.get_qs()
state1 = state0 * np.conj(state0[0])
state = state1/np.linalg.norm(state1, ord=2)
print(state)

[ 0.30151134+0.00000000e+00j -0.60302269+6.03022689e-01j
  0.30151134+2.30137075e-17j  0.30151134+2.30137075e-17j]

而回想我们之前的输入数据经归一化之后为

print(vec)

[ 0.30151134+0.j         -0.60302269+0.60302269j  0.30151134+0.j
  0.30151134+0.j        ]

可见，得到的结果和所期望的相一致。