Variational Quantum Circuits

Variational quantum circuits (VQCs) are parameterized quantum circuits whose gate parameters can be optimized classically. They form the backbone of variational quantum algorithms like VQE and QAOA.

What Are Variational Circuits?

A variational quantum circuit consists of:

Parameterized gates: Quantum gates with tunable parameters (e.g., rotation angles)
Cost function: An observable (typically a Hamiltonian) to minimize
Classical optimizer: Updates gate parameters based on measured cost

The optimization loop:

The VQCircuit Module

The vqcircuit module provides tools for building and optimizing parameterized circuits:

python

from pyqpanda3 import vqcircuit

Key Components

Component	Description
`VQCircuit`	Parameterized quantum circuit builder
`Parameter`	Parameter container for variational parameters
`DiffMethod`	Differentiation method (ADJOINT_DIFF)
`ResGradients`	Gradient values for single parameter set
`ResNGradients`	Gradient values for N parameter sets
`ResGradientsAndExpectation`	Combined gradients and expectation value

VQCircuit Class API

The VQCircuit class is the central interface for constructing, evaluating, and differentiating parameterized quantum circuits. This section documents every public method with usage examples.

Constructor

python

from pyqpanda3 import vqcircuit, core

# Create a VQCircuit (pre_split defaults to True)
vqc = vqcircuit.VQCircuit()

The constructor accepts an optional bool pre_split parameter (defaults to True), which controls whether the circuit is pre-split for gradient computation.

Adding Gates

Variational gates are appended to the circuit using the standard << operator, just like building an ordinary QProg. The difference is that parameter indices (of type MutableParamType) are supplied instead of numeric angles.

python

from pyqpanda3 import core, vqcircuit

# Build a single-layer variational circuit on 2 qubits
# Param() is a method on VQCircuit that returns parameter references
vqc = vqcircuit.VQCircuit()
param_theta_0 = vqc.Param(0)   # parameter index 0
param_theta_1 = vqc.Param(1)   # parameter index 1
param_phi_0   = vqc.Param(2)   # parameter index 2
param_phi_1   = vqc.Param(3)   # parameter index 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)

Setting Parameter Values

Before computing expectations or gradients, you assign concrete numeric values to the parameter slots:

python

from pyqpanda3 import vqcircuit
import numpy as np

# Assign values to all parameters at once via __call__
theta = np.array([0.5, 1.2, 0.8, 2.1])
result = vqc(theta)  # returns VQCircuitResult

Computing Expectation Values

The expectation value of a Hamiltonian with respect to the parameterized state is obtained through the simulation utilities. Using a PauliOperator as the observable:

python

from pyqpanda3 import core
from pyqpanda3.hamiltonian import PauliOperator

# Define a simple Z-Z Hamiltonian on 2 qubits
hamiltonian = PauliOperator({"Z0 Z1": 1.0, "Z0": 0.5, "Z1": 0.3})

# Build the quantum program
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)

# Compute the expectation value (statevector backend)
exp_val = core.expval_pauli_operator(prog, hamiltonian)
print(f"Expectation value: {exp_val:.6f}")

Computing Gradients

Gradients are computed with respect to every parameter using the adjoint differentiation method:

python

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

# Compute gradients using adjoint differentiation
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

# Assign concrete parameter values
theta = np.array([0.5, 1.2, 0.8, 2.1])

# The gradient API is an instance method on VQCircuit
grad_result = vqc.get_gradients(theta, hamiltonian, diff_method)

# Access individual gradient values
for i in range(len(grad_result)):
    g = grad_result.at(i)
    print(f"  dE/dtheta[{i}] = {g:.6f}")

Computing Gradients and Expectation Together

For efficiency, gradients and the expectation value can be computed in a single pass:

python

# Obtain both gradients and expectation in one evaluation
combined = vqc.get_gradients_and_expectation(theta, hamiltonian, diff_method)

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

Batch Evaluation with ResNGradients

When evaluating gradients for multiple parameter sets at once, the ResNGradients container is returned:

python

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

# Evaluate gradients for a batch of parameter configurations
# Flatten all parameter sets into a single 1D array
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,
])

# Use the batch overload with param_group_total
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 Method Summary

Method	Signature	Return Type	Description
Constructor	`VQCircuit(pre_split=True)`	`VQCircuit`	Create circuit (pre_split controls gradient pre-computation)
Add gate	`vqc << gate`	`VQCircuit`	Append a variational gate to the circuit
Param ref	`vqc.Param(idx)`	`MutableParamType`	Get parameter reference for gate construction
Gradients	`vqc.get_gradients(theta, H, method)`	`ResGradients`	Compute gradients of $⟨ H ⟩$
Batch gradients	`vqc.get_gradients(theta, H, n, method)`	`ResNGradients`	Gradients for `n` parameter sets
Combined	`vqc.get_gradients_and_expectation(theta, H, method)`	`ResGradientsAndExpectation`	Gradients and expectation together

Variational Gates (VQGate)

Variational gates are parameterized gate wrappers. They accept MutableParamType objects instead of fixed angles, allowing parameters to be updated during optimization.

VQGate functions are exported from the core module when called with MutableParamType parameters:

Single-Qubit Parameterized Gates

Gate	Signature	Description
`RX(qubit, param)`	`core.RX(idx, MutableParamType)`	RX rotation with parameter
`RY(qubit, param)`	`core.RY(idx, MutableParamType)`	RY rotation with parameter
`RZ(qubit, param)`	`core.RZ(idx, MutableParamType)`	RZ rotation with parameter
`P(qubit, param)`	`core.P(idx, MutableParamType)`	Phase gate with parameter
`U1(qubit, param)`	`core.U1(idx, MutableParamType)`	U1 gate with parameter
`U2(qubit, p1, p2)`	`core.U2(idx, param1, param2)`	U2 with 2 parameters
`U3(qubit, p1, p2, p3)`	`core.U3(idx, p1, p2, p3)`	U3 with 3 parameters
`U4(qubit, p1, p2, p3, p4)`	`core.U4(idx, p1, p2, p3, p4)`	U4 with 4 parameters
`RPHI(qubit, p1, p2)`	`core.RPhi(idx, param1, param2)`	RPhi with 2 parameters

Two-Qubit Parameterized Gates

Gate	Signature	Description
`CU(ctrl, targ, p1-p4)`	`core.CU(ctrl, targ, ...)`	Controlled-U with 4 params
`CP(ctrl, targ, param)`	`core.CP(ctrl, targ, param)`	Controlled-phase
`CRX(ctrl, targ, param)`	`core.CRX(ctrl, targ, param)`	Controlled-RX
`CRY(ctrl, targ, param)`	`core.CRY(ctrl, targ, param)`	Controlled-RY
`CRZ(ctrl, targ, param)`	`core.CRZ(ctrl, targ, param)`	Controlled-RZ
`RXX(ctrl, targ, param)`	`core.RXX(ctrl, targ, param)`	RXX rotation
`RYY(ctrl, targ, param)`	`core.RYY(ctrl, targ, param)`	RYY rotation
`RZZ(ctrl, targ, param)`	`core.RZZ(ctrl, targ, param)`	RZZ rotation
`RZX(ctrl, targ, param)`	`core.RZX(ctrl, targ, param)`	RZX rotation

Mixed Fixed and Variational Parameters

Multi-parameter gates support mixing fixed values with variational parameters:

python

from pyqpanda3 import core, vqcircuit

# Create a VQCircuit and get a variational parameter reference
vqc = vqcircuit.VQCircuit()
param_index = vqc.Param(0)

# U3 with first two params fixed, third variational
gate = core.U3(0, 1.57, 0.0, param_index)

# U4 with only last param variational
gate = core.U4(0, 0.0, 0.0, 0.0, param_index)

# U2 with first param variational, second fixed
gate = core.U2(0, param_index, 3.14)

Parameter Management

Parameter Class

The Parameter class manages the variational parameters:

python

from pyqpanda3 import vqcircuit

# The Parameter class is typically used through VQCircuit
# VQCircuit manages parameters internally, mutable_parameter_total() returns the total count
vqc = vqcircuit.VQCircuit()

# Query the number of mutable parameters in the circuit (initially 0)
n = vqc.mutable_parameter_total()
print(f"Number of mutable parameters: {n}")

MutableParamType

The MutableParamType specifies where a parameter value comes from:

Index-based: param = vqcircuit.Param([idx_dim0, idx_dim1, ...]) — references a multi-dimensional parameter
Scalar index: param = vqcircuit.Param(idx) — references a single parameter by index

python

from pyqpanda3 import core, vqcircuit

# Create a VQCircuit to hold parameters
vqc = vqcircuit.VQCircuit()

# Create parameterized RX gate using VQCircuit.Param()
# The MutableParamType specifies which parameter slot to use
param_0 = vqc.Param(0)  # References parameter at index 0
param_1 = vqc.Param(1)  # References parameter at index 1

# Append variational gates to the VQCircuit
vqc << core.RX(0, param_0)
vqc << core.RY(1, param_1)

Gradient Computation

Adjoint Differentiation

The adjoint differentiation method efficiently computes gradients by evolving the circuit forward and then backward (adjoint). This avoids the $2 n$ circuit evaluations required by the parameter shift rule.

Reference: Jones and Gacon, "Efficient calculation of gradients in classical simulations of variational quantum algorithms", arXiv:2009.02823.

python

from pyqpanda3 import vqcircuit

# Specify the differentiation method
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

ResGradients - Single Parameter Set

ResGradients holds gradient values for a single set of parameters:

python

# Assume gradients is a ResGradients object from VQCircuit computation
gradients = vqc_result.gradients()  # Get ResGradients

# Get number of gradients
n_grads = len(gradients)

# Get all gradient values as a list
all_grads = gradients.gradients()

