André-Louis Cholesky Created 2026-10-18 Updated 2026-11-05
André-Louis Cholesky was a mathematician and military surveyor whose computational work grew out of geodesy. During the revision of French triangulation surveys, he developed the method now known as the Cholesky decomposition.[ref]
https://www.polytechnique.edu/bibliotheque/actualites/cholesky-andre-x1895. École polytechnique historical archive, A la section géodésique de l'armée
For a real symmetric positive-definite matrix , the Cholesky decomposition writes with lower triangular matrix and positive diagonal entries. A linear system can then be solved by two triangular substitutions. The link with surveying is natural: reconciling redundant, imperfect measurements by least squares leads to such systems when the design matrix has full column rank.Identity source: www.polytechnique.edu/bibliotheque/actualites/cholesky-andre-x1895
Portrait of André-Louis Cholesky
. Unknown photographer (Commons field: no conegut). Public domain (PD-old as declared on Commons). Source and usage information on the linked Wikimedia file page. Source. Cholesky decomposition Created 2026-09-28 Updated 2026-10-25
Every real symmetric matrix that is a positive-definite matrix has a unique Cholesky decomposition , where is a lower triangular matrix with positive diagonal entries.[ref]
DPOTRF (LAPACK 3.12.1) by LAPACK project. DPOTRF, LAPACK 3.12.1: Purpose, lower-triangular factorization; UPLO and INFO argument documentation
For a complex Hermitian positive-definite matrix, the corresponding factorization is , with the conjugate transpose and the diagonal of positive and real.[ref]ZPOTRF (LAPACK 3.12.1) by LAPACK project. ZPOTRF, LAPACK 3.12.1: Purpose, Hermitian factorization with conjugate transpose
The factorization turns a linear system into two triangular systems.For example, matching entries in givesChoose the positive square roots in order: , and . ThereforeThe second radicand, , is the scalar Schur complement remaining after eliminating the first variable. A positive diagonal of alone would not guarantee that this quantity is positive; Sylvester's criterion explains the relevant leading principal minors.
For that right-hand side, first use forward substitution to solve : and , so . Then use backward substitution in a triangular system to solve : and , giving . No matrix inverse is required. For several right-hand sides with the same , keep the same factor and repeat only these triangular solves.
Computed three-by-three LU decomposition Created 2026-10-14 Updated 2026-10-31
Here is an original example of Gaussian elimination that can be checked entirely with integers. Start withThe first pivot element is . The elimination multipliers are and . Apply the elementary row operations and to obtainNow the second pivot element is , and . The elementary row operation givesImportantly, comes from the current reduced matrix, not the original entry . To reconstruct from , its second row is twice the first row of plus the second; its third row is minus the first row of plus three times the second plus the third. These are exactly the coefficients in the lower triangular matrix .
Gaussian elimination turns a three-by-three matrix into an upper triangular matrix, while elimination multipliers are stored in a lower triangular matrix
. Details. The unpivoted calculation here illustrates the structure of LU decomposition; applying partial pivoting to this particular would choose a different first row and different factors. Both procedures can represent the same linear system when their permutation matrices are handled consistently.
Diagonal matrix entry Created 2026-10-18 Updated 2026-10-25
A diagonal matrix entry uses the same row and column index. The diagonal entries of a lower triangular matrix determine its determinant by their product.
Forward substitution 2026-10-14
Forward substitution solves a linear system with a lower triangular matrix and nonzero diagonal, from the first row downward:The first sum is empty, so . Each subsequent row uses only components already found. For a unit triangular matrix , omit the division because . The complementary procedure, back substitution, solves an upper triangular matrix system from the last row upward.[ref]
Gauss-Seidel method Created 2026-10-05 Updated 2026-10-18
With , the Gauss-Seidel method solves . It uses each newly computed coordinate immediately in the remaining triangular solve. It is successive over-relaxation with parameter one, and converges from every start when is a symmetric positive-definite matrix.
To solve a linear system , split its matrix into the diagonal matrix and the strictly lower and upper parts , retaining their actual signs. Every diagonal entry must be nonzero. A sweep updates coordinates in order:The Jacobi method uses the old vector throughout; the Gauss-Seidel method reuses each fresh coordinate. Thus the update order is part of the method.[ref]
Notes for Ma221 Lecture 21, Nov 10 2009. Math 221, Matrix Computations, Fall 2009; Gauss-Seidel: component loop, splitting and update-order remark, web lines 372–380
The matrix splitting has , , so its iteration matrix isCompute the update by forward substitution with the lower triangular matrix , rather than forming a matrix inverse. For a sparse matrix, a sweep only visits stored entries and can overwrite the working vector. Dependencies on fresh coordinates constrain parallel execution. This makes the Gauss-Seidel method an inexpensive solver component, including in a multigrid method.[ref]
Notes for Ma221 Lecture 21, Nov 10 2009. Math 221, Matrix Computations, Fall 2009; Basic iterative methods, cheap splitting requirement; Gauss-Seidel loop and red-black ordering; opening list of solver building blocks
Laplace expansion Created 2026-10-18 Updated 2026-10-20
Expand the determinant along any fixed row using the displayed formula, where deletes row and column . The signed matrix minors are its cofactors. Applying the same formula to the transpose gives the column expansion. The row formula follows from the Leibniz formula for determinants by grouping permutations according to the selected column in row . For a lower triangular matrix, expanding along the first row reduces the determinant to times a smaller triangular determinant; induction gives the product of all diagonal entries, even when some vanish.
LU decomposition Created 2026-09-29 Updated 2026-10-14
An LU decomposition writes a square matrix as , with a lower triangular matrix and an upper triangular matrix. We use a unit triangular matrix , meaning . The factors record Gaussian elimination: contains the reduced equations, while stores the elimination multipliers needed to reconstruct the original equations. This makes Reusing LU factors for multiple right-hand sides inexpensive.[ref]
For an invertible matrix, unpivoted LU decomposition requires nonzero pivot elements at every elimination step; invertibility alone does not guarantee this. In general use Pivoted LU decomposition, with the explicit convention and a permutation matrix . Once the factors exist, solving a linear system needs forward substitution followed by back substitution, rather than forming a matrix inverse.
Matrix decomposition 2026-10-14
A matrix decomposition expresses a matrix as a matrix product of factors with simpler structure. An LU decomposition uses triangular factors, a QR decomposition uses an orthogonal matrix and a triangular factor, and a singular value decomposition separates input directions, stretching factors and output directions. The factorization is an identity; the chosen structure determines what it helps compute.
For an original two-by-two example,The outer factors are a lower triangular matrix and an upper triangular matrix, each with unit diagonal as in a unit triangular matrix; the middle factor is a diagonal matrix. Their matrix product is , and the final diagonal entry is the Schur complement of the first entry . This is the scalar version of the factorization under block elimination.
To solve , solve with the three factors in their displayed order from left to right: first the lower triangular factor, then the diagonal factor, then the upper triangular factor. The intermediate vectors are and , giving the concise answer . A new right-hand side uses the same factors. This is why Reusing LU factors for multiple right-hand sides is useful: factorization records work that otherwise would be repeated, while forward substitution and back substitution apply it cheaply.[ref]
Convex Optimization by Stephen Boyd and Lieven Vandenberghe. §C.4.1, printed pp. 673–674: block factor-solve interpretation
Different decompositions reveal different features. A Cholesky decomposition uses positive definiteness of a symmetric matrix, while a singular value decomposition also applies to rectangular matrices and exposes matrix rank. The existence and numerical behaviour of the factors depend on their own hypotheses; for example, an unpivoted LU decomposition can fail when a chosen pivot element is zero even though the matrix is an invertible matrix.
Matrix splitting Created 2026-09-24 Updated 2026-10-18
A matrix splitting chooses an easy part of a matrix and corrects for the remaining part . To solve the linear system , first rewrite it as , then use the old vector on the right and solve for the new vector on the left. The matrix inverse notation describes the update mathematically; each step should solve with rather than construct its matrix inverse.[ref]
Notes for Ma221 Lecture 21, Nov 10 2009. Math 221, Matrix Computations, Fall 2009; Basic iterative methods: definition A = M − K and recurrence (*), with the two conditions for a useful splitting
The same identity also yields a correction based on the linear-system residual:This connects a matrix splitting to matrix preconditioning: rescales the imbalance before adding a correction. Choosing gives and the exact solution in one step, but requires solving the original linear system; choosing cheaply matters as much as the convergence criterion for a stationary iteration.
With the sign convention , is the diagonal matrix of diagonal entries and contain the actual lower and upper entries, without reversing their signs. The Jacobi method uses , ; the Gauss-Seidel method uses , . The first requires only division by diagonal entries; the second uses a lower triangular matrix and forward substitution. These choices explain their different information flows.
We prove by mathematical induction on that every linear map on an -dimensional complex vector space admits both triangular forms. For or this is immediate. When , the fundamental theorem of algebra supplies an eigenvector ; its span is an invariant subspace.
By parts (ii) and (iii), has dimension and is an invariant subspace of the dual map . Apply the induction hypothesis to to choose a basis with an upper triangular matrix. Extend it by a linear functional satisfying . Such a functional exists by basis extension, and cannot belong to . In this basis of , the matrix of has the formwhere is an upper triangular matrix; here the last diagonal entry follows by evaluating at , where .
There is a basis of with , since the dual basis construction gives a finite-dimensional bidual evaluation map isomorphism. In particular . The matrix transpose relationship for a dual map givesa lower triangular matrix. Reversing this basis turns a lower triangular matrix into an upper triangular matrix. If is the permutation matrix reversing the order, we have provedThis completes the induction by triangularization using dual spaces.
Since every , the diagonal matrix is invertible. SetThe division is by the diagonal entry of row . Thus is an upper triangular matrix with diagonal entries one and .
To use symmetry without assuming the conclusion, multiply by the matrix inverse of and its transpose:The matrix is symmetric because is a symmetric matrix. It is also an upper triangular matrix, since every factor on the right is an upper triangular matrix. A symmetric upper triangular matrix is diagonal. Its diagonal is exactly , because and have diagonal entries one. Hence , givingThis is the LDL decomposition.
Take the positive square roots of the diagonal entries:The matrix is a lower triangular matrix with positive diagonal, andThis is a Cholesky decomposition. For a nonzero vector , the vector is nonzero because is invertible. ThereforeHere is the Euclidean norm. Thus is a positive-definite matrix.
Use Gaussian elimination with partial pivoting, choosing a largest absolute entry in the uneliminated part of each pivot column. The first pivot is two in row one, with multipliers for rows two, three and four. The second-column candidates are then , so row two stays as pivot and the remaining nonzero multiplier is for row four. The augmented matrix at this stage isThe third pivot is three, so swap rows three and four. The final multiplier is , leaving last pivot . Carrying that swap through the previously stored multipliers givesThe transformed right side is . Back substitution gives , , and , hence .
For the eigenvalue claim, put . It is invertible because a permutation matrix and a unit lower triangular matrix are invertible. The factorization gives , soSimilar matrices have the same characteristic polynomial and therefore the same eigenvalues with algebraic multiplicities. This proves the claim without changing the prescribed elimination method.
The elimination multipliers are , , and . The unit lower triangular matrix in the LU factorization is thereforeDirect multiplication checksThe positive signs in undo the row subtractions used during Gaussian elimination.
The following R code implements the given recursion for the Cholesky decomposition. The input is assumed to be a real symmetric matrix that is a positive-definite matrix.
mychol <- function(A) {
n <- nrow(A)
if (n == 1L) return(matrix(sqrt(A[1L, 1L]), 1L, 1L))
j <- 2L:n
L <- matrix(0, n, n)
L[1L, 1L] <- sqrt(A[1L, 1L])
b <- A[j, 1L, drop = FALSE] / L[1L, 1L]
L[j, 1L] <- b
S <- A[j, j, drop = FALSE] - tcrossprod(b)
L[j, j] <- mychol(S)
L
}Here
tcrossprod(b) is the outer product . The drop = FALSE arguments preserve matrix dimensions, including the final one-by-one subproblem, and the base case also returns a matrix.To verify correctness, partition . The code sets , and recursively factors . This Schur complement is a positive-definite matrix: for ,If is the recursive answer, then the returned lower triangular matrix satisfiesThe one-dimensional identity starts this mathematical induction proof. Each recursive call reduces the matrix size by one and therefore terminates.
An inverse matrix exists only if the lower triangular matrix has no zero diagonal entries; this necessary assumption is implicit in the request for . PartitionWriting its inverse in the same block form and multiplying gives the recursive inverse of a lower triangular matrix:Indeed the lower-left block of the product is , and the diagonal blocks are identities. This gives an implementation with no call to
solve:ltinv <- function(L) {
n <- nrow(L)
if (n == 1L) return(matrix(1 / L[1L, 1L], 1L, 1L))
j <- 2L:n
Cinv <- ltinv(L[j, j, drop = FALSE])
B <- matrix(0, n, n)
B[1L, 1L] <- 1 / L[1L, 1L]
B[j, j] <- Cinv
B[j, 1L] <- -as.vector(Cinv %*% L[j, 1L, drop = FALSE]) / L[1L, 1L]
B
}The preserved matrix dimensions also handle . A mathematical induction using the block identity proves that every returned matrix is the desired inverse. Applied to the Cholesky decomposition from (a), all diagonal entries are positive, so all divisions are valid.
For either an upper triangular matrix or a lower triangular matrix,In the permutation formula for the determinant, every contribution except the diagonal product includes an entry on the zero side of the diagonal. Equivalently, successive cofactor expansions along the first column of an upper triangular matrix prove the formula by mathematical induction.
Let column of be , so . Every diagonal entry of the invertible lower triangular matrix is nonzero, since its determinant is the product of those entries. In row one through row , the right-hand side is zero. Forward substitution first gives , then , and successively for all . These are precisely the entries above the diagonal in column of . Every column satisfies this, proving that is lower triangular. This is triangularity of an inverse matrix.
Apply Gaussian elimination to the supplied nonsymmetric matrix,The first pivot element is , with elimination multiplier . Perform , givingThe second pivot element is , with elimination multiplier . Perform , then use the third pivot element and elimination multiplier to perform . The LU factorization isHere is a unit lower triangular matrix, hence an invertible matrix, and is an upper triangular matrix. No row interchange is needed: all three pivot elements used in elimination are nonzero. The final diagonal entry is zero, with an entirely zero last row; there is no further division by it. The three nonzero rows of are linearly independent, so the matrix rank is and is a singular matrix. Direct matrix multiplication gives .
Choose a fixed time increment , let , and approximate by . The Backward Euler method evaluates the transport term at the new time:Rearrange to the linear systemAlthough is a singular matrix, is an invertible matrix for every . In each row, its diagonal entry is and the sum of the absolute values of its off-diagonal entries is ; hence it is a strictly diagonally dominant matrix. Its Gershgorin discs have centre and radius , excluding zero. In particular, the zero eigenvalue of becomes an eigenvalue of .
- Form once. Compute a Pivoted LU decomposition , storing the permutation matrix , the unit lower triangular matrix , and the upper triangular matrix .
- Set .
- At each step form . Solve by forward substitution, then solve by back substitution.
- Advance and repeat until the desired final time. Choose the fixed increment so the final time is an integer multiple of , or refactor if a shorter final step is used.
Both triangular solves use the factors of . The LU factorization of the singular transport matrix from part (a) is not a solver for this new linear system. No matrix inverse needs to be formed. Since and are fixed, the same LU factorization serves every step; only the right-hand side changes. A varying would give without changing these factors.
As a check, multiplying the update by givesexactly reproducing the inventory balance at the sampled times, apart from roundoff error. The Backward Euler method has temporal global error for smooth solutions, so its A-stability does not remove the need to select for accuracy.


