Benchmarking Noise-Resilient Optimizers for Variational Quantum Algorithms

Variational quantum algorithms (VQAs) [1] aim to extract useful results from noisy intermediate-scale quantum (NISQ) hardware [2]. In them, a quantum device evaluates expectation values of a parametrized circuit, and a classical optimizer uses those estimates to update the circuit parameters iteratively.

Because every expectation value is estimated from a finite number of measurements, or shots, the optimizer never sees the objective function exactly. The resulting shot noise can obscure small improvements and becomes increasingly costly to suppress. Bärligea et al. [3], for example, showed that this uncertainty creates serious scaling challenges for variational quantum optimization. The classical optimizer must therefore make progress using noisy information without requiring an impractical number of circuit measurements.

This tutorial investigates how several optimizer families behave as both the number of qubits and the measurement precision change. We first introduce two explicitly noise-aware optimizer proposals, HOPSO [4] and ANATRA [5], together with ROSALIN [6] and standard gradient-based and derivative-free baselines. We then benchmark them with Tequila on Sherrington-Kirkpatrick spin-glass instances.

Background

Harmonic Oscillator-Based Particle Swarm Optimization (HOPSO)

Harmonic oscillator-based particle swarm optimization (HOPSO) replaces the algebraic velocity update of conventional particle swarm optimization with the motion of a damped harmonic oscillator [7]. Its core idea is that each particle oscillates around an attractor \(a_{j,d}\) obtained from its personal best \(p_{j,d}\) and the swarm’s global best \(g_d\). In one dimension, the idealized trajectory has the form

\[ x(t)=A_0 e^{-\lambda t}\cos(\omega t+\theta)+a_{j,d}, \]

where \(A_0\) is the initial amplitude, \(\lambda\) is the damping rate, \(\omega\) is the angular frequency, and \(\theta\) is the phase. Damping gradually shifts the behavior from broad exploration to local refinement.

For variational circuits, the search coordinates are rotation angles and are therefore periodic. The VQA adaptation of HOPSO [4] treats personal and global best positions on a circle and prevents the oscillation amplitude from vanishing prematurely. Defining

\[ r_{j,d}=|p_{j,d}-g_d|\bmod 2\pi, \]

the amplitude floor is based on the shorter circular distance,

\[ A_{\mathrm{th},j,d} =\frac{1}{2}\min\!\left(r_{j,d},2\pi-r_{j,d}\right)m, \]

where \(m\) controls the minimum remaining oscillation amplitude.

Circular attractor computation. A linear average of \(10^\circ\) and \(350^\circ\) is \(180^\circ\), although the two angles are only \(20^\circ\) apart on the circle. HOPSO instead uses the shorter wrapped displacement before forming the weighted attractor. The following excerpt shows the central periodic update used in the implementation.

# HOPSO Velocity & Position Update Snippet
for i in range(population_size):
    # Circular/periodic shortest deltas
    d_pb = (P[i] - X[i] + np.pi) % (2 * np.pi) - np.pi
    d_gb = (G - X[i] + np.pi) % (2 * np.pi) - np.pi

    # Calculate attractor
    attractor_delta = (c1 * d_pb + c2 * d_gb) / (c1 + c2)
    attractor = X[i] + attractor_delta

    # Update velocity via Damped Harmonic Oscillator physics
    amp = amplitude_factor * np.abs(d_pb - d_gb)
    V[i] = inertia * damping * V[i] + amp * (attractor - X[i])
    X[i] = X[i] + V[i]

Adaptive Shot Allocation and Trust-Region Optimization

Noise-aware optimizers can intervene at different points in the optimization loop. ROSALIN changes how a finite measurement budget is distributed across Hamiltonian terms and gradient components. ANATRA instead changes how a derivative-free local model is constructed when objective values are uncertain.

ROSALIN

In VQE, the energy is the expectation value of a Hamiltonian written as a sum of Pauli operators,

\[ H=\sum_i c_iP_i. \]

