变分量子线路

变分量子线路（VQC）是参数化的量子线路，其量子门参数可以通过经典方法进行优化。它们构成了 VQE 和 QAOA 等变分量子算法的基础。

---

什么是变分线路？

变分量子线路由以下部分组成：

1. 参数化量子门：带有可调参数（如旋转角度）的量子门
2. 代价函数：需要最小化的可观测量（通常为哈密顿量）
3. 经典优化器：根据测量结果更新量子门参数

优化循环：

---

VQCircuit 模块

vqcircuit 模块提供了构建和优化参数化线路的工具：

from pyqpanda3 import vqcircuit

核心组件

组件	描述
VQCircuit	参数化量子线路构建器
Parameter	变分参数容器
DiffMethod	微分方法（ADJOINT_DIFF）
ResGradients	单参数集的梯度值
ResNGradients	N 个参数集的梯度值
ResGradientsAndExpectation	梯度与期望值的组合结果

---

VQCircuit 类 API

VQCircuit 类是构建、求值和微分参数化量子线路的核心接口。本节记录了每个公共方法及其使用示例。

构造函数

from pyqpanda3 import vqcircuit, core

# 创建 VQCircuit（pre_split 默认为 True）
vqc = vqcircuit.VQCircuit()

构造函数接受一个可选的 bool pre_split 参数（默认为 True），用于控制是否对线路进行预拆分以计算梯度。

添加量子门

变分量子门通过标准的 << 运算符追加到线路中，就像构建普通的 QProg 一样。不同之处在于，参数索引（类型为 MutableParamType）被提供，而不是数值角度。

from pyqpanda3 import core, vqcircuit

# 在 2 个量子比特上构建单层变分线路
# Param() 是 VQCircuit 上的方法，返回参数引用
vqc = vqcircuit.VQCircuit()
param_theta_0 = vqc.Param(0)   # 参数索引 0
param_theta_1 = vqc.Param(1)   # 参数索引 1
param_phi_0   = vqc.Param(2)   # 参数索引 2
param_phi_1   = vqc.Param(3)   # 参数索引 3

vqc << core.RY(0, param_theta_0)
vqc << core.RY(1, param_theta_1)
vqc << core.CNOT(0, 1)
vqc << core.RY(0, param_phi_0)
vqc << core.RY(1, param_phi_1)

设置参数值

在计算期望值或梯度之前，需要为参数槽分配具体的数值：

from pyqpanda3 import vqcircuit
import numpy as np

# 创建 VQCircuit 并定义参数（需先添加量子门才能设置参数）
vqc = vqcircuit.VQCircuit()
param_theta_0 = vqc.Param(0)
param_theta_1 = vqc.Param(1)
vqc << core.RY(0, param_theta_0)
vqc << core.RY(1, param_theta_1)
vqc << core.CNOT(0, 1)

# 通过 __call__ 一次性为所有参数赋值
theta = np.array([0.5, 1.2])
result = vqc(theta)  # 返回 VQCircuitResult

计算期望值

通过模拟工具获取哈密顿量相对于参数化态的期望值。使用 PauliOperator 作为可观测量：

from pyqpanda3 import core
from pyqpanda3.hamiltonian import PauliOperator

# 在 2 个量子比特上定义一个简单的 Z-Z 哈密顿量
hamiltonian = PauliOperator({"Z0 Z1": 1.0, "Z0": 0.5, "Z1": 0.3})

# 构建量子程序
prog = core.QProg()
prog << core.RY(0, 0.5) << core.RY(1, 1.2)
prog << core.CNOT(0, 1)
prog << core.RY(0, 0.8) << core.RY(1, 2.1)

# 计算期望值（态矢量后端）
exp_val = core.expval_pauli_operator(prog, hamiltonian)
print(f"Expectation value: {exp_val:.6f}")

计算梯度

使用伴随微分方法计算每个参数的梯度：

from pyqpanda3 import vqcircuit
from pyqpanda3.hamiltonian import PauliOperator
import numpy as np

# 使用伴随微分方法计算梯度
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

# 赋值具体参数值
theta = np.array([0.5, 1.2, 0.8, 2.1])

# 构建一个 2 量子比特变分线路用于梯度计算
vqc = vqcircuit.VQCircuit()
param_theta_0 = vqc.Param(0)
param_theta_1 = vqc.Param(1)
vqc << core.RY(0, param_theta_0)
vqc << core.RY(1, param_theta_1)
vqc << core.CNOT(0, 1)

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

# 梯度 API 是 VQCircuit 上的实例方法
grad_result = vqc.get_gradients(theta, hamiltonian, diff_method)

# 访问各个梯度值
for i in range(len(grad_result)):
    g = grad_result.at(i)
    print(f"  dE/dtheta[{i}] = {g:.6f}")

同时计算梯度和期望值

为了提高效率，梯度和期望值可以在一次计算中同时获得：

from pyqpanda3 import vqcircuit
import numpy as np

# 在一次求值中同时获得梯度和期望值
# 使用上面已定义的 vqc、theta、hamiltonian、diff_method
combined = vqc.get_gradients_and_expectation(theta, hamiltonian, diff_method)

print(f"Expectation: {combined.expectation_val():.6f}")
print(f"Gradients:   {combined.gradients()}")

使用 ResNGradients 进行批量求值

当需要同时对多组参数求梯度时，返回 ResNGradients 容器：

from pyqpanda3 import vqcircuit
from pyqpanda3.hamiltonian import PauliOperator
import numpy as np

# 对多组参数配置求梯度
# 将所有参数集展平为一维数组
batch_size = 3
theta_flat = np.array([
    0.1, 0.2, 0.3, 0.4,
    0.5, 0.6, 0.7, 0.8,
    1.0, 1.1, 1.2, 1.3,
])

# 使用带有 param_group_total 的批量重载
# 使用上面已定义的 vqc、hamiltonian、diff_method
n_grad_result = vqc.get_gradients(theta_flat, hamiltonian, batch_size, diff_method)

print(f"Number of parameter sets evaluated: {len(n_grad_result)}")
for i in range(len(n_grad_result)):
    grad_set = n_grad_result.at(i)
    print(f"  Set {i} gradients: {grad_set.gradients()}")

API 方法汇总

方法	签名	返回类型	描述
构造函数	VQCircuit(pre_split=True)	VQCircuit	创建线路（pre_split 控制梯度预计算）
添加量子门	vqc << gate	VQCircuit	向线路追加一个变分量子门
参数引用	vqc.Param(idx)	MutableParamType	获取用于量子门构建的参数引用
梯度	vqc.get_gradients(theta, H, method)	ResGradients	计算 $⟨H⟩$ 的梯度
批量梯度	vqc.get_gradients(theta, H, n, method)	ResNGradients	n 组参数的梯度
组合结果	vqc.get_gradients_and_expectation(theta, H, method)	ResGradientsAndExpectation	同时获取梯度和期望值

---

