hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

convergence sometimes fails for multiple eigenvalues; revert k>1 case, adjust tests and add notes

Browse source
This commit is contained in:
fredrik 2018-12-06 19:27:13 +01:00
parent 6d6049ede0
commit fe53e3f306
4 changed files with 141 additions and 7 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

28
acb_mat/eig_enclosure_rump.c
View file

@ -148,7 +148,7 @@ acb_approx_mag(mag_t res, const acb_t x)
/* Extract k largest rows to freeze */
static void
partition_X(slong * u, slong * v, const acb_mat_t X)
partition_X_sorted(slong * u, slong * v, const acb_mat_t X, slong prec)
{
slong i, j, n, k, c;
slong * row_idx;
@ -181,7 +181,9 @@ partition_X(slong * u, slong * v, const acb_mat_t X)
if (mag_cmp(row_mag + j, row_mag + j + 1) > 0)
{
mag_swap(row_mag + j, row_mag + j + 1);
c = row_idx[j]; row_idx[j] = row_idx[j + 1]; row_idx[j + 1] = c;
c = row_idx[j];
row_idx[j] = row_idx[j + 1];
row_idx[j + 1] = c;
}
}
}
@ -189,6 +191,7 @@ partition_X(slong * u, slong * v, const acb_mat_t X)
/* Not frozen rows of the approximation. */
for (i = 0; i < n - k; i++)
u[i] = row_idx[i];
/* Frozen rows of the approximation. */
for (i = 0; i < k; i++)
v[i] = row_idx[n - k + i];
@ -198,6 +201,22 @@ partition_X(slong * u, slong * v, const acb_mat_t X)
mag_clear(t);
}
static void
partition_X_trivial(slong * u, slong * v, const acb_mat_t X, slong prec)
{
slong n, k, i;
n = acb_mat_nrows(X);
k = acb_mat_ncols(X);
/* Not frozen rows of the approximation. */
for (i = 0; i < n - k; i++)
u[i] = i;
/* Frozen rows of the approximation. */
for (i = 0; i < k; i++)
v[i] = n - k + i;
}
void
acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t X, const acb_mat_t A,
const acb_t lambda_approx, const acb_mat_t X_approx, slong prec)
@ -221,7 +240,10 @@ acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t X, const acb_mat
/* Frozen rows of the approximation. */
v = flint_malloc(sizeof(slong) * k);
partition_X(u, v, X_approx);
if (k == 1)
partition_X_sorted(u, v, X_approx, prec);
else
partition_X_trivial(u, v, X_approx, prec);
mag_init(eps);
acb_mat_init(R, n, n);

89
acb_mat/test/t-eig_multiple.c
View file

