Approximation of stability radii for large-scale dissipative Hamiltonian systems

A linear time-invariant dissipative Hamiltonian (DH) system ẋ=(J−R)Qx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\dot x = (J-R)Q x$\end{document}, with a skew-Hermitian J, a Hermitian positive semidefinite R, and a Hermitian positive 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 perturbations of R in the form R + BΔCH for given matrices B, C, and another with respect to Hermitian perturbations in the form R + BΔBH,Δ = Δ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.


Introduction
Linear time-invariant Dissipative Hamiltonian (DH) systems are of the forṁ x = (J − R)Qx, (1.1) where Q = Q H ∈ C n×n is a Hermitian positive definite matrix denoted by Q > 0, whereas J = −J H ∈ C n×n , and the dissipation matrix R = R H ∈ C n×n is 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 engineering, 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 necessarily asymptotically stable, since A may 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 exponential growth for some initial conditions. Asymptotic stability in the presence of uncertainty can only be guaranteed when the stability radius of the system is reasonably 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 perturbations, 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.

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 (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 N of small norm, may result in instability. Here, we emphasize that the perturbations due to the matrix N are 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 uncertainties 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 −1 is sparse, whereas Q is 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 −1 rather than Q.

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 R of (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 iR refers to the imaginary axis in the complex plane, (A) to the spectrum of a matrix A, and A 2 to the spectral norm of A. Definition 1 Consider a DH system of the form (1.1) and suppose that B ∈ C n×m and C ∈ C p×n , n ≥ m, p, are given restriction matrices.
-The non-Hermitian restricted stability radius r(R; B, C) with respect to perturbations of R under the fixed restriction matrices B, C is defined by -The Hermitian restricted stability radius with respect to Hermitian perturbations of R under the fixed restriction matrix B is defined by 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 N is tolerable to preserve the asymptotic stability of the system in (1.4). The relevant stability radius for a specified is given by  (R; B). Indeed, these two stability radii can differ considerably, see also [21]. Our main concern is the approximation of these two stability radii for large-scale DH systems when n m, 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 dimension. 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 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 r Herm (R; B) poses more challenges. Unlike the characterization for r(R; B, C), the eigenvalue optimization characterization for r Herm (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 Hermite interpolation properties between the original and projected objective eigenvalue functions hold.
Outline The paper is organized as follows. Section 2 proposes subspace frameworks for approximating r(R; B, C), while Section 3 concerns r Herm (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 2 and 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].

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] for an asymptotically stable DH system of the form (1.1) and restriction matrices B ∈ C n×m , C ∈ C p×n . Above, we use the notation

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 V k and W k of equal dimension, as well as matrices V k and W k , respectively, whose columns form orthonormal bases for these subspaces. The state x of system (2.3) is restricted to V k , i.e., in (2.3) the state x is replaced by V k x k , and it is imposed that the residual of the differential part after this restriction is orthogonal to W k . This gives rise to the reduced order system , and C k := CV k . Then letting G k (s) := C k (s E k − A k ) −1 B k be the reduced order transfer function, and assuming that E k is invertible, a maximizer ω k+1 ∈ arg max ω∈R σ max (G k (iω)) can be computed efficiently and reliably if the dimensions of V k , W k are small by employing the globally convergent method in [8,9]. Once such an ω k+1 has been determined, then the subspaces V k and W k are expanded into larger subspaces V k+1 and W k+1 in such a way that the corresponding reduced transfer function G k+1 (s) satisfies the Hermite interpolation conditions Denoting the range space of a matrix by A by Ran(A), it is shown in [2] that the inclusions ensure that the conditions in (2.5) are satisfied. For instance when m = p, in practice, we set V k+1 := orth V k (iω k+1 I − A) −1 B , and W k+1 := orth W k (iω k+1 I − A) −H C H , where orth(M) denotes a matrix whose columns form an orthonormal basis for Ran(M). The definitions of the matrices V k+1 , W k+1 can be modified in a simple way when m = p so that they have equal number of columns; we refer to [2] for the details. The procedure is then repeated with the expanded subspaces V k+1 , W k+1 spanned by the columns of V k+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) for all k. The solution of the reduced problem max ω∈R σ max (G k (iω)) by the method in [8,9] requires extraction of all eigenvalues of 2d × 2d pencils, where d is the order of the system, i.e., the dimension of the subspaces V k , W k . Assuming that d is small, this is much cheaper than solving the linear systems (iω k+1 I − A) −1 B, (iω k+1 I − A) −H C H involving the full order system, which is needed to form V k+1 , W k+1 . Even though A = (J − R)Q has the DH structure, this is in general not true for the reduced system A k . In the next subsection, we modify the procedure of [2] to preserve the DH structure.

A structure-preserving subspace framework
In this subsection, we introduce an interpolating and DH structure-preserving version 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 interpolation result, see [16,Remark 3]. For this result and the subsequent arguments, we introduce the matrix-valued function depending on s ∈ C associated 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 G k (s) be the transfer function for the reduced system

7)
and W k is such that W H k V k = I , then, denoting the -th derivatives of G(s) and G k (s) at the point s with G ( ) ( s) and G ( ) k ( s), we have The characterization of r(R; B, C) involves the maximization of the largest singular 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 ∈ C and b ∈ C m , choose V k to satisfy (2.7). Then define i.e., W k V H k is an oblique projector onto Ran(QV k ). The reduced system matrices A k , B k , C k then satisfy where we define (2.12) and Moreover, E k = W H k V k = 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. (2.2). Furthermore, for a given point s ∈ C and a given tangent direction b ∈ C m , suppose that V k is a matrix with orthonormal columns such that (2.12), and B k , C k as in (2.13). Then the resulting reduced order systemẋ