变分量子门（VQGate）

变分量子门是参数化的量子门封装器。它们接受 MutableParamType 对象而不是固定角度，允许在优化过程中更新参数。

当使用 MutableParamType 参数调用时，VQGate 函数从 core 模块导出：

单量子比特参数化量子门

量子门	签名	描述
RX(qubit, param)	core.RX(idx, MutableParamType)	带参数的 RX 旋转
RY(qubit, param)	core.RY(idx, MutableParamType)	带参数的 RY 旋转
RZ(qubit, param)	core.RZ(idx, MutableParamType)	带参数的 RZ 旋转
P(qubit, param)	core.P(idx, MutableParamType)	带参数的相位门
U1(qubit, param)	core.U1(idx, MutableParamType)	带参数的 U1 门
U2(qubit, p1, p2)	core.U2(idx, param1, param2)	带 2 个参数的 U2
U3(qubit, p1, p2, p3)	core.U3(idx, p1, p2, p3)	带 3 个参数的 U3
U4(qubit, p1, p2, p3, p4)	core.U4(idx, p1, p2, p3, p4)	带 4 个参数的 U4
RPHI(qubit, p1, p2)	core.RPhi(idx, param1, param2)	带 2 个参数的 RPhi

双量子比特参数化量子门

量子门	签名	描述
CU(ctrl, targ, p1-p4)	core.CU(ctrl, targ, ...)	带 4 个参数的受控 U 门
CP(ctrl, targ, param)	core.CP(ctrl, targ, param)	受控相位门
CRX(ctrl, targ, param)	core.CRX(ctrl, targ, param)	受控 RX 门
CRY(ctrl, targ, param)	core.CRY(ctrl, targ, param)	受控 RY 门
CRZ(ctrl, targ, param)	core.CRZ(ctrl, targ, param)	受控 RZ 门
RXX(ctrl, targ, param)	core.RXX(ctrl, targ, param)	RXX 旋转
RYY(ctrl, targ, param)	core.RYY(ctrl, targ, param)	RYY 旋转
RZZ(ctrl, targ, param)	core.RZZ(ctrl, targ, param)	RZZ 旋转
RZX(ctrl, targ, param)	core.RZX(ctrl, targ, param)	RZX 旋转

混合固定参数与变分参数

多参数量子门支持将固定值与变分参数混合使用：

from pyqpanda3 import core, vqcircuit

# 创建 VQCircuit 并获取变分参数引用
vqc = vqcircuit.VQCircuit()
param_index = vqc.Param(0)

# U3，前两个参数固定，第三个为变分参数
gate = core.U3(0, 1.57, 0.0, param_index)

# U4，仅最后一个参数为变分参数
gate = core.U4(0, 0.0, 0.0, 0.0, param_index)

# U2，第一个参数为变分参数，第二个固定
gate = core.U2(0, param_index, 3.14)

---

参数管理

Parameter 类

Parameter 类管理变分参数：

from pyqpanda3 import vqcircuit

# Parameter 类通常通过 VQCircuit 使用
# VQCircuit 内部管理参数，mutable_parameter_total() 返回参数总数
vqc = vqcircuit.VQCircuit()

# 查询线路中的变分参数总数（初始为 0）
n = vqc.mutable_parameter_total()
print(f"Number of mutable parameters: {n}")

MutableParamType

MutableParamType 指定参数值的来源：

• 基于索引：param = vqc.Param([idx_dim0, idx_dim1, ...]) — 引用多维参数
• 标量索引：param = vqc.Param(idx) — 通过索引引用单个参数

from pyqpanda3 import core, vqcircuit

# 创建 VQCircuit 来持有参数
vqc = vqcircuit.VQCircuit()

# 使用 VQCircuit.Param() 创建参数化的 RX 门
# MutableParamType 指定使用哪个参数槽
param_0 = vqc.Param(0)  # 引用索引 0 处的参数
param_1 = vqc.Param(1)  # 引用索引 1 处的参数

# 将变分量子门追加到 VQCircuit
vqc << core.RX(0, param_0)
vqc << core.RY(1, param_1)

---

梯度计算

伴随微分

伴随微分方法通过先正向再反向（伴随）演化线路来高效计算梯度。这避免了参数偏移规则所需的 $2n$ 次线路求值。

参考文献：Jones 和 Gacon，"Efficient calculation of gradients in classical simulations of variational quantum algorithms"，arXiv:2009.02823。

from pyqpanda3 import vqcircuit

# 指定微分方法
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

ResGradients - 单参数集

ResGradients 保存单组参数的梯度值：

# 假设 gradients 是来自 VQCircuit 计算的 ResGradients 对象
gradients = vqc_result.gradients()  # 获取 ResGradients

# 获取梯度数量
n_grads = len(gradients)

# 获取所有梯度值列表
all_grads = gradients.gradients()

# 获取索引 i 处参数的梯度
grad_i = gradients.at(i)

# 获取多维参数的梯度
grad_md = gradients.at([idx_dim0, idx_dim1])

# 字符串表示
print(gradients)

ResNGradients - N 个参数集

ResNGradients 保存 N 组参数的梯度值（批量求值）：

from pyqpanda3 import vqcircuit
from pyqpanda3.hamiltonian import PauliOperator
import numpy as np

# 假设 n_gradients 是一个 ResNGradients 对象
n_gradients = vqc.get_gradients(theta_flat, hamiltonian, batch_size, diff_method)

# 获取参数集数量（N）
n_sets = len(n_gradients)

# 获取第 i 组参数的 ResGradients
grad_set_i = n_gradients.at(i)

# 获取所有梯度数据（嵌套列表）
all_data = n_gradients.data()

# 字符串表示
print(n_gradients)

ResGradientsAndExpectation - 组合结果

ResGradientsAndExpectation 同时保存梯度值和期望值：

# 假设 result 是一个 ResGradientsAndExpectation 对象
result = vqc.get_gradients_and_expectation(theta, hamiltonian, diff_method)

# 获取期望值
exp_val = result.expectation_val()

# 获取所有梯度值列表
grads = result.gradients()

# 获取特定参数的梯度
grad_i = result.at([idx])

# 获取组合数据：(期望值, [梯度])
data = result.data()

# 字符串表示
print(result)

---

数学基础

参数化量子态

将参数化线路 $U(\theta )$ 作用于 $|0{⟩}^{\otimes n}$ 产生：

$$|\psi (\theta )⟩=U(\theta )|0{⟩}^{\otimes n}$$

期望值

代价函数是哈密顿量 $H$ 的期望值：

$$⟨H{⟩}_{\theta }=⟨\psi (\theta )|H|\psi (\theta )⟩$$

伴随微分梯度

关于参数 ${\theta }_k$ 的梯度为：

$$\frac{\partial ⟨H⟩}{\partial {\theta }_k}=2\text{Re}[⟨{\psi }_L|{\psi }_R⟩]$$

其中：

