Decompositions

Matrix factorizations: QR, LU, Hessenberg, Householder. More...

Functions
template<typename T , Size R, Size C>
constexpr Matrix< T, R, R >	consteig::house (Matrix< T, R, C > a)
	Compute a Householder reflector for a square matrix.
template<typename T , Size R, Size C>
constexpr QRMatrix< T, R >	consteig::qr_hessenberg (const Matrix< T, R, C > a)
	QR decomposition optimized for upper Hessenberg matrices.
template<typename T , Size R, Size C>
constexpr QRMatrix< T, R >	consteig::qr (const Matrix< T, R, C > a)
	General QR decomposition using Givens rotations.
template<typename T , Size R, Size C, Size L = R>
constexpr PHMatrix< T, R >	consteig::hess (Matrix< T, R, C > a)
	Reduce a square matrix to upper Hessenberg form.
template<typename T , Size S>
constexpr LUMatrix< T, S >	consteig::lu (const Matrix< T, S, S > &a)
	LU decomposition with partial pivoting.
template<typename T , Size S>
constexpr Matrix< T, S, 1 >	consteig::lu_solve (const LUMatrix< T, S > &lu, const Matrix< T, S, 1 > &b)
	Solve the linear system Ax = b given the LU factorization of A.

Detailed Description

Matrix factorizations: QR, LU, Hessenberg, Householder.

Function Documentation

house()

template<typename T , Size R, Size C>

constexpr Matrix< T, R, R > consteig::house ( Matrix< T, R, C > a )

constexpr

Compute a Householder reflector for a square matrix.

Returns an orthogonal matrix H = I - 2*v*v^T that zeros the first column of a below the first subdiagonal. Row 0 of the first column is untouched; row 1 receives the reflected norm; rows 2 through R-1 become zero.

Used internally by hess for Hessenberg reduction.

Template Parameters

T	Floating-point scalar type.
R	Number of rows (must equal C).
C	Number of columns.

Parameters

a	Square input matrix (only the first column is used).

Returns

R×R Householder reflector matrix.

Precondition

R == C and T must be floating-point (both enforced by static_assert).

qr_hessenberg()

template<typename T , Size R, Size C>

constexpr QRMatrix< T, R > consteig::qr_hessenberg ( const Matrix< T, R, C > a )

constexpr

QR decomposition optimized for upper Hessenberg matrices.

Factors an upper Hessenberg matrix as A = Q * R using Givens rotations. Because A has at most one subdiagonal nonzero per column, each rotation only needs to update two adjacent rows of R and two adjacent columns of Q, giving better performance than the general qr() for Hessenberg inputs.

Used internally by the eigenvalue solver's QR iteration loop. Prefer the general qr for non-Hessenberg matrices.

Template Parameters

T	Floating-point scalar type.
R	Number of rows.
C	Number of columns.

Parameters

a	Upper Hessenberg square matrix.

Returns

QRMatrix containing _q (orthogonal) and _r (upper-triangular).

Precondition

R == C (enforced by static_assert).

qr()

template<typename T , Size R, Size C>

constexpr QRMatrix< T, R > consteig::qr ( const Matrix< T, R, C > a )

constexpr

General QR decomposition using Givens rotations.

Factors a square matrix as A = Q * R where Q is orthogonal and R is upper triangular. Uses Givens rotations rather than Gram-Schmidt, avoiding the loss of orthogonality that Gram-Schmidt can suffer from floating-point cancellation.

For matrices already in Hessenberg form, qr_hessenberg is faster.

Template Parameters

T	Floating-point scalar type.
R	Number of rows.
C	Number of columns.

Parameters

a	Square input matrix.

Returns

QRMatrix containing _q (orthogonal) and _r (upper-triangular).

Precondition

R == C (enforced by static_assert).

static constexpr consteig::Matrix<double, 3, 3> A{{
    {12.0, -51.0, 4.0},
    { 6.0, 167.0, -68.0},
    {-4.0,  24.0, -41.0}
}};
static constexpr auto qr_result = consteig::qr(A);
// qr_result._q is orthogonal, qr_result._r is upper triangular
// A ≈ qr_result._q * qr_result._r

hess()

template<typename T , Size R, Size C, Size L = R>

constexpr PHMatrix< T, R > consteig::hess ( Matrix< T, R, C > a )

constexpr

Reduce a square matrix to upper Hessenberg form.

Computes H = P^T * A * P where H is upper Hessenberg (zero below the first subdiagonal) and P is an orthogonal matrix accumulated from Householder reflectors. This is a similarity transformation: H and A have the same eigenvalues.

Reducing to Hessenberg form before QR iteration cuts the cost of each QR step from O(n^3) to O(n^2), making the overall eigenvalue solver O(n^3).

Used internally by eig, eig_shifted_qr, and eig_double_shifted_qr.

Template Parameters

T	Floating-point scalar type.
R	Number of rows (must equal C).
C	Number of columns.
L	Internal recursion parameter; do not specify (defaults to R).

Parameters

a	Square input matrix.

Returns

PHMatrix containing _p (orthogonal) and _h (Hessenberg).

Precondition

R == C and T must be floating-point (both enforced by static_assert).

lu()

template<typename T , Size S>

constexpr LUMatrix< T, S > consteig::lu ( const Matrix< T, S, S > & a )

constexpr

LU decomposition with partial pivoting.

Factors a square matrix as PA = LU where P is a row permutation, L is unit lower triangular, and U is upper triangular. Partial pivoting selects the largest magnitude element as the pivot to control rounding error growth.

Used internally by eigenvectors for inverse iteration. Can also be used directly for solving linear systems via lu_solve.

Template Parameters

T	Scalar type (works with both real and Complex types).
S	Matrix dimension.

Parameters

a	Square input matrix.

Returns

LUMatrix containing _l, _u, and permutation _p.

lu_solve()

template<typename T , Size S>

constexpr Matrix< T, S, 1 > consteig::lu_solve ( const LUMatrix< T, S > & lu, const Matrix< T, S, 1 > & b )

constexpr

Solve the linear system Ax = b given the LU factorization of A.

Uses the result of lu to solve in two triangular passes:

1. Forward substitution to solve Ly = Pb.
2. Backward substitution to solve Ux = y.

Nearly singular systems (diagonal element of U below 1e-30) are handled gracefully to support inverse iteration in eigenvectors.

Template Parameters

T	Scalar type.
S	System dimension.

Parameters

lu	LU factorization from lu.
b	Right-hand side vector (S×1).

Returns

Solution vector x (S×1) satisfying Ax = b.