Theorem 2 Consider a linear system as in (2.3) with the transfer function G(s) in
that satisfies According to Theorem 2, for a given s ∈ C, setting we obtain G( s) = G k ( s) and G ( s) = G k ( s), and thus the Hermite interpolation conditions are satisfied. This suggests the greedy subspace framework outlined in Algorithm 1 for the approximation of r(R; B, C).
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 imaginary 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 convergence 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 ω ∈ R one LU factorization of the matrix T (i ω) suffices. For large n, the computation 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].

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: 1. The relative distance between ω k and ω k−1 is less than a prescribed tolerance, i.e., |ω k − ω k−1 | < ε · |ω k |. 2. Letting f k := max ω∈R σ max (G k (iω)), two consecutive iterates f k , f k−1 are close enough in a relative sense, i.e., |f k − f k−1 | < ε · f k . 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 k max = 80 for both Algorithm 1 and the subspace framework from [2].
Both Algorithms seem to converge locally in our experiments. The initial interpolation 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.

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 B and C chosen 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, where f BB and f SF 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

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 at a given ω ∈ R in 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 M and the stiffness matrix K( ) and thus the matrix Q −1 are sparse, but Q turns out to be dense. Trying to invert Q −1 and/or solve a linear system with the coefficient matrix Q is computationally very expensive and would require full matrix storage. This difficulty can be avoided by exploiting that Hence, to compute X, Y as in (2.16), we proceed as follows.
A second observation that further speeds up the computations is the particular structure of the coefficient matrix 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, Table 3 for some values of .

