KERNEL METHODS FOR THE
APPROXIMATION OF NONLINEAR SYSTEMS

arXiv:1108.2903v3 [math.OC] 31 Mar 2016

JAKE BOUVRIE† ¶ AND BOUMEDIENE HAMZI‡ ¶
Abstract. We introduce a data-driven model approximation method for nonlinear control
systems, drawing on recent progress in machine learning and statistical dimensionality reduction. The
method is based on embedding the nonlinear system in a high (or infinite) dimensional reproducing
kernel Hilbert space (RKHS) where linear balanced truncation may be carried out implicitly. This
leads to a nonlinear reduction map which can be combined with a representation of the system
belonging to a RKHS to give a closed, reduced order dynamical system which captures the essential
input-output characteristics of the original model. Working in RKHS provides a convenient, general
functional-analytical framework for theoretical understanding. Empirical simulations illustrating the
approach are also provided.

1. Introduction. Data-based modelling of nonlinear dynamical systems has
been addressed by many authors. For example, several methods have been developed in Time Series Analysis ([19] for example) and System Identification ([52], [51]
for example). Coifman et al. discuss data-based modelling of a stochastic Langevin
system [9]. Archambeau et al. [2] proposed methods to approximate SDEs from data.
Smale and Zhou use kernel methods to approximate a hyperbolic dynamical system
[45].
In this paper we propose a scheme for the approximation of nonlinear systems
using balanced model-order reduction. A key, and to our knowledge, novel point
of departure from the literature on nonlinear model reduction is that our approach
marries approximation and dimensionality reduction methods known to the machine
learning and statistics communities with existing ideas in linear and nonlinear control. In particular, we apply a method similar to kernel PCA (Principal Component
Analysis) as well as function approximation in Reproducing Kernel Hilbert Spaces
(RKHSes) to the problem of balanced model reduction. Working in RKHSes provides
a convenient, general functional-analytical framework for theoretical understanding as
well as a ready source of existing results and error estimates. The approach presented
here is also strongly empirical, in that observability and controllability, and in some
cases the dynamics of the nonlinear system are estimated from simulated or measured
trajectories.
The approach we propose begins by viewing the controllability and observability
energies for nonlinear systems as Gramians in a high (or possibly infinite) dimensional RKHS. These Gramians are approximated empirically and then simultaneously
diagonalized in order to identify directions which, in the RKHS, are both the most
observable and the most controllable.
The assumption that it is possible to apply the method of linear balancing to
a nonlinear system when lifted to a RKHS is far more reasonable than applying
the linear theory in the original space hoping for the best. Working in the high
dimensional RKHS allows one to perform linear operations on a representation of
the system’s state and output which can capture strong nonlinearities. Moreover,
† Laboratory for Computational and Statistical Learning, Massachusetts Institute of Technology,
Cambridge, MA 02138, USA. (jvb@csail.mit.edu)
‡ Department of Mathematics,
Imperial College London, London SW7 2AZ, UK.
(b.hamzi@imperial.ac.uk)
¶ Parts of this work were done while both authors were at Department of Mathematics, Duke
University, Durham, NC 27708, USA.

1

working in an RKHS readily include linear spaces1 and, therefore, our approach not
only covers existing Linear Theory but also extends the range of available spaces
where it is reasonable to assume that one will obtain better results when dealing with
a nonlinear problem. Therefore, a system for which existing model reduction methods
fail, may be approximated by a lower dimensional system when mapped into a RKHS.
This situation closely parallels the problem of linear separability in data classification:
A dataset which is not linearly separable might be easily separated when mapped into
a nonlinear feature space. The decision boundary is linear in this feature space, but
is nonlinear in the original data space. Essentially, we are proposing to apply linear
methods to nonlinear systems once mapped into a high (possibly infinite-dimensional)
Hilbert space2 .
Nonlinear reduction of the state space already allows to the design of simpler
controllers, but is only half of the picture. One would also like to be able to write a
closed, reduced dynamical system whose input-output behavior closely captures that
of the original system. This problem is the focus of the second half of our paper,
where we again exploit helpful properties of RKHS in order to provide such a closed
system.
The paper is organized as follows. In the next section we provide the relevant
background for model reduction and balancing for linear and nonlinear control systems. We then adapt and extend balancing techniques described in the background
section to the current RKHS setting in Section 3. Section 4 then proposes a method
for determining a closed, reduced nonlinear control system in light of the reduction
map described in Section 3. Finally, Section 5 provides experiments illustrating an
application of the proposed methods to some nonlinear systems where the method of
linear balancing does not apply in Rn since the systems we simulated are not linearly
controllable and the origin is not asymptotically stable but the same method of linear
balancing applies to the nonlinear system after being lifter to a RKHS. Appendix A
contains a review of Learning Theory, Appendix B contains a description of kernel
PCA.
Preliminary results of this work can be found in [5, 6].
2. Background. Several approaches have been proposed for the reduction of linear control systems in view of control, but few exist for finite or infinite-dimensional
nonlinear control systems. For linear systems, the pioneering “Input- Output balancing” approach proposed by B.C. Moore [29] observes that the important states are the
ones that are both easy to reach and that generate a lot of energy at the output. If a
large amount of energy is required to reach a certain state but the same state yields
a small output energy, the state is unimportant for the input-output behavior of the
system. The goal is then to find the states that are both the most controllable and the
most observable. One way to determine such states is to find a change of coordinates
where the controllability and observability Gramians (which can be viewed as a measure of the controllability and the observability of the system) are equal and diagonal.
States that are difficult to reach and that don’t significantly affect the output are
then ignored or truncated. A system expressed in the coordinates where each state is
1 When using the embedding Φ : x 7→ x and kernel k(x, y) = hx, yi, cf. Appendix A for meaning
of these objects.
2 Let us note here that it is important to choose the right RKHS in order to perform this sort of
computations. In our work, we used existing Universal Kernels such as the Gaussian or polynomial
kernels (cf. Proposition A.4) but, in general, properly choosing the right RKHS is an open problem
even in classical Learning Theory.

2

equally controllable and observable is called its balanced realization.
A proposal for generalizing this approach to nonlinear control systems was advanced by J. Scherpen [36], where suitably defined controllability and observability energy functions reduce to Gramians in the linear case. In general, to find the
balanced realization of a system one needs to solve a set of Hamilton-Jacobi and
Lyapunov equations (as we will discuss below). Moore [29] proposed an alternative,
data-based approach for balancing in the linear case. This method uses samples of the
impulse response of a linear system to construct empirical controllability and observability Gramians which are then balanced and truncated using Principal Components
Analysis (PCA, or Proper Orthogonal Decomposition (POD) [17]). This data-driven
strategy was then extended to nonlinear control systems with a stable linear approximation by Lall et al. [25], by effectively applying Moore’s method to a nonlinear
system by way of the Galerkin projection. Despite the fact that the balancing theory
underpinning their approach assumes a linear system, Lall and colleagues were able
to effectively reduce some nonlinear systems.
Phillips et al. [32] has also studied reduction of nonlinear circuit models in the case
of linear but unbalanced coordinate transformations and found that approximation
using a polynomial RKHS could offer computational advantages. Gray and Verriest mention in [15] that studying algebraically defined Gramian operators in RKHS
may provide advantageous approximation properties, though the idea is not further
explored. Finally, Coifman et al. [9] discuss reduction of an uncontrolled stochastic Langevin system. There, eigenfunctions of a combinatorial Laplacian, built from
samples of trajectories, provide a set of reduction coordinates but does not provide
a reduced system. This method is related to kernel principal components (KPCA)
using a Gaussian kernel, however reduction in this study is carried out on a simplified
linear system outside the context of control.
In the following sections we review balancing of linear and nonlinear systems as
introduced in [29] and [36]. See also [37] for a good survey on balancing for linear and
nonlinear systems.
2.1. Balancing of Linear Systems. Consider a linear control system
ẋ =
y =

Ax + Bu,
,
Cx,

(2.1)

where (A, B) is controllable, (A, C) is observable and A is Hurwitz. We define the
controllability and the observability Gramians as, respectively,
R∞
R∞ >
>
Wc = 0 eAt BB >eA t dt, and Wo = 0 eA t C >CeAt dt.
These two matrices can be viewed as a measure of the controllability and the observability of the system [29]. For instance, consider the past energy [36, 35], Lc (x0 ),
defined as the minimal energy required to reach x0 from 0 in infinite time
Z
1 0
||u(t)||2 dt,
(2.2)
Lc (x0 ) =
inf
2 −∞
u∈L2 (−∞,0),
x(−∞)=0,x(0)=x0

and the future energy [36, 35], Lo (x0 ), defined as the output energy generated by
releasing the system from its initial state x(t0 ) = x0 , and zero input u(t) = 0 for
t ≥ 0, i.e.
Z
1 ∞
||y(t)||2 dt,
(2.3)
Lo (x0 ) =
2 0
3

for x(t0 ) = x0 and u(t) = 0, t ≥ 0. In the linear case, it can be shown that
−1
Lc (x0 ) = 21 x>
0 Wc x0 ,

and Lo (x0 ) = 12 x>
0 W o x0 .

(2.4)

The columns of Wc span the controllable subspace while the nullspace of Wo coincides
with the unobservable subspace. As such, Wc and Wo (or their estimates) are the key
ingredients in many model reduction techniques. It is also well known that Wc and
Wo satisfy the Lyapunov equations [29]
AWc + Wc A> = −BB >,

A>Wo + Wo A = −C >C.

(2.5)

Several methods have been developed to solve these equations (see [16, 26, 27] for
example). As mentioned at the beginning of this section, it is also possible to estimate
the Gramians from empirical data, cf. (3.4) and (3.6) below.
The idea behind balancing is to find a representation where the system’s observable and controllable subspaces are aligned so that reduction, if possible, consists
of eliminating the states that are least controllable and which are also the least observable. More formally, we would like to find a new coordinate system such that
Σ := Wc = Wo = diag{σ1 , · · · , σn } where σ1 ≥ σ2 ≥ · · · ≥ σn > 0.
Theorem 2.1. [11] If (A, B) is controllable and (A, C) is observable, then the
eigenvalues of Wo Wc are similarity invariants, i.e. they are independent of the choice
of the state-space representation of (2.1). Moreover, there exists a state-space representation where
Σ := Wc = Wo = diag{σ1 , · · · , σn }

(2.6)

with σ1 ≥ σ2 ≥ · · · ≥ σn > 0 are the square roots of the eigenvalues of Wo Wc . Such
representations are called balanced, and the system is in balanced form. The σi ’s, i =
1, ..., n are called the Hankel singular values. The state space expressed in the transformed coordinates (QAQ−1 , QB, CQ−1 ) is balanced and QWc Q> = Q−>Wo Q−1 = Σ
where Q ∈ Rn×n .
The Hankel singular values, σi |ni=1 , are the square roots of the eigenvalues of
Wo Wc and are the singular values of the Hankel operator
H = Ψo ◦ Ψc

(2.7)

that characterizes the input-output behaviour of the system (2.1) with Ψc , the controllability operator, which maps u ∈ L2 (−∞, 0] to x(0) and Ψo , the observability
operator, which maps x(0) to y(t), t ≥ 0 with no input applied for t ≥ 0 [11]. More
precisely,
Ψc : L2 (−∞, 0] → Cn
R0
u 7→ −∞ e−Aτ Bu(τ )dτ

(2.8)

and
Ψ o : Cn
x(0) = x0

→ L
2 [0, ∞)
CeAt x0 , for t ≥ 0,
7→
0, otherwise.

