1 Introduction↩︎

As a general principle for solving difficult problems, the Homotopy paradigm appears across diverse domains under different names, for example, Graduated Non-Convexity [1] and Gaussian homotopy [2] for optimization, homotopy continuation [3] for polynomial root-finding, and annealed Langevin dynamics [4] for sampling. Specifically, the Homotopy paradigm firstly construct an explicit homotopy interpolation from a simple, easily solved source problem to a complex target problem. Then, the solution of the complex problem is progressively approached by tracing the implicit trajectory along this interpolation path, effectively circumventing the challenges of direct solution.

Practical solvers for these problems often follow a predictor-corrector (PC) structure, where a predictor advances along the outer homotopy interpolation and a corrector iteratively refines the solution [5]. Despite their widespread use, these solvers rely on manually designed heuristics for step sizes and termination rules, which are typically suboptimal and task-specific. Furthermore, these methods have been independently developed in each domain, and no prior work has systematically unified these efforts under a single framework. We argue that this unification is crucial: it enables the design of a general solver that applies across problem instances, rather than requiring ad-hoc, per-problem solutions.

Building on this perspective, we propose Neural Predictor-Corrector (NPC), a plug-and-play framework that replaces heuristic rules with automatically learned policies. Instead of manually designed rules, NPC treats the choice of predictor and corrector strategies as a sequential decision-making process [6] and employs reinforcement learning (RL) [7] to adaptively learn effective policies. Crucially, we adopt an amortized training regime: a single offline training phase over a distribution of problem instances produces a policy that can be directly deployed on new instances from the same problem without per-instance fine-tuning.

We evaluate NPC on four representative homotopy tasks: Graduated Non-Convexity for robust optimization [1], Gaussian homotopy for global optimization [2], homotopy continuation for polynomial root-finding [3], and annealed Langevin dynamics for sampling  [4]. Through experiments on four representative problems, our approach is validated for strong generalization to previously unseen instances. Furthermore, the results reveal a dual advantage: our method not only consistently outperforms existing approaches in computational efficiency, but also demonstrates superior numerical stability, thereby underscoring the benefits of our proposed architecture.

In summary, our main contributions are as follows:

• To the best of our knowledge, we are the first to unify diverse problems, including robust optimization, global optimization, polynomial system root-finding, and sampling, under the homotopy paradigm, thereby revealing their common predictor-corrector structure across these problems. This enables a unified solver framework, rather than per-problem solutions.

• We introduce Neural Predictor-Corrector (NPC), the first reinforcement learning-based framework that automatically learns predictor and corrector policies, replacing hand-crafted heuristics with learned, adaptive strategies.

• Extensive experiments across multiple homotopy problems demonstrate that NPC significantly outperforms other methods in efficiency, while achieving higher stability and enabling efficient, training-free deployment on previously unseen instances.

2 Related Works↩︎

Although PC solvers appear across multiple domains, these lines of research have largely evolved independently. We review them here and highlight gaps that motivate our work. A full discussion of related works is provided in 11.

PC schemes trace solution trajectories along explicit homotopy interpolations. In robust optimization, Graduated Non-Convexity (GNC) gradually increases non-convexity to avoid poor local minima, with iterative solvers performing corrections [1], [8]. Gaussian homotopy methods construct progressively less smoothed objectives to track minimizers along bandwidth reduction [2], [9]–[11]. Polynomial system root-finding uses homotopy continuation with PC integration to trace roots from a simple start system [3], [12], [13]. In sampling, annealed Langevin dynamics and Sequential Monte Carlo define sequences of intermediate distributions with PC steps [4], [14], [15]. Across all these domains, predictor and corrector components are typically hand-designed, requiring per-instance tuning and limiting generalization.

Recent work has introduced learning into homotopy pipelines, showing efficient and effective improvements on Gaussian homotopy [16], sampling [17], combinatorial optimization [18], and polynomial root-finding [19], [20]. However, prior learning-based methods either focus on a single homotopy component or require specialized per-instance training.

RL has been applied to learn optimizers or adapt algorithmic parameters, showing benefits on some optimization and sampling tasks [21]–[26]. However, these works do not address the full predictor–corrector control problem across diverse homotopy classes, nor do they leverage amortized training to produce a single policy transferable across instances.

3 Homotopy Paradigm as a Unified Perspective↩︎

In this section, we introduce a unified perspective on diverse problems. We begin in 3.1 by introducing the homotopy paradigm, a general principle that underlies a wide range of problems. Next, in 3.2, we show that the corresponding practical solvers can all be instantiated within a common predictor-corrector (PC) framework. Finally, in 3.3, we discuss four representative problems together with their homotopy formulations and PC implementations, thereby illustrating the breadth and utility of this unified perspective.

3.1 Homotopy Paradigm↩︎

As shown in 1, the homotopy paradigm provides a general principle for solving complex problem \(g(\mathbf{x})\). Specifically, the homotopy paradigm defines a continuous interpolation \(H(\mathbf{x},t)\) from a simple source problem \(H(\mathbf{x},0)=f(\mathbf{x})\) with known solutions to a complex target problem \(H(\mathbf{x},1)=g(\mathbf{x})\). By tracing the implicit solution trajectory \(\mathbf{x}^*(t)\) along this interpolation as \(t\) varies from \(0\) to \(1\), one progressively transforms the source solution into the target solution. The source problem and interpolation are explicitly defined by the user, while the target solution is implicitly determined along the trajectory.

