模拟

本教程涵盖 pyqpanda3 中所有可用的量子模拟器(simulator)，解释何时以及如何使用每一种模拟器。

---

概述

pyqpanda3 提供多个模拟器后端(simulator backend)，每种针对不同的使用场景进行了优化：

模拟器	态表示	噪声支持	最佳用途
CPUQVM	态矢量(state vector)	是（基于采样）	通用模拟
GPUQVM	态矢量（GPU）	是	带 CUDA 的大型线路
DensityMatrixSimulator	密度矩阵(density matrix)	是（精确）	噪声分析、态层析
Stabilizer	稳定子态(stabilizer state)	是	Clifford 线路、快速基准测试
PartialAmplitudeQVM	部分态矢量	否	大型线路、特定振幅

---

CPUQVM

CPUQVM 是主要的态矢量模拟器(statevector simulator)。它通过维护量子系统的完整态矢量并直接应用酉操作来模拟量子线路。

基本用法

创建一个 CPUQVM 实例，并使用指定的测量采样次数(shots)运行量子程序：

from pyqpanda3 import core

# 构建贝尔态线路
prog = core.QProg()
prog << core.H(0)
prog << core.CNOT(0, 1)
prog << core.measure([0, 1], [0, 1])

# 创建模拟器并运行
machine = core.CPUQVM()
machine.run(prog, shots=10000)

# 获取测量结果
result = machine.result()
counts = result.get_counts()
print(f"Measurement counts: {counts}")
# 预期输出: {'00': ~5000, '11': ~5000}

获取结果

QResult 对象提供了多种提取模拟数据的方法：

# 以字典形式获取测量计数
counts = result.get_counts()
# 返回: {'00': 5023, '11': 4977}

# 获取特定量子比特的概率
prob_list = result.get_prob_list([0, 1])
# 返回: [0.5023, 0.0, 0.0, 0.4977]

# 以字典形式获取概率
prob_dict = result.get_prob_dict([0, 1])
# 返回: {'00': 0.5023, '01': 0.0, '10': 0.0, '11': 0.4977}

# 获取采样次数
num_shots = result.shots()
print(f"Shots: {num_shots}")

# 打印格式化结果
result.print_results()

获取态矢量

对于没有测量的线路，态矢量直接表示量子态：

# 构建不含测量的线路
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)

machine = core.CPUQVM()
machine.run(prog, shots=1)

# 获取态矢量
sv = machine.result().get_state_vector()
print(f"State vector: {sv}")
# 预期: [1/sqrt(2), 0, 0, 1/sqrt(2)]

$n$ 个量子比特的态矢量 $|\psi ⟩$ 有 $2^n$ 个复数振幅(amplitude)：

$$|\psi ⟩=\sum _{i=0}^{2^n-1}{\alpha }_i|i⟩$$

其中 $|i⟩$ 表示计算基态(computational basis state)，${\alpha }_i$ 是满足以下条件的复数振幅：

$$\sum _{i=0}^{2^n-1}|{\alpha }_i{|}^2=1$$

带噪声运行

CPUQVM 通过 NoiseModel 参数支持噪声模拟(noise simulation)：

from pyqpanda3 import core

# 创建噪声模型
noise = core.NoiseModel()
error = core.depolarizing_error(0.01)
noise.add_all_qubit_quantum_error(error, core.GateType.CNOT)

# 构建并运行含噪声线路
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
prog << core.measure([0, 1], [0, 1])

machine = core.CPUQVM()
machine.run(prog, shots=10000, model=noise)
print(machine.result().get_counts())
# 有噪声时，可能在 '01' 和 '10' 中看到少量计数

输出：

{'00': 4950, '11': 4920, '01': 70, '10': 60}

CPUQVM 高级功能

除了基本的线路执行，CPUQVM 还提供了多项高级功能，用于提取量子态的详细信息并优化线路结构。

展平量子程序

复杂线路通常包含嵌套子线路(nested subcircuit)、控制流操作或复合门。展平(flatten)程序会将这些解析为原始门(primitive gate)的线性序列，从而减少模拟开销并使线路结构更加明确。

from pyqpanda3 import core

# 构建具有嵌套结构的线路
prog = core.QProg()
prog << core.H(0)
prog << core.H(1)
prog << core.CNOT(0, 1)
prog << core.measure([0, 1], [0, 1])

# 将程序展平为线性门序列
flat_prog = prog.flatten()

# 运行展平后的程序
machine = core.CPUQVM()
machine.run(flat_prog, shots=10000)
print(machine.result().get_counts())

展平在性能分析(profiling)或转译(transpilation)之前特别有用，因为它确保每个操作都是原始门。线路与其展平形式之间的关系为：

$$|{\psi }_{\text{flat}}⟩=U_{\text{flat}}|0^n⟩=U_k\cdots U_2{U}_1|0^n⟩$$

其中每个 $U_i$ 是按顺序应用的单个原始门。

访问完整态矢量

运行不含测量的线路时，可以获取完整的态矢量以检查各个振幅和相位。

from pyqpanda3 import core
import cmath

# 在 3 个量子比特上制备叠加态
prog = core.QProg()
prog << core.H(0) << core.H(1) << core.H(2)

machine = core.CPUQVM()
machine.run(prog, shots=1)

# 获取完整态矢量
sv = machine.result().get_state_vector()

# 检查振幅的幅值和相位
for idx, amp in enumerate(sv):
    binary = format(idx, '03b')
    magnitude = abs(amp)
    phase = cmath.phase(amp)
    print(f"|{binary}>: amplitude = {amp:.4f}, "
          f"|amp| = {magnitude:.4f}, phase = {phase:.4f}")

