[eigen] Re: LU precision tuning

[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)
  {
    for(RealScalar precision = machine_epsilon<RealScalar>();
        precision < sqrt(machine_epsilon<RealScalar>());
        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;
    }
  }
}

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