Computing Unstructured and Structured Polynomial Pseudospectrum Approximations

TL;DR

This paper introduces a new efficient method for approximating the pseudospectra of matrix polynomials using rank-one perturbations, improving accuracy over random perturbation methods for both structured and unstructured cases.

Contribution

The paper presents a novel approach leveraging rank-one perturbations inspired by Wilkinson's analysis to compute polynomial pseudospectra more efficiently and accurately.

Findings

Method outperforms random perturbation approaches in accuracy

Effective for both structured and unstructured pseudospectra

Computational efficiency is significantly improved

Abstract

In many applications it is important to understand the sensitivity of eigenvalues of a matrix polynomial to perturbations of the polynomial. The sensitivity commonly is described by condition numbers or pseudospectra. However, the computation of pseudospectra of matrix polynomials is very demanding computationally. This paper describes a new approach to computing approximations of pseudospectra of matrix polynomials by using rank-one or projected rank-one perturbations. These perturbations are inspired by Wilkinson's analysis of eigenvalue sensitivity. This approach allows the approximation of both structured and unstructured pseudospectra. Computed examples show the method to perform much better than a method based on random rank-one perturbations both for the approximation of structured and unstructured (i.e., standard) polynomial pseudospectra.

Figures

Tables

Table 1. Table 1: Example 1: Eigenvalue condition numbers.

$i$	${\lambda }_i$	$\kappa \left({\lambda }_i\right)$	${\kappa }^{𝒮}\left({\lambda }_i\right)$
$1$	$-1.6907$	$23.2593$	$7.0577$
$2$	$-0.9225+1.1935i$	$5.9741$	$1.8875$
$3$	$-0.9225-1.1935i$	$5.9741$	$1.8875$
$4$	$0.5245+1.3668i$	$34.2042$	$11.5406$
$5$	$0.5245-1.3668i$	$34.2042$	$11.5406$
$6$	$0.4113+0.7192i$	$17.4605$	$8.3749$
$7$	$0.4113-0.7192i$	$17.4605$	$8.3749$
$8$	$0.6637$	$18.3210$	$9.8822$
$9$	$0.2045$	$7.4414$	$7.3777$
$10$	$-0.5701$	$6.2696$	$3.7923$

Table 2. Table 2: Example 2: Eigenvalue condition numbers.

$i$	${\lambda }_i$	$\kappa \left({\lambda }_i\right)$
$1$	$-0.8848+8.4415i$	$27.2147$
$2$	$-0.8848-8.4415i$	$27.2147$
$3$	$0.0947+2.5229i$	$0.9276$
$4$	$0.0947-2.5229i$	$0.9276$
$5$	$-0.9180+1.7606i$	$2.3301$
$6$	$-0.9180-1.7606i$	$2.3301$

Equations

$$P(\lambda )=A_m{\lambda }^m+A_{m-1}{\lambda }^{m-1}+\dots +A_1\lambda +A_0,$$

$$A_{0}\mbox{\boldmath{$x$}}=-\lambda A_{1}\mbox{\boldmath{$x$}},$$

$$\Lambda (P)=\left\{\lambda \in \mathbb{C}:\det (P(\lambda ))=0\right\}.$$

$$\mathcal{A}(P,\epsilon ,\Delta ,\omega )=\left\{\sum _{j=0}^m(A_j+\epsilon {\Delta }_j){\lambda }^j:\Vert {\Delta }_j{\Vert }_F\le {\omega }_j,j=0,\dots ,m\right\}.$$

$${\Lambda }_{\epsilon }(P)=\left\{z\in \Lambda (Q):Q\in \mathcal{A}(P,\epsilon ,\Delta ,\omega )\right\}.$$

$$\kappa(\lambda)=\frac{\omega(|\lambda|)}{|\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}|},$$

$$\Delta_{j}=\eta\omega_{j}\mathrm{e}^{-\mathrm{i}j\arg(\lambda)}\mbox{\boldmath{$y$}}\mbox{\boldmath{$x$}}^{H},\qquad j=0,\dots,m\,,$$

$$\sum_{j=0}^{m}\Delta_{j}\lambda^{j}(\varepsilon)\mbox{\boldmath{$x$}}(\varepsilon)+\sum_{j=1}^{m}(A_{j}+\epsilon\Delta_{j})j\lambda^{j-1}(\varepsilon)\lambda^{\prime}(\varepsilon)\mbox{\boldmath{$x$}}(\varepsilon)+\sum_{j=0}^{m}(A_{j}+\epsilon\Delta_{j})\lambda^{j}(\varepsilon)\mbox{\boldmath{$x$}}^{\prime}(\varepsilon)=\mbox{\bf 0}.$$

$$\sum_{j=0}^{m}\Delta_{j}\lambda^{j}\mbox{\boldmath{$x$}}+\sum_{j=1}^{m}A_{j}j\lambda^{j-1}\lambda^{\prime}(0)\mbox{\boldmath{$x$}}+\sum_{j=0}^{m}A_{j}\lambda^{j}\mbox{\boldmath{$x$}}^{\prime}(0)=\mbox{\bf 0},$$

$$P^{\prime}(\lambda)\lambda^{\prime}(0)\mbox{\boldmath{$x$}}=-P(\lambda)\mbox{\boldmath{$x$}}^{\prime}(0)-\sum_{j=0}^{m}\Delta_{j}\lambda^{j}\mbox{\boldmath{$x$}}.$$

