Optimized Integration Contours

• Numerically optimized continuous integration contours are paths whose geometry is determined by task-specific numerical objectives instead of traditional analytic forms.
• They improve numerical stability in tasks such as high-order derivative estimation, eigenvalue isolation, and variance reduction in complex integral evaluations.
• Optimization methods span discrete grid searches, branch-and-bound techniques, and deep learning, yielding significant speedups and accuracy improvements.

Numerically optimized continuous integration contours are integration paths, contour families, or contour-targeted design objects whose geometry is chosen by an explicit numerical criterion rather than by fixed analytic convention. In the literature considered here, optimization serves several distinct but related purposes: minimizing the condition number of Cauchy integrals for high-order derivatives, isolating targeted eigenspectra in contour-integral generalized eigenvalue solvers, reducing integrand variance under complex contour deformation for sector-decomposed Feynman integrals, mitigating the sign problem by maximizing the mean phase factor on deformed manifolds, and accelerating level-set estimation by maximizing contour-based expected improvement under Gaussian-process surrogates (Bornemann et al., 2011, Chen et al., 2 Nov 2025, Jones et al., 16 Jan 2026, Pasztor et al., 9 Sep 2025, Bursa et al., 2021, Franey et al., 2010).

1. Conceptual scope and problem classes

The term “contour” is used in two closely related senses. In numerical contour integration proper, a contour is a closed path $\Gamma \subset \mathbb{C}$ or a deformed manifold in $\mathbb{C}^N$ on which an integral is evaluated. In surrogate-based sequential design, a contour is instead a level set $S(c)=\{x\in\chi:f(x)=c\}$, so optimization concerns the next sample location for estimating that set rather than a complex integration path. This broader usage is important because several optimization principles recur across both settings, especially the use of task-specific objective functions and global search procedures (Franey et al., 2010).

A numerical contour must typically satisfy geometric or analytic constraints in addition to optimizing a performance metric. For generalized eigenvalue problems, the contour should enclose exactly the target eigenvalues, remain sufficiently far from the spectrum to avoid ill-conditioning, and be as small as possible to reduce the number of quadrature nodes and the size of projected problems (Chen et al., 2 Nov 2025). For sector-decomposed Feynman integrals, a deformation $x_i \mapsto z_i=x_i+i\tau_i(x;s)$ is required to satisfy $\operatorname{Im}F(z;s)<0$ along the path, thereby avoiding pole crossing while reducing variance (Jones et al., 16 Jan 2026). For sign-problem mitigation, one deforms the original real contour to reduce phase oscillations of $e^{-S}$ while keeping the Jacobian tractable (Pasztor et al., 9 Sep 2025, Bursa et al., 2021). For high-order derivatives computed by Cauchy’s formula, the contour must retain winding number $1$ around the expansion point (Bornemann et al., 2011).

Setting	Contour or deformation	Optimization target
High-order derivatives	Closed contour $\gamma$	Minimize $\kappa(\gamma,n)$ or weighted length
GEP contour-integral solvers	Circle, ellipse, or polygon $\Gamma$	Enclose target eigenvalues with low cost
Sector-decomposed Feynman integrals	$\mathbb{C}^N$0	Reduce variance subject to $\mathbb{C}^N$1
Fermionic sign problem	$\mathbb{C}^N$2	Maximize average phase factor
Bose-gas contour deformation	Local field deformation $\mathbb{C}^N$3	Maximize mean phase factor
GP contour estimation	Level set $\mathbb{C}^N$4 and next sample $\mathbb{C}^N$5	Maximize contour-based EI

This range of uses suggests that “numerically optimized contours” are best understood as a family of task-adapted geometric constructions rather than a single algorithmic paradigm.

2. Objective functions and mathematical criteria

For high-order derivatives of a holomorphic function, the central quantity is the condition number of Cauchy’s integral,

$\mathbb{C}^N$6

Because the denominator is contour-independent by Cauchy’s theorem, minimizing $\mathbb{C}^N$7 is equivalent to minimizing the weighted length

$\mathbb{C}^N$8

The optimization problem is therefore geometric but directly tied to backward stability (Bornemann et al., 2011).

In contour-integral generalized eigenvalue solvers for $\mathbb{C}^N$9, the contour enters through the spectral projector

$S(c)=\{x\in\chi:f(x)=c\}$0

Moment-based methods such as CIRR compute

$S(c)=\{x\in\chi:f(x)=c\}$1

whereas FEAST-type subspace iteration applies $S(c)=\{x\in\chi:f(x)=c\}$2 repeatedly to a trial subspace. After quadrature discretization,

