[eigen] precision

[Michael Meyer <mjhmeyer@xxxxxxxxx>]

Hello,

This is my first post to the list.
I am very impressed with the design of this library and could not believe the benchmarks when I saw them.
So the following comments are intended to be constructive:

I think the precision (MathFunctions.h) is generally set too low
(I am on version 2.01, haven't switched to 2.1 yet).
Of course it depends on the context what a reasonable value is.

I am thinking about the function ei_is_much_smaller_than when it is used to check that
a number is sufficiently close to zero (probably the major use).
In that case it is perfectly reasonable to set precision to a quantity smaller than unit roundoff (machine epsilon).

A sloppy argument against that would run as follows: if the number u is computed from larger numbers a,b,...
and is itself smaller than roundoff in a,b,... then it cannot possibly be meaningful.
This is not generally true, it depends on the nature of the computation, for example

double u=1e-30, a=1;

u=u/sqrt(1+a);      

the roundoff in the denominator is likely bigger than u but this does not matter for the computation,
the result will still be extremely precise.

I played around with the Hessenberg tridiagonalization applied to the Hilbert matrix H in dimension 400
(scalar type double, H=QTQ', T tridiagonal, Q orthogonal) and checked how close H-QTQ' is to the zero matrix
using the L1-norm.

The norm of the Hilbert matrix is || H || = 550

With the value of precision<double> as in the distribution I got

|| H - QTQ' || = 4.4e-4

I then set precision<double> to 1e-20 and got

|| H - QTQ' || = 2.8e-8.

Setting the pecision even higher did not improve the result.
The code is attached, sorry no includes since I renamed some the files:




/** \f$ L^1 \f$-norm: sum of absolute values of entries.
 */
template<typename MatrixType>
long double norm(const MatrixType& A)
{
    long double sum=0.0;
    for(int i=0; i<A.rows(); ++i)
    for(int j=0; j<A.cols(); ++j) sum+=ei_abs(A(i,j));

    return sum;
}


void testTridiagonal()
{
    int n=400;
    typedef MatrixXd MatrixType;
    MatrixType A(n,n);
   
    for(int i=0; i<n; ++i)
    for(int j=i; j<n; ++j){

        // Hilbert matrix
        A(i,j)=1/double(1+i+j);
        A(j,i)=A(i,j);
    }

    Tridiagonalization<MatrixType> Hess_A(A);
    MatrixType T(Hess_A.matrixT());
    MatrixType Q(Hess_A.matrixQ());

    T=(Q*T)*Q.adjoint();
   

    std::cerr << "`norm(A) = " << norm(A) << std::endl;
    std::cerr << "`norm(Q) = " << norm(Q) << std::endl;
    std::cerr << "`norm(QTQ') = " << norm(T) << std::endl;
    std::cerr << "`norm(A-QTQ') = " << norm(A-T) << std::endl;
    std::cerr << "Input integer to exit: ";

    int exit; cin >> exit;
}

Cheers,

Mike