$$\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}\cdot\lambda^{\prime}(0)=-\mbox{\boldmath{$y$}}^{H}P(\lambda)\mbox{\boldmath{$x$}}^{\prime}(0)-\mbox{\boldmath{$y$}}^{H}\sum_{j=0}^{m}\Delta_{j}\lambda^{j}\mbox{\boldmath{$x$}},$$

$$\lambda^{\prime}(0)=-\frac{\mbox{\boldmath{$y$}}^{H}\sum_{j=0}^{m}\Delta_{j}\lambda^{j}\mbox{\boldmath{$x$}}}{\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}}.$$

$$|\lambda^{\prime}(0)|=\frac{|\mbox{\boldmath{$y$}}^{H}(\sum_{j=0}^{m}\Delta_{j}\lambda^{j})\mbox{\boldmath{$x$}}|}{|\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}|}\leq\frac{\omega(|\lambda|)}{|\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}|},$$

$${\Lambda }_{{\omega }_m\epsilon }(A_m):=\left\{z\in \Lambda (A_m+E),E\in {\mathbb{C}}^{n\times n},\Vert E{\Vert }_F\le {\omega }_m\epsilon \right\}.$$

$$\epsilon <\underset{1\le i\le n}{\min }\frac{|{\lambda }_i(A_m)|}{\hat{\kappa }({\lambda }_i(A_m)){\omega }_m},$$

$${\epsilon }_{\ast }=\inf \{\Vert P(\lambda )-Q(\lambda ){\Vert }_F:Q(\lambda )\in {\mathbb{C}}^{n\times n}\text{mbox}isdefective\}.$$

$$\epsilon :=\underset{\begin{array}{c}1\le i\le mn\\ 1\le j\le mn\\ j\ne i\end{array}}{\min }\frac{|{\lambda }_i-{\lambda }_j|}{\kappa ({\lambda }_i)+\kappa ({\lambda }_j)}.$$

$$M{|}_{\hat{\mathcal{S}}}:=\frac{M{|}_{\mathcal{S}}}{\Vert M{|}_{\mathcal{S}}{\Vert }_F}$$

$${\mathcal{A}}^{\mathcal{S}}(P,\epsilon ,\omega ,\Delta )=\left\{\sum _{j=0}^m(A_j+\epsilon {\Delta }_j){\lambda }^j:{\Delta }_j\in {\mathcal{S}}_j,\Vert {\Delta }_j{\Vert }_F\le {\omega }_j,j=0,\dots ,m\right\}.$$

$$\kappa^{{\mathcal{S}}}(\lambda)=\frac{\omega^{{\mathcal{S}}}(|\lambda|)}{|\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}|},$$

$$\omega^{{\mathcal{S}}}(z)=\sum_{j=0}^{m}\|\mbox{\boldmath{$y$}}\mbox{\boldmath{$x$}}^{H}|_{{{\mathcal{S}}}_{j}}\|_{F}\,\omega_{j}z^{j}\,.$$

$$\Delta^{{\mathcal{S}}}_{j}=\eta\omega_{j}\mathrm{e}^{-\mathrm{i}j\arg(\lambda)}\mbox{\boldmath{$y$}}\mbox{\boldmath{$x$}}^{H}|_{\widehat{{\mathcal{S}}}_{j}},\qquad j=0,\dots,m,$$

$$|\lambda^{\prime}(0)|=\frac{|\mbox{\boldmath{$y$}}^{H}(\sum_{j=0}^{m}\Delta_{j}\lambda^{j})\mbox{\boldmath{$x$}}|}{|\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}|}=\frac{|\sum_{j=0}^{m}(\mbox{\boldmath{$y$}}^{H}\Delta_{j}\mbox{\boldmath{$x$}})\lambda^{j}|}{|\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}|},$$

$$|\lambda^{\prime}(0)|=\frac{\omega^{{\mathcal{S}}}(|\lambda|)}{|\mbox{\boldmath{$y$}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{$x$}}|}.$$

$$P(\lambda )=M{\lambda }^2+C\lambda +K,$$

$${\Lambda }_{\epsilon }^{\mathcal{S}}(P)=\left\{z\in \Lambda (Q):Q\in {\mathcal{A}}^{\mathcal{S}}(P,\epsilon ,\omega ,\Delta )\right\}.$$

$${\Lambda }_{{\omega }_m\epsilon }^{{\mathcal{S}}_m}(A_m):=\left\{z\in \Lambda (A_m+E),E\in {\mathcal{S}}_m,\Vert E{\Vert }_F\le {\omega }_m\epsilon \right\}.$$

$$\epsilon <\underset{1\le i\le n}{\min }\frac{|{\lambda }_i(A_m)|}{{\hat{\kappa }}_{{\mathcal{S}}_m}({\lambda }_i(A_m)){\omega }_m}.$$

$${\epsilon }_{\ast }^{\mathcal{S}}=\inf \{\Vert P(\lambda )-Q(\lambda ){\Vert }_F:Q(\lambda )\in \mathcal{S}\text{mbox}isdefective\}.$$

$${\epsilon }^{\mathcal{S}}:=\underset{\begin{array}{c}1\le i\le mn\\ 1\le j\le mn\\ j\ne i\end{array}}{\min }\frac{|{\lambda }_i-{\lambda }_j|}{{\kappa }^{\mathcal{S}}({\lambda }_i)+{\kappa }^{\mathcal{S}}({\lambda }_j)}\ge \epsilon .$$

Full text

1 Introduction