3.2 Predictor-Corrector Algorithm↩︎

While the homotopy paradigm specifies the abstract principle, an effective algorithm is needed to trace the implicit solution trajectory in practice. The PC method [5] provides such a concrete algorithmic framework. As shown in 2, PC decomposes trajectory tracking into two complementary steps:

• Predictor: Determines the next level of the homotopy interpolation and predicts the solution’s position at that level.

• Corrector: Iteratively refines the predicted solution to align it with the true solution trajectory, thereby preventing the accumulation of bias across levels.

The choice of predictor level schedule and corrector iteration count is often heuristic. Suboptimal settings can lead to inefficiency, instability, or failure to follow the trajectory accurately, motivating the development of adaptive or learning-based strategies for robust and efficient solution tracking.

3.3 Representative Homotopy Problems and Practical Solvers↩︎

To illustrate the breadth of homotopy paradigm applications, we describe four representative problems together with their corresponding homotopy interpolations and PC implementations.

Robust loss functions (e.g.., Geman–McClure [27]) mitigate the effect of outliers. However, they introduce strong non-convexity, increasing the risk of poor local minima. Graduated Non-Convexity (GNC) [1] addresses this challenge by defining a homotopy interpolation: \[\label{eq:GNC-GM} H(\mathbf{x},t) = \sum_i \frac{\bar{c}^2 \, {r}(\mathbf{x}, y_i)^2}{\bar{c}^2 + t\, {r}(\mathbf{x}, y_i)^2},\tag{1}\] where \(\bar{c}\) is a predefined parameter that controls the robustness of the GM loss, \(r(\cdot,\cdot)\) represents the residual function, and \(y_i\) denotes the measurements. This interpolation smoothly transitions from a convex quadratic loss (\(H(\mathbf{x},0)=\sum_{i=1}r(\mathbf{x},y_i)^2\)) to the original non-convex Geman–McClure loss (\(H(\mathbf{x},1)=g(\mathbf{x})=\sum_{i=1} \frac{\bar{c}^2 \, {r}(\mathbf{x}, y_i)^2}{\bar{c}^2 + {r}(\mathbf{x}, y_i)^2}\)). The predictor gradually increases non-convexity according to a predefined schedule, while the corrector refines the solution at each stage, often via a non-linear least squares optimizer (e.g.., Levenberg–Marquardt algorithm [28]). This homotopy strategy has proven highly effective in problems such as point cloud registration under severe outlier contamination [29]. Details are provided in 9.1.

Many optimization problems suffer from highly non-convex landscapes with narrow basins of attraction, making it difficult for solvers to converge to global or high-quality local minima. [10] address this challenge by progressively smoothing the target function through convolution with a Gaussian kernel \(\mathcal{N}(0, t \sigma^2)\): \[\label{eq:GH} H(\mathbf{x},t) = g(\mathbf{x}) \star \mathcal{N}(0, t \sigma^2),\tag{2}\] where \(\star\) denotes the convolution operator. This Gaussian smoothing enlarges the basin of attraction, allowing solvers to approach promising regions more reliably. The predictor progressively reduces the kernel bandwidth, while the corrector refines the solution at each stage. Details are provided in 9.2.

Root-finding for polynomial systems is challenging due to multiple solutions and computational complexity. [3] address this by starting from a source system \(f(\mathbf{x})=0\) with known roots and defining a linear homotopy: \[\label{eq:HC} H(\mathbf{x},t) = (1-t) f(\mathbf{x}) + t g(\mathbf{x}),\tag{3}\] tracing the solution trajectory from the source roots to the target roots. The predictor extrapolates the next solution along this path, while the corrector refines it using Gauss-Newton [30] iteration at each step, ensuring accuracy along the trajectory. Details are provided in 9.3.

Sampling from complex, high-dimensional distributions is challenging due to multi-modality and slow mixing. [4] address this by constructing a homotopy between a simple source distribution (e.g.., Gaussian) and the target distribution: \[\label{eq:ALD} H(\mathbf{x},t) \propto \exp\big(-(1-t) f(\mathbf{x}) - t g(\mathbf{x})\big).\tag{4}\] The predictor schedules the intermediate distributions, while Langevin dynamics acts as the corrector at each step, iteratively refining samples to match the current intermediate distribution. Details are provided in 9.4.

These examples collectively highlight the broad applicability of homotopy paradigm and the central role of predictor-corrector strategies, motivating the need for learning-based policy optimization.

4 Neural Predictor-Corrector with Reinforcement Learning↩︎

This section introduces the Neural Predictor-Corrector (NPC) framework, a general approach for homotopy problems that replaces heuristic step-size and termination rules with neural parameterizations learned via RL. As shown in 3, NPC reformulates the predictor-corrector process as a sequential decision problem: the predictor advances the homotopy level, while the corrector ensures accuracy, both guided by adaptive policies. We first present the NPC formulation (4.1), followed by its training with reinforcement learning (4.2).

4.1 Neural Predictor-Corrector↩︎

