# Computation of the nearest structured matrix triplet with common null space

Nicola Guglielmi111Gran Sasso Science Institute, Section of Mathematics, Via
Crispi 7, I-67010 L’ Aquila, Italy. Email<EMAIL_ADDRESS>Volker
Mehrmann33footnotemark: 3

###### Abstract

We study computational methods for computing the distance to singularity, the
distance to the nearest high index problem, and the distance to instability
for linear differential-algebraic systems (DAEs) with dissipative Hamiltonian
structure. While for general unstructured DAEs the characterization of these
distances is very difficult, and partially open, it has been shown in [32]
that for dissipative Hamiltonian systems and related matrix pencils there
exist explicit characterizations. We will use these characterizations for the
development of computational methods to compute these distances via methods
that follow the flow of a differential equation converging to the smallest
perturbation that destroys the property of regularity, index one or stability.

###### keywords:

dissipative Hamiltonian systems, structured distance to singularity,
structured distance to high index problem, structured distance to instability,
low-rank perturbation, differential-algebraic system

###### AMS:

15A18, 15A21, 65K05, 15A22

22footnotetext: Institut für Mathematik, MA 4-5, TU Berlin, Str. der 17 Juni
136, D–10623 Berlin, Germany. Email<EMAIL_ADDRESS>
##  1 Introduction

We derive computational methods for determining the distance to singularity,
the distance to the nearest high index problem, and the distance to
instability for _linear, time-invariant differential-algebraic systems (DAEs)
with dissipative Hamiltonian structure (dHDAEs)_. Such systems arise as
linearization of general dHDAEs along a stationary solution and have the form

$E\dot{x}=(J-R)x+f,$ (1)

with constant coefficient matrices $E,J,R\in{\mathbb{R}}^{n,n}$,
$J=-J^{\top}$, and $E=E^{\top},R=R^{\top}$ symmetric positive semidefinite, a
differentiable state function $x:{\mathbb{R}}\to{\mathbb{R}}^{n}$ and a right
hand side $f:{\mathbb{R}}\to{\mathbb{R}}^{n}$, see [5, 18, 24, 31, 32, 33, 37,
40, 38, 39] for slightly varying definitions and a detailed analysis of such
systems also in the context of the more general _port-Hamiltonian systems_.
The matrix $E$ is associated with the Hessian of the associated _Hamiltonian_
energy function, which in the quadratic case has the form
$\mathcal{H}(x)=\frac{1}{2}x^{\top}Ex$. It is well-known [5, 33, 40] that
pHDAEs satisfy a _dissipation inequality_
${\mathcal{H}}\big{(}x(t_{1})\big{)}-{\mathcal{H}}\big{(}x(t_{0})\big{)}\leq
0$ for $t_{1}\geq t_{0}$.

Such pHDAE systems arise in all areas of science and engineering [4, 5, 14,
33, 40] as linearizations, space discretization, or approximation of physical
systems and are usually model descriptions with uncertainties. It is therefore
important to know whether the model is close to an ill-posed or badly
formulated model, and this has been an important research topic recently, see
[1, 3, 6, 18, 19, 20, 28, 29, 32, 34]. Since the system properties of (1) are
characterized by investigating the corresponding _dissipative Hamiltonian (dH)
matrix pencil_

$L(\lambda):=\lambda E-(J-R),$ (2)

the discussed nearness problems can be characterized by determining the
distance to the nearest singular pencil, i.e., a pencil with a $\det(\lambda
E-J+R)$ identically zero, the distance to the nearest high-index problem,
i.e., a problem with Jordan blocks associated to the eigenvalue $\infty$ of
size bigger than one, or the nearest problem on the boundary of the unstable
region, i.e. a problem with purely imaginary eigenvalues. To compute these
distances is very difficult for general linear systems [7, 8, 10, 18, 19, 20,
22, 30]. However, if one restricts the perturbations to be structured, i.e.
one considers _structured distances_ within the class of linear time-invariant
dHDAEs, the situation changes completely, see [18, 19, 20, 31, 32], and one
obtains very elegant characterizations that can be used in numerical methods
to compute these distances.

These methods are usually based on non-convex optimization approaches. In
contrast to such approaches, we derive computational methods to compute these
structured distances by following the flow of a differential equation. This
approach has been shown to be extremely effective for computing the distance
to singularity for general matrix pencils [22] and we will show that this
holds even more so in the structured case.

Neither the methods based on non-convex optimization nor the methods based on
following a flow for general pencils or structured pencils are really feasible
for large scale problems. To treat the large sparse case they have to be
combined with projections on the sparsity structure and model reduction
methods, see [2, 3], which intertwine the optimization step with model
reduction via interpolation. Here we discuss only the small scale case, but
the combination with interpolation methods can be carried out in an analogous
way as in [3].

The paper is organized as follows. In Section 2 we recall a few basic results
about linear time-invariant dHDAE systems. In Section 3 we discuss
optimization methods that are based on gradient flow computations. Since the
cases of even and odd dimension are substantially different, in Section 4 we
specialize these methods for the optimization problems associated with the
three discussed distance problems for the case that the state dimension is
odd, while in Section 6 we discuss the case that the state dimension is even.
Since it is known that the optimal perturbations are rank two matrices, in
Section 5 for the odd size case we discuss the special situation that we
restrict the perturbation to be at most of rank two. In Section 8 we briefly
discuss the iterative procedure for computing the optimal $\varepsilon$ in the
upper level of the two level procedure. In all cases, we present numerical
examples.

##  2 Preliminaries

We use the following notation. The set of symmetric (positive semidefinite)
matrices in $\mathbb{R}^{n,n}$ is denoted by $\mathrm{Sym}^{n,n}$
($\mathrm{Sym}^{n,n}_{\geq 0}$), and the skew-symmetric matrices in
$\mathbb{R}^{n,n}$ by $\mathrm{Skew}^{n,n}$. By $\|X\|_{F}$ we denote the
Frobenius norm of a (possibly rectangular) matrix $X$, we extend this norm to
matrix tuples $\mathcal{X}=(X_{0},\dots,X_{k})$ via
$\|\mathcal{X}\|_{F}=\|[X_{0},\dots,X_{k}]\|_{F}$. For
$A,B\in{\mathbb{C}}^{n,n}$, we denote by

$\langle A,B\rangle=\text{tr}(B^{H}A)$

the Frobenius inner product on ${\mathbb{C}}^{n,n}$, where $B^{H}$ is the
conjugate transpose of $B$. The Euclidian norm in $\mathbb{R}^{n}$ is denoted
by $\|~{}\|$. By $\lambda_{\min}(X)$ we denote the smallest eigenvalue of
$X\in\mathrm{Sym}^{n,n}_{\geq 0}$. The real and imaginary part of a complex
matrix $A\in{\mathbb{C}}^{n,n}$ is denoted by ${\mbox{\rm Re}}(A)$,
${\mbox{\rm Im}}(A)$, respectively.

To characterize the properties of dHDAEs of the form (1), we exploit the
_Kronecker canonical form_ of the associated matrix pencil (2), see [17]. If
$\mathcal{J}_{n}(\lambda_{0})$ denotes the standard upper triangular Jordan
block of size $n\times n$ associated with an eigenvalue $\lambda_{0}$ and
$\mathcal{L}_{n}$ denotes the standard right Kronecker block of size
$n\times(n+1)$, i.e.,

$\mathcal{L}_{n}=\lambda\left[\begin{array}[]{cccc}1&0\\\ &\ddots&\ddots\\\
&&1&0\end{array}\right]-\left[\begin{array}[]{cccc}0&1\\\ &\ddots&\ddots\\\
&&0&1\end{array}\right],$

then for $E,A\in{\mathbb{C}}^{n,m}$ there exist nonsingular matrices
$S\in{\mathbb{C}}^{n,n}$ and $T\in{\mathbb{C}}^{m,m}$ that transform the
pencil to _Kronecker canonical form_ ,

$S(\lambda E-A)T=\mathop{\operator@font diag}\nolimits({\cal
L}_{\epsilon_{1}},\ldots,{\cal L}_{\epsilon_{p}},{\cal
L}^{\top}_{\eta_{1}},\ldots,{\cal L}^{\top}_{\eta_{q}},{\cal
J}_{\rho_{1}}^{\lambda_{1}},\ldots,{\cal J}_{\rho_{r}}^{\lambda_{r}},{\cal
N}_{\sigma_{1}},\ldots,{\cal N}_{\sigma_{s}}),$ (3)

where
$p,q,r,s,\epsilon_{1},\dots,\epsilon_{p},\eta_{1},\dots,\eta_{q},\rho_{1},\dots,\rho_{r},\sigma_{1},\dots,\sigma_{s}\in\mathbb{N}_{0}$
and $\lambda_{1},\dots,\lambda_{r}\in\mathbb{C}$, as well as ${\cal
J}_{\rho_{i}}^{\lambda_{i}}=I_{\rho_{i}}-\mathcal{J}_{\rho_{i}}(\lambda_{i})$
for $i=1,\dots,r$ and
$\mathcal{N}_{\sigma_{j}}=\mathcal{J}_{\sigma_{j}}(0)-I_{\sigma_{j}}$ for
$j=1,\dots,s$.

For real matrices and real transformation matrices $S,T$, the blocks ${\cal
J}_{\rho_{j}}^{\lambda_{j}}$ with
$\lambda_{j}\in\mathbb{C}\setminus\mathbb{R}$ are in _real Jordan canonical
form_ associated to the corresponding pair of conjugate complex eigenvalues,
the other blocks are the same. A real or complex eigenvalue is called
_semisimple_ if the largest associated Jordan block in the complex Jordan form
has size one and the sizes $\eta_{j}$ and $\epsilon_{i}$ are called the _left
and right minimal indices_ of $\lambda E-A$, respectively. A pencil $\lambda
E-A$, is called _regular_ if $n=m$ and $\det(\lambda_{0}E-A)\neq 0$ for some
$\lambda_{0}\in\mathbb{C}$, otherwise it is called _singular_ ;
$\lambda_{1},\dots,\lambda_{r}\in\mathbb{C}$ are called the finite eigenvalues
of $\lambda E-A$, and $\lambda_{0}=\infty$ is an eigenvalue of $\lambda E-A$
if zero is an eigenvalue of the $\lambda A-E$. The size of the largest block
${\cal N}_{\sigma_{j}}$ is called the _index_ $\nu$ of the pencil $\lambda
E-A$.

The definition of stability for differential-algebraic systems varies in the
literature. We call a pencil $\lambda E-A$ Lyapunov stable (asymptotically
stable) if it is regular, all finite eigenvalues are in the closed (open) left
half plane, and the ones lying on the imaginary axis (including $\infty$) are
semisimple [13]. Note that pencils with eigenvalues on the imaginary axis or
at $\infty$ are on the boundary of the set of asymptotically systems and those
with multiple, but semisimple, purely imaginary eigenvalues (including
$\infty$) lie on the boundary of the set of Lyapunov stable pencils.

The following theorem summarizes some results of [31, 32] for real dH pencils;
note that some of the results also hold in the complex case.

######  Theorem 1.

Let $E,R\in\mathbb{R}^{n,n}$ be symmetric and positive semidefinite, and
$J=-J^{\top}\in\mathbb{R}^{n,n}$. Then the following statements hold for the
pencil $L(\lambda)=\lambda E-J+R$.

  1. (i)

If $\lambda_{0}\in\mathbb{C}$ is an eigenvalue of $L(\lambda)$ then
$\operatorname{Re}(\lambda_{0})\leq 0$.

  2. (ii)

If $\omega\in\mathbb{R}$ and $\lambda_{0}=i\omega$ is an eigenvalue of
$L(\lambda)$, then $\lambda_{0}$ is semisimple. Moreover, if the columns of
$V\in\mathbb{C}^{m,k}$ form a basis of a regular deflating subspace of
$L(\lambda)$ associated with $\lambda_{0}$, then $RQV=0$.

  3. (iii)

