5.2 Error estimation for DEIM reduced systems
5.2.4 Efficient error estimation
The error estimator introduced in Theorem 5.2.9 makes explicit use of the local logarithmic Lipschitz constant LG[f](xr(t)), which is (as well as its global counterpart LG[f]) not readily available in most practical situations.
Therefore, letJ : Ω →Rd×d denote the Jacobian J(x)off atx ∈ Ω. With
the Taylor expansion off around xr(t) we obtain he(t),f(x(t))−f(xr(t))iG
||e(t)||2G = D
e(t),J(xr(t))e(t) +O
||e(t)||2GE
G
||e(t)||2G
= he(t),J(xr(t))e(t)iG
||e(t)||2G +O(||e(t)||G), (5.33) for anyt ∈ [0, T]. With Definition 5.2.5 and (5.25, p.178) this directly gives a first order approximation of the local logarithmic Lipschitz constant
LG[f](xr(t)) = LG[J(xr(t))] + O(||e(t)||G). (5.34) Using the approximation (5.34) avoids the need to obtain LG[f](xr(t)), but it comes with additional cost: The computation of the Jacobian logarithmic norm is expensive as it involves solving an eigenvalue problem of high di-mensionddue to the representation (5.26, p.178).
We will address this issue in the following discussion, where we propose to apply a suitable partial similarity transformation to the Jacobians which has been designed to preserve the largest eigenvalues of their symmetric parts.
The key ingredient is to perform a POD of their corresponding eigenvectors, which in turn allows us to bound the resulting eigenvalue approximation error in terms of the remaining eigenvalues of their covariance matrix. We explain this idea for general symmetric matrices in the following theorem, and refer to [97, 176] for details on POD and related error estimates. Recall here the MatLab style notation introduced in Chapter 1.
Theorem 5.2.11 (Approximation of eigenvalues for a family of symmet-ric matsymmet-rices). Let a continuous family of symmetric matrices H(t) ∈ Rd×d overt ∈ [a, b]be given and let [λ(t),q(t)] := λmax(H(t)) denote the largest eigenvalueλ(t) with corresponding normalized eigenvectorq(t) of H(t). Let
5.2. ERROR ESTIMATION FOR DEIM REDUCED SYSTEMS 187 supt∈[a,b]||H(t)|| ≤ CH hold. Further, define
R= Z b
a
q(t)q(t)Tdt,
and let QΣ2QT = Rbe the eigen-decomposition of Rwith
QTQ = I, Σ = diag (σ1, σ2, . . . , σd), σ1 ≥ σ2 ≥ · · · ≥σd > 0.
For k ≤ d letQk := Q[:,1:k]andλk(t) := λmax(QTkH(t)Qk). Then
b
Z
a
|λ(t)−λk(t)|dt ≤4CH
X
j>k
σ2j.
Proof. First, define the vector valued function w(t) := Σ−1QTq(t), t ∈ [a, b]. Note that
Z b a
w(t)w(t)Tdt = Z b
a
Σ−1QTq(t)qT(t)QΣ−1dt= Σ−1QTRQΣ−1 = I.
Thus, Rb
a wi(t)wj(t)dt = δij (the Kronecker delta), where wi(t) is the i-th component of w(t). Now, partition
Q = [Qk,Q˜k], Σ = diag
Σk,Σ˜k
, w(t) = [wTk(t),w˜Tk(t)]T, withQ˜k := Q(:,(k+1):d),Σk,Σ˜kdenoting the appropriate diagonal blocks ofΣandwk(t),w˜k(t)denoting the corresponding sub-vectors ofw(t). Put q(t) = qk(t) + ˜qk(t), with qk(t) := QkΣkwk(t), q˜k(t) := ˜QkΣ˜kw˜k(t), and note that qTk(t)˜qk(t) = 0for all t ∈ [a, b]. We observe that
Z b a
˜
qTk(t)˜qk(t)dt= Z b
a
˜
wTk(t) ˜Σ2kw˜k(t)dt = X
j>k
σ2j Z b
a
wj2(t)dt = X
j>k
σj2. (5.35)
In the following we omit the argument t and set q = q(t), qk = qk(t),
˜
qk = ˜qk(t), etc. Then
qTkHqk = (q−q˜k)TH(q−q˜k) =λ−2qTHq˜k + ˜qTkHq˜k
= λ−2λqTq˜k + ˜qTkHq˜k = λ−2λ˜qTkq˜k + ˜qTkHq˜k. (5.36)
Moreover, sinceµ ≤ ||H|| ≤ CH for any eigenvalueµofH, the definitions ofλ = λ(t) and λk = λk(t) imply
λ = sup
v6=0
vTHv
vTv ≥ sup
Qkvk6=0
vTkQTkHQkvk
vTkvk = λk ≥ qTkHqk
qTkqk . (5.37) Combining (5.36) and (5.37) provides with||qk||2 +||˜qk||2 = ||q||2 = 1 that
0≤ λ−λk ≤ λ− qTkHqk qTkqk
= λ− λ−2λ˜qTk˜qk + ˜qTkHq˜k
1−q˜Tkq˜k = λ˜qTkq˜k −q˜TkHq˜k 1− ||˜qk||2
=
λ− q˜TkHq˜k q˜Tkq˜k
||˜qk||2
1− ||˜qk||2 ≤2||H|| ||˜qk||2 1− ||˜qk||2, which is equivalent to
0 ≤ (λ−λk)(1− ||q˜k||2) ≤2||H|| ||q˜k||2. Hence,
0 ≤λ −λk ≤ (2||H|| +λ−λk)||˜qk||2 ≤ 4||H|| ||˜qk||2 ≤4CH ||˜qk||2, and therefore, using (5.35),
b
Z
a
|λ(t)−λk(t)|dt≤ 4CH
b
Z
a
||˜qk(t)||2dt = 4CH X
j>k
σj2
5.2. ERROR ESTIMATION FOR DEIM REDUCED SYSTEMS 189 as claimed.
Remark5.2.12. IfRis rank deficient, the above argument is still valid simply by replacingΣ−1 with the pseudo-inverse Σ+.
This result can be applied directly in the context of approximating the log-arithmic norms of the local Jacobians.
Proposition 5.2.13( Jacobian partial similarity transform). LetCCT denote the Cholesky decomposition of the weighting matrixG, and consider the family of symmetric matrices
H(t) := 1 2
CTJ(xr(t))C−T + CTJ(xr(t))C−TT .
Then, we have the corresponding values CH > 0, {σi}i=1...d and Q ∈ Rd×d from Theorem 5.2.11, that allow us to write
λ(t) = LG[J(xr(t))], λk(t) = LIk
QTkCTJ(xr(t))C−TQk
. (5.38)
Now we directly obtain the estimate
T
Z
0
LG[J(xr(t))]−LIk
QTkCTJ(xr(t))C−TQk
dt ≤CHX
j>k
σj2. (5.39)
Note here that in practice the matrix R in Theorem 5.2.11 is not available.
Instead,Qis obtained as the set of left singular vectors of the singular value decomposition (SVD) of a snapshot matrix
Sn =
rb−a n
h
q(t1) . . . q(tn)i .
Assuming equally spaced pointstj ∈ [a, b]with t1 = a, tn = b, we have
n→∞lim SnSTn =
b
Z
a
q(t)q(t)Tdt = R.
Thus, for a sufficiently large n, the sum of the neglected squared singular values ofSn will be arbitrarily close to the corresponding neglected eigen-values of R and can safely be used in the estimate. The detailed argument is similar to one given in [97]. However, for a given set of training data X = {xi | i = 1. . . n} ⊆ Ω, Algorithm 10 describes how to obtain Q in practice. Notationally, the method SV D performs the singular value de-compositionA= VΣWT for a matrix A.
Algorithm 10: Q = getTransformationMatrix(X)
1: Sn ←[]
2: for allxi ∈ X do
3: M ← CTJ(xi)C−T
4: [λi,qi] ← λmax 12(M +MT)
5: Sn ← [Sn qi]
6: end for
7: [Q,Σ,WT]← SV D q
b−a n Sn
8: returnQ
Before we can derive an efficient error estimator variant with the transfor-mation introduced above, there is one problem left to deal with. The reduced matrix from the right hand side of (5.38) is of small sizek ×k, however, its computation involves the transformed Jacobian CTJ(xr(t))C−T ∈ Rd×d which makes its computation infeasible during reduced simulations. Thus, we propose to apply a Matrix-DEIM approximation, which not only reduces evaluation costs for the Jacobian but also allows an efficient offline/online decomposition of (5.38), which we will discuss in Section 5.3. A similar idea named “Multi-Component EIM” has been formulated and successfully ap-plied in [171, §4.3.2]. Consequently, for anyA ∈ Rd×d we define the
trans-5.2. ERROR ESTIMATION FOR DEIM REDUCED SYSTEMS 191 formation
Υ : Rd×d →Rd
2
A 7→Υ[A] :=
AT1, AT2, . . . , ATd T
,
which maps the matrix entries of Acolumn-wise into a vector (equivalent to the MatLab operation A(:), also known asvec-operation).
Proposition 5.2.14 (Matrix DEIM). Choose MJ ≤ d and let UˆMJ,PˆMJ ∈ Rd
2 ×MJ denote the corresponding matrices for the DEIM basis and interpo-lation points obtained by application of the known DEIM approximation pro-cedure from Section 2.3/[23] to the vector-valued function Υ[J(x)]. Then, for mJ ≤ MJ, the mJ-th order MDEIM approximation ofJ is given via
J˜mJ(x) := Υ−1h
UˆmJ( ˆPmJUˆmJ)−1PˆTm
JΥ[J(x)]i , whereUˆmJ := ˆUMJ(:,1:mJ)and PˆmJ := ˆPMJ(:,1:mJ).
To this end, no extra assumptions on f need to be made that are not al-ready required by the standard DEIM, as the point-wise evaluations of the Jacobian entries can always be approximated via finite differences and point-wise evaluation of the underlying f. Of course, direct computation of Jaco-bian entries is preferred if possible. Moreover, the evaluation technique for reduced variables V z(t) carries over directly as proposed in the original work [23, §3.5]. We would also like to mention that any matrix approxima-tion technique using linear combinaapproxima-tions of basis matrices could be applied in this context, see [21] for an approach using a POD-basis with least squares weights and structure-preserving constraints.
Next we investigate the approximation quality for the logarithmic norm of the similarity transformed, MDEIM JacobianLIk
QTkCTJ˜mJ(yr(t))C−TQk against those of the true Jacobians LG[J(yr(t))]. Similar to Section 5.2.2, this is done outside the scope of the error estimation process, where we shall use the offline data from the 1D viscous Burgers model of Section 3.4.3
again. Consequently, Figure 5.10 shows the maximum (left) and mean (right)
Figure 5.10: Maximum (left) and mean (right) approximation error of the logarithmic norm overX
logarithmic norm approximation error overY for differentmJ andkvalues.
While values ofmJ = 13, k = 15are sufficient on average to have less than 1% relative error, in the worst casemJ = 28, k = 30 are needed to ensure the same tolerance over Y. But even for maximal mJ = k = 50, the rel-ative errors are already only 3.162 × 10−4 (worst case) and 8.222 × 10−6 (average). These results strongly indicate that both performing the simi-larity transformation and MDEIM approximation of Propositions 5.2.13 and 5.2.14, respectively, are suitable for a low-cost but accurate approximation of Jacobian logarithmic norms.
Now, using the above results from Propositions 5.2.13 and 5.2.14, we de-rive an efficient approximation of the estimator introduced in Theorem 5.2.9, where the computational complexity for obtaining β(t) (see (5.31c, p.184)) is independent of the high dimensiond.
Theorem 5.2.15(Estimator using local Jacobian logarithmic norms). Let the conditions from Theorem 5.2.8, Propositions 5.2.13 and 5.2.14 hold and replace (5.31c, p.184) by
β(t) := LIk
QTkCTJ˜mJ(xr(t))C−TQk
. (5.40)
Then the estimator (5.31a, p.184) can be approximated arbitrarily close with increasingmJ and k and is exactly reproduced formJ = d2, k = d.
5.2. ERROR ESTIMATION FOR DEIM REDUCED SYSTEMS 193 Proof. We have J ≡ J˜mJ for mJ = d2 as each entry of the Jacobian is interpolated. Together with (5.39, p.189) for k = d we obtain equality of (5.31c, p.184) and (5.40).
Remark 5.2.16. One other obvious alternative in order to obtain an approx-imation of the Jacobian logarithmic norm is to use the reduced projected Jacobian WTJ(xr(t)))V and thus only solve an eigenvalue problem of size r d. When also applying the MDEIM approximation of Proposi-tion 5.2.14, the only real difference lies in the transformaProposi-tion with V,W instead of Qk, with the advantage of not having to compute Q. However, as the projection subspace spanned by V is not designed to preserve any eigenvalues and the results from Proposition 5.2.13 are “optimal” in a sense, this approach is expected to have a lower approximation quality of the log-arithmic norms.
Addressing this problem by incorporating suitable information into V,W is certainly an interesting question for future work, but does not essentially avoid the necessity of a training sample set of eigenvectors during subspace computation.