对于给定量子纯态，求其中某一或某些量子比特的约化密度矩阵

对于给定量子纯态，求其中某一或某些量子比特的约化密度矩阵

在量子计算的结果中，我们可能只对其中几个量子比特感兴趣，这时就要求它们的约化密度矩阵。本文基于 MindQuantum 来实现这一点（撰写本文时，MindQuantum 的发行版为 0.8.0，尚不支持此功能，虽然后续版本中肯定会支持，但是考虑到网上目前很少有关于这方面内容的介绍，我还是稍微写得详细一些，便于读者参考。如果有地方表达不清，或者有错误，欢迎批评指正。）

MindQuantum 中量子态的表示

为了保证运算效率，MindQuantum 采用 little-endian 表示。我们先来说明一下 little-endian 表示。我们考虑下面一个实例：对第 0 个量子比特施加 X 门，而对第二个量子比特不施加作用。对结果量子态的传统描述方法（比如在 Nielsen 的 Quantum computation and quantum information 中）应为：

$$|q_0⟩\otimes |q_1⟩=|1⟩\otimes |0⟩=|10⟩=[0,0,1,0{]}^T.$$

但在 little-endian 表示下，序号小的量子比特在右边，而序号大的量子比特在左边，这样一来，量子态就成了

$$|q_1⟩\otimes |q_0⟩=|0⟩\otimes |1⟩=|01⟩=[0,1,0,0{]}^T.$$

如下：

import numpy as np
from mindquantum import *

circ = Circuit()
circ += X.on(0)
circ += I.on(1) # 可以直接省略该操作

sim = Simulator('mqvector', n_qubits=2)
sim.apply_circuit(circ)
print('量子态的狄拉克符号表示：')
print(sim.get_qs(True))

print('\n量子态的矢量表示：')
print(sim.get_qs())

量子态的狄拉克符号表示：
1¦01⟩

量子态的矢量表示：
[0.+0.j 1.+0.j 0.+0.j 0.+0.j]

约化密度矩阵：

密度矩阵最深刻的应用也许是作为描述复合量子系统子系统的工具。这种描述由约化密度矩阵提供。约化密度矩阵在分析复合量子系统时非常有用，几乎是不可缺少的。

假设有物理系统 A 和 B，其状态由密度矩阵 ${\rho }^{AB}$ 描述。对于系统 A，约化密度矩阵定义为：

$${\rho }^A\equiv {\mathrm{tr}}_B({\rho }^{AB}).$$

对于一个双量子比特系统的量子纯态 $|\psi ⟩=a_{00}|00⟩+a_{01}|01⟩+a_{10}|10⟩+a_{11}|11⟩$ 的密度矩阵 $\rho$:
$$\rho =|\psi ⟩⟨\psi |=\sum _{i,j,m,n=0}^{1,1,1,1}{a}_{ij}{a}_{mn}^{\ast }|ij⟩⟨mn|.$$
其中，矩阵元 $|ij⟩⟨mn|$ 在矩阵 $\rho$ 中的位置为二进制表示的第 $ij$ 行，第 $mn$ 列。比如矩阵元 $|11⟩⟨10|$ 在矩阵中的位置为第 3 行，第 2 列。（另一种角度认为 $|11⟩⟨10|$ 本身就是一个矩阵，这也无妨，因为不影响结果：矩阵 $|11⟩⟨10|$ 是一个第 3 行，第 2 列元素为 1，其他元素都为 0 的矩阵。）

而量子比特 $q_0$ 的约化密度矩阵
$${\rho }^0={\mathrm{tr}}_1\rho =\sum _{i,j,m,n=0}^{1,1,1,1}{a}_{ij}{a}_{mn}^{\ast }|i⟩⟨m|{\delta }_{j,n}.$$

也就是 ${\rho }^0$ 矩阵元 $|i⟩⟨m|$ 的值为 ${\sum }_{j,n=0}^{1,1}{a}_{ij}{a}_{mn}^{\ast }{\delta }_{j,n}$。接下来我们来代码实现。

import numpy as np
from mindquantum import *

circ = Circuit()
circ += X.on(0)
circ += H.on(1)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1)) # 用列矢量来表示量子态。

rho = psi@psi.T.conj() # 构建密度矩阵

rho_0 = np.zeros((2,2),dtype='complex') # 用于盛装第 0 个量子比特的约化密度矩阵