The index of $L(\lambda)$ is at most two.

  4. (iv)

All right and left minimal indices of $L(\lambda)$ are zero (if there are
any).

  5. (v)

The pencil $L(\lambda)$ is singular if and only if $\,\ker J\cap\ker E\cap\ker
R\neq\\{0\\}$.

Based on Theorem 1, in [32] the following distance problems were introduced
for dH pencils.

######  Definition 2.

Let $\mathcal{L}$ denote the class of square $n\times n$ real matrix pencils
of the form (2). Then

  1. 1.

the _structured distance to singularity_ is defined as

$d_{\operatorname{sing}}^{\mathcal{L}}\big{(}L(\lambda)):=\inf\big{\\{}\big{\|}\Delta_{L}(\lambda)\big{\|}_{F}\
\big{|}\ L(\lambda)+\Delta_{L}(\lambda)\in\mathcal{L}\mbox{ and is
singular}\big{\\}};$ (4)

  2. 2.

the _structured distance to the nearest high-index problem_ is defined as

$d_{\operatorname{hi}}^{\mathcal{L}}\big{(}L(\lambda)):=\inf\big{\\{}\big{\|}\Delta_{L}(\lambda)\big{\|}_{F}\
\big{|}\ L(\lambda)+\Delta_{L}(\lambda)\in\mathcal{L}\mbox{ and is of
index}\geq 2\big{\\}};$ (5)

  3. 3.

the _structured distance to instability_ is defined as

$d_{\operatorname{inst}}^{\mathcal{L}}\big{(}L(\lambda)\big{)}:=\inf\big{\\{}\big{\|}\Delta_{L}(\lambda)\big{\|}_{F}\
\big{|}\ L(\lambda)+\Delta_{L}(\lambda)\in\mathcal{L}\mbox{ and is
unstable}\big{\\}}.$ (6)

Here $\Delta_{L}(\lambda)=\lambda\Delta_{E}-\Delta_{J}+\Delta_{R}$, with
$\Delta_{J}\in\mathrm{Skew}^{n,n}$,
$E+\Delta_{E},R+\Delta_{R}\in\mathrm{Sym}^{n,n}_{\geq 0}$, and
$\big{\|}[\Delta_{J},\Delta_{R},\Delta_{E}]\big{\|}_{F}=\big{\|}[\Delta_{L}(\lambda)]\big{\|}_{F}$.

It has also been shown in [32] that these distances can be characterized as
follows.

######  Theorem 3.

Let $\lambda E-J+R\in\mathcal{L}$. Then the following statements hold.

  1. 1.

Define for a matrix $Y\in{\mathbb{R}}^{n,n}$, the matrix
$\Delta_{Y}^{u}=-uu^{\top}Y-Yuu^{\top}+uu^{\top}Yuu^{\top}$. The distance to
singularity (4) is attained with a perturbation $\Delta_{E}=\Delta_{E}^{u}$,
$\Delta_{J}=\Delta_{J}^{u}$, and $\Delta_{R}=\Delta_{R}^{u}$ for some
$u\in\mathbb{R}^{n}$ with $\|u\|_{2}=1$. It is given by

$d_{\operatorname{sing}}^{\mathcal{L}}\big{(}\lambda E-J+R\big{)}\\\
=\min_{u\in{\mathbb{R}}^{n}\atop\|u\|=1}\sqrt{2\|Ju\|^{2}+2\big{\|}(I-uu^{\top})Eu\big{\|}^{2}+(u^{\top}Eu)^{2}+2\big{\|}(I-uu^{\top})Ru\big{\|}^{2}+(u^{\top}Ru)^{2}}$

and is bounded as

$\sqrt{\lambda_{\min}(-J^{2}+R^{2}+E^{2})}\leq
d_{\operatorname{sing}}^{\mathcal{L}}\big{(}\lambda
E-J+R\big{)}\leq\sqrt{2\cdot\lambda_{\min}(-J^{2}+R^{2}+E^{2})}.$ (7)

  2. 2.

The structured distance to higher index (5) and the structured distance to
instability (6) coincide and satisfy

$d_{\operatorname{hi}}^{\mathcal{L}}\big{(}\lambda
E-J+R\big{)}=d_{\operatorname{inst}}^{\mathcal{L}}\big{(}\lambda
E-J+R\big{)}\\\
=\min_{u\in{\mathbb{R}}^{n}\atop\|u\|=1}\sqrt{2\big{\|}(I-uu^{\top})Eu\big{\|}^{2}+(u^{\top}Eu)^{2}+2\big{\|}(I-uu^{\top})Ru\big{\|}^{2}+(u^{\top}Ru)^{2}}$

and are bounded as

$\sqrt{\lambda_{\min}(E^{2}+R^{2})}\leq
d_{\operatorname{hi}}^{\mathcal{L}}\big{(}\lambda
E-J+R\big{)}=d_{\operatorname{inst}}^{\mathcal{L}}\big{(}\lambda
E-J+R\big{)}\leq\sqrt{2\cdot\lambda_{\min}(E^{2}+R^{2})}.$ (8)

With formulas and close upper and lower bounds available, these distances can
be computed by global constrained optimization methods such as [36]. Based on
our experience in computing the distance to instability for general matrix
pencils, where different computational methods were studied and it was shown
that gradient flow methods were extremely efficient, in the next section we
introduce such gradient methods to compute the discussed structured distances.

##  3 ODE-based gradient flow approaches

In the previous section we have seen that for dH pencils the distance to
singularity is characterized by the distance to the nearest common nullspace
of three structured matrices and the distance to high index and instability by
the distance to the nearest common nullspace of two symmetric positive
definite matrices, with perturbations that keep the structure.

The perturbation matrices that give the structured distance to singularity can
be alternatively expressed as

$\displaystyle(\Delta E_{*},\Delta R_{*},\Delta J_{*})$ $\displaystyle=$
$\displaystyle\arg\min\limits_{\Delta E,\Delta R,\Delta J}\|(\Delta E,\Delta
R,\Delta J)\|$ (9) subj. to $\displaystyle E+\Delta E,R+\Delta
R\in\mathrm{Sym}^{n,n}_{\geq 0},\Delta J\in\mathrm{Skew}^{n,n},$ and
$\displaystyle(E+\Delta E)x=0,\quad(R+\Delta R)x=0,\quad(J+\Delta J)x=0$
$\displaystyle\mbox{for some}\ x\in\mathbb{R}^{n},x\neq 0.$

Then $d_{\operatorname{sing}}^{\mathcal{L}}\big{(}\lambda
E-J+R\big{)}=\|(\Delta E_{*},\Delta R_{*},\Delta J_{*})\|$ and our algorithmic
approach to minimize this functional is based on this reformulation.

###  3.1 A two-level minimization

To determine the minimum in (9) we use a two-level minimization. As an inner
iteration, for a perturbation size $\varepsilon$, we consider perturbed
matrices $E+\varepsilon\Delta$, $R+\varepsilon\Theta$ and
$J+\varepsilon\Gamma$ with $\|(\Delta,\Theta,\Gamma)\|_{F}\leq 1$ satisfying
the constraints in (9). Let us denote

  * (i)

by $(\lambda,x)$ an eigenvalue/eigenvector pair of $E+\varepsilon\Delta$
associated with the smallest eigenvalue and $\|x\|=1$;

  * (ii)

by $(\nu,u)$ an eigenvalue/eigenvector pair of $R+\varepsilon\Theta$
associated with the smallest eigenvalue, and $\|u\|=1$;

  * (iii-a)

if $n$ is even, by $({\bf i}\mu,w)$ an eigenvalue/eigenvector pair of
$J+\varepsilon\Gamma$, with $\mu>0$ such that ${\bf i}\mu$ is the eigenvalue
with smallest positive imaginary part and $\|w\|=1$,

  * (iii-b)

if $n$ is odd, by $(0,w)$ an eigenvalue/eigenvector pair of
$J+\varepsilon\Gamma$ (this exists for all $\Gamma$).

In the inner iteration, for any fixed $\varepsilon$ we compute a (local)
minimizer of (9) that is, however, different for even or odd $n$.

### The case that $n$ is odd

In this case the skew-symmetric matrix always has a zero eigenvalue (with an
associated real eigenvector $w$) so that the only contribution to the
optimization is through the alignment of $w$ with $x$ and $u$. Hence, the
functional to be minimized in (9) can be expressed in the simplified form

$F_{\varepsilon}^{\mathrm{od}}(\Delta,\Theta,\Gamma)=\frac{1}{2}\Big{(}\lambda^{2}+\nu^{2}+1-(x^{\top}u)^{2}+1-(x^{\top}w)^{2}\Big{)}$
(10)

with $\|\left(\Delta,\Gamma,\Theta\right)\|_{F}\leq 1$. It is, however,
possible to include a further term $1-|u^{\top}w|^{2}$ in the functional,
which does not change the solution but may have an impact on the conditioning
of the problem and hence the numerical performance.

### The case that $n$ is even

In this case, when two eigenvalues $\pm{\bf i}\mu$ ($\mu>0$) coalesce at $0$,
they form a semi-simple double eigenvalue and the associated eigenvectors
$w=w_{1}+{\bf i}w_{2}$ and $\overline{w}=w_{1}-{\bf i}w_{2}$ form a two-
dimensional nullspace spanned by the two real vectors $w_{1}$ and $w_{2}$.
These can be assumed to be orthogonal to each other, i.e.
$w_{1}^{\top}w_{2}=0$ and have the same norm $1/\sqrt{2}$ so that still
$\|w\|=1$. Using $w_{1},w_{2}$, we define the real orthogonal matrix

$W=\sqrt{2}\left[w_{1},w_{2}\right],$

and to satisfy the constraint in (9), we require that

$Wz=x\qquad\mbox{for some}\ z\in\mathbb{R}^{2}.$

This leads to the minimization of

$\|Wz-x\|\qquad\mbox{for some}\ z\in\mathbb{R}^{2}.$

Since $W$ is orthogonal, the solution is $z=W^{\top}x$, and the functional to
be minimized takes the form

$1-x^{\top}WW^{\top}x=1-2\,(x^{\top}w_{1})^{2}-2\,(x^{\top}w_{2})^{2},$

which is positive if $x$ does not lie in the range of $W$ and zero otherwise.

In summary, the functional in the even case is given by

$F_{\varepsilon}^{\mathrm{ev}}(\Delta,\Theta,\Gamma)=\frac{1}{2}\Big{(}\lambda^{2}+\nu^{2}+\mu^{2}+1-(x^{\top}u)^{2}+1-2\,(x^{\top}{\mbox{\rm
Re}}(w))^{2}-2\,(x^{\top}{\mbox{\rm Im}}(w))^{2}\Big{)}$ (11)

with $\|\left(\Delta,\Gamma,\Theta\right)\|_{F}\leq 1$.

######  Remark 4.

In both the odd and the even case we have that

$\min_{\|\left(\Delta,\Gamma,\Theta\right)\|_{F}\leq
1}F_{\varepsilon}(\Delta,\Theta,\Gamma)=\min_{\|\left(\Delta,\Gamma,\Theta\right)\|_{F}=1}F_{\varepsilon}(\Delta,\Theta,\Gamma).$

