Computing a canonical form of a matrix pencil

1.   Introduction

Let $A, B \in \mathbb{C}^{n\times n}$, with $A-\lambda B$ being a regular matrix pencil, i.e., $\det(A-\lambda B)$ is not identically zero. If the pencil is regular with respect to the unit circle, i.e., $\det(A- \lambda B) \neq 0$ for all $\lambda$ with $|\lambda| = 1$, then there exists a decomposition (see, e.g., [4, chap. 10]),

in which $T, Q \in \mathbb{C}^{n\times n}$ are nonsingular, and the eigenvalues of $A_1$ lie inside the open unit disk.

We are using $I_n$ to denote the identity matrix of order $n$. Throughout the paper, we will also be denoting the 2-norm of a given matrix $X$ as $\|X\|_2$, while $X^{*}$ will stand for its conjugate transpose.

The decomposition (1.1) is referred to as the canonical form of the pencil $A-\lambda B$. It characterizes that the set of eigenvalues of the pencil $A-\lambda B$ is the union of the set of eigenvalues inside the open unit disk (all of which are eigenvalues of $A_1$) and the set of eigenvalues outside the closed unit disk (all of which are reciprocals of the eigenvalues of $B_2$). The decomposition (1.1) is not unique. However, any other decomposition must have the form

with

where $Y_1$ and $Y_2$ are nonsingular matrices. Therefore, in practice, a canonical form can be found only up to such a transformation.

By using a Möbius transform, it is clear that the unit circle can be replaced by a more general region of the complex plane. In particular, the matrices $A_1$ and $B_2$ can be chosen in Jordan canonical form with a $B_2$ nilpotent (i.e., having only zero eigenvalues), and this yields the Weierstrass canonical form (see, e.g., [14, chap. 6]).

The spectral projector $\mathcal{P}_r$ (resp., $\mathcal{P}_l$) onto the right (resp. left) deflating subspace corresponding to the eigenvalues inside the unit circle can be expressed by the complex contour integral (see, e.g., [4, chap. 10], [6, chap. 4]):

These spectral projectors are connected by the relations $\mathcal{P}_l A = A \mathcal{P}_r$, $\mathcal{P}_l B = B \mathcal{P}_r.$ They contain all required information on the spectral properties of the pencil $A-\lambda B$. They may be useful, in part or as a whole, for example in model reduction (see, e.g., ^[2] and references therein), or when analyzing linear differential algebraic systems (see, e.g., ^[11]).

Using (1.1), the projectors (1.3) simplify to

Hence, if a decomposition (1.2) is known, the projectors are computable through (1.4). However, from a numerical algorithmic viewpoint, there are other ways to efficiently compute $\mathcal{P}_r$ and $\mathcal{P}_l$ without resorting to the canonical form. For example by applying a spectral dichotomy method or an inverse-free spectral divide and conquer method ^[1,10] or methods devised in the context of linear differential algebraic equations ^[11], although a canonical form –or part of it– may still be required in some circumstances.

While the knowledge of a canonical form (1.2) immediately implies the knowledge of the projectors through (1.4), it turns out that the converse question has not received much attention. This is the subject of our work: determine the canonical form knowing only the projector $\mathcal{P}_r$ or $\mathcal{P}_l$. In other words, solve the reverse problem of going from (1.4) to (1.2). In this case, the factor forms shown in the equalities of (1.4) are not known but are merely identities that presuppose (1.2). Specifically, the newness of the present note is thus: From the projector $\mathcal{P}_r$ or $\mathcal{P}_l$, we determine the matrices $T$, $Q$, $A_1$ and $B_2$ that characterize the representation (1.1) or the matrices in the equivalent form (1.2).

As noted earlier, choosing $B_2$ nilpotent results in the Weierstrass canonical form (also called Kronecker canonical form) that generalizes the Jordan canonical form meant for single matrices to regular matrix pairs. It can thus give insights into the structure of a generalized eigenproblem with the matrix pair $(A, B)$ and its underlying applications. As an illustrative example, consider the analysis of a linear system of differential algebraic equations (DAEs) with constant coefficients,

where $B$ is a singular matrix so that the above cannot simply be multiplied by $B^{-1}$ to revert to ordinary differential equations (ODEs). A canonical form can help transform such DAEs to ease their analysis (e.g., to explore different model parameters or different input vectors $f$).