# Get gradient for parameter at index i
grad_i = gradients.at(i)

# Get gradient for multi-dimensional parameter
grad_md = gradients.at([idx_dim0, idx_dim1])

# String representation
print(gradients)

ResNGradients - N Parameter Sets

ResNGradients holds gradient values for N sets of parameters (batch evaluation):

python

# Assume n_gradients is a ResNGradients object from vqc.get_gradients(..., batch_size, ...)
n_gradients = vqc.get_gradients(theta_flat, hamiltonian, batch_size, diff_method)

# Get number of parameter sets (N)
n_sets = len(n_gradients)

# Get ResGradients for the i-th parameter set
grad_set_i = n_gradients.at(i)

# Get all gradient data as nested list
all_data = n_gradients.data()

# String representation
print(n_gradients)

ResGradientsAndExpectation - Combined Result

ResGradientsAndExpectation holds both gradient values and the expectation value:

python

# Assume result is a ResGradientsAndExpectation object from vqc.get_gradients_and_expectation(...)
result = vqc.get_gradients_and_expectation(theta, hamiltonian, diff_method)

# Get the expectation value
exp_val = result.expectation_val()

# Get all gradient values as a list
grads = result.gradients()

# Get gradient for specific parameter
grad_i = result.at([idx])

# Get combined data: (expectation, [gradients])
data = result.data()

# String representation
print(result)

Mathematical Foundation

Parameterized Quantum State

A parameterized circuit $U (θ)$ applied to $| 0 ⟩^{\otimes n}$ produces:

| ψ (θ) ⟩ = U (θ) | 0 ⟩^{\otimes n}

Expectation Value

The cost function is the expectation value of a Hamiltonian $H$ :

⟨ H ⟩_{θ} = ⟨ ψ (θ) | H | ψ (θ) ⟩

Gradient via Adjoint Differentiation

The gradient with respect to parameter $θ_{k}$ is:

\frac{\partial ⟨ H ⟩}{\partial θ_{k}} = 2 Re [⟨ ψ_{L} | ψ_{R} ⟩]

where:

| ψ_{L} ⟩ = \frac{\partial U (θ)}{\partial θ_{k}} | 0 ⟩, | ψ_{R} ⟩ = U^{†} (θ) H | ψ (θ) ⟩

The adjoint method avoids storing intermediate states by computing the inner product during a single forward-backward pass.

Parameter Shift Rule (Alternative)

For comparison, the parameter shift rule computes each gradient using two circuit evaluations:

\frac{\partial ⟨ H ⟩}{\partial θ_{k}} = \frac{1}{2} [⟨ H ⟩_{θ_{k} + π / 2} - ⟨ H ⟩_{θ_{k} - π / 2}]

This requires $2 n$ circuit evaluations for $n$ parameters, while adjoint differentiation requires only one forward-backward pass.

VQA Workflow

Building a VQCircuit Step by Step

This section walks through the complete construction of a variational quantum circuit from scratch, explaining each decision along the way. We will build a 3-qubit circuit suitable for a small optimization problem.

Step 1: Decide on the Ansatz Structure

Before writing any code, decide what ansatz (circuit structure) is appropriate for your problem. For this walkthrough we use a hardware-efficient ansatz with alternating rotation and entangling layers.

Step 2: Import Modules and Define Constants

python

from pyqpanda3 import core, vqcircuit
import numpy as np

# Circuit configuration
N_QUBITS = 3       # Number of qubits
DEPTH = 2          # Number of repeated layers
# Each layer has 2 rotations per qubit (RY then RZ)
N_PARAMS = N_QUBITS * 2 * DEPTH  # = 12 parameters total

Step 3: Create the Parameter Object

The Parameter object tracks all variational parameters and their current values.

python

# Initialize the parameter container
params = vqcircuit.Parameter()
print(f"Initial parameter count: {params.size()}")

Step 4: Build the Ansatz Circuit

Construct the circuit layer by layer. Each layer consists of RY and RZ rotations on every qubit, followed by CNOT entangling gates.

python

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  # running index for parameter allocation

    for layer in range(depth):
        # --- Rotation sub-layer ---
        for q in range(n_qubits):
            # RY rotation for qubit q (use 0.0 as placeholder)
            prog << core.RY(q, 0.0)
            param_idx += 1

            # RZ rotation for qubit q (use 0.0 as placeholder)
            prog << core.RZ(q, 0.0)
            param_idx += 1

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

        # Close the ring for full connectivity
        if n_qubits > 2:
            prog << core.CNOT(n_qubits - 1, 0)

    total_params = param_idx
    return prog, total_params

# Build the circuit
ansatz_prog, n_params_used = build_hardware_efficient_ansatz(N_QUBITS, DEPTH)
print(f"Circuit built with {n_params_used} parameters")

Step 5: Define the Observable (Hamiltonian)

Choose a Hamiltonian to serve as the cost function. For demonstration we use a transverse-field Ising model:

H = - \sum_{⟨ i, j ⟩} Z_{i} Z_{j} - h \sum_{i} X_{i}

