Principle:Haifengl Smile Matrix Decomposition

Overview

Matrix Decomposition (also called matrix factorization) is the process of expressing a matrix as a product of simpler, structured matrices. Each decomposition reveals different structural properties of the original matrix and enables different computational applications. The five principal decompositions in Smile are LU, QR, SVD, Eigenvalue (EVD), and Cholesky.

Matrix decomposition is the central stage of the Matrix_Decomposition_Pipeline workflow. It transforms a constructed matrix into a factored form that can then be used for solving linear systems, computing determinants, finding eigenvalues, performing dimensionality reduction, and more.

Theoretical Basis

LU Decomposition

The LU decomposition factors a matrix into a lower triangular matrix and an upper triangular matrix with row pivoting:

$PA=LU$

where $P$ is a permutation matrix, $L$ is unit lower triangular (ones on diagonal), and $U$ is upper triangular.

Properties:

• Exists for any square matrix (with pivoting)
• Computational cost: $\frac{2}{3}{n}^3$ flops for an $n\times n$ matrix
• The determinant is $\det (A)=\det (P{)}^{-1}\cdot \prod _{i=1}^n{u}_{ii}$
• If $u_{ii}=0$ for some $i$, the matrix is singular (the factorization still completes, but the matrix cannot be used for solving)

Applications: Solving $Ax=b$, computing determinants, matrix inversion.

QR Decomposition

The QR decomposition factors a matrix into an orthogonal matrix and an upper triangular matrix:

$A=QR$

where $Q\in {\mathbb{ℝ}}^{m\times n}$ has orthonormal columns ($Q^{T}Q=I$) and $R\in {\mathbb{ℝ}}^{n\times n}$ is upper triangular.

Properties:

• Always exists for any $m\times n$ matrix with $m\ge n$
• Computed via Householder reflections (LAPACK dgeqrf): $\frac{4}{3}{n}^{2}m-\frac{2}{3}{n}^3$ flops
• $Q$ is stored implicitly as a product of Householder reflectors $Q=H_1{H}_2\cdots H_k$ where $H_i=I-{\tau }_i{v}_i{v}_i^T$
• The explicit $Q$ can be reconstructed via dorgqr

Applications: Least squares problems $\min \Vert Ax-b{\Vert }_2$, eigenvalue algorithms (QR algorithm), orthogonalization.

Singular Value Decomposition (SVD)

The SVD factors any $m\times n$ matrix into three matrices:

$A=U\mathrm{\Sigma }{V}^T$

where $U\in {\mathbb{ℝ}}^{m\times k}$ and $V\in {\mathbb{ℝ}}^{n\times k}$ have orthonormal columns, $\mathrm{\Sigma }\in {\mathbb{ℝ}}^{k\times k}$ is diagonal with nonnegative entries ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_k\ge 0$, and $k=\min (m,n)$.

Properties:

• Always exists for any matrix
• Compact SVD: for $m>n$, only the first $n$ columns of $U$ are computed
• The singular values measure the "importance" of each component
• The matrix rank equals the number of nonzero singular values
• Condition number: $\kappa (A)={\sigma }_1/{\sigma }_k$

The Eckart-Young Theorem: The best rank-$r$ approximation to $A$ in the Frobenius or operator norm is:

$A_r=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^T$

Applications: Dimensionality reduction (PCA), pseudoinverse computation, matrix rank determination, data compression, latent semantic analysis.

Eigenvalue Decomposition (EVD)

For a square matrix $A\in {\mathbb{ℝ}}^{n\times n}$:

$Av=\lambda v$

where $\lambda$ is an eigenvalue and $v$ is the corresponding eigenvector.

Symmetric case: If $A=A^T$, then: $A=V\mathrm{\Lambda }{V}^T$

where $V$ is orthogonal and $\mathrm{\Lambda }=\text{diag}({\lambda }_1,\dots ,{\lambda }_n)$ with all real eigenvalues.

General case: Eigenvalues may be complex. LAPACK computes the real and imaginary parts separately: ${\lambda }_i=w_{r,i}+j\cdot w_{i,i}$. Left eigenvectors satisfy $v^{T}A=\lambda v^T$.

Properties:

• Symmetric matrices have real eigenvalues and orthogonal eigenvectors
• Smile uses LAPACK dsyevd for symmetric and dgeev for general matrices
• dsyevd uses divide-and-conquer, which is faster than dsyev for large matrices

Applications: PCA (eigendecomposition of covariance matrix), spectral clustering, stability analysis, Markov chain analysis.

Cholesky Decomposition

For a symmetric positive definite (SPD) matrix:

$A=L{L}^T$

where $L$ is lower triangular with positive diagonal elements.

Properties:

• Exists if and only if $A$ is SPD
• Computational cost: $\frac{1}{3}{n}^3$ flops -- half the cost of LU
• If dpotrf returns a nonzero info code, the matrix is not positive definite
• The determinant is $\det (A)={\left(\prod _{i=1}^n{l}_{ii}\right)}^2$

Applications: Efficient solving of SPD systems, Monte Carlo simulation with correlated variables, Kalman filters (unscented transform).

Choosing the Right Decomposition

Decomposition	Matrix Requirements	Best For	Cost
LU	Square	Solving $Ax=b$, determinants	$\frac{2}{3}{n}^3$
Cholesky	Symmetric positive definite	Solving SPD systems (2x faster than LU)	$\frac{1}{3}{n}^3$
QR	$m\ge n$ (overdetermined)	Least squares, numerical stability	$\frac{4}{3}{n}^{2}m$
SVD	Any	Rank, pseudoinverse, dimensionality reduction	$O(mn\min (m,n))$
EVD	Square (symmetric preferred)	Spectral analysis, PCA	$O(n^3)$

Destructive Nature of Decompositions

An important implementation detail: all decompositions in Smile overwrite the input matrix. The caller must make a copy before decomposing if the original matrix is still needed:

// The matrix A will be overwritten by lu()
DenseMatrix A_copy = A.copy();
LU lu = A_copy.lu();
// A_copy now contains the LU factors, not the original data

This design avoids unnecessary memory allocation when the original matrix is no longer needed.

Relationship to the Pipeline

Construction --> Arithmetic --> Decomposition --> Solving --> Result Extraction
                                     ^
                                     |
                              (this principle)

Related

• Implementation:Haifengl_Smile_DenseMatrix_Decomposition
• Principle:Haifengl_Smile_Matrix_Arithmetic
• Principle:Haifengl_Smile_Linear_System_Solving
• Principle:Haifengl_Smile_Decomposition_Result_Extraction

Knowledge Sources

• Smile
• LAPACK Users' Guide
• Golub & Van Loan, Matrix Computations

Domains

Linear_Algebra, Numerical_Computing, Matrix_Theory