Classical PC algorithms differ across homotopy problems in how they define prediction and correction, yet share a key limitation: their step-size schedules and termination criteria are governed by fixed heuristics. Such heuristics fail to adapt to varying solution trajectories, where small steps are needed for sharp transitions but larger steps improve efficiency when the trajectory is smooth.

• Homotopy Level: Current position along the interpolation path.

• Corrector Statistics: Iteration count and attained tolerance from the previous step, capturing both convergence efficiency and deviation from the predicted trajectory.

• Convergence Velocity: Relative change in an optimality metric between consecutive levels, reflecting the speed of convergence. For optimization and root-finding, this is the relative change in the objective value. For sampling, it is the change in a statistical distance such as Kernelized Stein Discrepancy (KSD) [31] between the empirical sample distribution and the target distribution across consecutive levels.

Given the state \(s\), NPC outputs a two-part action \(a\):

• Step Size \(\Delta t\): Controls the predictor’s advance along the homotopy path.

• Corrector Termination: Convergence threshold \(\epsilon\) or maximum number of updates, balancing accuracy and efficiency.

As shown in 5, the NPC solver operates in an iterative loop to trace the solution path of a given homotopy problem \(H\). Each iteration consists of three key stages. First, a neural network (the NPC module) dynamically determines next actions for both the predictor and corrector. Second, the predictor advances the homotopy level to \(t_n\) and predicts the solution \(\mathbf{x}_n\) at this level. Third, the corrector iteratively refines this prediction until the convergence criteria are met. Finally, performance statistics are collected and fed back to the NPC module to inform its decisions in the next iteration, creating an adaptive, closed-loop system.

4.2 Reinforcement Learning for NPC↩︎

Because the predictor-corrector procedure is non-differentiable and early decisions influence the entire trajectory, supervised or self-supervised training is inadequate. These approaches would require assuming that local geometric structures of the solution trajectory remain consistent across instances, which rarely holds in practice. We instead employ RL, which inherently evaluates sequential decisions by their cumulative effect and enables learning policies that generalize across instances within the same problem class. The reward function is designed to promote both accuracy and efficiency:

• Step-wise Accuracy (\(r_t^{\text{acc}}\)): Encourages faithful trajectory tracking, based on convergence velocity or relative error change in the target problem.

• Terminal Efficiency Bonus (\(r^{\text{eff}}\)): Rewards shorter corrector sequences, formulated as \(T_{\max} - T\), where \(T_{\max}\) is a predefined upper bound and \(T\) is the total corrector iterations.

Consequently, the cumulative reward \(R\) for an episode is defined as: \(R=(\sum_{t=1}^T\lambda_1r_t^{\text{acc}}) + \lambda_2r^{\text{eff}}\), where \(\lambda_1, \lambda_2\) are scaling coefficients detailed in 9. This reward design enables agent to balance accuracy and efficiency across the homotopy trajectory.

Sequential decision-making in homotopy problems entails that early step-size choices affect all subsequent levels. Self-supervised learning fails in this context because measuring the future contribution of a step size is infeasible: it depends on the local geometric properties of the trajectory at future homotopy levels, which are unknown in advance. Relying on such assumptions risks overfitting to the training landscapes. This challenge is analogous to the dilemma discussed in [21], where although the problem domains differ, the core issue of long-term dependency and overfitting is similar.

Reinforcement learning, by contrast, inherently evaluates actions based on cumulative outcomes, allowing NPC to adapt to diverse solution trajectories without assuming consistent local geometry. Amortized training further improves generalization: by training over a distribution of problem instances, NPC learns a policy that can be applied efficiently to unseen instances within the same problem class.

5 Experiments↩︎

5.1 Implementation Details↩︎

Following the RL formulation in 4, we employ Proximal Policy Optimization (PPO) [32], an on-policy algorithm well-suited for continuous state and action spaces. Implementation is based on the open-source Stable Baselines3 library [33]. The policy and value functions are parameterized as multi-layer perceptrons (MLPs) with two hidden layers of 16 units each and ReLU activations. All other hyperparameters use the default values provided by Stable Baselines3. To account for varying problem formulations and noise levels across tasks, reward signals are scaled appropriately to ensure stable learning and comparability across tasks. Details are provided in 9. All experiments are conducted on a 12-core 5.0 GHz Intel Core i7-12700KF CPU and an NVIDIA GeForce RTX 3060 GPU, unless otherwise specified.

In all tables, Iter \(\downarrow\) records the total number of corrector iterations (rather than predictor iterations, which are more commonly used to measure progress in homotopy problems), and Time \(\downarrow\) reports runtime in milliseconds. The best results are bolded and the second-best results in 7 are underlined. All results represent the average over 50 independent trials.

5.2 Problem 1 : Robust Optimization via GNC↩︎

We evaluate NPC in the context of robust optimization using the GNC framework, comparing it against the classical GNC (Classic GNC) approach and the iteratively reweighted least-squares (IRLS) version [8]. The evaluation covers two spatial perception tasks with high outlier ratios: point cloud registration [34] (95% outliers) and multi-view triangulation [35] (50% outliers). Our NPC model is trained solely on the Aquarius dataset from the EPFL Geometric Computing Laboratory, demonstrating its cross-instance generalization capabilities.