对于处于等概率叠加的 $n$ 个量子比特，每个振幅为：

$${\alpha }_i=\frac{1}{\sqrt{2^n}}\mathrm{\forall }i\in \{0,\dots ,2^n-1\}$$

提取概率分布

可以提取任意量子比特子集的概率分布，从而实现对多量子比特态的边缘分析(marginal analysis)。

from pyqpanda3 import core

# 构建纠缠的 3 量子比特态（GHZ 态）
prog = core.QProg()
prog << core.H(0)
prog << core.CNOT(0, 1)
prog << core.CNOT(1, 2)

machine = core.CPUQVM()
machine.run(prog, shots=1)

# 获取所有 3 个量子比特的完整概率分布
result = machine.result()
prob_dict = result.get_prob_dict([0, 1, 2])
print("Full probability distribution:")
for state, prob in sorted(prob_dict.items()):
    print(f"  P(|{state}>) = {prob:.6f}")

# 获取仅量子比特 0 和 2 的边缘概率
prob_marginal = result.get_prob_dict([0, 2])
print("\nMarginal distribution (qubits 0, 2):")
for state, prob in sorted(prob_marginal.items()):
    print(f"  P(|{state}>) = {prob:.6f}")

# 以有序列表形式获取概率
prob_list = result.get_prob_list([0, 1, 2])
print(f"\nProbability list: {[f'{p:.6f}' for p in prob_list]}")

对于 GHZ 态 $|\text{GHZ}⟩=\frac{1}{\sqrt{2}}(|000⟩+|111⟩)$，概率分布为：

$$P(x)=\{\begin{array}{ll}0.5& \text{if }x=000\text{ or }x=111\\ 0& \text{otherwise}\end{array}$$

---

GPUQVM

GPUQVM 为大型量子线路提供 CUDA 加速模拟。它使用与 CPUQVM 相同的接口，但利用 GPU 并行性进行态矢量操作。

注意：GPUQVM 仅在 pyqpanda3 以 USE_CUDA=ON 编译时可用。

用法

from pyqpanda3 import core

prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1) << core.measure([0, 1], [0, 1])

# 创建 GPU 模拟器
gpu_machine = core.GPUQVM()
gpu_machine.run(prog, shots=10000)

result = gpu_machine.result()
print(result.get_counts())

何时使用 GPUQVM

在以下情况下 GPU 模拟具有优势：

• 量子比特数超过 20（态矢量大小 > $2^{20}\times 16$ 字节 = 16 MB）
• 需要大量采样以确保统计精度
• 线路深度适中但门数量较多

$n$ 个量子比特的内存需求为 $2^n\times 16$ 字节（复数双精度）。例如：

量子比特数	态矢量大小	内存
20	$2^{20}$	16 MB
25	$2^{25}$	512 MB
30	$2^{30}$	16 GB

---

DensityMatrixSimulator

DensityMatrixSimulator 维护量子系统的完整密度矩阵 $\rho$，支持精确噪声模拟和混合态分析。

什么是密度矩阵？

密度矩阵(density matrix)推广了态矢量，可以描述纯态(pure state)和混合态(mixed state)：

$$\rho =\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|$$

对于纯态 $|\psi ⟩$，密度矩阵为：

$$\rho =|\psi ⟩⟨\psi |$$

基本用法

from pyqpanda3 import core

# 创建密度矩阵模拟器
dm_sim = core.DensityMatrixSimulator()

# 构建线路（不需要测量）
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)

# 运行模拟
dm_sim.run(prog)

# 获取密度矩阵
dm = dm_sim.density_matrix()
print(f"Density matrix shape: {dm.shape}")
print(f"Density matrix:\n{dm}")

概率

DensityMatrixSimulator 提供多种计算概率的方法：

# 按索引获取特定态的概率
prob_0 = dm_sim.state_prob(0)
print(f"P(|00⟩) = {prob_0}")

# 按二进制字符串获取概率
prob_01 = dm_sim.state_prob("01")
print(f"P(|01⟩) = {prob_01}")

# 所有概率
all_probs = dm_sim.state_probs()
print(f"All probabilities: {all_probs}")

# 特定量子比特的概率
qubit_probs = dm_sim.state_probs([0, 1])
print(f"Qubit [0,1] probabilities: {qubit_probs}")

测量态 $|i⟩$ 的概率为：

$$P(i)=⟨i|\rho |i⟩={\rho }_{ii}$$

约化密度矩阵

偏迹(partial trace)操作计算量子比特子集的约化密度矩阵(reduced density matrix)：

# 获取量子比特 0 的约化密度矩阵
reduced = dm_sim.reduced_density_matrix([0])
print(f"Reduced density matrix for qubit 0:\n{reduced}")

对二分系统(bipartite system)的量子比特 $B$ 的偏迹：

$${\rho }_A={\text{Tr}}_B({\rho }_{AB})=\sum _j⟨j_B|{\rho }_{AB}|j_B⟩$$

噪声模拟

DensityMatrixSimulator 在精确噪声模拟方面表现优异：

from pyqpanda3 import core

# 创建噪声模型
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.amplitude_damping_error(0.05),
    core.GateType.H
)
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.01),
    core.GateType.CNOT
)

# 带噪声运行
dm_sim = core.DensityMatrixSimulator()
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog, noise)

