Computing the Kreiss Constant of a Matrix

TL;DR

This paper introduces the first globally convergent algorithms for accurately computing the Kreiss constant of a matrix, with improved complexity and scalable optimization methods for large matrices.

Contribution

It presents novel algorithms for computing the Kreiss constant with proven convergence and reduced computational complexity, including continuous and discrete-time variants.

Findings

Algorithms achieve arbitrary accuracy in computing Kreiss constants.

Complexity reduced from $ ext{O}(n^6)$ to $ ext{O}(n^4)$ on average.

Efficient local optimization methods for large-scale matrices.

Abstract

We establish the first globally convergent algorithms for computing the Kreiss constant of a matrix to arbitrary accuracy. We propose three different iterations for continuous-time Kreiss constants and analogues for discrete-time Kreiss constants. With standard eigensolvers, the methods do $\mathcal{O}(n^6)$ work, but we show how this theoretical work complexity can be lowered to $\mathcal{O}(n^4)$ on average and $\mathcal{O}(n^5)$ in the worst case via divide-and-conquer variants. Finally, locally optimal Kreiss constant approximations can be efficiently obtained for large-scale matrices via optimization.

Tables

Table 1. Table 1: The upper half of the table shows data for our continuous-time example, while the lower half is for the discrete-time example. A dash indicates the column is not relevant for the given method. Chebfun is simply taking the max of the chebfuns produced to make Fig.   1 , where Chebfun was supplied a reasonably small interval already known to contain the maximizer attaining 𝒦 ​ ( A ) 𝒦 𝐴 \mathcal{K}(A) .

	$𝒦\left(A\right)$	$z_0$	level-set tests	restarts	seconds
Chebfun	$1.29186707005845\times {10}^5$	—	—	—	97.28
Algorithm 1	$1.29186707011257\times {10}^5$	$6+6𝐢$	14	2	1.95
Algorithm 1	$1.29186707015035\times {10}^5$	$6+6𝐢$	2	2	0.76
Algorithm 2	$1.29181067678395\times {10}^5$	$1$	89	—	21.48
Chebfun	1.89501339090609	—	—	—	74.04
Algorithm 1	1.89501339090580	$-1+1𝐢$	15	3	7.18
Algorithm 1	1.89501339090580	$-1+1𝐢$	3	3	1.62
Algorithm 2	1.89501305930067	$2$	81	—	41.04

Equations

$$x_{k+1}=A{x}_k$$

$$\mathcal{K}(A)\le \underset{k\ge 0}{\sup }\Vert A^k\Vert \le \mathrm{e}n\mathcal{K}(A),$$

$$\mathcal{K}(A)=\underset{z\in \mathbb{C},|z|>1}{\sup }(|z|-1)\Vert (zI-A{)}^{-1}\Vert .$$

$$\mathcal{K}(A)=\underset{\epsilon >0}{\sup }\frac{{\rho }_{\epsilon }(A)-1}{\epsilon },$$

$${\rho }_{\epsilon }(A)$$

$$=\max \{|z|:z\in \mathbb{C},\Vert (zI-A{)}^{-1}\Vert \ge {\epsilon }^{-1}\}.$$

$$\dot{x}=Ax,$$

$$\mathcal{K}(A)\le \underset{t\ge 0}{\sup }\Vert {\mathrm{e}}^{tA}\Vert \le \mathrm{e}n\mathcal{K}(A)$$

$$\mathcal{K}(A)=\underset{z\in \mathbb{C},\mathrm{Re}z>0}{\sup }(\mathrm{Re}z)\Vert (zI-A{)}^{-1}\Vert$$

$$\mathcal{K}(A)=\underset{\epsilon >0}{\sup }\frac{{\alpha }_{\epsilon }(A)}{\epsilon },$$

$${\alpha }_{\epsilon }(A)$$

$$=\max \{\mathrm{Re}z:z\in \mathbb{C},\Vert (zI-A{)}^{-1}\Vert \ge {\epsilon }^{-1}\}.$$

$$\tfrac{\partial^{2}}{\partial x\partial y}\lambda_{j}\bigg{|}_{x=\hat{x},y=\hat{y}}=q_{j}^{*}\tfrac{\partial^{2}A(\hat{x},\hat{y})}{\partial x\partial y}q_{j}+2\sum_{k\neq j}\frac{q_{j}^{*}\tfrac{\partial A(\hat{x},\hat{y})}{\partial x}q_{k}\cdot q_{j}^{*}\tfrac{\partial A(\hat{x},\hat{y})}{\partial y}q_{k}}{\lambda_{j}-\lambda_{k}}.$$

$$\dot{x}=Ax+Bu,$$

$$\tau (A,B)=\underset{z\in \mathbb{C}}{\min }{\sigma }_{\min }\left(\left[\begin{array}{cc}A-zI& B\end{array}\right]\right)=\underset{x,y\in \mathbb{R}}{\min }f(x,y),$$

$$f(x,y)=f(x+\eta ,y)=\gamma .$$

$$\tau (A,B)>\gamma -\frac{\eta }{2}.$$

$$\mathcal{K}(A{)}^{-1}=\underset{x>0,y\in \mathbb{R}}{\inf }{\sigma }_{\min }\left(\frac{(x+\mathbf{i}y)I-A}{x}\right)=\underset{x>0,y\in \mathbb{R}}{\inf }g(x,y),$$

$$g(x,y)={\sigma }_{\min }(G(x,y))\text{and}G(x,y)=\frac{(x+\mathbf{i}y)I-A}{x}.$$

$$\frac{\partial G(x,y)}{\partial x}=\frac{xI-((x+\mathbf{i}y)I-A)}{x^2}=\frac{A-\mathbf{i}yI}{x^2}\text{and}\frac{\partial G(x,y)}{\partial y}=\frac{\mathbf{i}I}{x}.$$

$$\nabla g(\hat{x},\hat{y})=\mathrm{Re}\left[\begin{array}{c}{u}^{\ast }\frac{\partial G(x,y)}{\partial x}v\\ u^{\ast }\frac{\partial G(x,y)}{\partial y}v\end{array}\right].$$

$$\left[\begin{array}{cc}0& G(x,y)\\ G(x,y{)}^{\ast }& 0\end{array}\right],$$

$$\frac{{\partial }^{2}G(x,y)}{\partial x^2}=\frac{-2(A-\mathbf{i}yI)}{x^3},\frac{{\partial }^{2}G(x,y)}{\partial y^2}=0,\text{and}\frac{{\partial }^{2}G(x,y)}{\partial x\partial y}=\frac{-\mathbf{i}I}{x^2}.$$

$$\mathcal{K}(A{)}^{-1}=\underset{r>1,\theta \in [0,2\pi )}{\inf }{\sigma }_{\min }\left(\frac{r{\mathrm{e}}^{\mathbf{i}\theta }I-A}{r-1}\right)=\underset{r>1,\theta \in [0,2\pi )}{\inf }h(r,\theta ),$$

$$h(r,\theta )={\sigma }_{\min }(H(r,\theta ))\text{and}H(r,\theta )=\frac{r{\mathrm{e}}^{\mathbf{i}\theta }I-A}{r-1}.$$

$$\frac{\partial H(r,\theta )}{\partial r}=\frac{(r-1){\mathrm{e}}^{\mathbf{i}\theta }I-(r{\mathrm{e}}^{\mathbf{i}\theta }I-A)}{(r-1{)}^2}=\frac{A-{\mathrm{e}}^{\mathbf{i}\theta }I}{(r-1{)}^2}\text{and}\frac{\partial H(r,\theta )}{\partial \theta }=\frac{\mathbf{i}r{\mathrm{e}}^{\mathbf{i}\theta }I}{r-1},$$

$$\frac{{\partial }^{2}H(r,\theta )}{\partial r^2}=\frac{-2(A-{\mathrm{e}}^{\mathbf{i}\theta }I)}{(r-1{)}^3},\frac{{\partial }^{2}H(r,\theta )}{\partial {\theta }^2}=\frac{-r{\mathrm{e}}^{\mathbf{i}\theta }I}{r-1},\text{and}\frac{{\partial }^{2}H(r,\theta )}{\partial r\partial \theta }=\frac{-\mathbf{i}{\mathrm{e}}^{\mathbf{i}\theta }I}{(r-1{)}^2}.$$

$$g(x,y)=g(x,y+\eta )=\gamma .$$

$$\mathcal{K}(A{)}^{-1}>\gamma -\frac{\eta }{2x_{\star }}.$$

$$b_1=(x_{\star },y_{\star }-{\eta }_1)\text{and}{b}_2=(x_{\star },y_{\star }+{\eta }_2)$$

$$g(x_{\star },y_{\star }-{\eta }_1)=g(x_{\star },y_{\star }+{\eta }_2)=\gamma ,$$

Full text

Keywords: discontinuity of Kreiss constants, inverses of Kronecker sums, distance to uncontrollability algorithms, transient growth, pseudospectra

Notation: $\|\cdot\|$ denotes the spectral norm, $\sigma_{\min}(\cdot)$ the smallest singular value, $\Lambda(\cdot)$ the spectrum, $J=\begin{bmatrix}0&I\\ -I&0\end{bmatrix}$, a matrix $A\in\mathbb{C}^{2n\times 2n}$ is Hamiltonian if $(JA)^{*}=JA$ and symplectic if $A^{*}JA=J$, a matrix pencil $A-\lambda B$ is symplectic if $A^{*}JA=B^{*}JB$, $\mathrm{e}$ is Euler’s number $2.71828\ldots$, and $I_{n}$ is the $n\times n$ identity, though we will often omit the subscript when the dimension is clear.

1 Introduction

Given a matrix $A\in\mathbb{C}^{n\times n}$, the ordinary difference equation

[TABLE]

is asymptotically stable if $A$ is Schur stable, i.e., if $\rho(A)<1$, where $\rho$ denotes the spectral radius. While $\rho(A)$ tells one about the asymptotic behavior of (1.1), it does not convey information about its transient behavior. For that, we can look at the Kreiss Matrix Theorem, which says for any matrix $A\in\mathbb{C}^{n\times n}$ [TE05, Eq. 18.2]

[TABLE]

where the Kreiss constant $\mathcal{K}(A)$ is given by [TE05, p. 143]

[TABLE]

As also noted in [TE05, p. 143], $\mathcal{K}(A)$ has an equivalent formulation of

[TABLE]

where the $\varepsilon$-pseudospectral radius $\rho_{\varepsilon}$ is defined by

[TABLE]

For any matrix $A\in\mathbb{C}^{n\times n}$, $\mathcal{K}(A)\geq 1$, and $\mathcal{K}(A)$ may be arbitrarily large. As is well known, a matrix $A$ is power-bounded, i.e., $\mathcal{K}(A)<\infty$, if and only if $\rho(A)\leq 1$ and all eigenvalues of $A$ with modulus 1 are nondefective. If $A$ is normal and $\rho(A)\leq 1$, then $\mathcal{K}(A)=1$.

As discussed on [TE05, p. 177], the original statement by Kreiss in 1962 [Kre62] actually had a far looser upper bound than (1.2): approximately $c^{n^{n}}\mathcal{K}(A)$. The reduction of the constant factor to its current form in fact occurred over nearly thirty years in at least nine separate steps, with Spijker proving the conjecture of [LT84, p. 590] to finally obtain the (in a certain sense) tight factor of $\mathrm{e}n$ in 1991 [Spi91].

The Kreiss Matrix Theorem also comes in a continuous-time variant for an ordinary differential equation

[TABLE]

which is asymptotically stable if $A$ is Hurwitz stable, i.e., if $\alpha(A)<0$, where $\alpha$ denotes the spectral abscissa. In this case, the Kreiss Matrix Theorem states [TE05, Eq. 18.8]

[TABLE]

where by [TE05, Eq. 14.7], $\mathcal{K}(A)$ is now equivalently given by either

[TABLE]

or

[TABLE]

and where the $\varepsilon$-pseudospectral abscissa $\alpha_{\varepsilon}$ is defined by

[TABLE]

Like the discrete-time case, $\mathcal{K}(A)\geq 1$ and can be arbitrary large. If $A$ is normal and $\alpha(A)\leq 0$, then $\mathcal{K}(A)=1$.

Despite the wealth of work done over decades towards making the upper bound of the Kreiss Matrix Theorem now tight, there has been no algorithm given to actually compute $\mathcal{K}(A)$ with guarantees. In the literature, $\mathcal{K}(A)$ is often just approximated by plotting (1.4) or (1.9); e.g., see [EK17] and [Men06, Chapter 3.4.1].

In this paper, we propose the first globally convergent algorithms to compute both continuous- and discrete-time Kreiss constants to arbitrary accuracy. We assume that $A$ is nonnormal, as otherwise computing its Kreiss constant just involves checking if $A$ is unstable. Furthermore, we assume that $\alpha(A)<0$ or $\rho(A)<1$ holds, respectively, in the continuous- or discrete-time case. With standard eigensolvers, the three different methods we propose all have $\mathcal{O}(n^{6})$ work complexities,111 We use the standard convention of treating dense computations of singular values, eigenvalues, solutions of Sylvester equations, etc., as atomic operations with cubic costs in the dimensions of the associated matrices. We additionally assume these costs become linear in the dimension of the matrices when corresponding sparse methods are available.

