Adaptive Error Budget Allocation for Fault-Tolerant Quantum Resource Estimation: A Metaheuristic ApproachPubID: pubid: Accepted for publication at IEEE Quantum Week 2026.

Fault-tolerant quantum computation requires encoding logical qubits into a larger number of physical qubits using a quantum error-correcting code, such that errors can be detected and corrected without destroying the encoded quantum information[9, 12]. The surface code[8, 4] has emerged as the leading candidate architecture for fault-tolerant quantum computation, primarily because it requires only nearest-neighbor interactions on a two-dimensional qubit lattice while maintaining a comparatively high error-correction threshold of approximately 1% per physical gate operation [24]. This combination of architectural simplicity and noise tolerance makes it particularly well-suited for practical implementations of quantum error correction and, consequently, for the error budget optimization framework developed in this work.

To execute a fault-tolerant quantum algorithm within an allowed overall error budget $\varepsilon$, the error budget must account for all sources of error throughout the computation. Following the framework established by Beverland et al. [2] and implemented in the AQRE [27], the total error budget is partitioned into three components:

where $\varepsilon_{\text{logical}}$ denotes the portion of the error budget reserved for logical qubit operations; $\varepsilon_{\text{tstates}}$ represents the portion allocated to producing reliable magic states through distillation; and $\varepsilon_{\text{rotations}}$ corresponds to the portion assigned to approximating arbitrary rotations through Clifford+T decompositions[27, 18]. These allocations are interdependent within the resource estimation pipeline, requiring careful balancing to minimize hardware overhead while maintaining the desired computational reliability[7].

The allocation of the error budget is based on a constrained optimization problem on the probability simplex. Given a total error budget $\varepsilon_{\text{total}}$ representing the maximum tolerable failure probability for the quantum algorithm, we seek an allocation vector $\mathbf{s}=(s_{L},s_{T},s_{R})$ that distributes error fractions among three functionally distinct subsystems: logical operations ($L$), distillation of the T-state ($T$) and rotation synthesis ($R$). The allocation is constrained to the $\varepsilon$-interior of the 2-simplex:

where $\varepsilon$ ensures that no subsystem is assigned a degenerate allocation that would cause the resource estimator to produce infeasible configurations. Under this constraint, $\mathbf{S}_{\varepsilon}$ forms a truncated version of the standard 2-simplex, resulting in a compact two-dimensional feasible region over which optimization is performed, with each subsystem receiving a corresponding fraction of the total error. The AQRE then acts as a black-box mapping $\mathbf{s}\mapsto(Q(\mathbf{s}),R(\mathbf{s}))$, where the resulting resource estimates vary in a piecewise manner due to the underlying discrete design choices.

The optimization problem formulated by the resource estimation model is inherently non-convex, discontinuous, and structurally disordered as the resource estimator takes continuous allocation input to physical resources through discrete internal decisions, such as integer-valued code distance choice and finite T-factory design catalogs. Minimal changes in the error allocation vector can induce sudden change of the regime, causing step-wise variations in the number of physical qubits and execution time. As a result, the objective function exhibits a rugged, multimodal surface with a large number of local minima and no cost separations. These properties have a negative impact on the use of gradient-based optimization algorithms, which assume smoothness and differentiability at the edges of regimes that is not met, and restrict the utility of deterministic local search algorithms that do not have escape mechanisms of suboptimal basins of attraction. Additionally, the simplex constraint of allocation vector also restricts possible search space and thus any optimization algorithm must ensure that its search remains budget feasible at each step.

The resource estimator returns two competing quantities for any allocation: physical qubit count $Q(\mathbf{s})$ and algorithmic runtime $R(\mathbf{s})$. Optimizing either one in isolation is insufficient, as reducing qubits often increases runtime and vice versa. We therefore define a composite cost function that jointly captures space–time complexity:

where $\alpha\in[0,1]$ controls the trade-off between qubit count and runtime. This formulation corresponds to a geometric mean objective, ensuring scale-invariant balancing between competing resource dimensions. In practice, we use $\alpha=0.5$, yielding:

which is monotonically equivalent to minimizing the space–time volume $Q\cdot R$, the standard metric in fault-tolerant resource estimation. The optimization problem is therefore:

where $\mathbf{s}^{*}$ denotes the allocation that achieves the minimum value of the estimator-induced space-time objective within the feasible simplex $\mathcal{S}_{\varepsilon}$.

