add acb_mat_randtest_eig, acb_mat_eig_global_enclosure

acb_mat.h

@@ -94,6 +94,8 @@ void acb_mat_set_round_arb_mat(acb_mat_t dest, const arb_mat_t src, slong prec);
 
 void acb_mat_randtest(acb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits);
 
+void acb_mat_randtest_eig(acb_mat_t A, flint_rand_t state, acb_srcptr E, slong prec);
+
 /* I/O */
 
 void acb_mat_fprintd(FILE * file, const acb_mat_t mat, slong digits);
@@ -406,6 +408,8 @@ void acb_mat_det(acb_t det, const acb_mat_t A, slong prec);
 
 int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec);
 
+void acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec);
+
 /* Special functions */
 
 void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec);

acb_mat/eig_global_enclosure.c

@@ -0,0 +1,63 @@
+/*
+    Copyright 2018 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_mat.h"
+
+void
+acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec)
+{
+    acb_mat_t Y, R1, R2;
+    slong i, j, n;
+    mag_t r1, r2;
+
+    n = acb_mat_nrows(A);
+
+    acb_mat_init(Y, n, n);
+    acb_mat_init(R1, n, n);
+    acb_mat_init(R2, n, n);
+    mag_init(r1);
+    mag_init(r2);
+
+    /* Y ~= inv(R) */
+    acb_mat_one(R1);
+    acb_mat_approx_solve(Y, R, R1, prec);
+
+    /* R2 = Y*R - I */
+    acb_mat_mul(R2, Y, R, prec);
+    for (i = 0; i < n; i++)
+        acb_sub_ui(acb_mat_entry(R2, i, i), acb_mat_entry(R2, i, i), 1, prec);
+
+    acb_mat_bound_inf_norm(r2, R2);
+
+    if (mag_cmp_2exp_si(r2, 0) < 0)
+    {
+        /* R1 = Y*(AR - RD) */
+        acb_mat_mul(R2, A, R, prec);
+        for (i = 0; i < n; i++)
+            for (j = 0; j < n; j++)
+                acb_submul(acb_mat_entry(R2, i, j), acb_mat_entry(R, i, j), E + j, prec);
+        acb_mat_mul(R1, Y, R2, prec);
+
+        acb_mat_bound_inf_norm(r1, R1);
+        mag_geom_series(r2, r2, 0);
+        mag_mul(eps, r1, r2);
+    }
+    else
+    {
+        mag_inf(eps);
+    }
+
+    acb_mat_clear(R1);
+    acb_mat_clear(R2);
+    acb_mat_clear(Y);
+    mag_clear(r1);
+    mag_clear(r2);
+}

acb_mat/randtest_eig.c

@@ -0,0 +1,58 @@
+/*
+    Copyright (C) 2018 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_mat.h"
+
+void
+acb_mat_randtest_eig(acb_mat_t A, flint_rand_t state, acb_srcptr E, slong prec)
+{
+    slong n, i, j, ebits;
+    acb_mat_t U, Q;
+
+    n = acb_mat_nrows(A);
+    ebits = 1 + n_randint(state, 5);
+
+    acb_mat_init(U, n, n);
+    acb_mat_init(Q, n, n);
+
+    /* Skew-Hermitian matrix */
+    acb_mat_randtest(Q, state, prec, 1);
+    if (n_randint(state, 2))
+        acb_mat_get_mid(Q, Q);
+    for (i = 0; i < n; i++)
+    {
+        for (j = i + 1; j < n; j++)
+        {
+            acb_neg(acb_mat_entry(Q, i, j), acb_mat_entry(Q, j, i));
+            acb_conj(acb_mat_entry(Q, i, j), acb_mat_entry(Q, i, j));
+        }
+        arb_zero(acb_realref(acb_mat_entry(Q, i, i)));
+    }
+
+    acb_mat_exp(Q, Q, prec);
+
+    acb_mat_randtest(U, state, prec, ebits);
+    if (n_randint(state, 2))
+        acb_mat_get_mid(U, U);
+    for (i = 0; i < n; i++)
+        for (j = 0; j < i; j++)
+            acb_zero(acb_mat_entry(U, i, j));
+
+    for (i = 0; i < n; i++)
+        acb_set(acb_mat_entry(U, i, i), E + i);
+
+    acb_mat_mul(U, Q, U, prec);
+    acb_mat_conjugate_transpose(Q, Q);
+    acb_mat_mul(A, U, Q, prec);
+
+    acb_mat_clear(U);
+    acb_mat_clear(Q);
+}

acb_mat/test/t-eig_global_enclosure.c

