Efficient a-posteriori error estimation for nonlinear ... - Semantic Scholar
D. Wirtz and B. Haasdonk
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems Stuttgart, December 2010 (rev. June 2011) Institute of Applied Analysis and Numerical Simulation, University of Stuttgart, Pfaffenwaldring 57 D-70569 Stuttgart, Germany {daniel.wirtz,bernard.haasdonk}@mathematik.uni-stuttgart.de www.agh.ians.uni-stuttgart.de
Abstract In this paper we consider the topic of model reduction for nonlinear dynamical systems based on kernel expansions. Our approach allows for a full offline/online decomposition and efficient online computation of the reduced model. In particular we derive an a-posteriori state-space error estimator for the reduction error. A key ingredient is a local Lipschitz constant estimation that enables rigorous a-posteriori error estimation. The computation of the error estimator is realized by solving an auxiliary differential equation during online simulations. Estimation iterations can be performed that allow a balancing between estimation sharpness and computation time. Numerical experiments demonstrate the rigour and effectiveness of the error bounds. Keywords nonlinear dynamical systems · model reduction · kernel methods · a-posteriori error estimates · offline/online decomposition · subspace projection
Preprint Series Issue No. 2010-83, rev. June 2011 Stuttgart Research Centre for Simulation Technology (SRC SimTech) SimTech – Cluster of Excellence Pfaffenwaldring 7a 70569 Stuttgart [email protected] www.simtech.uni-stuttgart.de
2
D. Wirtz and B. Haasdonk
1 Introduction Nowadays, modeling of real world processes like biochemical reactions or electric circuits naturally leads to a formulation as dynamical systems with inputs and outputs. Mathematically, they are described by systems of ordinary differential equations (ODEs) and can be roughly categorized into linear and nonlinear types. Even though computational power has significantly increased over the last years, high-resolution models often result in large-scale dynamical systems that are expensive to simulate. In this context the need for fast simulation is particularly evident in many-query scenarios. They comprise parameter studies or inverse problems where multiple simulations have to be performed for different inputs or initial state configurations. Additionally, some dynamical systems model processes in a real-time setting like control components and thus also require fast computation without strong hardware. Within all those settings fast and rigorous error estimation procedures are important in order to quantify the model errors introduced by the reduction process. Consequently, model reduction techniques for dynamical systems are nowadays subject to intensive research. Work has been done for various types of linear systems, e.g. time invariant systems [1], time-variant [2] or parameterized [3]. For an overview we refer the reader to [4–6]. Reduction techniques for nonlinear dynamical systems have been less investigated, not at least because of the arising difficulties. As with linear systems, most of them involve projection of the system into a lower dimensional subspace. One well known method for nonlinear model reduction is the trajectory piecewise linear (TPWL) approach [7], which is extended by moment matching techniques in [8] or to a piecewise-polynomial scheme in [9]. Various extensions of the balanced truncation procedure [10] to nonlinear systems are investigated in [11–13]. An “approximate reduction” method is introduced in [14] and model reduction for weakly nonlinear systems by bilinearization has been discussed in [15], for example. In addition, the special class of bilinear quadratic nonlinearities has recently been investigated in [16]. A further promising technique for nonlinear model reduction of dynamical systems is by means of discrete empirical interpolation [17]. Finally, computing subspaces by proper orthogonal decomposition (POD) of presumably statistically representative trajectories [18] is an expensive but well-established method. This method is also known as principal component analysis (PCA), see [19] and references therein for an overview. In this work we adopt principles of the reduction technique proposed in [20], extend the class of usable kernels and provide novel, efficient a-posteriori error estimators for the resulting reduced systems. The structure of our paper is as follows: In Section 2 we introduce our base dynamical system and discuss aspects of the reduction process. Section 3 is concerned with a-posteriori error estimates and introduces our new local Lipschitz constant estimation method. Next, Section 4 presents numerical experiments for synthetic dynamical systems and we conclude with Section 5. Some auxiliary mathematical details are presented in the appendix. 2 Reduction of kernel based systems The starting point for our investigations is the reduction method introduced in [20], which combines subspace projection with kernel methods. In model reduction, the latter are used to approximate the nonlinearities of dynamical systems in kernel spaces and promise a large potential in the field. Kernel methods comprise applications from machine learning like support vector regression [21] as well as kernel interpolation with corresponding theoretical foundation [22–24]. However, those approximation techniques are out of the scope of this article. Instead we focus on nonlinear systems whose inner dynamics are already given by a kernel expansion and investigate the projection behaviour and the resulting a-posteriori error estimators. 2.1 The base dynamical system Our central assumption is to have a nonlinear kernel expansion f (x) =
N X i=1
ci Φ(x, xi ),
(1)
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
3
where X := {x1 , . . . , xn } ⊂ Rd are the centers or support vectors of the expansion and ci ∈ Rd the coefficient vectors. The function Φ is a kernel, which is basically a symmetric function Φ : Rd ×Rd −→ R. We omit further characterizations of kernels and mention here that the key analytical property is that for some open Ω ⊂ Rd each symmetric positive definite kernel spans a unique Hilbert space NΦ (Ω) of functions. Those Hilbert spaces have a special reproducing property and are commonly referred to as reproducing kernel Hilbert spaces (RKHS). In our context those RKHS serve as a base class for the d dynamical systems functions f , so (1) implies f ∈ (NΦ (Ω)) . For more details on kernels, RKHS and their applications we refer to [22, 25, 26], for example. Now, the base system considered throughout this paper is of the form x0 (t) = f (x) + Bu(t), x(0) = x0 , y(t) = Cx(t),
(2) (3)
where x(t) ∈ Rd denotes the state of the system at times t ∈ [0, T ], 0 ≤ T < ∞, x0 ∈ Rd the initial condition, u : [0, T ] → Rm an input function with B ∈ Rd×m and output y(t) ∈ Rk with C ∈ Rk×d . Additionally, in some applications different norms than the standard Euclidean norm are taken for the state space. Thus, let G ∈ Rd×d be a symmetric positivepdefinite matrix. Then G defines a scalar product hx, yiG := xt Gy on Rd with induced norm ||x||G := hx, xiG . Choosing G := Id will yield the standard Euclidean inner product h·, ·i and 2-norm ||·||, where Id denotes the d-dimensional identity matrix. For example, G is usually chosen as the Gram-matrix of the finite element/ finite volume basis for dynamical systems obtained from partial differential equation discretizations. Having introduced the underlying dynamical system (2)-(3) we now detail the reduction approach using subspace projection. We assume to have two biorthogonal projection matrices V, W ∈ Rd×r , W t V = Ir with r d denoting the reduced system’s dimension. The reduced system is then obtained by applying a Petrov-Galerkin projection to the full system: z 0 (t) = W t f (V z(t)) + W t Bu(t), z(0) = W t x0 , y r (t) = CV z(t),
(4) (5)
where z(t) ∈ Rr now denotes the reduced system’s state and y r (t) the approximate output. Let U := hv1 , . . . , vr i ⊂ Rd be the space spanned by the columns vi ∈ Rd of V . Consequently, the reconstructed approximate solution and output are given by xr (t) := V z(t) ∈ U and y r (t) = Cxr (t). The requirements on V, W are now to yield a good approximation xr (t) ≈ x(t) and thus y r (t) ≈ y(t) for different initial values and inputs. For the remainder of this work we will assume the matrix computation method for V, W as a black-box since our method is applicable using any basis/projection matrix generation method. The rest of this section is concerned with the different projection aspects of the reduction technique for kernel-based dynamical systems. We will see that the structure of the kernel expansion allows for dramatic reduction of the computational costs using certain kernels.
2.2 System projection At first, the projection into Rr by W t can be applied directly to the coefficient vectors ci via cri := PN W t ci , i = 1 . . . N . This results in a reduced function f r (z) := i=1 cri Φ(V z, xi ), with new coefficient r r vectors ci ∈ R . Unfortunately, the reduced system (4) is still expensive to simulate given that f is evaluated at xr (t) ∈ Rd . In order to avoid input arguments of high dimension d, we consider two special classes of kernels. 2.2.1 Inner product kernels The first type of kernels that allow efficient argument evaluations are the inner product kernels also mentioned in [20]. It is assumed that Φ(x, y) = φ(hx, yiG ) for some scalar function φ : R → R. Then Φ(V z, xi ) = φ(hV z, xi iG ) = φ(hz, zi i) =: Φr (z, zi ) for zi := V t Gxi ∈ Rr , i = 1 . . . N . This way, it is sufficient to project the centers xi into the low-dimensional subspace and the evaluations of Φ can
4
D. Wirtz and B. Haasdonk
be computed efficiently and loss-less during the reduced simulation via Φr . Some examples for those kernels are the linear kernel Φ(x, y) = hx, yiG or polynomial kernels Φ(x, y) = (1 + hx, yiG )p of degree p ∈ N. 2.2.2 Translation- and rotation invariant kernels In extension to [20] we present here a further class of kernels that allow for efficient loss-less argument evaluations: The translation and rotation-invariant kernels Φ(x, y) := φ(||x − y||G ) for some scalar function φ : R+ 0 → R. We impose the additional requirement xi ∈ U, i.e. xi = V zi for some zi ∈ Rr , i = 1 . . . N . Then we obtain Φ(V z, xi ) = φ(||V z − V zi ||G ) = φ(||z − zi ||V t GV ) =: Φr (z, zi ) with V t GV being a small Rr×r matrix inducing a new norm on Rr . Note that the assumption xi ∈ U is of a technical nature. We either extend U by the span of the xi , or, if the kernel expansion is created with knowledge of U , one can choose xi ∈ U in the first place. Those kernels are also commonly referred to 2 as radial basis functions, of which the Gaussian kernel Φ(x, y) = exp(− ||x − y|| /β 2 ) is probably the most popular example. For further examples and characterizations of the above mentioned kernels we refer to [21, 26].
2.3 The reduced system In conclusion, for any of the specific kernels above we obtain from (4)-(5) an equivalent reduced system: z 0 (t) =
N X
cri Φr (z, zi ) + B r u(t),
(6)
i=1
z(0) = z0 ,
y r (t) = C r z(t),
(7)
with z0 := W t x0 , B r := W t B, C r := CV . The evaluation complexity is independent of d and the system has only low dimensional components. Note again, that this reduction technique (using inner-product kernels) has been investigated in [20]. Despite our new class of kernels, investigation of the reduction properties is out of scope of this article. Instead, our main focus is the error estimation.
3 A-posteriori error estimates When dealing with reduced dynamical systems it is essential to have means to compute or at least estimate the error that is made during simulations. Still, little work has been done so far in the field of error estimations for nonlinear dynamical systems. A-priori error analysis of POD reduction schemes in different applications are considered e.g. in [27, 28] and a-priori error bounds for certain nonlinear systems are derived in [29]. Error estimates using dual-weighted-residual methods can be found in [30] and [7] perform a-posteriori error estimation for their TPWL approach. Recent work on a-posteriori error estimation has also been done for nonlinear parametrized evolution equations [31] and a related scheme for error estimators of linear but parameterized dynamical systems can be found in [3, 32]. Note that in many cases differential inequalities such as the comparison lemma or Gronwall’s lemma [33, p.36] are used to bound ODE solutions/errors a-priori. Instead, we will use it to provide a-posteriori error bounds which can be computed along with the reduced system for any input u(t). Our estimators bound the error in the state-space for finite time intervals instead of giving frequency domain error bounds. For the reduced system (6)-(7) we derive rigorous and efficient error estimators. In Section 3.1 we introduce the estimator structure obtained using the comparison lemma [33, p.32] (also applied in [7]). Section 3.2 introduces the local Lipschitz constant estimations utilized in Section 3.3 to obtain our main estimator result. We complete our results with convergence analysis and computational aspects in Sections 3.4 and 3.5.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
5
3.1 Global Lipschitz constant a-posteriori error estimator The state space error e(t) := x(t) − xr (t) for any t ∈ [0, T ] is given by the solution of the error system e0 (t) = f (x(t)) + Bu(t) − V W t f (xr (t)) − V W t Bu(t), e(0) = x0 − V W t x0 .
(8) (9)
Estimating the norm of the derivative and application of the comparison lemma yields a first aposteriori error estimator. Let Ex0 := ||(Id − V W t )x0 ||G be the initial error for the remainder of this work. Theorem 1 (Global Lipschitz constant a-posteriori error Estimator (GLE)) Let Φ be a kernel. Further, let Φ be LΦ -Lipschitz continuous with respect to the first variable. Then the state space error is bounded via ||e(t)||G ≤ ∆GLE (t) ∀ t ∈ [0, T ] with Zt ∆GLE (t) :=
α(s)eβ(t−s) ds + Ex0 ,
0 α(t) := Id − V W t f (xr (t)) + Bu(t) G , β := LΦ
N X
(10)
||ci ||G .
i=1
Proof (Proof ) With (1) we directly obtain N X ||f (x(t)) − f (xr (t))||G = ci Φ(x(t), xi ) − Φ(xr (t), xi ) i=1
≤
N X
(11) G
||ci ||G Φ(x(t), xi ) − Φ(xr (t), xi )
i=1
≤ LΦ
N X
||ci ||G ||e(t)||G = β ||e(t)||G .
i=1
Now from (8) we get
||e0 (t)||G = f (x(t)) − f (xr (t)) + Id − V W t f (xr (t)) + Bu(t) G ≤ β ||e(t)||G + α(t).
Finally, application of the comparison lemma [33, p.32] for the ODE v 0 (t) = βv(t) + α(t), v(0) = Ex0 yields the error estimator as analytical solution. Note that the error estimator ∆GLE (t) is rigorous, continuous and monotonically increasing. Once a bound for the state error is available the output error ey (t) := y(t) − y r (t) can be bounded by ||ey (t)|| ≤ ||C||G ||e(t)||G ≤ ||C||G ∆GLE (t), t ∈ [0, T ], where ||C||G is the G-induced matrix norm. 3.2 Local Lipschitz constants Closer investigation of the error estimator from Theorem 1 shows that if the kernel Lipschitz constant LΦ is large, the estimator will grow fast and the estimation will be overly conservative. Therefore, we propose a local Lipschitz constant computation that utilizes the current position xr (t) of the reduced state variable and an a-priori error bound ||e(t)||G ≤ C(t) to improve the estimation procedure. The benefit of knowing xr (t) is that one argument of the estimation ||f (x(t)) − f (xr (t))||G ≤ β ||x(t) − xr (t)||G is fixed, which will be used to compute a local Lipschitz constant using secant gradients. Also, the a-priori bound C(t) will enable efficient computation of those secant gradients utilizing the kernel expansion center positions.
6
D. Wirtz and B. Haasdonk
The key to our local Lipschitz constant computations is to use the class of kernels introduced in Section 2.2.2. In particular, we use radial basis functions Φ(x, y) := φ(||x − y||), x, y ∈ Rd induced by a certain class of scalar functions φ: + Definition 1 (Bell functions) A function φ : R+ 0 → R0 is called a bell function, if
i) ii) iii) iv)
φ ∈ C 2 (R+ 0 ), ||φ||L∞ ≤ B, B > 0, φ0 (0) ≤ 0, ∃ r0 > 0 : φ00 (r)(r − r0 ) > 0 ∀ r 6= r0 .
(12) (13) (14)
Here, r0 denotes the unique turning point where φ is strictly concave/convex on [0, r0 ]/[r0 , ∞[, re2 spectively. For kernel is induced by a bell function φ(r) = e−(r/β) . Then we √ example, the Gaussian have r0 = β/ 2 and LΦ = |φ0 (r0 )| in the context of Theorem 1. Let φ denote a bell function for the remainder of this work. We will state results regarding bell functions first and connect these to the multidimensional case later. + For any s ∈ R+ 0 , we introduce a secant gradient function Ss : R0 −→ R with ( φ(r)−φ(s) , r 6= s r−s Ss (r) := (15) 0 φ (s), r = s. Well-definedness and continuous differentiability of Ss can easily be verified by elementary analysis. The following lemma characterizes the minimum secant gradient position. Lemma 1 (Minimum secant gradient) Let s ∈ R+ 0 . Then rs := arg min+ Ss (r) r∈R0
(16)
exists, is unique and exactly one of the following cases holds: s < r0 ≤ r s , rs ≤ r0 < s, s = r0 = rs .
(17) (18) (19)
For the proof we refer to SectionA.1 of the appendix. Figure 1 shows two examples for minimum secant gradients and the position of rs with each s < r0 (left) and s > r0 (right), using the Gaussian inducing bell function φ(r) = exp(−r2 /β 2 ) with β = 2. Now, further assume to restrict admissible r values to r ∈ BC (s) ∩ R+ 0 for a C > 0, where BC (s) denotes an open ball of radius C around s and BC (s) the closure. With this condition we can state our novel result regarding local Lipschitz constant computations. Proposition 1 (Local Lipschitz estimation using secant gradients) Let C > 0, s ∈ R+ 0 and Ω := BC (s) ∩ R+ . Then |φ(r) − φ(s)| ≤ |S (r )||r − s| ∀ r ∈ Ω with s C,s 0 ( s + sign (rs − s) C rs ∈ /Ω rC,s = (20) rs rs ∈ Ω. Proof (Proof ) Assume s 6= r0 . Then Lemma2 in SectionA.2 of the appendix implies that Ss is monotonically decreasing on [0, rs ] and increasing on [rs , ∞[. For rs ∈ / Ω, as Ss is negative, |Ss | takes its maximum value at the border of Ω that is closest to rs . Since arg minr∈R+ Ss (r) = arg maxr∈R+ |Ss (r)|, 0 0 case rs ∈ Ω follows by definition. The case s = r0 implies rs = r0 and hence case two, as of course r0 ∈ Ω = BC (r0 ) and C > 0. We finally conclude |φ(x) − φ(y)| ≤ ||Ss ||L∞ (Ω) |r − s| = |Ss (rC,s )||r − s|. Next, we generalize these results to the multidimensional kernel setting.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems Bell function demo, s=0.282843, rs=2.058702 1.2
Bell function demo, s=5.391323, rs=0.400439 1.2
φ φ(rs) + φ’(rs)(rs−s)
φ(s)
1
7
0.8
φ φ(rs) + φ’(rs)(rs−s)
φ(r ) s
1
0.8
φ(r0)
0.6
φ(rs)
0.4
0.4
0.2
0.2
φ(rm)
0
r0
s
φ(r0)
0.6
rs
φ(rm)
0
rm
−0.2
rs
r0
φ(s) s
rm
−0.2
−0.4
−0.4 0
1
2
3
4
5
Fig. 1 Gaussian φ with β = 2, r0 =
√
6
2 and rm =
7
√
0
2
1−e−1/2
1
2
3
4
5
6
7
≈ 3.5942
Corollary 1 (Local Lipschitz constants for bell function kernels) Let x, y, z ∈ Rd , C > 0, Φ(x, z) = φ(||x − z||G ) and s = ||y − z||G . Then we have |Φ(x, z)−Φ(y, z)| ≤ |Ss (rC,s )| ||x − y||G ∀ x ∈ BC (y). Proof (Proof ) Fix y, z ∈ Rd , C > 0, let R := x ∈ Rd ||x − z||G ∈ BC (s) . Then using Proposition 1 we estimate |Φ(x, z) − Φ(y, z)| = |φ(||x − z||G ) − φ(||y − z||G )| ≤ |Ss (rC,s )| ||x − z|| − ||y − z|| G
G
≤ |Ss (rC,s )| ||x − y||G ∀ x ∈ R.
The fact that BC (y) ⊂ R finishes the proof.
3.3 Improved error estimator Finally, the local Lipschitz constant estimations from Section 3.2 can be used to state our main result, which is an improved version of the error estimator derived in Theorem 1. As shortly mentioned in the beginning, we additionally assume to have an a-priori coarse error bound ||e(t)||G ≤ C(t) ∀ t ∈ [0, T ],
(21)
which can be C(t) = ∞ if no further knowledge is available. We introduce the notation di (t) := ||xr (t) − xi || ,
i = 1 . . . N,
(22)
for the distance of the reduced state xr (t) to each expansion center xi during the reduced simulation. In fact, the coarse error bound means nothing but x(t) ∈ BC(t) (xr (t)) ∀ t ∈ [0, T ]. Consequently, for any t ∈ [0, T ], i ∈ {1 . . . N } we identify x, y, z ∈ Rd and C > 0 from Corollary 1 with x(t), xr (t) of the full/reduced solution, the kernel expansion centers xi and C(t), respectively. This allows to obtain local Lipschitz constant estimations at each time t using the current reduced simulation’s state xr (t) as “center” of locality: Theorem 2 (Local Secant gradient Lipschitz a-posteriori error Estimator (LSLE)) Let the error system be given as in (8)-(9), where the kernel Φ of the expansion (1) is induced by a bell
8
D. Wirtz and B. Haasdonk
function φ. Further, let C(t) be an a-priori error bound. Then the state space error is bounded via ||e(t)||G ≤ ∆C LSLE (t) ∀ t ∈ [0, T ] with Zt ∆C LSLE (t)
:=
Rt
β(r)dr
α(s)es
ds + Ex0 ,
0
α(t) := Id − V W t f (xr (t)) + Bu(t) G , β(t) :=
N X
||ci ||G Sdi (t) (rC(t),di (t) ) .
(23) (24)
i=1
Proof (Proof ) Similar to (11), with (1) and Corollary 1 we deduce that ||f (x(t)) − f (xr (t))||G ≤ β(t) ||e(t)||G . Performing an analogous estimation as in Theorem 1 we obtain the comparison lemma ODE v 0 (t) = β(t)v(t) + α(t), v(0) = Ex0 , whose analytical solution yields the error estimator. We see in the α(t) term that if any input u(t) happens to maintain the full system in the subspace U , then the resulting zero approximation error will be verified a-posteriori by both the GLE and LSLE error estimators. Before we present numerical experiments we complete the estimator results with some theory regarding the coarse bound and computational aspects.
3.4 Utilization of a-priori bound C(t) We still need to justify the assumption (21) because it requires existence of an a-priori bound in order to obtain an a-posteriori error estimator. The core idea is to use the a-priori bound in an iterative scheme starting with C(t) ≡ ∞, which is trivially satisfied. Then the resulting a-posteriori error estimate is reused as a-priori bound for another iteration of the error estimation procedure. A nice feature of this procedure is that it allows to balance the error bound’s sharpness against the computational costs. The next theorem shows that actual progress can be made under reasonable assumptions and the iterative scheme reproduces the first estimate in the worst case. Theorem 3 (Convergence of error estimator iterations) Let the conditions of Theorem 2 hold Cn−1 and define C0 (t) ≡ ∞ and Cn (t) := ∆LSLE (t), n ≥ 1. Recall the distances di (t) from (22), rs from (16) and for n ∈ N let Γn := t ∈ [0, T ] ∃ i ∈ {1 . . . N } : rdi (t) ∈ / BCn (t) (di (t)) . (25) If Γ1 6= ∅, then Γn+1 = 6 ∅ ∀ n ∈ N, ( n < ∆C Cn+1 LSLE (t), ∆LSLE (t) n = ∆C LSLE (t),
t > inf Γn+1 t ≤ inf Γn+1
∀ n ∈ N.
(26)
n Else if Γ1 = ∅, then ∆C LSLE (t) = C1 (t) ∀ n ≥ 1, t ∈ [0, T ]. Further, the sequence of iterates converges Cn uniformly to a continuous lower error bound ∆∞ LSLE (t) := lim ∆LSLE (t), t ∈ [0, T ].
n→∞
For the proof we refer to Section A.3 of the appendix. The crucial assumption here is that the first estimation C1 is small enough, so that for some time t and expansion center xi the minimum secant gradient point rdi (t) is not contained in BC1 (t) (di (t)). This is reasonable in the context of kernel-based model reduction, since the expansion centers are naturally scattered over the high dimensional state space.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
9
3.5 Computation of error estimators The error estimators from Theorem 1 and 2 allow an efficient computation during the reduced simulation with a complexity independent from the full system’s (high) dimension. As already mentioned in the proofs of the respective theorems, the error estimators are given by the solution of the onedimensional ODE u0 (t) = β(t)v(t)+α(t) with initial condition u(0) = Ex0 (For the GLE from Theorem 1 we have a constant β). This solution can be computed “on the fly” by adding an extra dimension to the reduced system. The efficient computation is possible due to the full offline/online separability of the error estimators and an effective scheme for solving the minimum secant gradient problem, which we will present next. 3.5.1 Offline/online decomposition A vital ingredient to the efficient computation of the error estimators is to decompose the α(t) ODE parts (10),(23) in an offline/online fashion similar to [3]. In this setting, the offline phase comprises expensive computations with complexity depending on the high dimension d and has to be performed once. The online phase is performed for different inputs u(t) and its complexity is dependent only on the reduced dimension r. Here, we confine ourselves to detail the case for the LSLE presented in Theorem 2, as the GLE from Theorem 1 is decomposed the same way but with constant β(t) term. Proposition 2 (Offline/online decomposition of the LSLE) Let the conditions from Theorem 2 hold. Define M := (α1 . . . αN ) ∈ Rd×N and recall the norm inducing matrix G from Section 2.1. Then the offline computations are Ex0 = ||x0 − V W t x0 ||G , ||ci ||G , i = 1 . . . N for the β(t) term and ˜ := Id − V W t M, ˜ := Id − V W t B, M B ˜ t GM ˜ ∈ RN ×N , ˜ t GB ˜ ∈ RN ×m , ˜ t GB ˜ ∈ Rm×m M1 := M M2 := M M3 := B regarding α(t). With v(0) = Ex0 , the online part consists of computing q t t v 0 (t) = β(t)v(t) + ϕ(t) M1 ϕ(t) + 2ϕ(t) M2 u(t) + u(t)t M3 u(t), N ϕ(t) := Φr (z(t), zi ) i=1 ∈ RN .
(27)
The proof follows directly from the definitions. Notice that if there are no inputs u(t) for a given system only the matrix M1 has to be computed as M2 and M3 do not appear. However, all the offline matrices M1 , M2 , M3 are small matrices independent of the original systems dimension. 3.5.2 Efficient computation of secant gradient Lipschitz constants When using the improved error estimation from Theorem 2, another important aspect for fast online evaluation is efficient computation of the β(t) term. In the setting of a reduced simulation this term reads as N X β(t) = ||ci ||G Sdi (t) (rC(t),di (t) ) , i=1
rC(t),di (t)
( di (t) + sign di (t) − rdi (t) C(t), rdi (t) ∈ / BC(t) (di (t)) = , rdi (t) , rdi (t) ∈ BC(t) (di (t))
rdi (t) = arg min+ r∈R0
φ(r) − φ(di (t)) . r − di (t)
So, in theory, N small minimization problems need to be solved at any time t ∈ [0, T ], whereas in practice Lemma 1 allows to reduce the amount of problems to be solved. For example, say di (t) < r0 for a t ∈ [0, T ], i ∈ {1 . . . N }. Then condition (17) yields r0 ≤ rdi (t) . So, if di (t) + C(t) < r0 there is no need to compute rdi (t) as rC(t),di (t) = di (t) + C(t). However, we propose using a penalized Newton iteration scheme for efficient computation of the minimizer itself, which is detailed in Proposition 3 in the appendix. As the reduced state variable xr (t) changes continuously, so does rs , which keeps the number of Newton iterations for subsequent (next time-step) computations of rs small when using the old values as starting point for the new iterations.
10
D. Wirtz and B. Haasdonk
4 Numerical experiments In this part we present numerical experiments for synthetic dynamical systems using kernel expansions as nonlinearity. Before we state the system settings, we introduce two more error estimator modifications. First, the convergence result from Theorem 3 motivates a heuristic variant of our local estimators. This originates from the fact that the limit function ∆∞ LSLE reproduces itself when used as a-priori bound within the iterations. Numerically (up to the integration error) this can be exploited using the error estimate from the previous time step as bound at the current time step. We will refer to this time-discrete method by “LSLE TD”. Experiments indicate that this variant indeed seems to bound ||f (x(t))−f (xr (t))||G the iterated LSLE estimators from below. Second, when setting β(t) = we obtain ||x(t)−xr (t)||G the smallest possible estimation, since this is the best local Lipschitz constant at any t ∈ [0, T ]. As this version requires x(t) it is not considered a practical error estimator but rather a theoretical “Lower Bound”. Our test setting is given as follows: Let d = 240000 to represent a large-scale system, T = 20, G = Id . The kernel expansion (1) uses N = 20 and (xi )j := N50 −1 (i − 1), i = 1 . . . N, j = 1 . . . d. 2 2 The kernel used is a Gaussian Φ(x, y) = exp(− ||x − y|| /β ) with β = 224, which is chosen to have Φ(xi , xj ) < 10−5 ∀ |i − j| ≥ 2. This way a certain locality of the kernel expansion is ensured. Finally, the expansion coefficient vectors are given via (ci )j = exp(−(xi )j /15), i = 1 . . . N, j = 1 . . . d. Further we set x0 = (1 . . . 1)t ∈ Rd and B = (1 . . . 1)t ∈ Rd×1 . To emphasize the input-dependent behavior of 2 1 our estimators, we choose two system inputs u1 (t) = 25 sin( 3t ) and u2 (t) = 21 e−(12−t) which represent one oscillating and one localized input. We use an explicit Euler scheme with time-step ∆t = 0.05 as solver. In order to investigate the error estimator behavior, we decided to have means of controlling the projection subspace quality. Since the system setup is homogeneous in each √ component, we would ˆ = (1, . . . , 1)t / d as projection matrices. obtain zero approximation error for any input using Vˆ = W ˆ Then, for a given degree θ we set the rotated subspace projection matrices to V := RVˆ , W := RW ˆ for an orthogonal block matrix R := diag(R, Id−200 ) with cos(θ) − sin(θ) 200×200 ˆ . R := diag(R2 , . . . , R2 ) ∈ R , R2 = sin(θ) cos(θ) This way, θ continuously controls the quality of the projection subspace. The results for the different error estimators can be seen in Figure 2. The rigorosity of all estimators is verified as all are upper
True error GLE LSLE LSLE, 1 It LSLE, 2 It LSLE, 5 It LSLE TD Lower bound
−3
−4
10
0
2
4
6
8
10 Time
12
14
16
18
True error GLE LSLE LSLE, 1 It LSLE, 2 It LSLE, 5 It LSLE TD Lower bound
−3
10
Relative error estimates
Error estimates
10
20
−4
10
0
2
4
6
8
10 Time
12
14
16
18
20
Fig. 2 Absolute (left) and relative (right) error estimates using θ = .0005 and input u1
bounds of the true error. Iterations of the LSLE estimator improve the error bound. The fact that the LSLE TD variant behaves like a tight lower bound for their corresponding iterative schemes (LSLE, 5 It and LSLE TD are indistinguishable) supports the applicability of our heuristic LSLE TD approach.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
11
Furthermore, the LSLE TD variant grows with the same rate as the lower bound estimation up to about t = 16, which shows the effectiveness of the local Lipschitz constants together with the a-priori bounds. Compared to the GLE, the LSLEs have a significantly lower exponential increase rate. One can see that the first and second LSLE iterations fall back to their specific increase rate after a certain time at which the a-priori bound is too big to have a positive effect regarding the choice (20). Figure 3 displays the computation times versus the error estimations (left) and the output of the full and reduced simulations along with the output error bounds computed by the LSLE TD estimator (right) for θ = 0.05 and u2 . Note that identical conclusions can be drawn by using either input for Figures 2 and 3. The full error computation takes about 20s, and the GLE estimator is cheap but too conservative to have a significant meaning. Comparing identical symbols reveals that indeed iterating the error estimation process increases the computational time, but improves the accuracy of the estimators. Hence, the iteration number can be used as a balancing parameter between online run-time and error bound accuracy. The output of the system using input u2 is nicely bounded by the LSLE TD estimator.
True error GLE LSLE LSLE, 1 It LSLE, 2 It LSLE, 5 It LSLE TD
20
18
16
1.6
1.4
14
12
y(t) / yr(t)
Comp. time [s]
Full system Reduced system Lower bound Upper bound
1.5
10
1.3
1.2 8
1.1 6
1 4
0.9
2 0
10
10
10
20
10
30
10
40
∆(T)
10
50
10
60
10
70
10
0
2
4
6
8
10 t
12
14
16
18
20
Fig. 3 Computation times and output for u2 and θ = .05. The output error is bounded using the LSLE TD variant
Note that the reduced and full models’ outputs are indistinguishable in this plot. Table 1 shows the estimated errors at T = 20 for different θ values and both inputs. Albeit the true ||e(T )||G
GLE
LSLE
LSLE 2 It
LSLE 5 It
LSLE TD
Low. b.
θ = 0.0005, u1
1.3E−5
1.1E+73
3.2E+11
1.4E+2
6.6E−5
6.6E−5
2.2E−5
θ = 0.005, u1
1.3E−4
1.1E+74
3.2E+12
1.2E+5
6.7E−4
6.6E−4
2.2E−4
θ = 0.05, u1
1.3E−3
1.1E+75
3.2E+13
1.2E+8
2.3E+0
6.7E−3
2.2E−3
θ = 0.5, u1
1.3E−2
1.1E+76
3.1E+14
1.0E+11
9.2E+5
8.0E−2
2.2E−2
θ = 0.0005, u2
2.3E−5
1.1E+73
8.5E+11
2.3E+3
1.8E−4
1.8E−4
3.0E−5
θ = 0.005, u2
2.3E−4
1.1E+74
8.5E+12
1.0E+6
2.7E−3
1.8E−3
3.0E−4
θ = 0.05, u2
2.3E−3
1.1E+75
8.5E+13
5.6E+8
3.6E+2
1.9E−2
3.0E−3
θ = 0.5, u2
2.3E−2
1.1E+76
8.4E+14
3.4E+11
3.5E+7
3.6E+0
3.1E−2
Model
Table 1 Errors of estimation runs for different θ values and inputs u1 , u2 at T = 20
error is of the same magnitude for both inputs and any θ, the estimators perform differently for u1 and u2 . This table emphasizes the superiority of the LSLE compared to the GLE, as the non-iterated LSLE variant is about 62 orders of magnitude better than the GLE.
5 Conclusion and Outlook In this article we investigated the topic of model reduction for kernel-based nonlinear dynamical systems and presented a series of improving efficient a-posteriori error estimators. Special classes of kernels allow for the creation of a small reduced system (6)-(7) via loss-less low-dimensional argument evaluations.
12
D. Wirtz and B. Haasdonk
The error estimators are obtained simultaneously with the reduced simulation by solving an additional ODE. Due to a complete decomposition of the error estimators in an offline/online fashion, ultimately a rapid simulation including rigorous error bounds is ensured. Numerical experiments show the superiority of the local Lipschitz estimator using minimum secant gradients obtained via solving a small Newton problem. Furthermore, iterating the error estimation process allows to balance the computation costs against estimation sharpness. Future work will target at the improvement of the error estimators regarding the coefficient vector norms ||ci ||G in (24) as locality can also be used there to sharpen the error estimations. Finally, incorporation of techniques to obtain kernel expansion approximations (1) of arbitrary nonlinear dynamical systems and extension to parameterized systems will be in the focus of our research.
Acknowledgements The authors would like to thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart and the Baden-Wrttemberg Stiftung gGmbH.
References 1. P. Benner, Numerical linear algebra for model reduction in control and simulation, GAMM-Mitt. 29 (2) (2006) 275–296. 2. H. Sandberg, A. Rantzer, Balanced truncation of linear time-varying systems, IEEE T. Automat. Contr. 49 (2) (2004) 217 – 229. 3. B. Haasdonk, M. Ohlberger, Efficient reduced models and a-posteriori error estimation for parametrized dynamical systems by offline/online decomposition, Math. Comp. Model. Dyn. 17 (2) (2011) 145–161. 4. A. Antoulas, An overview of approximation methods for large-scale dynamical systems, Annu. Rev. Contr. 29 (2005) 181–190. 5. A. Antoulas, Approximation of Large–Scale Dynamical Systems, SIAM Publications, 2005. 6. K. Zhou, J. C. Doyle, K. Glover, Robust and Optimal Control, Prentice Hall PTR, 1996. 7. M. Rewienski, A trajectory piecewise-linear approach to model order reduction of nonlinear dynamical systems, Ph.D. thesis, MIT, Cambridge, USA (2003). 8. B. Bond, L. Daniel, Parameterized model order reduction of nonlinear dynamical systems, in: Proc. of ICCAD 2005, 2005, pp. 487–494. 9. N. Dong, J. Roychowdhury, Piecewise polynomial nonlinear model reduction, in: Proc. of DAC 2003, 2003, pp. 484–489. 10. B. Moore, Principal component analysis in linear systems: Controllability, observability, and model reduction, IEEE T. Automat. Contr. 26 (1) (1981) 17–32. 11. J. Scherpen, Balancing for nonlinear systems, Syst. Contr. Lett. 21 (2) (1993) 143–153. 12. S. Lall, J. Marsden, S. Glavaski, A subspace approach to balanced truncation for model reduction of nonlinear control systems, Int. J. Robust Nonlin. 12 (5) (2002) 519–535. 13. M. Condon, R. Ivanov, Empirical balanced truncation of nonlinear systems, J. Nonlinear Sci. 14 (2004) 405–414. 14. P. Tabuada, A. D. Ames, A. Julius, G. J. Pappas, Approximate reduction of dynamic systems, Syst. Contr. Lett. 57 (7) (2008) 538–545. 15. J. Phillips, Projection frameworks for model reduction of weakly nonlinear systems., in: Proc. of DAC 2000, 2000, pp. 184–189. 16. P. Benner, T. Breiten, Krylov-subspace based model reduction of nonlinear circuit models using bilinear and quadratic-linear approximations, Max-Planck-Institute for Dynamics of Complex Technical Systems, Magdeburg, Germany. Submitted (2010). 17. S. Chaturantabut, D. Sorensen, Discrete empirical interpolation for nonlinear model reduction, in: Proc. of CDC/CCC 2009, 2009, pp. 4316–4321. 18. K. Kunisch, V. Volkwein, Proper Orthogonal Decomposition for Optimality Systems, ESAIM, Math. Model. Numer. Anal. 42 (2008) 1–23.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
13
19. I. T. Jolliffe, Principal Component Analysis, Springer-Verlag, 2002. 20. J. Phillips, J. Afonso, A. Oliveira, L. Silveira, Analog macromodeling using kernel methods, in: Proc. of ICCAD 2003, 2003, pp. 446–453. 21. B. Schlkopf, A. J. Smola, Learning with Kernels, The MIT Press, 2002. 22. H. Wendland, Scattered data approximation, Cambridge University Press, 2005. 23. L. Ma, Z. Wu, Kernel based approximation in sobolev spaces with radial basis functions, Appl. Math. Comput. 215 (6) (2009) 2229–2237. 24. H. Wendland, Computational aspects of radial basis function approximation, in: Topics in Multivariate Approximation and Interpolation, Elsevier, 2006, pp. 231–256. 25. N. Aronszajn, Theory of reproducing kernels, Trans. Am. Math. Soc. 68 (3) (1950) 337–404. 26. C. Berg, J. Christensen, P. Ressel, Harmonic Analysis on Semigroups, Springer, 1984. 27. M. Hinze, S. Volkwein, Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: Error estimates and suboptimal control in dimension reduction of large-scale systems, in: Dimension Reduction of Large-Scale Systems, Lecture Notes in Computational and Applied Mathematics, Vol. 45, Springer Berlin Heidelberg, 2005, pp. 261–306. 28. S. Volkwein, Model Reduction using Proper Orthogonal Decomposition, Lecture notes, University of Graz (2008). 29. T. Reis, M. Heinkenschloss, Model reduction with a-priori error bounds for a class of nonlinear electrical circuits, in: Proc. of CDC/CCC 2009, 2009, pp. 5376–5383. 30. M. Meyer, H. Matthies, Efficient model reduction in non-linear dynamics using the Karhunen-Lo`eve expansion and dual-weighted-residual methods, Comput. Mech. 31 (2003) 179–191. 31. M. Drohmann, B. Haasdonk, M. Ohlberger, Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation, Preprint Angewandte Mathematik und Informatik 02/10-N, University of Mnster, Submitted to SISC (2010). 32. N. Jung, A. Patera, B. Haasdonk, B. Lohmann, Model order reduction and error estimation with an application to the parameter-dependent eddy current equation, Math. Comp. Model. Dyn., Accepted. 33. J. Hale, Ordinary differential equations, Wiley-Interscience, 1969.
Appendix A Detailed proofs A.1 Proof of Lemma 1 Proof (Proof ) It is clear that φ0 (r0 ) ≤ Ss (r) ≤ 0 ∀ r ∈ R+ 0 . Using the boundedness of φ we see that φ(r) − φ(s) ≤ B + φ(s) r→∞ |Ss (r)| ≤ −→ 0. r−s |r − s| As Ss is continuous and Ss (r0 ) < 0 ∀ s we see that rs < ∞ and hence obtain well-definedness of (16). 0 s (rs ) , so Next we show the position inequalities. Necessary condition for a minimizer is 0 = Ss0 (rs ) = φ (rsr)−S s −s 0 φ (rs ) = Ss (rs ). At first consider the case s < r0 . We continue with proof by contradiction and therefore assume φ(r)−φ(s) we had rs < r0 . Now choose a r ∈ R+ = φ0 (ξ) for 0 with max(rs , s) < r < r0 . Then we have Ss (r) = r−s some ξ ∈]s, r[ by the mean value theorem. Now as φ is strictly concave on [0, r0 ] by the bell function condition (14) and φ0 (r) ≤ 0 by (13) we obtain Ss (rs ) = φ0 (rs ) > φ0 (ξ) = Ss (r), which contradicts with the minimizing property of rs as obviously rs is not a minimizer of Ss . The case r0 < s is obtained by an analogous argument using the strict convexity of φ on [r0 , ∞[. Finally, the case r0 = s is clear when considering the limit process for both former cases. Now, uniqueness of the minimizers follows directly using the necessary condition for minimizers and the concavity/convexity ”on the other side“ as e.g. for two assumed minimizers rs 6= rs0 we always obtain a contradiction Ss (rs ) = φ0 (rs ) 6= φ0 (rs0 ) = Ss (rs0 ) to the minimal choice. Having s = r0 directly gives uniqueness by (19).
A.2 Auxiliary lemma for the proof of Proposition 1. + 0 Lemma 2 Let s ∈ R+ 0 \{r0 } and Ss be given as in (15). Then Ss (r)(r − rs ) > 0 ∀ r ∈ R0 \{rs }.
14
D. Wirtz and B. Haasdonk 0
φ (r)−Ss (r) 0 Proof (Proof ) Choose s ∈ R+ . First consider the case s < r0 . Then we 0 and note that Ss (r) = r−s differentiate locations of r: – r < s < r0 ≤ rs : We have r − rs < 0, r − s < 0 and Ss (r) > φ0 (r) by concavity. 0 0 φ00 (t)−Ss (t) s (t) – r = s < r0 ≤ rs : We have Ss0 (s) = lim Ss0 (t) = lim φ (t)−S = lim = φ00 (s) − lim Ss0 (t). This t−s 1 t→s
00
t→s
t→s
t→s
means Ss0 (s) = φ 2(s) < 0 and as r < rs we have Ss0 (r)(r − rs ) > 0. – s < r ≤ r0 ≤ rs : We have r − rs < 0, r − s > 0 and Ss (r) < φ0 (r) by concavity. – s < r0 ≤ r < rs : Since φ0 (r0 ) < Ss (r0 ) and φ0 (rs ) = Ss (rs ) we must have φ0 (r) < Ss (r) ∀ r ∈ [r0 , rs [. Otherwise, since φ0 (·) − Ss (·) is continuous, there would exist an r0 ∈ [r0 , rs [ with Ss (r0 ) = φ0 (r0 ) < φ0 (rs ) = Ss (rs ) as r0 < rs and φ is strictly convex on [r0 , ∞[. But this is a contradiction to the minimal choice of rs , and together with r − rs < 0, r − s > 0 we obtain Ss0 (r)(r − rs ) > 0. 0 00 φ00 (rs )−2Ss (rs ) s) – s < r0 ≤ rs < r: At rs we have Ss0 (rs ) = 0 and thus Ss00 (rs ) = = φrs(r > 0, as φ is rs −s −s 0 strictly convex for r > r0 . So there exists an > 0 with Ss (rs + ) > 0. In fact for any point r > r0 with Ss0 (r) = 0 we have Ss00 (r) > 0. But as Ss0 is continuous there cannot be two consecutive local minima without a local maximum. Hence, Ss0 (r) cannot equal zero for r > rs and since Ss0 (rs + ) > 0 we must have Ss0 (r) > 0 ∀ r > rs . The case r0 < s is obtained using an analogous argumentation.
A.3 Proof of Theorem 3 Proof (Proof ) At first we see that Γn is open (in fact a union of open intervals) since di (t), Cn (t) and rdi (t) are continuous in t. PN For n ∈ N let βn (t) := i=1 ||ci ||G |Sdi (t) (rCn (t),di (t) )| and si := di (t) as shorthand for any fixed t ∈ [0, T ]. We prove (26) by induction and n begin with n = 0. We haveoΓ1 6= ∅ by assumption. So, for t ∈ Γ1 condition (25) implies that M := i ∈ {1 . . . N } rsi ∈ / BC1 (t) (si ) is nonempty. Further, we know that |Ssi si + sign (rsi − si ) C1 (t) | < |Ssi (rsi )| ∀ i ∈ M as Lemma 2 holds since si 6= r0 ∀ i ∈ M . Then we deduce X X β1 (t) = ||ci ||G Ssi (rC1 (t),si ) + ||ci ||G Ssi (rC1 (t),si ) i∈M (20)
=
X
i∈M /
X ||ci ||G Ssi si + sign (rsi − si ) C1 (t) + ||ci ||G Ssi (rsi ) i∈M /
i∈M
inf Γ1 we know with (28) that Z
Rs
Z
β1 (r)dr
α(s)e0
ds
inf Γn , i.e. Cn+1 (t) < Cn (t) ∀ t > inf Γn by definition. This means Γn+1 6= ∅ as Γn ⊆ Γn+1 . So, for any t ∈ Γn+1 we have rCn+1 (t),si = si − Cn+1 (t) > si − Cn (t) = rCn (t),si which gives βn+1 (t) < βn (t) analogous to (28). Splitting the integral as above using Γn+1 directly gives (26), concluding the induction. n If Γ1 = ∅ we have βn ≡ β1 ∀ n ≥ 1 and thus ∆C LSLE ≡ C1 ∀ n ≥ 1. n Next, for each t ∈ [0, T ] we see from (26) that ∆C LSLE (t) n∈N is monotonically decreasing and bounded by Cn ∞ n zero. Thus, we have a pointwise convergence of ∆C LSLE to the limit function ∆LSLE (t) := limn→∞ ∆LSLE (t). From for ∀ n, m ∈ N, m ≥ n we have βm (t) ≤ βn (t) ∀ t ∈ [0, T ]. This the inductionC we know that Cn m n Cm implies ∆C LSLE (t) − ∆LSLE (t) ≤ ∆LSLE (T ) − ∆LSLE (T ) ∀ t ∈ [0, T ]. Taking the supremum on the left and considering m → ∞ we obtain n→∞ ∞ Cn Cn ≤ ∆∞ ∆LSLE − ∆LSLE LSLE (T ) − ∆LSLE (T ) −→ 0, L∞ (0,T )
which finally shows uniform convergence and continuity of ∆∞ LSLE (T ).
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
15
B Penalized Newton iteration We first state some preparatory lemmata before we propose our penalized Newton iteration in Proposition 3. Lemma 3 Let φ be a bell function. Then the mapping + η(s) : R+ 0 −→ R0 , s 7−→ rs = arg min Ss (r) is monotonically decreasing. r∈R+ 0
Proof (Proof ) Note here that η is well-defined by Lemma 1. Choose some s1 , s2 ∈ R+ 0 with s1 6= s2 . Without loss of generality assume s1 < s2 . Let ri := η(si ), i ∈ {1, 2}. At first consider the case s1 < s2 ≤ r0 . Define l(r) := φ(s1 ) + φ0 (r1 )(r − s1 ). Then we see that l(s2 ) ≤ φ(s2 ), as φ(s2 ) − φ(s1 ) φ(r1 ) − φ(s1 ) l(s2 ) − φ(s1 ) ≥ = φ0 (r1 ) = . s2 − s1 r1 − s1 s2 − s1 Using φ0 (r1 ) =
φ(r1 )−φ(s1 ) , r1 −s1
l can also be written as l(r) = φ(r1 ) + φ0 (r1 )(r − r1 ). So we conclude
φ(r1 ) − φ(s2 ) φ(r2 ) − φ(s2 ) ≤ r2 − s2 r1 − s2 φ(r1 ) − l(s2 ) φ(r1 ) − φ(s1 ) ≤ = = φ0 (r1 ), r1 − s2 r1 − s1 and since φ(x) is convex for x > r0 we must have r1 ≥ r2 . The proof for the case r0 ≤ s1 < s2 is analogous using the concavity on [0, r0 ]. φ0 (r2 ) =
Lemma 4 (Bounds of rs ) Let φ be a bell function. Then we have 0 ≤ rs < rm ∀ s ∈ R+ for rm :=
φ(0)r0 . φ(0) − φ(r0 )
Proof (Proof ) Fix an s ∈ R+ . We already know that 0 ≤ rs by definition. At first we see that rm = φ(0) r0 φ(0)−φ(r > r0 φ(0) = r0 . So, if s > r0 we have rs ≤ r0 < rm by Lemma 1, (18). Obviously, for s = r0 we φ(0) 0) have rs = r0 < rm . This leaves us with the case s < r0 , i.e. r0 ≤ rs by (17). Let ls (r) := φ(r0 ) + Ss (r0 )(r − r0 ) be the line through φ(r0 ) with gradient Ss (r0 ). Since Ss (r0 ) < 0 and ls (r0 ) = φ(r0 ) > 0 we know that 0 ∃! ηs ∈ R+ 0 : ls (ηs ) = 0. Also, as 0 > Ss (r0 ) > φ (r0 ) we have ∃ δ > 0 : φ(r0 + δ) < ls (r0 + δ). But as φ(ηs ) > 0, ls (ηs ) = 0 and φ − ls is continuous we know that ∃ t ∈ ]r0 , ηs [ : φ(t) = ls (t). Further we have Sr0 (t) =
φ(t) − φ(r0 ) ls (t) − φ(r0 ) φ(r0 ) + Ss (r0 )(t − r0 ) − φ(r0 ) = = = Ss (r0 ). t − r0 t − r0 t − r0
Since φ is convex on [r0 , ∞[, ∃! ξ ∈]r0 , t[: φ0 (ξ) = Sr0 (t). With the minimizing property of rs we see that φ0 (ξ) = Sr0 (t) = Ss (r0 ) > Ss (rs ) = φ0 (rs ), but this means rs < ξ by the strict convexity of φ on [r0 , ∞[. Finally, as s ∈ R+ ∀ s ∈ R+ 0 was arbitrary, we obtain rs < ξ < t < ηs 0 . Now Lemma 3 yields rs ≤ rs=0 = η(0) < η0 ∀ s ∈ R+ , and solving l (η ) = 0 gives the desired bound rm = η0 . 0 0 0 φ(0)r0 . Define the φ(0)−φ(r0 ) n(ν) 2 n0 (ν) − ν. Then rs =
Proposition 3 (Penalized Newton iteration) Let φ be a bell function and rm := unpenalized objective function as n(r) := φ0 (r) − Ss (r) and c(ν) := arg min φ(r)−φ(s) r−s r∈R+ 0
n0 (ν)2 , d(ν) 4n(ν)
:=
can be computed as root of the penalized objective function np given by
s ≤ r0 :
r0 < s :
2 r ≤ r0 c(r0 )(r + d(r0 )) r0 < r < rm np (r) := n(r), c(rm )(r + d(rm ))2 , r ≥ rm 2 r≤0 c(0)(r + d(0)) , 0 < r < r0 np (r) := n(r), c(r0 )(r + d(r0 ))2 , r ≥ r0
Proof (Proof ) At first consider the penalty polynomials and define pν (x) = c(ν)(x + d(ν))2 . Then it is easy to validate that pν (ν) = n(ν), p0ν (ν) = n0 (ν). Those polynomials will be used to extend n at appropriate points. Next, existence and uniqueness of rs have already been shown in Lemma 1. At first consider the case s < r0 . 0 s (r) , and since s < rs by Lemma As the necessary condition for a minimizer r ≥ r0 we have 0 = Ss0 (r) = φ (r)−S r−s 0 0 1, 0 = n(r) = φ (r) − Ss (r) is also sufficient for r ≥ r0 . Since n(s) = φ (s) − Ss (s) = 0, we replace n(r) for r < r0 by pr0 (r). Even though Ss0 (r) > 0 ∀ r > rs by Lemma 2 and hence n(r) 6= 0 ∀ r > rs we see that limr→∞ n(r) = 0. To avoid the Newton iteration being drawn to ∞ we replace n(r) by prm for r ≥ rm as we know that rs < rm ∀ s ∈ R+ 0 from Lemma 4. For s > r0 we have 0 ≤ rs and so we replace n(r) by p0 for r ≤ 0 and n(r) by pr0 for r0 ≤ r since again n(s) = 0. Of course, s = r0 means rs = r0 and no iterations have to be performed.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems Stuttgart, December 2010 (rev. June 2011) Institute of Applied Analysis and Numerical Simulation, University of Stuttgart, Pfaffenwaldring 57 D-70569 Stuttgart, Germany {daniel.wirtz,bernard.haasdonk}@mathematik.uni-stuttgart.de www.agh.ians.uni-stuttgart.de
Abstract In this paper we consider the topic of model reduction for nonlinear dynamical systems based on kernel expansions. Our approach allows for a full offline/online decomposition and efficient online computation of the reduced model. In particular we derive an a-posteriori state-space error estimator for the reduction error. A key ingredient is a local Lipschitz constant estimation that enables rigorous a-posteriori error estimation. The computation of the error estimator is realized by solving an auxiliary differential equation during online simulations. Estimation iterations can be performed that allow a balancing between estimation sharpness and computation time. Numerical experiments demonstrate the rigour and effectiveness of the error bounds. Keywords nonlinear dynamical systems · model reduction · kernel methods · a-posteriori error estimates · offline/online decomposition · subspace projection
Preprint Series Issue No. 2010-83, rev. June 2011 Stuttgart Research Centre for Simulation Technology (SRC SimTech) SimTech – Cluster of Excellence Pfaffenwaldring 7a 70569 Stuttgart [email protected] www.simtech.uni-stuttgart.de
2
D. Wirtz and B. Haasdonk
1 Introduction Nowadays, modeling of real world processes like biochemical reactions or electric circuits naturally leads to a formulation as dynamical systems with inputs and outputs. Mathematically, they are described by systems of ordinary differential equations (ODEs) and can be roughly categorized into linear and nonlinear types. Even though computational power has significantly increased over the last years, high-resolution models often result in large-scale dynamical systems that are expensive to simulate. In this context the need for fast simulation is particularly evident in many-query scenarios. They comprise parameter studies or inverse problems where multiple simulations have to be performed for different inputs or initial state configurations. Additionally, some dynamical systems model processes in a real-time setting like control components and thus also require fast computation without strong hardware. Within all those settings fast and rigorous error estimation procedures are important in order to quantify the model errors introduced by the reduction process. Consequently, model reduction techniques for dynamical systems are nowadays subject to intensive research. Work has been done for various types of linear systems, e.g. time invariant systems [1], time-variant [2] or parameterized [3]. For an overview we refer the reader to [4–6]. Reduction techniques for nonlinear dynamical systems have been less investigated, not at least because of the arising difficulties. As with linear systems, most of them involve projection of the system into a lower dimensional subspace. One well known method for nonlinear model reduction is the trajectory piecewise linear (TPWL) approach [7], which is extended by moment matching techniques in [8] or to a piecewise-polynomial scheme in [9]. Various extensions of the balanced truncation procedure [10] to nonlinear systems are investigated in [11–13]. An “approximate reduction” method is introduced in [14] and model reduction for weakly nonlinear systems by bilinearization has been discussed in [15], for example. In addition, the special class of bilinear quadratic nonlinearities has recently been investigated in [16]. A further promising technique for nonlinear model reduction of dynamical systems is by means of discrete empirical interpolation [17]. Finally, computing subspaces by proper orthogonal decomposition (POD) of presumably statistically representative trajectories [18] is an expensive but well-established method. This method is also known as principal component analysis (PCA), see [19] and references therein for an overview. In this work we adopt principles of the reduction technique proposed in [20], extend the class of usable kernels and provide novel, efficient a-posteriori error estimators for the resulting reduced systems. The structure of our paper is as follows: In Section 2 we introduce our base dynamical system and discuss aspects of the reduction process. Section 3 is concerned with a-posteriori error estimates and introduces our new local Lipschitz constant estimation method. Next, Section 4 presents numerical experiments for synthetic dynamical systems and we conclude with Section 5. Some auxiliary mathematical details are presented in the appendix. 2 Reduction of kernel based systems The starting point for our investigations is the reduction method introduced in [20], which combines subspace projection with kernel methods. In model reduction, the latter are used to approximate the nonlinearities of dynamical systems in kernel spaces and promise a large potential in the field. Kernel methods comprise applications from machine learning like support vector regression [21] as well as kernel interpolation with corresponding theoretical foundation [22–24]. However, those approximation techniques are out of the scope of this article. Instead we focus on nonlinear systems whose inner dynamics are already given by a kernel expansion and investigate the projection behaviour and the resulting a-posteriori error estimators. 2.1 The base dynamical system Our central assumption is to have a nonlinear kernel expansion f (x) =
N X i=1
ci Φ(x, xi ),
(1)
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
3
where X := {x1 , . . . , xn } ⊂ Rd are the centers or support vectors of the expansion and ci ∈ Rd the coefficient vectors. The function Φ is a kernel, which is basically a symmetric function Φ : Rd ×Rd −→ R. We omit further characterizations of kernels and mention here that the key analytical property is that for some open Ω ⊂ Rd each symmetric positive definite kernel spans a unique Hilbert space NΦ (Ω) of functions. Those Hilbert spaces have a special reproducing property and are commonly referred to as reproducing kernel Hilbert spaces (RKHS). In our context those RKHS serve as a base class for the d dynamical systems functions f , so (1) implies f ∈ (NΦ (Ω)) . For more details on kernels, RKHS and their applications we refer to [22, 25, 26], for example. Now, the base system considered throughout this paper is of the form x0 (t) = f (x) + Bu(t), x(0) = x0 , y(t) = Cx(t),
(2) (3)
where x(t) ∈ Rd denotes the state of the system at times t ∈ [0, T ], 0 ≤ T < ∞, x0 ∈ Rd the initial condition, u : [0, T ] → Rm an input function with B ∈ Rd×m and output y(t) ∈ Rk with C ∈ Rk×d . Additionally, in some applications different norms than the standard Euclidean norm are taken for the state space. Thus, let G ∈ Rd×d be a symmetric positivepdefinite matrix. Then G defines a scalar product hx, yiG := xt Gy on Rd with induced norm ||x||G := hx, xiG . Choosing G := Id will yield the standard Euclidean inner product h·, ·i and 2-norm ||·||, where Id denotes the d-dimensional identity matrix. For example, G is usually chosen as the Gram-matrix of the finite element/ finite volume basis for dynamical systems obtained from partial differential equation discretizations. Having introduced the underlying dynamical system (2)-(3) we now detail the reduction approach using subspace projection. We assume to have two biorthogonal projection matrices V, W ∈ Rd×r , W t V = Ir with r d denoting the reduced system’s dimension. The reduced system is then obtained by applying a Petrov-Galerkin projection to the full system: z 0 (t) = W t f (V z(t)) + W t Bu(t), z(0) = W t x0 , y r (t) = CV z(t),
(4) (5)
where z(t) ∈ Rr now denotes the reduced system’s state and y r (t) the approximate output. Let U := hv1 , . . . , vr i ⊂ Rd be the space spanned by the columns vi ∈ Rd of V . Consequently, the reconstructed approximate solution and output are given by xr (t) := V z(t) ∈ U and y r (t) = Cxr (t). The requirements on V, W are now to yield a good approximation xr (t) ≈ x(t) and thus y r (t) ≈ y(t) for different initial values and inputs. For the remainder of this work we will assume the matrix computation method for V, W as a black-box since our method is applicable using any basis/projection matrix generation method. The rest of this section is concerned with the different projection aspects of the reduction technique for kernel-based dynamical systems. We will see that the structure of the kernel expansion allows for dramatic reduction of the computational costs using certain kernels.
2.2 System projection At first, the projection into Rr by W t can be applied directly to the coefficient vectors ci via cri := PN W t ci , i = 1 . . . N . This results in a reduced function f r (z) := i=1 cri Φ(V z, xi ), with new coefficient r r vectors ci ∈ R . Unfortunately, the reduced system (4) is still expensive to simulate given that f is evaluated at xr (t) ∈ Rd . In order to avoid input arguments of high dimension d, we consider two special classes of kernels. 2.2.1 Inner product kernels The first type of kernels that allow efficient argument evaluations are the inner product kernels also mentioned in [20]. It is assumed that Φ(x, y) = φ(hx, yiG ) for some scalar function φ : R → R. Then Φ(V z, xi ) = φ(hV z, xi iG ) = φ(hz, zi i) =: Φr (z, zi ) for zi := V t Gxi ∈ Rr , i = 1 . . . N . This way, it is sufficient to project the centers xi into the low-dimensional subspace and the evaluations of Φ can
4
D. Wirtz and B. Haasdonk
be computed efficiently and loss-less during the reduced simulation via Φr . Some examples for those kernels are the linear kernel Φ(x, y) = hx, yiG or polynomial kernels Φ(x, y) = (1 + hx, yiG )p of degree p ∈ N. 2.2.2 Translation- and rotation invariant kernels In extension to [20] we present here a further class of kernels that allow for efficient loss-less argument evaluations: The translation and rotation-invariant kernels Φ(x, y) := φ(||x − y||G ) for some scalar function φ : R+ 0 → R. We impose the additional requirement xi ∈ U, i.e. xi = V zi for some zi ∈ Rr , i = 1 . . . N . Then we obtain Φ(V z, xi ) = φ(||V z − V zi ||G ) = φ(||z − zi ||V t GV ) =: Φr (z, zi ) with V t GV being a small Rr×r matrix inducing a new norm on Rr . Note that the assumption xi ∈ U is of a technical nature. We either extend U by the span of the xi , or, if the kernel expansion is created with knowledge of U , one can choose xi ∈ U in the first place. Those kernels are also commonly referred to 2 as radial basis functions, of which the Gaussian kernel Φ(x, y) = exp(− ||x − y|| /β 2 ) is probably the most popular example. For further examples and characterizations of the above mentioned kernels we refer to [21, 26].
2.3 The reduced system In conclusion, for any of the specific kernels above we obtain from (4)-(5) an equivalent reduced system: z 0 (t) =
N X
cri Φr (z, zi ) + B r u(t),
(6)
i=1
z(0) = z0 ,
y r (t) = C r z(t),
(7)
with z0 := W t x0 , B r := W t B, C r := CV . The evaluation complexity is independent of d and the system has only low dimensional components. Note again, that this reduction technique (using inner-product kernels) has been investigated in [20]. Despite our new class of kernels, investigation of the reduction properties is out of scope of this article. Instead, our main focus is the error estimation.
3 A-posteriori error estimates When dealing with reduced dynamical systems it is essential to have means to compute or at least estimate the error that is made during simulations. Still, little work has been done so far in the field of error estimations for nonlinear dynamical systems. A-priori error analysis of POD reduction schemes in different applications are considered e.g. in [27, 28] and a-priori error bounds for certain nonlinear systems are derived in [29]. Error estimates using dual-weighted-residual methods can be found in [30] and [7] perform a-posteriori error estimation for their TPWL approach. Recent work on a-posteriori error estimation has also been done for nonlinear parametrized evolution equations [31] and a related scheme for error estimators of linear but parameterized dynamical systems can be found in [3, 32]. Note that in many cases differential inequalities such as the comparison lemma or Gronwall’s lemma [33, p.36] are used to bound ODE solutions/errors a-priori. Instead, we will use it to provide a-posteriori error bounds which can be computed along with the reduced system for any input u(t). Our estimators bound the error in the state-space for finite time intervals instead of giving frequency domain error bounds. For the reduced system (6)-(7) we derive rigorous and efficient error estimators. In Section 3.1 we introduce the estimator structure obtained using the comparison lemma [33, p.32] (also applied in [7]). Section 3.2 introduces the local Lipschitz constant estimations utilized in Section 3.3 to obtain our main estimator result. We complete our results with convergence analysis and computational aspects in Sections 3.4 and 3.5.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
5
3.1 Global Lipschitz constant a-posteriori error estimator The state space error e(t) := x(t) − xr (t) for any t ∈ [0, T ] is given by the solution of the error system e0 (t) = f (x(t)) + Bu(t) − V W t f (xr (t)) − V W t Bu(t), e(0) = x0 − V W t x0 .
(8) (9)
Estimating the norm of the derivative and application of the comparison lemma yields a first aposteriori error estimator. Let Ex0 := ||(Id − V W t )x0 ||G be the initial error for the remainder of this work. Theorem 1 (Global Lipschitz constant a-posteriori error Estimator (GLE)) Let Φ be a kernel. Further, let Φ be LΦ -Lipschitz continuous with respect to the first variable. Then the state space error is bounded via ||e(t)||G ≤ ∆GLE (t) ∀ t ∈ [0, T ] with Zt ∆GLE (t) :=
α(s)eβ(t−s) ds + Ex0 ,
0 α(t) := Id − V W t f (xr (t)) + Bu(t) G , β := LΦ
N X
(10)
||ci ||G .
i=1
Proof (Proof ) With (1) we directly obtain N X ||f (x(t)) − f (xr (t))||G = ci Φ(x(t), xi ) − Φ(xr (t), xi ) i=1
≤
N X
(11) G
||ci ||G Φ(x(t), xi ) − Φ(xr (t), xi )
i=1
≤ LΦ
N X
||ci ||G ||e(t)||G = β ||e(t)||G .
i=1
Now from (8) we get
||e0 (t)||G = f (x(t)) − f (xr (t)) + Id − V W t f (xr (t)) + Bu(t) G ≤ β ||e(t)||G + α(t).
Finally, application of the comparison lemma [33, p.32] for the ODE v 0 (t) = βv(t) + α(t), v(0) = Ex0 yields the error estimator as analytical solution. Note that the error estimator ∆GLE (t) is rigorous, continuous and monotonically increasing. Once a bound for the state error is available the output error ey (t) := y(t) − y r (t) can be bounded by ||ey (t)|| ≤ ||C||G ||e(t)||G ≤ ||C||G ∆GLE (t), t ∈ [0, T ], where ||C||G is the G-induced matrix norm. 3.2 Local Lipschitz constants Closer investigation of the error estimator from Theorem 1 shows that if the kernel Lipschitz constant LΦ is large, the estimator will grow fast and the estimation will be overly conservative. Therefore, we propose a local Lipschitz constant computation that utilizes the current position xr (t) of the reduced state variable and an a-priori error bound ||e(t)||G ≤ C(t) to improve the estimation procedure. The benefit of knowing xr (t) is that one argument of the estimation ||f (x(t)) − f (xr (t))||G ≤ β ||x(t) − xr (t)||G is fixed, which will be used to compute a local Lipschitz constant using secant gradients. Also, the a-priori bound C(t) will enable efficient computation of those secant gradients utilizing the kernel expansion center positions.
6
D. Wirtz and B. Haasdonk
The key to our local Lipschitz constant computations is to use the class of kernels introduced in Section 2.2.2. In particular, we use radial basis functions Φ(x, y) := φ(||x − y||), x, y ∈ Rd induced by a certain class of scalar functions φ: + Definition 1 (Bell functions) A function φ : R+ 0 → R0 is called a bell function, if
i) ii) iii) iv)
φ ∈ C 2 (R+ 0 ), ||φ||L∞ ≤ B, B > 0, φ0 (0) ≤ 0, ∃ r0 > 0 : φ00 (r)(r − r0 ) > 0 ∀ r 6= r0 .
(12) (13) (14)
Here, r0 denotes the unique turning point where φ is strictly concave/convex on [0, r0 ]/[r0 , ∞[, re2 spectively. For kernel is induced by a bell function φ(r) = e−(r/β) . Then we √ example, the Gaussian have r0 = β/ 2 and LΦ = |φ0 (r0 )| in the context of Theorem 1. Let φ denote a bell function for the remainder of this work. We will state results regarding bell functions first and connect these to the multidimensional case later. + For any s ∈ R+ 0 , we introduce a secant gradient function Ss : R0 −→ R with ( φ(r)−φ(s) , r 6= s r−s Ss (r) := (15) 0 φ (s), r = s. Well-definedness and continuous differentiability of Ss can easily be verified by elementary analysis. The following lemma characterizes the minimum secant gradient position. Lemma 1 (Minimum secant gradient) Let s ∈ R+ 0 . Then rs := arg min+ Ss (r) r∈R0
(16)
exists, is unique and exactly one of the following cases holds: s < r0 ≤ r s , rs ≤ r0 < s, s = r0 = rs .
(17) (18) (19)
For the proof we refer to SectionA.1 of the appendix. Figure 1 shows two examples for minimum secant gradients and the position of rs with each s < r0 (left) and s > r0 (right), using the Gaussian inducing bell function φ(r) = exp(−r2 /β 2 ) with β = 2. Now, further assume to restrict admissible r values to r ∈ BC (s) ∩ R+ 0 for a C > 0, where BC (s) denotes an open ball of radius C around s and BC (s) the closure. With this condition we can state our novel result regarding local Lipschitz constant computations. Proposition 1 (Local Lipschitz estimation using secant gradients) Let C > 0, s ∈ R+ 0 and Ω := BC (s) ∩ R+ . Then |φ(r) − φ(s)| ≤ |S (r )||r − s| ∀ r ∈ Ω with s C,s 0 ( s + sign (rs − s) C rs ∈ /Ω rC,s = (20) rs rs ∈ Ω. Proof (Proof ) Assume s 6= r0 . Then Lemma2 in SectionA.2 of the appendix implies that Ss is monotonically decreasing on [0, rs ] and increasing on [rs , ∞[. For rs ∈ / Ω, as Ss is negative, |Ss | takes its maximum value at the border of Ω that is closest to rs . Since arg minr∈R+ Ss (r) = arg maxr∈R+ |Ss (r)|, 0 0 case rs ∈ Ω follows by definition. The case s = r0 implies rs = r0 and hence case two, as of course r0 ∈ Ω = BC (r0 ) and C > 0. We finally conclude |φ(x) − φ(y)| ≤ ||Ss ||L∞ (Ω) |r − s| = |Ss (rC,s )||r − s|. Next, we generalize these results to the multidimensional kernel setting.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems Bell function demo, s=0.282843, rs=2.058702 1.2
Bell function demo, s=5.391323, rs=0.400439 1.2
φ φ(rs) + φ’(rs)(rs−s)
φ(s)
1
7
0.8
φ φ(rs) + φ’(rs)(rs−s)
φ(r ) s
1
0.8
φ(r0)
0.6
φ(rs)
0.4
0.4
0.2
0.2
φ(rm)
0
r0
s
φ(r0)
0.6
rs
φ(rm)
0
rm
−0.2
rs
r0
φ(s) s
rm
−0.2
−0.4
−0.4 0
1
2
3
4
5
Fig. 1 Gaussian φ with β = 2, r0 =
√
6
2 and rm =
7
√
0
2
1−e−1/2
1
2
3
4
5
6
7
≈ 3.5942
Corollary 1 (Local Lipschitz constants for bell function kernels) Let x, y, z ∈ Rd , C > 0, Φ(x, z) = φ(||x − z||G ) and s = ||y − z||G . Then we have |Φ(x, z)−Φ(y, z)| ≤ |Ss (rC,s )| ||x − y||G ∀ x ∈ BC (y). Proof (Proof ) Fix y, z ∈ Rd , C > 0, let R := x ∈ Rd ||x − z||G ∈ BC (s) . Then using Proposition 1 we estimate |Φ(x, z) − Φ(y, z)| = |φ(||x − z||G ) − φ(||y − z||G )| ≤ |Ss (rC,s )| ||x − z|| − ||y − z|| G
G
≤ |Ss (rC,s )| ||x − y||G ∀ x ∈ R.
The fact that BC (y) ⊂ R finishes the proof.
3.3 Improved error estimator Finally, the local Lipschitz constant estimations from Section 3.2 can be used to state our main result, which is an improved version of the error estimator derived in Theorem 1. As shortly mentioned in the beginning, we additionally assume to have an a-priori coarse error bound ||e(t)||G ≤ C(t) ∀ t ∈ [0, T ],
(21)
which can be C(t) = ∞ if no further knowledge is available. We introduce the notation di (t) := ||xr (t) − xi || ,
i = 1 . . . N,
(22)
for the distance of the reduced state xr (t) to each expansion center xi during the reduced simulation. In fact, the coarse error bound means nothing but x(t) ∈ BC(t) (xr (t)) ∀ t ∈ [0, T ]. Consequently, for any t ∈ [0, T ], i ∈ {1 . . . N } we identify x, y, z ∈ Rd and C > 0 from Corollary 1 with x(t), xr (t) of the full/reduced solution, the kernel expansion centers xi and C(t), respectively. This allows to obtain local Lipschitz constant estimations at each time t using the current reduced simulation’s state xr (t) as “center” of locality: Theorem 2 (Local Secant gradient Lipschitz a-posteriori error Estimator (LSLE)) Let the error system be given as in (8)-(9), where the kernel Φ of the expansion (1) is induced by a bell
8
D. Wirtz and B. Haasdonk
function φ. Further, let C(t) be an a-priori error bound. Then the state space error is bounded via ||e(t)||G ≤ ∆C LSLE (t) ∀ t ∈ [0, T ] with Zt ∆C LSLE (t)
:=
Rt
β(r)dr
α(s)es
ds + Ex0 ,
0
α(t) := Id − V W t f (xr (t)) + Bu(t) G , β(t) :=
N X
||ci ||G Sdi (t) (rC(t),di (t) ) .
(23) (24)
i=1
Proof (Proof ) Similar to (11), with (1) and Corollary 1 we deduce that ||f (x(t)) − f (xr (t))||G ≤ β(t) ||e(t)||G . Performing an analogous estimation as in Theorem 1 we obtain the comparison lemma ODE v 0 (t) = β(t)v(t) + α(t), v(0) = Ex0 , whose analytical solution yields the error estimator. We see in the α(t) term that if any input u(t) happens to maintain the full system in the subspace U , then the resulting zero approximation error will be verified a-posteriori by both the GLE and LSLE error estimators. Before we present numerical experiments we complete the estimator results with some theory regarding the coarse bound and computational aspects.
3.4 Utilization of a-priori bound C(t) We still need to justify the assumption (21) because it requires existence of an a-priori bound in order to obtain an a-posteriori error estimator. The core idea is to use the a-priori bound in an iterative scheme starting with C(t) ≡ ∞, which is trivially satisfied. Then the resulting a-posteriori error estimate is reused as a-priori bound for another iteration of the error estimation procedure. A nice feature of this procedure is that it allows to balance the error bound’s sharpness against the computational costs. The next theorem shows that actual progress can be made under reasonable assumptions and the iterative scheme reproduces the first estimate in the worst case. Theorem 3 (Convergence of error estimator iterations) Let the conditions of Theorem 2 hold Cn−1 and define C0 (t) ≡ ∞ and Cn (t) := ∆LSLE (t), n ≥ 1. Recall the distances di (t) from (22), rs from (16) and for n ∈ N let Γn := t ∈ [0, T ] ∃ i ∈ {1 . . . N } : rdi (t) ∈ / BCn (t) (di (t)) . (25) If Γ1 6= ∅, then Γn+1 = 6 ∅ ∀ n ∈ N, ( n < ∆C Cn+1 LSLE (t), ∆LSLE (t) n = ∆C LSLE (t),
t > inf Γn+1 t ≤ inf Γn+1
∀ n ∈ N.
(26)
n Else if Γ1 = ∅, then ∆C LSLE (t) = C1 (t) ∀ n ≥ 1, t ∈ [0, T ]. Further, the sequence of iterates converges Cn uniformly to a continuous lower error bound ∆∞ LSLE (t) := lim ∆LSLE (t), t ∈ [0, T ].
n→∞
For the proof we refer to Section A.3 of the appendix. The crucial assumption here is that the first estimation C1 is small enough, so that for some time t and expansion center xi the minimum secant gradient point rdi (t) is not contained in BC1 (t) (di (t)). This is reasonable in the context of kernel-based model reduction, since the expansion centers are naturally scattered over the high dimensional state space.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
9
3.5 Computation of error estimators The error estimators from Theorem 1 and 2 allow an efficient computation during the reduced simulation with a complexity independent from the full system’s (high) dimension. As already mentioned in the proofs of the respective theorems, the error estimators are given by the solution of the onedimensional ODE u0 (t) = β(t)v(t)+α(t) with initial condition u(0) = Ex0 (For the GLE from Theorem 1 we have a constant β). This solution can be computed “on the fly” by adding an extra dimension to the reduced system. The efficient computation is possible due to the full offline/online separability of the error estimators and an effective scheme for solving the minimum secant gradient problem, which we will present next. 3.5.1 Offline/online decomposition A vital ingredient to the efficient computation of the error estimators is to decompose the α(t) ODE parts (10),(23) in an offline/online fashion similar to [3]. In this setting, the offline phase comprises expensive computations with complexity depending on the high dimension d and has to be performed once. The online phase is performed for different inputs u(t) and its complexity is dependent only on the reduced dimension r. Here, we confine ourselves to detail the case for the LSLE presented in Theorem 2, as the GLE from Theorem 1 is decomposed the same way but with constant β(t) term. Proposition 2 (Offline/online decomposition of the LSLE) Let the conditions from Theorem 2 hold. Define M := (α1 . . . αN ) ∈ Rd×N and recall the norm inducing matrix G from Section 2.1. Then the offline computations are Ex0 = ||x0 − V W t x0 ||G , ||ci ||G , i = 1 . . . N for the β(t) term and ˜ := Id − V W t M, ˜ := Id − V W t B, M B ˜ t GM ˜ ∈ RN ×N , ˜ t GB ˜ ∈ RN ×m , ˜ t GB ˜ ∈ Rm×m M1 := M M2 := M M3 := B regarding α(t). With v(0) = Ex0 , the online part consists of computing q t t v 0 (t) = β(t)v(t) + ϕ(t) M1 ϕ(t) + 2ϕ(t) M2 u(t) + u(t)t M3 u(t), N ϕ(t) := Φr (z(t), zi ) i=1 ∈ RN .
(27)
The proof follows directly from the definitions. Notice that if there are no inputs u(t) for a given system only the matrix M1 has to be computed as M2 and M3 do not appear. However, all the offline matrices M1 , M2 , M3 are small matrices independent of the original systems dimension. 3.5.2 Efficient computation of secant gradient Lipschitz constants When using the improved error estimation from Theorem 2, another important aspect for fast online evaluation is efficient computation of the β(t) term. In the setting of a reduced simulation this term reads as N X β(t) = ||ci ||G Sdi (t) (rC(t),di (t) ) , i=1
rC(t),di (t)
( di (t) + sign di (t) − rdi (t) C(t), rdi (t) ∈ / BC(t) (di (t)) = , rdi (t) , rdi (t) ∈ BC(t) (di (t))
rdi (t) = arg min+ r∈R0
φ(r) − φ(di (t)) . r − di (t)
So, in theory, N small minimization problems need to be solved at any time t ∈ [0, T ], whereas in practice Lemma 1 allows to reduce the amount of problems to be solved. For example, say di (t) < r0 for a t ∈ [0, T ], i ∈ {1 . . . N }. Then condition (17) yields r0 ≤ rdi (t) . So, if di (t) + C(t) < r0 there is no need to compute rdi (t) as rC(t),di (t) = di (t) + C(t). However, we propose using a penalized Newton iteration scheme for efficient computation of the minimizer itself, which is detailed in Proposition 3 in the appendix. As the reduced state variable xr (t) changes continuously, so does rs , which keeps the number of Newton iterations for subsequent (next time-step) computations of rs small when using the old values as starting point for the new iterations.
10
D. Wirtz and B. Haasdonk
4 Numerical experiments In this part we present numerical experiments for synthetic dynamical systems using kernel expansions as nonlinearity. Before we state the system settings, we introduce two more error estimator modifications. First, the convergence result from Theorem 3 motivates a heuristic variant of our local estimators. This originates from the fact that the limit function ∆∞ LSLE reproduces itself when used as a-priori bound within the iterations. Numerically (up to the integration error) this can be exploited using the error estimate from the previous time step as bound at the current time step. We will refer to this time-discrete method by “LSLE TD”. Experiments indicate that this variant indeed seems to bound ||f (x(t))−f (xr (t))||G the iterated LSLE estimators from below. Second, when setting β(t) = we obtain ||x(t)−xr (t)||G the smallest possible estimation, since this is the best local Lipschitz constant at any t ∈ [0, T ]. As this version requires x(t) it is not considered a practical error estimator but rather a theoretical “Lower Bound”. Our test setting is given as follows: Let d = 240000 to represent a large-scale system, T = 20, G = Id . The kernel expansion (1) uses N = 20 and (xi )j := N50 −1 (i − 1), i = 1 . . . N, j = 1 . . . d. 2 2 The kernel used is a Gaussian Φ(x, y) = exp(− ||x − y|| /β ) with β = 224, which is chosen to have Φ(xi , xj ) < 10−5 ∀ |i − j| ≥ 2. This way a certain locality of the kernel expansion is ensured. Finally, the expansion coefficient vectors are given via (ci )j = exp(−(xi )j /15), i = 1 . . . N, j = 1 . . . d. Further we set x0 = (1 . . . 1)t ∈ Rd and B = (1 . . . 1)t ∈ Rd×1 . To emphasize the input-dependent behavior of 2 1 our estimators, we choose two system inputs u1 (t) = 25 sin( 3t ) and u2 (t) = 21 e−(12−t) which represent one oscillating and one localized input. We use an explicit Euler scheme with time-step ∆t = 0.05 as solver. In order to investigate the error estimator behavior, we decided to have means of controlling the projection subspace quality. Since the system setup is homogeneous in each √ component, we would ˆ = (1, . . . , 1)t / d as projection matrices. obtain zero approximation error for any input using Vˆ = W ˆ Then, for a given degree θ we set the rotated subspace projection matrices to V := RVˆ , W := RW ˆ for an orthogonal block matrix R := diag(R, Id−200 ) with cos(θ) − sin(θ) 200×200 ˆ . R := diag(R2 , . . . , R2 ) ∈ R , R2 = sin(θ) cos(θ) This way, θ continuously controls the quality of the projection subspace. The results for the different error estimators can be seen in Figure 2. The rigorosity of all estimators is verified as all are upper
True error GLE LSLE LSLE, 1 It LSLE, 2 It LSLE, 5 It LSLE TD Lower bound
−3
−4
10
0
2
4
6
8
10 Time
12
14
16
18
True error GLE LSLE LSLE, 1 It LSLE, 2 It LSLE, 5 It LSLE TD Lower bound
−3
10
Relative error estimates
Error estimates
10
20
−4
10
0
2
4
6
8
10 Time
12
14
16
18
20
Fig. 2 Absolute (left) and relative (right) error estimates using θ = .0005 and input u1
bounds of the true error. Iterations of the LSLE estimator improve the error bound. The fact that the LSLE TD variant behaves like a tight lower bound for their corresponding iterative schemes (LSLE, 5 It and LSLE TD are indistinguishable) supports the applicability of our heuristic LSLE TD approach.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
11
Furthermore, the LSLE TD variant grows with the same rate as the lower bound estimation up to about t = 16, which shows the effectiveness of the local Lipschitz constants together with the a-priori bounds. Compared to the GLE, the LSLEs have a significantly lower exponential increase rate. One can see that the first and second LSLE iterations fall back to their specific increase rate after a certain time at which the a-priori bound is too big to have a positive effect regarding the choice (20). Figure 3 displays the computation times versus the error estimations (left) and the output of the full and reduced simulations along with the output error bounds computed by the LSLE TD estimator (right) for θ = 0.05 and u2 . Note that identical conclusions can be drawn by using either input for Figures 2 and 3. The full error computation takes about 20s, and the GLE estimator is cheap but too conservative to have a significant meaning. Comparing identical symbols reveals that indeed iterating the error estimation process increases the computational time, but improves the accuracy of the estimators. Hence, the iteration number can be used as a balancing parameter between online run-time and error bound accuracy. The output of the system using input u2 is nicely bounded by the LSLE TD estimator.
True error GLE LSLE LSLE, 1 It LSLE, 2 It LSLE, 5 It LSLE TD
20
18
16
1.6
1.4
14
12
y(t) / yr(t)
Comp. time [s]
Full system Reduced system Lower bound Upper bound
1.5
10
1.3
1.2 8
1.1 6
1 4
0.9
2 0
10
10
10
20
10
30
10
40
∆(T)
10
50
10
60
10
70
10
0
2
4
6
8
10 t
12
14
16
18
20
Fig. 3 Computation times and output for u2 and θ = .05. The output error is bounded using the LSLE TD variant
Note that the reduced and full models’ outputs are indistinguishable in this plot. Table 1 shows the estimated errors at T = 20 for different θ values and both inputs. Albeit the true ||e(T )||G
GLE
LSLE
LSLE 2 It
LSLE 5 It
LSLE TD
Low. b.
θ = 0.0005, u1
1.3E−5
1.1E+73
3.2E+11
1.4E+2
6.6E−5
6.6E−5
2.2E−5
θ = 0.005, u1
1.3E−4
1.1E+74
3.2E+12
1.2E+5
6.7E−4
6.6E−4
2.2E−4
θ = 0.05, u1
1.3E−3
1.1E+75
3.2E+13
1.2E+8
2.3E+0
6.7E−3
2.2E−3
θ = 0.5, u1
1.3E−2
1.1E+76
3.1E+14
1.0E+11
9.2E+5
8.0E−2
2.2E−2
θ = 0.0005, u2
2.3E−5
1.1E+73
8.5E+11
2.3E+3
1.8E−4
1.8E−4
3.0E−5
θ = 0.005, u2
2.3E−4
1.1E+74
8.5E+12
1.0E+6
2.7E−3
1.8E−3
3.0E−4
θ = 0.05, u2
2.3E−3
1.1E+75
8.5E+13
5.6E+8
3.6E+2
1.9E−2
3.0E−3
θ = 0.5, u2
2.3E−2
1.1E+76
8.4E+14
3.4E+11
3.5E+7
3.6E+0
3.1E−2
Model
Table 1 Errors of estimation runs for different θ values and inputs u1 , u2 at T = 20
error is of the same magnitude for both inputs and any θ, the estimators perform differently for u1 and u2 . This table emphasizes the superiority of the LSLE compared to the GLE, as the non-iterated LSLE variant is about 62 orders of magnitude better than the GLE.
5 Conclusion and Outlook In this article we investigated the topic of model reduction for kernel-based nonlinear dynamical systems and presented a series of improving efficient a-posteriori error estimators. Special classes of kernels allow for the creation of a small reduced system (6)-(7) via loss-less low-dimensional argument evaluations.
12
D. Wirtz and B. Haasdonk
The error estimators are obtained simultaneously with the reduced simulation by solving an additional ODE. Due to a complete decomposition of the error estimators in an offline/online fashion, ultimately a rapid simulation including rigorous error bounds is ensured. Numerical experiments show the superiority of the local Lipschitz estimator using minimum secant gradients obtained via solving a small Newton problem. Furthermore, iterating the error estimation process allows to balance the computation costs against estimation sharpness. Future work will target at the improvement of the error estimators regarding the coefficient vector norms ||ci ||G in (24) as locality can also be used there to sharpen the error estimations. Finally, incorporation of techniques to obtain kernel expansion approximations (1) of arbitrary nonlinear dynamical systems and extension to parameterized systems will be in the focus of our research.
Acknowledgements The authors would like to thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart and the Baden-Wrttemberg Stiftung gGmbH.
References 1. P. Benner, Numerical linear algebra for model reduction in control and simulation, GAMM-Mitt. 29 (2) (2006) 275–296. 2. H. Sandberg, A. Rantzer, Balanced truncation of linear time-varying systems, IEEE T. Automat. Contr. 49 (2) (2004) 217 – 229. 3. B. Haasdonk, M. Ohlberger, Efficient reduced models and a-posteriori error estimation for parametrized dynamical systems by offline/online decomposition, Math. Comp. Model. Dyn. 17 (2) (2011) 145–161. 4. A. Antoulas, An overview of approximation methods for large-scale dynamical systems, Annu. Rev. Contr. 29 (2005) 181–190. 5. A. Antoulas, Approximation of Large–Scale Dynamical Systems, SIAM Publications, 2005. 6. K. Zhou, J. C. Doyle, K. Glover, Robust and Optimal Control, Prentice Hall PTR, 1996. 7. M. Rewienski, A trajectory piecewise-linear approach to model order reduction of nonlinear dynamical systems, Ph.D. thesis, MIT, Cambridge, USA (2003). 8. B. Bond, L. Daniel, Parameterized model order reduction of nonlinear dynamical systems, in: Proc. of ICCAD 2005, 2005, pp. 487–494. 9. N. Dong, J. Roychowdhury, Piecewise polynomial nonlinear model reduction, in: Proc. of DAC 2003, 2003, pp. 484–489. 10. B. Moore, Principal component analysis in linear systems: Controllability, observability, and model reduction, IEEE T. Automat. Contr. 26 (1) (1981) 17–32. 11. J. Scherpen, Balancing for nonlinear systems, Syst. Contr. Lett. 21 (2) (1993) 143–153. 12. S. Lall, J. Marsden, S. Glavaski, A subspace approach to balanced truncation for model reduction of nonlinear control systems, Int. J. Robust Nonlin. 12 (5) (2002) 519–535. 13. M. Condon, R. Ivanov, Empirical balanced truncation of nonlinear systems, J. Nonlinear Sci. 14 (2004) 405–414. 14. P. Tabuada, A. D. Ames, A. Julius, G. J. Pappas, Approximate reduction of dynamic systems, Syst. Contr. Lett. 57 (7) (2008) 538–545. 15. J. Phillips, Projection frameworks for model reduction of weakly nonlinear systems., in: Proc. of DAC 2000, 2000, pp. 184–189. 16. P. Benner, T. Breiten, Krylov-subspace based model reduction of nonlinear circuit models using bilinear and quadratic-linear approximations, Max-Planck-Institute for Dynamics of Complex Technical Systems, Magdeburg, Germany. Submitted (2010). 17. S. Chaturantabut, D. Sorensen, Discrete empirical interpolation for nonlinear model reduction, in: Proc. of CDC/CCC 2009, 2009, pp. 4316–4321. 18. K. Kunisch, V. Volkwein, Proper Orthogonal Decomposition for Optimality Systems, ESAIM, Math. Model. Numer. Anal. 42 (2008) 1–23.
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
13
19. I. T. Jolliffe, Principal Component Analysis, Springer-Verlag, 2002. 20. J. Phillips, J. Afonso, A. Oliveira, L. Silveira, Analog macromodeling using kernel methods, in: Proc. of ICCAD 2003, 2003, pp. 446–453. 21. B. Schlkopf, A. J. Smola, Learning with Kernels, The MIT Press, 2002. 22. H. Wendland, Scattered data approximation, Cambridge University Press, 2005. 23. L. Ma, Z. Wu, Kernel based approximation in sobolev spaces with radial basis functions, Appl. Math. Comput. 215 (6) (2009) 2229–2237. 24. H. Wendland, Computational aspects of radial basis function approximation, in: Topics in Multivariate Approximation and Interpolation, Elsevier, 2006, pp. 231–256. 25. N. Aronszajn, Theory of reproducing kernels, Trans. Am. Math. Soc. 68 (3) (1950) 337–404. 26. C. Berg, J. Christensen, P. Ressel, Harmonic Analysis on Semigroups, Springer, 1984. 27. M. Hinze, S. Volkwein, Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: Error estimates and suboptimal control in dimension reduction of large-scale systems, in: Dimension Reduction of Large-Scale Systems, Lecture Notes in Computational and Applied Mathematics, Vol. 45, Springer Berlin Heidelberg, 2005, pp. 261–306. 28. S. Volkwein, Model Reduction using Proper Orthogonal Decomposition, Lecture notes, University of Graz (2008). 29. T. Reis, M. Heinkenschloss, Model reduction with a-priori error bounds for a class of nonlinear electrical circuits, in: Proc. of CDC/CCC 2009, 2009, pp. 5376–5383. 30. M. Meyer, H. Matthies, Efficient model reduction in non-linear dynamics using the Karhunen-Lo`eve expansion and dual-weighted-residual methods, Comput. Mech. 31 (2003) 179–191. 31. M. Drohmann, B. Haasdonk, M. Ohlberger, Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation, Preprint Angewandte Mathematik und Informatik 02/10-N, University of Mnster, Submitted to SISC (2010). 32. N. Jung, A. Patera, B. Haasdonk, B. Lohmann, Model order reduction and error estimation with an application to the parameter-dependent eddy current equation, Math. Comp. Model. Dyn., Accepted. 33. J. Hale, Ordinary differential equations, Wiley-Interscience, 1969.
Appendix A Detailed proofs A.1 Proof of Lemma 1 Proof (Proof ) It is clear that φ0 (r0 ) ≤ Ss (r) ≤ 0 ∀ r ∈ R+ 0 . Using the boundedness of φ we see that φ(r) − φ(s) ≤ B + φ(s) r→∞ |Ss (r)| ≤ −→ 0. r−s |r − s| As Ss is continuous and Ss (r0 ) < 0 ∀ s we see that rs < ∞ and hence obtain well-definedness of (16). 0 s (rs ) , so Next we show the position inequalities. Necessary condition for a minimizer is 0 = Ss0 (rs ) = φ (rsr)−S s −s 0 φ (rs ) = Ss (rs ). At first consider the case s < r0 . We continue with proof by contradiction and therefore assume φ(r)−φ(s) we had rs < r0 . Now choose a r ∈ R+ = φ0 (ξ) for 0 with max(rs , s) < r < r0 . Then we have Ss (r) = r−s some ξ ∈]s, r[ by the mean value theorem. Now as φ is strictly concave on [0, r0 ] by the bell function condition (14) and φ0 (r) ≤ 0 by (13) we obtain Ss (rs ) = φ0 (rs ) > φ0 (ξ) = Ss (r), which contradicts with the minimizing property of rs as obviously rs is not a minimizer of Ss . The case r0 < s is obtained by an analogous argument using the strict convexity of φ on [r0 , ∞[. Finally, the case r0 = s is clear when considering the limit process for both former cases. Now, uniqueness of the minimizers follows directly using the necessary condition for minimizers and the concavity/convexity ”on the other side“ as e.g. for two assumed minimizers rs 6= rs0 we always obtain a contradiction Ss (rs ) = φ0 (rs ) 6= φ0 (rs0 ) = Ss (rs0 ) to the minimal choice. Having s = r0 directly gives uniqueness by (19).
A.2 Auxiliary lemma for the proof of Proposition 1. + 0 Lemma 2 Let s ∈ R+ 0 \{r0 } and Ss be given as in (15). Then Ss (r)(r − rs ) > 0 ∀ r ∈ R0 \{rs }.
14
D. Wirtz and B. Haasdonk 0
φ (r)−Ss (r) 0 Proof (Proof ) Choose s ∈ R+ . First consider the case s < r0 . Then we 0 and note that Ss (r) = r−s differentiate locations of r: – r < s < r0 ≤ rs : We have r − rs < 0, r − s < 0 and Ss (r) > φ0 (r) by concavity. 0 0 φ00 (t)−Ss (t) s (t) – r = s < r0 ≤ rs : We have Ss0 (s) = lim Ss0 (t) = lim φ (t)−S = lim = φ00 (s) − lim Ss0 (t). This t−s 1 t→s
00
t→s
t→s
t→s
means Ss0 (s) = φ 2(s) < 0 and as r < rs we have Ss0 (r)(r − rs ) > 0. – s < r ≤ r0 ≤ rs : We have r − rs < 0, r − s > 0 and Ss (r) < φ0 (r) by concavity. – s < r0 ≤ r < rs : Since φ0 (r0 ) < Ss (r0 ) and φ0 (rs ) = Ss (rs ) we must have φ0 (r) < Ss (r) ∀ r ∈ [r0 , rs [. Otherwise, since φ0 (·) − Ss (·) is continuous, there would exist an r0 ∈ [r0 , rs [ with Ss (r0 ) = φ0 (r0 ) < φ0 (rs ) = Ss (rs ) as r0 < rs and φ is strictly convex on [r0 , ∞[. But this is a contradiction to the minimal choice of rs , and together with r − rs < 0, r − s > 0 we obtain Ss0 (r)(r − rs ) > 0. 0 00 φ00 (rs )−2Ss (rs ) s) – s < r0 ≤ rs < r: At rs we have Ss0 (rs ) = 0 and thus Ss00 (rs ) = = φrs(r > 0, as φ is rs −s −s 0 strictly convex for r > r0 . So there exists an > 0 with Ss (rs + ) > 0. In fact for any point r > r0 with Ss0 (r) = 0 we have Ss00 (r) > 0. But as Ss0 is continuous there cannot be two consecutive local minima without a local maximum. Hence, Ss0 (r) cannot equal zero for r > rs and since Ss0 (rs + ) > 0 we must have Ss0 (r) > 0 ∀ r > rs . The case r0 < s is obtained using an analogous argumentation.
A.3 Proof of Theorem 3 Proof (Proof ) At first we see that Γn is open (in fact a union of open intervals) since di (t), Cn (t) and rdi (t) are continuous in t. PN For n ∈ N let βn (t) := i=1 ||ci ||G |Sdi (t) (rCn (t),di (t) )| and si := di (t) as shorthand for any fixed t ∈ [0, T ]. We prove (26) by induction and n begin with n = 0. We haveoΓ1 6= ∅ by assumption. So, for t ∈ Γ1 condition (25) implies that M := i ∈ {1 . . . N } rsi ∈ / BC1 (t) (si ) is nonempty. Further, we know that |Ssi si + sign (rsi − si ) C1 (t) | < |Ssi (rsi )| ∀ i ∈ M as Lemma 2 holds since si 6= r0 ∀ i ∈ M . Then we deduce X X β1 (t) = ||ci ||G Ssi (rC1 (t),si ) + ||ci ||G Ssi (rC1 (t),si ) i∈M (20)
=
X
i∈M /
X ||ci ||G Ssi si + sign (rsi − si ) C1 (t) + ||ci ||G Ssi (rsi ) i∈M /
i∈M
inf Γ1 we know with (28) that Z
Rs
Z
β1 (r)dr
α(s)e0
ds
inf Γn , i.e. Cn+1 (t) < Cn (t) ∀ t > inf Γn by definition. This means Γn+1 6= ∅ as Γn ⊆ Γn+1 . So, for any t ∈ Γn+1 we have rCn+1 (t),si = si − Cn+1 (t) > si − Cn (t) = rCn (t),si which gives βn+1 (t) < βn (t) analogous to (28). Splitting the integral as above using Γn+1 directly gives (26), concluding the induction. n If Γ1 = ∅ we have βn ≡ β1 ∀ n ≥ 1 and thus ∆C LSLE ≡ C1 ∀ n ≥ 1. n Next, for each t ∈ [0, T ] we see from (26) that ∆C LSLE (t) n∈N is monotonically decreasing and bounded by Cn ∞ n zero. Thus, we have a pointwise convergence of ∆C LSLE to the limit function ∆LSLE (t) := limn→∞ ∆LSLE (t). From for ∀ n, m ∈ N, m ≥ n we have βm (t) ≤ βn (t) ∀ t ∈ [0, T ]. This the inductionC we know that Cn m n Cm implies ∆C LSLE (t) − ∆LSLE (t) ≤ ∆LSLE (T ) − ∆LSLE (T ) ∀ t ∈ [0, T ]. Taking the supremum on the left and considering m → ∞ we obtain n→∞ ∞ Cn Cn ≤ ∆∞ ∆LSLE − ∆LSLE LSLE (T ) − ∆LSLE (T ) −→ 0, L∞ (0,T )
which finally shows uniform convergence and continuity of ∆∞ LSLE (T ).
Efficient a-posteriori error estimation for nonlinear kernel-based reduced systems
15
B Penalized Newton iteration We first state some preparatory lemmata before we propose our penalized Newton iteration in Proposition 3. Lemma 3 Let φ be a bell function. Then the mapping + η(s) : R+ 0 −→ R0 , s 7−→ rs = arg min Ss (r) is monotonically decreasing. r∈R+ 0
Proof (Proof ) Note here that η is well-defined by Lemma 1. Choose some s1 , s2 ∈ R+ 0 with s1 6= s2 . Without loss of generality assume s1 < s2 . Let ri := η(si ), i ∈ {1, 2}. At first consider the case s1 < s2 ≤ r0 . Define l(r) := φ(s1 ) + φ0 (r1 )(r − s1 ). Then we see that l(s2 ) ≤ φ(s2 ), as φ(s2 ) − φ(s1 ) φ(r1 ) − φ(s1 ) l(s2 ) − φ(s1 ) ≥ = φ0 (r1 ) = . s2 − s1 r1 − s1 s2 − s1 Using φ0 (r1 ) =
φ(r1 )−φ(s1 ) , r1 −s1
l can also be written as l(r) = φ(r1 ) + φ0 (r1 )(r − r1 ). So we conclude
φ(r1 ) − φ(s2 ) φ(r2 ) − φ(s2 ) ≤ r2 − s2 r1 − s2 φ(r1 ) − l(s2 ) φ(r1 ) − φ(s1 ) ≤ = = φ0 (r1 ), r1 − s2 r1 − s1 and since φ(x) is convex for x > r0 we must have r1 ≥ r2 . The proof for the case r0 ≤ s1 < s2 is analogous using the concavity on [0, r0 ]. φ0 (r2 ) =
Lemma 4 (Bounds of rs ) Let φ be a bell function. Then we have 0 ≤ rs < rm ∀ s ∈ R+ for rm :=
φ(0)r0 . φ(0) − φ(r0 )
Proof (Proof ) Fix an s ∈ R+ . We already know that 0 ≤ rs by definition. At first we see that rm = φ(0) r0 φ(0)−φ(r > r0 φ(0) = r0 . So, if s > r0 we have rs ≤ r0 < rm by Lemma 1, (18). Obviously, for s = r0 we φ(0) 0) have rs = r0 < rm . This leaves us with the case s < r0 , i.e. r0 ≤ rs by (17). Let ls (r) := φ(r0 ) + Ss (r0 )(r − r0 ) be the line through φ(r0 ) with gradient Ss (r0 ). Since Ss (r0 ) < 0 and ls (r0 ) = φ(r0 ) > 0 we know that 0 ∃! ηs ∈ R+ 0 : ls (ηs ) = 0. Also, as 0 > Ss (r0 ) > φ (r0 ) we have ∃ δ > 0 : φ(r0 + δ) < ls (r0 + δ). But as φ(ηs ) > 0, ls (ηs ) = 0 and φ − ls is continuous we know that ∃ t ∈ ]r0 , ηs [ : φ(t) = ls (t). Further we have Sr0 (t) =
φ(t) − φ(r0 ) ls (t) − φ(r0 ) φ(r0 ) + Ss (r0 )(t − r0 ) − φ(r0 ) = = = Ss (r0 ). t − r0 t − r0 t − r0
Since φ is convex on [r0 , ∞[, ∃! ξ ∈]r0 , t[: φ0 (ξ) = Sr0 (t). With the minimizing property of rs we see that φ0 (ξ) = Sr0 (t) = Ss (r0 ) > Ss (rs ) = φ0 (rs ), but this means rs < ξ by the strict convexity of φ on [r0 , ∞[. Finally, as s ∈ R+ ∀ s ∈ R+ 0 was arbitrary, we obtain rs < ξ < t < ηs 0 . Now Lemma 3 yields rs ≤ rs=0 = η(0) < η0 ∀ s ∈ R+ , and solving l (η ) = 0 gives the desired bound rm = η0 . 0 0 0 φ(0)r0 . Define the φ(0)−φ(r0 ) n(ν) 2 n0 (ν) − ν. Then rs =
Proposition 3 (Penalized Newton iteration) Let φ be a bell function and rm := unpenalized objective function as n(r) := φ0 (r) − Ss (r) and c(ν) := arg min φ(r)−φ(s) r−s r∈R+ 0
n0 (ν)2 , d(ν) 4n(ν)
:=
can be computed as root of the penalized objective function np given by
s ≤ r0 :
r0 < s :
2 r ≤ r0 c(r0 )(r + d(r0 )) r0 < r < rm np (r) := n(r), c(rm )(r + d(rm ))2 , r ≥ rm 2 r≤0 c(0)(r + d(0)) , 0 < r < r0 np (r) := n(r), c(r0 )(r + d(r0 ))2 , r ≥ r0
Proof (Proof ) At first consider the penalty polynomials and define pν (x) = c(ν)(x + d(ν))2 . Then it is easy to validate that pν (ν) = n(ν), p0ν (ν) = n0 (ν). Those polynomials will be used to extend n at appropriate points. Next, existence and uniqueness of rs have already been shown in Lemma 1. At first consider the case s < r0 . 0 s (r) , and since s < rs by Lemma As the necessary condition for a minimizer r ≥ r0 we have 0 = Ss0 (r) = φ (r)−S r−s 0 0 1, 0 = n(r) = φ (r) − Ss (r) is also sufficient for r ≥ r0 . Since n(s) = φ (s) − Ss (s) = 0, we replace n(r) for r < r0 by pr0 (r). Even though Ss0 (r) > 0 ∀ r > rs by Lemma 2 and hence n(r) 6= 0 ∀ r > rs we see that limr→∞ n(r) = 0. To avoid the Newton iteration being drawn to ∞ we replace n(r) by prm for r ≥ rm as we know that rs < rm ∀ s ∈ R+ 0 from Lemma 4. For s > r0 we have 0 ≤ rs and so we replace n(r) by p0 for r ≤ 0 and n(r) by pr0 for r0 ≤ r since again n(s) = 0. Of course, s = r0 means rs = r0 and no iterations have to be performed.