Following the metrics defined in [36], we report the rotation error (\(E_R\)) and translation error (\(E_t\)) in ¿tbl:tab:GNC95point95cloud95simplified? for each method. Additionally, [tbl:tab:GNC95triangulation] presents the 3D point reconstruction error (\(E_p\)), defined as the Euclidean distance between reconstructed and ground-truth 3D points. As shown in ¿tbl:tab:GNC95point95cloud95simplified? [tbl:tab:GNC95triangulation], NPC achieves accuracy comparable to Classic GNC, whereas IRLS, tailored for a specific task, performs poorly on triangulation and lacks generalization. In terms of efficiency, NPC significantly boosts GNC’s performance: on point cloud registration, it reduces iterations by approximately 70-80% and runtime by 80-90% without compromising accuracy. These results demonstrate that NPC preserves the robustness of Classic GNC while substantially improving efficiency and generalization.

5.3 Problem 2 : Global Optimization via GH↩︎

We evaluate NPC in the GH setting for non-convex function minimization. We compare our method with two categories of baselines: (i) the single loop GH methods, SLGH\(_r\) (\(\gamma=0.995\)) and SLGH\(_d\) (\(\eta=10^{-4}\)) [10], (ii) the Gaussian smoothing method, PGS (\(N=20\)) [11],

5.4 Problem 3 : Polynomial Root-Finding via HC↩︎

We evaluate NPC in the context of polynomial system root-finding using HC. Experiments are conducted on two categories of tasks: polynomial system benchmarks [37]–[39] and a computer vision problem (UPnP [40]) for generalized camera pose estimation from 2D–3D correspondences. 9 lists the specific polynomial systems used, with the first entries as classical benchmarks and the last as computer vision task. We compare NPC-accelerated HC with Classic HC and Simulator HC [20]. Both Classic

5.5 Problem 4 : Sampling via Annealed Langevin Dynamics (ALD)↩︎

We evaluate NPC in the context of ALD for sampling from complex distributions. Target distributions include a 40-mode Gaussian mixture model (GMM), a 10-dimensional funnel distribution [41], and a 4-particle double-well (DW-4) potential [42]. The NPC agent is trained on the 10-mode GMM with randomly sampled coefficients to learn generalizable policies for accelerating ALD. We compare our method against classic ALD [4]

5.6 Ablation Study of RL State Components↩︎

To assess the contribution of each component in the RL state, we perform an ablation study on the six datasets used for the GNC point cloud registration task, retraining the NPC agent with one component removed at a time. As summarized in 13, removing any single state component causes the agent to adopt a more conservative strategy, resulting in an increased number of corrector iterations relative to the full state. This tendency typically manifests as the agent selecting smaller predictor step sizes or stricter corrector tolerances to ensure convergence in the absence of complete

5.7 Analysis of Efficiency-Precision Trade-off↩︎

We analyze the efficiency-precision trade-off by benchmarking our NPC-accelerated method against classical GNC and ALD. Classical approaches require manual tuning of homotopy parameters, resulting in a performance curve where higher precision typically demands more iterations. By contrast, our NPC-accelerated method bypasses this manual exploration by learning a policy that directly identifies an optimal operating point. This learned policy inherently balances the predictor step size and corrector tolerance to maximize efficiency at a given precision level. The practical benefit is visualized in 14. For both GNC and ALD tasks, the single point representing our method lies well below the classical trade-off curves, clearly illustrating a substantial reduction in iterations at comparable precision.

6 Conclusion↩︎

This paper introduces Neural Predictor–Corrector (NPC), a reinforcement learning framework for homotopy solvers. By unifying diverse problems, including robust optimization, global optimization, polynomial system root-finding, and sampling, under the homotopy paradigm, their solvers are shown to universally follow a PC structure. NPC replaces handcrafted heuristics with adaptive learned policies and employs an amortized training regime, enabling one-time offline training and efficient, training-free deployment on new instances. Extensive experiments demonstrate that NPC generalizes effectively to unseen instances, consistently outperforms existing approaches in computational efficiency, and exhibits superior numerical stability. These findings position learning-based policy search as a practical, generalizable, and efficient alternative to traditional heuristic strategies. Looking ahead, this paradigm opens promising avenues for extending homotopy methods to broader classes of optimization and sampling problems. Nonetheless, we also acknowledge its current limitation, which is discussed in 12.

7 Ethics statement↩︎

Our work unifies diverse problem domains governed by the homotopy paradigm into a single framework and, based on it, proposes a general, learning-based solver NPC. Our experiments are conducted on publicly available academic benchmarks and synthetic data, involving no human subjects or sensitive personal information. We do not foresee any direct negative societal impacts or dual-use concerns, as the primary application of our work is to provide a more efficient and robust tool for scientific inquiry.

8 Reproducibility statement↩︎

To ensure reproducibility, we specify the sources for all real-world datasets and the parameters used to generate synthetic data. In addition, 9 provides additional implementation details, covering the specific problem formulations and the hyperparameters used in our experiments. Our code and pretrained models will also be released to the public.

9 Implementation details↩︎

9.1 Details of Problem 1 : Robust Optimization via GNC ( 5.2)↩︎