Measuring every term equally can be wasteful when the coefficients \(c_i\) vary strongly. ROSALIN (Random Operator Sampling for Adaptive Learning with Individual Number of shots) samples Hamiltonian terms with probabilities proportional to \(|c_i|\), reweights their contributions to retain an unbiased estimator, and combines this operator sampling with adaptive shot allocation for stochastic-gradient updates [6].

The separate tutorial Testing Shot-Frugal Optimizers for Variational Quantum Algorithms studies adaptive shot-frugal optimization under a total-shot budget. ROSALIN improved steadily with the available measurements and remained competitive there.

ANATRA

The Adaptive Noisy Model-Based Trust-Region Algorithm (ANATRA) [5] is a derivative-free optimizer designed for noisy quantum loss functions. Instead of estimating gradients directly, it constructs a local quadratic surrogate from energy values evaluated near the current parameter iterate \(\theta_k\). The next step minimizes this surrogate within a trust region of radius \(\Delta_k\).

ANATRA achieves noise resilience through two main mechanisms:

• Surrogate smoothing: Fitting the quadratic model to several sampled points reduces the influence of any single noisy objective value.
• Adapative sampling radius: A standard trust-region method shrinks \(\Delta_k\) after unsuccessful steps. Under shot noise, points that are too close have energy differences that are indistinguishable from evaluation noise. ANATRA therefore decouples the sampling radius from the step radius and enforces a minimum sampling distance related to the noise level.

The excerpt below shows how the implementation computes this minimum radius and requests a better interpolation geometry when the sampled energy variation is dominated by noise.

# ANATRA Noise-Aware Sampling Radius & Geometry Check 


# 1. Compute Lipschitz constant & noise-aware minimum sampling radius
if valid_geometry:
    L_new = np.max(np.linalg.eigvals(H_hessian))
    L = np.maximum(L_new, 1.0)
min_sample_delta = np.sqrt(2.0 * r * epsilon / L)

# 2. Decouple sampling radius (delta_bar) from trust-region step radius (delta_k)
delta_bar = max(delta_k, min_sample_delta)

# 3. Test if function variation across sample set is dominated by evaluation noise
noise_dominated = np.all(np.abs(hF[Mind] - hF[xk_center]) < epsilon)
if noise_dominated:
    # Trigger geometry improvement subroutine (Improving poisedness of X)
    Mdir, mp, valid, Cres, Gres, Hres, Mind = formquad_lagrange(
        X[0:nf + 1, :], Res[0:nf + 1, :], delta_bar, xk_center, np_max, Par, 0, Mind
    )

2.3 Standard Baseline Optimizers

To benchmark the specialized noise-resilient algorithms discussed above, we compare them to well-established classical optimizers. ADAM and stochastic gradient descent (SGD) are available through Tequila, while COBYLA and Powell’s method use SciPy implementations.

• ADAM and SGD update parameters using estimated gradients. SGD follows the current gradient estimate directly, whereas ADAM combines moving averages of gradients and squared gradients to introduce momentum and coordinate-dependent learning rates.
• COBYLA and Powell’s method use objective values without explicit gradients. COBYLA builds local linear models inside a trust region, while Powell’s method performs sequential line searches along a changing set of directions.

As these methods were not designed specifically for VQAs, they provide useful reference points for asking whether a noise-aware construction yields a practical advantage in the benchmark below.

Experimental Setup

The experiment is designed to separate two effects that are easily conflated: a larger circuit creates a higher-dimensional optimization problem, while a larger shot setting makes each sampled objective value more precise. We therefore vary the system size from \(N=3\) to \(N=8\) qubits and the nominal shot setting from \(S=10\) to \(S=100{,}000\). Every optimizer is tested on random Sherrington-Kirkpatrick (SK) spin glass Hamiltonians using the same hardware-efficient ansatz, built from \(R_y\) and \(R_z\) rotations and CNOT entangling layers.

For each pair \((N,S)\) of system size and shot setting, we generate ten problem instances and optimize each from ten seeded initial parameter vectors. The resulting \(100\) runs sample both instance-to-instance variation and sensitivity to initialization. Reusing these seeds across optimizers keeps the comparison aligned, and the exact minimum of each small SK instance provides the reference energy.

