Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "Matrix.hpp"
       4             : 
       5             : /** 
       6             :  * @brief   QR factorization using Householder reflectors.
       7             :  * 
       8             :  * Factorizes an m×n matrix with m >= n into an m×m unitary factor Q and an m×n 
       9             :  * upper triangular factor R.
      10             :  * 
      11             :  * It can be used for solving square systems of equations or for finding a least
      12             :  * squares solution to an overdetermined system of equations.
      13             :  * 
      14             :  * This version does not use column pivoting, and is not rank-revealing.
      15             :  * 
      16             :  * @ingroup Factorizations
      17             :  */
      18             : class HouseholderQR {
      19             :   public:
      20             :     /// @name Constructors
      21             :     /// @{
      22             : 
      23             :     /// Default constructor.
      24             :     HouseholderQR() = default;
      25             :     /// Factorize the given matrix.
      26          11 :     HouseholderQR(const Matrix &matrix) { compute(matrix); }
      27             :     /// Factorize the given matrix.
      28           1 :     HouseholderQR(Matrix &&matrix) { compute(std::move(matrix)); }
      29             : 
      30             :     /// @}
      31             : 
      32             :   public:
      33             :     /// @name Factorization
      34             :     /// @{
      35             : 
      36             :     /// Perform the QR factorization of the given matrix.
      37             :     void compute(Matrix &&matrix);
      38             :     /// Perform the QR factorization of the given matrix.
      39             :     void compute(const Matrix &matrix);
      40             : 
      41             :     /// @}
      42             : 
      43             :   public:
      44             :     /// @name   Retrieving the Q factor
      45             :     /// @{
      46             : 
      47             :     /// Compute the product QᵀB, overwriting B with the result.
      48             :     void apply_QT_inplace(Matrix &B) const;
      49             :     /// Compute the product QᵀB.
      50             :     Matrix apply_QT(const Matrix &B) const;
      51             :     /// Compute the product QᵀB.
      52             :     Matrix &&apply_QT(Matrix &&B) const;
      53             : 
      54             :     /// Compute the product QB, overwriting B with the result.
      55             :     void apply_Q_inplace(Matrix &X) const;
      56             :     /// Compute the product QB.
      57             :     Matrix apply_Q(const Matrix &X) const;
      58             :     /// Compute the product QB.
      59             :     Matrix &&apply_Q(Matrix &&B) const;
      60             : 
      61             :     /// Compute the unitary matrix Q and copy it to the given matrix.
      62             :     void get_Q_inplace(SquareMatrix &Q) const;
      63             :     /// Compute the unitary matrix Q.
      64             :     SquareMatrix get_Q() const;
      65             : 
      66             :     /// @}
      67             : 
      68             :   public:
      69             :     /// @name   Retrieving the R factor
      70             :     /// @{
      71             : 
      72             :     /// Get the upper-triangular matrix R, reusing the internal storage.
      73             :     /// @warning    After calling this function, the HouseholderQR object is no
      74             :     ///             longer valid, because this function steals its storage.
      75             :     Matrix &&steal_R();
      76             : 
      77             :     /// Copy the upper-triangular matrix R to the given matrix.
      78             :     void get_R_inplace(Matrix &R) const;
      79             :     /// Get a copy of the upper-triangular matrix R.
      80             :     Matrix get_R() const &;
      81             :     /// Get the upper-triangular matrix R.
      82             :     Matrix &&get_R() && { return steal_R(); }
      83             : 
      84             :     /// @}
      85             : 
      86             :   public:
      87             :     /// @name   Solving systems of equations and least-squares problems
      88             :     /// @{
      89             : 
      90             :     /// Solve the system AX = B or QRX = B.
      91             :     /// Matrix B is overwritten with the result X. If the matrix A is square,
      92             :     /// no new allocations occur, and the storage of B is reused for X.
      93             :     /// If A is not square, new storage will be allocated for X.
      94             :     void solve_inplace(Matrix &B) const;
      95             :     /// Solve the system AX = B or QRX = B.
      96             :     Matrix solve(const Matrix &B) const;
      97             :     /// Solve the system AX = B or QRX = B.
      98             :     Matrix &&solve(Matrix &&B) const;
      99             :     /// Solve the system Ax = b or QRx = b.
     100             :     Vector solve(const Vector &B) const;
     101             :     /// Solve the system Ax = b or QRx = b.
     102             :     Vector &&solve(Vector &&B) const;
     103             : 
     104             :     /// @}
     105             : 
     106             :   public:
     107             :     /// @name   Access to internal representation
     108             :     /// @{
     109             : 
     110             :     /// Check if this object contains a valid factorization.
     111          19 :     bool is_factored() const { return state == Factored; }
     112             : 
     113             :     /// Get the internal storage of the strict upper-triangular part of R and
     114             :     /// the Householder reflector vectors W.
     115           1 :     const Matrix &get_RW() const & { return RW; }
     116             :     /// @copydoc    get_RW
     117             :     Matrix &&get_RW() && { return std::move(RW); }
     118             :     /// Get the internal storage of the diagonal elements of R.
     119           1 :     const Vector &get_R_diag() const & { return R_diag; }
     120             :     /// @copydoc    get_R_diag
     121             :     Vector &&get_R_diag() && { return std::move(R_diag); }
     122             : 
     123             :     /// @}
     124             : 
     125             :   private:
     126             :     /// The actual QR factorization algorithm.
     127             :     void compute_factorization();
     128             :     /// Back substitution algorithm for solving upper-triangular systems RX = B.
     129             :     void back_subs(const Matrix &B, Matrix &X) const;
     130             : 
     131             :   private:
     132             :     /// Result of a Householder QR factorization: stores the strict
     133             :     /// upper-triangular part of matrix R and the full matrix of scaled
     134             :     /// Householder reflection vectors W. The reflection vectors have norm √2.
     135             :     Matrix RW;
     136             :     /// Contains the diagonal elements of R.
     137             :     Vector R_diag;
     138             : 
     139             :     enum State {
     140             :         NotFactored = 0,
     141             :         Factored    = 1,
     142             :     } state = NotFactored;
     143             : };
     144             : 
     145             : /// Print the Q and R matrices of a HouseholderQR object.
     146             : /// @related    HouseholderQR
     147             : std::ostream &operator<<(std::ostream &os, const HouseholderQR &qr);