[eigen] Re: LU precision tuning

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

So, here are the first results for floats and doubles (attached file
biscuitlu.log).
Hopefully Hauke can get the results for complex numbers: my laptop is
on its knees.
Also attached, for reference, the latest version of biscuitlu.cpp.

The conclusion, for float and double (let's see tomorrow if complex
numbers are the same), is that the table look the same in both cases.
It looks roughly like

size  precision
32:   4 * machine_epsilon
64:   8 * machine_epsilon
128:   16 * machine_epsilon
256:   32 * machine_epsilon
512:   128 * machine_epsilon
1024: 512 * machine_epsilon
2048: 1024 * machine_epsilon

Given that it's harmless to overestimate a little but we should not
underestimate (would make the rank computation wrong), it's obvious
that we need look no further than this formula:

precision = size * machine_epsilon

A nice coincidence: this turns out to be exactly Higham's formula that
we already use in LDLt.

So I implemented it in LU.h. If the results look different for complex
numbers, we'll have to update it.

The next step was supposed to be to fine-tune 'precision' for optimal
precision of the LU decomposition and solver. But in the meanwhile, it
occured to me that it is obvious that the smaller the 'precision', the
more precise the LU and solver. So we don't need to look any further,
we just use the above formula.

So, unless complex numbers do a prank on us, this can be the end of the story.

Cheers,
Benoit

2009/5/7 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
>>>>
>>>
>>
>

repeating each test 50 times
starting void test_lu_rank() [with Scalar = float]
size 32, precision 4.768e-07 = (machine epsilon)*2^2
size 64, precision 9.536e-07 = (machine epsilon)*2^3
size 128, precision 1.9072e-06 = (machine epsilon)*2^4
size 256, precision 3.8144e-06 = (machine epsilon)*2^5
size 512, precision 7.6288e-06 = (machine epsilon)*2^6
size 1024, precision 3.05152e-05 = (machine epsilon)*2^8
size 2048, precision 0.000122061 = (machine epsilon)*2^10
starting void test_lu_rank() [with Scalar = double]
size 32, precision 8.88e-16 = (machine epsilon)*2^2
size 64, precision 1.776e-15 = (machine epsilon)*2^3
size 128, precision 3.552e-15 = (machine epsilon)*2^4
size 256, precision 7.104e-15 = (machine epsilon)*2^5
size 512, precision 2.8416e-14 = (machine epsilon)*2^7
size 1024, precision 5.6832e-14 = (machine epsilon)*2^8
size 2048, precision 2.27328e-13 = (machine epsilon)*2^10

#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>
bool test_lu_rank_one_size_one_precision(int size, typename NumTraits<Scalar>::Real precision)
{
  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);
    if(rank != lu.rank()) return false;
  }
  return true;
}

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 = 32; size <= 1024; size *= 2)
  {
    bool finished = false;
    int i = 0;
    for(RealScalar precision = machine_epsilon<RealScalar>();
        !finished;
        precision *= 2, i++)
    {
      finished = test_lu_rank_one_size_one_precision<Scalar>(size, precision);
      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> >();  
}