mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
Merge pull request #127 from argriffing/interval-frobenius-norm
add interval frobenius norm
This commit is contained in:
Fredrik Johansson
commit
508583b0f8
12 changed files with 640 additions and 139 deletions
|
|
@ -150,7 +150,9 @@ void acb_mat_transpose(acb_mat_t mat1, const acb_mat_t mat2);
|
|||
|
||||
void acb_mat_bound_inf_norm(mag_t b, const acb_mat_t A);
|
||||
|
||||
void acb_mat_bound_fro_norm(mag_t b, const acb_mat_t A);
|
||||
void acb_mat_frobenius_norm(arb_t res, const acb_mat_t A, slong prec);
|
||||
|
||||
void acb_mat_bound_frobenius_norm(mag_t b, const acb_mat_t A);
|
||||
|
||||
/* Arithmetic */
|
||||
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@
|
|||
#include "acb_mat.h"
|
||||
|
||||
void
|
||||
acb_mat_bound_fro_norm(mag_t b, const acb_mat_t A)
|
||||
acb_mat_bound_frobenius_norm(mag_t b, const acb_mat_t A)
|
||||
{
|
||||
slong i, j, r, c;
|
||||
mag_t t;
|
||||
|
|
@ -45,7 +45,7 @@ acb_mat_bound_fro_norm(mag_t b, const acb_mat_t A)
|
|||
{
|
||||
for (j = 0; j < c; j++)
|
||||
{
|
||||
acb_ptr z = acb_mat_entry(A, i, j);
|
||||
acb_srcptr z = acb_mat_entry(A, i, j);
|
||||
arb_get_mag(t, acb_realref(z));
|
||||
mag_addmul(b, t, t);
|
||||
arb_get_mag(t, acb_imagref(z));
|
||||
52
acb_mat/frobenius_norm.c
Normal file
52
acb_mat/frobenius_norm.c
Normal file
|
|
@ -0,0 +1,52 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2016 Arb authors
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb_mat.h"
|
||||
|
||||
void
|
||||
acb_mat_frobenius_norm(arb_t res, const acb_mat_t A, slong prec)
|
||||
{
|
||||
slong i, j, r, c;
|
||||
|
||||
r = acb_mat_nrows(A);
|
||||
c = acb_mat_ncols(A);
|
||||
|
||||
arb_zero(res);
|
||||
|
||||
if (r == 0 || c == 0)
|
||||
return;
|
||||
|
||||
for (i = 0; i < r; i++)
|
||||
{
|
||||
for (j = 0; j < c; j++)
|
||||
{
|
||||
acb_srcptr z = acb_mat_entry(A, i, j);
|
||||
arb_addmul(res, acb_realref(z), acb_realref(z), prec);
|
||||
arb_addmul(res, acb_imagref(z), acb_imagref(z), prec);
|
||||
}
|
||||
}
|
||||
|
||||
arb_sqrtpos(res, res, prec);
|
||||
}
|
||||
278
acb_mat/test/t-frobenius_norm.c
Normal file
278
acb_mat/test/t-frobenius_norm.c
Normal file
|
|
@ -0,0 +1,278 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2016 Arb authors
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb_mat.h"
|
||||
|
||||
static void
|
||||
_acb_mat_set_fmpq_mat_fmpq_mat(acb_mat_t C,
|
||||
const fmpq_mat_t rsrc, const fmpq_mat_t isrc, slong prec)
|
||||
{
|
||||
slong i, j;
|
||||
for (i = 0; i < acb_mat_nrows(C); i++)
|
||||
{
|
||||
for (j = 0; j < acb_mat_ncols(C); j++)
|
||||
{
|
||||
arb_set_fmpq(acb_realref(acb_mat_entry(C, i, j)),
|
||||
fmpq_mat_entry(rsrc, i, j), prec);
|
||||
arb_set_fmpq(acb_imagref(acb_mat_entry(C, i, j)),
|
||||
fmpq_mat_entry(isrc, i, j), prec);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static void
|
||||
_acb_mat_conjugate_transpose(acb_mat_t B, const acb_mat_t A)
|
||||
{
|
||||
slong i, j;
|
||||
acb_mat_transpose(B, A);
|
||||
for (i = 0; i < acb_mat_nrows(B); i++)
|
||||
for (j = 0; j < acb_mat_ncols(B); j++)
|
||||
acb_conj(acb_mat_entry(B, i, j), acb_mat_entry(B, i, j));
|
||||
}
|
||||
|
||||
static void
|
||||
_fmpq_mat_sum_of_squares(fmpq_t res, const fmpq_mat_t Q)
|
||||
{
|
||||
slong i, j;
|
||||
fmpq_zero(res);
|
||||
for (i = 0; i < fmpq_mat_nrows(Q); i++)
|
||||
{
|
||||
for (j = 0; j < fmpq_mat_ncols(Q); j++)
|
||||
{
|
||||
fmpq_addmul(res, fmpq_mat_entry(Q, i, j), fmpq_mat_entry(Q, i, j));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("frobenius_norm....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
/* compare to the exact rational norm */
|
||||
for (iter = 0; iter < 10000; iter++)
|
||||
{
|
||||
fmpq_mat_t Qr, Qi;
|
||||
fmpq_t q;
|
||||
acb_mat_t A;
|
||||
slong n, qbits, prec;
|
||||
|
||||
n = n_randint(state, 8);
|
||||
qbits = 1 + n_randint(state, 100);
|
||||
prec = 2 + n_randint(state, 200);
|
||||
|
||||
fmpq_mat_init(Qr, n, n);
|
||||
fmpq_mat_init(Qi, n, n);
|
||||
fmpq_init(q);
|
||||
|
||||
acb_mat_init(A, n, n);
|
||||
|
||||
fmpq_mat_randtest(Qr, state, qbits);
|
||||
fmpq_mat_randtest(Qi, state, qbits);
|
||||
|
||||
/* compute the square of the exact rational norm */
|
||||
{
|
||||
fmpq_t qr, qi;
|
||||
fmpq_init(qr);
|
||||
fmpq_init(qi);
|
||||
_fmpq_mat_sum_of_squares(qr, Qr);
|
||||
_fmpq_mat_sum_of_squares(qi, Qi);
|
||||
fmpq_add(q, qr, qi);
|
||||
fmpq_clear(qr);
|
||||
fmpq_clear(qi);
|
||||
}
|
||||
|
||||
_acb_mat_set_fmpq_mat_fmpq_mat(A, Qr, Qi, prec);
|
||||
|
||||
/* check that the arb interval contains the exact value */
|
||||
{
|
||||
arb_t a;
|
||||
arb_init(a);
|
||||
|
||||
acb_mat_frobenius_norm(a, A, prec);
|
||||
arb_mul(a, a, a, prec);
|
||||
|
||||
if (!arb_contains_fmpq(a, q))
|
||||
{
|
||||
flint_printf("FAIL (containment, iter = %wd)\n", iter);
|
||||
flint_printf("n = %wd, prec = %wd\n", n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("Qr = \n"); fmpq_mat_print(Qr); flint_printf("\n\n");
|
||||
flint_printf("Qi = \n"); fmpq_mat_print(Qi); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(Q)^2 = \n");
|
||||
fmpq_print(q); flint_printf("\n\n");
|
||||
|
||||
flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(A)^2 = \n");
|
||||
arb_printd(a, 15); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(A)^2 = \n");
|
||||
arb_print(a); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_clear(a);
|
||||
}
|
||||
|
||||
/* check that the upper bound is not less than the exact value */
|
||||
{
|
||||
mag_t b;
|
||||
fmpq_t y;
|
||||
|
||||
mag_init(b);
|
||||
fmpq_init(y);
|
||||
|
||||
acb_mat_bound_frobenius_norm(b, A);
|
||||
mag_mul(b, b, b);
|
||||
mag_get_fmpq(y, b);
|
||||
|
||||
if (fmpq_cmp(q, y) > 0)
|
||||
{
|
||||
flint_printf("FAIL (bound, iter = %wd)\n", iter);
|
||||
flint_printf("n = %wd, prec = %wd\n", n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("Qr = \n"); fmpq_mat_print(Qr); flint_printf("\n\n");
|
||||
flint_printf("Qi = \n"); fmpq_mat_print(Qi); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(Q)^2 = \n");
|
||||
fmpq_print(q); flint_printf("\n\n");
|
||||
|
||||
flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
flint_printf("bound_frobenius_norm(A)^2 = \n");
|
||||
mag_printd(b, 15); flint_printf("\n\n");
|
||||
flint_printf("bound_frobenius_norm(A)^2 = \n");
|
||||
mag_print(b); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
mag_clear(b);
|
||||
fmpq_clear(y);
|
||||
}
|
||||
|
||||
fmpq_mat_clear(Qr);
|
||||
fmpq_mat_clear(Qi);
|
||||
fmpq_clear(q);
|
||||
acb_mat_clear(A);
|
||||
}
|
||||
|
||||
/* check trace(A^H A) = frobenius_norm(A)^2 */
|
||||
for (iter = 0; iter < 10000; iter++)
|
||||
{
|
||||
slong m, n, prec;
|
||||
acb_mat_t A, AH, AHA;
|
||||
acb_t t;
|
||||
|
||||
prec = 2 + n_randint(state, 200);
|
||||
|
||||
m = n_randint(state, 10);
|
||||
n = n_randint(state, 10);
|
||||
|
||||
acb_mat_init(A, m, n);
|
||||
acb_mat_init(AH, n, m);
|
||||
acb_mat_init(AHA, n, n);
|
||||
acb_init(t);
|
||||
|
||||
acb_mat_randtest(A, state, 2 + n_randint(state, 100), 10);
|
||||
_acb_mat_conjugate_transpose(AH, A);
|
||||
acb_mat_mul(AHA, AH, A, prec);
|
||||
acb_mat_trace(t, AHA, prec);
|
||||
acb_sqrt(t, t, prec);
|
||||
|
||||
/* check the norm bound */
|
||||
{
|
||||
mag_t low, frobenius;
|
||||
|
||||
mag_init(low);
|
||||
acb_get_mag_lower(low, t);
|
||||
|
||||
mag_init(frobenius);
|
||||
acb_mat_bound_frobenius_norm(frobenius, A);
|
||||
|
||||
if (mag_cmp(low, frobenius) > 0)
|
||||
{
|
||||
flint_printf("FAIL (frobenius norm bound)\n", iter);
|
||||
flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("lower(sqrt(trace(A^H A))) = \n");
|
||||
mag_printd(low, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("bound_frobenius_norm(A) = \n");
|
||||
mag_printd(frobenius, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
mag_clear(low);
|
||||
mag_clear(frobenius);
|
||||
}
|
||||
|
||||
/* check the norm interval */
|
||||
{
|
||||
arb_t frobenius;
|
||||
|
||||
arb_init(frobenius);
|
||||
acb_mat_frobenius_norm(frobenius, A, prec);
|
||||
|
||||
if (!arb_overlaps(acb_realref(t), frobenius))
|
||||
{
|
||||
flint_printf("FAIL (frobenius norm overlap)\n", iter);
|
||||
flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("sqrt(trace(A^H A)) = \n");
|
||||
acb_printd(t, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("frobenius_norm(A) = \n");
|
||||
arb_printd(frobenius, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_clear(frobenius);
|
||||
}
|
||||
|
||||
acb_mat_clear(A);
|
||||
acb_mat_clear(AH);
|
||||
acb_mat_clear(AHA);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -25,16 +25,6 @@
|
|||
|
||||
#include "acb_mat.h"
|
||||
|
||||
static void
|
||||
_acb_mat_conjugate_transpose(acb_mat_t B, const acb_mat_t A)
|
||||
{
|
||||
slong i, j;
|
||||
acb_mat_transpose(B, A);
|
||||
for (i = 0; i < acb_mat_nrows(B); i++)
|
||||
for (j = 0; j < acb_mat_ncols(B); j++)
|
||||
acb_conj(acb_mat_entry(B, i, j), acb_mat_entry(B, i, j));
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
|
|
@ -134,6 +124,7 @@ int main()
|
|||
|
||||
flint_printf("trace(ab) = \n"); acb_printd(trab, 15); flint_printf("\n\n");
|
||||
flint_printf("trace(ba) = \n"); acb_printd(trba, 15); flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
acb_clear(trab);
|
||||
|
|
@ -145,65 +136,6 @@ int main()
|
|||
acb_mat_clear(ba);
|
||||
}
|
||||
|
||||
/* check trace(A^H A) = frobenius_norm(A)^2 */
|
||||
for (iter = 0; iter < 10000; iter++)
|
||||
{
|
||||
slong m, n, prec;
|
||||
acb_mat_t A, AH, AHA;
|
||||
acb_t t;
|
||||
mag_t low, fro;
|
||||
|
||||
prec = 2 + n_randint(state, 200);
|
||||
|
||||
m = n_randint(state, 10);
|
||||
n = n_randint(state, 10);
|
||||
|
||||
acb_mat_init(A, m, n);
|
||||
acb_mat_randtest(A, state, 2 + n_randint(state, 100), 10);
|
||||
|
||||
acb_mat_init(AH, n, m);
|
||||
_acb_mat_conjugate_transpose(AH, A);
|
||||
|
||||
acb_mat_init(AHA, n, n);
|
||||
acb_mat_mul(AHA, AH, A, prec);
|
||||
|
||||
acb_init(t);
|
||||
acb_mat_trace(t, AHA, prec);
|
||||
acb_sqrt(t, t, prec);
|
||||
|
||||
mag_init(low);
|
||||
acb_get_mag_lower(low, t);
|
||||
|
||||
mag_init(fro);
|
||||
acb_mat_bound_fro_norm(fro, A);
|
||||
|
||||
if (mag_cmp(low, fro) > 0)
|
||||
{
|
||||
flint_printf("FAIL (frobenius norm)\n", iter);
|
||||
flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("lower(sqrt(trace(A^H A))) = \n");
|
||||
mag_printd(low, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("upper(frobenius_norm(A)) = \n");
|
||||
mag_printd(fro, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
acb_clear(t);
|
||||
|
||||
mag_clear(low);
|
||||
mag_clear(fro);
|
||||
|
||||
acb_mat_clear(A);
|
||||
acb_mat_clear(AH);
|
||||
acb_mat_clear(AHA);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
|
|
|
|||
|
|
@ -142,7 +142,9 @@ void arb_mat_transpose(arb_mat_t mat1, const arb_mat_t mat2);
|
|||
|
||||
void arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A);
|
||||
|
||||
void arb_mat_bound_fro_norm(mag_t b, const arb_mat_t A);
|
||||
void arb_mat_frobenius_norm(arb_t res, const arb_mat_t A, slong prec);
|
||||
|
||||
void arb_mat_bound_frobenius_norm(mag_t b, const arb_mat_t A);
|
||||
|
||||
/* Arithmetic */
|
||||
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@
|
|||
#include "arb_mat.h"
|
||||
|
||||
void
|
||||
arb_mat_bound_fro_norm(mag_t b, const arb_mat_t A)
|
||||
arb_mat_bound_frobenius_norm(mag_t b, const arb_mat_t A)
|
||||
{
|
||||
slong i, j, r, c;
|
||||
mag_t t;
|
||||
51
arb_mat/frobenius_norm.c
Normal file
51
arb_mat/frobenius_norm.c
Normal file
|
|
@ -0,0 +1,51 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2016 Arb authors
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "arb_mat.h"
|
||||
|
||||
void
|
||||
arb_mat_frobenius_norm(arb_t res, const arb_mat_t A, slong prec)
|
||||
{
|
||||
slong i, j, r, c;
|
||||
|
||||
r = arb_mat_nrows(A);
|
||||
c = arb_mat_ncols(A);
|
||||
|
||||
arb_zero(res);
|
||||
|
||||
if (r == 0 || c == 0)
|
||||
return;
|
||||
|
||||
for (i = 0; i < r; i++)
|
||||
{
|
||||
for (j = 0; j < c; j++)
|
||||
{
|
||||
arb_srcptr x = arb_mat_entry(A, i, j);
|
||||
arb_addmul(res, x, x, prec);
|
||||
}
|
||||
}
|
||||
|
||||
arb_sqrtpos(res, res, prec);
|
||||
}
|
||||
235
arb_mat/test/t-frobenius_norm.c
Normal file
235
arb_mat/test/t-frobenius_norm.c
Normal file
|
|
@ -0,0 +1,235 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2016 Arb authors
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "arb_mat.h"
|
||||
|
||||
static void
|
||||
_fmpq_mat_sum_of_squares(fmpq_t res, const fmpq_mat_t Q)
|
||||
{
|
||||
slong i, j;
|
||||
fmpq_zero(res);
|
||||
for (i = 0; i < fmpq_mat_nrows(Q); i++)
|
||||
{
|
||||
for (j = 0; j < fmpq_mat_ncols(Q); j++)
|
||||
{
|
||||
fmpq_addmul(res, fmpq_mat_entry(Q, i, j), fmpq_mat_entry(Q, i, j));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("frobenius_norm....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
/* compare to the exact rational norm */
|
||||
for (iter = 0; iter < 10000; iter++)
|
||||
{
|
||||
fmpq_mat_t Q;
|
||||
fmpq_t q;
|
||||
arb_mat_t A;
|
||||
slong n, qbits, prec;
|
||||
|
||||
n = n_randint(state, 8);
|
||||
qbits = 1 + n_randint(state, 100);
|
||||
prec = 2 + n_randint(state, 200);
|
||||
|
||||
fmpq_mat_init(Q, n, n);
|
||||
fmpq_init(q);
|
||||
|
||||
arb_mat_init(A, n, n);
|
||||
|
||||
fmpq_mat_randtest(Q, state, qbits);
|
||||
_fmpq_mat_sum_of_squares(q, Q);
|
||||
|
||||
arb_mat_set_fmpq_mat(A, Q, prec);
|
||||
|
||||
/* check that the arb interval contains the exact value */
|
||||
{
|
||||
arb_t a;
|
||||
arb_init(a);
|
||||
|
||||
arb_mat_frobenius_norm(a, A, prec);
|
||||
arb_mul(a, a, a, prec);
|
||||
|
||||
if (!arb_contains_fmpq(a, q))
|
||||
{
|
||||
flint_printf("FAIL (containment, iter = %wd)\n", iter);
|
||||
flint_printf("n = %wd, prec = %wd\n", n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("Q = \n"); fmpq_mat_print(Q); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(Q)^2 = \n");
|
||||
fmpq_print(q); flint_printf("\n\n");
|
||||
|
||||
flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(A)^2 = \n");
|
||||
arb_printd(a, 15); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(A)^2 = \n");
|
||||
arb_print(a); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_clear(a);
|
||||
}
|
||||
|
||||
/* check that the upper bound is not less than the exact value */
|
||||
{
|
||||
mag_t b;
|
||||
fmpq_t y;
|
||||
|
||||
mag_init(b);
|
||||
fmpq_init(y);
|
||||
|
||||
arb_mat_bound_frobenius_norm(b, A);
|
||||
mag_mul(b, b, b);
|
||||
mag_get_fmpq(y, b);
|
||||
|
||||
if (fmpq_cmp(q, y) > 0)
|
||||
{
|
||||
flint_printf("FAIL (bound, iter = %wd)\n", iter);
|
||||
flint_printf("n = %wd, prec = %wd\n", n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("Q = \n"); fmpq_mat_print(Q); flint_printf("\n\n");
|
||||
flint_printf("frobenius_norm(Q)^2 = \n");
|
||||
fmpq_print(q); flint_printf("\n\n");
|
||||
|
||||
flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
flint_printf("bound_frobenius_norm(A)^2 = \n");
|
||||
mag_printd(b, 15); flint_printf("\n\n");
|
||||
flint_printf("bound_frobenius_norm(A)^2 = \n");
|
||||
mag_print(b); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
mag_clear(b);
|
||||
fmpq_clear(y);
|
||||
}
|
||||
|
||||
fmpq_mat_clear(Q);
|
||||
fmpq_clear(q);
|
||||
arb_mat_clear(A);
|
||||
}
|
||||
|
||||
/* check trace(A^T A) = frobenius_norm(A)^2 */
|
||||
for (iter = 0; iter < 10000; iter++)
|
||||
{
|
||||
slong m, n, prec;
|
||||
arb_mat_t A, AT, ATA;
|
||||
arb_t t;
|
||||
|
||||
prec = 2 + n_randint(state, 200);
|
||||
|
||||
m = n_randint(state, 10);
|
||||
n = n_randint(state, 10);
|
||||
|
||||
arb_mat_init(A, m, n);
|
||||
arb_mat_init(AT, n, m);
|
||||
arb_mat_init(ATA, n, n);
|
||||
arb_init(t);
|
||||
|
||||
arb_mat_randtest(A, state, 2 + n_randint(state, 100), 10);
|
||||
arb_mat_transpose(AT, A);
|
||||
arb_mat_mul(ATA, AT, A, prec);
|
||||
arb_mat_trace(t, ATA, prec);
|
||||
arb_sqrt(t, t, prec);
|
||||
|
||||
/* check the norm bound */
|
||||
{
|
||||
mag_t low, frobenius;
|
||||
|
||||
mag_init(low);
|
||||
arb_get_mag_lower(low, t);
|
||||
|
||||
mag_init(frobenius);
|
||||
arb_mat_bound_frobenius_norm(frobenius, A);
|
||||
|
||||
if (mag_cmp(low, frobenius) > 0)
|
||||
{
|
||||
flint_printf("FAIL (bound)\n", iter);
|
||||
flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("lower(sqrt(trace(A^T A))) = \n");
|
||||
mag_printd(low, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("bound_frobenius_norm(A) = \n");
|
||||
mag_printd(frobenius, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
mag_clear(low);
|
||||
mag_clear(frobenius);
|
||||
}
|
||||
|
||||
/* check the norm interval */
|
||||
{
|
||||
arb_t frobenius;
|
||||
|
||||
arb_init(frobenius);
|
||||
arb_mat_frobenius_norm(frobenius, A, prec);
|
||||
|
||||
if (!arb_overlaps(t, frobenius))
|
||||
{
|
||||
flint_printf("FAIL (overlap)\n", iter);
|
||||
flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("sqrt(trace(A^T A)) = \n");
|
||||
arb_printd(t, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("frobenius_norm(A) = \n");
|
||||
arb_printd(frobenius, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_clear(frobenius);
|
||||
}
|
||||
|
||||
arb_mat_clear(A);
|
||||
arb_mat_clear(AT);
|
||||
arb_mat_clear(ATA);
|
||||
arb_clear(t);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -136,65 +136,6 @@ int main()
|
|||
arb_mat_clear(ba);
|
||||
}
|
||||
|
||||
/* check trace(A^T A) = frobenius_norm(A)^2 */
|
||||
for (iter = 0; iter < 10000; iter++)
|
||||
{
|
||||
slong m, n, prec;
|
||||
arb_mat_t A, AT, ATA;
|
||||
arb_t t;
|
||||
mag_t low, fro;
|
||||
|
||||
prec = 2 + n_randint(state, 200);
|
||||
|
||||
m = n_randint(state, 10);
|
||||
n = n_randint(state, 10);
|
||||
|
||||
arb_mat_init(A, m, n);
|
||||
arb_mat_randtest(A, state, 2 + n_randint(state, 100), 10);
|
||||
|
||||
arb_mat_init(AT, n, m);
|
||||
arb_mat_transpose(AT, A);
|
||||
|
||||
arb_mat_init(ATA, n, n);
|
||||
arb_mat_mul(ATA, AT, A, prec);
|
||||
|
||||
arb_init(t);
|
||||
arb_mat_trace(t, ATA, prec);
|
||||
arb_sqrt(t, t, prec);
|
||||
|
||||
mag_init(low);
|
||||
arb_get_mag_lower(low, t);
|
||||
|
||||
mag_init(fro);
|
||||
arb_mat_bound_fro_norm(fro, A);
|
||||
|
||||
if (mag_cmp(low, fro) > 0)
|
||||
{
|
||||
flint_printf("FAIL (frobenius norm)\n", iter);
|
||||
flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec);
|
||||
flint_printf("\n");
|
||||
|
||||
flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("lower(sqrt(trace(A^T A))) = \n");
|
||||
mag_printd(low, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("upper(frobenius_norm(A)) = \n");
|
||||
mag_printd(fro, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_clear(t);
|
||||
|
||||
mag_clear(low);
|
||||
mag_clear(fro);
|
||||
|
||||
arb_mat_clear(A);
|
||||
arb_mat_clear(AT);
|
||||
arb_mat_clear(ATA);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
|
|
|
|||
|
|
@ -161,10 +161,14 @@ Norms
|
|||
Sets *b* to an upper bound for the infinity norm (i.e. the largest
|
||||
absolute value row sum) of *A*.
|
||||
|
||||
.. function:: void acb_mat_bound_fro_norm(mag_t b, const acb_mat_t A)
|
||||
.. function:: void acb_mat_frobenius_norm(acb_t res, const acb_mat_t A, slong prec)
|
||||
|
||||
Sets *b* to an upper bound for the Frobenius norm (i.e. the square root
|
||||
of the sum of squares of magnitudes of entries) of *A*.
|
||||
Sets *res* to the Frobenius norm (i.e. the square root of the sum
|
||||
of squares of entries) of *A*.
|
||||
|
||||
.. function:: void acb_mat_bound_frobenius_norm(mag_t res, const acb_mat_t A)
|
||||
|
||||
Sets *res* to an upper bound for the Frobenius norm of *A*.
|
||||
|
||||
|
||||
Arithmetic
|
||||
|
|
|
|||
|
|
@ -152,10 +152,14 @@ Norms
|
|||
Sets *b* to an upper bound for the infinity norm (i.e. the largest
|
||||
absolute value row sum) of *A*.
|
||||
|
||||
.. function:: void arb_mat_bound_fro_norm(mag_t b, const arb_mat_t A)
|
||||
.. function:: void arb_mat_frobenius_norm(arb_t res, const arb_mat_t A, slong prec)
|
||||
|
||||
Sets *b* to an upper bound for the Frobenius norm (i.e. the square root
|
||||
of the sum of squares of entries) of *A*.
|
||||
Sets *res* to the Frobenius norm (i.e. the square root of the sum
|
||||
of squares of entries) of *A*.
|
||||
|
||||
.. function:: void arb_mat_bound_frobenius_norm(mag_t res, const arb_mat_t A)
|
||||
|
||||
Sets *res* to an upper bound for the Frobenius norm of *A*.
|
||||
|
||||
Arithmetic
|
||||
-------------------------------------------------------------------------------
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue