{'fraction_non_alphanumeric': 0.04272444519166106, 'fraction_numerical': 0.01687542030934768, 'mean_word_length': 4.104043762737317, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kernel methods for center manifold approximation and a data-based version of the Center Manifold Theorem

For dynamical systems with a non hyperbolic equilibrium, it is possible to significantly simplify the study of stability by means of the center manifold theory. This theory allows to isolate the complicated asymptotic behavior of the system close to the equilibrium point and to obtain meaningful predictions of its behavior by analyzing a reduced order system on the so-called center manifold. Since the center manifold is usually not known, good approximation methods are important as the center manifold theorem states that the stability properties of the origin of the reduced order system are the same as those of the origin of the full order system. In this work, we establish a data-based version of the center manifold theorem that works by considering an approximation in place of an exact manifold. Also the error between the approximated and the original reduced dynamics are quantified. We then use an apposite data-based kernel method to construct a suitable approximation of the manifold close to the equilibrium, which is compatible with our general error theory. The data are collected by repeated numerical simulation of the full system by means of a high-accuracy solver, which generates sets of discrete trajectories that are then used as a training set. The method is tested on different examples which show promising performance and good accuracy.


Introduction
Center manifold theory plays an important role in the study of the stability of dynamical systems when the equilibrium point is not hyperbolic. It isolates the complicated asymptotic behavior by locating the center manifold which is an invariant manifold tangent to the subspace spanned by the eigenspace of eigenvalues on the imaginary axis. Then, the dynamics of the original system will be essentially determined by the restriction of this dynamics on the center manifold since the local dynamic behavior "transverse" to this invariant manifold is relatively simple as it corresponds to the flows in the local stable (and unstable) manifolds. In practice, one does not compute the center manifold and its dynamics exactly since this requires the resolution of a quasilinear partial differential equation which is not easily solvable. In most cases of interest, an approximation of degree two or three of the solution is sufficient. Then, the reduced dynamics on the center manifold can be determined, its stability can be studied and then conclusions about the stability of the original system can be obtained. For a detailed treatment of this topic, we refer e.g. to [8,14,15,16,22,28].
In this paper we consider data-based approximations of the center manifold, and our main result proves in great generality that any sufficiently accurate approximation can be used to draw conclusions on the asymptotic stability of the system, thus establishing one direction of a data-based center manifold theorem.
We then turn our attention to concrete procedures to construct approximants that provide that kind of error control, and to this end we revise and refine a kernel-based approximation method that we introduced in [12]. This approximation can be constructed based solely on the knowledge of points on the trajectories of the system close to the equilibrium point, that are obtained in practice by the numerical integration of the system. This method is formulated in the Reproducing Kernel Hilbert Spaces (RKHS) of a strictly positive kernel (see [2]), which have provided strong mathematical foundations for analyzing dynamical systems in other settings (see e.g. [1,3,4,5,6,7,9,11,12,13,18,19]). Our method explicitly incorporates into the approximant some prior knowledge on the structure of the manifold. This is possible using some RKHS theory, and we prove that this constrained problem admits a unique solution. In particular, we show how to control the error of this approximation, and the resulting bound guarantees that this kernel-based approximation is compatible with the general approximation discussed above. Moreover, these bounds are such that they can be applied even if the approximant is constructed using a subset of the dataset selected by means of certain greedy methods.
The paper is organized as follows. In Section 2 we recall the definition of the center manifold and its role in the study of the stability of certain equilibria. In Section 3, we prove a data-based version of the center manifold theorem. In section 4 we describe kernel methods to construct a data-based approximation of the center manifold. Some numerical experiments are given in Section 5.

The center manifold
We consider an open set D ⊂ R n with 0 ∈ D, and a dynamical systeṁ where f : D → R n is a continuously differentiable function and x := x(t) ∈ D for all t. We are interested in the study of the behavior of this possibly high-dimensional system around an equilibrium point, which we may assume w.l.o.g. to be the origin, i.e., f (0) = 0.
Here and in the following of this paper, we use the notation ∇ x g(x) := (∂ 1 g(x), . . . , ∂ n g(x)) T and D x f (x) := (∂ i f j (x)) 1≤i≤n,1≤j≤m to denote the gradient of a function g : R n → R and the Jacobian matrix of a vector-valued function f : R n → R m . Moreover, a norm · without subscript refers to the Euclidean 2-norm on R n for some n, or to the induced norm of a matrix in R n×n .
Letting L := D x f (x) |x=0 ∈ R n×n , we can rewrite (1) aṡ with the nonlinear component N(x) := f (x) − Lx, and denote as σ R (L) the set of real parts of the eigenvalues of L. A classical result relates the stability of the equilibrium with the spectrum of L, and in particular it is known that if all the eigenvalues of L have negative real part, i.e., σ R (L) ⊂ R <0 , then the origin is asymptotically stable, and if L has some eigenvalues with positive real parts, then the origin is unstable. If instead σ R (L) ⊂ R ≤0 then the linearization fails to determine the stability properties of the origin, and thus the analysis of this situation requires to employ additional tools.
In this case, we can first use a linear change of coordinates to separate the eigenvalues of L with zero and with negative real part, i.e., we can rewrite (1) as 1 Intuitively, we expect the stability of the equilibrium to only depend on the nonlinear term N 1 (x, y). This intuition turns out to be correct, and indeed it can be properly formalized by means of the Center Manifold Theorem.
Before stating this theorem, we start by recalling a sufficient condition for the existence of a center manifold. Observe that the following assumption on the smoothness of N 1 , N 2 require that f ∈ C 2 (D).
Theorem 2.1 (Theorem 8.1 in [17]). Assume that N 1 and N 2 are twice continuously differentiable on D and are such that Furthermore, assume that the eigenvalues of L 1 have zero real part and the eigenvalues of L 2 have negative real part. Then there exists a neighborhood B ⊂ R d of the origin and a center manifold h : B → R m for (2), i.e., y = h(x) is an invariant manifold 2 for (2), h is continuously differentiable, and Under the assumptions of this theorem, using (2) we can deduce that h satisfies the PDE and the following center manifold theorem ensures that there are continuously differentiable solutions to this PDE. Moreover, the theorem also allows to deduce the stability of the origin of the full order system (2) from the stability of the origin of a reduced order system called the center dynamics. . Under the assumptions of Theorem 2.1, the equilibrium x = 0, y = 0 of the full order system (2) is asymptotically stable (resp. unstable) if and only if the equilibrium x = 0 of the reduced dynamic (dynamics on the center manifold)ẋ is asymptotically stable (resp. unstable).
1 Please note that x will be a generic variable related to the state of a dynamical system and thus may have varying dimension and meaning throughout this paper. This result guarantees that, upon solving the PDE (5) for the unknown h : R d → R m , the problem of analyzing the stability properties of the system (2) reduces to analyzing the nonlinear stability of the lower dimensional system (6). This second problem is of smaller dimension and thus, provided the knowledge of h, the approach is attractive to obtain information on the system (1) via a reduced order model.
Observe that the dimension d of this reduced system is equal to the number of eigenvalues with zero real part (counted with multiplicity) of the original system (1), and thus it may take any value between 1 and the full dimension n, depending on the particular dynamical system under consideration. The actual dimensional reduction provided by the center manifold is thus highly varying depending on the problem.
We remark that this methodology is valid also for parameterized dynamical systems and is used to study the stability of dynamical systems with bifurcations.
Remark 2.1. An exact knowledge of the center manifold h is not required for this equilibrium analysis, since it is sufficient to have an approximated solution of the PDE (5). Indeed, it is frequently sufficient to compute only the low degree terms of the Taylor series expansion of h around x = 0, i.e., if (·) [k] is the degree k part of the Taylor series of h, the approximation h(x) ≈ h [1] x + h [2] (x) + h [3] is sufficient to obtain an approximation of the dynamics of order ε d as x ≤ ε.
Nevertheless, even this approximated knowledge of h can be difficult to obtain in practice for a general ODE.

A data-based version of the Center Manifold Theorem
In this section we show that one direction of Theorem 2.2 holds even when the center manifold is replaced by a sufficiently accurate approximation. Namely, we assume to have access only to an approximantĥ : R d → R m of the center manifold, instead of h itself, and we consider the approximated reduced dynamicsẋ This system is now both low dimensional and completely explicit, sinceĥ is known. Observe in particular that the change of coordinates required to obtain L 1 and N 1 from L and N is based solely on a diagonalization of L.
Depending on different conditions that the approximant may satisfy, we show in Theorem 3.2 that the asymptotic stability of the origin of the full order dynamics follows from the asymptotic stability of the approximated center dynamics. To formulate this data-based version of the Center Manifold Theorem, we first need to recall the following result. Then there exists a neighborhood B ⊂ D of zero and a Lyapunov function for the system, i.e., a continuously differentiable function V : We can now state and prove our main result. Actually, Theorem 3.2. can be understood in a more general way as a "robust" center manifold theorem asĥ does not necessarily need to be a data-based approximant. But as we only use such data-based approximants in the following, this terminology data-based Center Manifold Theorem seems appropriate.
Assume that the point x = 0 is an asymptotically stable equilibrium of the approximated center dynamicsẋ = L 1 x + N 1 (x,ĥ(x)).
Let V be the function obtained by applying Theorem 3.1 to (9), let α(·) and k be as in (8), and assume that there exists p > 0 such that Under these assumptions, define Then, if there exist q ≥ p and ε ≤ ε 0 such that the equilibrium x = 0, y = 0 of the full order dynamics (2) is asymptotically stable.
Proof. First we prove that if x = 0 is an equilibrium of (9) then (x, y) = (0, 0) is an equilibrium of the full order system. Indeed, since h(0) = 0 by definition, we have that (12) implies that h(0) = 0 whatever is the value of q > 0. Thus, since x = 0 is an equilibrium of (9), we have and thus x = 0 is an equilibrium of the reduced system (6). By Theorem 2.2 it follows that (0, 0) is an equilibrium of (2). In order to prove asymptotic stability we adapt the proof of Theorem 8.2 in [17]. We consider the change of variables which transforms (2) intȯ and this transformation is such that this new system has the same stability properties of (2) (see [17]). Observe that the difference y(t) − h(x(t)) represents the deviation of the trajectory from the center manifold at any time t. This may be identically zero if y(0) = h(x(0)) since in this case the solution (x(t), y(t)) will lie in the manifold for all t ≥ 0, and if so the second equation in (13) is trivially true.
Defining now it is proven in [17] that by manipulation of the previous equations it follows thaṫ and furthermore that for any radius ρ > 0 there exist constants k i := k i (ρ) > 0 such that for all (x, w) in the ball B ρ := B ρ (0, 0) with center (0, 0) and radius ρ it holds where the constants k 1 and k 2 can be made arbitrarily small by choosing ρ small enough. We now consider the reduced order system (9) and show that if the origin x = 0 is asymptotically stable, then the origin (x, w) = (0, 0) is asymptotically stable for the full system in the form (13).
Since x = 0 is asymptotically stable, by Theorem 3.1 there are a continuously differentiable and positive function V and a constant k > 0 such that in a neighborhood B of the origin is a strictly increasing continuous function with α(0) = 0. Furthermore, by choosing ρ small enough we may assume that x ∈ B for all (x, w) ∈ B ρ , and thus both the bounds (15) and (16) hold in B ρ . Moreover, in the proof of Theorem 8.2 in [17] it is shown that there is a matrix P which is positive definite and solves the Lyapunov equation L T 2 P + P L 2 = −I, where I is the identity matrix. We denote as λ min , λ max > 0 the minimal and maximal eigenvalues of this matrix.
Using V and P , we define a candidate Lyapunov function for the full order system. By construction this function is zero in (x, w) = (0, 0) and strictly positive in B ρ \ {(0, 0)}. If we show that its derivative along the trajectories of the full order system is strictly negative, it follows (see Theorem 4.1 in [17]) that the origin (x, w) = (0, 0) is asymptotically stable for the full system. To see this, we compute the derivative of V along the trajectories of the full order system asV We take from the proof of Theorem 8.2 in [17] the estimate of the second term, i.e., To bound instead the first term of (18), we add and subtract N 1 (x,ĥ(x)), we use the two inequalities (16) and the bound (15) with i = 1, and via the Cauchy-Schwartz inequality we obtain Moreover, since N 1 is continuously differentiable by assumption, for all x the function w → N 1 (x, w) is Lipschitz continuous for (x, w) ∈ B ρ , with a Lipschitz constant that is bounded by Inserting (19), (20) and (21) into (18) gives for all (x, w) ∈ B ρ the bounḋ The term − kk 1 is always strictly positive if k 1 , k 2 are made sufficiently small by choosing an appropriately small value of ρ. Thus, we are left with the inequalitẏ and using the error bound (12) we get thaṫ If x = 0 and w = 0, the last term ensures thatV(x, w) < 0. Moreover, the term i.e., when ε ≤ ε 0 , which is true by assumption. It follows thatV(x, w) < 0 for all (x, w) ∈ B ρ \ {(0, 0)}, and thus the origin (0, 0) is asymptotically stable for the full order system (13), and we are done.
We conclude this section with some remarks on the theorem and on its proof. (12) is sharp in q, since the assumption that the approximation error behaves as α when x → 0, i.e., that (12) holds with q ≥ p, is crucial for the proof. Indeed, if only an error of order x q with q < p can be obtained in (12), then (11) would give ε 0 = 0. Thus it would not be possible to conclude the proof by showing thatV(x, w) < 0 in a neighborhood of (0, 0), no matter how small we choose the neighborhood, unlessĥ = h.

Remark 3.1 (Weaker error conditions -Error rate). Observe that the error statement
This is the case also if one is trying to prove stability instead of asymptotic stability of the origin of the full order system. Indeed, in this case it would be sufficient to prove thaṫ , equality is possible), but also in this case the same problem as before would persist. (12) is instead not necessarily optimal with respect to the constant ε. Indeed, in (11) we defined for simplicity ε 0 as a function of ℓ D , which is the Lipschitz constant of N 2 in D with respect to its second variable. Inspecting the proof, it is clear that it would instead be sufficient to define ε 0 in terms of ℓ ρ (see (21)). Since ℓ ρ ≤ ℓ D , the new value of ε 0 would be larger, and thus the condition ε ≤ ε 0 could be made less restrictive. Remark 3.3 (Unstable equilibrium). The Center Manifold Theorem 2.2 guarantees also that if the origin is unstable in the reduced system, then (0, 0) is unstable in the full system. We do not see how to extend this result to our case.

