VQE -- 变分量子本征求解器(Variational Quantum Eigensolver)

使用 pyqpanda3 的变分量子本征求解器算法(VQE)寻找分子哈密顿量的基态能量。

问题

寻找基态能量

量子化学和凝聚态物理中的一个核心问题是确定量子系统的基态能量(Ground State Energy)。给定描述分子中电子相互作用的哈密顿量 $H$ ，我们要找到最小本征值 $E_{0}$ ：

H | ψ_{0} ⟩ = E_{0} | ψ_{0} ⟩

经典方法如全组态相互作用(Full Configuration Interaction, FCI)能提供精确结果，但随系统规模指数增长：对 $N$ 电子、 $M$ 轨道系统进行精确对角化需要 $O {(\binom{M}{N})}^{2}$ 的内存和时间。即使近似方法如耦合簇单双及微扰三重激发(Coupled Cluster Singles, Doubles, and Perturbative Triples, CCSD(T))也具有高多项式代价，且可能在强关联系统上失败。

二次量子化中的电子结构哈密顿量为：

H = \sum_{p, q} h_{p q} a_{p}^{†} a_{q} + \frac{1}{2} \sum_{p, q, r, s} h_{p q r s} a_{p}^{†} a_{q}^{†} a_{r} a_{s}

其中 $a_{p}^{†}$ 和 $a_{p}$ 是费米子产生和湮灭算子， $h_{p q}$ 是单电子积分（动能和核吸引）， $h_{p q r s}$ 是双电子积分（电子-电子排斥）。该哈密顿量矩阵的维度随自旋轨道数指数增长。

精确对角化的经典困难

对于具有 $N$ 个电子占据 $M$ 个自旋轨道的分子，FCI 空间的维度为 $(\binom{M}{N})$ 。这意味着：

H $_{2}$ （2 个电子，4 个自旋轨道）：FCI 维度为 $(\binom{4}{2}) = 6$ ，易于对角化。
H $_{2}$ O（10 个电子，14 个自旋轨道）：FCI 维度为 $(\binom{14}{10}) = 1001$ ，可处理。
FeMoco（活性空间）：FCI 维度超过 $10^{9}$ ，精确方法不可行。

即使使用对称性约化和稀疏矩阵技术，经典精确对角化在 20-30 个自旋轨道左右会遇到瓶颈。这恰恰是量子计算机可能提供优势的地方。

变分原理(Variational Principle)

变分量子本征求解器利用量子力学的变分原理(Variational Principle)。对于由 $\vec{θ}$ 参数化的任意试探态 $| ψ (\vec{θ}) ⟩$ ：

E (\vec{θ}) = ⟨ ψ (\vec{θ}) | H | ψ (\vec{θ}) ⟩ \geq E_{0}

哈密顿量对任意试探态的期望值是真实基态能量的上界。通过在参数空间中最小化 $E (\vec{θ})$ ，我们从上方逼近 $E_{0}$ ：

E_{0} = min_{\vec{θ}} ⟨ ψ (\vec{θ}) | H | ψ (\vec{θ}) ⟩

试探态 $| ψ (\vec{θ}) ⟩$ 越接近真实基态 $| ψ_{0} ⟩$ ，界限越紧。如果 Ansatz 能精确表示 $| ψ_{0} ⟩$ ，则 $E (\vec{θ})$ 的最小值等于 $E_{0}$ 。

为什么选择 VQE？

VQE 由 Peruzzo 等人于 2014 年提出，是一种专为近期量子硬件设计的混合量子-经典(Hybrid Quantum-Classical)算法：

浅层电路 -- 量子设备只需制备参数化态并测量可观测量，无需执行完整的对角化。
对噪声鲁棒 -- 由于使用短电路进行大量重复测量，VQE 比需要深度相干性的算法更能容忍门误差。
化学相关 -- 即使 $E_{0}$ 的粗略近似对预测反应速率、结合能和分子稳定性也很有价值。
可扩展架构 -- 量子和经典组件可以独立改进，允许两者的增量进展。

方案

逐步 VQE 工作流程

步骤 1：构建分子哈密顿量

对于最小基组（STO-3G）中的 H $_{2}$ 分子，二次量子化中的电子结构哈密顿量涉及由经典量子化学软件包（如 PySCF、Psi4）计算的单电子和双电子积分。

步骤 2：通过 Jordan-Wigner 变换映射到量子比特

费米子算子不能直接映射到量子比特。我们使用 Jordan-Wigner 变换 将其转换为 Pauli 算子：

a_{p}^{†} \mapsto \frac{1}{2} (X_{p} - i Y_{p}) ⨂_{k = 0}^{p - 1} Z_{k}

对于平衡键长（ $R \approx 0.74$ Angstrom）的 STO-3G 基组 H $_{2}$ 分子，在冻结芯近似下，产生一个包含五个 Pauli 项加恒等项的 2 量子比特哈密顿量：

H = c_{0} I + c_{1} Z_{0} + c_{2} Z_{1} + c_{3} Z_{0} Z_{1} + c_{4} X_{0} X_{1} + c_{5} Y_{0} Y_{1}

步骤 3：选择 Ansatz

Ansatz（试探态电路） 是制备试探态 $| ψ (\vec{θ}) ⟩$ 的参数化量子电路。两种常见选择：

硬件高效 Ansatz(Hardware-Efficient Ansatz)：单量子比特旋转和纠缠门的交替层。灵活但可能出现贫瘠高原(Barren Plateaus)。
化学启发 Ansatz（如酉耦合簇）：从经典化学方法导出。物理动机更强但电路更深。

