#include <include/linalg/HouseholderQR.hpp>

Detailed Description

QR factorization using Householder reflectors.

Factorizes an m×n matrix with m >= n into an m×m unitary factor Q and an m×n upper triangular factor R.

It can be used for solving square systems of equations or for finding a least squares solution to an overdetermined system of equations.

This version does not use column pivoting, and is not rank-revealing.

Examples: Basic-matrix-operations-no-iostream.ino, Basic-matrix-operations.ino, and QR-PerfTest.cpp.

Definition at line 18 of file HouseholderQR.hpp.

Collaboration diagram for HouseholderQR:

Constructors
	HouseholderQR ()=default
	Default constructor. More...

	HouseholderQR (const Matrix &matrix)
	Factorize the given matrix. More...

	HouseholderQR (Matrix &&matrix)
	Factorize the given matrix. More...

Factorization
void	compute (Matrix &&matrix)
	Perform the QR factorization of the given matrix. More...

void	compute (const Matrix &matrix)
	Perform the QR factorization of the given matrix. More...

Retrieving the Q factor
void	apply_QT_inplace (Matrix &B) const
	Compute the product QᵀB, overwriting B with the result. More...

Matrix	apply_QT (const Matrix &B) const
	Compute the product QᵀB. More...

Matrix &&	apply_QT (Matrix &&B) const
	Compute the product QᵀB. More...

void	apply_Q_inplace (Matrix &X) const
	Compute the product QB, overwriting B with the result. More...

Matrix	apply_Q (const Matrix &X) const
	Compute the product QB. More...

Matrix &&	apply_Q (Matrix &&B) const
	Compute the product QB. More...

void	get_Q_inplace (SquareMatrix &Q) const
	Compute the unitary matrix Q and copy it to the given matrix. More...

SquareMatrix	get_Q () const
	Compute the unitary matrix Q. More...

Retrieving the R factor
Matrix &&	steal_R ()
	Get the upper-triangular matrix R, reusing the internal storage. More...

void	get_R_inplace (Matrix &R) const
	Copy the upper-triangular matrix R to the given matrix. More...

Matrix	get_R () const &
	Get a copy of the upper-triangular matrix R. More...

Matrix &&	get_R () &&
	Get the upper-triangular matrix R. More...

Solving systems of equations and least-squares problems
void	solve_inplace (Matrix &B) const
	Solve the system AX = B or QRX = B. More...

Matrix	solve (const Matrix &B) const
	Solve the system AX = B or QRX = B. More...

Matrix &&	solve (Matrix &&B) const
	Solve the system AX = B or QRX = B. More...

Vector	solve (const Vector &B) const
	Solve the system Ax = b or QRx = b. More...

Vector &&	solve (Vector &&B) const
	Solve the system Ax = b or QRx = b. More...

Access to internal representation
bool	is_factored () const
	Check if this object contains a valid factorization. More...

const Matrix &	get_RW () const &
	Get the internal storage of the strict upper-triangular part of R and the Householder reflector vectors W. More...

Matrix &&	get_RW () &&
	Get the internal storage of the strict upper-triangular part of R and the Householder reflector vectors W. More...

const Vector &	get_R_diag () const &
	Get the internal storage of the diagonal elements of R. More...

Vector &&	get_R_diag () &&
	Get the internal storage of the diagonal elements of R. More...

Private Types
enum	State { NotFactored = 0 , Factored = 1 }

Private Member Functions
void	compute_factorization ()
	The actual QR factorization algorithm. More...

void	back_subs (const Matrix &B, Matrix &X) const
	Back substitution algorithm for solving upper-triangular systems RX = B. More...

Private Attributes
Matrix	RW
	Result of a Householder QR factorization: stores the strict upper-triangular part of matrix R and the full matrix of scaled Householder reflection vectors W. More...

Vector	R_diag
	Contains the diagonal elements of R. More...

enum HouseholderQR::State	state = NotFactored

Related Functions
(Note that these are not member functions.)
std::ostream &	operator<< (std::ostream &os, const HouseholderQR &qr)
	Print the Q and R matrices of a HouseholderQR object. More...

Member Enumeration Documentation

◆ State

enum State

