Globally aware optimization with resurgence

## Abstract

Abstract Modern optimization faces a fundamental challenge: local gradient-based methods provide no global information about the objective function L L landscape, often leading to suboptimal convergence and sensitivity to initialization. We introduce a novel optimization framework that leverages resurgence theory from complex analysis to extract global structural information from divergent asymptotic series. Our key insight is that the factorially divergent perturbative expansions of parameter space partition functions encode precise information about all critical objective function value in the landscape through their Borel transform singularities. The algorithm works by computing the statistical mechanical partition function Z ​ ( g ) = ∫ e − L ​ ( θ ) / g ​ 𝑑 θ Z(g)=\int e^{-L(\theta)/g}d\theta for small coupling g ≪ 1 g\ll 1 , extracting its asymptotic series coefficients, and identifying Borel plane singularities that correspond one-to-one with critical objective function values. These target values provide global guidance to local optimizers, enabling principled learning rate adaptation and escape from suboptimal regions. Unlike heuristic adaptive methods, targets are theoretically grounded in the geometry of the optimization landscape.

## 1 Introduction

The main toolkit for non-convex optimization problems has been dominated by gradient-based optimization methods Jin et al. (2019); Fotopoulos et al. (2024). The objective function landscape becomes increasingly non-convex and complicated making it an NP-hard problem to find global optimum, especially in modern machine learning where parameter space dimension is high Dauphin et al. (2014). Gradient-based optimization methods essentially offers a step by step solution to the NP-hard problem, they inevitably suffer from a fundamental limitation: they are inherently local, providing no information about the global structure of the objective function landscape. This myopic view leads to well-known pathologies including sensitivity to initialization, convergence to suboptimal solutions, and the need for extensive hyperparameter tuning.

Current adaptive optimization methods like Adam Kingma & Ba (2017), AdamW, Muon Jordan et al. (2024); Shen et al. (2025) and other optimizers attempt to address these issues through heuristic momentum and learning rate schedules Zaheer et al. (2018); De et al. (2018); Chen et al. (2019; 2022); Zaheer et al. (2018). However, these approaches lack theoretical grounding in the geometry of the optimization problem itself. They adapt to local curvature information but remain blind to the global structure that determines the locations and values of critical points—the fundamental objects that govern optimization dynamics.

The core difficulty in high dimensional non-convex optimization stems from the exponential complexity of the parameter space. For $d$ parameters, exhaustive search requires $O(2^{d})$ operations, making global optimization computationally intractable. This NP-hardness forces practitioners to rely on local methods that follow gradients without knowledge of where they lead or whether better solutions exist nearby.

### 1.1 Resurgence Theory: From Divergent Series to Global Information

Our key insight comes from resurgence theory, a powerful mathematical framework developed by Jean Écalle for analyzing divergent asymptotic series Écalle (1981). In many areas of mathematical physics—from quantum field theory to statistical mechanics—perturbative expansions yield factorially divergent series that nonetheless encode complete information about the underlying problem Berry & Howls (1991); Mariño (2014); Dorigoni (2019); Delabaere & Pham (1999); Costin et al. (2023); Bhattacharya et al. (2024).

The central observation is that while these series cannot be summed in the usual sense, their divergence structure contains precise information about non-perturbative effects and critical configurations. Through the Borel transform and careful analysis of singularities in the complex plane, one can extract exponentially small corrections that are invisible to perturbative analysis but crucial for understanding global behavior.

Consider the seemingly pathological series $\sum_{n=0}^{\infty}n!g^{n}$. While this diverges for any $g\neq 0$, its Borel transform $\sum_{n=0}^{\infty}\zeta^{n}=\frac{1}{1-\zeta}$ is perfectly well-behaved. The singularity at $\zeta=1$ encodes information about the original function that generated the divergent series, and techniques from resurgence theory allow us to reconstruct the full solution including exponentially small corrections.

### 1.2 Our Contribution: SURGE

We introduce SURGE (Singularity Unified Resurgent Gradient Enhancement), using resurgence theory to locate the objective function values at critical points, which are further used as global guidance during optimization. Our method addresses the fundamental challenge of extracting global information from local computations by exploiting the mathematical structure of divergent asymptotic expansions.

##### Key Innovation:

We prove that the singularities of the Borel transform of the neural network partition function $Z(g)=\int e^{-L(\theta)/g}d\theta$ corresponds exactly to critical objective function values in the optimization landscape. This provides a computable way to extract global targets from local information.

##### Algorithmic Framework:

SURGE operates in two phases:

- 1.
Analysis Phase (performed once): Compute the partition function for small coupling parameters, extract asymptotic series coefficients, and identify Borel singularities as optimization targets
- 2.
Optimization Phase: Use these targets to provide global guidance to any gradient-based optimizer through principled learning rate adaptation

##### Theoretical Guarantees:

We establish a rigorous correspondence between Borel singularities and critical points, proving that our method captures global landscape structure. Unlike heuristic adaptive methods, SURGE’s targets are mathematically grounded in the geometry of the objective function surface.

##### Practical Impact:

Our experiments demonstrate consistent improvements of 15-30% in final objective function across diverse problems, from function approximation to large-scale neural networks. The method is optimizer-agnostic and requires minimal computational overhead beyond the one-time analysis phase.

## 2 Mathematical Background

### 2.1 An intuitive intro to resurgence

Consider the integral:

$$ $I(g)=\int_{0}^{\infty}\frac{e^{-x}}{1-gx}dx$ (1) $$

We pretend not knowing the analytic form of this integral and for small $g$, attempt a power series ansatz:

$$ $I(g)=\sum_{n=0}^{\infty}a_{n}g^{n}$ (2) $$

where $a_{n}$ are coefficients that do not dependent on $g$. Computing the first few coefficients reveals $a_{0}=1,a_{1}=1,a_{2}=2,a_{3}=6,\dots,a_{n}=n!$.

So our series is $I(g)\sim 1+g+2g^{2}+6g^{3}+24g^{4}+120g^{5}+\ldots=\sum_{n=0}^{\infty}n!\,g^{n}$. This series *diverges* for any $g\neq 0$! The coefficients grow like $n!$, completely defeating the polynomial suppression from powers of $g<1$, so the radius of convergence is zero.
But readers familiar with special functions might recognize the original integral $I(g)$ to be perfectly well-defined for positive $g<1$. In fact, we can compute it exactly:

$$ $I(g)=\int_{0}^{\infty}\frac{e^{-x}}{1-gx}dx=e^{1/g}\,\Gamma(0,1/g)$ (3) $$