# 检查生成的混合态
dm = dm_sim.density_matrix()
print(f"Trace: {dm_sim.density_matrix().trace()}")
print(f"State probabilities: {dm_sim.state_probs()}")

引入噪声后，演化遵循：

$${\rho }^{\prime }=\sum _i{K}_{i}U\rho U^{†}{K}_i^{†}$$

其中 $K_i$ 是噪声信道的 Kraus 算子(Kraus operator)。

密度矩阵操作

DensityMatrixSimulator 提供了强大的工具，用于分析超越简单概率提取的量子态。这些操作对于量子信息分析、纠缠量化和噪声表征至关重要。

偏迹和约化密度矩阵

偏迹操作通过对剩余量子比特求迹（边缘化）来关注子系统。这是经典联合概率分布边缘化的量子类似。

from pyqpanda3 import core
import numpy as np

# 创建最大纠缠贝尔态
dm_sim = core.DensityMatrixSimulator()
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog)

# 完整的 2 量子比特密度矩阵
dm_full = dm_sim.density_matrix()
print(f"Full density matrix shape: {dm_full.shape}")
# 输出: (4, 4)

# 仅量子比特 0 的约化密度矩阵
reduced_q0 = dm_sim.reduced_density_matrix([0])
print(f"Reduced density matrix for qubit 0:\n{reduced_q0}")
# 对于贝尔态，每个单量子比特约化态都是最大混合态

# 仅量子比特 1 的约化密度矩阵
reduced_q1 = dm_sim.reduced_density_matrix([1])
print(f"Reduced density matrix for qubit 1:\n{reduced_q1}")

对于贝尔态 $|{\mathrm{\Phi }}^{+}⟩=\frac{1}{\sqrt{2}}(|00⟩+|11⟩)$，完整密度矩阵为：

$${\rho }_{AB}=|{\mathrm{\Phi }}^{+}⟩⟨{\mathrm{\Phi }}^{+}|=\frac{1}{2}\left(\begin{array}{cccc}1& 0& 0& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 1\end{array}\right)$$

对量子比特 B 求迹得到量子比特 A 的最大混合态(maximally mixed state)：

$${\rho }_A={\text{Tr}}_B({\rho }_{AB})=\frac{I}{2}=\frac{1}{2}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$

纯度计算

纯度(purity)量化量子态的"混合"程度。纯态的纯度为 1，而 $d$ 维最大混合态的纯度为 $1/d$。

from pyqpanda3 import core
import numpy as np

dm_sim = core.DensityMatrixSimulator()

# 情况 1：纯态（无噪声）
prog_pure = core.QProg()
prog_pure << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog_pure)
dm_pure = dm_sim.density_matrix()

# 纯度为 Tr(rho^2)
purity_pure = np.trace(dm_pure @ dm_pure).real
print(f"Pure state purity: {purity_pure:.6f}")
# 输出: 1.000000

# 情况 2：混合态（含去极化噪声）
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.1),
    core.GateType.CNOT
)
dm_sim2 = core.DensityMatrixSimulator()
dm_sim2.run(prog_pure, noise)
dm_noisy = dm_sim2.density_matrix()

purity_noisy = np.trace(dm_noisy @ dm_noisy).real
print(f"Noisy state purity: {purity_noisy:.6f}")
# 输出: < 1.0 (mixed state)

纯度定义为：

$$\text{Tr}({\rho }^2)=\sum _{i,j}|{\rho }_{ij}{|}^2$$

纯度满足 $1/d\le \text{Tr}({\rho }^2)\le 1$，其中 $d=2^n$ 是希尔伯特空间(Hilbert space)维度。

冯诺依曼熵

冯诺依曼熵(von Neumann entropy)量化量子态中的不确定性（或纠缠量）。纯态的熵为零，最大混合态的熵最大。

from pyqpanda3 import core
import numpy as np

def von_neumann_entropy(rho):
    """计算冯诺依曼熵 S(rho) = -Tr(rho * log2(rho))。"""
    eigenvalues = np.linalg.eigvalsh(rho.real)
    # 过滤掉零特征值以避免 log(0)
    eigenvalues = eigenvalues[eigenvalues > 1e-12]
    entropy = -np.sum(eigenvalues * np.log2(eigenvalues))
    return entropy

dm_sim = core.DensityMatrixSimulator()

# 纯态：熵 = 0
prog_pure = core.QProg()
prog_pure << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog_pure)
dm_pure = dm_sim.density_matrix()
print(f"Pure state entropy: {von_neumann_entropy(dm_pure):.6f}")
# 输出: 0.000000

# 噪声态：熵 > 0
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.1),
    core.GateType.H
)
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.05),
    core.GateType.CNOT
)
dm_sim2 = core.DensityMatrixSimulator()
dm_sim2.run(prog_pure, noise)
dm_noisy = dm_sim2.density_matrix()
print(f"Noisy state entropy: {von_neumann_entropy(dm_noisy):.6f}")
# 输出: > 0 (信息因噪声而丢失)

# 通过约化密度矩阵计算纠缠熵
reduced = dm_sim.reduced_density_matrix([0])
print(f"Entanglement entropy S(rho_A): {von_neumann_entropy(reduced):.6f}")
# 对于贝尔态，S(rho_A) = 1.0（1 量子比特的最大值）

输出：

Pure state entropy: 0.000000
Noisy state entropy: 0.057234
Entanglement entropy S(rho_A): 1.000000

冯诺依曼熵定义为：

$$S(\rho )=-\text{Tr}(\rho {\log }_2\rho )=-\sum _i{\lambda }_i{\log }_2{\lambda }_i$$

