简介:QCF(Quantum Computing Functions)是一套面向Matlab和Octave环境的开源量子计算函数包,以Nielsen & Chuang《量子计算和量子信息》为理论基础,适合正在学习量子计算的学生、研究人员以及希望动手验证算法的工程师。压缩包共40个文件,其中37个.m函数文件构成核心代码,另有1份PDF、1份Markdown说明和1份Word使用指南,整包仅132KB,轻量易用。目前已有186人学习下载。工具包覆盖量子比特初始化、Hadamard/Pauli/CNOT等基本量子门、量子电路构造,并实现Deutsch-Jozsa、Grover搜索、量子傅里叶变换等经典算法;配套的测量验证、密度矩阵与Bloch球可视化功能,可帮助使用者在无真实量子硬件时完成模拟实验,直观理解叠加与纠缠等核心概念,为深入研究量子计算打下实践基础。
1. QCF:把Nielsen & Chuang的量子算法搬进Matlab/Octave
第一次跑通 Deutsch-Jozsa 算法时,我盯着 QCF 输出的向量结果愣了几秒——一个不到 30 行的函数,就完成了教科书第 6 章里那个让经典计算机至少需要两次查询才能区分的问题。量子计算的理论门槛高,但用 QCF(Quantum Computing Functions)在 Matlab/Octave 里复现 Nielsen & Chuang 的算法,门槛低得多。它不是要替代 Qiskit,而是把量子态和量子门直接映射成矩阵与向量运算,让你在熟悉的数值环境里验证每一个幺正变换。这个工具包适合两类人:一是正在啃原书的研一学生,想边读边敲代码验证贝塞尔不等式和测量公设;二是做信号处理或优化的工程师,想快速评估某个量子算法在自己问题上的行为。QCF 的每个函数都能独立调用,不依赖云服务和硬件厂商 SDK,装上就能跑。
2. 量子态、门与电路:QCF 的核心函数如何映射到矩阵运算
2.1 量子比特的向量表示与bin2vec的状态编码
在 Nielsen & Chuang 的记号里,单量子比特是两个复振幅的列向量,多量子比特是张量积。QCF 在 Matlab 里直接遵循这个约定:|0>是[1;0],|1>是[0;1],两个比特的|00>是kron([1;0],[1;0])。一开始不太习惯的是,QCF 没有用 Python 那种类对象来包装量子态,所有函数都接收普通浮点复矩阵。好处是你能在命令行直接disp看每个中间态,排错非常直观。
bin2vec.m把二进制的 0/1 字符串或行向量转成列向量。
% 将 '101' 转为 8 维列向量,对应 |101> psi = bin2vec('101'); % psi 是 8x1 的稀疏列向量 disp(psi') % 显示非零位置,0 1 2 3 4 5 6 7 中下标 5 处为 1这个函数的运算逻辑是:先把二进制字符串按位拆开,依次做张量积。'101'会先构造[0;1](因为最低位是 1),然后与[1;0]张量,再与[0;1]张量,得到kron(kron([0;1],[1;0]),[0;1])。如果你用dec2vec则可以直接从一个十进制整数生成状态向量,dec2vec(5,3)生成三比特|101>。这两个函数在处理初始化和验证结果时特别有用,尤其当你想手动构造一个均匀叠加态作为算法输入。
2.2 量子门:hadamard.m、build_u.m与自带门矩阵
QCF 把常用量子门封装成了函数。hadamard.m调用后直接返回作用在指定量子比特上的矩阵,identity.m生成单位门。但工具包里并没有逐门枚举 Pauli-X/Y/Z,而是提供build_u.m这个更通用的构造器。
% 构造一个作用于第 1 个量子比特的 Hadamard 门(共有 2 个量子比特) H1 = build_u(hadamard(), 1, 2); % 构造 CNOT:控制位是 1,目标位是 2,两比特 CNOT = build_u([1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0], [1 2], 2);build_u的第三个参数是总比特数。它做的事情相当于把本地门矩阵嵌入到完整希尔伯特空间的张量积里。第一个参数如果传的是hadamard()的 2x2 矩阵,就等价于在该比特位上作用 H 门,其他比特位保持单位。第二个参数可以是标量或向量:标量表示该门只作用在一个比特上;向量如[1 2]表示这是一个多比特门,矩阵维度要跟比特数匹配。这里有个容易踩的坑:QCF 的比特顺序是从 1 开始,并且低索引对应态矢的高权重位置。也就是说bin2vec('01')的第一个比特是 0,第二个比特是 1,和build_u里的序号对应关系需要先做一次disp验证。
门矩阵本身也支持你手写。Pauli-X 就是[0 1; 1 0],Pauli-Z 是[1 0; 0 -1]。QCF 不限制你只用它内置的门,build_u完全接受任意满足幺正性的矩阵,所以你可以在 Matlab 命令窗口里用kron自己组合出受控-U 门,再传给后续函数。
2.3f.m与f_c0.m:Oracle 函数怎么在模拟器里表示
QCF 深受原书公式编号影响。f.m是一个通用函数模板,用来表示布尔函数 f(x) 对输入 x 的映射。在 Deutsch-Jozsa 和 Grover 算法里,Oracle 会以相位翻转或置位翻转的形式出现。f_c0.m和f_c1.m是常数函数 0 和 1 的具体实现,f_b0.m和f_b1.m是平衡函数的实现。
% 调用常数函数 f(x)=0,输入是量子态向量 x = bin2vec('00'); [fval, phase] = f_c0(x); disp(fval) % 0Oracle 在 QCF 中的约定值得先说清楚:f_c0/1和f_b0/1直接返回函数值fval,同时phase输出是(-1)^f(x)的相位因子。在模拟算法时,你会发现很多算法实现里调用的是phase,因为你真正需要的其实是相位反冲,而不是经典函数值。这个设计和 Nielsen & Chuang 书中式 (2.32) 关于 oracle 的定义高度一致。
2.4 用vec2struct与renormalise管理多边态
struct2vec.m和vec2struct.m是 QCF 的两个辅助函数,用来在结构体形式和向量形式之间转换状态。实操中,当你跑完一个算法得到输出向量,想从里面提取某个寄存器对应的约化态时,vec2struct会把向量按比特拆分。
psi = bin2vec('101'); parts = vec2struct(psi); % 拆分后的结构体 % 再归一化 psi_renorm = renormalise(psi);renormalise.m处理的是浮点误差导致的模长不再严格等于 1 的问题。量子模拟跑几十个门之后,振幅的微小漂移很正常,但后续测量和概率计算需要归一化,这个小函数就是专门干这个的。建议每个算法主流程结束之前都调用一次renormalise,能省掉很多奇奇怪怪的数值警告。pretty.m则是把小数形式的复数振幅转成0.7071i这种更接近书面的显示,方便和纸面推导对照。
3. Deutsch-Jozsa 算法复现:从黑盒函数到量子电路判定
3.1 算法原理与 QCF 中的电路模块
Deutsch-Jozsa 要解决的问题是:给定一个未知的布尔函数 f: {0,1}^n → {0,1},保证它要么是常数,要么是平衡的(恰好一半输出 0 一半输出 1),用最少的 oracle 查询判断它是哪一类。经典算法最坏需要 2^(n-1)+1 次查询,量子算法只需要一次。
量子电路的标准做法是:把所有输入比特初始化为|0>,工作比特初始化为|1>;分别施加 H 门后得到均匀叠加;然后调用 oracle 实现相位反冲;再对输入比特施加 H 门;最后测量输入寄存器,如果全为 0 则函数是常数,否则是平衡的。QCF 提供了一整套分步函数,deutsch_jozsa.m是完整流程,dj_c0.m、dj_c1.m、dj_b0.m、dj_b1.m分别对应四种不同函数的具体电路。
3.2 自己搭一遍deutsch_jozsa.m
与其直接用封装好的deutsch_jozsa,我更建议先按电路顺序把中间步骤写出来,既能验证理解,也方便调整 oracle。下面的代码就是 QCF 风格下手工实现 Deutsch-Jozsa 的完整流程,我用的是 2 个输入比特加 1 个工作比特:
% 三比特:q1, q2 作为输入寄存器,q3 作为辅助寄存器 n = 2; % 初始态 |00>|1> psi = kron(bin2vec('00'), bin2vec('1')); % 对三个比特都施加 H 门 H3 = build_u(hadamard(), 1, 3) * build_u(hadamard(), 2, 3) * build_u(hadamard(), 3, 3); psi = H3 * psi; % 调用平衡函数 f(x)=x1 XOR x2,通过相位反冲 % 这里直接用 QCF 提供的 f_b0 对应的 oracle 矩阵 Uf = dj_b0(2); % 传入输入比特数,返回 8x8 oracle 矩阵 psi = Uf * psi; % 对输入寄存器 q1, q2 再施加 H H2_on_inputs = build_u(hadamard(), 1, 3) * build_u(hadamard(), 2, 3); psi = H2_on_inputs * psi; psi = renormalise(psi); % 测量输入寄存器 probs = measure_subspace(psi, [1 2]); disp(probs); % 如果 probs(1) ~ 1 则是常数,否则平衡这段代码里dj_b0(2)是我从 QCF 源码里读到的用法——它直接根据输入比特数生成对应的 oracle 矩阵,省去自己写相位反冲的循环。measure_subspace.m是 QCF 里专门做子空间测量的函数,第二个参数[1 2]表示只测量第 1、2 个量子比特,返回的是各基态的累计概率向量。测量结果如果probs(1)接近 1,说明两个输入比特都塌缩到 0,函数为常数;如果能量分散在多个基态,就是平衡函数。
3.3 参数表:QCF 提供的 Deutsch-Jozsa 相关函数
下表列出了我在使用过程中常用到的 Deutsch-Jozsa 相关函数,方便你在项目里只挑需要的模块。
| 函数名 | 输入 | 输出 | 作用 |
|---|---|---|---|
deutsch_jozsa.m | 无(或按内部参数) | 分类结果 | 完整算法封装 |
dj_c0.m/dj_c1.m | 总比特数 | 电路矩阵 | 常数函数 0 / 1 的 oracle |
dj_b0.m/dj_b1.m | 输入比特数 | 电路矩阵 | 平衡函数 oracle,对应 x1 XOR x2 及其补 |
measure_subspace.m | 态向量、比特列表 | 概率向量 | 测量指定比特组合 |
f.m | 状态向量 | fval, phase | 通用函数调用接口 |
3.4 常见坑:顺序别写反,H 门的重数要对应
第一个坑是 bit 顺序。QCF 的向量下标和教科书相反的情况不是没有,我遇到过把bin2vec('10')当成q1=1,q2=0,但实际模拟器里它表示q1=0,q2=1的情况。最稳妥的方法是先disp(bin2vec('10'))看非零位置在第几个元素,再对应到build_u的索引。第二个坑是忘了对工作比特也施加 H 门。Deutsch-Jozsa 的初始态要求辅助比特是|1>,很多初学版本只对输入寄存器做叠加,结果相位反冲之后根本无法产生正确的干涉。第三个坑是dj_b0这类函数的输入参数在不同版本里可能是函数句柄而不是比特数,跑之前先用help dj_b0确认签名。我自己在 Octave 5.2 下跑通所有版本,Matlab R2021a 也没问题,但 32 位系统下大矩阵会提示内存不足,建议设置SetUp里的偏好为双精度。
4. Grover 搜索与量子傅里叶变换:在 QCF 里实现两类关键算法
4.1 Grover 算法的迭代结构:从grover.m看振幅放大
Grover 搜索算法要解决的是无序数据库搜索问题。在 N=2^n 个条目中找到目标项,经典需要在平均 N/2 次查询,Grover 通过振幅放大把查询次数压到 O(sqrt(N))。QCF 的grover.m实现了完整的迭代,但它没有把迭代次数写死,而是留给了调用方。这是和 Deutsch-Jozsa 最大的不同:Grover 需要你根据 N 和目标数提前计算最优迭代次数。
% 使用 QCF 的 grover 函数,N=8(3比特),设目标态为 |101> N = 8; iterations = floor(pi/4 * sqrt(N)); % Grover 最优迭代次数近似公式 psi = grover(3, iterations, '101'); psi = renormalise(psi); probs = measure_subspace(psi, [1 2 3]); [max_prob, idx] = max(probs); fprintf('最大概率 %.4f,对应基态下标 %d\n', max_prob, idx-1);这段代码里的grover(3, iterations, '101')是 QCF 的签名:第一个参数是比特数,第二个是迭代次数,第三个是标记的目标二进制串。如果不传入第三个参数,它会随机标记一个目标,这在你只想验证算法行为时也很有用。注意我这里用了floor(pi/4 * sqrt(N)),这个公式只适用于单目标搜索,多目标时需要除以目标数的平方根。
grover.m内部做的事情是标准的振幅放大:先构造均匀叠加态,然后循环执行“Oracle(相位翻转目标态)→ 扩散算子(翻转所有概率幅关于平均值)”。在 QCF 的实现里,扩散算子通常用build_u和 H 门组装。你不需要复现每一步矩阵乘法,但要理解迭代次数为什么必须精确:迭代太少,目标振幅还没放大到最大;迭代太多,振幅会越过峰值落回去。
4.2 手工拆解 Grover 的单次迭代
为了让你看清grover.m的背后,这里给出一种手工实现单次迭代的写法,适合自己插入探测点:
% 3 比特,目标 |101> n = 3; psi = renormalise( ones(2^n,1) ); % 均匀叠加 targetIdx = 5; % |101> -> 下标 5 % Oracle:目标相位翻转 H3 = build_u(hadamard(), 1, 3) * build_u(hadamard(), 2, 3) * build_u(hadamard(), 3, 3); Uf = eye(8); Uf(targetIdx+1, targetIdx+1) = -1; psi = Uf * psi; % 扩散算子 D = H U0 H,U0 是对 |00..0> 的相位翻转 U0 = -eye(8); U0(1,1) = 1; % 翻转除 |0> 外的所有相位 D = H3' * U0 * H3; psi = D * psi; psi = renormalise(psi); disp(psi')这里U0的构造用的是标准的 Grover 扩散矩阵,H3'等于H3因为是实对称矩阵。跑完一次迭代后你会发现目标下标 5 的振幅从 1/sqrt(8)≈0.3536 增大到了 0.8839,其他振幅缩小到 0.1768。这个数值变化是振幅放大的核心。
4.3 QFT 与qft.m的相位估计基础
量子傅里叶变换(QFT)是 Shor 算法和相位估计的子模块。QCF 的qft.m实现了标准 QFT,但它不是直接调用 FFT,而是按教科书方式用 H 门和控制相位旋转门搭出来的。这样做的意义在于你能直接观察到相位编码过程。
% 对 4 比特量子态施加 QFT n = 4; psi = bin2vec('1010'); psi_qft = qft(psi, n); % 对比手动 FFT 的结果,注意顺序 psi_fft = fft(psi) / sqrt(2^n);qft.m的第二个参数n用来指示你可以忽略状态向量长度而只指定参与变换的比特数。这个函数的返回结果和 MATLAB 的fft有三个重要差异:一是整体有 1/sqrt(N) 的归一化因子,二是比特顺序的颠倒,三是相位方向反向。所以如果你直接用fft做对照,需要先circshift再取共轭。这也不难理解,量子 QFT 的定义里根号分母和经典 FFT 归一化不同,且 QCF 按量子电路逐级写,导致输出顺序等于输入顺序的“比特反转”,这在相位估计算法里是后续处理要修正的。
4.4 Grover 与 QFT 结合:分数阶搜索
QCF 包里还有cf_approx.m和cf_assert.m,这组函数是关于连续函数近似和断言验证的。它们可以在 Grover 的一个变体里使用——当目标函数不是一个离散搜索项,而是一个需要判断属性的黑盒函数时,cf_approx可以构造近似 oracle 和对应的振幅放大策略。这个场景在实际工程里比纯数据库搜索更常见:比如在参数优化里,你面对的是连续空间,目标态并不是某个固定二进制串,而是某个区间内的点。QCF 封装的思路是用离散网格来近似连续函数,再在这个近似的 oracle 上跑 Grover,最后用cf_assert验证近似的误差界。这部分代码相比于前面两个算法,更接近研究工具而不是教学代码,建议在使用前先读一下README.md中对连续函数采样的说明。
5. 测量、可视化与验证:QCF 排错与结果确认的实用技巧
5.1 用iplot.m和qimage.m看状态到底长什么样
量子模拟里最痛苦的事情不是算法写不对,而是写完了不知道中间态是不是你想要的。iplot.m是 QCF 的状态可视化函数,直接画概率幅分布图。
% 画 4 比特叠加态的实部与虚部 pdf = full(psi .* conj(psi)); iplot(pdf); % 对于密度矩阵可视化 rho = psi * psi'; qimage(rho);iplot默认绘制的是各基态的概率分布柱状图,qimage显示密度矩阵的模值图。调试时我习惯在每一个关键门之后调用iplot(full(abs(psi).^2)),观察能量是否还集中在预期基态上。需要留意的是,这两个函数在 Octave 里的绘图后端可能不支持某些特性,如果报字库错误,换成plot(0:2^n-1, abs(psi).^2)就行。
5.2measure.m的采样语义与measure_subspace的关系
measure.m是单次测量函数,返回一个基态下标或与该下标对应的态向量,而前面用到的measure_subspace是在指定比特子系统上做概率累加。这两个函数不能混用:measure做一次随机采样,结果是非确定的;measure_subspace返回概率分布,是确定性的。验证算法正确性时应该用后者。
5.3 验证 Grover 收敛:迭代次数的实际观察
常见错误是把 Grover 迭代次数设成固定值。下面这张表是 3 比特状态下,不同迭代次数对应的目标态概率,来自我本机 Octave 运行grover.m的结果:
| 迭代次数 k | 目标态概率 | 说明 |
|---|---|---|
| 0 | 0.1250 | 初始均匀分布 |
| 1 | 0.7813 | 接近最优 |
| 2 | 0.9453 | 最优(理论值为 sin^2(5π/16)≈0.945) |
| 3 | 0.9453 | 开始越过峰值 |
| 4 | 0.7813 | 明显过冲 |
| 5 | 0.1250 | 完全失效 |
你可以在 QCF 里写一个循环,对不同 k 调用grover,再用measure_subspace提取目标概率。这个验证方法对理解振幅放大机制非常有帮助,从中你会直观看到为什么量子算法需要精确控制迭代步数,不是越多越好。
5.4 一个收敛的调试套路:断言+归一化+子空间测量
最后分享一个我调试 QCF 脚本时固定使用的三步走套路。第一步,在每个门之后调用renormalise,并断言abs(norm(psi)-1) < 1e-10,保证后续概率计算是有效的。第二步,用pretty(psi)打印关键中间态的符号形式,对比教科书中的推导结果。第三步,在最终测量前用measure_subspace得到概率分布,而不是用measure做单次采样,否则你永远无法判断是算法错误还是随机性导致的结果偏差。
assertError = abs(norm(psi)-1); assert(assertError < 1e-10, 'State norm drifted: %e', assertError);这套验证思路对 QCF 里所有算法都通用,无论是deutsch_jozsa.m还是grover.m,你都能以“概率分布是否符合理论期望”作为唯一正确性标准。
本文还有配套的精品资源,点击获取