python

from pyqpanda3.hamiltonian import PauliOperator

# Build the Hamiltonian as a PauliOperator
hamiltonian_terms = {}

# ZZ coupling terms (nearest neighbor on a ring)
for i in range(N_QUBITS):
    j = (i + 1) % N_QUBITS
    hamiltonian_terms[f"Z{i} Z{j}"] = -1.0

# Transverse field terms
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}")

Step 6: Initialize Parameters

Set initial parameter values. A common choice is random angles uniformly distributed in $[0, 2 π)$ :

python

# Random initialization
np.random.seed(42)
initial_theta = np.random.uniform(0, 2 * np.pi, n_params_used)
print(f"Initial parameters: {initial_theta}")

Step 7: Evaluate the Cost Function

python

def cost_function(theta_array):
    """Compute the expectation value of the Hamiltonian for given parameters."""
    # Build a new circuit with the given parameter values
    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)

    # Evaluate expectation of the full Hamiltonian
    return core.expval_pauli_operator(prog, hamiltonian)

# Test the cost function
energy_0 = cost_function(initial_theta)
print(f"Initial energy: {energy_0:.6f}")

Step 8: Run the Optimization Loop

python

from scipy.optimize import minimize

# Optimize with COBYLA (derivative-free)
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}")

Step 9: Verify the Result

After optimization, inspect the final quantum state and check convergence:

python

# Rebuild circuit at optimal parameters
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)

# Get the final state vector
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}")

# Compute measurement probabilities
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}")

Batch Parameter Evaluation

When exploring the parameter landscape or performing population-based optimization, it is efficient to evaluate gradients for multiple parameter sets simultaneously. The ResNGradients container is designed for this purpose.

Why Batch Evaluation?

Evaluating gradients for $N$ parameter sets individually requires $N$ separate simulation runs, each with its own circuit construction overhead. Batch evaluation amortizes setup costs and is particularly useful for:

Grid search over a discrete parameter mesh
Evolutionary algorithms that maintain a population of candidate solutions
Parallel gradient comparison to identify the most promising descent direction

Batch Evaluation Workflow

Example: Evaluating a Parameter Grid

python

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

# --- Build a simple 2-qubit variational circuit ---
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

# Define the Hamiltonian (Z-Z interaction)
H = PauliOperator({"Z0 Z1": 1.0, "Z0": 0.5})

# --- Create a grid of parameter values ---
grid_size = 5
theta_0_vals = np.linspace(0, np.pi, grid_size)
theta_1_vals = np.linspace(0, np.pi, grid_size)

# Collect all parameter combinations
batch = []
for t0 in theta_0_vals:
    for t1 in theta_1_vals:
        batch.append([t0, t1])

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

Computing Batch Gradients

python

# Compute gradients for each parameter set
# Using the adjoint differentiation method
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

gradient_results = []
energies = []

for theta in batch:
    prog = build_2q_circuit(theta)

    # Compute expectation value
    energy = core.expval_pauli_operator(prog, H)
    energies.append(energy)

    # For demonstration, compute parameter-shift gradients manually
    shift = np.pi / 2
    grads = []
    for i in range(2):
        # Forward shift
        theta_plus = theta.copy()
        theta_plus[i] += shift
        prog_plus = build_2q_circuit(theta_plus)
        e_plus = core.expval_pauli_operator(prog_plus, H)

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

        # Gradient via parameter shift rule
        grad_i = 0.5 * (e_plus - e_minus)
        grads.append(grad_i)

    gradient_results.append(grads)

# Find the parameter set with the steepest descent
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}")

Analyzing the Energy Landscape

python

# Reshape results into a 2D grid for visualization analysis
energy_grid = np.array(energies).reshape(grid_size, grid_size)

# Report the energy landscape statistics
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}")

# Identify local minima
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}")

Gradient Verification and Debugging

When a variational algorithm fails to converge, the gradient computation is often the culprit. This section describes strategies to verify gradient correctness and debug common issues.

Numerical Gradient Check

The most reliable verification method is to compare analytic gradients against a finite-difference approximation:

\frac{\partial f}{\partial θ_{k}} \approx \frac{f (θ_{k} + ϵ) - f (θ_{k} - ϵ)}{2 ϵ}

python

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):
        # Perturb parameter k forward
        theta_plus = theta.copy()
        theta_plus[k] += epsilon
        f_plus = cost_fn(theta_plus)

        # Perturb parameter k backward
        theta_minus = theta.copy()
        theta_minus[k] -= epsilon
        f_minus = cost_fn(theta_minus)

        # Central difference
        grads[k] = (f_plus - f_minus) / (2 * epsilon)

    return grads


# Example usage: verify adjoint gradients against numerical ones
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)

# Test at a specific parameter point
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))

Gradient Consistency Check

For a robust check, compare gradients at multiple parameter points:

python

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.
    """
    # Compute at two different epsilon values
    g1 = numerical_gradient(cost_fn, theta, epsilon=epsilon)
    g2 = numerical_gradient(cost_fn, theta, epsilon=epsilon * 0.1)

    # They should agree to within the tolerance
    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

# Verify at several random parameter points
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)

Common Debugging Checklist

When gradients appear incorrect, check the following:

Detecting Barren Plateaus

A barren plateau occurs when gradient magnitudes vanish exponentially with system size, making optimization impossible. Detect it by monitoring gradient statistics:

python

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

    # Compute variance of each gradient component
    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}")

    # Heuristic: if variance is below 1e-8, barren plateau is likely
    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

# Run the check
check_barren_plateau(demo_cost, n_params=3)

Logging Gradient Norms During Optimization

Monitoring gradient norms during optimization helps identify convergence issues:

python

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}")

        # Check if gradient is decreasing
        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!")

# Usage during optimization
monitor = GradientMonitor()

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

# Run with monitoring
result = minimize(monitored_cost, np.array([0.5, 1.0, 1.5]),
                  method='COBYLA', options={'maxiter': 50})
monitor.report()

Performance Optimization

This section provides practical tips for improving the performance of variational quantum circuit computations in pyqpanda3.

Choose the Right Differentiation Method

Method	Cost per gradient	Precision	Backend requirement
Adjoint	$O (n_{gates})$	Machine precision	Statevector
Parameter shift	$2 n_{params}$ circuits	Machine precision	Any backend
Finite difference	$2 n_{params}$ circuits	Depends on $ϵ$	Any backend

Recommendation: Always use adjoint differentiation when working with a statevector simulator. It computes all $n$ gradients in a single forward-backward pass.

python

# Preferred: adjoint differentiation (single pass)
diff_method = vqcircuit.DiffMethod.ADJOINT_DIFF

Reduce Circuit Depth

Shallower circuits lead to faster simulation and smaller numerical errors:

python

# Instead of deep hardware-efficient ansatz (many layers)
# Use problem-inspired ansatz with minimal depth

# BAD: Generic deep circuit
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)  # placeholder
        for q in range(n_qubits - 1):
            prog << core.CNOT(q, q + 1)
    return prog

# GOOD: Problem-tailored shallow circuit (e.g., for H2)
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

Exploit Hamiltonian Symmetry

Group commuting Pauli terms to reduce the number of circuit evaluations:

python

from pyqpanda3.hamiltonian import PauliOperator

# When measuring a Hamiltonian with many Pauli terms,
# group terms that can be measured simultaneously (same Pauli basis)

# For example, Z0 Z1 and Z0 can both be measured in the Z basis
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

        # Build a combined PauliOperator for this basis group
        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

Optimizer Selection Guide

python

# Gradient-free (robust, easy to use)
from scipy.optimize import minimize
result_cobyla = minimize(cost_function, theta0, method='COBYLA',
                         options={'maxiter': 200})

# Gradient-based (faster convergence when gradients are accurate)
result_lbfgsb = minimize(cost_function, theta0, method='L-BFGS-B',
                         jac=gradient_function,
                         options={'maxiter': 200})

Parameter Initialization Strategies

Good initialization can drastically reduce the number of optimization iterations:

python

import numpy as np

# Strategy 1: Random uniform
theta_rand = np.random.uniform(0, 2 * np.pi, n_params)

# Strategy 2: Small random perturbation around zero
theta_small = np.random.normal(0, 0.1, n_params)

# Strategy 3: Identity-block initialization (avoids barren plateaus)
# Initialize so that U(theta_0) is close to identity
theta_identity = np.zeros(n_params)

# Strategy 4: Structured initialization based on problem symmetry
# For QAOA, initialize gamma near 0.5 and beta near 0.25
theta_qaoa = np.concatenate([
    np.full(p, 0.5),    # gamma parameters
    np.full(p, 0.25),   # beta parameters
])

Caching Intermediate Results

When the cost function is called repeatedly with similar parameters, caching can save computation time:

python

from functools import lru_cache
from pyqpanda3.hamiltonian import PauliOperator

# Cache the expectation value for Pauli terms (when circuit is deterministic)
# Note: only works when parameters are hashable (tuples, not arrays)

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

# Convert theta to a hashable tuple before calling
theta_tuple = tuple(theta_test)
exp_zz = cached_expectation(theta_tuple, "Z0 Z1")

Complete Example: VQE for H₂

python

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

# Step 1: Define the H₂ Hamiltonian (simplified)
# H = c0*I + c1*Z0 + c2*Z1 + c3*Z0*Z1 + c4*X0*X1 + c5*Y0*Y1
# Coefficients for bond length 0.735 Å
h2_coeffs = {
    "": -1.0523792,
    "Z0": 0.39793742,
    "Z1": -0.39793742,
    "Z0 Z1": -0.01128010,
    "X0 X1": 0.18093119,
}

# Step 2: Define ansatz circuit (hardware-efficient)
def build_ansatz(params_array):
    """Build a parameterized ansatz for H₂."""
    prog = core.QProg()

    # Parameterized single-qubit rotations
    prog << core.RY(0, params_array[0])
    prog << core.RY(1, params_array[1])

    # Entangling layer
    prog << core.CNOT(0, 1)

    # Second rotation layer
    prog << core.RY(0, params_array[2])
    prog << core.RY(1, params_array[3])

    # Entangling layer
    prog << core.CNOT(0, 1)

    return prog

# Step 3: Define cost function
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

        # Create Pauli operator for this term
        op = PauliOperator({pauli_str: 1.0})
        exp_val = core.expval_pauli_operator(prog, op, shots=1000)
        total += coeff * exp_val

    return total

# Step 4: Optimize using scipy
from scipy.optimize import minimize

# Initial parameters
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")

Interpreting VQE Results

The output of the VQE optimization provides several pieces of information that require careful interpretation.

Energy Accuracy

The exact ground state energy of $H_{2}$ at bond length 0.735 A is approximately $- 1.137275$ Hartree. The VQE result should be compared against this value:

Δ E = E_{VQE} - E_{exact}

python

# Compare VQE result against the exact value
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 Parameter Analysis

The optimized parameters reveal the structure of the ground state circuit:

python

# Inspect optimal parameters
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)")

# Check if any parameters are near 0 or pi (indicates gate is effectively identity or Pauli)
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)")

State Composition

Reconstruct the optimal quantum state and examine which computational basis states contribute:

python

# Build circuit at optimal parameters
optimal_prog = build_ansatz(optimal_theta)

# Obtain the state vector
machine = core.CPUQVM()
machine.run(optimal_prog)
state = machine.result().get_state_vector()

# Print probability distribution over basis states
print("\nBasis state probabilities:")
for i, amp in enumerate(state):
    prob = abs(amp) ** 2
    if prob > 0.001:
        label = format(i, '02b')  # 2-bit binary string
        # Map qubit ordering (bit 0 = rightmost)
        print(f"  |{label}> : amplitude = {amp:+.6f}, probability = {prob:.6f}")

# For H2, the ground state is approximately:
# |psi> ~ 0.9939|00> - 0.1102|11>
# This corresponds to the Hartree-Fock state with small double-excitation

Convergence Behavior

Analyze the optimizer's convergence trajectory:

python

# Track energy at each function evaluation
energy_history = []

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

# Re-run optimization with tracking
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})

# Report convergence statistics
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}")

# Check if energy stabilized (last 10 evaluations)
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.")

Scaling to Larger Molecules

The $H_{2}$ example uses 2 qubits and 4 parameters. For larger molecules the resource requirements grow:

Molecule	Qubits	Parameters (UCCSD)	Pauli Terms
$H_{2}$ (STO-3G)	2	4	5
LiH (STO-3G)	4	~20	~100
${BeH}_{2}$ (STO-3G)	6	~60	~200
$H_{2} O$ (STO-3G)	8	~150	~600

Variational Circuit Construction Patterns

Hardware-Efficient Ansatz

python

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

    for d in range(depth):
        # Single-qubit rotation layer
        for q in range(n_qubits):
            prog << core.RY(q, params[idx])
            idx += 1
            prog << core.RZ(q, params[idx])
            idx += 1

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

    return prog

QAOA Ansatz

python

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 Reference

ADJOINT_DIFF

The adjoint differentiation method:

Reference: Jones and Gacon (2020)
Complexity: $O (n_{gates})$ per gradient computation (vs $O (n_{params})$ circuit evaluations for parameter shift)
Precision: Machine precision (no statistical noise)
Requirement: Statevector simulation backend

Gradient Result API Summary

ResGradients

Method	Return	Description
`.at(idx)`	`float`	Gradient for parameter at index `idx`
`.at([dims])`	`float`	Gradient for multi-dim parameter
`.gradients()`	`list[float]`	All gradient values
`.__len__()`	`int`	Number of gradients
`.__str__()`	`str`	String representation

ResNGradients

Method	Return	Description
`.at(idx)`	`ResGradients`	Gradients for idx-th parameter set
`.data()`	`list[list[float]]`	All gradient data
`.__len__()`	`int`	Number of parameter sets (N)
`.__str__()`	`str`	String representation

ResGradientsAndExpectation

Method	Return	Description
`.expectation_val()`	`float`	Expectation value
`.gradients()`	`list[float]`	All gradient values
`.at([dims])`	`float`	Gradient for specific parameter
`.data()`	`tuple[float, list[float]]`	(expectation, [gradients])
`.__str__()`	`str`	String representation

QAOA Complete Example: MaxCut

The Quantum Approximate Optimization Algorithm (QAOA) is a variational algorithm designed for combinatorial optimization problems. This section demonstrates QAOA applied to the MaxCut problem on a small graph.

Problem Definition

Given an undirected graph $G = (V, E)$ , MaxCut asks for a partition of the vertices $V = V_{0} \cup V_{1}$ that maximizes the number of edges crossing the cut:

MaxCut (G) = max_{x \in {0, 1}^{n}} \sum_{(i, j) \in E} 1 [x_{i} \neq x_{j}]

This can be encoded as minimizing a cost Hamiltonian:

H_{C} = - \sum_{(i, j) \in E} \frac{1}{2} (1 - Z_{i} Z_{j})

The negative sign ensures that minimizing $H_{C}$ corresponds to maximizing the cut size.

Example Graph

We use a 4-vertex ring graph with edges $(0, 1), (1, 2), (2, 3), (3, 0)$ :

The optimal MaxCut value for this graph is 4 (all edges can be cut by partitioning {0,2} vs {1,3}).

QAOA Circuit Structure

The QAOA circuit with $p$ layers alternates between the cost unitary and the mixer unitary:

| γ, β ⟩ = \prod_{l = 1}^{p} e^{- i β_{l} H_{M}} e^{- i γ_{l} H_{C}} | + ⟩^{\otimes n}

where $H_{M} = \sum_{i} X_{i}$ is the mixer Hamiltonian.

Step 1: Define the Problem Graph

python

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

# Define the 4-vertex ring graph
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}")

Step 2: Construct the Cost Hamiltonian

The cost Hamiltonian encodes the MaxCut objective. Each edge contributes a $- Z_{i} Z_{j}$ term (dropping the constant):

python

# Build the cost Hamiltonian from graph edges
# H_C = sum_{(i,j) in E} Z_i Z_j  (we minimize this)
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}")

Step 3: Build the QAOA Ansatz

python

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)  # Number of QAOA layers
    assert len(beta) == p, "gamma and beta must have the same length"

    prog = core.QProg()

    # Step A: Prepare uniform superposition
    for q in range(n_qubits):
        prog << core.H(q)

    # Step B: Apply p alternating layers
    for layer in range(p):
        # Cost unitary: e^{-i * gamma * H_C}
        # For each edge (i, j), apply RZZ(i, j, 2 * gamma)
        for i, j in edges:
            prog << core.RZZ(i, j, 2 * gamma[layer])

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

    return prog

# Test with a single layer
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")

Step 4: Define the Cost Function

python

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)

    # Compute expectation of cost Hamiltonian
    expectation = core.expval_pauli_operator(prog, cost_hamiltonian)

    return expectation

# Evaluate at initial parameters
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}")

Step 5: Optimize the QAOA Parameters

python

from scipy.optimize import minimize

# Run optimization with multiple random restarts
best_cost = float('inf')
best_params = None
n_restarts = 10

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

    # Minimize using 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:]}")

Step 6: Sample the Optimal Solution

python

# Build circuit at optimal parameters
optimal_gamma = best_params[:p].tolist()
optimal_beta = best_params[p:].tolist()
optimal_prog = build_qaoa_circuit(n_qubits, edges, optimal_gamma, optimal_beta)

# Sample from the circuit to get candidate solutions
n_shots = 1000
machine = core.CPUQVM()
machine.run(optimal_prog, shots=n_shots)
counts = machine.result().get_counts()

# Evaluate each measured bitstring's cut value
print("\nTop solutions by frequency:")
sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)

for bitstring, count in sorted_counts[:5]:
    # Compute the cut value for this bitstring
    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}")

# Find the best bitstring overall
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}")

Step 7: Increase QAOA Depth

Higher $p$ values improve the approximation ratio at the cost of more parameters:

python

# Run QAOA with p=2 layers for better approximation
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 Approximation Ratio

The approximation ratio measures how close the QAOA result is to the optimal solution:

α = \frac{⟨ H_{C} ⟩_{QAOA}}{MaxCut (G)}

python

# Compute approximation ratios
optimal_cut = n_edges  # For the 4-vertex ring, optimal cut = 4

# The cost Hamiltonian is sum of Z_i Z_j
# When all edges are cut, each Z_i Z_j = -1, so cost = -n_edges
# Approximation ratio: |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}")

# Theoretical guarantee: QAOA with p=1 achieves alpha >= 0.6924
# for 3-regular graphs (Farhi et al. 2014)
print(f"\nTheoretical p=1 lower bound (3-regular): 0.6924")

Extending to Larger Graphs

The same framework applies to larger graphs. Here is a helper for arbitrary graph definitions:

python

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.
    """
    # Build cost Hamiltonian
    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()

    # Sample the optimal circuit
    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()

    # Find best bitstring
    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,
    }

