Quantum Computing

前面的“量子因式分解系列”讨论了Shor算法和它的量子电路，今天切换到编程模式，用谷歌的Cirq框架来实现Shor算法，并用它来解决简单的因式分解问题。 导入依赖包 首先还是导入需要的包： import fractions import math import random import sympy import cirq from cirq.contrib.svg import SVGCircuit from typing import Callable, List, Optional, Sequence, Union %matplotlib inline 指数取模运算 还记得Shor算法的量子电路吗？ 首先，我们需要一个“指数取模”的量子变换 \(U_f\) : \begin{matrix} U_f |x\rangle|y\rangle \rightarrow |x\rangle|ya^x \mod N\rangle \end{matrix} 它的输入 \(|x\rangle|y\rangle\) ，输出是 \(|x\rangle|ya^x \mod N\rangle\) 。 在Cirq里面怎么实现这个变换呢？一种办法是按照我之前的文章《聊聊量子因式分解算法的实现》，从加法器和乘法器开始，逐层往上实现。这个工程量就大了，我们也没有真的量子计算机可以跑不是？另一种方法是，利用Cirq框架提供的ArithmeticOperation接口，用经典计算来模拟这个变换： def mod_pow(a:int, x:int, N:int, y:int=1) -> int: p = y for i in range(x): p = (p * a) % N return p class ModularExp(cirq.ArithmeticOperation): def __init__( self, y: Sequence[cirq.Qid], x: Union[int, Sequence[cirq.Qid]], a: int, N: int, ) -> None: if len(y) < N.bit_length(): raise ValueError( f'Register with {len(y)} qubits is too small for N {N}' ) self.y = y self.x = x self.a = a self.N = N def registers(self) -> Sequence[Union[int, Sequence[cirq.Qid]]]: return self.y, self.x, self.a, self.N def with_registers( self, *new_registers: Union[int, Sequence['cirq.Qid']], ) -> 'ModularExp': if len(new_registers) != 4: raise ValueError( f'Expected 4 registers (y, x, a, ' f'N), but got {len(new_registers)}' ) y, x, a, N = new_registers if not isinstance(y, Sequence): raise ValueError(f'y must be a qubit register, got {type(y)}') if not isinstance(a, int): raise ValueError(f'a must be a classical constant, got {type(a)}') if not isinstance(N, int): raise ValueError(f'N must be a classical constant, got {type(N)}') return ModularExp(y, x, a, N) def apply(self, *register_values: int) -> int: assert len(register_values) == 4 y, x, a, N = register_values if y >= N: return y # return (y * a ** x) % N return mod_pow(a=a, x=x, N=N, y=y) def _circuit_diagram_info_( self, args: cirq.CircuitDiagramInfoArgs, ) -> cirq.CircuitDiagramInfo: assert args.known_qubits is not None wire_symbols: List[str] = [] t, e = 0, 0 for qubit in args.known_qubits: if qubit in self.y: if t == 0: if isinstance(self.x, Sequence): e_str = 'e' else: e_str = str(self.x) wire_symbols.append(f'ModularExp(t*{self.a}**{e_str} % {self.N})') else: wire_symbols.append('t' + str(t)) t += 1 if isinstance(self.x, Sequence) and qubit in self.x: wire_symbols.append('e' + str(e)) e += 1 return cirq.CircuitDiagramInfo(wire_symbols=tuple(wire_symbols)) 利用cirq.ArithmeticOperation接口，模拟实现一个指数取模变换ModularExp，这其中的核心就是apply函数： ...