@ -142,6 +142,95 @@ int main()
acb_clear(b);
}
/* Test convergence for DFT matrices */
for (iter = 0; iter < 50 * arb_test_multiplier(); iter++)
{
acb_mat_t A, R, QC;
acb_ptr E;
acb_t t;
fmpq_mat_t Q, Qinv;
slong i, n, c0, c1, c2, c3;
slong prec;
int algorithm, result;
n = n_randint(state, 30);
algorithm = n_randint(state, 2);
acb_mat_init(A, n, n);
acb_mat_init(R, n, n);
E = _acb_vec_init(n);
acb_init(t);
acb_mat_init(QC, n, n);
fmpq_mat_init(Q, n, n);
fmpq_mat_init(Qinv, n, n);
/* The current algorithm is not robust enough. */
#if 0
do {
fmpq_mat_randtest(Q, state, 2 + n_randint(state, 10));
} while (!fmpq_mat_inv(Qinv, Q));
#else
fmpq_mat_one(Q);
fmpq_mat_one(Qinv);
#endif
for (prec = 32; ; prec *= 2)
{
if (prec > 10000)
{
flint_printf("FAIL: unsuccessful, prec > 10000\n\n");
flint_printf("algorithm = %d, iter %wd\n\n", algorithm, iter);
flint_abort();
}
acb_mat_dft(A, 0, prec);
#if 0
acb_mat_set_fmpq_mat(QC, Q, prec);
acb_mat_mul(A, A, QC, prec);
acb_mat_set_fmpq_mat(QC, Qinv, prec);
acb_mat_mul(A, QC, A, prec);
#endif
acb_mat_approx_eig_qr(E, NULL, R, A, NULL, 0, prec);
if (algorithm == 0)
result = acb_mat_eig_multiple_rump(E, A, E, R, prec);
else
result = acb_mat_eig_multiple(E, A, E, R, prec);
/* Verify the known eigenvalues + multiplicities */
if (result)
{
c0 = c1 = c2 = c3 = 0;
for (i = 0; i < n; i++)
{
acb_set_d_d(t, 1.0, 0.0);
c0 += acb_contains(E + i, t);
acb_set_d_d(t, -1.0, 0.0);
c1 += acb_contains(E + i, t);
acb_set_d_d(t, 0.0, 1.0);
c2 += acb_contains(E + i, t);
acb_set_d_d(t, 0.0, -1.0);
c3 += acb_contains(E + i, t);
}
result = (n == 0 || (c0 == (n+4)/4 && c1 == (n+2)/4 && c2 == (n-1)/4 && c3 == (n+1)/4));
}
if (result)
break;
}
acb_mat_clear(A);
acb_mat_clear(R);
acb_mat_clear(QC);
_acb_vec_clear(E, n);
acb_clear(t);
fmpq_mat_clear(Q);
fmpq_mat_clear(Qinv);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");

20
acb_mat/test/t-eig_simple.c
View file

@ -156,9 +156,10 @@ int main()
/* Test convergence, given companion matrices */
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_mat_t A, R;
acb_mat_t A, R, QC;
acb_ptr E;
acb_ptr roots;
fmpq_mat_t Q, Qinv;
acb_poly_t f;
slong i, j, n, prec, count, count2;
int algorithm, success;
@ -170,6 +171,9 @@ int main()
acb_poly_init(f);
acb_mat_init(A, n, n);
acb_mat_init(R, n, n);
fmpq_mat_init(Q, n, n);
fmpq_mat_init(Qinv, n, n);
acb_mat_init(QC, n, n);
for (i = 0; i < n; i++)
{
@ -182,6 +186,10 @@ int main()
goto new_root;
}
do {
fmpq_mat_randtest(Q, state, 2 + n_randint(state, 100));
} while (!fmpq_mat_inv(Qinv, Q));
success = 0;
for (prec = 32; !success; prec *= 2)
@ -189,7 +197,7 @@ int main()
if (prec > 10000)
{
flint_printf("FAIL: unsuccessful, prec > 10000\n\n");
flint_printf("algorithm = %d\n\n", algorithm);
flint_printf("algorithm = %d, iter %wd\n\n", algorithm, iter);
flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n");
flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n");
flint_printf("roots = \n");
@ -202,7 +210,12 @@ int main()
}
acb_poly_product_roots(f, roots, n, prec);
acb_mat_companion(A, f, prec);
acb_mat_set_fmpq_mat(QC, Q, prec);
acb_mat_mul(A, A, QC, prec);
acb_mat_set_fmpq_mat(QC, Qinv, prec);
acb_mat_mul(A, QC, A, prec);
acb_mat_approx_eig_qr(E, NULL, R, A, NULL, 0, prec);
@ -250,6 +263,9 @@ int main()
}
}
fmpq_mat_clear(Q);
fmpq_mat_clear(Qinv);
acb_mat_clear(QC);
acb_mat_clear(A);
acb_mat_clear(R);
acb_poly_clear(f);

11
doc/source/acb_mat.rst
View file

@ -581,10 +581,17 @@ Component and error operations
Eigenvalues and eigenvectors
-------------------------------------------------------------------------------
The functions in this section are experimental. There may be classes
The functions in this section are experimental. There are classes
of matrices where the algorithms fail to converge even as
*prec* is increased, or for which the error bounds are much worse
than necessary. Manually balancing badly scaled matrices may help.
than necessary. In some cases, it can help to manually precondition
the matrix *A* by applying a similarity transformation `T^{-1} A T`.
* If *A* is badly scaled, take `T` to be a matrix such that the entries
of `T^{-1} A T` are more uniform (this is known as balancing).
* Simply taking `T` to be a random invertible matrix can help if an
algorithm fails to converge despite `A` being well-scaled. (This
can be the case when dealing with multiple eigenvalues.)
.. function:: 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)