Re: [eigen] Blocked QR algorithm - lapack compatible ?

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Ok to be more precise LAPACK now provides two methods to generate Householder reflectors: dlarfg and dlarfp. The former works like ours makeHouseholder function. So that why I said our hcoeffs is the same than Lapack. But recently they added a second routine dlarfp which enforce positive diagonal entries. See this paper for the "why" and the "how": www.netlib.org/lapack/lawnspdf/lawn203.pdf. Nevertheless this discussion is orthogonal to the initial one about the way data are stored in a compact fashion.

gael

On Tue, Jun 8, 2010 at 10:58 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Really ??? When I look at Lapack documentation our storage should really be the same:
- the upper triangular part (including the diagonal) represents the R matrix
- the strictly lower triangular part of the matrix stores the householder vectors v_i, and the householder coefficients h_i are stored in a separate vector (called tau in Lapack) such that each Householder reflector H_i is defined as:

H_i = I - h_i w_i w_i^T

where w_i^T = [0 .. 0 1 v_i^T]

The above is true for both Eigen and Lapack.

Actually, after checking again, the only difference I found is in the way the householder reflectors are computed, so I was wrong. The difference is that Lapack makes sure the diagonal entries are all >0 (another difference is that their routines is more careful about numerical precision). But that difference should not make any difference in the way the compact matrixQR and the hCoeffs vector (or tau) are used. To be sure I've changed our makeHouseholder routine to make beta>0, and the following piece of code

   int s = 10;
  MatrixXd m(10,10);
  m.setRandom();
 
  HouseholderQR<MatrixXd> qr(m);
  std::cerr << "Eigen's QR:\n" << qr.matrixQR() << "\n\nh_coeffs = " << qr.hCoeffs().transpose() << "\n";
 
  // lapack
  int info;
  VectorXd tau(s), w(s);
  dgeqr2_(&s,&s,m.data(),&s,tau.data(),w.data(),&info);
  std::cerr << "Lapack's QR:\n" << m << "\n\ntau = " << tau.transpose() << "\n";


give me:



Eigen's QR:                                                                                       
    1.68335    0.299504    -1.02341    0..421776    0.489904    0.354943     0.16256    0.811929     1.36868    0.538918
   0.210607     1.54456    -1.12169    0.702928     0.49701   -0.305953  -0.0807894    0.354537    0.165679    0.326214
  -0.564518    0.382985     1.66599   -0.728678   -0.461996     1.42203  -0.0842296     1.18131    0.604993   -0.222919
  -0.595109    0.146875   -0.214137    0.955957    0.634162    -0.90593     1.03467   -0.307614    0.218605    0.211342
  -0.820851   -0.511857    0.129739   0.0766891     1.65791   -0.698759    0.528311   -0.300805    0.241369    0.365926
   0.603102   -0.856756   -0.107592     -0.2841  -0.0171067    0.880522   -0.788467   -0.607822 -0.00167755    -1.03731
   0.328576    -0.31682    0.220283   0.0138095   -0.342978   0.0463581    0.979142   -0.513038   -0.585443  -0.0788398
  -0.534867   -0.206075    0.243301   -0.393836   -0.247999    0.109516    0.403584    0.362374    0.271345    0.581348
   0.443132    0.464533   -0.269054   0.0523833    0.309469    0.497054    0.845007     -2.0881    0.284359    0.597509
   -0.10762   -0.145642    0.181877   -0.260819   -0.355856  -0.0033688    0.122173     1.31234    -1.22468     1.02892

h_coeffs = 0.595821 0.786126  1..55346  1.52362   1.4267  1.58577  1.05717  0..28239 0.800049        0
Lapack's QR:
    1.68335    0.299504    -1.02341    0.421776    0.489904    0.354943     0.16256    0.811929     1.36868    0.538918
   0.210607     1.54456    -1.12169    0.702928     0.49701   -0.305953  -0.0807894    0.354537    0.165679    0.326214
  -0.564518    0.382985     1.66599   -0.728678   -0.461996     1.42203  -0.0842296     1.18131    0.604993   -0.222919
  -0.595109    0.146875   -0.214137    0.955957    0.634162    -0.90593     1.03467   -0.307614    0.218605    0.211342
  -0.820851   -0.511857    0.129739   0.0766891     1.65791   -0.698759    0.528311   -0.300805    0.241369    0.365926
   0.603102   -0.856756   -0.107592     -0.2841  -0.0171067    0.880522   -0.788467   -0.607822 -0.00167755    -1.03731
   0.328576    -0.31682    0.220283   0.0138095   -0.342978   0.0463581    0.979142   -0.513038   -0.585443  -0.0788398
  -0.534867   -0.206075    0.243301   -0.393836   -0.247999    0.109516    0.403584    0.362374    0.271345    0.581348
   0.443132    0.464533   -0.269054   0.0523833    0.309469    0.497054    0.845007     -2.0881    0.284359    0.597509
   -0.10762   -0.145642    0.181877   -0.260819   -0.355856  -0.0033688    0.122173     1.31234    -1.22468     1.02892

