[eigen] Re: LU precision tuning

[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 = 10;

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 <= 1000; size += 100)
  {
    bool finished = false;
    for(RealScalar precision = machine_epsilon<RealScalar>();
        !finished;
        precision *= 2)
    {
      float result = test_lu_rank_one_size_one_precision<Scalar>(size, precision);
      cout << "size " << size << ", precision " << precision << ", average relative error in rank: " << result << endl;
      if(result < 1.f / (size * repeat)) finished = true;
    }
  }
}

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> >();  
}