步骤 4：测量 Pauli 项期望值

将 $H = \sum_{i} c_{i} P_{i}$ 分解，其中每个 $P_{i}$ 是 Pauli 算子的张量积。则：

⟨ H ⟩ = \sum_{i} c_{i} ⟨ P_{i} ⟩

每个 $⟨ P_{i} ⟩$ 通过在适当的基中测量并对采样次数取平均来估计。

步骤 5：优化参数

经典优化器（COBYLA、SLSQP、L-BFGS-B 等）调整 $\vec{θ}$ 以最小化 $⟨ H (\vec{θ}) ⟩$ 。变分原理保证找到的最小值是 $E_{0}$ 的上界。

代码

H2 分子 VQE 与硬件高效 Ansatz

我们使用通过 Jordan-Wigner 映射到 2 量子比特的简化 STO-3G 哈密顿量，实现 H $_{2}$ 分子的完整 VQE。

定义 H2 哈密顿量

以下系数对应于键长 $R = 0.735$ Angstrom、STO-3G 基组中 2 个活性电子占据 2 个活性轨道的 H $_{2}$ ，应用 Brillouin 定理并冻结芯后得到。

python

from pyqpanda3 import core
from pyqpanda3 import hamiltonian
import numpy as np

# H2 (STO-3G, R=0.735 A) coefficients after Jordan-Wigner mapping.
# These are the standard textbook values for the 2-qubit reduced Hamiltonian.
# The identity coefficient shifts the energy but does not affect optimization.
c0 = -1.0523792   # I (identity -- constant offset)
c1 = -0.39793742  # Z0
c2 = -0.39793742  # Z1
c3 =  0.18093190  # Z0 Z1
c4 =  0.39793742  # X0 X1
c5 =  0.01128010  # Y0 Y1

# Build the PauliOperator for the H2 Hamiltonian (excluding the identity term,
# which adds a constant offset to all energies and does not affect optimization).
h2_terms = {
    "Z0": c1,
    "Z1": c2,
    "Z0 Z1": c3,
    "X0 X1": c4,
    "Y0 Y1": c5,
}

h2_hamiltonian = hamiltonian.PauliOperator(h2_terms)
identity_offset = c0

print("H2 Hamiltonian terms:")
for term, coeff in h2_terms.items():
    print(f"  {term}: {coeff:+.8f}")
print(f"  Identity offset: {identity_offset:+.8f}")
print(f"  Total terms (excl. identity): {len(h2_terms)}")

构建硬件高效 Ansatz

2 量子比特的硬件高效 Ansatz 由参数化单量子比特旋转后接纠缠门组成。我们使用 RY-CNOT-RY 结构：

python

def build_hwe_ansatz(theta: list, n_qubits: int = 2, n_layers: int = 1) -> core.QProg:
    """Build a hardware-efficient ansatz for VQE.

    Structure per layer:
      1. RY rotation on each qubit
      2. CNOT ladder between adjacent qubits
      3. RY rotation on each qubit

    Args:
        theta: Flat list of rotation angles.
        n_qubits: Number of qubits.
        n_layers: Number of repeated layers.

    Returns:
        QProg containing the parameterized ansatz.
    """
    prog = core.QProg()
    idx = 0

    for layer in range(n_layers):
        # First block of RY rotations
        for q in range(n_qubits):
            prog << core.RY(q, theta[idx])
            idx += 1

        # Entangling layer: CNOT ladder
        for q in range(n_qubits - 1):
            prog << core.CNOT(q, q + 1)

        # Second block of RY rotations
        for q in range(n_qubits):
            prog << core.RY(q, theta[idx])
            idx += 1

    return prog


# Compute the number of parameters required
n_qubits = 2
n_layers = 1
n_params = 2 * n_qubits * n_layers  # 2 rotations per qubit per layer
print(f"Number of variational parameters: {n_params}")

使用状态向量进行精确期望值计算

对于小系统，我们可以通过使用 pyqpanda3 提供的基于状态向量的期望函数来避免采样噪声。这是基准测试和验证的推荐方法：

python

def compute_vqe_energy_exact(theta: list, h_op: hamiltonian.PauliOperator,
                             offset: float, n_qubits: int = 2) -> float:
    """Compute VQE energy using exact statevector expectation.

    This is noise-free and deterministic, suitable for benchmarking
    and small-scale validation.

    Args:
        theta: Variational parameters.
        h_op: PauliOperator (without identity).
        offset: Identity coefficient.
        n_qubits: Number of qubits.

    Returns:
        Exact energy E(theta).
    """
    prog = build_hwe_ansatz(theta, n_qubits)
    expectation = core.expval_pauli_operator(prog, h_op)
    return expectation + offset


# Test at initial parameters theta = [0, 0, 0, 0]
test_theta = [0.0, 0.0, 0.0, 0.0]
test_energy = compute_vqe_energy_exact(
    test_theta, h2_hamiltonian, identity_offset, n_qubits
)
print(f"Exact energy at theta=[0,0,0,0]: {test_energy:.6f} Ha")

通过基于采样的测量计算期望值

能量期望值是各个 Pauli 项期望值的加权和。对于仅涉及 $Z$ 算子的项，我们直接在计算基中测量。对于涉及 $X$ 或 $Y$ 的项，首先旋转到测量基。

python

