Matrix Decompositions — LU, Cholesky, QR, SVD¶

Algorithm Flow¶

LU Decomposition (Doolittle)¶

Goal: Factor square matrix \(A\) into \(A = P \cdot L \cdot U\), where \(L\) is unit lower triangular, \(U\) is upper triangular, and \(P\) is a permutation matrix.

Algorithm (with partial pivoting):

1. Copy \(A\) to workspace \(U\), initialize \(L = I\), \(P = I\)

2. For each column \(k = 0, \dots, n-2\):

3. Pivot search: find max |element| in \(U[k..n-1, k]\), record row \(p\)

4. Swap: swap rows in \(U\), \(P\), and partial rows in \(L\)

5. Eliminate: for each row \(i = k+1, \dots, n-1\):

   • \(L_{ik} = U_{ik} / U_{kk}\)
   • \(U_{i, k..n-1} \;\;{-}{=}\;\; L_{ik} \cdot U_{k, k..n-1}\)
   • Return \(L, U, P\)

Complexity: \(O(n^3)\) FLOPs.

Guidance: Always enable pivoting for general matrices; only skip if diagonal dominance is guaranteed.

Cholesky Decomposition¶

Goal: Factor symmetric positive definite \(A\) into \(A = L \cdot L^T\), where \(L\) is lower triangular with positive diagonal.

Algorithm (in-place):

For each column \(j = 0, \dots, n-1\):

1. Diagonal: \(L_{jj} = \sqrt{A_{jj} - \sum_{k=0}^{j-1} L_{jk}^2}\)

2. Sub-diagonal: for each \(i = j+1, \dots, n-1\):

3. \(L_{ij} = \left( A_{ij} - \sum_{k=0}^{j-1} L_{ik} L_{jk} \right) / L_{jj}\)

Key Properties:

• Speed: ~half the FLOPs of LU (exploits symmetry)
• Stability: numerically stable for SPD, no pivoting needed
• Precondition: verify \(A\) is SPD via is_symmetric() + is_positive_definite()

Complexity: \(O(n^3 / 3)\) FLOPs.

QR Decomposition (Householder)¶

Goal: Factor \(m \times n\) matrix \(A\) into \(A = Q \cdot R\), with \(Q\) orthogonal, \(R\) upper triangular.

Algorithm (for each column \(k = 0, \dots, n-1\), \(m \geq n\)):

1. Extract \(\mathbf{x} = A[k..m-1, k]\)

2. Build Householder vector \(\mathbf{v}\):

3. \(\alpha = -\text{sign}(x_0) \cdot \|\mathbf{x}\|_2\)

4. \(\mathbf{v} = \mathbf{x} - \alpha \cdot \mathbf{e}_0\)
5. \(\beta = 2 / \|\mathbf{v}\|_2^2\)
6. Apply reflection: \(A[k..m-1, k..n-1] \;\;{=}\;\; A - \beta \cdot \mathbf{v} \cdot (\mathbf{v}^T \cdot A)\)
7. Accumulate \(Q\): \(Q = Q \cdot (I - \beta \cdot \mathbf{v} \cdot \mathbf{v}^T)\)

MGS vs Householder QR:

SVD Decomposition (Bilateral Jacobi + QR Refinement)¶

Goal: Factor \(m \times n\) matrix \(A\) into \(A = U \cdot \Sigma \cdot V^T\), with \(U, V\) orthogonal, \(\Sigma\) diagonal (\(\sigma_1 \geq \dots \geq \sigma_r \geq 0\)).

Algorithm:

Phase 1 — Gram matrix:

1. Compute \(G = A^T \cdot A\) (size \(n \times n\))
2. Problem reduces to eigen-decomposition of \(G\)

Phase 2 — Jacobi rotations:

1. Initialize \(V = I_n\), workspace \(B = G\)

2. Sweep all off-diagonal pairs \((p, q)\):

3. If \(|B_{pq}|\) > tolerance, compute Givens angle \(\theta\):

   • \(\tau = (B_{qq} - B_{pp}) / (2 \cdot B_{pq})\)
   • \(t = \text{sign}(\tau) / (|\tau| + \sqrt{1 + \tau^2})\)
   • \(c = 1 / \sqrt{1 + t^2}\), \(s = t \cdot c\)
   • Apply: \(B = J^T \cdot B \cdot J\), \(V = V \cdot J\)
   • Repeat sweeps until \(\max |B_{pq}| < \text{tolerance}\)

