Cholesky Decomposition

Every so often in mathematics, a structure hides a shortcut. Cholesky decomposition is one of the cleanest examples: if you know a matrix is symmetric ($A = A^T$) and positive-definite (all eigenvalues positive, or equivalently $x^T A x > 0$ for all $x \neq 0$), you can factor it as

$$A = L L^T$$

where $L$ is a lower-triangular matrix. Because $L$ and $L^T$ carry the same entries — just mirrored — you only need to compute half as many numbers as you would for a general LU decomposition.

The method is named after André-Louis Cholesky, a French military officer and geodesist who devised it around 1910 while adjusting survey measurements for the map of Crete. He never published it himself; the algorithm reached the world through a posthumous note by his colleague Commandant Benoit in 1924.

Today Cholesky decomposition is one of the workhorses of numerical linear algebra: it solves systems $Ax = b$ efficiently, powers Kalman filters, Gaussian process regression, and the "square root" step inside countless optimization algorithms.

Below is a $3 \times 3$ symmetric positive-definite matrix $A$. The demo computes the Cholesky factor $L$ one column at a time, highlighting each entry as it is calculated using the recurrence formulas. You can edit the diagonal entries to try different matrices, then verify that $L L^T = A$.

<!-- {{c_html_intro}} -->
<div class="container">
  <div class="section">
    <h3>{{label_matrix_a}}</h3>
    <p class="hint">{{hint_edit}}</p>
    <div class="matrix-wrap">
      <table id="matA" class="matrix"></table>
    </div>
  </div>
  <div class="section">
    <h3>{{label_factor_l}}</h3>
    <p class="hint">{{hint_l}}</p>
    <div class="matrix-wrap">
      <table id="matL" class="matrix"></table>
    </div>
  </div>
</div>
<div class="btns">
  <button id="btnStep" type="button">{{btn_step}}</button>
  <button id="btnAll" type="button">{{btn_all}}</button>
  <button id="btnReset" type="button" class="ghost">{{btn_reset}}</button>
</div>
<div id="status" class="status"></div>
<div id="formula" class="formula"></div>