其中 ${\lambda }_i$ 是 $\rho$ 的特征值(eigenvalue)。对于 $n$ 个量子比特的最大混合态：

$$S\left(\frac{I}{2^n}\right)=n$$

下表总结了这些量与物理态的关系：

物理量	纯态	最大混合态（$n$ 量子比特）	物理含义
纯度 $\text{Tr}({\rho }^2)$	1	$2^{-n}$	混合度
冯诺依曼熵 $S(\rho )$	0	$n$	不确定性
迹 $\text{Tr}(\rho )$	1	1	归一化

---

Stabilizer

Stabilizer 模拟器使用稳定子形式(stabilizer formalism)高效模拟 Clifford 线路。它可以处理非常大量的量子比特，因为它仅跟踪泡利稳定子生成元(stabilizer generator)而非完整态矢量。

Clifford 门

Stabilizer 仅支持 Clifford 门，Clifford 门在共轭变换下将泡利算子映射为泡利算子：

• H（Hadamard）：映射 $X\leftrightarrow Z$，$Y\leftrightarrow -Y$
• S（相位）：映射 $X\to Y$，$Y\to -X$
• CNOT（受控非）：映射 $X\otimes I\to X\otimes X$，$I\otimes Z\to Z\otimes Z$

像 T、RX、RY（角度非 $\pi /2$）等门不是 Clifford 门。

用法

from pyqpanda3 import core

# 构建 Clifford 线路
prog = core.QProg()
prog << core.H(0)
prog << core.S(1)
prog << core.CNOT(0, 1)
prog << core.measure([0, 1], [0, 1])

# 在 Stabilizer 模拟器上运行
stab = core.Stabilizer()
stab.run(prog, shots=10000)

result = stab.result()
print(result.get_counts())

带噪声运行

# Stabilizer 也支持噪声模拟
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.pauli_x_error(0.01),
    core.GateType.H
)

stab = core.Stabilizer()
stab.run(prog, shots=10000, model=noise)

何时使用 Stabilizer

• 验证 Clifford 线路的正确性
• 量子纠错模拟（表面码(surface code)、Steane 码）
• 大量子比特数基准测试（可支持 1000+ 量子比特）
• 在完整模拟器上运行前的快速健全性检查

稳定子形式详解

稳定子形式不通过振幅而是通过使态不变的泡利算子集合（稳定子）来表示量子态。对于 $n$ 个量子比特，一个稳定子态由 $n$ 个独立的对易泡利算子唯一确定。

稳定子模拟的工作原理

稳定子态 $|\psi ⟩$ 由其稳定子群 $\mathcal{S}=\{P\in {\mathcal{P}}_n:P|\psi ⟩=|\psi ⟩\}$ 定义，其中 ${\mathcal{P}}_n$ 是 $n$ 量子比特泡利群。该群由 $n$ 个独立生成元 $g_1,g_2,\dots ,g_n$ 生成：

$$\mathcal{S}=⟨g_1,g_2,\dots ,g_n⟩$$

关键洞见在于 Clifford 门通过共轭将稳定子生成元变换为新的稳定子生成元：

$$g_k\to U{g}_k{U}^{†}$$

此更新的时间复杂度为每门 $O(n^2)$，使得在 $n$ 个量子比特上 $m$ 个门的线路的整体复杂度为 $O(n^{2}m)$，相比之下态矢量模拟为 $O(2^{n}m)$。

稳定子生成元示例

考虑计算基态 $|00⟩$。其稳定子生成元为 $Z\otimes I$ 和 $I\otimes Z$：

from pyqpanda3 import core

# 示例 1：|00> 态
# 稳定子：ZI, IZ（约定：第一个字母 = 量子比特 0）
# 含义：Z_0|00> = +|00> 且 Z_1|00> = +|00>

prog = core.QProg()
# 未施加任何门 -- 态为 |00>
prog << core.measure([0, 1], [0, 1])

stab = core.Stabilizer()
stab.run(prog, shots=1000)
print(f"|00> state: {stab.result().get_counts()}")
# 输出: {'00': 1000}  (确定性的)

# 示例 2：对量子比特 0 施加 H -> |+0>
# 稳定子变化：ZI -> XI, IZ 保持 IZ
# 新稳定子：XI, IZ
prog2 = core.QProg()
prog2 << core.H(0)
prog2 << core.measure([0, 1], [0, 1])

stab2 = core.Stabilizer()
stab2.run(prog2, shots=1000)
print(f"|+0> state: {stab2.result().get_counts()}")
# 输出: {'00': ~500, '10': ~500}  (量子比特 0 随机)

# 示例 3：贝尔态（H + CNOT）
# |00> --对 q0 施加 H--> |+0> --施加 CNOT--> (|00>+|11>)/sqrt(2)
# 稳定子：ZI -> XI -> XX
#              IZ -------> IZ -> ZZ
# 最终稳定子：XX, ZZ
prog3 = core.QProg()
prog3 << core.H(0) << core.CNOT(0, 1)
prog3 << core.measure([0, 1], [0, 1])

stab3 = core.Stabilizer()
stab3.run(prog3, shots=1000)
print(f"Bell state: {stab3.result().get_counts()}")
# 输出: {'00': ~500, '11': ~500}  (关联的)

下表追踪了贝尔态线路中的稳定子生成元变化：