In many problems in science and engineering it is important to know the sensitivity of the eigenvalues of a square matrix to perturbations. The pseudospectrum is an important aid for shedding light on the sensitivity. Many properties and applications of the pseudospectrum of a matrix are discussed by Trefethen and Embree [23]; see also [6, 7, 10, 13, 20]. However, the computation of pseudospectra is a computationally demanding task except for very small matrices. Therefore, the development of numerical methods for the efficient computation of pseudospectra of medium-sized matrices, or partial pseudospectra of large matrices, has received considerable attention; see [3, 15, 18, 19, 25, 26].

The present paper is concerned with the computation of pseudospectra of matrix polynomials of the form

[TABLE]

where $\lambda\in\mathbb{C}$ and $A_{j}\in\mathbb{C}^{n\times n}$, $j=0,\ldots,m$. We will assume that $\det(A_{m})\neq 0$. Then $P$ has $mn$ finite eigenvalues, i.e., there are no eigenvalues at infinity. Matrix polynomials of this kind arise in many applications in systems and control theory; see, e.g., [8, 9, 14]. The case $m=1$ corresponds to the generalized eigenvalue problem

[TABLE]

and the special case $A_{1}=-I_{n}$ yields a standard eigenvalue problem. Here and throughout this paper $I_{n}$ denotes the identity matrix of order $n$. In some applications the matrices $A_{j}$ in (1.1) have a structure that should be respected, such as being symmetric, skew-symmetric, banded, Toeplitz, or Hankel.

The sensitivity of the eigenvalues of a matrix polynomial (1.1) to perturbations in the matrices $A_{j}$ is important in applications. This question therefore has received considerable attention; see, e.g., [4, 10, 12, 21, 22] and references therein. When the matrices $A_{j}$ are structured, it is natural to only consider perturbations that are similarly structured.

Define the spectrum of $P$,

[TABLE]

Given a set of matrices $\Delta=\{\Delta_{0},\ldots,\Delta_{m}\}$, $\Delta_{j}\in\mathbb{C}^{n\times n}$, and a set of weights $\omega=\{\omega_{0},\ldots,\omega_{m}\}$, $\omega_{j}\geq 0$ for all $j$, we let the class of admissible perturbed matrix polynomials be

[TABLE]

The parameters $\omega_{j}\geq 0$, $j=0,\dots,m$, determine the maximum norm of the perturbation $\Delta_{j}$ of each matrix $A_{j}$, where $\|\cdot\|_{F}$ denotes the Frobenius norm. For instance, to keep $A_{j}$ unperturbed, we set $\omega_{j}=0$.

One approach to investigate the sensitivity of the spectrum of a matrix polynomial to admissible perturbations is to compute and plot the $\varepsilon$-pseudospectrum of $P$ for several $\varepsilon$-values, where the $\varepsilon$-pseudospectrum of $P(\lambda)$ for $\varepsilon>0$ is defined by

[TABLE]

The computation of a $\varepsilon$-pseudospectrum of a matrix polynomial generally is very computationally intensive, in fact, it is much more demanding than the computation of the $\varepsilon$-pseudospectrum of a single matrix; see Tisseur and Higham [22] for a discussion on several numerical methods including approaches based on using a transfer function, random perturbations, and projections to small-scale problems. The computations use the companion form of the matrix polynomial $P$. This requires working with matrices of order $mn$, whose generalized Schur factorization is computed. Therefore, the computational methods can be expensive to apply when $mn$ is fairly large and an approximation of the $\varepsilon$-pseudospectrum is determined on a mesh with many points. Details and counts of arithmetic floating point operations are provided in [22].

This paper describes a novel approach to approximate the $\varepsilon$-pseudospectra of $P$ by choosing particular rank-one perturbations of the matrices $A_{j}$ (or projected rank-one perturbations in case $A_{j}$ has a structure that is to be respected). The use of these rank-one perturbations yields approximations of the $\varepsilon$-pseudospectrum (1.3) for a lower computational cost than the computation of the $\varepsilon$-pseudospectrum. Our approach is inspired by Wilkinson’s analysis of eigenvalue perturbation of a single matrix; see [24]. It generalizes an approach recently developed in [18] for the efficient computation of structured or unstructured pseudospectra of a single matrix.

This paper is organized as follows. Section 2 reviews results on the sensitivity of a simple eigenvalue of a matrix polynomial, pseudospectra and the distance from defectivity for matrix polynomials is considered in Section 3, while the corresponding discussions for structured perturbations can be found in Sections 4 and 5. Algorithms for computing approximate structured and unstructured pseudospectra for matrix polynomials are described in Section 6, and a few computed examples are presented in Section 7. Finally, Section 8 contains concluding remarks.

2 The condition number of a simple eigenvalue of a matrix polynomial

Consider the matrix polynomial (1.1) and assume that the determinant of the leading coefficient matrix, $A_{m}$, is nonvanishing. Let $\lambda_{0}\in\mathbb{C}$ be an eigenvalue of $P$. Then the linear system of equations P(\lambda_{0})\mbox{\boldmath{x}}=\mbox{\bf 0} has a nonzero solution \mbox{\boldmath{x}}_{0}\in\mathbb{C}^{n} (a right eigenvector), and there is a nonzero vector \mbox{\boldmath{y}}_{0}\in\mathbb{C}^{n} such that \mbox{\boldmath{y}}_{0}^{H}P(\lambda_{0})=\mbox{\bf 0}^{H} (left eigenvector). Here the superscript H denotes transposition and complex conjugation. The algebraic multiplicity of $\lambda$ is its multiplicity as a zero of the scalar polynomial $\det(P(\lambda))$. The algebraic multiplicity is known to be larger than or equal to the geometric multiplicity of $\lambda_{0}$, which is the dimension of the null space of $P(\lambda_{0})$. The following result by Tisseur [21, Theorem 5] is important for the development of our numerical method. We therefore present a proof for completeness.