tau = 0.595821 0.786126  1.55346  1.52362   1.4267  1.58577  1.05717  0.28239 0.800049        0


and for the sake on completeness, with the original (i.e. current) makeHouseolder() routine, Eigen gives me:

Eigen's QR:
   -1.68335   -0.299504     1.02341   -0.421776   -0.489904   -0.354943    -0.16256   -0.811929    -1.36868   -0.538918
 -0.0893648    -1.54456     1.12169   -0.702928    -0.49701    0.305953   0.0807894   -0.354537   -0.165679   -0.326214
   0.239536   -0.181557     1.66599   -0.728678   -0.461996     1.42203  -0.0842296     1.18131    0.604993   -0.222919
   0.252516  -0.0204998   -0.478563   -0.955957   -0.634162     0.90593    -1.03467    0.307614   -0.218605   -0.211342
   0.348304    0.447056   -0.364374   -0.114426     1.65791   -0.698759    0.528311   -0.300805    0.241369    0.365926
  -0.255908    0.491757  -0.0821699   -0.448179   0.0034012    0.880522   -0.788467   -0.607822 -0.00167755    -1.03731
  -0.139422     0.16815    0.317734   -0.544321   -0.483593    0.328056    0.979142   -0.513038   -0.585443  -0.0788398
   0.226955    0.206507  -0.0196116    0.576286    0.138886  -0.0462263    0.617459   -0.362374   -0.271345   -0.581348
   -0.18803    -0.36656     0.01843    0.171011    0.248586    0.159354    0.218406  -0.0142162    0.284359    0.597509
  0.0456651    0..110862    0.120987    0.232734   -0.150844  -0.0333722    0.370857  -0.0213748    -0.81316     1.02892

h_coeffs = 1.40418 1.18158 1.34693 1.03855 1.49509 1.76015 1.27674 1.99868 1.20393       0


Gael

On Tue, Jun 8, 2010 at 3:23 AM, Thomas Capricelli <orzel@xxxxxxxxxxxxxxx> wrote:

Hi,

Indeed, the way lapack and eigen stores the householder transformations is different.

Here you can see the code i use to go from one to the other:

http://bitbucket.org/eigen/eigen/src/tip/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h#cl-470

only those lines are important:

------------------------------------------------------------------------

ColPivHouseholderQR<JacobianType> qrfac(fjac);

        fjac = qrfac.matrixQR();

        wa1 = fjac.diagonal();

        fjac.diagonal() = qrfac.hCoeffs();

        permutation = qrfac.colsPermutation();

        // TODO : avoid this:

        for(Index ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors

------------------------------------------------------------------------

I use the eigen qr decompozition, but the code later one is raw lapack and relies on the lapack 'packing'.

In both cases, the packing is done within the (strict) lower triangular part of a matrix, a vector stored on the diagonal, and another vector stored "outside" the matrix.

There are two differences, though:

* 'diagonal' and 'external vector' are inversed between eigen and lapack. That is, the lapack diagonal is what eigen calls hCoeffs, and the eigen diagonal is what lapack calls 'wa1'.

* vectors in the lower triangular parts are 'scaled' (check my code) in one implementation and not in the other one.

Hope this helps!

Thomas

--

Thomas Capricelli <orzel@xxxxxxxxxxxxxxx>

http://www.freehackers.org/thomas

In data lunedì 07 giugno 2010 17:35:19, Benoit Jacob ha scritto:

> Actually, pinging Thomas, because I remember he faced the exact same

> issue. The conversion is not hard to do at all, but there is something

> to do.

>

> Benoit

>

> 2010/6/7 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:

> > I'm pretty sure that out householder stuff isn't LAPACK compatible:

> > either our hCoeffs or our housholder vectors are different. Don't

> > remember which. Since Gael says the hCoeffs are compatible, I guess

> > it's the vectors then, that are different.