(2.9)

Clearly, x0 = Ψc u(t) for u(t) ∈ L2 (−∞, 0] is the system state at t = 0 due to the
past input and y(t) = Ψ0 x0 , t ≥ 0 is the future output due to the initial state
4

x0 with the input set to zero. In fact, H characterizes the system’s future output
y(t) = Hu(t), t ≥ 0 based on the past input u(t), t ≤ 0. More precisely, if x(−∞) = 0,
Z 0
Hu(t) = Ψo Ψc u(t) =
CeA(t−τ ) Bu(τ )dτ, for t ≥ 0
(2.10)
−∞

When H is known to be a compact operator, then its adjoint operator H∗ is also
compact and the composition H∗ H is a self-adjoint compact operator with the spectral
decomposition
H∗ H =

∞
X

σi2 h·, Ψi iL2 Ψi ,

σi ≥ 0,

(2.11)

i=1

hΨi , Ψj iL2 = δij ,

hΨi , (H∗ H)(Ψi )iL2 = σi2

(2.12)

where σi2 is an eigenvalue of H∗ H with the corresponding eigenvector Ψi , ordered as
σ1 ≥ · · · ≥ σn > 0 and σi≥n+1 = 0 are the Hankel singular values of the input-output
system Σ. For square linear systems, the nonzero eigenvalues of the Hankel operator
associated to the system are the nonzero eigenvalues of the cross Gramian defined as
the solution, Wx , of AWx + Wx A + BC = 0 [37].
We also have, for every x0 ∈ Cn ,
Ψc Ψ∗c x0 = Wc x0 ,

Ψ∗o Ψo x0 = Wo x0 .

(2.13)

Thus Wc and Wo are the matrix representations of the operators Ψc Ψ∗c and Ψ∗o Ψo
[11].
For model reduction, typically one looks for a gap in the singular values {σi }
for guidance as to where truncation should occur. If there is a k such that σk 
σk+1 , then the states most responsible for governing the input-output relationship
of the system are (x1 , · · · , xk ) while (xk+1 , . . . , xn ) are assumed to make negligible
contributions.
Theorem 2.2. [11] Consider a linear system (2.1) with its associated Hankel
operator H (2.7). Let Hk be the Hankel operator of the reduced order linear system of
order k. Then,
||H − Hk || = σk+1

(2.14)

If x(0) = 0, the error between y, the output of the full order system, and yr , the output
of the reduced order system with k state variables, satisfies
 X

n
||y − yr ||2 ≤ 2
σj ||u||2
(2.15)
j=k+1

If F is unstable then the controllability and observability quantities defined in (2.2)
are undefined since the integrals will be unbounded. There may, however, still exist
solutions to the Lyapunov equations (2.5) when F is unstable [48, 20]. Other approaches to balancing unstable linear systems exist (see [49, 50, 18, 54] for the method
of LQG balancing for example).
Although several methods also exist for computing Q [26, 27], it is common to
simply compute the Cholesky decomposition of Wo so that Wo = ZZ >, and form the
1
SVD U Σ2 U > of Z >Wc Z. Then Q is given by Q = Σ 2 U >Z −1 . We also note that the
problem of finding the coordinate change Q can be seen as an optimization problem [1]
of the form minQ trace[QWc Q∗ + Q−∗ Wo Q−1 ].
5

2.2. Balancing of Nonlinear Systems. In the nonlinear case, the energy functions Lc and Lo in (2.2) and (2.3) are obtained by solving both a Lyapunov and a
Hamilton-Jacobi equation. Here we follow the development of Scherpen [36, 37]. Consider the nonlinear system

Pm
ẋ = f (x) + i=1 gi (x)ui ,
(2.16)
Σ:
y = h(x),
with x ∈ Rn , u ∈ Rm , y ∈ Rp , f (0) = 0, gi (0) = 0 for 1 ≤ i ≤ m, and h(0) = 0.
Moreover, assume the following Hypothesis.
Assumption A: The linearization of (2.16) around the origin is controllable, observable
and A = ∂f
∂x |x=0 is asymptotically stable.
The controllability operator Ψc : U → X with X = Rn and U = Lm
2 [0, ∞), and
the observability operator Ψo : X → Y with Y = Lp2 [0, ∞) for this system are defined
by

ẋ = −f (x) − g(x)u, x(∞) = 0,
0
Ψc : u 7→ x :
(2.17)
x0 = x(0)

0

Ψo : x 7→ y :



ẋ = f (x),
y = h(x)

x0 = x(0),

(2.18)

As in the linear case, Ψc and Ψo represent the input-to-state behavior and the stateto-output behavior, respectively. The Hankel operator for the nonlinear system Σ in
(2.16) is given by the composition of Ψc and Ψo
H := Ψo ◦ Ψc

(2.19)

Consider the norm-minimizing inverse Ψ†c : X → U
Ψ†c : x0 7→ u := argminΨc (u)=x0 ||u||.

(2.20)

From this point of view Lc in (2.2) and Lo in (2.3) are
Lc (x0 ) :=

1 † 0 2
||Ψ (x )|| ,
2 c

Lo (x0 ) :=

1
||Ψo (x0 )||2
2

(2.21)

Theorem 2.3. [36, 35] Consider the nonlinear system Σ defined in (2.16). If
the origin is an asymptotically stable equilibrium of f (x) on a neighborhood W of the
origin, then for all x ∈ W , Lo (x) is the unique smooth solution of
1
∂Lo
(x)f (x) + h>(x)h(x) = 0,
∂x
2

Lo (0) = 0

(2.22)

under the assumption that (2.22) has a smooth solution on W . Furthermore for all
x ∈ W , Lc (x) is the unique smooth solution of
∂Lc
1 ∂Lc
∂ >Lc
(x)f (x) +
(x)g(x)g >(x)
(x) = 0,
∂x
2 ∂x
∂x

Lc (0) = 0

(2.23)

under the assumption that (2.23) has a smooth solution L̄c on W and that the origin
L̄c
is an asymptotically stable equilibrium of −(f (x) + g(x)g >(x) ∂∂x
(x)) on W .
6

With the controllability and the observability functions on hand, the inputnormal/output-diagonal realization of system (2.16) can be computed by way of a
coordinate transformation. More precisely,
Theorem 2.4. [36, 35] Consider system (2.16) under Assumption A and the
assumptions in Theorem 2.3. Then, there exists a neighborhood W of the origin and
coordinate transformation x = ϕ(z) on W P
converting the energy functions into the
n
form Lc (ϕ(z)) = 12 z >z and Lo (ϕ(z)) = 12 i=1 zi2 σi (zi )2 , where σ1 (x) ≥ σ2 (x) ≥
· · · ≥ σn (x). The functions σi (·) are called Hankel singular value functions.
Analogous to the linear case, the system’s states can be sorted in order of importance by sorting the singular value functions, and reduction proceeds by removing the
least important states.
In the above framework for balancing of nonlinear systems, one needs to solve
(or numerically evaluate) the PDEs (2.22), (2.23) and compute the coordinate change
x = ϕ(z), however there are no systematic methods or tools for solving these problems. Various approximate solutions based on Taylor series expansions have been
proposed [23, 22, 13]. Newman and Krishnaprasad [30] introduce a statistical approximation based on exciting the system with white Gaussian noise and then computing
the balancing transformation using an algorithm from differential topology. As mentioned earlier, an essentially linear empirical approach was proposed in [25]. In this
paper, we combine aspects of both data-driven approaches and analytic approaches
by carrying out linear balancing of nonlinear control systems in a suitable RKHS.
3. Empirical Balancing of Nonlinear Systems in RKHS. We consider a
general nonlinear system of the form

ẋ = f (x, u)
(3.1)
y = h(x)
with x ∈ Rn , u ∈ Rm , y ∈ Rp , f (0, 0) = 0, and h(0) = 0. Let R(x0 ) = {x0 ∈ Rn :
∃ u ∈ L∞ (R, Rm ) and ∃ T ∈ [0, ∞) such that x(0) = x0 and x(T ) = x0 } be the
reachable set from the initial condition x(0) = x0 .
Hypothesis H:3 The system (3.1) is zero-state observable, its linearization around
the origin is controllable, and the origin of ẋ = f (x, 0) is asymptotically stable.
We treat the problem of estimating the observability and controllability Gramians
as one of estimating an integral operator from data in a reproducing kernel Hilbert
space (RKHS) [3], cf. Appendix for definition and key results on RKHSes. Our
approach hinges on the key modeling assumption that the nonlinear dynamical system
can be embedded in an appropriate high (or possibly infinite) dimensional RKHS where
the method of linear balancing can be applied. More precisely, we will essentially
assume that there is an RKHS H and maps Φ, Ψ : Rn → H; x 7→ H such that
controllability and observability energies of the nonlinear system (3.1) are “linearized”,
i.e. that in H they have an expression similar to the one in the linear case (2.4) and,
therefore, can be written as
Lc (x) u

1 T
Φ (x)W−1
c Φ(x),
2

Lo (x) u

1 T
Ψ (x)Wo Ψ(x),
2

(3.2)

with Wc , Wo ∈ RN ×N , N  n, are very large dimensional matrices4 .
3 Let us note here that this assumption is similar to Assumption A above and that we made it
mainly out of convenience but is not necessary as illustrated in the 2D and 7D examples below.
4 If the system is affine in the input, we will use the PDEs (2.22) and (2.23) to find Φ, Ψ, W and
c
Wo . We leave such analysis for future work.

7

By “linearization” here, we mean mapping the nonlinear system in a higher dimensional Hilbert space where linear theory can be applied. To illustrate this point
[46], consider a polynomial in R, p(x) = α + βx + γx2 where α, β, and γ are reals. If we consider the map Φ : R → R3 defined as Φ(x) = [1 x x2 ]T then
p(x) = α · [1 x x2 ]T = α · Φ(x) is an affine polynomial in the variable Φ(x).
Another example to illustrate our thought process is the one of Support Vector
Machines (SVMs) that we referred to in the introduction. More precisely, consider
the problem of classifying points in a data set D = ((x1 , y1 ), · · · , (xn , yn )) with xi ∈
X with X a set and yi = ±1, i.e. trying to find w ∈ Rd with ||w||2 = 1 and
b ∈ R such that hw, xi i + b > 0 for all i with yi = +1, and hw, xi i + b < 0 for
all i with yi = −1, i.e. that the linear hyperplane characterized by (w, b) perfectly
separates the set D into two groups of data points, the ones with yi = +1 and the ones
with yi = −1. Sometimes, it is possible to find such an hyperplane but, in general,
finding a linear hyperplane that perfectly separates a given data set D is not always
possible and finding (w, b) will not be possible. To solve the classification problem, the
SVM algorithm maps the input data (x1 , · · · , xn ) into a (possibly infinite-dimensional)
Hilbert space H, the so-called feature space, by a typically nonlinear map Φ : X → H
called the feature map. Then one looks for a linear hyperplane that separates the data
((Φ(x1 ), y1 ), · · · , (Φ(xn ), yn )), i.e. one looks for (w, b) in H. When this is possible,
the data in D will be classified in two categories, the ones with yi = +1 and the ones
with yi = −1, but the separating curve for D will not be a linear hyperplane in the
original space (even if it is a linear hyperplane in H). An important property of SVMs
is that for every dataset D without contradicting points, i.e. (xi , yi ) and (xj , yj ) with
xi = xj and yi 6= yj , there exists a feature map that allows the perfect separation by
a hyperplane in the feature space [47].
Our aim is to generalize this way of thinking to nonlinear dynamical systems, i.e.
given a problem for a nonlinear dynamical system, we map the state variables by a
typically nonlinear map Φ : X → H where H is (possibly infinite-dimensional) Hilbert
space in which the computations become simpler5 . In our case, “simpler” means
applying Linear Theory. In this paper, we will focus on the problem of approximation
of nonlinear control systems by applying the method of linear balancing in RKHSes.
We leave other applications for future work.
Covariance operators in RKHSes and their empirical estimates are the objects of
primary importance and contain the information needed to perform model reduction.
In particular, the (linear) observability and controllability Gramians are estimated
and diagonalized in the RKHS, but capture nonlinearities in the original state space.
The reduction approach we propose adapts ideas from kernel PCA (KPCA) [41] and is
driven by a set of simulated or sampled system trajectories, extending and generalizing
the work of Moore [29] and Lall et al. [25].
Our method works quite well since the controllability and the observability energies in the linear case can be expressed as inner products (2.4) and working in
RKHSes allows to find nonlinear versions of linear algorithms that can be expressed
in terms of inner products (this is the so-called kernel trick in Learning Theory, see
Appendix A for more explanations). Hence the empirical Gramians defined below for
nonlinear systems can be viewed as reasonable approximations of the controllability
and observability energies for nonlinear systems.
In the development below we lift state vectors of the system (3.1) into a Hilbert
5 One could also think of the methods in Quantum Mechanics where one constructs a Hilbert
space from measurements in order to perform computations.