$$|{\psi }_L⟩=\frac{\partial U(\theta )}{\partial {\theta }_k}|0⟩,|{\psi }_R⟩=U^{†}(\theta )H|\psi (\theta )⟩$$

伴随方法通过在单次正向-反向演化过程中计算内积，避免了存储中间态。

参数偏移规则（替代方法）

作为对比，参数偏移规则通过两次线路求值来计算每个梯度：

$$\frac{\partial ⟨H⟩}{\partial {\theta }_k}=\frac{1}{2}[⟨H{⟩}_{{\theta }_k+\pi /2}-⟨H{⟩}_{{\theta }_k-\pi /2}]$$

对于 $n$ 个参数，这需要 $2n$ 次线路求值，而伴随微分仅需一次正向-反向演化。

---

变分量子算法工作流

---

逐步构建 VQCircuit

本节将从零开始完整构建一个变分量子线路，并逐步解释每个决策。我们将构建一个适用于小型优化问题的 3 量子比特线路。

第一步：确定拟设结构

在编写任何代码之前，首先确定适合您问题的拟设（线路结构）。本教程使用硬件高效拟设，包含交替的旋转层和纠缠层。

第二步：导入模块和定义常量

from pyqpanda3 import core, vqcircuit
import numpy as np

# 线路配置
N_QUBITS = 3       # 量子比特数
DEPTH = 2          # 重复层数
# 每层每个量子比特有 2 个旋转（先 RY 后 RZ）
N_PARAMS = N_QUBITS * 2 * DEPTH  # 共 12 个参数

第三步：创建参数对象

Parameter 对象跟踪所有变分参数及其当前值。

# 初始化参数容器
params = vqcircuit.Parameter()
print(f"Initial parameter count: {params.size()}")

第四步：构建拟设线路

逐层构建线路。每层由每个量子比特上的 RY 和 RZ 旋转组成，随后是 CNOT 纠缠门。

def build_hardware_efficient_ansatz(n_qubits, depth):
    """
    Construct a hardware-efficient ansatz with the given qubit count and depth.
    Returns a QProg and the total number of parameters used.
    """
    prog = core.QProg()
    param_idx = 0  # 参数分配的运行索引

    for layer in range(depth):
        # --- 旋转子层 ---
        for q in range(n_qubits):
            # 为量子比特 q 分配 RY 参数（使用 0.0 作为占位符）
            prog << core.RY(q, 0.0)
            param_idx += 1

            # 为量子比特 q 分配 RZ 参数（使用 0.0 作为占位符）
            prog << core.RZ(q, 0.0)
            param_idx += 1

        # --- 纠缠子层（线性最近邻） ---
        for q in range(n_qubits - 1):
            prog << core.CNOT(q, q + 1)

        # 闭合环路以实现全连接
        if n_qubits > 2:
            prog << core.CNOT(n_qubits - 1, 0)

    total_params = param_idx
    return prog, total_params

# 构建线路
ansatz_prog, n_params_used = build_hardware_efficient_ansatz(N_QUBITS, DEPTH)
print(f"Circuit built with {n_params_used} parameters")

第五步：定义可观测量（哈密顿量）

选择一个哈密顿量作为代价函数。演示中使用横向场 Ising 模型：

$$H=-\sum _{⟨i,j⟩}{Z}_i{Z}_j-h\sum _i{X}_i$$

from pyqpanda3.hamiltonian import PauliOperator

# 将哈密顿量构建为 PauliOperator
hamiltonian_terms = {}

# ZZ 耦合项（环上的最近邻）
for i in range(N_QUBITS):
    j = (i + 1) % N_QUBITS
    hamiltonian_terms[f"Z{i} Z{j}"] = -1.0

# 横向场项
h_field = 0.5
for i in range(N_QUBITS):
    hamiltonian_terms[f"X{i}"] = -h_field

hamiltonian = PauliOperator(hamiltonian_terms)
print(f"Hamiltonian terms: {hamiltonian_terms}")

第六步：初始化参数

设置初始参数值。常见的选择是在 $[0,2\pi )$ 范围内均匀分布的随机角度：

# 随机初始化
np.random.seed(42)
initial_theta = np.random.uniform(0, 2 * np.pi, n_params_used)
print(f"Initial parameters: {initial_theta}")

第七步：评估代价函数

def cost_function(theta_array):
    """Compute the expectation value of the Hamiltonian for given parameters."""
    # 使用给定参数值构建新线路
    prog = core.QProg()
    idx = 0
    for layer in range(DEPTH):
        for q in range(N_QUBITS):
            prog << core.RY(q, theta_array[idx]); idx += 1
            prog << core.RZ(q, theta_array[idx]); idx += 1
        for q in range(N_QUBITS - 1):
            prog << core.CNOT(q, q + 1)
        if N_QUBITS > 2:
            prog << core.CNOT(N_QUBITS - 1, 0)

    # 计算完整哈密顿量的期望值
    return core.expval_pauli_operator(prog, hamiltonian)

# 测试代价函数
energy_0 = cost_function(initial_theta)
print(f"Initial energy: {energy_0:.6f}")

第八步：运行优化循环

from scipy.optimize import minimize

# 使用 COBYLA（无导数）进行优化
result = minimize(
    cost_function,
    initial_theta,
    method='COBYLA',
    options={'maxiter': 300, 'rhobeg': 0.5}
)

print(f"Optimized energy: {result.fun:.6f}")
print(f"Optimal parameters: {result.x}")
print(f"Converged: {result.success}")
print(f"Function evaluations: {result.nfev}")

第九步：验证结果

优化完成后，检查最终量子态并验证收敛性：

# 在最优参数处重建线路
optimal_prog = core.QProg()
idx = 0
for layer in range(DEPTH):
    for q in range(N_QUBITS):
        optimal_prog << core.RY(q, result.x[idx]); idx += 1
        optimal_prog << core.RZ(q, result.x[idx]); idx += 1
    for q in range(N_QUBITS - 1):
        optimal_prog << core.CNOT(q, q + 1)
    if N_QUBITS > 2:
        optimal_prog << core.CNOT(N_QUBITS - 1, 0)

# 获取最终态矢量
machine = core.CPUQVM()
machine.run(optimal_prog)
final_state = machine.result().get_state_vector()
print(f"State vector norm: {np.linalg.norm(final_state):.6f}")

# 计算测量概率
probs = np.abs(final_state) ** 2
for i, p in enumerate(probs):
    if p > 0.01:
        bitstring = format(i, f'0{N_QUBITS}b')
        print(f"  |{bitstring}>: {p:.4f}")

---

批量参数求值

在探索参数空间或执行基于种群的优化时，同时对多组参数求梯度是高效的。ResNGradients 容器就是为此设计的。

为什么需要批量求值？

单独对 $N$ 组参数求梯度需要 $N$ 次独立的模拟运行，每次都有自己的线路构建开销。批量求值可以摊销设置成本，特别适用于：