To see this, consider $(\Delta_{*},\Theta_{*},\Gamma_{*})$ of Frobenius norm
less than or equal to $1$ giving a minimizer of the left-hand side and suppose
that $(x_{*},\lambda_{*})$ is the minimizing eigenvalue/eigenvector pair of
$E+\varepsilon\Delta_{*}$. Then choosing a matrix $\widetilde{\Delta}\neq 0$
such that $\widetilde{\Delta}x_{*}=0$ and
$\langle\Delta_{*},\widetilde{\Delta}\rangle=0$, for a suitable $\theta$ the
matrix $\Delta_{1}=\Delta_{*}+\theta\widetilde{\Delta}$ is of unit Frobenius
norm and has
$F_{\varepsilon}(\Delta_{1},\Theta_{*},\Gamma_{*})=F_{\varepsilon}(\Delta_{*},\Theta_{*},\Gamma_{*})$.

Using the functionals (10), respectively (11), in our approach the local
minimizer of
$\min_{\|\left(\Delta,\Gamma,\Theta\right)\|_{F}=1}F_{\varepsilon}(\Delta,\Theta,\Gamma)$
is determined as an equilibrium point of the associated gradient system. Note,
however, that in general this may not be a global minimizer.

For the outer iteration we consider a continuous branch, as a function of
$\varepsilon$, of the minimizers
$\left(\Delta(\varepsilon),\Gamma(\varepsilon),\Theta(\varepsilon)\right)$ and
vary $\varepsilon$ iteratively in order find the smallest solution of the
scalar equation

$f(\varepsilon)=F_{\varepsilon}\left(\left(\Delta(\varepsilon),\Gamma(\varepsilon),\Theta(\varepsilon)\right)\right)=0$

with respect to $\varepsilon$.

Fig. 1: The function $f(\varepsilon)$ in a neighbourhood of
$\varepsilon^{\star}$. For $\varepsilon\geq\varepsilon^{\star}$ it becomes
identically zero.

######  Remark 5.

Note that the techniques for the distance to higher index or instability
follow directly by setting $J=0$ and not perturbing it.

###  3.2 Derivatives of eigenvalues and eigenvectors

The considered minimization is an eigenvalue optimization problem. We will
solve this problem by integrating a differential(-algebraic) equation with
trajectories that follow the gradient descent and satisfy further constraints.
To develop such a method, we first recall a classical result, see e.g. [25],
for the derivative of a simple eigenvalue and an associated eigenvector of a
matrix $C(t)$ with respect to variations in a real parameter $t$ of the
entries. Here we use the notation $\dot{C}(t):=\frac{d}{dt}C(t)$ to denote the
derivative with respect to $t$.

######  Lemma 6.

[25, Section II.1.1] Consider a continuously differentiable matrix valued
function $C(t):\mathbb{\mathbb{R}}\to\mathbb{\mathbb{R}}^{n,n}$, with $C(t)$
normal (i.e., $C(t)C(t)^{\top}=C(t)^{\top}C(t)$ for all $t$). Let $\lambda(t)$
be a simple eigenvalue of $C(t)$ for all $t$ and let $x(t)$ with $\|x(t)\|=1$
be the associated (right and left) eigenvector. Then $\lambda(t)$ is
differentiable with

$\dot{\lambda}(t)=x(t)^{H}\dot{C}(t)x(t).$ (12)

For $A\in\mathrm{Sym}^{n,n}$ consider a perturbation matrix
$\varepsilon\Delta(t)\in\mathrm{Sym}^{n,n}$ that depends on a real parameter
$t$. By Lemma 6, for a simple eigenvalue $\lambda(t)\in\mathbb{R}$ of
$A+\varepsilon\Delta(t)$ with associated eigenvector $x(t)$, $\|x(t)\|=1$, we
have (omitting the dependence on $t$)

$\frac{1}{2}\frac{d}{dt}\lambda^{2}=\varepsilon\,\lambda\,x^{\top}\dot{\Delta}x=\varepsilon\,\lambda\,\langle
xx^{\top},\dot{\Delta}\rangle.$ (13)

Similarly, which is needed in the case that $n$ is even, for all $t$, if ${\bf
i}\mu(t)\in{\bf i}\mathbb{R}$ ($\mu(t)\geq 0$) is a simple eigenvalue of a
matrix-valued function $B+\varepsilon\Theta(t)\in\mathrm{Skew}^{n,n}$, with
associated eigenvector $w(t)$, $\|w(t)\|=1$, then we have

$\frac{1}{2}\frac{d}{dt}|\mu|^{2}=\varepsilon\,\mu\,\langle{\bf
i}ww^{H},\dot{\Theta}\rangle=-\varepsilon\,\mu\,\langle{\mbox{\rm
Im}}\left(ww^{H}\right),\dot{\Theta}\rangle.$ (14)

To derive the gradient system associated with our optimization problem, we
make use of the following definition.

######  Definition 7.

Let $M\in{\mathbb{C}}^{n,n}$ be a singular matrix with a simple zero
eigenvalue. The _group inverse_(reduced resolvent) of $M$ is the unique matrix
$G$ satisfying

$MG=GM,\qquad GMG=G,\quad\mbox{and}\quad MGM=M.$

It is well-known, see [35], that for a singular and normal matrix
$M\in{\mathbb{C}}^{n,n}$ with simple eigenvalue zero, its group inverse $G$ is
equal to the Moore-Penrose pseudoinverse $M^{+}$. We have the following Lemma.

######  Lemma 8.

[35, Theorem 2] Consider a sufficiently often differentiable matrix function

$C:\mathbb{R}\to\mathbb{C}^{n,n}.$

Let $\lambda(t)$ be a simple eigenvalue of $C(t)$ for all $t$ and let $x(t)$,
with $\|x(t)\|=1$ be the associated right eigenvector function. Moreover, let
$M(t)=C(t)-\lambda(t)I$ and let $G(t)$ be the group inverse of $M(t)$. Then
$x(t)$ satisfies the system of differential equations

$\dot{x}=x\,x^{H}G(t)\dot{M}(t)x-G(t)\dot{M}(t)x.$ (15)

Moreover, if $C(t)$ is pointwise normal, then

$\dot{x}(t)=-G(t)\dot{M}(t)x(t)$ (16)

After these preparations, in the following sections we determine the
associated gradient systems for the functionals (10), respectively (11).

##  4 Gradient flow, odd state dimension

In this section we consider the case that the state dimension is odd and
construct the gradient system optimization algorithm for the functional (10).

###  4.1 Computation of the gradient

The functional $F_{\varepsilon}^{\mathrm{od}}(\Delta,\Theta,\Gamma)$ in (10)
has several parts. Applying Lemma 6 for perturbations $\varepsilon\Delta(t)$
of $E$ and $\varepsilon\Theta(t)$ of $R$, the computation of the gradient of
the part $\frac{1}{2}\Big{(}\lambda^{2}+\nu^{2}\Big{)}$ is obtained from the
expressions

$\displaystyle\frac{1}{2}\frac{d}{dt}\lambda^{2}=\varepsilon\,\lambda\,\langle
xx^{\top},\dot{\Delta}\rangle,\quad\frac{1}{2}\frac{d}{dt}\nu^{2}=\varepsilon\,\nu\,\langle
uu^{\top},\dot{\Theta}\rangle.$

Considering orthogonal projections with respect to the Frobenius inner product
onto the matrix manifold $\mathrm{Sym}^{n,n}$, we identify the constrained
gradient directions of these terms as

$\dot{\Delta}\propto\lambda\,xx^{\top},\qquad\dot{\Theta}\propto\nu\,uu^{\top},$

respectively. (Here $\propto$ denotes proportionality.) In order to treat the
other terms, we observe that

$\frac{1}{2}\frac{d}{dt}\left(|x^{\top}u|^{2}\right)=\frac{1}{2}\frac{d}{dt}\left(x^{\top}uu^{\top}x\right)=x^{\top}uu^{\top}\dot{x}+u^{\top}xx^{\top}\dot{u},$

and thus

$\displaystyle\frac{1}{2}\frac{d}{dt}\left(1-|x^{\top}u|^{2}\right)$
$\displaystyle=$
$\displaystyle\varepsilon\left((x^{\top}u)u^{\top}G\dot{\Delta}x+(u^{\top}x)x^{\top}N\dot{\Theta}u\right)$
$\displaystyle=$
$\displaystyle\varepsilon\left(\Big{\langle}\theta\,G^{\top}ux^{\top},\dot{\Delta}\Big{\rangle}+\Big{\langle}\theta\,N^{\top}xu^{\top},\dot{\Theta}\Big{\rangle}\right),$

where $\theta=x^{\top}u$, $G$ is the pseudoinverse of
$E+\varepsilon\Delta-\lambda I$, and $N$ is the pseudoinverse of
$R+\varepsilon\Theta-\nu I$.

Since $n$ is odd, which means that (generically) $0$ is a simple eigenvalue of
$J+\varepsilon\Gamma$, for the last term of (10) we have

$\displaystyle\frac{1}{2}\frac{d}{dt}\left(1-(x^{\top}w)^{2}\right)$
$\displaystyle=$
$\displaystyle\varepsilon\left(\Big{\langle}\eta\,G^{\top}wx^{\top},\dot{\Delta}\Big{\rangle}+\Big{\langle}\eta\,P^{\top}xw^{\top},\dot{\Gamma}\Big{\rangle}\right),$

where $\eta=x^{\top}w$ and $P$ is the pseudoinverse of $J+\varepsilon\Gamma$.

###  4.2 The gradient system of ODEs for the flow in the odd case

In order to compute the steepest descent direction, we minimize the gradient
of $F_{\varepsilon}$ and collect the summands involving $\dot{\Delta}$,
$\dot{\Theta}$ and those involving $\dot{\Gamma}$. Letting

$\displaystyle p$ $\displaystyle=$ $\displaystyle\theta G^{\top}u\ +\eta
G^{\top}w,$ $\displaystyle q$ $\displaystyle=$ $\displaystyle\theta
N^{\top}x,$ (17) $\displaystyle r$ $\displaystyle=$ $\displaystyle\eta
P^{\top}x,$

we have

$\displaystyle\frac{d}{dt}F_{\varepsilon}(\Delta,\Theta,\Gamma)$
$\displaystyle=$ $\displaystyle\varepsilon\,\langle\left(\lambda
x+p\right)\,x^{\top},\dot{\Delta}\rangle+\varepsilon\,\langle\left(\nu
u+q\right)\,u^{\top},\dot{\Theta}\rangle+\varepsilon\,\langle
r\,w^{\top},\dot{\Gamma}\rangle$ $\displaystyle=$
$\displaystyle\varepsilon\left(\mathrm{Sym}\left(\langle\left(\lambda
x+p\right)\,x^{\top}\right),\dot{\Delta}\rangle+\langle\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right),\dot{\Theta}\rangle+\varepsilon\,\langle\mathrm{Skew}\left(r\,w^{\top}\right),\dot{\Gamma}\rangle\right),$

where we have used the structural properties of $\dot{\Delta},\dot{\Theta}$
(symmetric) and $\dot{\Gamma}$ (skew-symmetric) and the property that for real
matrices $A$ and $B$, $\langle\mathrm{Sym}(A),\mathrm{Skew}(B)\rangle=0$.
Equation (4.2) identifies the gradient of the functional,

${\mathbf{G}}=\left(\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right),\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right),\mathrm{Skew}\left(r\,w^{\top}\right)\right):=\left({\mathbf{G}}_{E},{\mathbf{G}}_{R},{\mathbf{G}}_{J}\right).$
(19)

Since we want to impose a norm constraint on the perturbation
$(\Delta,\Theta,\Gamma)$ we need the following result.

######  Lemma 9 (Direction of steepest admissible ascent).

Let
${\mathbf{G}}=\left({\mathbf{G}}_{E},{\mathbf{G}}_{R},{\mathbf{G}}_{J}\right)\in\mathbb{R}^{n,3n}$,
${\mathbf{Z}}=\left(Z_{1},Z_{2},Z_{3}\right),{\mathbf{M}}=\left(\Delta,\Theta,\Gamma\right)\in\mathbb{R}^{n,3n}$
with ${\|(\Delta,\Theta,\Gamma)\|_{F}=1}$. A solution of the optimization
problem