The horizontal axis of the following heatmaps should be interpreted carefully. For fixed-shot methods, \(S\) controls the measurements used for an expectation-value estimate; adaptive methods may distribute measurements differently. Moreover, optimizers use different numbers of objective or gradient evaluations. The figures therefore compare final solution quality at a nominal sampling setting, not total quantum measurement cost. A direct resource-efficiency comparison would instead place every method under the same total circuit-shot budget.

The compact worker below shows how each seeded trial constructs the problem, initializes the ansatz, dispatches the selected optimizer, and stores the optimized and exact energies.

def worker_task(args_tuple):
    """ Executes a single optimization run """
    n, nshots, p_idx, r_idx, optimizer, ansatz, n_layers = args_tuple
    
    # 1. Establish deterministic seeds
    seed_problem = p_idx
    seed_run = 1000 * p_idx + r_idx
    
    # 2. Build the random SK Model and Tequila Ansatz
    prob = IsingProblem.create_random_instance(n=n, which="SK", seed=seed_problem)
    ref_min = prob.get_min_cost_value() # Calculate the exact analytical minimum
    O, vars_ = build_ising_random(n=n, seed=seed_problem, n_layers=n_layers, ansatz=ansatz)

    # 3. Initialize parameters normally around 0
    rng = np.random.default_rng(seed_run)
    init = {v: float(rng.normal(0.0, 0.1)) for v in vars_}

    # 4. Dispatch to the corresponding Optimizer
    if optimizer.lower() == "hopso":
        kwargs = {"population_size": 10, "maxiter": 30, "seed": seed_run, "samples": nshots, "c1": 1.0, "c2": 2.0}
        res = run_hopso(O, vars_, init, **kwargs)
        
    elif optimizer.lower() in ["rosalin", "anatra"]:
        res = run_custom_optimizer(O, optimizer, vars_, init, extra_kwargs={"s_min": 10, "samples": nshots})
        
    return {
        "found_min": res["energy"],
        "ref_min": ref_min,
        "nfev": res.get("fes", 0)
    }

Results and Analysis

For every run, we compare the optimized energy \(E_{\mathrm{found}}\) with the exact ground-state energy \(E_{\mathrm{exact}}\) through the relative error

\[ \text{Relative Error} =\frac{|E_{\mathrm{found}}-E_{\mathrm{exact}}|}{|E_{\mathrm{exact}}|}. \]

Each heatmap reports the mean over the \(100\) runs for one optimizer, where the horizontal axis gives the nominal shot setting, the vertical axis the number of qubits, and the logarithmic color scale moves from low error in green to high error in red. These two-dimensional plots let us ask separately whether an optimizer benefits from more precise evaluations and how its error changes across the tested problem sizes.

Gradient-Based Optimizers (ADAM, SGD)

ADAM and SGD produced comparatively low mean errors throughout this experiment, but the dependence on shots is not the usual monotonic decrease one might expect from more precise objective estimates. ADAM is especially insensitive to the shot setting: for \(N=8\), its mean relative error stays between \(0.017\) and \(0.036\) over the entire range from \(10\) to \(100{,}000\) shots. The corresponding SGD data show a similarly weak dependence, with several system sizes even giving lower errors at the smallest shot settings.

These results show that the tested gradient updates continue to make progress under strong sampling noise.

Gradient-Free Direct Search (COBYLA, Powell)

COBYLA and Powell’s method display a different pattern. Their mean errors fall markedly as the shot setting increases, consistent with direct-search steps benefiting from more precise comparisons between nearby objective values. For COBYLA at \(N=3\), the mean relative error decreases from \(0.451\) at \(10\) shots to \(0.001\) at \(100{,}000\) shots. At \(N=8\), the corresponding decrease is from \(0.834\) to \(0.093\), so the larger instances remain evidently harder within the same optimizer configuration.

The heatmap reveals a strong association between sampling precision and final error for this method.

4.3 ROSALIN and HOPSO