Remark 3.4 (Data based Center Manifold Theorem -Part II).
We currently have no proof of the other implication of the Center Manifold Theorem, i.e., proving that if the origin of the full order dynamics is asymptotically stable, then also the origin of the approximated reduced dynamics is asymptotically stable.
This converse implication would be of practical interest in order to guarantee the wellposedness of the problem of studying the asymptotic stability of an equilibrium of the reduced and approximated system. Indeed, this result would ensure that if one starts from a full system with an asymptotically stable equilibrium and constructs a sufficiently accurate approximant, then also this approximant has an asymptotically stable equilibrium.
Remark 3.5 (Finding the correct value of q). We will see in the following sections that it is possible to obtain data-based approximantsĥ which satisfy the error bound (12) with arbitrarily large q. Nevertheless, this process requires that q is known before constructingĥ, but in order to apply Theorem 3.2 we need q ≥ p, and the value of p depends on the behavior of the Lyapunov function V , which in turns depends onĥ itself. This circular argument may rule out the possibility to design a systematic way to find the correct value of q.
To solve this issue, it would be interesting to know if the exponent p in (10) can be deduced from properties of the full system alone, independently ofĥ. Namely, one may expect that if the exact reduced system has a Lyapunov function V h which satisfies (10) with a given p, then any sufficiently accurate approximantĥ would define an approximated reduced system (9) with a Lyapunov function Vĥ which satisfies (10) with the same p.
This kind of results may be possibly derived from results on the perturbation of Lyapunov functions (see e.g. [20]).