$\displaystyle{\mathbf{Z}}_{\star}$ $\displaystyle=$
$\displaystyle\arg\min_{\|{\mathbf{Z}}\|_{F}=1,\,\,\langle{\mathbf{M}},{\mathbf{Z}}\rangle=0}\
\langle{\mathbf{G}},{\mathbf{Z}}\rangle$ subj. to
$\displaystyle\Delta,\Theta\in\mathrm{Sym}^{n,n},\Gamma\in\mathrm{Skew}^{n,n},$

is given by

$\displaystyle\mu{\mathbf{Z}}_{\star}$ $\displaystyle=$
$\displaystyle-{\mathbf{G}}+\varrho\,{\mathbf{M}},$ (22)
$\displaystyle\varrho$ $\displaystyle=$
$\displaystyle\left(\langle\Delta,\mathrm{Sym}\left((\lambda
x+p)\,x^{\top}\right)\rangle+\langle\Theta,\mathrm{Sym}\left((\nu
u+q)\,u^{\top}\right)\rangle+\langle\Gamma,\mathrm{Skew}\left(r\,w^{\top}\right)\rangle\right)$

where $\mu$ is the Frobenius norm of the matrix on the right-hand side. The
solution is unique if ${\mathbf{G}}$ is not a multiple of ${\mathbf{M}}$.

###### Proof.

The result follows on noting that the function to minimize is a real inner
product on $\mathbb{R}^{n,3n}$, and the real inner product with a given vector
(which here is a matrix) is minimized over a subspace by orthogonally
projecting the vector onto that subspace. The expression in (22) is the
orthogonal projection of $G$ to the tangent space at ${\mathbf{M}}$ of the
manifold of matrices of unit Frobenius norm. ∎

Taking into consideration projection with respect to the Frobenius inner
product of the vector field onto the manifolds of symmetric and skew-symmetric
matrices, this leads to the system of differential equations for the
perturbation matrices

$\displaystyle\dot{\Delta}$ $\displaystyle=$
$\displaystyle-\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)+\varrho\Delta,$ $\displaystyle\dot{\Theta}$
$\displaystyle=$ $\displaystyle-\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)+\varrho\Theta,$ (23) $\displaystyle\dot{\Gamma}$
$\displaystyle=$
$\displaystyle-\mathrm{Skew}\left(r\,w^{\top}\right)+\varrho\Gamma,$

where, for $X\in\mathbb{R}^{n,n}$, $\mathrm{Sym}(X)=\frac{X+X^{\top}}{2}$,
$\mathrm{Skew}(X)=\frac{X-X^{\top}}{2}$, and

$\varrho=\left(\langle\Delta,\mathrm{Sym}\left((\lambda
x+p)\,x^{\top}\right)\rangle+\langle\Theta,\mathrm{Sym}\left((\nu
u+q)\,u^{\top}\right)\rangle+\langle\Gamma,\mathrm{Skew}\left(r\,w^{\top}\right)\rangle\right)$

is used to ensure the norm conservation, i.e.
$\langle(\dot{\Delta},\dot{\Theta},\dot{\Gamma}),\left(\Delta,\Theta,\Gamma\right)\rangle=0$.

######  Theorem 10.

Let $(\Delta(t),\Theta(t),\Gamma(t))$ of unit Frobenius norm satisfy the
differential equation (23). If $\lambda(t)$ is a simple eigenvalue of
$A+\varepsilon\Delta(t)$, then

$\displaystyle\frac{d}{dt}F_{\varepsilon}\left(\Delta(t),\Theta(t),\Gamma(t)\right)$
$\displaystyle\leq$ $\displaystyle 0.$ (24)

###### Proof.

The result follows directly by the fact that (23) is a constrained gradient
system.

In this way we have preserved the symmetry of $E,R$ and the skew-symmetry of
$J$. It may happen however, that along the solution trajectory of (23), due to
the projection on the matrix manifolds the smallest eigenvalue $\nu$ of
$R+\varepsilon\Theta$ and/or the smallest eigenvalue $\lambda$ of
$E+\varepsilon\Delta$ become negative. In this case the perturbed system is
not a dissipative Hamiltonian system any longer. This, however, is in general
not an issue for the optimization algorithm, since the dynamical gradient
system leads to eigenvalues $\nu$, $\lambda$, with $|\nu|$ as small as
possible, for a given $\varepsilon$, and thus drives them to zero when
$\varepsilon=\varepsilon^{*}$, so that in the limiting situation also the
positive semidefiniteness of $E+\varepsilon\Delta$ and $R+\varepsilon\Theta$
holds.

###  4.3 Stationary points of (23) and low rank property

In this subsection we discuss the existence of stationary points of the
solution trajectory of (23).

######  Lemma 11.

Let $\varepsilon$ be fixed and
$F_{\varepsilon}^{\mathrm{od}}(\Delta,\Theta,\Gamma)>0$. Let $\lambda$ be a
simple eigenvalue of $E+\varepsilon\Delta$ with associate normalized
eigenvector $x$, let $\nu$ be a simple eigenvalue of $R+\varepsilon\Theta$
with associate normalized eigenvector $u$, and let $0$ be a simple eigenvalue
of $J+\varepsilon\Gamma$ with associate normalized eigenvector $w$. Then, in
the generic situation, i.e., if $\lambda,\nu\neq 0$, $\theta=x^{\top}u\neq 0$,
and $\eta=x^{\top}w\in(0,1)$, we have

$\lambda x+p\neq 0,\qquad\nu u+q\neq 0,\qquad\mbox{and}\qquad r\neq 0.$ (25)

###### Proof.

The proofs for the three cases are similar.

  * (i)

Exploiting the property that $Gx=0$, see [35], we obtain that $x^{\top}p=0$.
If we had $x^{\top}\left(\lambda x+p\right)=0$, then this would imply $\lambda
x^{\top}x=0$ and thus, since $\lambda\neq 0$, we get a contradiction, since
$\|x\|=1$.

  * (ii)

Exploiting the property $Nu=0$, we obtain that $u^{\top}q=0$. If we had
$u^{\top}\left(\nu u+q\right)=0$, then we get $\nu u^{\top}u=0$, and again we
have a contradiction.

  * (iii)

Having assumed $\eta=x^{\top}w\neq 1$ we have that $x$ and $w$ are not
aligned. As a consequence $\eta P^{\top}x\neq 0$.

Using Lemma 11, we have the following characterization of stationary points..

######  Theorem 12.

Let $(\Delta(t),\Theta(t),\Gamma(t))$ of unit Frobenius norm satisfy the
differential equation (23). Moreover, suppose that for all $t$

$F_{\varepsilon}^{\mathrm{od}}(\Delta(t),\Theta(t),\Gamma(t))>0.$

and that $0\neq\lambda(t)\in\mathbb{R}$ is a simple eigenvalue of
$E+\varepsilon\Delta(t)$ with normalized eigenvector $x(t)$, that
$0\neq\nu(t)\in\mathbb{R}$ is a simple eigenvalue of $R+\varepsilon\Theta(t)$
with associated eigenvector $u(t)$, and that $J+\varepsilon\Gamma(t)$ has a
null vector $w(t)$.

Then the following are equivalent (here we omit the argument $t$):

  * ($1$)

$\displaystyle{\frac{d}{dt}}F_{\varepsilon}^{\mathrm{od}}(\Delta,\Theta,\Gamma)=0$;

  * ($2$)

$\dot{\Delta}=0$, $\dot{\Theta}=0$, $\dot{\Gamma}=0$;

  * ($3$)

$\Delta$ is a multiple of the rank-$2$ matrix $\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)$; $\Theta$ is a multiple of the rank-$2$ matrix
$\mathrm{Sym}\left(\left(\nu u+q\right)\,u^{\top}\right)$; $\Gamma$ is a
multiple of the rank-$2$ matrix $\mathrm{Skew}\left(r\,w^{\top}\right)$ with
$p,q,r$ given by (17).

The proof follows directly by equating to zero the right hand side of (23) and
by Lemma 11, which prevents the matrices to be zero. ∎

We have also the following extremality property.

######  Theorem 13.

Consider the functional (10) and suppose that
$F_{\varepsilon}^{\mathrm{od}}(\Delta,\Theta,\Gamma)>0$. Let
$\Delta_{*}\in\mathrm{Sym}^{n,n},\Theta_{*}\in\mathrm{Sym}^{n,n}$ and
$\Gamma_{*}\in\mathrm{Skew}^{n,n}$ with
$\|\left(\Delta_{*},\Theta_{*},\Gamma_{*}\right)\|_{F}=1$. Let
$0\neq\lambda_{*}\in\mathbb{R}$ be a simple eigenvalue of
$E+\varepsilon\Delta_{*}$ with associated eigenvector $x$, let
$0\neq\nu_{*}\in\mathbb{R}$ be a simple eigenvalue of
$R+\varepsilon\Theta_{*}$ with associated eigenvector $u$, and let
$J+\varepsilon\Gamma_{*}$ have a null vector $w$. Then the following are
equivalent:

  * (i)

Every differentiable path $(\Delta(t),\Theta(t),\Gamma(t))$ (for small $t\geq
0$) with the properties that $\|(\Delta(t),\Theta(t),\Gamma(t))\|_{F}\leq 1$,
that both $\lambda(t)$ and $\nu(t)$ are simple eigenvalues of
$E+\varepsilon\Delta(t)$ and $R+\varepsilon\Theta(t)$, with associated
eigenvectors $x(t)$ and $u(t)$, respectively, and for which $w(t)$ is the null
vector of $J+\varepsilon\Gamma(t)$, so that
$\Delta(0)=\Delta_{*},\Theta(0)=\Theta_{*},\Gamma(0)=\Gamma_{*}$, satisfies

$\frac{d}{dt}F_{\varepsilon}^{\mathrm{od}}(\Delta(t),\Theta(t),\Gamma(t))\geq
0.$

  * (ii)

The matrix $\Delta_{*}$ is a multiple of the rank two matrix
$\mathrm{Sym}\left(\left(\lambda x+p\right)\,x^{\top}\right)$, $\Theta_{*}$ is
a multiple of the rank two matrix $\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)$, and $\Gamma_{*}$ is a multiple of the rank two
matrix $\mathrm{Skew}\left(r\,w^{\top}\right)$ with $p,q,r$ given by (17).

###### Proof.

Lemma 11 ensures that $\lambda x+p\neq 0$.

Assume that (i) does not hold. Then there exists a path
$(\Delta(t),\Theta(t),\Gamma(t))$ through $(\Delta_{*},\Theta_{*},\Gamma_{*})$
such that
$\frac{d}{dt}F_{\varepsilon}(\Delta(t),\Theta(t),\Gamma(t))\big{|}_{t=0}<0$.
The steepest descent gradient property shows that also the solution path of
(23) passing through $(\Delta_{*},\Theta_{*},\Gamma_{*})$ is such a path.

Hence $(\Delta_{*},\Theta_{*},\Gamma_{*})$ is not a stationary point of (23),
and Theorem 12 then yields that (ii) does not hold.

Conversely, if

$\left(\Delta_{*},\Theta_{*},\Gamma_{*}\right)\not\propto\left(\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right),\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)\mathrm{Skew}\left(r\,w^{\top}\right)\right)$

then $(\Delta_{*},\Theta_{*},\Gamma_{*})$ is not a stationary point of (23),
and Theorems 12 and 10 yield that
$\frac{d}{dt}F_{\varepsilon}(\Delta(t),\Theta(t),\Gamma(t))\big{|}_{t=0}<0$
along the solution path of (23). ∎

###  4.4 Sparsity preservation

