Computing nearest stable matrix pairs

TL;DR

This paper introduces a new approach to find the nearest stable matrix pair by reformulating the problem with dissipative Hamiltonian matrix pairs, enabling efficient computation of stable approximations.

Contribution

The paper proposes a reformulation using dissipative Hamiltonian matrix pairs, simplifying the feasible set and allowing fast gradient methods for stability approximation.

Findings

Reformulation with DH matrix pairs simplifies the feasible set.

The method enables fast gradient-based computation of stable matrix pairs.

The approach provides a practical solution for stability approximation in matrix pairs.

Abstract

In this paper, we study the nearest stable matrix pair problem: given a square matrix pair $(E,A)$, minimize the Frobenius norm of $(\Delta_E,\Delta_A)$ such that $(E+\Delta_E,A+\Delta_A)$ is a stable matrix pair. We propose a reformulation of the problem with a simpler feasible set by introducing dissipative Hamiltonian (DH) matrix pairs: A matrix pair $(E,A)$ is DH if $A=(J-R)Q$ with skew-symmetric $J$, positive semidefinite $R$, and an invertible $Q$ such that $Q^TE$ is positive semidefinite. This reformulation has a convex feasible domain onto which it is easy to project. This allows us to employ a fast gradient method to obtain a nearby stable approximation of a given matrix pair.

Figures

Equations

$$E\dot{x}=Ax+f,$$

$$F(t,x,\dot{x})=0,$$

$$x(t_0)=x_0.$$

$$E=W\left[\begin{array}[]{cc}I_{q}&0\\ 0&N\end{array}\right]T\quad\text{and}\quad A=W\left[\begin{array}[]{cc}J&0\\ 0&I_{n-q}\end{array}\right]T,$$

$$E\dot{x}=(A+BK)x$$

$$\underset{(M,X)\in \mathbb{S}}{\inf }\{{\Vert E-M\Vert }_F^2+{\Vert A-X\Vert }_F^2\},$$

$$(E,A)=\left(\left[\begin{array}[]{ccc}1&0&0\\ 0&0&0\\ 0&0&0\end{array}\right],\ \left[\begin{array}[]{ccc}-1&0&2\\ 0&1&0\\ 0&0&1\end{array}\right]\right).$$

$$\Delta_{E}=\left[\begin{array}[]{ccc}\epsilon_{1}&0&0\\ 0&\epsilon_{2}&\epsilon_{3}\\ 0&0&0\end{array}\right]\quad\text{and}\quad\Delta_{A}=\left[\begin{array}[]{ccc}\delta&0&0\\ 0&0&0\\ 0&0&0\end{array}\right],$$

$$E=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right],J=\left[\begin{array}{ccc}0& 2& 0\\ -2& 0& 0\\ 0& 0& 0\end{array}\right],R=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right],Q=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],$$

$$(J-R)Qx=\lambda Ex.$$

$$x^H{Q}^{T}JQx-x^H{Q}^{T}RQx=\lambda x^H{Q}^{T}Ex,$$

$$\mathrm{Re}(\lambda )=-\frac{x^H{Q}^{T}RQx}{x^{\ast }{Q}^{T}Ex}\le 0,$$

$$(J-R)Qx=\lambda Ex.$$

$$x^H{Q}^T(J-R)Qx=\lambda x^H{Q}^{T}Ex.$$

$$\mathrm{Re}(\lambda )=-\frac{x^H{Q}^{T}RQx}{x^H{Q}^{T}Ex}\le 0,$$

$$(\lambda E-(J-R)Q)x_0=0,(\lambda E-(J-R)Q)x_1$$

$$(\lambda E-(J-R)Q)x_2$$

$$⋮$$

$$(\lambda E-(J-R)Q)x_{k-1}$$

$$(\lambda E-JQ)x_0=0\text{mbox}\mathrm{and}RQ{x}_0=0.$$

$$(\lambda E-(J-R)Q)x_1=E{x}_0.$$

$$x_0^H(\lambda Q^{T}E-(Q^{T}JQ-Q^{T}RQ))x_1=x_0^H{Q}^{T}E{x}_0.$$

$$-x_1^H(\lambda Q^{T}E-Q^{T}JQ)x_0+x_1^H{Q}^{T}RQ{x}_0=x_0^H{Q}^{T}E{x}_0,$$

$$x^{H}E=0,x^H(J-R)QU=0.$$

$$Q=V,J=\frac{(A{V}^{-1})-(A{V}^{-1}{)}^T}{2},\text{and}R=-\frac{(A{V}^{-1})+(A{V}^{-1}{)}^T}{2},$$

$$V^{T}RV=-\frac{V^T((A{V}^{-1})+(A{V}^{-1}{)}^T)V}{2}=-\frac{V^{T}A+A^{T}V}{2}\succ 0,$$

$$(\tilde{E},\tilde{A})$$

$$\mathbb{S}=\left\{(E,(J-R)Q)\in {\mathbb{R}}^{n,n}\times {\mathbb{R}}^{n,n}:J^T=-J,R\succ 0,Q\text{invertible s.t.}{Q}^{T}E⪰0\right\}=:{\mathbb{S}}_{DH}^{\succ 0}.$$

$$\underset{(M,X)\in \mathbb{S}}{\inf }\{{\Vert E-M\Vert }_F^2+{\Vert A-X\Vert }_F^2\}=\underset{(M,(J-R)Q)\in {\mathbb{S}}_{DH}^{\succ 0}}{\inf }\{{\Vert E-M\Vert }_F^2+{\Vert A-(J-R)Q\Vert }_F^2\}.$$

$${\mathbb{S}}_{DH}^{⪰0}:=\left\{(E,(J-R)Q)\in {\mathbb{R}}^{n,n}\times {\mathbb{R}}^{n,n}:J^T=-J,R⪰0,Q^{T}E⪰0\right\}.$$

$$\underset{(M,(J-R)Q)\in {\mathbb{S}}_{DH}^{\succ 0}}{\inf }\{{\Vert E-M\Vert }_F^2+{\Vert A-(J-R)Q\Vert }_F^2\}=\underset{(M,(J-R)Q)\in {\mathbb{S}}_{DH}^{⪰0}}{\inf }\{{\Vert E-M\Vert }_F^2+{\Vert A-(J-R)Q\Vert }_F^2\}.$$

$$\underset{(M,X)\in \mathbb{S}}{\inf }\{{\Vert E-M\Vert }_F^2+{\Vert A-X\Vert }_F^2\}=\underset{(M,(J-R)Q)\in {\mathbb{S}}_{DH}^{⪰0}}{\inf }\{{\Vert E-M\Vert }_F^2+{\Vert A-(J-R)Q\Vert }_F^2\}.$$