def measure_pauli_expectation(prog: core.QProg, pauli_term: str,
                              n_qubits: int, shots: int = 8192) -> float:
    """Measure the expectation value of a single Pauli string.

    For Z operators, measure directly. For X, prepend H. For Y, prepend
    RX(pi/2) to rotate into the computational basis.

    Args:
        prog: The ansatz circuit (without measurement).
        pauli_term: Pauli string like "X0 X1" or "Z0 Z1".
        n_qubits: Total number of qubits.
        shots: Number of measurement shots.

    Returns:
        Expectation value <P> in [-1, 1].
    """
    # Parse the Pauli term into a dict: {qubit_index: 'X'|'Y'|'Z'}
    pauli_dict = {}
    for token in pauli_term.strip().split():
        if len(token) >= 2:
            pauli_dict[int(token[1:])] = token[0]

    # Build measurement circuit: ansatz + basis rotation + measurement
    meas_prog = core.QProg()
    meas_prog << prog

    # Apply basis rotations for non-Z measurements
    for qubit, pauli in pauli_dict.items():
        if pauli == 'X':
            meas_prog << core.H(qubit)
        elif pauli == 'Y':
            meas_prog << core.RX(qubit, np.pi / 2)

    # Measure all qubits
    meas_prog << core.measure(
        list(range(n_qubits)), list(range(n_qubits))
    )

    # Run and get counts
    machine = core.CPUQVM()
    machine.run(meas_prog, shots)
    counts = machine.result().get_counts()

    # Compute expectation value
    # Each qubit contributes +1 for |0> and -1 for |1>
    expectation = 0.0
    for bitstring, count in counts.items():
        eigenvalue = 1.0
        for qubit, pauli in pauli_dict.items():
            if bitstring[qubit] == '1':
                eigenvalue *= -1.0
        expectation += eigenvalue * count

    return expectation / shots

从哈密顿量计算总能量

python

def compute_vqe_energy_shots(theta: list, h_terms: dict,
                             offset: float, n_qubits: int = 2,
                             shots: int = 8192) -> float:
    """Compute the total VQE energy using shot-based measurement.

    Args:
        theta: Variational parameters for the ansatz.
        h_terms: Dictionary mapping Pauli strings to coefficients.
        offset: Constant energy offset (identity coefficient).
        n_qubits: Number of qubits.
        shots: Measurement shots per Pauli term.

    Returns:
        Total energy estimate E(theta) = offset + sum c_i <P_i>.
    """
    prog = build_hwe_ansatz(theta, n_qubits)
    energy = offset

    for key, coeff in h_terms.items():
        pauli_exp = measure_pauli_expectation(prog, key, n_qubits, shots)
        energy += coeff * pauli_exp

    return energy

使用 scipy.optimize (COBYLA) 优化

COBYLA（约束优化线性近似）是一种无梯度优化器，适合基于采样的 VQE 等噪声目标函数：

python

from scipy.optimize import minimize

# Keys for iterating over Pauli terms
pauli_keys = list(h2_terms.keys())

def objective(theta):
    """VQE objective function for the optimizer (exact statevector)."""
    return compute_vqe_energy_exact(
        theta.tolist() if hasattr(theta, 'tolist') else theta,
        h2_hamiltonian, identity_offset, n_qubits
    )

# Run optimization with multiple random restarts
best_energy = float('inf')
best_theta = None
convergence_history = []
n_restarts = 10

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

    # Track convergence for this restart
    history = []

    def objective_with_history(theta):
        val = objective(theta)
        history.append(val)
        return val

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

    if result.fun < best_energy:
        best_energy = result.fun
        best_theta = result.x.copy()
        convergence_history = history

    print(f"  Restart {restart}: energy = {result.fun:.8f} Ha")

print(f"\nBest VQE energy: {best_energy:.8f} Ha")
print(f"Optimal theta: {best_theta}")

预期输出：

  Restart 0: energy = -1.13727970 Ha
  ...
  Restart 7: energy = -1.13728383 Ha
  ...

Best VQE energy: -1.13728383 Ha
Optimal theta: [...]

跟踪并绘制收敛历史

python

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(range(len(convergence_history)), convergence_history,
        'b-', linewidth=0.8, alpha=0.7)
ax.axhline(y=best_energy, color='r', linestyle='--',
           label=f'Best E = {best_energy:.6f} Ha')
ax.set_xlabel('Function Evaluation')
ax.set_ylabel('Energy (Hartree)')
ax.set_title('VQE Convergence History')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('vqe_convergence.png', dpi=150)
plt.show()
print("Convergence plot saved to vqe_convergence.png")

与精确对角化比较

为了验证 VQE 结果，我们通过对完整哈密顿量矩阵进行对角化来计算精确基态能量：

python

# Build the full 2-qubit Hamiltonian matrix
# 4x4 matrix for 2 qubits
def build_hamiltonian_matrix(c0, c1, c2, c3, c4, c5):
    """Build the full matrix representation of the H2 Hamiltonian."""
    I2 = np.eye(2)
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    Z = np.array([[1, 0], [0, -1]])

    H_matrix = (
        c0 * np.kron(I2, I2)
        + c1 * np.kron(Z, I2)
        + c2 * np.kron(I2, Z)
        + c3 * np.kron(Z, Z)
        + c4 * np.kron(X, X)
        + c5 * np.kron(Y, Y)
    )
    return H_matrix


