Symmetric positive definite matrices constitute the most mathematically elegant and computationally well-behaved class of matrices in numerical linear algebra. They arise throughout computational science, including the discretization of elliptic partial differential equations, Hessian matrices in convex optimization, covariance and information matrices in statistics, and the normal equations of linear least squares.
When a matrix possesses both symmetry and positive definiteness, standard numerical algorithms simplify dramatically. Gaussian elimination requires no pivoting to ensure backward stability, memory requirements are halved because only one triangular factor must be stored, and the total computational work is reduced by fifty percent compared to the general nonsymmetric case.
Definitions and Characterizations
The algebraic characterization of positive definiteness is stated in terms of a quadratic form. For a symmetric matrix \(\bA\), the scalar mapping \(\bx \mapsto \bx^T\bA\bx\) assigns an energy value to every direction in \(\fR^n\).
A matrix \(\bA \in \fR^{n \times n}\) is symmetric positive definite if \(\bA = \bA^T\) and the quadratic form is strictly positive for every nonzero vector \(\bx \in \fR^n\): \[
\begin{align}
\bx^T\bA\bx > 0 \qquad \text{for all } \bx \neq \bzero.
\end{align}
\] If the inequality is non-strict, so that \(\bx^T\bA\bx \geq 0\) for all \(\bx \in \fR^n\), the matrix is termed symmetric positive semidefinite.
Let \(\bA \in \fR^{n \times n}\) be a real symmetric matrix. The following statements are mathematically equivalent.
The matrix \(\bA\) is positive definite (\(\bx^T\bA\bx > 0\) for all \(\bx \neq \bzero\)).
All eigenvalues of \(\bA\) are strictly positive: \(\lambda_i(\bA) > 0\) for each \(i = 1, ..., n\).
All leading principal minors are positive: \(\det(\bA[1:k, 1:k]) > 0\) for each \(k = 1, ..., n\).
There exists a unique lower triangular matrix \(\bL \in \fR^{n \times n}\) with strictly positive diagonal entries such that \(\bA = \bL\bL^T\).
There exists a nonsingular matrix \(\bW \in \fR^{n \times n}\) such that \(\bA = \bW^T\bW\).
We demonstrate the equivalence between the quadratic form, the eigenvalues, and the existence of a factorization. By the Spectral Theorem, any real symmetric matrix has an orthogonal eigendecomposition \(\bA = \bQ\boldsymbol{\Lambda}\bQ^T\), where \(\bQ\) is orthogonal and \(\boldsymbol{\Lambda} = \operatorname{diag}(\lambda_1, ..., \lambda_n)\) contains the real eigenvalues. For any vector \(\bx \neq \bzero\), let \(\by = \bQ^T\bx\). Since \(\bQ\) is nonsingular, \(\by \neq \bzero\). The quadratic form becomes \[
\begin{align}
\bx^T\bA\bx = \bx^T\bQ\boldsymbol{\Lambda}\bQ^T\bx = \by^T\boldsymbol{\Lambda}\by = \sum_{i=1}^n \lambda_i y_i^2.
\end{align}
\] If all \(\lambda_i > 0\), this sum is strictly positive for any \(\by \neq \bzero\), proving that positive eigenvalues imply positive definiteness. Conversely, if \(\lambda_j \leq 0\) for some index \(j\), choosing \(\bx = \bq_j\) (the \(j\)-th eigenvector) gives \(\bx^T\bA\bx = \lambda_j \|\bq_j\|_2^2 \leq 0\), showing that positive definiteness implies all eigenvalues are positive.
To establish the equivalence with the factorization \(\bA = \bW^T\bW\), let \(\bW = \boldsymbol{\Lambda}^{1/2}\bQ^T\), where \(\boldsymbol{\Lambda}^{1/2} = \operatorname{diag}(\sqrt{\lambda_1}, ..., \sqrt{\lambda_n})\). This matrix is well-defined and nonsingular because each \(\lambda_i > 0\). Then \(\bW^T\bW = \bQ\boldsymbol{\Lambda}^{1/2}\boldsymbol{\Lambda}^{1/2}\bQ^T = \bQ\boldsymbol{\Lambda}\bQ^T = \bA\). Conversely, if \(\bA = \bW^T\bW\) with \(\bW\) nonsingular, then for every \(\bx \neq \bzero\), we have \(\bx^T\bA\bx = (\bW\bx)^T(\bW\bx) = \|\bW\bx\|_2^2 > 0\).
(Practical tests for positive definiteness) Although Sylvester’s criterion on leading principal minors is theoretically elegant, computing \(n\) determinants is numerically impractical. In scientific software, the standard numerical test for positive definiteness is attempting to compute the Cholesky factorization. If the algorithm successfully completes with all positive square roots, the matrix is guaranteed to be positive definite; if a non-positive pivot is encountered, the matrix is indefinite or singular.
Geometry of the Energy Norm
Every symmetric positive definite matrix induces a natural inner product and a corresponding geometry on \(\fR^n\). The level sets of the quadratic form define ellipsoids whose principal axes align with the eigenvectors of \(\bA\).
Let \(\bA \in \fR^{n \times n}\) be symmetric positive definite. The \(\bA\)-inner product of two vectors \(\bx, \by \in \fR^n\) is defined by \[
\begin{align}
\langle \bx, \by \rangle_\bA = \bx^T\bA\by.
\end{align}
\] The induced norm, known as the \(\bA\)-norm or energy norm, is given by \[
\begin{align}
\|\bx\|_\bA = \sqrt{\bx^T\bA\bx}.
\end{align}
\] The unit sphere in this norm, \(\{\bx \in \fR^n : \|\bx\|_\bA = 1\}\), forms an \(n\)-dimensional ellipsoid whose semi-axis lengths are \(1/\sqrt{\lambda_i}\), oriented along the eigenvectors \(\bq_i\) of \(\bA\).
(Geometric interpretation of the energy ellipsoid) Consider the diagonal positive definite matrix \[
\begin{align}
\bA = \begin{pmatrix} 4 & 0 \\ 0 & 1 \end{pmatrix}.
\end{align}
\] The unit sphere in the \(\bA\)-norm satisfies the equation \(4x_1^2 + x_2^2 = 1\). This describes an ellipse in the plane with semi-major axis of length \(1\) along the \(x_2\)-direction and semi-minor axis of length \(1/2\) along the \(x_1\)-direction. The direction corresponding to the larger eigenvalue \(\lambda_1 = 4\) is geometrically stiffer, meaning that moving along \(x_1\) incurs a higher energy penalty than moving along \(x_2\).
The Cholesky Factorization Algorithm
Because an SPD matrix satisfies \(\bA = \bA^T\), the lower and upper triangular factors in Gaussian elimination must mirror each other: \(\bU = \bL^T\) after accounting for diagonal scalings. The Cholesky algorithm computes this single triangular factor directly.
Every symmetric positive definite matrix \(\bA \in \fR^{n \times n}\) has a unique factorization of the form \[
\begin{align}
\bA = \bL\bL^T,
\end{align}
\] where \(\bL \in \fR^{n \times n}\) is lower triangular with strictly positive diagonal entries \(\ell_{jj} > 0\).
We derive the recursive formulas for the entries of \(\bL\) by equating entry \((i,j)\) of \(\bA\) with the product \(\bL\bL^T\). Because \(\bL\) is lower triangular, the product yields \[
\begin{align}
a_{ij} = \sum_{k=1}^n \ell_{ik} (\bL^T)_{kj} = \sum_{k=1}^{\min(i,j)} \ell_{ik} \ell_{jk}.
\end{align}
\] For the diagonal entry where \(i = j\), the equation becomes \[
\begin{align}
a_{jj} = \sum_{k=1}^{j-1} \ell_{jk}^2 + \ell_{jj}^2 \implies \ell_{jj} = \sqrt{a_{jj} - \sum_{k=1}^{j-1} \ell_{jk}^2}.
\end{align}
\] Because \(\bA\) is positive definite, the quantity under the radical is strictly positive at every step. For the strictly subdiagonal entries in column \(j\) where \(i > j\), we obtain \[
\begin{align}
a_{ij} = \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} + \ell_{ij} \ell_{jj} \implies \ell_{ij} = \frac{1}{\ell_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right).
\end{align}
\] Evaluating these formulas column by column for \(j = 1, ..., n\) computes the factor \(\bL\) uniquely.
Let \(\bA \in \fR^{n \times n}\) be a symmetric positive definite matrix. The lower triangular factor \(\bL\) is computed via the following column-oriented algorithm.
For \(j = 1, 2, ..., n\): \begin{enumerate}
Compute the diagonal entry: \[
\begin{align}
v = a_{jj} - \sum_{k=1}^{j-1} \ell_{jk}^2.
\end{align}
\] If \(v \leq 0\), the matrix is not positive definite and the algorithm terminates. Otherwise, set \(\ell_{jj} = \sqrt{v}\).
For each row index \(i = j+1, ..., n\), compute the subdiagonal entries: \[
\begin{align}
\ell_{ij} = \frac{1}{\ell_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right).
\end{align}
\]
\end{enumerate}
(Step-by-step calculation of a \(3 \times 3\) Cholesky factor) Consider the symmetric positive definite matrix \[
\begin{align}
\bA = \begin{pmatrix}
4 & 2 & -2 \\
2 & 10 & 2 \\
-2 & 2 & 6
\end{pmatrix}.
\end{align}
\] We compute the entries of the lower triangular factor \(\bL\) column by column.
First column (\(j = 1\)): \[
\begin{align}
\ell_{11} &= \sqrt{a_{11}} = \sqrt{4} = 2, \\
\ell_{21} &= \frac{a_{21}}{\ell_{11}} = \frac{2}{2} = 1, \\
\ell_{31} &= \frac{a_{31}}{\ell_{11}} = \frac{-2}{2} = -1.
\end{align}
\]
Second column (\(j = 2\)): \[
\begin{align}
\ell_{22} &= \sqrt{a_{22} - \ell_{21}^2} = \sqrt{10 - 1^2} = \sqrt{9} = 3, \\
\ell_{32} &= \frac{1}{\ell_{22}} \left( a_{32} - \ell_{31}\ell_{21} \right) = \frac{1}{3} \left( 2 - (-1)(1) \right) = \frac{3}{3} = 1.
\end{align}
\]
Third column (\(j = 3\)): \[
\begin{align}
\ell_{33} &= \sqrt{a_{33} - \ell_{31}^2 - \ell_{32}^2} = \sqrt{6 - (-1)^2 - 1^2} = \sqrt{4} = 2.
\end{align}
\]
The complete Cholesky factor is \[
\begin{align}
\bL = \begin{pmatrix}
2 & 0 & 0 \\
1 & 3 & 0 \\
-1 & 1 & 2
\end{pmatrix}.
\end{align}
\] Direct matrix multiplication confirms that \(\bL\bL^T = \bA\).
Stability and Complexity Analysis
The Cholesky algorithm has outstanding numerical stability properties that distinguish it from general LU factorization.
For a symmetric positive definite matrix \(\bA \in \fR^{n \times n}\), the Cholesky factorization requires asymptotically \(\frac{1}{3}n^3\) floating-point operations and \(n\) square roots. This represents half the computational work and half the memory storage of a standard LU factorization.
At column \(j\), computing the diagonal entry requires \(j-1\) multiplications and subtractions followed by one square root. Computing the \(n-j\) subdiagonal entries requires \((n-j)(j-1)\) multiplications and subtractions, and \(n-j\) divisions. Summing the dominant multiplicative terms across all columns gives \[
\begin{align}
\sum_{j=1}^n 2(n-j)j = 2n \sum_{j=1}^n j - 2 \sum_{j=1}^n j^2
= 2n \left( \frac{n^2}{2} \right) - 2 \left( \frac{n^3}{3} \right)
= \frac{1}{3}n^3 \text{ flops}.
\end{align}
\] Solving \(\bA\bx = \bb\) using \(\bL\by = \bb\) and \(\bL^T\bx = \by\) requires \(2n^2\) flops, identical to standard triangular substitution.
(Inherent stability without pivoting) A critical property of the Cholesky decomposition is that row interchanges are entirely unnecessary. For every entry in the factor, the diagonal identity \(a_{ii} = \sum_{k=1}^i \ell_{ik}^2\) implies that \(\ell_{ik}^2 \leq a_{ii}\) for all \(k \leq i\). Consequently, the entries of the factor \(\bL\) are bounded directly by the diagonal entries of the original matrix: \[
\begin{align}
|\ell_{ik}| \leq \sqrt{a_{ii}} \leq \sqrt{\max_m a_{mm}}.
\end{align}
\] There is no possibility of element growth during factorization. The algorithm is unconditionally backward stable, satisfying \(\|\delta\bA\|_2 / \|\bA\|_2 \leq c n \varepsilon_{\text{mach}}\) for a small constant \(c\).
(Condition number of SPD matrices) For any symmetric positive definite matrix \(\bA\), the condition number in the Euclidean norm simplifies to the ratio of extreme eigenvalues: \[
\begin{align}
\kappa_2(\bA) = \frac{\lambda_{\max}(\bA)}{\lambda_{\min}(\bA)}.
\end{align}
\] Furthermore, the Cholesky factor satisfies \(\kappa_2(\bL) = \sqrt{\kappa_2(\bA)}\), showing that the triangular solves operate on matrices with significantly milder condition numbers than the original operator.
Compute the Cholesky factorization by hand for the matrix \[
\begin{align}
\bA = \begin{pmatrix} 9 & 3 & -6 \\ 3 & 26 & -7 \\ -6 & -7 & 9 \end{pmatrix}.
\end{align}
\]
Using the Cholesky factor computed above, solve \(\bA\bx = (6, 44, 2)^T\) via forward and backward substitution.
Let \(\bA \in \fR^{n \times n}\) be an SPD matrix. Prove that the largest entry of \(\bA\) in absolute value must lie on the main diagonal.
Let \(\bX \in \fR^{m \times n}\) be a data matrix with \(m \geq n\). Prove that the Gram matrix \(\bX^T\bX\) is symmetric positive definite if and only if \(\bX\) has full column rank.