New numerical algorithm for deflation of infinite and zero eigenvalues and full solution of quadratic eigenvalue problems

TL;DR

This paper introduces an improved numerical algorithm for accurately computing all eigenvalues and eigenvectors of quadratic matrix pencils, especially handling zero and infinite eigenvalues more robustly and efficiently.

Contribution

It enhances the existing quadeig algorithm by improving stability, robustness, and the ability to deflate additional infinite eigenvalues using a new canonical form approach.

Findings

Enhanced backward stability and numerical robustness.

Ability to deflate more infinite eigenvalues.

Superior performance demonstrated through empirical testing.

Abstract

This paper presents a new method for computing all eigenvalues and eigenvectors of quadratic matrix pencil. It is an upgrade of the quadeig algorithm by Hammarling, Munro and Tisseur, which attempts to reveal and remove by deflation certain number of zero and infinite eigenvalues before QZ iterations. Proposed modifications of the quadeig framework are designed to enhance backward stability and to make the process of deflating infinite and zero eigenvalues more numerically robust. In particular, careful preprocessing allows scaling invariant/component-wise backward error and thus better condition number. Further, using an upper triangular version of the Kronecker canonical form enables deflating additional infinite eigenvalues, in addition to those inferred from the rank of leading coefficient matrix. Theoretical analysis and empirical evidence from thorough testing of the software implementation confirm superior numerical performances of the proposed method.

Equations

$$Q(\lambda )x=({\lambda }^{2}M+\lambda C+K)x=0,M,C,K\in {\mathbb{C}}^{n\times n}.$$

$$C_2(\lambda )=\left(\begin{array}{cc}C& -{\mathbb{I}}_n\\ K& 0_n\end{array}\right)-\lambda \left(\begin{array}{cc}-M& 0_n\\ 0& -{\mathbb{I}}_n\end{array}\right)\equiv A-\lambda B.$$

$$\eta (\lambda \ne \infty ,x)=\frac{\Vert ({\lambda }^{2}M+\lambda C+K)x{\Vert }_2}{(|\lambda {|}^2\Vert M{\Vert }_2+|\lambda |\Vert C{\Vert }_2+\Vert K{\Vert }_2)\Vert x{\Vert }_2};\eta (\lambda =\infty ,x)=\frac{\Vert Mx{\Vert }_2}{\Vert M{\Vert }_2\Vert x{\Vert }_2},$$

$$\eta (\lambda ,x)=\min \{ϵ:({\lambda }^2\tilde{M}+\lambda \tilde{C}+\tilde{K})x=0,\Vert \delta M{\Vert }_2\le ϵ\Vert M{\Vert }_2,\Vert \delta C{\Vert }_2\le ϵ\Vert C{\Vert }_2,\Vert \delta K{\Vert }_2\le ϵ\Vert K{\Vert }_2\}.$$

$$M=\left(\begin{array}{cc}{M}_0& 0_{3\times 2}\\ 0_{2\times 3}& 0_{2\times 2}\end{array}\right),C=\left(\begin{array}{cc}{C}_0& 0_{3\times 2}\\ 0_{2\times 3}& 0_{2\times 2}\end{array}\right),K=\left(\begin{array}{cc}{K}_0& -F_0^T\\ F_0& 0_{2\times 2}\end{array}\right),$$

$$C_0=\left(\begin{array}{ccc}-1.52143& -1.55168& 1.55168\\ 3.22064& 3.28467& -3.28467\\ -3.22064& -3.28467& 3.28467\end{array}\right),K_0=\left(\begin{array}{ccc}67.4894& 69.2393& -69.2393\\ 69.8124& 1.68624& -1.68617\\ -69.8123& -1.68617& -68.2707\end{array}\right).$$

$$\begin{array}[]{ll}\lambda_{1}=\texttt{-5.1616e-002 -2.2435e-001i}&\lambda_{6}=\texttt{1.9573e+06 - 3.3993e+06i}\cr\lambda_{2}=\texttt{-5.1616e-002 +2.2435e-001i}&\lambda_{7}=\texttt{1.9573e+06 + 3.3993e+06i}\cr\lambda_{3}=\texttt{2.7700e+05 - 4.7988e+05i}&\lambda_{8}=\texttt{-3.9305e+06 + 0.0000e+00i}\cr\lambda_{4}=\texttt{2.7700e+05 + 4.7988e+05i}&\lambda_{9}=\texttt{Inf}\cr\lambda_{5}=\texttt{-5.5417e+05 + 0.0000e+00i }&\lambda_{10}=\texttt{Inf}\end{array}$$

$$\begin{array}[]{l}\lambda_{1}=\texttt{-5.161621336216381e-002 -2.243476109085836e-001i}\cr\lambda_{2}=\texttt{-5.161621336216381e-002 +2.243476109085836e-001i}\cr\lambda_{3}=\cdots=\lambda_{8}=\texttt{Inf},\;\;\;\lambda_{9}=\lambda_{10}=\texttt{-Inf}.\end{array}$$

$$\omega (\lambda ,x)$$

$$AP=Q\left(\begin{array}{c}R\\ 0\end{array}\right),\text{mbox}where|R_{ii}{|}^2\ge \sum _{k=i}^j|R_{kj}{|}^2,\text{mbox}forall1\le i\le j\le n.$$

$$R=\left(\begin{array}{cc}{R}_{[11]}& R_{[12]}\\ 0& R_{[22]}\end{array}\right),R_{[11]}\in {\mathbb{C}}^{k\times k},$$

$$\text{mbox}then(A+\Delta A)P\equiv (A-Q\left(\begin{array}{cc}0& 0\\ 0& R_{[22]}\\ 0& 0\end{array}\right)P^T)P=Q\left(\begin{array}{cc}{R}_{[11]}& R_{[12]}\\ 0& 0\\ 0& 0\end{array}\right),\Delta A\equiv -Q\left(\begin{array}{cc}0& 0\\ 0& R_{[22]}\\ 0& 0\end{array}\right)P^T.$$

$$\Vert \Delta A{\Vert }_F=\Vert R_{[22]}{\Vert }_F\le \sqrt{n-k}|R_{k+1,k+1}|.$$

$$\sqrt{n-k}|R_{k+1,k+1}|/\Vert A{\Vert }_F\le \gamma ,$$

$$(A+\delta A)\tilde{P}=\hat{Q}\left(\begin{array}{c}\tilde{R}\\ 0\end{array}\right),\Vert \delta A{\Vert }_F\le {ϵ}_{qr}\Vert A{\Vert }_F,\Vert \tilde{Q}-\hat{Q}{\Vert }_F\le {ϵ}_{qr}.$$

$$\Vert \delta A(:,i){\Vert }_2\le {ϵ}_{qr}\Vert A(:,i){\Vert }_2,i=1,\dots ,n.$$

$$(A+\delta A+\Delta A)\tilde{P}=\hat{Q}\left(\begin{array}{cc}{\tilde{R}}_{[11]}& {\tilde{R}}_{[12]}\\ 0& 0\\ 0& 0\end{array}\right),\Vert \delta A{\Vert }_F\le {ϵ}_{qr}\Vert A{\Vert }_F,\Vert \tilde{Q}-\hat{Q}{\Vert }_F\le {ϵ}_{qr},$$

$$\Delta A\equiv -\hat{Q}\left(\begin{array}{cc}0& 0\\ 0& {\tilde{R}}_{[22]}\\ 0& 0\end{array}\right){\tilde{P}}^T\text{mbox}(estimatedusing(\text{ref}eq:DA-R22)),(\Delta A\tilde{P})(:,1:k)=0.$$

$$\Vert A(1,:){\Vert }_{\infty }\ge \Vert A(2,:){\Vert }_{\infty }\ge \dots \ge \Vert A(m,:){\Vert }_{\infty },$$

$$\underset{i=1:m}{\max }\frac{\Vert \delta A(i,:){\Vert }_{\infty }}{\Vert A(i,:){\Vert }_{\infty }}\le {ϵ}_{qr}\underset{i=1:m}{\max }\frac{{\alpha }_i}{\Vert A(i,:){\Vert }_{\infty }}\equiv {ϵ}_{qr}^{\rightrightarrows },\text{mbox}where{\alpha }_i=\underset{j,k}{\max }|{\tilde{A}}_{ij}^{(k)}|,$$

$$\alpha_{i}\leq\left\{\begin{array}[]{ll}\sqrt{m-i+1}(1+\sqrt{2})^{i-1}\|A(i,:)\|_{\infty},&i\leq n\cr(1+\sqrt{2})^{n-1}\|A(i,:)\|_{\infty},&i>n\end{array}\right.,$$

$$\left(\begin{array}{c}{R}_{[11]}^{\ast }\\ R_{[12]}^{\ast }\end{array}\right)=Z_R^{\ast }\left(\begin{array}{c}{T}_{[11]}^{\ast }\\ 0\end{array}\right),$$

$$A=Q\left(\begin{smallmatrix}T_{[11]}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}\end{smallmatrix}\right)Z,\;\;\mbox{where $Z=Z_{R}P^{T}$.}$$

$$A+\delta A+\Delta A+\hat{Q}\left(\begin{array}{c}\delta \tilde{R}\\ 0\end{array}\right){\tilde{P}}^T=\hat{Q}\left(\begin{array}{cc}{\Pi }_1& 0\\ 0& {\mathbb{I}}_{m-k}\end{array}\right)\left(\begin{array}{cc}{\tilde{T}}_{[11]}& 0_{k,n-k}\\ 0_{m-k,k}& 0\end{array}\right){\hat{Z}}^{\ast }{\Pi }_2{\tilde{P}}^T,$$

$$(\tilde{R}+\delta \tilde{R}{)}^{\ast }{\Pi }_1={\Pi }_2^T\hat{Z}\left(\begin{array}{c}{\tilde{T}}_{[11]}^{\ast }\\ 0\end{array}\right),$$

$$(A+\delta A+\Delta A+\hat{Q}\left(\begin{array}{c}\delta \tilde{R}\\ 0\end{array}\right){\tilde{P}}^T)=\hat{Q}\left(\begin{array}{cc}{\Pi }_1& 0\\ 0& {\mathbb{I}}_{m-k}\end{array}\right)\left(\begin{array}{cc}{\tilde{T}}_{[11]}& 0\\ 0& 0\end{array}\right){\hat{Z}}^{\ast }{\Pi }_2{\tilde{P}}^T,$$

$$(P_{r,M}M){\Pi }_M=Q_M{R}_M,R_M=\left(\begin{array}{cccc}\ast & \ast & \ast & \ast \\ 0& \ast & \ast & \ast \\ 0& 0& \ast & \ast \\ 0& 0& 0& 0\end{array}\right)=\left(\begin{array}{c}{\hat{R}}_M\\ 0_{n-r_M,n}\end{array}\right),$$

$$\Vert \delta {\tilde{Q}}_M{\Vert }_F\equiv \Vert {\hat{Q}}_M-{\tilde{Q}}_M{\Vert }_F\le ϵ,\text{mbox}and\Vert \delta M(:,i){\Vert }_2\le {ϵ}_{qr}\Vert M(:,i){\Vert }_2,i=1,\dots ,n.$$

$$(M+\delta M){\tilde{\Pi }}_M={\hat{Q}}_M{\tilde{R}}_M$$

$$\widetilde{R}_{M}=\begin{pmatrix}(\widetilde{R}_{M})_{[11]}&(\widetilde{R}_{M})_{[12]}\cr\mathbf{0}_{n-k,k}&(\widetilde{R}_{M})_{[22]}\end{pmatrix},\;\;\mbox{where the $(n-k)\times(n-k)$ block $(\widetilde{R}_{M})_{[22]}$ "is small"}.$$

$$(M+\delta M+\underset{\Delta M}{\underbrace{{\hat{Q}}_M\left(\begin{array}{cc}0& 0\\ 0& -({\tilde{R}}_M{)}_{[22]}\end{array}\right){\tilde{\Pi }}_M^T}}){\tilde{\Pi }}_M={\hat{Q}}_M\left(\begin{array}{cc}({\tilde{R}}_M{)}_{[11]}& ({\tilde{R}}_M{)}_{[12]}\\ 0_{n-k,k}& 0_{n-k,n-k}\end{array}\right).$$

Full text

1 Introduction

We propose a new numerical method for computing all eigenvalues with the corresponding left and right eigenvectors of the $n\times n$ quadratic eigenvalue problem

[TABLE]

Such problem is usually associated with a damped system $M\ddot{q}(t)+C\dot{q}(t)+Kq(t)=f(t)$, which arises in many applications in science and engineering. For an excellent review of applications and numerical methods we refer to [33].

