Stabilising the Metzler matrices with applications to dynamical systems

TL;DR

This paper develops efficient algorithms to find the closest stable Metzler matrices to unstable ones, using spectral methods, with applications in dynamical systems and large-scale matrices up to dimension 2000.

Contribution

It introduces a generalized spectral optimization approach for stabilizing Metzler matrices, extending existing methods to larger matrices and diverse applications.

Findings

Explicit solutions and algorithms for closest stable/unstable Metzler matrices.

Demonstrated efficiency on matrices up to dimension 2000.

Applications to dynamical systems and switching systems.

Abstract

Metzler matrices play a crucial role in positive linear dynamical systems. Finding the closest stable Metzler matrix to an unstable one (and vice versa) is an important issue with many applications. The stability considered here is in the sense of Hurwitz, and the distance between matrices is measured in $l_\infty,\ l_1$, and in the max norms. We provide either explicit solutions or efficient algorithms for obtaining the closest (un)stable matrix. The procedure for finding the closest stable Metzler matrix is based on the recently introduced selective greedy spectral method for optimizing the Perron eigenvalue. Originally intended for non-negative matrices, here is generalized to Metzler matrices. The efficiency of the new algorithms is demonstrated in examples and by numerical experiments in the dimension of up to 2000. Applications to dynamical systems, linear switching systems, and sign-matrices are considered.

Tables

Table 1. Table 1a: Average number of iterations for maximization, for full Metzler families

$N\setminus d$	25	100	500	2000
50	3	3	3.1	3
100	3.2	3.1	3.1	3.1
250	3.2	3.2	3.1	3.1

Table 2. Table 1b: Average computing time for maximization, for full Metzler families

$N\setminus d$	25	100	500	2000
50	0.02s	0.09s	0.5s	2.77s
100	0.04s	0.19s	1.02s	6.5s
250	0.1s	0.44s	2.19s	136.77s

Table 3. Table 2a: Average number of iterations for minimization, for full Metzler families

$N\setminus d$	25	100	500	2000
50	3.1	3.1	3.1	3.1
100	3.1	3	3	3.1
250	3.3	3.2	3.1	3.1

Table 4. Table 2b: Average computing time for minimization, for full Metzler families

$N\setminus d$	25	100	500	2000
50	0.02s	0.09s	0.54s	2.5s
100	0.04s	0.19s	1s	4.6s
250	0.1s	0.45s	2.32s	129.35s

Table 5. Table 3a: Average number of iterations for maximization, for sparse Metzler families

$N\setminus d$	25	100	500	2000
50	5.4	4.4	4.2	4
100	5.7	4.8	4.2	4.2
250	6.1	4.8	4.2	4.2

Table 6. Table 3b: Average computing time for maximization, for sparse Metzler families

$N\setminus d$	25	100	500	2000
50	0.04s	0.13s	0.64s	3.16s
100	0.07s	0.26s	1.3s	6.03s
250	0.19s	0.69s	3.13s	151.39s

Table 7. Table 4a: Average number of iterations for minimization, for sparse Metzler families

$N\setminus d$	25	100	500	2000
50	7.2	5	4.1	4.1
100	6.9	5.2	4.3	4.1
250	12.2	5.4	4.6	4.2

Table 8. Table 4b: Average computing time for minimization, for sparse Metzler families

$N\setminus d$	25	100	500	2000
50	0.87s	0.14s	0.66s	3.15s
100	0.14s	0.28s	1.29s	6.9s
250	0.46s	0.85s	3.35s	169.89s

Table 9. Table 5a: Effects of sparsity, number of iterations

${\gamma }_i$	9-15	16-21	22-51	52-76
MAX	4	4	4	3.9
MIN	4.2	4	4.2	3.6

Table 10. Table 5b: Effects of sparsity, computing time

${\gamma }_i$	9-15	16-21	22-51	52-76
MAX	2.54s	2.5s	2.5s	2.49s
MIN	2.76s	2.66s	2.96s	2.33s

Table 11. Table 4: Results for full Metzler matrices

$d$	50	100	250	500	750	1000
$\#$	4.9	7.6	7.3	8.3	9.9	7.4
$t$	0.19s	1.08s	3.06s	20.78s	80.98s	115.03s

Table 12. Table 5: Results for sparse Metzler matrices, with 9-15% sparsity of each row

$d$	50	100	250	500	750	1000
$\#$	7.1	9.6	12.6	13.4	17.1	15.4
$t$	5.08s	19.5s	63.01s	152.54s	484.06s	773.18s

Table 13. Table 6: Effects of sparsity

${\gamma }_i$	3-8	9-15	16-21	22-51	52-76
$\#$	15.3	15.9	10.4	21	10.2
$t$	1503.58s	664.77s	427.3s	381.17s	140.63s

Table 14. Table 7: Results for totally positive matrices

$d$	50	100	250	500	750	1000
$\#$	8.4	10.6	13.5	22.8	21.7	23.6
$t$	0.45s	1.61s	5.32s	29.38s	77.69s	179.7s

Table 15. Table 8: Results for sparse non-negative matrices, with 9-15% sparsity of each row

$d$	50	100	250	500	750	1000
$\#$	12.7	13.1	12.5	17.5	17.1	21.7
$t$	0.97s	1.32s	4.84s	21.2s	53.25s	129.62s

Table 16. Table 9: Effects of sparsity, for non-negative matrices

${\gamma }_i$	3-8	9-15	16-21	22-51	52-76
$\#$	17.8	19.7	16.5	25.2	25.1
$t$	105.22s	84.25s	130.35s	154.02s	233.16s

Equations

$$\begin{array}[]{lcc}\|X\|_{\max}&=&\max\limits_{1\leqslant i,j\leqslant d}\ |x_{ij}|;\vspace*{0.05cm}\\ \|X\|_{1}&=&\max\limits_{1\leqslant j\leqslant d}\ \sum_{i=1}^{d}|x_{ij}|;\vspace*{0.05cm}\\ \|X\|_{\infty}&=&\max\limits_{1\leqslant i\leqslant d}\ \sum_{j=1}^{d}|x_{ij}|.\end{array}$$

$$\rho (A)=\max \{|\lambda ||\lambda \mathrm{is}\mathrm{an}\mathrm{eigenvalue}\mathrm{of}A\in \mathcal{M}\}.$$

$$\eta (A)=\max \{\mathrm{Re}(\lambda )|\lambda \mathrm{is}\mathrm{an}\mathrm{eigenvalue}\mathrm{of}A\in \mathcal{M}\}.$$

$$\rho (A)=\eta (A).$$

$$\rho (A+hI)=\eta (A)+h,$$