8

space H, i.e. we consider a mapping Φ : Rn → H and analyze the nonlinear system
whose state is Φ(x).
3.1. Empirical Gramians in RKHS. Following [29], we estimate the controllability Gramian by exciting each coordinate of the input with impulses6 while
setting x0 = 0. One can also further excite using rotations of impulses as suggested
in [25], however for simplicity we consider only the original signals proposed in [29].
Let ui (t) = δ(t)ei be the i-th excitation signal, and let xi (t) be the corresponding
response of the system. Form the matrix


X(t) = x1 (t) · · · xm (t) ∈ Rn×m ,

(3.3)

so that X(t) is seen as a data matrix with column observations given by the respective
responses xi (t). Then Wc ∈ Rn×n is given by
1
Wc =
m

Z

∞

X(t)X(t)>dt.

(3.4)

0

We can approximate this integral by sampling the matrix function X(t) within a finite
time interval [0, T ] assuming the regular partition {ti }N
i=1 , ti = (T /N )i. This leads to
the empirical controllability Gramian
N
X
cc = T
W
X(ti )X(ti )>.
mN i=1

(3.5)

As described in [29], the observability Gramian is estimated by fixing u(t) = 0,
setting x0 = ei for i = 1, . . . , n, and measuring the corresponding system output responses y i (t). As before, assemble the responses into a matrix Y (t) = [y 1 (t) · · · y n (t)] ∈
co are
Rp×n . The observability Gramian Wo ∈ Rn×n and its empirical counterpart W
given by
Z
1 ∞
Wo =
Y (t)>Y (t)dt
(3.6)
p 0
and
N
X
co = T
W
Ye (ti )Ye (ti )>
pN i=1

(3.7)

where Ye (t) = Y (t)>. The matrix Ye (ti ) ∈ Rn×p can be thought of as a data matrix
with column observations
>
dj (ti ) = yj1 (ti ), . . . , yjn (ti ) ∈ Rn ,

j = 1, . . . , p, i = 1, . . . , N

(3.8)

so that dj (ti ) corresponds to the response at time ti of the single output coordinate
j to each of the (separate) initial conditions x0 = ek , k = 1, . . . , n. This convention
will lead to greater clarity in the steps that follow.
6 This is not a limitation of our approach, other input signals can be used such as a white gaussian
noise, cf. [6] for preliminary results.

9

3.2. Model Order Reduction Map. The method we propose consists, in
essence, of collecting samples and then performing a process similar to “simultaneous principal components analysis” on the controllability and observability Gramian
estimates in the (same) RKHS.
As mentioned above, given a choice of the kernel K defining a RKHS H, principal
components in the feature space can be computed implicitly in the original input
space using K. It is worth emphasizing however that we will be co-diagonalizing two
Gramians in the feature space by way of a non-orthogonal transformation; the process
bears a resemblance to (K)PCA, and yet is distinct. Indeed the favorable properties
associated with an orthonormal basis are no longer available, the quantities we will in
practice diagonalize are different, and the issue of data-centering must be considered
with some additional care.
cc can be viewed as the
First note that the empirical controllability Gramian W
sample covariance of a collection of N · m vectors, scaled by T
N
N m
X
T XX j
cc = T
x (ti )xj (ti )>
X(ti )X(ti )> =
W
mN i=1
mN i=1 j=1

(3.9)

where X(t) is defined in (3.3) and the observability Gramian can be similarly viewed
as the sample covariance of a collection of N · p vectors
p
N
T XX
c
Wo =
dj (ti )dj (ti )>
pN i=1 j=1

(3.10)

where the dj are defined in Equation (3.8).
We can thus consider three quantities of interest:
• The controllability kernel matrix Kc ∈ RN m×N m of kernel products
(Kc )µν = K(xµ , xν ) = hΦ(xµ ), Φ(xν )iH

(3.11)

for µ, ν = 1, . . . , N m where we have re-indexed the set of vectors {xj (ti )}i,j =
{xµ }µ to use a single linear index.
• The observability kernel matrix Ko ∈ RN p×N p ,
(Ko )µν = K(dµ , dν ) = hΦ(dµ ), Φ(dν )iH

(3.12)

for µ, ν = 1, . . . , N p, where we have again re-indexed the set {dj (ti )}i,j =
{dµ }µ for simplicity.
• The Hankel kernel matrix Ko,c ∈ RN p×N m ,
(Ko,c )µν = K(dµ , xν ) = hΦ(dµ ), Φ(xν )iH

(3.13)

for µ = 1, . . . , N p, ν = 1, . . . , N m.
We have chosen the suggestive terminology “Hankel kernel matrix” above because
>
the square-roots of the nonzero eigenvalues of the matrix Ko,c Ko,c
are the empirical
Hankel singular values of the system mapped into feature space7 , where we assume
7 The relation between the singular values and the eigenfunctions of the Hankel operator for the
nonlinear system given by H in (2.19) and the extension of (2.11)-(2.12) to the nonlinear setting
> is an open problem that we
using the empirical Hankel singular values and eigenvectors of Ko,c Ko,c
leave for future work.

10

that the method of linear balancing can be applied. This assertion will be proved
immediately below. Note that ordinarily, N m, N p  n and Kc , Ko will be rank
deficient.
Before proceeding we consider the issue of data centering in feature space. PCA
and kernel PCA assume that the data have been centered in order to make the problem
translation invariant. In the setting considered here, we have two distinct sets of data:
the observability samples and the controllability samples. A reasonable centering
convention centers the data in each of these datasets separately. Let Ψ denote the
matrix whose columns are the observability samples mapped into feature space by
the feature map Φ, and let Φ be the matrix similarly built from the feature space
representation of the controllability samples. Then
Ko = Ψ>Ψ,

Kc = Φ>Φ,

and Ko,c = Ψ>Φ.

(3.14)

The above equation reduces to (3.5) and (3.7) in the linear case. In fact, when
Φ(x) = x, (Kc )µν = hxµ , xν iRn , (Ko )µν = hdµ , dν iRn , (Ko,c ) = hdµ , xν iRn .
Assume for the moment that there are M observability data samples and N
controllability samples, and let 1N , 1M denote the length N , M vectors of all ones,
respectively. We can define centered versions of the feature space data matrices Φ, Ψ
as
e = Φ − µc 1> ,
Φ
N

e = Ψ − µo 1>
Ψ
M

(3.15)

where µc := N −1 Φ1N and µo := M −1 Ψ1M . We will need two centered quantities in
the development below. Let us note here that, in practice, we do not need to compute
µc and µo as detailed above. Moreover, there is no need to compute the embedding
Φ since all computations are done in terms of the kernels.
The first centered quantity we consider is the centered version of Ko,c , namely
e >Φ.
e Although one cannot compute µc , µo explicitly from the data, we can
e
Ko,c = Ψ
e o,c by observing that
compute K
e o,c = Ψ − µo 1>
K
M
= Ko,c −

>

Φ − µc 1>
N

>
1
N Ko,c 1N 1N

−



>
1
M 1M 1M Ko,c

+

>
>
1
N M 1M 1M Ko,c 1N 1N .

(3.16)

The second quantity we’ll need is a centered version of the empirical observability
feature map
>
ko (x) := Ψ>Φ(x) = K(x, d1 ), . . . , K(x, dM )

(3.17)

where x ∈ Rd is the state variable and the observability samples {dj } are again indexed
by a single variable as in Equation (3.12). Centering follows reasoning similar to that
of the Hankel kernel matrix immediately above:
eo (x) = Ψ − µo 1>
k
M
= ko (x) −

>

Φ(x) − µc

1
N Ko,c 1N

−



>
1
M 1M 1M ko (x)

+

>
1
N M 1M 1M Ko,c 1N .

(3.18)

e o,c ,
Note: Throughout the remainder of this paper we will drop the special notation K
e
ko (x) and assume that Ko,c , ko (x) are centered appropriately.
With the quantities defined above, we can co-diagonalize the empirical Gramians
(balancing) and reduce the dimensionality of the state variable (truncation) in feature
11

space by carrying out calculations in the original data space. As we are assuming that
the method of linear balancing applies in the feature space, the order of the model
can be reduced by discarding small Hankel values {Σii }ni=q+1 , and projecting onto the
subspace associated with the first q < n largest eigenvalues. The following key result
describes this process:
Theorem 3.1 (Balanced Reduction in Feature Space). Consider the nonlinear
control system (3.1) and its responses (3.3) and (3.8) to impulses from the input and
initial condition, respectively. Let K be a Mercer kernel, Ko,c be the Hankel kernel
>
matrix defined in (3.13), and Ko,c Ko,c
= V Σ2 V > be its SVD decomposition with
Σ = diag{σ1 , · · · σN p } and σi ≥ σi+1 for i = 1, · · · N p. Also consider the controllability
kernel matrix Kc (3.11), the observability kernel matrix Ko (3.12).
If there is a spectral gap in Σ, i.e. there is q such that σq >> σq+1 then balanced
reduction in the RKHS can be accomplished by applying the state-space reduction map
Π : Rn → Rq given by
Π(x) = Tq>ko (x),

x ∈ Rn

(3.19)

−1/2

, Vq are the eigenvectors that correspond to the largest q Hankel
where Tq = Vq Σq
singular values. ko (x) is the empirical observability feature map (3.18).
Proof. We assume
 the data have been centered in feature space. Let Φ be a
matrix with columns Φ(xj (ti )) , i = 1, . . . , N, j = 1, . . . , m, so that
X = ΦΦ>,

(3.20)

is the feature space controllability Gramian counterpart to Equation (3.9). Similarly,
let Ψ be a matrix with columns Φ(dj (ti )) , i = 1, . . . , N, j = 1, . . . , p, so that
Y = ΨΨ>,

(3.21)

is the feature space observability Gramian counterpart to Equation (3.10). Since by
definition K(x, y) = hΦ(x), Φ(y)iF , we also have that Kc = Φ>Φ and Ko = Ψ>Ψ.
In general the Gramians X, Y are infinite dimensional whereas the kernel matrices
Kc , Ko are necessarily of finite dimension.
We now carry out linear balancing on (X, Y ) in the feature space (RKHS). First,
take the SVD of X 1/2 Ψ so that
U ΣV > = X 1/2 Ψ