but by also developing so-called divide-and-conquer variants, we show how this theoretical work complexity reduces to $\mathcal{O}(n^{4})$ on average and $\mathcal{O}(n^{5})$ in the worst case. Our work also shows that locally optimal approximations to $\mathcal{K}(A)$ can be efficiently and reliably obtained for large-scale matrices via standard optimization techniques. Furthermore, we establish some variational properties of Kreiss constants, including that the Kreiss constant of a matrix $A$ is not a continuous function with respect to the entries of $A$. Finally, as a side effect of our work, we also propose an important modification to the distance-to-uncontrollability algorithms of [Gu00, BLO04, GMO*+*06] to greatly improve their reliability in practice.

The contributions of the paper are structured as follows. In §2, we introduce the so-called distance to uncontrollability and present a theorem of Gu [Gu00, Theorem 3.1] and distance-to-uncontrollability algorithms that are based upon it. In §3, we establish variational properties of Kreiss constants that we will need and show that there is a potentially exploitable similarity to computing the distance to uncontrollability. Then, in §4, we develop theorems for continuous-time $\mathcal{K}(A)$ that are analogues of the aforementioned theorem of Gu, but nevertheless show that, due to key structural differences, existing distance-to-uncontrollability algorithms will not directly extend to Kreiss constants. By developing a so-called globality certificate in §5, we present our first algorithm for computing continuous-time Kreiss constants, which is an optimization-with-restarts method using backtracking, and then establish an asymptotically faster divide-and-conquer variant that is inspired by [GMO*+*06]. We then modify the premise of our globality certificate in a crucial way in §6 to develop a second certificate with significant structural differences and properties; this alternative certificate enables two other algorithms for continuous-time $\mathcal{K}(A)$, an optimization-with-restarts method without backtracking and a trisection iteration, which can be considered closer analogues of the distance-to-uncontrollability algorithms of [BLO04]. Faster divide-and-conquer versions of these two methods are also developed. In §7, we consider the case of discrete-time $\mathcal{K}(A)$ and present discrete-time analogues of all of our continuous-time algorithms and associated theoretical results; to the best of our knowledge, this is also the first extension of the 2D level-set ideas of [Gu00] to a discrete-time setting, and it turns out to have surprising differences. Numerical examples are presented in §8, with concluding remarks made in §9.

For this paper, the following general definition and two theorems will be needed. The theorems can be found in several places in various forms, such as [Lan64, OW95] and [Lau05, Theorem 13.16], respectively.

Definition 1.1**.**

Given a domain $\Omega\subseteq\mathbb{R}$, a function $f:\Omega\to\mathbb{R}$ has a global Lipschitz constant (GLC) of $c\geq 0$ if $|f(x)-f(y)|\leq c|x-y|$ for all $x,y\in\Omega$.

Theorem 1.2**.**

For $x,y\in\mathbb{R}$, let $A(x,y)$ be a twice-differentiable $n\times n$ Hermitian matrix family, and for a point $(\hat{x},\hat{y})$, let $\lambda_{1}\geq\ldots\geq\lambda_{n}$ be the eigenvalues of $A(\hat{x},\hat{y})$ with associated unit-norm eigenvectors $q_{1},\ldots,q_{n}$. Then assuming $\lambda_{j}$ is unique,

[TABLE]

Theorem 1.3**.**

Let $A\in\mathbb{C}^{n\times n}$ with $\lambda_{j}\in\Lambda(A)$ for $j=1,\ldots,n$ and $B\in\mathbb{C}^{m\times m}$ with $\mu_{k}\in\Lambda(B)$ for $k=1,\ldots,m$. Then the Kronecker sum $A\oplus B=I_{m}\otimes A+B\otimes I_{n}$ has eigenvalues $\lambda_{j}+\mu_{k}$, for all pairs of $j$ and $k$.

2 Computing the distance to uncontrollability

Given $A\in\mathbb{C}^{n\times n}$ and $B\in\mathbb{C}^{n\times m}$, consider the linear control system

[TABLE]

where the state $x\in\mathbb{C}^{n}$ and $u\in\mathbb{C}^{m}$, the control input, are both dependent on time. The system (2.1) is controllable if given respective initial and final states $x(0)$ and $x(T)$ with $T>0$, there exists a control $u(\cdot)$ that realizes some trajectory $x(\cdot)$ with endpoints $x(0)$ and $x(T)$. The distance to uncontrollability, which we denote as $\tau(A,B)$, can be computed via solving the nonconvex optimization problem [Eis84]

[TABLE]

where $f(x,y)=\sigma_{\min}(F(x,y))$ and $F(x,y)=\begin{bmatrix}A-(x+\mathbf{i}y)I&B\end{bmatrix}$. The first practical algorithm to address computing $\tau(A,B)$ is due to Gu [Gu00], based on the following result [Gu00, Theorem 3.1].

Theorem 2.1**.**

Let $\gamma,\eta\geq 0$ be given. If $\tau(A,B)\leq\gamma$ and $\eta\in[0,2(\gamma-\tau(A,B))]$, then there exists a pair $x,y\in\mathbb{R}$ such that

[TABLE]

Corollary 2.2**.**

Let $\gamma,\eta\geq 0$ be given. If there do not exist any pairs $x,y\in\mathbb{R}$ such that (2.3) holds, then

[TABLE]

The proof of Theorem 2.1 relies on the fact that $f(x,y)$ has a GLC of 1 with respect to either $x$ or $y$.

What [Gu00, Section 3.2] additionally devised was a sequence of computations to verify whether either (2.3) or (2.4) holds for a given choice of $\gamma$ and $\eta$. This verification procedure, when using exact arithmetic, is able to detect and find any points $(x,y)$ such that (2.3) is satisfied. If so, the test returns these points and $\tau(A,B)\leq\gamma$ is verified. Otherwise, the test asserts no pairs satisfy (2.3) and so (2.4) instead must hold. Using this procedure, Gu proposed a bisection-like scheme to estimate $\tau(A,B)$ to within a factor of two. For initialization, $\gamma\coloneqq f(0,0)$ and $\eta\coloneqq\gamma$. If the test verifies $\tau(A,B)\leq\gamma$, then $\gamma$ and $\eta$ are both halved (so $\eta=\gamma$ still holds) and the test is done with these smaller values. Otherwise, (2.4) holds and so $\gamma$ and $\tau(A,B)$ are within a factor of two and Gu’s method terminates.

As noted in [BLO04, p. 358], it is tempting to try to obtain $\tau(A,B)$ to higher precision via a true bisection method, i.e., one that would update both upper and lower bounds, unlike Gu’s method which only updates an upper bound. The problem with this approach is that in order to ascertain whether the current estimate $\gamma$ is essentially a lower bound to $\tau(A,B)$ via (2.4), one would have to perform the verification procedure for $\eta\approx 0$. Unfortunately, this is not tenable in the presence of rounding errors, as Gu’s procedure becomes more and more numerically unreliable as $\eta\to 0$, i.e., points satisfying (2.3) may not be detected. Consequently, in practice, the lower bound will generally be erroneously updated at some point, thus preventing convergence to $\tau(A,B)$.222In 6.3, we discuss the unreliability of Gu’s procedure in more detail and explain how our new tests avoid key numerical pitfalls. Besides being useful for our $\mathcal{K}(A)$ algorithms, our modifications can also improve the reliability of the $\tau(A,B)$ methods of [Gu00, BLO04, GMO*+*06].

Using Gu’s verification procedure (as Gu specified it, i.e., without modifications), Burke, Lewis, and Overton instead proposed a trisection algorithm [BLO04, Algorithm 5.2] that balances how much the lower bound is updated with how quickly the value of $\eta$ is decreased, precisely to postpone the numerical unreliability of Gu’s procedure as long as possible. Trisection works as follows. Let $L\coloneqq 0$ and $U\coloneqq f(0,0)$ be initial lower and upper bounds, respectively. Then on the $k$th iteration, $\eta_{k}\coloneqq\tfrac{2}{3}(U-L)$ and $\gamma_{k}\coloneqq L+\eta_{k}$ are set as the current values of $\eta$ and $\gamma$ for the verification test. If the test finds points satisfying (2.3), then the upper bound is updated $U\coloneqq\gamma_{k}$. Otherwise, (2.4) holds so we know that $\tau(A,B)\geq\gamma_{k}-\tfrac{\eta_{k}}{2}=L+\tfrac{\eta_{k}}{2}$, and so now the lower bound can be updated $L\coloneqq L+\tfrac{\eta_{k}}{2}$. Thus, the new interval has length $\tfrac{2}{3}(U-L)$. As the trisection algorithm is linearly convergent, Gu’s verification procedure will have to be invoked many times, which already is $\mathcal{O}(n^{6})$ work by itself using standard dense eigensolvers. Thus, the trisection algorithm also has a large constant factor term hidden away in its asymptotic work complexity. Also, trisection is not a panacea for the numerical issues of Gu’s procedure; although $\eta_{k}\to 0$ only in the limit, if $\tau(A,B)$ is small, then $\eta_{k}$ must become commensurately small in order for trisection to attain any digits of accuracy; see Lemma B.1 and Corollary B.2 in Appendix B.

In the same paper, Burke, Lewis, and Overton also proposed a second algorithm for $\tau(A,B)$ [BLO04, Algorithm 5.3] and advocated it as preferable to trisection. This optimization-with-restarts method also relies on Gu’s verification procedure, now as a globality certificate, and additionally, on the fact that $f(x,y)$ is semialgebraic and so $f(x,y)$ only has a finite number of locally minimal values; see [BLO04, p. 359]. This second method thus works by using optimization techniques to find a minimizer of (2.2) with function value $f_{k}$ and then uses Gu’s verification procedure with carefully chosen values of $\gamma$ and $\eta$ so that the test checks if $f_{k}$ is sufficiently close to $\tau(A,B)$; for some relative tolerance $\texttt{tol}>0$, the specific values are $\gamma\coloneqq f_{k}\cdot(1-0.5\cdot\texttt{tol})$ and $\eta\coloneqq f_{k}\cdot\texttt{tol}$.333 Note that [BLO04] writes these in an equivalent but different form using $\delta_{1}=\gamma$ and $\delta_{2}=\gamma-\tfrac{\eta}{2}$. Otherwise, if $f_{k}\not\approx\tau(A,B)$ to tolerance, the certificate provides one or more new starting points from which optimization can be restarted with the guarantee that a better (lower) minimum of (2.2) will be found, hence optimization is restarted in a loop until the certificate indeed asserts that the desired accuracy has been attained. By construction, $f_{k}$ is monotonically decreasing and optimization-with-restarts must terminate with $f_{k}\approx\tau(A,B)$ to tolerance in a finite number of restarts. Although it is not clear exactly how many restarts will occur, only a handful are typically needed in practice, if any. Furthermore, the optimization phases are relatively cheap, requiring $\mathcal{O}(n^{3})$ work with a relatively low constant factor, since minimizers of $f(x,y)$ can generally be found with superlinear or even quadratic convergence. As a result, optimization-with-restarts is almost always many times faster than trisection. However, optimization-with-restarts can still be susceptible to the numerical difficulties of Gu’s verification procedure, since $\eta$ itself may still become very small, e.g., if either high accuracy is desired or $\tau(A,B)$ is small.

Finally, to address the high cost of Gu’s verification procedure, though not necessarily its numerical issues, [GMO*+*06] proposed a divide-and-conquer strategy that lowers the asymptotic work complexity of Gu’s procedure to $\mathcal{O}(n^{4})$ on average and $\mathcal{O}(n^{5})$ in the worst case. This benefits all of the aforementioned algorithms.

3 Variational properties and the inverse of the Kreiss constant

We now establish some variational properties of Kreiss constants, which in turn show that locally optimal approximations to $\mathcal{K}(A)$ can be efficiently computed via optimization, even if $A$ is large. We also show how the problem of computing $\mathcal{K}(A)$ shares some similarity with computing $\tau(A,B)$. We begin with the following result.

Lemma 3.1**.**

The Kreiss constant $\mathcal{K}$ is not always continuous at $A$, as it may instantaneously jump to/from $\infty$. However, $\mathcal{K}$ is continuous at $A$ if $\alpha(A)<0$ holds (continuous-time case) or $\rho(A)<1$ holds (discrete-time case).

Proof.

We begin with the second part. If $\alpha(A)<0$, then $zI-A$ is invertible for all $z$ such that $\mathrm{Re}\,z>0$. By continuity of singular values, $\|(zI-A)^{-1}\|$ in (1.8) is continuous at $A$, and thus so is $\mathcal{K}(A)$. Via an analogous argument with (1.3), the continuity claim also holds for the discrete-time case.

We prove the first part by example. Let $A(\delta)\coloneqq\begin{bmatrix}-0.5&0\\ 0&\delta\end{bmatrix}$ for real scalar $\delta\geq 0$. As $A(\delta)$ is always normal and $\alpha(A(0))=0$, $\mathcal{K}(A(0))=1$ holds. However, as $\alpha(A(\delta))>0$ for any $\delta>0$, $\mathcal{K}(A(\delta))=\infty$ for any $\delta>0$. Using the same example with $\delta\geq 1$ for the discrete-time case, we have that $\mathcal{K}(A(1))=1$ and $\mathcal{K}(A(\delta))=\infty$ for all $\delta>1$. ∎

3.1 The continuous-time case

Identifying $\mathbb{C}$ with $\mathbb{R}^{2}$, consider the inverse of the continuous-time Kreiss constant (1.8), i.e.,

[TABLE]

where

[TABLE]