Full text

Keywords. dissipative Hamiltonian system, distance to stability, convex optimization,

1 Introduction

We study the stability and instability under perturbations for systems of linear time-invariant differential-algebraic equations (DAEs) of the form

[TABLE]

on the unbounded interval ${\mathbb{I}}=[t_{0},\infty)$, where $E,A\in\mathbb{R}^{n,n}$ and $f\in C^{0}({\mathbb{I}},\mathbb{R}^{n})$ is sufficiently smooth. Here $\mathbb{R}^{n,n}$ denotes the real $n\times n$ matrices and $C^{k}(D,\mathbb{R}^{n})$ denotes the $k$-times differentiable functions from a set $D$ to $\mathbb{R}^{n}$. Systems of the form (1.1) arise from linearization around stationary solutions of initial value problems for general implicit systems of DAEs

[TABLE]

with an initial condition

[TABLE]

Here we assume that $F\in C^{0}({\mathbb{I}}\times{\mathbb{D}}_{x}\times{\mathbb{D}}_{\dot{x}},{\mathbb{R}}^{n})$ is sufficiently smooth and ${\mathbb{D}}_{x},{\mathbb{D}}_{\dot{x}}\subseteq{\mathbb{R}}^{n}$ are open sets, and following [16], we use the following solution and stability concepts.

Definition 1**.**

Consider system (1.2) with sufficiently smooth $F$.

• i)

A function $x:{\mathbb{I}}\rightarrow\mathbb{R}^{n}$ is called a solution of (1.2) if $x\in C^{1}({\mathbb{I}},\mathbb{R}^{n})$ and $x$ satisfies (1.2) pointwise.

• ii)

It is called a solution of the initial value problem (1.2)–(1.3) if $x$ is a solution of (1.2) and satisfies (1.3).

• iii)

An initial condition (1.3) is called consistent if the corresponding initial value problem has at least one solution.

• iv)

System (1.2) is called regular if for every consistent initial condition $x_{0}$ it has a unique solution.

Definition 2**.**

Suppose that the DAE (1.2) is regular and (1.3) is a consistent initial condition. A solution $x:t\mapsto x(t;t_{0},x_{0})$ of (1.2)–(1.3) is called

1. i)

stable* if for every $\varepsilon>0$ there exists $\delta>0$ such that*

1. (a)

the initial value problem (1.2) with initial condition $x(t_{0})=\hat{x}_{0}$ is solvable on ${\mathbb{I}}$ for all consistent $\hat{x}_{0}\in{\mathbb{R}}^{n}$ with $\|\hat{x}_{0}-x_{0}\|<\delta$; 2. (b)

the solution $x(t;t_{0},\hat{x}_{0})$ satisfies $\|x(t;t_{0},\hat{x}_{0})-x(t;t_{0},x_{0})\|<\varepsilon$ on ${\mathbb{I}}$. 2. ii)

asymptotically stable* if it is stable and there exists $\varrho>0$ such that*

1. (a)

the initial value problem (1.2) with initial condition $x(t_{0})=\hat{x}_{0}$ is solvable on ${\mathbb{I}}$ for all consistent $\hat{x}_{0}\in{\mathbb{R}}^{n}$ with $\|\hat{x}_{0}-x_{0}\|<\varrho$; 2. (b)

the solution $x(t;t_{0},\hat{x}_{0})$ satisfies $\lim_{t\to\infty}\|x(t;t_{0},\hat{x}_{0})-x(t;t_{0},x_{0})\|=0$.

To characterize regularity and stability for the linear constant coefficient case we introduce the following notation. A square matrix pair $(E,A)$ with $E,A\in\mathbb{R}^{n,n}$ is called regular if the matrix pencil $zE-A$ is regular, i. e. if $\operatorname{det}(\lambda E-A)\neq 0$ for some $\lambda\in\mathbb{C}$, otherwise it is called singular. For a regular matrix pair $(E,A)$, the roots of the polynomial $\operatorname{det}(zE-A)$ are called finite eigenvalues of the pencil $zE-A$ or of the pair $(E,A)$, i. e. $\lambda\in\mathbb{C}$ is a finite eigenvalue of the pencil $zE-A$ if there exists a vector $x\in\mathbb{C}^{n}\setminus\{0\}$ such that $(\lambda E-A)x=0$, and $x$ is called an eigenvector of $zE-A$ corresponding to the eigenvalue $\lambda$. A regular pencil $zE-A$ has $\infty$ as an eigenvalue if $E$ is singular.

A regular real matrix pair $(E,A)$ (with $E,A\in{\mathbb{R}}^{n,n}$) can be transformed to Weierstraß canonical form [11], i. e. there exist nonsingular matrices $W,T\in{\mathbb{C}}^{n,n}$ such that

[TABLE]

where $J\in{\mathbb{C}}^{q,q}$ is a matrix in Jordan canonical form associated with the $q$ finite eigenvalues of the pencil $zE-A$ and $N\in{\mathbb{C}}^{n-q,n-q}$ is a nilpotent matrix in Jordan canonical form corresponding to $n-q$ times the eigenvalue $\infty$. If $q<n$ and $N$ has degree of nilpotency $\nu\in\{1,2,\ldots\}$, i. e. $N^{\nu}=0$ and $N^{i}\neq 0$ for $i=1,\ldots,\nu-1$, then $\nu$ is called the index of the pair $(E,A)$. If $E$ is nonsingular, then by convention the index is $\nu=0$. A pencil $zE-A$ is of index at most one if it is regular with exactly $r:=\text{rank}(E)$ finite eigenvalues, see e. g. [23, 28]. In this case the $n-r$ copies of the eigenvalue $\infty$ are semisimple (non-defective).

Then we have the following well-known results, see e. g. [8, 10, 17].

Proposition 1**.**

Consider the initial value problem (1.1).

• i)

The linear DAE (1.1) is regular if and only if the matrix pair $(E,A)$ is regular **[16]**.

• ii)

If the pair $(E,A)$ is regular and of index at most one, then the initial value problem (1.1), (1.3) with a consistent initial value $x_{0}$ is stable if all the finite eigenvalues of $zE-A$ are in the closed left half of the complex plane and those on the imaginary axis are semisimple.

• iii)

If the pair $(E,A)$ is regular and of index at most one then the initial value problem (1.1), (1.3) with a consistent initial value $x_{0}$ is asymptotically stable if all the finite eigenvalues of $zE-A$ are in the open left half of the complex plane.

