https://doi.org/10.1007/s10444-020-09763-5
Approximation of stability radii for large-scale
dissipative Hamiltonian systems
Nicat Aliyev1·Volker Mehrmann2·Emre Mengi3
Received: 30 August 2018 / Accepted: 17 January 2020 /
©The Author(s) 2020
Abstract
A linear time-invariant dissipative Hamiltonian (DH) system ˙x=(J −R)Qx, with
a skew-Hermitian J, a Hermitian positive semidefinite R, and a Hermitian posi-
tive definite Q, is always Lyapunov stable and under further weak conditions even
asymptotically stable. By exploiting the characterizations from Mehl et al. (SIAM
J. Matrix Anal. Appl. 37(4), 1625–1654, 2016), we focus on the estimation of two
stability radii for large-scale DH systems, one with respect to non-Hermitian per-
turbations of Rin the form R+BCHfor given matrices B,C, and another with
respect to Hermitian perturbations in the form R+BBH, =H. We propose
subspace frameworks for both stability radii that converge at a superlinear rate in
theory. The one for the non-Hermitian stability radius benefits from the DH structure-
preserving model order reduction techniques, whereas for the Hermitian stability
radius we derive subspaces yielding a Hermite interpolation property between the full
and projected problems. With the proposed frameworks, we are able to estimate the
two stability radii accurately and efficiently for large-scale systems which include a
finite-element model of an industrial disk brake.
Keywords Dissipative hamiltonian system ·Robust stability ·Stability radius ·
Eigenvalue optimization ·Subspace projection ·Structure-preserving subspace
framework ·Hermite interpolation
Mathematics Subject Classification (2010) 65F15 ·93D09 ·93A15 ·90C26
Communicated by: Anthony Nouy
Supported by Deutsche Forschungsgemeinschaft via Project ME 40-1 of Priority Program 1897
Calm, Smooth and Smart.
Supported by Deutsche Forschungsgemeinschaft via Project A02 of Sonderforschungsbereich 910.
This article belongs to the Topical Collection: Model reduction of parametrized Systems
Guest Editors: Anthony Nouy, Peter Benner, Mario Ohlberger, Gianluigi Rozza, Karsten Urban and
Karen Willcox
Extended author information available on the last page of the article.
Published online: 6 February 2020
Adv Comput Math (2020) 46: 6
1 Introduction
Linear time-invariant Dissipative Hamiltonian (DH) systems are of the form
˙x=(J −R)Qx, (1.1)
where Q=QH∈Cn×nis a Hermitian positive definite matrix denoted by Q>
0, whereas J=−JH∈Cn×n, and the dissipation matrix R=RH∈Cn×nis
positive semidefinite, denoted by R≥0. DH systems arise as homogeneous parts of
port-Hamiltonian (PH) systems without inputs and outputs.
PH and DH systems play an essential role in many areas of science and engineer-
ing, see e.g. [18,26]. They are automatically Lyapunov stable, i.e., all eigenvalues of
A=(J −R)Q are in the closed left half of the complex plane, and any eigenvalue on
the imaginary axis is semisimple, see [21]. However, DH systems are not necessar-
ily asymptotically stable, since Amay have purely imaginary eigenvalues. Even if a
DH system is asymptotically stable, small perturbations that do not preserve the DH
structure may cause the system to become unstable so that the solution exhibits expo-
nential growth for some initial conditions. Asymptotic stability in the presence of
uncertainty can only be guaranteed when the stability radius of the system is reason-
ably large, i.e., the system is not only asymptotically stable but also reasonably away
from the set of systems with purely imaginary eigenvalues. Hence, an estimation of
the stability radius, which involves an optimization problem over all admissible per-
turbations, is an important ingredient of a proper stability analysis. One concrete real
application where the knowledge of the stability radius is invaluable is the avoidance
of brake squeals, which is presented next.
1.1 A motivating real world application
Disk brake squeal is a well-known problem in mechanical engineering. It occurs
due to self-excited vibrations caused by instability at the pad-rotor interface [1]. In
[12], FE models for disk brakes are derived in the form of second order differential
equations
M¨x+D() ˙x+K()x =f, (1.2)
with large and sparse symmetric coefficient matrices M,D(),andK() ∈Rn×n
depending on the rotational speed >0suchthatM>0, K() > 0andD() ≥0.
The function frepresents a forcing term or control; hence, it may be set as f=
0 for the stability analysis. The incorporation of gyroscopic effects, modeled by a
term G() ˙x, with G() =−G()H, and circulatory effects, modeled by a non-
symmetric term Nx, gives rise to a system
M¨x+(D() +G()) ˙x+(K() +N)x =0. (1.3)
Expressing (1.3) (with f=0) as a first order ODE followed by a straightforward
manipulation yields the system
˙
z=(J −R)Q
z, (1.4)
Adv Comput Math (2020) 46:6
Page 2 of 28
6
where
J=−G() −K()
K() 0,R=D() N
00
,Q=M0
0K() −1
,
z=˙x
x,
z=Q−1z.
(1.5)
In the absence of the circulatory effects, i.e., when N=0, this is a DH system and, as
a result, it is Lyapunov stable and typically it is also asymptotically stable. However,
small circulatory effects, i.e., perturbations by a non-symmetric Nof small norm,
may result in instability. Here, we emphasize that the perturbations due to the matrix
Nare restricted to the finite element nodes on the brake pad and the associated nodes
on the disk. Moreover, the number of such nodes is several orders of magnitude
smaller than the system size.
In modern brakes, damping devices (shims) that increase the damping of the brake
pad are used to avoid resonances between the disk and the brake pad. The uncertain-
ties and the additional damping due to shims mainly affect the matrices D() and
N.
Note that for this brake squeal example leading to the DH system in (1.4)and(1.5),
the matrix Q−1is sparse, whereas Qis usually not. This however does not cause
any difficulty in our approaches; as elaborated on later when we perform numerical
experiments, the ideas we introduce throughout the paper can be modified so that
they are based on Q−1rather than Q.
1.2 Scope
In this paper, we study the stability radii of large-scale DH systems of the form (1.1)
in the presence of perturbations. We typically assume that the DH systems at hand are
sparse; however, the approaches we propose may also be applicable to dense systems
as long as it is feasible to solve a few linear systems of size equal to the order of the
system.
Motivated by the disk brake example, we focus on perturbations in the coefficient
matrix Rof (1.1). In particular, we consider the following two stability radii for a DH
system of the form (1.1) originally introduced in [21], where the notation iRrefers
to the imaginary axis in the complex plane, (A) to the spectrum of a matrix A,and
A2to the spectral norm of A.
Definition 1 Consider a DH system of the form (1.1) and suppose that B∈Cn×m
and C∈Cp×n,n≥m, p, are given restriction matrices.
– The non-Hermitian restricted stability radius r(R;B,C) with respect to pertur-
bations of Runder the fixed restriction matrices B,Cis defined by
r(R;B,C) := inf{2|∈Cm×p,
((J −(R+BC))Q)∩iR=∅}.
Adv Comput Math (2020) 46:6 Page 3 of 28 6
– The Hermitian restricted stability radius with respect to Hermitian perturbations
of Runder the fixed restriction matrix Bis defined by
rHerm(R;B) := inf{2|=H∈Cm×m,
(J −(R+BBH))Q∩iR=∅}.(1.6)
The significance of the two stability radii defined above in the context of the disk
brake example is discussed next.
Example 1 In the brake-squeal example in Section 1.1, it is of interest to know
whether, for given , the norm of the non-symmetric matrix Nis tolerable to pre-
serve the asymptotic stability of the system in (1.4). The relevant stability radius for
a specified is given by
inf {N2|(A(N))∩iR=∅
},(1.7)
where
A(N) := −G() −K()
K() 0−D() 0
00
+0N
00
M0
0K() −1
.
The perturbations in Nonly act on the FE nodes associated with the brake pad.
Hence, the stability radius in (1.7) corresponds to the non-Hermitian restricted stabil-
ity radius r(R;B,C) with the restriction matrices B=Bp0Hand C=0Cp
with blocks Bp,C
passociated with these finite element nodes. On the other hand, the
Hermitian stability radius rHerm(R;B) is associated with restricted Hermitian pertur-
bations with B=Bp0H,C =BHwhen the Hermitian matrix D() is subject
to Hermitian perturbations.
Note that the Hermitian stability radius rHerm(R;B)is finite (see Theorem 3 below
and [21, Corollary 4.10]) and can be arbitrarily small. Additionally, by definition, we
have r(R;B,BH)≤rHerm(R;B). Indeed, these two stability radii can differ con-
siderably, see also [21]. Our main concern is the approximation of these two stability
radii for large-scale DH systems when nm, p. We propose subspace frameworks
that replace the large-scale DH system with a small-scale system whose stability radii
approximate those of the large-scale system usually well.
Contributions The stability radius under non-Hermitian restricted perturbations
r(R;B,C) can be characterized as the reciprocal of the H∞norm of the transfer
function of a linear control system, see in particular (2.1) in the next section. It should
be noted that, due to the DH structure, the singular value function that needs to be
optimized in the H∞-norm computation often has many local maximizers.
For the estimation of r(R;B,C), we propose a DH structure–preserving subspace
framework that converges at a superlinear rate with respect to the subspace dimen-
sion. The idea of the new method is partly inspired by the subspace framework in
[2] applicable for the approximation of the H∞norm of an asymptotically stable
system, but we form the subspaces in a different way by combining it with DH
structure–preserving model order reduction techniques, see [11,15,16,24]. Both the
Adv Comput Math (2020) 46:6
Page 4 of 28
6
framework in [2] and the one proposed here seem to converge to local maximizers
that are not necessarily global maximizers. However, our numerical experience on
synthetic examples indicates that the DH structure–preserving framework converges
to the global maximizer more frequently compared with the unstructured framework
in [2]. There are other approaches [5–7,10,17,23] for approximating the H∞norm
of a large-scale asymptotically stable system, but none of these employ subspace
projections.
The Hermitian restricted stability radius rHerm(R;B) poses more challenges.
Unlike the characterization for r(R;B,C), the eigenvalue optimization characteri-
zation for rHerm(R;B) (given in Theorem 3) does not appear to be related to any
transfer function. We derive a new subspace projection framework in which Her-
mite interpolation properties between the original and projected objective eigenvalue
functions hold.
Outline The paper is organized as follows. Section 2proposes subspace frameworks
for approximating r(R;B,C), while Section 3concerns rHerm(R;B). For the latter,
even the solution of the small-scale problems is challenging. Hence, we first discuss
how small-scale problems can be solved in Section 3.1. The new structured subspace
framework is then presented in Section 3.2. In Sections 2and 3, we illustrate the
new frameworks by means of numerical experiments on the disk brake example and
several synthetic examples.
A longer version of this article discussing also perturbations in coefficients other
than R, and featuring additional numerical results is available in [4].
2 The non-Hermitian stability radius for large-scale problems
In this section, we study the estimation of the non-Hermitian restricted stability radius
r(R;B,C) for large-scale asymptotically stable DH systems. Our approach operates
on the characterization [21, Theorem 3.3]
r(R;B,C) =1/GH∞,
GH∞:= sup
s∈C+
σmax(G(s)) =max
ω∈Rσmax(G(iω)). (2.1)
for an asymptotically stable DH system of the form (1.1) and restriction matrices
B∈Cn×m,C∈Cp×n. Above, we use the notation
G(s) := CQ(sI −(J −R)Q)−1B, (2.2)
the notation σmax(G(s)) for the largest singular value of G(s),andC+:= {s∈
C|Re(s) > 0}for the open right-half of the complex plane. Note that the second
equality in the second line in (2.1) follows from the maximum modulus principle, as
the matrix (J −R)Q does not have any eigenvalue in C+.
The matrix-valued function G(s) in (2.2) corresponds to the transfer function of
the control system
˙x=
Ax +
Bu, y =
Cx, (2.3)
Adv Comput Math (2020) 46:6 Page 5 of 28 6
where we set
A:= (J −R)Q,
B:= Band
C:= CQ. Moreover, the quantity
GH∞is commonly referred as the H∞norm of the transfer function G(s).
2.1 The subspace framework of [2]
Recently, in [2], a subspace framework for the approximation of the H∞norm of a
large-scale asymptotically stable system has been proposed. At iteration k, it requires
two subspaces Vkand Wkof equal dimension, as well as matrices Vkand Wk, respec-
tively, whose columns form orthonormal bases for these subspaces. The state xof
system (2.3) is restricted to Vk, i.e., in (2.3)thestatexis replaced by Vkxk, and it is
imposed that the residual of the differential part after this restriction is orthogonal to
Wk.
This gives rise to the reduced order system
Ek˙xk=
Akxk+
Bku, yk=
Ckxk,(2.4)
with
Ek:= WH
kVk,
Ak:= WH
k
AVk,
Bk:= WH
k
B,and
Ck:=
CVk. Then letting
Gk(s) :=
Ck(s
Ek−
Ak)−1
Bkbe the reduced order transfer function, and assuming
that
Ekis invertible, a maximizer ωk+1∈arg maxω∈Rσmax(Gk(iω)) can be com-
puted efficiently and reliably if the dimensions of Vk,Wkare small by employing the
globally convergent method in [8,9]. Once such an ωk+1has been determined, then
the subspaces Vkand Wkare expanded into larger subspaces Vk+1and Wk+1in such
a way that the corresponding reduced transfer function Gk+1(s) satisfies the Hermite
interpolation conditions
σmax(G(iωk+1)) =σmax(Gk+1(iωk+1)) ,
d
dω σmax(G(iω))ω=ωk+1
=d
dω σmax(Gk+1(iω))ω=ωk+1
. (2.5)
Denoting the range space of a matrix by Aby Ran(A),itisshownin[2] that the
inclusions
Ran((iωk+1I−
A)−1
B) ⊆Vk+1,Ran((
C(iωk+1I−
A)−1)H)⊆Wk+1,
ensure that the conditions in (2.5) are satisfied. For instance when m=p,
in practice, we set Vk+1:= orth Vk(iωk+1I−
A)−1
B,andWk+1:=
orth Wk(iωk+1I−
A)−H
CH, where orth(M) denotes a matrix whose columns
form an orthonormal basis for Ran(M). The definitions of the matrices Vk+1,W
k+1
can be modified in a simple way when m= pso that they have equal num-
ber of columns; we refer to [2] for the details. The procedure is then repeated
with the expanded subspaces Vk+1,Wk+1spanned by the columns of Vk+1,W
k+1.
The sequence {ωk}always seems to converge to a local maximizer, but a formal
proof is open at the moment. In [2], assuming that the sequence {ωk}converges
to a local maximizer, it is proven that this sequence converges at a superlinear
rate due to the satisfaction of (2.5)forallk. The solution of the reduced problem
maxω∈Rσmax(Gk(iω)) by the method in [8,9] requires extraction of all eigenval-
ues of 2d×2dpencils, where dis the order of the system, i.e., the dimension of
the subspaces Vk,Wk. Assuming that dis small, this is much cheaper than solving
Adv Comput Math (2020) 46:6
Page 6 of 28
6
the linear systems (iωk+1I−
A)−1
B,(iωk+1I−
A)−H
CHinvolving the full order
system, which is needed to form Vk+1,W
k+1.
Even though
A=(J −R)Q has the DH structure, this is in general not true for
the reduced system
Ak. In the next subsection, we modify the procedure of [2]to
preserve the DH structure.
2.2 A structure-preserving subspace framework
In this subsection, we introduce an interpolating and DH structure–preserving ver-
sion of the subspace projection framework. Structure-preserving subspace projection
methods in the context of model order reduction of large-scale PH and DH systems
have been proposed in [15,16,24]. Our approach is inspired by the following inter-
polation result, see [16, Remark 3]. For this result and the subsequent arguments, we
introduce the matrix-valued function
T(s):= (J −R)Q −sI (2.6)
depending on s∈Cassociated with the DH system in (2.3).
Theorem 1 (Right Tangential Interpolation [16]) Let G(s) be the transfer function
for the system (2.3), and let Gk(s) be the transfer function for the reduced system
(2.4). For given
s∈C,
b∈Cmand a positive integer N,if
T(
s)−
B
b∈Vkfor =1,...,N, (2.7)
and Wkis such that WH
kVk=I, then, denoting the -th derivatives of G(s) and
Gk(s) at the point
swith G()(
s) and G()
k(
s), we have
G()(
s)
b=G()
k(
s)
bfor =0,...,N −1 (2.8)
provided that both
sI −
Aand
sI −
Akare invertible.
The characterization of r(R;B,C) involves the maximization of the largest sin-
gular value of the transfer function G(s) on the imaginary axis. In the following, we
make use of Theorem 1 to obtain a reduced order system satisfying the interpolation
conditions (2.8) while retaining the structure.
For a given
s∈Cand
b∈Cm, choose Vkto satisfy (2.7). Then define
Wk:= QVk(V H
kQVk)−1,(2.9)
so that
WH
kVk=Iand (WkVH
k)2=WkVH
k,(2.10)
i.e., WkVH
kis an oblique projector onto Ran(QVk). The reduced system matrices
Ak,
Bk,
Ckthen satisfy
Ak=WH
k
AVk=WH
k(J −R)QVk=WH
k(J −R)WkVH
kQVk
=(Jk−Rk)Qk,(2.11)
Adv Comput Math (2020) 46:6 Page 7 of 28 6
where we define
Jk:= WH
kJW
k=−JH
k,R
k:= WH
kRWk=RH
k≥0,
Qk:= VH
kQVk=QH
k>0,(2.12)
and
Ck=
CVk=CQVk=CWkVH
kQVk=CkQk,C
k:= CWk,
Bk=WH
kB=: Bk. (2.13)
Moreover,
Ek=WH
kVk=I. Using this construction and Theorem 1, we deduce
the following result that indicates particular projection subspaces to achieve Hermite
interpolation between the full and reduced order systems in a structure-preserving
way.
Theorem 2 Consider a linear system as in (2.3)with the transfer function G(s) in
(2.2). Furthermore, for a given point
s∈Cand a given tangent direction
b∈Cm,
suppose that Vkis a matrix with orthonormal columns such that
T(
s)−B
b∈Ran(Vk)for =1,...,N
for some positive integer N,whereT(s) is as in (2.6). Define Wkas in (2.9),
Jk,R
k,Q
kas in (2.12), and Bk,C
kas in (2.13). Then the resulting reduced order
system
˙xk=(Jk−Rk)Qkxk+Bku,
yk=CkQkxk
is a DH system with the transfer function
Gk(s) := CkQk(sI −(Jk−Rk)Qk)−1Bk(2.14)
that satisfies
G()(ˆs)ˆ
b=G()
k(ˆs)ˆ
bfor =0...,N −1,
where G()(ˆs),G()
k(ˆs) denote the -th derivatives of G(s), Gk(s) at ˆs.
According to Theorem 2, for a given
s∈C, setting
Vk:= orth T(
s)−1BT(
s)−2B,W
k:= QVk(V H
kQVk)−1,
we obtain G(
s) =Gk(
s) and G(
s) =G
k(
s), and thus the Hermite interpolation
conditions
σmax(G(
s)) =σmax(Gk(
s)), d
dsσmax(G(s))s=
s
=d
dsσmax(Gk(s))s=
s
(2.15)
are satisfied. This suggests the greedy subspace framework outlined in Algorithm 1
for the approximation of r(R;B,C).
Adv Comput Math (2020) 46:6
Page 8 of 28
6
In line 6 of Algorithm 1, at every iteration, the subspace framework computes the
H∞norm of a reduced system, in particular it computes the point iω∗on the imagi-
nary axis where this H∞norm is attained. Then the current left and right subspaces
are expanded in a way so that the resulting reduced system still has the DH structure
and its transfer function Hermite interpolates the original transfer function at iω∗.
Since the Hermite interpolation conditions (2.15) are satisfied at
s=iω1,...,iωk
at the beginning of iteration k, the analysis in [2] shows a superlinear rate of con-
vergence for the sequence {ωk}, assuming that this sequence converges to a local
maximizer. We remark that σk+1=σmax(G(iωk+1)) due to line 10 of the algorithm
and the Hermite interpolation property. Hence, if the sequence {ωk}converges to a
global maximizer, the reciprocal of the limit of the sequence {σk}yields r(R;B,C).
The most expensive part of Algorithm 1 occurs in lines 2 and 7, where several
linear systems have to be solved. If this is done with a direct solver, for each value
ω∈Rone LU factorization of the matrix T(iω) suffices. For large n, the computa-
tion time is usually dominated by these LU factorizations. In contrast, the reduced
problem in line 6 involves a small system; its solution can be obtained efficiently and
robustly by the algorithm in [8,9].
2.3 Numerical experiments
In this subsection, we illustrate the performance of MATLAB implementations of
Algorithm 1 and the subspace framework from [2] via several numerical experiments
carried out in MATLAB R2017b on a machine with Mac OS 10.12.6, an Intel®
Core™i5-4260U CPU and 4GB RAM.
Algorithm 1 and the subspace framework from [2] are terminated when at least
one of the following three conditions is fulfilled:
Adv Comput Math (2020) 46:6 Page 9 of 28 6
1. The relative distance between ωkand ωk−1is less than a prescribed tolerance,
i.e., |ωk−ωk−1|<ε·|ωk|.
2. Letting fk:= maxω∈Rσmax(Gk(iω)), two consecutive iterates fk,f
k−1are close
enough in a relative sense, i.e., |fk−fk−1|<ε·fk.
3. The number of iterations exceeds a specified integer, i.e., k>k
max.
In all of the presented numerical examples, we use ε=10−12 and kmax =80 for
both Algorithm 1 and the subspace framework from [2].
Both Algorithms seem to converge locally in our experiments. The initial interpo-
lation points affect whether convergence to a global maximizer occurs. A procedure
to choose the initial interpolation points is outlined below. It uses an interval [l,u],
ideally containing a global minimizer. A reasonable strategy is to set l,u equal to
the approximations of the imaginary parts of the eigenvalues of (J −R)Q with the
smallest and largest imaginary parts.
(i) First, we discretize [l,u]into ρequally spaced points, say ω0,1,...,ω
0,ρ,
including the end-points l,u,whereρis specified by the user.
(ii) Then we either approximate the eigenvalues z1,...,z
ρof (J −R)Q closest
to iω0,1,...,iω0,ρ (using eigs in MATLAB), respectively, or otherwise set
zj=iω0,j for j=1,...,ρ.
(iii) Next, we permute z1,...,z
ρinto zj1,...,z
jρwhere {j1,...,j
ρ}={1,...,ρ}
is such that σmax(G(izj1)) ≥···≥σmax(G(izjρ)).
(iv) The interpolation points ω1,...,ω
κemployed initially are chosen as the
imaginary parts of zj1,...,z
jκ, where again κ≤ρis specified by the user.
For the synthetic examples below, we use l=−2000, u=0, ρ=25, κ=10, and,
in part (ii), z1,...,z
ρare set equal to the eigenvalues closest to iω1,...,iωρ.Forthe
disk brake example, we use l=−2·105,u=−1.7 ·105,ρ=15, κ=15, and, in
part (ii), we set zj=iω0,j for j=1,...,ρ. The values of land uemployed in these
examples are chosen based on rough sketches of the largest singular value functions
to be maximized.
2.3.1 Results on synthetic examples
We now present results for 2000 dense linear DH systems with random coefficient
matrices J,R,Q of order 500, and the restriction matrices Band Cchosen as 500×2
and 2 ×500 random matrices; MATLAB code generating such random coefficient
matrices are openly available.1
The relative errors of the values returned by Algorithm 1 and [2, Algorithm
1] for the first 5 random examples are presented in Table 1,wherefBB and
fSF denote the results returned by the Boyd-Balakrishnan (BB) algorithm [8,9]
and the subspace framework (Algorithm 1 or [2, Algorithm 1]). Note that for the
1http://home.ku.edu.tr/∼emengi/software/DH-stabradii-files/randomDH.m
This function needs to be called as [J,R,Q,B,C] = randomDH(500,2,2) to generate such random
matrices.
Adv Comput Math (2020) 46:6
Page 10 of 28
6
Table 1 Relative errors of the values returned by Algorithm 1 and [2, Algorithm 1] on 5 dense random
DH examples of order 500. The number of subspace iterations, run-time (in seconds) are also listed
|fSF −fBB|/fBB # iterations run-time
Ex. Alg. 1 [2]Alg.1[2]Alg.1[2]
12.13·10−14 1.85 ·10−14 868.86.8
22.98·10−15 5.60 ·10−2446.86.7
35.88·10−15 5.51 ·10−16 637.26.6
43.26·10−15 1.63 ·10−14 9 53 7.3 33.3
51.23·10−14 2.77 ·10−15 7157.27.7
second example, [2, Algorithm 1] returns an estimate with a significant error (high-
lighted in bold face). The unstructured distance to instability [27]((J −R)Q) :=
inf 2|∃∈Cn×n, ((J −R)Q +) ∩iR=∅
for all of these 5 examples
is larger than r(R;B,C) by ten times or even more.
Considering all 2000 examples, the non-Hermitian stability radii are estimated
with relative errors less than the prescribed tolerance 10−12 for 1921 of these exam-
ples by Algorithm 1 and for 1821 by [2, Algorithm 1]. These numbers correspond
to 96.45% and 91.45% accuracy for Algorithms 1 and [2, Algorithm 1], respectively.
Algorithms 1 performs 10.92 subspace iterations on average, smaller than the average
number of iterations 14.95 by [2, Algorithm 1]. The average run-times on this data
set are quite similar, specifically 10.48 and 9.65 s for Algorithms 1 and [2, Algorithm
1], respectively.
Table 2gives more detailed insight into the performance of the algorithms. It
reports the accuracy of the results returned by Algorithm 1 and [2, Algorithm 1], as
well as the average number of iterations and the average run-times, when only those
examples in the data set for which the rank of Ris greater than 14,20,30,40,50,
respectively, are taken into account. (Note that for all examples we have Rank(R) >
14.) We call a result returned by one of the subspace frameworks accurate if the rela-
tive error |fSF −fBB|/fBB is less than the prescribed tolerance 10−12. According to
Table 2 The accuracy of Algorithm 1 and [2, Algorithm 1] with respect to Rank(R) on the DH systems of
order 500. The second column lists the number of examples in the data set that satisfy the rank constraint,
whereas the third column lists the number of examples among these for which the returned result has error
less than the prescribed tolerance. The average number of iterations and the average run-times (in seconds)
are listed in the fifth and sixth columns
Rank # accurate Accuracy % # iter. Run-time
of R> # Ex. Alg. 1 [2]Alg.1[2]Alg.1[2]Alg.1[2]
14 2000 1929 1829 96.45% 91.45% 10.92 14.95 10.48 9.65
20 1853 1807 1744 97.52% 94.12% 10.39 14.13 9.55 8.98
30 1619 1595 1562 98.52% 96.48% 9.75 13.08 8.96 8.43
40 1390 1382 1363 99.42% 98.06% 9.20 12.09 8.61 8.08
50 1162 1157 1141 99.57% 98.19% 8.85 11.39 8.47 7.89
Adv Comput Math (2020) 46:6 Page 11 of 28 6
the figures in the table, not only the accuracy of the algorithms improves as the rank
of Rincreases, but also the average number of iterations and the average run-time
decrease.
2.3.2 The FE model of a disk brake
The only large-scale computation required by Algorithm 1 (after the preprocessing
step that yield initial interpolation points) is the solution of the linear systems
T(iω) X =Band T(iω)2Y=B(2.16)
at a given ω∈Rin lines 2 and 7, where T(iω) is as in (2.6). For the DH system of
the disk brake in (1.4)and(1.5), the mass matrix Mand the stiffness matrix K()
and thus the matrix Q−1are sparse, but Qturns out to be dense. Trying to invert Q−1
and/or solve a linear system with the coefficient matrix Qis computationally very
expensive and would require full matrix storage.
This difficulty can be avoided by exploiting that
((J −R)Q −iωI)−1B=Q−1((J −R) −iωQ−1)−1B.
Hence, to compute X, Y as in (2.16), we proceed as follows.
1. We first solve ((J −R) −iωQ−1)
X=Bfor
X,andsetX=Q−1
X.
2. Then we solve ((J −R) −iωQ−1)
Y=Xfor
Y,andsetY=Q−1
Y.
A second observation that further speeds up the computations is the particular struc-
ture of the coefficient matrix ((J −R) −iωQ−1)with N=0. Setting
M(iω;) :=
iωM +D() +G(),wehave
(J −R) −iωQ−1=−
M(iω;) −K()
K() −iωK() .
Hence, to solve ((J −R) −iωQ−1)Z =Wfor a given W=WT
1WT
2Tand the
unknown Z=ZT
1ZT
2Twith W1,W
2,Z
1,Z
2having all equal number of rows,
we perform a column block permutation and then eliminate the lower left block to
obtain −K() −
M(iω;)
0K() +iω
M(iω;) Z2
Z1=W1
W2−iωW1.
At every subspace iteration, the highest costs arise from the computation of the LU
factorizations of the sparse matrices K() and K() +iω
M(iω;).
We have applied Algorithm 1 to estimate the non-Hermitian stability radius
r(R;B,BT)for a DH system of the form (1.4)and(1.5) resulting from a FE brake
model with N=0. In this example, G(), K(), D(), M ∈R4669×4669 yielding
J,R,Q ∈R9338×9338, whereas B∈R9338×3.
The plot of the computed estimate values for r(R;B,BT)with respect to the rota-
tion speed ∈[2.5,1800]is presented on the left in Fig. 1. It seems that r(R;B,BT)
as a function of exhibits nonsmoothness at several points, but such points in the
figure are usually local or global maximizers. Two notable maximizers where non-
smoothness is present are near =1120 and =1590. Denoting the former of
Adv Comput Math (2020) 46:6
Page 12 of 28
6
0 200 400 600 800 1000 1200 1400 1600 1800
2
4
6
8
10
12 10-3
1100 1105 1110 1115 1120 1125 1130
4
6
8
10
12 10-3
Fig. 1 (Left) Plot of the approximation to the stability radius r(R;B,BT)by Algorithm for the DH system
specified by (1.4)and(1.5) with N=0 as it varies with the rotation speed ∈[2.5,1800]. (Right) Plot
is zoomed for ∈[1100,1130]
these maximizers with ∗, the plot of r(R;B,BT)for ∈[1100,1130]containing
∗is zoomed on the right-hand plot in Fig. 1. The nonsmoothness here is due to the
fact that σmax(BTQ(iωI −(J −R)Q)−1B) as a function of ωhas multiple global
maximizers for the particular parameter value ∗. A similar remark applies to several
other points where nonsmoothess is present.
The estimated values for r(R;B,BT)are listed in Table 3for some values of
. For comparison purposes, we also estimate r(R;B,BT)by employing [2, Algo-
rithm 1]. (This requires the solution of T(iω)HY=CHinstead of the solution
of T(iω)2Y=Bin (2.16). It is straightforward to extend the approach discussed
above for the efficient solution of T(iω)X =Bto this setting.) The two algorithms
return the same values of r(R;B,BT)up to prescribed tolerances for all values in
the table. In all cases, 2 subspace iterations are sufficient to achieve the prescribed
accuracy for Algorithm 1. Both algorithms in the end replace the full problem of
order 9338 with considerably smaller reduced systems, in particular of order 102 for
Algorithm 1.
Table 3 Estimated values for the stability radii r(R;B,BT)in the second column by Algorithm 1 and
[2, Algorithm 1] for a few values for the system originating from the FE brake model. The other
columns display ω∗corresponding to arg maxωσmax(BTQ(iωI −(J −R)Q)−1B), the number of subspace
iterations, total run-time (in seconds) and subspace dimension at termination
r(R;B,BT)# iter. Run-time Dimension
Alg. 1 & [2]ω∗Alg. 1 [2]Alg.1[2]Alg.1[2]
5 0.00893 −1.741 ·1052 10 30.9 46.8 102 75
50 0.00904 −1.741 ·1052 13 31.4 55.1 102 84
300 0.00752 −1.742 ·1052 3 31.5 31.2 102 54
1100 0.00834 −1.789 ·1052 7 30.8 39.7 102 66
1200 0.00407 −1.742 ·1052 8 30.9 41.9 102 69
Adv Comput Math (2020) 46:6 Page 13 of 28 6
3 Estimation of the hermitian restricted stability radius
In this section, we focus on Hermitian structured perturbations in the dissipation
matrix, in particular approximation of rHerm(R;B) as in (1.6). The algorithms
we propose exploit the following characterization of rHerm(R;B) established in
[21, Theorem 4.9], which is originally presented in terms of an orthonormal basis
U(λ)for the kernel of (I −BB+)T (λ),whereT(λ)is as in (2.6)andB+denotes the
Moore-Penrose pseudoinverse of B. However, it turns out that U(λ) does not have to
be orthonormal, rather the theorem can be stated in terms of any basis for the kernel
of (I −BB+)T (λ). In our version, given in Theorem 3 below, we employ a particular
basis that simplifies the formulas and facilitates the computation.
Theorem 3 For an asymptotically stable DH system of the form (1.1), and a
restriction matrix B∈Cn×mof full column rank, let
1. T(λ)in (2.6)at a given λ∈Cbe such that T(λ)is invertible,
2. L(λ) be a lower triangular Cholesky factor of
H0(λ) := (BHQT (λ)−1B)H(BHQT (λ)−1B), (3.1)
i.e., L(λ) is a lower triangular matrix satisfying
H0(λ) =L(λ)L(λ)H,
3. H0(λ) := L(λ)−1L(λ)−H,
4. H1(λ) := i(
H1(λ) −
H1(λ)H), with
H1(λ) :=
L(λ)−1(BHQT (λ)−1B)HL(λ)−H.
Then rHerm(R;B) is finite, and given by
rHerm(R;B) =inf
ω∈Rsup
t∈R
λmin(H0(iω) +tH1(iω))1/2
,
where λmin(·)denotes the smallest eigenvalue of its Hermitian matrix argument, and
the inner supremum is finite if and only if H1(iω) is indefinite.
We first describe a numerical technique for small-scale problems in Section 3.1
and then develop a subspace framework in Section 3.2, both based on the eigenvalue
optimization characterization of rHerm(R;B) in Theorem 3. Note that in this section,
Bis always assumed to be of full column rank.
3.1 An algorithm for small-scale problems
The eigenvalue optimization characterization of rHerm(R;B) is a min-max prob-
lem, where the inner maximization problem is concave; indeed, it can alternatively
be expressed as a semidefinite program (SDP). Formally, for a given ω∈R,and
H0(iω), H1(iω) representing the Hermitian matrices defined in Theorem 3, let us
consider
ηHerm(R;B,iω) := sup
t∈R
λmin(H0(iω) +tH1(iω))
=sup{z|z, t ∈Rs.t. H0(iω) +tH1(iω) −zI ≥0}(3.2)
Adv Comput Math (2020) 46:6
Page 14 of 28
6
so that
rHerm(R;B) =inf
ω∈RηHerm(R;B,iω). (3.3)
In the following subsections, we solve the minimization problem inside the square
root in (3.3).
3.1.1 Inner maximization problems
The characterization in the second line of (3.2) is a linear convex SDP. Some of
the most widely used techniques to solve a linear convex SDP are different forms
of interior-point methods. Implementations of some of these by the interior-point
methods available through the package cvx [13,14].
An alternative is to employ the software package eigopt [22] for the eigenvalue
optimization problem in the first line of (3.2). This software package is based on
an approach that approximates the eigenvalue function with a piecewise quadratic
function that lies globally above the eigenvalue function. It requires a global upper
bound for the second derivative of the eigenvalue function, which can be chosen as
any slightly positive real number (e.g., 10−6) for the concave eigenvalue function in
the first line of (3.2).
In our experience, eigopt performs the computation of ηHerm(R;B,iω) signifi-
cantly faster than cvx. The only downside is that an interval containing the optimal
tfor (3.2) must be supplied to eigopt, whereas such an interval is not needed by
cvx. In practice, we employ eigopt due to its efficiency.
3.1.2 Outer minimization problems
The minimum of ηHerm(R;B,iω) with respect to ω∈Ryields the square of the
radius rHerm(R;B), and it turns out that the minimizing ω∈Rcorresponds to the
point iωthat first becomes an eigenvalue on the imaginary axis under the smallest
perturbation possible [21, Theorem 4.9]. This is a nonconvex optimization problem.
Indeed, the objective ηHerm(R;B,iω) mayevenblowupatsomeω.
We again resort to eigopt for the minimization of ηHerm(R;B,iω).Forthesake
of completeness, a formal description is provided in Algorithm 2 below, where we
use the abbreviations
ηj:= ηHerm(R;B,iωj)and η
j:= d
dω ηHerm(R;B,iω)ω=ωj
.
Every convergent subsequence of the sequence {ωk}in Algorithm 2 is guaranteed
to converge to a global minimizer of ηHerm(R;B,iω) provided γis chosen suffi-
ciently small. Moreover, by continuity, it also follows that ηHerm(R;B,iωk)in the
limit approaches rHerm(R;B) for all γvalues sufficiently small. In our experiments,
modestly small values of γ,suchasγ=−20, seem to lead to convergence.
Adv Comput Math (2020) 46:6 Page 15 of 28 6
At step kof the algorithm, ηHerm(R;B,iω) and its derivative need to be computed
at ωk. For instance, one can rely on one of the two approaches (cvx or eigopt)in
Section 3.1.1 for the computation of ηHerm(R;B,iωk), and employ a finite difference
formula to approximate its derivative, which would require an additional computation
of ηHerm(R;B,iω) at an ωclose to ωk.
3.2 An interpolation-based subspace framework for large-scale problems
In this section, we consider the Hermitian restricted stability radius for large DH
systems. The characterization in Theorem 3 is in terms of the matrix-valued functions
H0(iω) and H1(iω), which are of small size provided that Bhas few columns. The
large-scale nature of this characterization is hidden in the matrix-valued function
T(iω) =(J −R)Q −iωI, because the computation of H0(iω) and H1(iω) requires
the solution of the linear system T(iωk)Z =B.
3.2.1 Derivation of the reduced problems
To cope with the large-scale setting, we benefit from structure-preserving two-sided
projections similar to those described in Section 2.2. In particular, for a given matrix
Vkwith orthonormal columns, we define Wkas in (2.9), Jk,R
k,Q
kas in (2.12), and
Bkas in (2.13), respectively. Again, Vk,Wkrefer to the right, left subspaces spanned
by the columns of Vk,W
k.
However, we no longer have a tool such as Theorem 1 to establish interpolation
results, because there is no apparent transfer function, as there is indeed no apparent
linear PH system that can be tied to the eigenvalue optimization characterization. The
following simple observation turns out to be very useful. The subsequent arguments
make use of the matrix-valued function
Tk(λ) := (Jk−Rk)Qk−λI (3.4)
for λ∈Cassociated with the reduced problems.
Adv Comput Math (2020) 46:6
Page 16 of 28
6
Lemma 1 Consider a DH model (1.1), and a reduced model ˙xk=(Jk−Rk)Qkxk
with coefficients as in (2.12)defined in terms of a given Vkand Wkas in (2.9).
Letting T(λ)and Tk(λ) be as in (2.6)and (3.4), respectively, we then have Tk(λ) =
WH
kT(λ)V
kfor all λ∈C.
Proof From the definitions of Jk,R
k,Q
kwe obtain
Tk(λ) =(WH
k(J −R)Wk)V H
kQVk−λI
=WH
k(J −R)QVk−λWH
kVk
=WH
k((J −R)Q −λI)Vk=WH
kT(λ)V
k,
where we have employed (2.10) in the second equality.
In the following, we develop a subspace framework in which the Hermite inter-
polation properties between the original problem defined in terms of T(λ) and the
reduced problem defined in terms of Tk(λ) hold. For this, we first show an auxiliary
interpolation result between BHQT (λ)−1Band its reduced counterpart.
Lemma 2 Consider a DH model (1.1), and a reduced model ˙xk=(Jk−Rk)Qkxk
with coefficients as in (2.12)defined in terms of a given Vkand Wkas in (2.9).
Furthermore, let T(λ)and Tk(λ) be as in (2.6)and (3.4), respectively. For a given
λ∈Csuch that T(
λ) and Tk(
λ) are invertible, the following assertions hold:
(i) If Ran(T (
λ)−1B) ⊆Vk,thenBHQT (
λ)−1B=BH
kQkTk(
λ)−1Bk.
(ii) Additionally, if Ran(T (
λ)−2B) ⊆Vkand the orthonormal basis Vkfor Vkis of
the form Vk=
Vk
Vk, where the columns of
Vkform an orthonormal basis
for Ran(T (
λ)−1B),then BHQT (
λ)−2B=BH
kQkTk(
λ)−2Bk.
Proof (i) If Ran(T (
λ)−1B) ⊆Vk,then
BHQT (
λ)−1B=BHQVkVH
kT(
λ)−1B=BHWkVH
kQVkVH
kT(
λ)−1B
=BH
kQkVH
kT(
λ)−1B,
(3.5)
where the first equality is due to the fact that VkVH
kis the orthogonal projector
onto Vk, and the second follows, since WkVH
kis a projector onto Ran(QVk).
We complete the proof by showing that VH
kT(
λ)−1B=Tk(
λ)−1Bk.Tothis
end, let Z:= T(
λ)−1B,andZkbe such that VkZk=Z. (There exists a unique
Zkwith this property, because Ran(Z) ⊆Vk.) Then T(
λ)Z =Bimplies that
T(
λ)VkZk=B, and thus WH
kT(
λ)VkZk=WH
kB=Bk. Hence, by Lemma 1,
we see that Zk=Tk(
λ)−1Bk, which in turn implies
VH
kT(
λ)−1B=VH
kZ=VH
k(VkZk)=Tk(
λ)−1Bk. (3.6)
(ii) Following the steps at the beginning of the proof of part (i), in particular from
a line of reasoning similar to that in (3.5), we have
BHQT (
λ)−2B=BH
kQkVH
kT(
λ)−2B.
Adv Comput Math (2020) 46:6 Page 17 of 28 6
It remains to show that VH
kT(
λ)−2B=Tk(
λ)−2Bk. To this end, we exploit
VH
kT(
λ)−2B=(V H
kT(
λ)−1
Vk)(
VH
kT(
λ)−1B). (3.7)
Now define Z:= T(
λ)−1
Vkand Zksuch that VkZk=Z(once again such a
Zkexists uniquely, because Ran(Z) ⊆Vk)sothatT(
λ)VkZk=
Vk. Then, due
to the property WH
kVk=I
k,wehaveWH
kT(
λ)VkZk=WH
k
Vk=I
k,m,where
I
k,m is the matrix consisting of the first mcolumns of the
k×
kidentity matrix
I
kand
k:= dim Vk≥m. By Lemma 1, this implies that Zk=Tk(
λ)−1I
k,m,so,
for the term inside the first parenthesis on the right-hand side of (3.7), we have
VH
kT(
λ)−1
Vk=VH
kZ=VH
k(VkZk)=Tk(
λ)−1I
k,m. (3.8)
For the term inside the second parenthesis on the right-hand side of (3.7), we
make use of the following observation:
VH
kT(
λ)−1B=
VH
k
VH
kT(
λ)−1B=
VH
kT(
λ)−1B
0,(3.9)
where the last equality follows, since the columns of
Vkform an orthonormal
basis for Ran(T (
λ)−1B). Inserting these in (3.7), we obtain
VH
kTk(
λ)−2Bk=(V H
kT(
λ)−1
Vk)(
VH
kT(
λ)−1B)
=Tk(
λ)−1I
k,m(
VH
kT(
λ)−1B)
=Tk(
λ)−1I
k(V H
kT(
λ)−1B)
=Tk(
λ)−1I
k(Tk(
λ)−1Bk)=Tk(
λ)−2Bk,
where in the second and third equality we exploit (3.8)and(3.9), respectively,
and in the fourth equality we employ the equivalence in (3.6).
Our reduced problems are expressed in terms of the reduced versions of H0(λ),
L(λ),andH1(λ) defined via
Hk,0(λ) := Lk(λ)−1Lk(λ)−H,(3.10)
with Lk(λ) denoting a lower triangular Cholesky factor of
Hk,0(λ) := (BH
kQkTk(λ)−1Bk)H(BH
kQkTk(λ)−1Bk), (3.11)
and Hk,1(λ) := i(
Hk,1(λ) −
Hk,1(λ)H)with
Hk,1(λ) := Lk(λ)−1(BH
kQkTk(λ)−1Bk)HLk(λ)−H. (3.12)
To ensure the uniqueness of L(λ) and Lk(λ), we define them as the Cholesky factors
of
H0(λ) and
Hk,0(λ) with real and positive entries along the diagonal. We emphasize
here that the matrix functions H0(λ),L(λ),H1(λ) and their reduced counterparts are
of the same size. However, the reduced matrix-valued functions are defined in terms
of Tk(λ), which is small if the dimensions of Vk,Wkare small, in contrast to T(λ)in
the definitions of H0(λ),L(λ),H1(λ).
Adv Comput Math (2020) 46:6
Page 18 of 28
6
Our goal is to come up with a reduced counterpart of ηHerm(R;B,iω),that
Hermite-interpolates the full function at prescribed points. For a given Vk,andWkas
in (2.9), we introduce
ηHerm
k(R;B,iω) := sup
t∈R
λmin(Hk,0(iω) +tHk,1(iω)).
Next, we establish that the quantity ηHerm
k(R;B,iω) is independent of the choice
of the basis for the subspace Vk. For this result, we introduce the notation
ηHerm
Vk,Wk(R;B,iω)) to emphasize the particular choices of the bases Vk,W
kused in the
definition of ηHerm
k(R;B,iω). Similarly, we indicate the choices of the bases used in
the definitions of Jk,R
k,Q
k,andBkwith the help of the notations JWk,RWk,QVk,
and BWk.
Lemma 3 Let the columns of Vkand
Vkform orthonormal bases for the subspace
Vk, and let Wk:= QVk(V H
kQVk)−1,
Wk:= Q
Vk(
VH
kQ
Vk)−1.Then
ηHerm
Vk,Wk(R;B,iω) =ηHerm
Vk,
Wk(R;B,iω)
for all ω∈R.
Proof The proof makes use of the characterization
ηHerm(R;B,iω) =
inf{2|=H,iω∈(J −(R +BBH))Q}2(3.13)
shown in [21, Theorem 4.9].
Now since the columns of Vkand
Vkform bases for the same space, there exists a
unitary matrix Psuch that
Vk=VkP. Furthermore,
Wk=Q(VkP )(P HVH
kQVkP)
−1=QVkPPH(V H
kQVk)−1P=WkP.
The assertion then follows from (3.13) and the following set of equivalences:
iω∈(JWk−RWk)QVk−(BWkBH
Wk)QVk⇐⇒
det(WH
k(J −R)Q −BBHQ−iωIVk)=0⇐⇒
det(
WH
k(J −R)Q −BBHQ−iωI
Vk)=0⇐⇒
iω∈(J
Wk−R
Wk)Q
Vk−(B
WkBH
Wk)Q
Vk.
Above, the first and third equivalences are due to the identities in (2.10), whereas the
second equivalence is due to
Vk=VkPand
Wk=WkP.
We now prove our main interpolation result.
Theorem 4 (Hermite Interpolation of Hermitian Backward Errors) Consider a DH
model (1.1), and a reduced model ˙xk=(Jk−Rk)Qkxkwith coefficients as in (2.12)
defined in terms of a given Vkand Wkas in (2.9). Furthermore, let T(λ)and Tk(λ)
be as in (2.6)and (3.4), respectively. For a given ω∈R, suppose that the subspace
Vk:= Ran(Vk)is such that
Ran(T (iω)−1B) ⊆Vk,Ran(T (iω)−2B) ⊆Vk. (3.14)
Adv Comput Math (2020) 46:6 Page 19 of 28 6
(i) The quantity ηHerm(R;B,iω) is finite if and only if ηHerm
k(R;B,iω) is finite. If
ηHerm(R;B,iω) is finite, then
ηHerm(R;B,iω) =ηHerm
k(R;B,iω). (3.15)
(ii) Moreover, if ηHerm(R;B,iω) and ηHerm
k(R;B,iω) are differentiable at ω,then
we have
d
dω ηHerm(R;B,iω)ω=ω
=d
dω ηHerm
k(R;B,iω)ω=ω
. (3.16)
Proof Without loss of generality, we may assume that the matrix Vkis such
that Vk=
Vk
Vk,with the columns of
Vkforming an orthonormal basis for
Ran(T (iω)−1B). By Lemma 3, it suffices to prove the claims for such a particular
choice of orthonormal basis.
(i) By the definitions of
H0(λ),
Hk,0(λ) in (3.1), (3.11), respectively, and by part
(i) of Lemma 2, we have
H0(iω) =(BHQT (iω)−1B)H(BHQT (iω)−1B)
=(BH
kQkTk(iω)−1Bk)H(BH
kQkTk(iω)−1Bk)=
Hk,0(iω).
This also implies that L(iω) =Lk(iω) due to the uniqueness of the
Cholesky factors of
H0(iω),
Hk,0(iω),soH0(iω) =L(iω)−1L(iω)−H=
Lk(iω)−1Lk(iω)−H=Hk,0(iω). Furthermore, recalling the definitions of
H1(λ), H1(λ) in part 4 of Theorem 3, the definitions of
Hk,1(λ), Hk,1(λ) in
(3.12) and the line above (3.12), we have
H1(iω) =L(iω)−1(BHQT (iω)−1B)HL(iω)−H
=Lk(iω)−1(BH
kQkTk(iω)−1Bk)HLk(iω)−H=
Hk,1(iω)
and H1(iω) =Hk,1(iω). From the latter equality and Theorem 3, we deduce
that ηHerm(R;B,iω) is finite if and only if ηHerm
k(R;B,iω) is finite. Addition-
ally, if ηHerm(R;B,iω) is finite, then
ηHerm(R;B,iω) =maxt∈Rλmin(H0(iω) +tH1(iω))
=maxt∈Rλmin(Hk,0(iω) +tHk,1(iω)) =ηHerm
k(R;B,iω),
completing the proof of (3.15).
(ii) Now let us suppose that ηHerm(R;B,iω) and ηHerm
k(R;B,iω) are differentiable
at ω. To prove the interpolation property in the derivatives, we benefit from ana-
lytical expressions. It follows from the standard perturbation theory of simple
eigenvalues (e.g., see [20]) that
d
dω ηHerm(R;B,iω)ω=ω
=vHd
dω H0(iω) +
tH1(iω)ω=ωv,
d
dω ηHerm
k(R;B,iω)ω=ω
=vHd
dω Hk,0(iω) +
tHk,1(iω)ω=ωv,
(3.17)
Adv Comput Math (2020) 46:6
Page 20 of 28
6
where
tis such that
t∈arg max
tλmin(H0(iω) +tH1(iω)) =arg max
tλmin(Hk,0(iω) +tHk,1(iω))
and vis a unit eigenvector corresponding to λmin(H0(iω) +
tH1(iω)) =λmin(
Hk,0(iω) +
tHk,1(iω)). Thus, it suffices to prove
d
dωH0(iω)ω=ω
=d
dωHk,0(iω)ω=ω
,d
dωH1(iω)ω=ω
=d
dωHk,1(iω)ω=ω
.
We first establish the interpolation properties between the derivatives of H0(iω)
and Hk,0(iω) at ω=ω. To this end, observe that parts (i) and (ii) of Lemma 2
imply that
d
dω
H0(iω)ω=ω
=i(BHQT (iω)−1B)H(BHQT (iω)−2B)
−i(BHQT (iω)−2B)H(BHQT (iω)−1B)
=i(BH
kQkTk(iω)−1Bk)H(BH
kQkTk(iω)−2Bk)
−i(BH
kQkTk(iω)−2Bk)H(BH
kQkTk(iω)−1Bk)=d
dω
Hk,0(iω)ω=ω
.
(3.18)
Next, for a given ω, consider the Cholesky factorizations
H0(iω) =
L(iω)L(iω)Hand
Hk,0(iω) =Lk(iω)Lk(iω)H. By differentiating these two
equations and exploiting (3.18), we have ]selectfont
L(iω) d
dωL(iω)H+d
dωL(iω)L(iω)Hω=ω
=
Lk(iω) d
dωLk(iω)H+d
dωLk(iω)Lk(iω)Hω=ω
=
L(iω) d
dωLk(iω)H+d
dωLk(iω)L(iω)Hω=ω
,
wherewehaveusedL(iω) =Lk(iω), established in part (i), in the sec-
ond equality. Thus, both d
dωL(iω)|ω=ωand d
dωLk(iω)|ω=ωare lower triangular
solutions of the matrix equation
d
dω
H0(iω)ω=ω
=d
dω
Hk,0(iω)ω=ω
=L(iω)XH+XL(iω)H,
which has a unique lower triangular solution with real entries along the
diagonal, so d
dωL(iω)|ω=ω=d
dωLk(iω)|ω=ω. Furthermore, it follows from
the definitions of H0(iω) and Hk,0(iω) in part 3 of Theorem 3 and (3.10)
that L(iω)H0(iω)L(iω)H=I=Lk(iω)Hk,0(iω)Lk(iω)H, which we dif-
ferentiate at ω=ω, and exploit L(iω) =Lk(iω),d
dωL(iω)|ω=ω=
d
dωLk(iω)|ω=ω,H0(iω) =Hk,0(iω) to deduce that the equality d
dωH0(iω)|ω=ω
=d
dωHk,0(iω)|ω=ωholds.
Finally, by exploiting L(iω)
H1(iω)L(iω)H=(BHQT (iω)−1B)H,which
follows from the definition of H1(iω), Hk,1(iω) in part 4 of Theorem 3, we
Adv Comput Math (2020) 46:6 Page 21 of 28 6
obtain
d
dω L(iω)
H1(iω)L(iω)Hω=ω
=−i(BHQT (iω)−2B)H
=−i(BH
kQkTk(iω)−2Bk)H
=d
dω Lk(iω)
Hk,1(iω)Lk(iω)Hω=ω
,
where we exploit L(iω) =Lk(iω),d
dωL(iω)|ω=ω=d
dωLk(iω)|ω=ω,
H1(iω) =
Hk,1(iω) to deduce d
dω
H1(iω)|ω=ω=d
dω
Hk,1(iω)|ω=ω. It immediately follows
that d
dωH1(iω)|ω=ω=d
dωHk,1(iω)|ω=ω, and the proof of (3.16) is complete.
Remark 1 The function ηHerm(R;B,iω) is differentiable at ωwhenever
1. ηHerm(R;B,iω) is finite, equivalently H1(iω) is indefinite,
2. the global minimum of λmin(H0(iω) +tH1(iω)) is attained at a unique t,and
3. λmin(H0(iω) +
tH1(iω)) is simple, where
t:= arg maxt∈Rλmin(H0(iω) +
tH1(iω)).
These conditions guarantee also the differentiability of ηHerm
k(R;B,iω) at ωpro-
vided that the subspace inclusion in (3.14) holds. This latter differentiability property
is due to H0(iω) =Hk,0(iω) and H1(iω) =Hk,1(iω) from part (i) of Theorem 4.
Additionally, when the concave function g(t) := λmin(H0(iω) +tH1(iω)) attains
its maximum, the maximizer is nearly always unique. The function g(t) is the min-
imum of mreal analytic functions [19,25], each corresponding to an eigenvalue of
H0(iω)+tH1(iω).Ifg(t) does not have a unique maximizer, then at least one of these
real analytic functions must be constant and equal to ηHerm(R;B,iω) =maxtg(t).
Thus, a simple sufficient condition that ensures the uniqueness of the maximizer is
that H1(iω) has full rank, in which case all eigenvalues of H0(iω) +tH1(iω) blow
up, either to ∞or −∞,ast→∞implying that each of the real analytic functions
is non-constant.
3.2.2 The subspace framework for
r
Herm(
R
;
B
)
The Hermite interpolation result of Theorem 4 immediately suggests the subspace
framework in Algorithm 3 for the approximation of rHerm(R;B). This resembles
the structure-preserving subspace framework to estimate the non-Hermitian stability
radius r(R;B,C), especially in the way the subspaces Vk,Wkare built. At every iter-
ation, a reduced problem is solved in line 5 as discussed in Section 3.1, specifically
by employing Algorithm 2. Letting ωbe a global minimizer of the reduced problem,
the subspaces are expanded so that the original function ηHerm(R;B,iω) is Hermite
interpolated by its reduced counterpart at ω=ω.
Assuming that the sequence {ωk}converges to a minimizer ω∗of the function
ηHerm(R;B,iω) such that
Adv Comput Math (2020) 46:6
Page 22 of 28
6
1. H1(iω∗)is indefinite with full rank, and
2. λmin(H0(iω∗)+t∗H1(iω∗)) is simple, with t∗:= arg maxt∈Rλmin(H0(iω∗)+
tH1(iω∗)),
the analysis in [3] shows that the sequence {ωk}converges at a superlinear rate.
As discussed in Remark 1, the two conditions above ensure the differentiability of
ηHerm(R;B,iω) and ηHerm
k(R;B,iω) for all large kin a neighborhood of ω∗.This
differentiability property is essential for the convergence analysis in [3].
3.3 Numerical experiments
In this section, we present numerical results for the estimation of the Hermitian sta-
bility radius. Most of the examples are large-scale, and, hence, involve the application
of the proposed subspace framework (Algorithm 3). However, recall that Algorithm 3
uses Algorithm 2 for solving reduced subproblems. For Algorithm 3, we use the same
stopping criteria as in Section 2.3, but now fk:= minω∈RηHerm
k(R;B,iω),andwe
set kmax =80 and ε=10−6. The larger value for εis due to the fact that we compute
the derivatives of ηHerm
k(R;B,iω) with respect to ωby means of forward differences,
which yield approximate values accurate up to about half of the double precision.
We form the initial subspaces nearly as described in Section 2.3 for the non-
Hermitian stability radius. The only difference is that, in part (iii), the points
z1,...,z
ρare permuted into zj1,...,z
jρso that ηHerm(R;B,izj1)≤ ··· ≤
ηHerm(R;B,izjρ). For the synthetic examples in Section 3.3.1, the parameter val-
ues l=−200, u=0, ρ=κ=5 are used, and z1,...,z
ρare set equal to the
eigenvalues of (J −R)Q closest to iω1,...,iωρ. As for the disk brake example in
Section 3.3.2,weusel=−2·105,u=−1.7 ·105,ρ=κ=15, as well as zj=iωj
for j=1,...,ρ.
Adv Comput Math (2020) 46:6 Page 23 of 28 6
Fig. 2 A DH system of order 20 with random system matrices. Eigenvalues of (J −R)Q are marked with
red asterisks, whereas the black regions represent the eigenvalues of all matrices of the form (J −(R +
BBH))Q for 100000 randomly chosen real symmetric matrices with 2=rHerm(R;B). The red
circle marks 1.9794i, where 1.9794 is the computed global minimizer of ηHerm(R;B,iω) over ω∈R.
3.3.1 Synthetic examples
We first present a numerical result for a small dense random example, where
J,R,Q ∈R20×20 and B∈R20×2. Application of Algorithm 2 to this example yields
rHerm(R;B) =0.0501, and the point that is first reached on the imaginary axis under
the smallest perturbation is 1.9794i. Figure 2illustrates that there exist real symmet-
ric matrices with 2=0.0501 such that 1.9794i is contained (nearly) in the
spectrum of (J −(R +BBH))Q.
Next we experiment with larger dense random matrices, specifically with random
A∈Rn×n,B∈Rn×2. The results of Algorithm 3 are reported in Tables 4and 5.For
n=1000, as shown in Table 4, the direct application of Algorithm 2 and the subspace
framework (Algorithm 3) return exactly the same values for rHerm(R;B) up to the
prescribed tolerances, yet the subspace framework already needs less computing time
than Algorithm 2. Listed in this table is also the unstructured distance to instability
((J −R)Q) defined in Section 2.3.1. It appears that ((J −R)Q) overestimates
rHerm(R;B) considerably.
Table 4 Performance of Algorithm 3 to estimate rHerm(R;B) on dense random DH systems of order
1000. The results from Algorithm 2, the non-Hermitian stability radius r(R;B,BH)and the unstructured
distance to instability ((J −R)Q) are also included. The fifth column contains the number of subspace
iterations, the sixth the run-times (in seconds)
rHerm(R;B) r(R;B,BH) ((J −R)Q) # iter. Run-time
#Alg.2&3 Alg.3Alg.2 Alg.3
1 0.012922 0.011263 2.04 7 2090.0 463.5
2 0.010110 0.010104 0.36 1 1544.1 91.8
3 0.014067 0.012826 3.63 5 2840.2 335.6
4 0.009605 0.007189 0.94 11 2926.5 739.6
Adv Comput Math (2020) 46:6
Page 24 of 28
6
Table 5 Performance of Algorithm 3 to estimate rHerm(R;B) on dense random DH systems of order
n=2000 and n=4000. The fifth column contains the number of subspace iterations, and the sixth the
run-times (in seconds)
rHerm(R;B) r(R;B,BH)# iter. Run-time
n# Alg. 3 Alg. 3 Alg. 3
2000 1 0.010217 0.005231 6 408.6
2000 2 0.003835 0.003754 13 811.8
2000 3 0.012456 0.010876 5 369.2
2000 4 0.010096 0.010031 3 221.8
4000 1 0.008982 0.008705 1 213.1
4000 2 0.008376 0.006803 3 343.2
4000 3 0.010400 0.008646 2 267.8
4000 4 0.002068 0.002055 11 828.9
For larger systems, Algorithm 2 becomes computationally too expensive, so we
do not report the results of Algorithm 2. For the same reason, we also do not report
((J −R)Q) for larger systems. Rather, we report only the results of Algorithm 3
for n=2000 and n=4000 in Table 5, as well as the non-Hermitian stability radii.
We remark that in these experiments the number of subspace iterations to reach the
prescribed accuracy is usually small and seems independent of n.
All of these examples involve optimization of highly nonconvex functions.
Figure 3depicts ηHerm(R;B,iω) as a function of ωwith the solid curve for the
first example in Table 4with n=1000. The same figure also depicts the reduced
function ηHerm
k(R;B,iω) with the dashed curve for the same example at termination
after 7 subspace iterations. Even though the reduced function ηHerm
k(R;B,iω) uses
projected matrices onto 48 dimensional subspaces, it captures the original function
ηHerm(R;B,iω) remarkably well near the global minimizer ω∗=−70.9581.
Fig. 3 Application of
Algorithm 3 to estimate
rHerm(R;B) on a dense random
DH system of order 1000. The
full and reduced functions
ηHerm(R;B,iω) and
ηHerm
k(R;B,iω) at termination
after 7 subspace iterations are
plotted with the solid and
dashed curves, respectively. The
point (ω∗,η
Herm(R;B,iω∗))
with ω∗denoting the global
minimizer of ηHerm(R;B,iω) is
marked with a circle
-200 -150 -100 -50 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4 10-3
Adv Comput Math (2020) 46:6 Page 25 of 28 6
Table 6 Estimated values for the Hermitian stability radii rHerm(R;B) by Algorithm 3 in the second
column for the FE model of a disk brake of order 9338 for a few values of .Thecolumnω∗displays
the computed global minimizer of ηHerm(R;B,iω) by the algorithm, whereas the last columns depicts the
total run-time (in seconds) for Algorithm 3
r
Herm(R;B) r(R;B,BT)ω
∗Run-time
5 0.01038 0.00893 −1.938 ×105227.4
50 0.01001 0.00904 −1.938 ×105222.3
100 0.00988 0.00914 −1.938 ×105223.9
1100 0.00837 0.00834 −1.789 ×105228.5
1200 0.00408 0.00407 −1.742 ×105224.1
3.3.2 FE model of a disk brake
We also applied our implementation of Algorithm 3 to the FE model of the disk
brake described in Section 2.3.2. The estimated values for the Hermitian structured
radius rHerm(R;B) along with the associated minimizer ω∗of ηHerm(R;B,iω),as
well as the run-time (in seconds) are listed in Table 6for a few values of .Only
one subspace iteration is performed in these examples and the subspace dimension at
termination is 96. It is worth comparing the values of rHerm(R;B) in this table with
those for the non-Hermitian stability radius r(R;B,BT), also included in the table.
As expected, the listed values for the Hermitian structured stability radii are larger.
The estimates for the two stability radii differ by significant amounts for smaller
values of , but only slightly at larger .
The results show that the FE model of the disk brake is close to instability for
all . They indicate that further damping via shims would be helpful to avoid the
instabilities due to uncertainties.
4 Concluding remarks
We have proposed subspace frameworks to estimate the stability radii for large-
scale dissipative Hamiltonian systems. At every iteration, we apply DH structure-
preserving projections to small subspaces, then compute the corresponding stability
radii for the resulting reduced system by exploiting its eigenvalue optimization
characterization. We expand the subspaces used in the projections so that Hermite
interpolation properties between the objective eigenvalue function of the full and the
reduced problems are attained at the optimizer of the reduced problem.
The proposed frameworks seem to converge to local optimizers, but a formal proof
of this observation is open at the moment. They do converge at a superlinear rate in
theory under the local convergence assumption. With the purpose of improving the
likelihood of convergence to a global maximizer, we have initiated the subspaces to
attain Hermite interpolation at several points on the imaginary axis between the full
and initial reduced problems.
Adv Comput Math (2020) 46:6
Page 26 of 28
6
MATLAB implementations of the proposed algorithms and subspace frameworks,
as well as some of the data are made publicly available.2
A related research direction that is currently investigated is the maximization of
the stability radii when J,R,Q depend on parameters, e.g., the FE model of the disk
brake depends on the rotation speed .
Acknowledgments The authors thank two anonymous referees for helpful comments on initial versions
of this manuscript.
Funding Information Open access funding provided by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long
as you give appropriate credit to the original author(s) and the source, provide a link to the Creative
Commons licence, and indicate if changes were made. The images or other third party material in this
article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use
is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission
directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/
by/4.0/.
References
1. Akay, A.: Acoustics of friction. J. Acoust. Soc. Am. 111(4), 1525–1548 (2002)
2. Aliyev, A., Benner, P., Mengi, E., Schwerdtner, P., Voigt, M.: Large-scale computation of L∞-norms
by a greedy subspace method. SIAM J. Matrix Anal. Appl. 38(4), 1496–1516 (2017)
3. Aliyev, A., Benner, P., Mengi, E., Voigt, M.: A subspace framework for H∞-norm minimization.
arXiv:1508.04214v2 [math.NA]. Submitted to SIAM J. Matrix Anal. Appl (2019)
4. Aliyev, N., Mehrmann, V., Mengi, E.: Computation of stability radii for large-scale dissipative hamil-
tonian systems. Preprint 2018-09, Institut f¨ur Mathematik, TU Berlin, D-10623 Berlin, FRG. http://
www.math.tu-berlin.de/menue/forschung/veroeffentlichungen/preprints 2018 (2018)
5. Benner, P., Lowe, R., Voigt, M.: L∞norm computation for large-scale descriptor systems using
structured iterative eigensolvers. Numerical Algebra Control and Optimization 8(1), 119–133 (2018)
6. Benner, P., Mitchell, T.: Faster and more accurate computation of the H∞norm via optimization.
SIAM J. Sci. Comput. 40(5), A3609–A3635 (2018)
7. Benner, P., Voigt, M.: A structured pseudospectral method for H∞-norm computation of large-scale
descriptor systems. Math. Control Signals Syst. 26(2), 303–338 (2014)
8. Boyd, S., Balakrishnan, V.: A regularity result for the singular values of a transfer matrix and a
quadratically convergent algorithm for computing its L∞-norm. Syst. Control Lett. 15(1), 1–7 (1990)
9. Bruinsma, N.A., Steinbuch, M.: A fast algorithm to compute the H∞-norm of a transfer function
matrix. Syst. Control Lett. 14(4), 287–293 (1990)
10. Freitag, M.A., Spence, A., Van Dooren, P.: Calculating the H∞-norm using the implicit determinant
method. Linear Algebra Appl. 35(2), 619–635 (2014)
11. Gallivan, K., Vandendorpe, A., Van Dooren, P.: Model reduction of mIMO systems via tangential
interpolation. SIAM J. Matrix Anal. Appl. 26(2), 328–349 (2005)
12. Gr¨
abner, N., Mehrmann, V., Quraishi, S., Schr¨
oder, C., von Wagner, U.: Numerical methods for para-
metric model reduction in the simulation of disc brake squeal. Z. Angew. Math. Mech. 96, 1388–1405
(2016). https://doi.org/10.1002/zamm.201500217
2http://home.ku.edu.tr/∼emengi/software/DH-stabradii.html
Adv Comput Math (2020) 46:6 Page 27 of 28 6
13. Grant, M., Boyd, S.: Graph implementations for nonsmooth convex programs. In: Blondel, V., Boyd,
S., Kimura, H. (eds.) Recent Advances in Learning and Control, Lecture Notes in Control and
Information Sciences, pp. 95–110. Springer-Verlag Limited (2008)
14. Grant, M., Boyd, S.: CVX: Matlab Software for disciplined convex programming, version 2.1 http://
cvxr.com/cvx (2014)
15. Gugercin, S., Polyuga, R.V., Beattie, C., van der Schaft, A.J.: Interpolation-based model reduction for
port-Hamiltonian systems. In: Proceedings of the 48th IEEE Conference on Decision and Control,
2009 held jointly with the 2009 28th Chinese Control Conference. CDC/CCC 2009, pp. 5362–5369
(2009)
16. Gugercin, S., Polyuga, R.V., Beattie, C., van der Schaft, A.J.: Structure-preserving tangential
interpolation for model reduction of port-Hamiltonian systems. Automatica 48(9), 1963–1974 (2012)
17. Guglielmi, N., G¨urb¨uzbalaban, M., Overton, M.L.: Fast approximation of the H∞-norm via optimiza-
tion over spectral value sets. SIAM J. Matrix Anal. Appl. 34(2), 709–737 (2013)
18. Jacob, B., Zwart, H.: Linear Port-Hamiltonian Systems on Infinite-Dimensional Spaces Operator
Theory: Advances and Applications, vol. 223. Birkh¨
auser/Springer Basel AG, Basel (2012)
19. Kato, T.: Perturbation Theory for Linear Operators. Springer, Berlin (1995)
20. Lancaster, P.: On eigenvalues of matrices dependent on a parameter. Numer. Math. 6, 377–387 (1964)
21. Mehl, C., Mehrmann, V., Sharma, P.: Stability radii for linear Hamiltonian systems with dissipa-
tion under structure-preserving perturbations. SIAM J. Matrix Anal. Appl. 37(4), 1625–1654 (2016).
https://doi.org/10.1137/16M1067330
22. Mengi, E., Yildirim, E.A., Kilic¸, M.: Numerical optimization of eigenvalues of Hermitian matrix
functions. SIAM J. Matrix Anal. Appl. 35(2), 699–724 (2014). https://doi.org/10.1137/130933472
23. Mitchell, T., Overton, M.L.: Hybrid expansion-contraction: a robust scaleable method for approximat-
ing the H∞norm. IMA J. Numer. Anal. 36(3), 985–1014 (2016)
24. Polyuga, R.V., van der Schaft, A.J.: Structure preserving moment matching for port-Hamiltonian
systems: Arnoldi and Lanczos. IEEE Trans. Autom. Cont. 56(6), 1458–1462 (2011)
25. Rellich, F.: Perturbation Theory of Eigenvalue Problems. Notes on Mathematics and its Applications.
Gordon and Breach, New York (1969)
26. van der Schaft, A.J., Jeltsema, D.: Port-Hamiltonian systems theory: an introductory overview. Found
Trends Syst Control 1(2-3), 173–378 (2014)
27. Van Loan, C.F.: How near is a stable matrix to an unstable matrix? Contemporary Math. 47, 465–478
(1985)
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps
and institutional affiliations.
Affiliations
Nicat Aliyev1·Volker Mehrmann2·Emre Mengi3
Nicat Aliyev
naliye[email protected]
Volker Mehrmann
[email protected]berlin.de
Emre Mengi
1Institute of Mathematics and Mechanics, Azerbaijan National Academy of Sciences, B.Vahabzade 9,
AZ1141, Baku, Azerbaijan
2Institut fur Mathematik, Technische Universit¨
at Berlin,
MA 4-5. Str. des 17 Juni 136, 0623, Berlin, Germany
3Department of Mathematics, Koc¸ University, Rumeli Feneri Yolu 34450, Sarıyer, Istanbul, Turkey
Adv Comput Math (2020) 46:6
Page 28 of 28
6