Proposition 2.1**.**

Let $\lambda\in\Lambda(P)$ be a simple eigenvalue, i.e. $\lambda\notin\Lambda(P^{\prime})$, with corresponding right and left eigenvectors $x$ and $y$ of unit Euclidean norm. Here $P^{\prime}$ denotes the derivative of $\lambda\rightarrow P(\lambda)$. Then the condition number of $\lambda$ is given by

[TABLE]

where $\omega(z)=\omega_{m}z^{m}+\ldots+\omega_{0}$. The maximal perturbations are

[TABLE]

for any unimodular $\eta\in\mathbb{C}$.

Proof.

Differentiating \sum_{j=0}^{m}(A_{j}+\epsilon\Delta_{j})\lambda^{j}(\varepsilon)\mbox{\boldmath{x}}(\varepsilon)=0 with respect to $\varepsilon\in\mathbb{C}$ yields

[TABLE]

Setting $\varepsilon=0$, one obtains

[TABLE]

where $\lambda=\lambda(0)$. It follows that

[TABLE]

Applying \mbox{\boldmath{y}}^{H} to both the right-hand side and left-hand side of this equality yields

[TABLE]

where we note that \mbox{\boldmath{y}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{x}}\neq 0 because $\lambda$ is a simple eigenvalue; see [2, Theorem 3.2]. Observing that \mbox{\boldmath{y}}^{H}P(\lambda)=\mbox{\bf 0}^{H}, and dividing by \mbox{\boldmath{y}}^{H}P^{\prime}(\lambda)\mbox{\boldmath{x}}, one has

[TABLE]

Taking absolute values yields

[TABLE]

where the inequality follows from the bounds $\|\Delta_{j}\|_{F}\leq\omega_{j}$, $j=0,\ldots,m$. Finally, letting the matrix $\Delta_{j}$ be a rank-one matrix of the form \eta\omega_{j}\mathrm{e}^{-\mathrm{i}j\arg(\lambda)}\mbox{\boldmath{y}}\mbox{\boldmath{x}}^{H} with unimodular $\eta\in\mathbb{C}$ (and therefore of Frobenius norm $\omega_{j}$) for all $j=0,\ldots,m$ shows the proposition. ∎

Remark 2.2*.*

Consider the standard eigenvalue problem with $m=1$, $A_{0}=A$, and $A_{1}=-I_{n}$. Then $P(\lambda)=A-\lambda I_{n}$ and $P^{\prime}(\lambda)=-I_{n}$. Setting $\omega_{0}=1$ and $\omega_{1}=0$, Proposition 2.1 yields the standard eigenvalue condition number \kappa(\lambda)=1/|\mbox{\boldmath{y}}^{H}\mbox{\boldmath{x}}|. When instead $A_{0}=A$ and $A_{1}=-B$, we obtain $P(\lambda)=A-\lambda B$ (and $P^{\prime}(\lambda)=-B$), and the proposition gives the generalized eigenvalue condition number \kappa(\lambda)=(\omega_{0}+\omega_{1}|\lambda|)/|\mbox{\boldmath{y}}^{H}B\mbox{\boldmath{x}}|; see [11].

Remark 2.3*.*

If $n=1$, the polynomial is scalar-valued. Let $\lambda$ be a simple root of $P$. Then the condition number of $\lambda$ is $\omega(|\lambda|)/|P^{\prime}(\lambda)|$.

3 The $\varepsilon$-pseudospectrum of a matrix polynomial and the distance from

defectivity

The $\varepsilon$-pseudospectrum of $P(\lambda)$ given by (1.3) is bounded if and only if $\det(A_{m}+\varepsilon\Delta_{m})\neq 0$ for all $\Delta_{m}$ such that $\|\Delta_{m}\|_{F}\leq\omega_{m}$. Therefore the boundedness of $\Lambda_{\varepsilon}(P)$ is guaranteed if $\varepsilon$ is such that the origin does not belong to the $\omega_{m}\varepsilon$-pseudospectrum of $A_{m}\in\mathbb{C}^{n\times n}$, which is given by

[TABLE]

It is easy to see that, if $\varepsilon$ satisfies the constraint

[TABLE]

then a first order analysis suggests that no component of $\Lambda_{\omega_{m}\varepsilon}(A_{m})$, which is approximately a disk of radius $\widehat{\kappa}(\lambda_{i}(A_{m}))\omega_{m}\varepsilon$ centered at $\lambda_{i}(A_{m})$ for $\omega_{m}\varepsilon$ small enough, can contain the origin. The origin is on the border of the disk centered at $\lambda_{i}(A_{m})$ when $|\lambda_{i}(A_{m})|={\widehat{\kappa}(\lambda_{i}(A_{m}))\omega_{m}\varepsilon}$. Here $\widehat{\kappa}(\lambda(M))$ denotes the traditional condition number of the eigenvalue $\lambda$ of the matrix $M\in\mathbb{C}^{n\times n}$.