The organization of the paper is as follows. After recalling some past works aimed at computing spectral projectors in Section 2, we give details about our new method in Section 3, followed by examples in Section 4 and finally some concluding remarks in Section 5.

2.   On determining the spectral projectors

Given that this study assumes that the spectral projectors are known, we cite here for ease of reference some common approaches used to obtain them in practice. Several spectral dichotomy algorithms credit their roots to the early works of Bulgakov, Godunov and Malyshev, of which ^[5,10] are examples among many others. In their extensive paper ^[1], Bai, Demmel, and Ming Gu stated that they built and improved on such early efforts to offer two inverse-free, highly parallel, spectral divide and conquer algorithms: one for computing an invariant subspace of a nonsymmetric matrix and another one for computing the left and right deflating subspaces of a regular matrix pencil. Bai et al. ^[1] included specifics on how to use a Möbius transform to split the dichotomy along a more general region of the complex plane than the unit circle. Aside from their practical usefulness, the algorithms in ^[1] were deemed more stable than other parallel algorithms for the nonsymmetric eigenproblem based on the matrix sign function.

In the course of investigating linear constant coefficient DAEs for properties such as local solvability or asymptotical stability, März ^[11] indicated the need to decouple those DAEs in canonical form by means of spectral projectors. The work included an iterative algorithm, dubbed matrix and projector chain, producing a sequence of matrices that become stationary and from which the desired projector emerges.

There are other methods ^[12] that also retrieve the projector through an iterative process. Finally, we can cite works such as ^[9,13] which first performed some form of decomposition where eigenvalues come into play (a generalized Schur factorization and block-diagonalization), allowing to ultimately obtain the spectral projections by either solving a system of generalized Sylvester equations or via a special reordering of the Schur factorization.

3.   Determining a canonical form using a spectral projector

3.1.   Main result

We now turn to the newness of our contribution, which is the determination of the canonical form knowing only the projector $\mathcal{P}_r$ or $\mathcal{P}_l$. We suppose that $\mathcal{P}_r$ is known (a similar reasoning can be applied to $\mathcal{P}_l$ and, therefore, is omitted), and denote by

the results of a column orthonormalization process that computes the matrices $U_1$ of size $n \times n_1$ and $V_2$ of size $n \times n_2$ whose columns form orthonormal bases of the range spaces $\mathsf{ran}(\mathcal{P}_r)$ and $\mathsf{ran}(I_n - \mathcal{P}_r)$. We have the following:

Theorem 3.1. Assume that the pencil $A-\lambda B$ is regular with respect to the unit circle, and consider the canonical form (1.1) with unknown matrices $T$, $Q$, $A_1$, and $B_2$. Set $\mathbb{Q} = \begin{bmatrix} U_1 & V_2 \end{bmatrix}^{-1}$ and $\mathbb{T} = \left(A+(B-A) \mathcal{P}_r\right)\mathbb{Q}^{-1}$, then the matrices $\mathbb{Q}$ and $\mathbb{T}$ are nonsingular and we have the canonical form

where the eigenvalues of $\mathbb{A}_{1}$ are inside the open unit disk and the reciprocals of the eigenvalues of $\mathbb{B}_{2}$ are outside the closed unit disk; the matrices $\mathbb{A}_{1}$, $A_1$ and $\mathbb{B}_{1}$, $B_2$ are related by

where $Y_{1}$ and $Y_{2}$ are nonsingular matrices.

Proof. The orthonormalization process needed in (3.1) can be done (see, for example, ^[7]) via a QR, SVD (singular value decomposition), or other method. Using, for instance, the SVD, we obtain

where $U_1$, $W_1$, ${V}_1$, and $Z_1$ are $n \times n_1$ matrices, $U_2$, $W_2$, $V_2$, $Z_2$ are $n \times n_2$ matrices, and $\begin{bmatrix} U_1 & U_2 \end{bmatrix}$, $\begin{bmatrix} W_1 & W_2 \end{bmatrix}$, $\begin{bmatrix} V_2 & V_1 \end{bmatrix}$ and $\begin{bmatrix} Z_2 & Z_1 \end{bmatrix}$ are unitary matrices. The matrix $\Sigma_1$ is an $n_1 \times n_1$ diagonal matrix containing the nonzero singular values of $\mathcal{P}_r$ and $\Gamma_2$ is an $n_2 \times n_2$ diagonal matrix containing the nonzero singular values of $I_n - \mathcal{P}_r$. Note that $n_1 = \mathsf{trace}(\mathcal{P}_r) = \mathsf{rank}(\mathcal{P}_r)$. The properties $\mathcal{P}_r^2 = \mathcal{P}_r$ and $(I_n-\mathcal{P}_r)^2 = I_n- \mathcal{P}_r$ imply that

