用cuQuantum框架进行量子计算模拟

Simulate Quantum Computation Using cuQuantum Framework

摘要: 量子计算模拟是使用经典计算机和相应框架对量子计算进行模拟，是设计和实现量子算法的有力工具。cuQuantum作为模拟量子计算的框架，为使用GPU进行量子计算模拟提供了重要工具，但在使用该框架时，由于表示法与直观有一定差距，带来了使用上的一些问题，通过研究和学习该框架，以及使用该框架对量子线路的基本操作进行模拟，澄清了框架使用过程中的相关问题，为后续设计和实现量子算法提供基础。

Abstract: The simulation of quantum computation is a powerful tool for designing and implementing quantum algorithms by using classical computers and corresponding frameworks. The cuQuantum, as a framework for simulating quantum computing, provides an important tool for quantum computing simulation using GPU. However, since there is a certain gap between representation and intuition when using this framework, it brings some problems in the usage. By studying and learning the framework and using the framework to simulate the basic operation of quantum circuits, the related problems in the usage of the framework are clarified, which provides a basis for the subsequent design and implementation of quantum algorithms.

1. 引言

量子计算是利用量子力学原理对经典计算模型的一种改进，有望突破强Church-Turing论题，实现对复杂计算问题的超多项式级加速 [1]。量子计算的提出可以追溯到费曼，但证明量子计算优势的成果当数Shor算法 [2] ——应用此算法可以快速分解大整数，在此前后相继出现了Simon算法 [3] 和Grover算法 [4]，都表明量子算法比经典算法效率更高。

量子算法要量子计算机才能执行，一般情况下，获取量子计算机过于困难，虽然IBM等公司提供了量子计算机的云接口，但开放的接口只是不到10个量子比特的量子计算机，对于量子算法特别是量子智能算法的研究无疑是不够的。要想更好设计和实现量子算法，就需要使用模拟的方法，在经典机器上仿真、调试，达到要求后，再使用量子计算机运行，这样设计和实现算法的效率更高。NVIDIA公司于2021年年末推出了cuQuantum框架，使用该框架和相应的GPU进行量子计算模拟将大大提高效率。

2. 量子计算的基本操作

2.1. 量子比特

在电子计算机中，我们用0、1表示计算机的状态，信息被编码为0、1所组成的字符串，而相应的0、1被称为比特；在量子计算中我们使用叠加态

$|s\rangle =a|0\rangle +b|1\rangle$，其中a、b为复数，满足 ${\left|a\right|}^2+{\left|b\right|}^2=1$

表示量子计算机的状态，a、b表示概率幅，其模的平方 ${\left|a\right|}^2$ 、 ${\left|b\right|}^2$ 代表测量后得到状态 $|0\rangle$ 、 $|1\rangle$ 的概率。我们也可将 $|s\rangle$ 用基和坐标形式表示为：

$|s\rangle =\left(|0\rangle ,|1\rangle \right)(\; a\; b\; )$

在表示多位量子比特和表示对多位量子比特进行操作时涉及张量积，使用坐标表示法更简洁，我们现在来看看如何表示多位量子比特。

第一种情况，多位量子比特可由单个量子比特组合而成，我们以两位量子比特为例：

$|s_0{s}_1\rangle =|s_0\rangle |s_1\rangle =|s_0\rangle \otimes |s_1\rangle$

其中

$|s_0\rangle =\left(|0\rangle ,|1\rangle \right)\left(\begin{array}{c}{a}_0\\ b_0\end{array}\right)$， $|s_1\rangle =\left(|0\rangle ,|1\rangle \right)\left(\begin{array}{c}{a}_1\\ b_1\end{array}\right)$

可以得出

$\begin{array}{c}|s_0{s}_1\rangle =|s_0\rangle |s_1\rangle =|s_0\rangle \otimes |s_1\rangle =\left[\left(|0\rangle ,|1\rangle \right)\left(\begin{array}{c}{a}_0\\ b_0\end{array}\right)\right]\otimes \left[\left(|0\rangle ,|1\rangle \right)\left(\begin{array}{c}{a}_1\\ b_1\end{array}\right)\right]\\ =\left[\left(|0\rangle ,|1\rangle \right)\otimes \left(|0\rangle ,|1\rangle \right)\right]\left[\left(\begin{array}{c}{a}_0\\ b_0\end{array}\right)\otimes \left(\begin{array}{c}{a}_1\\ b_1\end{array}\right)\right]=\left(|00\rangle ,|01\rangle ,|10\rangle ,|11\rangle \right)\left(\begin{array}{c}\begin{array}{c}{a}_0{a}_1\\ a_0{b}_1\end{array}\\ \begin{array}{c}{b}_0{a}_1\\ b_0{b}_1\end{array}\end{array}\right).\end{array}$

