Re: [eigen] Specialized QR

[Andrea Arteaga <yo.eres@xxxxxxxxx>]

Ok, I worked around the Givens QR decomposition and I find I reached some good result. I implemented a QR_Givens<MatrixType> class, which looks like the "normal" QR<MatrixType> class. The algorithmus is based on fast Givens rotations, where the matrix is updated in each step only there, where is strictly needed. The result is: the QR_Givens class computes the QR decomposition about 4 times faster than the QR class. The relative errors are sometimes bigger than those of QR, but I never had errors bigger than 1e-13.
(I tested if Q^t * Q is the identity and if Q*R = A.)

Furthermore, I noticed that Eigen QR class has (run-time) troubles with mxn matrices, where m<n (I become this error: "../eigen2/Eigen/src/Core/MapBase.h:166: Eigen::MapBase<Derived>::MapBase(const typename Eigen::ei_traits<T>::Scalar*, int, int) [with Derived = Eigen::Block<Eigen::Block<Eigen::Matrix<double, 50, 100, 0, 50, 100>, 50, 1, 1, 32>, 33331, 1, 1, 32>]: Assertion `(data == 0) || ( rows > 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows) && cols > 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols))' failed.", here with Rows = 50, Cols = 100). Well, even if this is not a very important feature, my QR_Givens class works also well with this kind of matrices.

In this moment my class works only with fixed-sized matrices. As I have time, I will implement the dynamically-sized version and try with larger matrices.

I also implemented a rank() method, but this is not very stable yet. It works fine with "normal" matrices, but has problems for some "strange matrices". For example the Matrix

0, 0, 2, 3, 4
0, 1, 4, 1, 3
0, 0, 0, 3, 1
0, 0, 0, 0, 2
0, 0, 0, 0, 0

has rank 4, but QR::rank() will return 0, QR_Givens::rank() will return 3. I haven't implemented the other methods yet, inclusive the solver. Before that, i'd like to test more the decomposition.

Another question: the QR class returns a quadratic matrix R and a rectangular matrix Q; QR_Givens instead returns a rectangular R and a "really orthogonal" quadratic matrix Q. Which behaviour is better? We can implement both methods for each matrix (e.g.: matrixR() and totalMatrixR(), matrixQ() and totalMatrixQ()), if one finds useful to have them.

Well, I attach my code. Have a look and if you find it good, I will send a patch as soon as I implement the rest of the methods.

Cheers.
Andrea Arteaga

#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