模拟

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

---

概述

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

graph TD
    A[Choose a Simulator] --> B{Need noise?}
    B -->|Yes| C{Circuit type?}
    B -->|No| D{Circuit size?}
    C -->|Clifford only| E[Stabilizer]
    C -->|General| F[DensityMatrixSimulator]
    D -->|Small < 25 qubits| G[CPUQVM]
    D -->|Large state vector| H{CUDA available?}
    H -->|Yes| I[GPUQVM]
    H -->|No| J[PartialAmplitudeQVM]
    D -->|Only one amplitude| J2[SingleAmplitudeQVM]
    D -->|Weak entanglement| J3[MPSQVM]
    D -->|Only an expectation value| J4[SparseHeisenbergSimulator]
    E --> K[Fast Clifford simulation]
    F --> L[Full noise modeling]
    G --> M[Statevector simulation]
    I --> N[GPU-accelerated]
    J --> O[Partial amplitudes only]
    J2 --> O2[Single amplitude query]
    J3 --> O3[MPS simulation]
    J4 --> O4[Heisenberg-picture expectation]

模拟器	态表示	噪声支持	最佳用途
CPUQVM	态矢量(state vector)	是（基于采样）	通用模拟
GPUQVM	态矢量（GPU）	是	带 CUDA 的大型线路
DensityMatrixSimulator	密度矩阵(density matrix)	是（精确）	噪声分析、态层析
Stabilizer	稳定子态(stabilizer state)	是	Clifford 线路、快速基准测试
PartialAmplitudeQVM	部分态矢量	否	大型线路、特定振幅
SingleAmplitudeQVM	张量网络	否	大型线路中查询单个振幅
MPSQVM	矩阵乘积态	是	弱纠缠的大型线路
SparseHeisenbergSimulator	海森堡绘景（稀疏）	否	快速计算泡利哈密顿量期望值（非采样、无态矢量）

---

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}

运行结果：

Measurement counts: {'11': 4972, '00': 5028}

获取结果

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

from pyqpanda3 import core

# 先运行贝尔态线路
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1) << core.measure([0, 1], [0, 1])
machine = core.CPUQVM()
machine.run(prog, shots=10000)
result = machine.result()

# 以字典形式获取测量计数
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()

运行结果：

Shots: 10000

获取态矢量

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

from pyqpanda3 import core

# 构建不含测量的线路
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)]

运行结果：

State vector: [(0.7071067811865476+0j), 0j, 0j, (0.7071067811865476+0j)]

$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' 中看到少量计数

运行结果：

{'10': 49, '01': 58, '00': 4979, '11': 4914}

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())

运行结果：

{'11': 2551, '10': 2462, '01': 2461, '00': 2526}

展平在性能分析(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}")

运行结果：

|000>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000
|001>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000
|010>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000
|011>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000
|100>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000
|101>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000
|110>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000
|111>: amplitude = 0.3536+0.0000j, |amp| = 0.3536, phase = 0.0000

对于处于等概率叠加的 $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]}")

运行结果：

Full probability distribution:
  P(|000>) = 0.500000
  P(|001>) = 0.000000
  P(|010>) = 0.000000
  P(|011>) = 0.000000
  P(|100>) = 0.000000
  P(|101>) = 0.000000
  P(|110>) = 0.000000
  P(|111>) = 0.500000

Marginal distribution (qubits 0, 2):
  P(|00>) = 0.500000
  P(|01>) = 0.000000
  P(|10>) = 0.000000
  P(|11>) = 0.500000

Probability list: ['0.500000', '0.000000', '0.000000', '0.000000', '0.000000', '0.000000', '0.000000', '0.500000']

对于 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])

# Use CPUQVM (GPUQVM available with USE_CUDA=ON build)
machine = core.CPUQVM()
machine.run(prog, shots=10000)
result = machine.result()
print(result.get_counts())

运行结果：

{'11': 4938, '00': 5062}

何时使用 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}")

运行结果：

Density matrix shape: (4, 4)
Density matrix:
[[0.5+0.j 0. +0.j 0. +0.j 0.5+0.j]
 [0. +0.j 0. +0.j 0. +0.j 0. +0.j]
 [0. +0.j 0. +0.j 0. +0.j 0. +0.j]
 [0.5+0.j 0. +0.j 0. +0.j 0.5+0.j]]

概率

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

from pyqpanda3 import core

# 创建密度矩阵模拟器并制备贝尔态
dm_sim = core.DensityMatrixSimulator()
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog)

# 按索引获取特定态的概率
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}")

运行结果：

P(|00⟩) = 0.5000000000000001
P(|01⟩) = 0.0
All probabilities: [0.5000000000000001, 0.0, 0.0, 0.5000000000000001]
Qubit [0,1] probabilities: [0.5000000000000001, 0.0, 0.0, 0.5000000000000001]

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

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

约化密度矩阵

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

from pyqpanda3 import core

# 创建密度矩阵模拟器并制备贝尔态
dm_sim = core.DensityMatrixSimulator()
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog)

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

运行结果：

Reduced density matrix for qubit 0:
[[0.5+0.j 0. +0.j]
 [0. +0.j 0.5+0.j]]

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

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

噪声模拟

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

from pyqpanda3 import core

# Create noise model
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.amplitude_damping_error(0.05),
    core.GateType.H
)

# Run with noise
dm_sim = core.DensityMatrixSimulator()
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog, noise)

# Check resulting mixed state
dm = dm_sim.density_matrix()
print(f"Trace: {dm_sim.density_matrix().trace()}")
print(f"State probabilities: {dm_sim.state_probs()}")

运行结果：

Trace: (1.0000000000000002+0j)
State probabilities: [0.5250000000000001, 0.0, 0.0, 0.47500000000000003]

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

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

密度矩阵操作

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

偏迹和约化密度矩阵

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

from pyqpanda3 import core

# 创建最大纠缠贝尔态
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)$$

graph LR
    A[Full Density Matrix rho_AB] --> B[Partial Trace Tr_B]
    A --> C[Partial Trace Tr_A]
    B --> D["Reduced rho_A = I/2"]
    C --> E["Reduced rho_B = I/2"]
    D --> F[Maximally mixed]
    E --> G[Maximally mixed]
    F --> H["Purity = 0.5 < 1"]

纯度计算

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

from pyqpanda3 import core
import numpy as np

dm_sim = core.DensityMatrixSimulator()

# Case 1: Pure state (no noise)
prog_pure = core.QProg()
prog_pure << core.H(0) << core.CNOT(0, 1)
dm_sim.run(prog_pure)
dm_pure = dm_sim.density_matrix()

# Purity is Tr(rho^2)
purity_pure = np.trace(dm_pure @ dm_pure).real
print(f"Pure state purity: {purity_pure:.6f}")
# Output: 1.000000

# Case 2: Mixed state (with depolarizing noise on H gate)
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.1),
    core.GateType.H
)
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}")
# Output: < 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):
    """Compute von Neumann entropy S(rho) = -Tr(rho * log2(rho))."""
    eigenvalues = np.linalg.eigvalsh(rho.real)
    # Filter out zero eigenvalues to avoid log(0)
    eigenvalues = eigenvalues[eigenvalues > 1e-12]
    entropy = -np.sum(eigenvalues * np.log2(eigenvalues))
    return entropy

dm_sim = core.DensityMatrixSimulator()

# Pure state: entropy = 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}")
# Output: 0.000000

# Noisy state: entropy > 0
noise = core.NoiseModel()
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.1),
    core.GateType.H
)
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}")
# Output: > 0 (information lost to noise)

# Entanglement entropy via reduced density matrix
reduced = dm_sim.reduced_density_matrix([0])
print(f"Entanglement entropy S(rho_A): {von_neumann_entropy(reduced):.6f}")
# For a Bell state, S(rho_A) = 1.0 (maximum for 1 qubit)

Output:

Pure state entropy: 0.000000
Noisy state entropy: 0.122292
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())

带噪声运行

from pyqpanda3 import core

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

# 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⟩$$

graph TD
    A[Initial state: 0^n] --> B["Stabilizers: Z_1, Z_2, ..., Z_n"]
    B --> C{Apply Clifford gate}
    C -->|H on qubit i| D["Conjugate generators: X_i <-> Z_i"]
    C -->|S on qubit i| E["Conjugate generators: Z_i -> Y_i"]
    C -->|CNOT i,j| F["Conjugate generators: X_i -> X_i X_j, Z_j -> Z_i Z_j"]
    D --> G[Updated stabilizer tableau]
    E --> G
    F --> G
    G --> H{More gates?}
    H -->|Yes| C
    H -->|No| I[Read measurement outcomes]

关键洞见在于 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, '01': ~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)

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

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

何时使用 PartialAmplitudeQVM

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

---

期望值计算

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

expval_hamiltonian

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

from pyqpanda3 import core, hamiltonian

# 定义哈密顿量
H = hamiltonian.Hamiltonian({"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 后端：

from pyqpanda3 import core, hamiltonian

# 构建线路并定义可观测量
prog = core.QProg()
prog << core.H(0) << core.CNOT(0, 1)
H = hamiltonian.Hamiltonian({"Z0": 1.0, "Z0 Z1": 0.5})

# 使用 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"
)

---

模拟器选择指南

graph LR
    subgraph "Simulation Task"
        A[Start]
    end

    subgraph "Circuit Analysis"
        B{Clifford only?}
        C{Need noise?}
        D{How many qubits?}
    end

    subgraph "Simulator"
        E[Stabilizer]
        F[DensityMatrixSimulator]
        G[CPUQVM]
        H[GPUQVM]
        I[PartialAmplitudeQVM]
    end

    A --> B
    B -->|Yes| C
    B -->|No| C
    C -->|Yes, exact| F
    C -->|Yes, sampling| D
    C -->|No| D
    D -->|< 25| G
    D -->|25-30, CUDA| H
    D -->|> 30| I

    style E fill:#4CAF50,color:#fff
    style F fill:#2196F3,color:#fff
    style G fill:#FF9800,color:#fff
    style H fill:#9C27B0,color:#fff
    style I fill:#F44336,color:#fff

决策标准

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

---

精度选择与能力查询（0.5.0）

每个 CPU 系列模拟器在构造时接受算术 precision："float32" 大致将状态内存占用减半并加速内核，"float64"（默认）保持完整双精度。该选择在实例生命周期内固定，并通过只读属性暴露。

import numpy as np
from pyqpanda3.core import CPUQVM, H, QProg

cpu32 = CPUQVM(precision="float32")
cpu32.init_state(np.array([1, 0], dtype=np.complex64))

print(cpu32.precision, cpu32.supported_precisions)
print(cpu32.current_state().dtype, cpu32.current_state().shape)

运行结果：

float32 ('float32', 'float64')
complex64 (2,)

current_state() 以模拟器原生精度返回状态（此处为 complex64，float64 为 complex128），且不做 dtype 转换。在 CPUQVM 上它镜像的是 init_state(...) 装入的状态——run() 不会填充它，读取 run() 产生的状态请改用 result().get_state_vector()。结果 API 不受影响：result() 与 get_state_vector() 返回与 float64 运行相同的对象，现有分析代码无需修改。

未知的精度名与"合法但该后端不支持"是两类不同的错误，均为 ValueError 的子类：

from pyqpanda3.core import CPUQVM

try:
    CPUQVM(precision="fp32")
except ValueError as e:
    print(type(e).__name__, "-", e)

try:
    CPUQVM(precision="exact")
except ValueError as e:
    print(type(e).__name__, "-", e)

运行结果：

InvalidPrecisionError - unknown precision 'fp32'
UnsupportedPrecisionError - backend 'CPUQVM' does not support precision 'exact'; supported precisions are: float32, float64

哪个后端支持什么？

simulator_capabilities() 无需构造任何对象即可回答。注册表还会报告某后端在当前主机上能否真正构造（无 CUDA 的机器上 GPUQVM 报告 available=False，此时 pyqpanda3.core.GPUQVM 为 None，不能构造；用 simulator_capabilities("GPUQVM").available 判断）：

from pyqpanda3.core import simulator_capabilities

for cap in simulator_capabilities():
    print(f"{cap.name:26s} {','.join(cap.supported_precisions):16s}"
          f" available={str(cap.available):5s} domain={cap.simd_domain}")

运行结果：

CPUQVM                     float32,float64  available=True  domain=StateVector
DensityMatrixSimulator     float32,float64  available=True  domain=Density
GPUQVM                     float64          available=False domain=NotApplicable
MPSQVM                     float32,float64  available=True  domain=Tensor
PartialAmplitudeQVM        float32,float64  available=True  domain=Amplitude
SingleAmplitudeQVM         float32,float64  available=True  domain=Amplitude
SparseHeisenbergSimulator  float32,float64  available=True  domain=Sparse
Stabilizer                 exact            available=True  domain=NotApplicable

SIMD 后端诊断

态矢量内核为多种 SIMD 指令集编译，并在进程启动时选定。pyqpanda3.core.simd 展示该决定，也可以覆盖它（偏好作用于之后构造的模拟器）：

from pyqpanda3.core import simd

cap = simd.capabilities()
print(cap.architecture, list(cap.available_backends), cap.active_backend)

simd.set_preference("scalar")  # 强制标量内核，例如用于调试
print(simd.active_backend())
simd.set_preference("auto")    # 恢复探测到的默认值
print(simd.active_backend())

运行结果（在支持 AVX2 的 x86_64 主机上）：

x86_64 ['scalar', 'sse2', 'avx2', 'avx512'] avx2
scalar
avx2

MPS 键维截断与运行统计

MPSQVM 以矩阵乘积态存储状态，因此其内存随**键维（bond dimension）**而非 2**n 增长。这正是它在约 30 个量子比特以上成为实用选择的原因，但它同时也是一种近似：每个双比特门之后，模拟器执行一次奇异值分解并丢弃最不重要的 Schmidt 系数。两个构造关键字控制这一丢弃过程，stats() 则报告实际发生了什么。

import numpy as np
from pyqpanda3.core import MPSQVM, QProg, RY, CZ

rng = np.random.default_rng(7)
prog = QProg()
for q in range(6):
    prog << RY(q, rng.uniform(0, np.pi))
for layer in range(4):
    pairs = [(0, 1), (2, 3), (4, 5)] if layer % 2 == 0 else [(1, 2), (3, 4)]
    for a, b in pairs:
        prog << CZ(a, b) << RY(a, rng.uniform(0, np.pi)) << RY(b, rng.uniform(0, np.pi))

capped = MPSQVM(precision="float64", max_bond_dimension=2)
capped.run(prog, 1)
print("capped :", capped.stats().discarded_weight, capped.stats().max_bond_dimension_seen)

exact = MPSQVM()
exact.run(prog, 1)
print("exact  :", exact.stats().discarded_weight, exact.stats().max_bond_dimension_seen)

运行结果：

capped : 0.024353811677256556 2
exact  : 1.8613363972962274e-34 4

两次运行描述同一条线路，唯一区别是键维上限：

构造关键字	含义	默认值
max_bond_dimension	每次截断后保留的键维硬上限。	None（不设上限）
truncation_atol	低于该权重的奇异值将被丢弃。	float32 为 1e-6，float64 为 1e-12

stats() 返回 MPSRunStats 记录，含两个字段：discarded_weight 是迄今因截断而丢弃的概率权重总量，max_bond_dimension_seen 是运行期间实际达到过的最大键维。在上面的受限运行中，上限生效于 2，discarded_weight 升至约 2.4e-2——这个数字就是截断所付出的保真度代价，因而是选择上限时需要关注的量。未设限的运行达到键维 4，仅丢弃了约 1.9e-34 的浮点噪声底。

实践中，上限通过扫描确定：逐步提高 max_bond_dimension 运行同一线路，直到 discarded_weight 低于应用所需的精度为止。由于容差默认值与精度相关，在相同的 max_bond_dimension 下，float32 模拟器本就比 float64 模拟器截断得更激进。

完整参考（含精度异常层次与单振幅缓存统计）见精度与能力。

---

性能考虑

内存使用

对于 $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)       # Layer 1
small << core.H(1)       # Layer 1 (parallel with H(0))
small << core.CNOT(0, 1) # Layer 2 (depends on qubits 0, 1)
small << core.H(2)       # Layer 2 (parallel with CNOT, different qubit)
small << core.CNOT(1, 2) # Layer 3 (depends on qubit 1 from layer 2)
print("Small circuit: depth ~3, 5 gates")

graph TD
    subgraph "Layer 1"
        H0["H(0)"]
        H1["H(1)"]
    end
    subgraph "Layer 2"
        CX01["CNOT(0,1)"]
        H2["H(2)"]
    end
    subgraph "Layer 3"
        CX12["CNOT(1,2)"]
    end
    H0 --> CX01
    H1 --> CX01
    CX01 --> CX12
    H2 --> CX12

门计数和深度共同决定了模拟成本。对于在 $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.Hamiltonian({"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)$ 是同步开销。

graph LR
    A[Start Performance Tuning] --> B{Qubit count?}
    B -->|n < 15| C[1-2 threads]
    B -->|15 <= n < 25| D[4-8 threads]
    B -->|n >= 25| E[Max threads or GPU]
    C --> F[Threading overhead dominates]
    D --> G[Good parallel speedup]
    E --> H[Essential for feasibility]

线程配置遵循以下指南：

场景	推荐线程数	理由
小型线路（$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):
    """Run a benchmark on a single simulator and return timing data."""
    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

# Define a noise model with realistic error rates
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
)

# Build a 3-qubit circuit
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])

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

# Noisy simulation (CPUQVM, shot-based noise)
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}")

# Exact noise simulation (DensityMatrixSimulator, single-qubit gate noise only)
prog_no_meas = core.QProg()
prog_no_meas << core.H(0) << core.CNOT(0, 1) << core.CNOT(1, 2)
noise_dm = core.NoiseModel()
noise_dm.add_all_qubit_quantum_error(
    core.depolarizing_error(0.005), core.GateType.H
)
dm_sim = core.DensityMatrixSimulator()
dm_sim.run(prog_no_meas, noise_dm)
exact_probs = dm_sim.state_probs()
print(f"Exact probs:    {[f'{p:.4f}' for p in exact_probs]}")