# Example: solve for a 5-vertex path graph
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}")

Next Steps

Hamiltonian & Pauli Operators - Define observables for VQE cost functions
Simulation - Run variational circuits on simulators
Quantum Information - Analyze variational state properties

Variational Quantum Circuits ​

What Are Variational Circuits? ​

The VQCircuit Module ​

Key Components ​

VQCircuit Class API ​

Constructor ​

Adding Gates ​

Setting Parameter Values ​

Computing Expectation Values ​

Computing Gradients ​

Computing Gradients and Expectation Together ​

Batch Evaluation with ResNGradients ​

API Method Summary ​

Variational Gates (VQGate) ​

Single-Qubit Parameterized Gates ​

Two-Qubit Parameterized Gates ​

Mixed Fixed and Variational Parameters ​

Parameter Management ​

Parameter Class ​

MutableParamType ​

Gradient Computation ​

Adjoint Differentiation ​

ResGradients - Single Parameter Set ​

ResNGradients - N Parameter Sets ​

ResGradientsAndExpectation - Combined Result ​

Mathematical Foundation ​

Parameterized Quantum State ​

Expectation Value ​

Gradient via Adjoint Differentiation ​

Parameter Shift Rule (Alternative) ​

VQA Workflow ​

Building a VQCircuit Step by Step ​