If the matrices $E,R$ and $J$ have a given sparsity pattern, then we may
include as a constraint that the perturbations do not alter the sparsity
structure. In terms of the Frobenius norm, it is immediate to obtain the
constrained gradient system. Denoting by $\Pi_{E}$, $\Pi_{R}$, and $\Pi_{J}$,
respectively, projections onto the manifold of sparse matrices with the given
sparsity pattern and structure of $E$, $R$ and $J$, then we get

$\displaystyle\dot{\Delta}$ $\displaystyle=$
$\displaystyle-\Pi_{E}\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)+\varrho\Delta,$ $\displaystyle\dot{\Theta}$
$\displaystyle=$ $\displaystyle-\Pi_{R}\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)+\varrho\Theta,$ (26) $\displaystyle\dot{\Gamma}$
$\displaystyle=$
$\displaystyle-\Pi_{J}\mathrm{Skew}\left(r\,w^{\top}\right)+\varrho\Gamma,$

where

$\varrho=\left(\langle\Delta,\Pi_{E}\mathrm{Sym}\left((\lambda
x+p)\,x^{\top}\right)\rangle+\langle\Theta,\Pi_{R}\mathrm{Sym}\left((\nu
u+q)\,u^{\top}\right)\rangle+\langle\Gamma,\Pi_{J}\mathrm{Skew}\left(r\,w^{\top}\right)\rangle\right).$

After deriving formulas, in the next subsection we illustrate the properties
of the optimization procedure with a numerical example.

###  4.5 A numerical example

Let $n=5$ and consider the randomly generated matrices

$\displaystyle E$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rrrrr}0.15&0.02&-0.04&0.02&-0.04\\\
0.02&0.22&0&-0.01&-0.03\\\ -0.04&0&0.11&-0.07&-0.04\\\
0.02&-0.01&-0.07&0.01&0.10\\\ -0.04&-0.03&-0.04&0.10&0.39\end{array}\right],$
$\displaystyle R$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rrrrr}0.49&-0.13&0.05&-0.15&0.11\\\
-0.13&0.23&-0.05&-0.10&-0.19\\\ 0.05&-0.05&0.48&-0.06&0.02\\\
-0.15&-0.10&-0.06&0.55&0.16\\\ 0.11&-0.19&0.02&0.16&0.48\end{array}\right],$
$\displaystyle J$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rrrrr}0&-0.27&-0.03&-0.01&0.21\\\
0.27&0&-0.15&0.03&0.11\\\ 0.03&0.15&0&0.07&-0.07\\\ 0.01&-0.03&-0.07&0&0.05\\\
-0.21&-0.11&0.07&-0.05&0\end{array}\right].$

Running the two level iteration with an initial value of the functional
$F_{0}^{\mathrm{od}}\left(\cdot,\cdot,\cdot\right)=0.9181$, we find a
perturbation at a distance (rounded to four digits) $\varepsilon^{*}=0.35681$
with a common null space given by the vector

$c=\left[\begin{array}[]{rrrrr}0.2195\\\ -0.6664\\\ -0.0639\\\ 0.3187\\\
-0.6341\end{array}\right]$

and the computed perturbations are given by

$\displaystyle\Delta E$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rrrrr}-0.0385&0.0912&0.0089&0.0251&0.1408\\\
0.0912&-0.2114&-0.0218&-0.0903&-0.3482\\\
0.0089&-0.0218&-0.0009&0.0019&-0.0293\\\
0.0251&-0.0903&0.0019&0.1628&-0.0123\\\
0.1408&-0.3482&-0.0293&-0.0123&-0.4801\end{array}\right],$
$\displaystyle\Delta R$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rrrrr}-0.0689&0.1166&0.0148&-0.1044&0.1242\\\
0.1166&-0.1164&-0.0307&0.1766&-0.1433\\\
0.0148&-0.0307&-0.0049&0.0208&-0.0319\\\
-0.1044&0.1766&0.0208&-0.1592&0.1884\\\
0.1242&-0.1433&-0.0319&0.1884&-0.1679\end{array}\right],$ $\displaystyle\Delta
J$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rrrrr}0&0.0887&-0.0118&-0.0474&0.0852\\\
-0.0887&0&0.0496&0.0003&0.0027\\\ 0.0118&-0.0496&0&0.0257&-0.0488\\\
0.0474&-0.0003&-0.0257&0&-0.0030\\\
-0.0852&-0.0027&0.0488&0.0030&0\end{array}\right].$

##  5 Rank two dynamics

Theorem 12 motivates to search for a differential equation on the manifold of
rank two symmetric/skew-symmetric matrices, which still leads to a gradient
system for $F_{\varepsilon}^{\mathrm{od}}$, but in addition requires the
derivatives of the matrices $\Delta,\Theta,\Gamma$ lying in the respective
tangent spaces.

Let ${\mathcal{M}}_{2}^{n,n}:=\\{X\in\mathbb{R}^{n\times n}:{\rm
rank}(X)=2\\}$. Then we restrict the perturbations to the matrix manifolds

$\Delta,\Theta\in{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}},\qquad\Gamma\in{\mathcal{M}}_{2}^{\mathrm{Skew}},$
(33)

where
${\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}={\mathcal{M}}_{2}^{n,n}\cap\mathrm{Sym}^{n,n}$
and
${\mathcal{M}}_{2}^{\mathrm{Skew}^{n,n}}={\mathcal{M}}_{2}^{n,n}\cap\mathrm{Skew}^{n,n}$.

Following [26], every real symmetric rank two matrix $X$ of dimension $n\times
n$ can be written in the form

$X=USU^{\rm T},$ (34)

where $U\in\mathbb{R}^{n\times 2}$ has orthonormal columns, i.e., $U^{\rm
T}U=I_{2}$ and $S\in\mathrm{Sym}^{2\times 2}$. Here we will not assume that
$S$ is diagonal. Note that the representation (34) is not unique; indeed
replacing $U$ by $\widetilde{U}=UU_{1}$ with orthogonal
$U_{1}\in\mathbb{R}^{2\times 2}$ and correspondingly $S$ by
$\widetilde{S}=U_{1}^{\rm T}SU_{1}$, yields the same matrix $X=USU^{\rm
T}=\widetilde{U}\widetilde{S}\widetilde{U}^{\rm T}$.

As a compensation for the non-uniqueness in the decomposition (34), we will
use a unique decomposition in the tangent space. Let ${\mathcal{V}}_{n,2}$
denote the Stiefel manifold of real $n\times 2$ matrices with orthonormal
columns. The tangent space at $U\in{\mathcal{V}}_{n,2}$ is given by

$\displaystyle T_{U}{\mathcal{V}}_{n,2}$ $\displaystyle=$
$\displaystyle\\{\dot{U}\in\mathbb{R}^{n\times 2}:\ \dot{U}^{\rm T}U+U^{\rm
T}\dot{U}=0\\}=\\{\dot{U}\in\mathbb{R}^{n\times 2}:\ U^{\rm
T}\dot{U}\,\hbox{is skew-symmetric}\\}.$

Following [26], every tangent matrix $\dot{X}\in
T_{X}{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}$ is of the form

$\dot{X}=\dot{U}SU^{\rm T}+U\dot{S}U^{\rm T}+US\dot{U}^{\rm T},$

where $\dot{S}\in\mathrm{Sym}^{2\times 2}$, $\dot{U}\in
T_{U}{\mathcal{V}}_{n,2}$, and $\dot{S},\dot{U}$ are uniquely determined by
$\dot{X}$ and $U,S$, if we impose the orthogonality condition

$U^{\rm T}\dot{U}=0.$

We note the following lemma adapted from [26].

######  Lemma 14.

The orthogonal projection onto the tangent space
$T_{X}{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}$ at $X=USU^{\rm
T}\in{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}$ is given by

$P^{\mathrm{Sym}}_{X}(Z)=Z-(I-UU^{\rm T})Z(I-UU^{\rm T})$ (35)

for $Z\in\mathrm{Sym}^{n\times n}$.

Analogous results hold for $Y\in\mathrm{Skew}^{n,n}$, with
$S\in\mathrm{Skew}^{2,2}$.

###  5.1 A differential equation for rank two matrices

To derive the differential equation in the rank two case, we replace in (23)
the right-hand sides by the orthogonal projections to
$T_{\Delta}{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}$,
$T_{\Theta}{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}$, and
$T_{\Gamma}{\mathcal{M}}_{2}^{\mathrm{Skew}^{n,n}}$, respectively, so that
solutions starting with rank-two will retain rank two for all $t$. This gives
the differential equations

$\displaystyle\dot{\Delta}$ $\displaystyle=$
$\displaystyle-P^{\mathrm{Sym}}_{\Delta}\Bigl{(}\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)+\varrho\Delta\Bigr{)},$
$\displaystyle\dot{\Theta}$ $\displaystyle=$
$\displaystyle-P^{\mathrm{Sym}}_{\Theta}\Bigl{(}\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)+\varrho\Theta\Bigr{)},$ (36)
$\displaystyle\dot{\Gamma}$ $\displaystyle=$
$\displaystyle-P^{\mathrm{Skew}}_{\Gamma}\Bigl{(}\mathrm{Skew}\left(r\,w^{\top}\right)+\varrho\Gamma\Bigr{)},$

where again $p,q$ and $r$ are defined by (17) and

$\displaystyle\varrho$ $\displaystyle=$
$\displaystyle\Big{\langle}\Delta,P^{\mathrm{Sym}}_{\Delta}\Bigl{(}\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)\Bigr{)}\Big{\rangle}+$ $\displaystyle+$
$\displaystyle\Big{\langle}\Theta,P^{\mathrm{Sym}}_{\Theta}\Bigl{(}\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)\Bigr{)}\Big{\rangle}+\Big{\langle}\Gamma,P^{\mathrm{Skew}}_{\Gamma}\Bigl{(}\mathrm{Skew}\left(r\,w^{\top}\right)\Bigr{)}\Big{\rangle}.$

Since for $X\in{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}$ and
$Z\in\mathrm{Sym}^{n,n}$, we have $P_{X}(X)=X$ and $\langle X,Z\rangle=\langle
X,P_{X}(Z)\rangle$, (and analogous properties hold for
$X\in{\mathcal{M}}_{2}^{\mathrm{Skew}^{n,n}}$ and
$Y\in{\mathcal{M}}_{2}^{\mathrm{Skew}^{n,n}}$), the system of differential
equations can be rewritten as

$\displaystyle\dot{\Delta}$ $\displaystyle=$
$\displaystyle-P^{\mathrm{Sym}}_{\Delta}\Bigl{(}\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)\Bigr{)}+\varrho\Delta,$
$\displaystyle\dot{\Theta}$ $\displaystyle=$
$\displaystyle-P^{\mathrm{Sym}}_{\Theta}\Bigl{(}\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)\Bigr{)}+\varrho\Theta,$ (37)
$\displaystyle\dot{\Gamma}$ $\displaystyle=$
$\displaystyle-P^{\mathrm{Skew}}_{\Gamma}\Bigl{(}\mathrm{Skew}\left(r\,w^{\top}\right)\Bigr{)}+\varrho\Gamma,$

with

$\displaystyle\varrho$ $\displaystyle=$
$\displaystyle\left(\left\langle\Delta,P^{\mathrm{Sym}}_{\Delta}\left(\mathrm{Sym}\left((\lambda
x+p)\,x^{\top}\right)\right)\right\rangle+\left\langle\Theta,P^{\mathrm{Sym}}_{\Theta}\left(\mathrm{Sym}\left((\nu
u+q)\,u^{\top}\right)\right)\right\rangle\ \right.$ $\displaystyle+$
$\displaystyle\left.\left\langle\Gamma,P^{\mathrm{Skew}}_{\Gamma}\left(\mathrm{Skew}\left(r\,w^{\top}\right)\right)\right\rangle\right).$

This system differs from (23) in that the free gradient terms are replaced by
their orthogonal projections on the rank two manifold of the corresponding
structure.