$S(c)=\{x\in\chi:f(x)=c\}$3

so contour design directly controls both the conditioning of shifted linear systems and the number of required solves (Chen et al., 2 Nov 2025).

For sector-decomposed Feynman integrals, the objective is variance reduction under a deformation constrained by analyticity and pole avoidance. The deformed integrand is

$S(c)=\{x\in\chi:f(x)=c\}$4

and the loss is

$S(c)=\{x\in\chi:f(x)=c\}$5

with

$S(c)=\{x\in\chi:f(x)=c\}$6

The paper shows that minimizing the estimated error of a Quasi-Monte Carlo sample is ill-defined, and accordingly uses global variance as the robust proxy (Jones et al., 16 Jan 2026).

In sign-problem applications, the optimization criterion is the average phase factor

$S(c)=\{x\in\chi:f(x)=c\}$7

Direct contour optimization introduces

$S(c)=\{x\in\chi:f(x)=c\}$8

and minimizes $S(c)=\{x\in\chi:f(x)=c\}$9 over a finite-parameter deformation $x_i \mapsto z_i=x_i+i\tau_i(x;s)$0. In the Bose-gas formulation the same aim is expressed as maximizing the mean phase factor for $x_i \mapsto z_i=x_i+i\tau_i(x;s)$1, or equivalently minimizing $x_i \mapsto z_i=x_i+i\tau_i(x;s)$2 or $x_i \mapsto z_i=x_i+i\tau_i(x;s)$3 (Pasztor et al., 9 Sep 2025, Bursa et al., 2021).

For contour estimation under Gaussian-process surrogates, the improvement function is designed to reward both $x_i \mapsto z_i=x_i+i\tau_i(x;s)$4 being close to the target level $x_i \mapsto z_i=x_i+i\tau_i(x;s)$5 and the predictive uncertainty $x_i \mapsto z_i=x_i+i\tau_i(x;s)$6 being large. Ranjan et al. proposed

$x_i \mapsto z_i=x_i+i\tau_i(x;s)$7

with $x_i \mapsto z_i=x_i+i\tau_i(x;s)$8. Franey, Ranjan, and Chipman then used the modified criterion

$x_i \mapsto z_i=x_i+i\tau_i(x;s)$9

where $\operatorname{Im}F(z;s)<0$0. This modification is specifically motivated by simpler, piecewise-monotonic partial derivatives that support branch-and-bound (Franey et al., 2010).

3. Optimization mechanisms and algorithmic architectures

One major line of work converts contour selection into a discrete global optimization problem. For high-order derivatives, Bornemann and Wechslberger embed a square grid $\operatorname{Im}F(z;s)<0$1 into the complex domain, assign each graph edge $\operatorname{Im}F(z;s)<0$2 the weight

$\operatorname{Im}F(z;s)<0$3

and seek a shortest enclosing walk. Provan’s algorithm reduces the search to a finite family of walks built from shortest paths and a single edge. A constrained version chooses $\operatorname{Im}F(z;s)<0$4, runs Dijkstra once from $\operatorname{Im}F(z;s)<0$5, then checks enclosing walks of the form $\operatorname{Im}F(z;s)<0$6. The method is accelerated by diagonal edges, adaptive refinement in a tubular neighborhood of the current walk, and Clenshaw–Curtis quadrature along the final straight segments (Bornemann et al., 2011).

A second line uses branch-and-bound over hyperrectangles. In contour-based sequential design, one assumes bounds

$\operatorname{Im}F(z;s)<0$7

and exploits the facts that $\operatorname{Im}F(z;s)<0$8, $\operatorname{Im}F(z;s)<0$9 for $e^{-S}$0, and $e^{-S}$1 for $e^{-S}$2. This yields conservative lower and upper bounds from the $e^{-S}$3 corner evaluations in $e^{-S}$4-space. The branch-and-bound routine splits the active hyperrectangle along its longest edge, updates global lower and upper bounds, and prunes any cell whose upper bound falls below the incumbent lower bound. The implementation recipe explicitly recommends a priority queue of active rectangles keyed by $e^{-S}$5 (Franey et al., 2010).

A third line is learning-based contour prediction. DeepContour first applies a Fourier Neural Operator, the Eigen-Neural-Operator, to predict an approximate spectrum $e^{-S}$6 from discretized physical parameters. It then applies a one-dimensional KDE sparsity function

$e^{-S}$7

to locate a splitting point $e^{-S}$8 between spectral clusters. Intervals with more than $e^{-S}$9 predicted eigenvalues are split recursively; smaller intervals are wrapped by contours, for example a tight circle with center $1$0 and radius $1$1. The resulting contours are then passed to CIRR or FEAST and solved in parallel (Chen et al., 2 Nov 2025).