• 在离散参数网格上进行网格搜索
• 维护候选解种群的进化算法
• 并行梯度比较以确定最有前景的下降方向

批量求值工作流

示例：评估参数网格

from pyqpanda3 import core, vqcircuit
from pyqpanda3.hamiltonian import PauliOperator
import numpy as np

# --- 构建一个简单的 2 量子比特变分线路 ---
def build_2q_circuit(theta):
    """Build a 2-qubit circuit with 2 variational parameters."""
    prog = core.QProg()
    prog << core.RY(0, theta[0])
    prog << core.RY(1, theta[1])
    prog << core.CNOT(0, 1)
    prog << core.RY(0, theta[0])
    return prog

# 定义哈密顿量（Z-Z 相互作用）
H = PauliOperator({"Z0 Z1": 1.0, "Z0": 0.5})

# --- 创建参数值网格 ---
grid_size = 5
theta_0_vals = np.linspace(0, np.pi, grid_size)
theta_1_vals = np.linspace(0, np.pi, grid_size)

# 收集所有参数组合
batch = []
for t0 in theta_0_vals:
    for t1 in theta_1_vals:
        batch.append([t0, t1])

print(f"Total parameter sets: {len(batch)}")

计算批量梯度

# 计算每组参数的梯度
# 使用伴随微分方法
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

gradient_results = []
energies = []

for theta in batch:
    prog = build_2q_circuit(theta)

    # 计算期望值
    energy = core.expval_pauli_operator(prog, H)
    energies.append(energy)

    # 为演示目的，手动计算参数偏移梯度
    shift = np.pi / 2
    grads = []
    for i in range(2):
        # 正向偏移
        theta_plus = theta.copy()
        theta_plus[i] += shift
        prog_plus = build_2q_circuit(theta_plus)
        e_plus = core.expval_pauli_operator(prog_plus, H)

        # 反向偏移
        theta_minus = theta.copy()
        theta_minus[i] -= shift
        prog_minus = build_2q_circuit(theta_minus)
        e_minus = core.expval_pauli_operator(prog_minus, H)

        # 通过参数偏移规则计算梯度
        grad_i = 0.5 * (e_plus - e_minus)
        grads.append(grad_i)

    gradient_results.append(grads)

# 找到下降最陡的参数组
grad_norms = [np.linalg.norm(g) for g in gradient_results]
best_idx = np.argmin(energies)

print(f"Best energy: {energies[best_idx]:.6f} at params {batch[best_idx]}")
print(f"Gradient at best: {gradient_results[best_idx]}")
print(f"Grad norm: {grad_norms[best_idx]:.6f}")

分析能量景观

# 将结果重塑为二维网格以进行可视化分析
energy_grid = np.array(energies).reshape(grid_size, grid_size)

# 报告能量景观统计信息
print(f"Energy range: [{np.min(energy_grid):.6f}, {np.max(energy_grid):.6f}]")
print(f"Energy mean:  {np.mean(energy_grid):.6f}")
print(f"Energy std:   {np.std(energy_grid):.6f}")

# 识别局部极小值
flat_idx = np.argmin(energy_grid)
min_row, min_col = np.unravel_index(flat_idx, energy_grid.shape)
print(f"Grid minimum at theta_0={theta_0_vals[min_col]:.3f}, "
      f"theta_1={theta_1_vals[min_row]:.3f}")
print(f"Minimum energy: {energy_grid[min_row, min_col]:.6f}")

---

梯度验证与调试

当变分算法未能收敛时，梯度计算往往是问题所在。本节描述了验证梯度正确性和调试常见问题的策略。

数值梯度检查

最可靠的验证方法是将解析梯度与有限差分近似进行比较：

$$\frac{\partial f}{\partial {\theta }_k}\approx \frac{f({\theta }_k+ϵ)-f({\theta }_k-ϵ)}{2ϵ}$$

from pyqpanda3 import core, vqcircuit
from pyqpanda3.hamiltonian import PauliOperator
import numpy as np

def numerical_gradient(cost_fn, theta, epsilon=1e-5):
    """
    Compute numerical gradients via central finite differences.

    Args:
        cost_fn: Callable that takes a parameter array and returns a float.
        theta: Current parameter values (numpy array).
        epsilon: Finite difference step size.

    Returns:
        numpy array of gradient values, one per parameter.
    """
    n = len(theta)
    grads = np.zeros(n)

    for k in range(n):
        # 正向扰动参数 k
        theta_plus = theta.copy()
        theta_plus[k] += epsilon
        f_plus = cost_fn(theta_plus)

        # 反向扰动参数 k
        theta_minus = theta.copy()
        theta_minus[k] -= epsilon
        f_minus = cost_fn(theta_minus)

        # 中心差分
        grads[k] = (f_plus - f_minus) / (2 * epsilon)

    return grads


# 示例用法：将伴随梯度与数值梯度进行比较
def build_demo_circuit(theta):
    prog = core.QProg()
    prog << core.RY(0, theta[0])
    prog << core.RZ(0, theta[1])
    prog << core.RY(1, theta[2])
    prog << core.CNOT(0, 1)
    return prog

H_demo = PauliOperator({"Z0 Z1": 1.0, "X0": 0.3})

def demo_cost(theta):
    prog = build_demo_circuit(theta)
    return core.expval_pauli_operator(prog, H_demo)

# 在特定参数点进行测试
theta_test = np.array([0.5, 1.0, 1.5])
num_grads = numerical_gradient(demo_cost, theta_test)

print("Numerical gradients:", num_grads)
print("Gradient magnitudes:", np.abs(num_grads))

梯度一致性检查

进行更稳健的检查，在多个参数点比较梯度：

def verify_gradients_at_point(theta, cost_fn, epsilon=1e-5, tol=1e-4):
    """
    Verify that numerical gradients are self-consistent
    at multiple epsilon values.
    """
    # 在两个不同的 epsilon 值下计算
    g1 = numerical_gradient(cost_fn, theta, epsilon=epsilon)
    g2 = numerical_gradient(cost_fn, theta, epsilon=epsilon * 0.1)

    # 它们应该在容差范围内一致
    diff = np.max(np.abs(g1 - g2))
    consistent = diff < tol

    print(f"  Epsilon={epsilon:.1e} vs {epsilon*0.1:.1e}: max diff = {diff:.2e}")
    print(f"  Consistent: {consistent}")
    return consistent

# 在多个随机参数点进行验证
np.random.seed(0)
for trial in range(5):
    theta_random = np.random.uniform(0, 2 * np.pi, 3)
    print(f"Trial {trial}: theta = {theta_random}")
    verify_gradients_at_point(theta_random, demo_cost)

常见调试检查清单

当梯度出现异常时，请检查以下事项：

检测贫瘠高原

贫瘠高原（barren plateau）发生在梯度幅值随系统规模呈指数衰减时，导致优化无法进行。通过监控梯度统计信息来检测：