To obtain the differential equation in a form that uses the factors in
$X=USU^{\rm T}$ rather than the full $n\times n$ matrix $X$, we use the
following result, whose proof is similar to that given in [26, Prop. 2.1].

######  Lemma 15.

For $X=USU^{\rm T}\in{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}}$ with nonsingular
$S\in\mathrm{Sym}^{2\times 2}$ and with $U\in\mathbb{R}^{n\times 2}$ having
orthonormal columns, the equation $\dot{X}=P^{\mathrm{Sym}}_{X}(Z)$ with $Z$
symmetric is equivalent to $\dot{X}=\dot{U}SU^{\rm T}+U\dot{S}U^{\rm
T}+US\dot{U}^{\rm T},$ where

$\displaystyle\dot{S}$ $\displaystyle=$ $\displaystyle U^{\rm T}ZU,$
$\displaystyle\dot{U}$ $\displaystyle=$ $\displaystyle(I-UU^{\rm T})ZUS^{-1}.$

An analogous statement holds for $Y\in{\mathcal{M}}_{2}^{\mathrm{Skew}^{n,n}}$
and $Z$ skew-symmetric.

With the ansatz $Z={}-\lambda
xx^{\top}-\frac{1}{2}px^{\top}+\frac{1}{2}xp^{\top}$, and introducing for
$S_{1},U_{1}$, the quantities $g_{1}=U_{1}^{\top}x\in\mathbb{R}^{2}$ and
$h_{1}=U_{1}^{\top}p\in\mathbb{R}^{2}$, this yields that the differential
equation (37) for $\Delta=U_{1}S_{1}U_{1}^{T}$ is equivalent to the following
system of differential equations

$\displaystyle\dot{S}_{1}$ $\displaystyle=$ $\displaystyle{}-\lambda
g_{1}g_{1}^{\top}-\frac{1}{2}\left(g_{1}h_{1}^{\top}+h_{1}g_{1}^{\top}\right)+\varrho
S_{1},$ $\displaystyle\dot{U}_{1}$ $\displaystyle=$
$\displaystyle\left({}-\lambda
xg_{1}^{\top}-\frac{1}{2}\left(xh_{1}^{\top}+pg_{1}^{\top}\right)+U_{1}\,\left(\lambda
g_{1}g_{1}^{\top}+\frac{1}{2}\left(g_{1}h_{1}^{\top}+h_{1}g_{1}^{\top}\right)\right)\right)S_{1}^{-1}.$

Similarly, for $\Theta=U_{2}S_{2}U_{2}^{T}$, setting
$g_{2}=U_{2}^{\top}u\in\mathbb{R}^{2},h_{2}=U_{2}^{\top}q\in\mathbb{R}^{2}$,
we obtain the system of differential equations

$\displaystyle\dot{S}_{2}$ $\displaystyle=$ $\displaystyle{}-\nu
g_{2}g_{2}^{\top}-\frac{1}{2}\left(g_{2}h_{2}^{\top}+h_{2}g_{2}^{\top}\right)+\varrho
S_{2},$ $\displaystyle\dot{U}_{2}$ $\displaystyle=$ $\displaystyle\left({}-\nu
ug_{2}^{\top}-\frac{1}{2}\left(uh_{2}^{\top}+qg_{2}^{\top}\right)+U_{2}\,\left(\lambda
g_{2}g_{2}^{\top}+\frac{1}{2}\left(g_{2}h_{2}^{\top}+h_{2}g_{2}^{\top}\right)\right)\right)S_{2}^{-1}.$

Finally, for $\Gamma=U_{3}S_{3}U_{3}^{T}$ (with
$S_{3}\in\mathrm{Skew}^{2,2}$), setting
$g_{3}=U_{3}^{\top}w\in\mathbb{R}^{2},h_{3}=U_{3}^{\top}r\in\mathbb{R}^{2}$,
we obtain the system of differential equations

$\displaystyle\dot{S}_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left({}-g_{3}h_{3}^{\top}+h_{3}g_{3}^{\top}\right)+\varrho
S_{3},$ $\displaystyle\dot{U}_{3}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{2}\left({}-wh_{3}^{\top}+rg_{3}^{\top}\right)+U_{3}\,\left(\frac{1}{2}\left(g_{3}h_{3}^{\top}-h_{3}g_{3}^{\top}\right)\right)\right)S_{3}^{-1}.$

Having established differential equations for rank two factors, in the next
section we discuss the monotonicity of the functional.

###  5.2 Monotonicity of the functional

We have the following monotonicity result, which establishes that (37) is a
gradient system for
$F_{\varepsilon}^{\mathrm{od}}(\Delta(t),\Theta(t),\Gamma(t))$.

######  Theorem 16.

Let
$\Delta(t),\Theta(t)\in{\mathcal{M}}_{2}^{\mathrm{Sym}^{n,n}},\Gamma(t)\in{\mathcal{M}}_{2}^{\mathrm{Skew}^{n,n}}$
satisfy the differential equation (37), and suppose that

$F_{\varepsilon}^{\mathrm{od}}(\Delta(t),\Theta(t),\Gamma(t))>0.$

If $\lambda(t)$ is a simple eigenvalue of $E+\varepsilon\Delta(t)$, $\mu(t)$
is a simple eigenvalue of $R+\varepsilon\Theta(t)$, and $0$ a simple
eigenvalue of $J+\varepsilon\Gamma$, then

$\frac{d}{dt}F_{\varepsilon}^{\mathrm{od}}(\Delta(t),\Theta(t),\Gamma(t))<0.$
(41)

###### Proof.

We note that

$\displaystyle\frac{1}{2\varepsilon}\frac{d}{dt}F_{\varepsilon}^{\mathrm{od}}(\Delta,\Theta,\Gamma)$
$\displaystyle=$ $\displaystyle\lambda\,\langle
xx^{\top},\dot{\Delta}\rangle+\nu\,\langle
uu^{\top},\dot{\Theta}\rangle+\left(\Big{\langle}px^{\top},\dot{\Delta}\Big{\rangle}+\Big{\langle}qu^{\top},\dot{\Theta}\Big{\rangle}+\Big{\langle}rw^{\top},\dot{\Gamma}\Big{\rangle}\right)$
$\displaystyle=$ $\displaystyle{}-\Big{\langle}\left(\lambda
x+p\right)\,x^{\top},-P^{\mathrm{Sym}}_{\Delta}\Bigl{(}\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)\Bigr{)}\Big{\rangle}+\varrho\langle\left(\lambda
x+p\right)\,x^{\top},\Delta\rangle$ $\displaystyle{}-\Big{\langle}\left(\nu
u+q\right)\,u^{\top},-P^{\mathrm{Sym}}_{\Theta}\Bigl{(}\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)\Bigr{)}\Big{\rangle}+\varrho\langle\left(\nu
u+q\right)\,u^{\top},\Theta\rangle$
$\displaystyle{}-\Big{\langle}r\,w^{\top},-P^{\mathrm{Skew}}_{\Gamma}\Bigl{(}\mathrm{Skew}\left(r\,w^{\top}\right)\Bigr{)}\Big{\rangle}+\varrho\langle
r\,w^{\top},\Theta\rangle$ $\displaystyle=$
$\displaystyle\varrho^{2}-\Big{\|}P^{\mathrm{Sym}}_{\Delta}\Bigl{(}\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)\Bigr{)}\Big{\|}^{2}-\Big{\|}P^{\mathrm{Sym}}_{\Theta}\Bigl{(}\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)\Bigr{)}\Big{\|}^{2}$ $\displaystyle-$
$\displaystyle\Big{\|}P^{\mathrm{Skew}}_{\Gamma}\Bigl{(}\mathrm{Skew}\left(r\,w^{\top}\right)\Bigr{)}\Big{\|}^{2},$

with

$\displaystyle\varrho$ $\displaystyle=$
$\displaystyle\Big{\langle}\Delta,P^{\mathrm{Sym}}_{\Delta}\Bigl{(}\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)\Bigr{)}\Big{\rangle}+\Big{\langle}\Theta,P^{\mathrm{Sym}}_{\Theta}\Bigl{(}\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)\Bigr{)}\Big{\rangle}$ $\displaystyle+$
$\displaystyle\Big{\langle}\Gamma,P^{\mathrm{Skew}}_{\Gamma}\Bigl{(}\mathrm{Skew}\left(r\,w^{\top}\right)\Bigr{)}\Big{\rangle}.$

Applying the Cauchy-Schwarz inequality proves the assertion. ∎

###  5.3 Computational approach

We will use the explicit Euler method to carry out the numerical integration
of the gradient systems associated with the rank two perturbations. For this
we use Algorithm 1, an adaptation of the method proposed in [11], and solve
instead the system given by
$\dot{K}(t)=F\left(t,K(t)\,U_{0}^{\top}\right)\,U_{0}$ on the interval
$[t_{0},t_{1}]$ and to approximate the solution at time $t_{1}$ of the ODE
$\dot{S}(t)=U_{1}^{\top}\,F\left(t,U_{1}S(t)U_{1}^{\top}\right)$. If $n$ is
large, then the memory requirement and the computing time are significantly
reduced with respect to the integration of the full ODEs.

Data: Matrix $X_{0}=U_{0}S_{0}U_{0}^{\top}$, $F(t,X)$, $t_{0},t_{1}$, $h$

Result: Matrix $X_{1}=U_{1}S_{1}U_{1}^{\top}$

begin

1 Solve the $n\times 2$ ODE
$\dot{K}(t)=F\left(t,K(t)\,U_{0}^{\top}\right)\,U_{0},\ K(t_{0})=U_{0}S_{0}$.

2 Compute a $QR$-decomposition $K(t_{1})=U_{1}\,R$.

3 Integrate the $2\times 2$ ODE
$\dot{S}(t)=U_{1}^{\top}\,F\left(t,U_{1}S(t)U_{1}^{\top}\right)$,

with initial value
$S(t_{0})=U_{1}^{\top}X_{0}U_{1}=\left(U_{1}^{\top}U_{0}\right)X_{0}\left(U_{1}^{\top}U_{0}\right)^{\top}$.

4 Set $S_{1}=S(t_{1})/\|S(t_{1})\|$,(normalization).

5 Return $S_{1},U_{1}$

Algorithm 1 Low rank symmetry/skew-symmetry preserving integrator

Using the example of Section 4.5, integrating Equations (LABEL:ode-su1),
(LABEL:ode-su2), and (LABEL:ode-su3), we obtain the same distance
$\varepsilon^{*}$ and a common null vector $c$ of the same accuracy as when
integrating (23), i.e.

$\varepsilon^{*}S_{1}=\left[\begin{array}[]{rr}-0.2722&0\\\
0&0.0719\end{array}\right],\
\varepsilon^{*}S_{2}=\left[\begin{array}[]{rr}-0.2053&0\\\
0&0.0205\end{array}\right],\
\varepsilon^{*}S_{3}=\left[\begin{array}[]{rr}0&0.0529\\\
-0.0529&0\end{array}\right]$

and

$U_{1}=\left[\begin{array}[]{rr}-0.2291&0.0737\\\ 0.5521&-0.3270\\\
0.0449&0.0453\\\ 0.0682&0.9264\\\ 0.7975&0.1658\end{array}\right],\
U_{2}=\left[\begin{array}[]{rr}-0.3685&-0.3829\\\ 0.4904&-0.5782\\\
0.0689&0.1841\\\ -0.5685&-0.5870\\\ 0.5439&-0.3750\end{array}\right],\
U_{3}=\left[\begin{array}[]{rr}-0.8360&0.2633\\\ -0.1780&-0.6387\\\
0.4867&-0.0431\\\ 0.1176&0.3495\\\ -0.1367&-0.6315\end{array}\right].$

##  6 Gradient flow for even state dimension