Like $f(x,y)$, $g(x,y)$ is semialgebraic, and so in the open right half-plane, $g(x,y)$ must have only a finite number of locally minimal function values. We now derive the gradient and Hessian of $g(x,y)$ which will be useful for finding minimizers via quasi-Newton or Newton methods. Although singular values can vary nonsmoothly with respect to matrix entries, they are nevertheless locally Lipschitz, and so this nonsmoothness is confined to a set of measure zero. We first need the first partial derivatives of $G(x,y)$ for $x\neq 0$:

[TABLE]

Let $(\hat{x},\hat{y})$ be such that $g(\hat{x},\hat{y})\neq 0$ is a simple singular value of $G(x,y)$ with associated left and right singular vectors $u$ and $v$. Then, by standard perturbation theory for singular values, it follows that

[TABLE]

Now since $g(x,y)=\sigma_{\min}(G(x,y))$ is also the $n$th eigenvalue (in descending order) of the $2n\times 2n$ Hermitian matrix

[TABLE]

$\nabla^{2}g(\hat{x},\hat{y})$ can be computed by applying Theorem 1.2 to (3.5). Computationally, the necessary first and second partial derivatives of (3.5) can be obtained via the first partials given in (3.3) and the following second partial derivatives:

[TABLE]

Although the full eigendecomposition of (3.5) is needed, it can actually be constructed more or less for free given the full SVD of $G(\hat{x},\hat{y})$; see [BM19, Section 2.2] for details. If $(\hat{x},\hat{y})$ is additionally a distinct (up to conjugacy) global minimizer of $g(x,y)$, then the gradient and Hessian of $\mathcal{K}(A)^{-1}$ with respect to $A$ is equivalent to the gradient and Hessian of $g(\hat{x},\hat{y})$ with respect to $A$. The gradient and Hessian of $\mathcal{K}(A)$ is then simply obtained by applying the chain rule for the inverse.

The cost of obtaining $g(\hat{x},\hat{y})$, its gradient, and its Hessian is $\mathcal{O}(n^{3})$, as they can all be computed given the full SVD of $G(\hat{x},\hat{y})$. Although (3.1) is technically a constrained optimization problem, $g(x,y)\to\infty$ as $x$ approaches zero from the right, assuming that $\mathbf{i}y$ is not an eigenvalue of $A$. Thus, just returning $\infty$ as the value of $g(x,y)$ whenever $x\leq 0$ suffices for using unconstrained optimization solvers to find feasible local/global minimizers of (3.1). Provided $g(x,y)$ is sufficiently smooth about its stationary points, one can expect local quadratic convergence when using a Newton-based optimization method and superlinear convergence with a quasi-Newton method (forgoing the use of the Hessian). Note that scalable methods for computing smallest singular values, e.g., PROPACK [Lar], can also be used to compute $g(\hat{x},\hat{y})$ and its associated pair of left and right singular vectors in order to obtain $\nabla g(\hat{x},\hat{y})$. Thus, by combining such a sparse solver with a quasi-Newton method, one can efficiently obtain locally optimal approximations to Kreiss constants of large-scale matrices.

Remark 3.2**.**

One could also consider using optimization to find maximizers of (1.9), which has the benefit of working with only one optimization variable instead of two. However, computing $\alpha_{\varepsilon}(A)$ is substantially more expensive than the minimum singular value of a matrix; the quadratically-convergent criss-cross algorithm of [BLO03] to compute $\alpha_{\varepsilon}(A)$, as well as the faster method of [BM19], require computing all eigenvalues of $2n\times 2n$ matrices, often several times. Moreover, we have just shown how the explicit Hessian of $g(x,y)$ can easily be computed in order to obtain faster convergence of the iterates produced by optimization methods. Finally, for large-scale $A$ matrices, sparse methods for $\sigma_{\min}(A)$ are generally much faster and more reliable than those for approximating $\alpha_{\varepsilon}(A)$ [GO11, KV14].

That $g(x,y)$ is not so dissimilar to $f(x,y)$ for $\tau(A,B)$ indicates that it might be possible to adapt Gu’s verification procedure to develop globality certificates for $g(x,y)$. Combined with the optimization techniques discussed here, this would enable a globally convergent optimization-with-restarts method for Kreiss constants that terminates within a finite number of restarts.

3.2 The discrete-time case

Again identifying $\mathbb{C}$ with $\mathbb{R}^{2}$, but now using polar coordinates, consider the inverse of the discrete-time Kreiss constant (1.3), i.e.,

[TABLE]

where

[TABLE]

Naturally $h(r,\theta)$ has the same key properties as $g(x,y)$, i.e., it too is semialgebraic and locally Lipschitz. Thus, $h(r,\theta)$ has a finite number of locally minimal function values, and we can consider using optimization to find minimizers of $h(r,\theta)$. We will need the analogous gradients and Hessian of $h(r,\theta)$; for brevity, we just provide the first and second partial derivatives of $H(r,\theta)$ for $r\neq 1$ here, which are respectively

[TABLE]

and

[TABLE]

The costs to compute $h(r,\theta)$ along with its gradient and Hessian also remain as are described in §3.1 for $g(x,y)$, and the variational results above similarly allow the gradient and Hessian of $\mathcal{K}(A)^{-1}$ or $\mathcal{K}(A)$ to be computed in the discrete-time case. To ensure that optimization returns a feasible minimizer to (3.7), it suffices to return $\infty$ for the value of $h(r,\theta)$ whenever $r\leq 1$; this is because for $\mathrm{e}^{\mathbf{i}\theta}$ not an eigenvalue of $A$, $\lim_{r\to 1^{+}}h(r,\theta)\to\infty$. Thus, locally optimal approximations to discrete-time Kreiss constants can be computed efficiently, for small- or large-scale matrices. To develop a globally convergent algorithm, we will need to develop a discrete-time globality certificate.

4 Continuous-time Kreiss constant analogues of Gu’s theorem and their consequences

Before we present our globally convergent iterations for computing Kreiss constants, we first develop analogues of Theorems 2.1 and 2.2. For the time being, we consider continuous-time $\mathcal{K}(A)$ and begin by considering vertically oriented pairs of points on the $\gamma$-level set of $g(x,y)$. For $\tau(A,B)$, Gu considered pairs of level-set points oriented horizontally, but this choice was rather arbitrary. However, as we will soon see, for Kreiss constants the choice of orientation does have important consequences, both theoretically and for our new algorithms.

Theorem 4.1**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If $\mathcal{K}(A)^{-1}\leq\gamma$ and $\eta\in[0,2x_{\star}(\gamma-\mathcal{K}(A)^{-1})]$, then there exists a pair $x,y\in\mathbb{R}$ with $x>0$ such that

[TABLE]

Corollary 4.2**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If there do not exist any pairs $x,y\in\mathbb{R}$ with $x>0$ such that (4.1) holds, then

[TABLE]

To prove Theorem 4.1, we will use the following topology definition and result.

Definition 4.3**.**

Let $\mathcal{A}\subset\mathbb{C}$ be a bounded open (path) connected set and let $\mathcal{A}^{\mathrm{C}}$ be its complement. Furthermore, let $\mathcal{B}=\{z\in\mathbb{C}:z\text{ is in a bounded component of }\mathcal{A}^{\mathrm{C}}\}$. Then $\mathcal{A}\cup\mathcal{B}$ is the simply connected hull of $\mathcal{A}$, which we denote $\mathcal{A}^{\mathrm{H}}$.

Lemma 4.4**.**

Let $\mathcal{A},\mathcal{B}\subset\mathbb{C}$ both be bounded open (path) connected sets. If there exist points $b_{\mathrm{in}},b_{\mathrm{out}}\in\operatorname{\mathrm{bd}}\mathcal{B}^{\mathrm{H}}$ such that $b_{\mathrm{in}}\in\mathcal{A}$ and $b_{\mathrm{out}}\in\mathcal{A}^{\mathrm{C}}$, then $\operatorname{\mathrm{bd}}\mathcal{A}\cap\operatorname{\mathrm{bd}}\mathcal{B}\neq\varnothing$.

Proof.

Since $\operatorname{\mathrm{bd}}\mathcal{B}^{\mathrm{H}}\subset\operatorname{\mathrm{bd}}\mathcal{B}$, we can assume that $b_{\mathrm{out}}$ is in the interior of $\mathcal{A}^{\mathrm{C}}$, as otherwise the proof is done. As $\mathcal{B}^{\mathrm{H}}$ is a bounded open simply connected set, $\operatorname{\mathrm{bd}}\mathcal{B}^{\mathrm{H}}$ must be connected; see, e.g., [Bur79, p. 345]. By way of contradiction, suppose that $\operatorname{\mathrm{bd}}\mathcal{A}\cap\operatorname{\mathrm{bd}}\mathcal{B}=\varnothing$. Then $\operatorname{\mathrm{bd}}\mathcal{B}^{\mathrm{H}}\subset\operatorname{\mathrm{bd}}\mathcal{B}\subset\mathcal{A}\cup\mathrm{int}(\mathcal{A}^{\mathrm{C}})$. However, since $b_{\mathrm{in}}\in\mathcal{A}$, $b_{\mathrm{out}}\in\mathrm{int}(\mathcal{A}^{\mathrm{C}})$, and $\mathcal{A}$ and $\mathrm{int}(\mathcal{A}^{\mathrm{C}})$ are both open and nonempty disjoint sets, it follows by definition of connected that $\operatorname{\mathrm{bd}}\mathcal{B}^{\mathrm{H}}$ is in fact disconnected, a contradiction. Thus, $\operatorname{\mathrm{bd}}\mathcal{A}\cap\operatorname{\mathrm{bd}}\mathcal{B}\neq\varnothing$ holds. ∎

Proof of Theorem 4.1.

If $\gamma=\mathcal{K}(A)^{-1}$, the proof is trivially satisfied with $\eta=0$, so assume that $\gamma\in(\mathcal{K}(A)^{-1},1)$. Since $\alpha(A)<0$, it follows that $\lim_{x\to 0^{+}}g(x,y)=\infty$ for any $y\in\mathbb{R}$, and so $g(x_{\star},y_{\star})=\mathcal{K}(A)^{-1}$ with $x_{\star}>0$. Now consider the strict lower level set $\mathcal{L}_{\gamma}\coloneqq\{(x,y):g(x,y)<\gamma,x>0\}$, which is clearly open and also bounded; see [Mit19, Theorem 2.3]. Let $\mathcal{L}$ be the (open) connected component of $\mathcal{L}_{\gamma}$ such that $(x_{\star},y_{\star})\in\mathcal{L}$ and let $\mathcal{G}\coloneqq\mathcal{L}^{\mathrm{H}}$, i.e., the simply connected hull of $\mathcal{L}$.

Now following the proof of [Gu00, Theorem 3.1] a bit more closely, by continuity of $g(x,y)$ and the fact that $\mathcal{G}$ is bounded, there must exist points $b_{1},b_{2}\in\operatorname{\mathrm{bd}}\mathcal{G}$

[TABLE]

such that

[TABLE]

where $\eta_{1},\eta_{2}>0$. Furthermore, we can assume that $\eta_{1}$ and $\eta_{2}$ are the smallest positive values such that (4.3) holds with $b_{1},b_{2}\in\operatorname{\mathrm{bd}}\mathcal{G}$. Noting that

[TABLE]

whose numerator has a GLC of 1, it follows for any $y_{1},y_{2}\in\mathbb{R}$ that

[TABLE]

i.e., $g(x_{\star},y)$ with respect to $y$ has a GLC of $\tfrac{1}{x_{\star}}$. By applying this GLC to (4.3), it follows that

[TABLE]

and so $\eta_{1}+\eta_{2}\geq 2x_{\star}(\gamma-\mathcal{K}(A)^{-1})$.

Now suppose that $\eta\in(0,2x_{\star}(\gamma-\mathcal{K}(A)^{-1})]$ so $\eta\leq\eta_{1}+\eta_{2}$. Obviously (4.1) is satisfied if $\eta=\eta_{1}+\eta_{2}$, so assume that $\eta<\eta_{1}+\eta_{2}$. Let

[TABLE]

i.e., $\mathcal{G}$ shifted downward by the amount $\eta$, and consider the line segment joining $b_{1}$ and $b_{2}$ and the point $b_{\mathrm{in}}=(x_{\star},y_{\star}+\eta_{2}-\eta)$. As $b_{\mathrm{in}}$ must be on this line segment, but not at its endpoints, $b_{\mathrm{in}}\in\mathcal{G}$ and $b_{\mathrm{in}}\in\operatorname{\mathrm{bd}}\mathcal{G}_{\eta}$, since $b_{2}\in\operatorname{\mathrm{bd}}\mathcal{G}$. Let $b_{\mathrm{out}}\in\operatorname{\mathrm{bd}}\mathcal{G}_{\eta}$ be a lowermost point of $\operatorname{\mathrm{bd}}\mathcal{G}_{\eta}$. Then, as $\eta>0$, $b_{\mathrm{out}}\in\mathcal{G}^{\mathrm{C}}$, the complement of $G$. Since $\mathcal{G}$ and $\mathcal{G}_{\eta}$ are both bounded open connected sets in the plane and $\mathcal{G}_{\eta}=\mathcal{G}_{\eta}^{\mathrm{H}}$, Lemma 4.4 applies, and so $\operatorname{\mathrm{bd}}\mathcal{G}\cap\operatorname{\mathrm{bd}}\mathcal{G}_{\eta}\neq\varnothing$. Letting $(\tilde{x},\tilde{y})$ be any such point in $\operatorname{\mathrm{bd}}\mathcal{G}\cap\operatorname{\mathrm{bd}}\mathcal{G}_{\eta}$, it follows that $(\tilde{x},\tilde{y}+\eta)\in\operatorname{\mathrm{bd}}\mathcal{G}$, hence $(\tilde{x},\tilde{y})$ satisfies (4.1). ∎