第二种情况，多位量子比特处于量子纠缠态，这时的多位量子比特不能写为单个量子比特的张量积形式，我们以如下处于纠缠态的量子比特为例：

$|s_0{s}_1\rangle =\frac{1}{\sqrt{2}}|00\rangle +\frac{1}{\sqrt{2}}|11\rangle$

由于这种形式的量子比特无法写为两个独立的量子比特张量积的形式，用基与坐标方式可将其表示为：

$|s_0{s}_1\rangle =\left(|00\rangle ,|01\rangle ,|10\rangle ,|11\rangle \right)\left(\begin{array}{c}\begin{array}{c}\frac{1}{\sqrt{2}}\\ 0\end{array}\\ \begin{array}{c}0\\ \frac{1}{\sqrt{2}}\end{array}\end{array}\right)$.

2.2. 量子门

在经典计算中，逻辑门的作用是对比特进行逻辑运算，在量子计算中，我们要对量子比特进行操作就需要有量子门，将量子门作用于量子比特，即可完成相应的量子计算操作，而量子门对量子比特的作用一般都可视为线性变换。

对于单量子比特的量子操作，我们可使用线性代数中的表示法，假设我们有一个幺正算符U，将此算符作用于一个单量子比特，就相当于将此幺正算符对应的矩阵乘上概率幅向量，具体表示如下：

$U|s\rangle =\left(|0\rangle ,|1\rangle \right)U(\; a\; b\; )$

对于多量子比特构成系统的量子操作，通常情况下可仿照单量子比特量子操作的表示法表示，假设有一个幺正算符U，将其作用于一个双量子比特，可以表示为：

$U|s_0{s}_1\rangle =\left(|00\rangle ,|01\rangle ,|10\rangle ,|11\rangle \right)U\left(\begin{array}{c}\begin{array}{c}{a}_0{a}_1\\ a_0{b}_1\end{array}\\ \begin{array}{c}{b}_0{a}_1\\ b_0{b}_1\end{array}\end{array}\right)$

这种表示法适用于所有多量子比特系统，当多量子比特系统中量子操作相互独立，即可以将多量子操作拆分为几个单量子操作时，或是由多个单量子操作组合为一个多量子操作时，我们可以用张量积的方式表示，我们还是以双量子比特系统为例，在第一个量子比特 $|s_1\rangle$ 上施加了一个U₁操作，在第二个量子比特 $|s_2\rangle$ 上施加一个U₂操作，即

$\begin{array}{c}U|s_1{s}_2\rangle =\left(U_1|s_1\rangle \right)\left(U_2|s_2\rangle \right)=\left(U_1|s_1\rangle \right)\otimes \left(U_2|s_2\rangle \right)=\left[\left(|0\rangle ,|1\rangle \right)U_1\left(\begin{array}{c}{a}_1\\ b_1\end{array}\right)\right]\otimes \left[\left(|0\rangle ,|1\rangle \right)U_2\left(\begin{array}{c}{a}_2\\ b_2\end{array}\right)\right]\\ =\left[\left(|0\rangle ,|1\rangle \right)\otimes \left(|0\rangle ,|1\rangle \right)\right]\left(U_1\otimes U_2\right)\left[\left(\begin{array}{c}{a}_1\\ b_1\end{array}\right)\otimes \left(\begin{array}{c}{a}_2\\ b_2\end{array}\right)\right]\\ =\left(|00\rangle ,|01\rangle ,|10\rangle ,|11\rangle \right)\left(U_1\otimes U_2\right)\left(\begin{array}{c}\begin{array}{c}{a}_1{a}_2\\ a_1{b}_2\end{array}\\ \begin{array}{c}{b}_1{a}_2\\ b_1{b}_2\end{array}\end{array}\right)\end{array}$

故上述操作可写为两个操作对应矩阵的张量积表示。这种方式的操作合成可以推广到任意量子比特组合构成的量子线路的情况，这里不再讨论。

常用的量子门及其矩阵表示如表1和表2所示。

Table 1. Common single qubit quantum gates and their matrix representation [5]

表1. 常用单比特量子门及其矩阵表示 [5]

Table 2. Common multi-qubit quantum gates and their matrix representation [5]

表2. 常用多比特量子门及其矩阵表示 [5]

