mathematics//linear algebra//matrix factorization

A matrix factorization is the rewriting of a matrix as a product of simpler matrices (triangular, orthogonal, diagonal), and it is how every serious linear algebra routine actually works: to solve, to invert, to fit or to read a matrix, a library first breaks it into pieces whose own problems are trivial. A triangular system is solved by substitution, one unknown at a time; an orthogonal matrix is inverted by transposing it; a diagonal one by dividing. The cost goes into the factorization, of the order of \(n^3\) operations, and it is paid once: every later right-hand side costs only \(n^2\).

A matrix factorization is the rewriting of a matrix as a product of simpler matrices (triangular, orthogonal, diagonal), and it is how every serious linear algebra routine actually works: to solve, to invert, to fit or to read a matrix, a library first breaks it into pieces whose own problems are trivial. A triangular system is solved by substitution, one unknown at a time; an orthogonal matrix is inverted by transposing it; a diagonal one by dividing. The cost goes into the factorization, of the order of n3n^3n3 operations, and it is paid once: every later right-hand side costs only n2n^2n2.

The members differ in what they demand of the matrix and what they give back, and choosing among them is a short list of situations. A general square system calls for the LU decomposition, which is what solve runs and why it beats forming an inverse. A symmetric positive definite matrix (a covariance, a stiffness matrix, the normal matrix of a fit) takes the Cholesky decomposition, at half the cost of LU, and it is the factor behind the sigma points of an unscented Kalman filter and the square-root forms of the Kalman filter. An overdetermined fit takes the QR decomposition, the stable way to do least squares without squaring the condition number. When the question is how close the matrix is to losing rank, which directions it amplifies and which the data cannot see, the answer is the SVD, the most informative and most expensive of all. And when the matrix is a dynamics, a Markov transition or a covariance whose behaviour along its own directions is wanted, the eigendecomposition decouples it into independent scalar modes.

A factorization is chosen by structure, and structure is free speed and accuracy.

Telling the library that a matrix is symmetric, positive definite, banded or sparse (cho_solve, eigh, solve_banded, scipy.sparse) routinely buys a factor of two to a hundred in time and keeps results real and well behaved.

On embedded hardware the choice is also about survival in single precision. A Kalman covariance updated in float32 can lose its symmetry or its positive definiteness and make the filter diverge; the UD factorization, which stores the covariance as UDUTUDU^{\mathsf T}UDUT with UUU unit upper triangular and DDD diagonal, and its square-root relatives keep it positive definite by construction, at a modest extra cost per step.

Behind np.linalg and scipy.linalg sit LAPACK and an optimized BLAS (OpenBLAS, MKL), decades of tuning that no hand-written loop matches; on a microcontroller, CMSIS-DSP offers small dense routines in float32, enough for a Kalman filter and rarely for an SVD.