def check_barren_plateau(cost_fn, n_params, n_samples=100):
    """
    Estimate gradient variance across random parameter initializations.
    Exponentially small variance indicates a barren plateau.
    """
    grad_samples = np.zeros((n_samples, n_params))

    for s in range(n_samples):
        theta_rand = np.random.uniform(0, 2 * np.pi, n_params)
        grads = numerical_gradient(cost_fn, theta_rand)
        grad_samples[s] = grads

    # 计算每个梯度分量的方差
    grad_var = np.var(grad_samples, axis=0)
    mean_var = np.mean(grad_var)

    print(f"Mean gradient variance: {mean_var:.6e}")
    print(f"Per-parameter variances: {grad_var}")

    # 启发式判断：如果方差低于 1e-8，则可能存在贫瘠高原
    if mean_var < 1e-8:
        print("WARNING: Possible barren plateau detected!")
        print("  Consider: reducing circuit depth, using local cost functions,")
        print("  or initializing with identity-block structure.")
    else:
        print("Gradient variance appears healthy.")

    return mean_var

# 运行检查
check_barren_plateau(demo_cost, n_params=3)

优化过程中记录梯度范数

在优化过程中监控梯度范数有助于识别收敛问题：

class GradientMonitor:
    """Track gradient norms across optimization iterations."""

    def __init__(self):
        self.grad_norms = []
        self.energies = []

    def record(self, energy, gradients):
        """Record energy and gradient norm for this iteration."""
        self.energies.append(energy)
        norm = np.linalg.norm(gradients)
        self.grad_norms.append(norm)

    def report(self):
        """Print a summary of gradient behavior."""
        norms = np.array(self.grad_norms)
        print(f"Total iterations: {len(norms)}")
        print(f"Gradient norm range: [{np.min(norms):.6e}, {np.max(norms):.6e}]")
        print(f"Final gradient norm: {norms[-1]:.6e}")

        # 检查梯度是否在减小
        if len(norms) > 10:
            early = np.mean(norms[:10])
            late = np.mean(norms[-10:])
            print(f"Early avg norm: {early:.6e}, Late avg norm: {late:.6e}")
            if late < early:
                print("Gradient norm is decreasing (good convergence).")
            else:
                print("WARNING: Gradient norm not decreasing!")

# 在优化过程中使用
monitor = GradientMonitor()

def monitored_cost(theta):
    energy = demo_cost(theta)
    grads = numerical_gradient(demo_cost, theta)
    monitor.record(energy, grads)
    return energy

# 带监控地运行
result = minimize(monitored_cost, np.array([0.5, 1.0, 1.5]),
                  method='COBYLA', options={'maxiter': 50})
monitor.report()

---

性能优化

本节提供了在 pyqpanda3 中提高变分量子线路计算性能的实用建议。

选择合适的微分方法

方法	每个梯度的代价	精度	后端要求
伴随方法	$O(n_{\text{gates}})$	机器精度	态矢量
参数偏移	$2n_{\text{params}}$ 次线路求值	机器精度	任意后端
有限差分	$2n_{\text{params}}$ 次线路求值	取决于 $ϵ$	任意后端

建议：在使用态矢量模拟器时，始终使用伴随微分。它在单次正向-反向演化中计算所有 $n$ 个梯度。

# 推荐：伴随微分（单次通过）
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

减少线路深度

更浅的线路带来更快的模拟速度和更小的数值误差：

# 与其使用深层的硬件高效拟设（多层）
# 不如使用问题启发式的最小深度拟设

# 不佳：通用的深层线路
def deep_ansatz(n_qubits, depth=20):
    prog = core.QProg()
    for _ in range(depth):
        for q in range(n_qubits):
            prog << core.RY(q, 0.0)  # 占位符
        for q in range(n_qubits - 1):
            prog << core.CNOT(q, q + 1)
    return prog

# 推荐：问题定制的浅层线路（例如用于 H₂）
def shallow_h2_ansatz(theta):
    """Minimal 2-qubit UCCSD-inspired ansatz for H2."""
    prog = core.QProg()
    prog << core.X(1)                  # Initial Hartree-Fock state
    prog << core.RY(0, theta[0])
    prog << core.CNOT(0, 1)
    prog << core.RY(0, theta[1])
    return prog

利用哈密顿量对称性

将可交换的泡利项分组，以减少线路求值次数：

from pyqpanda3.hamiltonian import PauliOperator

# 在测量具有多个泡利项的哈密顿量时，
# 将可以在相同泡利基下同时测量的项分组

# 例如，Z0 Z1 和 Z0 都可以在 Z 基下测量
h_grouped = {
    "Z_basis": {"Z0 Z1": -0.01128010, "Z0": 0.39793742, "Z1": -0.39793742},
    "X_basis": {"X0 X1": 0.18093119},
    "identity": {"": -1.0523792},
}

def efficient_cost(theta, prog_fn, grouped_h):
    """Compute energy using grouped Pauli measurements."""
    prog = prog_fn(theta)
    total = 0.0

    for basis, terms in grouped_h.items():
        if basis == "identity":
            for pauli_str, coeff in terms.items():
                total += coeff
            continue

        # 为该基组构建组合的 PauliOperator
        pauli_dict = {k: v for k, v in terms.items()}
        op = PauliOperator(pauli_dict)
        exp_val = core.expval_pauli_operator(prog, op)
        total += exp_val

    return total

优化器选择指南

# 无导数（稳健，易于使用）
from scipy.optimize import minimize
result_cobyla = minimize(cost_function, theta0, method='COBYLA',
                         options={'maxiter': 200})

# 基于导数（当梯度准确时收敛更快）
result_lbfgsb = minimize(cost_function, theta0, method='L-BFGS-B',
                         jac=gradient_function,
                         options={'maxiter': 200})

参数初始化策略

良好的初始化可以大幅减少优化迭代次数：

import numpy as np

# 策略 1：均匀随机
theta_rand = np.random.uniform(0, 2 * np.pi, n_params)

# 策略 2：零附近的小随机扰动
theta_small = np.random.normal(0, 0.1, n_params)

# 策略 3：恒等块初始化（避免贫瘠高原）
# 初始化使 U(theta_0) 接近恒等操作
theta_identity = np.zeros(n_params)

# 策略 4：基于问题对称性的结构化初始化
# 对于 QAOA，gamma 初始化为接近 0.5，beta 初始化为接近 0.25
theta_qaoa = np.concatenate([
    np.full(p, 0.5),    # gamma 参数
    np.full(p, 0.25),   # beta 参数
])

缓存中间结果

当代价函数使用相似参数反复调用时，缓存可以节省计算时间：

from functools import lru_cache
from pyqpanda3.hamiltonian import PauliOperator

# 缓存泡利项的期望值（当线路是确定性的时）
# 注意：仅在参数可哈希时有效（元组，而非数组）