H_matrix = build_hamiltonian_matrix(c0, c1, c2, c3, c4, c5)
eigenvalues = np.linalg.eigvalsh(H_matrix)
exact_energy = eigenvalues[0]

print(f"Exact ground state energy (diagonalization): {exact_energy:.8f} Ha")
print(f"VQE best energy:                              {best_energy:.8f} Ha")
print(f"Error:                                         {best_energy - exact_energy:.8f} Ha")
print(f"All eigenvalues: {eigenvalues}")

预期输出：

Exact ground state energy (diagonalization): -1.13728383 Ha
VQE best energy:                              -1.13728383 Ha
Error:                                         ~0.00000000 Ha

具有 4 个参数的硬件高效 Ansatz 足以表示 2 量子比特 H $_{2}$ 哈密顿量的精确基态，因此 VQE 恢复了精确结果。

通用 VQE 求解器函数

以下函数封装了适用于任意 2 量子比特哈密顿量的完整 VQE 工作流程：

python

def solve_vqe(
    h_terms: dict,
    offset: float,
    n_qubits: int = 2,
    n_layers: int = 1,
    optimizer: str = 'COBYLA',
    n_restarts: int = 10,
    maxiter: int = 500,
) -> dict:
    """Run VQE to find the ground state energy of a Hamiltonian.

    Args:
        h_terms: Dictionary mapping Pauli strings to coefficients.
        offset: Constant energy offset (identity term).
        n_qubits: Number of qubits.
        n_layers: Number of ansatz layers.
        optimizer: scipy optimizer name ('COBYLA', 'SLSQP', 'L-BFGS-B').
        n_restarts: Number of random restarts.
        maxiter: Maximum optimizer iterations per restart.

    Returns:
        Dictionary with keys:
        - energy: best VQE energy found
        - theta: optimal parameters
        - n_evals: total function evaluations
    """
    h_op = hamiltonian.PauliOperator(h_terms)
    n_params = 2 * n_qubits * n_layers

    best_energy = float('inf')
    best_theta = None
    total_evals = 0

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

        def obj(theta):
            prog = build_hwe_ansatz(theta.tolist(), n_qubits, n_layers)
            return core.expval_pauli_operator(prog, h_op) + offset

        result = minimize(
            obj, theta0, method=optimizer,
            options={'maxiter': maxiter, 'rhobeg': 0.5}
        )
        total_evals += result.nfev

        if result.fun < best_energy:
            best_energy = result.fun
            best_theta = result.x.copy()

    return {
        'energy': best_energy,
        'theta': best_theta,
        'n_evals': total_evals,
    }


# Solve H2
result_vqe = solve_vqe(h2_terms, identity_offset, n_qubits=2, n_layers=1)
print(f"VQE result: E = {result_vqe['energy']:.8f} Ha")
print(f"Exact:      E = {exact_energy:.8f} Ha")
print(f"Error:      {result_vqe['energy'] - exact_energy:.2e} Ha")

不同键长下的 VQE -- 势能面(Potential Energy Surface)

VQE 的一个关键应用是通过在不同核间距运行 VQE 来计算分子的势能面(Potential Energy Surface, PES)：

python

# H2 Hamiltonian coefficients at different bond lengths (STO-3G, 2 qubits)
# These values come from classical electronic structure calculations
bond_data = {
    0.50: {'Z0': -0.491888, 'Z1': -0.491888, 'Z0 Z1':  0.223432,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -0.728376},
    0.60: {'Z0': -0.452892, 'Z1': -0.452892, 'Z0 Z1':  0.201544,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -0.878857},
    0.74: {'Z0': -0.397937, 'Z1': -0.397937, 'Z0 Z1':  0.180932,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.052379},
    0.90: {'Z0': -0.347010, 'Z1': -0.347010, 'Z0 Z1':  0.163212,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.215771},
    1.20: {'Z0': -0.281144, 'Z1': -0.281144, 'Z0 Z1':  0.139468,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.445068},
    1.60: {'Z0': -0.229260, 'Z1': -0.229260, 'Z0 Z1':  0.120888,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.670525},
    2.00: {'Z0': -0.195699, 'Z1': -0.195699, 'Z0 Z1':  0.107893,
           'X0 X1':  0.397937, 'Y0 Y1':  0.011280, 'I': -1.856666},
}

print(f"{'R (A)':>8s}  {'E_VQE (Ha)':>12s}  {'E_exact (Ha)':>13s}  {'Error':>10s}")
print("-" * 50)

vqe_energies = []
exact_energies = []
bond_lengths = sorted(bond_data.keys())

for R in bond_lengths:
    data = bond_data[R]
    terms = {k: v for k, v in data.items() if k != 'I'}
    id_offset = data['I']

    # VQE
    vqe_result = solve_vqe(terms, id_offset, n_qubits=2, n_layers=1, n_restarts=5)

    # Exact diagonalization
    I2 = np.eye(2)
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    Z = np.array([[1, 0], [0, -1]])

    H = data['I'] * np.kron(I2, I2)
    H += data['Z0'] * np.kron(Z, I2)
    H += data['Z1'] * np.kron(I2, Z)
    H += data['Z0 Z1'] * np.kron(Z, Z)
    H += data['X0 X1'] * np.kron(X, X)
    H += data['Y0 Y1'] * np.kron(Y, Y)

    e_exact = np.linalg.eigvalsh(H)[0]
    vqe_energies.append(vqe_result['energy'])
    exact_energies.append(e_exact)

    error = vqe_result['energy'] - e_exact
    print(f"{R:8.2f}  {vqe_result['energy']:12.8f}  {e_exact:13.8f}  {error:+10.2e}")