private

Enumerator
NotFactored
Factored

Definition at line 139 of file HouseholderQR.hpp.

Constructor & Destructor Documentation

◆ HouseholderQR() [1/3]

HouseholderQR ( )

default

Default constructor.

◆ HouseholderQR() [2/3]

HouseholderQR ( const Matrix & matrix )

inline

Factorize the given matrix.

Definition at line 26 of file HouseholderQR.hpp.

◆ HouseholderQR() [3/3]

HouseholderQR ( Matrix && matrix )

inline

Factorize the given matrix.

Definition at line 28 of file HouseholderQR.hpp.

Member Function Documentation

◆ compute() [1/2]

void compute ( Matrix && matrix )

Perform the QR factorization of the given matrix.

Examples: QR-PerfTest.cpp.

Definition at line 9 of file HouseholderQR.ipp.

◆ compute() [2/2]

void compute ( const Matrix & matrix )

Perform the QR factorization of the given matrix.

Definition at line 15 of file HouseholderQR.ipp.

◆ apply_QT_inplace()

void apply_QT_inplace ( Matrix & B ) const

Compute the product QᵀB, overwriting B with the result.

Implementation

void HouseholderQR::apply_QT_inplace(Matrix &B) const {
    assert(is_factored());
    assert(RW.rows() == B.rows());
    // Apply the Householder reflectors to each column of B.
    for (size_t i = 0; i < B.cols(); ++i) {
        // Recall that the Householder reflector H is applied as follows:
        //     bᵢ'[k+1:m] = H·bᵢ[k+1:m]
        //                = bᵢ[k+1:m] - wₖ·wₖᵀ·bᵢ[k+1:m]
        for (size_t k = 0; k < RW.cols(); ++k) {
            // Compute wₖᵀ·bᵢ
            double dot_product = 0;
            for (size_t r = k; r < RW.rows(); ++r)
                dot_product += RW(r, k) * B(r, i);
            // Subtract wₖ·wₖᵀ·bᵢ
            for (size_t r = k; r < RW.rows(); ++r)
                B(r, i) -= RW(r, k) * dot_product;
        }
    }
}

Definition at line 189 of file HouseholderQR.cpp.

◆ apply_QT() [1/2]

Matrix apply_QT ( const Matrix & B ) const

Compute the product QᵀB.

Definition at line 21 of file HouseholderQR.ipp.

◆ apply_QT() [2/2]

Matrix && apply_QT ( Matrix && B ) const

Compute the product QᵀB.

Definition at line 27 of file HouseholderQR.ipp.

◆ apply_Q_inplace()

void apply_Q_inplace ( Matrix & X ) const

Compute the product QB, overwriting B with the result.

Implementation

void HouseholderQR::apply_Q_inplace(Matrix &X) const {
    assert(is_factored());
    assert(RW.rows() == X.rows());
    // Apply the Householder reflectors in reverse order to each column of X.
    for (size_t i = 0; i < X.cols(); ++i) {
        // Recall that the Householder reflector H is applied as follows:
        //     xᵢ'[k+1:m] = H·xᵢ[k+1:m]
        //                = xᵢ[k+1:m] - wₖ·wₖᵀ·xᵢ[k+1:m]
        for (size_t k = RW.cols(); k-- > 0;) {
            // Compute wₖᵀ·xᵢ
            double dot_product = 0;
            for (size_t r = k; r < RW.rows(); ++r)
                dot_product += RW(r, k) * X(r, i);
            // Subtract wₖ·wₖᵀ·xᵢ
            for (size_t r = k; r < RW.rows(); ++r)
                X(r, i) -= RW(r, k) * dot_product;
        }
    }
}

Definition at line 215 of file HouseholderQR.cpp.

◆ apply_Q() [1/2]

Matrix apply_Q ( const Matrix & X ) const

Compute the product QB.

Examples: QR-PerfTest.cpp.

Definition at line 32 of file HouseholderQR.ipp.

◆ apply_Q() [2/2]

Matrix && apply_Q ( Matrix && B ) const

Compute the product QB.

Definition at line 38 of file HouseholderQR.ipp.

◆ get_Q_inplace()

