1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 #ifndef EIGEN_LLT_H
26 #define EIGEN_LLT_H
28 /** \ingroup cholesky_Module
29   *
30   * \class LLT
31   *
32   * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
33   *
34   * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
35   *
36   * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
37   * matrix A such that A = LL^* = U^*U, where L is lower triangular.
38   *
39   * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
40   * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
41   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
42   * situations like generalised eigen problems with hermitian matrices.
43   *
44   * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
45   * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
46   * has a solution.
47   *
48   * \sa MatrixBase::llt(), class LDLT
49   */
50  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
51   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
52   * the strict lower part does not have to store correct values.
53   */
54 template<typename MatrixType> class LLT
55 {
56   private:
57     typedef typename MatrixType::Scalar Scalar;
58     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
59     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
61     enum {
62       PacketSize = ei_packet_traits<Scalar>::size,
64     };
66   public:
68     /**
69     * \brief Default Constructor.
70     *
71     * The default constructor is useful in cases in which the user intends to
72     * perform decompositions via LLT::compute(const MatrixType&).
73     */
74     LLT() : m_matrix(), m_isInitialized(false) {}
76     LLT(const MatrixType& matrix)
77       : m_matrix(matrix.rows(), matrix.cols()),
78         m_isInitialized(false)
79     {
80       compute(matrix);
81     }
83     /** \returns the lower triangular matrix L */
84     inline Part<MatrixType, LowerTriangular> matrixL(void) const
85     {
86       ei_assert(m_isInitialized && "LLT is not initialized.");
87       return m_matrix;
88     }
90     /** \deprecated */
91     inline bool isPositiveDefinite(void) const { return m_isInitialized && m_isPositiveDefinite; }
93     template<typename RhsDerived, typename ResultType>
94     bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const;
96     template<typename Derived>
97     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
99     void compute(const MatrixType& matrix);
101   protected:
102     /** \internal
103       * Used to compute and store L
104       * The strict upper part is not used and even not initialized.
105       */
106     MatrixType m_matrix;
107     bool m_isInitialized;
108     bool m_isPositiveDefinite;
109 };
111 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
112   */
113 template<typename MatrixType>
114 void LLT<MatrixType>::compute(const MatrixType& a)
115 {
116   assert(a.rows()==a.cols());
117   m_isPositiveDefinite = true;
118   const int size = a.rows();
119   m_matrix.resize(size, size);
120   // The biggest overall is the point of reference to which further diagonals
121   // are compared; if any diagonal is negligible compared
122   // to the largest overall, the algorithm bails.  This cutoff is suggested
123   // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
124   // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
125   // Algorithms" page 217, also by Higham.
126   const RealScalar cutoff = machine_epsilon<Scalar>() * size * a.diagonal().cwise().abs().maxCoeff();
127   RealScalar x;
128   x = ei_real(a.coeff(0,0));
129   m_matrix.coeffRef(0,0) = ei_sqrt(x);
130   if(size==1)
131   {
132     m_isInitialized = true;
133     return;
134   }
135   m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
136   for (int j = 1; j < size; ++j)
137   {
138     x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
139     if (x < cutoff)
140     {
141       m_isPositiveDefinite = false;
142       continue;
143     }
145     m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
147     int endSize = size-j-1;
148     if (endSize>0) {
149       // Note that when all matrix columns have good alignment, then the following
150       // product is guaranteed to be optimal with respect to alignment.
151       m_matrix.col(j).end(endSize) =
152         (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy();
154       // FIXME could use a.col instead of a.row
156         - m_matrix.col(j).end(endSize) ) / x;
157     }
158   }
160   m_isInitialized = true;
161 }
163 /** Computes the solution x of \f$A x = b \f$ using the current decomposition of A.
164   * The result is stored in \a result
165   *
166   * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
167   *
168   * In other words, it computes \f$b = A^{-1} b \f$ with
169   * \f${L^{*}}^{-1} L^{-1} b \f$ from right to left.
170   *
171   * Example: \include LLT_solve.cpp
172   * Output: \verbinclude LLT_solve.out
173   *
174   * \sa LLT::solveInPlace(), MatrixBase::llt()
175   */
176 template<typename MatrixType>
177 template<typename RhsDerived, typename ResultType>
178 bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const
179 {
180   ei_assert(m_isInitialized && "LLT is not initialized.");
181   const int size = m_matrix.rows();
182   ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
183   return solveInPlace((*result) = b);
184 }
186 /** This is the \em in-place version of solve().
187   *
188   * \param bAndX represents both the right-hand side matrix b and result x.
189   *
190   * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
191   *
192   * This version avoids a copy when the right hand side matrix b is not
193   * needed anymore.
194   *
195   * \sa LLT::solve(), MatrixBase::llt()
196   */
197 template<typename MatrixType>
198 template<typename Derived>
199 bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
200 {
201   ei_assert(m_isInitialized && "LLT is not initialized.");
202   const int size = m_matrix.rows();
203   ei_assert(size==bAndX.rows());
204   matrixL().solveTriangularInPlace(bAndX);