在前面的Shor算法系列中，我讨论了用量子计算进行因式分解的思路。但是这些基本上都是理论上的推导，实际要实现量子因式分解，有哪些难点？解决的思路又是什么？今天就来聊聊我的实现方法。 Shor算法实现的难点：指数取模 在Shor算法中，我们用量子计算来寻找函数 \(f(x)=a^x \mod N\) 的周期（order），并且这一步也是Shor算法中唯一的『量子』部分。 为什么Shor算法这么设计？因为“order finding”这一步目前还没有找到高效的经典算法，而量子计算机则可以用多项式代价求解。 问题来了，这个电路看起来也太『简单』了吧！这么难的一个问题，这样几个电路块就解决了？ 实际上，为了方便讨论算法，上面的这个电路图作了大量的抽象，其中关键的一步，就是把函数 \(f(x)\) 包装成了一个计算『指数取模』的可逆变换 \(U_f\) 。有了这样一个可逆变换，我们才能给它输入 \(2^n\) 个基态的叠加态，实现对经典算法的指数级加速。 这个 \(U_f\) 变换有2个寄存器，分别都是n个量子位： 第一个是寄存器 \(|x\rangle\) ，就是要对a取的指数， 第二个是寄存器是辅助的 \(|1\rangle\) （注意，这也是个n量子位的寄存器）。 在另一端， \(U_f\) 输出 \(|x\rangle\) ，以及 \(|a^x \mod N\rangle\) 。 那这个『指数取模』可逆变换 \(U_f\) 怎么实现呢？ 一个偷懒的方法是，用经典电路实现这样一个变换，把它接入到量子电路里。但是，我们给它的输入是 \(2^n\) 个基态的叠加态，而这个经典电路一次只能计算一个f(x)值，遇到这种叠加态的输入就傻眼了：这不还是 \(O(2^n)\) 的时间复杂度么？运行这样的电路，等于只是用经典电路模拟了Shor算法而已。 所以，这个『指数取模』可逆变换 \(U_f\) 必须是『量子』的，不能用经典电路来实现。这也是Shor算法实现中的一个主要的挑战和难点。 接下来，我会用自顶向下的方式，逐个分解讨论。 量子指数取模 这个『量子指数取模』变换的实现，其实方法不止一种，这里结合我对相关论文的理解，提出我的一个实现方法。 首先， \(a^x \mod N\) 里面，N是给定的数，a是每次运行电路前随机选定的一个数，所以运行时候也是已知的，这两个都可以看作是已知的常量。 然后，变量x是n位二进制数，而我们知道n位二进制数可以表示成这样的形式： \begin{matrix} x = \left( x_{n-1} x_{n-2} \dots x_0 \right)_2 \\ = x_{n-1}2^{n-1} + x_{n-2}2^{n-2} \dots x_02^0 \\ \end{matrix} 其中的 \(x_{n-1}, x_{n-2}, \dots x_0\) 都是0或1。因为x是用在a的指数上， \(a^x \mod N\) 也可以写成： ...

在前文中我初步尝试了Lisp的数值计算和绘图功能，本着“拿着锤子找钉子”的精神，这次我们来搞个实际的应用，用Lisp来演示和可视化量子因式分解算法（Shor​算法）。 首先还是引用numcl和vgplot包： (ql:quickload :numcl) (ql:quickload :vgplot) 假设我们要分解的数是N（前文的例子中N=21），Shor算法需要先在1到N-1之间随机选一个与N互质的数a，然后用量子计算求解函数 \(f(x)=a^x \mod N\) 的周期r。要可视化这个求解的过程，我们得先知道r，所以先写个辅助函数： (defun a-x-mod-N (a x N) (mod (expt a x) N)) ;; Find order for function f(x)=a^x mod N ;; also returns the "modulo sequence" of the function. ;; E.g. a=11, N=21 -> returns order r=6 and sequence (1 11 16 8 4 2) (defun shor-find-mods (a N) (labels ((next-mod (a x N mods) (let ((r (a-x-mod-n a x N))) (if (equal r (car mods)) ; check if the mod sequence is repeating mods ; stop when mod sequence starts repeating (next-mod a ; otherwise continue with a^(x+1) mod N (+ 1 x) N (append mods (list r))))))) (let ((mod-seq (next-mod a 0 N nil))) ;; "values" returns multiple values ;; By default the first value is used, but caller can get other values using ;; multiple-value-bind: ;; (multiple-value-bind (a b) (foo)) (values (length mod-seq) mod-seq)))) a-x-mod-N 函数很简单，就是返回 \(a^x \mod N\) 的值。 ...