Phase 3 — Sort & truncate:

1. Sort eigenvalues (= \(\sigma^2\)) descending, permute \(V\) columns
2. Take square roots: \(\sigma_i = \sqrt{|\lambda_i|}\)
3. Recover \(U\): \(U = A \cdot V \cdot \Sigma^{-1}\)
4. Numerical rank: \(\#\{i \mid \sigma_i > \text{tolerance}\}\)

Complexity: ~\(O(n^3)\).

Caveats:

• Computing via \(A^T A\) squares the condition number; direct bilateral Jacobi on \(A\) is more stable
• Tolerance selection directly affects rank and pseudo-inverse result

MATRIX PROPERTIES & DECOMPOSITIONS¶

Matrix Property Checks¶

Check Symmetry¶

bool Mat::is_symmetric(float tolerance = 1e-6f) const;

Description:

Check whether a matrix is symmetric within the given tolerance. A matrix A is symmetric if A = A^T, i.e., A(i,j) = A(j,i) for all i, j.

Mathematical Principle:

For a symmetric matrix, all eigenvalues are real, and eigenvectors can be chosen to be orthogonal. Symmetric matrices are fundamental in many applications, especially in structural dynamics and optimization.

Parameters:

• float tolerance : Maximum allowed difference |A(i,j) - A(j,i)| (default: 1e-6).

Returns:

bool - true if approximately symmetric, false otherwise.

Usage Insights:

• Eigendecomposition: Symmetric matrices can use more efficient and stable eigendecomposition methods (e.g., Jacobi method).

• Cholesky Decomposition: Only symmetric positive definite matrices can be decomposed using Cholesky decomposition.

• Structural Dynamics: Stiffness and mass matrices in structural analysis are typically symmetric.

Check Positive Definiteness¶

bool Mat::is_positive_definite(float tolerance = 1e-6f, int max_minors_to_check = -1) const;

Description:

Check if a matrix is positive definite using Sylvester's criterion. A symmetric matrix A is positive definite if x^T A x > 0 for all non-zero vectors x, or equivalently, all eigenvalues are positive.

Mathematical Principle:

Sylvester's criterion states that a symmetric matrix is positive definite if and only if all leading principal minors are positive. The function checks the first few leading minors and diagonal elements for efficiency.

Parameters:

• float tolerance : Tolerance for numerical checks (default: 1e-6).
• int max_minors_to_check : Number of leading principal minors to check. -1 checks all; >0 checks only the first N.

Returns:

bool - true if matrix is positive definite, false otherwise.

Usage Insights:

• Cholesky Decomposition: Positive definite matrices can be decomposed using Cholesky decomposition, which is faster and more stable than LU decomposition.

• Optimization: Positive definite Hessian matrices indicate local minima in optimization problems.

• Stability Analysis: In control systems, positive definiteness of certain matrices ensures system stability.

Matrix Decomposition Structures¶

LU Decomposition Structure¶

struct Mat::LUDecomposition
{
    Mat L;                 // Lower triangular matrix (with unit diagonal)
    Mat U;                 // Upper triangular matrix
    Mat P;                 // Permutation matrix (if pivoting used)
    bool pivoted;          // Whether pivoting was used
    tiny_error_t status;   // Computation status

    LUDecomposition();
};

Description:

Container for LU decomposition results. The decomposition A = P * L * U (with pivoting) or A = L * U (without pivoting), where L is lower triangular with unit diagonal, U is upper triangular, and P is a permutation matrix.

Mathematical Principle:

LU decomposition factors a matrix into lower and upper triangular matrices, enabling efficient solution of linear systems. With pivoting, it handles near-singular matrices better.

Cholesky Decomposition Structure¶

struct Mat::CholeskyDecomposition
{
    Mat L;                 // Lower triangular matrix
    tiny_error_t status;   // Computation status

    CholeskyDecomposition();
};

Description:

Container for Cholesky decomposition results. For symmetric positive definite matrices, A = L * L^T, where L is lower triangular.

Mathematical Principle:

Cholesky decomposition is a specialized LU decomposition for symmetric positive definite matrices. It requires only half the storage and computation of LU decomposition.

QR Decomposition Structure¶

struct Mat::QRDecomposition
{
    Mat Q;                 // Orthogonal matrix (Q^T * Q = I)
    Mat R;                 // Upper triangular matrix
    tiny_error_t status;   // Computation status

    QRDecomposition();
};

Description:

Container for QR decomposition results. A = Q * R, where Q is orthogonal (Q^T * Q = I) and R is upper triangular.

Mathematical Principle:

QR decomposition expresses a matrix as the product of an orthogonal matrix and an upper triangular matrix. It's numerically stable and fundamental for least squares problems.

SVD Decomposition Structure¶

struct Mat::SVDDecomposition
{
    Mat U;                 // Left singular vectors (orthogonal matrix)
    Mat S;                 // Singular values (diagonal matrix or vector)
    Mat V;                 // Right singular vectors (orthogonal matrix, V^T)
    int rank;              // Numerical rank of the matrix
    int iterations;        // Number of iterations performed
    tiny_error_t status;   // Computation status

    SVDDecomposition();
};

Description:

Container for SVD decomposition results. A = U * S * V^T, where U and V are orthogonal matrices, and S contains singular values on the diagonal.

Mathematical Principle:

SVD is the most general matrix decomposition. The singular values reveal the matrix's rank, condition number, and enable computation of pseudo-inverse for rank-deficient matrices.

Matrix Decomposition Methods¶

LU Decomposition¶

Mat::LUDecomposition Mat::lu_decompose(bool use_pivoting = true) const;

Description:

Compute LU decomposition: A = P * L * U (with pivoting) or A = L * U (without pivoting). Efficient for solving multiple systems with the same coefficient matrix.

Mathematical Principle:

• Without pivoting: A = L * U, where L has unit diagonal

• With pivoting: P * A = L * U, where P is a permutation matrix

The decomposition enables solving Ax = b by solving Ly = Pb (forward substitution) then Ux = y (back substitution).

Parameters:

• bool use_pivoting : Whether to use partial pivoting for numerical stability (default: true).

Returns:

LUDecomposition containing L, U, P matrices and status.

Usage Insights:

• Multiple RHS: Once decomposed, solve multiple systems with different right-hand sides efficiently using solve_lu().

• Determinant: det(A) = det(P) * det(L) * det(U) = det(P) * det(U) (since det(L) = 1).

• Inverse: Can compute A^(-1) by solving LUx = eᵢ for each unit vector eᵢ.

• Performance: O(n³) for decomposition, O(n²) for each solve after decomposition.

Pitfalls:

• Pivoting is not optional in practice: Without pivoting, a perfectly valid matrix can still produce catastrophic round-off error.

• L is unit-diagonal: Do not expect diagonal values in L to carry scale information; the scale lives in U.

• Sign matters for determinants: Every row permutation changes the determinant sign through P.

Cholesky Decomposition¶

Mat::CholeskyDecomposition Mat::cholesky_decompose() const;

Description:

Compute Cholesky decomposition: A = L * L^T for symmetric positive definite matrices. Faster than LU for SPD matrices, commonly used in structural dynamics.

Mathematical Principle:

For a symmetric positive definite matrix A, there exists a unique lower triangular matrix L with positive diagonal elements such that A = L * L^T. This is essentially a specialized LU decomposition that takes advantage of symmetry.

Parameters:

None (matrix must be symmetric positive definite).

Returns:

CholeskyDecomposition containing L matrix and status.

Usage Insights:

• Efficiency: Requires approximately half the computation and storage of LU decomposition.

• Stability: More stable than LU for symmetric positive definite matrices.

• Applications:

• Structural dynamics: Mass and stiffness matrices are often SPD

• Optimization: Hessian matrices in Newton's method

• Statistics: Covariance matrices

• Error Handling: Returns error if matrix is not symmetric or not positive definite.

Pitfalls:

• SPD only: A matrix that is merely symmetric is not enough; positive definiteness is required.

• Symmetry tolerance: In floating-point code, tiny asymmetries can cause a failure unless the implementation tolerates them.

• Best when structure is known: If you already know the matrix is SPD, this is the right solver family; otherwise check first with is_symmetric() and is_positive_definite().

QR Decomposition¶

Mat::QRDecomposition Mat::qr_decompose() const;

Description:

Compute QR decomposition: A = Q * R, where Q is orthogonal and R is upper triangular. Numerically stable, used for least squares and orthogonalization.

Mathematical Principle:

QR decomposition expresses any matrix as the product of an orthogonal matrix Q (Q^T * Q = I) and an upper triangular matrix R. The decomposition is computed using the modified Gram-Schmidt process with re-orthogonalization.

Parameters:

None.

Returns:

QRDecomposition containing Q and R matrices and status.

Usage Insights:

• Least Squares: For overdetermined system Ax ≈ b, the solution minimizes ||Ax - b||₂ is x = R^(-1) * Q^T * b.

• Numerical Stability: QR decomposition is more stable than normal equations for least squares problems.

• Eigendecomposition: QR algorithm uses QR decomposition iteratively to find eigenvalues.

• Rank Revealing: The rank of A equals the number of non-zero diagonal elements of R.

Pitfalls:

• Orthogonality can drift: Gram-Schmidt based QR is convenient, but for badly conditioned matrices it is less robust than Householder QR.

• Not a general inverse: QR is excellent for least squares, but not the first choice for square exact solves unless you need its stability properties.

• Rank deficiency: Small diagonal values in R usually mean near-dependence. Treat them with tolerance, not exact comparisons.

SVD Decomposition¶

Mat::SVDDecomposition Mat::svd_decompose(int max_iter = 100, float tolerance = 1e-6f) const;

Description:

Compute Singular Value Decomposition: A = U * S * V^T. Most general decomposition, used for rank estimation, pseudo-inverse, dimension reduction. Uses iterative method based on eigendecomposition.

Mathematical Principle:

SVD decomposes any m × n matrix A into:

• U: m × m orthogonal matrix (left singular vectors)

• S: m × n diagonal matrix (singular values σ₁ ≥ σ₂ ≥ ... ≥ σᵣ ≥ 0)

• V: n × n orthogonal matrix (right singular vectors)

The singular values reveal the matrix's fundamental properties: rank, condition number, and numerical behavior.

Parameters:

• int max_iter : Maximum number of iterations (default: 100).

• float tolerance : Convergence tolerance (default: 1e-6).

Returns:

SVDDecomposition containing U, S, V matrices, rank, and status.

Usage Insights:

• Rank Estimation: The numerical rank is the number of singular values above the tolerance threshold.

• Pseudo-Inverse: A⁺ = V * S⁺ * U^T, where S⁺ has 1/σᵢ for non-zero σᵢ.

• Dimension Reduction: Truncated SVD (keeping only largest singular values) provides low-rank approximation.

• Condition Number: κ(A) = σ₁ / σᵣ, where σᵣ is the smallest non-zero singular value.

• Applications:

• Least squares for rank-deficient systems

• Principal Component Analysis (PCA)

• Image compression

• Noise reduction

Pitfalls:

• Cost: SVD is the most informative decomposition, but also one of the most expensive.

• Tolerance defines rank: Numerical rank is not absolute; it changes with the threshold you choose.

• Memory footprint: Full U and V storage can be heavy on embedded targets.

• Interpretation: Singular values tell you conditioning and energy concentration, but they do not by themselves solve a model-selection problem.

Solving Linear Systems Using Decompositions¶

Solve using LU Decomposition¶

static Mat Mat::solve_lu(const LUDecomposition &lu, const Mat &b);

Description:

Solve linear system Ax = b using precomputed LU decomposition. More efficient than solve() when solving multiple systems with the same coefficient matrix.

Mathematical Principle:

Given A = P * L * U, solve Ax = b by:

1. Solve Ly = Pb (forward substitution)
2. Solve Ux = y (back substitution)

Parameters:

• const LUDecomposition &lu : Precomputed LU decomposition.

• const Mat &b : Right-hand side vector (N×1).

Returns:

Mat - Solution vector (N×1).

Usage Insights:

• Multiple RHS: After computing LU decomposition once, solve multiple systems efficiently.

• Performance: O(n²) per solve vs O(n³) for full solve, significant savings for multiple RHS.

• Memory: Reuses the decomposition, avoiding repeated computation.

Pitfalls:

• Same matrix, many right-hand sides: This is where LU shines. If A changes every time, the gain shrinks quickly.

• Pivot information matters: The permutation matrix is part of the solve; ignoring it gives the wrong answer.

Solve using Cholesky Decomposition¶

static Mat Mat::solve_cholesky(const CholeskyDecomposition &chol, const Mat &b);

Description:

Solve linear system Ax = b using precomputed Cholesky decomposition. More efficient than LU for symmetric positive definite matrices.

Mathematical Principle:

Given A = L * L^T, solve Ax = b by:

1. Solve Ly = b (forward substitution)
2. Solve L^T x = y (back substitution)

Parameters:

• const CholeskyDecomposition &chol : Precomputed Cholesky decomposition.

• const Mat &b : Right-hand side vector (N×1).

Returns:

Mat - Solution vector (N×1).

Usage Insights:

• Efficiency: Faster than LU for SPD matrices, both in decomposition and solving.

• Stability: More numerically stable for SPD matrices.

• Applications: Structural dynamics, optimization, statistics.

Pitfalls:

• Only SPD: If the matrix is not SPD, this method is not just inefficient; it is invalid.

• Check before decomposing: A quick symmetry/positive-definite check is cheaper than debugging a failed factorization later.

Solve using QR Decomposition (Least Squares)¶

static Mat Mat::solve_qr(const QRDecomposition &qr, const Mat &b);

Description:

Solve linear system using QR decomposition. Provides least squares solution for overdetermined systems (more equations than unknowns).

Mathematical Principle:

For Ax ≈ b (overdetermined), the least squares solution minimizes ||Ax - b||₂. Using A = Q * R:

• x = R^(-1) * Q^T * b

This avoids the numerically unstable normal equations A^T * A * x = A^T * b.

Parameters:

• const QRDecomposition &qr : Precomputed QR decomposition.

• const Mat &b : Right-hand side vector (M×1, where M ≥ N).

Returns:

Mat - Least squares solution vector (N×1).

Usage Insights:

• Overdetermined Systems: Handles cases where there are more equations than unknowns.

• Numerical Stability: More stable than solving normal equations directly.

• Applications:

• Curve fitting

• Data regression

• Signal processing

Pitfalls:

• Least squares, not exact solve: If the system is inconsistent, QR gives the best fit, not a contradiction-free exact answer.

• Still needs rank awareness: Very small diagonal values in R mean the fit may be ill-conditioned or underdetermined.

Pseudo-Inverse¶

static Mat Mat::pseudo_inverse(const SVDDecomposition &svd, float tolerance = 1e-6f);

Description:

Compute the Moore-Penrose pseudo-inverse A⁺ using SVD decomposition. Works for rank-deficient or non-square matrices where the regular inverse doesn't exist.

Mathematical Principle:

For A = U * S * V^T, the pseudo-inverse is A⁺ = V * S⁺ * U^T, where S⁺ has 1/σᵢ for singular values σᵢ > tolerance, and 0 otherwise.

Properties of Pseudo-Inverse:

• A * A⁺ * A = A

• A⁺ * A * A⁺ = A⁺

• (A * A⁺)^T = A * A⁺

• (A⁺ * A)^T = A⁺ * A

Parameters:

• const SVDDecomposition &svd : Precomputed SVD decomposition.

• float tolerance : Threshold for singular values (default: 1e-6). Singular values below this are treated as zero.

Returns:

Mat - Pseudo-inverse matrix.

Usage Insights:

• Rank-Deficient Systems: Provides solution for systems where A is not full rank.

• Minimum Norm Solution: For underdetermined systems, gives the solution with minimum ||x||₂.

• Least Squares: For overdetermined systems, gives the least squares solution.

• Applications:

• Control systems

• Signal processing

• Machine learning (regularization)

Pitfalls:

• Tolerance is a modeling choice: The pseudo-inverse depends on the singular-value cutoff. Changing tolerance can change the effective solution.

• Not always sparse-friendly: A dense pseudo-inverse can destroy structure and inflate memory use.

• Prefer solve when possible: If the system is square and well-conditioned, direct solving is usually cheaper.