Neural contour deformation for Feynman integrals follows a related but distinct strategy. The guided ansatz

$1$2

is generalized either to $1$3 or to $1$4. The guided-deformation network uses kinematic invariants $1$5 as input and has 4 hidden layers with 40 units each, tanh activations, and linear output; the free-deformation network takes both $1$6 and $1$7, uses 3 hidden layers with 40 units each and GELU activations, and relies on automatic differentiation for the Jacobian determinant. Training proceeds in two phases: first with large $1$8 gradually decayed at fixed learning rate, and then with learning-rate reduction on plateau once $1$9 is small (Jones et al., 16 Jan 2026).

In sign-problem studies, direct continuous contour optimization typically uses low-dimensional analytic ansätze. One-dimensional deformations take the single-valued form $\gamma$0, where $\gamma$1 may be a Fourier series $\gamma$2 or a piecewise-linear function. Gradients of the negative log-sign are evaluated as phase-quenched expectation values, and optimization is performed with Adam or stochastic gradient descent. In the higher-dimensional Bose-gas setting, the deformation is designed so that the Jacobian matrix becomes upper-block-triangular, yielding an $\gamma$3 local Jacobian factor per site and therefore $\gamma$4 total determinant cost (Pasztor et al., 9 Sep 2025, Bursa et al., 2021).

4. Representative formulations across domains

The oldest formulation in this group is the use of optimized contours for Cauchy integrals. For $\gamma$5, the numerical problem is not the existence of a contour but the choice of one that does not amplify quadrature error through a large condition number. The grid-path shortest-enclosing-walk construction is particularly effective when singular geometry makes circles unattractive, such as branch cuts or nearby singularities (Bornemann et al., 2011).

In generalized eigenvalue problems, contour optimization serves as a preprocessing stage for projection methods. The contour itself determines the subset of eigenpairs retained by the spectral projector, the conditioning of the shifted linear systems $\gamma$6, and the size of the projected subproblems. DeepContour’s combination of FNO spectral prediction with KDE partitioning replaces “scouting” methods such as short Arnoldi runs with a one-shot surrogate, then formalizes contour construction from predicted spectral gaps (Chen et al., 2 Nov 2025).

For Feynman integrals after sector decomposition, contour deformation is both a validity and an efficiency problem. The deformation must preserve the correct $\gamma$7 prescription by maintaining $\gamma$8, but it is also optimized to reduce integrand variance over an entire phase-space region rather than at a single kinematic point. The resulting learned map $\gamma$9 is designed to be frozen and reused, including in neural-network-based integrators that require a fixed analytic contour (Jones et al., 16 Jan 2026).

For fermionic toy models and lattice Bose gas, the central issue is the sign problem. Here contour optimization is not aimed primarily at quadrature error or spectral isolation but at increasing the mean phase factor under reweighting. The literature distinguishes three geometric objects: Lefschetz thimbles defined by downward flow into critical points, finite flow-time manifolds $\kappa(\gamma,n)$0, and directly optimized single-valued contours $\kappa(\gamma,n)$1. In the Bose-gas model, explicit first-order and second-order deformations are derived in a small parameter $\kappa(\gamma,n)$2, then generalized to a flexible local ansatz with free nonnegative parameters $\kappa(\gamma,n)$3 (Pasztor et al., 9 Sep 2025, Bursa et al., 2021).

A distinct but conceptually related formulation appears in deterministic computer experiments. The contour $\kappa(\gamma,n)$4 is a level set of interest, and the optimization problem is to choose new simulator evaluations that most efficiently improve the estimate of that set. The Gaussian-process surrogate provides $\kappa(\gamma,n)$5 and $\kappa(\gamma,n)$6, and the contour-based expected improvement criterion explicitly trades proximity to $\kappa(\gamma,n)$7 against predictive uncertainty, making global mean-squared-error reduction unnecessary when only the $\kappa(\gamma,n)$8-level set matters (Franey et al., 2010).

5. Numerical behavior and empirical comparisons

The numerical literature emphasizes that optimized contours can yield large gains relative to standard heuristics. In the Cauchy-derivative setting, for

$\kappa(\gamma,n)$9

with derivative order $\Gamma$0 at $\Gamma$1, the shortest enclosing walk on a $\Gamma$2 grid with diagonals gives $\Gamma$3, whereas the optimal circle gives $\Gamma$4. The paper also reports nearly identical performance between optimized walks and optimal circles for benign entire functions such as $\Gamma$5 and $\Gamma$6, indicating that the advantage is strongest in singular or branch-cut geometries (Bornemann et al., 2011).