for j in range(2): # rho_0 的第 j 行
    for n in range(2): # rho_0 的第 n 列
        for im in range(2): # q_1 的状态有两种： 0 和 1
            index_row = int(f'{im}'+f'{j}', base=2) # 我们求量子比特 q_0 的约化密度矩阵，所以要对左边量子比特 q_1 的所有状态求和。
            index_col = int(f'{im}'+f'{n}',base=2)  # 注意两个量子比特的左右位置关系
            rho_0[j,n] += rho[index_row, index_col]
            
print('q_0 的约化密度矩阵为：\n')
print(rho_0)

# 检验一下
circ = Circuit()
circ += X.on(0)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1)) # 用列矢量来表示量子态。
rho_tem = psi@psi.T.conj() # 构建密度矩阵
print('\n正确答案为：\n')
print(rho_tem)

q_0 的约化密度矩阵为：

[[0.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

正确答案为：

[[0.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

同理，可以得到量子比特 $q_1$ 的约化密度矩阵 ${\rho }^1$

rho_1 = np.zeros((2,2),dtype='complex')    

for i in range(2):
    for m in range(2):
        for jn in range(2): # q_0 的状态有两种： 0 和 1
            index_row = int(f'{i}'+ f'{jn}', base=2) # 我们求量子比特 q_1 的约化密度矩阵，所以要对左边量子比特 q_0 的所有状态求和。
            index_col = int(f'{m}'+ f'{jn}', base=2) # 注意两个量子比特的左右位置关系
            rho_1[i,m] += rho[index_row, index_col]
            
print('q_1 的约化密度矩阵为：\n')
print(rho_1)

# 检验一下
circ = Circuit()
circ += H.on(0)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1)) # 用列矢量来表示量子态。
rho_tem = psi@psi.T.conj() # 构建密度矩阵
print('\n正确答案为：\n')
print(rho_tem)

q_1 的约化密度矩阵为：

[[0.5+0.j 0.5+0.j]
 [0.5+0.j 0.5+0.j]]

正确答案为：

[[0.5+0.j 0.5+0.j]
 [0.5+0.j 0.5+0.j]]

接下来求解多量子比特系统中，某一量子比特的约化密度矩阵。

对于多个量子比特系统的密度矩阵：

$$\rho =\sum _{i,...,j,...,k,l,...,m,...,n=0}^{1,...1,...1,1,...1,...1}{a}_{i...j...k}{a}_{l...m...n}^{\ast }|i...j...k⟩⟨l...m....n|.$$

对某个量子比特 $q_j$ 的约化密度矩阵为 ${\rho }^j={\sum }_{i,...,j,...,k,l,...,m,...,n=0}^{1,...1,...1,1,...1,...1}{a}_{i...j...k}{a}_{l...m...n}^{\ast }|i...j...k⟩⟨l...m....n|{\delta }_{(i...j-1,j+1,...k),(l...,m-1,m+1,...n)}.$

对给定的量子态，求解其中某个量子比特的约化密度矩阵。

比如在下面的例子中，我们设感兴趣的量子比特为 2。

import numpy as np
from mindquantum import *

circ = Circuit()
circ += H.on(0)
circ += X.on(1)
circ += Y.on(2)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1))
rho = psi@psi.T.conj()
# print(rho)

rho_i = np.zeros((2,2), dtype='complex')
n_qubits = int(np.log2(len(psi)))

index = 2 # 量子比特的序号

for i in range(2):
    for j in range(2):
        for k in range(2**(n_qubits-1)):
            str_res = '{0:b}'.format(k).rjust(n_qubits-1,'0')
            str_row = str_res[:len(str_res)-index] + '{0:b}'.format(i) + str_res[len(str_res)-index:] # 拼接字符串，注意比特顺序
            index_row = int(str_row, base=2)
            str_col = str_res[:len(str_res)-index] + '{0:b}'.format(j) + str_res[len(str_res)-index:]
            index_col = int(str_col, base=2)
            rho_i[i,j] += rho[index_row, index_col]
#             print(rho[index_row, index_col])
print(f'计算得到的第 {index} 个量子比特的约化密度矩阵为：\n')
print(rho_i)

# 检验一下
circ = Circuit()
circ += Y.on(0)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1)) # 用列矢量来表示量子态。
rho_tem = psi@psi.T.conj() # 构建密度矩阵
print('\n正确答案为：\n')
print(rho_tem)

计算得到的第 2 个量子比特的约化密度矩阵为：