量子计算机解决因式分解问题的Shor算法，前面讨论了许多背后的数学原理，未免有点抽象，今天用一个简单的例子来实际演示一下，相信可以对Shor算法有更直观的认识。 因式分解问题：N=21 人肉计算，数字不能用的太大 :-) 这里就假设我们要分解N=21这个数，以它为例子。 根据Shor算法，首先我们要在1到N之间选个随机数a，然后看a和N是否互质： 如果a和N不互质， \(GCD(a, N) \ne 1\) ，那么我们很走运， \(GCD(a, N)\) 就是N的因子之一，问题解决。 如果a和N互质，比如说a=11，继续下一步，用量子计算机寻找函数 \(f(x)=a^x \mod N\) 的周期r，也就是使得 \(1=a^r \mod N\) 的最小正整数r。 在这个例子里，假设我们挑选的随机数a=11，所以我们要找的是这个函数的周期： \begin{matrix} f(x)=11^x \mod 21 \end{matrix} 我们简单算一下就可以知道，这个函数的周期是6: \(x\) \(f(x)=11^x \mod 21\) 0 1 1 11 2 16 3 8 4 4 5 2 6 1 … … 但是在这个例子中，函数周期是未知的，需要用量子计算机来求解。 找函数周期 前文已经讨论过，量子计算寻找函数周期的电路长这样： 这里面的2个寄存器，分别需要多少个量子位呢？ 一般的规则是，寄存器能表达的最大二进制数q，应该不小于N的平方： \(N^2 \le q \le 2N^2\) （这个规则是从概率推导出来的，具体可参见相关教科书）。因此，量子位的数目n，应该不小于 \(\log N^2\) 。在这个例子里面，N=21 ，所以每个寄存器需要9个量子位， \(q=2^9=512\) 。 经过 \(U_f\) 后，我们其实就不关心第2个寄存器 \(|\beta\rangle\) 的状态了，换句话说我们『丢弃』了第2个寄存器。但在量子计算中，『丢弃』实际上也意味着一种『隐性的测量』（implicit measurement），所以，可以认为 \(|\beta\rangle\) 隐性的坍缩到了某个值z。 而这时候的第1个寄存器 \(|\alpha\rangle\) ，必然处于使得 \(11^x \mod 21 = z\) 的所有 \(|x\rangle\) 的叠加态。 ...

在前文『(5.2) Shor算法详解』中，我留了一个尾巴需要补充证明，今天就来填上这个坑，把补丁补上。 在前文讨论概率的推导过程中，有一步看起来不起眼，就是我们需要保证 \(\alpha\) 在0到 \(\pi/2\) 之间（其中 \(\alpha = \frac{A+1}{2} |\theta|\) ）。 但这一步其实很重要，为什么呢？因为只有保证 \(\alpha\) 在0到 \(\pi/2\) 之间，我们才能说 \(\sin \alpha \ge 2\alpha / \pi\) ： 如果 \(\alpha\) 超过了 \(\pi/2\) ，这个性质就不成立了。而这个不等式，是后面我们推导出 \(|\tilde f(y)|^2\) ，也就是Shor算法成功概率的前提。 回顾一下， \(\theta = 2\pi \frac{ry \mod q}{q}\) ，并且y满足不等式 \(|ry \mod q| \le r/2\) 。 我们想要的不等式是： \begin{matrix} \alpha = \frac{A+1}{2} \vert \theta \vert \le \frac{\pi}{2} \end{matrix} 怎样才能证明这个不等式呢？ 我自己研究了一下，应该可以这样证明。 首先来看一下A是怎么计算的。如果我们把 \(l, l+r, \dots, l+Ar\) 画在数轴上，大致是这样的： 所以A的意义是在0到q-1之间，去掉开头的offset和最后的余数t，中间能容纳A个完整的周期（A+1个数据点）。不难得出： \(A=\lfloor \frac{q-l-1}{r} \rfloor - 1\) 。 ...