Evolution of the system along the approximated reduced dynamics
An approximated center manifold can also be used to define approximated trajectories of the reduced system. In the following statement we show that the approximation error ofĥ can be used to bound the deviation between the true and the reduced center dynamics.
In the following we denote as diam(B) the diameter of the set B.
Furthermore, define the Lipschitz constants If x 0 ∈ B, and if x(t), z(t) are the solutions of the initial value problems then for all t > 0 it holds Proof. We define g : R n → R n as and set µ B := sup x∈B g(x) . Then a classical perturbation result (see e.g. Theorem 3.4 in [17]) states that Following the same steps as in the proof of Theorem 3.2 and using the bound (23) we derive the bound and this concludes the proof.

Kernel-based approximation of the center manifold
We discuss now a concrete approach for the construction of a data-based approximation of the center manifold, based on kernel methods. The algorithm that we are discussing has been introduced in [12]. We recall it here for completeness, and in addition we present a proof of the optimality of the approximant and we derive error bounds that may be combined with Theorem 3.2 and Theorem 3.3. The approximantĥ to the center manifold is constructed in a data-based way, i.e., we assume to know the values taken by the true manifold h on a finite set of input values. To collect these input-output pairs one can first derive the splitting (2) by analyzing the eigenvalues and eigenvectors of the linearization L := D x f (x) |x=0 ∈ R n×n , as discussed in Section 2. When this splitting is known, it is possible to consider several initial conditions (x 0 , y 0 ) and the corresponding Initial Value Problems (IVPs) given by the coupling of the initial values with the dynamics (2). These IVPs may be solved numerically with high accuracy to collect sequences of pairs (x i , y i ) that are numerical trajectories of the IVP at a given time instant, and which are collected into a dataset X N : for a certain N ∈ N. Thanks to Theorem 2.1, we know that if the point (x i , y i ) is sufficiently close to (0, 0), then y i ≈ h(x i ). These sets X n , Y n can then be regarded as collections of sufficiently accurate pairs of input-output values of the true manifold h.
Observe that a center manifold h needs not to be unique, as in general multiple solutions of the PDE (5) may exist. Nevertheless, the solutions may differ at most by a term exponentially decaying to zero, which may thus be neglected (in terms of approximation) if the points (x i , y i ) are sufficiently close to the origin. To be more precise, it may be that the points in X N , Y N of larger magnitude lie on different manifolds, but the distance between these different branches becomes negligible compared to the approximation error, which is assumed to be algebraic (see Theorem 3.2).