Step 1: Decide on the Ansatz Structure ​

Step 2: Import Modules and Define Constants ​

Step 3: Create the Parameter Object ​

Step 4: Build the Ansatz Circuit ​

Step 5: Define the Observable (Hamiltonian) ​

Step 6: Initialize Parameters ​

Step 7: Evaluate the Cost Function ​

Step 8: Run the Optimization Loop ​

Step 9: Verify the Result ​

Batch Parameter Evaluation ​

Why Batch Evaluation? ​

Batch Evaluation Workflow ​

Example: Evaluating a Parameter Grid ​

Computing Batch Gradients ​

Analyzing the Energy Landscape ​

Gradient Verification and Debugging ​

Numerical Gradient Check ​

Gradient Consistency Check ​

Common Debugging Checklist ​

Detecting Barren Plateaus ​

Logging Gradient Norms During Optimization ​

Performance Optimization ​

Choose the Right Differentiation Method ​

Reduce Circuit Depth ​

Exploit Hamiltonian Symmetry ​

Optimizer Selection Guide ​

Parameter Initialization Strategies ​

Caching Intermediate Results ​

Complete Example: VQE for H₂ ​

Interpreting VQE Results ​

Energy Accuracy ​

Optimal Parameter Analysis ​

State Composition ​

Convergence Behavior ​

