diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index a41a5f670..024590f3c 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -30,14 +30,32 @@ * * \class Tridiagonalization * - * \brief Trigiagonal decomposition of a selfadjoint matrix + * \brief Tridiagonal decomposition of a selfadjoint matrix * - * \param MatrixType the type of the matrix of which we are performing the tridiagonalization + * \tparam _MatrixType the type of the matrix of which we are computing the + * tridiagonal decomposition; this is expected to be an instantiation of the + * Matrix class template. * * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that: * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix. * - * \sa MatrixBase::tridiagonalize() + * A tridiagonal matrix is a matrix which has nonzero elements only on the + * main diagonal and the first diagonal below and above it. The Hessenberg + * decomposition of a selfadjoint matrix is in fact a tridiagonal + * decomposition. This class is used in SelfAdjointEigenSolver to compute the + * eigenvalues and eigenvectors of a selfadjoint matrix. + * + * Call the function compute() to compute the tridiagonal decomposition of a + * given matrix. Alternatively, you can use the Tridiagonalization(const MatrixType&) + * constructor which computes the tridiagonal Schur decomposition at + * construction time. Once the decomposition is computed, you can use the + * matrixQ() and matrixT() functions to retrieve the matrices Q and T in the + * decomposition. + * + * The documentation of Tridiagonalization(const MatrixType&) contains an + * example of the typical use of this class. + * + * \sa class HessenbergDecomposition, class SelfAdjointEigenSolver */ template class Tridiagonalization { @@ -46,21 +64,18 @@ template class Tridiagonalization typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef typename ei_packet_traits::type Packet; enum { Size = MatrixType::RowsAtCompileTime, SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, Options = MatrixType::Options, MaxSize = MatrixType::MaxRowsAtCompileTime, - MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1, - PacketSize = ei_packet_traits::size + MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 }; typedef Matrix CoeffVectorType; typedef typename ei_plain_col_type::type DiagonalType; typedef Matrix SubDiagonalType; - typedef typename ei_plain_row_type::type RowVectorType; typedef typename ei_meta_if::IsComplex, typename Diagonal::RealReturnType, @@ -74,22 +89,53 @@ template class Tridiagonalization Block,0 > >::ret SubDiagonalReturnType; - /** This constructor initializes a Tridiagonalization object for - * further use with Tridiagonalization::compute() + /** \brief Default constructor. + * + * \param [in] size Positive integer, size of the matrix whose tridiagonal + * decomposition will be computed. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via compute(). The \p size parameter is only + * used as a hint. It is not an error to give a wrong \p size, but it may + * impair performance. + * + * \sa compute() for an example. */ Tridiagonalization(int size = Size==Dynamic ? 2 : Size) : m_matrix(size,size), m_hCoeffs(size-1) {} + /** \brief Constructor; computes tridiagonal decomposition of given matrix. + * + * \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition + * is to be computed. + * + * This constructor calls compute() to compute the tridiagonal decomposition. + * + * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp + * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out + */ Tridiagonalization(const MatrixType& matrix) : m_matrix(matrix), m_hCoeffs(matrix.cols()-1) { _compute(m_matrix, m_hCoeffs); } - /** Computes or re-compute the tridiagonalization for the matrix \a matrix. + /** \brief Computes tridiagonal decomposition of given matrix. + * + * \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition + * is to be computed. * - * This method allows to re-use the allocated data. + * The tridiagonal decomposition is computed by bringing the columns of + * the matrix successively in the required form using Householder + * reflections. The cost is \f$ 4n^3/3 \f$ flops, where \f$ n \f$ denotes + * the size of the given matrix. + * + * This method reuses of the allocated data in the Tridiagonalization + * object, if the size of the matrix does not change. + * + * Example: \include Tridiagonalization_compute.cpp + * Output: \verbinclude Tridiagonalization_compute.out */ void compute(const MatrixType& matrix) { @@ -98,74 +144,191 @@ template class Tridiagonalization _compute(m_matrix, m_hCoeffs); } - /** \returns the householder coefficients allowing to - * reconstruct the matrix Q from the packed data. + /** \brief Returns the Householder coefficients. * - * \sa packedMatrix() + * \returns a const reference to the vector of Householder coefficients + * + * \pre Either the constructor Tridiagonalization(const MatrixType&) or + * the member function compute(const MatrixType&) has been called before + * to compute the tridiagonal decomposition of a matrix. + * + * The Householder coefficients allow the reconstruction of the matrix + * \f$ Q \f$ in the tridiagonal decomposition from the packed data. + * + * Example: \include Tridiagonalization_householderCoefficients.cpp + * Output: \verbinclude Tridiagonalization_householderCoefficients.out + * + * \sa packedMatrix(), \ref Householder_Module "Householder module" */ - inline CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; } + inline CoeffVectorType householderCoefficients() const { return m_hCoeffs; } - /** \returns the internal result of the decomposition. + /** \brief Returns the internal representation of the decomposition + * + * \returns a const reference to a matrix with the internal representation + * of the decomposition. + * + * \pre Either the constructor Tridiagonalization(const MatrixType&) or + * the member function compute(const MatrixType&) has been called before + * to compute the tridiagonal decomposition of a matrix. * * The returned matrix contains the following information: - * - the strict upper part is equal to the input matrix A - * - the diagonal and lower sub-diagonal represent the tridiagonal symmetric matrix (real). - * - the rest of the lower part contains the Householder vectors that, combined with - * Householder coefficients returned by householderCoefficients(), - * allows to reconstruct the matrix Q as follow: - * Q = H_{N-1} ... H_1 H_0 - * where the matrices H are the Householder transformations: - * H_i = (I - h_i * v_i * v_i') - * where h_i == householderCoefficients()[i] and v_i is a Householder vector: - * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] + * - the strict upper triangular part is equal to the input matrix A. + * - the diagonal and lower sub-diagonal represent the real tridiagonal + * symmetric matrix T. + * - the rest of the lower part contains the Householder vectors that, + * combined with Householder coefficients returned by + * householderCoefficients(), allows to reconstruct the matrix Q as + * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$. + * Here, the matrices \f$ H_i \f$ are the Householder transformations + * \f$ H_i = (I - h_i v_i v_i^T) \f$ + * where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and + * \f$ v_i \f$ is the Householder vector defined by + * \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$ + * with M the matrix returned by this function. * * See LAPACK for further details on this packed storage. + * + * Example: \include Tridiagonalization_packedMatrix.cpp + * Output: \verbinclude Tridiagonalization_packedMatrix.out + * + * \sa householderCoefficients() */ - inline const MatrixType& packedMatrix(void) const { return m_matrix; } + inline const MatrixType& packedMatrix() const { return m_matrix; } + /** \brief Reconstructs the unitary matrix Q in the decomposition + * + * \returns the matrix Q + * + * \pre Either the constructor Tridiagonalization(const MatrixType&) or + * the member function compute(const MatrixType&) has been called before + * to compute the tridiagonal decomposition of a matrix. + * + * This function reconstructs the matrix Q from the Householder + * coefficients and the packed matrix stored internally. This + * reconstruction requires \f$ 4n^3 / 3 \f$ flops. + * + * \sa Tridiagonalization(const MatrixType&) for an example, + * matrixT(), matrixQInPlace() + */ MatrixType matrixQ() const; + + /** \brief Reconstructs the unitary matrix Q in the decomposition + * + * This is an in-place variant of matrixQ() which avoids the copy. + * This function will probably be deleted soon. + */ template void matrixQInPlace(MatrixBase* q) const; + + /** \brief Constructs the tridiagonal matrix T in the decomposition + * + * \returns the matrix T + * + * \pre Either the constructor Tridiagonalization(const MatrixType&) or + * the member function compute(const MatrixType&) has been called before + * to compute the tridiagonal decomposition of a matrix. + * + * This function copies the matrix T from internal data. The diagonal and + * subdiagonal of the packed matrix as returned by packedMatrix() + * represents the matrix T. It may sometimes be sufficient to directly use + * the packed matrix or the vector expressions returned by diagonal() + * and subDiagonal() instead of creating a new matrix with this function. + * + * \sa Tridiagonalization(const MatrixType&) for an example, + * matrixQ(), packedMatrix(), diagonal(), subDiagonal() + */ MatrixType matrixT() const; - const DiagonalReturnType diagonal(void) const; - const SubDiagonalReturnType subDiagonal(void) const; + /** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition. + * + * \returns expression representing the diagonal of T + * + * \pre Either the constructor Tridiagonalization(const MatrixType&) or + * the member function compute(const MatrixType&) has been called before + * to compute the tridiagonal decomposition of a matrix. + * + * Example: \include Tridiagonalization_diagonal.cpp + * Output: \verbinclude Tridiagonalization_diagonal.out + * + * \sa matrixT(), subDiagonal() + */ + const DiagonalReturnType diagonal() const; + + /** \brief Returns the subdiagonal of the tridiagonal matrix T in the decomposition. + * + * \returns expression representing the subdiagonal of T + * + * \pre Either the constructor Tridiagonalization(const MatrixType&) or + * the member function compute(const MatrixType&) has been called before + * to compute the tridiagonal decomposition of a matrix. + * + * \sa diagonal() for an example, matrixT() + */ + const SubDiagonalReturnType subDiagonal() const; + + /** \brief Performs a full decomposition in place + * + * \param[in,out] mat On input, the selfadjoint matrix whose tridiagonal + * decomposition is to be computed. On output, the orthogonal matrix Q + * in the decomposition if \p extractQ is true. + * \param[out] diag The diagonal of the tridiagonal matrix T in the + * decomposition. + * \param[out] subdiag The subdiagonal of the tridiagonal matrix T in + * the decomposition. + * \param[in] extractQ If true, the orthogonal matrix Q in the + * decomposition is computed and stored in \p mat. + * + * Compute the tridiagonal matrix of \p mat in place. The tridiagonal + * matrix T is passed to the output parameters \p diag and \p subdiag. If + * \p extractQ is true, then the orthogonal matrix Q is passed to \p mat. + * + * The vectors \p diag and \p subdiag are not resized. The function + * assumes that they are already of the correct size. The length of the + * vector \p diag should equal the number of rows in \p mat, and the + * length of the vector \p subdiag should be one left. + * + * This implementation contains an optimized path for real 3-by-3 matrices + * which is especially useful for plane fitting. + * + * \note Notwithstanding the name, the current implementation copies + * \p mat to a temporary matrix and uses that matrix to compute the + * decomposition. + * + * Example (this uses the same matrix as the example in + * Tridiagonalization(const MatrixType&)): + * \include Tridiagonalization_decomposeInPlace.cpp + * Output: \verbinclude Tridiagonalization_decomposeInPlace.out + * + * \sa Tridiagonalization(const MatrixType&), compute() + */ static void decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true); - static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); - protected: + static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); static void _decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true); MatrixType m_matrix; CoeffVectorType m_hCoeffs; }; -/** \returns an expression of the diagonal vector */ template const typename Tridiagonalization::DiagonalReturnType -Tridiagonalization::diagonal(void) const +Tridiagonalization::diagonal() const { return m_matrix.diagonal(); } -/** \returns an expression of the sub-diagonal vector */ template const typename Tridiagonalization::SubDiagonalReturnType -Tridiagonalization::subDiagonal(void) const +Tridiagonalization::subDiagonal() const { int n = m_matrix.rows(); return Block(m_matrix, 1, 0, n-1,n-1).diagonal(); } -/** constructs and returns the tridiagonal matrix T. - * Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix. - * Therefore, it might be often sufficient to directly use the packed matrix, or the vector - * expressions returned by diagonal() and subDiagonal() instead of creating a new matrix. - */ template typename Tridiagonalization::MatrixType -Tridiagonalization::matrixT(void) const +Tridiagonalization::matrixT() const { // FIXME should this function (and other similar ones) rather take a matrix as argument // and fill it ? (to avoid temporaries) @@ -223,10 +386,9 @@ void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& } } -/** reconstructs and returns the matrix Q */ template typename Tridiagonalization::MatrixType -Tridiagonalization::matrixQ(void) const +Tridiagonalization::matrixQ() const { MatrixType matQ; matrixQInPlace(&matQ); @@ -240,6 +402,7 @@ void Tridiagonalization::matrixQInPlace(MatrixBase* q) con QDerived& matQ = q->derived(); int n = m_matrix.rows(); matQ = MatrixType::Identity(n,n); + typedef typename ei_plain_row_type::type RowVectorType; RowVectorType aux(n); for (int i = n-2; i>=0; i--) { @@ -248,7 +411,6 @@ void Tridiagonalization::matrixQInPlace(MatrixBase* q) con } } -/** Performs a full decomposition in place */ template void Tridiagonalization::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { diff --git a/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp b/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp new file mode 100644 index 000000000..109650041 --- /dev/null +++ b/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp @@ -0,0 +1,11 @@ +MatrixXd X = MatrixXd::Random(5,5); +MatrixXd A = X + X.transpose(); +cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl; + +Tridiagonalization triOfA(A); +cout << "The orthogonal matrix Q is:" << endl << triOfA.matrixQ() << endl; +cout << "The tridiagonal matrix T is:" << endl << triOfA.matrixT() << endl << endl; + +MatrixXd Q = triOfA.matrixQ(); +MatrixXd T = triOfA.matrixT(); +cout << "Q * T * Q^T = " << endl << Q * T * Q.transpose() << endl; diff --git a/doc/snippets/Tridiagonalization_compute.cpp b/doc/snippets/Tridiagonalization_compute.cpp new file mode 100644 index 000000000..0062a99e8 --- /dev/null +++ b/doc/snippets/Tridiagonalization_compute.cpp @@ -0,0 +1,9 @@ +Tridiagonalization tri; +MatrixXf X = MatrixXf::Random(4,4); +MatrixXf A = X + X.transpose(); +tri.compute(A); +cout << "The matrix T in the tridiagonal decomposition of A is: " << endl; +cout << tri.matrixT() << endl; +tri.compute(2*A); // re-use tri to compute eigenvalues of 2A +cout << "The matrix T in the tridiagonal decomposition of 2A is: " << endl; +cout << tri.matrixT() << endl; diff --git a/doc/snippets/Tridiagonalization_decomposeInPlace.cpp b/doc/snippets/Tridiagonalization_decomposeInPlace.cpp new file mode 100644 index 000000000..2e53af7df --- /dev/null +++ b/doc/snippets/Tridiagonalization_decomposeInPlace.cpp @@ -0,0 +1,10 @@ +MatrixXd X = MatrixXd::Random(5,5); +MatrixXd A = X + X.transpose(); +cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl; + +VectorXd diag(5); +VectorXd subdiag(4); +Tridiagonalization::decomposeInPlace(A, diag, subdiag); +cout << "The orthogonal matrix Q is:" << endl << A << endl; +cout << "The diagonal of the tridiagonal matrix T is:" << endl << diag << endl; +cout << "The subdiagonal of the tridiagonal matrix T is:" << endl << subdiag << endl; diff --git a/doc/snippets/Tridiagonalization_diagonal.cpp b/doc/snippets/Tridiagonalization_diagonal.cpp new file mode 100644 index 000000000..30ce77d13 --- /dev/null +++ b/doc/snippets/Tridiagonalization_diagonal.cpp @@ -0,0 +1,13 @@ +MatrixXcd X = MatrixXcd::Random(4,4); +MatrixXcd A = X + X.adjoint(); +cout << "Here is a random self-adjoint 4x4 matrix:" << endl << A << endl << endl; + +Tridiagonalization triOfA(A); +MatrixXcd T = triOfA.matrixT(); +cout << "The tridiagonal matrix T is:" << endl << triOfA.matrixT() << endl << endl; + +cout << "We can also extract the diagonals of T directly ..." << endl; +VectorXd diag = triOfA.diagonal(); +cout << "The diagonal is:" << endl << diag << endl; +VectorXd subdiag = triOfA.subDiagonal(); +cout << "The subdiagonal is:" << endl << subdiag << endl; diff --git a/doc/snippets/Tridiagonalization_householderCoefficients.cpp b/doc/snippets/Tridiagonalization_householderCoefficients.cpp new file mode 100644 index 000000000..e5d872880 --- /dev/null +++ b/doc/snippets/Tridiagonalization_householderCoefficients.cpp @@ -0,0 +1,6 @@ +Matrix4d X = Matrix4d::Random(4,4); +Matrix4d A = X + X.transpose(); +cout << "Here is a random symmetric 4x4 matrix:" << endl << A << endl; +Tridiagonalization triOfA(A); +Vector3d hc = triOfA.householderCoefficients(); +cout << "The vector of Householder coefficients is:" << endl << hc << endl; diff --git a/doc/snippets/Tridiagonalization_packedMatrix.cpp b/doc/snippets/Tridiagonalization_packedMatrix.cpp new file mode 100644 index 000000000..0f55d0c28 --- /dev/null +++ b/doc/snippets/Tridiagonalization_packedMatrix.cpp @@ -0,0 +1,8 @@ +Matrix4d X = Matrix4d::Random(4,4); +Matrix4d A = X + X.transpose(); +cout << "Here is a random symmetric 4x4 matrix:" << endl << A << endl; +Tridiagonalization triOfA(A); +Matrix4d pm = triOfA.packedMatrix(); +cout << "The packed matrix M is:" << endl << pm << endl; +cout << "The diagonal and subdiagonal corresponds to the matrix T, which is:" + << endl << triOfA.matrixT() << endl;