At each iteration, the algorithm generates a candidate allocation by perturbing the current solution. The perturbation must be large enough to escape local plateaus in the early stages of the search, yet small enough to refine promising allocations in the later stages. We achieve this by tying the step size to both the current temperature and the fraction of iterations completed.

At iteration $k$, a neighbor $\mathbf{s}^{\prime}$ is generated from the current allocation $\mathbf{s}^{(k)}$ by:

where each perturbation component is drawn independently as:

The step size $\eta_{k}$ starts large and shrinks as the search progresses, ensuring broad exploration early and focused refinement later. The projection operator $\text{Proj}_{\mathcal{S}_{\varepsilon}}$ maps the perturbed vector back onto the feasible simplex $\mathcal{S}_{\varepsilon}$ by clamping each component to a minimum of $\varepsilon$ and renormalizing to unit sum. This guarantees that every candidate allocation is feasible without requiring any auxiliary constraint-handling mechanism.

A deterministic acceptance of all candidate solutions is undesirable. Although cost reducing moves are always beneficial, the algorithm must occasionally permit higher cost transitions to escape local minima and ensure adequate exploration of the solution space. The standard Metropolis criterion $\exp(-\Delta/T)$ [17] controls this trade-off through temperature $T$. However, this standard form assumes that cost magnitudes are comparable across problem instances. In our setting, a small circuit may cost $10^{6}$, while a large circuit reaches $10^{15}$, making a fixed temperature meaningless across the benchmark suite.

We resolve this by normalizing the cost difference by the current cost value. Let $\Delta=\mathcal{J}(\mathbf{s}^{\prime})-\mathcal{J}(\mathbf{s}^{(k)})$ denote the cost difference. The acceptance rule is:

The denominator $T_{k}\cdot\mathcal{J}(\mathbf{s}^{(k)})$ means the algorithm evaluates the relative cost increase $\Delta/\mathcal{J}(\mathbf{s}^{(k)})$ rather than the absolute value of $\Delta$. A 5% cost increase is treated the same way whether the circuit costs $10^{6}$ or $10^{15}$, producing consistent exploration behavior across the entire benchmark suite.

The temperature $T$ controls how willing the algorithm is to accept worse solutions. High temperatures encourage exploration; low temperatures enforce exploitation. The rate at which $T$ decreases therefore determines how quickly the algorithm transitions between these two regimes.

A fixed cooling rate treats all phases identically, regardless of whether the algorithm is actively finding improvements or stagnating. Adaptive cooling schedules that adjust based on search progress have been proposed to improve efficiency [11]. We instead adapt the rate based on the recent improvement rate $\rho_{k}$, defined as the relative cost decrease of the most recent accepted improving move. The temperature update is:

where the cooling factor $\gamma$ responds to the optimization dynamics:

When the algorithm is finding large improvements ($\rho_{k}>\tau_{\text{high}}$), cooling slows down to give the search more time in the productive region. When improvement has stalled ($\rho_{k}<\tau_{\text{low}}$), cooling speeds up to push the algorithm toward convergence rather than wasting iterations. To prevent the temperature from dropping too low too early, periodic reheating is applied:

This brief temperature increase re-enables uphill moves at regular intervals, allowing the algorithm to escape plateaus it may have committed to prematurely.

A single SA run will converge to whatever local minimum lies in the basin of its starting point. If the starting point happens to be near a poor local minimum, the final solution will be poor regardless of how well the cooling schedule is tuned. To reduce this dependence on initialization, we run SA multiple times from different starting points and keep the best result.

The starting points are generated using Latin Hypercube Sampling (LHS)[16], which partitions each dimension of the unit hypercube into equal-probability intervals and places exactly one sample per interval:

The samples are then mapped onto the probability simplex through exponential normalization:

This produces starting points that are spread uniformly across the allocation space, unlike Dirichlet sampling which concentrates points near the simplex center. The globally best allocation found across all $n$ restarts is returned as the final optimized solution.

Every restart does not need to run for the full iteration budget. Once the search has stabilized and the cost is no longer changing, continuing to iterate wastes computation that could be used by other restarts. We detect this stabilization by monitoring the recent history of accepted cost values. Let $W$ denote a sliding window over the most recent accepted costs. Convergence is declared when the relative variation within this window falls below a threshold $\delta$:

When this condition is met, the accepted costs are effectively constant, indicating that the search has settled near a local minimum. The current restart terminates early, and the saved iterations are available for subsequent restarts that may still be exploring productive regions.

