Eigenvalue enclosures

Computing eigenvalues is a nonlinear problem, and ball arithmetic applied directly to an eigenvalue algorithm such as the QR iteration gives useless enclosures because the radii grow with every iteration. The functions of acb_mat.h therefore follow the approximate-and-certify approach of verification methods [Rum2010]: compute approximate eigenvalues \tilde\lambda_i and an approximate eigenvector matrix \tilde R by a floating-point algorithm without error control, and then prove, with a short computation in ball arithmetic, that the true eigenvalues lie within computed distances of the approximations. All statements below hold for every matrix in the balls of A.

Approximation#

acb_mat_approx_eig_qr reduces the midpoint matrix to upper Hessenberg form by Householder reflections and runs the implicitly shifted QR algorithm with deflation to obtain a Schur form A \approx Q T Q^* [GVL1996]; the diagonal of T gives the eigenvalues, and the eigenvectors are obtained by back substitution in T and multiplication by Q. The implementation is adapted from mpmath's eigen.py. An off-diagonal entry is treated as zero when it is smaller than the tolerance, by default 2^{-p}/(100n) relative to the neighbouring diagonal entries, and the iteration gives up if no deflation occurs within 14n + p/10 QR steps (by default); the function then returns 0, but still writes its current approximations, which may be usable. All arithmetic is done on midpoints at p bits.

A global enclosure#

This is the bound of Miyajima [Miy2010] in its simplest form. acb_mat_eig_global_enclosure computes Y as an approximate inverse of \tilde R and evaluates \delta and \varepsilon in ball arithmetic. The bound holds even if A has multiple eigenvalues or is not diagonalisable (if \tilde R is far from singular), but it does not say which disk contains which eigenvalue, or how many; it is used to obtain a quick inclusion region or as a starting point.

Simple eigenvalues#

If the eigenvalues are simple and well separated, the approximate diagonalisation can be certified for each eigenvalue individually. acb_mat_eig_simple_vdhoeven_mourrain (called by acb_mat_eig_simple) follows van der Hoeven and Mourrain [HM2017]:

The conditions on \eta in step 3 are those of the theorem of van der Hoeven and Mourrain [HM2017], which shows that there is an exact diagonalisation A = (T(I+E)) \Lambda (T(I+E))^{-1} with \|E\| \le 3\eta/\sigma; they are only needed for the eigenvectors. The cost is dominated by the ball matrix product AT and the preconditioned solve (Verified linear algebra). The enclosures have radius about \|H\|, which is of the order of the error of the approximation times the condition number of T; a slightly better enclosure could be obtained with the centre of the diagonal entries of D instead of the ball, as in [HM2017].

Clusters and multiple eigenvalues#

A multiple eigenvalue cannot be certified as simple, and an eigenvalue of a nontrivial Jordan block is perturbed by O(\epsilon^{1/k}) when A is perturbed by \epsilon, so its approximations are inaccurate. For such clusters, Rump's method [Rum2010, Section 13.4], also used in INTLAB's verifyeig, encloses a group of k eigenvalues together.

Let \tilde\lambda be an approximate eigenvalue of the cluster and \tilde X an n \times k matrix of approximations of a basis of the invariant subspace. Choose k row indices v (the k rows of \tilde X with the largest absolute value when k = 1, the last k rows otherwise); let U select the other rows. The unknowns are a correction Y (n \times k) whose rows in U correct \tilde X and whose rows in v form a k \times k matrix M with A X = X (\tilde\lambda I + M) for X = \tilde X + U U^T Y. With an approximate inverse R of (A - \tilde\lambda I) U U^T - \tilde X V^T, this equation becomes the fixed-point problem

acb_mat_eig_enclosure_rump implements this with \varepsilon-inflation: starting from \mathcal{Y}_0 = g(0), it enlarges \mathcal{Y} by 2^{-3+2j} \max|\mathcal{Y}| in iteration j and evaluates g(\mathcal{Y}) in ball arithmetic, for at most 5 + b iterations, where b is the bit length of p. The spectral radius in the eigenvalue bound is bounded by the Frobenius norm and by the Collatz–Wielandt bound \rho(N) \le \max_i (Nw)_i / w_i for positive vectors w, using the absolute values of approximate eigenvectors of N as w.

acb_mat_eig_multiple_rump groups the approximate eigenvalues into clusters of points within distance 2^{-3p/4} \max_i |\tilde\lambda_i| of the first point of the cluster (and on failure 2^{-p/2} \max_i |\tilde\lambda_i|), encloses each cluster of size k with acb_mat_eig_enclosure_rump, and succeeds if the enclosures are pairwise disjoint. Since each enclosure contains as many eigenvalues as its cluster has members and the cluster sizes add up to n, every eigenvalue of A is accounted for, and each enclosure is returned once per member. acb_mat_eig_multiple first tries the simple-eigenvalue algorithm, which is cheaper, and uses Rump's method only on failure.

As the theorem requires, acb_mat_eig_enclosure_rump accepts \mathcal{Y} only when every entry of g(\mathcal{Y}) lies in the interior of the corresponding entry of \mathcal{Y} (tested with acb_contains_interior); balls that touch the boundary are rejected.

Use in FLINT#

• acb_mat_approx_eig_qr (Approximation), acb_mat_eig_global_enclosure (Theorem 1), acb_mat_eig_simple, acb_mat_eig_simple_vdhoeven_mourrain, acb_mat_eig_simple_rump (Rump's method for each eigenvalue separately), acb_mat_eig_enclosure_rump, acb_mat_eig_multiple, acb_mat_eig_multiple_rump; see Eigenvalues and eigenvectors.
• There are no separate functions for real matrices; they are converted with acb_mat_set_arb_mat.
• gr_mat_diagonalization over gr_ctx_init_complex_acb calls acb_mat_approx_eig_qr and acb_mat_eig_simple (gr_mat.h).