The literature on (asymptotic) stability of constant coefficient DAEs is very ambiguous, see e. g. [4, 8, 28], and the review in [10]. This ambiguity arises from the fact that some authors consider only the finite eigenvalues in the stability analysis and allow the index of the pencil $(E,A)$ to be arbitrary, others consider regular high index pencils $zE-A$ as unstable by considering $\infty$ to be on the imaginary axis.

In this paper, we are interested in the problem of finding minimal perturbations to the system matrices $E+\Delta E$ and $A+\Delta A$ of an unstable DAE (1.1) that puts the system on the boundary of the stability region. This problem is closely related to the classical state-feedback-stabilization problem in descriptor control problems with $f=Bu$ for some control function $u$, see e. g. [3, 9, 23], where one chooses feedbacks $u=Kx$ so that the closed-loop descriptor system

[TABLE]

is asymptotically stable. However, then the perturbations are restricted only in $A$ while in derivative feedback one also allows feedbacks $E+BG$, see [5, 6]. Instead of the feedback stabilization problem, we extend the nearest stable matrix problem of [12, 26] to matrix pencils and we study the following problem.

Problem**.**

For a given pair $(E,A)\in{\mathbb{R}}^{n,n}\times{\mathbb{R}}^{n,n}$ find the nearest asymptotically stable matrix pair $(M,X)$. More precisely, let $\mathbb{S}$ be the set of matrix pairs $(M,X)\in{\mathbb{R}}^{n,n}\times{\mathbb{R}}^{n,n}$ that are regular, of index at most one, and have all finite eigenvalues in the open left half plane, then we wish to compute

[TABLE]

where ${\|\cdot\|}_{F}$ stands for the Frobenius norm of a matrix. **

In the following we use nearest stable matrix pair problem to refer to the above problem. This problem is the complementary problem to the distance to instability for matrix pairs, see [8] for complex pairs, and [10] for a survey on this problem. Since we require a stable pair to be regular it is also complementary to the distance to the nearest singular pencil, which is a long-standing open problem [7, 14, 22].

For ordinary differential equations, i. e. when $E=I_{n}$, stability solely depends on the eigenvalues of the matrix $A$. Recently, in [12] a new algorithm was proposed to compute a nearby stable perturbation to a given unstable matrix by reformulating this highly nonconvex problem into an equivalent optimization problem with a convex feasible region. This reformulation uses the fact that the set of stable matrices can be characterized as the set of dissipative Hamiltonian (DH) matrices. We recall from [12] that a matrix $A\in{\mathbb{R}}^{n,n}$ is called a DH matrix if $A$ can be factored as $A=(J-R)Q$ for some $J,R,Q\in{\mathbb{R}}^{n,n}$ such that $J^{T}=-J$, $R=R\succeq 0$, i. e. is positive semidefinite, and $Q=Q^{T}\succ 0$, i. e. is positive definite. In this way $A$ is stable if and only if $A$ is a DH matrix, see also [1, 20, 21].

In this paper we extend the ideas of [12] from the matrix case to the pencil case and discuss the problem of computing the nearest stable matrix pair. This is a challenging problem, since the set $\mathbb{S}$ of all stable matrix pairs is neither open nor closed with respect to the usual norm topology as the following example demonstrates.

Example 1**.**

Consider the matrix pair

[TABLE]

It is easy to check that $(E,A)$ is regular, of index one, and asymptotically stable with the only finite eigenvalue $\lambda_{1}=-1$, and thus $(E,A)\in\mathbb{S}$. Consider the perturbations $(\Delta_{E},\Delta_{A})$, where

[TABLE]

and the perturbed pair $(E+\Delta_{E},A+\Delta_{A})$. If we let $\delta=\epsilon_{1}=\epsilon_{2}=0$ and $\epsilon_{3}>0$, then the perturbed pair is still regular and all finite eigenvalues are in open left half plane, but it is of index two. If we let $\delta=\epsilon_{1}=\epsilon_{3}=0$ and $\epsilon_{2}>0$, then the perturbed pair is regular and index one, but has two finite eigenvalues $\lambda_{1}=-1$ and $\lambda_{2}=1/\epsilon_{2}>0$, and thus is not stable. This shows that $\mathbb{S}$ is not open. If we let $\epsilon_{1}=-\delta$, $\epsilon_{2}=\epsilon_{3}=0$, then as $\delta\rightarrow 1$ the perturbed pair becomes non-regular. This shows that $\mathbb{S}$ is not closed. **

It is also clear that the set $\mathbb{S}$ is highly nonconvex, see [26], and thus it is very difficult to find a globally optimal solution of the nearest stable matrix pair problem. To address this challenge, we follow the strategy suggested in [12] for matrices and reformulate the problem of computing the nearest stable matrix pair by extending the concept of dissipative Hamiltonian matrices to dissipative Hamiltonian matrix pairs (DH pairs).

The paper is organized as follows. In Section 2, we introduce and study DH matrix pairs. We provide several theoretical results and characterize the set of asymptotically stable matrix pairs in terms of DH matrix pairs. A reformulation of the nearest stable matrix pair problem using the DH characterization is obtained in Section 3. In Section 4, we propose a fast gradient method to solve the reformulated optimization problem. The effectiveness of the proposed algorithm is illustrated by several numerical examples in Section 5.

In the following, by $i{\mathbb{R}}$, we denote the imaginary axis of the complex plane, by $I_{n}$ the identity matrix of size $n\times n$, and by $\operatorname{null}(E)$ the nullspace of a matrix $E$.

2 Dissipative Hamiltonian matrix pairs

In this section, we construct the setup to approach the nearest stable matrix pair problem. As demonstrated in Example 1, the feasible set $\mathbb{S}$ is not open, not closed, non-bounded, and highly nonconvex, thus it is very difficult to work directly with the set $\mathbb{S}$. For this reason we reformulate the nearest stable matrix pair problem into an equivalent optimization problem with a simpler feasible set. For this, we extend the idea of a DH matrix from [12] to dissipative Hamiltonian matrix pairs.

Definition 3**.**

A matrix pair $(E,A)$, with $E,A\in\mathbb{R}^{n,n}$, is called a dissipative Hamiltonian (DH) matrix pair if there exists an invertible matrix $Q\in\mathbb{R}^{n,n}$ such that $Q^{T}E=E^{T}Q\succeq 0$, and $A$ can be expressed as $A=(J-R)Q$ with $J^{T}=-J$, $R^{T}=R\succeq 0$.