9.1.1 The Graduated Non-Convexity algorithm↩︎

Optimization problems that can be formulated as least-squares can utilize the robust kernel from 1 , which is represented as: \[\label{eq:GNC95obj95func} \mathbf{x}^* = \min_{\mathbf{x \in \mathcal{X}}, t \in \mathcal{T}} H(\mathbf{x},t).\tag{5}\] The GNC algorithm utilizes Black-Rangarajan Duality [27] to reformulate 5 into: \[\label{eq:GNC95obj95func95reformulate} \mathbf{x}^* = \min_{\mathbf{x \in \mathcal{X}}} \sum_{i=1}\left[ w_ir^2(\mathbf{y}_i,\mathbf{x}) + \Phi_{H_t}(w_i) \right],\tag{6}\] where \(w_i\) is the weight of the \(i^{th}\) measurement \(\mathbf{y}_i\), and the function \(\Phi_{H_t}(\cdot)\), whose expression depends on the choice of the robust cost function \(H_t\), defines a penalty on the weight \(w_i\). When \(H_t\) is defined by 1 , \(\Phi_{H_t}(w_i)\) is defined by \(\Phi_{H_t}(w_i)=\frac{1}{t}\bar c^2(\sqrt{w_i}-1)^2\). Moreover, the weight can be solved in closed form as a function of only \(t\) and residual \(r\).

Reformulating the problem as 6 simplifies the prediction step to updating the each weight \(w_i\) using 7 , rather than predicting the optimization variable \(\mathbf{x}\). \[\label{eq:GNC95weight95update} w_i=\left( \frac{\bar c^2}{tr^2(\mathbf{x},\mathbf{y}_i)+\bar c^2} \right)^2\tag{7}\]

We correct \(\mathbf{x}\) using a nonlinear optimization method defined by 8 . \[\label{eq:GNC95corrector} \mathbf{x}^*=\min_{\mathbf{x} \in \mathcal{X}} \sum_{i=1}w_ir^2(\mathbf{y}_i,\mathbf{x})\tag{8}\] In our experiments, point cloud registration employs a Gauss-Newton corrector, while multi-view triangulation uses a more robust Levenberg-Marquardt (LM) algorithm.

For the point cloud registration task, the reward scaling is set to \(\lambda_1=10^3\), and \(\lambda_2=10^{-3}\). For the multi-view triangulation task, the reward scaling is set to \(\lambda_1=10^{-1}\) and \(\lambda_2=10^{-3}\) due to its noise scale being significantly larger than that of point cloud registration.

9.2 Details of Problem 2 : Global Optimization via GH (5.3)↩︎

9.2.1 The Gaussian Homotopy algorithm↩︎

The equivalent expression for 2 is given by: \[H(\mathbf{x},t) = \int g(\mathbf{x}+t*\sigma)k(\sigma)d\sigma = \mathbb{E}_{\sigma \sim \mathcal{N}(0, \mathbf{I}_d)}[g(\mathbf{x}+t*\sigma)]\] where \(k(\sigma)=(2\pi)^{-\frac{d}{2}}e^{\frac{-\Vert\sigma \Vert^2}{2}}\) is referred to as the kernel.

The prediction process in Gaussian homotopy is implicit, as we only modify the shape of \(H(\mathbf{x},t)\) by varying the predictor’s homotopy level \(t\).

The correction for \(\mathbf{x}\) is performed using a momentum method [43], with the gradient update formulated as \[\label{eq:momentum95GH} \begin{align} \mathbf{v}_{t+1} &= \nabla_{\mathbf{x}} H(\mathbf{x}_t, t) + \beta \mathbf{v}_t \\ \mathbf{x}_{t+1} &= \mathbf{x}_t - \alpha \mathbf{v}_{t+1} \end{align}\tag{9}\] where \(\mathbf{v}_t\) is the velocity vector, with the initial velocity \(\mathbf{v}_0\) set to the zero vector, \(\beta\) is momentum coefficients, controlling the influence of past gradients, and \(\alpha\) is the learning rate, which determines the step size of the update. We set \(\alpha=0.01\) and \(\beta=0.8\) in our experiment. As the analytical computation of the gradient \(\nabla_{\mathbf{x}} H(\mathbf{x}_t, t)\) is not feasible for some Gaussian homotopy functions, we employ a zeroth-order method to obtain a numerical approximation. The calculation formula is as follows [44]: \[\label{eq:GH95dH47dx} \nabla_{\mathbf{x}} H(\mathbf{x}_t, t) =\nabla_{\mathbf{x}}\mathbb{E}_{\sigma \sim \mathcal{N}(0, \mathbf{I}_d)}[g(\mathbf{x}+t*\sigma)] =\frac{1}{t}\mathbb{E}_\sigma\left[(g(\mathbf{x}+t*\sigma)-g(\mathbf{x}))*\sigma\right]\tag{10}\]

The reward scaling is set to \(\lambda_1=1\), and \(\lambda_2=10^{-3}\).

9.2.2 The non-convex function minimization benchmarks↩︎

\[f(\mathbf{x}) = -20 \exp \left( -0.2 \sqrt{\frac{1}{n} \sum_{i=1}^{n} x_i^2} \right) - \exp \left( \frac{1}{n} \sum_{i=1}^{n} \cos(2\pi x_i) \right) + 20 + e.\]

\[f(x, y) = (x^2+y-11)^2 + (x+y^2-7)^2.\]

\[f(x, y) = 10 + x^2 + y^2 - 9\cos(2\pi x) - \cos(2\pi y)\]

9.3 Details of Problem 3 : Polynomial Root-Finding via HC (5.4)↩︎

9.3.1 The Homotopy-Continuation algorithm↩︎

The polynomial system root-finding problem is modeled in the form of 3 .

The prediction of \(\mathbf{x}(t+\Delta t)\) is performed using the Padé approximation. The Padé approximation polynomial \(R_{n,m}(x) = \frac{R_n(x)}{Q_m(x)}\) has the following form: \[\label{eq:pade95poly} R_{n,m}(x) = \frac{p_0 + p_1x + \dots + p_nx^n}{1 + q_1x + \dots + q_mx^m} = \frac{\sum_{j=0}^{n} p_jx^j}{1 + \sum_{k=1}^{m} q_kx^k}.\tag{11}\] It is equivalent to the power series given in 12 . \[\label{eq:pade95power95poly} \psi(x) := \sum_{k=1}^{\infty} c_k x^k.\tag{12}\] The basic Padé approximation principle is that, given two integers \(m,n \in \mathbb{N} \cup \{0\}\), we can find two polynomials \(P_n(x)\) of degree at most \(n\) and \(Q_m(x)\) of degree at most \(m\), such that the difference \(Q_m(x)f(x)-P_n(x)\) has an order of approximation of at least \(n+m+1\). In fact, this is mathematically equivalent to the requirement: \[Q_m(x)f(x) - P_n(x) = O(x^{n+m+1}), \label{eq:pade95power95transform}\tag{13}\] where \(O(x^N)\) denotes a power series of the form \(\sum_{n=N}^{\infty}c_nx^n\).