To quantify the benefit of optimized allocation, all results are compared against the uniform distribution baseline, which is the default setting of the AQRE:

Performance improvement is measured as the relative reduction in the space-time cost:

A positive value indicates that SA optimized allocation requires fewer physical resources than the uniform baseline, with larger values corresponding to greater resource savings.

A swarm of $N$ particles searches over generations $k=1,\dots,K$. Let $\mathbf{x}_{i}^{(k)}\in\mathbb{R}^{D}$ denote the position of particle $i$ at generation $k$, $\mathbf{p}_{i}$ its personal best position, and $\mathbf{g}$ the global best across the swarm. For each coordinate $d$, a stochastic attractor is formed as the random convex combination

The position update, derived from the inverse cumulative distribution of a particle bound in a Delta potential well [26], is

where the sign is chosen with equal probability and $\bar{p}_{d}=N^{-1}\sum_{i=1}^{N}p_{i,d}$ is the mean best position. The $\ln(1/u)$ term samples a displacement from a heavy-tailed distribution whose spatial reach is scaled by the distance from the particle to the mean best. This allows a particle to traverse an entire plateau and cross the adjacent code-distance boundary in a single step, without requiring gradient information from the piecewise-constant surface described in Section III-B.

The contraction–expansion coefficient $\beta_{k}$ is decreased linearly from $\beta_{\max}=1.0$ to $\beta_{\min}=0.5$ across the $K$ generations and plays a role analogous to the temperature in SA. Large values permit long-range tunneling during early exploration, while small values confine particles near their attractors during late refinement. Unlike the adaptive cooling in Eq. (9), this schedule is deterministic and requires no threshold parameters.

SA handles the simplex constraint through projection, where the perturbed vector in Eq. (6) is clamped to a minimum of $\varepsilon$ and renormalized. This introduces boundary bias because perturbations near the simplex edges are asymmetrically truncated. We instead map the open 2-simplex onto unconstrained $\mathbb{R}^{2}$ using the additive log-ratio (ALR) transform[1], taking $s_{L}$ as the reference component.

The inverse recovers the allocation via softmax over $\mathbf{z}=(0,\,y_{1},\,y_{2})$.

This mapping is a diffeomorphism between the open simplex and $\mathbb{R}^{2}$. The QPSO updates in Eqs. (17)–(18) operate entirely in log-ratio space, and the simplex constraint is satisfied exactly by construction at every evaluation. No projection, clamping, or renormalization is required at any step of the algorithm.

After the inverse transform, a minimum fraction guard ensures each component receives at least $\varepsilon_{i}=0.05$ of the total budget. The guard reserves $\varepsilon_{i}$ per component and distributes the remaining $1-3\varepsilon_{i}$ proportionally to the transformed fractions. This is not an algorithmic constraint, but a physical one imposed by the resource estimator’s code-distance ceiling of $d\leq 50$. It is applied as a post-processing step and does not affect the unconstrained geometry in which QPSO operates.

Each of $n$ independent restarts initializes particle positions by drawing from the symmetric Dirichlet distribution $\mathrm{Dir}(1,1,1)$, which is uniform on the simplex, and mapping the samples to the log-ratio space via Eq. (19). Within each restart, convergence is declared when the global best cost remains unchanged for 30 consecutive generations, after which the restart terminates and the saved evaluations are available to subsequent restarts. The best global allocation in all $n$ restarts is returned as the final solution.

Figure 2 compares the distribution of space-time improvement $(\mathcal{J}(\mathbf{s}_{\text{uniform}})-\mathcal{J}(\mathbf{s}^{*}))/\mathcal{J}(\mathbf{s}_{\text{uniform}})$ across 31 circuit families for SA, QPSO, and GT method. All per-family distributions and scatter plots in Figures 2–4 are computed from a single representative seed (seed 60) for SA and QPSO; multi-seed statistics are reported in Table I. The observed family-level trends are consistent across all three approaches, with the same circuit families consistently ranked at the top and bottom. This consistency suggests that the primary source of variation arises from the circuit-dependent error profile, rather than from optimizer-specific effects. Furthermore, both SA and QPSO increase the family-level averages relative to GT across most benchmarks, indicating that direct metaheuristic search can exploit additional opportunities for improvement while preserving the same underlying structural pattern.