This definition of a DH matrix pair is a natural generalization of that of a DH matrix, because for a standard DH matrix pair $(I_{n},A)$ (when $E=I_{n}$), $A$ is a DH matrix [12]. Note however, that this definition is slightly more restrictive than that of port-Hamiltonian descriptor systems in [2], where it is not required that the matrix $Q$ is invertible. In the following we call the matrix $R$ in a DH matrix pair $(E,(J-R)Q)$ the dissipation matrix as in the DH matrix case.

A DH-matrix pair is not necessarily regular as the following example shows.

Example 2**.**

The pair $(E,(J-R)Q)$ with

[TABLE]

is a DH matrix pair, but $(E,(J-R)Q)$ is singular, since $\operatorname{det}(zE-(J-R)Q)\equiv 0$. **

It is easy to see that if the matrices $E$ and $A$ have a common nullspace, then the pair $(E,A)$ is singular, but the converse is in general not true, see e. g. [7]. However, for a singular DH matrix pairs, the property that $Q$ is invertible guarantees that the converse also holds.

Lemma 1**.**

Let $(E,(J-R)Q)$ be a DH-matrix pair. Then $(E,(J-R)Q)$ is singular if and only if $\operatorname{null}(E)\cap\operatorname{null}((J-R)Q)\neq\emptyset$.

Proof.

The direction $(\Leftarrow)$ is immediate. For the other direction, let $(E,(J-R)Q)$ be singular, i. e. $\operatorname{det}(zE-(J-R)Q)\equiv 0$. Let $\lambda\in\mathbb{C}$ be such that $\mathop{\mathrm{Re}}{(\lambda)}>0$ and let $x\in\mathbb{C}^{n}\setminus\{0\}$ be such that

[TABLE]

Since $Q$ is nonsingular, we have that $Qx\neq 0$, and we can multiply with $(Qx)^{H}$ from the left to obtain

[TABLE]

where $x^{H}$ denotes the complex conjugate of a vector $x$. This implies that $x^{H}Q^{T}Ex=0$, because otherwise from (2.2) we would have

[TABLE]

since $Q^{T}RQ\succeq 0$ (as $R\succeq 0$), and $Q^{T}E\succeq 0$. But this is a contradiction to the fact that $\mathop{\mathrm{Re}}{(\lambda)}>0$. Therefore $x^{H}Q^{T}Ex=0$ and also $Q^{T}Ex=0$, and this implies that $Ex=0$ as $Q$ is invertible. Inserting this in (2.1), we get $(J-R)Qx=0$, i. e. $0\neq x\in\operatorname{null}(E)\cap\operatorname{null}((J-R)Q)$.

Lemma 1 gives a necessary and sufficient condition for a DH-matrix pair to be singular. However, if the dissipation matrix $R$ is positive definite, then regularity is assured, as shown in the following corollary.

Corollary 1**.**

Let $(E,(J-R)Q)$ be a DH matrix pair. If $R$ is positive definite then the pair is regular.

Proof.

By Lemma 1, a necessary condition for the pencil $\lambda E-(J-R)Q$ to be singular is that neither $E$ nor $(J-R)Q$ is invertible. Thus the result follows immediately by the fact that if $R\succ 0$ in a DH-matrix pair $(E,(J-R)Q)$, then $(J-R)Q$ is invertible. Indeed, suppose there exists $x\in\mathbb{C}\setminus\{0\}$ such that $(J-R)Qx=0$, then we have $x^{H}Q^{T}(J-R)Qx=0$. This implies $x^{H}Q^{T}RQx=0$, since $J^{T}=-J$. Since $Q$ is invertible and thus $Qx\neq 0$, this is a contradiction to the assumption that $R$ is positive definite.

In the following lemma, which is the matrix pair analogue of [12, Lemma 2], we localize the finite eigenvalues of a DH-matrix pair.

Lemma 2**.**

Let $(E,(J-R)Q)$ be a regular DH-matrix pair and let $L(z):=zE-(J-R)Q$. Then the following statements hold.

All finite eigenvalues of the pencil $L(z)$ are in the closed left half of the complex plane. 2. 2)

The pencil $L(z)$ has a finite eigenvalue $\lambda$ on the imaginary axis if and only if $RQx=0$ for some eigenvector $x$ of the pencil $zE-JQ$ associated with $\lambda$.

Proof.

Let $\lambda\in\mathbb{C}$ be an eigenvalue of $zE-(J-R)Q$ and let $x\in\mathbb{C}^{n}\setminus\{0\}$ be such that

[TABLE]

Multiplying (2.3) by $x^{H}Q^{T}$ from the left, we get

[TABLE]

Note that $x^{H}Q^{T}Ex\neq 0$, because if $x^{H}Q^{T}Ex=0$, then we have $Q^{T}Ex=0$ as $Q^{T}E\succeq 0$, and thus $Ex=0$ as $Q$ is invertible. Using this in (2.3), we have $(J-R)Qx=0$. This implies that $x\in\operatorname{null}(E)\cap\operatorname{null}((J-R)Q)$ which is, by Lemma 1, a contradiction to the regularity of the pair $(E,(J-R)Q)$.

Thus, by (2.4) we have

[TABLE]

because $Q^{T}RQ\succeq 0$ as $R\succeq 0$, and $x^{H}Q^{T}Ex>0$. This completes the proof of 1).

In the proof of 1), if $\lambda\in i\mathbb{R}$, then from (2.5) it follows that $x^{H}Q^{T}RQx=0$. This implies that $RQx=0$, since $R\succeq 0$. Using this in (2.3) implies that $(\lambda E-JQ)x=0$.

Conversely, let $\lambda\in i\mathbb{R}$ and $x\in\mathbb{C}^{n}\setminus\{0\}$ be such that $RQx=0$ and $(\lambda E-JQ)x=0$. Then this trivially implies that $\lambda$ is also an eigenvalue of the pencil $\lambda E-(J-R)Q$ with eigenvector $x$. This completes the proof of 2).

Making use of these preliminary results, we have the following stability characterization.

Theorem 1**.**

Every regular DH matrix pair $(E,(J-R)Q)$ of index at most one is stable.

Proof.

In view of Lemma 2, to prove the result it is sufficient to show that if $\lambda\in i\mathbb{R}$ is an eigenvalue of the pencil $zE-(J-R)Q$, then $\lambda$ is semisimple.

Let us suppose that $\lambda\in i\mathbb{R}$ is a defective eigenvalue of the pencil $zE-(J-R)Q$ and the set $\{x_{0},x_{1},\ldots,x_{k-1}\}$ forms a Jordan chain of length $k$ associated with $\lambda$, see e. g. [13], i. e. $x_{0}\neq 0$ and