The dimension $n$ can be large, e.g. in millions [26], and the matrices usually possess some structure such as e.g. sparsity that depends on the geometry of the underlying physical model, discretization method etc. In such cases, one is often interested in finding only those eigenvalues (with the corresponding eigenvectors) that meet some particular specifications; examples include the tasks of finding $k\ll n$ eigenvalues closest to some target value $\xi\in\mathbb{C}$, or all eigenvalues belonging to some prescribed $\Omega\subset\mathbb{C}$. In some other applications, the dimension $n$ may be moderate with non-sparse coefficient matrices and all eigenvalues, optionally with the corresponding eigenvectors, are requested.

Even in the large scale problems with sparse (or otherwise structured) coefficient matrices, an iterative method, such as e.g. the Second Order Arnoldi (SOAR [2]) uses the Rayleigh-Ritz projection and the full solution of the projected problem $(\lambda^{2}[X^{*}MX]+\lambda[X^{*}CX]+[X^{*}KX])v=\mathbf{0}$ of moderate size is needed. Here $X$ is an $n\times k$ orthonormal matrix that should span a subspace rich in the desired spectral information, and the dimension $k$ is moderate, $k\ll n$, so that the full solution with dense $k\times k$ coefficient matrices is feasible.

We focus to the complete solution of (1.1), i.e. finding all eigenvalues and eigenvectors; the dimension $n$ is assumed moderate. This is usually solved by linearization – transformation to an equivalent eigenvalue problem for linear pencils, such as e.g.

[TABLE]

From the space of linearizations, one is advised to use the one that preserves, as much as possible, structural (spectral) properties of the original problem. This is important because numerical methods can leverage particular structure when applied to an appropriate linearization. It should be stressed that linearizations algebraically equivalent to (1.1) are not numerically equivalent to the quadratic problem in the sense of sensitivity when solved in finite precision arithmetic.

For the linearized problem (1.2), the method of choice is the QZ method. In some cases this may run into difficulties, in particular if $M$ is exactly or nearly rank deficient, which leads to (numerically) infinite eigenvalues. Whether the QZ is troubled by the presence of the infinite eigenvalues is a tricky question. On one hand, from an expert point of view, [37] shows that the presence of infinite eigenvalues cannot derange the transmission of shifts during the QZ iterations. On the other hand, the filigran-work with plane rotations [37] that ensures this robustness is often replaced (in software implementations) with faster block-oriented procedures, and infinite eigenvalues may remain latent and come as spurious absolutely large values, with negative impact on the accuracy of the computed genuine finite eigenvalues. In the trade-off between accuracy and speed, the developers usually choose speed and flops seem more important than the number of correct digits in the output. However, the block oriented computation does provide considerably better performances, and this is an important aspect of software development.

So, in general, it is advantageous to deflate the infinite eigenvalues early in the computational scheme, using computational tools that can be efficiently implemented on the state of the art computers. Similarly, if $K$ is rank deficient, then its null space provides eigenvectors for the eigenvalue $\lambda=0$, and removing it in a preprocessing phase facilitates more efficient computation of the remaining eigenvalues. In both cases a nontrivial decision about the numerical rank has to be made.

This motivated [14] to develop a new deflation scheme that removes $n-r_{M}$ infinite and $n-r_{K}$ zero eigenvalues, where $r_{M}=\mathop{\rm rank}\nolimits{M}$ and $r_{K}=\mathop{\rm rank}\nolimits{K}$. The remaining generalized linear eigenvalue problem is of the dimension $r_{M}+r_{K}$; it may still have some infinite and zero eigenvalues, and their detection then depends on the performance of the QZ algorithm. The computation is done in the framework of the linearization (1.2), and the proposed approach has resulted in the quadeig algorithm.

For a pair $(\lambda,x)$ of an eigenvalue and its right eigenvector, the optimal backward error to test whether $(\lambda,x)$ is acceptable is defined as (see [34, §2.2])

[TABLE]

where, by definition, $\lambda=\infty$ is an eigenvalue of (1.1) if and only if $\mu=0$ is an eigenvalue of the reversed pencil $\mu^{2}K+\mu C+M$. The backward stability (or the residual) are usually measured by $\eta(\lambda,x)$ and the computation is deemed stable/acceptable if $\eta(\lambda,x)$ is, up to a moderate factor, of the order of the machine precision. This is justified by the fact that we have an exact eigenpair of a nearby pencil with coefficient matrices $\widetilde{M}=M+\delta M$, $\widetilde{C}=C+\delta C$, $\widetilde{K}=K+\delta K$, where

[TABLE]

Using this, essentially only necessary, condition is often used as sufficient for declaring the computed result reliable. The theoretical advances to identify condition numbers as measures of sensitivity of the eigenvalues [34, §2.3] are rarely used as part of numerical software solution. Unfortunately, such an approach is potentially prone to errors. The best way to illustrate the difficulty is by an example.

Example 1.1**.**

Consider the model of a two-dimensional three-link mobile manipulator (robot for e.g. facade cleansing, construction) [18].111This example is included in the NLEVP benchmark collection [3]. Linearisation around a particularly chosen working point yields a $5\times 5$ quadratic matrix polynomial (1.1) with the coefficient matrices

[TABLE]

where $F_{0}=\left(\begin{smallmatrix}1&0&0\\ 0&0&1\end{smallmatrix}\right)$ enforces certain constraints, $M_{0}=\left(\begin{smallmatrix}18.7532&-7.94493&7.94494\\ -7.94493&31.8182&-26.8182\\ 7.94494&-26.8182&26.8182\\ \end{smallmatrix}\right)$,

[TABLE]

The matrix $K$ has full rank, and the rank of $M$ is $r_{M}=3$. This means that there are at least $n-r_{M}=2$ infinite eigenvalues. We computed the eigenvalues using the quadeig algorithm [14] with default options as follows:222Al computations were done using Matlab 8.5.0.197613 (R2015a) under Windows 10. The roundoff unit in Matlab is $\texttt{eps}\approx 2.2\cdot^{-16}$.

[TABLE]

Next, we compute the eigenvalues by a direct application of the QZ algorithm to (1.2). The QZ algorithm found $8$ infinite and two finite eigenvalues:

[TABLE]

*If we check the backward errors, then both (1.5) and (1.6) (with the corresponding right eigenvectors) pass the backward stability test with $\eta(\lambda,x)$ at the level of machine roundoff for all eigenvalues. Yet, (1.5) has six finite eigenvalues more than (1.6), and those six values are $O(10^{6})$ in modulus. How can we tell whether those values correspond to genuine $O(10^{6})$ big finite eigenvalues, or whether they should be declared spurious? After all, the backward error (1.3) below $n\cdot\texttt{eps}$ seems reassuring, but the same holds for (1.6) as well. In a particular application, engineer’s intuition may help, but for a black-boxed numerical software this remains an important and difficult problem. This clearly shows that the problem is nontrivial; one can easily imagine similar ambiguity when solving a problem with, say, $n>10^{4}$, where different methods give completely different results, and all with small backward error (1.3). $\diamond$

An important message of Example 1.1 is that backward stability measured by (1.3) is only a necessary/desirable condition which may not be enough for reliable computation. More strengthen measure for backward stability is required, such as the component-wise backward error (residual)

[TABLE]

and the sensitivity of the eigenvalues must be estimated using appropriate condition number, see [15], [16], [17]. Using $\omega(\lambda,x)$ in Example 1.1 provides valuable information. Although quadeig is considered better that polyeig [14], in this example the QZ algorithm (including its preprocessing) seems more successful, see Figure 1.

Entry–wise backward stability (in the sense that the relative backward error in each matrix entry of $M$, $C$, $K$ is small) is in general not feasible; not even in the case $n=1$, see [25]. For practical usage, a mixed analysis (backward-forward) might be more appropriate, with more general reference values for measuring size of the entry-wise/column-wise backward errors, e.g. $|\delta M|\leq\epsilon\Xi_{M}$, where $\Xi_{M}\geq\mathbf{0}$.

To the best of our knowledge, quadeig [14] is the only published algorithm for the full solution of (1.1) that is supported by a backward error analysis and thoroughly tested, with publicly available software implementation. Our main contribution in this work is an upgrade of of the framework set by quadeig, and it has two main components. First, we treat the numerical rank issue and te backward error more carefully, so that our proposed method has better structured backward error in the sense of [15]. This is confirmed by theoretical analysis and numerical experiments. Second, we interpret the deflation strategy deployed in quadeig as the beginning step of computing the Kronecker canonical structure associated with the zero/infinite eigenvalue, which motivates taking further steps toward the Kronecker canonical form (KCF), thus opening possibility for removing additional infinite eigenvalues. Indeed, we show that our proposed method is capable of discovering infinite eigenvalues that are beyond the reach of quadeig. A preview of the performance of the new proposed method, designated as KVADeig333Pronounced as quadeig., is given in Figure 1.

1.1 Outline of the paper

The rest of the paper is organized as follows. In §2 we review, for the readers convenience, numerical computation of rank revealing decompositions. In particular, we review the pivoted QR factorization and the complete orthogonal (URV) decomposition, and in §2.2 we emphasize the importance of pivoting for more structured backward error bounds in case of unbalanced rows in data matrices. Further, in §3, we recall the Van Dooren’s algorithm [9] for computing the KCF of general matrix pencil $A-\lambda B$. We rewrite the Van Dooren’s algorithm in an upper triangular form, and use the rank revealing QR factorization instead of the SVD; this procedure is at the core of our new algorithm. Numerical equivalence of linearization and the original quadratic problem in the sense of backward stability is a delicate issue and special preprocessing by parameter scaling is recommended; see an enlightening discussion and examples in [17], and implementation details in quadeig [14]. Scaling is briefly reviewed in §4; in §4.2 we also propose balancing matrix entries by diagonal scalings of the coefficients of (1.1). In §5 we describe our modified approach to the quadeig type preprocessing. The goal is to adapt the rank revealing steps to the form of the KCF described in §3.2. In §5.3 we provide detailed backward error analysis of the first two steps of the reduction to the KCF, and prove backward stability in terms of the original data matrices. The proof is instructive because it shows how the modifications introduced in KVADeig contribute to backward stability. The global structure of the new algorithm is summarized in §6. Numerical experiments are presented in §7. Concluding remarks and comments on ongoing and future work are given in §8.

2 Rank revealing decompositions

Since detecting zero or infinite eigenvalues is based on numerical rank decision, we briefly discuss rank revealing decompositions (RRD, see [8]). For a general $m\times n$ matrix $A$, we say that $A=XDY^{*}$ is a rank revealing decomposition if both $X$ and $Y$ are of full column rank and well conditioned, and $D$ is diagonal nonsingular. In practical computation, such a decomposition is computed only approximately and we have $A+\delta A=\widetilde{X}\widetilde{D}\widetilde{Y}^{*}$, where $\delta A$ denotes initial uncertainty and/or the backward error that corresponds to the numerically computed $\widetilde{X}$, $\widetilde{D}$ and $\widetilde{Y}$. Hence, any decision on the rank actually applies to $A+\delta A$. We refer to [13] for a more in depth discussion and definition of a numerical rank.

Since the full rank matrices are open dense set in $\mathbb{C}^{m\times n}$ ($\mathbb{R}^{m\times n}$), it is unlikely that, in general, the rank deficiency will be determined correctly in finite precision computation. Furthermore, in many applications the matrix has been already contaminated by errors (initial uncertainty, previous computational steps etc.) and a firm statement about its rank is illusory. Since the zero and the infinite eigenvalues are rooted in rank deficiency, it is clear that revealing them numerically is a delicate issue.

In this section, we review the pivoted QR factorization, and the complete orthogonal decomposition (URV, ULV), which will be used as rank revealing decompositions. We discuss the backward stability, which will be an important ingredient in the analysis of the proposed algorithm.

2.1 QR factorization with column and complete pivoting

QR factorization with column pivoting is a tool of trade in a plethora of applications, in particular when the numerical rank of a matrix plays an important role. It is the main computational tool in the reduction procedure in the quadeig algorithm as well.

Particularly successful is the Businger–Golub [6] pivot strategy which, for $A\in\mathbb{C}^{m\times n}$, computes a permutation matrix $P$, a unitary $Q$ and an $\min(m,n)$ upper triangular (trapezoidal if $m<n$) matrix $R$ such that

[TABLE]

Here, for the sake of brevity, we consider only the case $m\geq n$. If $m<n$, then $R$ is $m\times n$ upper trapezoidal and the zero block in (2.1) is void. If $\rho=\mathrm{rank}(A)$, then $\prod_{i=1}^{\rho}R_{ii}\neq 0$ and $R(\rho+1:n,\rho+1:n)=\mathbf{0}$. If $k\in\{1,\ldots,n\}$, and if we introduce the block partition

[TABLE]

[TABLE]

If $k$ is such that $R_{[11]}$ is of full rank, then $A+\Delta A$ is of rank $k$ and

[TABLE]

Hence, if $\gamma>0$ is a given threshold, and if we can find an index $k$ ($1\leq k<n$) such that

[TABLE]