The derivation of the gradients in the case that the space dimension is even
is more complicated, since in this case the skew-symmetric matrix $J$ is not
guaranteed to have a zero eigenvalue.

###  6.1 Computation of the gradient of the functional (11)

Similarly to the odd case we have

$\displaystyle\frac{1}{2}\frac{d}{dt}\mu^{2}$ $\displaystyle=$
$\displaystyle{}-\varepsilon\,\mu\,\langle{\mbox{\rm
Im}}\left(ww^{H}\right),\dot{\Gamma}\rangle=\varepsilon\,\mu\,\langle\left({\mbox{\rm
Re}}(w)\,{\mbox{\rm Im}}(w)^{\top}-{\mbox{\rm Im}}(w)\,{\mbox{\rm
Re}}(w)^{\top}\right),\dot{\Gamma}\rangle.$

Considering orthogonal projections with respect to the Frobenius inner product
onto the respective matrix manifolds $\mathrm{Sym}^{n,n}$,
$\mathrm{Skew}^{n,n}$, we identify the constrained gradient directions of the
terms associated to eigenvalues as

$\dot{\Delta}\propto\lambda\,xx^{\top},\qquad\dot{\Theta}\propto\nu\,uu^{\top},\qquad\dot{\Gamma}\propto\mu\,\left({\mbox{\rm
Re}}(w)\,{\mbox{\rm Im}}(w)^{\top}-{\mbox{\rm Im}}(w)\,{\mbox{\rm
Re}}(w)^{\top}\right).$ (42)

Different to the odd case we have to consider