In GP contour estimation, direct one-step comparisons show that branch-and-bound finds larger expected improvement than a genetic algorithm under equal evaluation budgets. For the Branin contour at $\Gamma$7 with $\Gamma$8, average best $\Gamma$9 is $\mathbb{C}^N$00 for BNB versus $\mathbb{C}^N$01 for GA; for the 2D Levy contour at $\mathbb{C}^N$02 with $\mathbb{C}^N$03, $\mathbb{C}^N$04 versus $\mathbb{C}^N$05; and for simultaneous max/min on Branin with $\mathbb{C}^N$06, $\mathbb{C}^N$07 versus $\mathbb{C}^N$08. Over longer sequential runs, BNB drives Branin contour-divergence to near $\mathbb{C}^N$09 in about $\mathbb{C}^N$10 new runs, compared with about $\mathbb{C}^N$11 for GA, while static sampling remains about $\mathbb{C}^N$12. For the 4D Levy contour at $\mathbb{C}^N$13, BNB reduces contour-error by $\mathbb{C}^N$14 in one new run whereas GA needs about $\mathbb{C}^N$15 runs (Franey et al., 2010).

For generalized eigenvalue problems, DeepContour reports end-to-end speedups up to $\mathbb{C}^N$16 and CI-solver speedups up to $\mathbb{C}^N$17. In the Kirchhoff–Love Plate case at tolerance $\mathbb{C}^N$18, the Arnoldi-based baseline yields $\mathbb{C}^N$19 end-to-end and $\mathbb{C}^N$20 CI-solver speedup when replaced by DeepContour; GD and JD baselines are also outperformed, at $\mathbb{C}^N$21 and $\mathbb{C}^N$22, respectively. DeepContour’s total contour area is reported as $\mathbb{C}^N$23–$\mathbb{C}^N$24 smaller than Arnoldi-scout contours, and ENO spectral prediction reaches NMSE $\mathbb{C}^N$25 versus Arnoldi-scout $\mathbb{C}^N$26. In ablations, removing ENO and using an MLP misses $\mathbb{C}^N$27 eigenvalues on average, removing KDE increases solve time by a factor of $\mathbb{C}^N$28, and replacing ENO by a Krylov–Schur scout plus KDE misses about $\mathbb{C}^N$29 eigenvalues (Chen et al., 2 Nov 2025).

For sector-decomposed Feynman integrals, the one-loop bubble study shows that the guided network $\mathbb{C}^N$30 reproduces the analytic optimum $\mathbb{C}^N$31 that minimizes $\mathbb{C}^N$32, while the free network $\mathbb{C}^N$33 lowers $\mathbb{C}^N$34 by an additional $\mathbb{C}^N$35–$\mathbb{C}^N$36. In the two-loop elliptic box with $\mathbb{C}^N$37 sectors and $\mathbb{C}^N$38, the guided network reduces integrand variance by factors $\mathbb{C}^N$39–$\mathbb{C}^N$40 over SecDec default, and the free network achieves up to a further $\mathbb{C}^N$41 reduction. QMC error estimates scale like $\mathbb{C}^N$42 with similar $\mathbb{C}^N$43 across methods, but the $\mathbb{C}^N$44 minimizing the QMC error depends on the lattice size $\mathbb{C}^N$45, which is one reason the paper characterizes direct QMC-error minimization as ill-posed (Jones et al., 16 Jan 2026).

In the sign-problem literature, the gains are often measured through $\mathbb{C}^N$46, the average phase factor. For the Hubbard-like toy model at $\mathbb{C}^N$47, the reported values are $\mathbb{C}^N$48, $\mathbb{C}^N$49, $\mathbb{C}^N$50, and $\mathbb{C}^N$51. For the Gross–Neveu-like model, the corresponding values are $\mathbb{C}^N$52, $\mathbb{C}^N$53, $\mathbb{C}^N$54, and $\mathbb{C}^N$55; for the Thirring-like model, $\mathbb{C}^N$56, $\mathbb{C}^N$57, $\mathbb{C}^N$58, and $\mathbb{C}^N$59; and for the Chern–Simons-like model, $\mathbb{C}^N$60, $\mathbb{C}^N$61, $\mathbb{C}^N$62, and $\mathbb{C}^N$63. The same source notes that because the variance of a reweighted estimator grows like $\mathbb{C}^N$64, a change from $\mathbb{C}^N$65 to $\mathbb{C}^N$66 in the Hubbard model corresponds to a $\mathbb{C}^N$67-fold gain in effective statistics (Pasztor et al., 9 Sep 2025).

