Computing s and ${\rm sep}$

To explain s and ${\rm sep}$, we need to introduce the spectral projector P [94,76], and the separation of two matrices A and B, ${\rm sep}(A,B)$ [94,98].

We may assume the matrix A is in Schur form, because reducing it to this form does not change the values of s and ${\rm sep}$. Consider a cluster of $m \geq 1$ eigenvalues, counting multiplicities. Further assume the n-by-n matrix A is

$$A = \left( \begin{array}{cc} A_{11} & A_{12} \\ 0 & A_{22} \end{array} \right)$$ (4.1)

where the eigenvalues of the m-by-m matrix A₁₁ are exactly those in which we are interested. In practice, if the eigenvalues on the diagonal of A are in the wrong order, routine xTREXC can be used to put the desired ones in the upper left corner as shown.

We define the spectral projector, or simply projector P belonging to the eigenvalues of A₁₁ as

$$P = \left( \begin{array}{cc} I_m & R \\ 0 & 0 \end{array} \right)$$ (4.2)

where R satisfies the system of linear equations

A₁₁R - RA₂₂ = A₁₂. (4.3)

Equation (4.3) is called a Sylvester equation. Given the Schur form (4.1), we solve equation (4.3) for R using the subroutine xTRSYL.

We can now define s for the eigenvalues of A₁₁:

$$s = \frac{1}{\Vert P\Vert _2} = \frac{1}{\sqrt{1+\Vert R\Vert _2^2}}.$$ (4.4)

In practice we do not use this expression since |R|₂ is hard to compute. Instead we use the more easily computed underestimate

$$\frac{1}{\sqrt{1+\Vert R\Vert _F^2}}$$ (4.5)

which can underestimate the true value of s by no more than a factor $\sqrt { \min ( m,n-m ) }$. This underestimation makes our error bounds more conservative. This approximation of s is called RCONDE in xGEEVX and xGEESX.

The separation ${\rm sep}(A_{11},A_{22})$ of the matrices A₁₁ and A₂₂ is defined as the smallest singular value of the linear map in (4.3) which takes X to A₁₁X - XA₂₂, i.e.,

$${\rm sep}(A_{11},A_{22}) = \min_{X \neq 0} \frac{\Vert A_{11}X - XA_{22}\Vert _F} {\Vert X \Vert _F}.$$ (4.6)

This formulation lets us estimate ${\rm sep}(A_{11},A_{22})$ using the condition estimator xLACON [59,62,63], which estimates the norm of a linear operator $\Vert T \Vert _1$ given the ability to compute Tx and T^Tx quickly for arbitrary x. In our case, multiplying an arbitrary vector by T means solving the Sylvester equation (4.3) with an arbitrary right hand side using xTRSYL, and multiplying by T^T means solving the same equation with A₁₁ replaced by A₁₁^T and A₂₂ replaced by A₂₂^T. Solving either equation costs at most O(n³) operations, or as few as O(n²) if $m \ll n$. Since the true value of ${\rm sep}$ is |T|₂ but we use |T|₁, our estimate of ${\rm sep}$ may differ from the true value by as much as $\sqrt{m(n-m)}$. This approximation to ${\rm sep}$ is called RCONDV by xGEEVX and xGEESX.

Another formulation which in principle permits an exact evaluation of ${\rm sep}(A_{11},A_{22})$ is

$${\rm sep}(A_{11},A_{22}) = \sigma_{\min} ( I_{n-m} \otimes A_{11} - A_{22}^T \otimes I_m )$$ (4.7)

where $X \otimes Y \equiv [ x_{ij} Y]$ is the Kronecker product of X and Y. This method is generally impractical, however, because the matrix whose smallest singular value we need is m(n-m) dimensional, which can be as large as n²/4. Thus we would require as much as O( n⁴ ) extra workspace and O(n⁶) operations, much more than the estimation method of the last paragraph.

The expression ${\rm sep}(A_{11},A_{22})$ measures the ``separation'' of the spectra of A₁₁ and A₂₂ in the following sense. It is zero if and only if A₁₁ and A₂₂ have a common eigenvalue, and small if there is a small perturbation of either one that makes them have a common eigenvalue. If A₁₁ and A₂₂ are both Hermitian matrices, then ${\rm sep}(A_{11},A_{22})$ is just the gap, or minimum distance between an eigenvalue of A₁₁ and an eigenvalue of A₂₂. On the other hand, if A₁₁ and A₂₂ are non-Hermitian, ${\rm sep}(A_{11},A_{22})$ may be much smaller than this gap.