Home | History | Annotate | Download | only in SVD 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1 (at) gmail.com>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_BIDIAGONALIZATION_H
     11 #define EIGEN_BIDIAGONALIZATION_H
     12 
     13 namespace Eigen {
     14 
     15 namespace internal {
     16 // UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API.
     17 // At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class.
     18 
     19 template<typename _MatrixType> class UpperBidiagonalization
     20 {
     21   public:
     22 
     23     typedef _MatrixType MatrixType;
     24     enum {
     25       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
     26       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     27       ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
     28     };
     29     typedef typename MatrixType::Scalar Scalar;
     30     typedef typename MatrixType::RealScalar RealScalar;
     31     typedef typename MatrixType::Index Index;
     32     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
     33     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
     34     typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType;
     35     typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
     36     typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
     37     typedef HouseholderSequence<
     38               const MatrixType,
     39               CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, const Diagonal<const MatrixType,0> >
     40             > HouseholderUSequenceType;
     41     typedef HouseholderSequence<
     42               const MatrixType,
     43               Diagonal<const MatrixType,1>,
     44               OnTheRight
     45             > HouseholderVSequenceType;
     46 
     47     /**
     48     * \brief Default Constructor.
     49     *
     50     * The default constructor is useful in cases in which the user intends to
     51     * perform decompositions via Bidiagonalization::compute(const MatrixType&).
     52     */
     53     UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
     54 
     55     UpperBidiagonalization(const MatrixType& matrix)
     56       : m_householder(matrix.rows(), matrix.cols()),
     57         m_bidiagonal(matrix.cols(), matrix.cols()),
     58         m_isInitialized(false)
     59     {
     60       compute(matrix);
     61     }
     62 
     63     UpperBidiagonalization& compute(const MatrixType& matrix);
     64 
     65     const MatrixType& householder() const { return m_householder; }
     66     const BidiagonalType& bidiagonal() const { return m_bidiagonal; }
     67 
     68     const HouseholderUSequenceType householderU() const
     69     {
     70       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
     71       return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
     72     }
     73 
     74     const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
     75     {
     76       eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
     77       return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>())
     78              .setLength(m_householder.cols()-1)
     79              .setShift(1);
     80     }
     81 
     82   protected:
     83     MatrixType m_householder;
     84     BidiagonalType m_bidiagonal;
     85     bool m_isInitialized;
     86 };
     87 
     88 template<typename _MatrixType>
     89 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix)
     90 {
     91   Index rows = matrix.rows();
     92   Index cols = matrix.cols();
     93 
     94   eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols.");
     95 
     96   m_householder = matrix;
     97 
     98   ColVectorType temp(rows);
     99 
    100   for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k)
    101   {
    102     Index remainingRows = rows - k;
    103     Index remainingCols = cols - k - 1;
    104 
    105     // construct left householder transform in-place in m_householder
    106     m_householder.col(k).tail(remainingRows)
    107                  .makeHouseholderInPlace(m_householder.coeffRef(k,k),
    108                                          m_bidiagonal.template diagonal<0>().coeffRef(k));
    109     // apply householder transform to remaining part of m_householder on the left
    110     m_householder.bottomRightCorner(remainingRows, remainingCols)
    111                  .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1),
    112                                             m_householder.coeff(k,k),
    113                                             temp.data());
    114 
    115     if(k == cols-1) break;
    116 
    117     // construct right householder transform in-place in m_householder
    118     m_householder.row(k).tail(remainingCols)
    119                  .makeHouseholderInPlace(m_householder.coeffRef(k,k+1),
    120                                          m_bidiagonal.template diagonal<1>().coeffRef(k));
    121     // apply householder transform to remaining part of m_householder on the left
    122     m_householder.bottomRightCorner(remainingRows-1, remainingCols)
    123                  .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(),
    124                                              m_householder.coeff(k,k+1),
    125                                              temp.data());
    126   }
    127   m_isInitialized = true;
    128   return *this;
    129 }
    130 
    131 #if 0
    132 /** \return the Householder QR decomposition of \c *this.
    133   *
    134   * \sa class Bidiagonalization
    135   */
    136 template<typename Derived>
    137 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
    138 MatrixBase<Derived>::bidiagonalization() const
    139 {
    140   return UpperBidiagonalization<PlainObject>(eval());
    141 }
    142 #endif
    143 
    144 } // end namespace internal
    145 
    146 } // end namespace Eigen
    147 
    148 #endif // EIGEN_BIDIAGONALIZATION_H
    149