Definition, existence and optimality of the approximant
For the construction of the approximant we employ here matrix-valued kernels, which are suitable to approximate vector-valued functions. Details on kernel-based approximation can be found e.g. in [30], and the extension to the vectorial case is detailed e.g. in [21,35]. Here we recall only that a positive definite matrix-valued kernel on a set Ω is a function k : Ω × Ω → R m×m such that • k(x, y) = k(y, x) T for all x, y ∈ Ω, • the matrix (k(x i , x j )) N i,j=1 ∈ R mN ×mN is positive semidefinite for any set {x 1 , . . . , x N } ⊂ Ω of pairwise distinct points, for all N ∈ N.
Associated to a positive definite kernel there is a unique Hilbert space H of vector valued functions Ω → R m , named native space, where the kernel is reproducing, meaning that k(·, x)α is the Riesz representer of the directional point evaluation δ α x (f ) := α T f (x), for all α ∈ R m , x ∈ Ω. This is the functional setting where the following approximation process is defined.
In our case we aim at anĥ that provides an error bound like (12), and to achieve this we consider a kernel-based approximant that includes derivative information. For now we concentrate on a case that will provide an error bound of order q = 2 since it is somehow more natural given the definition of the center manifold, and we will discuss the extension to general q in the following.
To make this formal, we consider a twice continuously differentiable matrix-valued kernel k on Ω and define the approximant as the minimizer of a cost functional on H, which is defined according to the following three requirements.
Second, the approximant should provide a small distance between the measured output y i and the predictionĥ(x i ), for all (x i , y i ) in the dataset. To this end we introduce for each index i a positive definite weight matrix ω i ∈ R m×m which associates a possibly different importance to different data points and to different output dimensions. Using these matrices, we define Third, as discussed at the beginning of this section we do not know if a data pair (x i , y i ) lies exactly on the center manifold, i.e., if y i = h(x i ) holds, and thus an approximant which merely interpolates the data seems ill-suited for our purposes. For this reason, we introduce a regularization term in the functional to allow non-interpolatory solutions, as long as they have a small norm.
All together, we compute our approximant by minimizing over s ∈ H the functional J : under the constraint s(0) = 0, Ds(0) = 0.
for some λ > 0 independent of i. However, one might consider a more general approach, where the weight increases as the data tends to the origin, i.e. ω i ω j if x i ≤ x j , to give more importance to the samples that are expected to represent more closely the center manifold.
We now show that the problem that we just defined via (25) has a unique minimizerĥ. The proof is taken from [33], and we report it here for the sake of completeness.
has a unique solutionĥ ∈ H which can be represented aŝ where x N +1 := 0, and where the superscript ∂ (2) denotes that the derivative with regards to the second kernel component is taken. Moreover, the coefficient vectors α i , β i ∈ R m can be computed by solving the system with Proof. Letĥ ∈ H be given by (27). Then it follows thatĥ(X) = Y − W α,ĥ(0) = 0 and Dĥ(0) = 0, i.e.ĥ lies in the admissible set for the minimization problem. Let s ∈ H be an alternative admissible element. Then there exists g ∈ H with g(0) = 0, Dg(0) = 0 such that s =ĥ + g. It now holds Hence, it is sufficient to show that Since k is the reproducing kernel of H, and since Dg(0) = 0, it holds On the other hand we have Therefore, and equality holds if and only if g = 0. Hence,ĥ is the unique minimizer of J.
Observe in particular that the expression (27) allows easy manipulation of the approximant h. For example, first order derivatives are simply given by the formula

