简介压缩包内含量子计算算法的MATLAB实现示例面向量子计算初学者、课程设计学生或希望快速上手量子算法验证的研究者以5个轻量文件集中演示叠加态制备、量子门操作、测量与Grover搜索算法等核心内容。资源共5个文件包含4个.m脚本spadamard.m、measure_qubit.m、grover.m、cgate.m和1个README.md说明文档整体仅5KB结构清晰、便于直接阅读和运行。目前已有1870人学习下载适合快速入门量子计算编程实践。通过grover.m可理解无序数据库搜索的算法逻辑measure_qubit.m与cgate.m帮助掌握量子测量与受控非门操作spadamard.m演示Hadamard门构造叠加态的过程配合README运行说明读者无需真实量子硬件即可在MATLAB中观察量子态演化加深对量子并行与算法原理的认识为后续学习QFT、Shor等更复杂算法打下基础。1. 量子计算在MATLAB里“能跑”的真相不加速但要展开有一个反直觉的结论用MATLAB仿真量子算法本质上并不会获得任何“量子加速”反而要把指数级的状态空间老老实实地展开成矩阵和向量。量子计算之所以在MATLAB里可行是因为量子比特的叠加、干涉、纠缠在数学上全部坍缩成复数向量、酉矩阵和张量积运算——这几样恰好是MATLAB的看家本领。搞懂“量子算法的MATLAB实现”不是去学一套神秘硬件接口而是搞懂怎么把一个量子协议翻译成可运行、可验证的线性代数程序。这才是值得你花时间的地方它把量子算法从黑匣子变成了一段你能单步调试的代码也让“为什么能work”有了最直接的检验环境。2. 态矢量与酉演化MATLAB仿真量子算法的底层数学引擎2.1 量子比特的复数向量从布洛赫球到态矢量量子比特和经典比特最大的差别不是“0和1都能是”而是它活在复数概率幅上。一个量子比特的状态写作ψ α|0⟩ β|1⟩其中 α、β 是复数满足 |α|² |β|² 1测到 0 的概率是 |α|²测到 1 的概率是 |β|²。这个“概率幅是复数”的设定非常关键复数相位会带来干涉干涉才是量子算法比经典算法“多出来”的东西。布洛赫球只是把 α、β 的比值画在一个球面上方便你直觉地想象旋转但真到算的时候布洛赫球球面的几何直觉帮不上忙你手里真正拿得住的还是一组复数坐标。在MATLAB里一个量子比特就是一个 2×1 的复数列向量仅此而已qubit0 [1; 0]; % |0⟩ qubit1 [0; 1]; % |1⟩ psi [1/sqrt(2); 1i/sqrt(2)]; % 叠加态 (|0⟩ i|1⟩)/sqrt(2) % 验证归一化概率幅度平方和为1 prob_sum sum(abs(psi).^2); % 必须等于1代码里的1i是MATLAB的虚数单位别写成i后者可能被你之前循环里的变量覆盖掉这是第一处容易翻车的地方。归一化条件不是可选项是所有后续运算的不变量任何合法操作都必须让概率和保持为 1。如果你在仿真中看到某个门作用之后概率和变成了 1.2那说明你的门矩阵写错了它不是一个合法的量子门。当有 n 个量子比特时态空间不再按“每个比特各存一份”的方式扩展而是做张量积得到一个 2ⁿ 维的复向量。两个比特的组合态写成 ψ |q₁⟩ ⊗ |q₂⟩对应MATLAB里的kron(q1, q2)。这里的“⊗”不是可选的记号它就是态空间的乘法两个独立系统的组合自由度从 224 个“分量”变成 2×24 维空间维数相乘。所以 n 个比特的态向量长度是 2ⁿ这也是为什么经典仿真在 n 到 25 左右就会撞上内存墙。2.2 量子门就是稀疏矩阵与张量积Kron运算如何构造多比特门量子门在数学上是一个满足酉条件的矩阵U†U I其中 U† 是共轭转置。酉条件意味着操作是可逆的、概率守恒的。量子门不是任意矩阵随便写一个就能用这是仿真脚本里最常见的错误来源随手填了一个数值上不满足酉条件的矩阵后面的概率分布立刻变得不可解释。MATLAB里最常接触的三个单比特门X [0 1; 1 0]; % Pauli-X量子非门翻转 |0⟩ 和 |1⟩ I eye(2); H [1 1; 1 -1] / sqrt(2); % Hadamard门把计算基变成叠加基 % 验证酉条件 HdagH H * H; % 应为单位阵Hadamard门是几乎所有量子算法的起点它把 |0⟩ 变成 (|0⟩|1⟩)/√2等于是把一个确定的状态变成均匀叠加。在仿真里你不需要理解什么“微波脉冲”就是把向量乘一个矩阵。多比特门怎么来的答案是张量积。一个系统有两个比特时你想让 H 只作用在第 0 个比特上那么整体操作就是 H ⊗ I即第一个比特做 H第二个比特什么都不做。这个“什么都不做”也是要写进矩阵里的因为态空间是联合空间任何单个比特的操作都要提升到整个 2ⁿ 维空间中去。在MATLAB里I2 eye(2); H_on_bit0 kron(H, I2); % H作用在左比特右比特保持 H_on_bit1 kron(I2, H); % 反过来的组合受控门则是另一种构造思路。以最典型的 CNOT受控非门为例控制比特为 1 时目标比特翻转否则什么都不做。当控制比特是第 0 位、目标比特是第 1 位且约定第 0 位是二进制高位时CNOT 矩阵可以写成块对角形式CNOT_bit0_ctrl_bit1_target blkdiag(I2, X);这个矩阵左上 2×2 块是 I对应控制位为 0 的情况右下 2×2 块是 X对应控制位为 1 的情况。你看受控门本质上就是“按控制位的取值选择不同的块”在仿真里你不用去模拟什么物理机制矩阵的分块结构已经把条件逻辑写清楚了。2.3 测量与坍缩在仿真中的实现随机数与投影的等价物测量是量子力学里最“不线性”的一步。酉演化是可逆的、确定性的测量却是不可逆的、随机的。仿真代码里怎么处理这一步答案比你想的直白先算概率分布再按概率随机采样然后把态向量坍缩到测量到的那个分量上。function [measured_idx, state_after] measure_state(state) % state: 归一化态向量 probs abs(state).^2; % 各计算基的出现概率 r rand(); % 均匀随机数 cdf cumsum(probs); measured_idx find(cdf r, 1); % 按概率掷骰子 collapsed zeros(size(state)); collapsed(measured_idx) state(measured_idx) / abs(state(measured_idx)); state_after collapsed; % 坍缩后的态 end这里的关键不是代码本身而是“测量”在仿真中的语义它把连续的复数概率幅变成了一次离散的抽样。单次测量的结果没有意义你必须重复制备、测量很多次统计测量结果的频率才能反推概率分布。这也解释了学术论文里为什么所有测量结果都带“统计”二字——把测量实现成“取最大概率的那个基”是错误的仿真心智模型。3. 核心量子算法为什么能work相位、振幅与阶的三种玩法3.1 量子相位估计把本征相位写进二进制读数量子相位估计QPE是很多量子算法的“前置发动机”包括 Shor 分解、量子模拟都要调用它。解决的问题是这样给定一个酉算子 U 和它的本征态 |u⟩满足 U|u⟩ e^{2πiθ}|u⟩想求出相位 θ。θ 藏在复指数的相位里直接测量读不出来因为全局相位不可观测。QPE 的思路是绕道而行把相位“拷贝”到一个辅助寄存器上然后做逆量子傅里叶变换把相位变成二进制计数。受控 U 的幂次其实是在做一件事——当辅助寄存器处于计算基 |j⟩ 时整体产生相位 e^{2πiθj}。这个线性相位正是一个离散傅里叶变换的核。所以在仿真里QPE 的核心部分可以用一行fft直接验证theta 0.37; % 待估计的本征相位 t 4; % 相位寄存器比特数 M 2^t; % 特征态上的相位累积|j⟩ - exp(2πiθj)|j⟩ phase exp(2*pi*1i*theta*(0:M-1)).; psi ones(M,1)/sqrt(M); psi psi .* phase; % 逐元素乘积模拟受控U幂次的作用 % 逆QFT与MATLAB fft的关系逆QFT正是未归一化的DFT prob abs(fft(psi)).^2 / M; [~, k] max(prob); estimated_theta k / M; % 读出相位注意这里我用fft而不是ifftMATLAB 的ifft自带 1/M 归一化因子而量子逆傅里叶变换的归一化是 1/√M所以直接做一次ifft会多出一个 1/√M 的幅度偏差用fft再手工除以 √M才能让概率分布的峰值位置准确对应相位。这个细节做仿真时容易踩但对理解“QPE 为什么是傅里叶变换”反而有帮助量子傅里叶变换和经典 FFT 本质上是同一个线性变换只是基底的排列和归一化约定不同。测量结果 k/M 就是 θ 的 t 位二进制近似。误差来自离散化相位寄存器只有 2ᵗ 个格子θ 恰好在格子缝里时能量会泄漏到相邻的计数上。这就是分辨率与比特数的直接权衡。3.2 Grover搜索振幅放大为什么只需 O(√N) 次迭代Grover 算法解决的是无序搜索问题N 个条目里标记一个目标经典盲搜平均要试 N/2 次Grover 只要约 √N 次。它不是“并行遍历”而是用一个非常漂亮的几何过程——振幅放大。我一般这么理解把所有状态的均匀叠加记为 |s⟩目标态记为 |t⟩。Oracle 做的事是给目标态的振幅加一个负号即把 |s⟩ 中对应目标的那个分量翻转符号。这等于把整个态向量沿着“非目标子空间”做了一次反射。之后扩散算子 D 2|s⟩⟨s| − I 再把向量对 |s⟩ 方向反射一次。两次反射的复合效果就是让 |s⟩ 向 |t⟩ 方向旋转了 2θ 的角度其中 θ 是 |s⟩ 与 |t⟩ 的夹角。反复迭代每次转 2θ直到转到目标态附近。在MATLAB里这个几何过程可以非常直观地写成矩阵运算不需要按比特逐门展开——因为 Grover 迭代本身是一个全局操作n 4; N 2^n; target 13; % 要搜索的目标索引 % 均匀叠加态 psi ones(N,1) / sqrt(N); % Oracle: 对目标索引的振幅取负 oracle eye(N); oracle(target1, target1) -1; % MATLAB索引从1开始 % 扩散算子: D 2|s⟩⟨s| - I s ones(N,1) / sqrt(N); diffusion 2 * (s * s) - eye(N); % 迭代约 pi/4*sqrt(N) 次 r round(pi/4 * sqrt(N)); for k 1:r psi diffusion * (oracle * psi); end probs abs(psi).^2; [~, idx] max(probs); fprintf(估计目标索引: %d\n, idx - 1);这段代码最值得琢磨的是 Oracle 的写法。在这里我图方便把它写成了对角矩阵这在仿真里完全能跑通。但你要知道真实的 Grover 算法里 Oracle 是一个黑盒你并不知道目标是什么如果为了验证“搜索”效果而把 Oracle 写成对角矩阵那只是把答案告诉了程序。一个更诚实的仿真方式是让 Oracle 是一个函数内部把索引经过某种哈希或计算判断是否为目标。这样做才能让你体会到 Grover 的复杂度优势来自哪里迭代次数只与 √N 有关与目标的具体位置无关。3.3 Shor分解最出名的算法为什么最难仿真Shor 分解之所以出名是因为它理论上能攻破 RSA把大整数分解成质因数。它的核心步骤不是“直接找因数”而是把分解问题归结为求阶问题——找一个数 a求满足 aʳ ≡ 1 (mod N) 的最小整数 r。一旦得到 r就能用经典数论方法以高概率算出 N 的因子。求 r 这件事用的是量子相位估计。在仿真里复现 Shor 分解要冷静完整的 Shor 电路需要约 O(n³) 个门n 是被分解数的比特数。分解一个 15 位二进制4比特的数很简单但要分解一个 2048 位的 RSA 密钥对应的量子电路规模远超任何经典计算机的存储容量。这不是 MATLAB 不够好而是态矢量仿真本身面临指数墙量子寄存器 n 个比特仿真就要存 2ⁿ 个复数。所以 MATLAB 实现 Shor 通常有两种形态一是小规模整数分解的完整电路仿真二是把“阶估计”用算法级方式模拟跳过逐门电路直接构造相位估计的输入。理解 Shor 的关键不在那一堆门而在“阶估计为什么是相位估计”定义算子 U_a|x⟩ |ax mod N⟩它的本征态携带的相位恰好包含 a 的阶 r 的信息。相位估计测出的相位是 s/r 的形式其中 s 是某个整数对多个测量结果做连分数展开就能恢复出 r。整个过程里量子部分只负责“高效估计相位”经典部分负责从相位反推 r、再反推因子。这是一个很典型的“量子测量 经典后处理”分工也是你在仿真实验里最容易拆开的解耦点先用经典方法算好相位再单独验证后处理算法。4. 源码级实现细节决定脚本正确性与性能的关键约定4.1 索引映射比特顺序与二进制编码的约定量子仿真里的第一个约定问题是态向量下标 0..2ⁿ-1 对应的二进制串到底哪个比特是高位这个问题不统一后面写受控门、看测量结果时就会对不上。我建议一上来就用一个函数把全工程的比特顺序锁死function idx state_index(bits) % bits: 行向量每个元素是0或1bits(1)对应第0个量子比特 % 约定bits(1)是二进制最高位 weights 2.^(numel(bits)-1 : -1 : 0); idx bits * weights 1; % MATLAB索引从1起 end这个约定的含义是态向量第 idx 个分量对应比特串 bits而 bits(1) 是最高位。比如 n2 时索引 0 对应 |00⟩索引 1 对应 |01⟩索引 2 对应 |10⟩索引 3 对应 |11⟩。为什么强调这个因为kron的拼接顺序、受控门的块排列、测量结果的反推全部依赖它。你可以在脚本开头写一行注释甚至写个断言函数在所有入口检查比特串长度。我见过不少仿真脚本跑出来的概率分布完全正确、但解读结果时全反了就是因为测量后把二进制串的低位当成了高位。4.2 张量积顺序陷阱为什么同样的矩阵作用位置完全相反这是仿真脚本里翻车率最高的一个点kron(A, B)中 A 作用在左比特、B 作用在右比特这是固定的但很多人把“第0个比特”和“矩阵里第0位”搞混导致门作用错了目标。看这个对比psi kron(qubit0, qubit1); % |01⟩左比特是0右比特是1 % 想让H作用在右比特第1位上 psi1 kron(I2, H) * psi; % 正确右比特被H变换 % 错误写法kron(H, I2)会把H作用在左比特上 psi_wrong kron(H, I2) * psi;更隐蔽的是受控门的构造。同样是 CNOT控制比特在左还是右矩阵长相完全不同% 控制bit0目标bit1块对角形式 CNOT_ctrl0 blkdiag(I2, X); % 控制bit1目标bit0要按基底的排列重排 CNOT_ctrl1 [1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0];第二个矩阵不是简单地再套一个blkdiag因为二进制高位在左的约定下控制位变成低位之后块结构的位置也跟着变了。我的习惯是所有受控门都写成函数内部用state_index把控制位和目标的比特位置显式算出来再重排矩阵绝不手写这个排列。这种问题单步调试看不出来因为矩阵乘法不报错只有把测量结果的分布和预期对比时才暴露。4.3 全矩阵 vs 稀疏矩阵比特数的内存预算是硬约束态矢量仿真的天花板是内存不是时间。一个 n 比特的态向量要存 2ⁿ 个复数MATLAB 默认用 double 双精度一个复数占 16 字节。内存估算非常直白量子比特数态向量大小复数内存占用是否适合笔记本仿真10102416 KB轻松16655361 MB轻松20约 100 万16 MB轻松24约 1677 万268 MB勉强28约 2.68 亿4 GB困难30约 10.7 亿16 GB基本不可行注意到 24 比特是一个分水岭超过 24 比特光存一个态向量就吃掉几百 MB而一次 Grover 迭代还要创建 Oracle 和扩散矩阵内存瞬间翻几倍。此时可以用稀疏矩阵来存门n 20; Gate kron(X, speye(2^(n-1))); % 单比特门作用在第0位sparse存储稀疏矩阵对单比特门很友好因为它作用在 n 比特空间时非零元素只有 2ⁿ 的常数倍远少于全矩阵的 2²ⁿ。但要注意稀疏矩阵只是让“门”的存储变小了态向量本身永远是满的稠密向量而且矩阵乘向量之后的结果也是稠密的。所以 sparse 可以推迟内存墙的到来但救不了本质上的指数复杂度。真要做到 30 比特以上需要换张量网络态MPS之类的压缩表示那完全是另一套仿真范式。4.4 可运行的Grover核心骨架矩阵折叠与逐门展开写量子算法仿真时有两种风格一种是“逐比特展开电路”把每个门按比特位置kron到整个空间再逐个相乘另一种是“矩阵折叠”把整个算法迭代过程抽象成少数几个矩阵直接乘。前者更贴近真实硬件后者更贴近数学本质也更容易验证正确性。Grover 的骨架适合用矩阵折叠来写因为它的迭代操作非常规整function [est_idx, probs] grover_sim(n, target) N 2^n; psi ones(N,1) / sqrt(N); % 第1步Oracle目标振幅取负 O eye(N); O(target1, target1) -1; % 第2步扩散矩阵 s ones(N,1) / sqrt(N); D 2*s*s - eye(N); % 第3步迭代 r ≈ π/4·√N 次 r round(pi/4 * sqrt(N)); for k 1:r psi D * (O * psi); end probs abs(psi).^2; [~, est_idx] max(probs); est_idx est_idx - 1; % 转回从0开始的索引 end这段代码值得注意的机制点是Oracle 和扩散矩阵都是对全空间定义的所以每次迭代是两次全空间矩阵乘法。r 的取值不是越多越好过了 π/4·√N 这个峰值之后振幅会被“放大过头”目标概率反而下降。这就是 Grover 算法里最反直觉的地方——它不是单调收敛而是一个周期性的摆动。你可以把 r 从 1 一直跑到 2 倍建议值画一条成功概率曲线峰值位置会验证这个 π/4·√N 的公式。这种“先看到整体行为再理解单个步骤”的路径是仿真相对于真实硬件最大的优势。5. 避坑记录量子算法MATLAB仿真中最常见的五类翻车现场5.1 概率和不是1门矩阵写错了还是归一化丢了现象跑完一段演化后sum(abs(psi).^2)输出 1.3 或 0.7甚至出现复数。原因两种情况。一是门矩阵写得不满足酉条件比如把 Hadamard 门漏了 1/√2 的归一化因子二是某一步做矩阵乘法时误用了普通转置而不是共轭转置操作MATLAB 中A是共轭转置A.才是普通转置。对实矩阵两者一样一旦门含虚数相位用错就全盘皆输。解决在演化循环里每步断言概率和误差小于 1e-10写门矩阵时先用一行norm(U*U - eye(size(U)))验证酉条件。否则调试几天可能只是在追一个早已存在的符号错误。5.2 受控门结果和预期全反比特序约定没锁死现象CNOT 之后测量分布显示翻转的比特是控制位而不是目标位但矩阵明明是从教科书抄的。原因教科书里的 CNOT 矩阵默认控制位是高位你的态向量构建时却把高位赋给了另一个比特。kron本身不背锅锅在“高位”这个约定在不同参考系里指代不同。解决用state_index函数把控制位、目标位的比特位置显式算出来再通过行/列置换重排矩阵。凡是和“第几位比特”相关的门都走函数生成不手写矩阵。5.3 expm直接对高维哈密顿量求矩阵指数内存直接崩现象n20 时想用expm(-1i * H * t)模拟系统演化MATLAB 直接报内存不足或卡到天荒地老。原因expm 会把哈密顿量矩阵做特征分解或 Padé 近似这在 2ⁿ 维空间上是 O(2³ⁿ) 级别的高耗时操作而且中间临时矩阵可能比态向量大一个数量级。量子系统哈密顿量通常是稀疏的但矩阵指数稠密。解决不要对整个哈密顿量求指数。如果哈密顿量能分解成多项之和且各项不对易用 Trotter 分解把演化拆成多个稀疏门连乘如果只是验证算法逻辑则直接构造离散的酉门别走微分方程这条路。仿真量子算法时你应该把“矩阵指数”这个东西忘掉它不属于这个工具箱。5.4 测量结果每次跑都不一样以为算法不对其实是随机性现象Grover 仿真每次运行估计出的目标索引时对时错或者相位估计算出不同结果让人误以为代码有 bug。原因测量实现里用了rand()单次抽样天然有方差。量子测量本来就是一个随机过程正确的验证方式是多次重复取频率而不是看单次样本。解决在调试阶段用rng(固定种子)固定随机数保证可复现在验证阶段跑 1000 次测量画概率分布直方图看峰值是否落在目标位置。把“单次测量的结果”和“概率分布”分开看是量子仿真调试的重要心法。5.5 Oracle写成对角矩阵跑通了但什么都没证明现象Grover 仿真成功率高得惊人迭代几次就 100% 找到目标但你换一个目标索引成功率曲线完全不匹配理论值。原因把 Oracle 直接写成eye(N)然后对角改一个 -1等于把目标位置直接编码进了矩阵。Grover 的黑盒假设是“不知道目标在哪”这种实现把答案写进了程序复杂度分析自然失效。更隐蔽的是这个对角 Oracle 在 Grover 门之外还给整个态引入了相对相位影响扩散步骤的效果。解决把 Oracle 实现成函数内部通过计算或查找判断输入索引是否等于目标返回 ±1 的相位因子或者用一个只含经典逻辑的电路来模拟 Oracle。仿真是否正确判断标准是成功率随迭代次数的曲线是否贴合 π/4·√N 的规律而不是“最终是否找到了”。6. 设计取舍什么时候不该用MATLAB仿真量子算法MATLAB 做量子算法仿真优势不在规模在连通性。态矢量表示和线性代数教学天然无缝衔接量子门就是矩阵、测量就是抽样这对验证算法结构和教学演示极有价值。但如果你的目标不是验证思想而是逼近真实硬件或者想跑 30 比特以上MATLAB 的态矢量仿真就是错误工具。取舍的分界线大致可以画在 24 比特。24 比特以内的态矢量仿真普通笔记本能跑MATLAB 的开发效率和调试体验比任何框架都好超过 24 比特内存就开始成为主约束此时稀疏门只能延缓而无法解决指数膨胀。如果要做噪声模型、脉冲级模拟或用真实量子硬件后端跑实验那生态更丰富的专业量子计算框架更合适。如果你非要继续用 MATLAB 做大比特数仿真出路是张量网络态MPS把态向量压缩成一组低秩张量的乘积比特数到 50 以上也能近似仿真但它需要你放弃“所有概率幅都在内存里”的直观模型转而接受一套更抽象的压缩表示。用一个简单的判断表收束这个思考教学验证、算法结构调试、小规模数值实验选 MATLAB资源受限环境下的规模化仿真优先考虑稀疏矩阵加 MPS需要对接硬件的实验则应该选择支持后端切换的专业框架。这个选择不是谁好谁坏的问题而是你要理解每种工具的代价——MATLAB 的代价是显式矩阵乘法的 O(2ⁿ) 内存但它换来了最容易理解的调试路径。我自己的经历是一个教训有一次想验证某个变体 Grover 算法在 30 比特下的行为直接在 MATLAB 里初始化全态向量机器立刻内存耗尽我花了一晚上以为是自己代码哪里泄漏了。后来才意识到这是一个复杂度问题不是一个 bug 问题。从那以后我每次写量子仿真脚本都强制走一遍三件事先写清楚比特序约定再估算内存上限最后想清楚要验证的是算法数学还是物理过程。希望你也能用这三步把量子计算从黑匣子拆成可控的实验。希望帮到你。本文还有配套的精品资源点击获取