then $A$ is $\gamma$-close to the rank $k$ matrix $A+\Delta A$, whose pivoted QR factorization (2.3) is obtained from (2.1) by setting in the partition (2.2) the block $R_{[22]}$ to zero. Of course, we would take the smallest possible $k$ that satisfies (2.5). The essence of rank revealing capability of the factorization is in the fact that such a $k$ will very likely be visible on the diagonal of $R$ if $A$ is close to a rank $k$ matrix. This is due to the fact that the $|R_{ii}|$’s mimic the distribution of the singular values of $A$, [31], [32].

Remark 2.1**.**

The rank–$k$ approximation $A+\Delta A$ defined in (2.3) in general does not share the optimality property of the closest rank $k$ matrix $A_{k}$ from the SVD, but it has one distinctive feature: it always matches $A$ exactly at the selected $k$ columns, while $A_{k}$ in general does not match any part of $A$.

Remark 2.2**.**

An efficient and numerically reliable implementation of (2.1) is available e.g. in LAPACK [1] in the function xGEQP3, which is also under the hood of the Matlab’s function qr(); for the numerical and software details we refer to [10].

2.1.1 Backward stability

Computation of the QR factorization in finite precision arithmetic is backward stable: for the computed factors $\widetilde{P}$, $\widetilde{Q}$, $\widetilde{R}$, there exist a backward error $\delta A$ and a unitary matrix $\widehat{Q}$ such that

[TABLE]

Here $\epsilon_{qr}$ is bounded by a moderate function of the matrix dimensions times the machine roundoff $\varepsilon$. In fact, the backward stability is even stronger – the backward error in each column is small relative to its norm,

[TABLE]

This is an important feature if some columns of $A$ are, by its nature, much smaller than the largest ones (different weighting factors, different physical units); (2.7) assures that the computed factorization contains the information carried by small columns of $A$. While (2.6, 2.7) hold independent of pivoting, pivoting is important for the accuracy of the computed factorization, and for the rank revealing.

If, for a suitable partition of $\widetilde{R}$, analogous to (2.2), we can determine $k$ such that $\widetilde{R}_{[22]}$ can be chopped off, we have

[TABLE]

2.1.2 Complete pivoting

In some applications (e.g. weighted least squares) the rows of the data matrix may vary over several orders of magnitude, and it is desirable to have backward error that can be bounded row-wise analogously to (2.7). A pioneering work is done by Powell and Reid [30], who introduced QR factorization with complete pivoting. More precisely, in a $j$-th step, before deploying the Householder reflector to annihilate below-diagonal entries in the $j$-th column, row swapping is used to bring the absolutely largest entry to the diagonal position. As any pivoting, this precludes efficient blocking and using BLAS 3 level primitives.

[4] noted that the complete pivoting can be replaced with an initial sorting of the rows of $A$ in monotonically decreasing order with respect to the $\ell_{\infty}$ norm. If $P_{r}$ is the corresponding row permutation matrix, and if we set $A:=P_{r}A$, then

[TABLE]

and we proceed with the column pivoted factorization (2.1). An advantage of replacing the dynamic complete pivoting of Powell and Reid with the initial pre-sorting (2.10) followed by column pivoted QRF (2.1) is more efficient software implementation. An error analysis of this scheme and Householder reflector based QR factorization is given by Cox and Higham [7].

[TABLE]

and $\widetilde{A}^{(k)}$ is the $k$th computed (in finite precision arithmetic) intermediate matrix in the Householder QR factorization. As a result of initial row ordering and the column pivoting, [7] shows that

[TABLE]

where the factor $(1+\sqrt{2})^{n-1}$ is almost never experienced in practice.

Remark 2.3**.**

If we write the completely pivoted factorization as $P_{r}AP_{c}=Q\left(\begin{smallmatrix}R\cr\mathbf{0}\end{smallmatrix}\right)$, then $AP_{c}=(P_{r}^{T}Q)\left(\begin{smallmatrix}R\cr\mathbf{0}\end{smallmatrix}\right)$ is the column pivoted QR factorization (since $P_{r}$ is orthogonal); hence we can ease the notation of the complete pivoting and absorb the row permutation in the orthogonal factor, $Q\leftarrow P_{r}^{T}Q$. Note that the row pivoting brings nothing new to the rank revealing property that is encoded in the triangular factor (this is because of the essential uniqueness of the factorization). However, row pivoting has the important role in the numerical accuracy of the computed factorization.

2.2 The complete orthogonal factorization (ULV, URV)

Suppose that in the QR factorization (2.1), the matrix $A$ is of rank $k<n$, so that in the block partition (2.2) $R_{[22]}=\mathbf{0}$. In many instances, it is convenient to compress the trapezoidal matrix $(R_{[11]}\;R_{[12]})$ to (upper or lower) triangular form by an additional RQ or LQ factorization. The resulting complete orthogonal decomposition, often used in least squares computation [22], [19], is one of the computational subroutines in the quadeig algorithm.

The LQ factorization of $(R_{[11]}\;R_{[12]})$ is equivalent to computing the QR factorization

[TABLE]

where $T_{11}\in\mathbb{C}^{k\times k}$ is lower triangular and nonsingular. By a composition of these two steps we get the so called complete orthogonal decomposition of $A$

[TABLE]

In the QR factorization, the matrix $A$ is multiplied from the left by a sequence of unitary transformations. Hence, there is no mixing of the columns; we can analyse the process for each column separately; that is why the column-wise backward error bound (2.7) is natural and straightforward to derive. The transformations from the right in the pivoted QR factorization are the error free column interchanges.

On the other hand, (2.14) involves nontrivial two–sided transformations of $A$, and more careful implementation and error analysis are needed to obtain backward stability similar to the one described in §2.1. More precisely, the backward error in $A$ of the type (2.6) is a priori assured, because the transformations are unitary/orthogonal, but finer structure such as (2.7) requires some more work. A sharper backward error bound is attractive because the complete orthogonal decomposition is used in the reduction procedure in quadeig, and in our new proposed method.

Based on the discussion in §2.1.2, we compute with complete pivoting both the rank revealing QR factorization of $A$ and (2.13). The resulting scheme for computing the complete orthogonal decomposition is summarized in Algorithm 1.

2.2.1 Backward error analysis

The following theorem provides a backward error bound for the Algorithm 1. To our best knowledge, this approach to implementing and analyzing the URV/ULV is new.

Theorem 2.1**.**

Let $k$ be the numerical rank determined in step 2., with the computed trapezoidal (or triangular) factor $\widetilde{R}\in\mathbb{C}^{k\times n}$. Let $\widetilde{Q}$, $\widetilde{T}_{[11]}\in\mathbb{C}^{k\times k}$ and $\widetilde{Z}$ be the computed factors of complete orthogonal decomposition of $A$. Then there exist backward perturbations $\delta A$, $\Delta A$, $\delta\widetilde{R}$ and exactly unitary matrices $\widehat{Q}\approx\widetilde{Q}$, $\widehat{Z}\approx\widetilde{Z}$ that define the complete orthogonal decomposition

[TABLE]

where $\delta A$ is estimated as in (2.7), $\Delta A$ is as in (2.9), with an estimate as in (2.4), and $\|\delta\widetilde{R}(:,i)\|_{2}\leq\varepsilon_{3}(1+\varepsilon_{1})\|{A}(:,i)\|_{2}$ is bounded analogously to (2.11, 2.12).

Proof.

For the first step we have the relation (2.8), where the numerical rank $k$ is determined using a suitable thresholding. Here it is important to note that $\delta A$ satisfies (2.7), and that $\Delta A$ is zero at the $k$ leading pivotal columns selected by $\widetilde{P}$ (see (2.9)).

Set $\widetilde{R}=\begin{pmatrix}\widetilde{R}_{[11]}&\widetilde{R}_{[12]}\end{pmatrix}\in\mathbb{C}^{k\times n}$ and compute the QR factorization of $\widetilde{R}^{*}$. Here, we deploy complete pivoting for the following reason. We want the overall backward error to be small in each column of $A$ (or at least in the $k$ most important pivotal columns), relative to the norm of that particular column. To that end, first note that the computed factors $\widetilde{Z}$, $\widetilde{T}^{*}_{[11]}$ satisfy

[TABLE]

where $\widehat{Z}$ is unitary, $\|\widetilde{Z}-\widehat{Z}\|_{F}\leq\epsilon_{2}$ and, inserting (2.15) in (2.8), we obtain

[TABLE]

where $\delta A$ is from (2.6, 2.7), $\Delta A$ is as in (2.9). Now, the key observation is that by (2.11, 2.12), we can conclude that $\|\delta\widetilde{R}(:,i)\|\leq\epsilon_{3}\|\widetilde{R}(:,i)\|$, where, by (2.6) and (2.7), the column norms of $\widetilde{R}$ are bounded by $\|\widetilde{R}(:,i)\|\leq(1+\epsilon_{1})\|A(:,i)\|_{2}$.

Hence, the $k$ pivotal columns of $A$ (as selected by $\widetilde{P}$) have, individually, small backward errors of the type (2.7). Note that the complete pivoting in (2.15) is essential for column-wise small backward error in $\widetilde{R}$ and thus in $A$. ∎

2.3 Backward error in rank revealing QR factorizations of $M$ and $K$

The tool of trade in quadeig–type preprocessing is the pivoted QR factorization

[TABLE]

The initial row sorting before the column pivoted QR factorization is indicated by the matrices $P_{r,M}$, $P_{r,K}$. In the absence of row sorting $P_{r,M}=P_{r,K}=\mathbb{I}_{n}$.

In finite precision (see §2.1) the computed matrices $\widetilde{Q}_{M}$, $\widetilde{R}_{M}$, $\widetilde{\Pi}_{M}$ satisfy, independent of the choice of the permutation matrix $P_{r,M}$, $P_{r,M}(M+\delta M)\widetilde{\Pi}_{M}=\widehat{Q}_{M}\widetilde{R}_{M}$, where $\widehat{Q}_{M}=\widetilde{Q}_{M}+\delta\widetilde{Q}_{M}$ is exactly unitary with

[TABLE]

If $P_{r,M}$ is row sorting permutation, then, in addition, $\|\delta M(i,:)\|_{2}\leq\epsilon_{qr}^{\rightrightarrows}\|M(i,:)\|_{2}$, for $i=1,\ldots,n$. Since $P_{r,M}$ is unitary, we can absorb it into $\widetilde{Q}_{M}$ and $\widehat{Q}_{M}$ and redefine $\widetilde{Q}_{M}:=P_{r,M}^{T}\widetilde{Q}_{M}$, $\widehat{Q}_{M}:=P_{r,M}^{T}\widehat{Q}_{M}$ and write

[TABLE]

Analogous statements hold for the QRF $(P_{r,K}K)\Pi_{K}=Q_{K}R_{K}$, $R_{K}=\left(\begin{smallmatrix}\widehat{R}_{K}\cr\mathbf{0}_{n-r_{K},n}\end{smallmatrix}\right)$.

2.3.1 Backward error induced by truncation

If we truncate the triangular factor in an attempt to infer the numerical rank, we must push the truncated part into the backward error, as in (2.3). This changes the backward error structure, and the new error bounds depend on the truncation strategy and the threshold $\tau$.

Assume in (2.18) that we can partition $\widetilde{R}_{M}$ as

[TABLE]

Then we can write a backward perturbed rank revealing factorization

[TABLE]

Obviously, $\Delta M$ is zero at the $k$ "most linearly independent" columns of $M$ selected by the pivoting, $(\Delta M)\widetilde{\Pi}_{M}(:,1:k)=\mathbf{0}_{n,k}$. At the remaining $n-k$ columns we have

[TABLE]

where \zeta=1+O(n\mbox{\boldmath\varepsilon}) is a correction factor because the pivoting in floating point arithmetic cannot enforce the conditions (2.1) exactly. Consider the following choices of $k$:

$k$ is the first index for which $|((\widetilde{R}_{M})_{[22]})_{k+1,k+1}|\leq\tau|((\widetilde{R}_{M})_{[22]})_{k,k}|$. In that case

[TABLE]

• 2.

$k$ is the first index for which $|((\widetilde{R}_{M})_{[22]})_{k+1,k+1}|\leq\tau\cdot computed(\|M\|_{F})$. In that case

[TABLE]

• 3.

$k$ is the first index with $|((\widetilde{R}_{M})_{[22]})_{k+1,k+1}|\leq\tau\cdot computed(\max\{\|M\|_{F},\|C\|_{F},\|K\|_{F}\})$. In that case $\max_{j=1:n-k}\|(\Delta M)\widetilde{\Pi}_{M}(:,k+j)\|_{2}\leq\tau\cdot computed(\max\{\|M\|_{F},\|C\|_{F},\|K\|_{F}\})$. This strategy is used in quadeig with \tau=n\mbox{\boldmath\varepsilon}.

Remark 2.4**.**