From $\mathcal{P}_r U_1 = U_1$ and $(I_n- \mathcal{P}_r) V_2 = V_2$ we deduce that $\begin{bmatrix}I_{n_1} & 0 \\ 0 & 0 \end{bmatrix} QU_1 = QU_1$ and $\begin{bmatrix} 0 & 0 \\ 0 & I_{n_2} \end{bmatrix} QV_2 = QV_2$. Therefore $QU_1$ and $QV_2$ can be written as

where $Y_1$ and $Y_2$ are respectively $n_1 \times n_1$ and $n_2 \times n_2$ nonsingular matrices. As a consequence

This shows that the matrix $\mathbb{Q}$ is nonsingular and $\mathbb{Q}^{-1} = Q^{-1} \begin{bmatrix} Y_{1} & 0 \\ 0 & Y_{2} \end{bmatrix}.$ A simple calculation yields $A + (B-A) \mathcal{P}_r = TQ$, so $\mathbb{T}$ is nonsingular and $\mathbb{T}^{-1} = \mathbb{Q} (TQ)^{-1}$. Hence

□

3.2.   Perturbation results

Let the matrices $A$ and $B$ be slightly perturbed to $\tilde{A} = A + \Delta_{A}$ and $\tilde{B} = B + \Delta_{B}$, and let the spectral projection $\mathcal{P}_r$ be perturbed to ${\tilde{\mathcal{P}}}_r = \mathcal{P}_r + \Delta_{\mathcal{P}_r}$, accordingly. The projection ${\tilde{\mathcal{P}}}_r$ can be expressed as in (1.3) with $\tilde{A}$ and $\tilde{B}$ instead of $A$ and $B$. Our main objective is to study how far the matrices $\mathbb{T}$ and $\mathbb{Q}$ (in Theorem 3.1) will change when the matrices $A$ and $B$ are replaced by $\tilde{A}$ and $\tilde{B}$.

We assume that the perturbed pencil $\tilde{A} - \lambda \tilde{B}$ is regular and measure the size of perturbations by

Such a measure of perturbations appears when studying the $\epsilon$-pseudosepctrum of the pencil $\lambda B - A$ defined by (see, ^[15])

A related measure to $\epsilon$ is the norm of the resolvent $(\lambda B - A)^{-1}$ over the unit circle defined by

From (3.7) and (3.8), it is clear that $2\, r(A, B) \geq 1/\epsilon$ whenever $\lambda$ belongs to both the unit circle and the set $\Sigma_{\epsilon}(A, B)$. This entails certain restrictions on the existence of the perturbed projection, as the following proposition shows.

Proposition 1. Let the perturbations $\Delta_{A}$ and $\Delta_{B}$ be measured as in (3.6) and assume that $2 \, \epsilon\, r(A, B) < 1$, then the perturbed projection ${\tilde{P}}_r$ is well-defined and the perturbation $\Delta_{\mathcal{P}_r}$ satisfies $\|\Delta_{\mathcal{P}_r}\|_2 = \mathcal{O} \left(\epsilon\, r^2(A, B)\right).$

Proof. For $\lambda$ on the unit circle, we have

where

Since for $|\lambda| = 1$, we have

and we deduce that $C(\lambda)$ is nonsingular and

We also deduce that $\lambda \tilde{B} - \tilde{A}$ is nonsingular and, therefore, the projection ${\tilde{P}}_r$ is well-defined and is given by

A straightforward computation now gives

Therefore

where

We have

and using (3.6), (3.8), and (3.9), we deduce that

Since $1 + 2\, r(A, B) \|B\|_2 = \mathcal{O} \left(r(A, B)\right)$, we can write $\|\Delta_{\mathcal{P}_r}\|_2 = \mathcal{O} \left(\epsilon\, r^2(A, B)\right).$ □

Proposition 1 shows that the condition $2 \, \epsilon\, r(A, B) < 1$ ensures the existence of ${\tilde{P}}_r$, but in order to ensure the nearness of $\mathcal{P}_r$ and $\tilde{\mathcal{P}_r}$ and, hence, the smallness of $\|\Delta_{\mathcal{P}_r}\|_2$, we also need $\epsilon\, r^2(A, B)$ to be small. In particular, the resolvent should not be large on the unit circle. As we will see, the term $\epsilon\, r^2(A, B)$ essentially represents the error in the perturbed version of Theorem 3.1.