The higher-order Bose-gas contour study reports exponential decay rates $\mathbb{C}^N$68 for $\mathbb{C}^N$69, $\mathbb{C}^N$70, with $\mathbb{C}^N$71 for the simple first-order contour, $\mathbb{C}^N$72 for the first-order ansatz, $\mathbb{C}^N$73 for the simple second-order contour, and $\mathbb{C}^N$74 for the second-order ansatz. At $\mathbb{C}^N$75 and $\mathbb{C}^N$76, the mean phase rises from about $\mathbb{C}^N$77 for the simple first-order contour to about $\mathbb{C}^N$78 for the ansatz second-order contour. The decay rate also decreases with spatial dimension $\mathbb{C}^N$79, consistent with the interpretation $\mathbb{C}^N$80 (Bursa et al., 2021).

6. Misconceptions, limitations, and open directions

A recurrent misconception is that analytically distinguished contours are automatically numerically optimal. The sign-problem studies explicitly reject this. Lefschetz thimbles make $\mathbb{C}^N$81 constant on each thimble, but they do not in general maximize $\mathbb{C}^N$82 because multiple thimbles can contribute with different phases and the local parametrization Jacobian can itself be complex. The reported toy-model results further show that the holomorphic-flow value at large flow time, corresponding to thimbles, can be inferior to a finite optimal flow time, and that directly optimized continuous contours can outperform both (Pasztor et al., 9 Sep 2025).

A second misconception is that variance reduction and integration error minimization are interchangeable. The neural contour-deformation work for Feynman integrals shows that optimizing a contour to reduce the estimated error of a Quasi-Monte Carlo sample is ill-defined, because no single deformation parameter $\mathbb{C}^N$83 is best for all lattice sizes and sample shifts. The paper nevertheless treats variance reduction as a robust proxy for most integrators. This suggests a general principle: the numerical objective should be chosen to reflect algorithmic invariants rather than noisy or discretization-specific error surrogates (Jones et al., 16 Jan 2026).

Several methods retain nontrivial complexity barriers. Branch-and-bound for contour-targeted expected improvement is globally convergent under its bounding conditions, but the worst-case number of subregions grows exponentially in dimension $\mathbb{C}^N$84, even though pruning is often effective in practice (Franey et al., 2010). The grid-path shortest-enclosing-walk approach has a full Provan-search cost of $\mathbb{C}^N$85, with the constrained version reduced to $\mathbb{C}^N$86, but automatic choice of the bounding box and mesh resolution remains an open practical issue (Bornemann et al., 2011).

The eigenvalue-solver literature also makes the scope of current guarantees explicit. DeepContour’s stability discussion states that extending the method to non-Hermitian or complex spectra requires $\mathbb{C}^N$87D-KDE in $\mathbb{C}^N$88 and possibly polygonal contours. Its reported error heuristic is tied to standard analyses in which quadrature error decays like $\mathbb{C}^N$89 when the contour encloses all eigenvalues with margin $\mathbb{C}^N$90; in practice, $\mathbb{C}^N$91 or $\mathbb{C}^N$92 is stated to suffice for residuals down to $\mathbb{C}^N$93 on circles, but this does not eliminate the need for accurate spectral prediction (Chen et al., 2 Nov 2025).

Open directions are explicit across the corpus. For high-order derivatives, a rigorous asymptotic analysis of $\mathbb{C}^N$94 comparable to the existing theory for circles is stated to be open, and extension to more general contour-deformation problems such as Riemann–Hilbert integrals is said to require new theory (Bornemann et al., 2011). For Feynman integrals, the proposed extensions include higher-dimensional deformations, spline-based deformations, and coupling with normalizing-flow sampling (Jones et al., 16 Jan 2026). For sign-problem mitigation, higher-dimensional generalization proceeds by parameterizing $\mathbb{C}^N$95 with low-degree polynomials, normalizing flows, or neural networks with analytic Jacobians, then optimizing $\mathbb{C}^N$96 by Monte Carlo gradients (Pasztor et al., 9 Sep 2025).

Taken together, these results indicate that numerically optimized continuous integration contours form a mature but heterogeneous methodology. The common structure is the replacement of fixed contour heuristics by an explicit optimization problem combining geometric admissibility, task-specific numerical objectives, and computationally tractable search or learning procedures. The differences lie in what is being optimized—condition number, quadrature stability, variance, mean phase, projected subspace quality, or contour-estimation accuracy—and in whether the contour is a graph walk, a parametric curve, a deformed manifold, or a level set target.