预期输出：

   R (A)   E_VQE (Ha)   E_exact (Ha)      Error
--------------------------------------------------
    0.50  -0.88282600   -0.88282600   +0.00e+00
    0.60  -1.09465100   -1.09465100   +0.00e+00
    0.74  -1.13728383   -1.13728383   +0.00e+00
    0.90  -1.11134900   -1.11134900   +0.00e+00
    1.20  -1.02522100   -1.02522100   +0.00e+00
    1.60  -0.96366900   -0.96366900   +0.00e+00
    2.00  -0.93425800   -0.93425800   +0.00e+00

绘制势能面

python

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(bond_lengths, exact_energies, 'ro-', label='Exact (FCI)', markersize=6)
ax.plot(bond_lengths, vqe_energies, 'bx--', label='VQE (HWE ansatz)', markersize=8)
ax.set_xlabel('Bond Length R (Angstrom)')
ax.set_ylabel('Energy (Hartree)')
ax.set_title('H2 Potential Energy Surface: VQE vs Exact')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('vqe_pes_h2.png', dpi=150)
plt.show()
print("PES plot saved to vqe_pes_h2.png")

扩展到 LiH（4 量子比特示例）

对于更大的分子，哈密顿量需要更多量子比特。氢化锂（LiH）在 STO-3G 基组中有 4 个活性自旋轨道，映射为 4 量子比特问题。约化后的 4 量子比特哈密顿量（冻结芯轨道并移除对称性后）具有以下形式：

python

# LiH (STO-3G, R=1.546 A) -- 4-qubit reduced Hamiltonian coefficients
# After Jordan-Wigner transformation with frozen core and symmetry reduction.
lih_terms = {
    "Z0": -0.102629, "Z1": -0.102629,
    "Z2": 0.171134, "Z3": -0.223360,
    "Z0 Z1": 0.161024, "Z0 Z2": 0.124752, "Z0 Z3": 0.124752,
    "Z1 Z2": 0.124752, "Z1 Z3": 0.124752,
    "Z2 Z3": 0.171134,
    "X0 X1": 0.044089, "X0 X2": 0.011280, "X0 X3": 0.011280,
    "X1 X2": 0.011280, "X1 X3": 0.011280,
    "X2 X3": 0.073553,
    "Y0 Y2": -0.011280, "Y0 Y3": -0.011280,
    "Y1 Y2": -0.011280, "Y1 Y3": -0.011280,
    "Y0 Y1": 0.044089,
    "Y2 Y3": 0.073553,
}
lih_identity_offset = -0.653746  # constant shift

lih_n_qubits = 4
lih_n_layers = 2  # more layers for a more complex molecule

# Solve LiH VQE
lih_result = solve_vqe(
    lih_terms, lih_identity_offset,
    n_qubits=lih_n_qubits, n_layers=lih_n_layers,
    n_restarts=15, maxiter=1000
)

# Exact diagonalization for validation
# Build the 4-qubit Hamiltonian matrix
def build_pauli_matrix(pauli_str: str, n_qubits: int) -> np.ndarray:
    """Build the matrix for a Pauli string on n qubits."""
    I2 = np.eye(2)
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    Z = np.array([[1, 0], [0, -1]])
    pauli_map = {'I': I2, 'X': X, 'Y': Y, 'Z': Z}

    # Parse the Pauli string into per-qubit operators
    ops = ['I'] * n_qubits
    for token in pauli_str.strip().split():
        pauli_char = token[0]
        qubit_idx = int(token[1:])
        ops[qubit_idx] = pauli_char

    # Build the matrix via Kronecker products
    mat = pauli_map[ops[0]]
    for i in range(1, n_qubits):
        mat = np.kron(mat, pauli_map[ops[i]])
    return mat


H_lih = lih_identity_offset * np.eye(2**lih_n_qubits)
for term_str, coeff in lih_terms.items():
    H_lih += coeff * build_pauli_matrix(term_str, lih_n_qubits)

exact_lih = np.linalg.eigvalsh(H_lih)[0]

print(f"LiH VQE energy:  {lih_result['energy']:.8f} Ha")
print(f"LiH exact energy: {exact_lih:.8f} Ha")
print(f"Error:            {lih_result['energy'] - exact_lih:.8f} Ha")

对于 4 量子比特 LiH 问题，具有 2 层（16 个参数）的硬件高效 Ansatz 通常能达到化学精度（ $1.6 \times 10^{- 3}$ Ha）或更好的准确度。

解析

数学基础：变分原理深入理解

变分原理是线性代数应用于量子力学的基本结果。给定具有本征值 $E_{0} \leq E_{1} \leq \dots$ 的厄米算子 $H$ ，任意归一化态 $| ψ ⟩$ 可以在本征基中展开：

| ψ ⟩ = \sum_{k} c_{k} | ψ_{k} ⟩

能量期望值为：

⟨ ψ | H | ψ ⟩ = \sum_{k} | c_{k} |^{2} E_{k} \geq \sum_{k} | c_{k} |^{2} E_{0} = E_{0}

此不等式成立是因为所有本征值 $E_{k} \geq E_{0}$ 且 $\sum_{k} | c_{k} |^{2} = 1$ 。当且仅当 $k > 0$ 时 $c_{k} = 0$ （即 $| ψ ⟩ = | ψ_{0} ⟩$ ）时等式成立。