(3.22)

U Σ2 U > = (X 1/2 Ψ)(X 1/2 Ψ)> = X 1/2 Y X 1/2
2

VΣ V

>

= (X

1/2

>

Ψ) (X

1/2

>

>

Ψ) = Ψ ΦΦ Ψ =

(3.23)
>
Ko,c Ko,c
.

(3.24)

The last equality in Equation (3.23) follows since X is symmetric and therefore X 1/2
is too. The linear balancing transformation is then given by
T = Σ1/2 U >X −1/2 ,

(3.25)

and one can readily verify that T XT > = T −>Y T −1 = Σ. Here, inverses should be
interpreted as peudo-inverses when appropriate8 . From Equations (3.22)-(3.24), we
8 Such as in the case of X −1/2 when the number of data points is less than the dimension of the
RKHS.

12

see that U > = Σ−1 V >Ψ>X 1/2 and thus T = Σ−1/2 V >Ψ>. We can project an arbitrary mapped data point Φ(x) onto the (balanced) “principal” subspace of dimension
q spanned by the first q rows of T by computing
Tq Φ(x) = Σ−1/2
Vq>Ψ>Φ(x) = Σ−1/2
Vq>ko (x)
q
q

(3.26)

where ko (x) := Ψ>Φ(x) is the empirical observability feature map, recalling that Vq
>
is the matrix formed by taking the top q eigenvectors of Ko,c Ko,c
by Equation (3.24).
>
We note that square roots of the non-zero eigenvalues of Ko,c Ko,c
are exactly the
Hankel singular values of the system mapped into the feature space, under the assumption of linearity in the feature space. This can be seen by noting that λ+ (Y X) =
>
λ+ (X 1/2 Y X 1/2 ) = λ+ (Ko,c Ko,c
), where λ+ (·) refers to the non-zero eigenvalues of
>
its argument. In practice, we compute the largest eigenvalues of Ko,c Ko,c
instead of
performing its SVD.
In Section 4 below we show how to use the nonlinear reduction map (3.19) to
realize a closed, reduced order system which can approximate the original system to
a high degree of accuracy.

4. Closed Dynamics of the Reduced System. Given the nonlinear state
space reduction map Π : Rn → Rq , a remaining challenge is to construct a corresponding (reduced) dynamical system on the reduced state space which well approximates
the input-output behavior of the original system on the original state space. Setting
xr = Π(x) and applying the chain rule,

(4.1)
ẋr = JΠ (x)f (x, u) x=Π† (xr )
where Π† refers to an appropriate notion (to be defined) of the inverse of Π. However
we are faced with the difficulty that the map Π is not in general surjective (even if
q = n), and moreover one cannot guarantee that an arbitrary point in the RKHS
has a non-empty preimage under Φ [28]. We propose an approximation scheme to
get around this difficulty: The dynamics f will be approximated by an element of an
RKHS defined on the reduced state space. When f is assumed to be known explicitly
it can be approximated to a high degree of accuracy. An approximate, least-squares
notion of Π† will be given to first or second order via a Taylor series expansion,
but only where it is strictly needed – and at the last possible moment – so that
a first or second order approximation will not be as crude as one might suppose.
We will also consider, as an alternative, a direct approximation of JΠ (Π† (xr )) which
takes into account further properties of the reproducing kernel as well as the fact
that the Jacobian is to be evaluated at x = Π† (xr ) in particular. In both cases, the
important ability of the map Π to capture strong nonlinearities will not be significantly
diminished.
4.1. Representation of the dynamics in RKHS. The vector-valued map
f : Rn × Rm → Rn can be approximated by composing a set of n regression functions
(one for each coordinate) fˆi : Rq×m → R in an RKHS using the representer theorem
in Appendix A, with the reduction map Π. It is reasonable to expect that this
approximation will be better than directly computing f (Π† (xr ), u) using, for instance,
a Taylor expansion approximation for Π† which may ignore important nonlinearities
at a stage where crude approximations must be avoided. Let us note here that, for this
approximation part, we are using the representer theorem as described in Appendix
13

A and that the RKHS we use is not necessarily the same as the one where we have
been performing balancing as in the previous section.
Let x̃ = Π(x) denote a reduced state variable, and concatenate the input examples
x̃j := Π(xj ) ∈ Rq , uj ∈ Rm so that zj = (x̃j , uj ) ∈ Rq×m , and {(fi (xj , uj ), zj )}`j=1 is
a set of input-output training pairs describing the i-th coordinate of the map (x̃, u) 7→
f (x, u). The training examples should characterize “typical” behaviors of the system,
and can even re-use those trajectories simulated in response to impulses for estimating
the Gramians above. We will seek the function fˆi ∈ H which minimizes
`
X

2
fˆi (zj ) − fi (xj , uj ) + λi kfˆi k2H

j=1

where λi here is a regularization parameter. We have chosen the square loss, however
other suitable loss functions may be used (cf. [53, 12, 10, 34]). As mentioned in
section A, equation (A.12), fˆi takes the form
fˆi (z) = fˆi (Π(x), u) =

`
X

cij K f (z, zj ), i = 1, . . . , n,

(4.2)

j=1

where K f defines the RKHS Hf (and is unrelated to K used to estimate the Gramians). Note that although our notation takes the RKHS for each coordinate function
to be the same, in general this need not be true: different kernels may be chosen
(i,j)=(n,`)
for each function. Here the {cij }|(i,j)=(1,1) comprise a set of coefficients found using
the regularized least squares (RLS) algorithm and satisfying the algebraic equation
(A.13) with s = {(fi (xj , uj ), zj )}`j=1 as training examples. The kernel family and any
hyper-parameters can be chosen by cross-validation. For notational convenience we
will further define the vector-valued empirical feature map


kf (x̃, u) i := K f (x̃, u), zi

f
i
for i = 1, . . . , `. In this notation fˆi Π(x), u = c>
i k (x̃, u) where (ci )j = cj .
A broad class of systems seen in thePliterature [36] are also characterized by
m
separable dynamics of the form ẋ = f (x)+ i=1 gi (x)ui . In this case, one can estimate
the functions f and gi from examples {(Π(xj ), f (xj ))}j and {(Π(xj ), g(xj ))}j .
4.2. Approximation of the Jacobian Contribution. We turn to approximating the component JΠ Π† (xr ) appearing in Equation (4.1).
4.2.1. Inverse-Taylor Expansion. A simple solution is to compute a low-order
Taylor expansion of Π and then invert it using the Moore-Penrose pseudoinverse to
obtain the approximation. For example, consider the first order expansion Π(x) ≈
Π(a) + JΠ (a)(x − a). Then we can approximate Π† (xr ) (in the first-order, least-norm
sense) as

b † (xr ) := JΠ (a) † (xr − Π(a)) + a.
Π
(4.3)
We may start with a = x0 , but periodically update the expansion in different regions
of the dynamics if desired. A good expansion point could be the estimated preimage
eo (x) is the centered version
of xr (t) returned by the algorithm proposed in [24]. If k
eo (x), then
of the length M vector ko (x) defined by (3.17), and since Π(x) = TqT k
JΠ (x) =

∂Π(x)
= Tq> I −
∂x
14

>
1
M 1M 1M

 ∂ko (x)
∂x

(4.4)

where 1M is the length M vector of all ones. JΠ (x) is a matrix in Rq×n since Tq =
−1