[TABLE]

Note that by Lemma 2, we have that $(\lambda E-(J-R)Q)x_{0}=0$ implies that

[TABLE]

By (2.6), $x_{0}$ and $x_{1}$ satisfy

[TABLE]

Multiplying (2.8) by $x_{0}^{H}Q^{T}$ from the left, we obtain

[TABLE]

This implies that

[TABLE]

where the last equality follows by the fact that $Q^{T}E=E^{T}Q$, $J^{T}=-J$, and $R^{T}=R$. Thus, by using (2.7) in (2.9), we get $x_{0}^{H}Q^{T}Ex_{0}=0$. But this implies that $Q^{T}Ex_{0}=0$ as $Q^{T}E\succeq 0$ and $Ex_{0}=0$ as $Q$ is invertible. Since $x_{0}$ is an eigenvector of the pencil $zE-(J-R)Q$ to $\lambda$, $Ex_{0}=0$ implies that $(J-R)Qx_{0}=0$. This means that $0\neq x_{0}\in\operatorname{null}(E)\cap\operatorname{null}((J-R)Q)$, which contradicts the regularity of the pair $(E,(J-R)Q)$. Therefore there does not exist a vector $x_{1}\in\mathbb{C}^{n}$ satisfying (2.6). Hence $\lambda$ is semisimple.

We note that the proofs of Lemma 2 and Theorem 1 highly depend on the invertibility of $Q$. In fact, we have used the fact that $(Q^{T}E,Q^{T}(J-R)Q)$ is regular, because $Q$ invertible implies that $(E,(J-R)Q)$ is regular if and only if $(Q^{T}E,Q^{T}(J-R)Q)$ is regular. If $Q$ is singular and $Q^{T}E=E^{T}Q\succeq 0$, then the pair $(Q^{T}E,Q^{T}(J-R)Q)$ is always singular, but the pair $(E,(J-R)Q)$ may be regular.

However, parts of Lemma 2 and Theorem 1 also hold for singular $Q$ with an extra assumption on the eigenvectors of the pencil $zE-(J-R)Q$, as stated in the next result. This result is a generalization of [12, Lemma 2] for DH matrices with singular $Q$.

Theorem 2**.**

Let $(E,(J-R)Q)$ be a regular matrix pair with $J^{T}=-J$, $R\succeq 0$, and a singular $Q$ such that $E^{T}Q\succeq 0$. If no eigenvector $x$ of $L(z):=zE-(J-R)Q$ with respect to a nonzero eigenvalue satisfies $Q^{T}Ex=0$, then the following hold.

All finite eigenvalues of the pencil $L(z)$ are in the left half of the complex plane. 2. 2)

$L(z)$* has a finite eigenvalue $\lambda$ on the imaginary axis if and only if $RQx=0$ for some eigenvector $x$ of the pencil $zE-JQ$ with respect to eigenvalue $\lambda$.* 3. 3)

If $0\neq\lambda\in i\mathbb{R}$ is an eigenvalue of $zE-(J-R)Q$, then $\lambda$ is semisimple.

Note that if $E=I_{n}$ in Theorem 2, then $Q$ is positive semidefinite and singular, and the extra condition on the eigenvectors of pencil $zI_{n}-(J-R)Q$ (or, equivalently of matrix $(J-R)Q$) is trivially satisfied. In this case Theorem 2 coincides with the matrix case [12, Lemma 2].

The following result from [19] gives an equivalent condition for checking asymptotical stability of such a pair $(E,A)$ by using Lyapunov’s Theorem.

Theorem 3** ([19]).**

Consider a pair $(E,A)$ with $E,A\in\mathbb{R}^{n,n}$. The pair is regular, of index at most one and asymptotically stable if and only if there exists a nonsingular $V\in\mathbb{R}^{n,n}$ satisfying $V^{T}A+A^{T}V\prec 0$ and $E^{T}V=V^{T}E\succeq 0$.

Theorem 4**.**

Let $(E,A)$ be a matrix pair, where $E,A\in{\mathbb{R}}^{n,n}$. Then the following statements are equivalent.

$(E,A)$* is a DH-matrix pair with positive definite dissipation matrix.* 2. 2)

$(E,A)$* is regular, of index at most one, and asymptotically stable.*

Proof.

$1)\Rightarrow 2)$ Let $(E,A)$ be a DH-matrix pair with positive definite dissipation matrix, i. e., $A=(J-R)Q$ for some $R\succ 0$, $J^{T}=-J$, and nonsingular $Q$ with $Q^{T}E\succeq 0$. Clearly, by Corollary 1 $(E,(J-R)Q)$ is regular. Furthermore, $(E,(J-R)Q)$ has all its finite eigenvalues in the open left half plane. To see this, let $\lambda\in\mathbb{C}$ be a finite eigenvalue of the pencil $zE-(J-R)Q$. Then by Lemma 2 it follows that $\mathop{\mathrm{Re}}{(\lambda)}\leq 0$, and $\mathop{\mathrm{Re}}{(\lambda)}=0$ if and only if there exists $x\neq 0$ such that $(\lambda E-JQ)x=0$ and $0\neq Qx\in\operatorname{null}(R)$. But $\operatorname{null}(R)=\{0\}$ as $R\succ 0$.

To show that $(E,(J-R)Q)$ is of index at most one, we set $r:=\text{rank}(E)$ and assume that $U\in{\mathbb{R}}^{n,n-r}$ is an orthogonal matrix whose column spans $\text{null}(E)$. Then, see [15], $(E,(J-R)Q)$ is of index at most one if and only if \text{rank}(\left[\begin{array}[]{cc}E&(J-R)QU\end{array}\right])=n. Suppose that $x\in{\mathbb{C}}^{n}\in\setminus\{0\}$ is such that x^{H}\left[\begin{array}[]{cc}E&(J-R)QU\end{array}\right]=0. Then we have the two conditions

[TABLE]

Since $Q$ is invertible, we have $x^{H}EQ^{-1}=0$ and hence $(EQ^{-1})x=0$ because $EQ^{-1}\succeq 0$ as $E^{T}Q\succeq 0$. This shows that $Q^{-1}x\in\text{null}(E)$, and thus there exists $y\in{\mathbb{C}}^{n-r}$ such that $Q^{-1}x=Uy$, or, equivalently $x=QUy$. Using this in (2.10), we obtain that $x^{H}(J-R)x=0$. This implies that $x^{H}Jx=0$ and $x^{H}Rx=0$ as $J$ is skew-symmetric and $R$ is symmetric. But this is a contradiction to the assumption that $R\succ 0$. This completes the proof of $1)\Rightarrow 2)$.