Since by assumption $\det(A_{m})\neq 0$, the $\varepsilon$-pseudospectrum (1.3) has at most $mn$ bounded connected components. Any small connected component of the $\varepsilon$-pseudospectrum that contains exactly one simple eigenvalue $\lambda_{0}$ of the matrix polynomial $P$ is approximately a disk centered at $\lambda_{0}$ with radius $\kappa(\lambda_{0})\varepsilon$. A matrix polynomial $Q(\lambda)$ is said to be defective if it has an eigenvalue $\hat{\lambda}$, whose algebraic multiplicity is strictly larger than its geometric multiplicity; see [1]. Disjoint components of $\Lambda_{\varepsilon}(P)$ associated with distinct eigenvalues are, to a first order approximation, disjoint disks if $\varepsilon$ is strictly smaller than the distance $\varepsilon_{*}$ from defectivity of the matrix polynomial $P(\lambda)$, where

[TABLE]

A rough estimate of $\varepsilon_{*}$ is given by

[TABLE]

The disk centered at $\lambda_{i}$ is tangential to the disk centered at $\lambda_{j}$ when $|\lambda_{i}-\lambda_{j}|=(\kappa(\lambda_{i})+\kappa(\lambda_{j}))\,\varepsilon$. Let the index pair $\{\hat{\imath},\hat{\jmath}\}$ minimize the ratio (3.1) over all distinct eigenvalue pairs. We will refer to the eigenvalues $\lambda_{\hat{\imath}}$ and $\lambda_{\hat{\jmath}}$ as the most $\Lambda_{\varepsilon}$-sensitive pair of eigenvalues. We note that typically the most $\Lambda_{\varepsilon}$-sensitive pair of eigenvalues are not the eigenvalues with the largest condition numbers.

4 The structured condition number of a simple eigenvalue of a matrix polynomial

We briefly comment on structured eigenvalue condition numbers for a single matrix before turning to matrix polynomials. Consider the set ${{\mathcal{S}}}\,{\scriptscriptstyle{\begin{subarray}{c}\subset\\ \neq\end{subarray}}}\,\mathbb{C}^{n\times n}$ of structured matrices. For instance, the set may consist of symmetric, tridiagonal, or Toeplitz matrices. We are concerned with structured perturbations in ${{\mathcal{S}}}$. Let $M|_{{\mathcal{S}}}$ denote the matrix in ${\mathcal{S}}$ closest to $M\in\mathbb{C}^{n\times n}$ with respect to the Frobenius norm. This projection is used in definition of the eigenvalue condition number for structured perturbations, see [5, 16, 17, 18], where it is shown that the eigenvalue condition number for structured perturbations is smaller than the eigenvalue condition number for unstructured perturbations. We also will use the normalized projection

[TABLE]

in the definition of maximal structured perturbations in Proposition 4.1 below.

Matrix polynomials (1.1) are defined by $m+1$ matrices $A_{j}$, some or all of which may have a structure that is important for the application at hand. We refer to a matrix polynomial with at least one structured matrix $A_{j}$ as a structured matrix polynomial. To measure the sensitivity of the eigenvalues of a structured matrix polynomial to similarly structured perturbations, we proceed as follows. Let ${{\mathcal{S}}}_{j}$ be a set of structured matrices that the matrix $A_{j}$ of the matrix polynomial $P$ belongs to. If $A_{j}$ has no particular structure, then ${{\mathcal{S}}}_{j}={\mathbb{C}}^{n\times n}$. Introduce the set of sets of structured matrices ${{\mathcal{S}}}=\left\{{{\mathcal{S}}}_{0},{{\mathcal{S}}}_{1},\dots{\mathcal{S}}_{m}\right\}$ and let the class of admissible perturbed matrix polynomials be

[TABLE]

Proposition 4.1**.**

Let $\lambda\in\Lambda(P)$ be a simple eigenvalue with corresponding right and left eigenvectors $x$ and $y$ of unit Euclidean norm. Then the structured condition number of $\lambda$ is given by

[TABLE]

where

[TABLE]

The maximal perturbations are given by

[TABLE]

for any unimodular $\eta\in{\mathbb{C}}$.

Proof.

Differentiating \sum_{j=0}^{m}(A_{j}+\epsilon\Delta_{j})\lambda^{j}(\varepsilon)\mbox{\boldmath{x}}(\varepsilon)=\mbox{\bf 0} with respect to $\varepsilon$, as in the proof of Proposition 2.1, one obtains

[TABLE]

where $\Delta_{j}\in{{\mathcal{S}}}_{j}$ satisfies $\|\Delta_{j}\|_{F}\leq\omega_{j}$, $j=0,\ldots,m$. Substituting $\Delta_{j}$, for $j=0,\ldots,m$, by the structured matrix \eta\omega_{j}\mbox{\boldmath{y}}\mbox{\boldmath{x}}^{H}|_{\widehat{{\mathcal{S}}}_{j}}\in{{\mathcal{S}}}_{j} with Frobenius norm $\omega_{j}$, the upper bound \omega_{j}\,\|\mbox{\boldmath{y}}\mbox{\boldmath{x}}^{H}|_{{{\mathcal{S}}}_{j}}\|_{F} for |\mbox{\boldmath{y}}^{H}\Delta_{j}\mbox{\boldmath{x}}| is attained. Finally, letting \Delta_{j}=\eta\omega_{j}\mathrm{e}^{-\mathrm{i}j\arg(\lambda)}\mbox{\boldmath{y}}\mbox{\boldmath{x}}^{H}|_{\widehat{{\mathcal{S}}}_{j}} for all $j=0,\ldots,m$ gives

[TABLE]

This concludes the proof. ∎

Remark 4.2*.*