@lru_cache(maxsize=1024)
def cached_expectation(theta_tuple, pauli_str):
    """Cache expectation values keyed by parameter values and Pauli string."""
    theta = np.array(theta_tuple)
    prog = build_demo_circuit(theta)
    op = PauliOperator({pauli_str: 1.0})
    return core.expval_pauli_operator(prog, op)

# 在调用前将 theta 转换为可哈希的元组
theta_tuple = tuple(theta_test)
exp_zz = cached_expectation(theta_tuple, "Z0 Z1")

---

完整示例：H₂ 的 VQE

from pyqpanda3 import core, vqcircuit
from pyqpanda3.hamiltonian import PauliOperator
import numpy as np

# 第一步：定义 H₂ 哈密顿量（简化版）
# H = c0*I + c1*Z0 + c2*Z1 + c3*Z0*Z1 + c4*X0*X1 + c5*Y0*Y1
# 键长 0.735 Å 的系数
h2_coeffs = {
    "": -1.0523792,
    "Z0": 0.39793742,
    "Z1": -0.39793742,
    "Z0 Z1": -0.01128010,
    "X0 X1": 0.18093119,
}

# 第二步：定义拟设线路（硬件高效）
def build_ansatz(params_array):
    """Build a parameterized ansatz for H₂."""
    prog = core.QProg()

    # 参数化单量子比特旋转
    prog << core.RY(0, params_array[0])
    prog << core.RY(1, params_array[1])

    # 纠缠层
    prog << core.CNOT(0, 1)

    # 第二旋转层
    prog << core.RY(0, params_array[2])
    prog << core.RY(1, params_array[3])

    # 纠缠层
    prog << core.CNOT(0, 1)

    return prog

# 第三步：定义代价函数
def cost_function(params_array):
    """Compute expectation value of H₂ Hamiltonian."""
    prog = build_ansatz(params_array)
    total = 0.0

    for pauli_str, coeff in h2_coeffs.items():
        if pauli_str == "":
            total += coeff
            continue

        # 为该项创建泡利算符
        op = PauliOperator({pauli_str: 1.0})
        exp_val = core.expval_pauli_operator(prog, op, shots=1000)
        total += coeff * exp_val

    return total

# 第四步：使用 scipy 进行优化
from scipy.optimize import minimize

# 初始参数
theta0 = np.random.uniform(0, 2 * np.pi, 4)

result = minimize(cost_function, theta0, method='COBYLA',
                  options={'maxiter': 200, 'rhobeg': 0.5})

print(f"Optimal energy: {result.fun:.6f}")
print(f"Optimal parameters: {result.x}")
print(f"Expected ground state energy: -1.137275")

解读 VQE 结果

VQE 优化的输出提供了几个需要仔细解读的信息。

能量精度

${\text{H}}_2$ 在键长 0.735 A 处的精确基态能量约为 $-1.137275$ Hartree。VQE 结果应与此值进行比较：

$$\mathrm{\Delta }E=E_{\text{VQE}}-E_{\text{exact}}$$

# 将 VQE 结果与精确值进行比较
exact_energy = -1.137275
vqe_energy = result.fun
error = abs(vqe_energy - exact_energy)

print(f"VQE energy:    {vqe_energy:.6f} Hartree")
print(f"Exact energy:  {exact_energy:.6f} Hartree")
print(f"Absolute error: {error:.6f} Hartree")
print(f"Chemical accuracy threshold: 0.0016 Hartree (1 kcal/mol)")

if error < 0.0016:
    print("Result is within chemical accuracy!")
else:
    print("Result exceeds chemical accuracy; try more iterations or a better ansatz.")

最优参数分析

优化后的参数揭示了基态线路的结构：

# 检查最优参数
optimal_theta = result.x
print("Optimal parameters (radians):")
for i, t in enumerate(optimal_theta):
    print(f"  theta[{i}] = {t:.6f}  ({np.degrees(t):.2f} degrees)")

# 检查是否有参数接近 0 或 pi（表明该门实际上等价于恒等操作或泡利门）
for i, t in enumerate(optimal_theta):
    t_mod = t % (2 * np.pi)
    if t_mod < 0.1 or abs(t_mod - np.pi) < 0.1:
        print(f"  theta[{i}] is near a critical value (0 or pi)")

态组成

重建最优量子态并检查哪些计算基态有贡献：

# 在最优参数处构建线路
optimal_prog = build_ansatz(optimal_theta)

# 获取态矢量
machine = core.CPUQVM()
machine.run(optimal_prog)
state = machine.result().get_state_vector()

# 打印基态上的概率分布
print("\nBasis state probabilities:")
for i, amp in enumerate(state):
    prob = abs(amp) ** 2
    if prob > 0.001:
        label = format(i, '02b')  # 2 位二进制字符串
        # 映射量子比特排序（bit 0 = 最右边）
        print(f"  |{label}> : amplitude = {amp:+.6f}, probability = {prob:.6f}")

# 对于 H₂，基态约为：
# |psi> ~ 0.9939|00> - 0.1102|11>
# 这对应于带有小双激发的 Hartree-Fock 态

收敛行为

分析优化器的收敛轨迹：

# 在每次函数求值时记录能量
energy_history = []

def tracked_cost(params_array):
    energy = cost_function(params_array)
    energy_history.append(energy)
    return energy

# 带追踪地重新运行优化
theta0_retry = np.random.uniform(0, 2 * np.pi, 4)
tracked_result = minimize(tracked_cost, theta0_retry, method='COBYLA',
                          options={'maxiter': 200, 'rhobeg': 0.5})

# 报告收敛统计信息
print(f"Total function evaluations: {len(energy_history)}")
print(f"Initial energy: {energy_history[0]:.6f}")
print(f"Final energy:   {energy_history[-1]:.6f}")
print(f"Energy improvement: {energy_history[0] - energy_history[-1]:.6f}")

# 检查能量是否已稳定（最后 10 次求值）
last_10 = energy_history[-10:]
energy_std = np.std(last_10)
print(f"Energy std (last 10 evals): {energy_std:.6e}")
if energy_std < 1e-4:
    print("Energy has stabilized (good convergence).")
else:
    print("Energy is still fluctuating; consider more iterations.")

扩展到更大分子

${\text{H}}_2$ 示例使用了 2 个量子比特和 4 个参数。对于更大的分子，资源需求会增长：

分子	量子比特数	参数数（UCCSD）	泡利项数
${\text{H}}_2$ (STO-3G)	2	4	5
LiH (STO-3G)	4	~20	~100
${\text{BeH}}_2$ (STO-3G)	6	~60	~200
${\text{H}}_2\text{O}$ (STO-3G)	8	~150	~600

---

变分线路构建模式

硬件高效拟设

def hardware_efficient_ansatz(n_qubits, depth, params):
    """Build a hardware-efficient ansatz."""
    prog = core.QProg()
    idx = 0

    for d in range(depth):
        # 单量子比特旋转层
        for q in range(n_qubits):
            prog << core.RY(q, params[idx])
            idx += 1
            prog << core.RZ(q, params[idx])
            idx += 1

        # 纠缠层（线性连接）
        for q in range(n_qubits - 1):
            prog << core.CNOT(q, q + 1)

    return prog