$$a^{\prime}_{ij}=\left\{\begin{array}[]{cc}a_{ij},&a_{ij}\geqslant 0\ {\rm or}\ i=j\\ 0,&{\rm otherwise}\end{array}\right.$$

$$\left\{\begin{array}[]{l}\|X-A\|_{\max}\ \rightarrow\min\\ \eta(X)=0,\ X\in{\mathcal{M}}.\end{array}\right.$$

$$\eta (A+\alpha H)<0,0⩽\alpha <\frac{1}{\rho (-A^{-1}H)}$$

$$W(\alpha )=-A-\alpha H=-A-\alpha (-AB)=-A(I-\alpha B).\text{vspace}\ast 0.1cm$$

$$X<A+\frac{H}{\rho (-A^{-1}H)}$$

$$X<A-\frac{E}{⟨A^{-1}\mathbf{e},\mathbf{e}⟩}$$

$${\tau }_{\ast }=\frac{1}{\sum _{i,j}|a_{ij}^{\prime }|},$$

$$\rho (-A^{-1}E)=-⟨A^{-1}\mathbf{e},\mathbf{e}⟩.$$

$$\underset{X\in \mathcal{M}}{\min }\{\eta (X)|\Vert X-A{\Vert }_{\max }⩽\tau \}.$$

$$a_{ij}(\tau)=\left\{\begin{array}[]{l}a_{ii}-\tau,\ i=j\\ \max\{0,a_{ij}-\tau\},\ i\neq j.\end{array}\right.$$

$$\underset{X\in \mathcal{M}}{\min }\Vert X-A{\Vert }_{\max }.$$

$$H=\frac{A({\tau }_1)-A({\tau }_2)}{{\tau }_2-{\tau }_1}.$$

$${\tau }_{\ast }={\tau }_2-\frac{1}{\rho (-A^{-1}({\tau }_2)H)}.$$

$$\eta\Big{(}A\Big{(}\max_{ij}a_{ij}\Big{)}\Big{)}<0.$$

$$\{0\}\cup \{a_{ij}|a_{ij}>0\}$$

$$A(\tau )=A({\tau }_2)+\frac{{\tau }_2-\tau }{{\tau }_2-{\tau }_1}(A({\tau }_1)-A({\tau }_2))=A({\tau }_2)+({\tau }_2-\tau )H.$$

$${\mathcal{M}}_i=\{\mathbf{a}\in {\mathbb{R}}^d|a_{j\ne i}⩾0\},i=1,\dots ,d.$$

$$\left\{\begin{array}[]{l}\langle A^{(i)},{\boldsymbol{v}}_{k}\rangle\ \to\ \max\\ A^{(i)}\in{\mathcal{F}}_{i}\end{array}\right.$$

$$A=\left(\begin{array}[]{rrr}-2&\phantom{-}2&\phantom{-}0\\ \phantom{-}0&-6&\phantom{-}5\\ \phantom{-}2&\phantom{-}2&-9\\ \end{array}\right).$$

$${\mathcal{B}}_{\epsilon }(A)=\{X\in \mathcal{M}|\Vert A-X\Vert ⩽\epsilon \}$$

$${\mathcal{B}}_{\epsilon }^{(i)}=\{X^{(i)}\in {\mathcal{M}}_i|\Vert A^{(i)}-X^{(i)}\Vert ⩽\epsilon \},i=1,\dots ,d.$$

$$\left\{\begin{array}[]{l}\|X-A\|\ \rightarrow\min\\ \rho(X)=1.\end{array}\right.$$

$$X=A+{\tau }_{\ast }{E}_k.$$

$$\left\{\begin{array}[]{l}\|X-A\|\ \rightarrow\min\\ \rho(X)=h.\end{array}\right.$$

$$X=A+{\tau }_{\ast }{E}_k.$$

Full text

Keywords: Metzler matrix, spectral abscissa, spectral radius, Perron eigenvalue, sign-matrix, dynamical system, Hurwitz stability

AMS 2010 subject classification: 15A42, 15B35, 15B48, 90C26, 93C30

1. Introduction

Metzler matrices and the problem of finding their closest Hurwitz (un)stable counterparts arise in linear dynamical systems, differential equation analysis, electrodynamics, population dynamics, economics, etc. [3, 4, 5, 10, 22]. We refer to the problems of finding the closest (un)stable matrix as the (de)stabilization problems. By the closest matrix we mean a matrix $X$ which minimizes the distance from the starting matrix $A$, in some matrix norm. In this paper we are concerned with $l_{\infty},\ l_{1}$ and $\max$ norms. In other norms, such as Euclidean or Frobenius norms, the stabilization problem is hard and as a rule does not allow finding global solutions [10, 13].

The analogous problems for the non-negative matrices and Schur stability were considered in [11, 12, 13, 14]. In [14] an explicit solutions were given for the problems of destabilization. For the stabilization, an efficient algorithms were developed. Of special interest is the algorithm for finding the closest stable non-negative matrix in $l_{\infty}$ and $l_{1}$ norms, based on strikingly efficient greedy spectral simplex method. This method was devised in the same paper, and it is used for minimizing the spectral radius on the product families of (sparse) non-negative matrices. The greedy method can also be of independent interest, apart from the stabilization problem. Its modification, selective greedy method, was recently developed in [15]. It is significantly easier to implement and has a quadratic convergence.

In order to stabilize the Metzler matrix, we need to deal with the problem of optimizing spectral abscissa on families of Metzler families. In general, this problem is notoriously hard, since the objective function is neither convex nor concave, nor is Lipschitz. Because of this, there might be many points of local extrema which are, moreover, hard to identify. However, for matrix families with special structure, i.e. product (or uncertainty) structure, it is possible to find an efficient algorithm for optimizing the spectral abscissa. In this paper we present such algorithm, based on selective greedy method.

Furthermore, we discuss the application of stabilization problem for finding the closest stable Metzler sign-matrix. Metzler sign-matrices and sign-stability are interesting concepts with many applications [16, 20, 23, 24, 26], and here we present the algorithm for solving the sign-stabilization problem. Building upon this, we propose a procedure for stabilizing a positive linear switching system (LSS) [27, 28] having a certain structure.

The paper is structured as follows:

In Section 2 we deal with the Hurwitz (de)stabilization in the $\max$-norm, on the set of Metzler matrices. An explicit solution for the destabilization is given, while for the stabilization an iterative procedure is provided.

Section 3 generalizes both the greedy and selective greedy methods for optimizing Perron eigenvalue on Metzler product families. They remain as effective as when used on non-negative matrices. This is demonstrated by the numerical results: even for the dimension of 2000, the computational time does not exceed 3 minutes (on a standard laptop).

Section 4 is concerned with Hurwitz (de)stabilization in the $l_{\infty}$ norm111 See the beginning of the next subsection.. Again, an explicit solution for the closest unstable Metzler matrix is obtained. As for the stabilization, we provide an iterative algorithm. With slight modification, this algorithm can be applied for Schur stabilization of non-negative matrices. This modification significantly improves the computing time for Schur stabilization from [14].

The applications for Hurwitz stabilization (in $l_{\infty}$ norm) are presented in Section 5.

1.1. Some preliminaries

Throughout this paper we assume, unless otherwise is stated, that all our matrices are $d\times d$ square. The matrix norms considered are $l_{\infty},\ l_{1}$ and max norms, given by [2]:

[TABLE]

As it can be easily seen $\|X\|_{\infty}=\|X^{T}\|_{1}$. Consequently, all the results for the $l_{\infty}$ norm apply for $l_{1}$ norm as well: we just need to take the transpose of the matrix. Having this in mind, we shall only develop results for $l_{\infty}$ and max norms222 Norms without a subscript will be regarded as an $l_{\infty}$ norms (i.e $\|\cdot\|\ \equiv\|\cdot\|_{\infty}$).. We now define the following sets of matrices:

For two matrices $A$ and $B$ we write $A\geqslant B$ if $A-B\geqslant 0$, and $A\leqslant B$ if $B-A\geqslant 0$. Analogously, we define strict relations $A>B$ and $A<B$.

An apparent property of a strict Metzler matrix, which we will amply exploit, is that it can be translated to a non-negative matrix: i.e. if $A$ is a strict Metzler matrix, then there exists $h>0$ such that $A+hI$ is non-negative. Here (and throughout the paper) $I$ denotes the identity matrix.

The spectral radius and spectral abscissa of a Metzler matrix are fundamental concepts in the development of our work, which we define below:

[TABLE]

[TABLE]

From the Perron-Frobenius theory we know that a non-negative matrix has its spectral radius as a leading eigenvalue [1]. Since this implies that the leading eigenvalue is real and non-negative, we have:

[TABLE]

$\Box$**

Specially, a strictly positive matrix has a simple leading eigenvalue, and its corresponding leading eigenvector is unique and strictly positive. For the non-negative matrix, this might not always be the case: its leading eigenvalue might not be simple, its leading eigenvector can contain zero entries and not necessarily be unique.

Given the fact that any strict Metzler matrix $A$ can be translated to a non-negative one $A+hI$, for some $h>0$, and that both $A$ and $A+hI$ have the same sets of eigenvectors, from Proposition 1 we obtain:

[TABLE]

for any $h>0$ such that $A+hI$ is non-negative. $\Box$

It is well known that a spectral radius is monotone on the set of non-negative matrices [14]. Now, we prove that the same holds for the spectral abscissa on the set of Metzler matrices, that is:

Proof. Let $h>0$ be such that both $A+hI$ and $B+hI$ are non-negative. Since $B$ is entrywise bigger than $A$ we have $B+hI\geqslant A+hI$. By the monotonicity of spectral radius $\rho(B+hI)\geqslant\rho(A+hI)$ holds. Using Lemma 1, we obtain $\eta(B)\geqslant\eta(A)$. $\Box$

Even though the results of the Perron-Frobenius theory regarding the spectral radius do not hold for the strict Metlzer matrices, the following Theorem 1 reveals that, in this case, we should just change our focus to the spectral abscissa.

We denote by ${\mathcal{H}}$ the set of all weakly Hurwitz stable Metzler matrices, and by ${\mathcal{H}}_{s}$ the set of all strongly Hurwitz stable Metzler matrices.

Remark. When searching for the closest (un)stable Metzler matrix to the matrix $A$, the starting matrix need not necessarily be Metzler. It is easy to check that both the real matrix $A$ and matrix $A^{\prime}$, with entries

[TABLE]

have the same solution. Therefore, without the loss of generality, we shall always assume that our starting matrix is Metzler.

2. Problems of (de)stabilization in the max-norm

2.1. Closest Hurwitz unstable Metzler matrix in the max-norm

For a given Hurwitz stable Metzler matrix $A$ we consider the problem of finding its closest (weakly) Hurwitz unstable matrix $X$, with respect to the max-norm. In other words, for a matrix $A\in{\mathcal{H}}_{s}$ find a matrix $X$ that satisfies:

[TABLE]

The following set of results provides the answer to the posed problem. We start from

and its analogue for Metzler matrices:

[TABLE]

and \eta\Big{(}A+\frac{H}{\rho(-A^{-1}H)}\Big{)}=0.

Proof. Define matrix $B=-A^{-1}H$. Since $A$ is Hurwitz stable, from Lemma 4 we have $B\geqslant 0$. Define

[TABLE]

From Lemma 3 and Lemma 4, $W(\alpha)$ is invertible for every \alpha\in\Big{[}0,\frac{1}{\rho(B)}\Big{)}. Moreover, for every $\alpha$ from the given interval $W^{-1}(\alpha)$ is non-negative. Since $-W(\alpha)=A+\alpha H$ is Metzler, by Lemma 4, it is Hurwitz stable, so we have (1).

For the proof of the second part, we assume that matrix $B$ is strictly positive. Then its leading eigenvector ${\boldsymbol{v}}$ is also strictly positive, and we have -W\Big{(}\frac{1}{\rho(B)}\Big{)}{\boldsymbol{v}}=0. Hence, by the continuity of spectral abscissa on the set of Metzler matrices, we have \eta\Big{(}A+\frac{H}{\rho(B)}\Big{)}=0. $\Box$

The following result is a direct consequence of Lemma 5:

[TABLE]

are stable. $\Box$

Applying Corollary 1 to the special case of the non-negative matrix $H=E$, where $E$ is matrix of all ones, we arrive to the explicit formula for Hurwitz destabilization in the max-norm.

[TABLE]

are Hurwitz stable. Moreover, matrix $A-\frac{E}{\langle A^{-1}{\boldsymbol{e}},{\boldsymbol{e}}\rangle}$ is the closest Hurwitz unstable Metzler matrix with the distance to the starting matrix $A$ equal to

[TABLE]

where $a^{\prime}_{ij}$ are entries of the matrix $A^{-1}$.

Proof. The proof stems directly from the fact that

[TABLE]

$\Box$

2.2. Closest Hurwitz stable Metzler matrix in the max-norm

We now deal with a problem of Hurwitz stabilization in the max-norm. First, consider the following spectral abscissa minimization problem for a given matrix $A\in{\mathcal{M}}$ and parameter $\tau\geqslant 0$:

[TABLE]

[TABLE]

Proof. Clearly, a Metzler matrix $A(\tau)$ is feasible for (2). Moreover, for any other feasible solution $X$, we have $X\geqslant A(\tau)$. From the monotonicity of spectral abscissa we have that $A(\tau)$ is truly the optimal solution to our problem. $\Box$

Now, for a given Hurwitz unstable Metzler matrix $A$, consider the following minimization problem:

[TABLE]

Denote by $\tau_{*}$ its optimal value.

Proof. The statement follows directly from Lemma 6 and the fact the function $\eta(A(\tau))$ is monotonically decreasing in $\tau$. $\Box$

The matrix $A(\tau_{*})$ is the closest Hurwitz stable Metzler matrix to a given matrix $A\in{\mathcal{M}}$ with $\eta(A)>0$, and $\tau_{*}$ is the minimal distance. We now present the algorithm for computing the value $\tau_{*}$ and solving (3).

Algorithm 1: Finding the closest Hurwitz stable Metzler matrix in max-norm

Step 1. Sort all positive entries of matrix $A$ in an increasing order. Check if the highest entry is on the main diagonal. If so, take $\tau_{*}=\max_{i}a_{ii}$ and finish the procedure. Else, continue to the

Step 2. Using the monotonicity of the spectral abscissa, find by bisection in the entry number the value $\tau_{1}$, which is the largest between all positive entries $a_{ij}$ and zero, such that $\eta(A(a_{ij}))>0$. Find value $\tau_{2}$, which is the smallest entry of $A$ with $\eta(A(a_{ij}))<0$.

Step 3. Form the matrix

[TABLE]

Step 4. Compute

[TABLE]

$\diamond$

Proof. First, assume that the highest entry is on the main diagonal. Then, the matrix $A(\max_{i}a_{ii})$ will be non-positive diagonal matrix such that $\eta(A(\max_{i}a_{ii}))=0$. Since $\tau_{*}$, by Lemma 7, is the unique root of the equation $\eta(A(\tau))=0$, we can put $\tau_{*}=\max_{i}a_{ii}$.

Now, assume that the highest entry is off-diagonal. Then, the matrix $A(\max_{ij}a_{ij})$ is negative diagonal. Therefore, we have

[TABLE]

On the other hand, $\eta(A(0))=\eta(A)>0$. So, it is possible to find two values $\tau_{1}<\tau_{2}$ from the set

[TABLE]

such that $\eta(A(\tau_{1}))>0$ and $\eta(A(\tau_{2}))<0$. $A(\tau)$ is linear in $\tau$, and for every $\tau\in[\tau_{1},\tau_{2}]$ we have

[TABLE]

Applying Lemma 5 for Metzler matrix $A(\tau_{2})$ and non-negative matrix $H$ concludes the proof. $\Box$

3. Optimizing the spectral abscissa on product families of Metzler matrices

Stabilization in $l_{\infty}$ norm is much more delicate than in the max-norm. In order to properly address it, we first need to find a way to optimize spectral abscissa on Metzler product families333 For the reminder of the text, each time we mention Metzler product families or positive product families, or similar, it will be in the context of product families of Metzler matrices, product families of positive matrices, etc.. Optimizing spectral abscissa on such families can also be of independent interest, since product families occur quite often in applications [7, 8, 9, 19, 22, 29].

By ${\mathcal{M}}_{i}$ we denote the sets of vectors

[TABLE]

The matrices belonging to the product families are constructed by independently taking $i$th row from the uncertainty set ${\mathcal{F}}_{i}$. Moreover, product families can be topologically seen as product sets of their uncertainty sets: ${\mathcal{F}}={\mathcal{F}}_{1}\times\cdots\times{\mathcal{F}}_{d}$.

3.1. The greedy spectral method for the spectral abscissa

In [9] spectral simplex method for optimizing a spectral radius over non negative product families was presented. There, it was also observed that the theory behind this procedure can be generalized for the optimization of the spectral abscissa over Metzler product families.

[14] presents the greedy spectral simplex method (or shorter greedy method) for optimizing the spectral radius over positive product families, together with the modifications for the minimization on non-negative product families. Numerical experiments showed that this procedure is strikingly more efficient than the spectral simplex method. This efficiency was theoretically confirmed in [15], where it was proved that greedy method has local quadratic convergence. Here we present a version of this algorithm for maximizing spectral abscissa over Metzler product families.

Algorithm 2: Maximizing spectral abscissa over Metzler product families

Initialization. Let ${\mathcal{F}}={\mathcal{F}}_{1}\times\cdots\times{\mathcal{F}}_{d}$ be a product family of Metzler matrices. Taking arbitrary $A^{(i)}\in{\mathcal{F}}_{i}\,,\ i=1,\ldots,d$, form a Metzler matrix $A_{1}\in{\mathcal{F}}$ with rows $A^{(1)},\ldots,A^{(d)}$. Take its leading eigenvector ${\boldsymbol{v}}_{1}$. Denote $A_{1}^{(i)}\,=\,A^{(i)},\ i=1,\ldots,d$.

Main loop; the $k$th iteration. We have a matrix $A_{k}\in{\mathcal{F}}$ composed with rows $A_{k}^{(i)}\in{\mathcal{F}}_{i},i=1,\ldots,d$. Compute its leading eigenvector ${\boldsymbol{v}}_{k}$ (if it is not unique, take any of them) and for $i=1,\ldots d$, find a solution $\hat{A}^{(i)}\in{\mathcal{F}}_{i}$ of the problem

[TABLE]

For each $i=1,\ldots,d$, do:

If $\langle\hat{A}^{(i)},{\boldsymbol{v}}_{k}\rangle=\langle A_{k}^{(i)},{\boldsymbol{v}}_{k}\rangle$, then set $A_{k+1}^{(i)}=A_{k}^{(i)}$.

Otherwise, if $\langle\hat{A}^{(i)},{\boldsymbol{v}}_{k}\rangle>\langle A_{k}^{(i)},{\boldsymbol{v}}_{k}\rangle$, set $A_{k+1}^{(i)}=\hat{A}^{(i)}$.

Form the corresponding matrix $A_{k+1}$. If the first case took place for all $i$, i.e. if $A_{k+1}=A_{k}$, go to Termination. Otherwise, go to $(k+1)$st iteration.

Termination. If ${\boldsymbol{v}}_{k}>0$, the procedure is finished; $A_{k}$ is maximal in each row444 see Definition 11. and $\eta_{\max}=\eta(A_{k})$ is a solution. If ${\boldsymbol{v}}_{k}$ has some zero components, then the family ${\mathcal{F}}$ is reducible, and we need to stop the algorithm and factorize ${\mathcal{F}}$ (see section 3.4). $\diamond$

Remark. The procedure for minimizing the spectral abscissa is exactly the same, except that we change the row if $\langle\hat{A}^{(i)},{\boldsymbol{v}}_{k}\rangle<\langle A_{k}^{(i)},{\boldsymbol{v}}_{k}\rangle$, and we omit the requirement that ${\boldsymbol{v}}_{k}>0$.

3.2. Theoretical results

For general Metzler matrices, the described greedy algorithm may cycle. Hence, we one needs to implement some additional modifications. We will address this issue later on, but for now, we will focus on full Metzler families, since in this case the cycling does not occur. As with the positive and non-negative product families, the following theoretical results are the extension of the similar results from [8].

Proof. Let $A{\boldsymbol{u}}\geqslant\lambda{\boldsymbol{u}}$. Since $A$ is Metzler, there exists $h\geqslant 0$ such that $A+hI\geqslant 0$. We have $(A+hI){\boldsymbol{u}}\geqslant(\lambda+h){\boldsymbol{u}}$, and since the analogous of this Lemma for non-negative matrices and spectral radii holds [14], we obtain $\rho(A+hI)\geqslant\lambda+h$, and therefore $\eta(A)\geqslant\lambda$ (Lemma 1). The proof of the second statement follows by the same reasoning. $\Box$

Proof. Let $A{\boldsymbol{u}}>\lambda{\boldsymbol{u}}$. The statement $\eta(A)<\lambda$ is in direct contradiction with Lemma 8 and $\eta(A)=\lambda$ would imply that $\lambda$ is an eigenvalue of matrix $A$, which is not true. Hence, $\eta(A)>\lambda$ has to hold. The proof of the second statement follows the same reasoning. $\Box$

The following Proposition stems directly from Lemma 8.

1. if ${\boldsymbol{v}}>0$ and $A$ is maximal in each row with the respect to ${\boldsymbol{v}}$, then $\eta(A)=\max\limits_{X\in{\mathcal{F}}}\eta(X)$.

2. if $A$ is minimal in each row with the respect to ${\boldsymbol{v}}$, then $\eta(A)=\min\limits_{X\in{\mathcal{F}}}\eta(X)$. $\Box$*

Proof. Maximization. Let $A$ be maximal in each row w.t.r. to ${\boldsymbol{v}}$. Since for a full Metzler matrix its leading eigenvector is strictly positive, we can apply $1)$ of the previous Proposition, and therefore the maximal spectral abscissa is indeed achieved for matrix $A$.

Conversely, assume that $A$ maximizes the spectral abscissa on ${\mathcal{F}}$, but it is not maximal w.t.r. to its leading eigenvector ${\boldsymbol{v}}$ in each row. Construct matrix $A^{\prime}\in{\mathcal{F}}$ by changing all the necessary rows of the matrix $A$, so that $A^{\prime}$ becomes maximal w.t.r. to ${\boldsymbol{v}}$ in each row. Thus we have $A^{\prime}{\boldsymbol{v}}>A{\boldsymbol{v}}=\eta(A){\boldsymbol{v}}$. This, in turn, results in $\eta(A^{\prime})>\eta(A)$ (Corollary 2), which is incorrect. So, $A$ has to be maximal in each row w.t.r. to its eigenvector.

Minimization. Suppose now that $A$ is minimal in each row w.t.r. to ${\boldsymbol{v}}$. One direction of the equivalence is basically $2)$ of Proposition 2. The proof for the other direction goes the similar way as with the maximization, taking into account the strict positivity of the vector ${\boldsymbol{v}}$. $\Box$

Proof. Let $\varepsilon^{\prime}>0$ be such that the shifted product family ${\mathcal{F}}_{\varepsilon^{\prime}}={\mathcal{F}}+\varepsilon^{\prime}E$ is full Metzler family, and let $A_{\varepsilon^{\prime}}\in{\mathcal{F}}_{\varepsilon^{\prime}}$ be the matrix with maximal spectral radius. Since this matrix is full, by Proposition 3 it has to be maximal in each row. For all $0<\varepsilon\leqslant\varepsilon^{\prime}$, we associate one such a matrix $A_{\varepsilon}$. By compactness, there exists a sequence $\{\varepsilon_{k}\}_{k\in\mathbb{N}}$ such that $\varepsilon_{k}\rightarrow 0$ as $k\rightarrow+\infty$. Hence, matrices $A_{\varepsilon_{k}}$ converges to a matrix $A\in{\mathcal{F}}$ and their respective leading eigenvectors converge to a nonzero vector ${\boldsymbol{v}}$. By continuity, ${\boldsymbol{v}}$ is a leading eigenvector of $A$, and $A$ is maximal in each row w.t.r. to ${\boldsymbol{v}}$. $\Box$

Proof. Assume we use the greedy method for maximizing spectral abscissa (the same goes for minimization). Proposition 4 guarantees the existence of a matrix maximal in each row, and Proposition 3 that this matrix is the optimal one. Starting from matrix $A_{1}\in{\mathcal{F}}$ and iterating trough the algorithm we obtain the sequence of matrices $A_{2},A_{3},\ldots$ with the sequence of their respective spectral abscissas $\eta(A_{2}),\eta(A_{3}),\ldots$, which is increasing. Moreover, each $i$th row of of the given matrices is some vertex of the polyhedron ${\mathcal{F}}_{i}$, so the number of total states is finite. Therefore, we arrive at our solution in finite number of iterations, increasing the spectral abscissa in each one, until we reach the optimal matrix. $\Box$

3.3. Selective greedy method for spectral abscissa

To resolve the issue of cycling in [15] the selective greedy method was proposed. This method is as efficient as the greedy method, and it does not cycle even if the matrices are very sparse.

Notice that in this definition we have extended the notion of the selected leading eigenvector from [15] to include Metzler matrices as well, which will be fully justified further in the text.

Using the selective greedy method on Metzler matrices to avoid cycling in the sparse case seems like a logical move. However, it may happen that the power method applied to a strict Metzler matrix, starting with the initial vector ${\boldsymbol{v}}_{0}={\boldsymbol{e}}$, does not converge.

Example. We apply the power method to the matrix

[TABLE]

Starting from the vector ${\boldsymbol{v}}_{0}={\boldsymbol{e}}$, and obtain the following sequence of vectors:

${\boldsymbol{v}}_{1}=(0.242,\ 0,\ -0.97)\\ {\boldsymbol{v}}_{2}=(-0.025,\ -0.507,\ 0.862)\\ {\boldsymbol{v}}_{3}=(-0.094,\ -0.649,\ -0.755)\\ {\boldsymbol{v}}_{4}=(0.138,\ -0.694,\ 0.707)\\ \cdots$

It is clear that the sequence of vectors ${\boldsymbol{v}}_{k}$ is not going to converge, because the sign pattern keeps changing with every iteration.

Having in mind that a Metzler matrix can be translated to a corresponding non-negative matrix, we propose a translative power method.

Proof. If $A$ is non-negative, then translative power method is just the regular power method and we can use Proposition 5. Suppose now that $A$ is strictly Metzler. Applying a translative power method on it, taking ${\boldsymbol{e}}$ as an initial vector, we obtain a selected leading eigenvector ${\boldsymbol{v}}$ for a non-negative matrix $\tilde{A}=A+hI$. This eigenvector ${\boldsymbol{v}}=\lim\limits_{\varepsilon\to 0}{\boldsymbol{v}}_{\varepsilon}$ is a limit of a sequence of leading eigenvectors of perturbed positive matrices $\tilde{A}_{\varepsilon}=\tilde{A}+\varepsilon E$. Since the leading eigenvectors ${\boldsymbol{v}}_{\varepsilon}$ for matrices $\tilde{A}_{\varepsilon}$ are also the leading eigenvectors for the perturbed Metzler matrices $A_{\varepsilon}=A+\varepsilon E$, we have, by the definition, that ${\boldsymbol{v}}$ is a selected leading eigenvector of the matrix $A$. $\Box$

Proposition 6 allows us to use the Definition 13 when dealing with Metzler matrices as well, having in mind that in this case we need to resort to translative power method. Moreover,

Proof. Assume we have a product family ${\mathcal{F}}$ that contains strict Metzler matrices, since the case of non-negative families is already resolved by Theorem 5. Cycling of a selective greedy method on this family would also imply the cycling of the greedy spectral simplex method on the family ${\mathcal{F}}+\varepsilon E$ of full Metzler matrices, which is, given Theorem 4, impossible. $\Box$

Remark. As discussed in [15], when implementing selective greedy method we compute the selected leading eigenvectors by applying the power method not on the non-negative matrices $X_{k}$, obtained iterating through the procedure, but on the matrices $X_{k}+I$. This is due to the fact that, for any $k$, both matrices $X_{k}$ and $X_{k}+I$ have the same selected leading eigenvector, while for the matrix $X_{k}+I$ it will be unique, and maybe multiple. The uniqueness of this eigenvector is a guarantee that the power method will converge.

For the same reasons, if we obtain a strict Metzler matrix $X$ in some iteration, we will actually compute the selected leading eigenvector by applying power method on a non-negative matrix $X+(h+1)I$, where $h=|\min_{i}x_{ii}|$.

In addition, all remarks regarding the cycling due to the computational errors and how to avoid it, given in [15], apply for the Metzler matrices as well.

3.4. Optimizing spectral abscissa for reducible families

If we finish the (selective) greedy procedure with a leading eigenvector ${\boldsymbol{v}}_{k}$ that has zero entries, we might not obtain an optimal solution. In this case the product family ${\mathcal{F}}$ is reducible. This means that the linear space spanned by the vectors with the same support as ${\boldsymbol{v}}_{k}$ is invariant with the respect to all matrices of ${\mathcal{F}}$. One way to resolve this is to use the Frobenius factorization and run the maximization procedure on irreducible blocks ${\mathcal{F}}^{(j)}$, constructing the optimal solution as the block matrix. The optimal value $\eta_{\max}$ is the largest of all optimal values $\eta_{\max}^{(j)}$ among all blocks. The algorithm for the Frobenius factorization can be found in [31].

In practice, however, the Frobenius factorisation usually takes more time than the run of the greedy algorithm for the whole family. That is why it makes sense to avoid the factorisation procedure by making the family irreducible. Since having only one irreducible matrix is enough to make the whole family irreducible, a simple strategy is to simply include a matrix $H=\alpha P-\beta I$ to the family ${\mathcal{F}}$, where $P$ is a cyclic permutation matrix, $\alpha>0$, and $\beta\geqslant 0$. The matrix $P$ is irreducible, and so the family ${\mathcal{F}}\cup\{H\}$ is irreducible as well, hence the greedy method has to finish with a positive leading eigenvector. If the final matrix has no rows from $H$, then it is optimal for the original family ${\mathcal{F}}$. Otherwise we have to choose another pair $\alpha,\beta$, setting them to be smaller. However, if trying out various pairs of $\alpha,\beta$ produces the final matrix with some rows from $H$, we need to resort to Frobenius factorization.

3.5. Numerical results

We test the selective greedy method on both full Metzler product families and sparse Metzler product families with density parameters $\gamma_{i}$ (the percentage of nonzero entries of the vectors belonging to the uncertainty set $\mathcal{F}_{i}$).

The results of tests on full Metzler product families, as dimension $d$ and the size of the uncertainty sets $N$ vary, are given in Tables 1 and 2, for maximization and minimization, respectively. In Tables 3 and 4 we report the behaviour of selective greedy algorithm on sparse Metzler families. For each uncertainty set $\mathcal{F}_{i}$ the density parameter $\gamma_{i}$ is randomly chosen from the interval $9-15\%$ and the elements of the set $\mathcal{F}_{i}$ are randomly generated in accordance to it.

Table 5 shows how the number of iterations and computing time vary as the density parameter is changed. The dimension is kept fixed at $d=600$ and the cardinality of each product set at $|\mathcal{F}_{i}|=200$, while we vary the interval from which $\gamma_{i}$ takes value.

4. Problems of (de)stabilization in $l_{\infty}$ norm

Let $A$ be a Metzler matrix, and let

[TABLE]

denote an $\varepsilon$-ball of Metzler matrices around $A$. Optimization of the spectral abscissa on ${\mathcal{B}}_{\varepsilon}(A)$ is a crucial tool for the Hurwitz (de)stabilization. The greedy method developed in Section 3 can be applied on ${\mathcal{B}}_{\varepsilon}$, since it can be considered as a product set of the balls

[TABLE]

4.1. Closest Hurwitz unstable Metzler matrix in the $l_{\infty}$ norm

Let us first consider the problem of Schur destabilization of non-negative matrix in $l_{\infty}$ norm: if $A$ is a non-negative matrix with $\rho(A)<1$, find a matrix $X$ that satisfies

[TABLE]

An explicit solution to (4) was given in [14]. An important point is that the solution is always a non-negative matrix, entrywise bigger than $A$.

[TABLE]

$\Box$**

We can also solve a more general problem: for a given $h>0$ and non-negative matrix $A$ with $\rho(A)<h$, find the closest matrix $X$ having $\rho(X)=h$, i.e., find a solution to

[TABLE]

Here, as with the $h=1$ case, the solution $X$ is also non-negative and entrywise bigger than $A$. We can also prove a generalization of Theorem 7:

[TABLE]

Proof. The optimal matrix $X_{*}$ for (5) is also a solution to the maximization problem

[TABLE]

having the optimal value $\tau=\tau_{*}$. Let us characterize this matrix for arbitrary $\tau$. From Proposition 2 and Proposition 4, applied to the non-negative product families and spectral radii, we can conclude that $X$ is maximal in each row for the product family with the uncertainty sets:

[TABLE]

Conversely, every matrix $X$ with $\rho(X)=h$ which is maximal in each row w.r.t. a strictly positive leading eigenvector, solves (5).

Any matrix $X\in{\mathcal{B}}_{\tau}^{+}(A)$ with the leading eigenvector ${\boldsymbol{v}}$ is optimal in the $i$th row if and only if the scalar product $\langle X^{(i)}-A^{(i)},{\boldsymbol{v}}\rangle$ is maximal under the constraint $\langle X^{(i)}-A^{(i)},{\boldsymbol{e}}\rangle=\tau$. This maximum is equal to $r\tau$, where $r$ is the maximal component of the leading eigenvector ${\boldsymbol{v}}$. Denote the index of this component by $k$. Then

[TABLE]

Hence, if $X$ is maximal in each row, we have

[TABLE]

Furthermore, since each set ${\mathcal{B}}_{\tau}^{+}(A^{(i)})$ contains a strictly positive point, we have

[TABLE]

Therefore, the leading eigenvector ${\boldsymbol{v}}$ is strictly positive so, by Proposition 2, the matrix $X$ maximizes the spectral radius on the product family ${\mathcal{B}}_{\tau}^{+}(A)$.

Thus, the optimal matrix has the form (6) for some $k$, and $\|X-A\|=\tau_{*}$. It remains to find $k$ for which the value of $\tau$ is minimal.

Since $\rho(A+\tau_{*}E_{k})=h$, it follows that $\tau_{*}$ is the smallest positive root of the equation

[TABLE]

Since $\rho(\frac{A}{h})<1$ we have $(hI-A)^{-1}=\frac{1}{h}(I-\frac{A}{h})^{-1}=\frac{1}{h}\sum_{j=0}^{\infty}(\frac{A}{h})^{j}\geqslant 0$. Multiplying equation (7) by $\det(-(hI-A)^{-1})$ we obtain

[TABLE]

The matrix $\tau(hI-A)^{-1}E_{k}$ has only one nonzero column. This is the $k$th column, equal to $\tau(hI-A)^{-1}{\boldsymbol{e}}$. Hence,

[TABLE]

We conclude that the minimal $\tau$ corresponds to the biggest component of this vector. $\Box$

We can now move to a Hurwitz destabilization, or more formally put: for a given Metzler matrix $A$ with $\eta(A)<0$, find a solution to

[TABLE]

Proof. First, let us show that $X$ is Metzler. Assume the contrary, that it has some negative off-diagonal entries. We consider a matrix $X^{\prime}$ with entries

[TABLE]

We have $\|X^{\prime}-A\|<\|X-A\|$ and $\eta(X^{\prime})<0$. Since $X^{\prime}$ is Metzler, there exists $h>0$ such that $X^{\prime}_{h}=X^{\prime}+hI$ is non-negative. In addition, from Lemma 1, we obtain $\rho(X^{\prime}_{h})<h$.

Now, define the matrix $X_{h}=X+hI$ which contains off-diagonal negative entries. From the inequality $\|X_{h}^{k}\|\ <\|(X^{\prime}_{h})^{k}\|$, using Gelfand’s formula for spectral radius [2], we arrive at $\rho(X_{h})<\rho(X^{\prime}_{h})$. Since the largest real part of all the eigenvalues of $X$ is equal to zero, we have

[TABLE]

Adding $h$ on both sides of inequality and using $\rho(X)+h=\rho(X+hI)$, we obtain $h<\rho(X^{\prime}_{h})$, which is a contradiction. Therefore, matrix $X$ must not contain negative off-diagonal entries, i.e. it has to be Metzler.

Now let us prove that $X\geqslant A$. Assume the converse, that there exist some entries $x_{ij}$ of $X$ such that $x_{ij}<a_{ij}$. Define the matrix $\tilde{X}$ with entries

[TABLE]

$\tilde{X}$ is Metzler, $\tilde{X}\geqslant A$ and $\tilde{X}\geqslant X$. By the monotonicity of spectral abscissa we have $\eta(\tilde{X})\geqslant\eta(X)$. $X$ is the closest Hurwitz unstable Metzler matrix, so $\|\tilde{X}-A\|\ \geqslant\|X-A\|$. This is impossible, since $\|\tilde{X}-A\|\ <\|X-A\|$ holds. Hence, $X$ has to be entrywise bigger than $A$. $\Box$

[TABLE]

Proof. Lemma 4 ensures invertibility of the matrix $A$. We also remark that a Hurwitz stable Metzler matrix is necessary a strict Metzler matrix: if it were non-negative, than the biggest real part of its eigenvalues would also be non-negative, which is not the case.

Now, let $h>0$ be such that the matrix $\tilde{A}=A+hI$ is non-negative. We have $\rho(\tilde{A})<h$. Using Theorem 8 we obtain the closest non-negative matrix to $\tilde{A}$. Denote it by $\tilde{X}$. We have $\rho(\tilde{X})=h$ and

[TABLE]

where $\tau_{*}$ is the reciprocal of the biggest component of the vector $(hI-\tilde{A})^{-1}{\boldsymbol{e}}=-A^{-1}{\boldsymbol{e}}$. Let $X=\tilde{X}-hI$. $X$ is Metzler with $\eta(X)=0$, and satisfies (9).

Now we need to check that $X$ is really the closest Hurwitz unstable Metzler matrix for $A$. Assume that $Y\neq X$ is the true solution of (8). We have $\|Y-A\|\ <\|X-A\|=\tau_{*}$. Define $\tilde{Y}=Y+hI$. $\tilde{Y}$ is entrywise bigger than $\tilde{A}$, and $\rho(\tilde{Y})=h$. Since $\tilde{X}$ is the closest non-negative matrix to $\tilde{A}$ having spectral radius equal to $h$, the following is true:

[TABLE]

which is a contradiction. Therefore, our $X$ is really the optimal solution. $\Box$

Example. For a given Hurwitz stable Metzler matrix

[TABLE]

with $\eta(A)=-1$, the biggest component of the vector $-A^{-1}{\boldsymbol{e}}$ is the third one, and it is equal to $2.5$. Thus we get the matrix

[TABLE]

as the closest Hurwitz unstable.

4.2. Closest Hurwitz stable Metzler matrix in the $l_{\infty}$ norm

The problem of Schur stabilization for non-negative matrices in $l_{\infty}$ norm was also considered in [14]: if $A$ is a given non-negative matrix with $\rho(A)>1$, find a non-negative matrix $X$ that satisfies

[TABLE]

Notice that, in contrast to the problem of Schur destabilization (4), here we have to impose the non-negativity condition on our solution.

The story is the same with the Hurwitz (de)stabilization of Metzler matrix. For Hurwitz destabilization (8), requesting for the solution to be Metzler is redundant, a fact confirmed by Lemma 9. However, this is not the case with the Hurwitz stabilization; here, we have to explicitly impose the restriction on our solution to be Metzler. So, our problem can be written in the following manner:

[TABLE]

To obtain the solution, we solve the related problem

[TABLE]

and use the strategy that implements the bisection in $\tau$, together with the greedy procedure for minimizing the spectral abscissa described in Section 3.

First, we describe how to compute matrix $X$ explicitly in each iteration of the Algorithm 2 for minimization, thus finding the optimal solution on the ball of radius $\tau$.

Fix some $\tau>0$, and let $X^{\prime}$ be a matrix obtained in some iteration with ${\boldsymbol{v}}\geqslant 0$ as its leading eigenvector. We rearrange positive entries of ${\boldsymbol{v}}$: $v_{j_{1}}\geqslant\cdots\geqslant v_{j_{m}}$, where ${\mathcal{S}}$ is support of ${\boldsymbol{v}}$, and $|{\mathcal{S}}|\ =m$. In the next iteration of the greedy procedure we construct the matrix $X=(x_{ij})$:

For each $i\in{\mathcal{S}}$ we solve

[TABLE]

which we can rewrite as

[TABLE]

We can solve (12) explicitly:

[TABLE]

where

[TABLE]

For $i\in{\mathcal{S}}$, but $j\notin{\mathcal{S}}$, we take $x_{ij}=a_{ij}$; and if $i\notin{\mathcal{S}}$ we put $X^{(i)}=X^{{}^{\prime}(i)}$.

Alg. 3: Computing the closest Hurwitz stable Metzler matrix in $l_{\infty}$ norm

Step 0. Take $\frac{\|A\|}{2}$ as the starting value for $\tau$.

Step 1. Start the selective greedy method for minimizing the spectral abscissa on the ball of Metzler matrices ${\mathcal{B}}_{\tau}(A)$ (applying (13)). Iterate until a matrix $X$ is obtained with $\eta(X)<0$. When this is done, stop the greedy procedure, compute its leading eigenvector ${\boldsymbol{v}}$, and proceed to the next step.

If greedy procedure finishes finding the matrix with non-negative minimal spectral abscissa on ${\mathcal{B}}_{\tau}(A)$, keep implementing the bisection on $\tau$. Do this until the ball ${\mathcal{B}}_{\tau}(A)$ which contains a matrix with negative spectral abscissa is obtained.

Step 2. Construct matrices $C=(c_{ij})$ and $R$, as follows. For each $i\in{\mathcal{S}}$

[TABLE]

where $l_{i}$ is given by (14), while $R$ is a boolean matrix having ones on positions $(i,l_{i})$ and zeros in all other places. If some of the indices $i,j$ are not in the support, then we put $c_{ij}=x_{ij}$. We can write $X=C-\tau R$.

Denote by $\tau_{*}$ a potential optimal value of the problem (11). To determine this value we proceed to the next step.

Step 3. Since $\eta(X)<0$, we have by Lemma 4 that $X$ is invertible and

[TABLE]

Since $0=\eta(C-\tau_{*}R)=\eta(C-\tau R+(\tau-\tau_{*})R)$, it follows that $\det(-(C-\tau R)-(\tau-\tau_{*})R)=0$. From here we have

[TABLE]

Matrix $-(C-\tau R)^{-1}R$ is non-negative and $\frac{1}{\tau-\tau_{*}}$ is its (positive) leading eigenvalue. Therefore, in this step we find the potential optimal value $\tau_{*}$ by computing the leading eigenvalue $\lambda$ of the matrix $-(C-\tau R)^{-1}R$, and then calculating

[TABLE]

Step 4. To check if $\tau_{*}$ is really optimal, start iterating through the greedy procedure on the ball ${\mathcal{B}}_{\tau_{*}}(A)$, as in Step 1.

If during some iteration a matrix with negative spectral abscissa is obtained, $\tau_{*}$ is not the optiaml value. Stop the greedy procedure and return to Step 1, taking now $\tau_{*}$ as the starting value for $\tau$.

Else, if we finish the greedy procedure obtaining the matrix $X_{*}$ with minimal spectral abscissa $\eta(X_{*})=0$ on the ball ${\mathcal{B}}_{\tau_{*}}(A)$, we are done: $\tau_{*}$ is the optimal value for the problem (11), with $X_{*}$ as the corresponding optimal solution. $\diamond$

Remark. It is worth noting that in Step 1 (and Step 4) of Algorithm 3, we do not need to bring the greedy procedure to its completion. The moment we obtain the matrix with a negative spectral abscissa (not necessarily minimal), we can leave the greedy procedure, thus saving the computational time.

Example. For a given Hurwitz unstable Metzler matrix

[TABLE]

with $\eta(A)>0$, the described procedure computes the matrix

[TABLE]

as the closest stable Metzler with $\tau_{*}=10$.

4.3. Numerical results

In this subsection we report the numerical results for the implementation of Algorithm 3, as the dimension of the starting matrix $A$ is varied. $d$ represents the dimension of the problem, $\#$ the number of iterations required for the procedure to terminate, and $t$ is a computational time.

The next table shows how the sparsity affects the computations. The dimension is kept fixed at $d=850$ in all experiments, while the density parameter $\gamma_{i}$ for each row of the starting matrix is randomly chosen from the given interval.

4.4. Closest Schur stable non-negative matrix in the $l_{\infty}$ norm

Building upon the ideas of Algorithm 3, we propose the following modification of the algorithm for Schur stabilization, given in [14].

Alg. 4: Computing the closest Schur stable non-negative matrix in $l_{\infty}$ norm

Step 0. Take $\frac{\|A\|}{2}$ as the starting value for $\tau$.

Step 1. Start the selective greedy method procedure for minimizing the spectral radius on the ball of non-negative matrices ${\mathcal{B}}_{\tau}^{+}(A)$. Iterate until a matrix $X$ is obtained with $\rho(X)<1$. When this is done, stop the greedy procedure, compute its leading eigenvector ${\boldsymbol{v}}$, and proceed to the next step.

If greedy procedure finishes finding the Schur unstable matrix as the optimal one on ${\mathcal{B}}_{\tau}^{+}(A)$, keep implementing the bisection in $\tau$, until the ball ${\mathcal{B}}_{\tau}^{+}(A)$ containing strongly Schur stable matrix is reached.

Step 2. We construct matrices $C=(c_{ij})$ and $R$, as follows. For each $i\in{\mathcal{S}}$

[TABLE]

where $l_{i}$ the minimal index for which $\sum_{s=1}^{l_{i}}a_{ij_{s}}>\tau$. If this $l_{i}$ does not exist, we put $l_{i}=m$. $R$ is a boolean matrix with ones on positions $(i,l_{i})$ and zeros in all other places. If some of the indices $i,j$ are not in the support, then we put $c_{ij}=x_{ij}$. We can write $X=C-\tau R$.

Denote by $\tau_{*}$ a potential optimal value of the problem (10). To determine it, proceed to the next step.

Step 3. Since $\rho(X)<1$, we have by Lemma 3 that $I-X$ is invertible and

[TABLE]

Since $1=\rho(C-\tau_{*}R)=\rho(C-\tau R+(\tau-\tau_{*})R)$, it follows that $\det(I-(C-\tau R)-(\tau-\tau_{*})R)=0$. From here we have

[TABLE]

Matrix $[I-(C-\tau R)]^{-1}R$ is non-negative and $\frac{1}{\tau-\tau_{*}}$ is its (positive) leading eigenvalue.

Therefore, in this step we find the potential optimal value $\tau_{*}$ by computing the leading eigenvalue $\lambda$ of the matrix $[I-(C-\tau R)]^{-1}R$, and then calculating

[TABLE]

Step 4. To check if $\tau_{*}$ is really optimal, start iterating through the greedy procedure on the ball ${\mathcal{B}}_{\tau_{*}}^{+}(A)$, as in the Step 1.

If in some iteration a matrix $Y$ is obtained with $\rho(Y)<1$, then $\tau_{*}$ is not optimal. Stop the greedy procedure, return to Step 1, and continue doing the bisection taking now $\tau_{*}$ as the starting value.

Else, if we finish the greedy procedure obtaining the matrix $X_{*}$ with minimal spectral radius $\rho(X_{*})=1$ on the ball ${\mathcal{B}}_{\tau_{*}}^{+}(A)$, we are done: $\tau_{*}$ is the optimal value for the (10), with $X_{*}$ as the corresponding optimal solution. $\diamond$

We now present the numerical results for the implementation of Algorithm 4, as the dimension of the starting matrix $A$ is varied.

In the next table we shows how the sparsity affects the computations. The dimension is kept fixed $d=850$ in all experiments, while the density parameter for each row of the starting matrix is randomly chosen from the given interval.

If we compare these results with the numerical results for Schur stabilization from [14], we can observe a remarkable speed-up in the computational time. The first reason for this is because in the updated algorithm we do not need to conduct the greedy method until the end: we quit it as soon as the strongly Schur stable matrix is obtained. By quitting the greedy method before it finishes, we usually avoid obtaining matrices with zero spectral radius. As practical experiments showed, computation of the leading eigenvector for the zero spectral radius matrices can be drastically slow, especially for the very sparse big matrices. The big chunk of computational time in implementation of the algorithm from [14] in fact goes for iterating through zero spectral radius matrices and computing their leading eigenvectors. In the updated procedure this is effectively avoided, providing us with significantly faster computations.

Comparing the experimental results for Algorithm 3 and Algorithm 4, we see that Hurwitz stabilization is more time demanding than the Schur stabilization. In Algorithm 3 it is impossible to avoid matrices with spectral abscissa equal to zero, and the computation of their leading eigenvectors can also get very time consuming.

4.5. The importance of the solution set

The following example illustrates the importance of imposing the set of allowable solutions for the stabilization problems.

Example. Observe the non-negative matrix

[TABLE]

with spectral abscissa (i.e. spectral radius) greater than one. If we want to find a closest non-negative matrix having spectral abscissa (i.e. spectral radius) equal to one, we have a Schur stabilization problem. Solving it, we obtain the matrix

[TABLE]

as the closest Schur stable, with $\tau_{*}=5.764$. However if we expand our set of admissible solution to a set of Metzler matrices, we get

[TABLE]

as the solution. In this case the distance to the starting matrix will be $\tau_{*}=5.4$, which is smaller than for the closest non-negative matrix.

Remark In the previous example we found the closest Metzler matrix to a matrix $A$ having $\eta(A)=1$ by slightly modifying Step 4: we used greedy procedure on the balls of Metzler matrices ${\mathcal{B}}_{\tau}(A)$, instead on the balls of non-negative matrices ${\mathcal{B}}_{\tau}^{+}(A)$. Further, we changed spectral radius to spectral abscissa, and determined the indices $l_{i}$ using the formula (14).

5. Applications

5.1. Checking the reliability of a given data

One of the most prominent uses of Metzler matrices is for mathematical modelling [3, 4, 16]. Numerous phenomena can be modelled by positive linear dynamical systems. In most simple case these systems have the form $\dot{{\boldsymbol{x}}}=A{\boldsymbol{x}}$, where ${\boldsymbol{x}}={\boldsymbol{x}}(t)$ is a vector of time-dependent unknowns, and $A$ is a time-independent Metzler matrix. The coefficients of $A$ are determined from the gathered data. Often, researchers are interested in examining the stability of the model, which is reflected by the Hurwitz stability of matrix $A$. However, gathered data might contain some errors (i.e. due to the measurement imprecisions), which can lead us to the model with qualitative behaviour completely different than the real picture. To check if our data is reliable, we can optimize spectral abscissa on the ball ${\mathcal{B}}_{\varepsilon}(A)$.

For example, let us assume that the matrix of our model-system $A$ is Hurwitz stable and that each entry contains an error not bigger than $\varepsilon$. To check if our model is reliable, we need to maximize the spectral abscissa on ${\mathcal{B}}_{\varepsilon}(A)$.

Let $X_{*}$ be the optimal matrix. If $\eta(X_{*})<0$, we are safe and the qualitative behaviour of our model will agree with the real state of affair, even in the case of making the biggest predicted errors. On the other hand, if we obtain $\eta(X_{*})>0$, we cannot claim that our model reliably describes the phenomenon. In this case one needs to further refine our data gathering methods.

Similar reasoning can be applied if for our model we have $\eta(A)>0$, and we want to check if the real system is unstable as well. However, in this case we would actually need to minimize spectral abscissa on ${\mathcal{B}}_{\varepsilon}$.

5.2. Stabilization of 2D postitive linear switching systems

Positive linear switching systems (LSS) are an important tool of mathematical modelling. They are a point of an extensive research with many applications [27, 28, 29]. Let $A_{i},\ i=1,\ldots,N$ be family of Metzler matrices. A positive LSS is given with:

[TABLE]

where $x\in\mathbb{R}^{d}$ and $\sigma:\mathbb{R}_{\geqslant 0}\rightarrow\{1,\ldots,N\}$ is a piecewise constant switching signal.

A crucial question in the theory of LSS is the asymptotic stability under the arbitrary switching signal. The following theorem provides us with the necessary condition:

The converse of Theorem 10 is true only for the two-dimensional case [28]. Therefore, we can use Hurwitz stabilization on the unstable convex combinations to build an algorithm for the stabilization of the 2D positive LSS. By this we mean constructing a stable LSS from the original system, as shown below.

Assume that $2\times 2$ Metzler matrices $A_{i},\ i=1,\ldots,N$ are Hurwitz stable. If some of them is not, find and replace it with its closest stable. Suppose that there exists a matrix in $\textbf{{\rm co}}\{A_{1},\ldots,A_{N}\}$ that is not Hurwitz stable.

Let $A=\sum_{i=1}^{N}\alpha_{i}A_{i}$ be a convex combination in $\textbf{{\rm co}}\{A_{1},\ldots,A_{N}\}$ with the largest spectral abscissa $(\sum_{i=1}^{N}\alpha_{1}=1,\ \alpha_{i}\in[0,1],\ i=1,\ldots,N)$. Find its closest Hurwitz stable matrix $A^{\prime}$, and denote by $\tau_{*}$ the optimal distance.

We now need to decompose the matrix $A^{\prime}=\sum_{i=1}^{N}\alpha_{i}A^{\prime}_{i}$, while taking care about matrices $A_{i}$. We have

[TABLE]

Choose values $\tau_{i}$ and matrices $A^{\prime}_{i}\in{\mathcal{B}}_{\tau_{i}}(A_{i})$, so that (16) and $A^{\prime}=\sum_{i=1}^{N}\alpha_{i}fA^{\prime}_{i}$ are satisfied.

Let $A_{*}$ be a convex combination in $\textbf{{\rm co}}\{A^{\prime}_{1},\ldots,A^{\prime}_{N}\}$. If $\eta(A_{*})<0$, then we are done: the switching system built from matrices $A^{\prime}_{i}$ is asymptotically stable. Else, we should make a different choice of $\tau_{i}$, or restart everything, but now starting with matrices $A^{\prime}_{i}$.

5.3. Closest stable sign-matrix and its application to positive LSS

The notions of sign-matrices and sign-stability originated from the problems in economy [20, 21], ecology [24, 25], and chemistry [26]. Ever since, those concepts have been a point of interest in mathematical literature [16, 17, 18]. Metzler sign-matrices are very useful for the dynamical system modelling. They come in handy if we do not posses quantitative data, but just the information on the sign of coefficients of the observed system. By analysing the sign-matrix of the system, we can discern its qualitative behaviour. As we will see later on, Metzler sign-matrices can be used as a tool for stability analysis of linear switching systems.

Denote by ${\mathcal{M}}_{\rm sgn}$ set of all real Metzler matrices with entries from $\{-1,0,1\}$. Hurwitz stable matrices from ${\mathcal{M}}_{\rm sgn}$ have one peculiar property: replacing any entry by an arbitrary real number of the same sign does not influence the stability.

Example. The matrix $A\in{\mathcal{M}}_{\rm sgn}$

[TABLE]

is Hurwitz stable with $\eta(A)=-1$. Metzler matrix $A^{\prime}$ with the same sign pattern, given by

[TABLE]

has $\eta(A^{\prime})=-10^{-9}$, i.e., is also Hurwitz stable. It remains stable in spite of very big changes (${\mathcal{O}}(10^{9})$) of all non-negative entries!

This property can be derived by studying sign-matrices, which we introduce now.

In analogy with the real Metzler matrices, we can also consider Hurwitz stability of sign-matrices.

[TABLE]

[TABLE]

$M$* is (strongly) Hurwitz stable if and only if ${\rm sgn}(M)$ is (strongly) Hurwitz stable. Moreover, $M$ is strongly Hurwitz stable if and only if ${\rm sgn}(M)+I$ is the adjacency matrix of an acyclic graph. $\Box$*

From Theorem 11 we see why the Hurwitz stable matrices from ${\mathcal{M}}_{\rm sgn}$ do not change their stability, even as we change their non-zero entries.

The problem of Hurwitz stabilization of real Metzler matrices can be formulated for the Metzler sign-matrices as well. Theorem 11 gives a way to do so: for a given Hurwitz unstable Metzler sign-matrix $M$, we need to find closest Hurwitz stable Metzler sign-matrix $X$. In other words, $X$ should solve

[TABLE]

Remark. Notice that in Hurwitz stabilization of sign-matrix we allow our optimal solution to have negative spectral abscissa as well. We do this since in some cases a solution with zero spectral abscissa does not exist555 see Examples below..

Denote ${\mathcal{B}}_{k}^{\rm sgn}={\mathcal{B}}_{k}\cap{\mathcal{M}}_{\rm sgn}$. Using minimization of spectral abscissa on ${\mathcal{B}}_{k}^{\rm sgn}$ we can present a simple procedure for solving problem (17).

Algorithm 5: Finding the closest Hurwitz stable sign-matrix

Let $M$ be a Hurwitz unstable sign-matrix. Set $k=\lfloor{\frac{\|{\rm sgn}(A)\|}{2}}\rfloor$. By $X_{k}$ denote the sign-matrix such that $\eta_{k}=\eta({\rm sgn}(X_{k}))$ is minimal on ${\mathcal{B}}_{k}^{\rm sgn}$.

Doing a bisection in $k$, minimize the spectral abscissa on the ball ${\mathcal{B}}_{k}^{\rm sgn}$. Do this until sign-matrices $X_{k-1}$ and $X_{k}$ are obtained, with the (minimal) spectral abscissas $\eta_{k-1}>0$ and $\eta_{k}\leqslant 0$.

Take $k=k_{*}$ as the optimal value and $X_{k}=X_{*}$ as the optimal solution of (17). $\diamond$

Example. For the unstable Metzler sign-matrix

[TABLE]

we find

[TABLE]

with $\eta({\rm sgn}(M_{*}))=0$ as the closest sign-stable, and optimal distance $k_{*}=2$.

Example. For the unstable Metzler sign-matrix

[TABLE]

we find

[TABLE]

with $\eta({\rm sgn}(M_{*}))=-1$ as the closest sign-stable, and optimal distance $k_{*}=2$. The optimal solution with zero spectral abscissa is impossible to find, because at the distance $k=1$ we obtain the matrix

[TABLE]

with $\eta({\rm sgn}(M^{\prime}))=0.46$ as the one with the minimal spectral abscissa.

We have seen that Theorem 10 works both ways only in the two-dimensional case. However, there is a criterion of asymptotic stability, valid in any dimension, which involves sign-matrices and sign-stability.

Now, assume that the sign-matrices $M_{i},\ i=1,\ldots,N$ have negative diagonal entries, that they are Hurwitz stable, but their sum $M=\sum_{i=1}^{N}M_{i}$ is not. Using Algorithm 5, we can find closest Hurwitz sign-stable matrix $M^{\prime}$ to $M$. Let $k_{*}$ be the optimal value. Now we need to decompose $M^{\prime}$ into the sum $\sum_{i=1}^{N}M^{\prime}_{i}$, while taking into the consideration the structure of matrices $M_{i}$. We have

[TABLE]

We need to choose values $k_{i}$ and matrices $M^{\prime}_{i}\in{\mathcal{B}}_{k_{i}}^{\rm sgn}(M_{i})$, so that (18) and $M^{\prime}=\sum_{i=1}^{N}M^{\prime}_{i}$ are satisfied. Matrices $M^{\prime}_{i}$, including their sum, will have negative entries on the main diagonal and be Hurwitz stable. Therefore, by Theorem 12, any LSS we construct from matrices in ${\mathcal{Q}}(M^{\prime}_{i})$ will be asymptotically stable, while keeping the similar structure to the initial system.

Example. Observe a positive linear switching system, switching between the following three strongly Hurwitz stable matrices:

[TABLE]

Their corresponding sing-matrices are given by

[TABLE]

and their sum $M=M_{1}+M_{2}+M_{3}$ with

[TABLE]

Since sign-matrix $M$ is Hurwitz unstable (with $\eta({\rm sgn}(M))=1.303$), our positive LSS will not be asymptotically stable. Applying Algorithm 6 on it, we find its closest strongly Hurwitz stable sing-matrix

[TABLE]

with an optimal distance $k_{*}=1$. We select $k_{1}=k_{2}=k_{3}=1$ and decompose $M^{\prime}=M^{\prime}_{1}+M^{\prime}_{2}+M^{\prime}_{3}$ by choosing $M^{\prime}_{i}\in{\mathcal{B}}_{k_{i}}^{\rm sgn}(M_{i}),\ i=1,2,3$. We have:

[TABLE]

Following the sign pattern of matrices $M^{\prime}_{i}$ and correspondingly cutting the interdependencies in matrices $A_{i}$, we get

[TABLE]

A newly obtained LSS built from matrices $A^{\prime}_{i}$ will be asymptotically stable under the arbitrary switching.

Acknowledgement. Countless thanks to prof V. Yu. Protasov for the hours of inspiring discussions, all the valuable remarks and the support he was selflessly giving during the writing of this paper.