If the truncation strategy [math] is deployed, as in quadeig, it is necessary to assume that the coefficient matrices have been scaled so that their norms are nearly equal. Otherwise, such a truncation strategy may discard a block in $\widetilde{R}_{M}$ (and declare $M$ numerically rank deficient) because it is small as compared e.g. to $\|C\|_{F}$ or $\|K\|_{F}$. On the other hand, in quadeig, scaling the matrices is optional, and if the (also optional) deflation procedure is enabled, this opens a possibility for catastrophic error (severe underestimate of the numerical ranks) if the matrices are not scaled and if their norms differ by orders of magnitude.

3 Kronecker’s canonical form for general pencils

Canonical (spectral) structure of a matrix pencil $A-\lambda B$ is, through linearization, an extremely powerful tool for the analysis of quadratic pencils $Q(\lambda)=\lambda^{2}M+\lambda C+K$. In particular, since the second companion form (1.2) is strong linearization, the partial multiplicities, and thus the structure of all eigenvalues (including infinity), are preserved. In a numerical algorithm for the QEP, it is desirable to remove the zero and infinite eigenvalues as early as possible and, thus, canonical structure can be used to guide such a preprocessing step.

3.1 Van Dooren’s algorithm

In this section, we briefly review the numerical algorithm by Van Dooren [9], developed for the computation of the structure of eigenvalue $\lambda$ for a general $m\times n$ pencil $A-\lambda B$, i.e. the number and the orders of the Jordan blocks for $\lambda$. Here, however, we focus our attention only on computing the structure of the eigenvalue $\lambda=0$. Notice that for the infinite eigenvalue one can reverse the pencil. For an arbitrary finite eigenvalue, a suitably shifted pencil is used; see [9].

For the sake of completeness and later references, we briefly describe the main steps of the staircase reduction for the zero eigenvalue. The pencil $A-\lambda B$ is assumed regular (thus square, $n\times n$), and $\lambda=0$ is assumed to be among its eigenvalues.

Compute the singular value decomposition of $A$: $A=U_{A}\Sigma_{A}V^{*}_{A}$, and let $s_{1}=n-\mathop{\rm rank}\nolimits(A)$. (Since zero is assumed to be an eigenvalue of $A-\lambda B$, $A$ must be column rank deficient.) Note that AV_{A}=\left(\begin{array}[]{c | c}A_{2}&\mathbf{0}_{n\times s_{1}}\end{array}\right), where $A_{2}$ is of full column rank $n-s_{1}$. Partition BV_{A}=\left(\begin{array}[]{c | c}B_{2}&B_{1}\end{array}\right) in the compatible manner. If we multiply the pencil by $V_{A}$ from the right we get (A-\lambda B)V_{A}=\left(\begin{array}[]{c | c}A_{2}-\lambda B_{2}&-\lambda B_{1}\end{array}\right).

• 2.

Compute the SVD B_{1}=U_{B}\Sigma_{B}V^{*}_{B},\;\;\mbox{and set}\;\;U^{*}_{B}B_{1}=\left(\begin{array}[]{c}B_{1,1}\cr\hline\cr\mathbf{0}_{n-s_{1}\times s_{1}}\end{array}\right). Multiply the pencil $(A-\lambda B)V_{A}$ by $U^{*}_{B}$ from the left to get

[TABLE]

• 3.

Let $P_{B}$ be the permutation matrix that swaps the row blocks in the above partition. Thus, we have unitary matrices $P_{1}=P_{B}U^{*}_{B}$, $Q_{1}=V_{A}$ so that

[TABLE]

where $\left(\begin{smallmatrix}A_{2,2}\cr A_{2,1}\end{smallmatrix}\right)=P_{1}A_{2}\in\mathbb{C}^{n\times(n-s_{1})}$ is of full column rank.

This concludes the first step of the algorithm. Note that he rank of $B_{1}$ is $s_{1}$ (full column rank) since the initial matrix pencil is assumed regular, and

[TABLE]

$\mathop{\rm det}\nolimits B_{1,1}\neq 0$, which clearly exposes $s_{1}$ copies of zero in the spectrum, and reduces the problem to the pencil $A_{2,2}-\lambda B_{2,2}$ of lower dimension $n_{2}=n-s_{1}$. Also, we see how rank deficiency in $B_{1,1}$ implies singularity of the input pencil. Clearly, if $A_{22}$ is nonsingular, zero has been exhausted from the spectrum of $A-\lambda B$.

Otherwise, in the next step, we repeat the described procedure on the $n_{2}\times n_{2}$ pencil $A_{2,2}-\lambda B_{2,2}$ to obtain unitary matrices $\widehat{P}_{2}$, $\widehat{Q}_{2}$ so that444Here we abuse the notation slightly. In any block partitioned pencil, the pair of indices $(i,j)$ is interpreted as follows: $i$ is the column index in a block partition, counted from the right to the left; $j$ is the row index, counted bottom up, $(\leftarrow\mbox{column},\uparrow\mbox{row})$. Each such index pair applies to the current stage in the algorithm. So, $A_{2,1}-\lambda B_{2,1}$ in (3.3) is the $(\leftarrow 2,\uparrow 1)$ position in $P_{2}P_{1}(A-\lambda B)Q_{1}Q_{2}$, which is different from $A_{2,1}-\lambda B_{2,1}$ in (3.2). This notation reflects the direction of building the final result, from bottom right up.

[TABLE]

where $P_{2}=\mathop{\rm diag}\nolimits(\widehat{P}_{2},\mathbb{I}_{s_{1}}),Q_{2}=\mathop{\rm diag}\nolimits(\widehat{Q}_{2},\mathbb{I}_{s_{1}})$, and $s_{2}=n_{2}-\mathrm{rank}(A_{22})$, $n_{3}=n_{2}-s_{2}$. As in the first step, $B_{2,2}$ is nonsingular, and $\left(\begin{smallmatrix}A_{3,3}\cr A_{3,2}\end{smallmatrix}\right)$ is of full column rank. Furthermore, since $\left(\begin{smallmatrix}\mathbf{0}_{n_{3}+s_{2},s_{2}}\cr A_{2,1}\end{smallmatrix}\right)$ is a column block in the full column rank matrix, $A_{2,1}$ must have full column rank as well. This procedure is repeated until in an $\ell$th step we obtain

[TABLE]

with nonsingular $A_{\ell+1,\ell+1}$. Then, by design, $B_{i,i}$ has full rank $s_{i}$ for $i=1,\ldots,\ell$, and $A_{i,i-1}$ has full column rank $s_{i}$ for $i=2,\ldots,\ell$. See [9, Algorithm 3.1].

Proposition 3.1**.**

[9, Proposition 3.5]** The indicies $s_{j}$ given by Algorithm [9, Algorithm 3.1] completely determine the structure at $\lambda=0$ of the pencil $A-\lambda B$, i.e. $A-\lambda B$ has $s_{j}-s_{j+1}$ elementary divisors $\lambda^{j}$, $j=1,\ldots,\ell$.

Finally, we can conclude that this algorithm also determines the structure of zero eigenvalue for the quadratic eigenvalue problem via a (strong) linearization.

Theorem 3.1**.**

If [9, Algorithm 3.1] is applied to the pencil (1.2), then it completely determines the structure of the eigenvalue zero for the quadratic eigenvalue problem $Q(\lambda)x\equiv(\lambda^{2}M+\lambda C+K)x=\mathbf{0}$.

3.2 Computing Kronecker’s Canonical form using rank revealing QR factorization

It is clear that the key rank reveling in the Van Dooren’s algorithm can be done using the rank revealing QR factorization instead of the SVD. In §2.1 we argued that in practice rank revealing QR factorization is considered as reliable as the SVD, and that it may even have some advantages. If we replace the SVD $A=U_{A}\Sigma_{A}V_{A}^{*}$ with $A=Q_{A}R_{A}P_{A}^{T}$, then the post-multiplication by $V_{A}$ in the procedure described in §3.1 is replaced with column permutation, i.e. post-multiplication with the permutation matrix $P_{A}$, which is error free even in finite precision arithmetic. Also, we can rearrange the elimination process to obtain an upper block triangular matrix as follows.

Compute the rank revealing factorization of $A_{1,1}=A$ $A_{1,1}P_{A}=Q_{A}R_{A}$, and denote $s_{1}=n_{1}-\mathop{\rm rank}\nolimits(A)=n-\mathop{\rm rank}\nolimits(A)$. Now, Q_{A}^{*}A_{1,1}=\left(\begin{array}[]{c}A_{2}\cr\hline\cr\mathbf{0}_{s_{1}\times n}\end{array}\right). Partition Q^{*}_{A}B=\left(\begin{array}[]{c}B_{2}\cr\hline\cr B_{1}\end{array}\right) in compatible manner. Multiply the pencil $(A-\lambda B)$ by $Q^{*}_{A}$ on the left to get

[TABLE]

• 2.

Compute complete orthogonal decomposition of $B_{1}$, $B_{1}=U_{B}R_{B}V^{*}_{B}$. The column rank of $B_{1}$ is $s_{1}$, if the matrix pencil is regular, and B_{1}V_{B}=\left(\begin{array}[]{c| c}B_{1,1}&\mathbf{0}_{s_{1},n-s_{1}}\end{array}\right), where $B_{1,1}$ is upper triangular. Multiply the pencil (3.5) by $V_{B}$ on the right to get Q_{A}^{*}(A-\lambda B)V_{B}=\left(\begin{array}[]{c| c}A_{1,2}-\lambda B_{1,2}&A_{2,2}-\lambda B_{2,2}\cr\hline\cr\lambda B_{1,1}&\mathbf{0}\end{array}\right).

• 3.

Let $P_{B}$ be the permutation matrix for permuting the $s_{1}$ and $n-s_{1}$ column blocks. Define $P_{1}=Q^{*}_{A}$ and $Q_{1}=V_{B}P_{B}$. The first Jordan block for the eigenvalue $\lambda=0$ is deflated by the following orthogonal transformation:

[TABLE]

This reduction continues analogously to §3.1, but keeping the upper triangular structure; the details are in Algorithm 1.

The final form of $P(A-\lambda B)Q$ is given in (3.7). The deflation of infinite eigenvalues can be done by the same algorithm, but with the reversed pencil $B-\lambda A$.

[TABLE]

Remark 3.1**.**

For stronger backward stability, the rank revealing QR factorizations can be computed with full pivoting, e.g. $P_{A_{r}}A_{j,j}P_{A}=Q_{A}R_{A}$. The rest of the algorithm is changed simply by setting $Q_{A}\leftarrow P_{A,r}^{T}Q_{A}$, cf. Remark 2.3. The complete orthogonal decomposition is computed as described in §2.2.

4 Scaling and balancing

Backward stable computation of the eigenvalues of the linearized pencil (1.2) allows interpretation of the computed output as the exact spectral information of a pencil $A+\delta A-\lambda(B+\delta B)$, where $\delta A$ ($\delta B$) is small relative to $A$ ($B$) in an appropriate norm. Since $\delta A$, $\delta B$ are not structured (linear eigensolver that takes $A$ and $B$ on input is oblivious of the block structure and the original quadratic problem), the backward stability does not necessarily hold in terms of the original data $M$, $C$, $K$. In fact, it is easily seen that one potential source of the problem is when there is large variation in the norms of $M$, $C$, $K$, $\mathbb{I}_{n}$. This issue is addressed by parameter scaling.

4.1 Parameter scaling

The idea is to define a scaled problem $\widetilde{Q}(\mu)=\mu^{2}\widetilde{M}+\mu\widetilde{C}+\widetilde{K}$, which depends on two scalar parameters $\gamma$ and $\delta$, where $\lambda=\gamma\mu,\;\;\widetilde{Q}(\mu)=Q(\lambda)\delta=\mu^{2}(\gamma^{2}\delta M)+\mu(\gamma\delta C)+\delta K$. This scaling has no effect on backward error for QEP, but only on backward error for the linearization. In quadeig, [14] used two types of scaling.

Fan, Lin and Van Dooren scaling

[11] The scaling parameters $\gamma$ and $\delta$ are defined as the solution of minimization problem $\min_{\gamma,\delta}\max\{\|\widetilde{K}\|_{2}\!-\!1,\|\widetilde{C}\|_{2}\!-\!1,\|\widetilde{M}\|_{2}\!-\!1\}$, so that in the linearization (1.2) of $\widetilde{Q}(\mu)$ all nonzero blocks are of norm close to one. The solution of this minimization problem is $\gamma=\sqrt{{\|K\|_{2}}/{\|M\|_{2}}}$, $\delta={2}/{(\|K\|_{2}+\|C\|_{2}\gamma)}$.

Tropical scaling