We now consider horizontally oriented pairs of points on the $\gamma$-level set of $g(x,y)$, similar to Theorem 2.1 for $f(x,y)$. As the horizontal orientation is actually more complicated for $\mathcal{K}(A)$, we first establish the following intermediate general result.

Lemma 4.5**.**

Let $f:(0,\infty)\to[0,\infty)$ be continuous with a GLC of $c\geq 0$ and consider the function $q(x)\coloneqq\tfrac{f(x)}{x}$ on the same domain. For $a,b>0$ with $b-a=\eta>0$, if $q(a)=q(b)=\gamma$ and $x_{\star}=\operatorname{\arg\min}_{x\in[a,b]}q(x)$ with $\gamma_{\star}=q(x_{\star})$, then

[TABLE]

Proof.

We assume that $\gamma_{\star}<\gamma$ as otherwise the inequality clearly holds. Since $f(x)$ has a GLC of $c$, it follows that

[TABLE]

Meanwhile

[TABLE]

Combining the two yields $2x_{\star}(\gamma-\gamma_{\star})-\eta\gamma\leq c\eta$, thus completing the proof. ∎

Theorem 4.6**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If $\mathcal{K}(A)^{-1}\leq\gamma$ and \eta\in\big{[}0,\tfrac{2x_{\star}}{1+\gamma}(\gamma-\mathcal{K}(A)^{-1})\big{]}, then there exists a pair $x,y\in\mathbb{R}$ with $x>0$ such that

[TABLE]

Corollary 4.7**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If there do not exist any pairs $x,y\in\mathbb{R}$ with $x>0$ such that (4.6) holds, then

[TABLE]

Proof of Theorem 4.6.

The beginning of the proof is the same as the first paragraph of the proof of Theorem 4.1. Again consider the simply connected set $\mathcal{G}$ defined there that contains $(x_{\star},y_{\star})$, a global minimizer. By continuity of $g(x,y)$ and the fact that $\mathcal{G}$ is bounded, there must exist points $b_{1},b_{2}\in\operatorname{\mathrm{bd}}\mathcal{G}$

[TABLE]

such that

[TABLE]

where $\eta_{1},\eta_{2}>0$. We again assume that $\eta_{1}$ and $\eta_{2}$ are the smallest positive values such that (4.8) holds with $b_{1},b_{2}\in\operatorname{\mathrm{bd}}\mathcal{G}$. Applying Lemma 4.5 to $q(x)\coloneqq g(x,y_{\star})=\tfrac{f(x,y_{\star})}{x}$ with $a=x_{\star}-\eta_{1}$ and $b=x_{\star}+\eta_{2}$, we have that

[TABLE]

as the numerator of $g(x,y)$ as rewritten in (4.4) has a GLC of 1.

Now suppose that \eta\in\big{(}0,\tfrac{2x_{\star}}{1+\gamma}(\gamma-\mathcal{K}(A)^{-1})\big{]} so $\eta\leq\eta_{1}+\eta_{2}$. If $\eta=\eta_{1}+\eta_{2}$, (4.6) is clearly satisfied, so instead assume that $\eta<\eta_{1}+\eta_{2}$. Considering the set $\mathcal{\widehat{G}}\coloneqq\{(x-\eta,y):(x,y)\in\mathcal{G}\},$ i.e., $\mathcal{G}$ shifted to the left by the amount $\eta$, the rest of the proof follows analogously to the end of the proof of Theorem 4.1. ∎

Although we have derived Kreiss constant analogues of Theorems 2.1 and 2.2, Corollaries 4.2 and 4.7 in fact assert that the $\tau(A,B)$ algorithms of [Gu00, BLO04, GMO*+*06] will not directly extend to Kreiss constants. The crux of the problem is that these $\tau(A,B)$ methods all rely on the fact that (2.4) holds when there are no points satisfying (2.3), which recall, provided a way of computing $\tau(A,B)$ via Gu’s verification procedure to verify an upper or lower bound for $\tau(A,B)$. However, in the Kreiss constant setting, our lower bounds given in (4.2) and (4.7) are not as concrete, as they depend on $x_{\star}$, which is unknown. For a given $\gamma$ and $\eta$, we do not even know if the lower bounds provided (4.2) and (4.7) would be meaningful, as they might not even be positively valued. Thus, to develop algorithms for $\mathcal{K}(A)$, crucial departures must be made. We will do this via two different strategies.

5 A continuous-time $\mathcal{K}(A)$ algorithm based on fixed-distance pairs

Let $(\hat{x},\hat{y})$ be a local (but not global) minimizer of $g(x,y)$ and set $\gamma\coloneqq g(\hat{x},\hat{y})$. As $(\hat{x},\hat{y})$ is a local minimizer, we do not need to verify that $\gamma\geq\mathcal{K}(A)^{-1}$, as this is obviously true. To obtain an optimization-with-restarts algorithm, we do not necessarily need to verify a lower bound either. Instead, we can just aim to detect other (non-stationary) points on the $\gamma$-level set of $g(x,y)$. Using such level-set points to restart optimization, a better (lower) minimizer of $g(x,y)$ is guaranteed to be found. Assuming $\gamma\in[\mathcal{K}(A)^{-1},1)$, Theorems 4.1 and 4.6 assert that if $\eta>0$ is chosen sufficiently small, there must exist points satisfying either (4.1) or (4.6). Of course, we do not know a priori how small to choose $\eta$, so we propose to use backtracking, i.e., we can start with $\eta$ initially set to some large value and simply decrease it in a loop until level-set points for restarting optimization are found. As long as $(\hat{x},\hat{y})$ is not a global minimizer, this backtracking procedure must succeed in finding level-set points for restarting optimization. Meanwhile, when $(\hat{x},\hat{y})$ is a global minimizer, and so $\mathcal{K}(A)$ has been computed, the backtracking procedure can simply be terminated once $\eta$ falls below a tolerance. A high-level pseudocode using this backtracking-based globality certificate is given in Algorithm 1. To complete this algorithm, we now must develop a corresponding 2D level-set test for continuous-time Kreiss constants. We will do this by looking for points satisfying (4.1) or (4.6), i.e., level-set points that are a fixed distance $\eta$ apart, and develop a procedure inspired by Gu’s 2D level-set test for $\tau(A,B)$.

5.1 A 1D vertical level-set test

Before we develop our 2D level-set test for $g(x,y)$, we will need the following theorem which will allow us to obtain all the points on the $\gamma$-level set of $g(x,y)$ along a chosen vertical line.

Theorem 5.1**.**

Given $\gamma,x,y\in\mathbb{R}$, with $\gamma\geq 0$ and $x\neq 0$, $\gamma$ is a singular value of $G(x,y)$ defined in (3.2) if and only if $\mathbf{i}y$ is an eigenvalue of the Hamiltonian matrix

[TABLE]

Proof.

It is clear that the matrix is Hamiltonian. Suppose $\gamma$ is a singular value $G(x,y)$ with left and right singular vectors $u$ and $v$, and so

[TABLE]

Rearranging terms and multiplying the bottom block row by $-1$, this is equivalent to

[TABLE]

∎

In fact, for a given $\gamma\geq\mathcal{K}(A)^{-1}$ and $x\coloneqq\hat{x}\neq 0$, computing the imaginary eigenvalues of (5.1) may provide more than the $\gamma$-level set points of $g(x,y)$ along the vertical line $x=\hat{x}$. This is because if $\mathbf{i}\hat{y}$ is an eigenvalue of (5.1), Theorem 5.1 asserts that $\gamma$ is a singular value of $F(\hat{x},\hat{y})$, but not necessarily the minimum one; if this happens, $(\hat{x},\hat{y})$ would be on a level set of $g(x,y)$ lower than the $\gamma$-level set and thus be an even better point for restarting optimization. Finally, note that structure-preserving eigensolvers [BMX98, BMX99] can be used for more reliable detection of imaginary eigenvalues of Hamiltonian matrices like (5.1).

5.2 A 2D level-set test for fixed-distance pairs

We now derive a new 2D level-set test for $g(x,y)$. Per Theorems 4.1 and 4.6, the choice of orientation for pairs of points on the $\gamma$-level set of $g(x,y)$ has consequences, hence we will develop our new continuous-time Kreiss constant procedure for arbitrary orientation. Specifically, given $\eta>0$ and angle $\theta\in(-\tfrac{\pi}{2},\tfrac{\pi}{2}]$, we will look for points a fixed distance $\eta$ apart of the form $(\hat{x},\hat{y})$ and $(\hat{x}+\eta\cos\theta,\hat{y}+\eta\sin\theta)$ such that $g(x,y)=\gamma$ holds at both of them and $\hat{x}>0$.

Suppose that $\gamma$ is a singular value of both $G(x,y)$ and $G(x+\eta\cos\theta,y+\eta\sin\theta)$, with respective left and right singular vectors pairs $u$,$v$ and $\hat{u}$,$\hat{v}$. By Theorem 5.1 for $G(x,y)$ and following a similar argument as in its proof for $G(x+\eta\cos\theta,y+\eta\sin\theta)$, the following two Hamiltonian matrices must share an imaginary eigenvalue $\mathbf{i}y$, i.e.,

[TABLE]

Let $A_{1}$ and $A_{2}$ respectively denote the square matrices in (5.2a) and (5.2b) and let $W=\begin{bmatrix}v\\ u\end{bmatrix}\begin{bmatrix}\hat{v}^{*}&\hat{u}^{*}\end{bmatrix}\neq 0$ so that we have $A_{1}W=\mathbf{i}yW$ and $A_{2}W^{*}=\mathbf{i}yW^{*}$. To eliminate $y$, we take the conjugate transpose of the second and then add the two together to obtain the Sylvester equation:

[TABLE]

Thus, if (5.3) also has a nonzero solution $W\in\mathbb{C}^{2n\times 2n}$, $A_{1}$ and $-A_{2}^{*}$ must share an eigenvalue. As $A_{1}$ and $A_{2}$ are Hamiltonian matrices, their spectra have imaginary-axis symmetry, and so $A_{1}$ and $A_{2}$ must also have an eigenvalue in common. Now separating out all terms involving $x$, we get

[TABLE]

Rewriting both Sylvester forms using the vectorize operator, and letting $w=\operatorname{\mathrm{vec}}(W)$, results in the generalized eigenvalue problem

[TABLE]

Our 2D level-set test for fixed-distance pairs works as follows. Given $\gamma\in[\mathcal{K}(A)^{-1},1)$, $\eta\geq 0$, and $\theta\in(-\tfrac{\pi}{2},\tfrac{\pi}{2}]$, we first compute all the eigenvalues of (5.5). If there are no (finite) positive real eigenvalues of (5.5), then the test is finished and returns no level-set points. Otherwise, additional calculations are done to ascertain whether or not any level-set points have been detected. For each eigenvalue $\hat{x}>0$ of (5.5), we have that the vertical line specified by $\hat{x}$, and, if $|\theta|<\tfrac{\pi}{2}$, the vertical line $\hat{x}+\eta\cos\theta$, may contain points on the $\gamma$-level set of $g(x,y)$. To determine this, for each of these vertical lines, say, $x\coloneqq\hat{x}$, we then apply Theorem 5.1 and compute all the eigenvalues of the corresponding Hamiltonian matrix (5.1). If this matrix has no imaginary eigenvalues, then no level-set points have been detected on this vertical line. Otherwise, all the points $(\hat{x},\hat{y})$ such that $\mathbf{i}\hat{y}$ is an imaginary eigenvalue of this matrix are added to the list of detected level-set points of $g(x,y)$ to return. Optionally, for each of these points, one could additionally check whether or not $\gamma$ is the minimum singular value of $G(\hat{x},\hat{y})$, but this is not strictly necessary; as discussed previously, $(\hat{x},\hat{y})$ must be on the $\gamma$-level set of $g(x,y)$ or a lower one, either of which suffice for restarting optimization to obtain a better (lower) minimizer of $g(x,y)$. If no level-set points are detected, for any of the vertical lines, then the test returns no points. Otherwise, all the detected level-set points are returned.

Our new 2D level-set test differs from the procedure in [Gu00] (and [GMO*+*06]) in a significant way; in 6.3, we explain how our modifications here greatly improve the reliability of these 2D level-set tests.

In terms of cost, in extreme situations there may be up to $\mathcal{O}(n^{2})$ potential vertical lines detected, which means that $\mathcal{O}(n^{2})$ Hamiltonian eigenvalue problems of dimension $2n\times 2n$ must be solved. However, when using standard dense eigensolvers, the overall work complexity of our procedure is actually $\mathcal{O}(n^{6})$, as we must first compute the eigenvalues of (5.5), which is a matrix pencil with square matrices of dimension $4n^{2}$. In terms of constant factors, if $A$ is real, $\mathcal{A}_{1}$ is real if and only if $\theta=0$.

5.3 Properties of the eigensystem $\mathcal{A}_{1}w=x\mathcal{A}_{2}w$ and its solution

