Quadratic Forms and Positive Definiteness

A real symmetric matrix can be diagonalized with an orthonormal basis: A = QΛQᵀ, and in the rotated coordinates z = Qᵀx, the quadratic form becomes Σᵢ λᵢ zᵢ². This expression separates the contribution of each principal direction. Large positive eigenvalues penalize movement strongly, small values create shallow directions, and zeros create flat directions. The decomposition is not merely a proof device. It is the geometry behind condition numbers, principal components, Gaussian covariance, and second-order optimization.

Those signs are not an artifact of the chosen basis. The result that says so is shorter and more offhand than its standing suggests. James Joseph Sylvester published it in 1852, in a five-page note in the Philosophical Magazine. The business of that note is a proof that the roots of the characteristic equation of a symmetric matrix are real. The invariance is tossed off as a remark appended to it — “the number of positive and negative coefficients is invariable”. A line in the margin of someone else's proof became a law.

The modern statement is the same one. The Encyclopedia of Mathematics fixes the counts under any non-singular real change of variables — “The number p (respectively, n) of indices i for which b_i > 0 (or b_i < 0) is fixed” — records that “In its classical form, the law of inertia was established by J.J. Sylvester”, and names the invariant: “The pair (p, n) is called the signature of f”. Positive definiteness belongs to the form, not to the coordinates. It has belonged to it since 1852.

The quadratic form xᵀAx takes a vector, applies a matrix, and measures the alignment of the result with the original vector. Its value changes with both direction and magnitude. When A is symmetric positive semidefinite, the value is never negative. It can represent energy, curvature, variance, or a squared distance after a suitable transformation. That is why covariance ellipses, least-squares objectives, Gaussian exponents, and L2 regularization share the same geometry.

The standard test for positive definiteness carries a soldier's name, and it reached print only because a fellow officer carried it there. André-Louis Cholesky devised the factorization for adjusting geodetic grids by least squares. He “died from wounds received on the battle field on 31 August 1918 at 5 o'clock in the morning in the North of France”, as MacTutor's biography of him records. Commandant Benoit published the method in 1924, in the Bulletin géodésique, under a title whose last words are “Procédé du Commandant Cholesky”.

The death is recorded in the note itself, in a footnote: “Sur le Commandant Cholesky, tué à l’ennemi le 31 août 1918, voir la notice biographique insérée dans le volume du Bulletin géodésique de 1922”. The routine that every numerical library now reaches for first is a posthumous eleven-page item in a geodesy bulletin. Its author's fate is in the margin.

A symmetric matrix has orthogonal eigenvectors that define principal directions, and eigenvalues determine curvature along them. The four cases are not equally common once the matrix comes from data rather than from a textbook. Sagun and colleagues measured the Hessian of a trained deep network in 2016 and found neither positive definiteness nor full effective rank: “The eigenvalue distribution is seen to be composed of two parts, the bulk which is concentrated around zero, and the edges which are scattered away from zero.” Almost all the directions are flat. A handful carry the curvature.

The same measurement reached ImageNet scale in 2019. Ghorbani and colleagues computed the full Hessian spectrum rather than a few extreme eigenvalues, and found that the shape depends on the architecture: “In non-batch normalized networks, we observe the rapid appearance of large isolated eigenvalues in the spectrum ... In a batch normalized network, these two effects are almost absent.” The classification below is a vocabulary for reading such a spectrum. It is not a set of boxes real matrices fall into neatly.

Every covariance, Hessian, metric, or regularizer that appears can be read the same way. In production the reading is not done with an eigendecomposition. The test for positive definiteness is an attempted Cholesky factorization. LAPACK's reference documentation for dpotrf says that “DPOTRF computes the Cholesky factorization of a real symmetric positive definite matrix A”, and it localises the failure rather than merely reporting one: on exit, “> 0: if INFO = i, the leading principal minor of order i is not positive, and the factorization could not be completed”. Independent vendor implementations of the same interface repeat the convention. Intel oneMKL's ?potrf documents info > 0 as “The leading minor of order info is not positive definite, so the factorization could not be completed”. NAG documents the same routine as f07fdf, callable as dpotrf, and draws the conclusion out loud: “Hence A itself is not positive definite. This may indicate an error in forming the matrix A.” One factorization, and the error code is the order of the offending leading principal minor.

When the matrix is empirical, that failure is routine enough to have its own literature. Nicholas J. Higham posed the repair in 2002 as finding “the nearest symmetric positive semidefinite matrix with unit diagonal”, and reported where the damage comes from: “In the finance application the original matrix has many zero or negative eigenvalues; we show that for a certain class of weights the nearest correlation matrix has correspondingly many zero eigenvalues and that this fact can be exploited in the computation.” NAG ships that problem as routine g02aaf, which “computes the nearest correlation matrix, in the Frobenius norm, to a given square, input matrix”. Its definition of the target is plain: “A correlation matrix may be characterised as a real square matrix that is symmetric, has a unit diagonal and is positive semidefinite”. The algorithm inside g02aaf is not Higham's 2002 modified alternating projections. It is the quadratically convergent Newton method of Qi and Sun (2006), with the preconditioning of Borsdorf and Higham (2010) — the two works its documentation cites. The problem is Higham's. The solver has moved on.

A nonsingular covariance matrix Σ has precision matrix Σ⁻¹, and the zeros in that inverse are not a modelling convenience. Meinshausen and Buhlmann opened their 2006 paper in The Annals of Statistics with the fact itself: “The pattern of zero entries in the inverse covariance matrix of a multivariate normal distribution corresponds to conditional independence restrictions between variables.” A zero in the precision matrix is a conditional independence, not an approximation to one.

Finding those zeros is a solved computation. The graphical lasso, published in Biostatistics in 2008, states the move plainly: “We consider the problem of estimating sparse graphs by a lasso penalty applied to the inverse covariance matrix.” Then it reports the scale. A 1000-node problem, roughly 500,000 parameters, solved in at most a minute — 30 to 4000 times faster than the competing methods of the time.

Suppose height and arm span are strongly correlated in a population. The distance that respects that correlation is a quadratic form written with an inverse. NIST/SEMATECH's e-Handbook of Statistical Methods gives the statistic for p variables, in its section on “Hotelling's T squared”, as T^2 = n (xbar - mu_0)' S^{-1} (xbar - mu_0), and identifies S^{-1} as the inverse of the sample variance-covariance matrix.

Careful libraries never build that inverse. NAG's routine g03daf forms an upper triangular R_j by QR decomposition, so that S_j = R_j^T R_j. Then, per its description, “The distances are computed as z^T z, where R_j z = (x - xbar_j)” — a triangular solve where the formula shows a matrix inverse. That is also why a near-singular covariance degrades gradually instead of exploding. It is a quality-of-implementation choice rather than a universal practice: SciPy's scipy.spatial.distance.mahalanobis takes VI, “the inverse of the covariance matrix”, as an argument and evaluates the form from it directly.