Re: [eigen] GivensQR

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

2009/8/20 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
> On Wed, 19 Aug 2009, Benoit Jacob wrote:
>
>> I was intrigued by the option of computing/storing only the sequence
>> of Givens rotations, which doesn't seem like common practice at all; i
>> was expecting a GivensQR class to compute/store Q and/or R according
>> to what the user asked for.
>
> The HouseholderQR class does something similar. It doesn't store Q but only
> the sequence of Householder reflections. The matrix Q is only computed when
> the user asks for it.

Yes, for householder this is common practice. For Givens, I don't
know. It needs to be justified by a benchmark, since it is motivated
by hopefully good performance, that's why I tried to benchmark it.

>
>> So I wondered what was the use case for the approach of storing the
>> sequence of Givens rotations, and it seems to me that the best use
>> case for it should be when doing only 1 vector solve. Indeed, in that
>> case, avoiding the overhead of computing Q or R should be the most
>> beneficial.
>
> Another use case is solving an overdetermined system in the least-square
> sense. You don't need to compute Q; you only need to be able to multiply Q
> with a vector.

This is better done differently:

I understand what you mean -- store only the sequence of givens
rotations while computing the QR decompositions, then apply this
sequence of givens rotations to the vector.

However this isn't optimally efficient: the sequence of givens
rotations is traversed twice, it would be better to traverse it only
once, applying the givens rotations on the fly.

This is actually very easy to achieve with any QR decomposition algo
that allows to compute only R: just append your vector V to the right
of your matrix, as an additional column. Then in the resulting R
matrix, the last column is Q^t V. OK you wanted QV instead of Q^t V
but for that, just do instead a LQ decomposition of the transposed
matrix instead (yes, we want to have LQ, RQ etc... as soon as we find
time to code that). Also one might object that only the last column of
R is wanted so one should avoid computing the rest of R, but actually
R is workspace for the decomposition so one can't compute the
necessary (givens or householder) transformations without computing R.

----> actually I didn't think of that at once, but I could have made
my point much quicker: one can't compute the givens rotations without
computing R at the same time.

anyway, here's a benchmark (attached bench.cpp) for computing only R:

=== 12:51:40 ~$ g++ bench.cpp -o bench -O2 -DNDEBUG -I givensqr/ && ./bench
*** size 3 ***
GivensQR:  0.454789 sec
HouseholderQR:0.735767 sec
*** size 4 ***
GivensQR:  0.499081 sec
HouseholderQR:0.626632 sec
*** size 20 ***
GivensQR:  1.77154 sec
HouseholderQR:0.786832 sec
*** size 200 ***
GivensQR:  10.9403 sec
HouseholderQR:3.42356 sec


Indeed that does make GivensQR look good for small sizes.

Benoit

#include <Eigen/Eigen>
#include <bench/BenchTimer.h>

using namespace std;
using namespace Eigen;


template<int size, template<typename T> class decomp> double bench()
{
  enum { repeat = 10000000/(size*size),
         crows = size>10?Dynamic:size,
         ccols = size>10?Dynamic:size+1 };
  typedef Matrix<float,crows,ccols> mtype;
  typedef Matrix<float,crows,1> vtype;
  
  mtype m = mtype::Random(size,size+1);

  BenchTimer t;
  
  t.start();
  for(int i = 0; i < repeat; i++)
  {
     decomp<mtype> d(m);
     d.matrixR();
  }
  t.stop();
  return t.value();
}

int main()
{
  cout << "*** size 3 ***" << endl;
  cout << "GivensQR:  " << bench<3, GivensQR>() << " sec" << endl;
  cout << "HouseholderQR:" << bench<3, HouseholderQR>() << " sec" << endl;
  

  cout << "*** size 4 ***" << endl;
  cout << "GivensQR:  " << bench<4, GivensQR>() << " sec" << endl;
  cout << "HouseholderQR:" << bench<4, HouseholderQR>() << " sec" << endl;

  cout << "*** size 20 ***" << endl;
  cout << "GivensQR:  " << bench<20, GivensQR>() << " sec" << endl;
  cout << "HouseholderQR:" << bench<20, HouseholderQR>() << " sec" << endl;

  cout << "*** size 200 ***" << endl;
  cout << "GivensQR:  " << bench<200, GivensQR>() << " sec" << endl;
  cout << "HouseholderQR:" << bench<200, HouseholderQR>() << " sec" << endl;

}