$2)\Rightarrow 1)$ Consider a pair $(E,A)$, with $E,A\in\mathbb{R}^{n,n}$, that is regular, asymptotically stable, and of index at most one. Then by Theorem 3, there exist an nonsingular $V\in\mathbb{R}^{n,n}$ such that $V^{T}A+A^{T}V\prec 0$ and $\quad E^{T}V=V^{T}E\succeq 0$. Setting

[TABLE]

we have $J^{T}=-J$, $E^{T}Q=Q^{T}E\succeq 0$, and $R\succ 0$, as $V$ is invertible. Applying the Lyapunov inequality

[TABLE]

the assertion follows.

We conclude this section with the observation that the set of DH matrix pairs is invariant under orthogonal transformations of the matrix pair.

Lemma 3**.**

Let $(E,A)$ be a DH matrix pair, and let $U$ and $V$ be orthogonal matrices such that $\widetilde{E}=UEV$ and $\widetilde{A}=UAV$. Then $(\widetilde{E},\widetilde{A})$ is a DH pair.

Proof.

Since $(E,A)$ is DH, we have $A=(J-R)Q$ for some $J^{T}=-J$, $R\succeq 0$, and invertible $Q$ such that $Q^{T}E\succeq 0$. Using the orthogonality of $U$ and $V$, we have

[TABLE]

where $\widetilde{J}:=UJU^{T}$, $\widetilde{R}:=URU^{T}$, and $\widetilde{Q}:=UQV$. Since $U$ and $V$ are orthogonal, we have that $\widetilde{R}\succeq 0$ as $R\succeq 0$, and $\widetilde{Q}$ is invertible with $\widetilde{Q}^{T}\widetilde{E}\succeq 0$ as $Q$ is invertible with $Q^{T}E\succeq 0$. This shows that $(\widetilde{E},\widetilde{A})$ is again DH.

3 Reformulation of the nearest stable matrix pair problem

In this section, we exploit the results obtained in the previous section and derive a reformulation of the nearest stable matrix pair problem. By Theorem 4, the set $\mathbb{S}$ of all asymptotically stable matrix pairs can be expressed as the set of all DH matrix pairs with positive definite dissipation, i. e.,

[TABLE]

This characterization changes the feasible set and also the objective function in the nearest stable matrix pair problem as

[TABLE]

We have demonstrated in Example 1 that the set $\mathbb{S}$ of all stable matrix pairs is neither open nor closed and clearly the alternative characterization of $\mathbb{S}$ in terms of $\mathbb{S}_{DH}^{\succ 0}$ does not change this, since $\mathbb{S}_{DH}^{\succ 0}$ is not closed due to the constraint that $R\succ 0$ and $Q$ is invertible, and not open due to the constraint $Q^{T}E\succeq 0$. However, we can instead look at the set $\mathbb{S}_{DH}^{\succeq 0}$ containing all pairs of the form $(E,(J-R)Q)$ with $J^{T}=-J$, $R\succeq 0$ ($R$ can be singular), and $Q$ ($Q$ can be singular) such that $E^{T}Q\succeq 0$, i. e.

[TABLE]

Then $\mathbb{S}_{DH}^{\succeq 0}$ is the closure $\overline{\mathbb{S}_{DH}^{\succ 0}}$ of $\mathbb{S}_{DH}^{\succ 0}$. Following the arguments similar to that of Lemma 3 we also have that $\mathbb{S}_{DH}^{\succeq 0}$ is invariant under orthogonal transformations and, furthermore, we have that

[TABLE]

Note, however, that the set $\mathbb{S}_{DH}^{\succeq 0}$ is not bounded, and hence the infimum in the right hand side of (3.1) may not be attained.

Our observations lead to the following reformulation of the nearest stable matrix pair problem.

Theorem 5**.**

Let $(E,A)\in{\mathbb{R}}^{n,n}\times{\mathbb{R}}^{n,n}$. Then

[TABLE]

For the standard system (when $E$ is the identity matrix in (1.1)) stability solely depends on the eigenvalues of the matrix $A$. Thus making $A$ stable without perturbing the identity matrix gives an upper bound for the distance of $(I_{n},A)$ to the nearest stable matrix pair. This also follows from (3.2) because

[TABLE]

We note that the infimum on the right hand side of (3.3) is the distance of $A$ from the set of stable matrices [12]. We will demonstrate in our numerical experiments that (as expected) the inequality in (3.3) is strict. However, it is interesting to note that if $(J^{*},R^{*},Q^{*})$ is a stationary point of the right hand side of (3.3) then $(I_{n},J^{*},R^{*},Q^{*})$ is a stationary point of the left hand side of (3.3); see the discussion in Section 4.

A DAE with coefficient pair $(E,A)$ and nonsingular $E$ can be equivalently reformulated as a standard system $(I_{n},AE^{-1})$, and then stability can be determined by the eigenvalues of $AE^{-1}$. Thus again making $AE^{-1}$ stable, or equivalently, $(I_{n},AE^{-1})$ stable without perturbing the identity matrix, gives an upper bound for the distance of $(E,A)$ to the nearest stable matrix pair. Indeed, w.l.o.g. if we scale $(E,A)$ so that ${\|E\|}_{F}=1$, then

[TABLE]

where (3.4) follows by first using the fact that $Q^{T}E=E^{T}Q\succeq 0$ implies that $E^{-T}Q^{T}=QE^{-1}\succeq 0$, whereas in (3.5) we use the sub-multiplicativity of the Frobenius norm and ${\|E\|}_{F}=1$. Thus, for nonsingular $E$, an upper bound can be computed by the reduction to a standard system. However, the reduction process may be numerically unstable. Hence, it is advisable to work directly with the pair $(E,A)$ even for regular matrix pairs with invertible $E$.

Remark 1**.**

The eigenvalues of the regular DH pair $(E,(J-R)Q)$ can be highly sensitive to small perturbations. This happens if the finite or infinite eigenvalues are multiple and defective or if they are close to being multiple and the eigenvectors have a small angle. In fact, let $\lambda$ be a simple eigenvalue of $(E,(J-R)Q)$ and let $x$ be the corresponding eigenvector normalized to unit norm. By (2.5) we have

[TABLE]