# Compute total variation distance between ideal and noisy
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（完全不重叠的分布）。

跨模拟器验证工作流

graph TD
    A[Build Test Circuit] --> B[Run on Stabilizer]
    A --> C[Run on CPUQVM]
    A --> D[Run on DensityMatrixSimulator]
    B --> E{Results consistent?}
    C --> E
    D --> E
    E -->|Yes| F["Circuit verified (ideal case)"]
    E -->|No| G["Debug: check gate compatibility"]
    F --> H[Add noise model]
    H --> I[Compare noisy CPUQVM vs DensityMatrix]
    I --> J[Analyze noise impact with exact probs]

---

完整示例：模拟器比较

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

完整比较代码

from pyqpanda3 import core

# 构建 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 结果

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])

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 结果

from pyqpanda3 import core
import numpy as np

# 构建 3 量子比特 GHZ 态（不含测量）
prog_no_meas = core.QProg()
prog_no_meas << core.H(0) << core.CNOT(0, 1) << core.CNOT(1, 2)

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.j 0. +0.j]
 [0. +0.j 0.5+0.j]]
  Tr(rho_0^2) = 0.5000
  => Maximally mixed (entangled with other qubits)

含噪声模拟结果

from pyqpanda3 import core
import numpy as np

# Build a 3-qubit GHZ state (no measurement)
prog_no_meas = core.QProg()
prog_no_meas << core.H(0) << core.CNOT(0, 1) << core.CNOT(1, 2)

