Re: [eigen] Specialized QR |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Specialized QR*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Sun, 17 May 2009 03:55:45 +0200*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type; bh=/gc2F/sZZTTps3A1NqQz5FJoH+8X0LMQGqiikDFgRHU=; b=lMBaYt5rkHpWM/iRWwCsRr0RDX4FSC6kII4R4JyQNUo4ISenxWdHc+P+zGVVKh9dsC wHDHTVb0n7BM9nu2gtI6mVIa9jgg27gLvwUJ7sHBla3WYn6xIT2QY7ifkefNeLgbeZtf /uTrGf+fV18a/1Jmqoe+UN9wHotzCDSkb4quI=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; b=bFxU3CNRvT6G9tesuare4wVan0FJsaOhfORCnBM8gNR4mHPn6KR4qNZ675J+Fx+VZe W/tdKcRfmLViijcOGlKaQdvT0gRRRyYpwrPh01HC7JsJLweyMQYiFysX/5x8SQhcoWs+ G4DGftfQaxhxOGJrhqHweuLWeWJqew+SSsvD4=

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

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

**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:*Andrea Arteaga

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Migration to mercurial done** - 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/ |