Estimation of the hermitian restricted stability radius
In this section, we focus on Hermitian structured perturbations in the dissipation matrix, in particular approximation of r Herm (R; B) as in (1.6). The algorithms we propose exploit the following characterization of r Herm (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 (λ), where T (λ) is as in (2.6) and B + 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 ∈ C n×m of full column rank, let

·) denotes the smallest eigenvalue of its Hermitian matrix argument, and the inner supremum is finite if and only if H 1 (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 r Herm (R; B) in Theorem 3. Note that in this section, B is always assumed to be of full column rank.

An algorithm for small-scale problems
The eigenvalue optimization characterization of r Herm (R; B) is a min-max problem, where the inner maximization problem is concave; indeed, it can alternatively be expressed as a semidefinite program (SDP). Formally, for a given ω ∈ R, and H 0 (iω), H 1 (iω) representing the Hermitian matrices defined in Theorem 3, let us consider In the following subsections, we solve the minimization problem inside the square root in (3.3).

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ω) significantly faster than cvx. The only downside is that an interval containing the optimal t for (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.

Outer minimization problems
The minimum of η Herm (R; B, iω) with respect to ω ∈ R yields the square of the radius r Herm (R; B), and it turns out that the minimizing ω ∈ R corresponds 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ω) may even blow up at some ω.
We again resort to eigopt for the minimization of η Herm (R; B, iω). For the sake of completeness, a formal description is provided in Algorithm 2 below, where we use the abbreviations At step k of 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 .

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 H 0 (iω) and H 1 (iω), which are of small size provided that B has 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 H 0 (iω) and H 1 (iω) requires the solution of the linear system T (iω k )Z = B.

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 V k with orthonormal columns, we define W k as in (2.9), J k , R k , Q k as in (2.12), and B k as in (2.13), respectively. Again, V k , W k refer to the right, left subspaces spanned by the columns of V k , 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 for λ ∈ C associated with the reduced problems. (1.1), and a reduced modelẋ k = (J k − R k )Q k x k with coefficients as in (2.12) defined in terms of a given V k and W k as in (2.9).

Letting T (λ) and T k (λ) be as in (2.6) and (3.4), respectively, we then have T k (λ) = W H k T (λ) V k for all λ ∈ C.
Proof From the definitions of J k , R k , Q k we obtain where we have employed (2.10) in the second equality.
In the following, we develop a subspace framework in which the Hermite interpolation properties between the original problem defined in terms of T (λ) and the reduced problem defined in terms of T k (λ) hold. For this, we first show an auxiliary interpolation result between B H QT (λ) −1 B and its reduced counterpart. (1.1), and a reduced modelẋ k = (J k − R k )Q k x k with coefficients as in (2.12) defined in terms of a given V k and W k as in (2.9). Furthermore, let T (λ) and T k (λ) be as in (2.6) and (3.4), respectively. For a given λ ∈ C such that T ( λ) and T k ( λ) are invertible, the following assertions hold:  (3.5) where the first equality is due to the fact that V k V H k is the orthogonal projector onto V k , and the second follows, since W k V H k is a projector onto Ran(QV k ). We complete the proof by showing that