步骤	门	生成元 1	生成元 2	态
初始	--	$Z\otimes I$	$I\otimes Z$	$|00⟩$
H(0) 后	H	$X\otimes I$	$I\otimes Z$	$|+0⟩$
CNOT(0,1) 后	CNOT	$X\otimes X$	$Z\otimes Z$	$\frac{1}{\sqrt{2}}(|00⟩+|11⟩)$

扩展优势

稳定子表示使用大小为 $O(n^2)$ 的表(tableau)，而非大小为 $O(2^n)$ 的态矢量：

量子比特数 $n$	态矢量大小	稳定子表大小	加速倍数
10	16 KB	~100 字节	$1.6\times {10}^5$
20	16 MB	~400 字节	$4\times {10}^7$
30	16 GB	~900 字节	$1.8\times {10}^{10}$
100	16 GB $\times 2^{70}$	~10 KB	不可行 vs. 即时

---

PartialAmplitudeQVM

PartialAmplitudeQVM 计算态矢量的特定振幅，而无需维护完整的 $2^n$ 矢量。当你只需要某些概率振幅时，这对于大型线路非常有用。

用法

from pyqpanda3 import core

# 构建线路
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)

# 创建部分振幅模拟器
pavm = core.PartialAmplitudeQVM()
pavm.run(prog)

# 按基态字符串获取特定振幅
amplitudes = pavm.get_state_vector(["00", "01", "10", "11"])
print(f"Amplitude |00⟩: {amplitudes[0]}")
print(f"Amplitude |01⟩: {amplitudes[1]}")
print(f"Amplitude |10⟩: {amplitudes[2]}")
print(f"Amplitude |11⟩: {amplitudes[3]}")

何时使用 PartialAmplitudeQVM

• 只需要振幅的子集（例如，用于计算特定观测量）
• 线路太大，无法进行完整态矢量模拟
• 想要验证特定基态的概率

---

期望值计算

pyqpanda3 提供了用于计算观测量期望值(expectation value)的全局函数，这对于变分算法(variational algorithm)至关重要。

expval_hamiltonian

计算哈密顿量(Hamiltonian)的期望值 $⟨\psi |H|\psi ⟩$：

from pyqpanda3 import core, hamiltonian

# 定义哈密顿量
H = hamiltonian.PauliOperator({"Z0": 1.0, "Z0 Z1": 0.5})

# 构建态制备线路
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)

# 计算期望值
exp_val = core.expval_hamiltonian(
    prog,
    H,
    shots=1000,
    model=core.NoiseModel(),
    used_threads=4,
    backend="CPU"
)
print(f"⟨H⟩ = {exp_val}")

expval_pauli_operator

计算泡利算子的期望值：

from pyqpanda3 import core, hamiltonian

# 定义泡利算子
pauli_op = hamiltonian.PauliOperator({"X0 Y1": 0.5, "Z0": 0.3})

# 构建线路
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)

# 使用可选噪声模型计算期望值
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.001),
    core.GateType.CNOT
)

exp_val = core.expval_pauli_operator(
    prog,
    pauli_op,
    shots=1000,
    model=noise,
    used_threads=4,
    backend="CPU"
)
print(f"⟨P⟩ = {exp_val}")

多线程和 GPU 后端

两个期望值函数都支持并行化和 GPU 后端：

# 使用 8 个线程进行并行计算
exp_val = core.expval_hamiltonian(
    prog, H,
    shots=10000,
    used_threads=8,
    backend="CPU"
)

# 使用 GPU 后端（如果使用 CUDA 编译）
exp_val_gpu = core.expval_hamiltonian(
    prog, H,
    shots=10000,
    used_threads=4,
    backend="GPU"
)

---

模拟器选择指南

决策标准

场景	推荐模拟器	原因
通用量子线路	CPUQVM	完整态矢量、噪声支持
大型 Clifford 线路	Stabilizer	多项式扩展、快速
精确噪声分析	DensityMatrixSimulator	精确混合态演化
大型线路、GPU 可用	GPUQVM	硬件加速
大型线路、特定振幅	PartialAmplitudeQVM	减少内存
变分算法	CPUQVM + expval 函数	高效的期望值计算
纠错模拟	Stabilizer	可处理大量量子比特

---

性能考虑

内存使用

对于 $n$ 个量子比特，内存需求如下：

模拟器	内存复杂度	20 量子比特	25 量子比特	30 量子比特
CPUQVM	$O(2^n)$	16 MB	512 MB	16 GB
GPUQVM	$O(2^n)$（GPU 显存）	16 MB	512 MB	16 GB
DensityMatrixSimulator	$O(4^n)$	256 MB	256 GB	256 TB
Stabilizer	$O(n^2)$	~1 KB	~1 KB	~2 KB
PartialAmplitudeQVM	$O(2^n)$ 部分	可变	可变	可变

速度技巧

1. Clifford 线路使用 Stabilizer -- 快数个数量级
2. 在统计精度允许时减少采样次数（1000 次采样通常足够）
3. 期望值计算使用多线程（used_threads 参数）
4. 25+ 量子比特的线路使用 GPU 后端
5. 模拟前展平程序以减少开销：

from pyqpanda3 import core

prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
flat_prog = prog.flatten()

线路深度和门计数分析

在运行模拟之前，分析线路以估计资源需求是有用的。了解线路深度(depth)和门计数(gate count)有助于选择合适的模拟器并预测执行时间。

线路深度是当作用于不同量子比特的门并行执行时所需的时间步数。门计数是线路中门的总数。

from pyqpanda3 import core

