request: automatically find bounds in error-budgeting

[nelimee]

Motivation

At the moment, users are expected to provide sensible bounds for the noise parameters when doing error budgeting. See for example the presentation of the error budgeting feature on the Deltakit website:

from deltakit.explorer.analysis.error_budget import get_error_budget, SamplingParameters

# Among other things, the bounds for each noise parameters of the noise model
# are a required input.
bounds: list[tuple[float, float]] = [(1e-3, 2e-2)]
results = get_error_budget(
    simple_noise_model, P, num_rounds_per_distance, bounds,
    sampling_parameters=SamplingParameters(max_shots=500_000)
)

These bounds are used in error-budgeting to:

1. Discretise the interval formed by each bounds using a DiscretisationStrategy.
2. Compute $\Lambda$ and its standard deviation $\sigma_\Lambda$ on each of the points obtained in step 1.
3. Fit a degree $d \geqslant 1$ polynomial to the (noisy) points computed in step 2.
4. Compute the gradient of that degree $d$ polynomial at the point we are interested in (for the moment $\vec{p}$ or $\frac{\vec{p}}{2}$).

Ideally, this should be optional, with an automatic exploration triggering if the bounds are not provided by the user.

Proposed solution

Add a function to deltakit_explorer.analysis.error_budget to automatically compute bounds that make sense for each noise parameter for $\Lambda$ computation.

These bounds should be chosen such that:

1. the value of $\Lambda$ sufficiently changes when changing a single parameter of the noise model within the allowed bounds for that parameter and,
2. the resulting logical error probability is not too small or too high.

Point 1 here ensures that we are not restricting too much the values of $\Lambda$ with the bounds and that our gradient estimation will not be mostly noise due to sampling.

Point 2 ensures that we can compute $\Lambda$ using Monte-Carlo simulations without having to use billions of shots to see one error.

This new function could for example start with a small pre-defined initial step around the value of interest (e.g., $[p - \epsilon, p + \epsilon]$) and then exponentially grow $\epsilon$ to see if the estimated $\Lambda$ values are different enough (note that this generates values for $\Lambda$ which might then be used as a basis for fitting a polynomial, but these points will likely be pretty bad for polynomial fitting because they are quite far away from the Chebyshev nodes which are ideal for fitting).

Potential alternatives

Any other alternative is welcomed.

Additional context

• https://deltakit-docs.riverlane.com/en/stable/guide/error_budgeting.html
• https://en.wikipedia.org/wiki/Chebyshev_nodes

[nrdcid]

Hi, I'm Nathan. I'm interested in taking on this issue as part of the Deltakit Community Fund and have just submit my application. Here's an initial implementation approach I'm proposing:

Proposition: find_error_budget_bounds

The goal would be a function find_error_budget_bounds(noise_model, p, circuit_or_experiment, ...) in deltakit_explorer.analysis.error_budget that returns auto-computed bounds: list[tuple[float, float]] suitable for passing directly to get_error_budget

Implementation

For each noise parameter dimension independently:

• Start from $[p - \epsilon, p + \epsilon]$ for a small initial $\epsilon$ (e.g. $p \cdot 0.1$).
• Run a quick pilot simulation (small num_shots, configurable) at $p - \epsilon$, $p$, and $p + \epsilon$ to estimate $\Delta$ with error bars $\sigma_\Delta$
• Check two stopping criteria:
  • Sensitivity: $|\Delta (p + \epsilon) - \Delta(p - \epsilon)| > k \cdot \sigma_\Delta$ (signal-to-noise ratio exceeds threshold $k$, default $k=3$). Ensures the gradient estimate won't be dominated by sampling noise.
  • LEP feasibility: LEP at both endpoints falls within a configurable range (e.g., $[1e-5, 0.4]$). Ensures neither extreme is a simulation dead zone.

If either criterion fails, double $\epsilon$ and repeat up to a configurable max iterations. Once criteria are met, optionally refit the endpoint to the nearest Chebyshev nodes of the expanded interval, so the returned bounds are directly useful for polynomial fitting in step 3 of the existing pipeline.

Deliverables

• find_error_budget_bounds function with docstring, typed signature, and unit tests
• Integration with get_error_budget so that passing bounds=None auto-invokes the new function
• A short notebook example in docs/ demonstrating auto-bound discovery on a rotated planar code
• The bounds parameter on get_error_budget remains required by default, but can be made optional with a deprecation path

---

Open questions I'd like your input on:

• Should find_error_budget_bounds support multi-parameter noise models (i.e., axis-aligned search per parameter), or is a single-parameter implementation the right scope for a first PR?
• Is there a preferred pilot-shot budget that won't feel wasteful to users (e.g., max_shots=10_000 per probe point)?
• Any existing internal prototype or related code that you recommend studying?

Happy to adjust the scope based on your feedback.

[vishal1601-2005]

Hi @nelimee, @nrdcid — i am vishal . this is a great scoping question at the end of your proposal, and I think it's worth pushing on directly.

The axis-aligned per-parameter search is the right starting point, but I don't think it's safe to assume the joint feasible region for a multi-parameter noise model is just the Cartesian product of the per-parameter boxes. If two noise parameters interact — e.g. one channel partially compensating another's effect on Λ — a corner of the combined box (all parameters simultaneously at their individual extremes) could easily fail both the SNR and LEP feasibility checks even though each parameter passes them independently. Treating the marginal boxes as the final answer risks silently returning bounds that are infeasible in the region that actually matters for fitting.

Proposed Phase 2, building directly on find_error_budget_bounds:

Take the per-parameter box from Nathan's function as an initial candidate region.
Run pilot simulations at the box's vertices (not just each axis independently) to check joint feasibility — same SNR (k·σ_Δ) and LEP-band criteria, but evaluated where all parameters move together.
If any vertex fails, contract the box (e.g. simple trust-region-style shrinkage toward the center point p) and re-check, rather than assuming the marginal box is correct by construction.
For the polynomial fit itself (step 3/4 of the existing pipeline), generalize the Chebyshev-node refit to a tensor-product or sparse (Smolyak) grid rather than independent per-axis nodes — a full tensor grid blows up combinatorially as the number of noise parameters grows, so a sparse grid keeps the pilot-shot budget from exploding while still giving usable multivariate fit points.

This treats "the bounds are valid" as a hard constraint that has to hold jointly, not just marginally — which is the same shape of problem I've worked through before: deriving a bound from first principles and enforcing it as a hard constraint on a multi-dimensional space, rather than a per-axis heuristic that happens to usually work.

Questions for you both:

Does scoping this as a follow-on to find_error_budget_bounds (rather than folding it into the same PR) make sense, or would you rather see joint feasibility handled in v1?
Is there an existing noise model in Deltakit's test suite with known parameter coupling I could use to validate that the marginal-box assumption actually breaks down, before building the contraction logic?

[nelimee]

Hi both,

Should find_error_budget_bounds support multi-parameter noise models (i.e., axis-aligned search per parameter), or is a single-parameter implementation the right scope for a first PR?

Single parameter is fine, only one parameter varies at a time in error-budgeting.

Is there a preferred pilot-shot budget that won't feel wasteful to users (e.g., max_shots=10_000 per probe point)?

That will depend a lot on the circuit distance. In practice, we want the runtime to be low. Something like a few seconds would be nice, maximum a dozen of seconds in my opinion. That should be parametrisable by users.

Any existing internal prototype or related code that you recommend studying?

Not really. I think your plan is good. One note about

Once criteria are met, optionally refit the endpoint to the nearest Chebyshev nodes of the expanded interval, so the returned bounds are directly useful for polynomial fitting in step 3 of the existing pipeline.

It is not obvious to me that the Chebyshev nodes are the optimal ones as here we really don't care about what is happening at the boundaries, we only care that the gradient at one specific point (that will be around the "middle" of the interval) has a low variance. Because of that, I am not sure the Chebyshev nodes are optimal. But they are probably a good start!

Does scoping this as a follow-on to find_error_budget_bounds (rather than folding it into the same PR) make sense, or would you rather see joint feasibility handled in v1?
Is there an existing noise model in Deltakit's test suite with known parameter coupling I could use to validate that the marginal-box assumption actually breaks down, before building the contraction logic?

Note that one of the fundamental assumption of error budgeting is that your noise channels are independent. If they are not, you are outside of the contract for noise budgeting, and you should expect potentially bad results from the method. So I don't think that using such a model is a good counter example. Whether you could have a different error-budgeting method that does not perform this assumption is another question. I think it should be quite easy to do conceptually, but the computational cost might explode. I have not explored that path though.

So to me, this analysis is out of the scope of this issue.

[zkasuran]

On the Chebyshev question: the criterion you describe is already merged here, so the repo can settle it.

DiscretisationStrategy.C_OPTIMAL from #266 (merged in f59bb4c) minimises exactly "the gradient at one specific point has a low variance". _c_optimal_objective forms cov = (X.T @ X)^-1 for a candidate design and feeds it to _get_variance_of_gradient_estimation_at_point, the function of yours that #266 restored from _gradient.py, so the objective is the variance of the fitted slope at c. That is c-optimality in the classical sense of minimising the variance of the best linear unbiased estimator of a linear functional of the coefficients (Elfving 1952, https://doi.org/10.1214/aoms/1177729442, next to the page _discretisation.py already cites). Chebyshev nodes minimise a different thing, the worst case interpolation error across the whole interval, which is the boundary behaviour you say does not matter here. The fit is also a degree 3 least squares over 10 points rather than an interpolation, so nodes are not quite the right frame for it. Where Chebyshev nodes do bear on this is the fit bias that _approximate_derivative_at_point_from_values already warns about under Runge's phenomenon, not the variance.

Reproducing #266's gain on the docs example (P = [5e-3], bounds (1e-3, 2e-2), 10 points, degree 3, equal variance per point) I get a slope variance at p/2 of 0.43x linear for c-optimal against 0.85x for logarithmic. So a bounds finder can hand back the interval alone: the strategy re-picks its own points inside it anyway, so endpoints refitted to Chebyshev nodes get discarded.

Two things from the same read that bear on the stopping test.

get_error_budget evaluates the gradient at p/2, so the point the discretisation receives is p/2 and _check_interval demands a < p/2 < b. [p - eps, p + eps] therefore raises the Expected a < c < b ValueError for any eps <= p/2. For P = [5e-3] and eps = 0.1p that is a = 0.0045 against c = 0.0025. Doubling from there, the first interval the pipeline accepts is eps = 0.8p. The docstrings of get_error_budget and inverse_lambda_gradient_at both say the bounds must contain p, which holds only for the latter.

The quantity to stop on is then the variance of the gradient under the discretisation that will actually run, rather than an endpoint signal to noise ratio. An endpoint test cannot see the degree, the point count or the w=1/stddevs weights np.polyfit is given. Those weights matter. _c_optimal_objective assumes equal variance at every point. With a sigma rising toward the lower bound (I tried sigma proportional to x^-1/2 and to x^-1 as illustrations, not measured) the shipped design carries 1.35x and 4.0x the slope variance of one chosen with the true weights, so pilot standard deviations are the input it is missing. The feasibility band need not be constants either: get_error_budget already documents 100 / max_shots as the lower target with 1 / 2 above.

All of the above is from _discretisation.py, _gradient.py and _budget.py at f59bb4c, plus running those two discretisation functions directly for the numbers.