The following proposition gives an analogue of (3.1) for the perturbed projection.

Proposition 2. Let $U = \left(I_n + \Delta_{\mathcal{P}_r}\right)U_1$, $\tilde{U}_1 = U \left(U^*U\right)^{-\frac{1}{2}}$, $V = \left(I_n - \Delta_{\mathcal{P}_r}\right)V_1$, and $\tilde{V}_1 = V \left(V^*V\right)^{-\frac{1}{2}}.$ If $\epsilon$ is small enough so that Proposition 1 can be applied, then the columns of $\tilde{U}_1$ (resp., $\tilde{V}_1$) form an orthonormal basis of the range space $\mathsf{ran}(\mathcal{\tilde{P}}_r)$ (resp., $\mathsf{ran}(I_n - \mathcal{\tilde{P}}_r)$).

Proof. The assumption on $\epsilon$ ensures that Proposition 1 can be applied and, in particular, that $I_n \pm \Delta_{\mathcal{P}_r}$ is nonsingular. It follows that

From the properties of $\mathcal{P}_r$ (see (3.4) and (3.5)), we immediately deduce that

and these properties determine the null space of $\mathcal{\tilde{P}}_r$ and $I_n - \mathcal{\tilde{P}}_r$. □

Thus, if $A$ and $B$ are slightly perturbed to $\tilde{A} = A + \Delta_{A}$ and $\tilde{B} = B + \Delta_{B}$, as in Propositions 1 and 2, then the matrices $\mathbb{Q}$ and $\mathbb{T}$ in Theorem 3.1 can be replaced by

From Proposition 2, it is easy to check that $U = U_1 + \mathcal{O} \left(\epsilon\, r^2(A, B)\right)$, $U^*U = I_{n_1} + \mathcal{O} \left(\epsilon\, r^2(A, B)\right)$, and so $\tilde{U}_1 = U_1 + \mathcal{O} \left(\epsilon\, r^2(A, B)\right)$. Similarly, we have $\tilde{V}_1 = V_1 + \mathcal{O} \left(\epsilon\, r^2(A, B)\right)$. This leads to

4.   Examples

Two examples are given to illustrate Theorem 3.1. In the first example, the spectral projection is computed numerically by one of the methods referenced earlier, and in the second one, the spectral projection is obtained directly using the spectral properties of the pencil under consideration.

Example 1. Computations are carried out in MATLAB. Consider

The exact eigenvalues of the pencil $A - \lambda B$ are $-2$, $-0.5$, $0$, and $\infty$.

The spectal projection ${\mathcal P}_r$ onto the right deflating subspace corresponding to the eigenvalues inside the unit circle is computed using the algorithm developed in ^[1] to which the reader is referred for further details on convergence and stability. Below is a brief description.

Start with the matrices $A$ and $B$ as initializers for $A_0$ and $B_0$, then proceed iteratively: At each iteration $k = 0, 1, ...$, perform a $QR$ decomposition of the matrix $\begin{bmatrix} B_k \\-A_k \end{bmatrix}$, and replace the matrices $A_k$ and $B_k$ by, respectively, $Q_{12}^*A_k$ and $Q_{22}^*B_k$, where $Q_{12}$ (resp., $Q_{22}$) denotes the $n \times n$ submatrix of $Q$ ---that results from the QR decomposition--- formed from the first $n$ rows and last $n$ columns (resp., last $n$ rows and last $n$ columns). After a few iterations (see [1, Section 9] for the stopping criterion), a good approximation of the projection $\mathcal{P}_r$ is obtained from $(A_k+B_k)^{-1}B_k$. With the matrices of this example, after $10$ iterations, we obtain

Then Theorem 3.1 is used to compute the matrices $\mathbb{T}$, $\mathbb{Q}$, $\mathbb{A}_1$ and $\mathbb{B}_2$.

The 2-norm of the absolute error between $A$, $B$ and the computed canonical form is given by

and the eigenvalues of $\mathbb{A}_1$ are $-0.5$ and $6.9389\cdot 10^{-18}$ and those of $\mathbb{B}_2$ are $-0.5$ and $0$ (having reciprocals $-2$ and $\infty$). Hence, the exact values agree with the numerical predictions of Theorem 3.1.