# 构建中等复杂度的线路
prog = core.QProg()
prog << core.H(0)
prog << core.H(1)
prog << core.CNOT(0, 1)
prog << core.H(2)
prog << core.CNOT(1, 2)
prog << core.S(0)
prog << core.CNOT(0, 2)
prog << core.measure([0, 1, 2], [0, 1, 2])

# 展平以获取规范门列表
flat_prog = prog.flatten()

# 可以检查展平后的线路来计算门数和估计深度
# 序列中的每个门代表一个操作
print(f"Circuit contains the sequence of primitive gates after flattening")

# 构建较小的子线路来演示深度计数
small = core.QProg()
small << core.H(0)       # 层 1
small << core.H(1)       # 层 1（与 H(0) 并行）
small << core.CNOT(0, 1) # 层 2（依赖于量子比特 0, 1）
small << core.H(2)       # 层 2（与 CNOT 并行，不同量子比特）
small << core.CNOT(1, 2) # 层 3（依赖于层 2 的量子比特 1）
print("Small circuit: depth ~3, 5 gates")

门计数和深度共同决定了模拟成本。对于在 $n$ 个量子比特上有 $m$ 个门的线路：

• 态矢量模拟器（CPUQVM、GPUQVM）：时间复杂度 $O(m\cdot 2^n)$
• 密度矩阵模拟器：时间复杂度 $O(m\cdot 4^n)$
• 稳定子模拟器：时间复杂度 $O(m\cdot n^2)$

使用下表估计线路是否能放入内存：

门数	量子比特数	CPUQVM 时间估计	DensityMatrix 可行？
100	10	< 1 ms	是（~1 MB）
100	20	~10 ms	勉强（~256 MB）
1000	25	~1 s	否（~256 GB）
10000	30	~30 s	否
1000	100	不适用	否；使用 Stabilizer

线程配置和性能调优

pyqpanda3 提供多种机制来控制并行性并为你的硬件优化性能。了解这些选项可以显著减少大型线路的模拟时间。

控制线程数

used_threads 参数控制期望值计算和批量模拟中的并行性：

from pyqpanda3 import core

# 构建变分式线路（不含测量，用于期望值计算）
prog = core.QProg()
prog << core.H(0)
prog << core.RY(1, 0.5)
prog << core.CNOT(0, 1)

# 使用 1 个线程（单线程，适合调试）
prog_meas = core.QProg()
prog_meas << core.H(0)
prog_meas << core.RY(1, 0.5)
prog_meas << core.CNOT(0, 1)
prog_meas << core.measure([0, 1], [0, 1])
machine_1t = core.CPUQVM()
machine_1t.run(prog_meas, shots=10000)

# 对于期望值，显式控制线程
from pyqpanda3 import hamiltonian
H = hamiltonian.PauliOperator({"Z0": 1.0, "Z0 Z1": 0.5})

# 单线程计算
exp_val_1 = core.expval_hamiltonian(
    prog, H, shots=10000, used_threads=1, backend="CPU"
)

# 多线程计算（4 个线程）
exp_val_4 = core.expval_hamiltonian(
    prog, H, shots=10000, used_threads=4, backend="CPU"
)

# 最大并行性（使用所有可用核心）
import os
max_threads = os.cpu_count()
exp_val_max = core.expval_hamiltonian(
    prog, H, shots=10000, used_threads=max_threads, backend="CPU"
)
print(f"System has {max_threads} CPU cores available")

选择最优线程数

最优线程数取决于线路大小和线程开销之间的关系：

$$T_{\text{total}}(p)=\frac{T_{\text{serial}}}{p}+T_{\text{overhead}}(p)$$

其中 $p$ 是线程数，$T_{\text{serial}}$ 是单线程时间，$T_{\text{overhead}}(p)$ 是同步开销。

线程配置遵循以下指南：

场景	推荐线程数	理由
小型线路（$n<15$）	1-2	线程开销超过收益
中型线路（$15\le n<25$）	4-8	良好的并行加速
大型线路（$n\ge 25$）	所有可用核心	最大化吞吐量
GPU 可用	使用 GPU 后端	GPU 并行性远超 CPU
具有多个泡利项的期望值	所有核心	每项可独立求值

---

多模拟器对比工作流

在开发量子算法时，在多个模拟器之间比较结果以验证正确性并了解噪声影响是非常有价值的。本节提供系统的基准测试工作流。

基准测试设计

好的基准测试应衡量：

1. 正确性：所有模拟器对理想线路产生一致的结果
2. 性能：执行时间随线路大小按预期扩展
3. 噪声影响：含噪声模拟器显示出与理想结果的预期偏差

from pyqpanda3 import core
import time

def benchmark_simulator(simulator_cls, prog, shots, label, **kwargs):
    """在单个模拟器上运行基准测试并返回计时数据。"""
    start = time.perf_counter()
    sim = simulator_cls()
    sim.run(prog, shots=shots, **kwargs)
    elapsed = time.perf_counter() - start
    result = sim.result()
    return {
        "label": label,
        "time": elapsed,
        "counts": result.get_counts() if hasattr(result, 'get_counts') else None,
    }

# 构建 4 量子比特 GHZ 态用于基准测试
prog_ghz = core.QProg()
prog_ghz << core.H(0)
for i in range(3):
    prog_ghz << core.CNOT(i, i + 1)
prog_ghz << core.measure([0, 1, 2, 3], [0, 1, 2, 3])

# 不含测量的版本（用于密度矩阵模拟器）
prog_ghz_no_meas = core.QProg()
prog_ghz_no_meas << core.H(0)
for i in range(3):
    prog_ghz_no_meas << core.CNOT(i, i + 1)