Unlike the $4n^{2}\times 4n^{2}$ generalized eigenvalue problem (5.5), Gu derived a smaller $2n^{2}\times 2n^{2}$ generalized eigenvalue problem for $\tau(A,B)$ [Gu00, Eq. 3.13]. In [GMO*+*06, Section 3.1], this was then simplified further to a computationally easier $2n^{2}\times 2n^{2}$ standard eigenvalue problem [GMO*+*06, Eq. 3.7]. Via the following result, we show that the rank of $\mathcal{A}_{2}$ is $2n^{2}$, i.e., half its dimension.

Lemma 5.2**.**

Let $C\coloneqq\begin{bmatrix}aI&-bI\\ bI&-aI\end{bmatrix}\in\mathbb{C}^{2n\times 2n}$ with $a,b\in\mathbb{C}$ and $b\neq 0$. For $k\in\mathbb{N}$, the Kronecker sum $I_{2k}\otimes C+C\otimes I_{2k}=\mathcal{U}_{k}\mathcal{V}_{k}\in\mathbb{C}^{4kn\times 4kn}$ and $\mathcal{V}_{k}\mathcal{U}_{k}=2I_{k}\otimes C$, where

[TABLE]

Proof.

The factorization follows from the following if-and-only-if equivalences

[TABLE]

while $\mathcal{V}_{k}\mathcal{U}_{k}=I_{k}\otimes\begin{bmatrix}2aI&-bI\\ bI&0\end{bmatrix}+bI_{k}\otimes\begin{bmatrix}0&-I\\ I&-2ab^{-1}I\end{bmatrix}=2I_{k}\otimes C.$ ∎

Applying Lemma 5.2 to $\mathcal{A}_{2}$ with $a\coloneqq 1$, $b\coloneqq\gamma\neq 0$, and $k\coloneqq n$, $\mathcal{V}_{k}$ gives the reduced row echelon form of $\mathcal{A}_{2}$, and so the rank of $\mathcal{A}_{2}$ is $2n^{2}$. Thus, is it natural to ask if (5.5) can also be reduced to a $2n^{2}\times 2n^{2}$ generalized eigenvalue problem, and perhaps even a standard one. Unfortunately, the presence of nonzero off-diagonal blocks $\pm\gamma I$ in \big{[}\begin{smallmatrix}I&-\gamma I\\ \gamma I&-I\end{smallmatrix}\big{]} from $\mathcal{A}_{2}$ appears to prevent this. If one attempts to follow [Gu00, GMO*+*06] and similarly partition $W$ into four $n\times n$ blocks, multiplying out (5.4) for each block of $W$ results in four equations that all involve $x$. In contrast, in [Gu00, GMO*+*06], eigenvalue $x$ ($\alpha$ in their notation) only appears in the two corresponding equations for the diagonal blocks of $W$ ($X$ in their notation); for the off-diagonal blocks, eigenvalue $x$ (again $\alpha$ in their notation) does not appear in these other two equations, since as is seen in [GMO*+*06, Eq. 3.5], it ends up being multiplied by zero. Consequently, the reduction techniques of [Gu00, GMO*+*06] do not seem to be applicable to (5.5).

However, since $\mathcal{A}_{2}$ is singular, $\mathcal{A}_{1}-\lambda\mathcal{A}_{2}$ can at least be numerically deflated into a smaller pencil $\mathcal{\widetilde{A}}_{1}-\lambda\mathcal{\widetilde{A}}_{2}$ whose spectrum is the set of finite eigenvalues of (5.5); see the deflation routines of [BW20, Section 4.3] and [Köh21, Algorithm 3, Chapter 3.1], the former of which is implemented in the ml_ct_dss_adtf routine from MORLAB [BW19]. Although this iterative deflation technique is cubic work, and so also $\mathcal{O}(n^{6})$ work for (5.5), deflating matrix pencils this way can be faster than computing their eigenvalues with the QZ algorithm. Furthermore, for (5.5), deflation results in a matrix pencil of half the order, since it removes $2n^{2}$ infinite eigenvalues. As a result, deflating and then computing the eigenvalues of the resulting pencil may be even faster than computing the eigenvalues of (5.5) directly. That matrix $\mathcal{\widetilde{A}}_{2}$ in the reduced pencil is nonsingular also provides a second important benefit. Although we must compute the real eigenvalues of (5.5), we cannot expect its real eigenvalues to be exactly real in the presence of rounding errors. The key question then is how far away from the real axis can a computed eigenvalue be allowed to be while still being deemed a real eigenvalue of (5.5)? For the reduced pencil, a reliable tolerance is $\texttt{tol}\cdot\epsilon_{\mathrm{mach}}\cdot\|\mathcal{\widetilde{A}}_{2}^{-1}\mathcal{\widetilde{A}}_{1}\|_{\infty}$, where $\texttt{tol}>1$ is provided by the user and $\epsilon_{\mathrm{mach}}$ is the machine precision. Note that we still recommend computing the eigenvalues of $\mathcal{\widetilde{A}}_{1}-\lambda\mathcal{\widetilde{A}}_{2}$ as a generalized eigenvalue problem, instead of using the matrix $\mathcal{\widetilde{A}}_{2}^{-1}\mathcal{\widetilde{A}}_{1}$, as we have observed that the condition number of $\mathcal{\widetilde{A}}_{2}$ is generally very large in practice.

5.4 Faster computation of the real eigenvalues of $\mathcal{A}_{1}w=x\mathcal{A}_{2}w$

We now show how the $\mathcal{O}(n^{6})$ theoretical work complexity of Algorithm 1 can be reduced. To do this, we will work with the matrices in (5.5) and adapt the divide-and-conquer approach proposed in [GMO*+*06, Section 3.3.2] for faster computation of $\tau(A,B)$. At a high level, this efficiency improvement relies on two principles. First, although the matrices in (5.5) are $4n^{2}\times 4n^{2}$ in size, they arose from vectorizing the two corresponding $2n\times 2n$ Sylvester forms in (5.4). As such, applying $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$, or their inverses to a vector can actually be done with just $\mathcal{O}(n^{3})$ work. In turn, this means that for any shift $s\in\mathbb{C}$, $(\mathcal{A}_{1}-s\mathcal{A}_{2})^{-1}$ can be applied to a vector with $\mathcal{O}(n^{3})$ work (we will clarify how these computations are done in a moment). Consequently, a shift-and-invert eigenvalue solver, e.g., eigs in MATLAB, can be employed to find the eigenvalues of (5.5) that are closest to a shift $s$ with only $\mathcal{O}(n^{3})$ work. Second, given a matrix $\mathcal{X}\in\mathbb{C}^{q\times q}$, suppose one only wants its eigenvalues that are along a line segment, say, an interval $[0,D]$ on the real axis for some $D>0$. Then the recursive iteration given by [GMO*+*06, Algorithm 4], which uses a shift-and-invert eigensolver, can locate all eigenvalues of $\mathcal{X}$ in $[0,D]$ with at most $2q+1$ shifts in the worst case and $\mathcal{O}(\sqrt{q})$ shifts if the eigenvalues of $\mathcal{X}$ are distributed uniformly. For brevity, we forgo the details of describing [GMO*+*06, Algorithm 4], but note that the cost of choosing $D$ unnecessarily large only results in about four extra shifts [GMO*+*06, p. 490]. Thus, the theory says that by adapting this divide-and-conquer technique to compute the positive real eigenvalues of (5.5), the overall work complexity of Algorithm 1 will be reduced to $\mathcal{O}(n^{4})$ work on average and $\mathcal{O}(n^{5})$ in the worst case.

We now explain how $\mathcal{A}_{1}$, $\mathcal{A}_{2}$, and $(\mathcal{A}_{1}-s\mathcal{A}_{2})^{-1}$ can all be applied to a vector $w\in\mathbb{C}^{4n^{2}}$ with at most $\mathcal{O}(n^{3})$ work. As $\mathcal{A}_{2}$ has only $10n^{2}$ nonzero entries, it suffices to store it in a sparse matrix format. Computing $\mathcal{A}_{1}w$ can be done efficiently via vectorizing the first Sylvester form in (5.4), i.e., $\operatorname{\mathrm{vec}}\left(\begin{bmatrix}A&0\\ 0&-A^{*}\end{bmatrix}W+W\begin{bmatrix}A^{*}-\eta\mathrm{e}^{-\mathbf{i}\theta}I&-\gamma\eta\cos\theta I\\ \gamma\eta\cos\theta I&\eta\mathrm{e}^{\mathbf{i}\theta}I-A\end{bmatrix}\right),$ where $W\in\mathbb{C}^{2n\times 2n}$ and $w=\operatorname{\mathrm{vec}}(W)$. The dominant cost in obtaining $\mathcal{A}_{1}w$ is the two matrix multiplies with $W$, hence it too can be done in $\mathcal{O}(n^{3})$ work. For $y\in\mathbb{C}^{4n^{2}}$, we can efficiently obtain $w=(\mathcal{A}_{1}-s\mathcal{A}_{2})^{-1}y$ by considering $(\mathcal{A}_{1}-s\mathcal{A}_{2})w=y$. This “unvectorizes” into (5.3), provided that $x$ is replaced with $s$ and the zero in its right-hand side is replaced by $Y$, where $y=\operatorname{\mathrm{vec}}(Y)$. By solving the resulting Sylvester equation and vectorizing its solution $W$, $w=(\mathcal{A}_{1}-s\mathcal{A}_{2})^{-1}y$ is computed in $\mathcal{O}(n^{3})$ work.

6 Continuous-time $\mathcal{K}(A)$ algorithms based on variable-distance pairs

Having developed the first globally convergent iteration for continuous-time Kreiss constants, we now develop two more, namely Algorithms 1 and 2, which can be considered closer analogues of the two $\tau(A,B)$ methods of [BLO04] described in §2. As previously discussed at the end of §4, the lower bounds provided by Corollaries 4.2 and 4.7, due to their dependency on the unknown $x_{\star}$, prevent creating direct $\mathcal{K}(A)$ analogues of the trisection and optimization-with-restarts of [BLO04]. However, these lower bounds are the result of the assumption of looking for pairs of level-set points that are a fixed distance $\eta$ apart, as Gu originally used for Theorem 2.1. As we are about to show, if we change this assumption, then we can obtain different lower bound results than those given in (4.2) and (4.7).

Theorem 6.1**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If $\mathcal{K}(A)^{-1}\leq\gamma$ and $\eta\in[0,2(\gamma-\mathcal{K}(A)^{-1})]$, then there exists a pair $x,y\in\mathbb{R}$ with $x>0$ such that

[TABLE]

Corollary 6.2**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If there do not exist any pairs $x,y\in\mathbb{R}$ with $x>0$ such that (6.1) holds, then

[TABLE]

Proof of Theorem 6.1.

The proof follows similarly to the proof of Theorem 4.1 with the following modification. By using the variable distance $x\eta$, for $\gamma\in[\mathcal{K}(A)^{-1},1)$ we instead obtain that $\eta\in[0,2(\gamma-\mathcal{K}(A)^{-1})]$ is a sufficient condition for (6.1) to hold; this is because this choice cancels out the $x_{\star}$ in the proof, as $|y_{1}-y_{2}|=x_{\star}\eta$. ∎

Thus, by looking for vertical pairs of points that are this particular variable distance apart, i.e., $x\eta$, a corresponding certificate procedure would either assert that $\gamma\geq\mathcal{K}(A)^{-1}$ holds or that (6.2) does. Such a certificate would avoid the need for the backtracking procedure that was necessary for Algorithm 1. While this might seem to be obviously preferable to Algorithm 1, as a bit of foreshadowing, we note that the large eigenvalue problem that results for this variable-distance certificate is quite different from (5.5) and has its own downsides. We now describe this certificate to complete Algorithms 1 and 2.

6.1 A 2D level-set test for variable-distance pairs

Suppose $\gamma$ is both a singular value of $G(x,y)$ and $G(x,y+x\eta)$, with respective left and right singular vectors pairs $u$,$v$ and $\hat{u}$,$\hat{v}$. Via Theorem 5.1 for $G(x,y)$ and a similar argument as its proof for $G(x,y+x\eta)$, we have the following two Hamiltonian eigenvalue problems

[TABLE]

These two Hamiltonian matrices have an eigenvalue in common if

[TABLE]

has a nonzero solution $W\in\mathbb{C}^{2n\times 2n}$. Separating out the terms involving $x$, we have

[TABLE]

Rewriting both Sylvester forms using the vectorize operator, and letting $w=\operatorname{\mathrm{vec}}(W)$, we have the following generalized eigenvalue problem

[TABLE]

To check if (6.2) holds with this variable-distance certificate, we proceed similarly to our fixed-distance certificate from §5.2. Thus, after initially computing all the positive real eigenvalues of (6.6), we then apply Theorem 5.1 to each of these candidate vertical lines to see if we detect any points on a $\gamma$-level set (or lower) of $g(x,y)$.

Key Remark 6.3**.**