Remark 4.2 (Higher order derivatives).
Similarly to the case defined in Theorem 4.1, it is possible to define a modified optimization problem where the same cost functional J from (25) is used, but additionally also constraints on higher order derivatives are set. More specifically, for any p > 1, p ∈ N, one may consider the set of multi-indexes of the form a := (a 1 , . . . , a d ) ∈ N d with |a| := d i=1 a i , and require that where ∂ a f (x) := ∂ a 1 1 · . . . ∂ a d d f (x) for f : R d → R. Assuming that the kernel and the manifold h are smooth enough to impose these conditions, the optimization problem still has a unique solution of the form (27), but where now the second sum is a linear combination of all the p-th order derivatives of K. Moreover, the coefficients can still be found by solving a linear system like (28), augmented with blocks corresponding to higher order derivatives.
The main limitation of this approach is that imposing the conditions (29) requires the knowledge of ∂ a h i (0), 1 ≤ i ≤ m, for all |a| ≤ p and p > 1. These values are not known in general, unlike in the case p = 1 where they are defining properties of the center manifold.

Selection of the sampling points
If the technique of the previous section is used as it is, the approximant (27) is given by an expansion with N + 1 + m terms, where N is the number of points in the training set X N , Y N . Therefore, it might be inefficient to evaluate a model that is built using a large dataset. Furthermore, the computation of the coefficients in (27) requires the solution of the linear system (28) with m(N + 1) + m 2 rows and columns, and which can be severely illconditioned for non well-placed points X N .
To mitigate both problems, we employ an algorithm that aims at selecting small subsets XN , YN of X N , Y N such that the approximant computed with these sets is a sufficiently good approximation of the one which uses the full sets. The algorithm selects the points in a greedy way, i.e., one point at a time is selected and added to the current training set. In this way, it is possible to identify a good set without the need to solve a nearly infeasible combinatorial problem.
The selection is performed using the P -greedy method (see [10], and [34] for the matrixvalued case) applied to the kernel K, such that the set of points is selected before the computation of the approximant. The number of pointsN, and therefore the expansion size and evaluation speed, is depending on a prescribed target accuracy ε tol > 0. For details on the method implementation and its convergence properties we refer to [10,26,32,34,27].

