Agricultural Statistics & Data Science · Foundations · also Matrix Algebra, Matrix Theory, Mathematics for Statistics (part)
Vector spaces, rank, inverses, eigenvalues, quadratic forms and generalised inverses: the language of linear models, multivariate analysis and mixed models. Taught here with the statistician's applications (projection matrices, idempotence, the distribution of quadratic forms).
In the exams: GATE ST 'Matrix Theory' (heavy); GATE DA linear algebra (heavy); IIT JAM MS (heavy); AIEEA-PG Group 07 Unit II (present); AICE-JRF Ag Stat Unit 4 linear models (assumed)
Linear algebra is the grammar of modern statistics. A field trial with treatments and blocks, a regression of yield on rainfall and fertiliser, a mixed model for clones across sites, a principal component analysis of crop traits and a survey estimator with weights are all written and solved as matrix equations. The least-squares estimator is (X′X)⁻¹X′y; the fitted values come from a projection matrix; sums of squares in analysis of variance are quadratic forms whose distributions depend on whether their matrices are idempotent; principal components are eigenvectors; and when a design matrix is not of full rank, as in every one-way classification model written with an overall mean, generalised inverses take the place of ordinary inverses. A statistician who is fluent in matrices reads linear models, multivariate analysis and mixed models as natural extensions of one idea; one who is not finds them a collection of separate formulas.
This manual covers vector spaces, subspaces, linear independence, basis and dimension; rank, nullity and row-reduced echelon form; determinants, inverses and partitioned matrices; systems of linear equations and their consistency; inner products, orthogonality and Gram–Schmidt orthonormalisation; eigenvalues, eigenvectors and the Cayley–Hamilton theorem; diagonalisation and the spectral decomposition of symmetric matrices; quadratic forms and definiteness; idempotent and projection matrices; generalised inverses, including the Moore–Penrose inverse; the singular value decomposition; and Kronecker products and the vec operator. Each topic is linked to its statistical use.
Graduate entrance examinations in statistics, mathematical statistics and data science test matrix theory heavily, with short numerical questions on rank, determinants, eigenvalues, inverses, definiteness, traces of projection matrices and dimensions of solution spaces, and theoretical questions on idempotence, generalised inverses and quadratic forms. National eligibility tests in the mathematical sciences ask proof-style true-or-false items, and master's admission examinations in agricultural statistics include basic matrix algebra and determinants.
A vector space over the real numbers is a set of objects, called vectors, that can be added and multiplied by scalars, with the usual rules: addition is commutative and associative, there is a zero vector, every vector has a negative, and scalar multiplication distributes over addition. The standard example is Rⁿ, the set of n-tuples of real numbers, but the space of all m × n matrices, the space of polynomials of degree at most k, and the space of functions on an interval are also vector spaces. A subspace is a non-empty subset that is itself a vector space under the same operations; to check, it suffices that it contains the zero vector and is closed under addition and scalar multiplication. In R³, the subspaces are {0}, lines through the origin, planes through the origin and R³ itself; a plane that does not pass through the origin is not a subspace. The intersection of two subspaces is a subspace; their union usually is not. The column space of a matrix A, the set of all linear combinations Ax of its columns, and the null space, the set of solutions of Ax = 0, are the two subspaces that matter most in statistics: in regression, fitted values lie in the column space of the design matrix X, and the null space describes parameter combinations that the data cannot distinguish.
Vectors v₁, …, vₖ are linearly independent if the only solution of c₁v₁ + … + cₖvₖ = 0 is c₁ = … = cₖ = 0; otherwise one of them is a combination of the others. In Rⁿ, more than n vectors are always dependent. The vectors (1, 2, 3) and (2, 4, 6) are dependent because the second is twice the first. A set spans a space if every vector in the space is a linear combination of its members. A basis is a linearly independent spanning set, and every basis of a given space has the same number of vectors, called the dimension. Rⁿ has dimension n; the space of m × n matrices has dimension mn; the space of n × n symmetric matrices has dimension n(n + 1)/2, which is 6 for n = 3; and the space of n × n skew-symmetric matrices has dimension n(n − 1)/2, which is 3 for n = 3. If U and W are subspaces, dim(U + W) = dim U + dim W − dim(U ∩ W). Each vector has unique coordinates with respect to a given basis, and changing the basis changes the coordinates through an invertible matrix.
Elementary row operations, namely swapping two rows, multiplying a row by a non-zero scalar and adding a multiple of one row to another, do not change the solution set of a linear system or the row space of a matrix. Each corresponds to multiplying on the left by an invertible elementary matrix. Gaussian elimination uses them to reduce a matrix to row echelon form, in which each non-zero row starts with a leading entry (pivot) to the right of the pivot above it and zero rows are at the bottom; reduced row echelon form further makes each pivot 1 and the only non-zero entry in its column. The reduced form of a matrix is unique. For A with rows (1, 2, 3), (2, 4, 6) and (1, 1, 1), subtracting twice the first row from the second gives a zero row, and subtracting the first from the third gives (0, −1, −2), so there are two pivots.
The rank of a matrix is the number of pivots in its echelon form, which equals the dimension of its column space and of its row space: row rank equals column rank. The example above has rank 2. Rank properties used in statistics include rank(A) ≤ min(m, n) for an m × n matrix; rank(AB) ≤ min(rank A, rank B); rank(A′A) = rank(AA′) = rank(A); multiplying by an invertible matrix does not change rank; and rank(A + B) ≤ rank A + rank B. A square matrix is invertible if and only if it has full rank. The rank–nullity theorem states that for an m × n matrix, rank(A) + nullity(A) = n, where the nullity is the dimension of the null space; a 4 × 6 matrix of rank 3 therefore has a three-dimensional null space.
In linear models, the rank of the design matrix X determines how many parameters, or functions of parameters, can be estimated. The one-way classification model yᵢⱼ = μ + τᵢ + εᵢⱼ with t treatments has t + 1 parameters but its design matrix has rank t, because the column for μ is the sum of the treatment columns. Individual parameters are therefore not estimable, but contrasts such as τ₁ − τ₂ are; this rank deficiency is why generalised inverses are needed.
The determinant of a square matrix is a scalar that measures how the linear transformation scales volume, with sign indicating orientation. For a 2 × 2 matrix with rows (a, b) and (c, d), det = ad − bc; larger determinants are computed by cofactor expansion along a row or column or, more efficiently, by reducing to triangular form, whose determinant is the product of diagonal entries. Key properties: det(AB) = det A det B; det(A′) = det A; det(A⁻¹) = 1/det A; swapping two rows changes the sign; adding a multiple of one row to another leaves it unchanged; multiplying one row by k multiplies the determinant by k, so for an n × n matrix det(kA) = kⁿ det A (if det A = 2 for a 3 × 3 matrix, det 2A = 16); a matrix with two equal rows or a zero row has determinant 0; and det A equals the product of the eigenvalues. A square matrix is invertible if and only if its determinant is non-zero. In statistics, the determinant of a covariance matrix, the generalised variance, appears in the multivariate normal density and in likelihood ratio tests, and a determinant near zero signals near-collinearity.
The inverse A⁻¹ of a square matrix satisfies AA⁻¹ = A⁻¹A = I and exists only if det A ≠ 0. For a 2 × 2 matrix, A⁻¹ = (1/(ad − bc)) × the matrix with rows (d, −b) and (−c, a); the inverse of the matrix with rows (2, 1) and (1, 1) is the matrix with rows (1, −1) and (−1, 2). In general A⁻¹ = adj(A)/det A, where the adjugate is the transpose of the matrix of cofactors, though in practice inverses are computed by elimination. Rules: (AB)⁻¹ = B⁻¹A⁻¹; (A′)⁻¹ = (A⁻¹)′; the inverse of a symmetric matrix is symmetric; the inverse of a diagonal matrix inverts each diagonal entry; and the inverse of an orthogonal matrix is its transpose. The Sherman–Morrison–Woodbury formula gives the inverse of a matrix plus a low-rank update and underlies efficient updating of regression results when one observation is added or deleted.
A matrix can be partitioned into blocks, and blocks multiply like entries provided their dimensions conform. For a partitioned matrix with blocks A (top left), B, C and D, if A is invertible, det = det A × det(D − CA⁻¹B), where D − CA⁻¹B is the Schur complement of A. The inverse of a partitioned matrix can be written in terms of the inverse of the Schur complement. These formulas are central in statistics: the conditional distribution of one part of a multivariate normal vector given another has covariance equal to a Schur complement, Σ₂₂ − Σ₂₁Σ₁₁⁻¹Σ₁₂; partial correlations come from blocks of the inverse covariance matrix; and the variance of a subset of regression coefficients is the corresponding block of σ²(X′X)⁻¹.
A system Ax = b of m equations in n unknowns is consistent, having at least one solution, if and only if b lies in the column space of A, which is equivalent to rank(A) = rank([A | b]), the rank of the augmented matrix (the Rouché–Capelli or Kronecker–Capelli theorem). If the ranks are equal to n, the solution is unique; if they are equal but less than n, there are infinitely many solutions, forming an affine space of dimension n − rank(A), the general solution being a particular solution plus any vector in the null space of A. If rank([A | b]) > rank(A), the system is inconsistent. A homogeneous system Ax = 0 is always consistent and has a non-trivial solution if and only if rank(A) < n; for a square matrix, if and only if det A = 0. For three equations in three unknowns with a coefficient matrix of rank 2 and an augmented matrix of rank 2, the solutions form a line, one free parameter; if the augmented rank is 3, there is no solution.
Gaussian elimination with back substitution solves systems efficiently; Cramer's rule, which expresses each unknown as a ratio of determinants, is useful for small systems and theory but inefficient in practice. Numerical stability requires pivoting, choosing large pivots to limit rounding error, and ill-conditioned systems, whose coefficient matrices are nearly singular, amplify errors in the data.
In regression, the system Xβ = y is usually inconsistent, because y does not lie in the column space of X. Least squares instead solves the normal equations X′Xβ = X′y, which are always consistent, since X′y lies in the column space of X′, which equals the column space of X′X. When X has full column rank, X′X is invertible and β̂ = (X′X)⁻¹X′y is unique. When it does not, the normal equations have infinitely many solutions β̂ = (X′X)⁻X′y for any generalised inverse, but the fitted values Xβ̂ and estimable functions are the same for every solution.
The standard inner product of two vectors in Rⁿ is u′v = Σuᵢvᵢ, and the norm (length) is ‖u‖ = √(u′u). The angle θ between non-zero vectors satisfies cos θ = u′v/(‖u‖‖v‖), and the Cauchy–Schwarz inequality |u′v| ≤ ‖u‖‖v‖ guarantees that this cosine lies between −1 and 1; applied to centred data vectors it shows that a correlation coefficient lies between −1 and 1, and indeed the sample correlation is the cosine of the angle between the centred vectors. Vectors are orthogonal if u′v = 0. A set of mutually orthogonal unit vectors is orthonormal. More general inner products ⟨u, v⟩ = u′Wv with positive definite W define weighted distances, such as the Mahalanobis distance (x − μ)′Σ⁻¹(x − μ), which accounts for the variances and correlations of variables.
A square matrix Q is orthogonal if Q′Q = I, so Q⁻¹ = Q′; its columns form an orthonormal basis, it preserves lengths and angles, its determinant is ±1, and its eigenvalues have absolute value 1. Rotations and reflections are orthogonal transformations. Helmert matrices, orthogonal matrices whose first row is proportional to (1, 1, …, 1), are used to prove that the sample mean and sample variance of normal data are independent. The orthogonal complement of a subspace consists of all vectors orthogonal to every vector in it; Rⁿ is the direct sum of a subspace and its complement, so every vector splits uniquely into a component in the subspace and one orthogonal to it. In regression, y splits into fitted values in the column space of X and residuals orthogonal to it, which is why residuals sum to zero when the model has an intercept and why the total sum of squares splits into regression and residual parts.
The Gram–Schmidt process turns linearly independent vectors u₁, …, uₖ into an orthonormal set spanning the same space. Set e₁ = u₁/‖u₁‖; then for each later vector subtract its projections onto the earlier orthonormal vectors and normalise: v₂ = u₂ − (u₂′e₁)e₁, e₂ = v₂/‖v₂‖, and so on. For u₁ = (1, 1, 0) and u₂ = (1, 0, 1): e₁ = (1, 1, 0)/√2; u₂′e₁ = 1/√2, so v₂ = (1, 0, 1) − (1/2)(1, 1, 0) = (1/2, −1/2, 1), with norm √(3/2), and e₂ = v₂/√(3/2). Applied to the columns of a matrix, Gram–Schmidt gives the QR decomposition A = QR, with Q having orthonormal columns and R upper triangular. Least-squares software uses QR rather than inverting X′X because it is more accurate numerically; orthogonal polynomials for trend contrasts in designed experiments are also built by Gram–Schmidt.
A scalar λ is an eigenvalue of a square matrix A, with eigenvector x ≠ 0, if Ax = λx: the transformation stretches x by λ without changing its direction. Eigenvalues are the roots of the characteristic equation det(A − λI) = 0, a polynomial of degree n. For the matrix with rows (2, 1) and (1, 2), det(A − λI) = (2 − λ)² − 1 = 0 gives λ = 1 and 3, with eigenvectors (1, −1) and (1, 1). Basic facts: the sum of the eigenvalues equals the trace and their product equals the determinant; a triangular matrix has its diagonal entries as eigenvalues; A is singular if and only if 0 is an eigenvalue; if A has eigenvalue λ, then Aᵏ has λᵏ, A⁻¹ has 1/λ, A + cI has λ + c, and cA has cλ, with the same eigenvectors; A and A′ have the same eigenvalues; and AB and BA have the same non-zero eigenvalues. Eigenvectors belonging to distinct eigenvalues are linearly independent. If a 3 × 3 matrix has eigenvalues 1, 2 and 3, then its trace and determinant are both 6, the eigenvalues of A² are 1, 4 and 9, det(A⁻¹) = 1/6, and A − 2I is singular.
The Cayley–Hamilton theorem states that every square matrix satisfies its own characteristic equation. For a 2 × 2 matrix the characteristic polynomial is λ² − (tr A)λ + det A, so A² − (tr A)A + (det A)I = 0. For the matrix with rows (1, 2) and (3, 4), tr A = 5 and det A = −2, so A² = 5A + 2I, and multiplying by A⁻¹ gives A⁻¹ = (A − 5I)/2, the matrix with rows (−2, 1) and (1.5, −0.5). The theorem expresses any power of an n × n matrix, and its inverse, as a polynomial in A of degree at most n − 1, and is used in the theory of Markov chains and in computing matrix functions. The minimal polynomial, the monic polynomial of least degree that A satisfies, divides the characteristic polynomial and has the same roots.
Eigenvalues appear wherever variance is decomposed. In principal component analysis, the eigenvalues of a covariance or correlation matrix are the variances of the components and the eigenvectors their loadings. In canonical correlation and discriminant analysis, eigenvalues of products of covariance matrices measure association. The stationary distribution of a Markov chain is the left eigenvector of its transition matrix for eigenvalue 1. In population models such as the Leslie matrix, the dominant eigenvalue is the long-term growth rate. The ratio of the largest to the smallest eigenvalue of X′X, the condition number, diagnoses multicollinearity in regression.
A square matrix A is diagonalisable if A = PDP⁻¹ for an invertible P and diagonal D; the columns of P are eigenvectors and the diagonal of D the eigenvalues. A is diagonalisable if and only if it has n linearly independent eigenvectors, which is guaranteed when the n eigenvalues are distinct. Some matrices are not diagonalisable: the matrix with rows (1, 1) and (0, 1) has the repeated eigenvalue 1 but only one independent eigenvector. Diagonalisation makes powers easy, Aᵏ = PDᵏP⁻¹, which is how long-run behaviour of Markov chains and linear dynamic systems is studied.
Real symmetric matrices, which include every covariance matrix, every X′X and every projection matrix in statistics, have special properties: all eigenvalues are real; eigenvectors of distinct eigenvalues are orthogonal; and the matrix is always diagonalisable by an orthogonal matrix. The spectral decomposition is A = QΛQ′ = Σλᵢqᵢqᵢ′, where Q is orthogonal with the eigenvectors qᵢ as columns and Λ is diagonal with the eigenvalues. For the matrix with rows (2, 1) and (1, 2), q₁ = (1, 1)/√2 with λ = 3 and q₂ = (1, −1)/√2 with λ = 1, and A = 3q₁q₁′ + 1q₂q₂′. The decomposition gives matrix functions directly: A⁻¹ = QΛ⁻¹Q′, Aᵏ = QΛᵏQ′, and, for a non-negative definite A, a symmetric square root A^(1/2) = QΛ^(1/2)Q′. Symmetric square roots and the related Cholesky decomposition are used to generate correlated random variables, to standardise multivariate data and to transform a generalised least-squares problem with covariance Σ into an ordinary one.
If Σ is the covariance matrix of a random vector x, the spectral decomposition Σ = QΛQ′ means that the transformed variables y = Q′x are uncorrelated with variances λ₁ ≥ λ₂ ≥ …; these are the principal components. The first component direction q₁ maximises the variance a′Σa over unit vectors a, a result that follows from the spectral decomposition or from Lagrange multipliers, and the maximum is λ₁. The total variance, tr Σ, equals the sum of the eigenvalues, so λ₁/tr Σ is the proportion of variance explained by the first component.
A quadratic form in x ∈ Rⁿ is Q(x) = x′Ax = ΣΣaᵢⱼxᵢxⱼ, where A can always be taken symmetric. The form x² + 4xy + y² has matrix with rows (1, 2) and (2, 1): the off-diagonal entries are half the cross-product coefficient. A symmetric matrix (and its form) is positive definite if x′Ax > 0 for all x ≠ 0, positive semi-definite (non-negative definite) if x′Ax ≥ 0 for all x, negative definite or semi-definite with the inequalities reversed, and indefinite if the form takes both signs. Equivalent tests: A is positive definite if and only if all its eigenvalues are positive, if and only if all its leading principal minors are positive (Sylvester's criterion), and if and only if A = B′B for some non-singular B (the Cholesky factorisation A = LL′ with L lower triangular exists and is unique with positive diagonal). Positive semi-definite matrices have non-negative eigenvalues. The form 2x² + 2xy + 2y², with matrix rows (2, 1) and (1, 2) and eigenvalues 3 and 1, is positive definite; x² + 4xy + y², with eigenvalues 3 and −1, is indefinite. A positive definite matrix has positive diagonal entries, a positive determinant and a positive definite inverse.
Every covariance matrix is positive semi-definite, because Var(a′x) = a′Σa ≥ 0 for every a; it is positive definite unless some linear combination of the variables has zero variance, that is, unless the variables are exactly linearly dependent. X′X is positive semi-definite, and positive definite when X has full column rank. Sums of squares in analysis of variance are quadratic forms in the observations: the total corrected sum of squares is y′(I − J/n)y, where J is the matrix of ones, and the residual sum of squares is y′(I − H)y. The expected value of a quadratic form, E(x′Ax) = tr(AΣ) + μ′Aμ for x with mean μ and covariance Σ, gives the expected mean squares used to construct F-tests and to estimate variance components. Optimisation also relies on definiteness: a critical point of a function is a minimum where the Hessian is positive definite, which is why the log-likelihood's Hessian must be negative definite at a maximum.
A square matrix P is idempotent if P² = P. Its eigenvalues are 0 or 1, because Px = λx implies λ²x = λx; and its rank equals its trace, a fact used constantly to obtain degrees of freedom. If P is idempotent, so is I − P, and rank(I − P) = n − rank(P). The matrix J/n, all of whose entries are 1/n, is symmetric and idempotent with rank 1 and trace 1; it replaces every element of a vector by the mean. The centring matrix C = I − J/n is symmetric idempotent with rank n − 1 and trace n − 1, and Cy is the vector of deviations from the mean, so y′Cy = Σ(yᵢ − ȳ)², the corrected sum of squares with its n − 1 degrees of freedom.
A symmetric idempotent matrix is an orthogonal projection matrix: it projects vectors onto its column space along the orthogonal complement. In regression with an n × p design matrix X of full column rank, the hat matrix H = X(X′X)⁻¹X′ projects y onto the column space of X, giving fitted values ŷ = Hy; I − H projects onto the orthogonal complement, giving residuals e = (I − H)y. H is symmetric and idempotent with trace and rank p, and I − H has trace and rank n − p, the residual degrees of freedom. HX = X and (I − H)X = 0, so residuals are orthogonal to every column of X. The diagonal elements hᵢᵢ of H, called leverages, lie between 0 and 1, average p/n and identify observations with unusual combinations of explanatory variables; observations with leverage well above the average, often taken as above 2p/n, deserve checking. When X is not of full rank, H = X(X′X)⁻X′ is the same for every generalised inverse.
The link between idempotence and the chi-square distribution underlies all normal-theory tests. If x is distributed as Nₙ(0, I) and A is symmetric, x′Ax has a chi-square distribution if and only if A is idempotent, and its degrees of freedom are rank(A) = tr(A). More generally, for x ~ Nₙ(μ, Σ), x′Ax is non-central chi-square with r degrees of freedom and non-centrality μ′Aμ when AΣ is idempotent with rank r. Two quadratic forms x′Ax and x′Bx with x ~ N(0, I) are independent if and only if AB = 0 (Craig's theorem), and a linear form Bx is independent of x′Ax if BA = 0. Cochran's theorem states that if I = A₁ + … + Aₖ with ranks rᵢ summing to n, the quadratic forms x′Aᵢx are independent chi-squares with rᵢ degrees of freedom. This is the theorem that makes treatment, block and error sums of squares in a designed experiment independent, and their ratios F-distributed.
A matrix that is singular or not square has no ordinary inverse, but it has generalised inverses. A g-inverse of an m × n matrix A is any n × m matrix G satisfying AGA = A. G is not unique in general. Its key property is that when Ax = b is consistent, x = Gb is a solution, and the general solution is x = Gb + (I − GA)z for arbitrary z. A simple way to construct a g-inverse: find a non-singular r × r submatrix M of A, where r is the rank, invert it, transpose, put the result in the positions of M and zeros elsewhere, and transpose back; for a symmetric matrix this is simply inverting the non-singular block in place. If A has full column rank, (A′A)⁻¹A′ is a g-inverse (a left inverse), and if A is square and non-singular, the only g-inverse is A⁻¹. AG and GA are idempotent, and rank(AG) = rank(A).
Among all g-inverses, the Moore–Penrose inverse A⁺ is unique, defined by four conditions: AA⁺A = A; A⁺AA⁺ = A⁺; (AA⁺)′ = AA⁺; and (A⁺A)′ = A⁺A. It can be computed from the singular value decomposition: if A = UΣV′, then A⁺ = VΣ⁺U′, where Σ⁺ inverts the non-zero singular values and transposes. For a non-zero column vector a, a⁺ = a′/(a′a); for a = (1, 2, 2)′, a′a = 9, so a⁺ = (1, 2, 2)/9. For a diagonal matrix, the Moore–Penrose inverse inverts the non-zero diagonal entries and leaves zeros as zeros. When A has full column rank, A⁺ = (A′A)⁻¹A′; when it has full row rank, A⁺ = A′(AA′)⁻¹. The Moore–Penrose solution x = A⁺b of a consistent system is the solution of minimum length, and for an inconsistent system it is the least-squares solution of minimum length.
In the less-than-full-rank model y = Xβ + ε, the normal equations X′Xβ = X′y have solutions β° = (X′X)⁻X′y for any g-inverse. Different g-inverses give different β°, but a linear function λ′β is estimable, meaning it has a linear unbiased estimator, if and only if λ′ lies in the row space of X, equivalently λ′(X′X)⁻X′X = λ′; for estimable functions λ′β° is the same for every solution and is the best linear unbiased estimator (Gauss–Markov). In the one-way model, τ₁ − τ₂ is estimable and estimated by ȳ₁ − ȳ₂ whatever g-inverse is used, while τ₁ alone is not. Software handles rank deficiency either by g-inverses or by constraints such as setting one effect to zero; the estimable quantities, fitted values and sums of squares do not depend on the choice.
Every real m × n matrix A, square or not, of any rank, can be written A = UΣV′, where U (m × m) and V (n × n) are orthogonal and Σ is an m × n matrix with non-negative entries σ₁ ≥ σ₂ ≥ … ≥ 0, the singular values, on its diagonal and zeros elsewhere. The singular values are the square roots of the eigenvalues of A′A (or of AA′); the columns of V are eigenvectors of A′A and the columns of U eigenvectors of AA′. The number of non-zero singular values is the rank of A. For the matrix with rows (1, 1) and (1, 1), A′A has rows (2, 2) and (2, 2) with eigenvalues 4 and 0, so the singular values are 2 and 0 and the rank is 1. For a symmetric positive semi-definite matrix the SVD coincides with the spectral decomposition, and for a general symmetric matrix the singular values are the absolute values of the eigenvalues. The compact form A = Σσᵢuᵢvᵢ′ writes A as a sum of rank-one matrices.
The Eckart–Young theorem states that the best rank-k approximation to A, in the sense of least sum of squared differences (Frobenius norm), is obtained by keeping the k largest singular values and their vectors: Aₖ = Σᵢ₌₁ᵏ σᵢuᵢvᵢ′, with squared error equal to the sum of the squares of the discarded singular values. This is the mathematics behind principal component analysis computed from a centred data matrix (the right singular vectors are the loadings and the squared singular values divided by n − 1 are the component variances), correspondence analysis, the biplot, latent semantic analysis and the AMMI model for genotype-by-environment interaction in multi-environment crop trials, in which the interaction matrix is approximated by a few multiplicative terms. The condition number σ₁/σₙ of a matrix measures how much relative errors in data can be amplified in solving linear systems; large condition numbers of X indicate multicollinearity and numerically unstable regression. Numerical software computes least squares, ranks, pseudo-inverses and principal components through the SVD because it is the most stable decomposition.
The Kronecker product of an m × n matrix A and a p × q matrix B is the mp × nq block matrix A ⊗ B whose (i, j) block is aᵢⱼB. It is not commutative in general, but it obeys useful rules: (A ⊗ B)(C ⊗ D) = AC ⊗ BD when the products are defined; (A ⊗ B)′ = A′ ⊗ B′; (A ⊗ B)⁻¹ = A⁻¹ ⊗ B⁻¹ for non-singular A and B; tr(A ⊗ B) = tr A × tr B; rank(A ⊗ B) = rank A × rank B; and, for A of order m × m and B of order n × n, det(A ⊗ B) = (det A)ⁿ(det B)ᵐ. The eigenvalues of A ⊗ B are all products λᵢμⱼ of eigenvalues of A and B. If A is 2 × 2 with determinant 3 and B is 3 × 3 with determinant 2, det(A ⊗ B) = 3³ × 2² = 108.
The vec operator stacks the columns of a matrix into a single column vector: for an m × n matrix A, vec(A) has length mn. The fundamental identity is vec(ABC) = (C′ ⊗ A)vec(B), which converts matrix equations into ordinary linear systems; for example, the matrix equation AXB = C becomes (B′ ⊗ A)vec(X) = vec(C). Other identities include tr(A′B) = vec(A)′vec(B) and, for symmetric matrices, the vech operator, which stacks only the elements on and below the diagonal, removing duplicated entries; it is used when differentiating with respect to covariance matrices.
Designed experiments and multivariate models have natural Kronecker structure. In a two-way layout with a levels of A, b levels of B and r replicates, the design and contrast matrices can be written as Kronecker products of simpler one-factor matrices, and sums of squares for main effects and interaction correspond to projection matrices of the form P_A ⊗ J_b/b ⊗ J_r/r and so on. The covariance of a multivariate regression with n independent observations each with p correlated responses is Σ ⊗ Iₙ (or Iₙ ⊗ Σ, depending on the order of stacking). Separable spatial and spatio-temporal covariance models, and the covariance of observations in multi-environment trials where genotypic effects are correlated across environments, are Kronecker products, and the rules above let large inverses and determinants be computed from small ones, which makes such models feasible to fit.
A statistics student uses matrix theory to make sense of linear models and multivariate analysis: writing an experiment or regression as y = Xβ + ε, finding the rank of X to see which effects are estimable, computing (X′X)⁻¹X′y, showing that the hat matrix is idempotent with trace p, deriving the chi-square distribution of sums of squares, and obtaining principal components from a spectral decomposition. In examinations, the student must find ranks, determinants, eigenvalues and inverses quickly, decide definiteness from minors or eigenvalues, and recognise idempotent and projection matrices.
A statistician five years into practice relies on linear algebra when software output needs interpretation: a warning that X′X is singular means aliased effects in an unbalanced trial; large leverages identify plots or farms with unusual covariate combinations; a high condition number warns of collinearity among weather variables; and Kronecker-structured covariance models make multi-environment trials tractable. Understanding which quantities are invariant to the choice of constraint or g-inverse prevents misreporting of treatment effects from rank-deficient models.
A senior statistician with twenty-five years of experience reviews analyses for agencies and research programmes, decides whether a proposed model is identifiable before data are collected, designs experiments whose information matrices give the required precision for key contrasts, and guides colleagues in choosing decompositions that are numerically stable for large data. Such a person explains to scientists, without matrices, why some effects cannot be estimated from a design and what extra observations would make them estimable.
A researcher with a PhD and long experience uses matrix theory as the core tool of methodological work: deriving properties of new estimators in linear and mixed models, constructing optimal designs through properties of information matrices, developing algorithms based on QR, Cholesky and singular value decompositions for very large genomic and spatial data, and studying random matrices in high-dimensional statistics. Such researchers teach advanced linear models, examine theses and set the matrix-theory questions of fellowship examinations.
This manual covers matrix theory to the level used in linear models, multivariate analysis and graduate entrance examinations. It does not cover abstract linear maps between general vector spaces in depth, Jordan canonical form, matrix calculus (differentiation with respect to matrices) beyond a mention, matrix norms in general, or numerical algorithms in detail. Proofs are sketched rather than given in full.
Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions.
These notes are free to read here, with no account.
What the Laisha app adds:
Free for ten days, then a subscription. We would rather write that here than have you find it later.