对 VQE 的影响：

任意试探态给出 $E_{0}$ 的有效上界。
随着试探态逼近 $| ψ_{0} ⟩$ ，界限变紧。
Ansatz 的质量决定了我们能多接近基态。
即使较差的 Ansatz 也给出有用的上界。

变分原理还提供了自然的收敛判据：如果 $E ({\vec{θ}}_{new}) < E ({\vec{θ}}_{old})$ ，新参数保证更接近（在能量上）基态。

Ansatz 设计考量

Ansatz 的选择对 VQE 性能至关重要。关键考量：

硬件高效 Ansatz

本教程中使用的硬件高效 Ansatz 形式为：

U (\vec{θ}) = \prod_{l = 1}^{L} [\prod_{q} R Y_{q} (θ_{q, l}) \cdot CNOT ladder \cdot \prod_{q} R Y_{q} (θ_{q, l}^{'})]

性质：

每层为 $n$ 个量子比特添加 $2 n$ 个参数。
CNOT 阶梯提供纠缠。
$R Y$ 旋转在 Bloch 球的 $Y Z$ 平面上探索。
足够多的层可以逼近任意酉操作（通用性）。

化学启发 Ansatz (UCCSD)

酉耦合簇单双激发(Unitary Coupled Cluster Singles and Doubles, UCCSD) Ansatz 从经典耦合簇理论导出：

| ψ (\vec{θ}) ⟩ = e^{T (\vec{θ}) - T^{†} (\vec{θ})} | ψ_{HF} ⟩

其中 $T = \sum_{i a} θ_{i}^{a} a_{a}^{†} a_{i} + \sum_{i < j, a < b} θ_{i j}^{a b} a_{a}^{†} a_{b}^{†} a_{j} a_{i}$ 是簇算子， $i, j$ 标记占据轨道， $a, b$ 标记虚轨道。Hartree-Fock 态 $| ψ_{HF} ⟩$ 用作参考态。

UCCSD 产生更深的电路但编码了物理关联，使其更可靠地到达真实基态。然而，参数数量按 $O (N_{occ}^{2} N_{virt}^{2})$ 增长，对于较大分子可能代价高昂。

表达能力与可训练性的权衡

表达能力(Expressibility)：Ansatz 能达到的希尔伯特空间比例。更多参数通常增加表达能力。
可训练性(Trainability)：优化器是否能找到好的参数。高表达性 Ansatz 可能遇到贫瘠高原(Barren Plateaus)，梯度指数衰减。

理想的 Ansatz 有足够的表达能力表示基态，但又有足够结构避免贫瘠高原。对于化学问题，问题启发式 Ansatz 如 UCCSD 通过编码已知物理结构实现这一点。

Pauli 项分组以提高测量效率

哈密顿量分解为 Pauli 字符串，每个字符串必须单独测量（或在兼容组中）。关键洞察是相同基中的 Pauli 算子对易，可以同时测量。

对于我们的 H $_{2}$ 哈密顿量项 ${Z_{0}, Z_{1}, Z_{0} Z_{1}, X_{0} X_{1}, Y_{0} Y_{1}}$ ：

项	测量基	所需旋转
$Z_{0}$	计算（ $Z$ ）基	无
$Z_{1}$	计算（ $Z$ ）基	无
$Z_{0} Z_{1}$	计算（ $Z$ ）基	无（与 $Z_{0}$ 、 $Z_{1}$ 对易）
$X_{0} X_{1}$	$X$ 基	两个量子比特施加 $H$
$Y_{0} Y_{1}$	$Y$ 基	两个量子比特施加 $R X (π / 2)$

$Z$ 类项可以全部从同一组测量中估计，将不同的电路评估次数从 5 减少到 3（一次 $Z$ 基、一次 $X$ 基、一次 $Y$ 基）。这种技术称为 Pauli 分组(Pauli Grouping)或逐量子比特对易分组(Qubit-Wise Commutativity Grouping)。

对于更大的分子，存在更复杂的分组策略：

逐量子比特对易(Qubit-Wise Commutativity, QWC)：将每个量子比特独立对易的项分组。将组数从 $O (N^{4})$ 减少到约 $O (N^{2})$ 。
一般对易性：将作为矩阵对易的项分组。需要同时将组对角化的基旋转。可进一步减少组数。
经典阴影层析(Classical Shadow Tomography)：使用随机测量同时估计所有 Pauli 期望值，总测量次数更少。

优化景观与贫瘠高原

贫瘠高原(Barren Plateaus) 是变分量子算法的基本挑战。在贫瘠高原中，代价函数的梯度随量子比特数量指数衰减：

Var [\frac{\partial E}{\partial θ_{i}}] \sim O (\frac{1}{2^{n}})

这意味着对于大系统，能量景观变得指数平坦，基于梯度的优化器无法辨别正确方向。

贫瘠高原的原因

表达能力：均匀覆盖整个希尔伯特空间的高表达性 Ansatz 产生平坦景观。能量期望逼近 $H$ 的迹除以 $2^{n}$ ，方差指数级小。
纠缠：产生体积律纠缠的深层纠缠电路导致贫瘠高原。梯度信号在指数增长的希尔伯特空间中被稀释。
全局代价函数：使用全局代价函数（如高度非局域 $H$ 的 $⟨ H ⟩$ ）加剧了问题。测量少量量子比特上可观测量的局部代价函数更可训练。
噪声：硬件噪声进一步平坦化景观，将能量推向最大混合态。

恒等块初始化(Identity-Block Initialization)策略将每层起始点设在恒等附近（小角度），因此初始电路几乎不做任何操作，梯度景观表现良好。优化器随后根据需要逐渐增加电路表达能力。

对于 2 量子比特和 4 参数 Ansatz 的 H $_{2}$ ，贫瘠高原不是问题。它们在 10+ 量子比特和深层 Ansatz 的系统中变得关键。

误差缓解考量

真实量子硬件引入误差，使估计能量偏移。VQE 的常见误差缓解技术包括：

零噪声外推(Zero-Noise Extrapolation, ZNE)

在多个噪声水平下运行电路（通过拉伸门时间或插入恒等等效门对），然后外推到零噪声：

E (0) = \sum_{i} c_{i} E (λ_{i})

其中 $λ_{i}$ 是噪声缩放因子， $c_{i}$ 是外推系数。

读出误差缓解

校准测量误差矩阵 $M$ （读出的 $2^{n} \times 2^{n}$ 混淆矩阵）并校正观测概率：

{\vec{p}}_{corrected} = M^{- 1} {\vec{p}}_{observed}

对于 $n$ 个量子比特，完整校准矩阵需要 $2^{n}$ 个校准电路，但张量积近似将其减少到 $O (n)$ 个电路。

对称性验证

许多分子哈密顿量具有对称性（粒子数、自旋、宇称）。测量后，丢弃违反这些对称性的比特串。这种"后选择"以减少有效采样次数为代价，移除了一部分噪声结果。

概率误差消除

表征每个门的噪声信道并概率性地施加逆操作以消除误差。这产生无偏估计，但采样开销随电路深度指数增长。

标度分析

VQE 的计算代价随量子比特数（轨道数）和哈密顿量项数标度：

分子	量子比特	Pauli 项	每项采样次数	总采样次数
H $_{2}$ （STO-3G）	2	5	10,000	50,000
LiH（STO-3G）	4	~100	10,000	1,000,000
BeH $_{2}$ （STO-3G）	6	~200	10,000	2,000,000
H $_{2}$ O（STO-3G）	8	~600	10,000	6,000,000
NH $_{3}$ （STO-3G）	10	~2,000	10,000	20,000,000

Pauli 项数按 $N$ 个自旋轨道的 $O (N^{4})$ 增长，尽管对称性和约化可以将其减少到 $O (N^{3})$ 或更优。

与经典方法的比较

方法	标度	精度	适用性
FCI（精确）	$O (e^{N})$	精确	仅小系统（ $< 30$ 轨道）
CCSD(T)	$O (N^{7})$	弱关联下约 1 mHa	金标准，强关联下失败
DFT	$O (N^{3})$	约 10 mHa（依赖于泛函）	广泛使用，但精度不可预测
VQE (UCCSD)	$O (N^{4})$ 项	约 1 mHa（依赖 Ansatz）	强关联在 QPU 上有前景
VQE（硬件高效）	$O (N^{4})$ 项	可变（依赖深度）	近期可行，可预测性较低

VQE 并非要取代小规模、弱关联分子的经典方法。其价值在于解决经典方法失效的强关联(Strongly Correlated)系统：过渡金属络合物、键断裂过程和具有多参考特征的激发态。

经典优化器比较

不同优化器在 VQE 语境中有不同权衡：

优化器	类型	需要梯度	噪声容忍度	典型迭代次数
COBYLA	无导数	否	高	200~1000
Nelder-Mead	无导数	否	中	200~500
SLSQP	基于梯度	估计	低	50~200
L-BFGS-B	基于梯度	估计	低	30~100
SPSA	随机	近似	高	500~5000
Adam（量子）	基于梯度	参数移位	中	100~500

对于基于采样的 VQE，推荐无梯度优化器（COBYLA、Nelder-Mead），因为采样噪声使梯度估计不可靠。对于精确（状态向量）VQE，基于梯度的方法收敛更快。

SPSA（同时扰动随机逼近）特别适合硬件实验，因为无论参数数量多少，每次迭代仅需 2 次电路评估，而有限差分梯度方法需要 $O (n)$ 次。

常见陷阱与调试技巧

陷阱 1：错误的哈密顿量系数

系数对分子几何构型、基组和活性空间选择非常敏感。运行 VQE 前务必对照已知经典计算（FCI 或 CCSD）进行验证。

调试：通过对角化哈密顿量矩阵计算精确基态能量，并与文献值比较。对于 R=0.74 A 的 H $_{2}$ ， $E_{0} \approx - 1.1373$ Ha。

陷阱 2：Ansatz 表达能力不足

如果 Ansatz 无法表示基态，VQE 将收敛到严格高于 $E_{0}$ 的值。差距 $E_{VQE} - E_{0}$ 量化了 Ansatz 误差。

调试：增加层数并观察能量是否继续下降。如果在远高于 $E_{0}$ 处平台化，可能需要更改 Ansatz 结构。

陷阱 3：优化景观中的局部极小值

VQE 能量景观非凸，可能有多个局部极小值。

调试：使用多次随机重启（10~50 次）并比较最佳能量。如果不同重启收敛到显著不同的能量，景观有多个极小值，需要更多重启。

陷阱 4：错误的测量基旋转

忘记在测量前施加基旋转会导致错误的 Pauli 期望值。 $X$ 算子需要 Hadamard， $Y$ 算子需要 $R X (π / 2)$ 。

调试：在运行完整优化前，逐个 Pauli 项对照状态向量结果进行验证。

陷阱 5：测量采样次数不足

采样次数太少导致噪声能量估计，妨碍收敛。

调试：先用大量采样（10,000+）进行验证，然后减少到找到仍可接受精度的最小值。统计误差按 $1 / \sqrt{S}$ 标度，其中 $S$ 是采样次数。

陷阱 6：恒等项处理

恒等项 $c_{0} I$ 贡献常数偏移，通常从 PauliOperator 中排除以减少电路评估。忘记加回会给出不正确的绝对能量（尽管优化景观形状不变）。

调试：报告绝对能量时，始终计算并显示包含恒等偏移的总能量。

总结

组件	在 pyqpanda3 中的实现
哈密顿量	`hamiltonian.PauliOperator({"Z0": -0.398, "X0 X1": 0.398, ...})`
Ansatz 电路	`core.RY`、`core.CNOT`、`core.QProg`
期望值（精确）	`core.expval_pauli_operator(prog, hamiltonian)`
基于采样的测量	`core.CPUQVM().run(prog, shots)` + `result.get_counts()`
经典优化	`scipy.optimize.minimize` 配合 COBYLA/SLSQP
精确对角化	`numpy.linalg.eigvalsh(H_matrix)`

后续步骤

QAOA教程 -- 用于组合优化的变分算法
噪声表征 -- 为硬件上的 VQE 构建真实噪声模型
自定义量子门 -- 实现化学启发 Ansatz 的专用门

VQE -- 变分量子本征求解器(Variational Quantum Eigensolver) ​

问题 ​

寻找基态能量 ​

精确对角化的经典困难 ​

变分原理(Variational Principle) ​

为什么选择 VQE？ ​

方案 ​

逐步 VQE 工作流程 ​

步骤 1：构建分子哈密顿量 ​

步骤 2：通过 Jordan-Wigner 变换映射到量子比特 ​

步骤 3：选择 Ansatz ​

步骤 4：测量 Pauli 项期望值 ​

步骤 5：优化参数 ​

代码 ​

H2 分子 VQE 与硬件高效 Ansatz ​

定义 H2 哈密顿量 ​

构建硬件高效 Ansatz ​

使用状态向量进行精确期望值计算 ​

通过基于采样的测量计算期望值 ​

从哈密顿量计算总能量 ​

使用 scipy.optimize (COBYLA) 优化 ​

跟踪并绘制收敛历史 ​

与精确对角化比较 ​

通用 VQE 求解器函数 ​

不同键长下的 VQE -- 势能面(Potential Energy Surface) ​

绘制势能面 ​

扩展到 LiH（4 量子比特示例） ​

解析 ​

数学基础：变分原理深入理解 ​

Ansatz 设计考量 ​

硬件高效 Ansatz ​

化学启发 Ansatz (UCCSD) ​

表达能力与可训练性的权衡 ​

Pauli 项分组以提高测量效率 ​

优化景观与贫瘠高原 ​

贫瘠高原的原因 ​

误差缓解考量 ​

零噪声外推(Zero-Noise Extrapolation, ZNE) ​

读出误差缓解 ​

对称性验证 ​

概率误差消除 ​

标度分析 ​

与经典方法的比较 ​

经典优化器比较 ​

常见陷阱与调试技巧 ​

陷阱 1：错误的哈密顿量系数 ​

陷阱 2：Ansatz 表达能力不足 ​

陷阱 3：优化景观中的局部极小值 ​

陷阱 4：错误的测量基旋转 ​

陷阱 5：测量采样次数不足 ​

陷阱 6：恒等项处理 ​

总结 ​

后续步骤 ​

VQE -- 变分量子本征求解器(Variational Quantum Eigensolver)

问题

寻找基态能量

精确对角化的经典困难

变分原理(Variational Principle)

为什么选择 VQE？

方案

逐步 VQE 工作流程

步骤 1：构建分子哈密顿量

步骤 2：通过 Jordan-Wigner 变换映射到量子比特

步骤 3：选择 Ansatz

步骤 4：测量 Pauli 项期望值

步骤 5：优化参数

代码

H2 分子 VQE 与硬件高效 Ansatz

定义 H2 哈密顿量

构建硬件高效 Ansatz

使用状态向量进行精确期望值计算

通过基于采样的测量计算期望值

从哈密顿量计算总能量

使用 scipy.optimize (COBYLA) 优化

跟踪并绘制收敛历史

与精确对角化比较

通用 VQE 求解器函数

不同键长下的 VQE -- 势能面(Potential Energy Surface)

绘制势能面

扩展到 LiH（4 量子比特示例）

解析

数学基础：变分原理深入理解

Ansatz 设计考量

硬件高效 Ansatz

化学启发 Ansatz (UCCSD)

表达能力与可训练性的权衡

Pauli 项分组以提高测量效率

优化景观与贫瘠高原

贫瘠高原的原因

误差缓解考量

零噪声外推(Zero-Noise Extrapolation, ZNE)

读出误差缓解

对称性验证

概率误差消除

标度分析

与经典方法的比较

经典优化器比较

常见陷阱与调试技巧

陷阱 1：错误的哈密顿量系数

陷阱 2：Ansatz 表达能力不足

陷阱 3：优化景观中的局部极小值

陷阱 4：错误的测量基旋转

陷阱 5：测量采样次数不足

陷阱 6：恒等项处理

总结

后续步骤