We now explain how the designs of our 2D level-set tests intentionally deviate from the $\tau(A,B)$ procedure Gu proposed in [Gu00, pp. 995–997] and how this results in much better reliability. If one were to more closely follow Gu’s procedure as it is written, for each positive real eigenvalue $\hat{x}>0$ of (6.6), one would instead compute the eigenvalues of the two matrices in (6.3) and then check whether these two spectra have any imaginary eigenvalues in common. If a shared imaginary eigenvalue is detected, then a pair of level-set points a distance $x\eta$ apart has been detected444Note that while the bottom of [Gu00, p. 996] says that detected points would be on the $\gamma$-level set, technically that only holds if $\gamma$ is also the minimum singular value at both of these points, which may or may not be true. and $\gamma\geq\mathcal{K}(A)^{-1}$ must hold. Otherwise, if none of the pairs of eigenvalue problems from (6.3) share imaginary eigenvalues, then following Gu would mean asserting that lower bound (6.2) must hold. While Gu’s procedure is sound in exact arithmetic, trying to assert whether two matrices share an (imaginary) eigenvalue is exceptionally difficult to do reliably in the presence of rounding errors. Furthermore, this matching also assumes that the real eigenvalues $\hat{x}$ have been computed accurately enough such that the matrices (6.3) would indeed share an eigenvalue $\mathbf{i}\hat{y}$, assuming $(\hat{x},\hat{y})$ satisfies (6.1), which is another source of numerical uncertainty; see [GMO*+*06, Section 5]. Suppose that there are indeed points a distance $x\eta$ apart vertically on the $\gamma$-level set of $g(x,y)$, but rounding errors prevent their detection. In this case, asserting that (6.2) holds may be erroneous. This can be a critical failure because the lower bound is only true if no such pair of level-set points exist, not if the procedure fails to detect them due to numerical problems! This is a major reason why Gu’s procedure can have such numerical difficulties, particularly when $\eta$ is small. Though [GMO*+*06] improves Gu’s procedure to make it faster, it too follows the same the idea of checking whether or not two matrices share imaginary eigenvalues and thus also inherits these numerical problems. However, in the course of our work here, we have realized that Gu’s $\tau(A,B)$ procedure ironically follows Theorems 2.1 and 2.2 too closely, because as it turns outs, checking whether or not level-set points are a distance $x\eta$ (or $\eta$) apart is entirely unnecessary. While indeed (6.2) must hold if no points satisfy (6.1), another sufficient condition for (6.2) to hold is that the procedure itself does not generate any level-set points whatsoever, on the $\gamma$-level set or lower and/or as pairs or single points. If any level-set points are detected, clearly $\gamma\geq\mathcal{K}(A)^{-1}$ holds, while none being generated implies (6.1) cannot hold, which in turn asserts (6.2) must hold. Note that our new way of performing these 2D level-set tests is not only a theoretical improvement; we initially designed our tests to check for common imaginary eigenvalues, like Gu’s procedure, but discovered that the aforementioned numerical issues were impeding the reliability of the codes. Applying our modifications to the $\tau(A,B)$ algorithms of [Gu00, BLO04, GMO*+*06] should improve their reliability as well.

6.2 Properties of the eigensystem $\mathcal{B}_{1}w=x\mathcal{B}_{2}w$ and its solution

Unlike (5.5), which can at least be numerically deflated to an order $2n^{2}$ generalized eigenvalue problem, we now show that (6.6) cannot be similarly reduced, as $\mathcal{B}_{2}$ is generally nonsingular. In fact, $\mathcal{B}_{2}$ has an explicit inverse.

Theorem 6.4**.**

Let $C\coloneqq\begin{bmatrix}aI&-bI\\ bI&-aI\end{bmatrix}\in\mathbb{C}^{2n\times 2n}$ with $a,b\in\mathbb{C}$ and $b\neq 0$, and let $\mathcal{M}\coloneqq I_{2k}\otimes C+C\otimes I_{2k}$. For $k\in\mathbb{N}$ and $s\in\mathbb{C}$, define matrix $\mathcal{D}\coloneqq sI_{4kn}+\mathcal{M}$ and scalar $\beta\coloneqq s^{2}+4(b^{2}-a^{2})$. Then $\mathcal{D}$ is invertible with inverse

[TABLE]

if and only if $s$ and $\beta$ are both nonzero, where $\mathcal{U}_{k}$ and $\mathcal{V}_{k}$ are defined in (5.6). Moreover, if $k=n$, then the following simpler formula holds

[TABLE]

Proof.

Via Lemma 5.2 applied to $\mathcal{M}$, we have that $\mathcal{D}=sI_{4kn}+\mathcal{U}_{k}\mathcal{V}_{k}$. As the eigenvalues of $C$ are $\pm\sqrt{a^{2}-b^{2}}$, by Theorem 1.3, the eigenvalues of $\mathcal{M}$ are $\pm 2\sqrt{a^{2}-b^{2}}$ and zero. Thus, $\mathcal{D}$ is invertible if and only if $-s$ is not equal to any of these eigenvalues, which is equivalent to $s\neq 0$ and $\beta\neq 0$. For (6.7), we apply the Sherman-Morrison-Woodbury formula to $(sI_{4kn}+\mathcal{U}_{k}\mathcal{V}_{k})^{-1}$ and use $\mathcal{V}_{k}\mathcal{U}_{k}=2I_{k}\otimes C$ from Lemma 5.2, yielding

[TABLE]

For (6.8), now with $k=n$, first note that $C^{2}=\phi I_{2n}$ with $\phi=(a^{2}-b^{2})$, while

[TABLE]

where we have used the mixed-product property of $\otimes$. Then by multiplying (6.8) with $sI_{4n^{2}}+\mathcal{M}$ and noting that $\beta-s^{2}+4\phi=0$, it follows that

[TABLE]

∎

Applying Theorem 6.4 to $\mathcal{B}_{2}$, with $a\coloneqq 1$, $b\coloneqq\gamma\neq 0$, $s\coloneqq-\mathbf{i}\eta\neq 0$, we see that $\beta\neq 0$ holds if $\eta\neq\pm 2\sqrt{\gamma^{2}-1}$ $(\not\in\mathbb{R}$ if $|\gamma|<1$), hence $\mathcal{B}_{2}$ is generically invertible with

[TABLE]

where $\beta=4(\gamma^{2}-1)-\eta^{2}$. On the upside, one could thus perform our variable-distance certificate by computing the eigenvalues of $\mathcal{B}_{2}^{-1}\mathcal{B}_{1}$, while our fixed-distance certificate requires solving the general eigenvalue problem (5.5). However, using $\mathcal{B}_{2}^{-1}\mathcal{B}_{1}$ can introduce numerical issues and the order of this eigenvalue problem is $4n^{2}$, whereas for the fixed-distance case, we could instead solve the numerically deflated order $2n^{2}$ problem discussed in §5.3. Furthermore, from (6.9), it is clear that $\|\mathcal{B}_{2}^{-1}\|\to\infty$ as $\eta\to 0$; in the presence of rounding errors, this may make it quite difficult to reliably ascertain which eigenvalues of $\mathcal{B}_{2}^{-1}\mathcal{B}_{1}$ should be considered real-valued.

At this point, one might ask if it would have been better to consider a variable-distance certificate using a horizontal orientation instead of a vertical one. We consider this in Appendix A and note here that it results in the matrix pencil $\mathcal{B}_{1}-\lambda\mathcal{\widetilde{B}}_{2}$, where, like $\mathcal{B}_{2}$, the matrix $\mathcal{\widetilde{B}}_{2}$ is nonsingular but becomes singular as $\eta\to 0$.

6.3 Adapting divide-and-conquer for $\mathcal{B}_{1}w=x\mathcal{B}_{2}w$

We now describe the computations needed for a divide-and-conquer version of our variable-distance certificate. Note that the numerical reliability of this approach might be different than our fixed-distance certificate, as (6.6) only has finite eigenvalues, while (5.5) has finite and infinite eigenvalues. For $w\in\mathbb{C}^{4n^{2}}$, $\mathcal{B}_{2}w$ can be efficiently obtained by storing $\mathcal{B}_{2}$ in a sparse format. Computing $\mathcal{B}_{1}w$ is $\mathcal{O}(n^{3})$ work, which is done by vectorizing the first matrix in (6.5), i.e., $\operatorname{\mathrm{vec}}\left(\begin{bmatrix}A&0\\ 0&-A^{*}\end{bmatrix}W+W\begin{bmatrix}A^{*}&0\\ 0&-A\end{bmatrix}\right)$, where $W\in\mathbb{C}^{2n\times 2n}$ and $w=\operatorname{\mathrm{vec}}(W)$. To obtain $w=(\mathcal{B}_{1}-s\mathcal{B}_{2})^{-1}y$ for $y\in\mathbb{C}^{4n^{2}}$, consider $(\mathcal{B}_{1}-s\mathcal{B}_{2})w=y$. This “unvectorizes” into (6.4) provided that $x$ is replaced by $s$ and the zero on its right-hand side is replaced by $Y$, where $y=\operatorname{\mathrm{vec}}(Y)$. Vectorizing the solution of the resulting Sylvester equation yields $w=(\mathcal{B}_{1}-s\mathcal{B}_{2})^{-1}y$ in $\mathcal{O}(n^{3})$ work.

7 Algorithms for discrete-time $\mathcal{K}(A)$

To adapt Algorithms 1, 1, and 2 to compute discrete-time Kreiss constants, we need to develop discrete-time versions of the 2D level-set tests from §5.2 and §6.1; the other necessary components for the two optimization-with-restarts algorithms have already been discussed in §3.2. We begin by developing discrete-time analogues of Theorems 4.1 and 4.2. For the discrete-time case, we now additionally assume that $0\not\in\Lambda(A)$.

Theorem 7.1**.**

For $A\in\mathbb{C}^{n\times n}$ with $\rho(A)<1$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(r_{\star},\theta_{\star})$ be a global minimizer of (3.7). If $\mathcal{K}(A)^{-1}\leq\gamma$ and \eta\in\big{(}0,\tfrac{2(r_{\star}-1)}{1+\gamma}(\gamma-\mathcal{K}(A)^{-1})\big{]}, then there exists an $r>1$ and $\theta\in[0,2\pi)$ such that

[TABLE]

Corollary 7.2**.**

For $A\in\mathbb{C}^{n\times n}$ with $\rho(A)<1$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(r_{\star},\theta_{\star})$ be a global minimizer of (3.7). If there do not exist any pairs $r,\theta\in\mathbb{R}$ with $r>1$ such that (7.1) holds, then

[TABLE]

Proof of Theorem 7.1.

If $\gamma=\mathcal{K}(A)^{-1}$, the proof is trivially satisfied with $\eta=0$, so assume that $\gamma\in(\mathcal{K}(A)^{-1},1)$. Since $\rho(A)<1$, it follows that $\lim_{r\to 1^{+}}h(r,\theta)=\infty$ for any $\theta\in\mathbb{R}$, and so $h(r_{\star},\theta_{\star})=\mathcal{K}(A)^{-1}$ with $r_{\star}>1$. Now consider the strict lower level set $\mathcal{L}_{\gamma}\coloneqq\{(r,\theta):h(r,\theta)<\gamma,r>1\}$, which is clearly open and also bounded; see [Mit19, Theorem 3.2]. Let $\mathcal{L}$ be the (open) connected component of $\mathcal{L}_{\gamma}$ such that $(r_{\star},\theta_{\star})\in\mathcal{L}$ and let $\mathcal{G}\coloneqq\mathcal{L}^{\mathrm{H}}$, i.e., the simply connected hull of $\mathcal{L}$.

By continuity of $h(r,\theta)$ and the boundedness of $\mathcal{G}$, there must exist $b_{1},b_{2}\in\operatorname{\mathrm{bd}}\mathcal{G}$

[TABLE]

such that

[TABLE]

where $\eta_{1},\eta_{2}>0$. Furthermore, we can assume that $\eta_{1}$ and $\eta_{2}$ are the smallest positive values such that (7.3) holds with $b_{1},b_{2}\in\operatorname{\mathrm{bd}}\mathcal{G}$. Defining

[TABLE]

$\tilde{h}(r)$ has a global minimizer $r_{\star}-1$ on domain $(0,\infty)$ and its numerator has a GLC of 1. Applying Lemma 4.5 to $q(r)\coloneqq\tilde{h}(r)$ with $a=r_{\star}-\eta_{1}-1$ and $b=r_{\star}+\eta_{2}-1$ yields

[TABLE]

Now suppose that \eta\in\big{(}0,\tfrac{2(r_{\star}-1)}{1+\gamma}(\gamma-\mathcal{K}(A)^{-1})\big{]} so $\eta\leq\eta_{1}+\eta_{2}$. If $\eta=\eta_{1}+\eta_{2}$, (7.1) is clearly satisfied, so instead assume that $\eta<\eta_{1}+\eta_{2}$. Considering the set $\mathcal{\widehat{G}}\coloneqq\{(r-\eta,\theta):(r,\theta)\in\mathcal{G}\}$, without loss of generality, we can assume $\theta_{\star}=0$, and so $\mathcal{\widehat{G}}$ is simply $\mathcal{G}$ shifted left by the amount $\eta$. Thus, as in the end of the proof of Theorem 4.6, the rest of the argument follows analogously to the last paragraph of the proof of Theorem 4.1. ∎

7.1 A 1D circular level-set test

For all of our discrete-time $\mathcal{K}(A)$ algorithms, we will also need the following theorem.

Theorem 7.3**.**

Given $\gamma,r,\theta\in\mathbb{R}$, with $\gamma\geq 0$ and $r\neq 1$, $\gamma$ is a singular value of $H(r,\theta)$ defined in (3.8) if and only if $\mathrm{e}^{\mathbf{i}\theta}$ is an eigenvalue of the symplectic matrix pencil

[TABLE]

Furthermore, if $A$ is invertible and $r\neq 0$, zero is not an eigenvalue of (7.5) and the matrix pencil is regular.