/* {{c_css_intro}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: .4rem; }
h3 { margin: 0 0 .3rem; font-size: 1rem; color: #1d3557; }
.container { display: flex; gap: 1.5rem; flex-wrap: wrap; margin-bottom: .6rem; }
.section { flex: 1 1 140px; }
.hint { font-size: .82rem; color: #555; margin: 0 0 .4rem; line-height: 1.4; }
.matrix-wrap { display: inline-block; border-left: 2px solid #1d3557; border-right: 2px solid #1d3557;
               padding: 4px 8px; }
.matrix { border-collapse: separate; border-spacing: 4px; }
.matrix td { width: 44px; height: 36px; text-align: center; font: 600 14px ui-monospace, monospace;
             background: #e8eef3; border-radius: 6px; border: 1px solid #cdd9e3; color: #1d3557;
             transition: background .25s, color .25s; }
.matrix td.editable { background: #f0f4f8; cursor: text; }
.matrix td.zero { color: #aaa; }
.matrix td.active { background: #ffd166; border-color: #e0a800; color: #333; }
.matrix td.done { background: #c8f0d4; border-color: #5cb87a; }
.matrix td.empty { background: #f7f7f7; color: #bbb; border-color: #ddd; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 14px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button:disabled { opacity: .45; cursor: default; }
button.ghost { background: #fff; color: #1d3557; }
.status { font-size: .95rem; font-weight: 600; min-height: 1.3em; margin-bottom: .3rem; }
.status.ok { color: #0a7d33; }
.status.bad { color: #c92f3c; }
.status.info { color: #1d3557; }
.formula { font-size: .83rem; color: #555; min-height: 1.2em; font-family: ui-monospace, monospace; }

// Code not found

Notice the pattern: the diagonal entry $l_{jj}$ is a square root (it must be real because $A$ is positive-definite), and each off-diagonal entry $l_{ij}$ is just a dot product divided by $l_{jj}$. This column-by-column sweep visits each entry of $L$ exactly once — roughly half the work of a full LU decomposition on the same matrix.

How fast is Cholesky, and why does symmetry help so much?

• Operation count. A general LU decomposition of an $n \times n$ matrix costs about $\frac{2}{3}n^3$ floating-point operations. Cholesky costs roughly $\frac{1}{3}n^3$ — exactly half — because it only touches the lower triangle.
• Numerical stability. For symmetric positive-definite matrices, Cholesky is unconditionally stable: no pivoting is required. The diagonal entries $l_{jj} = \sqrt{a_{jj} - \sum_{k=1}^{j-1} l_{jk}^2}$ are always real and positive, so the algorithm never divides by zero or amplifies rounding errors the way partial pivoting must guard against in general LU.
• Existence and uniqueness. A matrix has a Cholesky factorization if and only if it is symmetric positive-definite — and when it exists, the factorization with positive diagonal is unique. This equivalence makes "does Cholesky succeed?" a practical test for positive-definiteness.
• Blocked and parallel variants. Modern LAPACK implementations use blocked Cholesky (DPOTRF), processing panels of columns at once to exploit cache locality and BLAS-3 matrix-matrix operations. On multi-core and GPU hardware, the triangular solves that follow (DPOTRS) are also highly parallel.

Cholesky is a solved, polished algorithm — not an open problem. Its complexity class is the same as matrix multiplication in general, and the $O(n^3)$ bound has been the practical standard since the 1960s. Current research focuses on sparse Cholesky (exploiting sparsity patterns in finite-element meshes) and randomized Cholesky for approximate factorizations of very large matrices.

Symmetric positive-definite matrices appear wherever you measure covariance, minimize a quadratic objective, or discretize a self-adjoint differential operator — and Cholesky is almost always the right tool:

• Statistics and machine learning. Covariance matrices are symmetric positive-semidefinite by construction. Cholesky gives you the "square root" of a covariance matrix: if $\Sigma = L L^T$ and $z \sim \mathcal{N}(0, I)$, then $Lz \sim \mathcal{N}(0, \Sigma)$. Gaussian process regression relies on this at every prediction step.
• Kalman filtering. The "square-root Kalman filter" propagates $L$ (the Cholesky factor of the covariance matrix) rather than the matrix itself, giving better numerical conditioning for navigation, tracking, and control systems.
• Convex optimization. Interior-point methods for semidefinite programming and quadratic programs solve systems $Ax = b$ with symmetric positive-definite $A$ at every iteration — Cholesky is the inner loop.
• Finite-element analysis. Discretizing an elliptic PDE produces a large, sparse, symmetric positive-definite stiffness matrix. Sparse Cholesky (with fill-reducing reordering like AMD or nested dissection) is the standard direct solver.
• Monte Carlo simulation. Generating correlated random variables from a multivariate distribution requires factoring the correlation matrix — a textbook Cholesky application.

In short: wherever symmetry and positive-definiteness arise together — which is most of applied mathematics — Cholesky delivers the fastest, numerically cleanest linear solve available.

Cholesky decomposition is a beautiful case of structure meeting efficiency: when you know a matrix is symmetric and positive-definite, you don't have to pretend it isn't. Exploit that structure and the work literally halves, pivoting disappears, and numerical stability comes for free.

André-Louis Cholesky discovered this while fitting survey measurements together on a Cretan hillside. A century later, his algorithm lives inside every statistical package, every Kalman filter, every interior-point optimizer — quietly doing its job in $\frac{1}{3}n^3$ flops with no fuss.

The lesson generalizes: the best algorithms are not always the most general ones. A method that knows something extra about its input — symmetry, sparsity, positivity — can be twice as fast, twice as stable, and far more elegant. Understanding when those assumptions hold is as important as understanding the algorithm itself.