QAOA 拟设

def qaoa_ansatz(n_qubits, p, params, cost_hamiltonian):
    """Build a QAOA ansatz with p layers."""
    gamma = params[:p]      # Cost parameters
    beta = params[p:2*p]    # Mixer parameters

    prog = core.QProg()

    # Initial superposition
    for q in range(n_qubits):
        prog << core.H(q)

    for layer in range(p):
        # Cost layer (apply e^{-iγH_C})
        # For Ising-type Hamiltonian, use RZZ gates
        prog << core.RZZ(0, 1, 2 * gamma[layer])
        # Add more problem-specific gates

        # Mixer layer (apply e^{-iβH_M})
        for q in range(n_qubits):
            prog << core.RX(q, 2 * beta[layer])

    return prog

---

DiffMethod 参考

ADJOINT_DIFF

伴随微分方法：

• 参考文献：Jones 和 Gacon (2020)
• 复杂度：每次梯度计算 $O(n_{\text{gates}})$（对比参数偏移规则的 $O(n_{\text{params}})$ 次线路求值）
• 精度：机器精度（无统计噪声）
• 要求：态矢量模拟后端

---

梯度结果 API 汇总

ResGradients

方法	返回	描述
.at(idx)	float	索引 idx 处参数的梯度
.at([dims])	float	多维参数的梯度
.gradients()	list[float]	所有梯度值
.__len__()	int	梯度数量
.__str__()	str	字符串表示

ResNGradients

方法	返回	描述
.at(idx)	ResGradients	第 idx 组参数的梯度
.data()	list[list[float]]	所有梯度数据
.__len__()	int	参数集数量（N）
.__str__()	str	字符串表示

ResGradientsAndExpectation

方法	返回	描述
.expectation_val()	float	期望值
.gradients()	list[float]	所有梯度值
.at([dims])	float	特定参数的梯度
.data()	tuple[float, list[float]]	（期望值，[梯度]）
.__str__()	str	字符串表示

---

QAOA 完整示例：最大割

量子近似优化算法（QAOA）是一种专为组合优化问题设计的变分算法。本节演示了将 QAOA 应用于小图上的最大割问题。

问题定义

给定无向图 $G=(V,E)$，最大割要求将顶点 $V=V_0\cup V_1$ 划分为两部分，使得跨越割的边数最大化：

$$\text{MaxCut}(G)=\underset{x\in \{0,1{\}}^n}{max}\sum _{(i,j)\in E}\mathbb{1}[x_i\ne x_j]$$

这可以编码为最小化代价哈密顿量：

$$H_C=-\sum _{(i,j)\in E}\frac{1}{2}(1-Z_i{Z}_j)$$

负号确保最小化 $H_C$ 对应于最大化割的大小。

示例图

我们使用一个 4 顶点环图，边为 $(0,1),(1,2),(2,3),(3,0)$：

该图的最优最大割值为 4（通过划分 {0,2} vs {1,3} 可以切断所有边）。

QAOA 线路结构

具有 $p$ 层的 QAOA 线路在代价酉算符和混合酉算符之间交替：

$$|\gamma ,\beta ⟩=\prod _{l=1}^p{e}^{-i{\beta }_l{H}_M}{e}^{-i{\gamma }_l{H}_C}|+{⟩}^{\otimes n}$$

其中 $H_M=\sum _i{X}_i$ 是混合哈密顿量。

第一步：定义问题图

from pyqpanda3 import core, vqcircuit
from pyqpanda3.hamiltonian import PauliOperator
import numpy as np

# 定义 4 顶点环图
edges = [(0, 1), (1, 2), (2, 3), (3, 0)]
n_qubits = 4
n_edges = len(edges)

print(f"Graph: {n_qubits} vertices, {n_edges} edges")
print(f"Edges: {edges}")
print(f"Optimal MaxCut value: {n_edges}")

第二步：构建代价哈密顿量

代价哈密顿量编码了最大割目标。每条边贡献一个 $-Z_i{Z}_j$ 项（去掉常数）：

# 从图的边构建代价哈密顿量
# H_C = sum_{(i,j) in E} Z_i Z_j  （我们最小化这个量）
cost_terms = {}
for i, j in edges:
    cost_terms[f"Z{i} Z{j}"] = 1.0

cost_hamiltonian = PauliOperator(cost_terms)
print(f"Cost Hamiltonian terms: {cost_terms}")

第三步：构建 QAOA 拟设

def build_qaoa_circuit(n_qubits, edges, gamma, beta):
    """
    Build a QAOA circuit with parameters gamma (cost) and beta (mixer).

    Args:
        n_qubits: Number of qubits (vertices in the graph).
        edges: List of (i, j) tuples defining the graph edges.
        gamma: List of cost parameters, one per QAOA layer.
        beta: List of mixer parameters, one per QAOA layer.

    Returns:
        A QProg containing the QAOA circuit.
    """
    p = len(gamma)  # QAOA 层数
    assert len(beta) == p, "gamma and beta must have the same length"

    prog = core.QProg()

    # 步骤 A：制备均匀叠加态
    for q in range(n_qubits):
        prog << core.H(q)

    # 步骤 B：应用 p 个交替层
    for layer in range(p):
        # 代价酉算符：e^{-i * gamma * H_C}
        # 对于每条边 (i, j)，应用 RZZ(i, j, 2 * gamma)
        for i, j in edges:
            prog << core.RZZ(i, j, 2 * gamma[layer])

        # 混合酉算符：e^{-i * beta * H_M}
        # 对于每个量子比特，应用 RX(q, 2 * beta)
        for q in range(n_qubits):
            prog << core.RX(q, 2 * beta[layer])

    return prog

# 使用单层进行测试
p = 1
gamma_init = [0.5]
beta_init = [0.25]

test_circuit = build_qaoa_circuit(n_qubits, edges, gamma_init, beta_init)
print("QAOA circuit built successfully")

第四步：定义代价函数

def qaoa_cost(params, n_qubits, edges, p):
    """
    Compute the QAOA cost function (expectation of H_C).

    Args:
        params: Flat array of [gamma_0, ..., gamma_{p-1}, beta_0, ..., beta_{p-1}].
        n_qubits: Number of qubits.
        edges: Graph edges.
        p: Number of QAOA layers.

    Returns:
        Expectation value of the cost Hamiltonian.
    """
    gamma = params[:p].tolist()
    beta = params[p:2*p].tolist()

    prog = build_qaoa_circuit(n_qubits, edges, gamma, beta)

    # 计算代价哈密顿量的期望值
    expectation = core.expval_pauli_operator(prog, cost_hamiltonian)

    return expectation