1
1M 1M T ) ∈
Vq Σq 2 ∈ RM ×q with Vq ∈ RM ×q , Koc ∈ RN p×N m , Σq ∈ Rq×q , (I − M

∂k
(x)
o
RM ×M , ko (x) ∈ RM ×1 , ∂x
∈ RM ×n . An example calculation of ∂x ko (x) i =
∂x K(x, di ) in the case of a polynomial kernel is given in the section immediately
below.

4.2.2. Exploiting Kernel Properties. For certain choices of the kernel K
defining the Gramian feature space H, one can exploit the fact that Kx and its derivative bear a special relationship, and potentially improve the estimate for JΠ (Π† (xr )).
Perhaps the most commonly used off-the-shelf kernel families are the polynomial and
Gaussian families. For any two kernels with hyperparameters p and q (respectively)
in one of these classes, we have that Kp = (Kq )p/q . We’ll consider the polynomial
kernel of degree d, Kd (x, y) := (1 + hx, yi)d in particular; the Gaussian case can be
derived using similar reasoning. For a polynomial kernel we have that
 d−1
∂Kd (x, y)
= dKd−1 (x, y)y > = d Kd (x, y) d y >.
∂x
Recalling that Kd (x, y) = hΦ(x), Φ(y)iH and xr = Π(x) given by (3.19), if Π were
invertible then we would have
∂Kd (x, y)
∂x

= d (Φ ◦ Π−1 )(xr ), Φ(y)
x=Π−1 (x

d−1
d

y >.

r)

The map Π is not injective however, and in addition the fibers of Φ may be potentially
empty, so we must settle for an approximation. It is reasonable then to define (Φ ◦
Π† )(xr ) as the solution to the convex optimization problem
min

kzkH

subj. to

kMq z − xr kRk = 0

z∈H

(4.5)

where Mq : H → Rk is defined as in Equation (3.26). If a point z ∈ H has a preimage in Rn this definition is consistent with composing Φ with the formal definition
Φ−1 (z) = {x ∈ Rn | Φ(x) = z} and noting that in this case Π ◦ Φ−1 = Mq (Φ ◦ Φ−1 ) =
Mq z. Furthermore, a trajectory xr (t) of the closed dynamical system on the reduced
statespace need not (and may not) have a counterpart in the original statespace by
virtue of the way in which Π† is used in our formulation of the reduction map and
corresponding reduced dynamical system.
One will recognize that the solution z ∗ to (4.5) is just the Moore-Penrose pseudoinverse z ∗ = Mq† xr . Inserting this solution into the feature map representation of
a kernel K gives the following definition for K(Π† (xr ), y):
K(Π† (xr ), y) = (Φ ◦ Π† )(xr ), Φ(y)
= Mq† xr , Φ(y)
=
=
=

H

= xr , (Mq>)† Φ(y)

H
> −1
xr , (Mq Mq ) Mq Φ(y)
xr , (Mq Mq>)−1 Π(y)
xr , (Tq>Ko Tq )−1 Π(y)

Rk

where the final equality follows applying Equations (3.22)-(3.24) and Tq is defined as
in Theorem 3.1. Substituting into the derivative for a polynomial kernel K = Kd
15

gives
∂Kd (x, y)
∂x

= d xr , (Tq>Ko Tq )−1 Π(y)

d−1
d y>

x=Π† (xr )

which immediately gives an expression for JΠ (Π† (xr )) using (4.4):
JΠ (x)|x=Π† (xr ) = Tq> I −

>
1
M 1M 1M



d xr , (Tq>Ko Tq )−1 Π(di )

d−1
d d>
i

(4.6)

Note that this approximation is global in the sense that the q × q matrix inverse
(Tq>Ko Tq )−1 need only be computed once9 ; no updating is required during simulation
of the closed system.

4.3. Reduced System Dynamics. Given an estimate fˆ Π(x), u of f (x, u) in
the RKHS Hf as in (4.2) and a notion of JΠ Π† (xr ) from (4.4) and the section above,
we can write down a closed dynamical system on the reduced statespace. We have

ẋr ≈ JΠ (x)fˆ(Π(x), u)
x=Π† (xr )

≈ JΠ (x) x=Π† (x ) C>kf (xr , u)
r
 ∂ko (x)
C>kf (xr , u)
∂x x=Π† (xr )

 > f
†
1
I−M
1M 1>
M Jk Π (xr ) C k (xr , u)

= Tq> I −
:= Tq>

>
1
M 1M 1M

(4.7)


o (x)
where C is a matrix with the vectors ci as its rows, and Jk Π† (xr ) := ∂k∂x
x=Π† (xr )
is the Jacobian of the empirical feature map defined in Equation (3.17). Here the
expression Jk Π† (xr ) should be interpreted as notation for either of the Jacobian
approximations suggested in Section 4.2.
Equation (4.7) is seen to give a closed nonlinear control system expressed solely
in terms of the reduced variable xr ∈ Rq :
(
ẋr = JΠ (x)|x=Π† (xr ) C>kf (xr , u)
(4.8)
ŷ = ĥ(xr )
where the map ĥ ◦ Π modeling the output function h : Rn → Rp is estimated as
described immediately below. Although the “true” reduced system does not actually
exist due to non-surjectivity of the feature map Φ, in many situations one can expect
that the above system will capture the essential input-output behavior of the original
system. We leave a precise analysis of the error in the approximations appearing
in (4.7) to future work10 .
9 We use the word “inverse” loosely. In practice one would use a numerically stable method, such
as an LU-factorization, which can be used to rapidly compute A−1 b for fixed A but many different
b.
10 There are different sources of error in this process of approximation. A first one due to truncation
in the RKHS which could be characterized by extending (2.15) using the eigenvalues of the Hankel
kernel matrix (3.13). There is also
 the error infinding the pseudo inverse of Φ in (4.5), i.e. that it
PN p
is tempting to write ||y − ŷ|| ≤ 2
j=q+1 σj ||u|| modulo the fact that Φ may not be surjective

and that ŷ in (4.8) is not what we will end up observing as output for the reduced order system in
Rq and the fact that there is also the error in approximating f (x, u)|x=Π† (xr ) by C>kf (xr , u) that
might be computed using results from [42] for a given u as illustrated in section 4.5 below.
16

4.4. Outputs of the Reduced System. Analogous to the case of the dynamics
f , we are faced with two possibilities for approximating

y = h Π† (xr ) .
(4.9)
We can apply
a crude Taylor series approximation to estimate Π† and therefore

†
h Π (xr ) , or as in Section 4.1, we can estimate a map
(ĥ ◦ Π) :

Rn
xr

→ Rp ,
7→ y

(4.10)

from the reduced state space to the output space directly, using
p RKHS methods.
Given samples s = {Π(xj ), yj }`j=1 , each coordinate function ĥi i=1 is given in the
familiar form (A.12)
ĥi (Π(x)) =

`
X


bij K h Π(x), Π(xj ) ,

(4.11)

j=1

where K h is the kernel chosen to define the RKHS, and may be different for each
coordinate and the bj satisfy the algebraic set of equations (A.13).
It should be noted that just given the state space reduction map Π, one can
immediately compare the output of the system defined by ĥ(xr ) to the original system
without defining a closed dynamics as above. In fact with Π and ĥ one can design a
simpler controller which takes as input the reduced state variable xr , but controls the
original system.
4.5. On Error Estimates. The problems of approximating the functions f
and h from time series can be viewed as learning problem as described in Appendix A
and a theoretical justification of our algorithm is guaranteed by the error estimates in
Theorem A.12. In fact, for the dynamical system x(k +1) = f (x(k)) (resp. the control
system x(k + 1) = f (x(k), u(k))), we have that f ∗ in (A.5) is the map f ∗ (x) = fi (x)
0
(resp. f ∗ (z) = f
i (x, u) with
 z = (x, u) ) and the samples s in (A.7) are (x(k), xi (k +
x(k)
1) + ηi ) (resp. (
, xi (k + 1) + ηi )).
u(k)


x(k)
Here fi is the unknown map x(k) → xi (k + 1) (resp.
→ xi (k + 1) ) and
u(k)
it plays the role of the unknown function (A.5).
The initial condition x(0) (resp. the input) is known and ηi are distributed
according to a probability measure ρx that satisfies the following condition (this is
the Special Assumption in [43]).
Assumption The measure ρx is the marginal on X = Rn of a Borel measure ρ
on X × R with zero mean supported on [−Mx , Mx ], Mx > 0.
Let’s note that in [43], the authors do not consider time series and that we apply
their results to time series. In the case of learning the map (4.9) ĥ ◦ Π : Rn → Rp , we
are looking at learning the map xr 7→ y from the samples s = {Π(xj ), yj }`j=1 while
in the case of learning the vector-valued map f : Rn × Rm → Rn , this corresponds to
learning the i-th coordinate of the map (x̃, u) 7→ f (x, u) for i = 1, · · · , n.
In the case of the dynamics x(k + 1) = f ∗ (x(k)) (resp. the control system x(k +
1) = f ∗ (x(k), u(k))), the error estimate (A.20) has the form
||x̂i (k + 1) − xi (k + 1)||2 ≤ 2Cx̄ Esamp + 2||x(k + 1)||2K (γ + 8Cx̄ ∆),
17

(4.12)

where ||xi (k + 1)||HK =

c2i,j
j=1 λj .

P∞


In the case of the unknown map (4.9) y = h Π† (xr ) = h̃(xr ), the error estimate
(A.20) has the form
||ŷ − y||2 ≤ 2Cx̄ Esamp + 2||h̃||2K (γ + 8Cx̄ ∆),
where ||h̃||HK =

(4.13)

c̃2i,j
j=1 λ̃j .

P∞

The first term in the right hand side of inequalities (4.12)-(4.13) represents the
error due to the noise (sampling error) and the second term represents the error due
to regularization (regularization error) and the finite-number of samples (integration
error).
4.6. Algorithm Summary. To summarize, the approach we have proposed
proceeds as follows
1. Given a nonlinear control system (2.16), let ui (t) = δ(t)ei be the i-th excitation signal for i = 1, . . . , m, and let xi (t) : t ∈ [0, ∞) 7→ xi (t) ∈ Rn
be the corresponding response of the system. Run the system and sample
the trajectories at times {tj }N
j=1 to generate a collection of N · m vectors
{xi (tj ) ∈ Rn }.
2. Fixing u(t) = 0 and setting x0 = ei for i = 1, . . . , n (separately), measure the
corresponding system output responses y i (t) : t ∈ [0, ∞) 7→ y i (t) ∈ Rp . As
before, sample the responses at times {tj }N
j=1 and save the collection of N · p
vectors {dk (ti )} defined as
>
dk (tj ) = yk1 (tj ), . . . , ykn (tj ) ∈ Rn , k = 1, . . . , p, j = 1, . . . , N
(4.14)
3. Choose a kernel K defining a RKHS H, and form the Hankel kernel matrix
Ko,c ∈ RN p×N m ,
(Ko,c )µν = K(dµ , xν ) µ = 1, . . . , N p, ν = 1, . . . , N m

(4.15)

where we have re-indexed the sets {dk (ti )} = {dµ }, {xi (tj )} = {xν } to use
single indices.
>
4. Compute the eigendecomposition Ko,c Ko,c
= V Σ2 V > assuming Ko,c has been
centered according to Equation (3.16).
5. The order of the model is reduced by discarding small eigenvalues {Σii }ni=q+1 ,
and projecting onto the subspace associated with the first q < n largest
eigenvalues. This leads to the state-space reduction map
Π : Rn →
x 7→

Rq
xr = Π(x) = Tq>ko (x)

(4.16)

−1/2

where Tq = Vq Σq
and ko (x) is the centered empirical observability feature
map given by Equation (3.18).
eo (x) = Ψ − µo 1>
k
M
= ko (x) −

>

Φ(x) − µc

1
N Ko,c 1N

−



>
1
M 1M 1M ko (x)

+

>
1
N M 1M 1M Ko,c 1N .

(4.17)

and
ko (x) = ΨT Φ(x) = (K(x, d1 ), · · · , K(x, dM ))T
18

(4.18)

6. From input/output pairs or simulated/measured trajectories, learn approximations of the dynamics and output function defined on the reduced state
space using, for instance, the representation theorem in section §A (equations
(A.12)-(A.13). The RKHS used to approximate these functions need not be
the same as the RKHS in which balanced truncation was carried out.
For the approximation of the dynamics, the representer theorem will be applied to learn the i − th coordinate of the map (Π(x), u) 7→ f (x, u) using the
samples s = (fi (xj , uj ), (Π(xj ), uj ))|`j=1 .
For the approximation of the output map, the representer theorem will be applied to learn the map xr = Π(x) 7→ y using the samples s = (Π(xj ), yj )|`j=1 .
7. Approximate the Jacobian contribution as described in section §4.2. If the
chosen kernels are polynomials, directly use (4.6).
8. Combine the approximations to determine an expression for a closed, reduced,
nonlinear dynamical system (4.8) as described in sections §. 4.3 and §. 4.4.
5. Experiments. We demonstrate an application of our method on two examples appearing in [31] (Examples 3.1. and 3.2, pgs. 52-54).
5.1. Two-Dimensional Exactly Reducible System. Consider the nonlinear
system
ẋ1
ẋ2
y

= −3x31 + x21 x2 + 2x1 x22 − x32
= 2x31 − 10x21 x2 + 10x1 x22 − 3x32 − u
= 2x1 − x2 .

(5.1)

It can be shown that this system has the same input-output relationship as the system
ẏ = −y 3 + u by rearranging terms so that
ẋ1 = −(2x1 − x2 )2 x1 + (x1 − x2 )3
ẋ2 = −(2x1 − x2 )2 x2 + 2(x1 − x2 )3 − u
y = 2x1 − x2 .
Defining the new variables z1 = 2x1 − x2 and z2 = x1 − x2 , the system can then be
re-written
ż1
ż2
y

= −z13 + u
= −z12 z2 − z23 + u
= z1 .

(5.2)

It can be seen that the variable z2 may be truncated because it doesn’t appear in the
expression of the output and thus doesn’t affect z1 . Let’s note here that the change
of variables x 7→ z is a linear one and, importantly, that the method of nonlinear
balancing as developed by Scherpen doesn’t apply for this example since the system
is not linearly controllable and the jacobian of the linearization has zero eigenvalues.
5.2. Seven-Dimensional System. We will also consider a 7-dimensional nonlinear system with one dimensional input and output:
ẋ1
ẋ2
ẋ3
ẋ4
ẋ5
ẋ6
ẋ7

=
=
=
=
=
=
=

−x31 + u
−x32 − x21 x2 + 3x1 x22 − u
−x33 + x5 + u
−x34 + x1 − x2 + x3 + 2u
x1 x2 x3 − x35 + u
x5 − x36 − x35 + 2u
−2x36 + 2x5 − x7 − x35 + 4u
19

(5.3)

with y = x1 − x22 + x3 + x4 x3 + x5 − 2x6 + 2x7 .
Here also the method of nonlinear balancing as developed by Scherpen doesn’t
apply for this example since the system is not linearly controllable and the jacobian
of the linearization has zero eigenvalues.
5.3. Experimental Setup. For both systems impulse and initial-condition responses of the system were simulated as described above, and 800 samples equally
spaced in the time interval [0, 5s] were sampled to build the Hankel kernel matrix
Ko,c given by the third degree polynomial kernel K(x, y) = (1 + hx, yi)3 . For the
2-D system we retained one component, and for the 7-D system we retained two for
the sake of variety. Thus the reduction map Π was defined by taking the top one or
two eigenvectors (scaled columns of T ) corresponding to the largest Hankel singular
values, giving a reduced state space of dimension one or two for the 2-D and 7-D
systems, respectively.
Next, a map from the reduced variable xr to ẋ was estimated following Section 4.1.
The same procedure was followed in both experiments. The control input was chosen
to be a 10hz square wave with peaks at ±1 at 50% duty cycle, and 1000 samples from
the simulated system in the interval [0, 5s] were mapped down using Π and then used
to solve the RLS regression problems, one for each state variable, again using a third
degree polynomial kernel. All initial conditions were set to zero. The desired outputs
(dependent variable examples) used to learn fˆ were taken to be the true function
f evaluated at the samples from the simulated state trajectory. We also added a
bias dimension of 1’s to the data to account for an offset, and used a fast leave-oneout cross-validation (LOOCV) computation [33] to select the optimal regularization
parameter.
We followed a similar process to learn the output function y = ĥ(xr ) for both
systems. Here we used a 10Hz square wave control input (peaks at ±2, 50% duty
cycle), zero initial conditions and 700 samples in the interval [0, 5s]. For this function
the Gaussian kernel K(x, y) = exp(−γkx − yk22 ) was used to demonstrate that our
method does not rely on any particular match between the form of the dynamics and
the type of kernel. The scale hyperparameter γ was chosen to be the reciprocal of the
average squared-distance between the training examples. We again used LOOCV to
select the RLS regularization parameter.
Finally, the comparisons between the exact and approximate reduced order systems were done using x0 = 0 and a control input different from those used
to learn

the dynamics and output functions: u(t) = 41 sin(2π3t) + sq(2π5t − π/2) where sq(·)
denotes the square wave function.
The Taylor series approximation for Π was done once, about x0 , and was not
updated further.
5.4. Results. For the 2D example, there is a clear gap in the singular valT
T
T
ues of Koc
Koc , σ1 (Koc
Koc ) = 89.3419, σ2 (Koc
Koc ) = 0.2574 while the other ones
are negligible. For the 7D example, there is also a clear gap in the singular values
T
T
T
T
of Koc
Koc , σ1 (Koc
Koc ) = 197.7821, σ2 (Koc
Koc ) = 46.9314, σ3 (Koc
Koc ) = 2.7132,
T
T
T
T
σ4 (Koc Koc ) = 0.3293, σ5 (Koc Koc ) = 0.0718, σ6 (Koc Koc ) = 0.0028, σ7 (Koc
Koc ) =
T
0.0010, σ6 (Koc Koc ) = 0.0001 while the other ones are negligible. We plot the first
100 singular values of the Hankel kernel matrix for both examples in logarithmic scale
in figure 5.3.
The simulated outputs ŷ(t) of the closed reduced systems as well as the output
y(t) of the original system together with a comparison with Lall et al. are plotted in
Figures 5.1 and 5.2, respectively. One can see that, even for a significantly different
20

0.04
True (Full order system)
Lall
True (exact reduced system)
RKHS

0.03

0.02

0.01

0

−0.01

−0.02

−0.03

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 5.1. Simulated output trajectories for the original and reduced 2-dimensional system.

input, the reduced systems closely capture the original systems. The main source
of error is seen to be over- and under-shoot near the square wave transients. This
error can be further reduced by simulating the system for different sorts of inputs
(and/or frequencies) and including the collected samples in the training sets used to
learn Π, fˆ and ĥ. Indeed, we have had some success driving example systems with
random uniform input in some cases. Finally, we note that our method is as good as
Lall et al.’s method especially that we assumed that the dynamics (3.1) is unknown
and the simulation results are obtained using (4.8) while in Lall et al. the dynamics
is assumed to be known.
6. Conclusion. We have introduced a new, empirical model reduction method
for nonlinear control systems. The method assumes that the method of linear balancing applies to nonlinear systems in a high dimensional feature space. This leads to a
nonlinear reduction map, which we suggest can be combined with representations of
the dynamics and output functions by elements of an RKHS to give a closed reduced
order dynamical system which captures the input-output characteristics of the original
system. We then demonstrated an application of our technique to a pair of nonlinear
systems and simulated the original and reduced models for comparison, showing that
the approach proposed here can yield good low-order nonlinear reductions of strongly
nonlinear control systems. We believe that techniques well known to the machine
learning and statistics communities can offer much to control and dynamical systems
research, and many further directions remain, including computing error estimates,
reduction of unstable systems, structure preserving systems, stochastic differential
equations (SDEs), and finding easily verifiable conditions of model reducibility of
nonlinear systems. For instance, we conjecture that a nonlinear system is model reducible if there exists an RKHS where there is a spectral gap in the gramians. Finally,
our results allows us to argue that working in RKHSes allows to develop methods for
21

0.3
Lall
True (Full order system)
RKHS

0.25

0.2
0.15

0.1

0.05
0

−0.05
−0.1

−0.15

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 5.2. Simulated output trajectories for the original and reduced 7-dimensional system.

a Data-Based Theory for (Nonlinear) Dynamical Systems.
7. Acknowledgements. BH thanks the European Commission and the Scientific and the Technological Research Council of Turkey (Tubitak) for financial support
received through a Marie Curie Fellowship.
Appendix A. Elements of Learning Theory.
In this section, we give a brief overview of reproducing kernel Hilbert spaces as
used in statistical learning theory. The discussion here borrows heavily from [10, 53,
46]. Early work developing the theory of RKHS was undertaken by I.J. Schoenberg
[38, 39, 40] and then N. Aronszajn [3]. Historically, RKHSes came from the question:
when is it possible to embed a metric space into a Hilbert space ?11
Definition A.1. Let H be a Hilbert space of functions on a set X . Denote by
hf, gi the inner product on H and let ||f || = hf, f i1/2 be the norm in H, for f and
g ∈ H. We say that H is a reproducing kernel Hilbert space (RKHS) if there exists
K : X × X → R such that
i. K has the reproducing property, i.e. ∀f ∈ H, f (x) = hf (·), K(·, x)i.
ii. K spans H, i.e. H = span{K(x, ·)|x ∈ X }.
K will be called a reproducing kernel of H. HK (X) will denote the RKHS H with
reproducing kernel K.
Definition A.2. Given a kernel K and inputs x1 , · · · , xn ∈ X , the n × n matrix
k := (K(xi , xj ))ij ,

(A.1)

11 A quasi-metric space (X, d) is embeddable in a Hilbert space H if there exists a mapping
Φ : X → H such that d(x, y) = ||Φ(x) − Φ(y)|| for all x and y in X. Schoenberg [38, 39, 40] proved
that a finite metric space (X, d) whose points are x0 · · · , xn can be embedded into a Hilbert space if
and only if the matrix A, whose entries are Aij = d(xi , x0 )2 + d(xj , x0 )2 − d(xi , xj )2 , is nonnegative
definite, i.e. A is a Gram matrix [8].

22

4

10

2d example
7d example
2

10

0

log(singular values of the Hankel kernel matrix

10

−2

10

−4

10

−6

10

−8

10

−10

10

−12

10

−14

10

−16

10

0

1

10

2

10

10

Fig. 5.3. First hundred singular values of the Hankel kernel matrix, in logarithmic scale, for
the 2D and 7D examples.

is called the Gram Matrix of k with respect to x1 , · · · , xn . The kernel K : X × X → R
for which for all n ∈ IN and distinct xi ∈ X gives rise to a strictly positive definite
Gram matrix is called a strictly positive definite kernel.
Definition A.3. (Mercer kernel map) A function K : X × X → R is called a
Mercer kernel if it is continuous, symmetric and positive definite.
The important properties of reproducing kernels are summarized in the following
proposition
Proposition A.4. If K is a reproducing kernel of a Hilbert space H, then
i. K(x, y) is unique.
ii. P
∀x, y ∈ X , K(x, y) = K(y, x) (symmetry).
m
iii.
i,j=1 αi αj K(xi , xj ) ≥ 0 for αi ∈ R and xi ∈ X (positive definitness).
iv. hK(x, ·), K(y, ·)iH = K(x, y).
v. Let c 6= 0. The following kernels, defined on a compact domain X ⊂ Rn , are
Mercer kernels: K(x, y) = x · y 0 (Linear), K(x, y) = (1 + x · y 0 )d , d ∈ IN
||x−y||2

(Polynomial), K(x, y) = e− σ2 , σ > 0 (Gaussian) .
Theorem A.5. Let K : X × X → R be a symmetric and positive definite
function. Then, there exists a Hilbert space of functions H defined on X admitting K
as a reproducing Kernel. Moreover, there exists a function Φ : X → H such that12
K(x, y) = hΦ(x), Φ(y)iH

for

x, y ∈ X .

Φ is called a feature map13 .
12 This

decomposition shows that kernels can be viewed as generalized dot products.
dimension of the RKHS can be infinite and corresponds to the Rdimension of the eigenspace
of the integral operator LK : L2ν (X ) → C(X ) defined as (LK f )(x) = K(x, t)f (t)dν(t) if K is a
Mercer kernel, for f ∈ L2ν (X ) and ν is a Borel measure on X .
13 The

23

Conversely, let H be a Hilbert space of functions f : X → R, with X compact, satisfying ∀x ∈ X , ∃κx > 0, such that |f (x)| ≤ κx ||f ||H . Then, H has a reproducing
kernel K.
Theorem A.6. Every sequence of functions (fn )n≥1 which converges strongly to
a function f in HK (X), converges also in the pointwise sense, that is, limn→∞ fn (x) =
f (x), for any point x ∈ X. Further, this convergence is uniform on every subset of
X on which x 7→ K(x, x) is bounded.
Remarks.
i. In theorem A.5, and using property [iv.] in Proposition A.4, we can take
Φ(x) := Kx := K(x, ·) in which case F = H – the “feature space” is the
RKHS. This is called the canonical feature map.
ii. The fact that Mercer kernels are positive definite and symmetric reminds us
of similar properties of Gramians and covariance matrices. This is an essential
fact that we are going to use in the following.
iii. In practice, we choose a Mercer kernel, such as the ones in [v.] in Proposition
A.4, and theorem A.5 guarantees the existence of a Hilbert space admitting
such a function as a reproducing kernel.
iv. Working in RKHSes allows to find nonlinear version of algorithms expressed
in terms of inner products [46]. In fact, if an algorithm contains the quantity
hx, x0 i then a nonlinear version of it, i.e. when x is replaced by φ(x) with
φ : Rn → H, would contain the quantity hφ(x), φ(x0 )i where φ represents the
nonlinearity. If we are working in an RKHS then hφ(x), φ(x0 )i := K(x, x0 )
and therefore we can replace all the quantities involving hx, x0 i in the original
algorithm by K(x, x0 ) in its nonlinear version.
/
Example. The following example is taken from [4]. Let V be the collection of
functions f with f 00 ∈ L2 [0, 1] and consider the subspace
W20 = {f (x) ∈ V : f, f 0 absolutely continuous and f (0) = f 0 (0) = 0}.
Define an inner product on W20 as
Z
hf, gi =

1

f 00 (t)g 00 (t)dt.

(A.2)

0

Using integration by parts and the Fundamental Theorem of Calculus, it can be shown
that for f ∈ W20 and any s ∈ [0, 1], f (s) can be written as
Z
f (s) =

1

(s − u)+ f 00 (u)du,

(A.3)

0

Since the reproducing kernel of the space W20 must satisfy f (s) = hf (·), R(·, s)i. From
2
= (s − u)+ .
(A.2) and (A.3), we deduce that K(·, s) is a function such that d K(u,s)
d2 u
Moreover, since K(·, s) ∈ W20 and using the property K(s, t) = hK(·, t), K(·, s)i, we
deduce that
Z 1
max(s, t) min2 (s, t) min3 (s, t)
K(s, t) = hK(·, t), K(·, s)i =
(t−u)+ (s−u)+ du =
−
2
6
0
(A.4)
4
24

RKHS play an important role in learning theory whose objective is to find an
unknown function
f∗ : X → Y

(A.5)

s = (xi , yi )|m
i=1 ,

(A.6)

from random samples

In the following we review results from [43] (for a more general setting, cf. [10])
about the special case when the data samples s are such that
Assumption 1: The samples in (A.6) have the special form
S : s = (x, yx )|x∈x̄ ,

(A.7)

∗
where x̄ = {xi }|d+1
i=1 and yx is drawn at random from f (x) + ηx , where ηx is drawn
from a probability measure ρx .
Here for each xP
∈ X, ρx is a probability measure with zero mean, and its variance
σx2 satisfies σ 2 := x∈x̄ σx2 < ∞. Let X be a closed subset of Rn and t̄ ⊂ X is a
discrete subset. Now, consider a kernel K : X × X → R and define a matrix (possibly
infinite) Kt̄,t̄ : `2 (t̄) → `2 (t̄) as
X
(Kt̄,t̄ a)s =
K(s, t)at , s ∈ t̄, a ∈ `2 (t̄),
(A.8)
t∈t̄

P
where `2 (t̄) is the set of sequences a = (at )t∈t̄ : t̄ → R with ha, bi = t∈t̄ at bt defining
an inner product. For example, we can take X = R and t̄ = {0, 1, · · · , d}.
In the case of dynamical systems such as the ones we are studying in this paper, we
are interested in learning the map F : x(k) 7→ x(k + 1) that characterises a dynamical
system x(k + 1) = F(x(k)) or the map F : (x(k), u(k)) 7→ x(k + 1) that characterises a
controlled dynamical system x(k +1) = F(x(k), u(k)). In order to get error estimates,
we could easily apply the following results to the case of the dynamical systems we
are interested in.
The case of approximating a function f ∗ ∈ HK from samples of the form (A.6) has
been studied in [43, 44]. The problem we are interested in is the one of reconstructing
f ∗ from s which can be expressed as the minimisation problem
X

f¯s,γ := argminf ∈HK,t̄
(f (x) − yx )2 + γ||f ||2K ,
(A.9)
x∈x̄

where γ ≥ 0. Moreover, in order to consider the case where x̄ is not defined by a
uniform grid on X, the authors of [43] introduced a weighting w := {wx }x∈x̄ on x̄
with wx > 014 . Let Dw be the diagonal matrix with ,aim diagonal entries {wx }x∈x̄ .
Then, ||Dw || ≤ ||w||∞ .
In this case, the regularisation scheme (A.9) becomes
X

2
2
¯
fs,γ := argminf ∈HK,t̄
wx (f (x) − yx ) + γ||f ||K ,
(A.10)
x∈x̄
14 A suggestion proposed in [43] is to consider the ρ −volume of the Voronoi associated with x̄.
X
1
Another example is w = 1 or if |x̄| = m < ∞, w = m
.

25

In learning theory, the minimization is taken over functions from a hypothesis
space often taken to be a ball of a RKHS HK associated to Mercer kernel K, and the
function fs that minimizes the empirical error Es is given by the representer theorem
(cf. [53, 12, 10] for derivation and more general forms of this theorem)
Theorem A.7. Assume f ∗ ∈ HK,t̄ and the standing hypotheses with X, K, t̄, ρ
as above, y as in (A.7). Suppose Kt̄,x̄ Dw Kx̄,t̄ + γKt̄,t̄ is invertible. Define L to be the
linear operator L = (Kt̄,x̄ Dw Kx̄,t̄ + γKt̄,t̄ )−1 Kt̄,x̄ Dw . Then the problem (A.10) has a
unique solution
X
(Ly)t Kt
(A.11)
fs,γ =
t∈t̄

The above theorem can be reformulated as follows
Theorem A.8. Let s ∈ Zm and λ ∈ R, λ > 0. The empirical target, i.e. the
function fλ,s = fs minimizing the regularized empirical error (A.9) over f ∈ HK ,
may be expressed as
fs (x) =

m
X

cj K(x, xj ),

(A.12)

j=1

where c = (c1 , · · · , cm ) is the unique solution of the well-posed linear system in Rm
λ m ci +

m
X

K(xi , xj )cj = yi ,

i = 1, · · · m,

(A.13)

j=1

Assumption 2 : For each x ∈ X,P
ρx is a probability measure with zero mean
1
supported on [−Mx , Mx ] with Bw := ( x∈x̄ wx Mx2 ) 2 < ∞.
Definition A.9. We say that x̄ is ∆−dense in X if for each y ∈ X there is
some x ∈ x̄ satisfying ||x − y||`∞ (Rn ) ≤ ∆.
Theorem A.10. (Sample Error) [Theorem 4, Propositions 2 and 3 in[43]]
Suppose
Kt̄,x̄ Dw Kx̄,t̄ + γKt̄,t̄ is invertible. Under the assumption (A.7), let fs,γ =
P
t∈t̄ ct Kt be the solution of (A.10) given in Theorem A.7 by c = Ly. Let Lw and κ
be
1/2
Lw = (Kt̄,x̄ Dw Kx̄,t̄ + γKt̄,t̄ )−1 Kt̄,x̄ Dw

(A.14)

κ := ||Kt̄,t̄ || ||(Kt̄,x̄ Dw Kx̄,t̄ + γKt̄,t̄ )−1 ||2

(A.15)

Then for every 0 < δ < 1, with confidence 1 − δ we have the sample error estimate


2
1
2 −1 2||Kt̄,t̄ Lw || ||Lw || Bw
||fs,γ − fx̄,γ ||2K ≤ Esamp := κσw
α
log
,
(A.16)
2
κσw
δ
where α is the increasing function defined for u > 1 as α(u) = (u − 1) log u. In
2
particular, Esamp → 0 when γ → ∞ or σw
→ 0.
Theorem A.11. (Regularization Error and Integration Error, Proposition 4 and
Theorem 5 in[43])
Under Assumptions 1 and 2. Let X̄ = (Xx )x∈x̄ be the Voronoi of X associated with x̄ and wx = ρX (Xx ). Define the Lipschitz norm on a subset X 0 ⊂ X as
26

||f ||Lip(X 0 ) := ||f ||L∞ (X 0 ) + sups,u∈X

|f (s)−f (u)|
||s−u||`∞ (Rn )
15

and assume that the inclusion map

of HK,t̄ into the Lipschitz space satisfies
P
Cx̄ := sup

x∈x̄

f ∈HK,t̄

wx ||f ||2Lip(X

x)

< ∞.

(A.17)

γ||Kt̄,t̄ || ||f ∗ ||2K
λ2x̄,w

(A.18)

||f ||2K

i. If f ∗ ∈ HK,t̄ and λx̄,w > 0, then
||fx̄,γ − f ∗ ||2K ≤

ii. If x̄ is ∆−dense, Cx̄ < ∞, and f ∗ ∈ HK,t̄ , then
||fx̄,γ − f ∗ ||2 ≤ ||f ∗ ||2K (γ + 8Cx̄ ∆)

(A.19)

Theorem A.12. (Sample, Regularization and Integration Errors) (Corollary 5
in [43])
Under Assumptions 1 and 2. Let X̄ = (Xx )x∈x̄ be the Voronoi of X associated
with x̄ and wx = ρX (Xx ). If x̄ is ∆−dense, Cx̄ < ∞, and f ∗ ∈ HK,t̄ , then, for every
0 < δ < 1, with probability at least 1 − δ there holds
||fs,γ − f ∗ ||2 ≤ 2Cx̄ Esamp + 2||f ∗ ||2K (γ + 8Cx̄ ∆),

(A.20)

where Esamp is given in (A.16).
If f ∗ is not an element of HK,t̄ then one also needs an estimate for the approximation error [42],[10].
2
Theorem A.13. Define fs,γ by (A.11). If L−r
K fρ ∈ LρX , then
||fs,γ − fρ ||L2ρ
When

1
2

X

≤ λr ||L−r
K fρ ||L2ρ ,

if

0 < r ≤ 1.

(A.21)

X

< r ≤ 1, we have
1

||fs,γ − fρ ||K ≤ λr− 2 ||L−r
K fρ ||L2ρ

(A.22)
X

Here fs,γ is taken as an approximation of the regression function fρ . Hence, minimizing over the (possibly infinite dimensional) Hilbert space, reduces to minimizing
over Rm . The series (A.11) converges absolutely and uniformly to f . We call learning
the process of approximating the unknown function f from random samples on Z.
In the following,
we assume that the kernels K are continuous and bounded by
p
κ = supx∈X K(x, x) < ∞.
Appendix B. Kernel PCA. Kernel PCA [41] will be a helpful starting point
for understanding the approach to balanced reduction introduced in this paper. We
briefly review the relevant background here.
Kernel PCA (KPCA) is a generalization of linear PCA that allows to take into
account nonlinear versions of observations. This is done by carrying out PCA in a
high dimensional RKHS through an injective, not necessarily surjective, map
Φ : Rn → H, x 7→ Φ(x).
15 This

(B.1)

assumption is true if X is compact and the inclusion map of HK,t̄ into the space of
Lipschitz functions on X is bounded which is the case when K is a C 2 Mercer kernel [55]. In fact,
if ||f ||Lip(X) ≤ C0 ||f ||K for each f ∈ HK,t̄ , then Cx̄ ≤ C02 ρX (X).
27

Φ is possibly nonlinear and H is a high-dimensional, possibly infinite-dimensional,
n
n
n
RKHS. Given the set of data x := {xi }N
i=1 ⊂ R , the covariance matrix C ∈ R × R
defined as
C=

N
1 X
xi xTi
N i=1

(B.2)

becomes a covariance matrix in H
CH =

N
1 X
Φ(xi )Φ(xi )T
N i=1

(B.3)

If H is infinite-dimensional, we think of Φ(xi )Φ(xi )T as a linear operator on H,
mapping x 7→PΦ(xj )hΦ(xj ), xi. We will assume the data are centered in the feature
space so that i Φ(xi ) = 0. If not, data may be centered according to the prescription
in [41].
Taking the feature map16 Φ : Rn → RN with Φi (x) = K(x, xi ) and given the set
n
N
of data x := {xi }N
by simply working with the
i=1 ⊂ R , we can consider PCA in R
covariance of the mapped vectors (B.3).
The principal subspaces are computed by diagonalizing CH through solving the
eigenvalue problem
λv = CH v,

(B.4)

for eigenvalues λ ≥ 0 and nonzero eigenvectors v ∈ RN . This problem can be expressed
in term of an inner-product by plugging (B.3) into (B.4) and getting
λv = CH v =

N
1 X
hΦ(xi ), viH Φ(xi ),
N i=1

(B.5)

However as is shown in [41], we can perform the computations directly in terms
of the kernel without explicitly knowing Φ. One can equivalently form the matrix K
of kernel products whose entries are (K)ij = K(xi , xj ) for i, j = 1, . . . , N , and solve
the eigenproblem
Kα = N λα,

(B.6)

in RN . Moreover the eigenvectors v in (B.4) and the eigenvectors α in (B.6) are
related through
vi = Ψαi ,


(B.7)

where Ψ := Φ(x1 ) · · · Φ(xM ) , and the non-zero eigenvalues of K and CH coincide.
The eigenvectors αi of K are subsequently normalized so that the eigenvectors vi
2
of CH have unit norm in the RKHS, leading to the condition kαi k = λ−1
i . Assuming
this normalization convention, sort the eigenvectors according to the magnitudes of
the corresponding eigenvalues in descending order, and form the matrix


Aq = α1 · · · αq , 1 ≤ q ≤ min(n, M ).
(B.8)
16 This

feature map is valid given property iv. in Proposition A.3.
28



Similarly, form the matrix Vq = v1 · · · vq , 1 ≤ q ≤ min(n, M ) of sorted eigenvectors of CH . The first q principal components of a vector x = Φ(x̃) in the feature
space are then given by Vq>x. It can be shown however (see [41]) that principal components in the feature space can be computed in the original space with kernels using
the map Π : RM → Rq
Π(x) := A>
(B.9)
q k(x),
>
where k(x) = K(x, x1 ), . . . , K(x, xM ) .
Kernel methods can be used to develop nonlinear generalizations of any algorithm
that can be expressed in terms of inner products [41] and KPCA is an illustration of
this approach. KPCA is viewed as a nonlinear version of PCA since PCA in Rn can
be reformulated as an eigenvalue problem in terms of inner products as in (B.5) but
with CH replaced by C and Φ(x) = x, i.e.
λv = Cv =

M
1 X
hxi , vixi ,
M i=1

(B.10)

given the expression of C in (B.2). If one wants to perform PCA on a nonlinear
version of the data (xi )|M
i=1 through a nonlinear map Φ, it is enough to replace x by
Φ(x) in the eigenvalue problem (B.10) to get (B.5).
Our goal in this paper is to extend this method to linear balancing in view of
applying it to nonlinear control systems.
REFERENCES
[1] Antoulas, A. C. (2005). Approximation of Large-Scale Dynamical Systems, SIAM Publications.
[2] C Archambeau, D Cornford, M Opper, J Shawe-Taylor, Gaussian process approximations of
stochastic differential equations, Journal of machine learning research, pp. 1-16, 64, 2007.
[3] Aronszajn, N. (1950). Theory of Reproducing Kernels. Trans. Amer. Math. Soc., 68:337-404.
[4] Berlinet, A. and C. Thomas-Agnan (2004). Reproducing Kernel Hilbert Spaces in Probability
and Statistics, Kluwer Academic Publishers, Norwell, MA.
[5] Bouvrie, J. and B. Hamzi (2010). Balanced Reduction of Nonlinear Control Systems in Reproducing Kernel Hilbert Space, in Proc. 48th Annual Allerton Conference on Communication, Control, and Computing, pp. 294-301. http://arxiv.org/abs/1011.2952.
[6] Bouvrie, J. and B. Hamzi (2012), Empirical Estimators for the Controllability Energy and Invariant Measure of Stochastically Forced Nonlinear Systems,in Proc. of the 2012 American
Control Conference (long version at http://arxiv.org/abs/1204.0563).
[7] Bouvrie, J. and B. Hamzi (2014), Embedology for Control and Random Dynamical Systems in
Reproducing Kernel Hilbert Spaces, in preparation.
[8] Cheney, W. and W. Light (2009). A Course in Approximation Theory, Graduate Studies in
Mathematics, vol. 101, AMS.
[9] Coifman, R. R., I. G. Kevrekidis, S. Lafon, M. Maggioni and B. Nadler (2008). Diffusion
Maps, Reduction Coordinates, and Low Dimensional Representation of Stochastic Systems,
Multiscale Model. Simul., 7(2):842-864.
[10] Cucker, F. and S. Smale (2001). On the mathematical foundations of learning. Bulletin of AMS,
39:1-49.
[11] Dullerud, G. E., and F. Paganini (2000). A Course in Robust Control Theory: a Convex
Approach, Springer.
[12] Evgeniou, T., M. Pontil, and T. Poggio (2000). Regularization networks and support vector
machines, Advances in Computational Mathematics, vol. 13, no. 1, pp. 1-50.
[13] Fujimoto, K. and D. Tsubakino (2008). Computation of nonlinear balanced realization and
model reduction based on Taylor series expansion, Systems and Control Letters, 57, 4, pp.
283-289.
[14] Fujimoto, K. and J. Scherpen (2010). Balanced realization and model order reduction for nonlinear systems based on singular value analysis, SIAM Journal on Control and Optimization,
Vol. 48, 7, pp. 180-194.
29

[15] Gray, W. S. and E. I. Verriest (2006). Algebraically Defined Gramians for Nonlinear Systems,
Proc. of the 45th IEEE CDC.
[16] Hammarling, S. J. (1982). Numerical Solution of the Stable, Non-negative Definite Lyapunov
Equation, IMA Journal of Numerical Analysis, vol. 2, pp. 303-323.
[17] Jolliffe, I.T. (2002). Principal Component Analysis, Springer.
[18] Jonckheere, E.A., and L. M. Silverman (1983). A New Set of Invariants for Linear Systems
- Application to Reduced Order Compensator Design. IEEE Transactions on Automatic
Control, AC-28, 10, 953-964.
[19] Kantz, H. and T. Schreiber. Nonlinear Time Series Analysis, Cambridge University Press, 2004.
[20] Kenney, C. and G. Hewer (1987). Necessary and Sufficient Conditions for Balancing Unstable
Systems, IEEE Transactions on Automatic Control, 32, 2, pp. 157-160.
[21] Krener, A. (2006). Model reduction for linear and nonlinear control systems. Bode Lecture,
45th IEEE Conference on Decision and Control.
[22] Krener, A. J. (2007). The Important State Coordinates of a Nonlinear System. In “Advances
in control theory and applications”, C. Bonivento, A. Isidori, L. Marconi, C. Rossi, editors,
pp. 161-170. Springer.
[23] Krener, A. J. (2008). Reduced order modeling of nonlinear control systems. In “Analysis and
Design of Nonlinear Control Systems”, A. Astolfi and L. Marconi, editors, pp. 41-62.
Springer.
[24] Kwok, J. T. and I.W. Tsang (2003). “The Pre-Image Problem in Kernel Methods”. In Proceedings of the Twentieth International Conference on Machine Learning (ICML).
[25] Lall, S., J. Marsden, and S. Glavaski (2002). A subspace approach to balanced truncation
for model reduction of nonlinear control systems, International Journal on Robust and
Nonlinear Control, 12, 5, pp. 519-535.
[26] Laub, A.J. (1980). On Computing “balancing” transformations, Proc. of the 1980 Joint Automatic Control Conference (ACC).
[27] Li, J.-R. (2000). Model Reduction of Large Linear Systems via Low Rank System Grammians.
Ph.D. thesis, Massachusetts Institute of Technology.
[28] Mika, S., B. Schölkopf, A. Smola, K. R. Müller, M. Scholz, and G. Rätsch (1998). Kernel PCA
and de-noising in feature spaces, In Proc. Advances in Neural Information Processing
Systems (NIPS) 11, pp. 536–542, MIT Press.
[29] Moore, B. (1981). Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduction, IEEE Tran. Automat. Control, 26, 1, pp. 17-32.
[30] Newman, A.J., and P. S. Krishnaprasad (2000). Computing balanced realizations for nonlinear
systems, Proc. of the Math. Theory of Networks and Systems (MTNS).
[31] Nilsson, O. (2009). On Modeling and Nonlinear Model Reduction in Automotive Systems, Ph.D.
thesis, Lund University.
[32] Phillips, J., J. Afonso, A. Oliveira and L. M. Silveira (2003). Analog Macromodeling using
Kernel Methods. In Proceedings of the IEEE/ACM International Conference on Computeraided Design.
[33] Rifkin, R., and R.A. Lippert. Notes on Regularized Least-Squares, CBCL Paper 268/AI Technical Report 2007-019, Massachusetts Institute of Technology, Cambridge, MA, May, 2007.
[34] Schaback, R. and H. Wendland. Kernel techniques: From machine learning to meshless methods, Acta Numerica, vol.15, pp. 543-639, 2006.
[35] Scherpen, J. (1993). Balancing for nonlinear systems. Syst. & Contr. Let., 21, pp 143-153.
[36] Scherpen, J. (1994). Balancing for Nonlinear Systems, Ph.D. thesis, University of Twente,
http://www.dcsc.tudelft.nl/~jscherpen/thesis.html.
[37] Scherpen, J. (2011). Balanced Realizations, Model Order Reduction, and the Hankel Operator,
the Control Systems Handbook, 2nd ed., Advanced Methods, Eds. William Levine, CRC
Press, Taylor & Francis Group, Chap. 4, pp. 1-24.
[38] Schoenberg, I. J. (1935). Remarks to Maurice Fréchet’s article ”Sur la définition axiomatique
d’une classe d’espace distanciés vectoriellement applicable sur l’espace de Hilbert”, Annals
of Mathematics 36 (1935), 724–732.
[39] Schoenberg, I. J. (1937) On certain metric spaces arising from euclidean spaces by a change of
metric and their imbedding in Hilbert space, Annals of Mathematics, (2), vol. 38 (1937),
pp. 787-793.
[40] Schoenberg, I. J. (1938). Metric spaces and positive definite functions, Trans. Amer. Math. Soc.
44 (1938), pp. 522-536.
[41] Schölkopf, B., Smola, A., and Müller, K (1998). Nonlinear component analysis as a kernel
eigenvalue problem. Neural Computation, 10(5):1299–1319.
[42] Smale, S. and D.,-X. Zhou (2003). Estimating the Approximation Error in Learning Theory,
Analysis and Applications, Volume 01, Issue 01.
30

[43] Smale, S. and D.,-X. Zhou, Shannon sampling and function reconstruction from point values,
Bull. Amer. Math. Soc. 41, 2004, pp. 279-305.
[44] Smale, S. and D.,-X. Zhou, Shannon sampling II: Connections to learning theory, Applied and
Computational Harmonic Analysis, Volume 19, Issue 3, November 2005, pp. 285-302.
[45] S. Smale and D.-X. Zhou, Online Learning with Markov Sampling, Anal. Appl., 7, pp. 87-113,
2009.
[46] Schölkopf, B., and A. J. Smola (2002). Learning with Kernels, The MIT Press.
[47] Steinwart, I. and A. Christmann (2008). Support Vector Machines, Springer.
[48] Therapos, C. P. (1989). Balancing Transformations for Unstable Nonminimal Linear Systems,
IEEE Transactions on Automatic Control, 34, 4, pp. 455-457.
[49] Verriest, E. (1981). Suboptimal LQG-Design via Balanced Realizations. Proc. of the 20th IEEE
CDC, pp. 686-687.
[50] Verriest, E. (1984). Approximation and Order Reduction in Nonlinear models using an RKHS
Approach, Proc. of the 18th Annual Conference on Information Sciences and Systems, pp.
197-201.
[51] Box, G. E. P. and G. M. Jenkins, Time Series Analysis: Forecasting and Control, Wiley, 2008.
[52] L. Ljung, System Identification: Theory for the User, Prentice Hall, 1999.
[53] Wahba, G. (1990). Spline Models for Observational Data, SIAM CBMS-NSF Regional Conference Series in Applied Mathematics 59.
[54] Weiland, S. (1991). Theory of approximation and disturbance attenuation of linear systems,
Doctoral dissertation, University of Groningen.
[55] Zhou, D.-X. (2003). Capacity of reproducing kernel spaces in learning theory, IEEE Transactions on Information Theory, Volume:49 , Issue: 7, pp. 1743 - 1752.

31

