对于给定量子纯态,求其中某一或某些量子比特的约化密度矩阵
在量子计算的结果中,我们可能只对其中几个量子比特感兴趣,这时就要求它们的约化密度矩阵。本文基于 MindQuantum 来实现这一点(撰写本文时,MindQuantum 的发行版为 0.8.0,尚不支持此功能,虽然后续版本中肯定会支持,但是考虑到网上目前很少有关于这方面内容的介绍,我还是稍微写得详细一些,便于读者参考。如果有地方表达不清,或者有错误,欢迎批评指正。)
MindQuantum 中量子态的表示
为了保证运算效率,MindQuantum 采用 little-endian 表示。我们先来说明一下 little-endian 表示。我们考虑下面一个实例:对第 0 个量子比特施加 X 门,而对第二个量子比特不施加作用。对结果量子态的传统描述方法(比如在 Nielsen 的 Quantum computation and quantum information 中)应为:
∣ q 0 ⟩ ⊗ ∣ q 1 ⟩ = ∣ 1 ⟩ ⊗ ∣ 0 ⟩ = ∣ 10 ⟩ = [ 0 , 0 , 1 , 0 ] T . |q_0\rangle \otimes |q_1\rangle = |1\rangle\otimes|0\rangle = |10\rangle = [ 0, 0, 1, 0]^T. ∣q0⟩⊗∣q1⟩=∣1⟩⊗∣0⟩=∣10⟩=[0,0,1,0]T.
但在 little-endian 表示下,序号小的量子比特在右边,而序号大的量子比特在左边,这样一来,量子态就成了
∣ q 1 ⟩ ⊗ ∣ q 0 ⟩ = ∣ 0 ⟩ ⊗ ∣ 1 ⟩ = ∣ 01 ⟩ = [ 0 , 1 , 0 , 0 ] T . |q_1\rangle \otimes |q_0\rangle = |0\rangle\otimes|1\rangle=|01\rangle= [0, 1, 0, 0]^T. ∣q1⟩⊗∣q0⟩=∣0⟩⊗∣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,其状态由密度矩阵 ρ A B \rho^{AB} ρAB 描述。对于系统 A,约化密度矩阵定义为:
ρ A ≡ t r B ( ρ A B ) . \rho^A\equiv\mathrm{tr}_B(\rho^{AB}). ρA≡trB(ρAB).
对于一个双量子比特系统的量子纯态
∣
ψ
⟩
=
a
00
∣
00
⟩
+
a
01
∣
01
⟩
+
a
10
∣
10
⟩
+
a
11
∣
11
⟩
|\psi\rangle=a_{00}|00\rangle+a_{01}|01\rangle+a_{10}|10\rangle+a_{11}|11\rangle
∣ψ⟩=a00∣00⟩+a01∣01⟩+a10∣10⟩+a11∣11⟩ 的密度矩阵
ρ
\rho
ρ:
ρ
=
∣
ψ
⟩
⟨
ψ
∣
=
∑
i
,
j
,
m
,
n
=
0
1
,
1
,
1
,
1
a
i
j
a
m
n
∗
∣
i
j
⟩
⟨
m
n
∣
.
\rho = |\psi\rangle\langle\psi| = \sum_{i,j,m,n=0}^{1,1,1,1} a_{ij}a^{*}_{mn}|ij\rangle\langle mn|.
ρ=∣ψ⟩⟨ψ∣=i,j,m,n=0∑1,1,1,1aijamn∗∣ij⟩⟨mn∣.
其中,矩阵元
∣
i
j
⟩
⟨
m
n
∣
|ij\rangle\langle mn|
∣ij⟩⟨mn∣ 在矩阵
ρ
\rho
ρ 中的位置为二进制表示的第
i
j
ij
ij 行,第
m
n
mn
mn 列。比如矩阵元
∣
11
⟩
⟨
10
∣
|11\rangle\langle 10|
∣11⟩⟨10∣ 在矩阵中的位置为第 3 行,第 2 列。(另一种角度认为
∣
11
⟩
⟨
10
∣
|11\rangle\langle 10|
∣11⟩⟨10∣ 本身就是一个矩阵,这也无妨,因为不影响结果:矩阵
∣
11
⟩
⟨
10
∣
|11\rangle\langle 10|
∣11⟩⟨10∣ 是一个第 3 行,第 2 列元素为 1,其他元素都为 0 的矩阵。)
而量子比特
q
0
q_0
q0 的约化密度矩阵
ρ
0
=
t
r
1
ρ
=
∑
i
,
j
,
m
,
n
=
0
1
,
1
,
1
,
1
a
i
j
a
m
n
∗
∣
i
⟩
⟨
m
∣
δ
j
,
n
.
\rho^0=\mathrm{tr}_1\rho=\sum_{i,j,m,n=0}^{1,1,1,1} a_{ij}a^{*}_{mn}|i\rangle\langle m|\delta_{j,n}.
ρ0=tr1ρ=i,j,m,n=0∑1,1,1,1aijamn∗∣i⟩⟨m∣δj,n.
也就是 ρ 0 \rho^0 ρ0 矩阵元 ∣ i ⟩ ⟨ m ∣ |i\rangle\langle m| ∣i⟩⟨m∣ 的值为 ∑ j , n = 0 1 , 1 a i j a m n ∗ δ j , n \sum_{j,n=0}^{1,1} a_{ij}a^{*}_{mn}\delta_{j,n} ∑j,n=01,1aijamn∗δ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 q_1 q1 的约化密度矩阵 ρ 1 \rho^1 ρ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]]
接下来求解多量子比特系统中,某一量子比特的约化密度矩阵。
对于多个量子比特系统的密度矩阵:
ρ = ∑ i , . . . , j , . . . , k , l , . . . , m , . . . , n = 0 1 , . . . 1 , . . . 1 , 1 , . . . 1 , . . . 1 a i . . . j . . . k a l . . . m . . . n ∗ ∣ i . . . j . . . k ⟩ ⟨ l . . . m . . . . n ∣ . \rho=\sum_{i,...,j,...,k,l,...,m,...,n=0}^{1,...1,...1,1,...1,...1}a_{i...j...k}a^*_{l...m...n}|i...j...k\rangle\langle l...m....n|. ρ=i,...,j,...,k,l,...,m,...,n=0∑1,...1,...1,1,...1,...1ai...j...kal...m...n∗∣i...j...k⟩⟨l...m....n∣.
对某个量子比特 q j q_j qj 的约化密度矩阵为 ρ j = ∑ i , . . . , j , . . . , k , l , . . . , m , . . . , n = 0 1 , . . . 1 , . . . 1 , 1 , . . . 1 , . . . 1 a i . . . j . . . k a l . . . m . . . n ∗ ∣ i . . . j . . . k ⟩ ⟨ l . . . m . . . . n ∣ δ ( i . . . j − 1 , j + 1 , . . . k ) , ( l . . . , m − 1 , m + 1 , . . . n ) . \rho^j = \sum_{i,...,j,...,k,l,...,m,...,n=0}^{1,...1,...1,1,...1,...1}a_{i...j...k}a^*_{l...m...n}|i...j...k\rangle\langle l...m....n|\delta_{(i...j-1,j+1,...k),(l...,m-1,m+1,...n)}. ρj=∑i,...,j,...,k,l,...,m,...,n=01,...1,...1,1,...1,...1ai...j...kal...m...n∗∣i...j...k⟩⟨l...m....n∣δ(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_0 q0 和 q 1 q_1 q1 构成的子系统的约化密度矩阵
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]]
当然,上面的方法比较笨,当量子比特数很大时,构造整体密度矩阵会有压力。更好的方案是直接对量子态矢量进行对应操作,这项优化任务就等后面有时间了再来完成吧。