From e640056d7fed74655a47c74024b3260904baeecf Mon Sep 17 00:00:00 2001 From: fredrik Date: Sat, 24 Nov 2018 23:41:05 +0100 Subject: [PATCH] add acb_mat_randtest_eig, acb_mat_eig_global_enclosure --- acb_mat.h | 4 + acb_mat/eig_global_enclosure.c | 63 +++++++++++++++ acb_mat/randtest_eig.c | 58 ++++++++++++++ acb_mat/test/t-eig_global_enclosure.c | 107 ++++++++++++++++++++++++++ doc/source/acb_mat.rst | 50 ++++++++++++ doc/source/credits.rst | 3 +- 6 files changed, 284 insertions(+), 1 deletion(-) create mode 100644 acb_mat/eig_global_enclosure.c create mode 100644 acb_mat/randtest_eig.c create mode 100644 acb_mat/test/t-eig_global_enclosure.c diff --git a/acb_mat.h b/acb_mat.h index dce02deb..af3c6db0 100644 --- a/acb_mat.h +++ b/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); diff --git a/acb_mat/eig_global_enclosure.c b/acb_mat/eig_global_enclosure.c new file mode 100644 index 00000000..c61c8308 --- /dev/null +++ b/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 . +*/ + +#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); +} diff --git a/acb_mat/randtest_eig.c b/acb_mat/randtest_eig.c new file mode 100644 index 00000000..c87d1ab5 --- /dev/null +++ b/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 . +*/ + +#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); +} diff --git a/acb_mat/test/t-eig_global_enclosure.c b/acb_mat/test/t-eig_global_enclosure.c new file mode 100644 index 00000000..dfe85007 --- /dev/null +++ b/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 . +*/ + +#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; +} diff --git a/doc/source/acb_mat.rst b/doc/source/acb_mat.rst index 4108e1a3..9db75a21 100644 --- a/doc/source/acb_mat.rst +++ b/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. diff --git a/doc/source/credits.rst b/doc/source/credits.rst index 11de7fa3..52a4d995 100644 --- a/doc/source/credits.rst +++ b/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 -