2.3. 测量

要得到最终结果，我们需要将量子比特变换为比特，从而沟通量子系统和经典系统，这样才能实现量子计算的目的。在量子力学中，我们对处于量子态的粒子进行测量，其量子态就会按照概率幅所确定的概率坍缩为某一种状态，而之后都以这种状态存在。在量子计算中，对量子比特进行测量，其量子比特就将以概率幅模的平方为概率变换为对应的比特。要对测量操作进行模拟，在模拟过程中就需要满足：以概率幅确定的概率来随机决定测量后的状态；测量后的状态一直保持不变，直到下次变换。在量子线路中，我们常常将测量操作作为量子计算的最终步骤。

在模拟测量操作时，我们先要获得量子比特 $|s\rangle =a|0\rangle +b|1\rangle$ 坍缩为某一状态 $|m\rangle$ 的概率(m取0、1两值)：

$p\left(m\right)=\langle s\left|\frac{1}{2}\left[I+{\left(-1\right)}^{m}Z\right]\right|s\rangle$

接着用随机数生成器产生一个[0, 1)区间的数q，比较q与 $p\left(m\right)$ 的关系：若 $q<p\left(m\right)$，则取状态 $|m\rangle$ ；否则，取状态 $|1-m\rangle$。测量之后量子比特 $|s\rangle =a|0\rangle +b|1\rangle$ 将坍缩为

$\frac{\frac{1}{2}\left[I+{\left(-1\right)}^{m}Z\right]|s\rangle }{\sqrt{p(\; m\; )}}$

3. cuQuantum框架 [6]

从上面的描述可以看出，量子计算主要使用矩阵，而这正是GPU处理的强项，使用GPU模拟量子计算是较好的选择。NVIDIA公司于2021年年底推出的cuQuantum，正是基于GPU的量子计算模拟框架，它分为cuStateVec和cuTensornet两个部分，cuStateVec是处理量子门运算的库，相对比较底层，也较为灵活，我们以下的量子计算模拟主要使用这个库；而cuTensornet主要负责张量网络计算，本文不涉及此库的相关内容。

由于cuQuantum框架是基于GPU进行量子计算模拟的框架，它建立在已有的CUDA基础上，其CPU、内存与GPU的数据传输与管理操作都由CUDA提供，cuQuantum仅提供模拟量子计算部分的API，本文第4节将使用cuStateVec库的部分API进行量子计算模拟。

4. 量子计算模拟

量子计算模拟是指使用经典计算设备对量子计算机进行模拟，主要模拟量子线路和量子操作。现阶段量子算法的描述以量子线路为主，量子线路是由量子门组成，为设计和实现量子算法提供便利，故在cuQuantum中所有操作都以量子门形式给出。本文的量子计算模拟环境如表3所示。

Table 3. Parameters of the simulation environment in this paper

表3. 本文模拟环境参数

4.1. 量子比特的cuQuantum框架表示

首先我们来看看如何用cuQuantum框架表示一个量子比特，假设有一个量子比特 $|s\rangle =2.3|0\rangle +1.4i|1\rangle$，用框架可表示为：

cuDoubleComplex h_sv[] = {{2.3, 0.0}, {0.0, 1.4}};

这里注意系数，根据定义一个量子比特在各个分量的系数都为复数，在cuQuantum中我们用{x, y}来表示复数x + yi，故系数2.3表示为{2.3, 0.0}，系数1.4i表示为{0.0, 1.4}。

接下来我们来看如何表示多个量子比特，一般来说，多个量子比特可表示为单量子比特的张量积，如第2节中的描述，但在cuQuantum中并非以张量积的方式来表示多个量子比特，而是用一个系数向量来表示，并以小端序进行编码，譬如：我们有一个2位的量子比特 $|s_0{s}_1\rangle =0.6|00\rangle +|01\rangle +0.3|10\rangle +0.5|11\rangle$，用框架表示为：

cuDoubleComplex h_sv[] = {{0.6, 0.0}, {0.3, 0.0}, {1, 0.0}, {0.5, 0.0}};

虽然在框架中，系数可以是任意的复数，但在实际表述时，应当将其归一化，一种常用的做法是仅对最终结果进行归一化处理。为了方便起见，以下论述中所用示例我们均没有进行归一化处理(cuQuantum框架没有强行要求归一化)。

4.2. 对单量子比特操作的模拟

在量子计算中最常用的量子门是Hadamard门，Pauli门，Phase门以及π/8门，如第1节所述，Pauli门包含3个分别用X、Y、Z来表示。(以下描述中，*表示乘，sqrt表示平方根)