$\gamma$ and $\delta$ are defined as tropical roots of max-times scalar quadratic polynomial $q_{\text{trop}}(x)=\max(\|M\|_{2}x^{2},\|C\|_{2}x,\|K\|_{2}),\;\;x\in[0,\infty\rangle.$ Define $\tau_{Q}=\frac{\|C\|_{2}}{\sqrt{\|M\|_{2}\|K\|_{2}}}$. If $\tau_{Q}\leq 1$, there is double root $\gamma^{+}=\gamma^{-}=\sqrt{\|K\|_{2}/\|M\|_{2}}$, and if $\tau_{Q}>1$ there are two distinct roots $\gamma^{+}={\|C\|_{2}}/{\|M\|_{2}}>\gamma^{-}={\|K\|_{2}}/{\|C\|_{2}}$. Hence, when $\tau_{Q}>1$, scaling with $\gamma=\gamma^{+},\;\;\delta=(q_{\text{trop}}(\gamma^{+}))^{-1}$ is used to compute the eigenvalues outside of the unit circle, and scaling with $\gamma=\gamma^{-},\;\;\delta=(q_{\text{trop}}(\gamma^{-}))^{-1}$ is used to compute the eigenvalue inside the unit circle. For further details on parameter scaling see [12].

4.2 Balancing

Backward error $\delta A$, such that $\|\delta A\|/\|A\|$ is small in some appropriate matrix norm may not guarantee high accuracy of the result (e.g. of the computed eigenvalues) if the corresponding condition number is too big. Also, the interpretation of the computed result in the backward error analysis framework may be difficult because small $\|\delta A\|/\|A\|$ does not guarantee that the backward relative error in the individual matrix entries is correspondingly small. This problem in particular escalates if the matrix entries vary over several orders of magnitude, and the smallest entries carry relevant information. Well known technique to cope with the unbalanced matrix entries is the technique of balancing: replace $A$ with a similar matrix $D^{-1}AD$, where $D$ is a diagonal matrix (or a permutation of diagonal matrix) that in some sense balances the matrix entries. For instance, for the non-symmetric eigenvalue problem, the balancing strategy proposed by [28], [29] ensures that for all $i$, the $i$th row and the $i$th column of $D^{-1}AD$ have (nearly) the same norm; usually the scaling factors $D_{ii}$ are replaced with closest power of the base of the machine arithmetic, so that the balancing step only shifts the exponents of the matrix entries, leaving the mantissas unchanged. For a software implementation see [20]. Another interesting strategy of is to make $D^{-1}AD$ closer to normal matrices [23]. It should be noted that balancing may be harmful, see [38].

Balancing has been recognized as an important preprocessing technique for eigensolvers for general pencils $A-\lambda B$, see e.g. [36], [23]. In some applications, it is convenient to apply it to matrix triplets. For instance, the frequency response matrix $C(\sigma E-A)^{-1}B$ of a MIMO descriptor linear system is efficiently computed for many values of $\sigma$ by using precomputed $m$-Hessenberg-triangular form of the triplet $(E,A,B)$. It is shown in [5] that the numerical robustness can be improved by pre-scaling $D_{l}AD_{r}$, $D_{l}ED_{r}$ and $D_{l}B$, where the diagonal matrices $D_{l}=\mathrm{diag}(10^{l_{i}})_{i=1}^{n}$ and $D_{r}=\mathrm{diag}(10^{r_{i}})_{i=1}^{n}$ are such that the range of magnitude orders of all elements in the matrices $D_{l}AD_{r}$, $D_{l}ED_{r}$ and $D_{l}B$ is minimal.

We will show that balancing can substantially improve numerical solution of the QEP, and that it should be considered as standard preprocessing technique, together with the parameter scaling. Balancing with $D_{l}$ and $D_{r}$ yields an equivalent QEP:

[TABLE]

For a computed eigenpair $(\lambda,x)$ of (4.1), the corresponding right eigenpair for the original problem is $(\lambda,D_{r}x)$. Note that this corresponds to the transformation of $A-\lambda B$ of the linearization (1.2) to $(D_{l}\oplus D_{l})(A-\lambda B)(D_{r}\oplus D_{l}^{-1}).$

The scaling changes the matrix entries so that e.g. $m_{ij}\mapsto 10^{l_{i}}m_{ij}10^{r_{j}}$, and logarithmic magnitude of the scaled element is $\log_{10}10^{l_{i}}|m_{ij}|10^{r_{j}}=l_{i}+r_{j}+\log_{10}|m_{ij}|$. The idea of balancing is to make these as close to zero as possible over all entries of all matrices, where the importance of each matrix is weighted with a corresponding factor $\alpha_{M}\geq 0$, $\alpha_{C}\geq 0$, $\alpha_{K}\geq 0$, respectively. This yields the minimization problem $\min_{l=(l_{i}),r=(r_{i})\in\mathbb{R}^{n}}\varphi(l,r)$, where

[TABLE]

which can be written in the generic form $\|Ax-b\|_{2}\rightarrow\min$ as follows.

Proposition 4.1**.**

For each index $i$ define row vectors $\mathbbm{1}^{\alpha_{M}}_{i}=\alpha_{M}\mathcal{I}(M(i,:))$, $\mathbbm{1}^{\alpha_{C}}_{i}=\alpha_{C}\mathcal{I}(C(i,:))$, $\mathbbm{1}^{\alpha_{K}}_{i}=\alpha_{K}\mathcal{I}(K(i,:))$, where $\mathcal{I}(\cdot)$ is the indicator vector, e.g. $(\mathcal{I}(M(i,:)))_{j}=1$ if $m_{ij}\neq 0$ and [math] otherwise. Further, define diagonal matrices $\mathbb{I}^{\alpha_{M}}_{i}=\mathrm{diag}(\mathbbm{1}^{\alpha_{M}}_{i})$, $\mathbb{I}^{\alpha_{C}}_{i}=\mathrm{diag}(\mathbbm{1}^{\alpha_{C}}_{i})$, $\mathbb{I}^{\alpha_{K}}_{i}=\mathrm{diag}(\mathbbm{1}^{\alpha_{K}}_{i})$ and column vectors in $\mathbb{R}^{n}$ $L^{\alpha_{M}}_{i}=\alpha_{M}[\log|M(i,:)|\circ\mathcal{I}(M(i,:))]^{T}$, $L^{\alpha_{M}}_{i}=\alpha_{M}[\log|M(i,:)|\circ\mathcal{I}(M(i,:))]^{T}$, $L^{\alpha_{M}}_{i}=\alpha_{M}[\log|M(i,:)|\circ\mathcal{I}(M(i,:))]^{T}$, where $\circ$ denotes the Hadamard product and $0\log 0$ is set to [math]. Then (4.2) is a linear least squares problem $\|Ax-b\|_{2}\to\text{min}$ with

[TABLE]

The corresponding system of normal equations, $Lx=p$, reads

[TABLE]

*with $n_{r_{i}}\!\!=\!\!\sum^{n}_{\begin{subarray}{c}j=1\\ m_{ij}\neq 0\end{subarray}}\alpha_{M}+\sum^{n}_{\begin{subarray}{c}j=1\\ c_{ij}\neq 0\end{subarray}}\alpha_{C}+\sum^{n}_{\begin{subarray}{c}j=1\\ k_{ij}\neq 0\end{subarray}}\alpha_{K},n_{c_{j}}\!\!=\!\!\sum^{n}_{\begin{subarray}{c}i=1\\ m_{ij}\neq 0\end{subarray}}\alpha_{M}+\sum^{n}_{\begin{subarray}{c}i=1\\ c_{ij}\neq 0\end{subarray}}\alpha_{C}+\sum^{n}_{\begin{subarray}{c}i=1\\ k_{ij}\neq 0\end{subarray}}\alpha_{K};$ $G\in\mathbb{R}^{n\times n}$ is the sum of incidence matrices of $M$, $C$ and $K$ scaled by $\alpha_{M},\alpha_{C}$ and $\alpha_{K}$: $G=\alpha_{M}\mathcal{I}(M)+\alpha_{C}\mathcal{I}(C)+\mathcal{I}(K)$, where e.g. $\mathcal{I}(M)_{ij}=1$ if $m_{ij}\neq 0$ and [math] otherwise. The entries of $c\in\mathbb{R}^{n}$ are $c_{i}=\alpha_{M}\sum^{n}_{\begin{subarray}{c}j=1\\ m_{ij}\neq 0\end{subarray}}\log|m_{ij}|+\alpha_{C}\sum^{n}_{\begin{subarray}{c}j=1\\ c_{ij}\neq 0\end{subarray}}\log|c_{ij}|+\alpha_{K}\sum^{n}_{\begin{subarray}{c}j=1\\ k_{ij}\neq 0\end{subarray}}\log|k_{ij}|$, and of $d\in\mathbb{R}^{n}$ are $d_{j}=\alpha_{M}\sum^{n}_{\begin{subarray}{c}i=1\\ m_{ij}\neq 0\end{subarray}}\log|m_{ij}|+\alpha_{C}\sum^{n}_{\begin{subarray}{c}i=1\\ c_{ij}\neq 0\end{subarray}}\log|c_{ij}|+\alpha_{K}\sum^{n}_{\begin{subarray}{c}i=1\\ k_{ij}\neq 0\end{subarray}}\log|k_{ij}|$. *

In the case of large scale sparse problem, the normal equations can be solved using the preconditioned conjugate gradients method, see e.g. [5].

In Figure 2, we illustrate the impact of balancing to the backward error in the case of one NLEVP benchmark example (damped beam). In general, this kind of preprocessing should also be used in the iterative methods.

5 Improved deflation process

In this section, we propose an extension of the quadeig reduction scheme toward the KCF outlined in §3.2. To motivate, we use an example from the NLEVP collection.

5.1 A case study example

This is a $10\times 10$ quadratic eigenvalue problem for the pencil $\mathcal{I}(\lambda)=\lambda^{2}M+\lambda C+K$, whose real eigenvalues and the corresponding eigenvectors give the intersection points of a sphere, a cylinder and a plane. Although of small dimension and very simple structure, this example is an excellent illustration of difficulties in solving nonlinear eigenvalue problems.

It has been shown in [27], [24] that $\mathcal{I}(\lambda)$ has only four finite eigenvalues: two real ones and a complex conjugate pair. We take this example as a case study and compute the spectrum by several mathematically equivalent methods. If one plainly applies the QZ to a linearization of $\mathcal{I}(\lambda)$, such as the first or the second companion form with the Fan-Lin-Van Dooren scaling, the spectrum appears as

[TABLE]

If we use the same method, but with the reversed pencil $\mu^{2}K+\mu C+M$, ($\lambda=1/\mu$) then from the first companion form QZ has computed $12$ finite eigenvalues, and from the second $10$. If we run the Matlab’s solver polyeig(), we obtain

[TABLE]

and if we scale the coefficient matrices then

[TABLE]

Almost perfect match in $\lambda_{1}$ and $\lambda_{2}$ is reassuring, but there is an obvious disagreement in the total number and the nature (real or complex) of finite eigenvalues. With an earlier version of Matlab, the results that correspond to (5.5), (5.8), (5.9), and (5.10) coincide in the numbers of finite eigenvalues; $\lambda_{1}$ and $\lambda_{2}$ are close up to machine precision, but the remaining computed finite eigenvalues are substantially different.

The rank of the matrix $M$ is exactly $3$, and it will be correctly determined numerically due to a particularly simple sparsity structure of $M$. The matrix $K$ is also sparse with $\kappa_{2}(K)\approx\texttt{4.09+03}$, so there is no numerical rank issue. In this situation, a preprocessing procedure such as in quadeig will reverse the pencil and deflate $7$ zero eigenvalues (infinite eigenvalues of the original problem) at the very beginning. The remaining eigenvalues are then computed (e.g. using quadeig) as555Here, to save the space, we display the computed values only to five digits.

[TABLE]

After the deflation of the $7$ zero eigenvalues, in the thus obtained linearization $A-\lambda B$, the rank of the matrix $A$ is $7$, and it can be determined exactly because of sparsity ($A$ has $6$ zero columns, and the remaining $7$ ones build a well conditioned $13\times 7$ submatrix of $A$). The matrix $B$ is well conditioned. This means that at least $6$ more zero eigenvalues are present in the reversed problem (infinities in the original problem); those are not detected by the QZ algorithm running on $A-\lambda B$.

Remark 5.1**.**

It should be noted that the successful removal of many infinite eigenvalues in (5.5), (5.8), (5.9), and (5.10) is due to the sparsity of $M$, $C$, $K$ (and the linearization $A-\lambda B$) that is successfully exploited by the preprocessing in the QZ algorithm. Recall, before the reduction to the triangular - Hessenberg form the matrices are scaled and permuted, as described in [36] in order to get an equivalent pencil $\hat{A}-\lambda\hat{B}$ with the block structure

[TABLE]

where $A_{11},A_{33},B_{11},B_{33}$ are upper triangular, and the matrices $D_{1}$ and $D_{2}$ are computed so that the elements of $D_{1}A_{22}D_{2}$ and $D_{1}B_{22}D_{2}$ have magnitudes close to one. Further, $G_{2}$ is the permutation matrix determined so that the ratios of the column norms of $D_{1}A_{22}D_{2}G_{2}$ to the corresponding column norms of $D_{1}B_{22}D_{2}G_{2}$ appear in decreasing order; $G_{1}$ is determined so that the ratios of row norms of $G_{1}D_{1}A_{22}D_{2}G_{2}$ to those of $G_{1}D_{1}B_{22}D_{2}G_{2}$ appear in decreasing order. See also [21].

5.2 Deflation process revisited – a modified scheme

Let us now get back to the initial idea of quadeig – removing the zero and the infinite eigenvalues in the preprocessing phase of the computation. To introduce a new scheme, consider removing the zeros from the spectrum and assume that a rank revealing QR factorization $KP_{K}=Q_{K}\left(\begin{smallmatrix}R_{K,1}\cr\mathbf{0}_{n-r_{K},n}\end{smallmatrix}\right)=Q_{K,1}R_{K,1},$ $Q_{K}=(Q_{K,1},Q_{K,2})$, has detected that $K$ is numerically rank deficient. Consider the transformation

[TABLE]

Note that $\mathrm{rank}(A)=n+\mathrm{rank}(K)$, so $A$ and $K$ have null spaces of equal dimensions. In essence, multiplication from the left with $Q_{K}^{*}\oplus Q_{K}^{*}$ (or with $Q_{M}^{*}\oplus Q_{K}^{*}$, or $\mathbb{I}_{n}\oplus Q_{K}^{*}$) is a rank revealing transformation of $A$. We now truncate the $s_{1}=n-r_{K}$ copies of the eigenvalue $\lambda=0$ and proceed with the truncated $(n+r_{K})\times(n+r_{K})$ pencil

[TABLE]

Note that using the definition (5.24) of $A_{22}-\lambda B_{22}$ in (5.23) yields

[TABLE]

At this point, quadeig removes the $n-r_{K}$ zero eigenvalues and leaves to the QZ to find the remaining ones, if present in the reduced problem. However, note that with $A_{11}:=A$ and $B_{11}:=B$, this procedure can be understood as the first step of the Van Dooren’s algorithm (actually its transposed version described in §3.2; see (3.6)), which removes the first Jordan block of the eigenvalue zero. The existence of further Jordan blocks of $\lambda=0$ can be established immediately by checking the rank of $A_{22}$.

From the definition of $A_{22}$ in (5.24), it follows that it suffices to compute a rank revealing QR factorization of its $n\times n$ submatrix666This is why we have not followed the strategy (6.20). $A_{22}(r_{K}+1:n+r_{K},1:n)$,

[TABLE]

This can be used to transform the pencil $A_{22}-\lambda B_{22}$ to

[TABLE]

If the factorization (5.26) shows no rank deficiency, there are no zeros in the spectrum of $A_{22}-\lambda B_{22}$. Equivalently, the reversed pencil $B_{22}-\mu A_{22}$ has no infinite eigenvalues.

Otherwise, $A_{22}(r_{K}+1:n+r_{K},1:n)$ is rank deficient; assume its rank to be $r_{22}$, $r_{22}<n$, and let $s_{2}=n+r_{K}-r_{22}$. Then $R_{A_{22}}=\left(\begin{smallmatrix}{R}_{A_{22},1}\cr\mathbf{0}_{n-r_{22},n}\end{smallmatrix}\right),\;\;{R}_{A_{22},1}\in\mathbb{C}^{r_{22}\times n},$ and

[TABLE]

To deflate the additional $m-r_{22}$ zeros, we need $(n-r_{22})\times(r_{K}+r_{22})$ zero block in the lower left corner of the matrix $\widehat{\mathbf{P}}_{2}B_{22}$. This is achieved using the complete orthogonal decomposition \left(\begin{array}[]{c|||c}\square&\triangle\end{array}\right)=U_{B}R_{B}V^{*}_{B} of the $(n-r_{22})\times(n+r_{k})$ block \left(\begin{array}[]{c|||c}\square&\triangle\end{array}\right) in (5.28).

If we write \left(\begin{array}[]{c|||c}\square&\triangle\end{array}\right)V_{B}=\left(\begin{array}[]{c c}\Diamond_{n-r_{22},n-r_{22}}&\mathbf{0}_{n-r_{22},r_{K}+r_{22}}\end{array}\right), $U_{B}R_{B}=\Diamond_{n-r_{22},n-r_{22}}$, and define the permutation $P_{B}$ to swap the first $(n-r_{22})$ and the remaining $(r_{K}+r_{22})$ columns, then deploying the transformation matrix

[TABLE]

from the right yields the targeted structure

[TABLE]

In (5.30), we slightly abuse the notation and denote the block $U_{B}R_{B}$ by $B_{22}$, to emphasize that this procedure follows the reduction from §3.2. If we premultiply (5.30) with $\mathbb{I}_{t_{K}+r_{22}}\oplus U_{B}^{*}$, then we can set $B_{22}=R_{B}$. Note that in (5.30) nonsingularity of $B_{22}$ (i.e. full row rank of the matrix block \left(\begin{array}[]{c|||c}\square&\triangle\end{array}\right) in (5.28)) is necessary for the regularity of the pencil. If the overall algorithm is designed to consider singular pencil an invalid input, then instead of the complete orthogonal decomposition we can use rank revealing LQ factorization with row (or complete) pivoting; in that case $U_{B}$ is a permutation matrix, and $R_{B}$ is lower triangular.777In this case, by additional permutations $R_{B}$ can be made upper triangular.

After this step, the structure of the linearization pencil is lost. To determine the existence of more zero eigenvalues, the rank of the $(r_{K}+r_{22})\times(r_{K}+r_{22})$ matrix $A_{33}$ is needed. For the further deflation we simply switch to the KCF computation as in §3.2.

5.3 Backward error

Since all involved transformations are numerically unitary, the classical backward error analysis of the transformations of $A-\lambda B$ in a unitary invariant matrix norm is straightforward. It is important, however, to keep the backward errors in terms of the original matrices $M$, $C$, $K$ (instead of the matrices $A$, $B$ of the linearization) and to estimate them in finer resolution – instead of measuring the backward error in terms of matrix (operator) norm, the goal is to obtain the size of relative error in individual columns, rows or matrix entries.

Such an analysis is usually technically involved and tedious. Here we provide the key elements of the analysis of the first two steps of the removal of zero eigenvalues. The idea is to provide insight, and to show how to keep small backward errors in the matrices of the quadratic problem, and to ensure it is column–wise small.

Theorem 5.1**.**

The reduction described in §5.2 computes (5.30) with small backward errors in the initial data matrices $M$, $C$ and $K$.

We provide the proof with detailed discussion in the following two subsections, which also contain the bounds on the backward errors.

5.3.1 The first step

The backward error of the rank revealing factorization of $K$ is as explained in §2.3; it contains both the floating point error $\delta K$ and the truncation error $\Delta K$ analogous to (2.19), i.e. $(K+\delta K+\Delta K)\widetilde{P}_{K}=\widehat{Q}_{K}\widetilde{R}_{K}=\widehat{Q}_{K,1}\widetilde{R}_{K,1}$. Set $\Delta_{\Sigma}K=\delta K+\Delta K$, and introduce block partitions $\widehat{Q}_{K}=(\widehat{Q}_{K,1},\widehat{Q}_{K,2})$, where $\widehat{Q}_{K,1}$ has $\widetilde{r}_{K}$ columns, and $\widetilde{r}_{K}$ is the computed numerical rank of $K$.

Let $\widetilde{X}_{11}=computed(\widetilde{Q}_{K}^{*}C)$. By the standard error analysis of floating point matrix multiplication we know hat $computed(\widetilde{Q}_{K}^{*}C)=\widetilde{Q}_{K}^{*}C+\mathfrak{G}_{C},\;\;|\mathfrak{G}_{C}|\leq\epsilon_{*}|\widetilde{Q}_{K}^{*}||C|$, where 0\leq\epsilon_{*}\leq O(n)\mbox{\boldmath\varepsilon}. Here we assume that $\widetilde{Q}_{K}$ is explicitly computed888As in LAPACK’s subroutines xORGQR. and then applied. If $\widetilde{Q}_{K}$ is stored and applied in factored form (using Householder vectors),999See e.g. LAPACK’s subroutines xORMQR. similar relation holds, but with more technical details, which we omit for the sake of brevity. Since we can represent the numerically unitary $\widetilde{Q}_{K}$ as $\widetilde{Q}_{K}=(\mathbb{I}+\mathfrak{E}_{C})\widehat{Q}_{K}$, $\|\mathfrak{E}_{C}\|_{2}\leq\epsilon_{qr}$,

[TABLE]

with column-wise estimates $\|\delta C(:,i)\|_{2}\leq(\|\mathfrak{E}_{C}^{*}\|_{2}+\epsilon_{*}\sqrt{n}(1+\|\mathfrak{E}_{C}^{*}\|_{2}))\|C(:,i)\|_{2}$.

By the same reasoning we get $\widetilde{Y}_{11}\equiv computed(\widetilde{Q}_{K}^{*}M)=\widehat{Q}_{K}(M+\delta M)$, where for each $i=1,\ldots,n$, $\|\delta M(:,i)\|_{2}\leq\epsilon_{M}\|M(:,i)\|_{2}$, and $\epsilon_{M}=(\epsilon_{qr}+\epsilon_{*}\sqrt{n}(1+\epsilon_{qr})).$

Hence, we can represent the computed equivalent of (5.23) as

[TABLE]

5.3.2 The second step

In the next step, we compute the rank revealing factorization (5.26) of the computed matrix, or, equivalently,

[TABLE]

For the computed factors $\widetilde{\Pi}_{A},\widetilde{Q}_{A_{22}},\widetilde{R}_{A_{22}}$ of floating point implementation it holds that

[TABLE]

where $\widehat{Q}_{A_{22}}$ is exactly unitary and close to the computed unitary factor $\widetilde{Q}_{A_{22}}$, and $\left(\begin{smallmatrix}\mathfrak{C}\\ \mathfrak{K}\end{smallmatrix}\right)$ is the backward error. The last relation in (5.35) illustrates how to push the backward errors $\mathfrak{C}$, $\mathfrak{K}$ into the original data.

[TABLE]

This allows a mixed stability interpretation indicated by $"\approx"$ above: if the computed part of the matrix that contains $\widetilde{Q}_{A_{22}}$ is changed by using the corresponding submatrices of $\widehat{Q}_{A_{22}}$ instead, then we have backward stability.101010This mixed formulation is preferred because it allows leaving the identity block unchanged, which leaves the underlying structure of the quadratic pencil intact. Besides, $\widehat{Q}_{A_{22}}$ and $\widetilde{Q}_{A_{22}}$ are $O(\varepsilon_{qr})$ close. Further, the backward stability holds for the reduced matrix (the leading $(n+\widetilde{r}_{K})\times(n+\widetilde{r}_{K})$) block.

Consider now the $B$–part of the pencil. The effective change is

[TABLE]

where $\delta\mathfrak{Q}=\widetilde{Q}_{A_{22}}^{*}\!\!-\!\widehat{Q}_{A_{22}}^{*}$. Note that $\widetilde{Y}_{11}(1:\widetilde{r}_{K},1:n)=\widehat{Q}_{K,1}^{*}(M+\delta M)=\widehat{Q}_{K,1}^{*}(M+\delta M+\widehat{Q}_{K,2}\mathfrak{E})$. Hence, we can push $\widehat{Q}_{K,2}\mathfrak{E}$ backward into $M$ without affecting the first $\widetilde{r}_{K}$ rows in the $B$–part. Here too we have a mixed stability interpretation

[TABLE]

Hence, up to an $O(\varepsilon_{qr})$ forward error (5.87), the computed $B$–part satisfies the backward perturbation relation (5.58) with exactly unitary transformations. This, together with (5.36), (5.56) establishes a mixed stability interpretation of the second reduction step.

Mixed stability can be replaced by pure backward stability as follows. In (5.35), write $\widetilde{Q}_{A_{22}}\widetilde{R}_{A_{22}}$ instead of $\widehat{Q}_{A_{22}}\widetilde{R}_{A_{22}}$; this only changes the backward error $\left(\begin{smallmatrix}\mathfrak{C}\cr\mathfrak{K}\end{smallmatrix}\right)$. Then repeat (5.36)–(5.56) and (5.58)–(5.87) (with $\widehat{Q}_{A_{22}}\widetilde{R}_{A_{22}}$ replaced by $\widetilde{Q}_{A_{22}}\widetilde{R}_{A_{22}}$) with updated $\mathfrak{E}$ and $\delta\widehat{Y}_{21}=\mathbf{0}$.

5.3.3 Comments on the sizes of the backward errors $\mathfrak{C}$ and $\mathfrak{K}$

The backward error of the rank revealing QR factorization (5.35) is column-wise bounded by

[TABLE]

Here, for simplicity, we assume full rank case, so the backward error above does not contain the truncation part. The backward errors in $C$ and $K$ are then bounded by

[TABLE]

Since both $\|\mathfrak{C}(:,i)\|_{2}$, $\|\mathfrak{K}(:,i)\|_{2}$ contain mixing of $C$ and $K$, we see the benefits of scaling and balancing which then help keeping the error in each matrix small relative to its corresponding part/column.

6 Putting it all together: the global structure of KVADeig

Here, we describe the global structure of the new procedure. We assume that the initial scaling and balancing are done as requested by an expert user.

6.1 Full reduction scheme

We will always attempt to deflate zero eigenvalues. This means that if there are infinite eigenvalues, we will switch to the reversed problem. If both zero and infinite eigenvalues are present, we work with the problem in which the number of zero eigenvalues is higher, and the structure is used in the deflation process of this eigenvalue. However, the information about the known number of infinite eigenvalues is transmitted in the next steps of the Algorithm 1. The first step are the rank revealing decompositions of the matrices $M$ and $K$. Then we have the following cases.

1. Both matrices $M$ and $K$ are regular

In this case, there is no deflation. However, the rank revealing decomposition of $M$ is used to reduce the matrix $B$ to upper triangular form, as in quadeig, to facilitate more efficient run of the QZ algorithm.

[TABLE]

2. One of the matrices $M$ or $K$ is singular

Assume that $K$ is singular. Otherwise, switch to the reversed problem. Before we start the deflation process of the $n-\mathop{\rm rank}\nolimits(K)=n-r_{K}$ zero eigenvalues, we check whether there is more than one Jordan block for this eigenvalue, that is, we compute the numerical rank of the $n\times n$ block matrix $\left(\begin{smallmatrix}Q^{*}_{K,2}C\cr\hline\cr{R}_{K,1}P^{T}_{K}\end{smallmatrix}\right)$, see (5.26). Depending on the computed rank, we have two cases.

2.1. Full rank

There are exactly $n\!-\!r_{K}$ zero eigenvalues. The deflation is done as in quadeig, i.e. along with the deflation, $B$ is reduced to upper triangular form.

[TABLE]

The reduced $(n+r_{K})\times(n+r_{K})$ pencil is

[TABLE]

2.2. Rank deficient

This means that there is more than one Jordan block for the zero eigenvalue, more precisely there are at least $n-\mathop{\rm rank}\nolimits({{\left(\begin{smallmatrix}Q^{*}_{K,2}C\cr\hline\cr{R}_{K,1}P^{T}_{K}\end{smallmatrix}\right)}})+n-r_{K}=2n-r_{22}-r_{K}$ zero eigenvalues. The first two blocks are deflated as described in §5. Let (5.30) be the transformed pencil, after the two steps of the deflation. Now, the rank of the matrix $A_{33}$ determines whether there are more zero eigenvalues. The full deflation of zero eigenvalues is finished by calling the Algorithm 1 on the truncated pencil $A_{33}-\lambda B_{33}$. As an output of the algorithm, we get the reduced pencil $A_{\ell+1,\ell+1}-\lambda B_{\ell+1,\ell+1}$ with $A_{\ell+1,\ell+1}$ regular, and the corresponding transformation matrices $Q_{p}$ and $P_{p}$ of order $r_{K}+r_{22}$. The final transformation matrices $Q$ and $P$ are

[TABLE]

Now, the remaining nonzero eigenvalues are computed using the QZ algorithm on the truncated pencil $A_{\ell+1,\ell+1}-\lambda B_{\ell+1,\ell+1}$.

3. Both matrices $M$ and $K$ are singular

In this case, we first check the ranks of both block matrices $\left(\begin{smallmatrix}Q^{*}_{K,2}C\cr\hline\cr\widehat{R}_{K}P^{T}_{K}\end{smallmatrix}\right),\;\;\left(\begin{smallmatrix}Q^{*}_{M,2}C\cr\hline\cr\widehat{R}_{M}P^{T}_{M}\end{smallmatrix}\right).$ The ranks of these matrices determine whether there exist more than one Jordan block for the zero and the infinite eigenvalue, respectively. Depending on their regularity, we distinguish the following three cases.

3.1. Both matrices are regular

This means that there are exactly $n-r_{M}$ infinite, and $n-r_{K}$ zero eigenvalues. The deflation of these eigenvalues is done analogously to the quadeig algorithm. For the reader’s convenience, we provide the details. Let

[TABLE]

Note the difference in the transformation from the right with (6.19): instead of $\Pi_{M}$, we now have $\mathbb{I}_{n}$, so that $Q_{M}^{*}M=R_{M}\Pi_{M}^{T}$ is not upper triangular. Preserving the triangular form in this moment does not seem important because it is likely that it will be destroyed in subsequent steps. In the next step, the complete orthogonal decomposition (i.e. URV decomposition, using unitary matrices $Q_{X}$ and $Z_{X}$) is computed as

[TABLE]

Then (6.59) can be further transformed as follows:

[TABLE]

The $(1,1)$ diagonal block in the new partition ($=\!\parallel\!=$) is $(r_{M}+r_{K})\times(r_{M}+r_{K})$, and

• •

$n-r_{M}=r_{K}$: $\widehat{X}_{ij}=\widetilde{X}_{ij}$, $i=1,2$, $j=2,3$;

• •

$n-r_{M}>r_{K}$: $\widehat{X}_{12}=\widetilde{X}_{12}(:,1:r_{K})$, $\widehat{X}_{13}=(\widetilde{X}_{12}(:,r_{K}+1:n-r_{M}),\;\widetilde{X}_{13})$,

$\widehat{X}_{22}=\widetilde{X}_{22}(:,1:r_{K})$, $\widehat{X}_{23}=(\widetilde{X}_{22}(:,r_{K}+1:n-r_{M}),\;\widetilde{X}_{23})$

• •

$n-r_{M}<r_{K}$: $\widehat{X}_{12}=(\widetilde{X}_{12},\;\widetilde{X}_{13}(:,1:r_{K}+r_{M}-n))$, $\widehat{X}_{13}=\widetilde{X}_{13}(:,r_{K}+r_{M}-n+1:r_{K})$

$\widehat{X}_{22}=(\widetilde{X}_{22},\;\widetilde{X}_{23}(:,1:r_{K}+r_{M}-n))$, $\widehat{X}_{23}=\widetilde{X}_{23}(:,r_{K}+r_{M}-n+1:r_{K})$

On the right hand side, the transformation reads, analogously,

[TABLE]

• •

$n-r_{M}=r_{K}$: $\widehat{Y}_{ij}=\widetilde{Y}_{ij}$, $i=1,2$, $j=2,3$;

• •

$n-r_{M}>r_{K}$: $\widehat{Y}_{12}=\widetilde{Y}_{12}(:,1:r_{K})$, $\widehat{Y}_{13}=(\widetilde{Y}_{12}(:,r_{K}+1:n-r_{M}),\;\widetilde{Y}_{13})$,

$\widehat{Y}_{22}=\widetilde{Y}_{22}(:,1:r_{K})$, $\widehat{Y}_{23}=(\widetilde{Y}_{22}(:,r_{K}+1:n-r_{M}),\;\widetilde{Y}_{23})$

• •

$n-r_{M}<r_{K}$: $\widehat{Y}_{12}=(\widetilde{Y}_{12},\;\widetilde{Y}_{13}(:,1:r_{K}+r_{M}-n))$, $\widehat{Y}_{13}=\widetilde{Y}_{13}(:,r_{K}+r_{M}-n+1:r_{K})$

$\widehat{Y}_{22}=(\widetilde{Y}_{22},\;\widetilde{Y}_{23}(:,1:r_{K}+r_{M}-n))$, $\widehat{Y}_{23}=\widetilde{Y}_{23}(:,r_{K}+r_{M}-n+1:r_{K})$

[TABLE]

and it immediately reveals that the original quadratic pencil is singular if $\mathrm{det}(R_{X})=0$. Otherwise, we first identify the $n-r_{K}$ zero eigenvalues and the $n-r_{M}$ infinite ones, and the remaining ones are computed from the linear generalized eigenvalue problem of the $(r_{M}+r_{K})\times(r_{M}+r_{K})$ pencil

[TABLE]

3.2. Only one of the matrices is singular

In any case, we use the structure of the linearization to deflate two Jordan blocks of the eigenvalue zero, meaning that if $\scriptscriptstyle\left(\begin{smallmatrix}Q^{*}_{M,2}C\cr\hline\cr\widehat{R}_{M}P^{T}_{M}\end{smallmatrix}\right)$ is singular, i.e if there are at least two Jordan block for infinite eigenvalues, we switch to the reversed problem. Thus, we have the following: there are at least $n-\mathop{\rm rank}\nolimits(K)+n-\mathop{\rm rank}\nolimits(\left(\begin{smallmatrix}Q^{*}_{K,2}C\cr\hline\cr\widehat{R}_{K}P^{T}_{K}\end{smallmatrix}\right))=n-r_{K}+n-r_{22}$ zero eigenvalues and exactly $n-r_{M}$ infinite eigenvalues. First, two know Jordan blocks of zero eigenvalue are deflated as described in section §5. As a result we get the transformed pencil as in (5.30).

As in the case 2.2, the existence of the additional zero eigenvalues depends on the rank of the matrix $A_{33}$, and deflation of zeros is completed by calling the Algorithm 1 on the pencil $A_{33}-\lambda B_{33}$. The output arguments are the transformation matrices $Q_{p},P_{p}$ of order $r_{K}+r_{22}$ and pencil $A_{\ell+1,\ell+1}-\lambda B_{\ell+1,\ell+1}$ with $A_{\ell+1,\ell+1}$ being regular.

At this point, we have removed all zero eigenvalues. It remains to remove $n-r_{M}$ infinite ones. This is done by using Algorithm 1 on the reversed pencil $B_{\ell+1,\ell+1}-\lambda A_{\ell+1,\ell+1}$. Since we already determined that there is just one Jordan block for infinite eigenvalue, by computing the rank of $\left(\begin{smallmatrix}Q^{*}_{M,2}C\cr\hline\cr\widehat{R}_{M}P^{T}_{M}\end{smallmatrix}\right)$, we send, as input data of the Algorithm 1, the information about $\mathop{\rm rank}\nolimits(B_{\ell+1,\ell+1})=\mathop{\rm rank}\nolimits(\left(\begin{smallmatrix}Q^{*}_{M,2}C\cr\hline\cr\widehat{R}_{M}P^{T}_{M}\end{smallmatrix}\right))$ and the number of the deflation steps, which is one. Finally, as the output we get a pencil $A_{\ell+2,\ell+2}-\lambda B_{\ell+2,\ell+2}$ with both $A_{\ell+2,\ell+2}$ and $B_{\ell+2,\ell+2}$ regular, and the transformation matrices $Q_{P1},P_{P1}$.

Let $n_{\ell+1}$ be the dimension of $A_{\ell+1,\ell+1}\!-\!\lambda B_{\ell+1,\ell+1}$.The final transformations $Q,P$ are

[TABLE]

The remaining finite nonzero eigenvalues are computed using the QZ algorithm for the truncated pencil $A_{\ell+2,\ell+2}-\lambda B_{\ell+2,\ell+2}$.

3.3. Both matrices are singular

The decision whether we consider original or the reversed problem depends on the total number of the detected zero and infinite eigenvalues. The goal is to use the structure of the linearization pencil to deflate the zeros. Hence, we will consider the problem with the higher number of the zeros.

Therefore, at the beginning of the deflation process have the following setting: there are at least $n-\mathop{\rm rank}\nolimits(K)+n-\mathop{\rm rank}\nolimits(\left(\begin{smallmatrix}Q^{*}_{K,2}C\cr\hline\cr\widehat{R}_{K}P^{T}_{K}\end{smallmatrix}\right))=n-r_{K}+n-r_{22}$ zeros and $n-\mathop{\rm rank}\nolimits(M)+n-\mathop{\rm rank}\nolimits(\left(\begin{smallmatrix}Q^{*}_{M,2}C\cr\hline\cr\widehat{R}_{M}P^{T}_{M}\end{smallmatrix}\right))$ infinite eigenvalues. First step is to deflate all detected zeros using the procedure described in the §5. The resulting pencil is (5.30).

Possible extra zeros are deflated using the Algorithm 1 on the $A_{33}-\lambda B_{33}$. As the output we get the pencil $A_{\ell+1,\ell+1}-\lambda B_{\ell+1,\ell+1}$ with $A_{\ell+1,\ell+1}$ regular, and $Q_{p}$ and $P_{p}$ transformation matrices. Denote with $n_{\ell+1}$ the dimension of the resulted pencil.

At this point, all the zero eigenvalues are deflated. Now, we want to deflate infinite eigenvalues. We already know the dimension of first two Jordan blocks, namely $n-r_{M}$ and $n-\mathop{\rm rank}\nolimits(\left(\begin{smallmatrix}Q^{*}_{M,2}C\cr\hline\cr\widehat{R}_{M}P^{T}_{M}\end{smallmatrix}\right))$. This information is forwarded to the Algorithm 1 along with the reversed pencil $B_{\ell+1,\ell+1}-\lambda A_{\ell+1,\ell+1}$. In this case, the algorithm will not compute the numerical rank in the first two steps, but it will use the forwarded information. The output is the pencil $A_{\ell+\ell_{1},\ell+\ell_{1}}-\lambda B_{\ell+1,\ell+1}$, with both $A_{\ell+\ell_{1},\ell+\ell_{1}}$ and $B_{\ell+1,\ell+1}$ regular, and the transformation matrices $Q_{p1}$, $Q_{p_{2}}$. The final transformation matrices $Q$, $P$ are

[TABLE]

The remaining (finite nonzero) eigenvalues are computed using the QZ algorithm on $A_{\ell+\ell_{1},\ell+\ell_{1}}-\lambda B_{\ell+1,\ell+1}$.

6.2 Eigenvector recovery

The process of extracting the eigenvectors of the quadratic eigenvalue problem from those of the linearization is as follows (see also [14]). The right eigenvectors $z$ and the left eigenvectors $w$ of the linearization $C_{2}(\lambda)$ are of the form

[TABLE]

where $x,y$ are the right and the left eigenvector for the quadratic problem. From this relation we see that, when the matrix $K$ is nonsingular, we have two choices for the right eigenvector, namely $z_{1}$ and $K^{-1}z_{2}$. Otherwise, we only have one choice,111111Up to a nonzero scaling parameter. $z_{1}$.

The eigenvectors of $C_{2}(\lambda)$ are obtained from the eigenvectors $\widetilde{z}$, $\widetilde{w}$ of the transformed pencil $PC_{2}(\lambda)Q$ as $z=Q\widetilde{z}$, $w=P^{*}\widetilde{w}$ where $P$, $Q$ contain all transformations from the reduction process described in §6.1. For the readers convenience, we provide the details of assembling the eigenvectors.

6.2.1 The right eigenvectors

Case 1

The matrix $K$ is nonsingular. We have two choices for $x$. Let $\widetilde{z}$ be the right eigenvector for the transformed linearization. The corresponding right eigenvector for $C_{2}(\lambda)$ is z=\left(\begin{array}[]{c}z_{1}\cr\hline\cr\hline\cr\hline\cr z_{2}\end{array}\right)=\begin{pmatrix}\Pi_{M}&\mathbf{0}\\ \mathbf{0}&\mathbb{I}_{n}\end{pmatrix}\left(\begin{array}[]{c}\widetilde{z}_{1}\cr\hline\cr\hline\cr\hline\cr\widetilde{z}_{2}\end{array}\right)=\left(\begin{array}[]{c}\Pi_{M}\widetilde{z}_{1}\cr\hline\cr\hline\cr\hline\cr\widetilde{z}_{2}\end{array}\right). Hence, the two candidates for the eigenvector $x$ are $\Pi_{M}\widetilde{z}_{1}$ and $K^{-1}\widetilde{z}_{2}$. We choose the one with smaller backward error.

Case 2.1

The matrix $K$ is singular, and $n-r_{K}$ zero eigenvalues are deflated. Eigenvectors corresponding to those eigenvalues span the nullspace of matrix $K$. The basis for this nullspace is computed as orthogonal complement of the range of $K^{*}$ using the QR decomposition of the upper triangular matrix $\widehat{R}^{*}_{K}$, $\widehat{R}^{*}_{K}=Q_{\widehat{R}^{*}_{K}}R_{\widehat{R}^{*}_{K}}$. The wanted vector is determined by the last $n-r_{K}$ columns of the orthogonal matrix $Q_{\widehat{R}^{*}_{K}}$.

The remaining eigenvalues and eigenvectors $\widetilde{z}=\left(\begin{smallmatrix}\widetilde{z}_{1}\cr\widetilde{z}_{2}\end{smallmatrix}\right)\in\mathbb{C}^{n+r_{K}}$ are computed from the transformed linearization (6.20). The corresponding eigenvector for $C_{2}(\lambda)$ is

[TABLE]

Case 3.1

Both matrices $M$ and $K$ are singular, and $n-r_{M}$ infinite and $n-r_{K}$ zero eigenvalues are deflated. The eigenvectors for zero eigenvalues are obtained as in the case (2.1), whilst the eigenvectors for the infinite eigenvalues form a basis for the nullspace of the matrix $M$. The basis is obtained as the orthogonal complement of the range of $M^{*}$, represented by the last $n-r_{M}$ columns of the orthogonal matrix $Q_{\widehat{R}^{*}_{M}}$ in the QR factorization $\widehat{R}^{*}_{M}=Q_{\widehat{R}^{*}_{M}}R_{\widehat{R}^{*}_{M}}$. The remaining eigenvalues and eigenvectors $\widetilde{z}\!\in\!\mathbb{C}^{r_{K}+r_{M}}\!$ are obtained from the pencil (6.90). The corresponding eigenvector of $C_{2}(\lambda)$ is

[TABLE]

The eigenvector $x$ is represented by the first $n\!\times\!1$ block of the $(n\!+\!r_{K})\!\times\!1$ vector $Z^{*}_{X}\begin{pmatrix}\mathbf{0}\\ \widetilde{z}\end{pmatrix}$.

Cases 2.2., 3.2. and 3.3

When Algorithm 1 is used for the removal of additional zero and/or infinite eigenvalues, the eigenvector recovery procedure is as follows. Let $n_{\ell+j}$ be the dimension of the truncated pencil $A_{\ell+j,\ell+j}-\lambda B_{\ell+j,\ell+j}$, $j=1,2,\ell_{1}$. Let $\widetilde{z}\in\mathbb{R}^{n_{\ell+j}}$ be the computed right eigenvector of $A_{\ell+j,\ell+j}-\lambda B_{\ell+j,\ell+j}$. If $n_{\ell+j}>n$, the right eigenvector is recovered as $x=Q(1:n,1:n)\widetilde{z}(1:n)$, and if $n_{\ell+j}<n$ then $x=Q(1:n,1:n_{\ell+\ell_{1}})\widetilde{z}$.

6.2.2 The left eigenvectors

Case 1

Let $\widetilde{w}$ be the left eigenvector for the transformed linearization $PC_{2}(\lambda)Q$. The corresponding left eigenvector for the linearization $C_{2}(\lambda)$ is

[TABLE]

The two candidates for the left eigenvector $y$ of quadratic problem are $Q_{M}\widetilde{w}_{1}$ and $\widetilde{w}_{2}$. We choose the one with smaller backward error.

Case 2.1

The left eigenvectors for zero eigenvalues are represented by the last $n\!-\!r_{K}$ columns of the matrix $Q_{K}$. Let \left(\begin{array}[]{c}\widetilde{w}_{1}\cr\hline\cr\hline\cr\hline\cr\widetilde{w}_{2}\end{array}\right)\!\!\in\!\!\mathbb{C}^{n+r_{K}} be the eigenvector for the deflated pencil (6.20). The corresponding eigenvector for the $2n\times 2n$ pencil, before truncation satisfies

[TABLE]

The vector $z$ for $C_{2}(\lambda)$ is {\displaystyle\left(\begin{array}[]{c}w_{1}\cr\hline\cr\hline\cr\hline\cr w_{2}\end{array}\right)=\begin{pmatrix}Q_{M}&\mathbf{0}\\ \mathbf{0}&Q_{K}\end{pmatrix}\left(\begin{array}[]{c}\widetilde{w}_{1}\cr\hline\cr\hline\cr\hline\cr\widetilde{w}_{2}\cr\widetilde{w}_{3}\end{array}\right)=\left(\begin{array}[]{c}Q_{M}\widetilde{w}_{1}\cr\hline\cr\hline\cr\hline\cr Q_{K}\begin{pmatrix}\widetilde{w}_{2}\\ \widetilde{w}_{3}\end{pmatrix}\end{array}\right)}. The left eigenvector $y$ for the original problem is now chosen between $Q_{M}\widetilde{w}_{1}$ and $Q_{K}\begin{pmatrix}\widetilde{w}_{2}\\ \widetilde{w}_{3}\end{pmatrix}$. The one with smaller backward error is preferred.

Case 3.1

The left eigenvectors for zero eigenvalues are represented by the last $n-r_{K}$ columns of $Q_{K}$, and for infinite eigenvalues by the last $n-r_{M}$ columns of $Q_{M}$. Let $\begin{pmatrix}\widetilde{w}_{1}\\ \widetilde{w}_{2}\end{pmatrix}\in\mathbb{C}^{r_{K}+r_{M}}$ be a left eigenvector for the truncated $(r_{K}+r_{M})\times(r_{K}+r_{M})$ pencil (6.90). The corresponding eigenvector for the pencil $QC_{2}(\lambda)V$ satisfies

[TABLE]

The components $\widetilde{w}^{*}_{3},\widetilde{w}^{*}_{4}$ are thus computed as

[TABLE]

[TABLE]

and the candidates for the left eigenvector $y$ are $Q_{M}\left(\begin{smallmatrix}\widetilde{w}_{1}\\ Q_{X}\widetilde{w}_{3}\end{smallmatrix}\right)$ and $Q_{K}\left(\begin{smallmatrix}\widetilde{w}_{2}\\ \widetilde{w}_{4}\end{smallmatrix}\right)$. Again, we choose the one with smaller backward error.

Cases 2.2, 3.2 and 3.3

As for the right eigenvectors, in these cases we have the same procedure for the recovery of the left eigenvectors. Write the transformed linearization pencil as $P(A-\lambda B)Q=\begin{pmatrix}A_{\ell+j,\ell+j}-\lambda B_{\ell+j,\ell+j}&X\\ \mathbf{0}&Y\end{pmatrix},\;\;{j=1,2,\ell_{1}}$. Now, the left eigenvector $w\in\mathbb{R}^{2n}$ for the transformed pencil $P(A-\lambda B)Q$ is $w=\left(\begin{smallmatrix}w_{1}\\ w_{2}\end{smallmatrix}\right)$, $w_{2}=-w_{1}^{*}XY^{-1}$, where $w_{1}$ is the computed left eigenvector of $A_{\ell+j,\ell+j}-\lambda B_{\ell+j,\ell+j}$.

For the left eigenvector of the original quadratic problem we choose between $P^{*}w(1:n)$ and $P^{*}w(n+1:2n)$, whichever has lower backward error.

Finally, note that we can choose to measure the backward error norm-wise or as more structured as discussed before.

7 Numerical examples

The potential benefits of the proposed modifications of the quadeig framework have been already illustrated in Figure 1 in Example 1.1, and in Figure 2 in §4.2. We now present a summary of testing the proposed method on the NLEVP collection of benchmark examples. With respect to the norm-wise backward error (1.3), the comparison of quadeig and KVADeig shows no substantial difference, as expected – both use parameter scaling as a default. As noted in [14], polyeig is weaker than quadeig in this respect. Furthermore, we noted that, as shown on the left panel in Figure 2, quadeig performs better than polyeig also in terms of the structured backward error. A similar conclusion is obtained with the power_plant example, see Figure 3.

8 Concluding remarks

We have proposed a modification of the quadeig algorithm [14]. In addition to the parameter scaling, we use diagonal balancing and allow the rank revealing QR factorizations with full pivoting; this is shown to enhance numerical/backward stability as measured by structured backward error. Our backward error analysis is instructive – it shows the mechanism of keeping the backward error small in the coefficient matrices of the original problem in the first two steps of the deflation process (thus including quadeig), and it also shows how the proposed modifications contribute to keeping the structured backward error small. Further we extend the deflation procedure of quadeig toward an upper triangular variation of the Kronecker Canonical Form, which opens the possibility of removing additional infinite eigenvalues that were beyond the reach of quadeig. We believe that our contributions to the reduction framework introduced by quadeig will be useful for both the full solution of medium size problems and the iterative methods (e.g. SOAR) where full solution of the projected problem is needed and where notrivial techniques of deflation, purging and locking may be critically dependent on it. These issues are the subject of our ongoing and future work. A LAPACK–style implementation of the method, with further fine details of numerical software development is in progress.

Acknowledgments

This research was part of the second authors thesis [35] and of the project IP-11-2013-9345 (Mathematical Modelling, Analysis and Computing With Applications to Complex Mechanical Systems) supported by the Croatian Science Foundation. The authors thank Serkan Gugercin (Virginia Tech), Luka Grubišić and Zvonimir Bujanović (University of Zagreb) for valuable comments.