Proof.

It is easy to verify (7.5) is symplectic, and under the additional assumptions, also regular and that zero cannot be an eigenvalue. Now suppose $\gamma$ is a singular value $H(r,\theta)$ with left and right singular vectors $u$ and $v$, and so

[TABLE]

Rearranging terms, this is equivalent to

[TABLE]

∎

Similar to Theorem 5.1, note that the unimodular eigenvalues of (7.5) correspond to points that are either on the $\gamma$-level set of $h(r,\theta)$ or on lower level sets. The structure-preserving eigensolvers of [BMX98, BMX99] can also be used for more reliable detection of unimodular eigenvalues of symplectic pencils like (7.5).

7.2 Adapting Algorithm 1 for discrete-time $\mathcal{K}(A)$

To create our first discrete-time 2D level-set test, we will again look for pairs of points a fixed distance $\eta\geq 0$ apart, but now we will do this along rays from the origin, i.e., for $\hat{\theta}\in[0,2\pi)$, we check if $h(\hat{r},\hat{\theta})=h(\hat{r}+\eta,\hat{\theta})=\gamma$ holds for some $\hat{r}>1$. Suppose $\gamma$ is a singular value of both $H(r,\theta)$ and $H(r+\eta,\theta)$ with respective left and right singular vector pairs $u,v$ and $\hat{u},\hat{v}$. Then by Theorem 7.3, we have that

[TABLE]

which respectively we denote as $M-\lambda N$ and $\widetilde{M}-\lambda\widetilde{N}$. Now define $W=\begin{bmatrix}v\\ u\end{bmatrix}\begin{bmatrix}\hat{v}^{*}&\hat{u}^{*}\end{bmatrix}\neq 0$. Multiplying the two equations above from the right side, respectively by $\begin{bmatrix}\hat{v}^{*}&\hat{u}^{*}\end{bmatrix}\widetilde{M}^{*}$ and $\begin{bmatrix}v^{*}&u^{*}\end{bmatrix}N^{*}$, yields

[TABLE]

If we take the conjugate transpose of (7.7b) and then multiply it by $\mathrm{e}^{\mathbf{i}\theta}$, we obtain

[TABLE]

Subtracting (7.8) from (7.7a) yields

[TABLE]

so as $M-\lambda N$ and $\widetilde{N}^{*}-\lambda\widetilde{M}^{*}$ are both regular, [Chu87, Theorem 1] states that these two pencils must share an eigenvalue if $W\neq 0$ solves the equation above. Moreover if this shared eigenvalue is $\mathrm{e}^{\mathbf{i}\theta}$, then $\mathrm{e}^{-\mathbf{i}\theta}$ is an eigenvalue of $\widetilde{N}-\lambda\widetilde{M}$, which in turn implies that $\mathrm{e}^{\mathbf{i}\theta}$ is an eigenvalue of $\widetilde{M}-\lambda\widetilde{N}$. Thus, (7.9) having a nonzero solution $W$ is a necessary condition for the two pencils in (7.6) to have eigenvalue $\mathrm{e}^{\mathbf{i}\theta}$ in common.

We now want to separate out the $r$ terms of $MW\widetilde{M}^{*}$ and $NW\widetilde{N}^{*}$. First, we have:

[TABLE]

Then $MW\widetilde{M}^{*}$ is

[TABLE]

which is equal to

[TABLE]

Vectorizing the above equation, with $w=\operatorname{\mathrm{vec}}(W)$, yields

[TABLE]

which we will abbreviate as

[TABLE]

Likewise, $NW\widetilde{N}^{*}$ is

[TABLE]

which is equal to

[TABLE]

Similarly vectorizing this gives

[TABLE]

which we will abbreviate as

[TABLE]

Therefore, we finally have the following quadratic eigenvalue problem:

[TABLE]

which we will abbreviate as:

[TABLE]

Thus, to perform our fixed-distance discrete-time $\mathcal{K}(A)$ certificate for adapting Algorithm 1, we first compute all the real eigenvalues $r>1$ of (7.20), which specify a set of concentric circles centered at the origin on which we may find level set points. Then, for each of these candidate radii, we apply Theorem 7.3 to see if we indeed detect any points on a $\gamma$-level set (or lower) of $h(r,\theta)$. Like our continuous-time algorithms, this procedure is also $\mathcal{O}(n^{6})$ work (albeit with a larger constant term) when using dense (quadratic) eigensolvers.

Note that for a generic quadratic eigenvalue problem of the form (7.20), if either $\mathcal{Q}_{0}$ or $\mathcal{Q}_{2}$ are nonsingular, then the problem is well posed, i.e., it has at least one solution and not infinitely many; see [BDD*+*00, p. 283]. While for our quadratic eigenvalue problem $\mathcal{Q}_{2}$ is singular (see Lemma B.4 in Appendix B), we now show that $\mathcal{Q}_{0}$ is generically nonsingular, and thus (7.20) is guaranteed to be well posed.

Theorem 7.4**.**

Let $\gamma,\eta\in\mathbb{R}$ both be positive. Then $\mathcal{Q}_{0}$ from (7.20) is nonsingular if and only if $A$ is nonsingular and $\gamma$ is not a singular value of $A$.

Proof.

For any vector $w\in\mathbb{C}^{4n^{2}}$, suppose $\mathcal{Q}_{0}w=0$, i.e.,

[TABLE]

which holds if and only if

[TABLE]

where $w=\operatorname{\mathrm{vec}}(W)$. By [Chu87, Theorem 1], the generalized Sylvester equation above has a unique solution if and only if the two matrix pencils

[TABLE]

are both regular and have no eigenvalues in common. Clearly $W=0$ satisfies (7.21), so as long as zero is the only solution, $\mathcal{Q}_{0}$ is nonsingular for the assumptions. To prove the forward direction, we thus show that these three conditions hold.

We begin with the first pencil, which is regular if for at least one value of $\lambda$, the resulting matrix is nonsingular. Using $\lambda=-1$ results in \big{[}\begin{smallmatrix}A&-\gamma I\\ -\gamma I&A^{*}\end{smallmatrix}\big{]}. Since $A$ is invertible, this matrix is invertible if and only if its Schur complement with respect to $A$ is, i.e.,

[TABLE]

which by definition of singular values, holds if and only if $\gamma$ is not a singular value of $A$. Thus, the first pencil is regular. For the second pencil, choosing $\lambda=0$ results in matrix $\begin{bmatrix}\eta I&\gamma(\eta-1)I\\ 0&A\end{bmatrix}$, which is invertible since $\eta\neq 0$ and $A$ is invertible. Thus the second pencil is also regular. Furthermore, this argument also establishes that zero cannot be an eigenvalue of the second pencil.

We now show that the two pencils do not have any eigenvalues in common. We begin by noting that some of the eigenvalues of the first pencil are infinity, while all the eigenvalues of the second pencil are finite. Let $\lambda\neq 0$ be a finite eigenvalue of the first pencil with eigenvector $\begin{bmatrix}v\\ u\end{bmatrix}\neq 0$, hence

[TABLE]

Since $\gamma\neq 0$ and $A$ is invertible, $u=0$ if and only if $v=0$, hence neither are zero. This means they can be rescaled to each have unit norm and so $\gamma$ is a singular value of $A$, a contradiction, hence $\lambda=0$ must hold. As zero cannot be an eigenvalue of the second pencil, this part of the proof is complete.

For the reverse direction, first suppose that $A$ is singular. Then zero must be an eigenvalue of the second pencil, hence the two pencils share zero as an eigenvalue and so $\mathcal{Q}_{0}$ is singular. Now suppose that $\gamma$ is a singular value of $A$ with left and right singular vectors $\tilde{u}$ and $\tilde{v}$ and again consider (7.22). Since $\begin{bmatrix}v\\ u\end{bmatrix}\coloneqq\begin{bmatrix}\tilde{v}\\ \tilde{u}\end{bmatrix}\neq 0$ is in the nullspace of the matrix given in (7.22) for any $\lambda\in\mathbb{C}$, the first pencil matrix is not regular and so $\mathcal{Q}_{0}$ is singular. ∎

Remark 7.5**.**

It is also possible to derive a globality certificate based on arcs instead of radial segments, where for $\hat{r}>1$, angle(s) $\hat{\theta}$ satisfying $h(\hat{r},\hat{\theta})=h(\hat{r},\hat{\theta}+\eta)=\gamma$ are sought. We also considered this, but it resulted in a quadratic eigenvalue problem where both of the corresponding $\mathcal{Q}_{0}$ and $\mathcal{Q}_{2}$ matrices were singular, and so it was unclear if this alternative quadratic eigenvalue problem was well posed or not.

7.3 Adapting Algorithms 1 and 2 for discrete-time $\mathcal{K}(A)$

Now following the variable-distance certificate idea from §6, we derive a new discrete-time version to extend Algorithms 1 and 2 so that they can compute discrete-time Kreiss constants.

Theorem 7.6**.**

For $A\in\mathbb{C}^{n\times n}$ with $\rho(A)<1$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(r_{\star},\theta_{\star})$ be a global minimizer of (3.7). If $\mathcal{K}(A)^{-1}\leq\gamma$ and $\eta\in[0,2(\gamma-\mathcal{K}(A)^{-1})]$, then there exists an $r>1$ and $\theta\in[0,2\pi)$ such that

[TABLE]

where $\beta\coloneqq 1-\delta$ and $\delta\coloneqq-\tfrac{\eta}{1+\gamma}$.

Corollary 7.7**.**

For $A\in\mathbb{C}^{n\times n}$ with $\rho(A)<1$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(r_{\star},\theta_{\star})$ be a global minimizer of (3.7). If there do not exist any pairs $r,\theta\in\mathbb{R}$ with $r>1$ such that (7.23) holds, then

[TABLE]

Proof of Theorem 7.6.

The proof follows similarly to the proof of Theorem 7.1, except that (7.23) corresponds to a distance of $\tilde{\eta}\coloneqq(\beta r+\delta)-r=\tfrac{\eta(r-1)}{1+\gamma}$ between the two level-set points, which leads to cancellation, and so $\eta\in[0,2(\gamma-\mathcal{K}(A)^{-1})]$. ∎

For brevity, we skip showing the lengthy derivation of the resulting discrete-time variable-distance certificate (it follows similarly to §7.2), and instead just give the key parts necessary to perform the computation. Suppose $\gamma$ is both a singular value of $H(r,\theta)$ and $H(\beta r+\delta,y)$ with respective left and right singular vectors pairs $u$,$v$ and $\hat{u}$,$\hat{v}$. Applying Theorem 7.3 to $H(r,\theta)$ and $H(\beta r+\delta,\theta)$ yields (7.6a) and (7.6b) but now with $r$ replaced by $\beta r+\delta$. For this modified pair of sympletic eigenvalue problems, the certificate derivation results in the large quadratic eigenvalue problem

[TABLE]

Although $\mathcal{\widetilde{Q}}_{2}$ is singular, $\mathcal{\widetilde{Q}}_{0}$ is nonsingular for the assumptions of Theorem 7.4, since $\eta\coloneqq\delta\neq 0$. Thus, (7.25) is well posed under the same assumptions.

7.4 Adapting divide-and-conquer for discrete-time $\mathcal{K}(A)$

We now show how the real eigenvalues of (7.20) and (7.25) may be computed using divide-and-conquer. We begin with (7.20) and form its (companion) linearization

[TABLE]

where $z=\begin{bmatrix}rw\\ w\end{bmatrix}$. Assuming an a priori upper bound $D>0$ is known for all real eigenvalues of (7.20), divide-and-conquer will sweep the interval $[1,D]$ to find all the real-valued eigenvalues in this range. We now detail how the necessary operations with the matrices in (7.26) can all be done in at most $\mathcal{O}(n^{3})$ work.

Consider doing matrix-vector products with either matrix in (7.26) and a vector $w=\begin{bmatrix}w_{1}\\ w_{2}\end{bmatrix}$, where $w_{1},w_{2}\in\mathbb{C}^{2n^{2}}$. The nontrivial parts of these products are:

[TABLE]

which are equal to the following respective vectorizations:

[TABLE]

where $w_{1}=\operatorname{\mathrm{vec}}(W_{1})$ and $w_{2}=\operatorname{\mathrm{vec}}(W_{2})$. The first two of these can be obtained in $\mathcal{O}(n^{3})$ work since they only involve matrix-matrix products with $2n\times 2n$ matrices. The third can be obtained in $\mathcal{O}(n^{2})$ work as the number of nonzero entries in $\mathcal{Q}_{2}$ is simply $8n^{2}$, hence one should just store $\mathcal{Q}_{2}$ in a sparse format. This is also fortunate as when applying shift-and-invert to a generalized eigenvalue problem $Ax=\lambda Bx$, solvers such as eigs in MATLAB require that $B$ is provided explicitly, even when the operator $(A-sB)^{-1}$ is given implicitly as a function handle.

Given a shift $s\in\mathbb{C}$ and a vector $y=\begin{bmatrix}y_{1}\\ y_{2}\end{bmatrix}$, where $y_{1},y_{2}\in\mathbb{C}^{2n^{2}}$, we now focus on

[TABLE]

which is equivalent to

[TABLE]

The bottom block row provides:

[TABLE]

Multiplying the top block row of (7.30) by $s$ and then substituting in (7.31), we get

[TABLE]