which implies that $(-\mathop{\mathrm{Re}}{(\lambda)},x)$ is an eigenvalue/eigenvector pair of the symmetric pencil $L(z):=zQ^{T}E-Q^{T}RQ$, since $Q^{T}E\succeq 0$ and $Q^{T}RQ\succeq 0$. Then, see e. g. [27], the normwise condition number of $\lambda$ of $L(z)$ is given by

[TABLE]

which implies that

[TABLE]

Therefore, if $\kappa(-\mathop{\mathrm{Re}}{(\lambda)},L)$ is large, then a small perturbation can significantly perturb the eigenvalues of $(E,(J-R)Q)$. **

4 Algorithmic solution to the nearest stable matrix pair problem

In this section, we propose an algorithm to solve the nearest stable matrix pair problem using the reformulation of Theorem 5, i. e., to solve

[TABLE]

One of the first algorithms that comes to mind to solve (4.1) is a block-coordinate descent (BCD) method (see for example [30]) with three blocks of variables $(J,R)$, $Q$ and $M$. In this case the subproblems for each of the three blocks are convex when the others are fixed (these are least-squares problems with linear or positive semidefinite constraints). However, we have observed that this method is not very efficient in practice, since convergence is slow when one gets close to a stationary point. Moreover, the BCD method can get easily stuck in saddle points. In particular, consider a matrix pair of the form $(I_{n},A)$ and let $(J,R,Q)$ be a fixpoint of the BCD method (with two blocks of variables $(J,R)$ and $Q$) for the nearest stable matrix problem

[TABLE]

Then $(I_{n},J,R,Q)$ is a fixpoint of the BCD method (with three blocks of variables $(J,R)$, $Q$ and $M$) for (4.1). By construction, then the update of $(J,R)$ and $Q$ cannot be improved for $M=I_{n}$ fixed, while $M=I_{n}$ is optimal, since it is feasible and minimizes ${\|M-I_{n}\|}_{F}$. This behavior of the BCD method is illustrated in the following example.

Example 3**.**

Consider a matrix pair $(E,A)$, where

[TABLE]

with eigenvalues $1$, $1\pm\sqrt{2}i$. The representation $(J,R,Q,I_{n})$ with

[TABLE]

is a fixpoint of the BCD method for (4.1) with error ${\|A-(J-R)Q\|}_{F}^{2}=3$ and we strongly believe that $(J-R)Q$ is the nearest stable matrix to $A$.

Using our fast projected gradient method that will be introduced below, initialized with this solution, we obtain a nearest stable matrix pair with much lower distance

[TABLE]

One possible reason why the formulation (4.1) does not seem to lead to good solutions is that the constraint $Q^{T}M\succeq 0$ couples the variables $Q$ and $M$ rather strongly. This motivated us to introduce a further reformulation of (4.1), where the feasible set is free of any constraint that involves coupling of some of its variables and onto which it is easy to project.

4.1 Reformulation of (4.1)

The formulation (4.1) is not convenient for standard optimization schemes, such as a projected gradient method, because the feasible set is highly non-convex with the constraints $Q^{T}M\succeq 0$, and this constraint seems to prevent standard optimization schemes to converge to good solutions, as demonstrated in Example 3. To overcome this dilemma we introduce a new variable $H=Q^{T}M=M^{T}Q\succeq 0$ in (4.1), which corresponds to the Hamiltonian for port-Hamiltonian descriptor systems, see [2]. This results in another reformulation of (4.1) with a modified feasible set and objective function

[TABLE]

The feasible set of (4.2) is rather simple, with no coupling of the variables, and it is relatively easy to project onto it.

Remark 2**.**

If one want to obtain a DH pair with $Q$ invertible and $R\succ 0$ so that $(M,(J-R)Q)$ is regular, asymptotically stable and is of index at most one, see Theorem 4, then the constraints $R\succeq 0$ and $H\succeq 0$ can be replaced with $R\succeq\delta I_{n}$ and $H\succeq\delta I_{n}$ for some small parameter $\delta>0$ that prevents $R$ and $H=Q^{T}M$ from being singular. As a consequence then $Q$ will be invertible, since otherwise $Q^{-T}H$ is unbounded. Numerically, it does not make the problem much different, the projection is still rather straightforward. Furthermore, in view of Remark 1, the constraints $R\succeq\delta I_{n}$ and $H\succeq\delta I_{n}$ will assure that every finite eigenvalue $\lambda$ of the DH pair $(M,(J-R)Q)$ satisfies

[TABLE]

where $\kappa{(-\mathop{\mathrm{Re}}{(\lambda)},L)}$ is the condition number of $-\mathop{\mathrm{Re}}{(\lambda)}$ being a simple eigenvalue of the symmetric pencil $L(z)=zQ^{T}M-Q^{T}RQ$. **

4.2 Fast projected gradient algorithm

Following the spirit of [12], we use a fast projected gradient algorithm to solve (4.2) despite the fact that it is a non-convex problem. The only non-trivial part in implementing this method is to compute the gradient of the objective function with respect to $Q$. We show in Appendix A that

[TABLE]

where $f(J,R,Q,H):={\|A-(J-R)Q\|}_{F}^{2}+{\|E-Q^{-T}H\|}_{F}^{2}$. A pseudocode for our fast projected gradient methods is presented in Algorithm 1.

Initialization

To initialize the fast projected gradient approach we use a similar strategy as in the matrix case, see [12]. For $Q=I_{n}$, the optimal solutions for $J$, $R$ and $H$ in (4.2) can be computed easily, and are given by

[TABLE]

where $\mathcal{P}_{\succeq}$ is the projection onto the set of positive semidefinite matrices.

Because of the non-convex nature of the problem, the solution obtained by our algorithm is highly sensitive to the initial point. We observed that the above initial point provides in general good solutions. Future work will require the development or more sophisticated initialization schemes and possibly globalization approaches to obtain better solutions.

Line-search

For $Q$ fixed, the Lipschitz constant of the gradient of $f$ with respect to $H$ (resp. $(J,R)$) is given by $L=\lambda_{\min}(QQ^{T})^{-1}$ (resp. $L=\lambda_{\max}(QQ^{T})$), where $\lambda_{\min}(X)$ (resp. $\lambda_{\max}(X)$) denotes the smallest (resp. largest) eigenvalue of matrix $X$. Hence, performing a gradient step with step length $\gamma=1$ in the variable $H$ (resp. $(J,R)$) will guarantee the decrease of the objective function. Since we initialize with $Q=I_{n}$, we choose an initial step length of $\gamma=1$. If it does not lead to a decrease of the objective function, then we decrease $\gamma$ by a constant factor (we used 2/3). This is a standard backtracking line-search strategy. If no decrease of the objective function is possible, i. e. if the step length is smaller than some threshold $\underline{\gamma}$, then we restart the fast gradient scheme with a standard gradient step, which guarantees convergence. At the next step, we use a slightly larger step length than the previous step to avoid the step length to go to zero, e. g. by multiplying $\gamma$ by a constant factor; we used 2.