void get_Q_inplace ( SquareMatrix & Q ) const

Compute the unitary matrix Q and copy it to the given matrix.

Definition at line 80 of file HouseholderQR.ipp.

◆ get_Q()

SquareMatrix get_Q ( ) const

Compute the unitary matrix Q.

Examples: Basic-matrix-operations-no-iostream.ino.

Definition at line 87 of file HouseholderQR.ipp.

◆ steal_R()

Matrix && steal_R ( )

Get the upper-triangular matrix R, reusing the internal storage.

Warning: After calling this function, the HouseholderQR object is no longer valid, because this function steals its storage.

Definition at line 66 of file HouseholderQR.ipp.

◆ get_R_inplace()

void get_R_inplace ( Matrix & R ) const

Copy the upper-triangular matrix R to the given matrix.

Definition at line 43 of file HouseholderQR.ipp.

◆ get_R() [1/2]

Matrix get_R ( ) const &

Get a copy of the upper-triangular matrix R.

Examples: Basic-matrix-operations-no-iostream.ino, and QR-PerfTest.cpp.

Definition at line 60 of file HouseholderQR.ipp.

◆ get_R() [2/2]

Matrix&& get_R ( ) &&

inline

Get the upper-triangular matrix R.

Definition at line 82 of file HouseholderQR.hpp.

◆ solve_inplace()

void solve_inplace ( Matrix & B ) const

Solve the system AX = B or QRX = B.

Matrix B is overwritten with the result X. If the matrix A is square, no new allocations occur, and the storage of B is reused for X. If A is not square, new storage will be allocated for X.

Implementation

void HouseholderQR::solve_inplace(Matrix &B) const {
    // If AX = B, then QRX = B, or RX = QᵀB, so first apply Qᵀ to B:
    apply_QT_inplace(B);
 
    // To solve RX = QᵀB, use back substitution:
 
    // If the matrix is square, B and X have the same size, so we can reuse B's
    // storage if we operate on B directly:
    if (RW.cols() == RW.rows()) {
        back_subs(B, B);
    }
    // If the matrix is rectangular, the sizes of B and X differ, so use a
    // separate result variable:
    else {
        Matrix X(RW.cols(), B.cols());
        back_subs(B, X);
        B = std::move(X);
    }
}

Definition at line 274 of file HouseholderQR.cpp.

◆ solve() [1/4]

Matrix solve ( const Matrix & B ) const

Solve the system AX = B or QRX = B.

Examples: Basic-matrix-operations-no-iostream.ino, and Basic-matrix-operations.ino.

Definition at line 93 of file HouseholderQR.ipp.

◆ solve() [2/4]

Matrix && solve ( Matrix && B ) const

Solve the system AX = B or QRX = B.

Definition at line 100 of file HouseholderQR.ipp.

◆ solve() [3/4]

Vector solve ( const Vector & B ) const

Solve the system Ax = b or QRx = b.

Definition at line 105 of file HouseholderQR.ipp.

◆ solve() [4/4]

Vector && solve ( Vector && B ) const

Solve the system Ax = b or QRx = b.

Definition at line 109 of file HouseholderQR.ipp.

◆ is_factored()

bool is_factored ( ) const

inline

Check if this object contains a valid factorization.

Definition at line 111 of file HouseholderQR.hpp.

◆ get_RW() [1/2]

const Matrix& get_RW ( ) const &

inline

Get the internal storage of the strict upper-triangular part of R and the Householder reflector vectors W.

Definition at line 115 of file HouseholderQR.hpp.

◆ get_RW() [2/2]

Matrix&& get_RW ( ) &&

inline

Get the internal storage of the strict upper-triangular part of R and the Householder reflector vectors W.

Definition at line 117 of file HouseholderQR.hpp.

◆ get_R_diag() [1/2]

const Vector& get_R_diag ( ) const &

inline

Get the internal storage of the diagonal elements of R.

Definition at line 119 of file HouseholderQR.hpp.

◆ get_R_diag() [2/2]

Vector&& get_R_diag ( ) &&

inline

Get the internal storage of the diagonal elements of R.

Definition at line 121 of file HouseholderQR.hpp.

◆ compute_factorization()