前文『5.1 Shor算法详解』讨论了在简化条件下，QFT如何提取函数周期信息，但是如果简化条件不满足呢？今天接着这个话题继续讨论。推导过程较长，涉及的知识点也比较多，我会根据我的学习经验，尽可能写的简明易懂。 回顾：简化条件 在前文中，我们构造周期函数 \(f(x)\) 的可逆变换 \(U_f\) ，给它制备适当的输入： 在输出端测量 \(|\beta\rangle\) ，使得它坍缩到某个值z，这时 \(\alpha\) 就会坍缩到 \(|l\rangle, |l+r\rangle, |l+2r\rangle, \dots, |l+Ar\rangle\) 的叠加态： \begin{matrix} |\alpha\rangle = \frac{1}{\sqrt{A+1}} \sum_{j=0}^{A}|jr+l\rangle \end{matrix} 其中\(l\) 称为offset，是满足 \(a^l \mod N = z\) 的最小整数。假设 \(|\alpha\rangle\) 有n个量子位，那么A就是在0到 \(2^n-1\) 之间有多少个周期。 接下来的关键一步，是对 \(|\alpha\rangle\) 作量子傅立叶变换，从中提取函数周期信息。 在前文中，我们假设了一个简化条件，即 \(q=2^n\) 能被函数周期 \(r\) 整除。【注：这里勘误一下，前文中把简化条件写成了“ \(2^n-1\) 能被 \(r\) 整除”，这是笔误】 在这个简化条件下， \(|\alpha\rangle\) 经过量子傅里叶变换(QFT)后再测量，会得到某个值 \(y=k\frac{q}{r}\) ，而这个y值与 \(q=2^n\) 之间的比值，等于某个整数k与周期r的比值，即 \(\frac{y}{q} = \frac{k}{r}\) 。这里y和q都是已知的，化简后得到的分数，其分母就是可能的 \(r\) 值。 那么，如果这个简化条件不成立呢？ 简化条件不成立？ 还是从QFT前的状态开始推导，QFT的作用是什么？ 和前文一样，我们知道把QFT作用在任意态 \(|\psi\rangle\) 上： ...

在前文“回到Shor算法”中我们用QFT提取函数周期信息，但是​这个过程的原理是什么？本文接着详解…… Shor算法关键：如何用QFT提取函数周期信息？ 用量子计算机求解Order Finding问题，实质上是要找到函数 \(f(x)=a^x \mod N\) 的周期。 Shor算法的第一步，是构造函数f(x)对应的可逆变换 \(U_f\) （这里的 \(|x\rangle, |y\rangle\) 都是多个量子位组成的寄存器）： 然后给它制备输入： 在输出端测量 \(|\beta\rangle\) ，使得它坍缩到某个值z，这时 \(\alpha\) 就会坍缩到 \(|l\rangle, |l+r\rangle, |l+2r\rangle, \dots, |l+Ar\rangle\) 的叠加态： \begin{matrix} |\alpha\rangle = \frac{1}{\sqrt{A+1}} \sum_{j=0}^{A}|jr+l\rangle \end{matrix} 其中\(l\) 称为offset，是满足 \(a^l \mod N = z\) 的最小整数。假设 \(|\alpha\rangle\) 有n个量子位，那么A就是在0到 \(2^n-1\) 之间有多少个周期。 这里面， \(l, j\) 都是未知的，并且每次运行这个电路都有可能不同，不能直接从中导出周期r。所以下一步我们要做的，就是用量子傅立叶变换从中提取出有用的周期信息。 接下来关键一步： 如果我们对 \(|\alpha\rangle\) 做量子傅立叶变换，会发生什么？ 先考虑一个简化情况 先考虑一个简化情况： 如果 \(2^n\) 能被 \(r\) 整除，那么我们定义 \(m=2^n/r\) ，显然 \(A=m-1\) ，并且 \(A+1=m=\frac{2^n}{r}\) 。 那么这时的 \(|\alpha\rangle\) 可以写成： ...