Using (7.28) to obtain the four matrix-vector products above, we can then solve for $w_{1}$ via solving the following generalized continuous-time algebraic Sylvester equation:

[TABLE]

where $w_{1}=\operatorname{\mathrm{vec}}(W_{1})$, $\hat{y}=\operatorname{\mathrm{vec}}(\widehat{Y})$, and the matrix pairs $M,N$ and $\widetilde{M}^{*},\widetilde{N}^{*}$ are respectively given in (7.6a) and (7.10), all with $r$ replaced with $s$. Per [KW89], solving (7.33) can be done in $\mathcal{O}(n^{3})$ work. Finally, $w_{2}$ is obtained via (7.31).

For (7.25), only a few minor modifications to the divide-and-conquer variant we have just explained for (7.20) are necessary. As $\mathcal{\widetilde{Q}}_{0}$ is equal to $\mathcal{Q}_{0}$ with $\eta\coloneqq\delta$ and $\mathcal{\widetilde{Q}}_{2}=\beta\mathcal{Q}_{2}$, the first and third equations in (7.28) can be used to do the corresponding matrix-vector products. For $w_{1}\in\mathbb{C}^{2n^{2}}$, we have

[TABLE]

8 Numerical experiments

To validate our methods for computing Kreiss constants, we implemented proof-of-concepts in MATLAB. The supplementary material includes the code, test examples, and a detailed description of the setup in order to reproduce the experiments in this paper. We plan to add “production-ready” implementations of our methods to a future release of ROSTAPACK [Mit].

We considered two $10\times 10$ stable matrices based on demos in EigTool [Wri02], one for continuous-time $\mathcal{K}(A)$ and a second for discrete-time $\mathcal{K}(A)$. For the former, we used $A=B-\kappa I$, where $B=\texttt{companion\_demo(10)}$ and $\kappa=1.001\alpha(B)$; this matrix has a large Kreiss constant, and as shown in Fig. 1, (1.9) has two local maximizers for this example. For the latter, we chose $A=\tfrac{1}{13}B+\tfrac{11}{10}I$, where $B=\texttt{convdiff\_demo(11)}$; while this matrix a small Kreiss constant, it is interesting for testing as, for this example, $h(r,\theta)$ has several local minimizers (see the supplementary material) and, per Fig. 1, (1.4) appears to be nonsmooth (though not at the maximizer). To verify restarting, we intentionally chose starting points such that global minimizers of $g(x,y)$ and $h(r,\theta)$ would not be found in the first round of optimization. In Table 1, we provide detailed metrics on our three algorithms and see in practice that the optimization-based Algorithms 1 and 1 are much faster and more accurate than the trisection-based Algorithm 2. The much higher numerical accuracy of our two optimization methods is also verified by comparing to Chebfun, with Chebfun being much slower, even though it was given a small interval containing a global maximizer attaining $\mathcal{K}(A)$. In fact, Chebfun is also much slower than Algorithm 2, which is the slowest of our three new methods, due to the linear convergence of trisection. As can be seen, Algorithm 1 appears to be the best of the algorithms. See the supplementary material for additional figures showing how Algorithms 1 and 1 progress one from minimizer to the next to converge to $\mathcal{K}(A)$.

We also evaluated our theoretically faster divide-and-conquer approaches. The supplementary material includes a demo for reproducing one key experiment of the many we performed. This demo uses the computations described in §5.4 and §6.3 to respectively compute eigenvalues of (5.5) and (6.6) closest to a given shift, where $\gamma$ and $\eta$ were set to the values used in the first 2D level-set test computed by Algorithm 1 when computing $\mathcal{K}(A)$ of our continuous-time example. However, we observed that these shift-and-invert closest eigenvalue computations were somewhat unreliable when compared to using eig in MATLAB on (5.5) and (6.6) directly. For space reasons, we forgo many of the details and instead make a few key observations. First, sometimes the eigenvalues computed by the shift-and-invert approach did not have even a single digit of agreement to those computed by eig, but as these eigenvalues were often very close to the origin, this is not so surprising. Second, we also observed that eigs in MATLAB would sometimes fail to return one of the closest eigenvalues to a given shift; interestingly, this phenomenon can even be observed in the paper where divide-and-conquer for $\tau(A,B)$ was proposed [GMO*+*06, top left plot of Fig. 3.1]. Third, we also observed that the large eigenvalue problems sometimes have very close conjugate pair eigenvalues, which may explain some of the aforementioned difficulties with extracting eigenvalues on the real axis accurately. As divide-and-conquer is built on the assumption that shift-and-invert eigensolvers are reliable, addressing the numerical issues above may just be a matter of picking the right sparse eigensolver. That being said, our initial testing of divide-and-conquer and its numerical reliability is merely a starting point for future work.

9 Concluding remarks

We have presented the first globally convergent algorithms for computing continuous- and discrete-time Kreiss constants to arbitrary accuracy. Algorithm 1, which uses our variable-distance 2D level-set tests, appears to be the fastest of the three different iterations proposed here, since it requires performing the least number of these tests. However, for continuous-time Kreiss constants, upcoming fast deflation methods [Köh21, Chapter 3.1] may enable our fixed-distance 2D level-set test of §5.2 to be computed with significantly less work, which could make Algorithm 1 more competitive, perhaps even faster than Algorithm 1. While we have outlined the key theoretical differences between the various large eigenvalue problems underlying our fixed- and variable-distance certificates, a systematic numerical evaluation could help illuminate whether there are any notable differences in terms of numerical reliability when using dense eigensolvers. Similarly, further investigation into the numerics of our asymptotically faster divide-and-conquer variants using sparse shift-and-invert eigensolvers could prove useful.

Acknowledgments

Many thanks to Mark Embree for suggesting to look at the problem of computing Kreiss constants and to Michael L. Overton for pointing out the possible connection to computing the distance to uncontrollability and hosting the author at the Courant Institute in New York for several visits. The author is also grateful for the many useful suggestions from the referees, one of whom deserves special mention for pointing out both a more concise proof of Lemma 4.5 and the simpler explicit inverse formulas for $\mathcal{B}_{2}$ and $\mathcal{\widetilde{B}}_{2}$ respectively given by (6.9) and (A.6).

Appendix A A 2D level-set test for variable-distance pairs with horizontal orientation

If we modify Theorem 4.6 to consider horizontal pairs of points a variable distance $\tfrac{x\eta}{1+\gamma}$ apart (instead a fixed distance $\eta$ apart), to induce cancellation in the proof, we obtain the following results.

Theorem A.1**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If $\mathcal{K}(A)^{-1}\leq\gamma$ and $\eta\in[0,2(\gamma-\mathcal{K}(A)^{-1})]$, then there exists a pair $x,y\in\mathbb{R}$ with $x>0$ such that

[TABLE]

where $\beta\coloneqq 1+\tfrac{\eta}{1+\gamma}$.

Corollary A.2**.**

For $A\in\mathbb{C}^{n\times n}$ with $\alpha(A)<0$, let $\gamma\in[0,1)$, $\eta\geq 0$, and $(x_{\star},y_{\star})$ be a global minimizer of (3.1). If there do not exist any pairs $x,y\in\mathbb{R}$ with $x>0$ such that (A.1) holds, then

[TABLE]

Proof of Theorem A.1.

The proof follows similarly to the proof of Theorem 4.6, except that (A.1) corresponds to a distance of $\tilde{\eta}\coloneqq\beta x-x=\tfrac{\eta x}{1+\gamma}$ between the two level-set points, which leads to cancellation, and so $\eta\in[0,2(\gamma-\mathcal{K}(A)^{-1})]$. ∎

The derivation of the corresponding verification procedure, to enable versions of Algorithms 1 and 2 using (A.1), is as follows. Suppose $\gamma$ is both a singular value of $G(x,y)$ and $G(\beta x,y)$ with respective left and right singular vectors pairs $u$,$v$ and $\hat{u}$,$\hat{v}$. Applying Theorem 5.1 to $G(x,y)$ and $G(\beta x,y)$ yields the following two Hamiltonian standard eigenvalue problems

[TABLE]

which have a common eigenvalue if

[TABLE]

has a nonzero solution $W\in\mathbb{C}^{2n\times 2n}$. Separating this into two Sylvester forms to isolate $x$ and then vectorizing yields the generalized eigenvalue problem

[TABLE]

$\mathcal{B}_{1}$ is defined in (6.6), and $w=\operatorname{\mathrm{vec}}(W)$. Applying Theorem 1.3 to $\mathcal{\widetilde{B}}_{2}$, we see that its eigenvalues are the pairwise sums of $\pm\sqrt{1-\gamma^{2}}$ and $\pm\beta\sqrt{1-\gamma^{2}}$, hence it is nonsingular if and only if $\beta\neq\pm 1$ and $\gamma\neq\pm 1$. Furthermore, via Lemma B.3,

[TABLE]

where $\phi\coloneqq(1+\beta^{2})(1-\gamma^{2})$ and $\omega\coloneqq(1-\beta^{2})^{-2}(1-\gamma^{2})^{-2}$. However, as $\eta\to 0$, $\beta\to 1$, and so $\mathcal{\widetilde{B}}_{2}$ will become closer and closer to being singular as the algorithms progress. This suggests that it would be better numerically to solve (A.5) as a generalized eigenvalue problem, rather than solving the standard eigenvalue problem $\mathcal{\widetilde{B}}_{2}^{-1}\mathcal{B}_{1}$. In terms of the number of arithmetic operations, if $A$ is real, using divide-and-conquer to compute the positive real eigenvalues of (A.5) only involves matrix problems with real coefficients. The two eigenvalue problems in (A.3) also have real coefficients when $A$ is real. Given $s\in\mathbb{R}$, $y=(\mathcal{B}_{1}-s\mathcal{\widetilde{B}}_{2})^{-1}w$ can be computed by solving the Sylvester equation (A.4) with $x$ replaced by $s$ and its zero right hand side replaced by $Y$, where $w=\operatorname{\mathrm{vec}}(W)$ and $y=\operatorname{\mathrm{vec}}(Y)$. In contrast, the computations for the divide-and-conquer variant from §6.3 requires complex arithmetic.

Appendix B Some technical results

Lemma B.1**.**

Let $\{\gamma_{j}\}$ be the iterates produced by trisection with $\{\gamma_{j}\}\to\gamma_{\star}$, and suppose that trisection terminates at the $k$th iterate with $|\gamma_{\star}-\gamma_{k}|\leq\psi\gamma_{\star}$ holding for some given relative error tolerance $\psi>0$. Then $\eta_{k}\leq(1+\psi)\gamma_{\star}$.

Proof.

By construction, $\gamma_{k}=L+\eta_{k}$ and $L\geq 0$ so $\eta_{k}\leq\gamma_{k}$ always holds. If $\gamma_{k}\geq\gamma_{\star}$, then $\psi\gamma_{\star}\geq|\gamma_{\star}-\gamma_{k}|=\gamma_{k}-\gamma_{\star}$, which implies $\eta_{k}\leq(1+\psi)\gamma_{\star}$. Otherwise, $\gamma_{k}<\gamma_{\star}$ must hold, and so it follows that $\psi\gamma_{\star}\geq|\gamma_{\star}-\gamma_{k}|=\gamma_{\star}-\gamma_{k}>\gamma_{k}-\gamma_{\star}$. ∎

Corollary B.2**.**

Given a fixed $\psi>0$, if $\eta_{k}>(1+\psi)\gamma_{\star}$, then the relative error of trisection at the $k$th iterate is bounded below by $\psi$, specifically $|\gamma_{\star}-\gamma_{k}|>\psi\gamma_{\star}$.

Lemma B.3**.**

Let $C\coloneqq\begin{bmatrix}aI&-bI\\ bI&-aI\end{bmatrix}\in\mathbb{C}^{2n\times 2n}$ with $a,b\in\mathbb{C}$ and let $\beta\in\mathbb{C}$. Then the matrix $\mathcal{D}\coloneqq I_{2n}\otimes C+\beta C\otimes I_{2n}$ is invertible if and only if $a^{2}\neq b^{2}$ and $\beta^{2}\neq 1$, where

[TABLE]

$\phi\coloneqq(1+\beta^{2})(a^{2}-b^{2})$, and $\omega\coloneqq(1-\beta^{2})^{-2}(a^{2}-b^{2})^{-2}$.

Proof.

First note that $C^{2}=(a^{2}-b^{2})I_{2n}$, while using the mixed-product property of $\otimes$ yields

[TABLE]

Using this equivalence for $\mathcal{D}^{2}$, the formula for $\mathcal{D}^{-1}$ is verified via

[TABLE]

which, by noting $\phi^{2}-4\beta^{2}(a^{2}-b^{2})^{2}=\omega^{-1}$, is equivalent to $I_{4n^{2}}$. ∎

Lemma B.4**.**

Matrix $\mathcal{Q}_{2}$ from (7.20) is singular.

Proof.

Suppose that $\mathcal{Q}_{2}w=0$ if and only if $w=0$, and so

[TABLE]

where the equivalence holds by unvectorizing the first equation and $w=\operatorname{\mathrm{vec}}(W)$. This generalized Sylvester equation has a unique solution if and only if

[TABLE]

are both regular matrix pencils and have no eigenvalues in common. However, as $\begin{bmatrix}0&\gamma I\\ 0&I\end{bmatrix}$ and $\begin{bmatrix}I&\gamma I\\ 0&0\end{bmatrix}$ are both singular, the two pencils share zero as an eigenvalue. ∎