The structured condition number (4.1) is bounded above by the (unstructured) condition number (2.1). In fact, the former can be much smaller than the latter. For instance, let us consider the quadratic eigenvalue problem P(\lambda)\mbox{\boldmath{x}}=\mathbf{0}, with \mbox{\boldmath{x}}\neq\mathbf{0}, where

[TABLE]

with the same structured mass matrix $M$, damping matrix $C$ and stiffness matrix $K$ as in [22, Section 4.2], i.e., $M:=I_{n}$, $C:=10\,\mathrm{tridiag}(-1,3,-1)$, and $K:=5\,\mathrm{tridiag}(-1,3,-1)$. The $2n$ eigenvalues of the polynomial matrix are real and negative. In more detail, the spectrum is split into two sets: $n$ eigenvalues are spread approximately uniformly in the interval $[-50,-10]$ and $n$ eigenvalues are clustered at $-0.5$. Figure 1 shows the unstructured (i.e., standard) condition numbers (top graph) and structured condiition numbers (bottom graph) for each eigenvalue. The unstructured condition numbers are seen to be much larger than the structured condition numbers.

5 The structured $\varepsilon$-pseudospectrum of a matrix polynomial and the

structured distance from defectivity

The ${{\mathcal{S}}}$-structured $\varepsilon$-pseudospectrum of $P(\lambda)$ is for $\varepsilon>0$ defined by

[TABLE]

One has that $\Lambda_{\varepsilon}^{{\mathcal{S}}}(P)$ is bounded if and only if $\det(A_{m}+\varepsilon\Delta_{m})\neq 0$ for all $\Delta_{m}\in{{\mathcal{S}}}_{j}$ such that $\|\Delta_{m}\|_{F}\leq\omega_{m}$. Thus, the boundedness of $\Lambda_{\varepsilon}^{{\mathcal{S}}}(P)$ is guaranteed if $\varepsilon$ is such that $0\notin\Lambda_{\omega_{m}\varepsilon}^{{\mathcal{S}}_{m}}(A_{m})$, where $\Lambda_{\omega_{m}\varepsilon}^{{\mathcal{S}}_{m}}(A_{m})$ denotes the structured $\omega_{m}\varepsilon$-pseudospectrum of $A_{m}\in{{\mathcal{S}}_{m}}$, which is defined by

[TABLE]

We will assume that $\varepsilon$ satisfies the constraint

[TABLE]

Then a first order analysis suggests that no component of $\Lambda_{\omega_{m}\varepsilon}^{{\mathcal{S}}_{m}}(A_{m})$ contains the origin. In fact, when $\epsilon>0$ is small, the component that contains the eigenvalue $\lambda_{i}(A_{m})$ of $A_{m}$ is approximately a disk of radius $\widehat{\kappa}_{{\mathcal{S}}_{m}}(\lambda_{i}(A_{m}))\omega_{m}\varepsilon$ centered at $\lambda_{i}(A_{m})$. Here $\widehat{\kappa}_{{\mathcal{S}}_{m}}(\lambda)$ denotes the ${{\mathcal{S}}_{m}}$-structured condition number of the eigenvalue $\lambda$ in $\Lambda(M)$, where $M$ belongs to the set ${{\mathcal{S}}_{m}}$ of structured matrices in $\mathbb{C}^{n\times n}$.

Any small connected component of $\Lambda_{\varepsilon}^{{\mathcal{S}}}(P)$ that contains exactly one simple eigenvalue $\lambda_{0}\in\Lambda(P)$ is approximately a disk centered at $\lambda_{0}$ with radius $\kappa^{{\mathcal{S}}}(\lambda_{0})\varepsilon$. Such disks of $\Lambda_{\varepsilon}^{{\mathcal{S}}}(P)$ for distinct eigenvalues are, to a first order approximation, disjoint if $\varepsilon$ is strictly smaller than the structured distance $\varepsilon^{{\mathcal{S}}}_{*}$ from defectivity of the matrix polynomial $P(\lambda)$. This distance is given by

[TABLE]

A rough estimate of $\varepsilon_{*}^{{{\mathcal{S}}}}$ is provided by

[TABLE]

Similarly as in Section 3, the disk centered at $\lambda_{i}$ is tangential to the disk centered at $\lambda_{j}$ when $|\lambda_{i}-\lambda_{j}|=(\kappa^{{\mathcal{S}}}(\lambda_{i})+\kappa^{{\mathcal{S}}}(\lambda_{j}))\,\varepsilon$. Let the index pair $\{\hat{\imath},\hat{\jmath}\}$ minimize the ratio (5.2) over all distinct eigenvalue pairs. We will refer to the eigenvalues $\lambda_{\hat{\imath}}$ and $\lambda_{\hat{\jmath}}$ as the most $\Lambda_{\varepsilon}^{{\mathcal{S}}}$-sensitive pair of eigenvalues. We note that usually the most $\Lambda_{\varepsilon}^{{\mathcal{S}}}$-sensitive pair of eigenvalues is not made up of the worst conditioned eigenvalues with respect to structured perturbations.

6 Algorithms

This section describes algorithms based on Propositions 2.1 and 4.1 for computing approximations of unstructured and structured pseudospectra of matrix polynomials.

Let \{\lambda_{i},\mbox{\boldmath{x}}_{i},\mbox{\boldmath{y}}_{i}\}_{i=1}^{mn} denote eigen-triplets made up of the eigenvalues $\lambda_{i}$ and associated left and right unit eigenvectors, \mbox{\boldmath{x}}_{i} and \mbox{\boldmath{y}}_{i}, respectively, of the matrix polynomial $P$ defined by (1.1). We will assume the eigenvalues to be distinct. If a matrix polynomial has multiple eigenvalues, then we can apply the algorithms to the ones of algebraic multiplicity one. Throughout this section $\mathrm{i}=\sqrt{-1}$.

