什么是量子态矢量？

在前面一篇量子系统模拟的博客中，我们介绍了使用python去模拟一个量子系统演化的过程。当我们尝试理解量子态和量子门操作时，可以通过其矩阵形式的运算来描述量子态演化的过程：

\[\left|\psi_t\right>=e^{-iHt}\left|\psi_0\right> \]

这里的狄拉克标记符号和矩阵指数运算，在这篇博客中同样进行了介绍。我们可以简单的将该过程理解为：一个矩阵和一个矢量进行了一个点乘操作，得到了一个更新后的态矢量：

\[\overrightarrow{x’}=\textbf{A}\overrightarrow{x} \]

在量子计算的框架下，由于通用的量子门操作都是酉矩阵(Unitary)，因此不论是更新前还是更新后的操作，得到的态矢量总是归一化的：

\[\left<\psi|\psi\right>=1 \]

因此，在本文中考虑采样时，为了方便计算，我们将态矢量先转换为概率幅矢量，再进行采样。概率幅矢量的特征表现为：

\[\sum_{i=0}^{2^n-1}p_i=1 \]

这里的\(n\)就表示该量子系统的比特数，一个量子系统的量子态元素个数，或者是概率幅的元素个数是比特数的指数倍数(跟量子比特所占用的能级数有关，最常用的是两能级系统，因此为\(2^n\)个元素个数)。

如何实现采样？

在上一个章节中所表述的是量子态的形式，在转换为概率幅矢量之后，其每一个元素都代表获取到当前二进制量子态的概率。这样我们获得一个量子态的态矢量或者概率幅矢量时，其实就是获得了该系统的概率分布。通过该概率分布，我们可以进行蒙特卡罗模拟：先在\([0,1)\)上面进行均匀随机撒点，同时将概率幅矢量转换为其对应的累积分布图，最后计算随机撒的点对应的累积分布图的位置，即可获得当前概率下的模拟采样，具体实现请参考如下示例。

采样示例一

我们先假设一个概率幅的分布，再对其进行采样。

给定一个指数下降的概率幅分布

这里我们先给定一个\(e^{-x}\)的概率分布函数，注意我们采取的是概率幅，因此要对其进行归一化的话只需要计算\(y_i=\frac{y_i}{\sum_jy_j}\)即可。

import numpy as np
import matplotlib.pyplot as plt

plt.figure()
x=[i for i in range(32)]
y=[np.exp(-i) for i in range(32)]
y/=sum(y)
plt.title('Distribution of Quantum States')
plt.xlabel('Quantum States')
plt.ylabel('Probabilities')
plt.xticks(x,('0'*(7-len(str(bin(i))))+str(bin(i)).split('b')[1] for i in range(32)),rotation=70)
plt.plot(x,y,color='black')
plt.show()

Python实现量子态采样
在这个案例中我们还使用了一些matplotlib的特殊绘图技巧，这里我们不展开介绍，后续会单独写一篇文章来分析常用的matploblib画图姿势。

均匀随机数

这里我们直接使用python的random函数，就可以生成\([0,1)\)之间的均匀随机数，撒点数量越多，呈现的均匀分布的结果就越明显。

import random
import matplotlib.pyplot as plt

x = [random.random() for i in range(1000)]
y = [random.random() for i in range(1000)]
plt.figure()
plt.plot(x,y,'.',color='red')
plt.show()

Python实现量子态采样

累积分布函数

所谓的累积分布函数，其实就是将前面获取到的概率幅矢量做一个累积叠加的操作，对应的计算方法如下：

\[y_i=\sum_{j<=i}y_j \]

我们很容易可以预测，在累积分布函数的终点一定是1,这是因为前面所定义的\(\sum_{i=0}^{2^n-1}p_i=1\)。

import numpy as np
import matplotlib.pyplot as plt

plt.figure()
x=[i for i in range(32)]
y=[np.exp(-i) for i in range(32)]
y/=sum(y)
for i in range(len(y)-1):
    y[i+1]+=y[i]
