用 Python 写个玩具级量子计算机 (1)
华中师范大学 hahakity
最近量子计算炒的很火，作为物理专业的学生，学了量子力学却不懂量子计算太可惜。这一节使用 Python 编程，构造一个玩具版本的量子计算机，辅助大家对量子力学和量子计算的理解，从而快速使用量子计算做些有意思的工作。
学习目标
- 单量子位 (Qubit)，多Qubit，叠加态 (superposition state), 纠缠态（entanglement)
- 几个常见的量子门 (Gate) 及其作用
- 单Qubit操作，多Qubit操作，态的直积，矩阵的直积
预备知识
- 线性代数，量子力学
极简介绍
量子计算像一本武林秘籍，这本秘籍的心法是量子力学的多世界诠释。
量子位 Qubit
传统计算机使用比特 0 和 1 以及它们之间的逻辑操作 and，or，not，xor 等位运算构造加减乘除，并搭建整个经典计算机体系。
在量子计算机上，基本的存储单元是量子位，即 Qubit。只要搞明白了 Qubit，以及它们之间的逻辑操作，原则上就可以搭建量子计算机。所以本虚拟机仅实现这些基本功能。
Qubit 可以是任意两个本征态 和 ，或者是它们的线性叠加 ，其中 是复数。
用矩阵表示，一个 Qubit 有两个基向量，
,
最简单可对应为自旋向上 与向下 。与经典比特不同的是，一个量子比特可以是自旋 与 的叠加态，在一个世界里粒子自旋向上，另一个世界里自旋向下，就像薛定谔的猫态（既死又活态），
当然，波函数要满足幺正条件: 。
更多的比特可以通过直积 生成，这个直积完成如下操作，
比如两量子比特，
两个 Qubit 生成了4个基向量，每个基向量只有一位为1，其余位为0，称为 one-hot 表示。
一个一般的两 Qubit 态可以展开为，
我们先用 Python 构造单个 Qubit，然后通过直积构造 n 个 Qubit 的 维基向量。
构造基向量
在 jupyter notebook 开头，导入库 numpy，scipy 以及输出数学公式的 Markdown库。
下面这段程序还定义了自旋向上，自旋向下态以及它们的集合 bit。
import numpy as npfrom scipy.linalg import kron
from IPython.display import Markdown as md
spin_up = np.array([[1, 0]]).T
spin_down = np.array([[0, 1]]).T
# bit[0] = |0>, bit[1] = |1>
bit = [spin_up, spin_down]
单个 Qubit 没什么可讲，我们来看一看两个或多个 Qubit 的基向量如何构造
下面定义一个函数，给定 Qubits 序列，返回它的 one hot 表示，
def basis(string='00010'):'''string: the qubits sequence'''
res = np.array([[1]])
# 从最后一位开始往前数，做直积
for idx in string[::-1]:
res = kron(bit[int(idx)], res)
return np.matrix(res)
构造基向量 ，输入 basis('10').A1, 返回,
array([0, 0, 1, 0])
注意：直积操作 用 scipy.linalg 中的 kron() 函数。
basis() 返回了 np.matrix 数据类型，它的 .A1 属性将矩阵转化为 1 维 numpy 数组。
两个Qubit，基向量是4维，3个Qubit是8维，10个Qubit基向量是1024维，在多体量子计算中，大家一般要计算超大矩阵的基态本征值，如果用 100 个 Qubit 进行 个态的多体量子系统计算，可远超经典计算机能够存储和处理的数据极限。自从费曼提出这个想法，大家就一直为量子计算机而疯魔。
这里我们随便试一下 20 个量子比特的基向量的维数，
len(basis('00100000001000011100'))输出为：1048576
即原则上 20 个 Qubit 可以解 1048576 x 1048576 维哈密顿矩阵的本征值问题。
如果是50 个 Qubit，那么单方向维数就是 1125899906842624 。
我们可以看一看 n 个 Qubit 支撑起来的希尔伯特空间，
def hilbert_space(nbit=5):nspace = 2**nbit
for i in range(nspace):
#bin(7) = 0b100
binary = bin(i)[2:]
nzeros = nbit - len(binary)
yield '0'*nzeros + binary
这个函数使用了 python 内置 bin() 函数，将数字转化为二进制序列。希尔伯特空间维数 随 Qubit 个数 n 呈指数增长，所以我们使用 yield 返回一个 generator，不保存数组，但可循环取值。举例，对于 3 个Qubit，希尔伯特空间基向量有8个，
for mi in hilbert_space(nbit=3):print(mi, end=',')
输出：000,001,010,011,100,101,110,111,
构造叠加态波函数
波函数是各个基向量的线性叠加，系数为复数，称为振幅。
如果给定系数和基向量，可以定义如下函数计算波函数，
def wave_func(coef=[], seqs=[]):'''返回由振幅和几个Qubit序列表示的叠加态波函数，
sum_i coef_i |psi_i> '''
res = 0
for i, a in enumerate(coef):
res += a * basis(seqs[i])
return np.matrix(res)
比如，对于 3个 Qubit 的叠加态波函数，
调用 wave_function 函数，
coef = [1j/np.sqrt(3), np.sqrt(2/3)]seqs = ['000', '100']
s = wave_func(coef, seqs)
s
输出是波函数 在计算机内部的存储格式，其中 1j 是复数中虚数的基本单位 。
matrix([[0. +0.57735027j],
[0. +0.j ],
[0. +0.j ],
[0. +0.j ],
[0.81649658+0.j],
[0. +0.j ],
[0. +0.j ],
[0. +0.j ]])
可以简单的验证，
s.H * s输出：matrix([[1.+0.j]])。
其中 s.H 计算了 的转置共轭 ，也称的厄密特。
s 为 numpy.matrix 数据类型, 内置了 H 属性。
波函数的展开形式
对于一个经过多次操作，已经看不出原形的波函数，可以将其朝各个基向量上投影，找出振幅，并分解。如果某个基为 , 则 朝这个方向投影得到的振幅系数为 。我们来定义一个投影函数，
def project(wave_func, direction):'''
return wave_func.H * direction
通过将波函数 朝各个基向量投影，可以写出人类好识别的形式，
def decompose(wave_func):'''将叠加态波函数分解'''
nbit = int(np.log2(len(wave_func)))
amplitudes = []
direct_str = []
for seq in hilbert_space(nbit):
direct = basis(seq)
amp = project(wave_func, direct).A1[0]
if np.linalg.norm(amp) != 0:
amplitudes.append(amp)
direct_str.append(seq)
return amplitudes, direct_str
decompose(wave_func) 将波函数分解，输出振幅序列和以Qubits序列表示的基向量列表，使用 jupyter notebook 的 Markdown 函数，可以输出漂亮的数学符号格式，
def print_wf(wf):coef, seqs = decompose(wf)
latex = ""
for i, seq in enumerate(seqs):
latex += r"%s$|%s\rangle$"%(coef[i], seq)
if i != len(seqs) - 1:
latex += "+"
return md(latex)
在 print_wf 里面，先将 wf 分解成振幅序列和 Qubit 序列，然后输出，
print_wf(s)输出：
纠缠态
上面例子中的 是叠加态但不是纠缠态。所有纠缠态都不能写成直积形式。再举个非纠缠态的例子，
而下面这个态不能展开为直积形式，为纠缠态，
这是两个费米子的纠缠态，一个自旋向上spin up，一个自旋向下spin down，总自旋为0。但是，没有人知道哪个 spin up，哪个 spin down，系统是两种可能的叠加态。同时，只要测量知道了一个费米子 spin up, 那另一个费米子一定是 spin down。无论对一方测量的时候，另一方离得有多远。它们是纠缠在一起的，总自旋为0。这就是爱因斯坦说的鬼魅般的超距相互作用。
量子门 Gate ：单 Qubit 操作
量子门对单个 Qubit 或两个 Qubit 操作。就像经典计算机上的逻辑门。最简单的量子门是全等门，
, 可以简单的验证，
, 以及
最重要的单量子比特操作的量子门是 NOT 门，简写为 X，
可以简单验证 。
即 NOT 门完成对量子比特的翻转。它有另外一个名字，Pauli 矩阵中的 ，
NOT 门命名为 X 的原因是除了 ，Pauli 矩阵中还有 和 。算上 ，任何 2X2 的厄米特矩阵都可以用它们4个展开。
其实，任何 维厄米特矩阵，都可以用 n 个 pauli 矩阵直积形成的基矩阵展开。
下面使用 np.matrix 函数构造一些常用的量子门，
# Define the one-qubit and 2-qubit operation gatesI = np.matrix("1 0; 0 1")
X = np.matrix("0 1; 1 0")
Y = np.matrix("0 -1j; 1j 0")
Z = np.matrix("1 0; 0 -1")
H = np.matrix("1 1; 1 -1") / np.sqrt(2)
CNOT = np.matrix("1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0")
SWAP = np.matrix("1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1")
gates = {'I':I, 'X':X, 'Y':Y, 'Z':Z, 'H':H, 'CNOT':CNOT, 'SWAP':SWAP}
门翻转 Qubit
举例：计算
print_wf(X * spin_up)返回：
计算
print_wf(X * spin_down)返回：
计算
print_wf(X * (0.6 * spin_up + 0.8 * spin_down))返回：
门带相位翻转
- 计算
返回
- 计算
返回：
门只引起 的相位改变
计算
print_wf(Z*spin_up)返回：
计算
print_wf(Z*spin_down)返回：
Pauli 矩阵做的事情就是将态矢量从 Bloch 球上一个方位旋转到另一个方位，比如 X 门将态矢量绕 x 轴转 180度，Y 门将态矢量绕 y 轴转 180 度，这两个旋转操作都将 z 翻转到了 -z，只不过翻转的路径不一样。如果直接做测量，得到的翻转后的概率是一致的。但是如果把它们翻转的操作当作一个复杂的量子计算机其中的一步，其导致的后续输出是不一样的，因为 X 和 Y 操作输出结果差一个相位 （ 1 + 0j 对比 0 + 1j ，注意python里面复数表示不是i，而是j ）。绕z轴旋转，对 |0> 不会起任何效果，对 |1> 也只有一个相位的变换。
下面是 Bloch 球的示意图，
注意，可以对单个量子比特连续进行多次操作。
Hadamard 门生成叠加态
单 Qubit 除了 NOT 门，还有一个非常常用的 Hadamard 门，用来构造叠加态。
这个门的作用就是把 或 变成叠加态，
返回：0.7071067811865475 +0.7071067811865475
print_wf(H * bit[1])返回：0.7071067811865475 +- 0.7071067811865475
注意这里 bit[0] 跟 spin_up 等价，bit[1] 跟 spin_down等价，H 是 Hadamard 门, 那一长串数字为：
对多个 Qubit 中的一位操作
如果想将 中的 1 翻转为 0，量子门该如何构造？此时虽然看上去是单量子比特操作，而实际上它是多量子比特操作，只不过前两个 0 上的操作为 不变操作 。
是一个 8 维列向量。所以门电路现在是 8 × 8 的矩阵，由 构造，
IIX
array(
[[0, 1, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 1, 0]])
我们将此矩阵作用在 上检查一下，
s001 = basis('001')IIX = kron(I, kron(I, X))
print_wf(IIX * s001)
输出： 。说明 IIX 操作正确，符合预期。
思考：直积求态，再直积求矩阵，最后做一个大矩阵与大向量的点乘非常浪费时间。在玩具模型里，我们能否直接对相应的 Qubit 做操作，跳过 这3个步骤呢？何时可以这样做？
两 Qubit 操作门
经典计算机里有对两个比特的逻辑操作，比如 AND，OR，XOR 操作，那么量子计算机里面有没有对多个量子比特的操作呢？事实证明是有的，比如 CNOT (Controled Not)操作, 就是作用在两个量子位上，它的矩阵形式是，
，两个量子比特的 one-hot 表示是，
CNOT操作示例
- 计算
print_wf(CNOT * s00)
输出：
- 计算
print_wf(CNOT * s01)
输出：
- 计算
print_wf(CNOT * s10)
输出：： ， 此时第一位为 1，第二位翻转。
- 计算
print_wf(CNOT * s11)
输出： ， 此时第一位为 1，第二位也翻转。
可以看出 CNOT 只有作用在 和 上时，第二位才为 。这等价于经典逻辑里的 XOR 操作，也就是单层人工神经网络搞不定，至少需要一个隐藏层才能描述的那个操作。
SWAP门
另一个对两量子比特做操作的门是 SWAP gate，它的作用是将 |01> 翻转为 |10>, 或者将 |10> 翻转成 |01>。其矩阵形式为，
，将SWAP矩阵左乘之前定义的两量子比特的one-hot 矢量表示，
就可以自然的看到SWAP门的效果。
计算
print_wf(SWAP * basis('10'))输出:
纠缠态的制备
对 的左边那位做 Hadamard 门操作变成
然后再接一个 CNOT 门操作，就会制备出纠缠态：
下面两行代码演示了如何用我们自己的玩具量子计算机制备上面这个纠缠态，
# HI 等于 H 与 I 的直积HI = kron(H, I)
print_wf(CNOT * HI * basis('00'))
输出：0.7071067811865475 +0.7071067811865475
记住 0.707... 等于 。
同样的 Hadamard+CNOT 门操作，选择不同的输入，|00>, |01>, |10> 以及 |11>, 将会产生如下4个纠缠态(同时也是贝尔不等式的四个贝尔态），
第四个纠缠态
也称 EPR 态，来自于爱因斯坦和另外两个科学家一起提出的用来推翻量子力学的著名假想实验，EPR佯谬。有了叠加态和纠缠态，就可以试着在量子计算机上做一些有意思的实验，比如 Bell 不等式实验。
更多的门
量子计算门算符是一组复数矩阵，其转制共轭等于其矩阵的逆，即要保证矩阵操作可逆，并且要保证输出结果的各种态的总几率为1 （幺正)。我们可以定义更多标准的量子门算符(参考文献 A Practical Quantum Instruction Set Architecture)，
- Pauli gates:
- Hadamard gate:
- Phase gates:
- Controlled phase gates:
- Cartesian rotation gates:
- Controlled XX gates:
- Swap gates:
总结
用简短的代码制作了一个玩具量子计算机，演示了单个 Qubit，多个 Qubit的制备与门电路操作。介绍了叠加态与纠缠态的概念与制备。
参考文献
- QuTip 文档
- 很早以前写的 hahakity：免费云量子计算机试用指南
- Quantum Computing, From Linear Algebra to Physical Realizations
- M. A. Neilsen and I. L. Chuang, Quantum Computation and QuantumInformation, Cambridge University Press
- https://github.com/adamisntdead/QuSimPy/blob/master/QuSim.py
- corbett/QuantumComputing
- Getting Started with Qiskit
下一篇：TFQ量子计算体验