Figure 3 evaluates the extent to which the magnitude of improvement can be anticipated from three pre-optimization circuit descriptors: the number of logical qubits $n_{q}$, the T-gate count $N_{T}$, and the rotation count $N_{\mathrm{rot}}$. For each descriptor, improvement is plotted on a logarithmic scale for SA, QPSO, and GT, with both Pearson correlation coefficient $r$ and the Spearman rank correlation $\rho$ reported. The results indicate that the features of the structural circuit are only weakly negatively correlated with the improvement in all methods. Specifically, for logical qubit count, the correlations range from $(r,\rho)=(-0.09,-0.08)$ for GT to $(-0.15,-0.15)$ for QPSO. For T-gate count, the correlations are similarly small, ranging from $(-0.18,-0.15)$ to $(-0.19,-0.18)$. Rotation count exhibits the strongest association among the three descriptors, but the effect remains limited, with correlations spanning approximately $(-0.24,-0.23)$ to $(-0.28,-0.21)$. These findings suggest that while larger or more gate-intensive circuits may experience marginally reduced improvement from budget reallocation, the overall association is consistent across SA, QPSO, and GT.

The observed weak dependence is itself a significant finding. If simple structural descriptors are strong predictors of optimization gain, one would expect to observe pronounced monotonic trends and a tighter clustering of data points around the regression lines in each panel. However, the data for all three methods exhibit substantial dispersion. Both small and large circuits are represented across the full range of improvement, as are circuits with moderate and high T-gate or rotation counts. The marginally stronger negative correlation for rotation count indicates that rotation synthesis exerts a somewhat greater influence on the optimization landscape compared to qubit count or T-gate count, but this effect remains insufficient to serve as a robust screening criterion.

Across all three panels, the most prominent feature is the pronounced bimodal distribution of improvement values, rather than any regression trend. For each of the descriptors, namely logical qubit count, T-gate count, and rotation count, the data are separated into two distinct bands: a high-improvement band spanning approximately 40% to 99%, and a low-improvement band from 0% to 15%, with few instances in the intermediate range. This bimodality is consistent across SA, QPSO, and GT. Both small circuits (fewer than ten logical qubits) and circuits with $10^{4}$ or more rotations appear in both bands. This separation reflects the discrete nature of the resource-estimation process: a reallocation either crosses a threshold that reduces code distance or synthesis overhead, leading to a substantial reduction in resource cost, or it does not, resulting in only marginal improvement. Intermediate cases are rare because code distance is integer-valued, and changes to it induce stepwise shifts in resource requirements.

As a result, no simple screening criterion based solely on circuit dimensions can reliably identify which instances will achieve substantial improvement. For instance, a circuit with $10^{5}$ rotations may belong to either the high or low improvement band, depending on the interplay among the three error mechanisms that contribute to the total resource cost. This variability is not captured by any single scalar measure of circuit size. The data across all three panels indicate that structural circuit features alone provide limited predictive power for optimization gain. The consistency of this result across SA, QPSO, and GT supports the conclusion that the primary factor governing improvement is the interaction between the circuit’s error profile and the estimator, rather than the circuit’s absolute scale. This observation motivates the development of predictors that directly quantify allocation imbalance, rather than relying on pre-optimization circuit counts.

The preceding analysis establishes that improvement is governed by the asymmetry of the error profile rather than by circuit scale, but does not provide a single quantitative measure of this asymmetry. To this end, we compute the Gini coefficient of the optimized allocation $\mathbf{s}^{*}$, a standard measure of distributional inequality that equals zero for the uniform allocation and approaches $(d-1)/d$ when the budget is concentrated on a single mechanism.

Figure 4 presents the relationship between the Gini coefficient of the optimized error budget allocation and the resulting metric improvement. Circuits within the same MQT benchmark family share a common gate topology and differ only in qubit count. Consequently, their transpiled gate structures, error budget compositions, and optimal allocations are inherently correlated. Treating these correlated instances as independent observations would artificially inflate the effective sample size and exaggerate the statistical significance. Both quantities are therefore aggregated to family-level means ($k=31$ families), ensuring that each data point represents a structurally distinct circuit architecture and that the reported correlations reflect genuine between family variation. The structural analysis in Subsection VII-B retains instance-level reporting because the primary interest there is the scaling of improvement with circuit size within each family, a within family effect that aggregation would obscure. The Gini analysis targets a between family effect, namely whether architecturally distinct circuit families with more asymmetric error budgets achieve systematically greater improvement, and family-level aggregation is the appropriate unit of analysis for this comparison.