# 在初始参数处求值
params_init = np.array(gamma_init + beta_init)
cost_0 = qaoa_cost(params_init, n_qubits, edges, p)
print(f"Initial cost: {cost_0:.6f}")

第五步：优化 QAOA 参数

from scipy.optimize import minimize

# 使用多次随机重启运行优化
best_cost = float('inf')
best_params = None
n_restarts = 10

for restart in range(n_restarts):
    # 随机初始参数
    theta0 = np.random.uniform(0, np.pi, 2 * p)

    # 使用 COBYLA 进行最小化
    result = minimize(
        qaoa_cost,
        theta0,
        args=(n_qubits, edges, p),
        method='COBYLA',
        options={'maxiter': 300, 'rhobeg': 0.5}
    )

    if result.fun < best_cost:
        best_cost = result.fun
        best_params = result.x.copy()

    print(f"  Restart {restart}: cost = {result.fun:.6f}, "
          f"converged = {result.success}")

print(f"\nBest cost found: {best_cost:.6f}")
print(f"Best gamma: {best_params[:p]}")
print(f"Best beta:  {best_params[p:]}")

第六步：采样最优解

# 在最优参数处构建线路
optimal_gamma = best_params[:p].tolist()
optimal_beta = best_params[p:].tolist()
optimal_prog = build_qaoa_circuit(n_qubits, edges, optimal_gamma, optimal_beta)

# 从线路中采样以获取候选解
n_shots = 1000
machine = core.CPUQVM()
machine.run(optimal_prog, shots=n_shots)
counts = machine.result().get_counts()

# 评估每个测量比特串的割值
print("\nTop solutions by frequency:")
sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)

for bitstring, count in sorted_counts[:5]:
    # 计算该比特串的割值
    cut_value = 0
    for i, j in edges:
        if bitstring[i] != bitstring[j]:
            cut_value += 1

    probability = count / n_shots
    print(f"  {bitstring}: cut={cut_value}, "
          f"count={count}, probability={probability:.3f}")

# 找到整体最优比特串
best_cut = 0
best_bitstring = None
for bitstring, count in counts.items():
    cut_value = 0
    for i, j in edges:
        if bitstring[i] != bitstring[j]:
            cut_value += 1
    if cut_value > best_cut:
        best_cut = cut_value
        best_bitstring = bitstring

print(f"\nBest solution found: {best_bitstring} with cut value {best_cut}")
print(f"Optimal cut value: {n_edges}")
print(f"Approximation ratio: {best_cut / n_edges:.4f}")

第七步：增加 QAOA 深度

更高的 $p$ 值以更多参数为代价改善近似比：

# 使用 p=2 层运行 QAOA 以获得更好的近似
p = 2

best_cost_p2 = float('inf')
best_params_p2 = None

for restart in range(10):
    theta0 = np.random.uniform(0, np.pi, 2 * p)

    result = minimize(
        qaoa_cost,
        theta0,
        args=(n_qubits, edges, p),
        method='COBYLA',
        options={'maxiter': 300, 'rhobeg': 0.5}
    )

    if result.fun < best_cost_p2:
        best_cost_p2 = result.fun
        best_params_p2 = result.x.copy()

print(f"\np=1 best cost: {best_cost:.6f}")
print(f"p=2 best cost: {best_cost_p2:.6f}")
print(f"Expected optimal: {-n_edges:.1f} (all edges cut)")

QAOA 近似比

近似比衡量 QAOA 结果与最优解的接近程度：

$$\alpha =\frac{⟨H_C{⟩}_{\text{QAOA}}}{\text{MaxCut}(G)}$$

# 计算近似比
optimal_cut = n_edges  # 对于 4 顶点环图，最优割 = 4

# 代价哈密顿量是 Z_i Z_j 的求和
# 当所有边都被切断时，每个 Z_i Z_j = -1，因此 cost = -n_edges
# 近似比：|cost| / n_edges
alpha_p1 = abs(best_cost) / optimal_cut
alpha_p2 = abs(best_cost_p2) / optimal_cut

print(f"Approximation ratio (p=1): {alpha_p1:.4f}")
print(f"Approximation ratio (p=2): {alpha_p2:.4f}")

# 理论保证：对于 3 正则图，p=1 的 QAOA 达到 alpha >= 0.6924
# (Farhi et al. 2014)
print(f"\nTheoretical p=1 lower bound (3-regular): 0.6924")

扩展到更大的图

同样的框架适用于更大的图。以下是任意图定义的辅助函数：

from pyqpanda3.hamiltonian import PauliOperator

def solve_maxcut_qaoa(edges, n_qubits, p=1, n_restarts=20):
    """
    Solve MaxCut using QAOA for an arbitrary graph.

    Args:
        edges: List of (i, j) tuples defining the graph.
        n_qubits: Number of vertices.
        p: QAOA depth parameter.
        n_restarts: Number of random restarts for optimization.

    Returns:
        Dictionary with optimization results.
    """
    # 构建代价哈密顿量
    cost_terms = {f"Z{i} Z{j}": 1.0 for i, j in edges}
    cost_H = PauliOperator(cost_terms)

    best_cost = float('inf')
    best_params = None

    for restart in range(n_restarts):
        theta0 = np.random.uniform(0, np.pi, 2 * p)

        def cost_fn(params):
            gamma = params[:p].tolist()
            beta = params[p:2*p].tolist()
            prog = build_qaoa_circuit(n_qubits, edges, gamma, beta)
            return core.expval_pauli_operator(prog, cost_H)

        result = minimize(cost_fn, theta0, method='COBYLA',
                          options={'maxiter': 300})

        if result.fun < best_cost:
            best_cost = result.fun
            best_params = result.x.copy()

    # 采样最优线路
    gamma = best_params[:p].tolist()
    beta = best_params[p:].tolist()
    optimal_prog = build_qaoa_circuit(n_qubits, edges, gamma, beta)
    machine = core.CPUQVM()
    machine.run(optimal_prog, shots=2000)
    counts = machine.result().get_counts()

    # 找到最优比特串
    best_cut = 0
    for bitstring in counts:
        cut = sum(1 for i, j in edges if bitstring[i] != bitstring[j])
        best_cut = max(best_cut, cut)

    return {
        'best_cost': best_cost,
        'best_params': best_params,
        'best_cut': best_cut,
        'optimal_cut': len(edges),
        'approximation_ratio': best_cut / len(edges),
        'counts': counts,
    }

# 示例：求解 5 顶点路径图
path_edges = [(0, 1), (1, 2), (2, 3), (3, 4)]
result_path = solve_maxcut_qaoa(path_edges, n_qubits=5, p=2)
print(f"Path graph result: cut={result_path['best_cut']}/{result_path['optimal_cut']}, "
      f"ratio={result_path['approximation_ratio']:.4f}")

---

下一步

• 哈密顿量与泡利算符 - 定义 VQE 代价函数的可观测量
• 模拟 - 在模拟器上运行变分线路
• 量子信息 - 分析变分态的性质