diff --git a/Eigen/QR b/Eigen/QR index 151c0b31b..a0575040c 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -36,6 +36,8 @@ namespace Eigen { */ #include "src/QR/QR.h" +#include "src/QR/FullPivotingHouseholderQR.h" +#include "src/QR/ColPivotingHouseholderQR.h" #include "src/QR/Tridiagonalization.h" #include "src/QR/EigenSolver.h" #include "src/QR/SelfAdjointEigenSolver.h" diff --git a/Eigen/StdVector b/Eigen/StdVector index a7f39ea37..f37de5ff6 100644 --- a/Eigen/StdVector +++ b/Eigen/StdVector @@ -136,10 +136,8 @@ class vector > { return vector_base::insert(position,x); } void insert(const_iterator position, size_type new_size, const value_type& x) { vector_base::insert(position, new_size, x); } -#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,1) +#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,2) // workaround GCC std::vector implementation - // Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&), - // no no need to workaround ! void resize(size_type new_size, const value_type& x) { if (new_size < vector_base::size()) @@ -147,9 +145,12 @@ class vector > else vector_base::insert(vector_base::end(), new_size - vector_base::size(), x); } -#elif defined(_GLIBCXX_VECTOR) +#elif defined(_GLIBCXX_VECTOR) && (!EIGEN_GNUC_AT_LEAST(4,1)) + // Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&), + // no no need to workaround ! using vector_base::resize; #else + // either GCC 4.1 or non-GCC // default implementation which should always work. void resize(size_type new_size, const value_type& x) { diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 1f4c6bf7a..7cb8853e6 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -745,6 +745,8 @@ template class MatrixBase /////////// QR module /////////// const HouseholderQR householderQr() const; + const ColPivotingHouseholderQR colPivotingHouseholderQr() const; + const FullPivotingHouseholderQR fullPivotingHouseholderQr() const; EigenvaluesReturnType eigenvalues() const; RealScalar operatorNorm() const; diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index e1e106b80..d70344deb 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -61,11 +61,22 @@ template struct ei_product_type enum { Rows = Lhs::RowsAtCompileTime, Cols = Rhs::ColsAtCompileTime, - Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime), - - value = ei_product_type_selector<(Rows >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Rows==1 ? 1 : Small)), - (Cols >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Cols==1 ? 1 : Small)), - (Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small))>::ret + Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime) + }; + + // the splitting into different lines of code here, introducing the _select enums and the typedef below, + // is to work around an internal compiler error with gcc 4.1 and 4.2. +private: + enum { + rows_select = Rows >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Rows==1 ? 1 : Small), + cols_select = Cols >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Cols==1 ? 1 : Small), + depth_select = Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small) + }; + typedef ei_product_type_selector product_type_selector; + +public: + enum { + value = product_type_selector::ret }; }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 6f3e6a788..414aaceca 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -117,6 +117,8 @@ template class Reverse; template class LU; template class PartialLU; template class HouseholderQR; +template class ColPivotingHouseholderQR; +template class FullPivotingHouseholderQR; template class SVD; template class LLT; template class LDLT; diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 769ba3d43..36f02d7ce 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2009 Benoit Jacob +// Copyright (C) 2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -73,10 +74,10 @@ void MatrixBase::makeHouseholder( EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) VectorBlock tail(derived(), 1, size()-1); - RealScalar tailSqNorm; + RealScalar tailSqNorm = size()==1 ? 0 : tail.squaredNorm(); Scalar c0 = coeff(0); - if( (size()==1 || (tailSqNorm=tail.squaredNorm()) == RealScalar(0)) && ei_imag(c0)==RealScalar(0)) + if(tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0)) { *tau = 0; *beta = ei_real(c0); diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index b5940c74b..a6670a49f 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -26,17 +26,32 @@ #ifndef EIGEN_JACOBI_H #define EIGEN_JACOBI_H +/** Applies the counter clock wise 2D rotation of angle \c theta given by its + * cosine \a c and sine \a s to the set of 2D vectors of cordinates \a x and \a y: + * \f$ x = c x - s' y \f$ + * \f$ y = s x + c y \f$ + * + * \sa MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() + */ template void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s); +/** Applies a rotation in the plane defined by \a c, \a s to the rows \a p and \a q of \c *this. + * More precisely, it computes B = J' * B, with J = [c s ; -s' c] and B = [ *this.row(p) ; *this.row(q) ] + * \sa MatrixBase::applyJacobiOnTheRight(), ei_apply_rotation_in_the_plane() + */ template inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) { RowXpr x(row(p)); RowXpr y(row(q)); - ei_apply_rotation_in_the_plane(x, y, c, s); + ei_apply_rotation_in_the_plane(x, y, ei_conj(c), ei_conj(s)); } +/** Applies a rotation in the plane defined by \a c, \a s to the columns \a p and \a q of \c *this. + * More precisely, it computes B = B * J, with J = [c s ; -s' c] and B = [ *this.col(p) ; *this.col(q) ] + * \sa MatrixBase::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane() + */ template inline void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) { @@ -45,6 +60,12 @@ inline void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, S ei_apply_rotation_in_the_plane(x, y, c, s); } +/** Computes the cosine-sine pair (\a c, \a s) such that its associated + * rotation \f$ J = ( \begin{array}{cc} c & s \\ -s' c \end{array} )\f$ + * applied to both the right and left of the 2x2 matrix + * \f$ B = ( \begin{array}{cc} x & y \\ _ & z \end{array} )\f$ yields + * a diagonal matrix A: \f$ A = J' B J \f$ + */ template bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) { @@ -128,12 +149,13 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& const Packet pc = ei_pset1(c); const Packet ps = ei_pset1(s); + ei_conj_helper cj; for(int i=0; i class LU inline void computeInverse(MatrixType *result) const { ei_assert(m_originalMatrix != 0 && "LU is not initialized."); + ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result); } @@ -456,6 +456,7 @@ template typename ei_traits::Scalar LU::determinant() const { ei_assert(m_originalMatrix != 0 && "LU is not initialized."); + ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!"); return Scalar(m_det_pq) * m_lu.diagonal().prod(); } @@ -533,7 +534,8 @@ bool LU::solve( ) const { ei_assert(m_originalMatrix != 0 && "LU is not initialized."); - + if(m_rank==0) return false; + /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. * So we proceed as follows: * Step 1: compute c = Pb. @@ -552,7 +554,8 @@ bool LU::solve( for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i); // Step 2 - m_lu.corner(Eigen::TopLeft,smalldim,smalldim).template triangularView() + m_lu.corner(Eigen::TopLeft,smalldim,smalldim) + .template triangularView() .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols())); if(rows>cols) { @@ -564,11 +567,10 @@ bool LU::solve( if(!isSurjective()) { // is c is in the image of U ? - RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : 0; - for(int col = 0; col < c.cols(); ++col) - for(int row = m_rank; row < c.rows(); ++row) - if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision)) - return false; + RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff(); + if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision)) + return false; } m_lu.corner(TopLeft, m_rank, m_rank) .template triangularView() diff --git a/Eigen/src/QR/ColPivotingHouseholderQR.h b/Eigen/src/QR/ColPivotingHouseholderQR.h new file mode 100644 index 000000000..ed4b84f63 --- /dev/null +++ b/Eigen/src/QR/ColPivotingHouseholderQR.h @@ -0,0 +1,266 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2009 Benoit Jacob +// +// 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 . + +#ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H +#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H + +/** \ingroup QR_Module + * \nonstableyet + * + * \class ColPivotingHouseholderQR + * + * \brief Householder rank-revealing QR decomposition of a matrix + * + * \param MatrixType the type of the matrix of which we are computing the QR decomposition + * + * This class performs a rank-revealing QR decomposition using Householder transformations. + * + * This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal + * numerical stability. + * + * \sa MatrixBase::colPivotingHouseholderQr() + */ +template class ColPivotingHouseholderQR +{ + public: + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) + }; + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix MatrixQType; + typedef Matrix HCoeffsType; + typedef Matrix IntRowVectorType; + typedef Matrix IntColVectorType; + typedef Matrix RowVectorType; + typedef Matrix ColVectorType; + typedef Matrix RealRowVectorType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ColPivotingHouseholderQR::compute(const MatrixType&). + */ + ColPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + + ColPivotingHouseholderQR(const MatrixType& matrix) + : m_qr(matrix.rows(), matrix.cols()), + m_hCoeffs(std::min(matrix.rows(),matrix.cols())), + m_isInitialized(false) + { + compute(matrix); + } + + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which + * *this is the QR decomposition, if any exists. + * + * \param b the right-hand-side of the equation to solve. + * + * \param result a pointer to the vector/matrix in which to store the solution, if any exists. + * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). + * If no solution exists, *result is left with undefined coefficients. + * + * \note The case where b is a matrix is not yet implemented. Also, this + * code is space inefficient. + * + * Example: \include ColPivotingHouseholderQR_solve.cpp + * Output: \verbinclude ColPivotingHouseholderQR_solve.out + */ + template + void solve(const MatrixBase& b, ResultType *result) const; + + MatrixType matrixQ(void) const; + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + */ + const MatrixType& matrixQR() const { return m_qr; } + + ColPivotingHouseholderQR& compute(const MatrixType& matrix); + + const IntRowVectorType& colsPermutation() const + { + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + return m_cols_permutation; + } + + inline int rank() const + { + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + return m_rank; + } + + protected: + MatrixType m_qr; + HCoeffsType m_hCoeffs; + IntRowVectorType m_cols_permutation; + bool m_isInitialized; + RealScalar m_precision; + int m_rank; + int m_det_pq; +}; + +#ifndef EIGEN_HIDE_HEAVY_CODE + +template +ColPivotingHouseholderQR& ColPivotingHouseholderQR::compute(const MatrixType& matrix) +{ + int rows = matrix.rows(); + int cols = matrix.cols(); + int size = std::min(rows,cols); + m_rank = size; + + m_qr = matrix; + m_hCoeffs.resize(size); + + RowVectorType temp(cols); + + m_precision = epsilon() * size; + + IntRowVectorType cols_transpositions(matrix.cols()); + m_cols_permutation.resize(matrix.cols()); + int number_of_transpositions = 0; + + RealRowVectorType colSqNorms(cols); + for(int k = 0; k < cols; ++k) + colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); + RealScalar biggestColSqNorm = colSqNorms.maxCoeff(); + + for (int k = 0; k < size; ++k) + { + int biggest_col_in_corner; + RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner); + biggest_col_in_corner += k; + + // if the corner is negligible, then we have less than full rank, and we can finish early + if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision)) + { + m_rank = k; + for(int i = k; i < size; i++) + { + cols_transpositions.coeffRef(i) = i; + m_hCoeffs.coeffRef(i) = Scalar(0); + } + break; + } + + cols_transpositions.coeffRef(k) = biggest_col_in_corner; + if(k != biggest_col_in_corner) { + m_qr.col(k).swap(m_qr.col(biggest_col_in_corner)); + ++number_of_transpositions; + } + + RealScalar beta; + m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); + m_qr.coeffRef(k,k) = beta; + + m_qr.corner(BottomRight, rows-k, cols-k-1) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + + colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2(); + } + + for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; + for(int k = 0; k < size; ++k) + std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); + + m_det_pq = (number_of_transpositions%2) ? -1 : 1; + m_isInitialized = true; + + return *this; +} + +template +template +void ColPivotingHouseholderQR::solve( + const MatrixBase& b, + ResultType *result +) const +{ + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + const int rows = m_qr.rows(); + const int cols = b.cols(); + ei_assert(b.rows() == rows); + + typename OtherDerived::PlainMatrixType c(b); + + Matrix temp(cols); + for (int k = 0; k < m_rank; ++k) + { + int remainingSize = rows-k; + c.corner(BottomRight, remainingSize, cols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); + } + + m_qr.corner(TopLeft, m_rank, m_rank) + .template triangularView() + .solveInPlace(c.corner(TopLeft, m_rank, c.cols())); + + result->resize(m_qr.cols(), b.cols()); + for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i); + for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero(); +} + +/** \returns the matrix Q */ +template +MatrixType ColPivotingHouseholderQR::matrixQ() const +{ + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + // compute the product H'_0 H'_1 ... H'_n-1, + // where H_k is the k-th Householder transformation I - h_k v_k v_k' + // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] + int rows = m_qr.rows(); + int cols = m_qr.cols(); + int size = std::min(rows,cols); + MatrixType res = MatrixType::Identity(rows, rows); + Matrix temp(rows); + for (int k = size-1; k >= 0; k--) + { + res.block(k, k, rows-k, rows-k) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); + } + return res; +} + +#endif // EIGEN_HIDE_HEAVY_CODE + +/** \return the column-pivoting Householder QR decomposition of \c *this. + * + * \sa class ColPivotingHouseholderQR + */ +template +const ColPivotingHouseholderQR::PlainMatrixType> +MatrixBase::colPivotingHouseholderQr() const +{ + return ColPivotingHouseholderQR(eval()); +} + + +#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h new file mode 100644 index 000000000..0ffcfe88c --- /dev/null +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -0,0 +1,395 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2009 Benoit Jacob +// +// 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 . + +#ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H +#define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H + +/** \ingroup QR_Module + * \nonstableyet + * + * \class FullPivotingHouseholderQR + * + * \brief Householder rank-revealing QR decomposition of a matrix + * + * \param MatrixType the type of the matrix of which we are computing the QR decomposition + * + * This class performs a rank-revealing QR decomposition using Householder transformations. + * + * This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal + * numerical stability. + * + * \sa MatrixBase::householderRrqr() + */ +template class FullPivotingHouseholderQR +{ + public: + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) + }; + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix MatrixQType; + typedef Matrix HCoeffsType; + typedef Matrix IntRowVectorType; + typedef Matrix IntColVectorType; + typedef Matrix RowVectorType; + typedef Matrix ColVectorType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&). + */ + FullPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + + FullPivotingHouseholderQR(const MatrixType& matrix) + : m_qr(matrix.rows(), matrix.cols()), + m_hCoeffs(std::min(matrix.rows(),matrix.cols())), + m_isInitialized(false) + { + compute(matrix); + } + + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which + * *this is the QR decomposition, if any exists. + * + * \param b the right-hand-side of the equation to solve. + * + * \param result a pointer to the vector/matrix in which to store the solution, if any exists. + * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). + * If no solution exists, *result is left with undefined coefficients. + * + * \note The case where b is a matrix is not yet implemented. Also, this + * code is space inefficient. + * + * Example: \include FullPivotingHouseholderQR_solve.cpp + * Output: \verbinclude FullPivotingHouseholderQR_solve.out + */ + template + bool solve(const MatrixBase& b, ResultType *result) const; + + MatrixType matrixQ(void) const; + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + */ + const MatrixType& matrixQR() const { return m_qr; } + + FullPivotingHouseholderQR& compute(const MatrixType& matrix); + + const IntRowVectorType& colsPermutation() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_cols_permutation; + } + + const IntColVectorType& rowsTranspositions() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_rows_transpositions; + } + + /** \returns the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \warning a determinant can be very big or small, so for matrices + * of large enough dimension, there is a risk of overflow/underflow. + * + * \sa MatrixBase::determinant() + */ + typename MatrixType::RealScalar absDeterminant() const; + + /** \returns the rank of the matrix of which *this is the QR decomposition. + * + * \note This is computed at the time of the construction of the QR decomposition. This + * method does not perform any further computation. + */ + inline int rank() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_rank; + } + + /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. + * + * \note Since the rank is computed at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline int dimensionOfKernel() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_qr.cols() - m_rank; + } + + /** \returns true if the matrix of which *this is the QR decomposition represents an injective + * linear map, i.e. has trivial kernel; false otherwise. + * + * \note Since the rank is computed at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline bool isInjective() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_rank == m_qr.cols(); + } + + /** \returns true if the matrix of which *this is the QR decomposition represents a surjective + * linear map; false otherwise. + * + * \note Since the rank is computed at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline bool isSurjective() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_rank == m_qr.rows(); + } + + /** \returns true if the matrix of which *this is the QR decomposition is invertible. + * + * \note Since the rank is computed at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline bool isInvertible() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return isInjective() && isSurjective(); + } + + /** Computes the inverse of the matrix of which *this is the QR decomposition. + * + * \param result a pointer to the matrix into which to store the inverse. Resized if needed. + * + * \note If this matrix is not invertible, *result is left with undefined coefficients. + * Use isInvertible() to first determine whether this matrix is invertible. + * + * \sa inverse() + */ + inline void computeInverse(MatrixType *result) const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!"); + solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result); + } + + /** \returns the inverse of the matrix of which *this is the QR decomposition. + * + * \note If this matrix is not invertible, the returned matrix has undefined coefficients. + * Use isInvertible() to first determine whether this matrix is invertible. + * + * \sa computeInverse() + */ + inline MatrixType inverse() const + { + MatrixType result; + computeInverse(&result); + return result; + } + + protected: + MatrixType m_qr; + HCoeffsType m_hCoeffs; + IntColVectorType m_rows_transpositions; + IntRowVectorType m_cols_permutation; + bool m_isInitialized; + RealScalar m_precision; + int m_rank; + int m_det_pq; +}; + +#ifndef EIGEN_HIDE_HEAVY_CODE + +template +typename MatrixType::RealScalar FullPivotingHouseholderQR::absDeterminant() const +{ + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); + return ei_abs(m_qr.diagonal().prod()); +} + +template +FullPivotingHouseholderQR& FullPivotingHouseholderQR::compute(const MatrixType& matrix) +{ + int rows = matrix.rows(); + int cols = matrix.cols(); + int size = std::min(rows,cols); + m_rank = size; + + m_qr = matrix; + m_hCoeffs.resize(size); + + RowVectorType temp(cols); + + m_precision = epsilon() * size; + + m_rows_transpositions.resize(matrix.rows()); + IntRowVectorType cols_transpositions(matrix.cols()); + m_cols_permutation.resize(matrix.cols()); + int number_of_transpositions = 0; + + RealScalar biggest; + + for (int k = 0; k < size; ++k) + { + int row_of_biggest_in_corner, col_of_biggest_in_corner; + RealScalar biggest_in_corner; + + biggest_in_corner = m_qr.corner(Eigen::BottomRight, rows-k, cols-k) + .cwise().abs() + .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); + row_of_biggest_in_corner += k; + col_of_biggest_in_corner += k; + if(k==0) biggest = biggest_in_corner; + + // if the corner is negligible, then we have less than full rank, and we can finish early + if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) + { + m_rank = k; + for(int i = k; i < size; i++) + { + m_rows_transpositions.coeffRef(i) = i; + cols_transpositions.coeffRef(i) = i; + m_hCoeffs.coeffRef(i) = Scalar(0); + } + break; + } + + m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner; + cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; + if(k != row_of_biggest_in_corner) { + m_qr.row(k).end(cols-k).swap(m_qr.row(row_of_biggest_in_corner).end(cols-k)); + ++number_of_transpositions; + } + if(k != col_of_biggest_in_corner) { + m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner)); + ++number_of_transpositions; + } + + RealScalar beta; + m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); + m_qr.coeffRef(k,k) = beta; + + m_qr.corner(BottomRight, rows-k, cols-k-1) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + } + + for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; + for(int k = 0; k < size; ++k) + std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); + + m_det_pq = (number_of_transpositions%2) ? -1 : 1; + m_isInitialized = true; + + return *this; +} + +template +template +bool FullPivotingHouseholderQR::solve( + const MatrixBase& b, + ResultType *result +) const +{ + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + if(m_rank==0) return false; + + const int rows = m_qr.rows(); + const int cols = b.cols(); + ei_assert(b.rows() == rows); + + typename OtherDerived::PlainMatrixType c(b); + + Matrix temp(cols); + for (int k = 0; k < m_rank; ++k) + { + int remainingSize = rows-k; + c.row(k).swap(c.row(m_rows_transpositions.coeff(k))); + c.corner(BottomRight, remainingSize, cols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); + } + + if(!isSurjective()) + { + // is c is in the image of R ? + RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff(); + if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision)) + return false; + } + m_qr.corner(TopLeft, m_rank, m_rank) + .template triangularView() + .solveInPlace(c.corner(TopLeft, m_rank, c.cols())); + + result->resize(m_qr.cols(), b.cols()); + for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i); + for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero(); + return true; +} + +/** \returns the matrix Q */ +template +MatrixType FullPivotingHouseholderQR::matrixQ() const +{ + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + // compute the product H'_0 H'_1 ... H'_n-1, + // where H_k is the k-th Householder transformation I - h_k v_k v_k' + // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] + int rows = m_qr.rows(); + int cols = m_qr.cols(); + int size = std::min(rows,cols); + MatrixType res = MatrixType::Identity(rows, rows); + Matrix temp(rows); + for (int k = size-1; k >= 0; k--) + { + res.block(k, k, rows-k, rows-k) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); + res.row(k).swap(res.row(m_rows_transpositions.coeff(k))); + } + return res; +} + +#endif // EIGEN_HIDE_HEAVY_CODE + +/** \return the full-pivoting Householder QR decomposition of \c *this. + * + * \sa class FullPivotingHouseholderQR + */ +template +const FullPivotingHouseholderQR::PlainMatrixType> +MatrixBase::fullPivotingHouseholderQr() const +{ + return FullPivotingHouseholderQR(eval()); +} + +#endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 90e6f8132..e5da6d691 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2008-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -39,7 +39,7 @@ * * Note that no pivoting is performed. This is \b not a rank-revealing decomposition. * - * \sa MatrixBase::qr() + * \sa MatrixBase::householderQr() */ template class HouseholderQR { @@ -54,6 +54,7 @@ template class HouseholderQR typedef Block MatrixRBlockType; typedef Matrix MatrixTypeR; typedef Matrix HCoeffsType; + typedef Matrix RowVectorType; /** * \brief Default Constructor. @@ -125,12 +126,12 @@ HouseholderQR& HouseholderQR::compute(const MatrixType& m_qr = matrix; m_hCoeffs.resize(size); - Matrix temp(cols); + RowVectorType temp(cols); for (int k = 0; k < size; ++k) { int remainingRows = rows - k; - int remainingCols = cols - k -1; + int remainingCols = cols - k - 1; RealScalar beta; m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake index ef36d2067..571774787 100644 --- a/cmake/EigenTesting.cmake +++ b/cmake/EigenTesting.cmake @@ -74,7 +74,7 @@ macro(ei_add_test testname) string(STRIP "${ARGV2}" ARGV2_stripped) string(LENGTH "${ARGV2_stripped}" ARGV2_stripped_length) if(${ARGV2_stripped_length} GREATER 0) - target_link_libraries(${targetname} "${ARGV2}") + target_link_libraries(${targetname} ${ARGV2}) endif(${ARGV2_stripped_length} GREATER 0) endif(${ARGC} GREATER 2) diff --git a/doc/C01_QuickStartGuide.dox b/doc/C01_QuickStartGuide.dox index 3d98e14f5..d43dbd72d 100644 --- a/doc/C01_QuickStartGuide.dox +++ b/doc/C01_QuickStartGuide.dox @@ -1,6 +1,6 @@ namespace Eigen { -/** \page TutorialCore Tutorial 1/3 - Core features +/** \page TutorialCore Tutorial 1/4 - Core features \ingroup Tutorial
\ref index "Overview" diff --git a/doc/C03_TutorialGeometry.dox b/doc/C03_TutorialGeometry.dox index 2d3d742c9..e349c68ce 100644 --- a/doc/C03_TutorialGeometry.dox +++ b/doc/C03_TutorialGeometry.dox @@ -1,6 +1,6 @@ namespace Eigen { -/** \page TutorialGeometry Tutorial 2/3 - Geometry +/** \page TutorialGeometry Tutorial 2/4 - Geometry \ingroup Tutorial
\ref index "Overview" diff --git a/doc/C05_TutorialLinearAlgebra.dox b/doc/C05_TutorialLinearAlgebra.dox index 46333cb58..e70298b47 100644 --- a/doc/C05_TutorialLinearAlgebra.dox +++ b/doc/C05_TutorialLinearAlgebra.dox @@ -1,7 +1,6 @@ - namespace Eigen { -/** \page TutorialAdvancedLinearAlgebra Tutorial 3/3 - Advanced linear algebra +/** \page TutorialAdvancedLinearAlgebra Tutorial 3/4 - Advanced linear algebra \ingroup Tutorial
\ref index "Overview" @@ -11,6 +10,9 @@ namespace Eigen { | \ref TutorialSparse "Sparse matrix"
+This tutorial chapter explains how you can use Eigen to tackle various problems involving matrices: +solving systems of linear equations, finding eigenvalues and eigenvectors, and so on. + \b Table \b of \b contents - \ref TutorialAdvSolvers - \ref TutorialAdvLU @@ -18,53 +20,129 @@ namespace Eigen { - \ref TutorialAdvQR - \ref TutorialAdvEigenProblems + \section TutorialAdvSolvers Solving linear problems -This part of the tutorial focuses on solving linear problem of the form \f$ A \mathbf{x} = b \f$, -where both \f$ A \f$ and \f$ b \f$ are known, and \f$ x \f$ is the unknown. Moreover, \f$ A \f$ -assumed to be a squared matrix. Of course, the methods described here can be used whenever an expression -involve the product of an inverse matrix with a vector or another matrix: \f$ A^{-1} B \f$. -Eigen offers various algorithms to this problem, and its choice mainly depends on the nature of -the matrix \f$ A \f$, such as its shape, size and numerical properties. +This part of the tutorial focuses on solving systems of linear equations. Such statems can be +written in the form \f$ A \mathbf{x} = \mathbf{b} \f$, where both \f$ A \f$ and \f$ \mathbf{b} \f$ +are known, and \f$ \mathbf{x} \f$ is the unknown. Moreover, \f$ A \f$ is assumed to be a square +matrix. + +The equation \f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution if \f$ A \f$ is invertible. If +the matrix is not invertible, then the equation may have no or infinitely many solutions. All +solvers assume that \f$ A \f$ is invertible, unless noted otherwise. + +Eigen offers various algorithms to this problem. The choice of algorithm mainly depends on the +nature of the matrix \f$ A \f$, such as its shape, size and numerical properties. + - The \ref TutorialAdvSolvers_LU "LU decomposition" (with partial pivoting) is a general-purpose + algorithm which works for most problems. + - Use the \ref TutorialAdvSolvers_Cholesky "Cholesky decomposition" if the matrix \f$ A \f$ is + positive definite. + - Use a special \ref TutorialAdvSolvers_Triangular "triangular solver" if the matrix \f$ A \f$ is + upper or lower triangular. + - Use of the \ref TutorialAdvSolvers_Inverse "matrix inverse" is not recommended in general, but + may be appropriate in special cases, for instance if you want to solve several systems with the + same matrix \f$ A \f$ and that matrix is small. + - \ref TutorialAdvSolvers_Misc "Other solvers" (%LU decomposition with full pivoting, the singular + value decomposition) are provided for special cases, such as when \f$ A \f$ is not invertible. + +The methods described here can be used whenever an expression involve the product of an inverse +matrix with a vector or another matrix: \f$ A^{-1} \mathbf{v} \f$ or \f$ A^{-1} B \f$. + + +\subsection TutorialAdvSolvers_LU LU decomposition (with partial pivoting) + +This is a general-purpose algorithm which performs well in most cases (provided the matrix \f$ A \f$ +is invertible), so if you are unsure about which algorithm to pick, choose this. The method proceeds +in two steps. First, the %LU decomposition with partial pivoting is computed using the +MatrixBase::partialLu() function. This yields an object of the class PartialLU. Then, the +PartialLU::solve() method is called to compute a solution. + +As an example, suppose we want to solve the following system of linear equations: + +\f[ \begin{aligned} + x + 2y + 3z &= 3 \\ + 4x + 5y + 6z &= 3 \\ + 7x + 8y + 10z &= 4. +\end{aligned} \f] + +The following program solves this system: -\subsection TutorialAdvSolvers_Triangular Triangular solver -If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the diagonal -are all not zero), then the problem can be solved directly using MatrixBase::solveTriangular(), or better, -MatrixBase::solveTriangularInPlace(). Here is an example: -
-\include MatrixBase_marked.cpp - -output: -\include MatrixBase_marked.out +\include Tutorial_PartialLU_solve.cpp + +output: \include Tutorial_PartialLU_solve.out
-See MatrixBase::solveTriangular() for more details. +There are many situations in which we want to solve the same system of equations with different +right-hand sides. One possibility is to put the right-hand sides as columns in a matrix \f$ B \f$ +and then solve the equation \f$ A X = B \f$. For instance, suppose that we want to solve the same +system as before, but now we also need the solution of the same equations with 1 on the right-hand +side. The following code computes the required solutions: + +
+\include Tutorial_solve_multiple_rhs.cpp + +output: \include Tutorial_solve_multiple_rhs.out +
+ +However, this is not always possible. Often, you only know the right-hand side of the second +problem, and whether you want to solve it at all, after you solved the first problem. In such a +case, it's best to save the %LU decomposition and reuse it to solve the second problem. This is +worth the effort because computing the %LU decomposition is much more expensive than using it to +solve the equation. Here is some code to illustrate the procedure. It uses the constructor +PartialLU::PartialLU(const MatrixType&) to compute the %LU decomposition. + +
+\include Tutorial_solve_reuse_decomposition.cpp + +output: \include Tutorial_solve_reuse_decomposition.out +
+ +\b Warning: All this code presumes that the matrix \f$ A \f$ is invertible, so that the system +\f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution. If the matrix \f$ A \f$ is not invertible, +then the system \f$ A \mathbf{x} = \mathbf{b} \f$ has either zero or infinitely many solutions. In +both cases, PartialLU::solve() will give nonsense results. For example, suppose that we want to +solve the same system as above, but with the 10 in the last equation replaced by 9. Then the system +of equations is inconsistent: adding the first and the third equation gives \f$ 8x + 10y + 12z = 7 \f$, +which implies \f$ 4x + 5y + 6z = 3\frac12 \f$, in contradiction with the seocond equation. If we try +to solve this inconsistent system with Eigen, we find: + +
+\include Tutorial_solve_singular.cpp + +output: \include Tutorial_solve_singular.out +
+ +The %LU decomposition with \b full pivoting (class LU) and the singular value decomposition (class +SVD) may be helpful in this case, as explained in the section \ref TutorialAdvSolvers_Misc below. + +\sa LU_Module, MatrixBase::partialLu(), PartialLU::solve(), class PartialLU. -\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices) -If the matrix \f$ A \f$ is small (\f$ \leq 4 \f$) and invertible, then a good approach is to directly compute -the inverse of the matrix \f$ A \f$, and then obtain the solution \f$ x \f$ by \f$ \mathbf{x} = A^{-1} b \f$. With Eigen, -this can be implemented like this: +\subsection TutorialAdvSolvers_Cholesky Cholesky decomposition -\code -#include -Matrix4f A = Matrix4f::Random(); -Vector4f b = Vector4f::Random(); -Vector4f x = A.inverse() * b; -\endcode +If the matrix \f$ A \f$ is \b symmetric \b positive \b definite, then the best method is to use a +Cholesky decomposition: it is both faster and more accurate than the %LU decomposition. Such +positive definite matrices often arise when solving overdetermined problems. These are linear +systems \f$ A \mathbf{x} = \mathbf{b} \f$ in which the matrix \f$ A \f$ has more rows than columns, +so that there are more equations than unknowns. Typically, there is no vector \f$ \mathbf{x} \f$ +which satisfies all the equation. Instead, we look for the least-square solution, that is, the +vector \f$ \mathbf{x} \f$ for which \f$ \| A \mathbf{x} - \mathbf{b} \|_2 \f$ is minimal. You can +find this vector by solving the equation \f$ A^T \! A \mathbf{x} = A^T \mathbf{b} \f$. If the matrix +\f$ A \f$ has full rank, then \f$ A^T \! A \f$ is positive definite and thus you can use the +Cholesky decomposition to solve this system and find the least-square solution to the original +system \f$ A \mathbf{x} = \mathbf{b} \f$. -Note that the function inverse() is defined in the LU module. -See MatrixBase::inverse() for more details. +Eigen offers two different Cholesky decompositions: the LLT class provides a \f$ LL^T \f$ +decomposition where L is a lower triangular matrix, and the LDLT class provides a \f$ LDL^T \f$ +decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. The latter +includes pivoting and avoids square roots; this makes the %LDLT decomposition slightly more stable +than the %LLT decomposition. The LDLT class is able to handle both positive- and negative-definite +matrices, but not indefinite matrices. +The API is the same as when using the %LU decomposition. -\subsection TutorialAdvSolvers_Symmetric Cholesky (for positive definite matrices) -If the matrix \f$ A \f$ is \b positive \b definite, then -the best method is to use a Cholesky decomposition. -Such positive definite matrices often arise when solving overdetermined problems in a least square sense (see below). -Eigen offers two different Cholesky decompositions: a \f$ LL^T \f$ decomposition where L is a lower triangular matrix, -and a \f$ LDL^T \f$ decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. -The latter avoids square roots and is therefore slightly more stable than the former one. \code #include MatrixXf D = MatrixXf::Random(8,4); @@ -74,15 +152,19 @@ VectorXf x; A.llt().solve(b,&x); // using a LLT factorization A.ldlt().solve(b,&x); // using a LDLT factorization \endcode -when the value of the right hand side \f$ b \f$ is not needed anymore, then it is faster to use -the \em in \em place API, e.g.: + +The LLT and LDLT classes also provide an \em in \em place API for the case where the value of the +right hand-side \f$ b \f$ is not needed anymore. + \code A.llt().solveInPlace(b); \endcode -In that case the value of \f$ b \f$ is replaced by the result \f$ x \f$. -If the linear problem has to solved for various right hand side \f$ b_i \f$, but same matrix \f$ A \f$, -then it is highly recommended to perform the decomposition of \$ A \$ only once, e.g.: +This code replaces the vector \f$ b \f$ by the result \f$ x \f$. + +As before, you can reuse the factorization if you have to solve the same linear problem with +different right-hand sides, e.g.: + \code // ... LLT lltOfA(A); @@ -91,40 +173,69 @@ lltOfA.solveInPlace(b1); // ... \endcode -\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT. +\sa Cholesky_Module, MatrixBase::llt(), MatrixBase::ldlt(), LLT::solve(), LLT::solveInPlace(), +LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT. -\subsection TutorialAdvSolvers_LU LU decomposition (for most cases) -If the matrix \f$ A \f$ does not fit in any of the previous categories, or if you are unsure about the numerical -stability of your problem, then you can use the LU solver based on a decomposition of the same name : see the section \ref TutorialAdvLU below. Actually, Eigen's LU module does not implement a standard LU decomposition, but rather a so-called LU decomposition -with full pivoting and rank update which has much better numerical stability. -The API of the LU solver is the same than the Cholesky one, except that there is no \em in \em place variant: -\code -#include -MatrixXf A = MatrixXf::Random(20,20); -VectorXf b = VectorXf::Random(20); -VectorXf x; -A.lu().solve(b, &x); -\endcode +\subsection TutorialAdvSolvers_Triangular Triangular solver -Again, the LU decomposition can be stored to be reused or to perform other kernel operations: -\code -// ... -LU luOfA(A); -luOfA.solve(b, &x); -// ... -\endcode +If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the +diagonal are all not zero), then the problem can be solved directly using the TriangularView +class. This is much faster than using an %LU or Cholesky decomposition (in fact, the triangular +solver is used when you solve a system using the %LU or Cholesky decomposition). Here is an example: -See the section \ref TutorialAdvLU below. +
+\include Tutorial_solve_triangular.cpp + +output: \include Tutorial_solve_triangular.out +
-\sa class LU, LU::solve(), LU_Module +The MatrixBase::triangularView() function constructs an object of the class TriangularView, and +TriangularView::solve() then solves the system. There is also an \e in \e place variant: + +
+\include Tutorial_solve_triangular_inplace.cpp + +output: \include Tutorial_solve_triangular_inplace.out +
+ +\sa MatrixBase::triangularView(), TriangularView::solve(), TriangularView::solveInPlace(), +TriangularView class. -\subsection TutorialAdvSolvers_SVD SVD solver (for singular matrices and special cases) -Finally, Eigen also offer a solver based on a singular value decomposition (SVD). Again, the API is the -same than with Cholesky or LU: +\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices) + +The solution of the system \f$ A \mathbf{x} = \mathbf{b} \f$ is given by \f$ \mathbf{x} = A^{-1} +\mathbf{b} \f$. This suggests the following approach for solving the system: compute the matrix +inverse and multiply that with the right-hand side. This is often not a good approach: using the %LU +decomposition with partial pivoting yields a more accurate algorithm that is usually just as fast or +even faster. However, using the matrix inverse can be faster if the matrix \f$ A \f$ is small +(≤4) and fixed size, though numerical stability problems may still remain. Here is an example of +how you would write this in Eigen: + +
+\include Tutorial_solve_matrix_inverse.cpp + +output: \include Tutorial_solve_matrix_inverse.out +
+ +Note that the function inverse() is defined in the \ref LU_Module. + +\sa MatrixBase::inverse(). + + +\subsection TutorialAdvSolvers_Misc Other solvers (for singular matrices and special cases) + +Finally, Eigen also offer solvers based on a singular value decomposition (%SVD) or the %LU +decomposition with full pivoting. These have the same API as the solvers based on the %LU +decomposition with partial pivoting (PartialLU). + +The solver based on the %SVD uses the class SVD. It can handle singular matrices. Here is an example +of its use: + \code #include +// ... MatrixXf A = MatrixXf::Random(20,20); VectorXf b = VectorXf::Random(20); VectorXf x; @@ -133,8 +244,23 @@ SVD svdOfA(A); svdOfA.solve(b, &x); \endcode -\sa class SVD, SVD::solve(), SVD_Module +%LU decomposition with full pivoting has better numerical stability than %LU decomposition with +partial pivoting. It is defined in the class LU. The solver can also handle singular matrices. +\code +#include +// ... +MatrixXf A = MatrixXf::Random(20,20); +VectorXf b = VectorXf::Random(20); +VectorXf x; +A.lu().solve(b, &x); +LU luOfA(A); +luOfA.solve(b, &x); +\endcode + +See the section \ref TutorialAdvLU below. + +\sa class SVD, SVD::solve(), SVD_Module, class LU, LU::solve(), LU_Module. diff --git a/doc/C07_TutorialSparse.dox b/doc/C07_TutorialSparse.dox index a8bfe006e..3a7182883 100644 --- a/doc/C07_TutorialSparse.dox +++ b/doc/C07_TutorialSparse.dox @@ -1,6 +1,6 @@ namespace Eigen { -/** \page TutorialSparse Tutorial - Getting started with the sparse module +/** \page TutorialSparse Tutorial 4/4 - Getting started with the sparse module \ingroup Tutorial
\ref index "Overview" diff --git a/doc/examples/Tutorial_PartialLU_solve.cpp b/doc/examples/Tutorial_PartialLU_solve.cpp new file mode 100644 index 000000000..80c393f9a --- /dev/null +++ b/doc/examples/Tutorial_PartialLU_solve.cpp @@ -0,0 +1,18 @@ +#include +#include + +using namespace std; +using namespace Eigen; + +int main(int, char *[]) +{ + Matrix3f A; + Vector3f b; + A << 1,2,3, 4,5,6, 7,8,10; + b << 3, 3, 4; + cout << "Here is the matrix A:" << endl << A << endl; + cout << "Here is the vector b:" << endl << b << endl; + Vector3f x; + A.partialLu().solve(b, &x); + cout << "The solution is:" << endl << x << endl; +} diff --git a/doc/snippets/Tutorial_solve_matrix_inverse.cpp b/doc/snippets/Tutorial_solve_matrix_inverse.cpp new file mode 100644 index 000000000..fff324446 --- /dev/null +++ b/doc/snippets/Tutorial_solve_matrix_inverse.cpp @@ -0,0 +1,6 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 4,5,6, 7,8,10; +b << 3, 3, 4; +Vector3f x = A.inverse() * b; +cout << "The solution is:" << endl << x << endl; diff --git a/doc/snippets/Tutorial_solve_multiple_rhs.cpp b/doc/snippets/Tutorial_solve_multiple_rhs.cpp new file mode 100644 index 000000000..fbb15165a --- /dev/null +++ b/doc/snippets/Tutorial_solve_multiple_rhs.cpp @@ -0,0 +1,10 @@ +Matrix3f A(3,3); +A << 1,2,3, 4,5,6, 7,8,10; +Matrix B; +B << 3,1, 3,1, 4,1; +Matrix X; +A.partialLu().solve(B, &X); +cout << "The solution with right-hand side (3,3,4) is:" << endl; +cout << X.col(0) << endl; +cout << "The solution with right-hand side (1,1,1) is:" << endl; +cout << X.col(1) << endl; diff --git a/doc/snippets/Tutorial_solve_reuse_decomposition.cpp b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp new file mode 100644 index 000000000..b4112adc4 --- /dev/null +++ b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp @@ -0,0 +1,13 @@ +Matrix3f A(3,3); +A << 1,2,3, 4,5,6, 7,8,10; +PartialLU luOfA(A); // compute LU decomposition of A +Vector3f b; +b << 3,3,4; +Vector3f x; +luOfA.solve(b, &x); +cout << "The solution with right-hand side (3,3,4) is:" << endl; +cout << x << endl; +b << 1,1,1; +luOfA.solve(b, &x); +cout << "The solution with right-hand side (1,1,1) is:" << endl; +cout << x << endl; diff --git a/doc/snippets/Tutorial_solve_singular.cpp b/doc/snippets/Tutorial_solve_singular.cpp new file mode 100644 index 000000000..da94ad445 --- /dev/null +++ b/doc/snippets/Tutorial_solve_singular.cpp @@ -0,0 +1,9 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 4,5,6, 7,8,9; +b << 3, 3, 4; +cout << "Here is the matrix A:" << endl << A << endl; +cout << "Here is the vector b:" << endl << b << endl; +Vector3f x; +A.partialLu().solve(b, &x); +cout << "The solution is:" << endl << x << endl; diff --git a/doc/snippets/Tutorial_solve_triangular.cpp b/doc/snippets/Tutorial_solve_triangular.cpp new file mode 100644 index 000000000..ff876f687 --- /dev/null +++ b/doc/snippets/Tutorial_solve_triangular.cpp @@ -0,0 +1,8 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 0,5,6, 0,0,10; +b << 3, 3, 4; +cout << "Here is the matrix A:" << endl << A << endl; +cout << "Here is the vector b:" << endl << b << endl; +Vector3f x = A.triangularView().solve(b); +cout << "The solution is:" << endl << x << endl; diff --git a/doc/snippets/Tutorial_solve_triangular_inplace.cpp b/doc/snippets/Tutorial_solve_triangular_inplace.cpp new file mode 100644 index 000000000..64ae4eea1 --- /dev/null +++ b/doc/snippets/Tutorial_solve_triangular_inplace.cpp @@ -0,0 +1,6 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 0,5,6, 0,0,10; +b << 3, 3, 4; +A.triangularView().solveInPlace(b); +cout << "The solution is:" << endl << b << endl; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ac4eb7633..fc43bbb1d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -78,7 +78,6 @@ else(QT4_FOUND) ei_add_property(EIGEN_MISSING_BACKENDS "Qt4 support, ") endif(QT4_FOUND) - if(TEST_LIB) add_definitions("-DEIGEN_EXTERN_INSTANTIATIONS=1") endif(TEST_LIB) @@ -122,6 +121,8 @@ ei_add_test(lu ${EI_OFLAG}) ei_add_test(determinant) ei_add_test(inverse ${EI_OFLAG}) ei_add_test(qr) +ei_add_test(qr_colpivoting) +ei_add_test(qr_fullpivoting) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") ei_add_test(svd) @@ -137,11 +138,6 @@ ei_add_test(regression) ei_add_test(stdvector) ei_add_test(resize) if(QT4_FOUND) - if(QT_QTCORE_LIBRARY_DEBUG) - set(QT_QTCORE_LIBRARY ${QT_QTCORE_LIBRARY_DEBUG}) - else(QT_QTCORE_LIBRARY_DEBUG) - set(QT_QTCORE_LIBRARY ${QT_QTCORE_LIBRARY_RELEASE}) - endif(QT_QTCORE_LIBRARY_DEBUG) ei_add_test(qtvector " " "${QT_QTCORE_LIBRARY}") endif(QT4_FOUND) ei_add_test(sparse_vector) diff --git a/test/lu.cpp b/test/lu.cpp index 4ad92bb11..75680b96b 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -44,7 +44,7 @@ template void lu_non_invertible() VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); VERIFY(!lu.isInjective()); VERIFY(!lu.isInvertible()); - VERIFY(lu.isSurjective() == (lu.rank() == rows)); + VERIFY(!lu.isSurjective()); VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); VERIFY(m1image.lu().rank() == rank); MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); @@ -53,7 +53,7 @@ template void lu_non_invertible() m2 = MatrixType::Random(cols,cols2); m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); - lu.solve(m3, &m2); + VERIFY(lu.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); m3 = MatrixType::Random(rows,cols2); VERIFY(!lu.solve(m3, &m2)); diff --git a/test/main.h b/test/main.h index 9ad014b64..e3866c68b 100644 --- a/test/main.h +++ b/test/main.h @@ -227,12 +227,11 @@ inline bool test_ei_isMuchSmallerThan(const MatrixBase& m, return m.isMuchSmallerThan(s, test_precision::Scalar>()); } -template -void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::MatrixBase& m) +template +void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m) { - typedef Derived MatrixType; typedef typename ei_traits::Scalar Scalar; - typedef Matrix VectorType; + typedef Matrix VectorType; MatrixType a = MatrixType::Random(rows,rows); MatrixType d = MatrixType::Identity(rows,cols); @@ -244,7 +243,7 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri HouseholderQR qra(a); HouseholderQR qrb(b); - m = (qra.matrixQ() * d * qrb.matrixQ()).lazy(); + m = qra.matrixQ() * d * qrb.matrixQ(); } } // end namespace Eigen diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp new file mode 100644 index 000000000..3a07c7131 --- /dev/null +++ b/test/qr_colpivoting.cpp @@ -0,0 +1,116 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2009 Benoit Jacob +// +// 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 . + +#include "main.h" +#include + +template void qr() +{ + /* this test covers the following files: QR.h */ + int rows = ei_random(20,200), cols = ei_random(20,200), cols2 = ei_random(20,200); + int rank = ei_random(1, std::min(rows, cols)-1); + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix SquareMatrixType; + typedef Matrix VectorType; + MatrixType m1; + createRandomMatrixOfRank(rank,rows,cols,m1); + ColPivotingHouseholderQR qr(m1); + VERIFY_IS_APPROX(rank, qr.rank()); + + MatrixType r = qr.matrixQR(); + // FIXME need better way to construct trapezoid + for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); + + MatrixType b = qr.matrixQ() * r; + + MatrixType c = MatrixType::Zero(rows,cols); + + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); + VERIFY_IS_APPROX(m1, c); + + MatrixType m2 = MatrixType::Random(cols,cols2); + MatrixType m3 = m1*m2; + m2 = MatrixType::Random(cols,cols2); + qr.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); +} + +template void qr_invertible() +{ + /* this test covers the following files: RRQR.h */ + typedef typename NumTraits::Real RealScalar; + int size = ei_random(10,50); + + MatrixType m1(size, size), m2(size, size), m3(size, size); + m1 = MatrixType::Random(size,size); + + if (ei_is_same_type::ret) + { + // let's build a matrix more stable to inverse + MatrixType a = MatrixType::Random(size,size*2); + m1 += a * a.adjoint(); + } + + ColPivotingHouseholderQR qr(m1); + m3 = MatrixType::Random(size,size); + qr.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); +} + +template void qr_verify_assert() +{ + MatrixType tmp; + + ColPivotingHouseholderQR qr; + VERIFY_RAISES_ASSERT(qr.matrixR()) + VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(qr.matrixQ()) +} + +void test_qr_colpivoting() +{ + for(int i = 0; i < 1; i++) { + // FIXME : very weird bug here +// CALL_SUBTEST( qr(Matrix2f()) ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + } + + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + } + + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); +} diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp new file mode 100644 index 000000000..d784e0d43 --- /dev/null +++ b/test/qr_fullpivoting.cpp @@ -0,0 +1,137 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2009 Benoit Jacob +// +// 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 . + +#include "main.h" +#include + +template void qr() +{ + /* this test covers the following files: QR.h */ + int rows = ei_random(20,200), cols = ei_random(20,200), cols2 = ei_random(20,200); + int rank = ei_random(1, std::min(rows, cols)-1); + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix SquareMatrixType; + typedef Matrix VectorType; + MatrixType m1; + createRandomMatrixOfRank(rank,rows,cols,m1); + FullPivotingHouseholderQR qr(m1); + VERIFY_IS_APPROX(rank, qr.rank()); + VERIFY(cols - qr.rank() == qr.dimensionOfKernel()); + VERIFY(!qr.isInjective()); + VERIFY(!qr.isInvertible()); + VERIFY(!qr.isSurjective()); + + + MatrixType r = qr.matrixQR(); + // FIXME need better way to construct trapezoid + for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); + + MatrixType b = qr.matrixQ() * r; + + MatrixType c = MatrixType::Zero(rows,cols); + + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); + VERIFY_IS_APPROX(m1, c); + + MatrixType m2 = MatrixType::Random(cols,cols2); + MatrixType m3 = m1*m2; + m2 = MatrixType::Random(cols,cols2); + VERIFY(qr.solve(m3, &m2)); + VERIFY_IS_APPROX(m3, m1*m2); + m3 = MatrixType::Random(rows,cols2); + VERIFY(!qr.solve(m3, &m2)); +} + +template void qr_invertible() +{ + typedef typename NumTraits::Real RealScalar; + typedef typename MatrixType::Scalar Scalar; + + int size = ei_random(10,50); + + MatrixType m1(size, size), m2(size, size), m3(size, size); + m1 = MatrixType::Random(size,size); + + if (ei_is_same_type::ret) + { + // let's build a matrix more stable to inverse + MatrixType a = MatrixType::Random(size,size*2); + m1 += a * a.adjoint(); + } + + FullPivotingHouseholderQR qr(m1); + VERIFY(qr.isInjective()); + VERIFY(qr.isInvertible()); + VERIFY(qr.isSurjective()); + + m3 = MatrixType::Random(size,size); + VERIFY(qr.solve(m3, &m2)); + VERIFY_IS_APPROX(m3, m1*m2); + + // now construct a matrix with prescribed determinant + m1.setZero(); + for(int i = 0; i < size; i++) m1(i,i) = ei_random(); + RealScalar absdet = ei_abs(m1.diagonal().prod()); + m3 = qr.matrixQ(); // get a unitary + m1 = m3 * m1 * m3; + qr.compute(m1); + VERIFY_IS_APPROX(absdet, qr.absDeterminant()); +} + +template void qr_verify_assert() +{ + MatrixType tmp; + + FullPivotingHouseholderQR qr; + VERIFY_RAISES_ASSERT(qr.matrixR()) + VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(qr.matrixQ()) +} + +void test_qr_fullpivoting() +{ + for(int i = 0; i < 1; i++) { + // FIXME : very weird bug here +// CALL_SUBTEST( qr(Matrix2f()) ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + } + + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + } + + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); +} diff --git a/test/stdvector.cpp b/test/stdvector.cpp index 6779feae5..a976f0cf1 100644 --- a/test/stdvector.cpp +++ b/test/stdvector.cpp @@ -159,5 +159,5 @@ void test_stdvector() // some Quaternion CALL_SUBTEST(check_stdvector_quaternion(Quaternionf())); - CALL_SUBTEST(check_stdvector_quaternion(Quaternionf())); + CALL_SUBTEST(check_stdvector_quaternion(Quaterniond())); }