Lemma 2 Consider a DH model
To this end, let Z := T ( λ) −1 B, and Z k be such that V k Z k = Z. (There exists a unique Z k with this property, because Ran(Z) ⊆ V k .) Then T ( λ)Z = B implies that T ( λ)V k Z k = B, and thus W H k T ( λ)V k Z k = W H k B = B k . Hence, by Lemma 1, we see that Z k = T k ( λ) −1 B k , which in turn implies (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

It remains to show that V H k T ( λ)
To this end, we exploit Now define Z := T ( λ) −1 V k and Z k such that V k Z k = Z (once again such a Z k exists uniquely, because Ran(Z) where I k,m is the matrix consisting of the first m columns of the k × k identity matrix I k and k := dim V k ≥ m. By Lemma 1, this implies that Z k = T k ( λ) −1 I k,m , so, for the term inside the first parenthesis on the right-hand side of (3.7), we have For the term inside the second parenthesis on the right-hand side of (3.7), we make use of the following observation: where the last equality follows, since the columns of V k form an orthonormal basis for Ran(T ( λ) −1 B). Inserting these in (3.7), we obtain 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 H 0 (λ), L(λ), and H 1 (λ) defined via with L k (λ) denoting a lower triangular Cholesky factor of and H k,1 (λ) To ensure the uniqueness of L(λ) and L k (λ), we define them as the Cholesky factors of H 0 (λ) and H k,0 (λ) with real and positive entries along the diagonal. We emphasize here that the matrix functions H 0 (λ), L(λ), H 1 (λ) and their reduced counterparts are of the same size. However, the reduced matrix-valued functions are defined in terms of T k (λ), which is small if the dimensions of V k , W k are small, in contrast to T (λ) in the definitions of H 0 (λ), L(λ), H 1 (λ).
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 V k , and W k as in (2.9), we introduce η Herm k (R; B, iω) := sup t∈R λ min (H k,0 (iω) + tH k,1 (iω)).
Next, we establish that the quantity η Herm k (R; B, iω) is independent of the choice of the basis for the subspace V k . For this result, we introduce the notation η Herm V k ,W k (R; B, iω)) to emphasize the particular choices of the bases V k , W k used in the definition of η Herm k (R; B, iω). Similarly, we indicate the choices of the bases used in the definitions of J k , R k , Q k , and B k with the help of the notations J W k , R W k , Q V k , and B W k .

Lemma 3 Let the columns of V k and V k form orthonormal bases for the subspace
for all ω ∈ R.
Proof The proof makes use of the characterization [21,Theorem 4.9]. Now since the columns of V k and V k form bases for the same space, there exists a unitary matrix P such that V k = V k P . Furthermore, The assertion then follows from (3.13) and the following set of equivalences: Above, the first and third equivalences are due to the identities in (2.10), whereas the second equivalence is due to V k = V k P and W k = W k P .
We now prove our main interpolation result. (1.1), and a reduced modelẋ k = (J k − R k )Q k x k with coefficients as in (2.12) defined in terms of a given V k and W k as in (2.9). Furthermore, let T (λ) and T k (λ) be as in (2.6) and (3.4), respectively. For a given ω ∈ R, suppose that the subspace V k := Ran(V k ) is such that

Theorem 4 (Hermite Interpolation of Hermitian Backward Errors) Consider a DH model
(3.14) Proof Without loss of generality, we may assume that the matrix V k is such that V k = V k V k , with the columns of V k forming an orthonormal basis for Ran(T (i ω) −1 B). By Lemma 3, it suffices to prove the claims for such a particular choice of orthonormal basis.
(i) By the definitions of H 0 (λ), H k,0 (λ) in (3.1), (3.11), respectively, and by part (i) of Lemma 2, we have . This also implies that L(i ω) = L k (i ω) due to the uniqueness of the Cholesky factors of H 0 (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 analytical expressions. It follows from the standard perturbation theory of simple eigenvalues (e.g., see [20] (3.17) where t is such that t ∈ arg max t λ min (H 0 (i ω) + tH 1 (i ω)) = arg max t λ min (H k,0 (i ω) + tH k, 1

Numerical experiments
In this section, we present numerical results for the estimation of the Hermitian stability 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 f k := min ω∈R η Herm k (R; B, iω), and we set k max = 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.

Synthetic examples
We first present a numerical result for a small dense random example, where J, R, Q ∈ R 20×20 and B ∈ R 20×2 . Application of Algorithm 2 to this example yields r Herm (R; B) = 0.0501, and the point that is first reached on the imaginary axis under the smallest perturbation is 1.9794i. Figure 2 illustrates that there exist real symmetric matrices with 2 = 0.0501 such that 1.9794i is contained (nearly) in the spectrum of (J − (R + B B H ))Q.
Next we experiment with larger dense random matrices, specifically with random A ∈ R n×n , B ∈ R n×2 . The results of Algorithm 3 are reported in Tables 4 and 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 r Herm (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 r Herm (R; B) considerably.  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 3 depicts η Herm (R; B, iω) as a function of ω with the solid curve for the first example in Table 4 with 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.

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 r Herm (R; B) along with the associated minimizer ω * of η Herm (R; B, iω), as well as the run-time (in seconds) are listed in Table 6 for 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 r Herm (R; B) in this table with those for the non-Hermitian stability radius r(R; B, B T ), 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.

Concluding remarks
We have proposed subspace frameworks to estimate the stability radii for largescale dissipative Hamiltonian systems. At every iteration, we apply DH structurepreserving 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.
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 .