Scaling to Larger Molecules ​

Variational Circuit Construction Patterns ​

Hardware-Efficient Ansatz ​

QAOA Ansatz ​

DiffMethod Reference ​

ADJOINT_DIFF ​

Gradient Result API Summary ​

ResGradients ​

ResNGradients ​

ResGradientsAndExpectation ​

QAOA Complete Example: MaxCut ​

Problem Definition ​

Example Graph ​

QAOA Circuit Structure ​

Variational Quantum Circuits

What Are Variational Circuits?

The VQCircuit Module

Key Components

VQCircuit Class API

Constructor

Adding Gates

Setting Parameter Values

Computing Expectation Values

Computing Gradients

Computing Gradients and Expectation Together

Batch Evaluation with ResNGradients

API Method Summary

Variational Gates (VQGate)

Single-Qubit Parameterized Gates

Two-Qubit Parameterized Gates

Mixed Fixed and Variational Parameters

Parameter Management

Parameter Class

MutableParamType

Gradient Computation

Adjoint Differentiation

ResGradients - Single Parameter Set

ResNGradients - N Parameter Sets

ResGradientsAndExpectation - Combined Result

Mathematical Foundation

Parameterized Quantum State

Expectation Value

Gradient via Adjoint Differentiation

Parameter Shift Rule (Alternative)

VQA Workflow

Building a VQCircuit Step by Step

Step 1: Decide on the Ansatz Structure

Step 2: Import Modules and Define Constants

Step 3: Create the Parameter Object

Step 4: Build the Ansatz Circuit

Step 5: Define the Observable (Hamiltonian)

Step 6: Initialize Parameters

Step 7: Evaluate the Cost Function

Step 8: Run the Optimization Loop

Step 9: Verify the Result

Batch Parameter Evaluation

Why Batch Evaluation?

Batch Evaluation Workflow

Example: Evaluating a Parameter Grid

Computing Batch Gradients

Analyzing the Energy Landscape

Gradient Verification and Debugging

Numerical Gradient Check

Gradient Consistency Check

Common Debugging Checklist

Detecting Barren Plateaus

Logging Gradient Norms During Optimization

Performance Optimization

Choose the Right Differentiation Method

Reduce Circuit Depth

Exploit Hamiltonian Symmetry

Optimizer Selection Guide

Parameter Initialization Strategies

Caching Intermediate Results

Complete Example: VQE for H₂

Interpreting VQE Results

Energy Accuracy

Optimal Parameter Analysis

State Composition

Convergence Behavior

Scaling to Larger Molecules

Variational Circuit Construction Patterns

Hardware-Efficient Ansatz

QAOA Ansatz

DiffMethod Reference

ADJOINT_DIFF

Gradient Result API Summary

ResGradients

ResNGradients

ResGradientsAndExpectation

QAOA Complete Example: MaxCut

Problem Definition

Example Graph

QAOA Circuit Structure