A kernel-based center-manifold theorem
Our goal is now to derive error bounds for the approximation process of Theorem 4.1. We aim to use these results in combination with Theorem 3.2 and Theorem 3.3, so we need bounds which contain terms that depend on x q , with a suitable power q.
Comprehensive approximation results are available if plain kernel-based interpolation or regularized interpolation are used instead of the approximant of Theorem 4.1, i.e., if no constraint is enforced. For scalar-valued target functions we refer to [30] for approximation results on interpolation and to [31] to account additionally for a constant regularization term. The extension to vector-valued outputs and general regularization functions has been very recently presented in [33].
In general terms, to obtain convergence we first assume that the sets XN are increasingly fine. A way to measure this is to require that the fill distance θN of XN in Ω, defined as is decreasing asN increases. Observe that in our case it is sufficient that Ω is a neighborhood of zero and not the full space where the ODE is defined.
Under the assumption that h is in the Hilbert space H, it holds where F (θN ) is a value that converges to zero as θN converges to zero, with a rate that depends on the particular kernel. We remark that both exponential and algebraic rates of convergence can be achieved (see e.g. [30] for a detailed presentation of this kind of error bounds in the scalar-valued case, and [35] for an extension to the vector-valued case). As a special case, one can consider sets of points such that with a dimension-dependent constant c d . In this case, the bound on the right hand-side of (30) becomes a function ofN only. We remark that the P -greedy algorithm of Section 4.2 selects sets of points that have this property at a very convenient computational cost (see [26]). The condition h ∈ H may be difficult to verify or to achieve in general. Nevertheless for some kernels one gets the notable situation where H is norm equivalent to the Sobolev space H t (Ω; R m ) with some t > d/2. We refer to [30] for a discussion of this property for m = 1, and to [33] in the general case.
In this case the convergence results require only the smoothness assumptions h ∈ H t (Ω; R m ), which seems reasonable in this context, and the error bound (30) reads This result is clearly not suitable to be used in combination with Theorem 3.2 and Theorem 3.3, since the upper bounds on the error are uniform in Ω, while the cited theorems require the error to decrease towards x = 0 with a given speed. The situation is different with the approximant defined in Theorem 4.1 and in this case we can prove the following theorem, which actually provides suitable bounds with q = 2.
Theorem 4.2. Let XN ⊂ Ω \ {0} be a set of pairwise distinct points and let θN be its corresponding fill distance. Assume that K is such that H is norm equivalent to the Sobolev space H t (Ω; R m ) with t > d/2 + 2. Then, if h ∈ H t (Ω; R m ), and ifĥ is the approximant of Section 3, it holds where C > 0 is a constant independent ofN , h and x, and where i.e., r max = λ if the constant regularization function (26) is used. In particular, if r is chosen so that r max ≤ θ t−d/2 N , then we have the bound Proof. By assumption onĥ we havê and therefore, by using a Taylor expansion with suitable remainder, there exists ξ ∈ B x (0) : Now, Theorem 3.2.6 in [33] for |a| = 2 provides the bound and inserting this result in (33) proves (32). The final bound is obtained by substituting the upper bound on r max .
We avoid reporting the full proof of this fact since it requires no technical modification over the proof of Theorem 4.2, but rather just the use of a more involved index set.
for a suitable constant c ′ . We recall that the condition (31) is satisfied by the points selected by the P -greedy algorithm (see [32]).

Numerical Examples
We test now our method on three different examples. In each example we start by identifying the splitting in (2) manually, and then by solving numerically the IVPs obtained by coupling the ODE system (2) with the 2 n different initial values (x 0 , y 0 ) ∈ {±0.8} n , where n is the dimension of the full system in each example. As numerical solver we use in all cases the implicit Euler scheme for initial time t 0 = 0, final time T = 1000 and with the time step ∆t = 0.1. The points (x i , y i ) lying on the computed numerical trajectories are first filtered to keep only those values that are sufficiently close to the equilibrium 0 ∈ R n , and this is done by removing all points outside of the cube [−0.1, 0.1] n . These first steps generate a dataset of N pairs (x i , y i ), on which the greedy algorithm of Section 4.2 is run to possibly further reduce the set of points to a final number ofN pairs. The tolerance ε tol which determines the stopping of the greedy algorithm (see Section 4.2) is set to values that are specified in each example. The resulting datasets (XN , YN ) are finally used to construct the approximants by solving the linear system (28), using in each case a constant regularization function with suitable values of λ (see Remark 26). Moreover, all the following center manifolds are scalar-valued, and thus plain positive definite kernels (i.e., matrix valued kernels with m = 1) are used. Further details on the setting and the parameters used to build and visualize our approximation to the center manifold are specified in each single example.
As a comparison, we also compute an algebraic approximation of the center manifold with the toolbox [25], which provides an approximated expansion of h on the d-dimensional monomial basis. The toolbox allows to compute approximants up to a given monomial degree, and in the following we always use the maximal possible value, which is 5. The toolbox requires as inputs the right hand side of the ODE, the value of the eigenvalues with zero real part, and a basis of the corresponding eigenspaces both for the linearization matrix and for its transpose.
To compute errors, we use an equally spaced test grid X test made of 1001 points in [−0.1, 0.1] for the first two examples, and of 101 2 points in [−0.1, 0.1] 2 in the third one. Additionally, using the PDE (5) we define the pointwise residual and we use it as a further way to quantify how well the approximated manifold models the true one.
Remark 5.1 (Polynomial kernel). Since an algebraic approximation of the manifold is sufficient in many cases, the polynomial kernel k(x, y) := (1 + γx T y) p , γ > 0, p > 0, seems particularly promising for the task. This kernel is positive definite, but not strictly positive definite. Moreover, the approximant (27) in this case takes the particularly simple form i.e., the terms corresponding to the interpolation conditions in zero are adding a linear term to the approximant.

