Two‑Stage Learning of Stabilizing Neural Controllers via Zubov Sampling and Iterative Domain Expansion
Authors: Haoyu Li
Xiangru Zhong
Bin Hu
Huan Zhang
Venue: NeurIPS 2025 spotlight ()
Keywords: Learning based Control
Formal Verification
Stability Analysis
Primary area: general_machine_learning

Abstract: Learning-based neural network (NN) control policies have shown impressive empirical performance. However, obtaining stability guarantees and estimates of the region of attraction of these learned neural controllers is challenging due to the lack of stable and scalable training and verification algorithms. Although previous works in this area have achieved great success, much conservatism remains in their frameworks. In this work, we propose a novel two-stage training framework to jointly synthesize a controller and a Lyapunov function for continuous-time systems. By leveraging a Zubov‑inspired region of attraction characterization to directly estimate stability boundaries, we propose a novel training-data sampling strategy and a domain-updating mechanism that significantly reduces the conservatism in training. Moreover, unlike existing works on continuous-time systems that rely on an SMT solver to formally verify the Lyapunov condition, we extend state-of-the-art neural network verifier $\alpha,\beta$-CROWN with the capability of performing automatic bound propagation through the Jacobian of dynamical systems and a novel verification scheme that avoids expensive bisection. To demonstrate the effectiveness of our approach, we conduct numerical experiments by synthesizing and verifying controllers on several challenging nonlinear systems across multiple dimensions. We show that our training can yield region of attractions with volume $5 - 1.5\cdot 10^{5}$ times larger compared to the baselines, and our verification on continuous systems can be up to $40-10{,}000$ times faster compared to the traditional SMT solver dReal. Our code is available at https://github.com/Verified-Intelligence/Two-Stage_Neural_Controller_Training.

Two-Stage Learning of Stabilizing Neural Controllers via Zubov Sampling and Iterative Domain Expansion
Haoyu Li Xiangru Zhong* Bin Hu Huan Zhang
University of Illinois Urbana- Chobaa-Champpign haoyuli50:11inoos com, xiaagruu401linnois.. edu binhu7@illinois com, huanChuan- com
Abstract
Learning-base basen neural network (NN) control policies have shown impressive em- pirical performance. However, obtaining stability guarantees and estimates of the region of attraction of these learned neural controllers is challenging due to the lack of stable and scalable training and verification algorithms. Although previous works in this area have achieved great success, much conservatism remains in their frameworks. In this work, we propose a novel two-stage training framework to jointly synthesize controller and a Lyapunov function for continuous-time systems. By leveraging a Zubov-inspired region of attraction characterization to directly estimate stability boundaries, propose a novel training-datt sampling strategy and a Nomain-updating mechanism that significantly reduces the conser- vatism in training. Moreover, unlike existing works continuous-time systems that rely on an SMT solver to formally verify the Lyapunov condition, we extend state-of-the -f--of---tt--tt--- art neural network verifier e.B-CROWN with the capability of per- forming automatic bound propagation through the Jacobian of dynamical systems and novel verification scheme that avoids expensive bisection. To demonstrate the effectiveness of our approach, we conduct numerical experiments by synthe- several challenging nonlinear systems across
sizing and verifying controllers multiple dimensions, We show that our training can yield region of attractions with volume 5 1.5 10t times larger compared to the basclines, and our verification on continuous systems can be up to 40 10,000 times faster compared to the traditional SMT solver dReal. Our code is available at htttts://gttiieii com/ StageSttge_Neae Lecreeeeeural__ee Controlllonrrrlerrr. Corrraorrroolllrrraiiiinnng
Verified-
Introduction
Learning- based neural network control policies have demonstrated great potential in complex systems due to their high expressiveness [19 53.. However, the applications of these neural controllers in safety-critical real-world physical systems raise concerns since they typically lack stability guarantees. A promising method for quantifying safety is estimating the region of attraction (ROA), which refers to the set of initial states from which a system is guaranteed to converge to a desired equilibrium point. The computation of the ROA is often formulated as the search for appropriate Lyapunov function satisfying some algebraic conditions over its sublevel sets (26] (1). Due to the highly complex and nonlinear nature of neural networks, searching for Lyapunov function is not easy. Many existing works that are able to make stability guarantees about a system have some restrictions on the problem structure, for example, requiring linear polynomial parameterizations of the system, where -fuuarrs technique [37 47 30 7 can applied. To enable certification for general nonlinear systems, many recent works explore data-driven methods synthesize neural certificate [8) such as
* "quual contribution.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Lyapunov functions 4
55
44
48
35
24
34
29 23 barrier Functions B
B2 7 10
25
16
and contraction metrics 38
39
36
22
to certify the system's 's stability.
In this work, we aim to synthesize neural controllers of continuous-time systems together with a Lyapunov certificate. Unlike typical machine learning tasks, this problem requires that the Lyapunov and boundary conditions hold over the entire continuous domain of interest, rather than statistically achieving high accuracy on discrete test set. However, since at each iteration only small number of points can be selected, the problem forces us to adopt specialized strategies both to pick the training domain and sample training data from the domain. We argue that these two are the key factors for the conservativeness issue observed in many previous works. For selecting the training domain, many previous works simply select fixed domain of interest and enforce the Lyapunov condition over it. While some previous works have proposed a new formulation to overcome the conservativeness by introducing an extra term to verify forward invariance and allow the ROA to touch the training domain nontrivially 48 35,, we argue that this does not fully solve the conservativeness issue and leaves the training domain selection a hard hyperparameter tune.
In addition to the training domain, training data selection also plays an important role. Besides simple random sampling that typically does not suffice for certification, previous works dealing with similar problems [4] 6,
55 48] typically achieve co-learning with a counterexample-guidee (CEGIS) approach, which generates counterexamples of the Lyapunov condition as the training data, either using verifier [ such as an SMT or MILP solver or gradient- based approach [44] 48) While these approaches have led to great success, such approaches can typically only refine upon already worked controller/I_yapo Lyapuuuo function pair and thus rely on carefully tuned
initializations 4 48 This limitation makes these works both conservative in their estimation of the ROA and also very difficult to apply to new systems without much domain -spaiin-peccfic knowledge.
In our work, we propose two-stage training framework with carefully designed curriculum domain expansion and training sample selection that can yield significantly less conservative ROA estimation inspired by Zubov theorem [56 24 18 28], which provides a characterization of the true ROA with a PDE. During the first stage, we propose Zubov-inspire training data sampling strategy by combining points from both inside and outside of the current ROA estimation better guide the learning. Moreover, we eliminate the conservativeness introduced by the training domain by dynamically expanding it using trajectory information. This stage can provide us an almost- working controller and Lyapunov function that we can refine upon. The second stage then involves counterexample-guiie learning (CEGIS) to eliminate all the counterexamples inside the ROA estimation and thus obtain verifiable Lyapunov certificate. Importantly, our training also fully utilizes the physics-informed Zubov loss and can eliminate the conservatism introduced by the inner approximation nature of the normal Lyapunov function. As we shall demonstrate in our numerical experiments, our approach can robustly produce controllers with significantly larger ROAs compared to previous sfate-of-the art approaches in all the benchmarks with various control input limits.
As we deal with continuous- time systems, verification over Jacobian vector products also presents a challenge. Previous works addressing the Lyapunov stability verification for continuous-tine systems typically rely on SMT-, MIP-, or SDP-based solvers
24
49
such as dReal 15] or
we also utilize the of-he-t neural network 5T for formally verifying Lyapunov stability. The key challenge is to handle the Jacobian-related operators in dynamical systems, which have not been discussed or implemented in existing neural network verifiers. To ensure verification tightness, we design new linear relaxations for these Jacobian operators in continuous time system learning such as Tanh and Sigmoid. Furthermore, we introduce a novel verification scheme to dynamically adjust the specifications when noticing counterexamples without incurring much computational cost to avoid expensive bisections as in [8 24).
Z3 11 Inspired by the success in 48
35
verification tool a,3 R-R-CROWN 52
3
45
46
To summarize, our main contributions include:
We introduce a novel two-stage Zubov-based training framework with carefully designed curriculum training domain expansion and training sample selection that significantly reduces the conservatism during training.
We customize and strengthen the advanced neural network verification tool a,8-CROWN with the ability of handling many commonly used Jacobian-related operators in continuuus-time systems. Moreover, we propose a novel verification algorithm that can effectively avoid expensive bisections without sacrificing on soundness.

---