print("\n--- DensityMatrixSimulator (Noisy) ---")
noise = core.NoiseModel()
# Apply depolarizing noise on idle qubits to model decoherence
noise.add_all_qubit_quantum_error(
    core.depolarizing_error(0.1), core.GateType.I
)
prog_noisy = core.QProg()
prog_noisy << core.H(0) << core.CNOT(0, 1) << core.CNOT(1, 2)
prog_noisy << core.I(0) << core.I(1) << core.I(2)
dm_sim2 = core.DensityMatrixSimulator()
dm_sim2.run(prog_noisy, noise)
noisy_probs = dm_sim2.state_probs()
print(f"Noisy probabilities: {[f'{p:.6f}' for p in noisy_probs]}")

# The noise leaks probability into previously zero states
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)")

# Quantify noise impact
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)")

With 10% depolarizing noise on idle qubits (modeled via I gates), the output shows:

--- DensityMatrixSimulator (Noisy) ---
Noisy probabilities: ['0.475417', '0.008194', '0.008194', '0.008194', '0.008194', '0.008194', '0.008194', '0.475417']
Purity: 0.860417 (< 1.0 indicates mixed state)
Probability leaked to other states: 0.049167

Stabilizer 结果

from pyqpanda3 import core

# 构建 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])

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.4754   | 0.4754   | 0.860   | 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 ⟩$$

---

MPSQVM —— 矩阵乘积态模拟

MPSQVM 用矩阵乘积态（Matrix Product State）后端运行量子程序。它不维护完整的 $2^n$ 维态矢量，而是把态表示为一组张量的收缩（MPS 形式），因此对纠缠较弱的线路在比特数较大时仍能保持较低内存占用——这正是 MPS 表示的优势所在。当线路纠缠增长（例如随机线路），MPS 的内部键维（bond dimension）会膨胀，内存优势随之消失，此时应改用 CPUQVM/GPUQVM。

MPSQVM.run(prog, shots, model) 与 CPUQVM 一样接受程序、采样次数和可选的 NoiseModel，运行后用 result() 取 QResult：

from pyqpanda3.core import MPSQVM, QProg, H, CNOT, measure, NoiseModel, depolarizing_error, GateType

# 4 比特 GHZ：H(0) + CNOT 阶梯
prog = QProg() << H(0)
for i in range(3):
    prog << CNOT(i, i + 1)
prog << measure(0, 0) << measure(1, 1) << measure(2, 2) << measure(3, 3)

# 无噪声 MPS 模拟：理想 GHZ 只有 0000/1111
qvm = MPSQVM()
qvm.run(prog, 2000)
ideal = qvm.result().get_counts()
print("only 0000/1111:", set(ideal.keys()) <= {"0000", "1111"})

# 含噪声 MPS 模拟：在 H 和 CNOT 上叠加去极化噪声
nm = NoiseModel()
nm.add_all_qubit_quantum_error(depolarizing_error(0.05), [GateType.H, GateType.CNOT])
noisy_qvm = MPSQVM()
noisy_qvm.run(prog, 2000, nm)
noisy = noisy_qvm.result().get_counts()
leaked = sum(v for k, v in noisy.items() if k not in ("0000", "1111"))
print("noisy leaked shots:", leaked, "of", sum(noisy.values()))

运行结果（计数随随机种子波动）：

only 0000/1111: True
noisy leaked shots: 294 of 2000