Example 1
We consider the 2-dimensional systeṁ Observe that in this case an explicit form of the true center manifold is not known.
Using the algebraic approximation toolbox [25], we get the approximation We then generate the dataset for the construction of the kernel approximants, which contains N = 38248 data pairs. Starting from these we run the greedy algorithm with ε tol = 10 −15 and with several setups: We use the polynomial kernels which result in the sets X 4 pol , X 5 pol , X 6 pol , X 1 Gauss and X 5 Gauss , containing 14, 15, 15, 7 and 8 points, respectively. The sets of target values are selected from Y N accordingly.
As an example, for the kernels k 4 pol and k 1 Gauss the approximant and the corresponding residual are depicted in Figure 1. In the figure we show also the approximationĥ alg and the residual r alg (x) = 12x 6 + 16x 8 , which is obtained by manually computing the value in (34) using the first two terms in (36). It is clear from the left panel in Figure 1 that the three approximants are very close. Moreover, the magnitude of the residual (right panel in Figure  1) is in all the three cases of the order of 10 −5 , with the algebraic approximation giving the smallest one. For each approximant, we first check how closeĥ is to meet the constraints of the optimization problem of Theorem 4.1. The first two columns of Table 1 show the valuesĥ(0) and Dĥ(0), which should be zero and indeed are in all cases at least smaller than 10 −9 in magnitude, and exactly zero for k 5 poly and k 5 gauss . We argue that the non-zero values are due to numerical inaccuracies when solving the system (28), and they should be considered numerical zeros. Indeed, the interpolation error max x i ∈XN |ĥ(x i ) − y i | on the training set, which is shown in the third column of the table, is several order of magnitudes larger than these non-zero values.
We further verify if the approximated manifolds can be used to derive information on the asymptotic stability of the reduced system. The reduced system in this case readṡ and thus x = 0 is an equilibrium independently of the values ofĥ(x). Since the algebraic approximation (36) is correct up to an error of O(x 5 ), we may assume that this equilibrium is asymptotically stable, since the right hand side is negative on the left of x = 0 and positive on the right. Under this assumption, we may conclude that an approximantĥ provides the correct characterization of the stability of the equilibrium of the reduced system if it is strictly negative in a neighborhood of zero, and zero only in zero. Looking at the fourth and first columns of Table 1, it is the case that max x∈Xtest\{0}ĥ (x) < 0 andĥ(0) ≤ 0 in all cases except forĥ 6 poly , and thus these four approximants are sufficiently accurate to predict the stability of the equilibrium. In the case ofĥ 6 poly we have still that max x∈Xtest\{0}ĥ (x) < 0, but the approximant is positive in zero. Since it is continuous, a finer test grid would find non negative values also for x = 0, and thus no stability is found in this case.
Finally, we expand the approximants as Taylor series of order 4 to draw a further comparison to the results given by [25]. These expansions are obtained by computing explicitly the ĥ (0) expansions of the kernel terms in (27) and then numerically computing the corresponding linear combinations. The coefficients of the Taylor expansions are given in Table 2. The accuracy also in this case is very good, even if it is degrading as the degree of the monomials increases.
For the first three terms the best approximation is provided by the approximants using the Gaussian kernel.

Example 2
We consider the 2-dimensional systeṁ This example is of relevance for this paper especially because an exact center manifold h(x) = x 2 is explicitly known (see [23], and [24] for a recent analysis of this derivation and a study of other properties of the system). Using the toolbox in [25], we get that h(x) = x 2 + O(x 5 ). For the kernel approximation this time we use the compactly supported Wendland kernel for space dimension p d = 1 and smoothness 2p k , with p k = 1 (see [29]) given by whose RKHS is the Sobolev space H p d /2+1/2+p k (Ω) = H 2 (Ω). Therefore, we have h ∈ H k wendland and hence we are able to apply the error theory of Section 4.3 and Theorem 4.2.
For this purpose we again run the P -Greedy algorithm with a tolerance of ε tol = 10 −10 which results in a point set X wendland = {x 1 , . . . , x 200 } containingN = 200 elements. Observe that these points are selected with an ordering corresponding to a sequentially improving covering of the full dataset (see [10]). To test the behavior of the error, we construct the full sequence of approximantsĥN corresponding to each of the incremental point sets XN := {x 1 , . . . , xN }, 1 ≤N ≤ 200. We use the regularization value λ := 10 −13 , and for each approximant we compute the error This discrete error serves as a proxy of the maximum error in (32). The values of eN are shown in the second row in Figure 2, and they clearly have a quadratic rate of decay to zero as a function ofN. Observe that this rate of convergence is significantly faster than the one predicted by Theorem 4.2 and Remark 4.4, which is equal to −s/d + 1/d + 1/2 = −2 + 3/2 = −1/2.
The pointwise error for the final approximantĥ 200 is depicted instead in the first row in Figure 2, left, and it is below 2 · 10 7 in the entire interval. A plot of the error restricted to the smaller interval [−0.01, 0.01], and computed by scaling X test to this interval, is shown in the right part of the the first row, and it is compared with the decay to zero of the functions |x| 3 and |x| 4 . On a qualitative level, this figure seems to suggests that the error around zero (see Theorem 3.2 and Theorem 4.2) is of order q = 3, but not q = 4.

Example 3
We consider now the (2 + 1)-dimensional systeṁ In this case the toolbox [25] returns an implicit parametrization of the manifold depending on the variables t, s 1 , s 2 , and given by the equations This system requires some manipulation to be in the form that we are interested in, since it is not directly providing an expression of y in terms of x 1 , x 2 . Observe that this is a limitation to the use of this methodology as a model order reduction technique.
To derive the expression for the manifold we invert the first two equations in (39), which are linear and invertible, to obtain and these imply that . Inserting this relation into the third equation in (39) gives finally where also in this case the approximation is limited to terms of degree at most 4. The initial dataset to compute the kernel approximants in this case consists of N = 78796 data pairs. We use the kernels k 4 poly (x, y) := (1 + x T y/2) 4 and k 1/2 Gauss (x, y) := e − x−y 2 2 /2 , and the greedy-selected sets with ε tol := 10 −10 have sizeN = 10 in both cases.