Example 2. Consider the matrix pencil $A-\lambda B$ with

where $K$ is an $n \times n$ matrix and $C$ is an $n \times m$ matrix of full rank with $n > m$. Such a pencil arises in the stability analysis of the steady-state solutions of equations modeling viscous incompressible flows (see, e.g., ^[3]). The case when $B$ has the form $\begin{bmatrix} M & 0 \\ 0 & 0 \end{bmatrix}$, where $M$ is a Hermitian positive definite matrix, can be brought back to the matrix $B$ in (4.1) by changing $K$ to $L^{-1}K(L^{-1})^*$ and $C$ to $L^{-1}C$, where $L$ is the Choleski factor of $M$.

Consider the $QR$ decomposition of $C$:

where $R_1$ is $m\times m$ nonsingular and upper triagular, $Q$ is $n \times n$ unitary, $Q_1$ is $n \times m$, and $Q_2$ is $n \times (n-m)$. From [3, Theorem 2.1], the pencil $A-\lambda B$ has $n-m$ finite eigenvalues and $2m$ infinite eigenvalues of algebraic multiplicity $2m$. The finite eigenvalues are given by the eigenvalues of the matrix $Q_2^*KQ_2$ and the corresponding eigenvectors are spanned by the range of $\begin{bmatrix} Q_2Q_2^* \\ -R_{1}^{-1}KQ_2 Q_2^* \end{bmatrix}.$

Let us find the spectral projector $\mathcal{P}_r$ corresponding to all finite eigenvalues (instead of the unit circle, we consider here any circle that encloses all finite eigenvalues of the pencil) and the resulting canonical form.

Using the properties ${\mathcal{P}}_r^2 = \mathcal{P}_r$, $B \mathcal{P}_r = \mathcal{P}_l B$, $\mathsf{ran} (\mathcal{P}_r) = \mathsf{ran} \left(\begin{bmatrix} Q_2Q_2^* \\ -R_{1}^{-1}KQ_2 Q_2^* \end{bmatrix}\right)$, and $\mathsf{trace}(\mathcal{P}_r) = \mathsf{rank}(\mathcal{P}_r) = n-m$, we deduce that

and in low rank forms

where $K_{ij} = Q_i^*KQ_j, i, j = 1, 2$. Note that $\begin{bmatrix} Q_2 \\ -R_{1}^{-1}K_{12}\end{bmatrix}$ and $\begin{bmatrix} 0 & Q_1 \\ R_{1}^{-1} & - R_{1}^{-1} (K_{11} - I_m) \end{bmatrix}$ have full rank and therefore their columns provide bases for the range spaces of $\mathcal{P}_r$ and $I_{n+m}-\mathcal{P}_r$.

Consider the $QR$ decompositions:

where $U_1$ and $V_2$ are matrices whose columns provide orthonormal bases for $\mathsf{ran}(\mathcal{P}_r)$ and $\mathsf{ran}(I_n - \mathcal{P}_r)$, and $Y_1$ and $Y_2$ are, respectively, $n \times n$ and $2m \times 2m$ (upper triangular) nonsingular matrices. Using (4.3), the matrices $\mathbb{Q}$ and $\mathbb{T}$ in Therorem 3.1 are given by

and

A straightforward computation gives the canonical form stated in Theorem 3.1 with

thus providing analytic solutions to a class of matrix pencil problems that have the form (4.1).

5.   Conclusions

It is often convenient in numerical computing to be able to retrieve useful quantities as byproducts from others, as is the case with estimating the residual or condition number when solving a linear system, or obtaining optional factors upon performing a matrix factorization. Here we have described a new method with which to determine the canonical form of a pencil by taking advantage of existing efficient numerical algorithms already developed to determine spectral projections. We performed a perturbation analysis of the proposed method and illustrated it with two examples. The first example was a walk-through that numerically went over the steps of the method, while the second example showed how to further build on the method to theoretically establish an analytic representation of the canonical form of a pencil that arises in the stability analysis of the steady-state solutions of equations modeling viscous incompressible flows. Given that much efforts have been devoted to determining the spectral projections compared to determining the canonical form, our proposed method is a potent bridge from the former to the latter, and it also helps fill a gap in the literature. Application areas where this newfound ability may be useful include model order reduction (see, e.g., ^[2,8]).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors have declared no conflicts of interest.