[[0.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

正确答案为：

[[0.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

对给定的量子态，求解其中某些量子比特的约化密度矩阵。

特例，对一个由 3 个量子比特构成的系统中，求由 $q_0$ 和 $q_1$ 构成的子系统的约化密度矩阵

import numpy as np
from mindquantum import *

circ = Circuit()
circ += X.on(0)
circ += H.on(1)
circ += I.on(2)
psi = circ.get_qs()

psi = np.reshape(psi, (8,1))
rho = psi@psi.T.conj()
rho_01 = np.zeros((4,4), dtype='complex')

n_qubits = int(np.log2(len(psi)))

for i in range(4):
    for j in range(4):
        for k in range(2):
            index_row = int(f'{k}'+  '{0:b}'.format(i).rjust(n_qubits-1,'0'), base=2)
            index_col = int(f'{k}'+  '{0:b}'.format(j).rjust(n_qubits-1,'0'), base=2)
            
            rho_01[i,j] += rho[index_row, index_col]
print('q0 和 q1 构成的子系统的约化密度矩阵为：\n')
print(rho_01)

# 检验一下
circ = Circuit()
circ += X.on(0)
circ += H.on(1)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1)) # 用列矢量来表示量子态。
rho_tem = psi@psi.T.conj() # 构建密度矩阵
print('\n正确答案为：\n')
print(rho_tem)

q0 和 q1 构成的子系统的约化密度矩阵为：

[[0. +0.j 0. +0.j 0. +0.j 0. +0.j]
 [0. +0.j 0.5+0.j 0. +0.j 0.5+0.j]
 [0. +0.j 0. +0.j 0. +0.j 0. +0.j]
 [0. +0.j 0.5+0.j 0. +0.j 0.5+0.j]]

正确答案为：

[[0. +0.j 0. +0.j 0. +0.j 0. +0.j]
 [0. +0.j 0.5+0.j 0. +0.j 0.5+0.j]
 [0. +0.j 0. +0.j 0. +0.j 0. +0.j]
 [0. +0.j 0.5+0.j 0. +0.j 0.5+0.j]]

普适，比如在下面的例子中，我们设感兴趣的量子比特为 [1,2]。

import numpy as np
from mindquantum import *
import copy

circ = Circuit()
circ += H.on(0)
circ += X.on(1)
circ += Y.on(2)
circ += Z.on(3)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1))
rho = psi@psi.T.conj()
# print(rho)

index_list = [1,2] # 量子比特的序号 注意大小顺序，不然后面构造字符串时可能会出错
index_list.sort()
rho_s = np.zeros((2**len(index_list), 2**len(index_list)), dtype='complex')
n_qubits = int(np.log2(len(psi)))

for i in range(2**len(index_list)):
    for j in range(2**len(index_list)):
        for k in range(2**(n_qubits-len(index_list))):
            str_res = '{0:b}'.format(k).rjust(n_qubits-len(index_list),'0')
            str_row = copy.deepcopy(str_res) # 注意要 deepcopy，不然 str_res 会动态变化
            str_col = copy.deepcopy(str_res)
            for i_index, index in enumerate(index_list):
                str_row = str_row[:len(str_row)-index] + '{0:b}'.format(i).rjust(len(index_list),'0')[-(i_index+1)] + str_row[len(str_row)-index:]
                str_col = str_col[:len(str_col)-index] + '{0:b}'.format(j).rjust(len(index_list),'0')[-(i_index+1)] + str_col[len(str_col)-index:]
            index_row = int(str_row, base=2)
            index_col = int(str_col, base=2)
            rho_s[i,j] += rho[index_row, index_col]
            
print(f'计算得到的第 {index_list} 个量子比特的约化密度矩阵为：\n')
print(rho_s)

# 检验一下
circ = Circuit()
circ += X.on(0)
circ += Y.on(1)
psi = circ.get_qs()

psi = np.reshape(psi, (len(psi),1)) # 用列矢量来表示量子态。
rho_tem = psi@psi.T.conj() # 构建密度矩阵
print('\n正确答案为：\n')
print(rho_tem)

计算得到的第 [1, 2] 个量子比特的约化密度矩阵为：

[[0.+0.j 0.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 0.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 0.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 0.+0.j 0.+0.j 1.+0.j]]

正确答案为：

[[0.+0.j 0.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 0.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 0.+0.j 0.+0.j 0.+0.j]
 [0.+0.j 0.+0.j 0.+0.j 1.+0.j]]

当然，上面的方法比较笨，当量子比特数很大时，构造整体密度矩阵会有压力。更好的方案是直接对量子态矢量进行对应操作，这项优化任务就等后面有时间了再来完成吧。