@@ -0,0 +1,107 @@
+/*
+    Copyright (C) 2018 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_mat.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("eig_global_enclosure....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 250 * arb_test_multiplier(); iter++)
+    {
+        acb_mat_t A, X;
+        acb_ptr E, F, B;
+        mag_t eps;
+        slong i, j, n, prec;
+        int success;
+
+        n = n_randint(state, 8);
+        prec = 2 + n_randint(state, 100);
+
+        acb_mat_init(A, n, n);
+        acb_mat_init(X, n, n);
+        E = _acb_vec_init(n);
+        F = _acb_vec_init(n);
+        B = _acb_vec_init(n);
+        mag_init(eps);
+
+        acb_mat_randtest_eig(A, state, E, prec);
+        acb_mat_approx_eig_qr(F, NULL, X, A, NULL, 0, prec);
+
+        /* perturb F further */
+        if (n_randint(state, 2))
+        {
+            for (i = 0; i < n; i++)
+            {
+                acb_randtest(B, state, prec, 1);
+                acb_mul_2exp_si(B, B, -n_randint(state, prec));
+                acb_add(F + i, F + i, B, prec);
+            }
+        }
+
+        acb_mat_eig_global_enclosure(eps, A, F, X, prec);
+
+        _acb_vec_set(B, F, n);
+        for (i = 0; i < n; i++)
+            acb_add_error_mag(B + i, eps);
+
+        for (i = 0; i < n; i++)
+        {
+            success = 0;
+
+            for (j = 0; j < n; j++)
+            {
+                if (acb_contains(B + j, E + i))
+                {
+                    success = 1;
+                    break;
+                }
+            }
+
+            if (!success)
+            {
+                flint_printf("FAIL\n\n");
+                flint_printf("A =\n");
+                acb_mat_printd(A, 20);
+                flint_printf("\n\ni = %wd\n\n", i);
+                flint_printf("\nE =\n");
+                for (j = 0; j < n; j++)
+                {
+                    acb_printn(E + j, 20, 0); flint_printf("\n");
+                }
+                flint_printf("\nB =\n");
+                for (j = 0; j < n; j++)
+                {
+                    acb_printn(B + j, 20, 0); flint_printf("\n");
+                }
+                flint_abort();
+            }
+        }
+
+        acb_mat_clear(A);
+        acb_mat_clear(X);
+        _acb_vec_clear(E, n);
+        _acb_vec_clear(F, n);
+        _acb_vec_clear(B, n);
+        mag_clear(eps);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/acb_mat.rst

@@ -104,6 +104,13 @@ Random generation
     Sets *mat* to a random matrix with up to *prec* bits of precision
     and with exponents of width up to *mag_bits*.
 
+.. function:: void acb_mat_randtest_eig(acb_mat_t mat, flint_rand_t state, acb_srcptr E, slong prec)
+
+    Sets *mat* to a random matrix with the prescribed eigenvalues
+    supplied as the vector *E*. The output matrix is required to be
+    square. We generate a random unitary matrix via a matrix
+    exponential, and then evaluate an inverse Schur decomposition.
+
 Input and output
 -------------------------------------------------------------------------------
 
@@ -545,3 +552,46 @@ Eigenvalues and eigenvectors
     return value indicates that the QR iteration converged numerically,
     but this is only a heuristic termination test and does not imply
     any statement whatsoever about error bounds.
+
+.. function:: void acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec)
+
+    Given an *n* by *n* matrix *A*, a length-*n* vector *E*
+    containing approximations of the eigenvalues of *A*,
+    and an *n* by *n* matrix *R* containing approximations of
+    the corresponding eigenvectors, computes a rigorous bound
+    `\varepsilon` such that every eigenvalue `\lambda` of *A* satisfies
+    `|\lambda - \hat \lambda_k| \le \varepsilon`
+    for some `\hat \lambda_k` in *E*.
+
+    In other words, the union of the balls
+    `B_k = \{z : |z - \hat \lambda_k| \le \varepsilon\}` is guaranteed to
+    be an enclosure of all eigenvalues of *A*.
+
+    Note that there is no guarantee that each ball `B_k` can be
+    identified with a single eigenvalue: it is possible that some
+    balls contain several eigenvalues while other balls contain
+    no eigenvalues. In other words, this method is not powerful enough
+    to compute isolating balls for the individual eigenvalues (or even
+    for clusters of eigenvalues other than the whole spectrum).
+    Nevertheless, in practice the balls `B_k` will represent
+    eigenvalues one-to-one with high probability if the
+    given approximations are good.
+
+    The output can be used to certify
+    that all eigenvalues of *A* lie in some region of the complex
+    plane (such as a specific half-plane, strip, disk, or annulus)
+    without the need to certify the individual eigenvalues.
+    The output is easily converted into lower or upper bounds
+    for the absolute values or real or complex parts
+    of the spectrum, and with high probability these bounds will be tight.
+    Using :func:`acb_add_error_mag` and :func:`acb_union`, the output
+    can also be converted to a single :type:`acb_t` enclosing
+    the whole spectrum of *A* in a rectangle, but note that to
+    test whether a condition holds for all eigenvalues of *A*, it
+    is typically better to iterate over the individual balls `B_k`.
+
+    This function implements the fast algorithm in Theorem 1 in
+    [Miy2010]_ which extends the Bauer-Fike theorem. Approximations
+    can, for instance, be computed using :func:`acb_mat_approx_eig_qr`.
+    No assumptions are made about the structure of *A* or the
+    quality of the given approximations.

doc/source/credits.rst

@@ -225,6 +225,8 @@ Bibliography
 
 .. [Kri2013] \A. Krishnamoorthy and D. Menon, "Matrix Inversion Using Cholesky Decomposition" Proc. of the International Conference on Signal Processing Algorithms, Architectures, Arrangements, and Applications (SPA-2013), pp. 70-72, 2013.
 
+.. [Miy2010] \S. Miyajima, "Fast enclosure for all eigenvalues in generalized eigenvalue problems", Journal of Computational and Applied Mathematics, 233 (2010), 2994-3004, https://dx.doi.org/10.1016/j.cam.2009.11.048
+
 .. [MPFR2012] The MPFR team, "MPFR Algorithms" (2012), http://www.mpfr.org/algo.html
 
 .. [NIST2012] National Institute of Standards and Technology, *Digital Library of Mathematical Functions* (2012), http://dlmf.nist.gov/
@@ -246,4 +248,3 @@ Bibliography
 .. [Tak2000] \D. Takahashi, "A fast algorithm for computing large Fibonacci numbers", Information Processing Letters 75 (2000) 243-246, http://www.ii.uni.wroc.pl/~lorys/IPL/article75-6-1.pdf
 
 .. [Tre2008] \L. N. Trefethen, "Is Gauss Quadrature Better than Clenshaw-Curtis?", SIAM Review, 50:1 (2008), 67-87, https://doi.org/10.1137/060659831
-