In this benchmark, ROSALIN remains close to relative error \(1\) across all tested system sizes and shot settings. In particular, increasing the nominal shot setting produces no systematic improvement in the heatmap.

This result should not be read as a general failure of adaptive operator sampling. In the related shot-frugal tutorial linked above, ROSALIN improved with the total measurement budget and remained competitive with the other adaptive optimizers. The contrast indicates that the present ROSALIN result needs to be understood together with its implementation, the commuting SK Hamiltonian, and the different meaning of the shot parameter before drawing a broader conclusion.

HOPSO also yields substantially larger errors than the gradient-based methods in this setup. Across the full grid, its mean relative errors range from \(0.162\) to \(0.714\). Increasing the shot setting does not reduce the error; for most system sizes, the error first rises and then plateaus.

The absence of improvement at high shot counts means that measurement imprecision alone cannot explain HOPSO’s performance here. A follow-up study would need to vary its population size, iteration budget, and dynamical hyperparameters before attributing the observed plateau to the optimizer itself.

Model-Based Trust-Region Optimization (ANATRA)

ANATRA shows the clearest improvement with increasing sampling precision. At \(N=8\), its mean relative error decreases from \(0.477\) at \(10\) shots to \(0.027\) at \(100{,}000\) shots. Similar transitions from high to low error appear for the other system sizes, although the shot setting required to reach the low-error region grows with \(N\).

Within the tested range, ANATRA therefore approaches the low errors of ADAM once sufficiently precise objective values are available. At low shot settings, however, its local surrogate models do not achieve the same solution quality.

Summary and Limitations

For these SK instances, ansatzes, hyperparameters, and evaluation limits, ADAM has the lowest errors with the weakest dependence on the nominal shot setting. ANATRA benefits strongly from additional sampling and becomes competitive at high shot counts. COBYLA and Powell also improve with more precise evaluations, while the tested ROSALIN and HOPSO configurations show little or no benefit from increasing the shot parameter.

Lastly note that these observations describe a finite benchmark over \(N=3\) to \(8\) qubits, so they should not be interpreted as asymptotic scaling results or evidence that one optimizer is generally superior. The heatmaps report means without displaying the run-to-run distributions, and the nominal shot setting does not equal the total measurement cost of an optimization. A stronger comparison would report uncertainty across instances, tune each method under a common protocol, and compare final error at equal total circuit-shot budgets.

References

[1]

M. Cerezo et al., “Variational quantum algorithms,” Nature Reviews Physics, vol. 3, no. 9, pp. 625–644, Aug. 2021, doi: 10.1038/s42254-021-00348-9.

[2]

J. Preskill, “Quantum computing in the NISQ era and beyond,” Quantum, vol. 2, p. 79, Aug. 2018, doi: 10.22331/q-2018-08-06-79.

[3]

A. Bärligea, B. Poggel, and J. M. Lorenz, “Scalability challenges in variational quantum optimization under stochastic noise,” Phys. Rev. A, vol. 112, no. 3, p. 032407, Sept. 2025, doi: 10.1103/rgyh-8xw8.

[4]

I. A. Mohammad, Y. Chernyak, and M. Plesch, “HOPSO: A robust classical optimizer for VQE.” arXiv, 2025. doi: 10.48550/ARXIV.2508.13651.

[5]

J. Larson, M. Menickelly, and J. Shi, “A novel noise-aware classical optimizer for variational quantum algorithms,” INFORMS Journal on Computing, vol. 37, no. 1, pp. 63–85, Jan. 2025, doi: 10.1287/ijoc.2024.0578.

[6]

A. Arrasmith, L. Cincio, R. D. Somma, and P. J. Coles, “Operator sampling for shot-frugal optimization in variational algorithms.” 2020. Available: https://arxiv.org/abs/2004.06252

[7]

Y. Chernyak, I. A. Mohammad, N. Masnicak, M. Pivoluska, and M. Plesch, “Harmonic oscillator based particle swarm optimization,” PLoS One, vol. 20, no. 6, p. e0326173, June 2025, doi: 10.1371/journal.pone.0326173.