Algorithm 1 describes our numerical method for the approximation of the $\varepsilon$-pseudospectrum of a matrix polynomial $P$ defined by matrices $A_{j}$, $j=0,\ldots,m$, without particular structure. The algorithm first determines an estimate $\varepsilon$ of the distance to defectivity (3.1) of the matrix polynomial and the indices $\hat{\imath}$ and $\hat{\jmath}$ of the most $\Lambda_{\varepsilon}$-sensitive pair of eigenvalues of $P$. It then computes the rank-one matrices $\Delta_{\hat{\imath}}$ and $\Delta_{\hat{\jmath}}$ defined in Proposition 2.1 for all (simple) eigenvalues $\lambda$ of $P$ and for equidistant values on the unit circle in the complex plane. This defines the perturbations of the polynomial $P$ at the eigenvalues $\lambda$. The spectra of the perturbations of $P$ so obtained are displayed. This simple approach typically provides valuable insight into properties of the $\varepsilon$-pseudospectrum of $P$.

Algorithm 2 is an analogue of Algorithm 1 for the approximation of the structured $\varepsilon$-pseudospectrum of a matrix polynomial. The algorithm differs from Algorithm 1 in that the distance to defectivity in the latter algorithm is replaced by the structured distance to defectivity (5.2) and the rank-one perturbations are replaced by structured rank-one perturbations defined in Proposition 4.1.

Both Algorithms 1 and 2 are easy to implement. The algorithms require the computation of the $mn$ eigenvalues of $n\times n$ polynomial matrices. Evaluating the spectrum of $2N$ perturbed polynomial matrices is the main computational burden and easily can be implemented efficiently on a parallel computer. However, a laptop computer was sufficient for the computed examples reported in the following section.

7 Numerical examples

The computations were performed on a MacBook Air laptop computer with a 1.8Ghz CPU and 4 Gbytes of RAM. All computations were carried out in MATLAB with about $16$ significant decimal digits.

Example 1. Consider the matrix polynomial $P(\lambda)=A_{2}\lambda^{2}+A_{1}\lambda+A_{0}$, where $A_{0}$ and $A_{1}$ are real $5\times 5$ matrices with normally distributed random entries with zero mean and variance, and $A_{2}$ is a real tridiagonal Toeplitz matrix of the same order with similarly distributed random diagonal, superdiagonal, and subdiagonal entries. We choose the weights $\omega_{i}=\|A_{i}\|_{F}$, $i=0:2$. The eigenvalues of $P$ and their standard and structured condition numbers are shown in Table 1. The structured condition numbers can be seen to be smaller than the standard condition numbers.

The estimate (3.1) of the (unstructured) distance from defectivity $\varepsilon_{*}$ is $\varepsilon_{1}=0.0127$. It is achieved for the indices $5$ and $7$, as well as for the indices $4$ and $6$, of the most $\Lambda_{\varepsilon}$-sensitive pairs of eigenvalues. The left plot in Figure 2 displays the spectrum of matrix polynomials of the form $P(\lambda)+\varepsilon_{1}\mathrm{e}^{\mathrm{i}\theta_{k}}W_{5}(\lambda)$ and $P(\lambda)+\varepsilon_{1}\mathrm{e}^{\mathrm{i}\theta_{k}}W_{7}(\lambda)$ for $\theta_{k}=2\pi(k-1)/N$, $k=1:N$, and $N=5\cdot 10^{2}$. Thus, the spectrum of $10^{3}$ matrix polynomials are determined. Details of the computations are described by Algorithm 1. We recall that the “curves” surrounding the eigenvalues $\lambda_{j}$ lie inside the $\varepsilon_{1}$-pseudospectrum of $P$. The figure illustrates that the eigenvalues $\lambda_{5}$ and $\lambda_{7}$ might coalesce already for a small perturbation of $P$.

We remark that since the matrices $A_{i}$, $i=0:2$, that define the matrix polynomial $P$ are real, the eigenvalues of $P$ appear in complex conjugate pairs. The pseudospectrum of matrix polynomials determined by real matrices is known to be symmetric with respect to the real axis in the complex plane. The fact that the left plot of Figure 2 is not symmetric with respect to the imaginary axis depends on that it only shows the spectra of the matrix polynomials $P(\lambda)+\varepsilon_{1}\mathrm{e}^{\mathrm{i}\theta_{k}}W_{5}(\lambda)$ and $P(\lambda)+\varepsilon_{1}\mathrm{e}^{\mathrm{i}\theta_{k}}W_{7}(\lambda)$ associated with the eigenvalues $\lambda_{5}$ and $\lambda_{7}$ of $P$, but not of the polynomials $P(\lambda)+\varepsilon_{1}\mathrm{e}^{\mathrm{i}\theta_{k}}W_{4}(\lambda)$ and $P(\lambda)+\varepsilon_{1}\mathrm{e}^{\mathrm{i}\theta_{k}}W_{6}(\lambda)$ associated with the eigenvalues $\lambda_{4}$ and $\lambda_{6}$. A plot of eigenvalues of all these polynomials is symmetric with respect to the real axis in the complex plane.