1/2
Gauss can be seen in Figure 3. The plots show a good accuracy for both kernels and among the two the polynomial kernel provides the best results. Observe that the error is small especially in a neighborhood of zero, as expected.
Following Theorem 3.3, we now try to simulate the evolution of the system (38) using the approximated reduced dynamics instead of the exact reduced dynamics. Since we do not know the center manifold explicitly, the trajectories of the exact reduced system are obtained by solving the full system, and then considering only the first two components of each trajectory point.
In more details, we fix an initial condition (x 1 , x 2 , y) 0 (and its coordinate projection (x 1 , x 2 ) 0 ) and we solve numerically the exact system (38) and the approximated systeṁ with the same setting used in all the experiments. For t i = t 0 , t 0 + ∆t, . . . , T we obtain a trajectory {(x 1 , x 2 , y) t i } T t i =t 0 of the full system, from which we consider the first two coordinates {(x 1 , x 2 ) t i } T t i =t 0 , and a trajectory {(x 1 ,x 2 ) t i } T t i =t 0 of the approximated and reduced system. We consider the two initial conditions x repr 0 := (x 1 , x 2 , y) 0 := (−0.8, −0.8, −0.8) and x gen 0 := (x 1 , x 2 , y) 0 := (−0.4, 0, −0.8). The first one has been used for the training ofĥ 4 poly while the second is new, and thus these are intended as a reproduction and a generalization experiment.
The corresponding exact and approximated numerical trajectories are shown in the first row of Figure 4. It is clear that for x repr 0 the evolution of the approximated system remains closer to the exact one than for x gen 0 , which is the expected behavior. Nevertheless, it is remarkable that in both cases the approximated trajectory converges to the same final state, which is the equilibrium (0, 0), and thus the reduced accuracy of the generalization case is not relevant for the long term evolution of the system. This fact may be explainable by looking at the second row of Figure 4, where we report the absolute and relative distance between the exact and reduced numerical trajectories in the two cases over the entire time interval. Indeed, the absolute error is decaying exponentially for both x repr 0 and x gen 0 with a rate which is numerically estimated to be e −0.05t i . This decay stops at a plateau value at around t i = 400. This exponential decay of the error is nevertheless driven by the fact that the norm x 2 1 + x 2 2 of the true numerical solution is exponentially decaying with roughly the same rate, until it reaches a final value (the approximation of the equilibrium) at t i ≈ 400. Indeed the relative error, for t i ≤ 350, is rapidly stabilizing to a constant value of 2.97% (reproduction) and 20.65% (generalization). Thus, a pretty large error in terms of approximation is in fact small in terms of solutions of the IVP.
This observation suggests once again that our approximation scheme may be used as an effective model order reduction technique, since an approximant computed using a few initial values is able to provide meaningful predictions for unseen parameters. This numerical observation would deserve further theoretical investigation.  Figure 4: The first row shows numerical trajectories computed with the full system (blue curves) and with the reduced and approximated one (red curves) for an initial value x repr 0 (left) and for an initial value x gen 0 (right). The second row shows the evolution over time of the absolute error compared with the rate e −0.05t i (left, with logarithmic scaled) and of the relative error (right, with linear scale and restricted to t i ≤ 350) between the full and reduced trajectories.
As a final analysis, we compute the Taylor expansion in (0, 0) of h 4 poly . Table 3 reports the coefficients of the expansion for the monomials up to degree four, computed as in Section 5.1. The terms of degree up to 2 have good agreement withĥ alg (which is exact up to degree four). The terms of degree three should be zero, and they are of magnitude 10 −3 , while the accuracy is degrading for the highest order terms, where no clear pattern nor resemblance with the exact values can be seen.

Conclusions
In this paper we proved a general data-based version of the main implication of the center manifold theorem, which allows to draw the same conclusions as the original theorem but examining a sufficiently accurate approximation in place of the true manifold. Additionally,  Table 3: Coefficients of the Taylor expansion in x = 0 of the approximated center manifold h 4 poly up to degree four.
such approximation can be used to evolve the system along approximated trajectories, and we proved bounds on their deviation from the exact ones. Moreover, we revisited an algorithm that we introduced in a previous work for the kernelbased approximation of the center manifold. We added details on its definition and optimality, and especially we proved rigorous bounds on the error between this approximant and the true manifold, and discussed how to use these bounds in combination with the theorems mentioned above.
Extensive numerical experiments suggested that the present method can reach a significant accuracy, and that it has the potential to be used as an effective technique to analyze the stability of dynamical systems.
Several points remain open and would deserve further investigation. Mainly, establishing a connection between the exponent q in Theorem 3.2 and the full order dynamical system would allow a more systematic application of the results of this paper on different problems. Moreover, results on the generalization power of the method for the simulation of IVPs with initial values not used during training are missing.