$\displaystyle\frac{1}{2}\frac{d}{dt}\left(|x^{\top}{\mbox{\rm
Re}}(w)|^{2}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{d}{dt}\left(x^{\top}{\mbox{\rm
Re}}(w){\mbox{\rm Re}}(w)^{\top}x\right)$ $\displaystyle=$ $\displaystyle
x^{\top}{\mbox{\rm Re}}(w){\mbox{\rm Re}}(w)^{\top}\dot{x}+{\mbox{\rm
Re}}(w)^{\top}xx^{\top}{\mbox{\rm Re}}(\dot{w}),$

and thus

$\displaystyle\frac{1}{2}\frac{d}{dt}\left({}-|x^{\top}{\mbox{\rm
Re}}(w)|^{2}\right)$ $\displaystyle=$
$\displaystyle\varepsilon\left((x^{\top}{\mbox{\rm Re}}(w)){\mbox{\rm
Re}}(w)^{\top}G\dot{\Delta}x+({\mbox{\rm
Re}}(w)^{\top}x)x^{\top}\left(\frac{P\dot{\Gamma}w+\overline{P}\dot{\Gamma}\overline{w}}{2}\right)\right)$
(43) $\displaystyle=$
$\displaystyle\varepsilon\left(\Big{\langle}\eta\,G^{\top}{\mbox{\rm
Re}}(w)x^{\top},\dot{\Delta}\Big{\rangle}+\Big{\langle}\eta\,{\mbox{\rm
Re}}(P^{H}xw^{H}),\dot{\Gamma}\Big{\rangle}\right),$

where $\eta=x^{\top}{\mbox{\rm Re}}(w)$, and $P$ is the pseudoinverse of
$J+\varepsilon\Gamma-{\bf i}\mu I$. Analogously

$\displaystyle\frac{1}{2}\frac{d}{dt}\left({}-|x^{\top}{\mbox{\rm
Im}}(w)|^{2}\right)$ $\displaystyle=$
$\displaystyle\varepsilon\left(\Big{\langle}\zeta\,G^{\top}{\mbox{\rm
Im}}(w)x^{\top},\dot{\Delta}\Big{\rangle}+\Big{\langle}\zeta\,{\mbox{\rm
Im}}(P^{H}xw^{H}),\dot{\Gamma}\Big{\rangle}\right),$

where $\zeta=x^{\top}{\mbox{\rm Im}}(w)$. Introduce $w_{\rm r}:={\mbox{\rm
Re}}(w)$ and $w_{\rm i}:={\mbox{\rm Im}}(w)$. In order to compute the steepest
descent direction, we minimize the gradient of $F_{\varepsilon}$ and collect
the summands involving $\dot{\Delta}$, $\dot{\Theta}$ and those involving
$\dot{\Gamma}$. This yields

$\frac{d}{dt}F_{\varepsilon}(\Delta,\Theta)=\varepsilon\,\langle\left(\lambda
x+p\right)\,x^{\top},\dot{\Delta}\rangle+\varepsilon\,\langle\left(\nu
u+q\right)\,u^{\top},\dot{\Theta}\rangle+\varepsilon\,\langle
W+\eta\,{\mbox{\rm Re}}(H)+\zeta\,{\mbox{\rm Im}}(H),\dot{\Gamma}\rangle,$

with

$\displaystyle p$ $\displaystyle=$ $\displaystyle\theta G^{\top}u\
+G^{\top}\left(\eta w_{\rm r}+\zeta w_{\rm i}\right),$ $\displaystyle q$
$\displaystyle=$ $\displaystyle\theta N^{\top}x,$ $\displaystyle W$
$\displaystyle=$ $\displaystyle w_{\rm r}w_{\rm i}^{\top}-w_{\rm i}w_{\rm
r}^{\top},$ $\displaystyle H$ $\displaystyle=$ $\displaystyle P^{H}xw^{H}.$

Taking into consideration projection with respect to the Frobenius inner
product of the vector field onto the manifolds of symmetric and skew-symmetric
matrices, this leads to the system of differential equations,

$\displaystyle\dot{\Delta}$ $\displaystyle=$
$\displaystyle-\mathrm{Sym}\left(\left(\lambda
x+p\right)\,x^{\top}\right)+\varrho\Delta,$ $\displaystyle\dot{\Theta}$
$\displaystyle=$ $\displaystyle-\mathrm{Sym}\left(\left(\nu
u+q\right)\,u^{\top}\right)+\varrho\Theta,$ (44) $\displaystyle\dot{\Gamma}$
$\displaystyle=$ $\displaystyle-\mathrm{Skew}\left(W+\eta\,{\mbox{\rm
Re}}(H)+\zeta\,{\mbox{\rm Im}}(H)\right)+\varrho\Gamma,$

where

$\displaystyle\varrho$ $\displaystyle=$
$\displaystyle\left(\langle\Delta,\mathrm{Sym}\left((\lambda
x+p)\,x^{\top}\right)\rangle+\langle\Theta,\mathrm{Sym}\left((\nu
u+q)\,u^{\top}\right)\rangle\right.$ $\displaystyle+$
$\displaystyle\left.\langle\Gamma,\mathrm{Skew}\left(W+\eta\,{\mbox{\rm
Re}}(H)+\zeta\,{\mbox{\rm Im}}(H)\right)\rangle\right)$

is again used to ensure the norm conservation. In this way we have again
obtained a structured flow with matrices in $\mathrm{Sym}^{n,n}$ and
$\mathrm{Skew}^{n,n}$, respectively.

######  Remark 17.

Similarly to the odd case it is possible to derive a rank two gradient system
and obtain a more effective numerical integration.

##  7 A unifying functional

One may also try to construct a unifying function that treats the odd and even
dimension case together. For this we denote by $x$ the eigenvector associated
to $\lambda$, the smallest eigenvalue of $E+\varepsilon\Delta$, with the goal
to make this the common null vector in the end. Introduce the alternative
functional

$\displaystyle\widetilde{F}_{\varepsilon}(\Delta,\Theta,\Gamma)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\Big{(}\lambda^{2}+\nu^{2}+\|\left(R+\varepsilon\Theta\right)x\|_{2}^{2}+\|\left(J+\varepsilon\Gamma\right)x\|_{2}^{2}\Big{)}$
(45) $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{(}\lambda^{2}+\nu^{2}+x^{\top}\left(R+\varepsilon\Theta\right)^{2}x-x^{\top}\left(J+\varepsilon\Gamma\right)^{2}x\Big{)},$

with $\|\left(\Delta,\Gamma,\Theta\right)\|_{F}=1$.

We observe that

$\displaystyle\frac{1}{2\varepsilon}\frac{d}{dt}\left(x^{\top}\left(R+\varepsilon\Theta\right)^{2}x\right)$
$\displaystyle=$
$\displaystyle-x^{\top}\left(R+\varepsilon\Theta\right)^{2}G\dot{\Delta}x+x^{\top}\left(R+\varepsilon\Theta\right)\dot{\Theta}x$
$\displaystyle=$ $\displaystyle{}-\langle
G\left(R+\varepsilon\Theta\right)^{2}xx^{\top},\dot{\Delta}\rangle+\langle\left(R+\varepsilon\Theta\right)xx^{\top},\dot{\Theta}\rangle,$

and similarly,

$\displaystyle\frac{1}{2\varepsilon}\frac{d}{dt}\left(x^{\top}\left(J+\varepsilon\Gamma\right)^{2}x\right)$
$\displaystyle=$
$\displaystyle{}-x^{\top}\left(J+\varepsilon\Gamma\right)^{2}G\dot{\Delta}x+x^{\top}\left(J+\varepsilon\Gamma\right)\dot{\Gamma}x$
$\displaystyle=$ $\displaystyle{}-\langle
G\left(J+\varepsilon\Gamma\right)^{2}xx^{\top},\dot{\Delta}\rangle+\langle\left(J+\varepsilon\Gamma\right)xx^{\top},\dot{\Gamma}\rangle.$

This leads to the system of ODEs

$\displaystyle\dot{\Delta}$ $\displaystyle=$
$\displaystyle-\mathrm{Sym}\left(\left(\lambda
x+s\right)\,x^{\top}\right)+\varrho\Delta,$ $\displaystyle\dot{\Theta}$
$\displaystyle=$ $\displaystyle-\mathrm{Sym}\left(\nu
u\,u^{\top}+t\,x^{\top}\right)+\varrho\Theta,$ (46)
$\displaystyle\dot{\Gamma}$ $\displaystyle=$
$\displaystyle\mathrm{Skew}\left(z\,x^{\top}\right)+\varrho\Gamma,$

where

$\varrho=\left(\langle\Delta,\mathrm{Sym}\left(\left(\lambda
x+s\right)\,x^{\top}\right)\rangle+\langle\Theta,\left(\nu
u\,u^{\top}+\mathrm{Sym}\left(t\,x^{\top}\right)\right)\rangle-\langle\Gamma,\mathrm{Skew}\left(z\,x^{\top}\right)\rangle\right)$

and

$\displaystyle s$ $\displaystyle=$ $\displaystyle
G\Big{(}\left(J+\varepsilon\Gamma\right)^{2}-\left(R+\varepsilon\Theta\right)^{2}\Big{)}x,$
$\displaystyle t$ $\displaystyle=$
$\displaystyle\left(R+\varepsilon\Theta\right)x,$ (47) $\displaystyle z$
$\displaystyle=$ $\displaystyle\left(J+\varepsilon\Gamma\right)x.$

Although this functional appears simpler to manage and does not require the
computation of two pseudo-inverses, our experiments seem to indicate that with
the previously considered functionals a higher accuracy can be reached.

##  8 The outer iteration for $\varepsilon$

In this section we discuss the outer iteration to compute the optimizing
$\varepsilon$. The simplest way to do this is by means of a bisection
technique.

Data: Matrices $E,R,J$, $k_{\max}$ (max number of iterations), $\delta$, and
tolerance tol

$\varepsilon_{0}$, $\varepsilon_{\rm lb}$ and $\varepsilon_{\rm ub}$ (starting
values for the lower and upper bounds for $\varepsilon^{\ast}$)

Result: $\varepsilon_{\delta}$ (upper bound for the distance),
$\Delta(\varepsilon^{\ast}),\Theta(\varepsilon^{\ast}),\Gamma(\varepsilon^{\ast})$

begin

1 Compute
$\Delta(\varepsilon_{0}),\Theta(\varepsilon_{0}),\Gamma(\varepsilon_{0})$

2 Compute $g(\varepsilon_{0})$

3 Set $k=0$

while  _$k\leq 1$ or $|\varepsilon_{\rm ub}-\varepsilon_{\rm lb}|>{\rm tol}$_
do

4 if  _$f({\varepsilon_{k}}) <{\rm tol}$_ then

Set $\varepsilon_{\rm ub}=\min(\varepsilon_{\rm ub},\varepsilon_{k})$

else

Set $\varepsilon_{\rm lb}=\max(\varepsilon_{\rm lb},\varepsilon_{k})$

5Compute $\varepsilon_{k+1}=(\varepsilon_{\rm lb}+\varepsilon_{\rm ub})/2$
(bisection step)

if  _$k=k_{\max}$_ then

Return interval $[\varepsilon_{\rm lb},\varepsilon_{\rm ub}]$

Halt

else

Set $k=k+1$

6Compute
$\Delta(\varepsilon_{k}),\Theta(\varepsilon_{k}),\Gamma(\varepsilon_{k})$

7 Compute $f(\varepsilon_{k})$

8Return $\varepsilon^{\ast}=\varepsilon_{k}$

Algorithm 2 Bisection method for distance approximation

###  8.1 An illustrative example

To illustrate the performance of the described algorithm, we consider a
scalable linear mass-spring-damper system, that has been used as a model
reduction test case in [23].

Fig. 2: The function $f(\varepsilon)$ for the mass-spring-damper example.

This test case generates matrices $M=M^{T}\geq 0,G=-G^{T},D=D^{T}\geq
0,K=K^{T}>0$ and, after first order formulation, leads to a dH pencil

$\lambda\underbrace{\begin{bmatrix}K&0&0\\\ 0&M&0\\\
0&0&0\end{bmatrix}}_{E}-\left(\underbrace{\begin{bmatrix}0&K&0\\\
-K&0&-G^{T}\\\ 0&G&0\end{bmatrix}}_{J}-\underbrace{\begin{bmatrix}0&0&0\\\
0&D&0\\\ 0&0&0\end{bmatrix}}_{R}\right)$

This pencil is regular and of index two. If one puts $\gamma I$ in the
$(3,3)$-block of the matrix $R$, then the distance to index $2$ and
instability is $\gamma$. Choosing the dimension $N=100$ we obtain matrices
$E,R,J\in\mathbb{R}^{3N+1}$.

We fix for $\gamma=10^{-1}$. The plot of the function $f(\varepsilon)$
obtained by integrating (37) for increasing $\varepsilon$, is given in Figure
2.

## Conclusions and further work

We have investigated a structured distance problem related to the study of
port-Hamiltonian systems, that is determining the closest triplet of matrices
to a given one, sharing a common null-space.

## Acknowledgments

The first author acknowledges that his research was supported by funds from
the Italian MUR (Ministero dell’Università e della Ricerca) within the PRIN
2017 Project “Discontinuous dynamical systems: theory, numerics and
applications” and by the INdAM Research group GNCS (Gruppo Nazionale di
Calcolo Scientifico).

The second author thanks Deutsche Forschungsgemeinschaft (DFG) for support
within the project B03 in CRC TRR 154.

## References

  * [1] F. Achleitner, A. Arnold, and V. Mehrmann. Hypocoercivity and controllability in linear semi-dissipative ODEs and DAEs. http://arxiv.org/abs/2104.07619 Submitted for publication, 2021.
  * [2] N. Aliyev, P. Benner, E. Mengi, P. Schwerdtner, and M. Voigt.  Large-scale computation of $L_{\infty}$-norms by a greedy subspace method. SIAM J. Matrix Anal. Appl., 38, 1496–1516, 2017.
  * [3] N. Aliyev, V. Mehrmann, and E. Mengi.  Computation of stability radii for large-scale dissipative Hamiltonian systems.  Advances Comp. Math., 46:6, 2020.
  * [4] R. Altmann, V. Mehrmann, and B. Unger. Port-Hamiltonian formulations of poroelastic network models. http://arxiv.org/abs/2012.01949. Submitted for publication, 2020.
  * [5] C. Beattie, V. Mehrmann, H. Xu, and H. Zwart.  Port-Hamiltonian descriptor systems.  Math. Control, Signals, Sys., 30:17, 2018.  https://doi.org/10.1007/s00498-018-0223-3.
  * [6] C. Beattie, V. Mehrmann, and P. Van Dooren.  Robust port-Hamiltonian representations of passive systems.  Automatica, 100:182–186, 2019.
  * [7] T. Berger, H. Gernandt, C. Trunk, H. Winkler, and M. Wojtylak.  A new bound for the distance to singularity of a regular matrix pencil.  In Proc. Appl. Mathematics and Mechanics, volume 17.1, pages 863–864. Wiley Online Library, 2017.
  * [8] T. Berger, H. Gernandt, C. Trunk, H. Winkler, and M. Wojtylak.  The gap distance to the set of singular matrix pencils.  Linear Alg. Appl., 564:28–57, 2019.
  * [9] K. E. Brenan, S. L. Campbell, and L. R. Petzold.  Numerical Solution of Initial-Value Problems in Differential Algebraic Equations.  SIAM Publications, Philadelphia, PA, 2nd edition, 1996.
  * [10] R. Byers, C. He, and V. Mehrmann.  Where is the nearest non-regular pencil.  Linear Algebra Appl., 285:81–105, 1998.
  * [11] G. Ceruti and C. Lubich.  Time integration of symmetric and anti-symmetric low-rank matrices and Tucker tensors.  BIT Numerical Mathematics, 60:591–614, 2020.
  * [12] L. Dai.  Singular Control Systems.  Number 118 in Lecture Notes in Control and Information Sciences. Springer-Verlag, Berlin, 1989.
  * [13] N.H. Du, V.H. Linh and V. Mehrmann,  Robust stability of differential-algebraic equations. In _Differential Algebraic Equation Forum._ Surveys in Differential-Algebraic Equations I, A. Ilchmann and T. Reis Edtrs. pp. 63–96, 2013.
  * [14] H. Egger, T. Kugler, B. Liljegren-Sailer, N. Marheineke, and V. Mehrmann,  On structure preserving model reduction for damped wave propagation in transport networks.  SIAM J. Sci. Comp., 40:A331–A365, 2018.
  * [15] R. W. Freund.  Structure-preserving model order reduction of rcl circuit equations.  In Model Order Reduction: Theory, Research Aspects and Applications, pages 49–73. Springer, 2008.
  * [16] R. W. Freund.  The SPRIM algorithm for structure-preserving order reduction of general rcl circuits.  In Model reduction for circuit simulation, pages 25–52. Springer, 2011.
  * [17] F. R. Gantmacher.  Theory of Matrices, volume 1.  Chelsea, New York, 1959.
  * [18] N. Gillis, V. Mehrmann, and P. Sharma.  Computing nearest stable matrix pairs.  Numer. Lin. Alg. Appl., 25:e2153, 2018.
  * [19] N. Gillis and P. Sharma.  On computing the distance to stability for matrices using linear dissipative hamiltonian systems.  Automatica, 85:113–121, 2017.
  * [20] N. Gillis and P. Sharma.  Finding the nearest positive-real system.  SIAM J. Matrix Anal. Appl., 56(2):1022–1047, 2018.
  * [21] N. Gräbner, V. Mehrmann, S. Quraishi, C. Schröder, and U. von Wagner.  Numerical methods for parametric model reduction in the simulation of disc brake squeal.  Z. Angew. Math. Mech., 96:1388–1405, 2016.
  * [22] N. Guglielmi, C. Lubich, and Volker Mehrmann.  On the nearest singular matrix pencil.  SIAM J. Matrix Anal. Appl., 38(3):776–806, 2017.
  * [23] S-A. Hauschild, N. Marheineke, and Volker Mehrmann.  Model reduction techniques for port‐Hamiltonian differential‐algebraic systems.  Control and Cybernetics, 48(1):1–19, 2019.
  * [24] B. Jacob and H. Zwart.  Linear port-Hamiltonian systems on infinite-dimensional spaces.  Operator Theory: Advances and Applications, 223. Birkhäuser/Springer Basel AG, Basel CH, 2012.
  * [25] T. Kato.  _Perturbation theory for linear operators_. Springer-Verlag, 1995.
  * [26] O. Koch and Ch. Lubich.  Dynamical low-rank approximation.  SIAM J. Matrix Anal. Appl., 29(2):434–454, 2007.
  * [27] P. Kunkel and V. Mehrmann.  Differential-Algebraic Equations. Analysis and Numerical Solution.  EMS Publishing House, Zürich, Switzerland, 2006.
  * [28] C. Mehl, V. Mehrmann, and P. Sharma.  Stability radii for linear hamiltonian systems with dissipation under structure-preserving perturbations.  SIAM J. Matrix Anal. Appl., 37:1625–1654, 2016.
  * [29] C. Mehl, V. Mehrmann, and P. Sharma.  Stability radii for real linear Hamiltonian systems with perturbed dissipation.  BIT Numerical Mathematics, 57:811–843, 2017.
  * [30] C. Mehl, V. Mehrmann, and M. Wojtylak.  On the distance to singularity via low rank perturbations.  Operators and Matrices, 9:733–772, 2015.
  * [31] C. Mehl, V. Mehrmann, and M. Wojtylak.  Linear algebra properties of dissipative Hamiltonian descriptor systems.  SIAM J. Matrix Anal. Appl., 39(3):1489–1519, 2018.
  * [32] C. Mehl, V. Mehrmann, and M. Wojtylak,  Distance problems for dissipative Hamiltonian systems and related matrix polynomials. http://arxiv.org/abs/2001.08902.  Linear Alg. Appl., https://doi.org/10.1016/j.laa.2020.05.026, 2020.
  * [33] V. Mehrmann and R. Morandin.  Structure-preserving discretization for port-hamiltonian descriptor systems.  In 58th IEEE Conf. Decision and Control (CDC), Nice, pages 6863–6868, 2019.  https://arXiv:1903.10451.
  * [34] V. Mehrmann and P. Van Dooren,  Optimal robustness of port-Hamiltonian systems,  SIAM J. Matrix Anal. Appl., 41:134–151, 2020.
  * [35] C.D. Meyer and G.W. Stewart.  Derivatives and perturbations of eigenvectors.  SIAM J. Numer. Anal., 25:679–691, 1988.
  * [36] M. L. Overton.  HANSO: Hybrid Algorithm for Non-Smooth Optimization.  http://www.cs.nyu.edu/overton/software/hanso. Date 23.11.20
  * [37] A. J. van der Schaft.  Port-Hamiltonian differential-algebraic systems.  In Surveys in Differential-Algebraic Equations I, pages 173–226. Springer-Verlag, 2013.
  * [38] A. J. van der Schaft and B. M. Maschke.  Hamiltonian formulation of distributed-parameter systems with boundary energy flow.  J. Geom. Phys., 42:166–194, 2002.
  * [39] A. van der Schaft and B. Maschke.  Generalized port-Hamiltonian dae systems.  Systems & Control Letters, 121:31–37, 2018.
  * [40] A. J. van der Schaft and D. Jeltsema.  Port-Hamiltonian systems theory: An introductory overview.  Found. and Trends in Systems and Control, 1(2-3):173–378, 2014\.
  * [41] W.C. Waterhouse.  The codimension of singular matrix pairs.  Linear Algebra Appl., 57:227–245, 1984.
  * [42] M.A. Woodbury  Inverting modified matrices.  Statistical Research Group, Memo. Rep. no. 42,  Princeton University, Princeton, N.J., 4 pp, 1950.