无噪声时 GHZ 态的测量结果严格落在 0000/1111 上（与 CPUQVM 一致）；叠加去极化噪声后，概率泄漏到其它比特串，泄漏比例反映了噪声强度。MPS 后端在物理上正确地传播了噪声对纠缠态的影响。

何时选 MPSQVM：线路纠缠较弱（如浅层 Clifford、GHZ 类、局部线路）且比特数较大时，MPS 的内存占用远低于完整态矢量；当线路接近最大纠缠（如随机量子线路）时，MPS 的键维增长会抵消这一优势，应改用态矢量模拟器。完整的 run/result 签名见 模拟器 API 参考。

---

SparseHeisenbergSimulator —— 无需采样的期望值计算

当你只需要泡利哈密顿量的期望值 <H>，且线路门集受支持、初态是直积态时，SparseHeisenbergSimulator 可以直接在海森堡绘景中计算期望值。它不是基于采样的模拟器，也不返回态矢量，因此针对这个特定问题远比完整态矢量模拟便宜。

from pyqpanda3.core import QProg, H, CNOT, SparseHeisenbergSimulator, sparse_heisenberg_expectation
from pyqpanda3.hamiltonian import Hamiltonian

zero = (1 + 0j, 0 + 0j)  # 每个比特为 |0>

# 制备贝尔态 (|00>+|11>)/sqrt(2)，再测量 <ZZ>：H(0) 把 qubit 0 置为 |+>，CNOT 关联两个比特
prog = QProg() << H(0) << CNOT(0, 1)

sim = SparseHeisenbergSimulator(num_qubits=2)
sim.set_initial_product_state([zero, zero])
sim.add_observable("ZZ", 1.0)
sim.evolve_backward(prog)

print(round(sim.expectation(), 6))

# 等价的一行写法（ZZ 也可写作带量子比特下标的 "Z0 Z1"）
ham = Hamiltonian({"Z0 Z1": 1.0})
print(round(sparse_heisenberg_expectation(prog, ham, initial_product_state=[zero, zero]), 6))

运行结果：

1.0
1.0

下面演示截断与 stats()：RX 旋转会把算符项数从 1 膨胀到 8，设 max_terms=4 后多余的项被按 l1 范数丢弃，stats() 给出丢弃项的统计。

from pyqpanda3.core import QProg, H, CNOT, RX, SparseHeisenbergSimulator

zero = (1 + 0j, 0 + 0j)
prog = QProg() << H(0) << CNOT(0, 1) << RX(0, 0.7) << RX(1, 1.3) << RX(2, 2.1)

# 不截断：观察算符项数膨胀到多少
full = SparseHeisenbergSimulator(num_qubits=3)
full.set_initial_product_state([zero, zero, zero])
full.add_observable("ZZZ", 1.0)
full.evolve_backward(prog)
print("no truncation terms:", full.stats().active_terms)

# 截断：最多保留 4 项，多余项按 l1 范数丢弃
trunc = SparseHeisenbergSimulator(num_qubits=3, max_terms=4, truncation_norm="l1")
trunc.set_initial_product_state([zero, zero, zero])
trunc.add_observable("ZZZ", 1.0)
trunc.evolve_backward(prog)
s = trunc.stats()
print("expectation:", round(trunc.expectation(), 6))
print("active_terms:", s.active_terms)
print("truncated_terms:", s.truncated_terms)
print("truncated_l1:", round(s.truncated_l1, 6))
print("coefficient_l1:", round(s.coefficient_l1, 6))

运行结果：

no truncation terms: 8
expectation: 0.400378
active_terms: 4
truncated_terms: 4
truncated_l1: 0.51565
coefficient_l1: 2.373074

当算符项数膨胀（例如经过非 Clifford 门）时，可用 max_terms 与 truncation_norm="l1" 截断，并通过 stats() 查看被丢弃的项（参见 API 参考）。

---

下一步

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