[eigen] Re: LU precision tuning

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

ok, last last version. After running for a while it turns out this is
really what's useful.

Still open TODO item: write similar program optimizing the precision of solve().

Benoit

2009/5/7 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> ok, last version: test exactly what's needed (previous one ended too
> soon or continued too long)
>
> 2009/5/7 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
>> oops, a little bug. find attached updated program (runs faster, more
>> significant output)
>>
>>
>> 2009/5/7 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
>>> Hi,
>>>
>>> This is mostly for Hauke but the attached file biscuitlu is worthy of
>>> being archived here, so here goes...
>>>
>>> In decompositions such as LU, we need to make decisions when to
>>> consider something as zero, and this depends on an arbitrary cutoff
>>> that can only be tuned by experiment. In LLt, we were lucky to have
>>> Keir who found a research paper telling us the optimal cutoff, but
>>> there doesn't seem to be anything like that for LU.
>>>
>>> So:
>>> - attached patch difflu allows LU to override the default precision level
>>> - attached patch biscuitlu tests the precision of LU::rank() for
>>> matrix sizes from 100 to 1000, for various precisions.
>>>
>>> This should allow us to infer the formula for the maximum precision we
>>> can use to allow for exact rank determination. Once this range of
>>> values is known, the next step will be to find the optimal 'precision"
>>> value to get optimally precise LU::solve(). So if you want to help,
>>> this 2nd program remains to be written !
>>>
>>> Benoit
>>>
>>
>

#include <Eigen/LU>
#include <Eigen/Array>
#include <ctime>

using namespace Eigen;
using namespace std;

const int repeat = 50;

template<typename Derived>
void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m)
{
  typedef typename Derived::RealScalar RealScalar;
  for(int a = 0; a < 3*(m.rows()+m.cols()); a++)
  {
    RealScalar d = Eigen::ei_random<RealScalar>(-1,1);
    int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number
    int j;
    do {
      j = Eigen::ei_random<int>(0,m.rows()-1);
    } while (i==j); // j is another one (must be different)
    m.row(i) += d * m.row(j);

    i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number
    do {
      j = Eigen::ei_random<int>(0,m.cols()-1);
    } while (i==j); // j is another one (must be different)
    m.col(i) += d * m.col(j);
  }
}

template<typename Scalar>
float test_lu_rank_one_size_one_precision(int size, typename NumTraits<Scalar>::Real precision)
{
  float result = 0;
  for(int i = 0; i < repeat; i++)
  {
    int rank = ei_random<int>(1, size-1);
    typedef Matrix<Scalar,Dynamic,Dynamic> Mat;
    Mat m = Mat::Random(size,size);
    for(int i = rank; i < size; i++) m.row(i).setZero();
    doSomeRankPreservingOperations(m);
    LU<Mat> lu(m, precision);
    result += float(abs(rank - lu.rank())) / size;  
  }
  return result / repeat;
}

template<typename Scalar>
void test_lu_rank()
{
  typedef typename NumTraits<Scalar>::Real RealScalar;
#ifdef __GNUC__
  cout << "starting " << __PRETTY_FUNCTION__ << endl;
#else
  cout << "starting " << __func__ << endl;
#endif
  for(int size = 100; size <= 1200; size += 100)
  {
    bool finished = false;
    int i = 0;
    for(RealScalar precision = machine_epsilon<RealScalar>();
        !finished;
        precision *= 2, i++)
    {
      float result = test_lu_rank_one_size_one_precision<Scalar>(size, precision);
      if(result < .99f / (size * repeat)) finished = true;
      if(finished)
      {
        cout << "size " << size << ", precision " << precision << " = (machine epsilon)*2^" << i << endl;
      }
    }
  }
}

int main()
{
  srand((unsigned int) time(NULL));
  cout << "repeating each test " << repeat << " times" << endl;
  test_lu_rank<float>();
  test_lu_rank<double>();
  test_lu_rank<complex<float> >();
  test_lu_rank<complex<double> >();  
}