We demonstrate through numerical experiments that our training can yield region of attraction with volume 1.5 10 times larger than the baselines. Moreover, across the systems under evaluation which we can formally verify, our verification achieves 40 - 10,000 times speedup compared to the previous SMT solver dReal.
2 Background and Problem Statement
We consider a nonlinear continuous- time state-feedback control system
j g(r, x(x)) f(x)
(1)
where is the system state, y denotes the open-loop dynamics, denotes the controller, and denotes the closed-loop nonlinear dynamics. We denote the solution of the system starting with initial condition
2o by e(t, ro), and denote the equilibrium state and the control input as x* x*/u" respectively. We denote the sublevelset {z f(z) & c} of some function f as f and correspondingly f=o for the levelset. Our objective is to jointly search for neural network control policy ug together with a Lyapunov function Vo that can certify the stability of the closed-loop dynamics ① according the Lyapunov theorem. The goal is to learn controller that maximizes the region of attraction (ROA), which has the formal definition below.
Definition 2.1 (Region of Attraction (ROA) (2]]). Region of attraction for the system I is the set R such that lim,-z (tt, zо) x* for all Fo R.
Theorem 2.2 (Lyapunov Theorem [26] 21)). Given forward invariant S containing the equilib-
rium, if there exists a function V S + such that we have
V(**)= = -0
V(x) > (Vx (M=*))
(2)
VV(2) f(x) < 0 (VIES)
then S is subset of the region of attraction.
Therefore, the objective of this work can be formalized as the optimization problem
max Vol(S), subject to Theorem 2.2 Va
(3)
To maximize the possible ROA, we adopt the Zubov theorem, which states that function that satisfies stricter set of constraints than the Lyapunov condition can express the real underlying ROA R.
Theorem 2.3 (Zubov Theorem [56, 24 18). Consider the nonlinear dynamical system in ① and
we assume that x* = without loss of generality. Let A c IK" and assume that there exists continuously differentiable function W A - E and positive definite function W R R satisfying the following conditions.
W(0) = O and 0 < W(x) < 1, Vx AM10 W(x) as 2.4 |
(4)
(5)
(= f(x) (=) W(z)), Vz
(6)
Then x* is asymptotically stable with the region of attraction is exactly A.
Here, equation O is called the Zubov equation and can be utilized as the physics- innfore- loss during our training. It is noticed that there exists a ground truth value for the Zubov function that can be obtained through the following theorem.
Theorem 2.4 (Liu et al.
(24). We adopt the convention that divergence of integration yields the value op. Then for any a 0, W (xx) tanh(a 5h x)]l dt) solves the Zubov equation equation
(w f(x) a(1 ((=)) W [x)
We refer the readers to 24 for more theoretical details. This theorem enables us to obtain the ground truth value of W during the training without solely depending on the PDE residual loss, making the training process more robust, and can provide better learning dynamics for the Lyapunov function.

---

3 Methods
3.1 Learning of Neural Controller and Lyapunov Function
In this section, we introduce our training pipeline for synthesizing the neural controller and a corresponding Lyapunov function. As we have discussed, the selection of the training domain and training samples are the most crucial factors for successful training. We shall see that the Zubov theorem provides us not only the physics-inspired loss, but also principled way of data selection. During training, we parameterize the controller and Lyapunov function as
uo(2) (r) = c. tanh((N.(x)) NN,(x*) + tanh ttahh'( (("n) Vo(x) (NN,) (7)
This parametrization ensures that we have 0 & V8 < and the crucial property ug(x*) u²
which ensures that the origin is also equilibrium for the closed-loop system. Moreover, this parametrization constrains the control input lie in the range [-c, c.. Sometimes, if non-symmetric limit [0, is more natural, append another ReLU layer ug- This operation can also be applied to each coordinate separately to support mixed control-input limits.
3.1.1 ROA Estimation Stage
The first stage of colearning which we call ROA estimation stage, aims to set up a mostly work- ing controller and an estimate of the Lyapunov function that can be further refined. As we have discussed, in this stage the most important aspect is the selection of training data and training do- main.
Eq. (B) Outside Eq. (9)
V=c
Zuboy Guided Sampling: Intuitively, a good training data selection scheme should mix the data points inside and outside of the running region of attraction. Selecting only "outside" points risks collapsing the sublevel set back toward the origin, while sampling exclusively "inside" offers no sig- nal to enlarge it. In low dimensions, a well-chosen training box plus random sampling can approxi- mate this mix, but as dimensionality grows, the ROA occupies vanishing fraction of the domain, and naive draws will almost always only select the points outside of the ROA. For example, in the cartpole system, the ROA we found at the end of training only occupies less than 1% of the full training domain, and this will only be more pronounced in even higher-dimennsonal systems, and therefore naive sampling is infeasible.
Figure 1: An illustration of the key idea used in the ROA estimation stage. Notice the balanced training sample selection from both in and out- side of the ROA estimation V The training domain is also dynamically expanded to include all the states along a convergent trajectory.
In the traditional Lyapunov function framework, is difficult to obtain this running estimation the ROA. Zubov theorem, on the contrary, characterizes the exact region of attraction by the sublevel set y Intuitively, our training pipeline, since our parametrization already deems that V cannot exactly equal 1, we propose picking the training data by selecting batch of points lying in V and another batch with (ax) as close to 1 as possible, with & 1. To achieve this, we select data by performing projected gradient descent (PGD) [27] with the following two objectives
Vinterior (0) Rel.U(Vo(x)- -c)
(8)
Lounside (0) (o(xx) 1/.
The first interior loss will push points to the sublevel set V and the outside loss will push points further away from the region of attraction as shown in Figure I A detailed discussion on how to pick is presented in Appendix A..
Dynamic Training Domain Expansion: While previous works attempt to resolve the conserva- tiveness introduced by fixed training domain by allowing the ROA to touch non-trivially with it, argue that the best way is instead to dynamically expand the domain during training. We thus

---

propose to utilize the trajectory information to dynamically update the training box. We start from relatively small training region and periodically perform trajectory simulation starting from the points inside the current ROA estimation ve with the same in the data selection. We then expand the training region to include the farthest states that the converged trajectories can reach. To also benefit simple systems where all trajectories converge to zero in monotone fashion, we further enlarge the box uniformly for the final training domain. If there is no convergent trajectory for many iterations, manually enlarge the box for further exploration. This approach simultancously resolves conservativeness in ROA posed by fixed domain and also resolves the challenging forward invariance issue. A visualization of our approach can be found in Figure II where the original blue training domain is expanded further to incorporate the states along the convergent trajectories. The domain updating algorithm can be found in Algorithm 2 in appendix.
Loss Function:
The training loss can be easily adapted from the Zubov theorem. The objective of the training is to train the Lyapunov function and controller to encourage Vo to behave more like Zubov function. Therefore, naturally we could minimize (0) and train on the PDE residual of (G) and the data loss (2.4). In the bandling of data loss, previous works such as 42] often simulates the trajectory till convergence to compute So | 1: I dl. This, however, would lead to very lengthy training and add difficulty in the implementation of efficient batched operations since different trajectories have different convergence speeds. We thus propose to replace this term using a Bellman-type dynamic programming approach. We note that the following rewrite is possible:
W (a) tanh(a
lo(t, x)
(10)
=
tanh(a a
| 1 dt + tanh (W(((,))))
With this approach, only simulation time T is required. In practice. we found that this term can benefit training even if a very small T like T' = 0.01s is adopted. We also utilized a similar strategy as in (42) to decouple the learning of the Lyapunov function and controller using the actor-critic framework and impose the boundary condition Vo(g) ~ 1 for E 0(2() to guide the early phase learning. Moreover, to ensure that the second CEGIS stage is numerically stable, we minimize (0) only with a hinge-type loss. Let > O be a small positive number. The final loss function is weighted combination of the following components:
(0) = ReLU (0) e),
(11)
Lnde(()) (0)
1 ki
(VVo(a;) (Mi, detach((((x:)) + a(1 Vo(a:)) Vo (2x 1m)
(12)
k2
Vo(o.) - tanh(a - ['
Tedata (0) =
(0(1, (t)- - tanh tant-aannaaaaaaaaa (Vo(((((T
(13)
122
Leontroller (0)
detach( (onliVVVioe VVo(<s)) [ u(x:)),
(14)
k
ka Ve(w))-1)
Looundary (0) =
(15)
ka
where ki, K2 Kg represents the number of data points. Here, the PDE loss minimizes the PDE residual (6) by optimizing the Lyapunov function, whereas the controller loss optimizes the controller akin to a Sontag-style formula. Here &; are training data sampled by our novel sampling scheme, and y: are the boundary points from ((2)). More details about the loss function design can be found in
appendix A.3
3.1.2 CaGeS-based Refinement
Many previous works 48 start their training with a counterexample-guidde approach (CEGIS (2). However, without a good initialization, the counterexamples generated cannot be casily eliminated and will accumulate through the training and cause failure eventually, which is
5

---

shown in the ablation studies Table We found that our first stage can yield great initialization for CEGIS- based training. The controller after the pretraining stage is already almost working, meaning that almost all trajectories starting within Voee will successfully converge to the equilibrium. However, as we have discussed, no certification can be given this time since random samples typically are not sufficient to guarantee the condition holds over the entire continuous domain. Therefore, another stage of CEGIS-basek finetuning is performed to ensure the certifiability of the system. To find counterexamples, we follow the previous works in 48 44 attack with the objective
to use cheap gradient-baset PGD
(x) = min( (VVo(z) Vo(x) - c Vo(x))
(16)
Maximizing this objective leads to potential Teex that violates the Lyapunov condition in the band V,5 (V, Here require a lower e' since there is no guarantee that V²(0) and near the origin there must exist some small violations. However, we can choose this c' to be just slightly larger than Vo(0). These counterexamples then serve as the dataset for the CEGIS training. We then eliminate these counterexamples with the loss function
Leegis (0; Feex) (eLLU(VVoTraxs))(aaaa
(17)
where it enforces the Lyapunov condition. Moreover, the overall shaping of the Lyapunov function should mostly be formed at the ROA estimation stage due to the presence of Zubov-inspired loss and should thus be only slightly refined at the CEGIS stage. Therefore, we denote vo as the value of (0) prior to the CEGIS finetuning, and apply the regularization term
Log ()) = max(ReLU(r)) Vo(0)), ReLU(Ve(() - -d)) + Σ (V((r) - 1|
(18)
ZERO
with the aim of preventing the levelset of Vg from changing drastically during CEGIS. The first term thus stabilizes the levelset structure near the origin by restricting Vo(0) to stay close to the original value, whereas the second term stabilizes the levelset near the boundary of the training domain and encourages the ROA to stay away from this boundary to prevent trajectories from escaping from it.
Algorithm 1 Training of Neural Controller and Lyapunov Function
Require: Network Parameters Starting Domain Q Maximum Iterations M1, M2, Integration Time T and Discretizations Timesteps dt,, Learning rate Number of Training Points N, Domain Update Frequency y, Number of Trajectories Nr. Lyapunov Upper Threshold c, Domain Expansion Factor B. PGD steps P and step size a,, Finetuning learning rate a2 and epochs K >> Stage I: ROA Estimation Stage
1: for = 1 to M do then
2: if 7
3:
Ω UPDATEDOMAIN Ω2 T, dt, Nr, c, 3)
t See algorithm
4: I - N data points from 88) and (
5: Tboundary N points from (2
6:
- 0 - a LL x, Tooundar. T, df)
p See loss function section.
7: Dataset + Ø
t Stage 2: CEGIS Refinement
8: for 1 to .M do
9:
+ P points from f
10:
for j = to P do
11:
< + V Leee (x)
t- See equation (16)
Dataset - (Dataset, x
Buffer has length limit. See appendix
13:
for j = 1 to K do
14:
0 a: V L Dataset) Le (0))
See section 3.1.2
15:
if No CEX in dataset for at least n steps then
break
17:
return Trained parameters 6, Final Training Domain Ω
3.2
Verification
Criterion: Since there no guarantee that (0) = 0, counterexamples near the origin are unavoid- able. As the Lyapunov condition enforces the Lyapunov function value to decrease along trajectories,

---

we propose to exclude an unverifiable region defined as a sublevel set of the Lyapunov function. If the excluded set is small enough, this verification suffices for most practical needs. Formally, we propose the verification criterion
Theorem 3.1. Let 0- < c1 < cz. Suppose we have
x E (y n St - VV(x) g(z,)((a) < O, A(x) < 0,
(19)
2s V
(20)
where Ta(z) is the outer normal vector of in at z. Then both V and V are forward invariant, and every trajectory initialized in V enters the smaller set V in finite time.
The proof of this theorem can be found in appendix section B.I Here equation (19) states the normal Lyapunov condition has hold in the band {æ : Ci & Vax) E 02) whereas equation (20) states that the trajectory cannot escape through the boundary. Compared to previous works that simply exclude a rectangle near the origin 488, our verification scheme more principled since it also provides some guarantees inside the unverifiable region. In systems where formal certification over the full domain is challenging, verifying the condition within thin band where C2 - C1 for arbitrarily small 6 still allows one to certify the forward invariance of the sublevel set V
Formal Verification: For formal verification, instead of the commonly used SMT solver dReal in previous works [15] 4 24), we hope to use a state-of-the-at neural network verification framework, 51 [4], rigorously verify the Lyapunov condition through 43
such as a.8-CROWN 52 33 45
46
the tool's efficient linear bound propagation and branch-and-bound procedures. However, several key challenges remain for the application of a,3-CROWN. First, since the continuous-time Lyapunov condition involves Jacobian- veecorr products, it is essential to support efficient and tight bound propa- gation for the Jacobian of various dynamical system operators, which were previously not supported in any NN verifiers. Second, even after applying a CEGIS procedure, small counterexamples can still arise near the levelset boundaries at 6j and C2 To identify the tightest verifiable values of e and C2, a lengthy and expensive bisection is typically required, as commonly done in prior work [48] 24
In our work, we extend a,,3- D.33-CROWN with the capability to efficiently and tightly bound the Jacobian of many commonly used operators such as Tanh, Sigmoid. We implemented new linear relaxations that are much tighter compared to naively treating their Jacobian as composite functions of already supported operators, where the details can be found in the appendix section B.3 Furthermore, we enable our verifier to dynamically update c1, C2 during verification to avoid expensive bisection while remain sound. Intuitively, during the branch-and- bound process, we pick out those domains currently with the worst bound and try to find counterexamples in those domains. If counterexamples are found, we update ei and 0o to exclude the counterexamples. Notice that after this update, all the verified subproblems during branch and bound still remain sound and do not need to be verified again. This effectively saves the time needed for the repetitive verification needed if c₁ and 02 are to be found by bisection. The full verification algorithm can be found in appendix section B.4
Evaluation Schemes:
Compared to discrete-time systems, where the final verification is performed on through domain. Continuous-
the discrete Lyapunov condition (v V(g(a, u(x)) V(x)
time systems require formally verifying the Jacobian vector product (19). The large number of multiplications is challenging for current verifiers since it hard to bound multiplication well with a linear relaxation. Therefore, even with the strongest verification tool, verification is hard to scale to higher dimensions due to bound tightness issues. We thus define three evaluation schemes we considered. As system dimensionality increases, we must rely increasingly more on empirical checks.
1. Formal Verification:
This scheme implies the condition in theorem 3.1 is formally verified throughout the domain with a verifier.
PGD Evaluation: In this scheme, we use strong PGD attack to attempt finding counterexamples of the condition in theorem (3.1). If counterexamples are not found, we say the empirical PGD evaluation succeeds.
3. Trajectory Evaluation: In this scheme, we randomly sample trajectories from the sublevel set e and determine whether the trajectories converge to the desired equilibrium.
We can notice that these schemes are more and more empirical but are also more and more scalable to higher dimensional systems. The formal verification scheme, even though it is the strongest verification and can imply the other two schemes, can be hard to scale up. The empirical PGD
7

---

evaluation can scale up better, but still fails on our hardest 12-dimensional system. However, the trajectory based evaluation can be easily scaled up to any dimension. While this sampling approach provides no formal guarantees, several previous works utilize this technique to validate the learned controller (0
4
Experiments
Setups: We demonstrate the effectiveness of our methods in several challenging benchmarks including four 2D systems (Van-Der-Pol. Double- Integraor,, Inverted Pendulum, Path Tracking), one 4D system (Cartpole), three 6D systems (PVTOL, 2D Quadrotor, Ducted Fan), and 12D system (3D Quadrotor).
System details follow prior works 44 48 351 and are also provided in appendix C.2 We compare
our training framework against baselines from DITI. [44], Yang et al. 48] and Shi et al. 35 We
note that all of them work with discrete-time systems. We focus on comparing to these baselines
due to limitations in the existing continuous-time methods. The works 4
42 13
incorrectly handle
the condition u(x*) = u², so that the origin fails to be the equilibrium for most systems. Further discussions and examples are included in the appendix C However, to ensure fair comparison, we still run and present the comparison with Fossil 2.0 [3 for those systems that their method handle correctly. For each method, we run the veriiinn/veiiiiatioon pipeline five times with different seeds. We compare training success rates and ROA volumes, where the latter is estimated following previous works 44 (35] as the fraction of points within the ROA sublevel set among a dense sample from the full domain. All models undergo the three evaluation steps from section 3.2 in reverse order. More experimental details can be found in appendix C
OLD
1%
30
20
X0.0
mml
(a) Double Integrator.
(b) Cartpole (x, 0) slice.
(c) Cartpole (2,, .i) slice.
Figure 2: Region-of-Attractin (ROA) for Double Integrator and Cartpole compared to baselines.
The dotted box shows the baseline's training box being used. For the double integrator system, is clear that the initial training box is limiting the baseline's 's potential but is not constraining us.
Two and Four Dimensional Systems: On all the 2D systems, our methods robustly produce much larger ROAs compared to existing baselines. On both the path tracking and the inverted pendulum system, our training framework yields 5 times bigger ROA compared to the baselines as shown in table Moreover, on the more challenging systems such as inverted pendulum with small torque, our method is the only one that can get perfect success rate. Visualizations of the learned ROA for Double Integrator can be seen in the left of Figure 2 We also test our method on more challenging Cartpole system with four dimensions. Compared to the baselines, we obtain an ROA with more than 330x larger than the previous state of the art, as can be seen in Table [] The visualizations for different projections of the ROA for Cartpole can be found on the right of Figure 2. The rest visualizations and details can be found in the Appendix C.6
Six Dimensional Systems: We also implement our method on three higher dimensional systems: PVTOL, 2D Quadrotor, and Ducted Fan. For each system, our training pipeline yields much larger ROA robustly. Compared to the baselines, we obtain an ROA with more than 150000x larger than the previous state of the art on average for 2D quadrotor and more than 380x larger for the PVTOL system as can be seen Table the ducted fan case, we are the only method that can achieve an empirically verifiable ROA. We choose one projection to visualize for each of the system given in Figure More visualizations can be found in the Appendix C.6 Due to the difficulty of continuous-time system verification, formal verification for these three systems is challenging.
8

---

2,5
-25
20
-10
10
-6
(a) PVTOL
(b) 2D Quadrotor
(c) Ducted Fan
Figure Slice of the ROA for the 6D benchmark systems compared baselines.
Table Success rates (Succ) and ROA volume (ROA) across methods. Dashes (-) denote systems that are not already covered by the baseline and we were unable to obtain successful runs. The "Scheme" column indicates the evaluation scheme for each system as described in section 3.2 Formal
formal verification; PGD empirical PGD evaluation; Trajectory trajectory evaluation. "B" in system implies sig-torque setting, and "S" implies Small-torque.
Schene
RIM
Stece
Stcee
State
Pol
Forual
3.00 0.58
0075
8075
50%
3.4
100%
Dogble
Fomal
268. 1 20.3
9076
15
16.19
6275
1300
4.33
00%
362.5
10.7
100%
Pandalom B
Foral
106.
23.5
9075
100%
70.6
10075
487
80%
2846-6
140
100%
217.34
Tracking
Formal
100%
0.20
19075
15.3
00%
122.0+3..7
100%
Puth Tracking
Formal
14.86 9 0.18
100%
12.5 6.5
100%
78.9 8 12.0
100%
Canpole
Forual
0.021 0.012
1007
0.5256
20%
306.1
64.2
100%
PVTUL
PGD
80%
(1:01
0.28) -10"
100%
Quodiotex
PGD
100%
80%
(6.64
-10"
100%
Docted
(4.31
(0) 10"
100%
-.64)-10"
However, empirical verification is confirmed for these systems, and in 2D Quadrotor, we can formally verify the forward invariance of the sublevelset 5/0022 using Theorem 3.1
dum
15
-30
-20
-2.0
2.5
0.0
(a) Projection on (z, 2)
(b) Projection on (Ø.p)
(c) Projections on (0, (.
Figure 4: ROA estimation for 3D quadrotor. As
from the flow map, all the trajectories in
our ROA estimation converges to the origin.
Twelye Dimensional System: Finally, we our training pipeline in 12 dimension 3D quadrotor system, in which none of the previous methods successfully give any ROA estimation. In this high dimensional case, we found that our pretraining framework works smoothly and can consistently yield working controller. However, we found that CEGIS in higher dimensions does not work perfectly and we oftentimes found that counterexamples can still be located using longer PGD with more restarts even if prolonged CEGIS is already performed. However, the empirical success of our approach confirmed by sampling over 109 trajectories inside the levelset V and confirm that they all converge. The flow map in Figure also confirms this.

---

Table 2: Success rates of different training schemes under three verification regimes.
Evaluation scheme (% success)
System
Training scheme
Formal
Empirical-PGD
Trajectory
Two-stage with Random sampling NOA-Fstimation stage only (Zabov Domain Expansion)
100%
Caripole
CEGIS only
0
0 100%
0 100%
Ours (full pipeline)
100%
Two-stage with Random sampling
0
0 100%
2-D Quadrotor
ROA Estimation stage only (Zubov + Domain Expansion)
0
CEGIS only
Ours (full pipeline)
100%
100%
Verification Time: As we have discussed, our verifier struggles to verify higher-dimensional systems due to the bound tightness issue over Jacobian vector product. However, we shall note that our approach is already much faster compared to the previous state-of-the-an verifier dReal [T5] used by many previous work 4 [4] for continuous time systems. In the 2D systems verification, our method yields 40 10000 x faster verification compared to dReal. And for the cartpole system, the dReal
fails to terminate. The full verification time comparison can be found in Table 3
We also investigate the effectiveness of our training design choices. We mainly focus on the effectiveness of our training sample selection scheme and the two-stage training by considering three training settings:
A Ablation Studies:
1. Two Stage training with domain expansion but with random sampling.
2. Only the first ROA estimation stage with Zubov inspired sampling and domain expansion.
3. Only CEGIS and domain expansion, no ROA estimation stage.
We then investigate the success rate of these three training schemes under different evaluations with two higher dimensional systems Cartpole and 2D quadrotor. The overall result can be found in Table We can found that without our careful design choices, neither training Scheme I or Scheme 3 can yield any verifiable result. However, with our designed ROA estimation stage, even without CEGIS our trained controller and I. Lyapunov function can robustly pass the trajectory-based verification. Then together with CEGIS stage, our full training pipeline consistently yields controller/Lyapunov function that can be either formally verified or pass the PGD attack evaluation.
Table 3: Example verification time for each solver across our benchmark systems with formal verification. For 2D quadrotor, is marked with a star due to only forward invariance of the set V g022 is verified by involving theorem 3.1 with C1 0.2, = 0.21.
Method
Van
Double Integrator
Pendulum
Pendulum
Tracking
Tracking S Cartpole 2DQundd*
dReal
359.27s 3.00
1479.55 3.64s
1709.19s 3.94s
113.72s 3.77
142.20s 3.67
a.9-CROWN
3,94s
76443s
104614s
Conclusion
In this work, we propose a two stage training pipeline with novel training domain and training data selection for obtaining stabilizing controller with large region of attraction. Through extensive numerical experiments across various dimensions, we demonstrate that our training pipeline can achieve ROA volumes that are 5 - 1.5- 10 bigger than the baselines. We also extend the state-of- the-art verifier a, R.3-CROWN with the capability of performing bound propagation through many Jacobian operators and introduce novel verification schemes to avoid expensive bisection. We demonstrated that our verification can yield a 40- 10,000 times speedup compared the previously most commonly used verifier dReal. Although our method achieves promising empirical performance regarding the ROA volume and dimensionality compared existing baselines, formal verification is still challenging for higher dimensional systems. As future works, it will be interesting to explore more
training strategies and to further tailor and improve the verifier itself for this specific problem. Moreover, it would be interesting to extend the framework to more general control system settings, such as discrete-time or hybrid systems, slowly time-varying systems, or to synthesize and verify robust controller for systems with perturbations.
10

---

Acknowledgment
H. Zhang is supported in part by the AI2050 program at Schmidt Sciences (AI2050 Early Career Fellowship) and NSF (IIS-2331967). Hu generously supported by the AFOSR award FA9550- 23-1-0732.
References
[1] Alessandro Abate, Daniele Ahmed. Mirco Giacobbe, and Andrea Peruffo. Formal synthesis of Iyapunov neural networks. IEEE Control Systems Letters, 5(3):773-778, 2020.
[2] Alessandro Abate, Cristina David, Pascal Kesseli, Daniel Kroening, and Elizabeth Polgreen. Counterexample guided inductive synthesis modulo theories. In International Conference on Computer Aided Verification, pages 770-288. Springer, 2018.
[3] Aaron D Ames, Xiangru Xu, Jessy W Grizzle, and Paulo Tabuada. Control barrier function based quadratic programs for safety critical systems. IEEE Transactions on Automatic Control, 6288)7686,-3876 2016.
[4] Ya-Chien Chang, Nima Roohi, and Sicun Gao. Neural Iyарunov control. Advances in neural information processing systems, 32, 2019.
[5] Shaoru Chen, Mahyar Fazlyab, Manfred Morari, George J Pappas, and Victor M Preciado. Learning Yyарunov functions for hybrid systems. In Proceedings of the 24th International Conference on Hybrid Systems: Computation and Control, pages 1-11, 2021.
[6] Hongkai Dai, Benoit Landry, Lujie Yang, Marco Pavone, and Russ Tedrake. Lyapunov neural-network control. arXiv preprint
2021.
[7] Hongkai Dai and Frank Permenter. Convex synthesis and verification of control-lyapunov and barrier functions with input constraints. In 2023 American Control Conference (ACC), pages 4116-4123, IEEE, 2023.
[8] Charles Dawson, Sieun Gao, and Chuchu Fan. Safe control with learned certificates: A survey of neural łyаpunov, barrier, and contraction methods for robotics and control. IEEE Transactions on Robotics, 39(3)1749-1176,, 2023.
[9] Charles Dawson, Bethany Lowenkamp, Dylan Goff, and Chuchu Fan. Learning safe, generaliz- able perception-based hybrid control with certificates. IEEE Robotics and Automation Letters, 1944 2022.
[10] Charles Dawson, Zengyi Qin, Sicun Gao, and Chuchu Fan. Safe nonlinear control using robust neural lyapunov-barrie functions. In Conference on Robot Learning, pages 1724-174 1735. PMLR, 2022.
[11] Leonardo De Moura and Nikolaj Bjørner. Z3: An efficient smt solver. In International conference on Tools and Algorithms for the Construction and Analysis of Systems, pages 337-340. Springer, 2008
[12] Alec Edwards, Andrea Peruffo, and Alessandro Abate. A general framework for verification and control of dynamical models via certificate synthesis. arXiv preprint arXiv:2309.00 060000 2023.
[13] Alec Edwards, Andrea Peruffo, and Alessandro Abate. Fossil 2.0: Formal certificate synthe- sis for the verification and control of dynamical models. In Proceedings of the 27th ACM International Conference on Hybrid Systems: Computation and Control, pages 1-10, 2024.
[14] Mahyar Fazlyab, Manfred Morari, and George Pappas.
Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming. IFEE Transactions Automatic Control, 67(1):1-15, 2020.

---

[15] Sicun Gao, Soonho Kong, and Edmund M Clarke. dreal: An smi solver for nonlinear theories over the reals. In International conference on automated deduction, pages 208 214..4. Springer, 2013.
[16] Pushpak Jagtap, George J Pappas, and Majid Zamani. Control barrier functions for unknown nonlinear systems using gaussian processes. In 2020 59th IEEE Conference on Decision and Control (CDC), pages 3699-3704. IEEE, 2020.
[17] Pushpak Jagtap, Sadegh Soudjani, and Majid Zamani. Formal synthesis of stochastic systems via control barrier certificates. IEEE Transactions on Automatic Control, 66(7):3097-3110 2020.
[18] Wei Kang, Kai Sun, and Liang Xu. Data-driven computational methods for the domain of attraction and zubov's 's equation. IEEE Transactions on Automatic Control, 6903):600-1111,, 2023.
[19] Elia Kaufmann, Leonard Bauersfeld, Antonio Loquercio, Matthias Müllrr, Vladlen Koltun, and Davide Scaramuzza. Champion-leved drone racing using deep reinforcement learning, Nature, 620(7976):982-95) 2023.
[20] Alex Kendall, Yarin Gal, and Roberto Cipolla. Multi-tash learning using uncertainty weigh losses for scene geometry and semantics. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7482-7491, 2018.
[21] Hassan Khalil and Jessy W Grizzle. Nonlinear systems, volume Prentice hall Upper Saddle River, NJ, 2002.
[22] Haoyu Li, Xiangru Zhong, Bin Hu, and Huan Zhang. Neural contraction metrics with formal guarantees for discrete-ting time nonlinear dynamical systems. In 7th Annual Learning for Dynamics & Control Conference, pages 1447-1459. PMLR, 2025.
[23] Jun Liu, Yiming Meng, Maxwell Fitzsimmons, and Ruikun Zhou. Tool lyznet: A lightweight python tool for learning and verifying neural Iyapunov functions and regions of attraction. In Proceedings of the 27th ACM International Conference on Hybrid Systems: Computation and Control, pages 1-8, 2024.
[24] Jun Liu, Yiming Meng, Maxwell Fitzsimmons, and Ruikun Zhou. Physics-informet neural network Yyapunov functions: Pde characterization, learning, and verification. Automatica, 175:112193, 2025.
[25] Simin Liu, Changliu Liu, and John Dolan. Safe control under input limits with neural control barrier functions. In Conference on Robot Learning, pages 1970-1980. PMLR, 2023.
[26] Aleksandr Mikhailovich Lyapunov. The general problem of the stability of motion. International journal of control, 55(3):331-533 1992.
[27] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu. Towards deep learning models resistant to adversarial attacks. arXiv preprint orXiv:1706.0600 2017.
[28] Yiming Meng, Ruikun Zhou, and Jun Liu. Yubov-koopman learning of maximal Iyapunov functions. In 2024 American Control Conference (ACC), pages 4020-4025. IEEE, 2024.
[29] Yiming Meng, Ruikun Zhou, and Jun Liu. Learning regions of attraction in unknown dynamical systems via zubov-koopman lifting: Regularities and convergence. IEEE Transactions on Automatic Control, 2025.
[30] Antonis Papachristodoulou and Stephen Prajna. On the construction of Iyapunov functions using the sum of squares decomposition. In Proceedings of the 4 IEEE Conference on Decision and Control, 2002., volume 3, pages 3482-3482 3487. IEEE, 2002.
[31] Stefan Ratschan and Zhikun She. Providing a basin of attraction to target region of polynomial systems by computation of lyapunov-like functions. SIAM Journal on Control and Optimization, 4877):437744394, 2010.
12

---

[32] Alexander Robey, Haimin Hu, Lars Lindemann, Hanwen Zhang, Dimos V Dimarogonas, Stephen Tu, and Nikolai Matni. Learning control barrier functions from expert demonstrations. In 2020 59th IEEE Conference on Decision and Control (CDC), pages 3717-3724. leee, 2020.
[33] Hadi Salman, Greg Yang, Huan Zhang, Cho-Jui Hsieh, and Pengchuan Zhang. A convex relaxation barrier to tight robustness verification of neural networks. Advances in Neural Information Processing Systems, $2:9835-98466 2019.
[34] Mohamed Serry, Haoyu Li, Ruikun Zhou, Huan Zhang, and Jun Liu. Safe domains attraction for discrete- time nonlinear systems: Characterization and verifiable neural network estimation. arXiv preprint arXiv:2506. 13961, 2025.
[35] Zhouxing Shi, Cho-Jui Hsieh, and Huan Zhang. Certified training with suncchandd-bound A case study on lyapunov-star stable neural control. arXiv preprint arXiv:2411.18215, 2024.
[36] Dawei Sun, Susmit Jha, and Chuchu Fan. Learning certified control using contraction metric. PMLR, 2021.
In Conference on Robot Learning, pages
[37] Russ Tedrake, Ian R Manchester, Mark Tobenkin, and John Roberts. Lqr-trees: Feedback motion planning via sums-of-squares verification. The International Journal of Robotics Research, 1938):008--105, 2010.
[38] Hiroyasu Tsukamoto and Soon-Jo Chung. Neural contraction metrics for robust estimation and control: convex optimization approach. IEEE Control Systems Letters, 5(1::211-216, 2020.
[39] Hiroyasu Tsukamoto, Soon-Jo Chung, Jean-Jacques Slotine, and Chuchu Fan. theoretical overview of neural contraction metrics for learning-basee control with guaranteed stability. In 2021 60th JEEE Conference on Decision and Control (CDC), pages 1949-2955 IEEE, 2021.
[40] Hiroyasu T'ukamoto, Soon-Jo Chung, and Jean-Jacques E Slotine. Neural stochastic contraction metrics for learning- baased control and estimation. IEEE Control Systems Letters, 5(5) 1830, 2020.
[41] Giorgio Valmorbida and James Anderson. Region of attraction estimation using invariant sets and rational дyарunov functions. Automatica, 75:37-45, 2017.
[42] Jiarui Wang and Mahyar Fazlyab. Actor-critic physics-informi neural Iyapunov control. IEEE Control Systems Letters, 2024.
[43] Shiqi Wang, Huan Zhang. Kaidi Xu, Xue Lin, Suman Jana, Cho-Jui Hsieh, and J Zico Kolter. Beta-CROWN: Efficient bound propagation with per-neuron split constraints for complete and incomplete neural network verification. Advances in Neural Information Processing Systems, 34, 2021.
[44] Junlin Wu, Andrew Clark, Yiannis Kantaros, and Yevgeniy Vorobeychik. Neural Iyapunov control for discrete-time systems. Advances in neural information processing systems, 36:2939 2955, 2023.
[45] Kaidi Xu, Zhouxing Shi, Huan Zhang, Yihan Wang, Kai- Wei Chang, Minlie Huang, Bhavya Kailkhura, Xue Lin, and Cho-Jui Hsich. Automatic perturbation analysis for scalable certified robustness and beyond. Advances in Neural Information Processing Systems, 33, 2020.
[46] Kaidi Xu, Huan Zhang. Shiqi Wang. Yihan Wang. Suman Jana, Xue Lin, and Cho-Jui Hsich. Fast and Complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers. In International Conference on Learning Representations, 2021.
[47] Lujie Yang, Hongkai Dai, Alexandre Amice, and Russ Tedrake. Approximate optimal controller synthesis for cart-poles and quadrotors via sums-of-squares. IEEE Robotics and Automation Letters, 8(11 811):7376-7383 2023.
[48] Lujie Yang. Hongkai Dai, Zhouxing Shi, Cho-Jui Hsieh, Russ Tedrake, and Huan Zhang. Lyapunov-stable neural control for state and output feedback: novel formulation. arXiv preprint arXiv:2404.0795. 2024.
13

---

[49] He Yin, Peter Seiler, and Murat Arcak. Stability analysis using quadratic constraints for systems with neural network controllers. IEEE Transactions on Automatic Control, 6784):1980-1987, 2021.
[50] Zhexuan Zeng, Ruikun Zhou, Yiming Meng, and Jun Liu. Data-driven optimal control of unknown nonlinear dynamical systems using the koopman operator. arXiv preprint arXiv:2412.0100 2024.
[51] Huan Zhang, Shiqi Wang. Kaidi Xu, Linyi Li, Bo Li, Suman Jana, Cho-Jui Hsieh, and J Zico Kolter. General cutting planes for bound- pordrd-rrppgatoo--baeee based neural network verification. Ad- vances in Neural Information Processing Systems, 2022.
[52] Huan Zhang, Tsui-Wei Weng, Pin-Yu Yu Chen, Cho-Jui Hsich, and Luca Daniel. Efficient neu- ral network robustness certification with general activation functions. Advances in Neural Information Processing Systems, 31:4939 4948, 2018.
[53] Xiaoning Zhao, Yougang Sun, Yanmin Li, Ning Jia, and Junqi Xu. Applications of machine learning in real-time control systems: review. Measurement Science and Technology, 2024.
[54] Duo Zhou, Jorge Chavez, Hesun Chen, Grani A Hanasusanto, and Huan Zhang. Clip-and-verift Linear constraint-driven domain clipping for accelerating neural network verification. In The l'hirty-nitth Annual Conference Neural Information Processing Systems, 2025.
[55] Ruikun Zhou, Thanin Quartz, Hans De Sterck, and Jun Liu, Neural Цуаршпоv control of unknown nonlinear systems with stability guarantees. Advances in Neural Information Processing Systems, 15:29113-29111 2022.
[56] Vladimir Ivanovich Zubov. Methods of AM Lyapunov and their application, volume 4439. US Atomic Energy Commission, 1961.
14

---

A Training Details
A.1 Domain Update Algorithm
The full domain update algorithm can be found in algorithm 2 It first samples a batch of points from the sublevelset V by PGD with the objective (), and then performs trajectory simulation starting from these points. The algorithm then pick out those convergent trajectories and update the training domain to include the farthest state possible these convergent trajectories can reach.
Algorithm 2 UpdateDDmain
Require: Controller ug. Lyapunov Function Vo. Domain Ω2 Integration Time T and Discretizations Timesteps dt,, Number of Trajectories Nr, Lyapunov Upper Threshold c, Domain Expansion
Factor /
1: Finit Nr samples from the sublevel 0 Vo
Data from (8)
2: CeStet of indices where trajectory pit,, min) converges
3: for each dimension d do
4: dmin minj x/ miminaa
5: dmax + ma (o,t, i
6: Og - [dmin dmax]
7: Ω - 3 x f2
To handle monotonic trajectories.
8: return
A.2 CEGIS Dataset
The dataset from line 12 of algorithm [ will not grow without limitations. In practice, we fix maximal size k allowed for this dataset. After each iteration, sort the dataset in descending order according to the violation of Lyapumov condition, and only keep the top samples.
A.3 Loss Function
As we have discussed in the main text, the training loss is simply an adaptation of the Zubov theorem. The goal is to guide Vo to behave more like the Zubov function through training. We adopt a similar loss function as in [42] and modify the data loss using the novel derivation (TOp. Moreover, to ensure that the second CEGIS stage is numerically stable, we minimize Vo(0) only with hinge-type loss. In the loss functions, Iy are training data from our novel data selection (8) and ( Y; are boundary points from ((22) and this boundary loss would help during the initial phase of learning by ensuring that the controller does not drive the trajectory to escape from at least two times the training domain. To ensure balanced weight of the loss term and the hyperparameter tuning, adopt loss weights that are learnable using the techniques presented (20).
A.4 Hyperparameters
Compared to previous works, our training framework significantly reduces the difficulty of hyper- parameter tuning when adapting to new systems. For example, the training domain now does not require too much tuning, and the training almost always works when starting from a small domain. However, several key hyperparameters still require careful adjustment.
Firstly, the weights for each loss term must be chosen with care. practice, we have found that the boundary loss and the zero loss benefit from higher weights compared to the other three loss terms. Conversely, assigning relatively small weight to the PDE residual loss tends yield better results. We hypothesize that this is because an overly dominant PDE loss introduces complex optimization landscape, which can impede the learning of other critical components. Secondly, the integration involved in the data loss must be approximated using numerical scheme. Both discretization timestep and the number of simulation steps play an important role in ensuring stable and efficient learning. We found that simple forward Fuler scheme, with small timestep such as dt 0.001 and approximately 50 simulation steps, often sufficient. this context, as long dt is not excessively small, smaller dt generally performs better than a larger one, since
15

---

the dynamics can exhibit stiffness when the controller is only partially learned. Moreover, during trajectory simulations used for domain updates, smaller dt and longer steps are recommended. Using large timestep in this setting can cause otherwise convergent trajectories to be recognized as non-convergent. Finally, the training data selection, the choice of e crucial. We typically fix to be very close to 1, like 0.95. However, as the dimension of the system goes up, the Zuboy equation might not be fitted well, so we reduce to be around 0.9 for better training stability.
B Verification Details
B.1 Proof of theorem 3.1
Proof. Notice that if the trajectory is in (V*c)) [ (/ Ω the Lyapunov function value will keep decreasing along the trajectory. Indeed, by the mean value theorem, we have
V(x(1)) - ((=(0) (<((T) < <
(21)
First, we shall show that V n 92 is forward invariant. Assume by contradiction that there exists V Ω such that o(T, x)) & V Ω for some T O. We define the first exit time as T = inf(t: x(t,x) & V c 22. By continuity, we have ((T*, E 0(2 o Q2), which satisfies
"(T x) 0( o (2) = (22 n V U (2 o V=O)
If o(T", e (A V condition (0) says that it cannot exit through the boundary of Ω. Otherwise, suppose TTT*, x) interval (T* - 5, T we have V(x( 02, which contradicts the decrease Lyapunov function. Next, we show that on finite time the trajectory will enter Ω V Indeed, suppose for the contrary that the trajectory stays within (V V Ω2 the Lyapunov value will decrease indefinitely as we know by the above derivation (xx(4)) < ((x((2)) for 4 > 12, which cannot be possible by the definition of the set. Finally, we show that Sj is forward invariant. With the same proof, it is easy to show that for any 0 < E <<< C2 the set V is forward invariant. Now if S1 is not forward invariant, we know by definition that there exists E & and &
E f2 V-ra We know by definition and continuity that on a small time
such that (ol(,, f)) 0 + e This violates the forward invariance of / for any < c', leading to contradiction.
B.2 The G-B-CROWN verifier
a,, R.B8-CRWNN has been developed to be the state-of-the-ar neural network verifier, where formally, the toolbox aims to certify
F(z) < 0 for all B
(22)
with F being a general computation graph and B being hyper-rectangle.
High Level Intuition: a,33CROWN couples fast symbolic linear bound propagation with input- space branch-and-bound For box B, it constructs an affine upper bound of F and provide certification based on the maximum value of this affine bound. If inconclusive, it splits B into sub-boxes and repeats. As boxes shrink, local linear relaxations tighten and nonlinearities stabilize, so on small sub-boxes the affine bound is often already exact. The full process is very efficient, as each bound computation essentially costs a backward pass; thus many sub-boxes can be processed in parallel (tens of thousands), enabling rapid pruning of large regions and certification with only coarse bounds tiny boxes. Moreover, thanks to the tensorized nature of bound propagation and the flexibility of the branch-and bound procedure, a,, 3-COWN can naturally handle general specifications with "or" and "and". In principle, a,3-CROWN can handle all cases traditional SMT solver, such as dReal, can handle, and in much more scalable fashion. We refer the readers to the
following papers for more details [52 45 45 46
51].
For our application, since the condition in theorem 3.1 can be effectively rewritten as
(VV(s) g(,, (x)) < 0) v (V(x) > 02) v (V(x) < 02) (o(x,u())) কl az) < 0) (x AB) (V(x) c2),
(23)
(24)
and each term is a computation graph, a, 3 CROWN can efficiently verify this condition.
16

---

B.3 Linear Relaxation of de tanh(æ) and d sigmoid(2)
dx
To compute the Jacobian of neural Lyapunov function involving tanh(z) and sigmoid(z), by the chain rule, we need to propagate the gradient through the operators d tanh and sigmoid. Naturally, to use CROWN [52) to compute bounds for the Jacobian of the network, we need to linearly relax these operators so that we can propagate linear bounds through them. Given
di
d tanh(æ) dæ
tanh tanh?(x),
(25)
d sigmoid(()) sigmoid(z) - sigmold()))
dx
4 sigimoid can be reduced to computing the composition Finding linear relaxations of d tanh and de
dz
of the linear relaxations of tanh(z), (z), x², sigmoid(z)) and x y. The linear relaxations of these operators are already supported by auto LiRPA [45). However, in practice, this results in overly conservative relaxations. Thus, in our work, we directly derive a novel linear relaxation for these two gradient operators.
d sigmoid(rz)) are "bell-shaped" functions having similar monotonicity and con- d (a))( and da dz vexity/coneavity properties, as shown in figure [ Thus, we apply similar strategies finding linear relaxations for both of them. Suppose the function to be bounded is denoted oxx), with the bounds of æ as x E [1,u]. u]. We aim to find lower bound a' + 3 and an upper bound a x + B such that
a 3r & ((()) & a x zu, Vo E [1, u).
(26)
where parameters at, B4, a", 3B are dependent on and u. In our following discussion, for a "bell- shaped" functions (xx), we denote = as the two inflection points (i.e. the solution o² 0), and we assume that + 2 0 without loss of generality. In many cases, the linear bound is chosen to be a tangent line of a,, and thus is solely determined by the tangent point d' or d"
1.0
0.25
0.8
0.20
ox sigmoid(x)
tanh(x)
0.6
0.15
0.4
0.10
0.2
0.05
0.0
0.00
-5.0
-2.5
0.0
2.5
5.0
-5.0
-2.5
0.0
2.5
5.0
X
X
Figure 5:
de tanh(x) (left) and
4.
sigmoid(()) (right) have similar monotonicity and convex- de
d
ity/concavity properties.
The segment to be bounded is concave if -2, < Zg. In this case, the lower bound is chosen to be the line directly connecting two endpoints (l)) and (u, (u)) For the upper bound, d can be any point between I and u (e.g. the middle point). See (a) figure (a).
The segment to be bounded is convex if 2. < 1 & In this case, the upper bound is chosen to be the line directly connecting two endpoints. For the lower bound, d= can be any point between / and (e.g. the middle point). See figure (b).
The segment bounded similarly an "S-shaped" function if zo u. See figure (c).
For the lower bound, let de be the tangent point of the tangent line that passes through (I)). If u di, the lower bound is chosen to be the line directly connecting two endpoints. Otherwise, dl can be any point between dj and u (e.g. the middle point).
For the upper bound, let d, be the tangent point of the tangent line that passes through (u, ou)). If du. the upper bound is chosen to be the line directly connecting two
17

---

endpoints. Otherwise, dl can be any point between 0 and dy (e.g. linearly interpolated given ranging from u to d,).
0.25
0.20
0.20
0.20
[.iimoidto)
(igmmodlx)
0.15
0.15
0.15
0.05
0.05
0.05
-2.5
2.5
-2.5
(a)
(b)
(c)
Figure 6: Illustration of linear relaxations "bell-shaped" function (sigmoid). (a) -2, < l
Ng. The segment to be bounded is concave. (b) 2g & & The segment to be bounded is convex.
(c) l < 20 < u. The segment bounded similarly to an "S-shaped" function.
B.4 Verification Algorithm
Although a CEGIS procedure has been applied during training, small counterexamples can still exist. Previous works use bisection identify the optimal boundaries of the region where the Lyapunov condition can be verified. This requires test running the entire verification procedure every time a new candidate boundary is set, which is very inefficient. In contrast, our work, thanks to the formulation of the region to verify as levelset defined by Ci and c in [9] and the input-split branch-and- pipeline, can adaptively adjust the values of C1 and 02 in single verification run, avoiding the tedious and time-consuming bisection process.
Specifically, during the branch-and-bound procedure, for each batch of subdomains, apart from identifying each of them as VERIFIED or UNKNOWN based on its certified bounds, we additionally empirically search for counterexamples on some of them (usually the one with the worst certified bound). We shall note that searching for counterexamples in a subdomain by either random or PGD attack is much casier than doing in the whole initial domain, and thus is more likely to find counterexample that were missed out the CEGIS based training stage. We then shrink the levelset by increasing e₁ decreasing C2 to exclude all the counterexamples. The key intuition here is although the specification to verify updated, the subdomains that has been verified remain verified, don't 't need to rerun the whole algorithm from the beginning again. The full verification algorithm is shown in algorithm 3
C Experiment Details
C.1 Continuous Time Baselines
critical condition for theoretical guarantees is that the controller satisfies u(=**) u"; otherwise, x is no longer an equilibrium point, rendering the Lyapunov-based guarantee vacuous. In NLC [4] observed that the provided code mistakenly uses a bias-free controller for both the inverted pendulum and the path-following dynamics. While this is correct for the inverted pendulum, where x* incorrect for the path-following dynamics, as the chosen parametrization still enforces u(() 0, which does not hold in this case. This exact issue has also been noted an open issue NLC's of- ficial GitHub repository: https: //ttthub//uubb com Y'aChienChaar inChangg////// Neural Cyyapunov-Contre issues/14 For the method proposed by Wang al. [42]. we found no special treatment enforce this condition. In fact, when testing the publicly released checkpoints from their GitHub repository, we found that none of them satisfied ux=") u². For example, in the Double Integrator system, 2(0)= 0.3733, and in the Van der Pol system, u())= = 0.0633, despite both systems having =* 0 as their correct equilibrium input. Moreover, we believe there are some subtle issues formulation of Fossil 2.0. In [2). which underpins Fossil 2.0, we notice that their Certificate misses the boundary condition, which should supposedly ensure either {a V(x) & () X is empty or any trajectory starting from this intersection will remain in X. We have observed that their method can
In the main text, discussed the limitations of the continuous baselines 4 12 42
verification
18

---

Algorithm 3 Adaptive Branch-and- -ooundd Verification with Levelset Adjustment
Require: Neural network f. specification functions V and V = VV S. initial domain Ω initial levelset thresholds c1 C2, batch size B, boundary update margin I> stack of unverified subdomains
1: S (Sac
2: while S 2 and C1 < C2 do
3:
B - BS)
pop up B subdomains > unverified subdomains in this batch
4:
2
5:
for each D in B do
6:
Compute certified bounds of (xx) and V(x) over D if condition 19 is not verified to satisfy then add D to
7:
C1 and 02 used here
8:
9:
U + Ø then
10:
D Dorrst WORSTDOMAIN(41) Empirically search D worsst for counterexample Zee if countereeample found then
greatest positive upper bound of V(x)
11:
13:
UPDDATETTRESSIOLDSSCCC.. C1, C2, 6)
14:
for each D in 21 do Split D into D1, D2; PUSH(S, D1, D2)
branching
16:
S = Ø then
return VERIFIED
all subdomains have been verified
18: else
return FALSIFIED
C1 C22 the levelset is completely excluded
20:
function UPDNTETTIRESSOOLDSSC.C. C2,
21:
x- V(rce)
22:
if |e < 02 then
23:
C1 - E
counterexample near lower bound
24:
else
25:
2 - 1- -
t counterexample near upper bound
occasionally return false positives, i.e., a controller and Lyapunov function are claimed to be found while the controller cannot stabilize any trajectory. Therefore, on top of their returned results, we ran an additional check with trajectory simulations, and recorded negative for a run if the controller does not work.
C.2 Dynamical Systems
Van Der Pol: The dynamics is given by
d d2
a 21 + p(1 2,)22 +
with = 1. We choose the control input limit as |u] < for all the basclines.
Double Integrator: The dynamics is given by
32
and constrain the torque limit as |u| & 1. Inverted Pendulum: The dynamics is given as
81 = 8₂
8 m|2
B
9 sin((1) +
1
ml2
consider two set of torque limits as in small torgue case presents the more challenging limit |u] 1.02mgl.
35) The large torque case sets |u] & 8.15mgl and the
19

---

Path Tracking; The dynamics is given as
d sin)
(0s)
0.
(0.)
We also consider set of torque limits for this dynamics. We set |u] 1.68 in the large torque case and |u| < in the small torque case following [35]
Cartpole: The dynamics is given by
1
8 My sin A (u mp sin 0(10² cos 0))
0 [m + mp sin 0)
(-u cos m21 cos 0 sin (me + mp)gsin 0).
Following DITL 44). we set 1.0, m, = 0.1, = 1.0,9 = 9.81 and |u| & 30.
2D Quadrotor: The dynamics is given by
i 1
(in( + u2)
1 cos(f) (u1 + w22 g
0 [1 U2)
with 0.486,1 0.25, 7 0.00383, 9 9.81. We also set the control input RR have
the constraints |u| & 1.25mg as in 35
PVTOL: The dynamics is given by
Up COS Φ Ug sin Ø Ug sin U2 cos ó
0
0
0
0
0
<(t) =
U sin の
+
0
0
(27)
0 cos
where g -98,, m ||u| c < 39.2 as in 44.
4.0.7 0.25, J 0.0475. We set the control input u E RR2 to have the limit
Ducted Fan: The dynamics is given by
(- dā + ug cos u1 sin 0)
m
1
Ú (-di un sin 0 un cos 0) g
m
0 7 u0
with m 11.2, 0.156, 0.0462, 0.1, and g 0.28. We set the control input
Mop #1]
]_
TR2 to be un E 10, 10), E [0,10]] 10].
3D Quadrotor:
1 m
p
Reol (,,00,,,, 0,+)) -
ò we + tan ()) (sin cos())):)
0 cos(o)w, sin d)w;
sin( cos ( w)
cos(f) ())
w - (-w x (Iw) + T)
20

---

where
sin s sin ø + cos vi sin 0 cos Φ (,, = cos d sin Φ + sin sin cos ø cos cos Ø
The control input u IR* represents the four rotor thrusts. These are mapped to generalized forces via
Tx Ty
0
1
0
-1
1
-c
The hyperparameters specified with m 0.486,1
0.225, g
9.81, diag( Tr., Ty. Iz)
diag(0.0049, 0.0049, 0.0088), c torque limit of R. is given by |u| 3.6.
1.1
0.03793, n
mg 0.486x9.81
1.190325. The
29
4
C.3 Examples of Domain Update
We present some examples on the effect of the dynamic domain expansion mechanism and demon- strate its effectiveness. The initial domain for Quadrotor picked by 0.4 multiplied the original domain from 48 351, and the other ones are designed to be small but overall are picked quite arbitrarily without much tuning. A comparison between the initial training domain and the one after the first stage training is presented in Table 4 We can notice that
domain update mechanism is able to enlarge the training domain and thus the ROA by a large amount, resulting in much bigger and much more less conservative ROA estimation compared to the baselines.
Table 4: Example start and end domain for each dynamical system. Lower/Upper bound of the domain is given.
System
Start Domain
End Domain
Double Integrator
# 9.6]
Van der Pol
+ [1, 1
[ 8.5]
Pendulum (small torque)
11.2
= 12.66 69.6]
Pendulum (large uoque)
11
[30, 120]
Path Tracking (small torque)
+ 22
+ [9.6, 7.2]
Parh Tracking (large torque)
22
+[12, 12]
Carppoee
0.1.1..011,0.1.... 0.4]
=|48, 3.6, 13.2, 13.2]
Quadrotor
10033
1.6,
112
19.2, 20.4,
PVTOL
= 0.1.0............000
13.6. 6. 3.4,
Ducted Fan
0.44, 0.4, 0.4.
0...0.4.04440 0.4.
6. 7. 10.8,
4.8, 25.21
C.4 Examples of verifiable e₁ and 02
Table gives some examples of the final verifiable C1, thresholds. To demonstrate that the unverifible hole is small, we also report the volume of the sublevelset / and the volume occupies the full ROA. We can see that for all the systems the unverifiable region occupies only less than 0.1% of the full ROA. For the dither-dimensionll (6-D) systems, none of the 5 x 107 samples fell into the small unverifiable region.
C.5 Training Time
provides the training time for method. The time needed for the ROA estimation stage is relatively stable since it is trained to fixed epoch. However, the time taken by CEGIS will depend the ei chosen and the hardness of counterexamples near C1. In practice, found that though the CEGIS time varies, all training processes complete within reasonable time.
C.6 More ROA Visualizations
We also provide more visualizations on the final ROA obtained through our methods figure 7 Across all visualized slices, our method yield much larger ROA estimation compared to the baselines.
21

---

Table 5: Examples of verifiable thresholds C1, C2, and the volume of the unverifiable sublevel set {x V(x) ai f with its proportion of the full ROA. The ROA is estimated in Monte-Carlo fashion with 50000000 points.
System
e1
<2
Val(V(z) s (1)
Proportion
Double Iniegrator
0.01
0.99
0.038
0.014%:
Pol
0.0106
0.015
0.025%
Pendulum (small loeque)
0.0103
0.98
0.17
0.016%
Pendulum (large lonque)
Path Tracking (small torque)
0.0108
0.99
0.057
0.08%
Path Tracking (large torque)
0.0105
0.040
0.03%
Cartpole
0.0068
0.95
0.0154
0.003%
2D Quairotor
0.096
0.89
0%:
Ducted Fan.
0.0058
0.89
(%%
Table 6: Training time (in seconds) for each system. We report separately for the ROA Estimation stage, CEGIS stage, and a combined overall time in seconds.
System
Estimation
CEGIS
Full
Double Integrator
1173.20 30.83
162.39 + 97.61
1335.6 123.80
Van der Pol
1201.10 25.35
187.53 28.14
1388.64 44.38
Pendulum (small torque)
1169.27 + 139.2 29
225.14 + 64.17
1394.41 + 134.53
Peadulum (large torque)
2946.59
153.15
Path Tracking (smaal lorque)
789.30 1. 29.64
131.55 1: 54.61
920.85 1 62.67
Tracking (large torque)
48.19
982.35 1
Cartpole
4421.06 06 + 78.97
3182.58 01.558+4 879.25
7603.65 + 881.36
Quadrotor
5704.92 251.83
12030.59
7258.85 85
17755.52 + 7110.10
PVTOL
7100.86 192.40
10240.46
1554.67
17341.33 1 1570.56
Ducted Fan
6132.89 1932.89-11 + 108.35
12129. 13
+ 3119.57
18262 03 + 3187,87
Quadrotor
C.7
Computation Resource
The training was performed on NVIDIA V100 gpus, each with 16GB memory. The verification was performed NVIDIA 5090 gpus with 32GB memory.
22

---

202
16.0
74
7.5
Torgeta
1.2
5.0
23
255
8 2.5
0.0
-1.
-1.
-105
-30.0
2500
0.0
30
(a) Van Der Pol.
(b) Path Tracking large.
(c) Path Tracking small.
100
20
15
lag Kcal
10-
50
40-
0
-5
-20
-200
20
20
20
20
(d) Pendulum big.
(e) Pendulum small.
(f) Cartpole (2, .)) slice.
V000
2
-10
10
y
PVTOL (I, y) slice.
(h) PVTOL (y. g) slice.
(i) 2D Quadrotor (I, y) slice.
15
(mo8)
(m/0)
-10
-10
10'
2
-6
2
Lm/
(m)
(j) 2D Quadrotor (x, .)) slice.
(k) Ducted fan (I, 3) slice.
(1) Ducted fan (1r .)) slice.
Figure 7: More ROA visualizations.
23

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Please see method and experiments section.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We clearly discussed the limitations in the conclusion section.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g. if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-texx system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess..
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]

---

Justification:
One theorem about verification condition is included in theorem 3.1 Assump- tion is fully discussed, proof can be found in appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Complete algorithms are exactly discussed. Hyperparameters are discussed in appendix, and GitHub link containing codes we used can be found in abstract.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
25

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification:
The codes open-sourced. All data used is generated along the way (details in code), no additional data source is needed.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://niis.cc.
Pholicc/guddes/Coddsubbisss caapssionoolyyyy for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ wwtps/cccppplggguiiieeeeeec
Stteee//CoSSSSbmsssS.Si. wsw.sssoiccgyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All parameters can be found in codes. Hyperparameters are discussed also in appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Each method is trained 5 times with different seeds, and success rate is reported. Variance of the result is also reported.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
26

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The machine being used is discussed in appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http://tp/e.ea.e...cooo////aaaaa
Answer: [Yes]
Justification:
The research conforms the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- cration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There no particular socictal impact of the work that the authors can foresee that are important to discuss.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
27

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
No such risks are presented as far as the authors concern.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: Code is original, no other data is used.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
28

---

If assets are released, the license, copyright information, and terms of use in the thcode has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssit
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: All codes are clearly documented.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: There are no such experiments.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
There are no such experiments.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
29

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLMs are only used for polishing the writing.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: /ttpar////eupp.... http//rrppsscccoooffeereeeee//// Cc/Confeerecces/2225//II. for what should or should be described.
30


Reviews:
[{"review_id": "mYlXXwiayY", "reviewer": "Reviewer_Pjvg", "summary": "This paper considers the problem of neural control for stabilization tasks with formal guarantees on correctness. It proposes a method for jointly synthesizing a neural control policy and a neural Lyapunov function, where the goal is also to maximize the region of attraction (ROA) within which the resulting neural control policy is guaranteed to stabilize the system. Following several recent works, the present method follows a counterexample-guided inductive synthesis (CEGIS) based approach where a learner and a verifier module are composed in a loop. The novelty behind the present method is that it proposes a novel loss function for the learning module, which is inspired by Zubov's theorem characterization of the ROA. This loss function allows more explicit maximization of the ROA compared to prior methods. On the verification side, the method utilizes three verification schemes. The first verification scheme is a formal verification method that uses interval bound propagation combined with linear relaxations of involved non-linear functions which can then be achieved using $\\alpha$-$\\beta$-CROWN. The second and the third verification scheme are empirical in nature and do not provide formal guarantees. The method is implemented and experimentally evaluated on a number of benchmarks collected from the literature, and the experiments show that the method is able to learn neural control policies and neural Lyapunov functions that give rise to larger ROAs compared to the baselines.", "questions": "1. Questions about the verification algorithm (see Weakness 2 above): To what function exactly is interval bound propagation applied? How is this achieved when the dynamics function $f$ is non-linear and involves non-linearities beyond sigmoid or tanh functions?\n\n2. In the experiments, which benchmarks have been successfully verified by the first verification scheme that provides formal guarantees?", "limitations": "The paper should be precise about which benchmarks in the experiments were formally verified, given that formal guarantees constitute one of the main motivations behind the paper.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "**Strengths**\n\n1. *Problem relevance and novelty.* The paper considers an important problem of neural control policy synthesis with formal stabilization guarantees. While this problem has received considerable attention, the proposed loss function that is motivated by Zubov's theorem is a new idea that more explicitly reasons about the ROA, hence allowing its maximization. This is a novel feature that distinguishes the present method from most previous works on this topic.\n\n2. *Experiments (learning).* Experiments demonstrate that the method is able to learn neural control policies and neural Lyapunov functions that provide larger ROAs compared to the baselines. It should be noted that, apart from lower dimensional benchmarks, these ROAs are not formally verified (see weaknesses below).\n\n**Weaknesses**\n\n1. *Presentation and the lack of clarity.* The technical sections of the paper are very hard to follow and only provide a high level overview. To actually understand any details of the learning procedure one needs to read the appendix, and even after reading the paper + appendix there are several aspects of the method that are unclear to me. This makes the paper non-self contained and ambiguous, and I believe that the presentation requires a major revision for the paper to be publishable : \\\n** For selecting training points inside and outside ROA, in l.159-165 it is proposed to take points with V^{<=c} as inside ROA and points at which V is close to 1 as outside ROA. But how is the value of this threshold c determined? Also, wouldn't this result in points outside ROA being very far away from the ROA boundary? \\\n** The loss function used for training is never properly defined in the paper and is deferred to the appendix. In l.181, it is simply said that the loss function is easily adapted from Zubov's theorem. This is insufficient for a paper to be self contained. While I eventually understood the loss function, it took me a while to decipher it since the appendix provides no intuition on how each loss term corresponds to conditions in Theorem 2.3 and 2.4. For instance, it is never properly mentioned that the positive definite function $\\omega$ in Theorem 2.4 is set to $\\omega(x) = ||x||^p$. \\\n** I could not understand the part on the CEGIS loop in Section 3.1.2 at all, especially how counterexamples are incorporated into the loss function, even after several reads.\n\n2. *Ambiguity about the verification procedure.* The verification procedure in Section 3.2 is described by basically saying that the method uses interval bound propagation and linear relaxation of non-linear functions such as sigmoid or tanh, which can then be delegated to $\\alpha$-$\\beta$-CROWN. But no further details are provided. To what function exactly is interval bound propagation applied? How is this achieved when the dynamics function $f$ is non-linear and involves non-linearities beyond sigmoid or tanh functions? At this point, I cannot say that I understand the verification procedure at all.\n\n3. *No guarantees around the origin.* The verification method does not provide formal guarantees on stability around the origin, see l.218-222. This is a disadvantage compared to some SMT-based methods that do address this issue but are not discussed in the paper, e.g. \\\n[1] Abate et al. \"Formal Synthesis of Lyapunov Neural Networks\". IEEE Control Systems Letters 5 (3), 2020\n\n4. *Experiments not formally verified for larger benchmarks.* The paper is motivated by the need for formal guarantees and achieving them without reliance on SMT solvers. However, when it comes to designing the verification module, the paper considers three verification schemes of which only the first one provides formal guarantees. The second and the third verification scheme are empirical and do not actually provide guarantees. While Section 4 is ambiguous about which verification scheme was applied to which benchmark, from the text my understanding is that only 2D and 4D benchmarks were verified by the first verification scheme. But this is in line with what existing SMT-based methods can formally verify as well, e.g. reference [8] in the paper which is not considered in the experiments. This should be properly clarified, also in the abstract and into.\n\n**Minor comments**\n- l.29 searching the Lyapunov function -> searching for the Lyapunov function\n- l.75 the the -> the\n- l.100 levelset -> level set\n- l.110 function satisfies -> function that satisfies\n- l.117 can obtained -> can be obtained\n- eq. (14) and l.255 is $g$ supposed to be $f$?", "was_revised": true, "final_justification": "The author response addresses and clarifies all my questions and I do feel more positive about the paper now, especially with the question about which experiments have been formally verified being clarified. I have increased my score accordingly.", "rebuttal": "We appreciate the reviewer’s careful reading, which uncovered several places where our exposition is unclear. While the current draft could undoubtedly be clearer, we believe the main concerns arise from presentation gaps rather than from substantive flaws. Below, we address each point to clarify misunderstandings, and we promise that these clarifications will appear in the revision. In addition, we also conducted **additional experiments** to compare with a new SMT baseline, fossil.\n\n### **W1: Clarity on training**\n\n***Value c & Outside ROA data selection.*** In our real experiments, this $c$ can be chosen to be 0.99 for all 2D systems and needs to be slightly more conservative for higher-dimensional systems (0.95 for 4D and 0.9 for even higher-dimensional systems). For the selection of data points outside the ROA, we shall notice that they are still constrained **within the training domain** and cannot diverge to infinity since we adopted **projected gradient descent**. Moreover, the exact location of these points is not crucial as long as they stay off the ROA. Since the Zubov equation provides a global condition (see loss function $L_{pde}(\\theta)$ and $L_{data}(\\theta))$, the data points will still provide information even though they are outside of ROA.\n\n***$\\omega$ in Theorem 2.4.***\nWe indeed made the assumption $\\omega(x) = ||x||^p$, but we mistakenly deleted it. We will make sure to add it back. \n\n***Loss function.*** The loss function is a weighted combination of the loss terms as in Eq. (16) - (20) as follows:\n\n$L(\\theta) = a_1 L_{zero}(\\theta) + a_2 L_{pde}(\\theta) + a_3 L_{data}(\\theta) + a_4 L_{controller}(\\theta) + a_5L_{boundary}(\\theta)$\n\nwhere we have\n\n* $L_{zero}(\\theta) = ReLU(V_\\theta^2(0)-0.01)$. This term **minimizes $V_\\theta(0)$ akin to Theorem 2.3**. We use a hinge loss to avoid potential numerical issues in the following CEGIS loop.\n* $L_{pde}(\\theta) = \\frac{1}{k_1} \\sum (\\nabla V_\\theta(x_i)^\\top f(x_i, detach(u(x_i))) + a(1+V_\\theta(x_i))(1-V_\\theta(x_i))||x_i||^p)^2$. This term **minimizes the PDE residual according to the Zubov equation in Theorem 2.4**. We detach the controller to let the Zubov Theorem only shape the Lyapunov function.\n* $L_{data}(\\theta) = \\frac{1}{k_2} \\sum (V_\\theta(x_i) - \\tanh(a\\int_0^T ||\\phi(t,x)||^p dt) - \\tanh^{-1}(V_\\theta(\\phi(T,x))))^2$. This term further **trains the Lyapunov function according to equation (10) in line 188**.\n* $L_{controller}(\\theta) = \\frac{1}{k_1}\\sum detach(\\nabla V_\\theta (x_i))^\\top f(x_i, u(x_i))$. This term **encourages the learned controller $u$ to follow a Sontag-style formula** with the current Lyapunov function.\n* $L_{boundary}(\\theta) = \\frac{1}{k_3} \\sum (V_\\theta(y_i) -1)^2$. Here the $y$'s are chosen from $\\partial(2\\Omega)$ where $\\Omega$ is our training domain. This term **regularizes early-stage training** by driving $V_\\theta(y_i)$ to $1$, and thus prevents the trajectory from diverging too early as $L_{controller}(\\theta)$ trains the system to have decreasing Lyapunov function value along trajectories.\n\nWe will move the loss function discussion to the main text and provide these intuitions to ensure the paper is self-contained.\n\n\n***CEGIS Loop.*** In Algorithm 1 line 13, we construct a **counterexample dataset**, which is then used in line 15, where the loss function $L_{cegis}$ is evaluated on these counterexamples. We will revise the notation for the loss functions $L_{cegis}(\\theta)$ as $L_{cegis}(\\theta;x)$ and uses $L_{cegis}(\\theta;\\text{Dataset})$ in Algorithm 1 line 15.\n\n### **W2 & Q1: Clarifications on α,β-CROWN**\n\nWe would like to clarify that we never mentioned that we used Interval Bound Propagation (IBP) in our paper. Instead, we used tighter linear bound propagation. We did not provide a detailed discussion of α,β-CROWN in our paper as it is a well-developed tool, and not our main contribution. α,β-CROWN is a GPU-accelerated neural network verifier capable of handling complex computation graphs with general nonlinearities. Intuitively, CROWN can be seen as a generalization of backpropagation. Just like how gradients are propagated through a network, CROWN propagates **linear upper/lower bounds** across the entire computation graph over an input region. Instead of computing derivatives, it computes **linear relaxations** for each nonlinear operation. These local relaxations are then composed layer by layer (analogous to how backpropagation works), resulting in global affine bounds of output with respect to inputs. This approach, known as **linear bound propagation**, is **much tighter than IBP**, which only tracks independent input-output intervals and ignores inter-variable dependencies. Moreover, CROWN can do branch-and-bound to further tighten the bounds. By applying CROWN to our computation graph $\\nabla V(x)f(x)$, we can provide formal verification in a more scalable manner compared to traditional SMT solvers like dReal.\n\nWhile the linear bound propagations have been presented in α,β-CROWN’s original paper [1] and extended to general computation graphs [2], we understand that it can be challenging to grasp these algorithms quickly, so we will **provide a toy example** on how it applies to our verification setting in the appendix.\n\nAlthough not our primary focus, our work does make new technical extensions to α,β-CROWN, including enabling support on computation graphs with Jacobian. Obtaining the **tightest possible** relaxations for nonlinearities with Jacobian operators is also challenging - the differentiation of many operators, such as tanh or sigmoid, can introduce some tricky nonlinearities like $1 - \\tanh^2(x)$. In principle, α,β-CROWN could handle this by treating it as a combination of operators, like treating $\\tanh^2(x)$ as the combination of $\\tanh(x)$ and $x^2$, and relaxing it using the combination of linear bounds of these more elementary nonlinearities. However, this bound will be loose. In our work, we implement a **tight relaxation** of the function $h(x) = 1 - \\tanh^2(x)$ (see Appendix B.2) directly to significantly tighten the bound of this operator (**over 6 times tighter**), and apply the same idea to sigmoid’s derivative. We discussed these in Appendix B.2.\n\n### **W3: Guarantees near origin**\n\nThe paper mentioned by the reviewer addresses only **polynomial systems and activations**, whereas our approach needs to handle **general nonlinearities**. In our setting, achieving exact guarantees is impossible due to verifiers’ numerical issues near the origin when handling general nonlinearities. To our knowledge, no prior work can achieve hole-free verification in this setting. Therefore, we instead provide alternative **theoretical guarantees** ensuring the excluded region is **forward-invariant**. In comparison, many previous works do not provide any guarantee at all, as discussed in lines 228-230. Moreover, the actual region excluded is **tiny** and thus our approach remains **practically effective**. For example, $\\lbrace V(x) \\le 0.0068\\rbrace$ is removed for the Cartpole system, which occupies **only 0.003% of the ROA**. Appendix C.4 shows that our hole excluded is **negligible** for all tested systems.\n\n### **W4 & Q2: High-dimensional Verification**\n\nWe will revise our paper to be clear about which systems are verified under which schemes. Moreover, we will change our terminology used in section 3.2. Namely, we will change \"Verification Schemes\" into \"Evaluation Schemes\", and change \"Empirical PGD Verification\" and \"Monte-Carlo Trajectory Verification\" to \"PGD Evaluation\" and \"Trajectory Evaluation\" to avoid misunderstandings. We note that although we could only verify forward invariance for 6D systems, this is **already the SOTA result for continuous-time systems**. Moreover, we achieve **significantly larger ROA than all the baselines**, which inherently makes the verification problem harder. Below, we list the evaluation schemes we apply for each system and the ROA volume improvement compared to the previous SOTA.\n\n|Dimensions|Systems|Traj.|PGD.|Formal Verification|ROA (compare to SOTA)|\n|:-:|:-:|:-:|:-:|:-:|:-:|\n|2D|Van Der Pol|✓|✓|Full|2.85x|\n||Double Int|✓|✓|Full|2.32x|\n||Pendulum B|✓|✓|Full|6.05x|\n||Pendulum S|✓|✓|Full|3.82x|\n||Path Tracking B|✓|✓|Full|5.07x|\n||Path Tracking S|✓|✓|Full|4.97x|\n|4D|Cartpole|✓|✓|Full|330.3x|\n|6D|2D Quadrotor|✓|✓|Forward Invariance ($\\lbrace 0.2 \\le V \\le 0.21\\rbrace$)|150000x|\n||PVTOL|✓|✓|-|383x|\n||Ducted Fan|✓|✓|-|N/A|\n|12D|3D Quadrotor|✓|-|-|N/A|\n\nThe reviewer mentions that some SMT baselines, such as NLC and Fossil, could scale on par with our method. We respectfully disagree and believe our method is much better than SMT baselines. We list our reasons below.\n\n- NLC and Fossil [3] have **soundness issues**. They mistakenly always use a bias-free network for controllers, leaving $u(0) = 0$. However, this is not true for Path Tracking and all >6D systems we considered. Further discussions can be found in Appendix C.1.\n- NLC and Fossil **do not scale to Cartpole system** that we excel in. Moreover, we have compared our results to Fossil on the systems that do not have soundness issues, which have been shown to be much better than NLC [3]. Our method achieves greater scalability and much larger ROA volumes compared to Fossil.  \n\n|Systems|Fossil 2.0||Ours||\n|:-:|:-:|:-:|:-:|:-:|\n||ROA|Succ|ROA|Succ|\n|Van Der Pol|3.05±0.58|60%|**57.6±3.4**|100%|\n| Double Int|258.1±36.3|80%|**302.5±10.7**|100%|\n| Pendulum S|106.5 ± 23.5|80%|**1169.2±124.5**|100%|\n| Pendulum B|285.2±53.3|80%|**2946.5±149.1**|100%|\n| Cartpole|—|0%|**306.1±54.2**|100%|\n\n### **References**\n[1] Zhang, Huan, et al. \"Efficient neural network robustness certification with general activation functions.\" NeurIPS 31 (2018).  \n[2] Xu, Kaidi, et al. \"Automatic perturbation analysis for scalable certified robustness and beyond.\" NeurIPS 33 (2020).  \n[3] Edwards, Alec, et al. \"Fossil 2.0: Formal certificate synthesis for the verification and control of dynamical models.\"​​HSCC. 2024."}, {"review_id": "ewQ5J2jASK", "reviewer": "Reviewer_QzKe", "summary": "This paper presents a novel approach for learning stabilizing neural controllers and corresponding Lyapunov functions, grounded in Zubov’s partial differential equation (PDE) framework. The motivation for using Zubov’s theorem lies in its natural connection to estimating the region of attraction (ROA). The authors introduce a Zubov-inspired sampling strategy that combines points from both inside and outside the current ROA estimate, allowing for more informed learning. This is followed by a counterexample-guided inductive synthesis (CEGIS) loop that iteratively refines the candidate Lyapunov function until a verifiable certificate is obtained. The framework employs the \nalpha-beta-CROWN verifier for neural network verification and is evaluated across a diverse set of benchmarks, from simple 2D systems like the inverted pendulum to a 12-dimensional 3D quadrotor. Experimental results demonstrate that the proposed approach produces significantly larger ROA estimates compared to baselines and achieves superior scalability in neural network verification relative to dReal.", "questions": "The FOSSIL framework appears to be a state-of-the-art tool for synthesizing neural Lyapunov certificates. How does your approach compare with FOSSIL in terms of ROA estimation, verification strength, and scalability? Could you clarify the reason for not including a direct comparison?", "limitations": "As noted earlier, the paper does not include comparisons with closely related tools like FOSSIL.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n\n1. The paper studies a timely and well-motivated problem. It is well-written and technically sound, with a non-trivial and original approach.\n2. The experimental benchmarks span a wide range of complexity and dimensionality, clearly demonstrating the benefits of the proposed approach relative to existing baselines.\n3. The authors implement several practical enhancements to alpha-beta CROWN, including linear relaxations of Jacobian operators and a branch-and-bound heuristic, improving the scalability and efficiency of verification.\n\nWeaknesses\n1. The introduction would benefit from a clearer articulation of the novelty of leveraging Zubov’s theorem in contrast to more traditional Lyapunov-based approaches.\n2. The paper omits several relevant citations related to CEGIS-based synthesis of Lyapunov functions, including:\n- Formal Synthesis of Lyapunov Neural Networks, Abate et al.\n- Counterexample-Guided Inductive Synthesis of Control Lyapunov Functions for Uncertain Systems, Masti et al., IEEE\n- Counterexample-Guided Synthesis of Control Lyapunov Functions for Switched Systems, Ravanbakhsh and Sankaranarayanan, IEEE\n- FOSSIL: A Software Tool for the Formal Synthesis of Lyapunov Functions and Barrier Certificates Using Neural Networks, Abate et al., ACM", "was_revised": true, "final_justification": "The authors have adequately addressed my concerns. I will maintain my original positive score.", "rebuttal": "We thank the reviewer for recognizing our contribution. We here provide a point-to-point response to the reviewer's question.\n\n### **Weakness 1: More introduction on Zubov approach**\nThank you for your suggestion, we will add a more detailed explanation of Zubov theorem in the introduction to clearly distinguish our work from the works that only utilize traditional Lyapunov theory.\n\n### **Weakness 2 & Question 1: Fossil baseline and discussion**\nThank you for pointing out the missing citations and baselines. We will add the citations of the mentioned papers, and we have since compared our method against Fossil 2.0 and present the results here.\n\nWe notice that Fossil has the same issue as NLC mentioned in Appendix C.1, as it **falsely assumes that $u(0) = 0$ on the equilibrium**. Thus, we here only present their results on the systems that satisfy the property $u(0) = 0$. It can be noticed that it **does not scale to the 4D** cartpole system, and **achieves smaller ROA volumes** than our method in 2D systems. \n\n|Systems|Fossil 2.0||Ours||\n|:-:|:-:|:-:|:-:|:-:|\n||ROA|Succ|ROA|Succ|\n|Van Der Pol|3.05±0.58|60%|**57.6±3.4**|100%|\n| Double Int|258.1±36.3|80%|**302.5±10.7**|100%|\n| Pendulum S|106.5 ± 23.5|80%|**1169.2±124.5**|100%|\n| Pendulum B|285.2±53.3|80%|**2946.5±149.1**|100%|\n| Cartpole|—|0%|**306.1±54.2**|100%|\n\nRegarding scalability and verification strength, Fossil is more similar to NLC [1] and uses SMT solvers like dReal or Z3 for verification and counterexample generation. Since we are using the more scalable and powerful verifier α,β-CROWN and PGD-based quick counterexample sampling, **our method is better in both verification and training scalability**, which is also proved by the experimental results.\n\n***Side note about Fossil.***\nWe believe there are some subtle issues in their verification formulation. In reference [2], which underpins Fossil 2.0, we notice that their Certificate 2 **misses the boundary condition**, which should supposedly ensure either $\\lbrace x: V(x) \\leq \\beta \\rbrace \\cap \\mathcal{X}$ is empty or any trajectory starting from this intersection will remain in $\\mathcal{X}$. We have observed that their method can occasionally return false positives even though $\\mathcal{X}$ is already chosen to be much bigger than $\\mathcal{X}_I$, i.e., a controller and Lyapunov function are claimed to be found while the controller cannot stabilize any trajectory. Therefore, on top of their returned results, we ran an additional check with trajectory simulations, and recorded a negative for a run if the controller does not work. \n\n### **References**\n[1] Chang, Ya-Chien, Nima Roohi, and Sicun Gao. \"Neural lyapunov control.\" _Advances in neural information processing systems_ 32 (2019).  \n[2] Edwards, Alec, Andrea Peruffo, and Alessandro Abate. \"A general framework for verification and control of dynamical models via certificate synthesis.\" _arXiv preprint arXiv:2309.06090_ (2023)."}, {"review_id": "41QlTxfzSn", "reviewer": "Reviewer_N7NM", "summary": "This paper presents a novel two-stage learning framework for synthesizing stabilizing neural controllers with formal stability guarantees. The authors propose a Zubov-inspired sampling strategy and iterative domain expansion mechanism to significantly reduce conservatism in training neural Lyapunov functions, while extending the α,β-CROWN verifier to handle Jacobian-related operators for efficient formal verification. The method demonstrates remarkable improvements, achieving 5-1.5×10^5 times larger regions of attraction compared to baselines and 40-10000× faster verification than SMT solvers across various nonlinear systems from 2D to 12D. The work makes important contributions in bridging the gap between empirical neural control performance and formal stability guarantees through its innovative training and verification pipeline.", "questions": "How does the method's computational complexity scale with system dimensionality, particularly for the verification stage which still faces challenges in higher dimensions?\n\nCould you elaborate on the robustness of the approach when dealing with model uncertainties or disturbances not accounted for during training?\n\nWhat modifications would be needed to apply this framework to hybrid or discrete-time systems commonly encountered in robotics applications?\n\nHave you considered combining the approach with online adaptation mechanisms to maintain stability guarantees under slowly varying system dynamics?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "While the methodology shows impressive empirical results across multiple nonlinear systems (achieving ROA volumes 5-1.5·10^5 times larger than baselines), its formal verification capabilities remain limited for higher-dimensional systems due to computational constraints in Jacobian-bound propagation. The theoretical contributions are well-grounded with proper assumptions and proofs, though the practical implementation requires careful hyperparameter tuning that may challenge real-world adoption.", "was_revised": true, "final_justification": "The proposed pipeline achieves impressive verified regions of attraction, yet its reliance on careful hyper-parameter tuning and its limited treatment of uncertainty, hybrid dynamics, and time-varying systems leave notable gaps.", "rebuttal": "We thank the reviewer for acknowledging our contribution and for highlighting several promising directions for future work. Below, we will have a point-to-point response to each of the questions raised.\n### **Weakness 1: Verification in High Dimension**\nWe share the reviewer’s concern that scaling verification to high-dimensional systems remains a key bottleneck, even though our method currently sets the state-of-the-art. We want to mention that we are actively pushing the frontiers of this research direction and have demonstrated results that are not possible with existing baselines. We would love to motivate more people to study this problem and further improve the verification pipeline for better real-world applications.\n\n### **Weakness 2: Hyperparameter Tuning Difficulty**\nWe thank the reviewer for bringing this up. We do agree that learning a stabilizing controller is quite a challenging setting, and our method is working towards the SOTA on making the tuning easier and training more stable. We have simplified the two most difficult hyperparameters, training domains and initializations, to tune in baselines such as [2,3], thanks to our novel training pipeline, and we achieved a much better success rate in Table 1.\n\n* One difficult hyperparameter is the **training domain**; selecting it not optimally may cause complete failure of the training. However, we have made this **hyperparameter selection mostly automatic** with our novel dynamic domain expansion. As Table 4 in Appendix C.3 shows, we begin with a loosely chosen initial domain and rely on dynamic expansion to refine it on the fly, yet the method still achieves very strong performance. \n* Another difficult hyperparameter is the **initialization**. Previous work depends on RL/LQR initializations can fail easily with a badly chosen initialization parameter. However, using our novel two-stage training framework, we eliminate the need for special initializations and achieve strong performance by training from scratch.\n\n### **Question 1: Verification Complexity**\nUnfortunately, there is no simple or uniform answer. In practice, verification complexity can **depend on many factors**, such as the system dimensions, the range of the verification domain, and the model architecture. For example, alpha,beta-CROWN was originally designed for robustness verification, where a typical problem is handling a perturbed image input through a CNN. In this setting, the input dimension is huge (e.g., 64 * 64 * 3); however, the verifier is still able to handle it due to the allowed perturbation on each dimension being small (around 0.01), and the computation graph is relatively simple. In our setting, each dimension of the **input domain can be huge** (around 10), and the **computation graph has a very complicated topology** due to the complexity of the dynamical system and the Lyapunov condition. Moreover, the computation graph involves **many unique nonlinear operators** such as tanh, sigmoid, their derivatives, and multiplication with both operands perturbed, which are hard to bound tightly. This issue presents a **common challenge** for all existing verification approaches to this problem. Despite this, **our training and verification method currently achieves state-of-the-art performance** for this class of problems, and we believe it advances the boundary of what is possible in high-dimensional nonlinear system verification. \n\n### **Question 2: Robustness under uncertainties/disturbances**\nThis is a very interesting direction that well deserves further exploration, although it falls outside the scope of our current paper. We believe verifying the trained model’s robustness under uncertainties or disturbances is **non-trivial**, as it would require a novel verification formulation of robust control theories like dissipativity theory. It is also possible that we need to reformulate our training framework under the robust control framework, which can be more non-trivial since it requires the combination of Zubov theory and dissipativity theory. \n\n### **Question 3: Extension to discrete/hybrid systems**\nOur current training framework’s success depends on the existence of a Zubov-type theory that can characterize the underlying ROA. For discrete-time systems, we notice that there is some recent work that discusses the discrete-time Zubov theorem [1], and our work should be readily **extendable to discrete-time settings**. However, there are still some **non-trivial challenges** in the application of dynamic domain expansion due to some connectedness issues. It is also highly non-trivial to extend our approach to hybrid systems, since a Zubov-type theory can be **more difficult to obtain** in this setting, and would require novel combinations of Zubov theory, hybrid system theory, and nonsmooth analysis. Both of these settings are beyond the scope of our current paper, but we believe it is an interesting and promising future direction. \n\n### **Question 4: Extension to slowly time-varying systems**\nCombining our approach with online adaptation mechanisms to maintain stability guarantees under slowly varying system dynamics is also not trivial. With prescribed bounds on the slow-changing parameter and its changing rate, we believe it is possible to address the slowly varying parameter in the robust control framework and combine our approach with dissipativity theory to train a robust controller. However, as commented above, combining our approach with dissipativity theory is highly non-trivial in the first place, and hence, **more future work will be needed** to pursue this interesting direction. \n\n### **References**\n[1] Serry, Mohamed, et al. \"Safe Domains of Attraction for Discrete-Time Nonlinear Systems: Characterization and Verifiable Neural Network Estimation.\" _arXiv preprint arXiv:2506.13961_ (2025).  \n[2] Yang, Lujie, et al. \"Lyapunov-stable neural control for state and output feedback: A novel formulation.\" arXiv preprint arXiv:2404.07956 (2024).  \n[3] Shi, Zhouxing, Cho-Jui Hsieh, and Huan Zhang. \"Certified training with branch-and-bound: A case study on lyapunov-stable neural control.\" arXiv preprint arXiv:2411.18235 (2024)."}, {"review_id": "FOnY1eBCj4", "reviewer": "Reviewer_NRaM", "summary": "This paper proposes a new two-stage Zubov-based training framework for synthesizing stabilizing neural network controllers, it contains a curriculum training domain expansion and training sample selection that significantly reduces the conservatism during training. The author also improves the advanced neural network verification tool and propose a novel verification algorithm. The empirical results demonstrate large improvements in ROA estimation and achieves significant verification speedups.", "questions": "- The paper overall looks good to me, but I have a general question reguarding to this paper -- why this paper is send to NeurIps rather than ICCAD or DAC? What's it contribution reguarding to ML side?", "limitations": "Yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength\n\n- The paper is logically organized and clearly written. The problem setup/motivation/methods are introduced with appropriate theoretical grounding.\n\n- The approach shows strong empirical scalability. And the speedups in verification and the increased ROA volumes have practical significance.\n\nWeakness\n\n- The lack of formal guarantees in higher-dimensional systems limits the strength of the paper's claims in safety-critical applications. The dependence on a verifier that still struggles with high-dimensional tight bounds means that many results are still empirical.\n\n- This is more likely a control-theoretic contribution, its broader impact within ML remain to be discussed.", "was_revised": true, "final_justification": "I have read the author's response and it has resolved all my mentioned concerns. So I keep the score as borderline accepted.", "rebuttal": "We appreciate the reviewer’s positive feedback. Below, we address each question the reviewer raised for our work.\n\n### **Weakness 1: Verification in High Dimensions**\nWe agree with the reviewer that scaling verification to high-dimensional systems remains a key bottleneck, even though our approach achieves state-of-the-art performance. We want to mention that we are actively pushing the frontiers of this research direction and have demonstrated results that are unattainable with existing baselines. We would love to motivate more people to study this problem and further improve the verification pipeline for better real-world applications.\n\n### **Weakness 2 & Question 1: Why NeurIPS and Significance to the ML community**\nAlthough ICCAD and DAC are great venues with papers on relevant topics, we believe our work presents its unique value to the machine learning community due to the following three reasons:\n\n* Our work lies at the **intersection of machine learning and control theory**, addressing the important problem of learning-based control. It demonstrates how machine learning techniques can be applied to a fundamental challenge in engineering and science. **Similar lines of work have been published in leading ML venues.**. For example, [1] is published at NeurIPS 2019, [2] at NeurIPS 2023, and [3] at ICML 2024. \n* More broadly, our work contributes to the problem of **Lyapunov function synthesis**, which has recently received increasing attention in the ML community, with related papers such as [4] at ICML 2025 and [5] at CVPR 2025. We believe our work shares similar merit and is suitable for publishing at a machine learning conference. \n* Finally, our work extends the **neural network verification** toolbox α,β-CROWN ([6-11, NeurIPS 2018, 2020, 2021, 2022, 2023, 2024]) in some nontrivial ways, such as **enabling automatic bound computation over Jacobian operators and providing new tight relaxations** for many commonly used derivatives. \n\nWe believe that publishing this work at a comprehensive machine learning venue like NeurIPS will engage both the neural network verification and control communities, and help make progress on the broader challenge of scalable, verifiable learning based control.\n\n### **References**\n[1] Chang, Ya-Chien, Nima Roohi, and Sicun Gao. \"Neural Lyapunov control.\" _Advances in neural information processing systems_ 32 (2019).  \n[2] Wu, Junlin, et al. \"Neural Lyapunov control for discrete-time systems.\" _Advances in neural information processing systems_ 36 (2023): 2939-2955.  \n[3] Yang, Lujie, et al. \"Lyapunov-stable neural control for state and output feedback: A novel formulation.\" _arXiv preprint arXiv:2404.07956_ (2024).  \n[4] Zou, Haohan, et al. \"Analytical Lyapunov Function Discovery: An RL-based Generative Approach.\" _arXiv preprint arXiv:2502.02014_ (2025).  \n[5] Zhao, Hanrui, et al. \"Learning-enabled Polynomial Lyapunov Function Synthesis via High-Accuracy Counterexample-Guided Framework.\" _Proceedings of the Computer Vision and Pattern Recognition Conference_. 2025.  \n[6] Zhang, Huan, et al. \"Efficient neural network robustness certification with general activation functions.\" Advances in neural information processing systems 31 (2018).  \n[7] Xu, Kaidi, et al. \"Automatic perturbation analysis for scalable certified robustness and beyond.\" Advances in Neural Information Processing Systems 33 (2020): 1129-1141.  \n[8] Wang, Shiqi, et al. \"Beta-crown: Efficient bound propagation with per-neuron split constraints for neural network robustness verification.\" Advances in neural information processing systems 34 (2021): 29909-29921.  \n[9] Zhang, Huan, et al. \"General cutting planes for bound-propagation-based neural network verification.\" Advances in neural information processing systems 35 (2022): 1656-1670.  \n[10] Kotha, Suhas, et al. \"Provably bounding neural network preimages.\" Advances in Neural Information Processing Systems 36 (2023): 80270-80290.  \n[11] Zhou, Duo, et al. \"Scalable neural network verification with branch-and-bound inferred cutting planes.\" Advances in Neural Information Processing Systems 37 (2024): 29324-29353."}]

Author rebuttal:

Meta-review: The paper proposes a method that jointly learns a neural control policy and a neural Lyapunov function, using a novel Zubov-inspired loss to explicitly maximize the region of attraction (ROA). Experiments show that this approach achieves larger ROAs than prior methods while ensuring formal correctness. Authors present several practical enhancements to alpha-beta CROWN, including linear relaxations of Jacobian operators and a branch-and-bound heuristic, improving the scalability and efficiency of verification. The experimental benchmarks span a wide range of complexity and dimensionality, clearly demonstrating the benefits of the proposed approach relative to existing baselines. One weakness at the moment is scaling to higher dimensional systems. Overall, the paper is a strong contribution. It is strongly suggested that the authors consider the suggestions for improving the presented by the reviewers.
Decision: Accept (spotlight)