首先我们来看Hadamard门的cuQuantum框架表示：

cuDoubleComplex matrix_H[] = {{1.0*sqrt(2)*0.5, 0.0}, {1.0*sqrt(2)*0.5, 0.0},

{1.0*sqrt(2)*0.5, 0.0}, {−1.0*sqrt(2)*0.5, 0.0}};

现在我们将此Hadamard门应用于一个单量子比特 $|s\rangle =2.3|0\rangle +1.4|1\rangle$，用cuQuantum框架对此量子比特进行操作，可表示为：

custatevecApplyMatrix(handle, d_sv, CUDA_C_64F, 1, matrix_H, CUDA_C_64F,

CUSTATEVEC_MATRIX_LAYOUT_ROW, 0, {0}, 1, {},

0, nullptr, CUSTATEVEC_COMPUTE_64F,

extraWorkspace, extraWorkspaceSizeInBytes);

d_sv是显存中的h_sv表示，可以直接用相应函数将h_sv复制到d_sv中，由于本文不讨论显卡的设备管理问题(如设备空间的分配、GPU与CPU的数据传输等问题)，故具体细节不在这里描述。上述函数使用后，d_sv数组就是计算结果(上例中d_sv = {{2.616295, 0.0}, {0.636396, 0.0}}，按照cuQuantum的表示法，可以看出最终结果为 $\text{2}.\text{616295}|0\rangle +0.\text{636396}|1\rangle$)。

其他几个单量子比特操作我们总结如表4所示。

Table 4. Operation simulation of single qubit quantum gate

表4. 单比特量子门的操作模拟

4.3. 对多量子比特操作的模拟

一般情况下，多量子比特是多个量子比特的组合，其表示法可以被看作多个量子比特的张量积，其对应操作可以表示为多个矩阵的张量积，第1节中有详细描述，用cuQuantum框架中对单量子比特操作模拟的做法，对大多数多量子比特的操作模拟而言都不存在问题，只需要记住多量子比特的cuQuantum表示法，下面以对两个量子比特的操作为例：

假设有 $|s_0{s}_1\rangle =0.6|00\rangle +|01\rangle +0.3|10\rangle +0.5|11\rangle$，我们要对第一个量子比特位进行X操作，对第二个量子比特位进行Y操作，则可以计算X与Y的张量积，将此张量积作用于上述量子比特，即：

$\begin{array}{c}\left(X\otimes Y\right)|s_0{s}_1\rangle =\left(|00\rangle |01\rangle |10\rangle |11\rangle \right)\left[\begin{array}{cccc}0& 0& 0& -i\\ 0& 0& i& 0\\ 0& -i& 0& 0\\ i& 0& 0& 0\end{array}\right]\left(\begin{array}{c}\begin{array}{c}0.6\\ 1\end{array}\\ \begin{array}{c}0.3\\ 0.5\end{array}\end{array}\right)\\ =-0.5i|00\rangle +0.3i|01\rangle -i|10\rangle +0.6i|11\rangle \end{array}$

在使用cuQuantum框架时，要将X、Y的位置进行交换，得到

$Y\otimes X=\left[\begin{array}{cccc}0& 0& 0& -i\\ 0& 0& i& 0\\ 0& -i& 0& 0\\ i& 0& 0& 0\end{array}\right]$

再将其作用于 $|s_0{s}_1\rangle$ 的cuQuantum表示法，即此时的量子操作将表示为

cuDoubleComplex matrix[] = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, −1.0},

{0.0, 0.0}, {0.0, 0.0}, {0.0, −1.0}, {0.0, 0.0},

{0.0, 0.0}, {0.0, 1.0}, {0.0, 0.0}, {0.0, 0.0},

{0.0, 1.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};

运用custatevecApplyMatrix函数后的结果为{{0.0, −0.5}, {0.0, −1.0}, {0.0, 0.3}, {0.0, 0.6}}，这里的结果是cuQuantum表示法，它表示 $-0.5i|00\rangle +0.3i|01\rangle -1.0i|10\rangle +0.6i|11\rangle$。

对于Swap门，在使用框架时，不需要作改变即可使用，即Swap门的变换矩阵可直接表示为

cuDoubleComplex matrix[] = {{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},

{0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0},

{0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},

{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}};

在量子计算中，受控操作是很重要的一类多量子比特操作，cuQuantum框架中通过设置控制位的方式，简化受控操作的矩阵表示，这里以Toffoli门为例，我们知道Toffoli门表示的量子操作是“如果两个控制位同时为1，则对第三个位取反，否则不改变第三个位”，对一个量子位进行取反的操作用cuQuantum框架可表示为