plt.title('Cumulative Distribution Function')
plt.xlabel('Quantum States')
plt.ylabel('Cumulative Probabilities')
plt.xticks(x,('0'*(7-len(str(bin(i))))+str(bin(i)).split('b')[1] for i in range(32)),rotation=70)
plt.plot(x,y,color='black')
plt.show()

Python实现量子态采样

量子态采样

这里我们将概率分布函数和模拟采样结果直接放在一起进行对比：

import numpy as np
import matplotlib.pyplot as plt
import random

plt.figure()
x=[i for i in range(32)]
y=[np.exp(-i) for i in range(32)]
y/=sum(y)
plt.plot(x,y,color='black',label='Real Distribution')
for i in range(len(y)-1):
    y[i+1]+=y[i]

sp=np.zeros(32)
for i in range(50):
    r = random.random()
    for j in range(32):
        if y[j]>r:
            sp[j]+=1/50
            break

plt.title('Sampling Results with 50 times')
plt.xlabel('Quantum States')
plt.ylabel('Probabilities')
plt.xticks(x,('0'*(7-len(str(bin(i))))+str(bin(i)).split('b')[1] for i in range(32)),rotation=70)
plt.plot(x,sp,color='red',label='Simulated Sampling')
plt.legend()
plt.show()

Python实现量子态采样
这里我们可以看到结果已经是非常的接近了，如果我们继续提高采样次数，结果当然会更加接近真实分布：

因此，在获得概率幅之后，我们可以根据场景对精度的要求，对该概率幅进行采样，到这里就完成了所有的功能实现。

采样示例二

除了单调函数外，这里我们再考虑另外一种形式的分布：正弦概率分布函数：

import numpy as np
import matplotlib.pyplot as plt

plt.figure()
x=[i for i in range(32)]
y=[np.sin(i/2)+1 for i in range(32)]
y/=sum(y)
plt.title('Distribution of Quantum States')
plt.xlabel('Quantum States')
plt.ylabel('Probabilities')
plt.xticks(x,('0'*(7-len(str(bin(i))))+str(bin(i)).split('b')[1] for i in range(32)),rotation=70)
plt.plot(x,y,color='black')
plt.show()

Python实现量子态采样
由于上面一个示例我们已经介绍完成了基本的操作流程和原理，这里我们就不过多的赘述，直接展示累积分布函数和最终的模拟采样效果：

累积分布函数

import numpy as np
import matplotlib.pyplot as plt

plt.figure()
x=[i for i in range(32)]
y=[np.sin(i/2)+1 for i in range(32)]
y/=sum(y)
for i in range(len(y)-1):
    y[i+1]+=y[i]
plt.title('Cumulative Distribution Function')
plt.xlabel('Quantum States')
plt.ylabel('Cumulative Probabilities')
plt.xticks(x,('0'*(7-len(str(bin(i))))+str(bin(i)).split('b')[1] for i in range(32)),rotation=70)
plt.plot(x,y,color='black')
plt.show()

Python实现量子态采样

模拟采样结果

import numpy as np
import matplotlib.pyplot as plt
import random

plt.figure()
x=[i for i in range(32)]
y=[np.sin(i/2)+1 for i in range(32)]
y/=sum(y)
plt.plot(x,y,color='black',label='Real Distribution')
for i in range(len(y)-1):
    y[i+1]+=y[i]

sp=np.zeros(32)
for i in range(3000):
    r = random.random()
    for j in range(32):
        if y[j]>r:
            sp[j]+=1/3000
            break

plt.title('Sampling Results with 3000 times')
plt.xlabel('Quantum States')
plt.ylabel('Probabilities')
plt.xticks(x,('0'*(7-len(str(bin(i))))+str(bin(i)).split('b')[1] for i in range(32)),rotation=70)
plt.plot(x,sp,color='red',label='Simulated Sampling')
plt.legend()
plt.show()

Python实现量子态采样