Re: [eigen] Specialized QR |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Specialized QR*From*: Andrea Arteaga <yo.eres@xxxxxxxxx>*Date*: Wed, 20 May 2009 00:25:49 +0200*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:received:from:to:subject:date :user-agent:references:in-reply-to:mime-version:content-type :message-id; bh=/KmzEEEWz1hy/zxQAthIWOf7hMU0xc/fBwzjna/3LpM=; b=s7qjQQakFTe73XRlOq1hSlYIP1EUkBZtKm5qCKLoNyQvWVEHLfFVfdrxPXvExhiQdk yrBNw/M6IJXp8jbeWo0zab5EjLKIxYs8w+qAvyGH+zfCOaMVWSH0n5lk1mlWi9MT1mgi RXBEMdEPc+EWBBzOZVDKsNJ+iN6rnKewSNkVY=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=from:to:subject:date:user-agent:references:in-reply-to:mime-version :content-type:message-id; b=gtkTwlw3EMGIbEcrlHJpkOEZpdxzNAyvxrSrGFcwhewoFfNNlq6IGJOcU4BHi+Wr1e 95dttOyilaT2EQy7GOr/evovNkbs28G5tQC2BEkVLAoJRXmGp835FoAy2lWh1XN9p1oc WVaM0Q6kCrI2DUXenRw4RqcRtnXeThMdxoQVk=

I did a mistake with my benchmark. My class isn't 4 times faster as the QR class! In facts, the performance of the two classes was almost the same. I used only dynamically-sized matriced, but in theory thereis no difference between dynamic and fixed (the function is exactly the same). I made four graphics, which explain where QR_Givens is faster than QR. I noticed that for almost square matrices Givens is faster, but for MxN matrices, where M>2*N, Householder is far better. For complex ones, Givens is only faster with sizes up to 10x10; then it becomes far slower. In facts this topic isn't very well-named. I wrote a new class, I didn't any specialization. The benchmark shows that Givens is actually faster for small-sized (and real) matrices, but in the general case Householder is better. I think that we should try with Gram-Schmidt for small-sized matrices with many rows. I think that this behaviour is more general. Do you agree? Andrea Arteaga |

**Attachment:
PerformanceComplex.png**

**Attachment:
PerformanceReal.png**

**Attachment:
TimeComplex.png**

**Attachment:
TimeReal.png**

// 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; typedef typename NumTraits<Scalar>::Real RealScalar; enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime, isDynamic = (MatrixType::RowsAtCompileTime == Dynamic) ? 1 : 0, isComplex = NumTraits<Scalar>::IsComplex, 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) { m_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, Rows, Cols> m_R; MatrixQType m_Q; }; template<class MatrixType> inline void QR_Givens<MatrixType>::compute_() { Scalar a, b, o1, o2, tmp; RealScalar r; const int rows = m_R.rows(), cols = m_R.cols(); const int q_cols = cols < rows ? cols : rows-1; m_Q.resize(rows, rows); m_Q.setIdentity(); for(int k1 = 0; k1 < q_cols; ++k1) { for (int k2 = rows - 1; k2 > k1; --k2) { b = m_R.coeff(k2, k1); a = m_R.coeff(k1, k1); r = ei_sqrt(ei_real(ei_conj(a)*a + ei_conj(b)*b)); o1 = a/r; // FIXME: o2 = b/r; // what if r == 0 <=> (a==0 && b==0) ? for (int i = k1; i < cols; ++i){ tmp = m_R.coeff(k1, i); m_R.coeffRef(k1, i) = ei_conj(o1)*m_R.coeff(k1, i) + ei_conj(o2)*m_R.coeff(k2, i); m_R.coeffRef(k2, i) = o1*m_R.coeff(k2, i) - o2*tmp; } for (int i = 0; i < rows; ++i){ tmp = m_Q.coeff(i, k1); m_Q.coeffRef(i, k1) = o1*m_Q.coeff(i, k1) + o2*m_Q.coeff(i, k2); m_Q.coeffRef(i, k2) = ei_conj(o1)*m_Q.coeff(i, k2) - ei_conj(o2)*tmp; } } } } template<class MatrixType> typename QR_Givens<MatrixType>::MatrixRType QR_Givens<MatrixType>::matrixR() const { return Part<MatrixRType, UpperTriangular>(m_R); } template<class MatrixType> typename QR_Givens<MatrixType>::MatrixQType QR_Givens<MatrixType>::matrixQ() const { return m_Q; } #endif // QR_GIVENS_H

**Follow-Ups**:**Re: [eigen] Specialized QR***From:*Benoit Jacob

**References**:**Re: [eigen] Specialized QR***From:*Andrea Arteaga

**Re: [eigen] Specialized QR***From:*Benoit Jacob

**Re: [eigen] Specialized QR***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: std::complex standard packing requirements was Re: [eigen] FFT for Eigen** - Next by Date:
**Re: [eigen] Specialized QR** - Previous by thread:
**Re: [eigen] Specialized QR** - Next by thread:
**Re: [eigen] Specialized QR**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |