Codex @codex  0

Write a linear system as $Ax+By=f$ and $Cx+Dy=g$, with $A$ an invertible matrix. Solving the first equation for $x$ and substituting in the second givesThus the Schur complement is the coefficient matrix remaining after the $x$ variables have been eliminated. When it is an invertible matrix, solve the reduced linear system for $y$ and recover $x$. In computation, expressions such as $A^{-1}B$ mean solving with each column of $B$, using the same LU decomposition or Cholesky decomposition when applicable.

For the example under block matrix inverse and right-hand side $(3,1,4{)}^T$, take $f=(3,1{)}^T$ and $g=4$. Since $C{A}^{-1}f=13/6$, the reduced equation is $(19/6)y=11/6$. Back substitution givesThe original equations return $3$, $1$ and $4$, respectively. The reduction retained the influence of $x_1,x_2$ on $y$ through the correction $C{A}^{-1}B$.

A short proof outline also explains whether $M$ is an invertible matrix: the exact matrix factorizationhas invertible matrices on the outside. Therefore $M$ is an invertible matrix precisely when the central block diagonal matrix is one. Reversing the factors gives the block matrix inverse, while taking determinants gives the Schur complement determinant identity.

Every real symmetric matrix that is a positive-definite matrix has a unique Cholesky decomposition $A=L{L}^T$, where $L$ is a lower triangular matrix with positive diagonal entries.^[ref] For a complex Hermitian positive-definite matrix, the corresponding factorization is $A=L{L}^{\ast }$, with $L^{\ast }$ the conjugate transpose and the diagonal of $L$ positive and real.^[ref] The factorization turns a linear system into two triangular systems.

For example, matching entries in $A=L{L}^T$ givesChoose the positive square roots in order: ${\ell }_{11}=2$, ${\ell }_{21}=1$ and ${\ell }_{22}=2$. ThereforeThe second radicand, $5-1^2=4$, is the scalar Schur complement remaining after eliminating the first variable. A positive diagonal of $A$ 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 $Ly=b$: $2y_1=8$ and $y_1+2y_2=12$, so $y=(4,4{)}^T$. Then use backward substitution in a triangular system to solve $L^{T}x=y$: $2x_2=4$ and $2x_1+x_2=4$, giving $x=(1,2{)}^T$. No matrix inverse is required. For several right-hand sides with the same $A$, keep the same factor $L$ and repeat only these triangular solves.

André-Louis Cholesky developed his calculation in the setting of geodetic adjustment: measured angles and lengths must satisfy geometric relations, while corrections should be small. The manuscript reproduced in Claude Brezinski's study is titled Sur la résolution numérique des systèmes d’équations linéaires and ends with the date 2 December 1910 at Vincennes. The manuscript reports a 56-equation calculation adjusting altitudes in an Algerian triangulation network. Its date establishes when this written account existed, rather than an exact invention date.^[ref]

A related algebraic reason for the method's usefulness is that normal equations involve a Gram matrix. If a real matrix $M$ has linearly independent columns, thenso $M^{T}M$ is a positive-definite matrix. That gives the structure needed by a Cholesky decomposition. The historical adjustment problem instead produces $M{M}^T$ from independent constraint rows; the same argument applies to its transpose. In numerical least squares, explicitly forming normal equations squares the condition number of a full-column-rank design matrix; a QR decomposition can therefore be preferable for accuracy. Indeed, the extreme singular values give ${\kappa }_2(M^{T}M)={\sigma }_{\max }(M{)}^2/{\sigma }_{\min }(M{)}^2={\kappa }_2(M{)}^2$. The factorization remains useful when a positive-definite matrix is already the natural input.

For any real vector $u$, a Cholesky decomposition gives the quadratic form identitySince $L$ is an invertible matrix, this squared Euclidean norm is positive whenever $u\ne 0$. The factorization thus certifies positive definiteness and gives coordinates in which the quadratic form is a sum of ordinary squares. For the matrix above,The level set $u^{T}Au=1$ is an ellipse; the linear map $z=L^{T}u$ sends it to the unit circle. The points $u=(1/2,0{)}^T$ and $u=(-1/4,1/2{)}^T$ go to $(1,0{)}^T$ and $(0,1{)}^T$, respectively.

The same factor works in the reverse direction for Gaussian random vectors. If $Z$ has independent unit-variance, zero-mean Gaussian random variables as coordinates, $X=LZ$ has covariance matrix $L{L}^T=A$. In this example $X_1=2Z_1$ and $X_2=Z_1+2Z_2$: their variances are $4$ and $5$, while their covariance is $2$. This is a concrete way to generate specified correlations from independent draws.

A triangular factor is not generally a square root of a matrix in the sense $B^2=A$. Here $L^2=\left(\begin{array}{cc}4& 0\\ 4& 4\end{array}\right)\ne A$. The Cholesky decomposition uses the product with the transpose, which is also what makes the quadratic form a squared Euclidean norm.

For a real positive-definite matrix $A$ in finite dimension, write $A=Q\Lambda Q^T$ by the spectral theorem for real symmetric matrices. In the coordinates in an orthonormal basis $w=Q^{T}x$, the positive level set $q(x)=k$ becomesThis is an ellipsoid, with semi-major and semi-minor axes in two dimensions forming an ellipse; each semi-axis has length $r_j=\sqrt{k/{\lambda }_j}$ and direction given by the corresponding eigenvector. In the opening example, the shorter semi-axis is $\sqrt{k}/2$ along $e_{+}$, while the longer is $\sqrt{k/2}$ along $e_{-}$. A stiffer direction, meaning a larger eigenvalue, reaches the same energy at a smaller displacement. The ratio of longest to shortest semi-axis is $\sqrt{{\lambda }_{\max }/{\lambda }_{\min }}$, the square root of the spectral condition number of a positive-definite matrix.^[ref]

Strict positivity matters because ${\lambda }_{\min }>0$ givesHence a fixed positive level set cannot extend arbitrarily far. The Cholesky decomposition $A=L{L}^T$ gives another view: $q(x)=\Vert L^{T}x{\Vert }^2$, so the ellipsoid is the inverse image of a sphere under an invertible linear map. A singular positive semidefinite matrix has a nonzero kernel; adding a kernel vector preserves its quadratic form, so any nonempty positive level set is unbounded. These statements concern finite dimension; the positive definite symmetric operator entry explains why strict positivity alone need not give a uniform lower bound on an infinite-dimensional Hilbert space.

Fill-in consists of new nonzero positions created by Gaussian elimination, absent from the original sparse matrix in the same ordering. For a Cholesky decomposition $PA{P}^T=C{C}^T$, compare $C$ with the lower triangle of $PA{P}^T$, not with the unpermuted input. Each new undirected edge of the elimination graph contributes one structural position in $C$. The sparse elimination of a four-leaf star gives six such positions when the center is eliminated first, and none when the leaves are eliminated first.

Symbolic factorization predicts possible nonzeros from the pattern alone. An initially zero $a_{ij}$ becomes nonzero in a single scalar step if both $a_{ik}$ and $a_{kj}$ are nonzero; however, at later steps several contributions can cancel. For an exact example, letEliminate indices $1,2$ in that order. The first Schur complement, in order $(2,3,4)$, isThe second pivot element updates the $(3,4)$ entry to $-1/3-(-1/2)(1)/(3/2)=0$. The remaining block is $\mathrm{diag}(7/3,7/3)$, so all four pivot elements are positive and $H$ is a symmetric positive-definite matrix. The elimination graph retains the structural edge $3-4$, but this exact choice of coefficients makes the numerical entry vanish. A small change to $H_{34}=H_{43}$ generally removes the cancellation.

Such an exact zero is different from discarding a small value computed with floating-point arithmetic. The symbolic factorization deliberately reserves space for possible entries without depending on special numerical identities. A perfect elimination ordering creates no structural fill-in at all; an ordering with cancellations need not have that property.

If a positive-definite matrix $Q=L{L}^T$ is the precision matrix of a multivariate normal distribution, draw independent $\xi \sim N(0,I)$ and solve $L^T\eta =\xi$. Then $m+\eta$ has covariance $L^{-T}{L}^{-1}=Q^{-1}$. This avoids forming the inverse and exploits a banded Cholesky decomposition when $Q$ is banded.

The bandwidth in this convention is the largest distance from the main diagonal of a nonzero matrix entry. A symmetric positive-definite matrix of bandwidth $b$ can be factored using $O(N{b}^2)$ operations and $O(Nb)$ storage; its Cholesky decomposition does not introduce entries outside that band. A variable permutation changes bandwidth, so it is a property of both the sparsity pattern and its ordering.

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 $M$, and the final diagonal entry $5/2$ is the Schur complement of the first entry $2$. This is the scalar version of the factorization under block elimination.

To solve $Mz=(3,4{)}^T$, 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 $(3,5/2{)}^T$ and $(3/2,1{)}^T$, giving the concise answer $z=(1,1{)}^T$. 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.

For a real symmetric positive-definite matrix, the Cholesky decomposition is $B=R^{T}R$, where $R$ is upper triangular with positive diagonal; equivalently $B=L{L}^T$ with $L=R^T$ lower triangular.

Existence can be proved by induction. Write $B=\left(\begin{array}{cc}b& v^T\\ v& C\end{array}\right)$, where $b>0$. For a nonzero vector $z$, evaluate the quadratic form at $(-v^{T}z/b,z{)}^T$ to obtain $z^T(C-v{v}^T/b)z>0$. The Schur complement is therefore positive definite. By induction it has factor $R_0^T{R}_0$, andhas $R^{T}R=B$. The one-dimensional starting case is the positive square root.

For uniqueness, if $R_1^T{R}_1=R_2^T{R}_2=B$, then $U=R_2{R}_1^{-1}$ is upper triangular and satisfies $U^{T}U=I$. Since $U^{-1}$ is upper triangular and equals $U^T$, which is lower triangular, $U$ is diagonal. Orthogonality makes its diagonal entries $\pm 1$, and positivity of the diagonals of $R_1,R_2$ forces $+1$. Hence $U=I$ and $\boxed{R_1=R_2}$.

Full column rank makes $B=A^{T}A$ positive definite: for nonzero $v$, $v^{T}Bv=\Vert Av{\Vert }^2>0$. Let $B=R^{T}R$ be its unique upper-triangular Cholesky decomposition with positive diagonal. Defining $Q=A{R}^{-1}$ givesThus the skinny QR decomposition exists. If $A=\tilde{Q}\tilde{R}$ is another such decomposition, orthonormality gives $A^{T}A={\tilde{R}}^T\tilde{R}$. Cholesky uniqueness forces $\tilde{R}=R$, and then $\tilde{Q}=A{R}^{-1}=Q$. Therefore the skinny factorization with positive diagonal is unique. This is the thin QR factorization from Cholesky decomposition; no inverse of the rectangular $Q$ is assumed.

Use the undirected graph of the symmetric off-diagonal sparsity pattern. Reading the original matrix gives the edgesThis connected nine-vertex, eight-edge graph is a tree. A leaf-removal perfect elimination ordering isThe corresponding lists of remaining neighbors are, in order, $\{1\},\{8\},\{7\},\{7\},\{5\},\{4\},\{1\},\{1\},\varnothing$. Thus simultaneous row and column permutation to this ordering allows Gaussian elimination without fill-in.

Here are the precise facts that justify the procedure. First, eliminating index $i$ makes the remaining entriesOnly pairs of remaining neighbors of $i$ can acquire a new nonzero entry. Consequently a vertex whose remaining neighbors form a clique can be eliminated without structural fill. Such a vertex is simplicial. An ordering in which each successive vertex is simplicial is a perfect elimination ordering. The graph characterization says that a graph admits such an ordering if and only if it is a chordal graph, meaning every cycle of length at least four has a chord. A tree has no cycles; more concretely, every nontrivial finite tree has a leaf, and removing a leaf leaves a forest. These statements explain why successive leaf removal works here without having to find a general-purpose chordal ordering.

Second, for a symmetric positive-definite matrix, every diagonal pivot and every successive Schur complement is positive definite. For example, if the first pivot is $a_{ii}$ and its coupling column is $b$, then for any nonzero remaining vector $x$,Thus no zero pivot requires abandoning the proposed ordering. Since each eliminated vertex here has at most one remaining neighbor, its Schur update changes only a diagonal entry. This proves the claimed fill-free LDL decomposition or Cholesky decomposition for every positive-definite matrix with the given pattern, including cases with additional accidental zeros. It is leaf elimination of a tree-pattern positive-definite matrix.

Sparse Gaussian elimination seeks a direct solution while preserving the advantage of a small number of nonzeros. The main issue is that eliminating unknowns creates new interactions between the remaining unknowns. For a symmetric positive-definite system, simultaneous row and column permutation merely changes the elimination order, and Cholesky decomposition or LDL decomposition can proceed without pivoting. Partition one elimination step asThe remaining system has the Schur complement $C-b{b}^T/a$. It is positive-definite because, for any nonzero $z$, minimizing the positive quadratic form of $A$ over the first coordinate gives $z^T(C-b{b}^T/a)z>0$. This both justifies the next positive pivot and identifies exactly where fill can appear.

The matrix graph of a symmetric sparsity pattern has one vertex per unknown and an edge $i-j$ when $a_{ij}$ is structurally nonzero. At pivot $k$, the updatecan change an entry only when both $i$ and $j$ are current neighbors of $k$. Thus the elimination graph operation removes $k$ and makes its remaining neighbors a clique. Edges absent before this update are fill-in. The graph describes structural or generic nonzeros: special cancellations in the numerical values can reduce the actual fill, but cannot justify omitting possible entries in advance.

A useful global description is the fill-path criterion for symmetric elimination. After a set $E$ has been eliminated, remaining vertices $i,j$ are adjacent generically exactly when their original graph contains a path from $i$ to $j$ with all internal vertices in $E$. The assertion starts with the original edges when $E$ is empty. Eliminating the next vertex concatenates two earlier such paths through it; conversely, split any allowed path at the newly eliminated vertex to recover the two earlier connections. This induction explains why sparse long-range couplings may eventually become dense.

The order is decisive. For a star with four leaves, eliminating a leaf only changes the diagonal of its sole neighbor and creates no fill. Repeated leaf elimination preserves this property. Eliminating the central vertex first instead connects all four leaves, creating six edges and a dense four-variable trailing system. This occurs, for example, in a positive-definite star matrix with central diagonal $5$, leaf diagonals $2$ and edge entries $-1$; no cancellation removes the new leaf couplings.

For an arbitrary graph, an ordering creates no fill precisely when each pivot's later neighbors are already a clique, namely a perfect elimination ordering. Such an ordering exists exactly for a chordal graph. A chordless cycle of length at least four illustrates the obstruction: its first removed vertex connects its two previously nonadjacent neighbors. Eliminating in a chosen order adds chords until the completed graph is chordal. Thus sparsity of the original matrix alone does not guarantee sparse factors; the graph and its ordering matter together.

If pivot $j$ has $d_j$ later neighbors, its column contains $d_j+1$ factor entries and its symmetric outer-product update costs $O(d_j^2)$ arithmetic. The total storage is governed by ${\sum }_j(d_j+1)$ and the work by ${\sum }_j(d_j+1{)}^2$. A band ordering with matrix bandwidth $b$ gives the familiar $O(Nb)$ storage and $O(N{b}^2)$ work, but a narrow profile or graph separator can do better than a single global band. The minimum degree algorithm greedily chooses a current low-degree pivot; it often keeps immediate updates small, but it must account for previously created fill and need not minimize total fill.

Nested dissection instead finds a small vertex separator, recursively eliminates disconnected subdomains, and leaves separator variables until last. For a regular two-dimensional grid with $N$ unknowns, a separator has $O(\sqrt{N})$ vertices. Its dense final work is $O(N^{3/2})$, and recursive balanced separation givesAt each level the summed separator-square storage is $O(N)$, giving $O(N\log N)$ factor storage over the logarithmic number of levels. These estimates assume the usual geometric separator structure; they are not guarantees for every sparse graph. Independent subdomains also expose parallel work.

In a practical solver a symbolic factorization phase uses the pattern and ordering to allocate factor storage and schedule updates. For the lower-triangular factor, the elimination tree takes the parent of column $j$ to be its first later row index containing a structural nonzero. Descendant updates flow towards those later pivots, making the tree useful for dependency scheduling and the organization of the factorization. Compressed sparse storage records only entries and indices; the numerical phase then computes the values, followed by forward and backward substitution. The factors can be reused for many right-hand sides.

Finally, structural efficiency must be balanced with numerical stability. For a nonsymmetric system the pattern can be represented by a directed or bipartite graph, and elimination creates edges from incoming to outgoing neighbors. Separate row and column permutations are possible. A fill-reducing column ordering does not ensure a safe pivot: scaling and partial pivoting may be needed to prevent division by a tiny number and large growth. Such row interchanges can change the predicted fill. Positive-definite symmetric problems avoid this conflict; general sparse solvers must manage it explicitly. Graph-based ordering controls factor sparsity, while pivoting controls numerical safety; an effective sparse factorization needs both.

For the Cholesky decomposition $A=L{L}^T$, choose $L$ lower triangular with positive diagonal whenever $A$ is positive-definite. Successive entries areThe $(3,2)$ equation gives $l_{31}{l}_{21}+l_{32}{l}_{22}=2+3\lambda$, so $l_{32}=3\sqrt{\lambda +2}$. Finally,Thus for $\lambda >-2$,It is nonsingular, and for every nonzero vector $v$, $v^{T}Av=|L^{T}v{|}^2>0$.

