Approximation of Lyapunov Functions from Noisy Data

0.29 MB·

by  Peter Giesl

Approximation of Lyapunov Functions from Noisy Data
P. Giesl, B. Hamzi, M. Rasmussen, K. N. Webster
∗ † ‡ §
January 8, 2016
6 Abstract
1
0 Methods havepreviously beendeveloped fortheapproximation ofLyapunov functions using
2
radial basis functions. However these methods assume that the evolution equations are known.

n
We consider the problem of approximating a given Lyapunov function using radial basis func-
a
J tions where the evolution equations are not known, but we instead have sampled data which is

7 contaminated withnoise. Wepropose analgorithm inwhichwefirstapproximate theunderlying

vector field, and use this approximation to then approximate the Lyapunov function. Our ap-
]
S proachcombineselementsofmachinelearning/statistical learningtheorywiththeexistingtheory
D
ofLyapunovfunction approximation. Errorestimatesareprovidedforouralgorithm.
.
h
t
a
m
1 Introduction
[

1
Ordinary differential equations model large classes of applications such as planetary motion, chemi-
v
8 cal reactions, population dynamics or consumer behaviour. A breakthrough in the understanding of
6
ordinary differential equations was initiated by Poincare´ and Lyapunovin the late 19th century, who
5
1 developedanapproachthatembracedtheuseoftopologicalandgeometricaltechniquesforthestudy
0
of dynamical systems. A key component of this theory is Lyapunov functions, which can be used to
.
1 determinethebasinofattractionofan asymptoticallystableequilibrium.
0
6 In general, it is not possible to find an explicit analytical expression for a Lyapunov function
1 associated to a nonlinear differential equation. Many methods have been proposed to numerically
:
v construct Lyapunov functions, see [13] for a recent review. These methods include the SOS (sums
i
X of squares) method, which constructs a polynomial Lyapunov function by semidefinite optimization
r [26]. Anothermethodconstructsacontinuouspiecewiseaffine(CPA)Lyapunovfunctionusinglinear
a
optimization [16]. A further method is based on Zubov’s equation and computes a solution of this
partial differentialequation[6]. Lyapunovfunctionscan also beconstructedusingset orientedmeth-
ods [15]. Themethod that is also used in this paper is based on approximatingthe solutionofa PDE
∗DepartmentofMathematics,UniversityofSussex,Falmer,BN19QH,UK([email protected])
†Department of Mathematics, Koc¸ University, Istanbul, Turkey & Departmentof Mathematics, Alfaisal University,
Riyadh,SaudiArabia([email protected])
‡DepartmentofMathematics,ImperialCollegeLondon,SW72AZ,UK([email protected])
§Department of Mathematics, Imperial College London, SW7 2AZ, UK & Potsdam Institute for Climate Impact
Research,POBox601203,14412Potsdam,Germany([email protected]).
1
ApproximationofLyapunovfunctionsfromnoisydata 2
using radial basis functions [11]. All these methods to approximate Lyapunov functions rely on the
knowledgeoftherighthand sideofthedifferentialequation.
In this paper, we develop a method to approximate Lyapunov functions where the right hand
side is unknown, but we have sampled data of the system, which is contaminated by noise. We will
first approximate the right hand side of the differential equation, and then use this approximation
to approximate the Lyapunov function. Our approach combines and develops previous results from
statistical learning theory [28, 29, 30] together with existing methods using radial basis functions
[11, 14], which usetheframeworkofreproducing kernelHilbert spaces (RKHS).
TheuseofRKHSspacestoapproximateimportantquantitiesindynamicalsystemshaspreviously
been exploitedby Smale and Zhou to approximatea hyperbolicdynamicalsystem[31]. Bouvrie and
Hamzi also use RKHS spaces to approximate some key quantities in control and random dynamical
systems[4, 5].
2 Setting of the Problem and Main Result
Weconsiderordinarydifferentialequationsoftheform
x˙ = f∗(x), (2.1)
wheref∗ : Rd Rd isasmoothvectorfield anddotdenotesdifferentiationwithrespect totime. We
→
definetheflow ϕ : Rd R Rd byϕ (ψ,t) := x(t), wherex(t) solves(2.1)withx(0) = ψ.
f∗ f∗
× →
We assume that (2.1) has a fixed point x that is exponentially asymptotically stable. Define the
basin of attraction as A(x) := ψ Rd lim ϕ (ψ,t) = x . Note that A(x) = and A(x) is
t→∞ f∗
{ ∈ | } 6 ∅
open. Subsets of the basin of attraction can be determined by the use of Lyapunov functions, which
are functions decreasing along solutions of (2.1). We consider two types of Lyapunov functions V
and T,as described inTheorems 3.2and 3.3below. TheseLyapunovfunctionssatisfy
V(x),f∗(x) = p(x), x A(x),
Rd
h∇ i − ∈
T(x),f∗(x) = c, x A(x) x ,
Rd
h∇ i − ∈ \{ }
where p is a smooth function with p(x) > 0 for x = x¯ and p(x¯) = 0, and c is a positive constant.
6
ThescalarproductsonthelefthandsidesarecalledtheorbitalderivativesofV andT withrespectto
(2.1), which are the derivativesof V and T along solutionsof (2.1). The orbital derivativesof V and
T arenegative,which impliesthat V and T are decreasing alongsolutions.
We assume that the function f∗ is unknown, but we have sampled data of the form (x ,y ) in
i i
X Rd,i = 1,...,m,withy = f∗(x )+η . Weassumethattheone-dimensionalrandomvariables
× i i xi
ηk Rd, where i = 1,...,m and k = 1,...,d, are independent random variables drawn from a
xi ∈
probability distribution with zero mean and variance (σk )2 bounded by σ2. Here X is a nonempty
xi
and compactsubsetofRd withC1 boundary.
In 5 we provide an algorithm to approximately reconstruct the Lyapunov functions V and T
by func§tions Vˆ and Tˆ. The following main theorem provides error estimates in a compact set
D ⊂
A(x) X,whichdependonthedensityofthedata,measuredbytwokeyquantities: thefilldistanceof
∩
the data h (see Definition 6.3) and the norm of the volumeweights w corresponding to theVoronoi
x
tessellationofthedata(seeDefinition5.1).
ApproximationofLyapunovfunctionsfromnoisydata 3
Theorem 2.1 Consider (2.1) such that f∗ Cν1(Rd,Rd) with ν (3d + 7)/2 if d is odd, or
1
∈ ≥
ν (3d + 12)/2 if d is even. Let τ ,τ R and k ,k N be such that τ = k + (d + 1)/2
1 1 2 1 2 1 1
≥ ∈ ∈
with τ = ν , and k = k (d + 2) (if d is odd) or k = k (d + 3) (if d is even). Define
1 1 2 1 2 1
⌈ ⌉ − −
τ := k +(d+1)/2.
2 2
Let Ω A(x) be a compact set and := Ω B (x) X, with ε > 0 small enough so that
ε
⊂ D \ ⊂
= . Forhx, w Rm and hq sufficientlysmall,thefollowingholds:
D 6 ∅ || ||
1. For every 0 < δ < 1, the reconstructionVˆ of theLyapunov functionV defined in Theorem 3.2
satisfiesthefollowingestimatewith probability1 δ:
−
h∇Vˆ,f∗iRd −h∇V,f∗iRd L∞(D) ≤ C||V||W2τ2(ΩV) hqk2−21
(cid:12)(cid:12) (cid:12)(cid:12) (cid:16)
(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) +||w||Rm +λr−23hx +λr−21 , (2.2)
λ√δ
(cid:19)
where Ω is acertaincompact subsetofA(x), and 1 < r 1.
V ⊃ D 2 ≤
2. For every 0 < δ < 1, the reconstruction Tˆ of the Lyapunov function T defined in Theorem 3.3
satisfiesthefollowingestimateswith probability1 δ:
−
h∇Tˆ,f∗iRd −h∇T,f∗iRd L∞(D) ≤ C||T||W2τ2(ΩT) hqk2−21
(cid:12)(cid:12) (cid:12)(cid:12) (cid:16)
(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) +||w||Rm +λr−32hx +λr−21 , (2.3)
λ√δ
(cid:19)
Tˆ −T L∞(Γ) ≤ Chqk˜2+21||T||W2τ2(ΩT), (2.4)
(cid:12)(cid:12) (cid:12)(cid:12)
where Γ is a non-characterist(cid:12)ic(cid:12) hyper(cid:12)s(cid:12)urface on which T has defined values (see Definition
(cid:12)(cid:12) (cid:12)(cid:12)
3.4), Ω is acertaincompact subsetofA(x),and 1 < r 1.
T ⊃ D 2 ≤
Themainpointisthattheexpressionsontherighthandsideof(2.2)–(2.4)canbemadearbitrarily
small as the data density increases and for suitably chosen λ (see equation (6.13)). Therefore the
orbital derivative of our Lyapunov function approximations Vˆ and Tˆ become arbitrarily close in the
infinity norm to those of V and T respectively. Estimate (2.2) implies that the orbital derivative of
Vˆ will be negativein (which does not contain a small neighbourhood of the equilibrium x), since
D
V(x),f∗(x) = p(x)wherepisapositivedefinitefunction(seeTheorem3.2). Theanalogous
Rd
h∇ i −
statement is true for Tˆ, since T(x),f∗(x) = c < 0. In principle the neighbourhood B (x)
Rd ε
h∇ i −
can shrinkas thedatadensityincreases (as hx and w Rm tend tozero).
|| ||
The above estimate contains λ > 0 as a regularisation parameter of our algorithm, and h as the
q
fill distance of a set of sampled points in Ω (resp. Ω ) of our choosing. Similarly, h is the fill
V T q˜
distance of a set of sampled points on Γ, which we are able to choose. The constants in the above
estimatesdependon d,σ, thechoiceoffunctionspaces forapproximationandthevectorfield f∗.
The rest of the paper is organised as follows. In 3 we provide the converse theorems for the
§
LyapunovfunctionsV and T. In 4 we setout theframeworkforthefunctionspaces that are usedto
§
approximate the Lyapunov functions, as well as previous results on the approximation of Lyapunov
functions when the right hand side of (2.1) is known. The algorithms themselves that are used to
computeVˆ and Tˆ are detailed in 5. In 6 we providean estimatefor our approximationof the right
§ §
hand sideof(2.1), whichisthen usedin theproofofTheorem2.1 in 7.
§
ApproximationofLyapunovfunctionsfromnoisydata 4
ΩV/ΩT
X
Γ
Γ˜
Ω
B (x)
ε x
Figure 1: Domains and sets used in the statement and proof of Theorem 2.1. The dotted lines show
the boundary of the set = Ω B (x) where the Lyapunov functions are approximated. We also
ε
D \
haveΩ ,Ω A(x).
V T
⊂
3 Converse theorems for Lyapunov functions
The concept of a Lyapunov function dates back to 1893, where Lyapunov introduced these func-
tions for the stabilityanalysis of an equilibriumfor a given differential equation, without theexplicit
knowledge of the solutions [19]. Many converse theorems have been proved that guarantee the ex-
istence of a Lyapunov function under certain conditions, see [11, 18, 21] for an overview. Massera
[22] provided the first main converse theorem for C1 vector fields where A(x) = Rd, with further
developments by several authors to prove the existence of smooth Lyapunov functions under weak
smoothnessassumptionson therighthandside(seee.g. [17, 20, 33]).
The existence of a Lyapunov function for system (2.1) with given values of the orbital derivative
has been shown by Bhatia [2, 3], as stated in thefollowingtheorems (see also [11]). We also refer to
[12] for a proof that the conditionson the function p givenhere are sufficient to define the Lyapunov
functionV, in contrasttotheconditionsgivenin [11]. Firstwemakethefollowingdefinition.
Definition 3.1 A continuousfunctionα : [0, ) [0, ) is a class functionifα(0) = 0 and α is
∞ → ∞ K
strictlymonotonicallyincreasing.
Theorem 3.2 Consider the autonomous system of differential equations x˙ = f∗(x), where f∗
∈
Cν1(Rd,Rd), ν 1, d N. We assume the system to have an exponentially asymptotically stable
1
≥ ∈
equilibriumx. Let p Cν1(Rd,R)bea functionwith thefollowingproperties:
∈
1. p(x) > 0 forx = x, andp(x) = 0.
6
2. Thereis a class functionα suchthatp(x x) α( x x )for allx Rd.
2
K − ≥ || − || ∈
Then there exists a Lyapunov function V Cν1(A(x),R) (where A(x) is the basin of attraction
∈
ofx), suchthat
V(x),f∗(x) = p(x) (3.1)
Rd
h∇ i −
holdsforallx A(x). The LyapunovfunctionV isuniquelydefined upto aconstant.
∈
ApproximationofLyapunovfunctionsfromnoisydata 5
Wemayalsochoosep(x)tobeapositiveconstantinequation(3.1)toobtainaLyapunovfunction
T,defined onA(x) x ,forwhich lim T(x) = .
x→x
\{ } −∞
Theorem 3.3 Consider the autonomous system of differential equations x˙ = f∗(x), where f∗
∈
Cν1(Rd,Rd), ν 1, d N. We assume the system to have an exponentially asymptotically stable
1
≥ ∈
equilibrium x. Then for all c R+, there exists a Lyapunov function T Cν1(A(x) x ,R) such
∈ ∈ \{ }
that
T(x),f∗(x) = c. (3.2)
Rd
h∇ i −
Moreover, lim T(x) = .
x→x
−∞
The Lyapunov function T will be uniquely defined if its values are given on a non-characteristic
hypersurfaceΓ A(x) [11]bya functionξ Cν1(Γ,R);that is,T(x) = ξ (x) forx Γ.
T T
⊂ ∈ ∈
Definition 3.4(Non-characteristic hypersurface) Consider x˙ = f∗(x), where f∗ Cν1(Rd,Rd),
∈
ν 1,d N. Leth Cν1(Rd,R)andrecallϕ : Rd R Rd denotestheflowmapping. Theset
1 f∗
≥ ∈ ∈ × →
Γ Rd iscalled a non-characteristichypersurfaceif
⊂
1. Γ iscompact,
2. h(x) = 0 ifand onlyif x Γ,
∈
3. h′(x) = h(x),f∗(x) < 0 holdsforallx Γ, and
Rd
h i ∈
4. foreachx A(x) x thereisa timeθ(x) Rsuch thatϕ (x,θ(x)) Γ.
f∗
∈ \{ } ∈ ∈
An example of a non-characteristic hypersurface is the level set of a (local) Lyapunov function,
see[11]. Inwhatfollowsweassumethatwehavechosenanon-characteristichypersurfaceΓtogether
withafunctionξ Cν1(Γ,R)so thatT isuniquelydefined.
T
∈
Remark3.5 The Lyapunov function V is smooth and defined on the entire basin of attractionA(x).
However, its orbital derivative vanishes at the equilibrium x, and therefore estimates that bound the
error of the orbital derivative for numerical approximations of V cannot guarantee negative orbital
derivative arbitrarilyclose to the equilibrium x. On the other hand, the Lyapunov function T is not
even defined at x and is unbounded near the equilibrium. However, its definition has the advantage
thatitis notrequiredthatwe knowwhere theequilibriumis.
InourapproachtoapproximateLyapunovfunctionsdirectlyfromdata,wewillprovideestimates
fortheapproximationofbothV andT. Thestrategyistofirst approximatef∗ fromthedata, anduse
thisinturn toapproximatetheLyapunovfunction.
4 Function spaces and approximation theorems for Lyapunov
functions
4.1 Reproducing kernel Hilbert spaces
The function spaces that we use to search for our approximationsto both f∗ and the Lyapunovfunc-
tionsV andT willbereproducingkernelHilbertspaces(RKHS). Forasurveyofthemainproperties
ofRKHSspaces mentionedinthissection,werefer to [7].
ApproximationofLyapunovfunctionsfromnoisydata 6
In order to define an RKHS function space we first fix a continuous, symmetric, positivedefinite
function (a “kernel”) K : X X R, and set K := K( ,x). Define theHilbert space by first
x K
× → · H
considering all finite linear combinations of functions K , that is a K with finitely many
a R nonzero. Aninnerproduct , on thisspace isdxefined by xi∈X i xi
i K
∈ h· ·i P
K ,K := K(x ,x )
h xi xjiK i j
and extendinglinearly. Onetakesthecompletionto obtain .
K
H
Alternatively, an equivalent definition of an RKHS is as a Hilbert space of real-valued functions
onX forwhich theevaluationfunctionalδ (f) := f(x) iscontinuousforallx X.
x
∈
Finite dimensional subspaces of can also be naturally defined by taking a finite number of
K
H
pointsx := x ,...,x X and consideringthelinearspan
1 m
{ } ⊂
:= sp K : x x .
K,x x
H { ∈ }
In practicewewillseek functionsinthesefinitedimensionalsubspaces as approximationsfor f∗.
WithintheseHilbertspaces thereproducing propertyholds:
K ,f = f(x), f . (4.1)
x K K
h i ∀ ∈ H
Ifwedenoteκ := sup K(x,x), then C(X)and itfollowsthat
x∈X HK ⊂
p f κ f , f . (4.2)
L∞(X) K K
|| || ≤ || || ∀ ∈ H
TheRKHS can alsobe defined by meansofan integraloperator. Let ρbeany (finite)strictly
K
H
positive Borel measure on X (e.g. Lebesgue measure) and L2(X) be the Hilbert space of square
ρ
integrablefunctionson X. Thendefine thelinearoperatorL : L2(X) C(X) by
K ρ →
(L f)(x) = K(x,y)f(y)dρ(y). (4.3)
K
ZX
When composed with the inclusion C(X) ֒ L2(X) we obtain an operator from L2(X) to L2(X),
→ ρ ρ ρ
which wealso denoteby L . L is then a self-adjointcompact operator, and it is also positiveifthe
K K
kernel K ispositivedefinite. Alsothemap
L1/2 : L2(X)
K ρ → HK
defines an isomorphismof Hilbert spaces. L1/2 is well defined as an operator on L2(X) in the sense
K ρ
thatL = L1/2 L1/2.
K K ◦ K
4.2 Sobolev space RKHS
In thispaperwewillwork withreproducingkernel Hilbertspacesthat areSobolevspaces. Giventhe
open domain Rd, for k N , 1 p < , the Sobolev space Wk( ) consists of all functions
B ⊂ ∈ 0 ≤ ∞ p B
f with weak derivatives Dαf Lp( ), α k. We also use the following notation to define the
∈ B | | ≤
(semi-)norms
1/p 1/p
f = Dαf p , f = Dαf p .
| |Wpk(B)  || ||Lp(B) || ||Wpk(B)  || ||Lp(B)
|α|=k |α|≤k
X X
   
ApproximationofLyapunovfunctionsfromnoisydata 7
Withp = thenormsare defined in thenatural way:
∞
f = sup Dαf and f = sup Dαf .
Wk(B) L∞(B) Wk(B) L∞(B)
| | ∞ || || || || ∞ || ||
α=k α≤k
Wewillalsousefractional orderSobolevspaces. Foradetaileddiscussionseee.g. [1].
The Sobolev embedding theorem states that for τ > d/2, Wτ(Rd) can be embedded into C(Rd),
2
andthereforeitfollowsthatWτ(Rd)isaRKHS,usingthefactthatthepointwiseevaluationfunctional
2
is then continuous. Several kernel functions are known to generate RKHS spaces that are norm-
equivalent to Sobolev spaces [24, 25]. We will choose to work with the Wendland functions [34].
These are positivedefinite, compactly supported radial basis function kernels that are represented by
aunivariatepolynomialon theirsupport.
Definition 4.1(Wendland function) Let l N, k N . Wedefinebyrecursion
0
∈ ∈
ψ (r) = (1 r)l
l,0 − +
1
andψ (r) = tψ (t)dt
l,k+1 l,k
Zr
forr R+. Here, x = x forx 0 andx = 0for x < 0.
∈ 0 + ≥ +
Settingl := d +k+1,theWendlandfunctionsarecharacterised byasmoothnessindexk N,
⌊2⌋ ∈
and belong to C2k(Rd). For a domain D Rd with a Lipschitz boundary, the Wendland function
⊂
kernel is given by K(x,y) := ψ (c x y ), c > 0, for x,y D. The Wendland function
l,k Rd
|| − || ∈
kernel generates an RKHS consisting of the same functions as the Sobolev space Wτ(D) with τ =
2
k+(d+1)/2,withanequivalentnorm[35,Corollary10.48]. ThereforethegeneratedSobolevspace
isofintegerorderwhen d isodd,and integerplusonehalfwhend iseven.
From now on we shall use RKHS spaces generated by Wendland function kernels. We will use
two such RKHS spaces for the two parts of our algorithm: to approximate the vector field f∗ in
(2.1) we use the space defined on X, corresponding to the Wendland function kernel K1 with
K1
H
smoothnessindexk , suchthatτ = k +(d+1)/2with τ = ν . Then isnorm-equivalentto
1 1 1 1 1 K1
⌈ ⌉ H
Wτ1(X). Inthiscase,whendisoddwehavethat τ = τ ,andtheassumptionthatν (3d+7)/2
2 ⌈ 1⌉ 1 1 ≥
impliesthatk d+3. Whendiseven, τ = τ +1/2,and theassumptionthatν (3d+12)/2
1 1 1 1
≥ ⌈ ⌉ ≥
givesk d+5 (cf. Theorem 2.1).
1
≥
In the second part of our algorithm we approximate the Lyapunov function. For this, we use the
RKHSspace defined on Ω (resp. Ω )correspondingtotheWendland functionkernelK2 with
K2 V T
H
smoothness index k , such that k = k (d + 2) if d is odd, or k = k (d + 3) if d is even.
2 2 1 2 1
− −
Correspondingly, this implies that k 1 when d is odd, and k 2 when d is even. Here, is
2 2 K2
≥ ≥ H
norm-equivalentto Wτ2(Ω ) (resp. Wτ2(Ω )), whereτ := k +(d+1)/2.
2 V 2 T 2 2
These function spaces consistof smoothfunctions as a consequence of thefollowinggeneralised
Sobolevinequality(for aproof,seee.g. [9, Chapter5.7, Theorem6]).
Lemma 4.2 Let Rd beaboundedopensetwithC1 boundary. Foru Wm( )wherem > d/2,
B ⊂ ∈ 2 B
we have
u C u , (4.4)
|| ||Cm−⌊d2⌋−1,γ(B) ≤ || ||W2m(B)
where γ = 1 ifd isodd, andγ isanyelement in(0,1)ifd iseven.
2
ApproximationofLyapunovfunctionsfromnoisydata 8
Corollary4.3 Forf andg , we have
K1 K2
∈ H ∈ H
f C f , (when d isodd), (4.5)
|| ||Ck1(X) ≤ || ||K1
f C f , (when d iseven), (4.6)
|| ||Ck1−1(X) ≤ || ||K1
g C g . (4.7)
|| ||C1(ΩV/ΩT) ≤ || ||K2
PROOF: Following the arguments given in [23], there exists a bounded extension operator E that
extendsf Wτ1(X)to Ef Wτ1(Rd) such thatEu = uon X. Then,from Lemma4.2wehave
∈ 2 ∈ 2
f = Ef Ef
|| ||Cτ1−⌊d2⌋−1,γ(X) || ||Cτ1−⌊d2⌋−1,γ(X) ≤ || ||Cτ1−⌊d2⌋−1,γ(Rd)
˜
≤ C||Ef||W2τ1(Rd) ≤ C||f||W2τ1(X).
Then (4.5) and (4.6) follow from using the norm equivalence of to Wτ1(X), and that τ =
HK1 2 1
k +(d+1)/2. Inequality(4.7)followssimilarly,alsousingk 1whendisodd,andk 2when
1 2 2
≥ ≥
d iseven. (cid:3)
4.3 Generalised interpolant and approximation theorems
In this section we introduce the generalised interpolant that is used to approximate the Lyapunov
functions V and T. Consider a general interpolation setting where Ω Rd is a bounded domain
⊂
having a Lipschitz boundary. Let L be a linear differential operator and q := q ,...,q Ω be
1 M
{ } ⊂
a set of pairwise distinct points which do not contain any singular points of L. (A point q Rd is a
∈
singularpointofL ifδ L = 0, seealso [14].) Wedefinelinearfunctionals
q
◦
λj(u) := δ L(u) = Lu(q ).
qj ◦ j
Definition 4.4(Generalized interpolant) SupposewehavevaluesLu(q ) = β fori = 1,...,M for
i i
a function u : Ω R. Given a sufficiently smooth kernel K( , ), define the generalised interpolant
→ · ·
as
M
s = α (δ L)yK( ,y),
j qj ◦ ·
j=1
X
where(δ L)y isthelinearfunctionappliedtooneargumentofthekernel. Thecoefficient vectorα
qj ◦
isthesolutionofAα = β = (β ) with theinterpolationmatrixA RM×M given by
i
∈
(A) = (δ L)x(δ L)yK(x,y).
ij qi ◦ qj ◦
Remark4.5 According to the above definition, we have Ls(q ) = β . In addition, it can be shown
i i
that the generalised interpolant above is the unique norm-minimal interpolant in , see [11, 14].
K
H
ThematrixAisguaranteedtobeinvertibleduetoourchoiceoftheWendlandfunctionkernelK [11,
Section3.2.2],andsinceqdoes notcontainanysingularpoints.
We concludethissectionby citingthefollowingtheorem from [14], which providesconvergence
estimates for approximating Lyapunov functions V (resp. T ) with the generalised interpolants s
g g 1
(resp. s )as above,foraknowndynamicalsystemx˙ = g(x).
2
ApproximationofLyapunovfunctionsfromnoisydata 9
Theorem 4.6 Let ν := τ with τ = k + (d + 1)/2, where k is the smoothness index of the
2 2 2 2 2
⌈ ⌉
compactly supported Wendland function. Let k > 1/2 if d is odd or k > 1 if d is even. Consider
2 2
thedynamicalsystemdefinedbytheordinarydifferentialequationx˙ = g(x),whereg Cν2(Rd,Rd).
∈
Let x Rd be an equilibrium such that the real parts of all eigenvalues of Dg(x) are negative, and
∈
supposeg isboundedin A(x).
Let Γ be a non-characteristic hypersurface as in Definition 3.4, with h Cν2(Rd,R) and ξ
T
∈ ∈
Cν2(Γ,R). LetV Wτ2(A(x),R),T Wτ2(A(x) x ,R)betheLyapunovfunctionsofTheorems
g ∈ 2 g ∈ 2 \{ }
3.2and 3.3forthesystemx˙ = g(x), withT (x) = ξ (x) forx Γ.
g T
∈
Given pairwisedistinctsetsofpointsq := (q )M inA(x) (notcontainingx)andq˜ := (q )M+N
i i=1 i i=M+1
in Γ, and let Ω A(x) be a boundeddomain with Lipschitz boundary,with q Ω. Let s and s be
1 2
⊂ ⊂
thegeneralisedinterpolantssatisfying
λi(s ) := s (q ),g(q ) = p(q ), i = 1,...,M,
g 1 h∇ 1 i i iRd − i
λi(s ) := s (q ),g(q ) = c, i = 1,...,M,
g 2 h∇ 1 i i iRd −
λi,0(s ) := s (q ) = ξ (q ), i = M +1,...,M +N,
2 2 i T i
and let the fill distance of the set of points q in Ω and q˜ in Γ be h and h respectively. Then for
q q˜
h ,h sufficientlysmall,thefollowingestimateshold:
q q˜
k−1
||h∇s1,giRd −h∇Vg,giRd||L∞(Ω) ≤ Chq 2||Vg||W2τ2(Ω), (4.8)
k−1
||h∇s2,giRd −h∇Tg,giRd||L∞(Ω) ≤ Chq 2||Tg||W2τ2(Ω), (4.9)
k+1
||s2 −Tg||L∞(Γ) ≤ Chq˜ 2||Tg||W2τ2(Ω). (4.10)
5 The Algorithm
Here we present the algorithm for which the estimate given in Theorem 2.1 holds. The algorithm is
actuallysplitintotwoparts. Thefirstpartcomputesf asanapproximationtof∗ (Algorithm1),and
z,λ
thesecondpartcomputesVˆ orTˆ asanapproximationtotheLyapunovfunctionsV orT respectively
given in Theorems 3.2 and 3.3 (Algorithm 2A or 2B). As discussed in 4.2, we will use two RKHS
§
spaces and for the two parts of the algorithm, corresponding to the Wendland function
K1 K2
H H
kernels K1 and K2 with smoothness indices k and k respectively. We recall that the smoothness
1 2
indicesarechosensuchthatτ = k +(d+1)/2with τ = ν ,andk = k (d+3)/2ifdisodd,
1 1 1 1 2 1
⌈ ⌉ −
ork = k (d+4)/2ifdis even.
2 1
−
Recallthatoursampleddatavaluesz := (x ,y )m (X Rd)m taketheformy = f∗(x )+η ,
i i i=1 ∈ × i i xi
with η Rd a random variable drawn from a probability distribution with zero mean and variance
x
∈
σ2.
x
Our approximation scheme for f∗ employs a regularised least squares algorithm (see e.g. [10]
anditsreferences) toapproximateeachcomponentf∗,k,k = 1,...,d. Wealsointroduceaweighting
w = w m correspondingto theVoronoi tessellationassociatedwiththepoints x m [32].
{ xi}i=1 { i}i=1
Definition 5.1(Voronoi tessellation) Let X Rd be compact. For a set of pairwise distinct points
⊂
x := x m Xm, the Voronoi tessellation is the collection of pairwise disjoint open sets (x),
{ i}i=1 ∈ Vi
i = 1,...,mdefined by
(x) = y X x y < x y , if i = j .
i i Rd j Rd
V { ∈ ||| − || || − || 6 }
ApproximationofLyapunovfunctionsfromnoisydata 10
Theweightingw = w m isthendefinedbyw = ρ( (x)),whereρisthestrictlypositiveBorel
{ xi}i=1 xi Vi
measurefrom 4.1.
§
Algorithm 1 Fix a regularisationparameter λ > 0, and define D Rm×m as the diagonal matrix
w
∈
with diagonal elements w , i = 1,...,m. The approximationf for f∗ is constructed component-
xi z,λ
wise. That is,for eachk = 1,...,d, weapproximatethek-th componentf∗,k byfk , defined
z,λ ∈ HK1
by fk = m a K1 , where the coefficients a := a m may be calculated as the solution to the
z,λ i=1 i xi { i}i=1
matrixequation
P
(A D A +λA )a = A D yk.
x w x x x w
Hereyk = (yk)m whereyk issimplythek-thcomponentofy Rd,andA Rm×m isasymmetric
i i=1 i i ∈ x ∈
matrixdefined by(A ) = K1(x ,x ).
x i,j i j
We note here that due to our choice of RKHS the matrix A is positive definite [14, Proposition
x
3.3],and thereforethematrix(A D A +λA ) isinvertible,as itisthesumoftwo positivedefinite
x w x x
matrices. The error in our approximation of f∗ by f is studied in 6. We will show that this error
z,λ
§
may bebounded in thesupremumnormon thedomain X, dependingon thedensityof thedatain X
and thechoiceofregularisationparameter λ.
Once we have our approximation f , then we construct our Lyapunov function approximation
z,λ
withthegeneralisedinterpolantas inTheorem 4.6, whereweset g = f .
z,λ
Therefore we have the following Algorithm 2A for Vˆ, or Algorithm 2B for Tˆ. These algorithms
involvesamplingourapproximationf at adiscreteset ofpoints.
z,λ
Algorithm 2A(ApproximationofV) First run Algorithm 1 on the sampled data set (x ,y )m to
i i i=1
computef ( )d.
z,λ K1
∈ H
Define a set of pairwise distinct points q := (q )M , with f (q ) = 0 for all i = 1,...,M.
i i=1 z,λ i 6
The approximation Vˆ for the Lyapunov function V (see Theorem 3.2) is given by Vˆ =
K2
∈ H
M b λi,y K2( ,y),where λi,y is thelinearfunctionalλi appliedto oneargument of thekernel.
i=1 i fz,λ · fz,λ fz,λ
Thecoefficients b := b M aregiven by
P { i}i=1
B b = p.
q
−
HereB RM×M isa symmetricmatrixdefined by(B ) = λi,x λj,y K2(x,y)andp = (p(q ))M .
q ∈ q i,j fz,λ fz,λ i i=1
As before, the choice of RKHS guarantees that the matrix B will be positive definite, provided
q
thatf (q ) = 0 foralli = 1,...,M.
z,λ i
6
To approximate T, we assume that a non-characteristic hypersurface for f∗ has been defined as
Γ = x A(x) x h(x) = 0 according to Definition 3.4, for which T(x) = ξ (x) on Γ,
T
{ ∈ \{ } | }
ξ Cν1(Γ,R).
T
∈
Algorithm 2B(ApproximationofT) First run Algorithm 1 on the sampled data set (x ,y )m to
i i i=1
computef ( )d.
z,λ K1
∈ H
Define a set of pairwise distinct points q := (q )M+N, with f (q ) = 0 for all i = 1,...,M,
i i=1 z,λ i 6
and q Γ for i = M +1,...,M +N. The approximation Tˆ for the Lyapunov function T
i K2
∈ ∈ H