In our implementation, we set \(n=2\) and \(m=1\). we can derive the following coefficients based on 13 : \[\begin{align} q_{1} &= -\frac{c_{3}}{c_{2}}, \\ p_{0} &= c_{0}, \\ p_{1} &= c_{1} + q_{1} c_{0}, \\ p_{2} &= c_{2} + q_{1} c_{1}, \end{align}\] where \(c_{0}=\mathbf{x}(t)\), \(c_1=\mathbf{x}'(t)\), \(c_2=\frac{1}{2}\mathbf{x}''(t)\), \(c_3=\frac{1}{6}\mathbf{x}'''(t)\). We obtain the derivative of \(\mathbf{x}\) by differentiating \(H(\mathbf{x}(t), t)\) with respect to t: \[\begin{align} \frac{\partial H}{\partial \mathbf{x}}\mathbf{x}'(t) &= -\frac{\partial H}{\partial t}, \\ \frac{\partial H}{\partial \mathbf{x}}\mathbf{x}''(t) &= -\left( \frac{\partial^2 H}{\partial \mathbf{x}\partial t}\mathbf{x}'(t) + \frac{\partial^2 H}{\partial t^2} \right), \\ \frac{\partial H}{\partial \mathbf{x}}\mathbf{x}'''(t) &= -\left( 2\frac{\partial^2H}{\partial \mathbf{x}\partial t}\mathbf{x}''(t) + \frac{\partial^3H}{\partial \mathbf{x}\partial t^2}\mathbf{x}''(t) + \frac{\partial^3 H}{\partial t^3} \right). \end{align}\] Consequently, \(\mathbf{x}(t+\Delta t)\) is according to the following equation: \[\label{eq:HC95x95pred} \mathbf{x}(t+\Delta t) = \frac{p_0 + p_1\Delta t + p_2\Delta t^2}{1 + q_1\Delta t}.\tag{14}\] If the denominator in 14 approaches zero, we revert to using a power series to predict \(\mathbf{x}(t+\Delta t)\). In this case, the prediction takes the form \(\mathbf{x}(t+\Delta t)=c_0 + c_1\Delta t + c_2\Delta t^2 + c_3\Delta t^3\).

We employ a Newton corrector in our experimental setup. At each iteration, \(\mathbf{x}\) is updated according to the following equation until the convergence criterion \(\Delta \mathbf{x} < \epsilon\) is met. \[\label{eq:HC95correct} \begin{align} \frac{\partial H(\mathbf{x},t+\Delta t)}{\partial \mathbf{x}}\Delta \mathbf{x} &= -H(\mathbf{x},t+\Delta t), \\ \mathbf{x} &= \mathbf{x} + \Delta \mathbf{x}. \end{align}\tag{15}\]

The reward scaling is set to \(\lambda_1=10^{-3}\), and \(\lambda_2=10^{-1}\).

9.3.2 The polynomial system benchmarks↩︎

\[\begin{align} f_0 &: \quad \left( \sum_{i=-n}^{n} x_i \right) - 1 = 0 \\ f_{k+1} &: \quad x_{-n}x_{n} + \left( \sum_{i=-n+1}^{n} x_i x_{k-i} \right) - x_k = 0 \quad (\text{for } k = 0, 1, \dots, n-1) \end{align}\]

\[\begin{align} f_0 &: \quad \sum_{j=0}^{n-1} x_j = 0 \\ f_1 &: \quad \sum_{j=0}^{n-1} x_j x_{(j+1) \pmod n} = 0 \\ \vdots \quad & \qquad \qquad \vdots \\ f_{n-2} &: \quad \sum_{j=0}^{n-1} \left( \prod_{k=0}^{n-2} x_{(j+k) \pmod n} \right) = 0 \\ f_{n-1} &: \quad \left( \prod_{j=0}^{n-1} x_j \right) - 1 = 0 \end{align}\]

\[\begin{align} x_i \left( \sum_{\substack{j=1 \\ j \neq i}}^{n} x_j^2 \right) - c x_i + 1 = 0 \quad \text{for } i = 1, \dots, n \end{align}\] In our implementation, we set \(c=1.1\).

\[\begin{align} 2n x_k - c x_k \left(1 + \sum_{i=1}^{n-1} \frac{k}{i+k} x_i\right) - 2n = 0 \quad \text{for } k = 1, \dots, n \end{align}\] In our implementation, we set \(c=0.51234\).

9.4 Details of Problem 4 : Sampling via ALD (5.5)↩︎

9.4.1 The annealed Langevin dynamic sampling algorithm↩︎

Annealed Langevin dynamics sampling obtains initial sample points from a simple distribution and uses a series of time-dependent potentials to control the update of the samples, as shown in 4 . Let the time-dependent potentials be \(H(\mathbf{x},t) \propto \exp\big(-(1-t) f(\mathbf{x}) - t g(\mathbf{x})\big)\).

The prediction process in ALD sampling is implicit, as we only modify the shape of \(H(\mathbf{x},t)\) by varying the predictor’s homotopy level \(t\).

In each iteration of the corrector, the positions of the samples are updated using the following formula: \[\label{eq:ALD95corrector} \mathbf{x}=\mathbf{x}-\frac{\xi}{2}\nabla_{\mathbf{x}} H(\mathbf{x}_t, t) + \sqrt{\xi}\sigma,\tag{16}\] where \(\xi\) is a Langevin step size, and \(\sigma \sim \mathcal{N}(0, \mathbf{I}_d)\) is a Gaussian noise vector.

The reward scaling is set to \(\lambda_1=10\), and \(\lambda_2=10^{-3}\).

9.4.2 Distributions↩︎

The 10 dimensions Funnel distribution defined as \[\label{eq:funnel} \begin{align} x_0 &\sim \mathcal{N}(0, \sigma^2) \\ x_i | x_0 &\sim \mathcal{N}(0, e^{x_0}), \quad \text{for } i=1, \dots, 9. \end{align}\tag{17}\] The funnel potential given as \[g(\mathbf{x}) = \frac{x_0^2}{2\sigma^2} + \frac{1}{2}\sum_{i=1}^9e^{-x_0}x_i^2\]

The DW-4 potential defined as \[g(\mathbf{x}) = \frac{1}{2\tau}\sum_{ij}a(d_{ij}-d_0)+b(d_{ij}-d_0)^2+c(d_{ij}-d_0)^4,\] where \(d_{ij}=\Vert x_i-x_j\Vert_2\) is the Euclidean distance between particles \(i\) and \(j\). In our implementation, we set \(\tau=1\), \(a=0\), \(b=-4\), and \(c=0.9\).

9.4.3 Metric↩︎

The Wasserstein-2 distance [45] is given by \[\mathcal{W}_2(\mu, \nu) = \left( \inf_{\pi} \int \pi(x, y)d(x, y)^2 \, dxdy \right)^{\frac{1}{2}},\] where \(\pi\) is the transport plan with marginals constrained to \(\mu\) and \(\nu\) respectively. In our implementation, we we use the Python Optimal Transport (POT) package [46] to compute this metric.

The Kernelized Stein Discrepancy [31] is defined as \[\begin{align} u_{q}(x,x^{\prime})=\;& s_{q}(x)^{\top}k(x,x^{\prime})s_{q}(x^{\prime}) + s_{q}(x)^{\top}\nabla_{x^{\prime}}k(x,x^{\prime}) \\ & + \nabla_{x}k(x,x^{\prime})^{\top}s_{q}(x^{\prime}) + \text{trace}(\nabla_{x,x^{\prime}}k(x,x^{\prime})), \\ \mathbb{S}(p,q) =\;& \mathbb{E}_{x,x'\sim p}[u_q(x,x')], \end{align}\] where \(s_q=-\nabla_\mathbf{x}g(\mathbf{x})\), and \(k(x,x')\) is a positive definite kernel. Specifically, we use the standard RBF kernel for KSD computation in this work.

10 Background on Reinforcement Learning↩︎

Reinforcement learning (RL) [7] provides a natural framework for learning adaptive strategies. It formulates sequential decision-making as an Markov Decision Process (MDP) with state space \(\mathcal{S}\), action space \(\mathcal{A}\), transition dynamics \(p(s_{t+1}|s_t,a_t)\), initial state distribution \(p_0(s_0)\), reward function \(r(s_t,a_t)\), and discount factor \(\gamma \in (0,1]\). The goal is to find an optimal policy \(\pi^*: \mathcal{S} \rightarrow \mathcal{A}\) that maximizes the expected cumulative reward along a trajectory \((s_0,a_0,\dots,s_T)\): \[\mathbb{E}\left[\sum_{t=0}^{T} \gamma^t r(s_t,a_t)\right].\]

11 Full Related Work↩︎

Although we mentioned in previous sections that methods from different fields essentially share the same predictor-corrector spirit, they have long evolved independently of each other. Our work is the first to unify these methods. In this section, we will review 1) classical predictor-corrector methods; 2) learning-based improvements on predictor-corrector methods; 3) efficient optimization and sampling methods via reinforcement learning.

1) Robust optimization: A core related technique is Graduated Non-Convexity (GNC), first proposed by [1]. This method employs a predictor-corrector approach with non-linear least-squares solvers to compute robust solutions. However, it relies on a hand-crafted, fixed iteration schedule, making it unsuitable for real-time robotics applications. Building on this work, [8] established a connection between GNC and the iteratively reweighted least-squares (IRLS) framework, based on which they designed a novel iteration strategy that achieved faster speeds in point cloud registration tasks [47]–[50]. Nevertheless, this strategy’s lack of generalizability to other problems remains its primary limitation. 2) Gaussian homotopy optimization: The underlying principle of this area was first introduced in [9]. More recently, [10] proposed a novel single-loop framework for the Gaussian homotopy method that simultaneously performs prediction and correction. Subsequently, [11] improved the algorithm’s convergence rate by adding an exponential power-N transformation prior to the Gaussian homotopy process. 3) Polynomial root-finding: Homotopy continuation [3], a numerical method for finding the roots of polynomial systems, uses a predictor-corrector scheme to track solution paths. Subsequent methods by [12] and [13] analyzed the properties of polynomials to introduce various improvements, enhancing the algorithm’s speed. 4) Sampling: In generative modeling, annealed Langevin dynamics [4], [14] utilizes a predictor-corrector method to sample from image probability distributions, where the correction step uses Langevin dynamics to restore samples to an equilibrium state. Similarly, Sequential Monte Carlo (SMC) methods [15] also apply a predictor-corrector approach to sample from posterior probability distributions, with a correction step that employs importance sampling to re-weight the samples.

1) Gaussian homotopy optimization: [16] is a novel model-based method that learns the entire continuation path for Gaussian homotopy. However, this approach requires specialized training for each problem. 2) Polynomial root-finding: Focusing on the more specific sub-problem of polynomial system root-finding, both [19] and [20] propose learning-based methods to determine the optimal starting system for homotopy continuation. 3) Combinatorial optimization: [18] proposes the Continuous Relaxation Annealing strategy, aiming to enhance unsupervised learning solvers for combinatorial optimization problems. 4) Sampling: [17] establishes a unifying framework based on path space measures and time-reversals, and proposes a novel log-variance loss that avoids differentiation through the SDE solver.

1) Optimization: [21] proposes a general framework by formulating an optimization algorithm as a reinforcement learning problem, where the optimizer is represented as a policy that learns to generate update steps directly, aiming to converge faster and find better optima than hand-engineered method. [22] utilizes reinforcement learning [51] to train an agent that dynamically selects the damping factor in the Levenberg-Marquardt algorithm [28] to accelerate convergence by reducing the number of iterations. 2) Sampling: [23] employs reinforcement learning to adaptively predict the denoising schedule via optimizing a reward function that encourages high image quality while penalizing an excessive number of denoising steps. [26] proposes a general framework named Reinforcement Learning Metropolis-Hastings, which aims to automatically design and optimize Markov Chain Monte Carlo (MCMC) samplers.

12 limitation and future work↩︎

One limitation of our work is that the NPC agent’s reward scale currently requires manual tuning for each problem instance based on its noise level to ensure stable and efficient training. The scale of step-wise rewards influences the training process’s convergence time, while an oversized terminal reward can nullify the guidance from step-wise rewards. This imbalance can prevent the agent from correctly tracking the solution trajectory, causing it to adopt myopic strategies to prematurely reach a terminal state. We conduct experiments on the point cloud registration task with different reward scaling factors. The results are shown in ¿tbl:tab:reward95scaling?.

To address this, two avenues for future work are promising. The first is to develop a mechanism that automatically adapts the reward scale. A more fundamental solution would be to investigate adaptive normalization techniques for the reward function, making the learning process inherently robust and eliminating manual tuning.

13 Full Experimental Results↩︎

This section provides complete experimental results, which are summarized in [subsec:GNC] [subsec:GH] [subsec:HC] due to limited space. ¿tbl:tab:GNC95point95cloud? shows the full results for the point cloud registration experiments via GNC, ¿tbl:tab:GH95non95convex95optimize95full? shows the full results for the non-convex function minimization experiments via GH, and ¿tbl:tab:HC95poly95benchmarks? presents the detailed results for the root-finding experiments on polynomial systems via HC. In addition, we present box plots for a subset of the experimental results in 15 to visually compare the different methods.

14 The Use of Large Language Models (LLMs)↩︎

We use LLMs to polish writing.

References↩︎

---

1. Equal contribution. \(^\dagger\) Project lead. \(\textsuperscript{\Letter}\)Corresponding authors.↩︎