void compute_factorization ( )

private

The actual QR factorization algorithm.

Precondition: RW contains the matrix A to be factorized; RW.rows() >= RW.cols(); R_diag.size() == RW.cols()

Postcondition: The strict upper-triangular part of RW contains the strict upper-triangular part of factor R. R_diag contains the diagonal of R.; The lower-triangular part (including the diagonal) of RW contains the Householder reflectors wₖ, with ‖wₖ‖ = √2.; apply_Q(get_R()) == A == get_Q() * get_R() (up to rounding errors)

Implementation

void HouseholderQR::compute_factorization() {
    // For the intermediate calculations, we'll be working with RW.
    // It is initialized to the rectangular matrix to be factored.
    // At the end of this function, RW will contain the strict
    // upper-triangular part of the matrix R (without the diagonal),
    // and the complete scaled matrix of reflection vectors W, which is a
    // lower-triangular matrix. The diagonal of R is stored separately in
    // R_diag.
 
    assert(RW.rows() >= RW.cols());
    assert(R_diag.size() == RW.cols());
 
    // Helper function to square a number
    auto sq = [](double x) { return x * x; };
 
    for (size_t k = 0; k < RW.cols(); ++k) {
        // Introduce a column vector x = A[k:M,k], it's the lower part of the
        // k-th column of the matrix.
        // First compute the norm of x:
 
        double sq_norm_x = 0;
        for (size_t i = k; i < RW.rows(); ++i)
            sq_norm_x += sq(RW(i, k));
        double norm_x = std::sqrt(sq_norm_x);
 
        // x consists of two parts: its first element, x₀, and the rest, xₛ
        //     x = (x₀, xₛ)
        // You can express the norm of x in terms of the norms of the two parts:
        //     ‖x‖² = x₀² + ‖xₛ‖²
        double &x_0 = RW(k, k);
 
        // The goal of QR factorization is to introduce zeros below the diagonal
        // in the R factor by transforming the vector x to a new vector that is
        // all zero, except for the first component.
        // Call this vector xₕ, the Householder reflection of x.
        // Since the transformation has to be unitary (Q is a unitary matrix),
        // this operation has to preserve the 2-norm. This means that the
        // nonzero component of xₕ has to have the same 2-norm (energy) as x:
        //     xₕ = (±‖x‖, 0, ..., 0) = ±‖x‖·̅e₁
        // where ̅e₁ is the first standard basis vector (1, 0, ..., 0).
        //
        // There are two (real) vectors that have the same energy as x in their
        // first component, because the sign doesn't affect the energy.
        // For numerical reasons, it's best to pick the sign opposite to the
        // sign of the first component of x, so that x and xₕ are far apart.
        //
        // The reflector vector vₖ is the difference between the original
        // vector x and its Householder reflection xₕ:
        //     vₖ = x - xₕ
        //        = x - (-sign(x₀)·‖x‖·̅e₁)
        //        = x + sign(x₀)·‖x‖·̅e₁
        //        = (x₀ + sign(x₀)·‖x‖, xₛ)
        //
        // Since vₖ will later be used to construct a projection matrix, it
        // should be normalized.
        // Computing the full norm is not necessary, since
        //     ‖vₖ‖² = (x₀ + sign(x₀)·‖x‖)² + ‖xₛ‖²
        //           = x₀² + 2·x₀·sign(x₀)·‖x‖ + ‖x‖² + ‖xₛ‖²
        //           = 2·x₀·sign(x₀)·‖x‖ + 2·‖x‖²
        //           = 2·|x₀|·‖x‖ + 2·‖x‖²
        //     ‖vₖ‖  = √2·√(|x₀|·‖x‖ + ‖x‖²)
        //
        // Normalize vₖ and call this uₖ:
        //     uₖ = vₖ / ‖vₖ‖
        //
        // For reasons that will become apparent in a moment, we'll keep the
        // factor of √2 in our normalization of vₖ, call this vector wₖ:
        //     wₖ = √2·uₖ
        //        = √2·vₖ / ‖vₖ‖
        //        = vₖ / √(|x₀|·‖x‖ + ‖x‖²)
        // Note how the sum only adds up numbers with the same sign. This
        // prevents catastrophic cancelations.
        //
        // It is clear that if ‖x‖ = 0, this normalization will fail. In that
        // case, we set uₖ = ̅e₁ or wₖ = √2·̅e₁.
        //
        // x will be overwritten by wₖ. The vector xₕ only has a single nonzero
        // component. It is saved in the R_diag vector.
 
        if (norm_x >= std::numeric_limits<double>::min() * 2) {
            double x_p = -std::copysign(norm_x, x_0); // -sign(x₀)·‖x‖
            double v_0 = x_0 - x_p;
            double norm_v_sq2 = std::sqrt(std::abs(x_0) * norm_x + sq_norm_x);
 
            // Overwrite x with vₖ:
            x_0 = v_0;
            // the other components of x (xₛ) are already equal to the bottom
            // part of vₖ, so they don't have to be overwritten explicitly.
 
            // Then normalize x (= vₖ) to obtain wₖ:
            for (size_t i = k; i < RW.rows(); ++i)
                RW(i, k) /= norm_v_sq2;
 
            // Save the first component of xₕ:
            R_diag(k) = x_p;
        } else {
            // Overwrite x with wₖ = √2·̅e₁:
            x_0 = std::sqrt(2);
            // the other components of x (xₛ) are already equal to zero, since
            // ‖x‖ = 0.
 
            // Save the first component of xₕ:
            R_diag(k) = 0;
        }
 
        // Now that the reflection vector vₖ (wₖ) is known, the rest of the
        // matrix A can be updated, that is A[k:m,k:n].
        //
        // The reflection was defined in terms of the reflection vector
        // vₖ = x - xₕ. To reflect x onto xₕ, you can first project x onto the
        // orthogonal complement of the space spanned by vₖ, call this
        // projection xₚ, and then add the difference between xₚ and x to xₚ
        // again, in order to reflect about the orthogonal complement of vₖ,
        // rather than projecting onto it:
        //     xₕ = xₚ + (xₚ - x)
        //        = 2·xₚ - x
        // The projection matrix onto the orthogonal complement of the space
        // spanned by vₖ is:
        //     P⟂ = I - uₖ·uₖᵀ
        // where I is the identity matrix, and uₖ the unit vector in the
        // direction of vₖ, as defined above.
        // This allows us to compute xₚ:
        //     xₚ = P⟂·x
        //        = I·x - uₖ·uₖᵀ·x
        //        = x - uₖ·uₖᵀ·x
        //     xₕ = 2·xₚ - x
        //        = 2·x - 2·uₖ·uₖᵀ·x - x
        //        = x - 2·uₖ·uₖᵀ·x
        // Because of our choice of wₖ = √2·uₖ earlier, this simplifies to
        //     xₕ = x - wₖ·wₖᵀ·x
        // The Householder reflector is thus the linear transformation given by
        //     H = I - wₖ·wₖᵀ
        //
        // The final step is to apply this reflector to the remaining part of
        // the matrix A, i.e. A[k:m,k:n].
        // The first column x = A[k:m,k] has already been updated, that's how
        // the reflector was computed, so only the submatrix A[k+1:m,k:n] has
        // to be updated.
        //     A'[k+1:m,k:n] = H·A[k+1:m,k:n]
        //                   = A[k+1:m,k:n] - wₖ·wₖᵀ·A[k+1:m,k:n]
        // This can be computed column-wise:
        //     aᵢ' = aᵢ - wₖ·wₖᵀ·aᵢ
        // where aᵢ is the i-th column of A.
 
        for (size_t c = k + 1; c < RW.cols(); ++c) {
            // Compute wₖᵀ·aᵢ
            double dot_product = 0;
            for (size_t r = k; r < RW.rows(); ++r)
                dot_product += RW(r, k) * RW(r, c);
            // Subtract wₖ·wₖᵀ·aᵢ
            for (size_t r = k; r < RW.rows(); ++r)
                RW(r, c) -= RW(r, k) * dot_product;
        }
    }
    state = Factored;
}