The resulting correlations are strong and consistent across all three methods. SA yields $r=0.90$ and $\rho=0.84$, while QPSO yields $r=0.91$ and $\rho=0.78$, and GT yields $r=0.78$ and $\rho=0.74$, all with $p\ll 0.001$. These values substantially exceed every structural predictor examined in Subsection VII-B and hold regardless of the underlying optimization strategy. The regression lines for SA and QPSO nearly overlap across the entire Gini range, with closely aligned slopes and broadly overlapping confidence intervals, while GT follows the same positive trend with a somewhat shallower slope. This indicates that the relationship between budget inequality and improvement is primarily determined by circuit architecture, while SA and QPSO recover somewhat larger gains from the same underlying asymmetry.

The attainable improvement increases systematically over the observed Gini range of approximately $0.2$ to $0.55$, rising from $20$–$30\%$ for families such as Grover, Shor, and the integer multiplier to $70$–$80\%$ for highly asymmetric families such as the CDKM ripple carry adder, full adder, and modular adder. This progression is visible for SA, QPSO, and GT alike, although the SA and QPSO trends lie consistently above GT across most of the observed range. The trend is approximately linear throughout, with no evidence of saturation. This progression has a direct interpretation. When the optimal allocation is relatively balanced, the uniform baseline already lies close to the cost minimum on the allocation simplex $\mathcal{S}_{\varepsilon}$, leaving limited room for reduction in $\mathcal{J}(\mathbf{s})$. When the optimal allocation is strongly unequal, the uniform baseline assigns substantial budget to inactive pathways, and the optimizer eliminates this misallocation by redirecting budget toward the active mechanisms, triggering code distance reductions that produce multiplicative savings in the space-time cost.

The Gini coefficient is computed from the optimized allocation $\mathbf{s}^{*}$ and therefore serves as a post hoc diagnostic rather than a predictive feature. It cannot be evaluated before optimization and does not replace the optimizer. Its value lies in providing a compact and interpretable explanation for the large cross family variation in improvement observed under SA, QPSO and GT, and in confirming that the gains reported in Table I are systematic consequences of error profile asymmetry rather than artifacts of any particular search procedure.

Table I compares the present results with the two most relevant prior approaches to error budget optimization for fault-tolerant resource estimation, namely the supervised method of Forster et al. [7] and the GT method of Asif et al. [23]. The comparison with the GT method is direct, as both approaches are evaluated on the same MQT Bench suite using the same uniform allocation baseline. In contrast, the comparison with the supervised method is not a circuit-by-circuit match, since the benchmark suite has expanded since that study and now includes additional circuit families. However, this comparison remains a meaningful point of reference because both studies address the same error budget allocation task within the same benchmark framework and report improvements relative to the same default baseline.

Within this comparative framework, Table I demonstrates both an increase in average improvement and a significant enhancement in computational efficiency. The supervised method achieves an average improvement of 15.6%, whereas the GT method increases this figure to 30.22% on the 433-circuit suite. The present SA and QPSO frameworks further advance these results, attaining average improvements of $35.18\pm 0.11\%$ and $33.84\pm 0.08\%$, respectively. Notably, these gains are accompanied by a reduction in total runtime for the full benchmark suite from approximately 287 minutes with the GT method to approximately 20-25 minutes. These findings indicate that the observed resource savings are achieved without incurring additional optimization overhead.

Timing results play a pivotal role in software workflows for fault-tolerant resource estimation, as the practical value of allocation quality depends on the timely delivery of solutions within iterative design and compilation flows. Our empirical evaluation demonstrates that the proposed framework advances the trade-off between solution quality and computational efficiency, yielding greater average improvements than the same-suite GT approach while reducing runtime by more than an order of magnitude. In contrast to supervised methods, our framework obviates the need for offline dataset construction, model training, or prior circuit exposure, thereby streamlining the optimization process and increasing its suitability for estimator-driven workflows.

The close correspondence between SA and QPSO further substantiates this interpretation. Their mean improvements differ by only 1.34 percentage points, and both achieve the same maximum reduction of 97.81%, completing the full benchmark suite within a narrow runtime interval. This consistency indicates that the observed gains are attributable to the underlying structure of the allocation problem rather than reliance on a particular optimizer, and supports the use of adaptive error-budget allocation as a practical optimization layer within fault-tolerant resource estimation workflows.