We compare the approximation of the $\varepsilon_{1}$-pseudospectrum shown in the left plot of Figure 2 with an approximation of the $\varepsilon_{1}$-pseudospectrum obtained by perturbing $P$ by random rank-one matrices. Specifically, the right plot of Figure 2 shows the spectrum of matrix polynomials of the form $P(\lambda)+\varepsilon_{1}\mathrm{e}^{\mathrm{i}\theta_{k}}E(\lambda)$ with $\theta_{k}=2\pi(k-1)/N$, $k=1:N$, where $N=10^{6}$, and $E(\lambda)=\sum_{h=0}^{m}\omega_{h}\lambda^{h}R_{h}$. Here $R_{h}$ is a rank-one random matrix scaled to have unit Frobenius norm. Despite using $10^{6}$ perturbations of $P$, which are many more perturbations than used for producing the left plot, the right plot of Figure 2 does not indicate that any eigenvalue of $P$ might coalesce under small perturbations of the matrix polynomial. This important property clearly is difficult to detect by using random rank-one perturbations.

Next we turn to structured pseudospectra and perturbations. We obtain from (5.2) the estimate $\varepsilon_{2}=0.0266$ of the structured distance from defectivity $\varepsilon^{{\mathcal{S}}}_{*}$. It is achieved for the eigenvalues $\lambda_{8}$ and $\lambda_{9}$. The left plot in Figure 3 displays the spectra of matrix polynomials of the form $P(\lambda)+\varepsilon_{2}\mathrm{e}^{\mathrm{i}\theta_{k}}W^{{\mathcal{S}}}_{8}(\lambda)$ and $P(\lambda)+\varepsilon_{2}\mathrm{e}^{\mathrm{i}\theta_{k}}W^{{\mathcal{S}}}_{9}(\lambda)$ with $\theta_{k}=2\pi(k-1)/N$, $k=1:N$, for $N=5\cdot 10^{2}$. The computations are described by Algorithm 2. The plot shows that the eigenvalues $\lambda_{8}$ and $\lambda_{9}$ might coalesce under small perturbations of $P$.

The right plot of Figure 2 shows the spectrum of matrix polynomials of the form $P(\lambda)+\varepsilon_{2}\mathrm{e}^{\mathrm{i}\theta_{k}}E^{{\mathcal{S}}}(\lambda)$ with $\theta_{k}=2\pi(k-1)/N$, $k=1:N$, where $N=10^{6}$, $E^{{\mathcal{S}}}(\lambda)=\sum_{h=0}^{m}\omega_{h}\lambda^{h}R^{{\mathcal{S}}}_{h}$, and $R^{{\mathcal{S}}}_{h}:=R_{h}|_{\widehat{{\mathcal{S}}}_{h}}$ is a unit-norm rank-one random matrix projected into ${{\mathcal{S}}}_{h}$. Despite using $10^{6}$ perturbations, the plot does not indicate that any eigenvalues of $P$ might coalesce under small structured perturbations. $\blacksquare$

Example 2. Consider the matrix polynomial $P(\lambda)=A_{2}\lambda^{2}+A_{1}\lambda+A_{0}$ defined by

[TABLE]

[TABLE]

This polynomial is discussed in [22, Section 4.1]. We choose $\omega=\{1,1,1\}$ similarly as in [22]. The eigenvalues and their condition numbers are shown in Table 2.

Figure 4 displays an approximation of the $\varepsilon$-pseudospectrum of $P(\lambda)$ obtained by letting $\varepsilon=10^{-0.8}$ (like in [22]) and computing the eigenvalues of matrix polynomials of the form $P(\lambda)+\varepsilon\mathrm{e}^{\mathrm{i}\theta_{k}}W_{1}(\lambda)$ for $\theta_{k}=2\pi(k-1)/10^{2}$, $k=1:10^{2}$, where $W_{1}(\lambda)$ is a maximal perturbation associated with the eigenvalue $\lambda_{1}$ (marked by red square) with the largest condition number. Details of the computations are described by Algorithm 1. $\blacksquare$

Example 3. We consider the matrix polynomial $P(\lambda)=M\lambda^{2}+C\lambda+K$ with the structure $\cal{S}$ defined in Remark 4.2. This polynomial is considered in [22, Section 4.2]. We choose the weights $\omega=\{\|K\|_{F},\|C\|_{F},\|M\|_{F}\}$ and obtain from (5.2) the estimate $\varepsilon_{2}=3.5709\cdot 10^{-7}$ of the structured distance from defectivity $\varepsilon^{{\mathcal{S}}}_{*}$. It is achieved for the eigenvalues $\lambda_{493}$ and $\lambda_{494}$. These eigenvalues are the most $\Lambda_{\varepsilon_{2}}^{{\mathcal{S}}}$-sensitive pair, but they are not the most ill-conditioned eigenvalues, despite that their relative distance is only $10^{-6}$.

Figure 5 displays the spectra of matrix polynomials of the form $P(\lambda)+\varepsilon_{2}\mathrm{e}^{\mathrm{i}\theta_{k}}W^{{\mathcal{S}}}_{493}(\lambda)$ and $P(\lambda)+\varepsilon_{2}\mathrm{e}^{\mathrm{i}\theta_{k}}W^{{\mathcal{S}}}_{494}(\lambda)$ with $\theta_{k}=2\pi(k-1)/10^{2}$, $k=1:10^{2}$. The computations are described by Algorithm 2. $\blacksquare$

8 Conclusions

This paper describes a novel and fairly inexpensive approach to determine the sensitivity of eigenvalues of a matrix polynomial. Eigenvalues of perturbed matrix polynomials are computed, where the perturbations are chosen to shed light on whether eigenvalues of the given matrix polynomial may coalesce under small perturbations.