既然有『量子傅立叶变换』，那么自然也有它的逆变换，今天就来聊聊这个话题。 量子傅立叶变换和逆变换 前文说过，对于任意状态 \(|\psi\rangle\) ： \begin{matrix} |\psi\rangle = \sum_{j=0}^{N-1}f(j)|j\rangle & (N=2^n) \end{matrix} 对它进行量子傅立叶变换QFT后，会得到一个新状态： \begin{matrix} |\tilde\psi\rangle = QFT|\psi\rangle = \sum_{k=0}^{N-1}\tilde f(k)|k\rangle, \\\ \\\ \tilde f(k) = \frac{1}{\sqrt{N}} \sum_{j=0}^{N-1} e^{2\pi ijk/N} f(j) \\\ = \frac{1}{\sqrt{N}} \sum_{j=0}^{N-1} \omega^{jk} f(j) & (\omega=e^{2\pi i/N}) \\\ \end{matrix} 那么量子傅里叶逆变换 \(QFT^{-1}\) 是什么样的？ 很简单，其实就是把e的指数加个负号： \begin{matrix} |\phi\rangle = \sum_{k=0}^{N-1} \tilde f(k) |k\rangle \\\ \\\ QFT^{-1} |\phi\rangle = \sum_{j=0}^{N-1} g(j) |j\rangle \\\ \\\ g(j) = \frac{1}{\sqrt{N}} \sum_{k=0}^{N} \omega^{-jk} \tilde f(k) \\\ \end{matrix} 可以看到，量子傅立叶逆变换和QFT非常相似，唯一的区别就是e的指数上多了个负号。 ...

最近中科院物理所发布了一个量子计算云平台ScQ，是一个10量子位的超导量子计算机，可以在云端免费用。看到这个新闻很兴奋，忍不住马上去试用了一下。 点击云平台的链接（http://q.iphy.ac.cn/），进去以后是这样的： 点击"Try ScQ Cloud"，进入量子电路页面： 先来试试精度，用2个量子比特，测量它们的初始状态（应该都是0）： 点击Submit，让我耐心等待结果： 然后结果是这样： 测量100次，96%情况下符合预期（00 ），有3%的情况下第2个量子位测到1… 试试H门产生叠加态： 预期是4个状态各25%，这个误差有点高……感觉可能我的取样次数（100）太少，增加到1000次，结果就好了很多： 再来试试Deutsch算法：f(x)=x, 预期输出1 (balanced) 取样1000次，92%的情况下得到正确的结果。 取样10000次，正确率达到接近97%： 第一次试用国产量子计算云平台的感受： 用户体验不错，界面直观易用，一键就可开始搭建电路 目前只支持10个量子位和基本的量子逻辑门，功能还比较有限 从我运行的结果看，误差还是有点高 提交电路后，有时候会返回运行错误，需要重试 为国产量子计算平台的进展感到高兴，希望将来能看到更多的消息！

之前的Shor算法系列，讲了利用量子傅立叶变换QFT来提取函数周期信息。这次歪个楼，聊聊量子傅立叶变换。 假如我们用n个量子位搭了一个QFT电路（后文会讲到如何实现），如何解读它的输入和输出呢？ 假设我们的电路输入状态是某个基态 \(|j\rangle\) ，它是由n个量子位组成的，我们把它展开成二进制表示的话： \begin{matrix} |j\rangle = |j_{n-1}j_{n-2}\dots j_0\rangle \\\ j = j_{n-1}2^{n-1} + j_{n-2}2^{n-2} + \dots + j_0 2^0 \end{matrix} 经过QFT，我们得到一个新的状态： \begin{matrix} QFT|j\rangle = \frac{1}{\sqrt{2^n-1}} \sum_{k=0}^{2^n-1} \omega^{jk} |k\rangle \\\ = \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{2\pi ijk/2^n} |k\rangle \end{matrix} 其中每个分量k的相位 \(e^{2\pi ijk/2^n}\) 。而因为其中的 \(k\) 也是一个n位二进制数，也就是说： \begin{matrix} k = k_{n-1}2^{n-1} + k_{n-2}2^{n-2} + \dots + k_0 2^0 \end{matrix} 那么这个相位里，是把k除以 \(2^n\) ，我们把k的二进制展开代入进去就变成了这样： \begin{matrix} k / 2^n = \frac{k_{n-1}}{2^1} + \frac{k_{n-2}}{2^2} + \dots + \frac{k_0}{2^n} \end{matrix} 等一下，这个式子是不是似曾相识？ ...