5 Numerical Examples

Our optimization code is available from https://sites.google.com/site/nicolasgillis/ and the numerical examples presented below can be directly run from this online code. The algorithm runs in $O(n^{3})$ operations, including projections on the set of positive semidefinite matrices, inversion of the matrix $Q$ and all necessary matrix-matrix products. Our algorithm can be applied on a standard desktop computer with $n$ up to a thousand. As far as we know, no other algorithm has been proposed for the computation of the nearest stable matrix pair. Therefore, we cannot compare to an existing method, so for illustration of our results we will only compare our fast gradient method with the projected gradient method (GM), which is simply FGM where restart is used at every step. For all experiments, we limit the CPU time to 10 seconds.

5.1 Case 1: $(I_{n},A)$

Let us take as $A$ the Grcar matrix from [12] and $E=I_{n}$. The Grcar matrix of order $k$ is a banded Toeplitz matrix with its subdiagonal set to $-1$, and both its main and $k$ superdiagonals set to 1.

We use $n=20$ and $k=3$. The nearest stable matrix $\tilde{A}$ to $A$ found in [12] satisfies ${\|A-\tilde{A}\|}_{F}^{2}=23.51$. Figure 5.1 shows the evolution of the objective function of (4.1). Our FGM achieves a objective function value ${\|A-(J-R)Q\|}_{F}^{2}+{\|E-Q^{-T}H\|}_{F}^{2}=6.28$. Note that FGM converges much faster than GM.

5.2 Case 2: Random $A$ and rank-deficient $E$

In our second example, we generate each entry of $A$ using the normal distribution of mean [math] and standard deviation $1$ (randn(n) in Matlab). We generate $E$ so as it has rank $r$, by taking $E$ as the best rank-$r$ approximation of a matrix generated randomly in the same way as $A$. We use $n=20$ and $r=3$.

For these types of matrices, usually the algorithm converges to a rank-deficient solution, i. e. $R,H,Q$ are not of full rank, hence it is useful to use the lower bound for the eigenvalues to obtain better numerical stability; we use $\delta=10^{-6}$.

Figure 5.2 shows the evolution of the objective function of (4.1) for a particular example. On all the examples that we have run, FGM always converged much faster and generated a significantly better solution than GM, similarly as shown on Figure 5.2.

5.3 Case 3: Mass-spring-damper model

Let us consider a lumped parameter, mass-spring-damper dynamical system, see, e. g., [29] with $n$ point masses and $n$ spring-damper pairs. The equations of motion can be written in the form of second order system of differential equations as

[TABLE]

where $M\succ 0$ is the mass matrix, $D\succeq 0$ is the damping matrix and $K\succ 0$ is the stiffness matrix, $f,u,\dot{u},\ddot{u}$ are the forces, displacement, velocity and acceleration vectors, respectively.

An appropriate first order formulation of (5.1) is associated with the DH pencil $\lambda E-A$, where

[TABLE]

Since $E$ is invertible, the pair $(E,A)$ is regular and of index zero, and thus from Theorem 1 is also stable. In order to make the pair unstable, we perturb the dissipation matrix to become indefinite via

[TABLE]

We use $n=10$ and $\epsilon=0.1$ which moves some eigenvalues to the open right half of the complex plane (see Figure 5.4 below).

The corresponding mass vector $m$ containing the point masses, spring vector $k$ containing the spring constants, and damping vector $c$ containing the damping constants are all equal to the vector

[TABLE]

Figure 5.3 shows the evolution of the objective function (4.2). We used the initialization from (5.2), i. e. we used the true $J$, $R$ and $Q$ from the known stable system to see whether an approximation better than the original stable pencil can be found. (Note that if we use our standard initialization, see Section 4.2, FGM converges to another local minimum with higher approximation error $32.70$.)

We observe that FGM converges much faster than GM, as in the two other examples, while it is able to generate a better approximation than the FGM from [12] applied on the nearest stable matrix problem $(I,E^{-1}A)$, as already noted in Section 5.1. We also observe that the algorithms are able, rather surprisingly, to significantly reduce the value of the objective compared to the initialization (which is the original stable pencil), from $21.97$ to $6.53$ for FGM (DH matrix), to $12.70$ for GM (DH pair), and to $4.09$ for FGM (DH pair). Actually, if we run FGM for $5$ more seconds, then it converges to a solution with a value $3.81$.

Figure 5.4 shows the location of the eigenvalues of the pencils $(E,A)$ generated by the FGM using a DH matrix and a DH pair, as well as those of the original and perturbed pencils.

Conclusion

We have considered the problem of computing the nearest stable matrix pair to a given unstable one. We have presented a new characterization of the set of stable matrix pairs using DH matrix pairs, defined as $(E,(J-R)Q)$ where $J=-J^{T}$, $R\succeq 0$, and $Q^{T}E\succeq 0$ with $Q$ invertible. This has allowed us to reformulate the nearest stable matrix pair problem in a more convenient way for which we develop an efficient fast gradient method.

Possible future work is related to develop more sophisticated ways to solve (4.2) and to apply globalization approaches to obtain better solutions.

Appendix A Gradient of $f$ with respect to $Q$

In this section, we explain how to compute the gradient of

[TABLE]

with respect to $Q$, denoted by $\nabla_{Q}f(Q)$. For the first term, we have

[TABLE]

For the second term, we need to compute $\nabla_{Q}{\|E-Q^{-T}H\|}_{F}^{2}$. Using the basic rules for matrix differentiation, see e. g. [18, 25],

1. (i)

Linear:

[TABLE] 2. (ii)

Product:

[TABLE] 3. (iii)

Inverse:

[TABLE]

we obtain (considering the transpose matrix $E^{T}-HQ^{-1}$),

[TABLE]

Using (ii) with $F(Q)=HQ^{-1}$ and $G(Q)=F(Q)^{T}$, we get for the first term

[TABLE]

Using (ii) again, now with $G(Q)=HQ^{-1}$ and $F(Q)=Y^{-T}H^{T}$ which is a constant, we obtain

[TABLE]

using (iii), $W=Y^{-T}H^{T}H$, and (i),

[TABLE]

For the second term we use (iii) and obtain

[TABLE]