aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen/src/Householder/HouseholderSequence.h
blob: 1e71e16a785d3e5269c6fff3a724c2aff05d1dc6 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
#define EIGEN_HOUSEHOLDER_SEQUENCE_H

namespace Eigen { 

/** \ingroup Householder_Module
  * \householder_module
  * \class HouseholderSequence
  * \brief Sequence of Householder reflections acting on subspaces with decreasing size
  * \tparam VectorsType type of matrix containing the Householder vectors
  * \tparam CoeffsType  type of vector containing the Householder coefficients
  * \tparam Side        either OnTheLeft (the default) or OnTheRight
  *
  * This class represents a product sequence of Householder reflections where the first Householder reflection
  * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by
  * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace
  * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but
  * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections
  * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods
  * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(),
  * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence.
  *
  * More precisely, the class %HouseholderSequence represents an \f$ n \times n \f$ matrix \f$ H \f$ of the
  * form \f$ H = \prod_{i=0}^{n-1} H_i \f$ where the i-th Householder reflection is \f$ H_i = I - h_i v_i
  * v_i^* \f$. The i-th Householder coefficient \f$ h_i \f$ is a scalar and the i-th Householder vector \f$
  * v_i \f$ is a vector of the form
  * \f[ 
  * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 
  * \f]
  * The last \f$ n-i \f$ entries of \f$ v_i \f$ are called the essential part of the Householder vector.
  *
  * Typical usages are listed below, where H is a HouseholderSequence:
  * \code
  * A.applyOnTheRight(H);             // A = A * H
  * A.applyOnTheLeft(H);              // A = H * A
  * A.applyOnTheRight(H.adjoint());   // A = A * H^*
  * A.applyOnTheLeft(H.adjoint());    // A = H^* * A
  * MatrixXd Q = H;                   // conversion to a dense matrix
  * \endcode
  * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators.
  *
  * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example.
  *
  * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
  */

namespace internal {

template<typename VectorsType, typename CoeffsType, int Side>
struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
{
  typedef typename VectorsType::Scalar Scalar;
  typedef typename VectorsType::Index Index;
  typedef typename VectorsType::StorageKind StorageKind;
  enum {
    RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
                                        : traits<VectorsType>::ColsAtCompileTime,
    ColsAtCompileTime = RowsAtCompileTime,
    MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
                                           : traits<VectorsType>::MaxColsAtCompileTime,
    MaxColsAtCompileTime = MaxRowsAtCompileTime,
    Flags = 0
  };
};

template<typename VectorsType, typename CoeffsType, int Side>
struct hseq_side_dependent_impl
{
  typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
  typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
  typedef typename VectorsType::Index Index;
  static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
  {
    Index start = k+1+h.m_shift;
    return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
  }
};

template<typename VectorsType, typename CoeffsType>
struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
{
  typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
  typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
  typedef typename VectorsType::Index Index;
  static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
  {
    Index start = k+1+h.m_shift;
    return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
  }
};

template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
{
  typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
    ResultScalar;
  typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
                 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
};

} // end namespace internal

template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
  : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
{
    enum {
      RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
      ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
      MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
      MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
    };
    typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
    typedef typename VectorsType::Index Index;

    typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType
            EssentialVectorType;

  public:

    typedef HouseholderSequence<
      VectorsType,
      typename internal::conditional<NumTraits<Scalar>::IsComplex,
        typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
        CoeffsType>::type,
      Side
    > ConjugateReturnType;

    /** \brief Constructor.
      * \param[in]  v      %Matrix containing the essential parts of the Householder vectors
      * \param[in]  h      Vector containing the Householder coefficients
      *
      * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The
      * i-th Householder coefficient \f$ h_i \f$ is given by \p h(i) and the essential part of the i-th
      * Householder vector \f$ v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the
      * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many
      * Householder reflections as there are columns.
      *
      * \note The %HouseholderSequence object stores \p v and \p h by reference.
      *
      * Example: \include HouseholderSequence_HouseholderSequence.cpp
      * Output: \verbinclude HouseholderSequence_HouseholderSequence.out
      *
      * \sa setLength(), setShift()
      */
    HouseholderSequence(const VectorsType& v, const CoeffsType& h)
      : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
        m_shift(0)
    {
    }

    /** \brief Copy constructor. */
    HouseholderSequence(const HouseholderSequence& other)
      : m_vectors(other.m_vectors),
        m_coeffs(other.m_coeffs),
        m_trans(other.m_trans),
        m_length(other.m_length),
        m_shift(other.m_shift)
    {
    }

    /** \brief Number of rows of transformation viewed as a matrix.
      * \returns Number of rows 
      * \details This equals the dimension of the space that the transformation acts on.
      */
    Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }

    /** \brief Number of columns of transformation viewed as a matrix.
      * \returns Number of columns
      * \details This equals the dimension of the space that the transformation acts on.
      */
    Index cols() const { return rows(); }

    /** \brief Essential part of a Householder vector.
      * \param[in]  k  Index of Householder reflection
      * \returns    Vector containing non-trivial entries of k-th Householder vector
      *
      * This function returns the essential part of the Householder vector \f$ v_i \f$. This is a vector of
      * length \f$ n-i \f$ containing the last \f$ n-i \f$ entries of the vector
      * \f[ 
      * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 
      * \f]
      * The index \f$ i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v
      * passed to the constructor.
      *
      * \sa setShift(), shift()
      */
    const EssentialVectorType essentialVector(Index k) const
    {
      eigen_assert(k >= 0 && k < m_length);
      return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
    }

    /** \brief %Transpose of the Householder sequence. */
    HouseholderSequence transpose() const
    {
      return HouseholderSequence(*this).setTrans(!m_trans);
    }

    /** \brief Complex conjugate of the Householder sequence. */
    ConjugateReturnType conjugate() const
    {
      return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
             .setTrans(m_trans)
             .setLength(m_length)
             .setShift(m_shift);
    }

    /** \brief Adjoint (conjugate transpose) of the Householder sequence. */
    ConjugateReturnType adjoint() const
    {
      return conjugate().setTrans(!m_trans);
    }

    /** \brief Inverse of the Householder sequence (equals the adjoint). */
    ConjugateReturnType inverse() const { return adjoint(); }

    /** \internal */
    template<typename DestType> inline void evalTo(DestType& dst) const
    {
      Matrix<Scalar, DestType::RowsAtCompileTime, 1,
             AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
      evalTo(dst, workspace);
    }