Definition at line 26 of file HouseholderQR.cpp.

◆ back_subs()

void back_subs	(	const Matrix &	B,
		Matrix &	X
	)		const

private

Back substitution algorithm for solving upper-triangular systems RX = B.

Implementation

void HouseholderQR::back_subs(const Matrix &B, Matrix &X) const {
    // Solve upper triangular system RX = B by solving each column of B as a
    // vector system Rxᵢ = bᵢ
    //
    //     ┌                 ┐┌     ┐   ┌     ┐
    //     │ r₁₁ r₁₂ r₁₃ r₁₄ ││ x₁ᵢ │   │ b₁ᵢ │
    //     │     r₂₂ r₂₃ r₂₄ ││ x₂ᵢ │ = │ b₂ᵢ │
    //     │         r₃₃ r₃₄ ││ x₃ᵢ │   │ b₃ᵢ │
    //     │             r₄₄ ││ x₄ᵢ │   │ b₄ᵢ │
    //     └                 ┘└     ┘   └     ┘
    //
    // b₄ᵢ = r₄₄·x₄ᵢ                     ⟺ x₄ᵢ = b₄ᵢ/r₄₄
    // b₃ᵢ = r₃₃·x₃ᵢ + r₃₄·x₄ᵢ           ⟺ x₃ᵢ = (b₃ᵢ - r₃₄·x₄ᵢ)/r₃₃
    // b₂ᵢ = r₂₂·x₂ᵢ + r₂₃·x₃ᵢ + r₂₄·x₄ᵢ ⟺ x₂ᵢ = (b₂ᵢ - r₂₃·x₃ᵢ + r₂₄·x₄ᵢ)/r₂₂
    // ...
 
    for (size_t i = 0; i < B.cols(); ++i) {
        for (size_t k = RW.cols(); k-- > 0;) {
            X(k, i) = B(k, i);
            for (size_t j = k + 1; j < RW.cols(); ++j) {
                X(k, i) -= RW(k, j) * X(j, i);
            }
            X(k, i) /= R_diag(k);
        }
    }
}

