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] 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.
Figure 1. . 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] For a complex Hermitian positive-definite matrix, the corresponding factorization is , with the conjugate transpose and the diagonal of positive and real.[ref] The factorization turns a linear system into two triangular systems.
For example, matching entries in gives
Choose the positive square roots in order: , and . Therefore
The 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.
Figure 1. . Original exact two-by-two example. Multiplying the displayed lower triangular matrix by its transpose recovers ; forward substitution and backward substitution in a triangular system give for . Details.
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 with
The first pivot element is . The elimination multipliers are and . Apply the elementary row operations and to obtain
Now the second pivot element is , and . The elementary row operation gives
Importantly, 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 .
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]
The matrix splitting has , , so its iteration matrix is
Compute 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]
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]
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 splitting with invertible produces
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]
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 form
where 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 gives
a 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 proved
This completes the induction by triangularization using dual spaces.
Since every , the diagonal matrix is invertible. Set
The 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 , giving
This is the LDL decomposition.
Take the positive square roots of the diagonal entries:
The matrix is a lower triangular matrix with positive diagonal, and
This is a Cholesky decomposition. For a nonzero vector , the vector is nonzero because is invertible. Therefore
Here 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 is
The 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 gives
The 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 , so
Similar 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 therefore
Direct multiplication checks
The 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 satisfies
The 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 . Partition
Writing 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.
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 , giving
The second pivot element is , with elimination multiplier . Perform , then use the third pivot element and elimination multiplier to perform . The LU factorization is
Here 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 system
Although 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 .
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 gives
exactly 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.