Re: [eigen] Specialized QR

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

OK, I finally had time to look at your (good) code.

Find attached an updated version where i correct 2 issues:
- we use 2 spaces to indent, no tabs. i changed it for you as it's
sometimes nontrivial to do automatically and my text editor (kate)
happens to make that easy.
- i added copyright+license info.

There are more issues:
- your class members are:
  Matrix<Scalar, 2, m_cols> m_qr;
  Matrix<Scalar, Rows, Cols> R_;
first, all members should start with a m_ prefix.
Second, is the member m_qr really well named if you have to store R separately?

- you use standard cmath functions like
      r = sqrt(a*a + b*b);
this should use ei_sqrt so it gives the user a chance to extend eigen
to work with a custom scalar type and using your QR with it. Btw, if
you wanted to use a standard function here, that would be std::sqrt.
The plain sqrt you're using is a C function defined only for double's
so you'd get warnings with floats.

- you have much room for optimization here in this inner loop:
      for (int i = k1; i < Cols; ++i){
        tmp = R_.coeff(k1, i);
        R_.coeffRef(k1, i) = o1*R_.coeff(k1, i) + o2*R_.coeff(k2, i);
        R_.coeffRef(k2, i) = o1*R_.coeff(k2, i) - o2*tmp;
      }
This should use Eigen expressions (to benefit from vectorization and
unrolling), but I agree that it's nontrivial because of the
interdependence (which forced you to introduce tmp). You can leave it
like that and let us deal with this, I agree that it's a rather tough
one to "Eigenify" (I'd say that it is best handled by a custom
unary-expression functor applied to matrix R_  -- and that matrix
should then be made row-major).
We have a similar issue below:
        tmprow = Q.row(k1);
        Q.row(k1) = o1 * tmprow + o2 * Q.row(k2);
        Q.row(k2) = o1 * Q.row(k2) - o2 * tmprow;
    }
Here you chose to use Eigen expressions but, the way you do it, you
traverse each row twice instead of once. Again it's the same problem,
it's OK to leave it to us to optimize.

2009/5/13 Andrea Arteaga <yo.eres@xxxxxxxxx>:
>> For what sizes?
> 3x3 up to 150x150, inclusive rectangular matrices (e.g.: 100x50, 20x10,...).
> All results are similars.

OK, interesting, now something occured to me: the most common use case
of QR is to get the matrixQ() and in both QR decompositions, this
involves a big amount of extra computation in addition to the QR
decomposition itself. So if you have time could you benchmark the
whole (QR dec + matrixQ()) ? instead of only benchmarking the QR dec
itself ?

>> So: R rectangular and Q "really unitary".
>> Notice that i say unitary instead of orthogonal because i want support
>> for complex matrices.
> Ok.
> With complex I haven't done much tests, but the code should work.

The one and only tricky thing to keep in mind, is that for complex
numbers, the dot product of vectors x and y becomes

x1 * conj(y1) + ...

i.e. complex conjugation is involved. This reflects in the fact that
your Givens rotations are going to be of the form
conj(a)     conj(b)
-b           a

Forgetting a complex conjugation is the only thing that can make your
code fail to work with complex numbers.

Cheers,
Benoit

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Andrea Arteaga <yo.eres@xxxxxxxxx>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.


#ifndef QR_GIVENS_H
#define QR_GIVENS_H

/**
  * QR_Givens
  *
  * QR decomposition of a matrix
  *
  * MatrixType is the type of the matrix of which we are computing the QR decomposition
  *
  * This class performs a QR decomposition using Givens rotations.
  */
template<class MatrixType> class QR_Givens
{
public:
  typedef typename MatrixType::Scalar Scalar;

  enum {
    Rows = MatrixType::RowsAtCompileTime,
    Cols = MatrixType::ColsAtCompileTime,
    q_cols = Cols < Rows ? Cols : Rows-1,
    m_cols = Cols < Rows ? (Cols*Rows - Cols*(Cols+1)/2) : (Rows*(Rows-1)/2)
  };

  typedef Matrix<Scalar, Rows, Rows> MatrixQType;
  typedef Matrix<Scalar, Rows, Cols> MatrixRType;

  QR_Givens(const MatrixType& matrix) {
    R_ = matrix;
    compute_();
  }

  /** Returns the upper triangular matrix.
    *
    * Dimensions are the same as the matrix of which *this is the QR decomposition.
    */
  MatrixRType matrixR() const;


  /** Returns the orthogonal matrix.
    *
    * This is a quadratic matrix with the same number of rows as the matrix of which *this
    *  is the QR decomposition.
    */
  MatrixQType matrixQ() const;

  /** Returns the rank of the matrix of which *this is the QR decomposition.
    *
    * This method is not stable yet.
    */
  int rank() const;

private:
  inline void compute_();

  Matrix<Scalar, 2, m_cols> m_qr;
  Matrix<Scalar, Rows, Cols> R_;
};


template<class MatrixType>
inline void QR_Givens<MatrixType>::compute_()
{
  Scalar a, b, r, o1, o2, tmp;
  int j = 0;
  for(int k1 = 0; k1 < q_cols; ++k1) {
    for (int k2 = Rows - 1; k2 > k1; --k2) {
      b = R_.coeff(k2, k1);
      if (b == Scalar(0)) {
        m_qr.coeffRef(0, j) = Scalar(1);
        m_qr.coeffRef(1, j++) = Scalar(0);
        continue;
      }
      a = R_.coeff(k1, k1);
      r = sqrt(a*a + b*b);
      o1 = a/r;
      o2 = b/r;

      m_qr.coeffRef(0, j) = o1;
      m_qr.coeffRef(1, j++) = o2;

      for (int i = k1; i < Cols; ++i){
        tmp = R_.coeff(k1, i);
        R_.coeffRef(k1, i) = o1*R_.coeff(k1, i) + o2*R_.coeff(k2, i);
        R_.coeffRef(k2, i) = o1*R_.coeff(k2, i) - o2*tmp;
      }
    }
  }
}

template<class MatrixType>
typename QR_Givens<MatrixType>::MatrixRType
QR_Givens<MatrixType>::matrixR() const
{
  return Part<MatrixRType, UpperTriangular>(R_);
}

template<class MatrixType>
typename QR_Givens<MatrixType>::MatrixQType
QR_Givens<MatrixType>::matrixQ() const
{
  MatrixQType Q;
  Q.setIdentity();
  int j = 0;
  Matrix<Scalar, 1, Rows> tmprow;
  Scalar o1, o2;
  for(int k1 = 0; k1 < q_cols; ++k1) {
    for(int k2 = Rows-1; k2 > k1; --k2){
        o1 = m_qr.coeff(0, j);
        o2 = m_qr.coeff(1, j++);
        tmprow = Q.row(k1);
        Q.row(k1) = o1 * tmprow + o2 * Q.row(k2);
        Q.row(k2) = o1 * Q.row(k2) - o2 * tmprow;
    }
  }

  return Q.transpose();
}

#endif // QR_GIVENS_H