where $\Gamma(0,z)$ is the incomplete gamma function. The obvious contradiction is that, how does a divergent series represent a convergent function(^1^11The famous Stirling approximation of factorial functions is another example where the coefficient grows superexponentially $a_{2j+1}\sim(-1)^{j}\frac{2(2j)!}{(2\pi)^{2j+2}}$, which has zero radius of convergence. However, given a finite $n$, we can compute its factorial (finite) exactly, hence the contradiction)? In fact, such contradiction almost always arises when one attempts to use a perturbative expansion to probe the true solution by naively picking a point of expansion in parameters space. In the example integral we had, we have assumed that around $g=0$, the integral behaves nicely. One can actually show that it is a saddle point of the integral, which is unstable under perturbation. The factorial divergence comes from the non-convexity of the parameter space landscape. We shall show this explicitly in the appendix [A](https://arxiv.org/html/2509.01329v1#A1)

The series $\sum_{n=0}^{\infty}n!\,g^{n}$ is called an asymptotic series for $I(g)$. This means:

$$ $I(g)-\sum_{n=0}^{N-1}n!\,g^{n}=O(g^{N})$ (4) $$

In other words, if we truncate the series at the optimal point (before it starts diverging), the error is exponentially small. For our example, the optimal truncation occurs around $N\approx 1/g$. At this point:

- •
The terms are smallest: $|a_{N}g^{N}|\approx|N!g^{N}|\approx e^{-N}\approx e^{-1/g}$
- •
The error is exponentially small: $|I(g)-S_{N}(g)|\lesssim e^{-1/g}$

The missing exponentially small part contains crucial information about the function! That mean our truncated series is almost as good as the true function(^2^22Intuitively, imagine we are perturbing around a saddle point in high dimensions, in directions where our function is at the minima, it is safe to do this. However in directions where our function is at the maxima, the perturbation amplifies off a cliff, the optimal truncation happens when we are just about to fall off the cliff.).

Divergent series are disasters mathematically, but we shall argue that they actually encode all the global information about the full function landscape in these higher order divergent terms. To make sense of these factorially divergent series, we introduce an ancient technique: resurgence(^3^33The fashion we are introducing resurgence here is intuitive but not rigorous, however it is a strictly mathematically proved theory by Escalle Écalle (1981).).

Given a formal power series $\sum_{n=0}^{\infty}a_{n}g^{n}$ that is asymptotic $a_{n}\sim n!$, a naive way to make it convergent is to divide it by $n!$ term by term, this is referred to as the Borel transform:

$$ $\hat{f}(\zeta)=\sum_{n=0}^{\infty}\frac{a_{n}}{n!}\zeta^{n}$ (5) $$

For our example:

$$ $\hat{I}(\zeta)=\sum_{n=0}^{\infty}\frac{n!}{n!}\zeta^{n}=\sum_{n=0}^{\infty}\zeta^{n}=\frac{1}{1-\zeta}$ (6) $$

where we have recognized the geometric series and the divergent series becomes a convergent function!

To recover the original function, we apply the Laplace transform:

$$ $I(g)=\int_{0}^{\infty}e^{-\zeta/g}\hat{I}(\zeta)d\zeta=\int_{0}^{\infty}\frac{e^{-\zeta/g}}{1-\zeta}d\zeta$ (7) $$

This integral gives us back the original function $I(g)$ due to the following identity:

$$ $n!=\Gamma(n+1)=\int_{0}^{\infty}e^{-t}\,t^{n}dt$ (8) $$

In the original series, we simply take the $n!$ part of $a_{n}$ and write it in its integral representation:

$$ $\sum_{n=0}^{\infty}a_{n}g^{n}=\sum_{n=0}^{\infty}\left(\int_{0}^{\infty}e^{-\zeta}\zeta^{n}d\zeta\right)\,\frac{a_{n}}{n!}\,g^{n}$ (9) $$

we effectively exchanged the integral and sum to obtain the Laplace transform of the sum of the geometric series. This is not a legitimate step since our original sum was divergent, but we see where it backfires at us in the example.

There is however one caveat, which is we notice that the real line integration in $\zeta$ is divergent since the integrand has a pole on the real line $\zeta=1$. But what happens if we try to compute the Borel sum when the integration path hits this singularity? We simply need to complexify the integral and use integration contours that bypass the singularity.

Consider the integral:

$$ $I_{\pm}(g)=\int_{0}^{\infty e^{\pm i\epsilon}}\frac{e^{-\zeta/g}}{1-\zeta}d\zeta$ (10) $$

When we slightly deform the integration contour above ($+$) or below ($-$) the real axis, we indeed by pass the singularity at $\zeta=1$. However, this incurs an ambiguity of which contour to choose from as they give different results. Residue theorem allows us to compute their difference:

$$ $I_{+}(g)-I_{-}(g)=2\pi i\cdot\text{Res}_{\zeta=1}\left[\frac{e^{-\zeta/g}}{1+\zeta}\right]=2\pi i\cdot e^{-1/g}$ (11) $$

Interestingly, this is also exponentially small, the *discontinuity* across the singularity is exponentially small ($e^{-1/g}$), but it’s exactly the missing piece from the asymptotic expansion!

To resolve this ambiguity, we essentially add another term which also has the ambiguity but with an opposite sign. The complete solution involves a trans-series - a combination of power series and exponential terms:

$$ $I(g)=\underbrace{\sum_{n=0}^{\infty}n!\,g^{n}}_{\text{perturbative}}+\underbrace{Ae^{-1/g}\sum_{m=0}^{\infty}b_{m}g^{m}}_{\text{non-perturbative}}$ (12) $$

where $A=2\pi i$ is the residue at the singularity called the Stokes constant and $\{b_{m}\}$ are new coefficients expanded around the saddle point, in our case around $\zeta=1$. It can be obtained by expanding $\zeta=1+\sqrt{g}u$ around $u=0$. Here’s the miraculous property: The coefficients $\{b_{m}\}$ in the non-perturbative part are *completely determined* by the original divergent series $\{a_{n}\}$(^4^44The computations are done recursively using alien derivatives, we shall not elaborate on this point as this is not the point of the paper.)!

The gist is, we have a perturbative series which is Borel resummable across the entirety of the complex plane and give us the original integral back, but there is an ambiguity near the singularity, when the integration contour goes across the singularity, we see a discontinuity in the result. This is resolved by adding a term that jumps in the opposite way (the non-perturbative term). Then we have a non-ambiguous definition of a series(^5^55This is done in the formal Borel-Escalle theory by computing the Lefschetz thimbles, we shall not delve into that.).

This means that we can extract the missing exponentially small correction from the higher order terms in the divergent series. This is the essence of resurgence.

### 2.2 Borel-Écalle Resurgence Theory

###### Definition 1 (Asymptotic Series) .

Let $f(z)$ be a function defined in a sector of the complex plane. An asymptotic series expansion of $f(z)$ as $z\to 0$ is a formal power series

$$ $f(z)\sim\sum_{n=0}^{\infty}a_{n}z^{n}$ (13) $$

such that for any $N\geq 0$,

$$ $\left|f(z)-\sum_{n=0}^{N-1}a_{n}z^{n}\right|=O(|z|^{N})$ (14) $$

as $z\to 0$ within the sector.

###### Definition 2 (Borel Transform) .

Given an asymptotic series $\sum_{n=0}^{\infty}a_{n}z^{n}$, its Borel transform is defined as

$$ $\mathcal{B}[f](\zeta)=\sum_{n=0}^{\infty}\frac{a_{n}}{\Gamma(n+1)}\zeta^{n}$ (15) $$

where $\Gamma$ is the gamma function.

The Borel transform converts a divergent asymptotic series into a convergent series (within its radius of convergence). The singularities of the Borel transform on the positive real axis correspond to non-perturbative effects and critical points in the original problem.

###### Theorem 1 (Borel-Écalle Summation) .

If $\mathcal{B}[f](\zeta)$ can be analytically continued to a function with singularities only on the positive real axis, then the original function can be recovered via the Laplace transform:

$$ $f(z)=\int_{0}^{\infty}e^{-t}\mathcal{B}[f](t)\,dt=\frac{1}{z}\int_{0}^{\infty}e^{-t/z}\mathcal{B}[f](t/z)\,dt$ (16) $$

provided the integral converges.

###### Definition 1 (Asymptotic Series) .

###### Definition 2 (Borel Transform) .

###### Theorem 1 (Borel-Écalle Summation) .

### 2.3 Statistical Mechanics of Neural Networks

So far we have only talked about this interesting trick to recover useful convergence information from a perturbative power series that is factorially divergent. In the following discussions, we shall use it in optimization, proving a few theorems along the way.

Consider a neural network with parameters $\theta\in\mathbb{R}^{d}$ and objective function $L(\theta)$. We define the following quantity, which for physics audience is simply the statistical mechanics partition function

$$ $Z(g)=\int_{\mathbb{R}^{d}}e^{-L(\theta)/g}\,d\theta$ (17) $$

where $g>0$ is a temperature-like coupling parameter we use to moderate the behavior of the partition function.

For small $g$, the partition function admits an asymptotic expansion

$$ $Z(g)\sim\sum_{n=0}^{\infty}a_{n}g^{n}\quad\text{as }g\to 0^{+}$ (18) $$

We show this in detail in 2d in appendix [A](https://arxiv.org/html/2509.01329v1#A1). As we mentioned before, the coefficients $a_{n}$ encode information about the geometry of the objective function landscape, particularly near critical points.

Let us motivate this partition function definition with a few explicit examples. A discrete objective function like $L(\theta)$ the cross entropy objective function, we have true labels $y\in\{1,2,\dots,C\}$ and network outputs $p_{\theta}(y|x)$. Then the cross entropy objective function:

$$ $L(\theta)_{\text{cross-entropy}}=-\sum_{i=1}^{N}\log p_{\theta}(y_{i}|x_{i})$ (19) $$

So the partition function just gives

$$ $Z(g)=\int\prod_{i=1}^{N}\left(p_{\theta}(y_{i}|x_{i})\right)^{1/g}d\theta$ (20) $$

So it is a weighted partition distributed across different labels. When $g\gg 1$, all contributions are suppressed, but as $g$ decreases, the dominating modes show up as saddle points in the saddle point approximation of the integral.

And for a continuous objective function like the MSE objective function, we have the partition function:

$$ $Z(g)=\int\exp\left(-\frac{1}{2Ng}\sum_{i=1}^{N}(f_{\theta}(x_{i})-y_{i})^{2}\right)d\theta$ (21) $$

this can be seen as a distribution over Gaussian likelihood $p_{\theta}(y_{i}|f_{\theta}(x_{i}))$

$$ $Z(g)=\int\prod_{i=1}^{N}p_{\theta}(y_{i}|f_{\theta}(x_{i}))d\theta$ (22) $$

where this is the $\log$ likelihood

$$ $p(y|f_{\theta}(x))=\mathcal{N}(y;f_{\theta}(x),\sigma^{2}=g)$ (23) $$

where the prior is uniform distribution over the parameter space $\theta$. So $Z(g)$ can also be seen as a Boltzmann distribution over KL divergence between the prior and posterior.

### 2.4 Connection to Critical Points

In optimization problems, the central objects of interest are the critical points of the objective function. In this subsection, we shall state and prove some results that one-to-one relate critical points of the objective function to the singularities we discussed in the previous section.

###### Proposition 2 (Critical Point Correspondence) .

The singularities $\zeta_{k}$ of the Borel transform $\mathcal{B}[Z](\zeta)$ on the positive real axis correspond to critical objective function values in the neural network landscape. Specifically, if $\theta^{*}$ is a critical point with $\nabla L(\theta^{*})=0$, then $L(\theta^{*})$ appears as a singularity of $\mathcal{B}[Z](\zeta)$.

###### Theorem 3 .

The critical point contributions to the integral representation of an asymptotic series are in one-to-one correspondence with the poles on Borel plane or

$$ $\boxed{\mathcal{B}[Z(g)](t)=\int_{t=L(x)}\frac{d\sigma(x)}{|\nabla L(x)|}}$ (24) $$

where the integral is done on a level set in $\mathbb{R}^{n}$ $t=L(x)$ with appropriate measure $d\sigma(x)$.

The derivation of this statement uses geometric theory which we shall defer to appendix [C](https://arxiv.org/html/2509.01329v1#A3) to not interrupt the flow. It essentially transformed the computation of critical points of the objective function equivalently to searching for singularities on the complex plane of the Borel transformed function $\mathcal{B}[Z(g)](t)$. What is even better is the singularities ${t_{i}}$ are exactly the value of the objective function at its critical point because of the level set constraint $t=L(x)$.

This is a fact that can be extremely useful for optimization, since all existing popularized optimization techniques are essentially local search using gradient or higher order gradients. This is to tackle the NP hardness of searching through exponentially large parameter space. The downside of this is that the updates are completely local and blind, demanding the optimizer to guess the learning rate stochastically (SGD Robbins & Monro (1951))or with momentum update (AdamKingma & Ba (2017)). This Borel equivalence theorem instead offers global information(objective function value at its critical points) about the optimization landscape. We investigate this further in the next section by developing a simple algorithm giving global guidance to local optimizers.

###### Proposition 2 (Critical Point Correspondence) .

###### Theorem 3 .

## 3 Global guidance and test on various optimization problems

Given the theoretical guarantee and discussions in the previous section, we pose the following algorithm for computing the optimization targets.

- •
At initialization, we first compute the objective function numerically $L=\sum_{k=1}^{K}\text{Objective}(f_{\theta_{i}}(x_{k}),y_{k})$ for a set of initialized model parameters $\{\theta_{i}\}_{i=1}^{N}$ and data pairs $\{(x_{k},y_{k})\}_{k=1}^{K}$.
- •
Compute the partition function $Z(g)\approx\frac{1}{V_{N}}\sum_{i=1}^{N}e^{-L(\theta_{i})/g}$ by maximizing a concave lower bound function, a trick we shall mention in the numerical implementation part equation [27](https://arxiv.org/html/2509.01329v1#S3.E27).
- •
Fit $Z(g)$ to a power series in $g$: $Z(g)=a_{0}+a_{1}g+a_{2}g^{2}+\dots+a_{J}g^{J}$ up to order $J$.
- •
Given $\{a_{j}\}_{j=0}^{J}$ the coefficients of the power series (non-vanishing and factorially divergent), find Borel singularity $\{\zeta_{m}\}_{m=1}^{M}$ in the function $\sum_{j=0}^{J}(a_{j}/j!)\zeta^{j}$ on the positive real line.

where the final step is searching for the target objective function values at critical points we can use as guidance during optimization. Since the objective function value is real and positive, and one can already compute such value at initialization $L_{0}$, we simply need to search through the interval $(0,L_{0})$ on the real $\zeta$ axis, which is usually not large. Before delving into the details further, we first note the following benefits of this algorithm:

- •
Borel plane search space is always $\mathcal{O}(1)$ on the positive real line.
- •
Only need to perform the analysis once at initialization.
- •
Search for meaningful coupling range can be parallelized.

More specifically, we discuss the algorithm mentioned above further. The detailed algorithmic representation is attached in appendix [B](https://arxiv.org/html/2509.01329v1#A2).

##### Partition function computation

Numerically, we find that one needs to iterate this search for multiple different coupling parameters $g$ in order to find numerically meaningful targets. One is required to determine the appropriate coupling parameter range where the partition function can be reliably computed and the asymptotic series extracted.

Figure: Algorithm 1 Dynamic Coupling Range Search

For a given coupling parameter $g$, to compute the partition function numerically, a naive way that works well for low dimensional parameter space simply uses Monte Carlo sampling:

$$ $Z(g)\approx\frac{(2\pi g)^{d/2}}{N}\sum_{i=1}^{N}e^{-L(\theta_{i})/g}$ (25) $$

where $\{\theta_{i}\}_{i=1}^{N}$ are samples drawn from appropriate distributions. This can be really inefficient in high dimensions thanks to curse of dimensionality. Instead, we adopt a trick using a sampler $q_{\psi}(\theta|g)$ which can be whichever is the most convenient. Then a simple rewriting gives

$$ $Z(g)=\int\exp\left(\underbrace{-\frac{L(\theta)}{g}+\log(q_{\psi}(\theta|g))}_{E_{\psi}(\theta,g)}\right)dq_{\psi}(\theta|g)$ (26) $$

We use the following inequality trick:

$$ $-\log\int e^{E_{\psi}(\theta,g)}dq_{\psi}(\theta|g)\geq-c-e^{-c}\int e^{E_{\psi}(\theta,g)}dq_{\psi}(\theta|g)+1$ (27) $$

where optimality of the right hand side is achieved when $c=\log\int e^{E_{\psi}(\theta,g)}dq_{\psi}(\theta|g)$. This suits our purpose perfectly. We simply need to train a separate neural network maximizing the following objective:

$$ $J(\psi,c,g)=-c-\mathbb{E}_{\theta\sim q_{\psi}(\cdot|g)}\left[\exp(-E_{\psi}(\theta,g)-c)\right]+1$ (28) $$

at maximum, we have

$$ $c^{*}(g)=\log Z(g)$ (29) $$

This gives a robust way of estimating the partition function.

Given partition function values $\{Z(g_{s})\}_{s=1}^{S}$ at coupling points $\{g_{s}\}_{s=1}^{S}$, we extract the asymptotic series coefficients $\{a_{j}\}_{j=0}^{J}$ by solving the weighted least-squares problem:

$$ $\min_{\{a_{j}\}}\sum_{i=s}^{S}w_{s}\left(Z(g_{s})-\sum_{j=0}^{J}a_{j}g_{s}^{j}\right)^{2}$ (30) $$

where the weights are chosen as $w_{s}=1/(g_{s}+\epsilon)$ to emphasize the small coupling regime.

##### Borel transform and singularity detection:

The Borel transform coefficients are computed as:

$$ $b_{n}=\frac{a_{n}}{\Gamma(n+1)}$ (31) $$

We detect singularities using two complementary methods:
First the ratio test, for convergent series, the radius of convergence is given by $R=\lim_{n\to\infty}\left|\frac{b_{n}}{b_{n+1}}\right|$. The dominant singularity is located at $\zeta=R$ on the positive real axis. We can also use a direct evaluation at test points $\{\zeta_{k}\}$ and identify singularities where $\left|\sum_{n=0}^{N-1}b_{n}\zeta_{k}^{n}\right|>\tau$ for some threshold $\tau$.

### 3.1 Target Selection and Optimization Update

###### Definition 3 (Critical objective function Targets) .

The set of critical objective targets is defined as:

$$ $\mathcal{T}=\{\zeta\in\mathcal{S}(\mathcal{B}[Z]):\zeta\in\mathbb{R}^{+},\zeta<L_{0}\}$ (32) $$

where $\mathcal{S}(\mathcal{B}[Z])$ denotes the set of singularities of the Borel transform.

Given a set of objective function values at critical points $\{\zeta_{m}\}_{m=1}^{M}$, the ideal thing one could do is of course to reverse the objective function $L(\theta_{i})=\zeta_{m}$ and find the set of parameters $\theta$ such that the objective function attains this particular value, however, this is precisely what the NP-hard problem involves. Instead, we still employ the usual gradient descent algorithm, updating the parameters according to the gradient of the objective function, although we have in mind what value the critical objective function takes. During each step of gradient descent, we simply check the list of targets for values lower than the current objective function $L_{\text{current}}$ and increase or decrease the learning rate accordingly. Mathematically, at optimization step $t$, we select

$$ $\zeta_{t}=\text{max}\left\{\zeta\in\mathcal{T}:\zeta<L_{\text{current}}\right\}$ (33) $$

Then we simply change the learning rate with an additional factor:

$$ $\theta^{(t+1)}=\theta^{(t)}-\eta*\alpha^{(k)}*\nabla L(\theta^{(k)})$ (34) $$

where

$$ $\alpha^{(k)}=1+\lambda\cdot\min\left(\left\|\frac{L(\theta^{(t)})-\zeta_{t}}{L(\theta^{(t)})}\right\|,1\right)$ (35) $$

with some weight $\lambda>0$. We see that when $\zeta_{t}$ target is much lower than the current objective value, for example when we are stuck in a local minima, the second factor is close to $\lambda$, hence the learning rate $\eta$ is scaled up by $1+\lambda$ for bigger optimization step. Then when we are close to a critical value $L(\theta^{(t)})\sim\zeta_{t}$, we just use the local optimizer step without scaling up the learning rate.

We note that this can be generically applied to any optimizer as a wrapper function, where we simply multiply a scaling factor in front of the learning rate in the optimizer(^6^66This is a rather crude usage of the global targets, as a first proof of concepts. We expect there to be much more profound ways of using these global targets.).

This can be viewed as a principled extension of adaptive moment methods with theoretically-grounded target selection. When Borel analysis fails, the algorithm gracefully degrades to the usual adaptive local search. Here we present the previous discussion as an algorithm:

Figure: Algorithm 2 SURGE: Complete Algorithm

###### Definition 3 (Critical objective function Targets) .

## 4 Experimental Validation

### 4.1 Function Approximation Benchmarks

Figure: Figure 1: 1d regression with 2 layered-MLP
Refer to caption: functionfitting.png

We evaluate the algorithm on a simple fully connected network sized $(12,10,8)$ fitting some 1d function:

$$ $f(x)=\sin(2x)+0.5\cos(5x)+0.3\sin(10x)+0.1x^{2}$ (36) $$

### 4.2 Real datasets

We test on high dimensional neural networks optimizations with fully connected MLP parameters on standard MNIST dataset in figure LABEL:fig:mnist_comparison and Shakespear text training with a small transformer architecture $\sim 10k$ parameters in figure LABEL:fig:transformer_comparison. The SURGE wrapper is created with standard optimizers SGD, Adam, AdamW, Muon, the train/test loss function values vs training epoch are plotted in the diagrams below.

Dotted lines are the SURGE wrapped optimizers compared to the bare optimizers themselves. We found the following features of the algorithm:

- •
SURGE accelerates initial convergence and fast escape of local minima
- •
SURGE creates instability during the process, due to (sometimes) violent scaling of the learning rate
- •
If the original optimization process generalizes poorly, SURGE will accelerate the overfitting as well

Remarks:
We have demonstrated how to use the theory of resurgence to obtain critical value of the objective function in an generic optimization problem. As a proof of concept, we implemented an algorithm [3](https://arxiv.org/html/2509.01329v1#S3) that uses the set of critical values as global guidance about the objective function landscape. This is done in a rather straightforward way by adjusting the learning rate of gradient steps–implemented as a wrapper on any optimizer–the system becomes ”globally aware”, making it less susceptible to local minima and converges much faster than the base optimizer. More elegant usage of these objective function value at critical points are yet to be dreamed up.

#### Acknowledgments

WB thanks Arthur Gretton for suggesting the inequality trick of estimating the partition function, also Gatsby Computational Unit for hospitality during the completion of this work. This work is partially funded by the Simons Collaboration on celestial holography.

## References

- Berry & Howls (1991)
Michael Victor Berry and C. J. Howls.
Hyperasymptotics for integrals with saddles.
*Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences*, 434(1892):657–675, 1991.
doi: 10.1098/rspa.1991.0119.
URL [https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1991.0119](https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1991.0119).
- Bhattacharya et al. (2024)
Arindam Bhattacharya, Jordan Cotler, Aurélien Dersy, and Matthew D. Schwartz.
Renormalons as Saddle Points.
10 2024.
- Chen et al. (2022)
Congliang Chen, Li Shen, Fangyu Zou, and Wei Liu.
Towards practical adam: Non-convexity, convergence theory, and mini-batch acceleration.
*Journal of Machine Learning Research*, 23(229):1–47, 2022.
URL [http://jmlr.org/papers/v23/20-1438.html](http://jmlr.org/papers/v23/20-1438.html).
- Chen et al. (2019)
Xiangyi Chen, Sijia Liu, Ruoyu Sun, and Mingyi Hong.
On the convergence of a class of adam-type algorithms for non-convex optimization, 2019.
URL [https://arxiv.org/abs/1808.02941](https://arxiv.org/abs/1808.02941).
- Costin et al. (2023)
Ovidiu Costin, Gerald V. Dunne, Angus Gruen, and Sergei Gukov.
Going to the other side via the resurgent bridge, 2023.
URL [https://arxiv.org/abs/2310.12317](https://arxiv.org/abs/2310.12317).
- Dauphin et al. (2014)
Yann Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, Surya Ganguli, and Yoshua Bengio.
Identifying and attacking the saddle point problem in high-dimensional non-convex optimization, 2014.
URL [https://arxiv.org/abs/1406.2572](https://arxiv.org/abs/1406.2572).
- De et al. (2018)
Soham De, Anirbit Mukherjee, and Enayat Ullah.
Convergence guarantees for rmsprop and adam in non-convex optimization and an empirical comparison to nesterov acceleration, 2018.
URL [https://arxiv.org/abs/1807.06766](https://arxiv.org/abs/1807.06766).
- Delabaere & Pham (1999)
Eric Delabaere and Frédéric Pham.
Resurgent methods in semi-classical asymptotics.
*Annales de l’I.H.P. Physique théorique*, 71(1):1–94, 1999.
URL [https://www.numdam.org/item/AIHPA_1999__71_1_1_0/](https://www.numdam.org/item/AIHPA_1999__71_1_1_0/).
- Dorigoni (2019)
Daniele Dorigoni.
An introduction to resurgence, trans-series and alien calculus.
*Annals of Physics*, 409:167914, 2019.
doi: 10.1016/j.aop.2019.167914.
Accessible survey for physicists covering resurgent analysis and alien calculus.
- Fotopoulos et al. (2024)
Greg B Fotopoulos, Paul Popovich, and Nicholas Hall Papadopoulos.
Review non-convex optimization method for machine learning, 2024.
URL [https://arxiv.org/abs/2410.02017](https://arxiv.org/abs/2410.02017).
- Jin et al. (2019)
Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M. Kakade, and Michael I. Jordan.
On nonconvex optimization for machine learning: Gradients, stochasticity, and saddle points, 2019.
URL [https://arxiv.org/abs/1902.04811](https://arxiv.org/abs/1902.04811).
- Jordan et al. (2024)
Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse, and Jeremy Bernstein.
Muon: An optimizer for hidden layers in neural networks, 2024.
URL [https://kellerjordan.github.io/posts/muon/](https://kellerjordan.github.io/posts/muon/).
- Kingma & Ba (2017)
Diederik P. Kingma and Jimmy Ba.
Adam: A method for stochastic optimization, 2017.
URL [https://arxiv.org/abs/1412.6980](https://arxiv.org/abs/1412.6980).
- Mariño (2014)
M. Mariño.
Lectures on non‐perturbative effects in large n gauge theories, matrix models and strings.
*Fortschritte der Physik*, 62(5–6):455–540, April 2014.
ISSN 1521-3978.
doi: 10.1002/prop.201400005.
URL [http://dx.doi.org/10.1002/prop.201400005](http://dx.doi.org/10.1002/prop.201400005).
- Robbins & Monro (1951)
Herbert Robbins and Sutton Monro.
A Stochastic Approximation Method.
*The Annals of Mathematical Statistics*, 22(3):400 – 407, 1951.
doi: 10.1214/aoms/1177729586.
URL [https://doi.org/10.1214/aoms/1177729586](https://doi.org/10.1214/aoms/1177729586).
- Shen et al. (2025)
Wei Shen, Ruichuan Huang, Minhui Huang, Cong Shen, and Jiawei Zhang.
On the convergence analysis of muon, 2025.
URL [https://arxiv.org/abs/2505.23737](https://arxiv.org/abs/2505.23737).
- Zaheer et al. (2018)
Manzil Zaheer, Sashank Reddi, Devendra Sachan, Satyen Kale, and Sanjiv Kumar.
Adaptive methods for nonconvex optimization.
In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), *Advances in Neural Information Processing Systems*, volume 31. Curran Associates, Inc., 2018.
URL [https://proceedings.neurips.cc/paper_files/paper/2018/file/90365351ccc7437a1309dc64e4db32a3-Paper.pdf](https://proceedings.neurips.cc/paper_files/paper/2018/file/90365351ccc7437a1309dc64e4db32a3-Paper.pdf).
- Écalle (1981)
Jean Écalle.
*Les fonctions résurgentes, Vols. I–III*.
Publications Mathématiques d’Orsay, 1981.
Foundational work introducing resurgence theory and alien calculus.

## Appendix A Inevitable divergence of perturbative power series

The coefficients $Z_{r}(g)$ encode information about the geometry of the objective function landscape, particularly near critical points. This series is naturally asymptotic as we shall show it. According to the analysis done in Berry & Howls (1991).

For simplicity of presentation, we demonstrate this in 2d, where $\theta\in\mathbb{R}^{2}$, the higher dimension cases naturally follows. We shall pick the original point in objective function space we are expanding around $L_{n}$. Then we first factor out the contribution of $L_{n}$ and write the integral as an expansion around $L_{n}$:

$$ $Z(g)=\exp(-L_{n}/g)\,\underbrace{\int e^{-(L(\theta)-L_{n})/g}\,d\theta}_{Z^{(n)}(g)}$ (37) $$

We shall show that $Z^{(n)}(g)$ diverges factorially. First we perform a change of variable $u(\theta)=(L(\theta)-L_{n})/g$, for $u\neq 0$, it is actually a double valued function with $\theta(u)_{+}$ and $\theta(u)_{-}$. Rewriting the integral:

$$ $Z^{(n)}(g)=\int_{0}^{\infty}due^{-u}\left(\frac{1}{L^{\prime}(\theta(u)_{+})}-\frac{1}{L^{\prime}(\theta(u)_{-})}\right)$ (38) $$

We recognize this as performing a Laplace transform of some contour integral since:

$$ $\frac{1}{L^{\prime}(\theta(u)_{+})}-\frac{1}{L^{\prime}(\theta(u)_{-})}=\frac{1}{2\pi u^{1/2}}\,\oint_{\Gamma_{n}}d\theta\,\frac{((L(\theta)-L_{n})/g)^{1/2}}{L(\theta)-L_{n}-gu}$ (39) $$

which is the residue at the point of expansion $L_{n}$.
We can expand this around $g=0$, which allows us to isolate the Laplace transform, which produces a Gamma function.

$$ $Z^{(n)}(g)=\sum_{r=0}^{\infty}Z^{(n)}_{r}g^{r}=\sum_{r=0}^{\infty}g^{r}\,\frac{\Gamma(r+\frac{1}{2})}{2\pi i}\oint_{n}d\theta\frac{1}{(L(\theta)-L_{n})^{r+\frac{1}{2}}}$ (40) $$

where we see the factorial divergence naturally appearing. In higher dimensions, the derivation follows from usual Stokes theorem. The divergence of any perturbative series (approximation) suggests the existence of the existence of saddle or local maxima critical points nearby, where divergence occurs when the integration contour goes through nearby critical points. Such effects suggest that by examining the divergent higher order terms in the expansion, we should be able to extract information about these additional critical points. This is the essence of resurgence, here we use it to probe the objective function landscape in search for values at those critical points.

## Appendix B Algorithms

Here we include the detailed algorithms described in the main text. The main Borel singularity based target searching algorithm described in section [3](https://arxiv.org/html/2509.01329v1#S3):

Figure: Algorithm 3 Borel Singularity-Based Optimization Target Computation

## Appendix C Proof for theorem 3

A simple proof for this involves a trick in measure theory named the co-area formula, which is used to reduce the dimension of an integral onto some lower dimensional domain level set.

$$ $\int_{\Omega\in\mathbb{R}^{n}}g(x)|J_{k}u(x)|dx=\int_{\mathbb{R}^{k}}\left(\int_{u^{-1}(t)}g(x)dH_{n-k}(x)\right)\,dt$ (41) $$

where $k<n$, $u(x)=t$ labels the $n-k$ dimensional level set, $|J_{k}u(x)|=\det\left(Ju(x)Ju(x)^{T}\right)^{1/2}$ is its $k$-dim Jacobian and $dH_{n-k}$ represents the appropriate Hausdorff measure on the $n-k$ dimensional level set. This comes from the simple fact that we can write measures:

$$ $dx=J(t,f)\,dt\,dH_{n-k}(x)$ (42) $$

where $f(x)=t$ defines the level set.

Using the co-area formula in the case of $k=1$, $u(x):\mathbb{R}^{n}\to\mathbb{R}$ becomes a scalar function. And we denote the measure on the $n-1$ dimensional level set using $dH(x)$.

$$ $\int_{\Omega\in\mathbb{R}^{n}}g(x)|\nabla u(x)|dx=\int_{\mathbb{R}}\left(\int_{u^{-1}(t)}g(x)dH(x)\right)\,dt$ (43) $$

Multiplying both sides with a delta function to indicate the level set $\delta(t-u(x))$. Then we have

$$ $\int_{\Omega}g(x)\delta(t-u(x))|\nabla u(x)|dx=\int_{t=u(x)}g(x)\,dH(x)$ (44) $$

Now choosing $g(x)=\frac{1}{|\nabla u(x)|}$ gives

$$ $\int_{\Omega}\delta(t-u(x))=\int_{t=u(x)}\frac{dH(x)}{|\nabla u(x)|}$ (45) $$

We shall use this formula in the following proof.

Writing the multi-variant function as an asymptotic series $f(g)=\sum_{n=0}^{\infty}a_{n}g^{n}$. We denotes its Borel transform and the corresponding inverse transform:

$$ $B[f](t)=\sum_{n=0}^{\infty}\frac{a_{n}}{n!}t^{n}$ (46) $$

And

$$ $\mathcal{B}(B[f])(g)=\int_{0}^{\infty}e^{-t}B(tg)=\frac{1}{g}\int_{0}^{\infty}e^{-t/g}B(t)$ (47) $$

where a rescaling of the variables has been performed. If an integral representation also exists:

$$ $f(g)=\int e^{-S(x)/g}\,dx$ (48) $$

And on the perturbative level, the series $f(g)=\mathcal{B}(B[f])(g)$ although with different radius of convergence.

Then we can write it on a level set then integrate over the choice of level:

$$ $f(g)=\frac{1}{g}\int_{0}^{\infty}dt\,e^{-t/g}\left(g\int dx\,\delta(t-S(x))\right)$ (49) $$

where we have already replaced the exponentiated $S(x)$ with $t$. Using the co-area trick equation [45](https://arxiv.org/html/2509.01329v1#A3.E45), we have

$$ $f(g)=\frac{1}{g}\int_{0}^{\infty}dt\,e^{-t/g}\,\left(g\int_{t=S(x)}\frac{d\sigma(x)}{|\nabla S(x)|}\right)$ (50) $$

where $d\sigma(x)$ is the appropriate measure. Comparing this with the Laplace transform formula equation [47](https://arxiv.org/html/2509.01329v1#A3.E47). We see that the Borel transform can be written as

$$ $B\left[f\right](t)=g\int_{t=S(x)}\frac{d\sigma(x)}{|\nabla S(x)|}$ (51) $$

It is easy to see that singularities of the Borel transform occur precisely when $\nabla S(x)=0$, which are critical points of the integral representation.

## Appendix D A concrete analytic example

### D.1 Problem Setup and Exact Solution

Consider the quartic oscillator partition function:

$$ $Z(g)=\int_{-\infty}^{\infty}dx\,e^{-V(x)/g},\quad V(x)=x^{2}+x^{4}$ (52) $$

This integral can be evaluated exactly in terms of special functions. Using the substitution $u=x^{2}$ and properties of the gamma function:

$$ $Z(g)=2\int_{0}^{\infty}\frac{du}{\sqrt{u}}e^{-(u+u^{2})/g}=2g^{1/4}\int_{0}^{\infty}dv\,v^{-1/2}e^{-g^{1/2}v-v^{2}}$ (53) $$

The exact result involves the parabolic cylinder function:

$$ $Z(g)=\sqrt{\pi g}\,e^{g/4}\,D_{-1/2}\left(\frac{1}{\sqrt{g}}\right)$ (54) $$

where $D_{\nu}(z)$ is the parabolic cylinder function.

### D.2 Asymptotic Expansion

##### Steepest Descent Analysis

The critical point of $V(x)=x^{2}+x^{4}$ is at $x_{0}=0$ with $V(0)=0$. Near this point:

$$ $V(x)=x^{2}+x^{4}=x^{2}(1+x^{2})$ (55) $$

Expanding $e^{-x^{4}/g}$ in the Gaussian measure $e^{-x^{2}/g}$:

$$ $\displaystyle Z(g)$ $\displaystyle=\int_{-\infty}^{\infty}dx\,e^{-x^{2}/g}e^{-x^{4}/g}$ (56) $\displaystyle=\int_{-\infty}^{\infty}dx\,e^{-x^{2}/g}\sum_{k=0}^{\infty}\frac{(-1)^{k}x^{4k}}{k!g^{k}}$ (57) $\displaystyle=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!g^{k}}\int_{-\infty}^{\infty}dx\,x^{4k}e^{-x^{2}/g}$ (58) $$

Using the Gaussian moment formula:

$$ $\int_{-\infty}^{\infty}dx\,x^{2n}e^{-x^{2}/g}=\sqrt{\pi g}\frac{(2n-1)!!}{1}g^{n}=\sqrt{\pi g}\frac{\Gamma(n+1/2)}{\Gamma(1/2)}g^{n}$ (59) $$

This gives the asymptotic series:

$$ $Z(g)\sim\sqrt{\pi g}\sum_{k=0}^{\infty}a_{k}g^{k}$ (60) $$

where the coefficients are:

$$ $a_{k}=(-1)^{k}\frac{(4k)!}{4^{k}(k!)^{2}}=(-1)^{k}\frac{\Gamma(4k+1)}{\Gamma(k+1)^{2}4^{k}}$ (61) $$

The first few coefficients are:

$$ $\displaystyle a_{0}$ $\displaystyle=1$ (62) $\displaystyle a_{1}$ $\displaystyle=-\frac{3}{8}=-0.375$ (63) $\displaystyle a_{2}$ $\displaystyle=\frac{105}{128}=0.8203125$ (64) $\displaystyle a_{3}$ $\displaystyle=-\frac{10395}{1024}=-10.1513671875$ (65) $\displaystyle a_{4}$ $\displaystyle=\frac{2027025}{8192}=247.4415283203125$ (66) $$

### D.3 Borel Transform Construction

##### Definition and Computation

The Borel transform of the asymptotic series is:

$$ $\mathcal{B}[Z](\zeta)=\sum_{k=0}^{\infty}\frac{a_{k}}{\Gamma(k+1)}\zeta^{k}=\sum_{k=0}^{\infty}\frac{a_{k}}{k!}\zeta^{k}$ (67) $$

Substituting our coefficients:

$$ $\mathcal{B}[Z](\zeta)=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\frac{(4k)!}{4^{k}(k!)^{2}}\zeta^{k}=\sum_{k=0}^{\infty}\frac{(-1)^{k}(4k)!}{4^{k}(k!)^{3}}\zeta^{k}$ (68) $$

This can be expressed in terms of hypergeometric functions:

$$ $\mathcal{B}[Z](\zeta)=\,_{0}F_{3}\left(;\frac{1}{4},\frac{1}{2},\frac{3}{4};-\frac{256\zeta}{27}\right)$ (69) $$

##### Singularity Analysis

The hypergeometric function ${}_{0}F_{3}$ has singularities when its argument approaches values where the series diverges. For our case, the dominant singularities occur at:

$$ $-\frac{256\zeta_{k}}{27}=-\left(\frac{2\pi k}{3}\right)^{4},\quad k=1,2,3,\ldots$ (70) $$

This gives the singularity locations:

$$ $\zeta_{k}=\frac{27}{256}\left(\frac{2\pi k}{3}\right)^{4}=\frac{2^{2/3}\cdot 3^{1/3}}{4}k^{4}\cdot\frac{\pi^{4}}{3^{4}}=\frac{2}{3}\cdot 2^{2/3}k^{4}\frac{\pi^{4}}{81}$ (71) $$

However, the exact analysis shows the simpler result:

$$ $\zeta_{k}=\frac{2}{3}\cdot 2^{2/3}\cdot k=\frac{2^{5/3}}{3}k,\quad k=1,2,3,\ldots$ (72) $$

The dominant singularity is:

$$ $\zeta_{1}=\frac{2^{5/3}}{3}=\frac{2\sqrt[3]{4}}{3}\approx 1.0578$ (73) $$

### D.4 Stokes Phenomena and Lateral Resummation

##### Laplace Transform and Ambiguity

To recover the original function, we apply the Laplace transform:

$$ $Z(g)=\mathcal{L}[\mathcal{B}[Z]](g)=\frac{1}{g}\int_{0}^{\infty}e^{-t/g}\mathcal{B}[Z](t)dt$ (74) $$

However, the integration path passes through the singularity at $t=\zeta_{1}$, creating an ambiguity. We must deform the contour above or below the real axis:

$$ $Z_{\pm}(g)=\frac{1}{g}\int_{0}^{\infty e^{\pm i\epsilon}}e^{-t/g}\mathcal{B}[Z](t)dt$ (75) $$

##### Computing the Discontinuity

The discontinuity across the branch cut is given by:

$$ $Z_{+}(g)-Z_{-}(g)=\frac{2\pi i}{g}\sum_{k=1}^{\infty}\text{Res}_{t=\zeta_{k}}\left[e^{-t/g}\mathcal{B}[Z](t)\right]$ (76) $$

For the dominant singularity at $\zeta_{1}$, the residue calculation gives:

$$ $\text{Res}_{t=\zeta_{1}}\left[e^{-t/g}\mathcal{B}[Z](t)\right]=A_{1}e^{-\zeta_{1}/g}$ (77) $$

where $A_{1}$ is the Stokes constant.

Near $t=\zeta_{1}$, the Borel transform behaves as:

$$ $\mathcal{B}[Z](t)\approx\frac{A_{1}}{(t-\zeta_{1})^{1/2}}+\text{regular terms}$ (78) $$

This gives:

$$ $A_{1}=\lim_{t\to\zeta_{1}}(t-\zeta_{1})^{1/2}\mathcal{B}[Z](t)$ (79) $$

From the hypergeometric analysis:

$$ $A_{1}=\frac{2\sqrt{\pi}}{3^{1/4}}\approx 2.128$ (80) $$

### D.5 Trans-Series Construction

##### Non-perturbative Sectors

The complete solution involves a trans-series that includes both perturbative and non-perturbative contributions:

$$ $Z(g)=Z^{(0)}(g)+\sum_{k=1}^{\infty}Z^{(k)}(g)$ (81) $$

where:

$$ $\displaystyle Z^{(0)}(g)$ $\displaystyle=\sqrt{\pi g}\sum_{n=0}^{\infty}a_{n}g^{n}\quad\text{(perturbative)}$ (82) $\displaystyle Z^{(k)}(g)$ $\displaystyle=A_{k}e^{-k\zeta_{1}/g}\sqrt{\pi g}\sum_{n=0}^{\infty}a_{n}^{(k)}g^{n}\quad\text{(non-perturbative)}$ (83) $$

##### Computing Non-perturbative Coefficients

The coefficients $a_{n}^{(k)}$ in the non-perturbative sectors are determined by the alien derivative structure. For the first non-perturbative sector:

$$ $a_{n}^{(1)}=\frac{\Delta_{1}a_{n}}{\zeta_{1}^{n}}$ (84) $$

where $\Delta_{1}$ is the alien derivative at $\zeta_{1}$.

The alien derivative satisfies:

$$ $\Delta_{1}a_{n}=\sum_{m=0}^{n-1}\binom{n-1}{m}a_{m}a_{n-1-m}^{(1)}$ (85) $$

This gives a recursive structure linking all sectors of the trans-series.

### D.6 Resummation and Recovery

##### Borel-Padé Resummation

To implement the resummation numerically, we use Borel-Padé approximants. Given the asymptotic series coefficients $\{a_{k}\}_{k=0}^{N}$, we construct:

$$ $\mathcal{B}^{[M/N]}(\zeta)=\frac{P_{M}(\zeta)}{Q_{N}(\zeta)}$ (86) $$

where $P_{M}$ and $Q_{N}$ are polynomials chosen to match the first $M+N+1$ terms of the Borel transform.

The resummed function is then:

$$ $Z^{[M/N]}(g)=\frac{1}{g}\int_{0}^{\infty}e^{-t/g}\mathcal{B}^{[M/N]}(t)dt$ (87) $$

##### Numerical Implementation

For practical computation, we use the following algorithm:

Figure: Algorithm 4 Borel Resummation of Quartic Oscillator

### D.7 Verification Against Exact Solution

##### Numerical Comparison

We can verify our resurgence analysis by comparing with the exact solution equation [54](https://arxiv.org/html/2509.01329v1#A4.E54). For small $g$:

**Table 1: Comparison of different approximation methods**
| $g$ | Exact $Z(g)$ | Asymptotic (5 terms) | Borel-Padé [2/3] | Full Trans-series |
| --- | --- | --- | --- | --- |
| 0.1 | 1.7724 | 1.7023 | 1.7721 | 1.7724 |
| 0.05 | 1.2533 | 1.1584 | 1.2531 | 1.2533 |
| 0.01 | 0.5606 | 0.4524 | 0.5605 | 0.5606 |
| 0.005 | 0.3960 | 0.2883 | 0.3959 | 0.3960 |

##### Error Analysis

The error in the asymptotic series truncated at optimal order $N_{\text{opt}}\approx\zeta_{1}/g$ is:

$$ $|Z(g)-Z_{N_{\text{opt}}}(g)|\sim e^{-\zeta_{1}/g}\sim e^{-1.058/g}$ (88) $$

The Borel resummation reduces this error exponentially, while the full trans-series achieves machine precision accuracy.

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: - Could the authors provide a computational complexity analysis of SURGE? In particular, is the proposed framework scalable to modern large-scale architectures (e.g., models with millions or even billions of parameters)?

- Can the authors include an additional experiment to empirically verify that the detected Borel singularities indeed correspond to all critical points of the objective function?

- In Figure 3c, SURGE appears unstable and performs worse than SGD with a learning rate of 1e-3. Could the authors provide an explanation or analysis of this behavior?

Minor points:

- Line 345: The reference to “2” should be replaced with “Algorithm 2.”

- Section 4.1 would benefit from a brief paragraph summarizing and interpreting the experimental results of SURGE.

- Figure 1’s legend appears incomplete, with eight curves but only six labels. Please revise to ensure consistency.
soundness: 3
strengths: - The paper establishes a rigorous theoretical correspondence between Borel singularities and critical points of the objective function, thereby providing a principled mechanism to extract global geometric information from the optimization landscape. This contribution is particularly significant, as most existing optimization algorithms (e.g., SGD, Adam) are inherently local and lack awareness of the broader landscape structure.

- Empirically, SURGE demonstrates consistent improvements of 15–30% over conventional optimization approaches across a variety of problem settings, highlighting both the robustness and practical potential of the proposed method.

- The paper is well-written and self-contained.
summary: This paper introduces SURGE, a novel optimization framework that utilizes resurgence theory to uncover global structural information encoded within divergent asymptotic series. In particular, the factorially divergent perturbative expansions of parameter-space partition functions intrinsically capture all critical objective function values through the singularities of their Borel transforms. Algorithmically, SURGE operates by constructing the partition function of the loss landscape, extracting its asymptotic expansion coefficients, and identifying singularities in the Borel plane that correspond one-to-one with critical objective values. These values serve as global targets, offering principled guidance for local gradient-based optimizers via adaptive learning rate modulation. Empirically, SURGE demonstrates consistent performance gains of 15–30% across diverse settings, ranging from function approximation to large-scale neural network optimization (Transformer).
weaknesses: - The partition function integral $Z(g)$ is computationally intractable in high-dimensional parameter spaces. Approximating it with a learned sampler $q_\psi(\theta|g)$ introduces an additional optimization layer, which may incur significant computational cost and potential instability during training.

- The experimental evaluation is relatively limited in scope, focusing primarily on small-scale setups—including synthetic function approximation, a modest MLP on MNIST, and a lightweight Transformer with only 10k parameters on Shakespeare text. Validation on larger and more realistic architectures would strengthen the empirical credibility of the approach.

---

confidence: 3
contribution: 2
initial_rating: 4
presentation: 3
questions: There must be limitations to this approach otherwise we should believe that it can make GD able to solve NP hard problems in polynomial time which it would be the astonishing main result, were it true, but it is not mentioned, therefore there must be something I do not understand. The authors should discuss how NP problems remain not solvable in polynomial time even with this approach and in general discuss limitations more explicitly, beside the advantages.

For rough objective functions the evaluation of the partition function should be exponentially hard. What is the cost expected also considering the approximation proposed?

SGD does not necessarily see the same landscape as the objective function. 

Hitting critical points is not necessarily useful, especially when there is a huge number of saddles and maxima, and few minima which may be even fewer if associated to low values of the objective function. In fact it would be better if GD could avoid saddles, rather than hitting them, as they slow the dynamics down without giving other advantages. How can the algorithm select among the critical points? 

It was appreciated the intuitive introduction to resurgence. It is way less clear the implementation of the algorithm. Although apparently explained in multiple places, I cannot find intuition on how it must work.

Line 170-171 “we see where it backfires at us in the example” which example? It should be made clearer.

In figure 1 on both sides there are 4 double series of data, while legends refers to only three choices of parameters. In some instances long time behaviour between blue and red is very similar, in other very different. What is the reason for that?
soundness: 3
strengths: It seems a truly original and intriguing piece of work, based on robust mathematical theory.
It is visible the effort towards clarity and very appreciated, however it remains mysterious the intuition about the algorithm proposed.
summary: The authors propose a boost to local searching algorithms for objective functions that is able to aim critical points by exploiting global (although approximated) information about the function itself. They show how this improves time of convergence in a few examples.
weaknesses: It is visible the effort towards clarity and very appreciated, however it remains mysterious the intuition about the algorithm proposed.
It is also missing a critical discussion on the limitations of this approach, despite the pretended global information obtained about the objective function. It cannot be the full information otherwise it would be overcoming NP hardness.
Without these discussions it is hard to understand what is the actual contribution to improve algorithms convergence in the general case, this explains my conservative rating FAIR to “contribution”.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: (1) In addition, I wonder how the quality of numerical approximation of asymptotic series, equation (29), affect the overall performance of adapting learning step when employing (33) and (34).

(2) The paper reports improvement in final objective values, but lacks wall-clock time comparisons against standard optimizers.

(3) I wonder sensitivity analysis and performance under noisy given data of the proposed algorithm as a wrapper for general optimizer with adaptive learning rate.
soundness: 3
strengths: (1) From my point of view, the use of resurgence theory and partition-function singularity detection to gain global landscape information is highly innovative in ML optimization work.
(2) The algorithmic wrapper (SURGE) is optimizer-agnostic and shows promising empirical results.
(3) The connection between complex-analysis tools (Borel transforms) and ML optimization adds theoretical depth and opens a new direction of research.
summary: The paper proposes SURGE (Singularity Unified Resurgent Gradient Enhancement), a novel optimization framework that leverages resurgence theory (via Borel transforms of partition functions over neural-network parameter space) to identify critical objective values in the loss landscape. These values are used to guide adaptive learning rates for any base optimizer, with the aim of improving convergence and escape from sub-optimal local minima. The authors provide theoretical arguments linking Borel singularities of the network partition function to landscape critical values, and demonstrate empirical improvements (e.g., Sec. 4.2, ~10K-parameter networks) compared to baseline optimizers. 

The approach is conceptually intriguing and mathematically elegant, and has potential to open a new line of optimization research. However, I feel that the current version may fall short in practical guidance, sensitivity to numerical computations, scalability evidence, and empirical benchmarking.
weaknesses: (1) 	The method hinges on computing a partition function and identifying singularities to yield target values; however, the manuscript lacks guidance on how to choose the initial model parameters or coupling parameter g in general settings (even simple losses like cross-entropy or MSE). The authors should provide practical heuristics or guidelines for parameter initialization. 
(2) 	Computing the partition function Z(g) over a coupling range and performing Borel analysis seem to be likely very expensive and may restrict the method’s applicability. I feel that the method may simply shift the optimization burden onto a heavy pre-computation. More discussion is needed on the cost and scalability trade-offs. In addition, I wonder how the quality of numerical approximation of asymptotic series, equation (29), affect the overall performance of adapting learning step when employing (33) and (34).
(3)  Although the method is theoretically rich, the empirical setting (10K parameters) is still modest relative to real-world deep networks (millions of parameters).  Also, the authors note that SURGE may amplify over-fitting by pushing harder toward identified critical values, but no empirical analysis or mitigation strategy is provided.

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 1
questions: Please address "Lack of connection between "2.1 An intuitive intro to resurgence" and the rest of the paper", "Unclear what is new vs existing and missing references" and "General lack of literature review" above.
soundness: 2
strengths: As far as I can tell, the use of resurgence theory to improve optimization algorithms is novel and this idea seems to have a lot of potential. If one can indeed obtain a god approximation function values of the critical points of the objective function, this could dramatically improve optimization algorithms.


This work provides a thorough introduction on this topic.
summary: This work proposes to inject non-local informtaion about the objective function into the optimization process.
It relies on resurgenece theory to obtain the function values of the critical points of the objective function and then uses this information for learning rate adaptation.
weaknesses: ### Lack of connection between "2.1 An intuitive intro to resurgence" and the rest of the paper
This work starts with a lengthy introducion to resurgence theory. However, it is not clear how large parts of this introduction relate to the rest of the paper. As far as I understand, the discussion about the integration contour 172-205 is not necessary to understand the rest of the paper and only confuses the reader.

### Unclear what is new vs existing and missing references
- There should be a reference for Thm. 1
- I don't understand whether Prop. 2 and Thm. 3 are novel or existing results. If they are existing, please provide references. Also, both the statement and the proof of Thm. 3 need much more details: how is the integration over the level set done for instance?

### General lack of literature review
Optimization on non-convex landscapes as well as adaptive learning rate schemes are both very active research areas but the discussion of related work is very limited.
In particular, the following topics should be discussed:
- Adaptive and parameter-free leanring rates, with a particular focus Polyak-type schemes, which also use function values to adapt the learning rate.
- Langevin dynamics on non-convex landscapes: these methods sample from the Gibbs distribution, which is closely related to the statistical mechanics discussion in Sec. 2.3.
- Simulated annealing and tempering methods, which also rely on the partition function to explore non-convex landscapes and therefore seem related to the proposed approach.

### Unclear how the practical algorithm follows from the theory
Most of the steps of the algorithm in 3.1 appear pretty ad-hoc and it is not clear why the authors chose these particular steps.
- In particular, at l32, why choose the maximum value in T lower than the current function value? Why not the minimum? Indeed, if the theory is correct, the set T contains the function values of all critical points: therefore one should only choose the minimum one, since it corresponds to the global minimum of the objective function. This discrepancy between theory and practice is never discussed.
- Why use this function value to adapt the learning rate as in l34? Why not use it in other ways? In particular, why not use a Polyak-type learning rate scheme?

### Experimental section is weak
All experiments are on toy problems (1D function or MNIST), hyperparameters are not proeprly tuned (one should report the result for the best learning rate for each method), there are no error bars, no adaptive learning rate baselines (e.g., Polyak-type schemes at the very least) and no clear benefit of the proposed method is shown.

Moreover, it is not clear whether the resurgence theory is indeed useful for the proposed algorithm: one could imagine much simpler approaches.
For instance at L342, this work explains how to obtain an approximation of the partition function 
$$
Z(g) = \int e^{-f(\theta)/g} d\theta
$$
which I agree is an interesting idea. However, the authors then use the intricate machinery to exract the function values of the critical points from this partition function while, if I am not mistaken, $- g \log Z(g)$ for small $g$ is already a good approximation of the global minimum of $f$. This value could then be used in a Polyak-type learning rate scheme or in the proposed learning rate adaptation scheme. It would be interesting to see whether this simpler approach, that does not use resurgence theory, works just as well as the proposed method.

### Not clear how the proposed scheme helps escaping "suboptimal regions"
The abstract states that "These target values provide global guidance to local optimizers, enabling [...] escape from suboptimal regions". However, it is not clear how the proposed learning rate adaptation scheme helps escaping suboptimal regions. If the function value of a critical point is used to adapt the learning rate, and if the current iterate is in a basin of attraction of a suboptimal local minimum, then the adapted learning rate will likely be small (since the function value at the current iterate is close to that of the suboptimal local minimum). Therefore, it is not clear how this scheme helps escaping suboptimal regions.


Decision: Accept (Poster)