cuDoubleComplex matrix[] = {{0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}};

这时我们把控制位定为0和1，目标位定为2，运用custatevecApplyMatrix函数，这个操作可写为：

custatevecApplyMatrix(handle, d_sv, CUDA_C_64F, 3, matrix, CUDA_C_64F,

CUSTATEVEC_MATRIX_LAYOUT_ROW, 0, {2}, 1, {0, 1},

2, nullptr, CUSTATEVEC_COMPUTE_64F,extraWorkspace, extraWorkspaceSizeInBytes);

d_sv在函数调用时传入

cuDoubleComplex h_sv[] = {{0.0, 0.0}, {0.0, 0.1}, {0.1, 0.1}, {0.1, 0.2},

{0.2, 0.2}, {0.3, 0.3}, {0.3, 0.4}, {0.4, 0.5}};

传入的参数用标准法表示为

$\begin{array}{c}\left(0.1i\right)|100\rangle +\left(0.1+0.1i\right)|010\rangle +\left(0.1+0.2i\right)|110\rangle +\left(0.2+0.2i\right)|001\rangle \\ +\left(0.3+0.3i\right)|101\rangle +\left(0.3+0.4i\right)|011\rangle +\left(0.4+0.5i\right)|111\rangle ,\end{array}$

运算后的结果为

{{0.0, 0.0}, {0.0, 0.1}, {0.1, 0.1}, {0.4, 0.5}, {0.2, 0.2}, {0.3, 0.3}, {0.3, 0.4}, {0.1, 0.2}}，

可用标准形式写为：

$\begin{array}{c}\left(0.1i\right)|100\rangle +\left(0.1+0.1i\right)|010\rangle +\left(0.4+0.5i\right)|110\rangle +\left(0.2+0.2i\right)|001\rangle \\ +\left(0.3+0.3i\right)|101\rangle +\left(0.3+0.4i\right)|011\rangle +\left(0.1+0.2i\right)|111\rangle .\end{array}$

4.4. 对测量操作的模拟

在量子计算中，要得到最终结果，就要对量子线路进行测量，在第1节中我们从概念上描述了测量操作，现在来看用cuQuantum框架如何实现测量操作。

cuQuantum框架提供了多个关于测量的函数，最常用的是custatevecMeasureOnZBasis函数，它使用Pauli-Z矩阵为基，可以达到模拟量子计算中测量操作的效果，下面以 $|s\rangle =\left(0.5+0.1i\right)|0\rangle +\left(0.4+0.3i\right)|1\rangle$ 为例，我们首先准备一个数组h_sv：

cuDoubleComplex h_sv[]= {{0.5, 0.1}, {0.4, 0.3}};

接下来产生一个[0, 1)中随机数randnum：

srand(time(0));

const double randnum = rand()/(RAND_MAX+1.0);

调用custatevecMeasureOnZBasis函数：

custatevecMeasureOnZBasis(

handle, d_sv, CUDA_C_64F, nIndexBits, &parity, basisBits, nBasisBits,

randnum, CUSTATEVEC_COLLAPSE_NORMALIZE_AND_ZERO);

就完成了对 $|s\rangle$ 的测量，测量结果在partiy变量中，测量之后的 $|s\rangle$ 用d_sv数组表示。要进一步使用测量后的结果直接使用d_sv数组即可，这也就模拟了测量后量子态的变化。表5显示了对 $|s\rangle$ 独立测量5次，每次测量的结果。(在程序设计时，不能连续测量，若连续测量——调用多次上述的测量函数，其结果在第一次测量后就已经确定，后4次其结果都为第1次的测量后结果，这也就是测量操作的特点。

Table 5. Measure $|s\rangle =\left(0.5+0.1i\right)|0\rangle +\left(0.4+0.3i\right)|1\rangle$

表5. 对 $|s\rangle =\left(0.5+0.1i\right)|0\rangle +\left(0.4+0.3i\right)|1\rangle$ 进行测量

5. 结论与展望

本文研究了使用cuQuantum框架进行量子计算模拟的方法，运用框架对常用的量子操作进行了模拟。下一步研究主要集中在使用cuQuantum框架实现和仿真量子计算智能和量子强化学习算法上。

基金项目

本文系四川开放大学科研课题“基于量子计算智能的组合优化问题研究”(课题编号：KTGCJS2021002Y)阶段性成果。