Definition at line 241 of file HouseholderQR.cpp.

Friends And Related Function Documentation

◆ operator<<()

std::ostream & operator<<	(	std::ostream &	os,
		const HouseholderQR &	qr
	)

Definition at line 121 of file HouseholderQR.ipp.

Member Data Documentation

◆ RW

Matrix RW

private

Result of a Householder QR factorization: stores the strict upper-triangular part of matrix R and the full matrix of scaled Householder reflection vectors W.

The reflection vectors have norm √2.

Definition at line 135 of file HouseholderQR.hpp.

◆ R_diag

Vector R_diag

private

Contains the diagonal elements of R.

Definition at line 137 of file HouseholderQR.hpp.

◆ state

enum HouseholderQR::State state = NotFactored

private

The documentation for this class was generated from the following files:

Detailed Description

Constructors

Factorization

Retrieving the Q factor

Retrieving the R factor

Solving systems of equations and least-squares problems

Access to internal representation

Private Types

Private Member Functions

Private Attributes

Related Functions

Member Enumeration Documentation

◆ State

Constructor & Destructor Documentation

◆ HouseholderQR() [1/3]

◆ HouseholderQR() [2/3]

◆ HouseholderQR() [3/3]

Member Function Documentation

◆ compute() [1/2]

◆ compute() [2/2]

◆ apply_QT_inplace()

Implementation

◆ apply_QT() [1/2]

◆ apply_QT() [2/2]

◆ apply_Q_inplace()

Implementation

◆ apply_Q() [1/2]

◆ apply_Q() [2/2]

◆ get_Q_inplace()

◆ get_Q()

◆ steal_R()

◆ get_R_inplace()

◆ get_R() [1/2]

◆ get_R() [2/2]

◆ solve_inplace()

Implementation

◆ solve() [1/4]

◆ solve() [2/4]

◆ solve() [3/4]

◆ solve() [4/4]

◆ is_factored()

◆ get_RW() [1/2]

◆ get_RW() [2/2]

◆ get_R_diag() [1/2]

◆ get_R_diag() [2/2]

◆ compute_factorization()

Implementation

◆ back_subs()

Implementation

Friends And Related Function Documentation

◆ operator<<()

Member Data Documentation

◆ RW

◆ R_diag

◆ state