    /** \internal */
    template<typename Dest, typename Workspace>
    void evalTo(Dest& dst, Workspace& workspace) const
    {
      workspace.resize(rows());
      Index vecs = m_length;
      if(    internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
          && internal::extract_data(dst) == internal::extract_data(m_vectors))
      {
        // in-place
        dst.diagonal().setOnes();
        dst.template triangularView<StrictlyUpper>().setZero();
        for(Index k = vecs-1; k >= 0; --k)
        {
          Index cornerSize = rows() - k - m_shift;
          if(m_trans)
            dst.bottomRightCorner(cornerSize, cornerSize)
               .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
          else
            dst.bottomRightCorner(cornerSize, cornerSize)
               .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());

          // clear the off diagonal vector
          dst.col(k).tail(rows()-k-1).setZero();
        }
        // clear the remaining columns if needed
        for(Index k = 0; k<cols()-vecs ; ++k)
          dst.col(k).tail(rows()-k-1).setZero();
      }
      else
      {
        dst.setIdentity(rows(), rows());
        for(Index k = vecs-1; k >= 0; --k)
        {
          Index cornerSize = rows() - k - m_shift;
          if(m_trans)
            dst.bottomRightCorner(cornerSize, cornerSize)
               .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
          else
            dst.bottomRightCorner(cornerSize, cornerSize)
               .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
        }
      }
    }

    /** \internal */
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
    {
      Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
      applyThisOnTheRight(dst, workspace);
    }

    /** \internal */
    template<typename Dest, typename Workspace>
    inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
    {
      workspace.resize(dst.rows());
      for(Index k = 0; k < m_length; ++k)
      {
        Index actual_k = m_trans ? m_length-k-1 : k;
        dst.rightCols(rows()-m_shift-actual_k)
           .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
      }
    }

    /** \internal */
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
    {
      Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
      applyThisOnTheLeft(dst, workspace);
    }

    /** \internal */
    template<typename Dest, typename Workspace>
    inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
    {
      workspace.resize(dst.cols());
      for(Index k = 0; k < m_length; ++k)
      {
        Index actual_k = m_trans ? k : m_length-k-1;
        dst.bottomRows(rows()-m_shift-actual_k)
           .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
      }
    }

    /** \brief Computes the product of a Householder sequence with a matrix.
      * \param[in]  other  %Matrix being multiplied.
      * \returns    Expression object representing the product.
      *
      * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this
      * and \f$ M \f$ is the matrix \p other.
      */
    template<typename OtherDerived>
    typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
    {
      typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
        res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
      applyThisOnTheLeft(res);
      return res;
    }

    template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;

    /** \brief Sets the length of the Householder sequence.
      * \param [in]  length  New value for the length.
      *
      * By default, the length \f$ n \f$ of the Householder sequence \f$ H = H_0 H_1 \ldots H_{n-1} \f$ is set
      * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that
      * is smaller. After this function is called, the length equals \p length.
      *
      * \sa length()
      */
    HouseholderSequence& setLength(Index length)
    {
      m_length = length;
      return *this;
    }

    /** \brief Sets the shift of the Householder sequence.
      * \param [in]  shift  New value for the shift.
      *
      * By default, a %HouseholderSequence object represents \f$ H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th
      * column of the matrix \p v passed to the constructor corresponds to the i-th Householder
      * reflection. After this function is called, the object represents \f$ H = H_{\mathrm{shift}}
      * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th
      * Householder reflection.
      *
      * \sa shift()
      */
    HouseholderSequence& setShift(Index shift)
    {
      m_shift = shift;
      return *this;
    }

    Index length() const { return m_length; }  /**< \brief Returns the length of the Householder sequence. */
    Index shift() const { return m_shift; }    /**< \brief Returns the shift of the Householder sequence. */

    /* Necessary for .adjoint() and .conjugate() */
    template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;

  protected:

    /** \brief Sets the transpose flag.
      * \param [in]  trans  New value of the transpose flag.
      *
      * By default, the transpose flag is not set. If the transpose flag is set, then this object represents 
      * \f$ H^T = H_{n-1}^T \ldots H_1^T H_0^T \f$ instead of \f$ H = H_0 H_1 \ldots H_{n-1} \f$.
      *
      * \sa trans()
      */
    HouseholderSequence& setTrans(bool trans)
    {
      m_trans = trans;
      return *this;
    }

    bool trans() const { return m_trans; }     /**< \brief Returns the transpose flag. */

    typename VectorsType::Nested m_vectors;
    typename CoeffsType::Nested m_coeffs;
    bool m_trans;
    Index m_length;
    Index m_shift;
};

/** \brief Computes the product of a matrix with a Householder sequence.
  * \param[in]  other  %Matrix being multiplied.
  * \param[in]  h      %HouseholderSequence being multiplied.
  * \returns    Expression object representing the product.
  *
  * This function computes \f$ MH \f$ where \f$ M \f$ is the matrix \p other and \f$ H \f$ is the
  * Householder sequence represented by \p h.
  */
template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
{
  typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
    res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
  h.applyThisOnTheRight(res);
  return res;
}

/** \ingroup Householder_Module \householder_module
  * \brief Convenience function for constructing a Householder sequence. 
  * \returns A HouseholderSequence constructed from the specified arguments.
  */
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
{
  return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
}

/** \ingroup Householder_Module \householder_module
  * \brief Convenience function for constructing a Householder sequence. 
  * \returns A HouseholderSequence constructed from the specified arguments.
  * \details This function differs from householderSequence() in that the template argument \p OnTheSide of
  * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft.
  */
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
{
  return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
}

} // end namespace Eigen

#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H