print("=== Ideal Circuit Benchmark ===")
shots = 10000

# 基准测试 CPUQVM
cpu_result = benchmark_simulator(core.CPUQVM, prog_ghz, shots, "CPUQVM")
print(f"CPUQVM: {cpu_result['time']:.4f}s, counts: {cpu_result['counts']}")

# 基准测试 Stabilizer
stab_result = benchmark_simulator(core.Stabilizer, prog_ghz, shots, "Stabilizer")
print(f"Stabilizer: {stab_result['time']:.4f}s, counts: {stab_result['counts']}")

# 基准测试 DensityMatrixSimulator
start_dm = time.perf_counter()
dm_sim = core.DensityMatrixSimulator()
dm_sim.run(prog_ghz_no_meas)
dm_time = time.perf_counter() - start_dm
dm_probs = dm_sim.state_probs()
print(f"DensityMatrix: {dm_time:.4f}s, probs: {[f'{p:.4f}' for p in dm_probs]}")

比较理想与含噪声结果

from pyqpanda3 import core

# 定义具有现实错误率的噪声模型
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.005), core.GateType.H
)
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.01), core.GateType.CNOT
)

# 构建 3 量子比特线路
prog = core.QProg()
prog << core.H(0)
prog << core.CNOT(0, 1)
prog << core.CNOT(1, 2)
prog << core.measure([0, 1, 2], [0, 1, 2])

# 理想模拟
cpu = core.CPUQVM()
cpu.run(prog, shots=10000)
ideal_counts = cpu.result().get_counts()
print(f"Ideal counts:   {ideal_counts}")

# 含噪声模拟（CPUQVM，基于采样的噪声）
cpu_noisy = core.CPUQVM()
cpu_noisy.run(prog, shots=10000, model=noise)
noisy_counts = cpu_noisy.result().get_counts()
print(f"Noisy counts:   {noisy_counts}")

# 精确噪声模拟（DensityMatrixSimulator）
prog_no_meas = core.QProg()
prog_no_meas << core.H(0) << core.CNOT(0, 1) << core.CNOT(1, 2)
dm_sim = core.DensityMatrixSimulator()
dm_sim.run(prog_no_meas, noise)
exact_probs = dm_sim.state_probs()
print(f"Exact probs:    {[f'{p:.4f}' for p in exact_probs]}")

# 计算理想与含噪声之间的总变差距离
ideal_p = {k: v / 10000 for k, v in ideal_counts.items()}
noisy_p = {k: v / 10000 for k, v in noisy_counts.items()}
all_keys = set(list(ideal_p.keys()) + list(noisy_p.keys()))
tvd = sum(abs(ideal_p.get(k, 0) - noisy_p.get(k, 0)) for k in all_keys) / 2
print(f"Total variation distance: {tvd:.4f}")

两个分布 $P$ 和 $Q$ 之间的总变差距离(total variation distance)为：

$$\text{TVD}(P,Q)=\frac{1}{2}\sum _x|P(x)-Q(x)|$$

此指标的值域为 0（完全相同的分布）到 1（完全不重叠的分布）。

跨模拟器验证工作流

---

完整示例：模拟器比较

本节提供了一个综合示例，在相同线路上使用所有模拟器，并对结果及其物理解释进行详细分析。

完整比较代码

from pyqpanda3 import core
import numpy as np

# 构建 3 量子比特 GHZ 态
prog = core.QProg()
prog << core.H(0)
prog << core.CNOT(0, 1)
prog << core.CNOT(1, 2)
prog << core.measure([0, 1, 2], [0, 1, 2])

# 不含测量的版本（用于密度矩阵分析）
prog_no_meas = core.QProg()
prog_no_meas << core.H(0) << core.CNOT(0, 1) << core.CNOT(1, 2)

print("=" * 60)
print("  COMPREHENSIVE SIMULATOR COMPARISON: 3-Qubit GHZ State")
print("=" * 60)
print(f"\nTheoretical state: |GHZ> = (|000> + |111>) / sqrt(2)")
print(f"Expected: P(|000>) = 0.5, P(|111>) = 0.5, all others = 0.0")

CPUQVM 结果

print("\n--- CPUQVM (State Vector, Shot-Based) ---")
cpu = core.CPUQVM()
cpu.run(prog, shots=10000)
cpu_counts = cpu.result().get_counts()
print(f"Raw counts: {cpu_counts}")

# 转换为概率以便比较
total = sum(cpu_counts.values())
cpu_probs = {k: v / total for k, v in sorted(cpu_counts.items())}
print("Empirical probabilities:")
for state, prob in cpu_probs.items():
    print(f"  P(|{state}>) = {prob:.4f}")

# 计算统计误差
for state in ['000', '111']:
    p = cpu_probs.get(state, 0)
    std_err = np.sqrt(p * (1 - p) / total)
    print(f"  {state}: {p:.4f} +/- {std_err:.4f} (1-sigma)")

预期输出显示围绕理想值的统计波动：

--- CPUQVM (State Vector, Shot-Based) ---
Raw counts: {'000': 5023, '111': 4977}
Empirical probabilities:
  P(|000>) = 0.5023
  P(|111>) = 0.4977
  000: 0.5023 +/- 0.0050 (1-sigma)
  111: 0.4977 +/- 0.0050 (1-sigma)

DensityMatrixSimulator 结果

print("\n--- DensityMatrixSimulator (Ideal, No Noise) ---")
dm_sim = core.DensityMatrixSimulator()
dm_sim.run(prog_no_meas)
ideal_probs = dm_sim.state_probs()
print(f"Exact probabilities: {[f'{p:.6f}' for p in ideal_probs]}")
print(f"P(|000>) = {dm_sim.state_prob('000'):.6f}")
print(f"P(|111>) = {dm_sim.state_prob('111'):.6f}")

# 验证纯度（纯态应为 1.0）
dm = dm_sim.density_matrix()
purity = np.trace(dm @ dm).real
print(f"Purity: {purity:.6f} (pure state = 1.0)")

# 约化密度矩阵（对 2 个量子比特求迹留下最大混合态）
reduced_0 = dm_sim.reduced_density_matrix([0])
print(f"\nReduced state for qubit 0:")
print(f"  rho_0 = {reduced_0}")
print(f"  Tr(rho_0^2) = {np.trace(reduced_0 @ reduced_0).real:.4f}")
print(f"  => Maximally mixed (entangled with other qubits)")

理想密度矩阵的预期输出：

--- DensityMatrixSimulator (Ideal, No Noise) ---
Exact probabilities: ['0.500000', '0.000000', '0.000000', '0.000000', '0.000000', '0.000000', '0.000000', '0.500000']
P(|000>) = 0.500000
P(|111>) = 0.500000
Purity: 1.000000 (pure state = 1.0)

Reduced state for qubit 0:
  rho_0 = [[0.5 0. ], [0.  0.5]]
  Tr(rho_0^2) = 0.5000
  => Maximally mixed (entangled with other qubits)

含噪声模拟结果

print("\n--- DensityMatrixSimulator (Noisy) ---")
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.02), core.GateType.CNOT
)
dm_sim2 = core.DensityMatrixSimulator()
dm_sim2.run(prog_no_meas, noise)
noisy_probs = dm_sim2.state_probs()
print(f"Noisy probabilities: {[f'{p:.6f}' for p in noisy_probs]}")

# 噪声将概率泄漏到先前为零的态
dm_noisy = dm_sim2.density_matrix()
purity_noisy = np.trace(dm_noisy @ dm_noisy).real
print(f"Purity: {purity_noisy:.6f} (< 1.0 indicates mixed state)")

# 量化噪声影响
leaked = sum(noisy_probs[i] for i in range(len(noisy_probs))
             if i not in [0, 7])
print(f"Probability leaked to other states: {leaked:.6f}")
print(f"  (0.0 = no noise, 0.5 = maximally noisy)")

在每个 CNOT 门上有 2% 去极化噪声的情况下，输出显示：

--- DensityMatrixSimulator (Noisy) ---
Noisy probabilities: ['0.4838', '0.0016', '0.0016', '0.0050', '0.0016', '0.0050', '0.0050', '0.4964']
Purity: 0.9871 (< 1.0 indicates mixed state)
Probability leaked to other states: 0.019800

Stabilizer 结果

print("\n--- Stabilizer (Clifford Simulation) ---")
stab = core.Stabilizer()
stab.run(prog, shots=10000)
stab_counts = stab.result().get_counts()
print(f"Stabilizer counts: {stab_counts}")

# Stabilizer 应该给出与 CPUQVM 相同的统计
# 因为 GHZ 线路是 Clifford 线路
total_stab = sum(stab_counts.values())
stab_probs = {k: v / total_stab for k, v in sorted(stab_counts.items())}
print("Stabilizer probabilities:")
for state, prob in stab_probs.items():
    print(f"  P(|{state}>) = {prob:.4f}")

--- Stabilizer (Clifford Simulation) ---
Stabilizer counts: {'000': 4989, '111': 5011}
Stabilizer probabilities:
  P(|000>) = 0.4989
  P(|111>) = 0.5011

结果总结

print("\n" + "=" * 60)
print("  SUMMARY OF RESULTS")
print("=" * 60)
print("""
| Simulator               | P(|000>) | P(|111>) | Purity  | Noise? |
|-------------------------|----------|----------|---------|--------|
| CPUQVM (ideal)          | 0.5023   | 0.4977   | N/A     | No     |
| DensityMatrix (ideal)   | 0.500000 | 0.500000 | 1.000   | No     |
| DensityMatrix (noisy)   | 0.4838   | 0.4964   | 0.987   | Yes    |
| Stabilizer              | 0.4989   | 0.5011   | N/A     | No     |

Key observations:
1. All ideal simulators agree on the 50/50 distribution within
   statistical error (CPUQVM, Stabilizer) or exactly (DensityMatrix).
2. The noisy DensityMatrix shows probability leaking into non-GHZ states,
   with purity dropping below 1.0.
3. The Stabilizer is the fastest for this Clifford circuit but cannot
   simulate non-Clifford gates.
4. Reduced density matrices confirm entanglement: each single-qubit
   reduced state is maximally mixed (purity = 0.5).
""")

理想态与含噪声态之间的保真度(fidelity)量化了整体噪声影响：

$$F({\rho }_{\text{ideal}},{\rho }_{\text{noisy}})={\left(\text{Tr}\sqrt{\sqrt{{\rho }_{\text{ideal}}}{\rho }_{\text{noisy}}\sqrt{{\rho }_{\text{ideal}}}}\right)}^2$$

对于纯态与混合态，简化为：

$$F(|\psi ⟩,\rho )=⟨\psi |\rho |\psi ⟩$$

---

下一步

• 噪声模拟 - 深入了解噪声模型和错误信道
• 哈密顿量与泡利算符 - 在变分算法中使用期望值
• 变分线路 - 构建带梯度计算的参数化线路