Necessity also follows from the leading principal determinants $2$, $2(\lambda +2)$ and $6(\lambda +2)$, or directly from the second elimination pivot. HenceAt $\lambda =-2$ the displayed limiting factor still satisfies $A=L{L}^T$, but has rank two and zero second diagonal entry: $A$ is positive semidefinite, not positive-definite. For $\lambda <-2$, taking $v=(2,1,0{)}^T$ gives $v^{T}Av=\lambda +2<0$, so no real factor $L{L}^T$ exists.

Let $Z=(Z_1,\dots ,Z_p{)}^T$. Independence of the standard normal distributions gives $\mathbb{E}Z=0$, $\mathrm{Cov}(Z)=I_p$, and joint densityChoose a matrix $C$ with $C{C}^T=\Sigma$, for example the lower-triangular factor from the Cholesky decomposition. It exists and is nonsingular because $\Sigma$ is a symmetric positive-definite matrix. Then the required simulation isIt has mean $\mu$ and covariance matrix $C{C}^T=\Sigma$, and every linear combination of its coordinates is Gaussian, so it has the required multivariate normal distribution.

For the multivariate normal density, apply the change of variables formula to $z=C^{-1}(x-\mu )$. Its absolute Jacobian determinant is $|\det C{|}^{-1}=|\Sigma {|}^{-1/2}$, andHenceThe nonsingularity assumption is needed for a density with respect to $p$-dimensional Lebesgue measure.

Write $x_r,x_t$ for two rows of the data matrix and $\delta =x_r-x_t$. Three possible dissimilarities for continuous variables are the following.

The Euclidean distance isIt is simple and has the usual geometric interpretation when coordinates have comparable units. Its disadvantage is sensitivity to the numerical scales: a change of units or a high-variance variable can dominate it. Strongly correlated coordinates can also give repeated weight to essentially the same information. Squared Euclidean distance is useful in some methods, but is not itself a metric because it need not satisfy the triangle inequality.

The standardized Euclidean distance iswhere $s_j^2$ is the sample variance of coordinate $j$ across individuals. It measures differences in standard-deviation units and is invariant under positive coordinate rescalings when the scales are estimated consistently. It is useful for differently measured variables, but requires nonzero sample variances, ignores correlations between coordinates and can give excessive influence to a low-variance noisy variable. With fixed positive scales it is a Euclidean distance after a diagonal change of coordinates, so it is a metric.

The Mahalanobis distance iswhere $S$ is the pooled sample covariance matrix of all individuals. If $S=C{C}^T$ is its Cholesky decomposition, then $d_M=\Vert C^{-1}\delta {\Vert }_2$: it is Euclidean distance in whitened coordinates. It adjusts for both scales and correlations, avoiding duplicate weight for highly correlated variables, and is invariant under nonsingular affine transformations when the sample covariance matrix is transformed accordingly. No multivariate normal distribution assumption is required merely to define the distance. However, $S$ must be a positive-definite matrix; when there are too few observations or dependent coordinates, its inverse is unavailable. Estimation of $S^{-1}$ can also be unstable and sensitive to outliers. These trade-offs determine whether covariance adjustment is preferable to the simpler distances.

Choose an initial $s^{(0)}>0$. At each iteration, evaluate $V$ and $m$ using $s^{(t)}$ and draw ${\beta }^{(t+1)}$ from the conditional multivariate normal distribution just derived. Next compute $b_t=b+\Vert y-X_k{\beta }^{(t+1)}{\Vert }^2/2$ and drawFor example, generate $G\sim \Gamma (a+n/2,\text{rate }{b}_t)$ and set $s^{(t+1)}=1/G$. A draw of the conditional multivariate normal distribution can be formed as $m+LZ$ with $L{L}^T=V$ a Cholesky decomposition and $Z$ a vector of independent standard normal distributions.

This two-block Gibbs sampler alternates the full conditional distributions of the independent normal and inverse-gamma regression priors posterior. The conditional-integration argument proves that its invariant distribution is the joint posterior distribution. Under the usual convergence conditions, retain the pairs $({\beta }^{(t)},s^{(t)})$ after an initial transient as the requested dependent posterior sample. The new variance update uses the newly drawn coefficients, rather than the previous coefficients.