| /* ----------------------------------------------------------------------- *//** |
| * |
| * @file SymmetricPositiveDefiniteEigenDecomposition_impl.hpp |
| * |
| *//* ----------------------------------------------------------------------- */ |
| |
| #ifndef MADLIB_DBAL_EIGEN_INTEGRATION_SPDEIGENDECOMPOSITION_IMPL_HPP |
| #define MADLIB_DBAL_EIGEN_INTEGRATION_SPDEIGENDECOMPOSITION_IMPL_HPP |
| |
| namespace madlib { |
| |
| namespace dbal { |
| |
| namespace eigen_integration { |
| |
| /** |
| * @brief Constructor that invokes the computation |
| * |
| * @param inMatrix Matrix to operate on. Note that the type is |
| * <tt>const MatrixType&</tt>, meaning that a temporary object will be |
| * created if the actual parameter is not of a subclass type. This means |
| * that memory will be copied -- on the positive side, this ensures memory |
| * will be aligned! |
| * @param inOptions A combination of DecompositionOptions |
| * @param inExtras A combination of SPDDecompositionExtras |
| */ |
| template <class MatrixType> |
| inline |
| SymmetricPositiveDefiniteEigenDecomposition<MatrixType> |
| ::SymmetricPositiveDefiniteEigenDecomposition( |
| const MatrixType &inMatrix, int inOptions, int inExtras) |
| : Base(inMatrix, inOptions) { |
| |
| computeExtras(inMatrix, inExtras); |
| } |
| |
| /** |
| * @brief Return the condition number of the matrix |
| * |
| * In general, the condition number of a matrix is the absolute value of the |
| * largest singular value divided by the smallest singular value. When a matrix |
| * is symmetric positive semi-definite, all eigenvalues are also singular |
| * values. Moreover, all eigenvalues are non-negative. |
| */ |
| template <class MatrixType> |
| inline |
| double |
| SymmetricPositiveDefiniteEigenDecomposition<MatrixType>::conditionNo() |
| const { |
| |
| const RealVectorType& ev = eigenvalues(); |
| |
| double numerator = ev(ev.size() - 1); |
| double denominator = ev(0); |
| |
| // All eigenvalues of a positive semi-definite matrix are |
| // non-negative, so in theory no need to take absolute values. |
| // Unfortunately, numerical instabilities can cause eigenvalues to |
| // be slightly negative. We should interpret that as 0. |
| if (denominator < 0) |
| denominator = 0; |
| |
| return numerator <= 0 ? std::numeric_limits<double>::infinity() |
| : numerator / denominator; |
| } |
| |
| /** |
| * @brief Return the pseudo inverse previously computed using computeExtras(). |
| * |
| * The result of this function is undefined if computeExtras() has not been |
| * called or the pseudo-inverse was not set to be computed. |
| */ |
| template <class MatrixType> |
| inline |
| const MatrixType& |
| SymmetricPositiveDefiniteEigenDecomposition<MatrixType>::pseudoInverse() const { |
| return mPinv; |
| } |
| |
| /** |
| * @brief Perform extra computations after the decomposition |
| * |
| * If the matrix has a condition number of less than 1e20 (currently |
| * this is hard-coded), it necessarily has full rank and is invertible. |
| * The Moore-Penrose pseudo-inverse coincides with the inverse and we |
| * compute it directly, using a \f$ L D L^T \f$ Cholesky decomposition. |
| * |
| * If the matrix has a condition number of more than 1e20, we are on the |
| * safe side and use the eigen decomposition for computing the |
| * pseudo-inverse. |
| * |
| * Since the eigenvectors of a symmtric positive semi-definite matrix |
| * are orthogonal, and Eigen moreover scales them to have norm 1 (i.e., |
| * the eigenvectors returned by Eigen are orthonormal), the Eigen |
| * decomposition |
| * |
| * \f$ M = V * D * V^T \f$ |
| * |
| * is also a singular value decomposition (where M is the |
| * original symmetric positive semi-definite matrix, D is the |
| * diagonal matrix with eigenvectors, and V is the unitary |
| * matrix containing normalized eigenvectors). In particular, |
| * V is unitary, so the inverse can be computed as |
| * |
| * \f$ M^{-1} = V * D^{-1} * V^T \f$. |
| * |
| * Only the <b>lower triangular part</b> of the input matrix |
| * is referenced. |
| */ |
| template <class MatrixType> |
| inline |
| void |
| SymmetricPositiveDefiniteEigenDecomposition<MatrixType>::computeExtras( |
| const MatrixType &inMatrix, int inExtras) { |
| |
| if (inExtras & ComputePseudoInverse) { |
| mPinv.resize(inMatrix.rows(), inMatrix.cols()); |
| |
| // FIXME: No hard-coded constant here |
| if (conditionNo() < 1e20) { |
| // We are doing a Cholesky decomposition of a matrix with |
| // pivoting. This is faster than the PartialPivLU that |
| // Eigen's inverse() method would use |
| mPinv = inMatrix.template selfadjointView<Eigen::Lower>().ldlt() |
| .solve(MatrixType::Identity(inMatrix.rows(), inMatrix.cols())); |
| } else { |
| if (!Base::m_eigenvectorsOk) |
| Base::compute(inMatrix, Eigen::ComputeEigenvectors); |
| |
| const RealVectorType& ev = eigenvalues(); |
| |
| // The eigenvalue are sorted in increasing order |
| Scalar epsilon = static_cast<double>(inMatrix.rows()) |
| * ev(ev.size() - 1) |
| * std::numeric_limits<Scalar>::epsilon(); |
| |
| RealVectorType eigenvectorsInverted(ev.size()); |
| for (Index i = 0; i < static_cast<Index>(ev.size()); ++i) { |
| eigenvectorsInverted(i) = ev(i) < epsilon |
| ? Scalar(0) |
| : Scalar(1) / ev(i); |
| } |
| mPinv = Base::eigenvectors() |
| * eigenvectorsInverted.asDiagonal() |
| * Base::eigenvectors().transpose(); |
| } |
| } |
| } |
| |
| } // namespace eigen_integration |
| |
| } // namespace dbal |
| |
| } // namespace madlib |
| |
| #endif // defined(MADLIB_DBAL_EIGEN_INTEGRATION_SPDEIGENDECOMPOSITION_IMPL_HPP) |