Linear algebra refresher

Bold lowercase letters are column vectors, bold capitals are matrices or random vectors, a prime means transpose and everything is real. The determinant of $\mathbf A$ is written $\lvert \mathbf A \rvert$ and $\mathbf 1$ is a vector of ones. The covariance matrix of a random vector $\mathbf X$ is $\boldsymbol\Sigma,$ and for observations $\mathbf x_1, \dots, \mathbf x_n$ with mean $\bar{\mathbf x}$ the sample covariance matrix is

The multivariate statistics reference has the results built on these facts.

Vectors

For ${\mathbf x = (x_1, \dots, x_n)'}$ and ${\mathbf y = (y_1, \dots, y_n)'}$ in $\mathbb R^n,$

Projection

The projection of $\mathbf x$ on a nonzero vector $\mathbf y$ is

Its length is $\lvert \mathbf x'\mathbf y \rvert / \lVert \mathbf y \rVert = \lVert \mathbf x \rVert\,\lvert \cos\theta \rvert,$ and what is left over, ${\mathbf x - \operatorname{proj}_{\mathbf y}\mathbf x},$ is orthogonal to $\mathbf y.$

Projecting a column of data $\mathbf y$ on $\mathbf 1$ gives $\bar y\,\mathbf 1$ and leaves the deviation vector ${\mathbf d = \mathbf y - \bar y\,\mathbf 1}.$ Its squared length $\mathbf d'\mathbf d = \sum_j (y_j -\nobreak \bar y)^2$ is ${n - 1}$ times the sample variance, and the sample correlation of two variables is the cosine of the angle between their deviation vectors.

Cauchy-Schwarz inequality

Equality holds exactly when one vector is a multiple of the other. This is why ${-1 \le \cos\theta \le 1},$ and so why correlations lie between $-1$ and $1.$ The extended version puts a positive definite matrix in the middle.

Linear independence

Vectors $\mathbf x_1, \dots, \mathbf x_k$ are linearly independent when $c_1\mathbf x_1 + \dots + c_k\mathbf x_k =\nobreak \mathbf 0$ only for ${c_1 = \dots = c_k = 0}.$ Otherwise they are dependent, and at least one of them is a linear combination of the others.

• More than $n$ vectors in $\mathbb R^n$ are always dependent.
• Nonzero vectors that are mutually orthogonal are independent.
• The columns of a matrix $\mathbf A$ are independent exactly when ${\mathbf A\mathbf c = \mathbf 0}$ only for ${\mathbf c = \mathbf 0}.$

Matrices

Products and special matrices

If $\mathbf A$ is ${m \times k}$ and $\mathbf B$ is ${k \times n},$ then $\mathbf A\mathbf B$ is ${m \times n},$ and its entry in row $i$ and column $j$ is

Order matters, and ${\mathbf A\mathbf B \ne \mathbf B\mathbf A}$ in general, even when both products exist.

• The identity $\mathbf I$ has ones on the diagonal and zeros elsewhere, and ${\mathbf I\mathbf A = \mathbf A\mathbf I = \mathbf A}.$
• A diagonal matrix $\operatorname{diag}(d_1, \dots, d_n)$ has zeros everywhere off the diagonal.
• A symmetric matrix has ${\mathbf A' = \mathbf A},$ so ${a_{ij} = a_{ji}}.$ The matrices $\boldsymbol\Sigma,$ $\mathbf S$ and every correlation matrix are symmetric, and so are $\mathbf X'\mathbf X$ and $\mathbf X\mathbf X'$ for any matrix $\mathbf X,$ and $\mathbf C\mathbf A\mathbf C'$ whenever $\mathbf A$ is symmetric.

Transpose and inverse

A square matrix $\mathbf A$ is invertible, or nonsingular, when some matrix $\mathbf A^{-1}$ has $\mathbf A\mathbf A^{-1} = \mathbf A^{-1}\mathbf A =\nobreak \mathbf I.$ The inverse rules below are for square matrices of the same size that are invertible, and a number ${c \ne 0}.$

In products the order reverses. A number is its own transpose, so ${\mathbf x'\mathbf A\mathbf y = \mathbf y'\mathbf A'\mathbf x}.$ The inverse of a symmetric matrix is symmetric, and for invertible $\mathbf A$ the system ${\mathbf A\mathbf x = \mathbf b}$ has exactly one solution, ${\mathbf x = \mathbf A^{-1}\mathbf b}.$ For a ${2 \times 2}$ matrix with ${ad - bc \ne 0},$

Orthogonal matrices

A square matrix $\mathbf Q$ is orthogonal when ${\mathbf Q'\mathbf Q = \mathbf Q\mathbf Q' = \mathbf I},$ so ${\mathbf Q^{-1} = \mathbf Q'}.$ Its columns have length 1 and are mutually orthogonal, and so are its rows. It keeps lengths and angles, since $(\mathbf Q\mathbf x)'(\mathbf Q\mathbf y) = \mathbf x'\mathbf y,$ so geometrically it is a rotation, possibly combined with a reflection, and ${\lvert \mathbf Q \rvert = \pm 1}.$

Principal components use an orthogonal matrix whose columns are the eigenvectors of $\boldsymbol\Sigma,$ so the components are the original variables seen along rotated axes.

Determinant, trace and rank

Determinant

For a ${2 \times 2}$ matrix

For larger matrices, expanding along any row $i$ gives

where $\mathbf M_{ij}$ is $\mathbf A$ with row $i$ and column $j$ deleted. Columns work the same way.

For ${n \times n}$ matrices $\mathbf A$ and $\mathbf B$ and a number $c,$

Watch for two traps. The scalar rule has $c^n,$ not $c,$ because each of the $n$ rows is multiplied by $c.$ And ${\lvert \mathbf A + \mathbf B \rvert}$ is not ${\lvert \mathbf A \rvert + \lvert \mathbf B \rvert}$ in general.

Geometrically, $\lvert \mathbf A \rvert$ is the factor by which $\mathbf A$ scales volume, with a sign for orientation. A zero determinant means some direction gets squashed flat.

Trace

The trace of a square matrix is the sum of its diagonal entries, $\operatorname{tr}(\mathbf A) = a_{11} + a_{22} + \dots + a_{nn}.$

The cyclic rule only lets you move a factor from one end of the product to the other. In general $\operatorname{tr}(\mathbf A\mathbf B\mathbf C)$ and $\operatorname{tr}(\mathbf A\mathbf C\mathbf B)$ differ, and $\operatorname{tr}(\mathbf A\mathbf B)$ is not $\operatorname{tr}(\mathbf A)\operatorname{tr}(\mathbf B).$

A trick that comes up again and again is writing a quadratic form, which is a single number, as a trace.

The first step holds because a ${1 \times 1}$ matrix is its own trace, and the second is the cyclic rule. It collapses a sum of quadratic forms into one trace, for example

Rank

The rank of a matrix is the largest number of linearly independent columns it has, which always equals the largest number of linearly independent rows. For an ${m \times n}$ matrix $\mathbf A$ the rank is at most $\min(m, n),$ and

Whenever the product $\mathbf A\mathbf B$ is defined,

Multiplying by an invertible matrix leaves the rank unchanged, and for nonzero vectors the outer product $\mathbf x\mathbf y'$ has rank 1.

For a square ${n \times n}$ matrix $\mathbf A$ these all say the same thing. The matrix $\mathbf A$ is invertible, ${\lvert \mathbf A \rvert \ne 0},$ ${\operatorname{rank}(\mathbf A) = n}$ (full rank), the columns are linearly independent, ${\mathbf A\mathbf x = \mathbf 0}$ only for ${\mathbf x = \mathbf 0},$ and no eigenvalue is 0.

A sample covariance matrix built from $n$ observations has rank at most ${n - 1},$ because the deviations ${\mathbf x_j - \bar{\mathbf x}}$ add up to $\mathbf 0.$ So with $p$ variables, $\mathbf S$ is singular whenever ${n \le p}.$

Eigenvalues

Definition

A number $\lambda$ is an eigenvalue of a square matrix $\mathbf A,$ with eigenvector $\mathbf x,$ when

So $\mathbf A$ only stretches $\mathbf x$ by the factor $\lambda$ and keeps it on the same line. Any nonzero multiple of an eigenvector is again an eigenvector, so it is common to scale to length 1 and write ${\mathbf e = \mathbf x / \lVert \mathbf x \rVert}.$ The sign is still a free choice, which is why two programs can report eigenvectors with opposite signs.

Characteristic equation

The eigenvalues are the solutions of

For an ${n \times n}$ matrix this is a polynomial equation of degree $n,$ so there are $n$ eigenvalues counting repeats. Some may be complex numbers when $\mathbf A$ is not symmetric. For each eigenvalue, the eigenvectors are the nonzero solutions of ${(\mathbf A - \lambda\mathbf I)\mathbf x = \mathbf 0}.$ For a ${2 \times 2}$ matrix the equation is

Useful facts

For any ${n \times n}$ matrix $\mathbf A$ with eigenvalues $\lambda_1, \dots, \lambda_n,$ repeats and complex values included,

The last row is why $\mathbf W^{-1}\mathbf B$ in MANOVA, with $\mathbf W$ positive definite and $\mathbf B$ symmetric, has the same eigenvalues as the symmetric matrix $\mathbf W^{-1/2}\mathbf B\mathbf W^{-1/2}$ from functions of a matrix, so they are real.

Symmetric matrices

When $\mathbf A$ is symmetric, as every covariance matrix is, the picture is much cleaner.

• Every eigenvalue is real.
• Eigenvectors for different eigenvalues are orthogonal.
• There is always a full set of $n$ orthonormal eigenvectors, meaning length 1 and mutually orthogonal, even when eigenvalues repeat. With a repeated eigenvalue they are not unique, since any orthonormal basis of the space spanned by that eigenvalue’s eigenvectors works.
• The rank is the number of nonzero eigenvalues, counted with multiplicity.

The orthogonality takes one line. If ${\mathbf A\mathbf e_i = \lambda_i\mathbf e_i}$ and ${\mathbf A\mathbf e_j = \lambda_j\mathbf e_j},$ then

using ${\mathbf A' = \mathbf A}$ in the middle, so ${\mathbf e_j'\mathbf e_i = 0}$ whenever ${\lambda_i \ne \lambda_j}.$

Worked example

Take

Here ${\operatorname{tr}(\mathbf A) = 7}$ and ${\lvert \mathbf A \rvert = 10 - 4 = 6},$ so the characteristic equation is

and the eigenvalues are ${\lambda_1 = 6}$ and ${\lambda_2 = 1}.$

For ${\lambda_1 = 6}.$ The system ${(\mathbf A - 6\,\mathbf I)\mathbf x = \mathbf 0}$ is

Both rows say ${x_1 = 2x_2},$ so ${\mathbf x = (2, 1)'}$ works and ${\mathbf e_1 = \frac{1}{\sqrt 5}(2, 1)'}.$

For ${\lambda_2 = 1}.$ Now ${(\mathbf A - \mathbf I)\mathbf x = \mathbf 0}$ has rows ${4x_1 + 2x_2 = 0}$ and ${2x_1 + x_2 = 0}.$ Both say ${x_2 = -2x_1},$ so ${\mathbf x = (1, -2)'}$ works and ${\mathbf e_2 = \frac{1}{\sqrt 5}(1, -2)'}.$

Checks.

The example continues in the next section.

Spectral decomposition

The decomposition

Every symmetric ${p \times p}$ matrix $\mathbf A$ can be written as

where $\lambda_1, \dots, \lambda_p$ are its eigenvalues, $\mathbf e_1, \dots, \mathbf e_p$ are matching orthonormal eigenvectors, ${\mathbf P = [\mathbf e_1 \ \cdots \ \mathbf e_p]}$ has them as columns and $\boldsymbol\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_p).$ Since $\mathbf P$ is orthogonal, the same fact reads ${\mathbf P'\mathbf A\mathbf P = \boldsymbol\Lambda}.$

The sum builds $\mathbf A$ from $p$ rank one pieces, one per eigenvector, each weighted by its eigenvalue.

Functions of a matrix

Apply the function to the eigenvalues and keep the eigenvectors.

For a positive definite $\mathbf A,$ which for a symmetric matrix means every eigenvalue is positive (see quadratic forms), all four exist and the square roots behave the way the notation suggests. The root $\mathbf A^{1/2}$ is symmetric and positive definite, and

For a positive definite covariance matrix this gives $\boldsymbol\Sigma^{-1/2},$ the matrix that standardizes a random vector.

Back to the example

For the example matrix $\mathbf A$ with eigenvalues ${\lambda_1 = 6}$ and ${\lambda_2 = 1},$ the rank one pieces are

and they rebuild $\mathbf A.$

The inverse comes out the same way and agrees with the ${2 \times 2}$ formula.

So does the square root.

Squaring it gives back $\mathbf A,$ because $(\mathbf e_i\mathbf e_i')(\mathbf e_j\mathbf e_j')$ is $\mathbf e_i\mathbf e_i'$ when ${i = j}$ and $\mathbf 0$ otherwise.

Projection matrices

A symmetric matrix $\mathbf H$ with ${\mathbf H\mathbf H = \mathbf H}$ is a projection matrix. Its eigenvalues are all 0 or 1, so its rank equals its trace. Two come up often.

• The centering matrix ${\mathbf C = \mathbf I - \frac{1}{n}\mathbf 1\mathbf 1'}$ turns a data column into its deviation vector, ${\mathbf C\mathbf y = \mathbf y - \bar y\,\mathbf 1}.$ It has rank ${n - 1},$ and $(n-1)\,\mathbf S = \mathbf X'\mathbf C\mathbf X$ for an ${n \times p}$ data matrix $\mathbf X.$
• For $\mathbf Z$ with $r$ linearly independent columns, ${\mathbf H = \mathbf Z(\mathbf Z'\mathbf Z)^{-1}\mathbf Z'}$ projects onto the span of those columns and gives the fitted values ${\hat{\mathbf y} = \mathbf H\mathbf y}$ in regression. It has rank $r.$

Quadratic forms

The form

For a symmetric ${p \times p}$ matrix $\mathbf A,$ the quadratic form is the number

For ${p = 2}$ it is $a_{11}x_1^2 + 2a_{12}x_1x_2 + a_{22}x_2^2.$ A nonsymmetric $\mathbf A$ gives the same form as its symmetric part ${\frac{1}{2}(\mathbf A + \mathbf A')},$ so assuming symmetry costs nothing.

With the spectral decomposition and ${\mathbf y = \mathbf P'\mathbf x},$

In the eigenvector coordinates a quadratic form is just a weighted sum of squares, so the signs of the eigenvalues decide its sign.

Definite matrices

Positive semidefinite is also called nonnegative definite. For a symmetric $\mathbf A,$ each of these is equivalent to being positive definite.

• every eigenvalue is positive
• it can be written as ${\mathbf A = \mathbf B'\mathbf B}$ for some matrix $\mathbf B$ with linearly independent columns
• every leading principal minor is positive, meaning the determinants of the top left blocks of sizes ${1 \times 1}, {2 \times 2}, \dots, {p \times p}$

For a ${2 \times 2}$ matrix the last test reads ${a_{11} > 0}$ and ${a_{11}a_{22} - a_{12}^2 > 0}.$ The example matrix passes with ${5 > 0}$ and ${10 - 4 = 6 > 0},$ matching its eigenvalues 6 and 1.

For positive semidefinite, the first two carry over with ${\ge 0}$ and any matrix $\mathbf B,$ but the minor test does not. The matrix $\begin{bmatrix} 0 & 0 \\ 0 & -1 \end{bmatrix}$ has leading minors 0 and 0 but is not positive semidefinite.

A few more facts.

• A positive definite matrix is invertible, has ${\lvert \mathbf A \rvert > 0}$ and positive diagonal entries, and its inverse is positive definite.
• A positive semidefinite matrix is positive definite exactly when it is invertible.
• If $\mathbf A$ is positive semidefinite, so is $\mathbf C\mathbf A\mathbf C'.$ If $\mathbf A$ is positive definite and the rows of $\mathbf C$ are linearly independent, $\mathbf C\mathbf A\mathbf C'$ is positive definite. This is why $\mathbf C\mathbf S\mathbf C'$ in a test of contrasts can be inverted when $\mathbf S$ is positive definite.

Distance and ellipsoids

For a positive definite $\mathbf A,$ the number

is a distance between the points $\mathbf x$ and $\mathbf y.$ The points with

form an ellipsoid centered at $\boldsymbol\mu,$ with axes along the eigenvectors $\mathbf e_i$ of $\mathbf A$ and half lengths $c/\sqrt{\lambda_i}.$ When ${\mathbf A = \boldsymbol\Sigma^{-1}}$ the eigenvectors are those of $\boldsymbol\Sigma,$ the eigenvalues are the reciprocals of those of $\boldsymbol\Sigma,$ and the half lengths become $c\sqrt{\lambda_i}$ with $\lambda_i$ now the eigenvalues of $\boldsymbol\Sigma.$ That is the picture behind density contours and Mahalanobis distance.

Covariance matrices

For any constant vector $\mathbf a,$

so every covariance matrix is positive semidefinite. It fails to be positive definite exactly when ${\operatorname{Var}(\mathbf a'\mathbf X) = 0}$ for some ${\mathbf a \ne \mathbf 0}.$ Then that combination of the variables is constant with probability 1, so any variable with a nonzero coefficient in $\mathbf a$ is determined exactly by the others through a linear formula. The same argument covers $\mathbf S,$ because $\mathbf a'\mathbf S\mathbf a$ is the sample variance of the numbers $\mathbf a'\mathbf x_1, \dots, \mathbf a'\mathbf x_n.$ Correlation matrices are positive semidefinite for the same reason.

Maximizing a quadratic form

Let $\mathbf B$ be symmetric ${p \times p}$ with eigenvalues $\lambda_1 \ge \lambda_2 \ge \dots \ge\nobreak \lambda_p$ and orthonormal eigenvectors $\mathbf e_1, \dots, \mathbf e_p.$ Then

and among nonzero $\mathbf x$ perpendicular to the first $k$ eigenvectors, for $k = 1, 2, \dots, p - 1,$

The ratio does not change when $\mathbf x$ is rescaled, so these are the same as maximizing $\mathbf x'\mathbf B\mathbf x$ over unit vectors ${\mathbf x'\mathbf x = 1},$ and any nonzero multiple of $\mathbf e_1$ also gives the maximum.

Why it works. Write ${\mathbf x = \mathbf P\mathbf y}.$ Then ${\mathbf x'\mathbf x = \mathbf y'\mathbf y}$ and

with equality at ${\mathbf y = (1, 0, \dots, 0)'},$ which is ${\mathbf x = \mathbf e_1}.$ Being perpendicular to $\mathbf e_1, \dots, \mathbf e_k$ forces ${y_1 = \dots = y_k = 0},$ and the same argument then gives $\lambda_{k+1}.$

Principal components. Take ${\mathbf B = \boldsymbol\Sigma}.$ Since $\operatorname{Var}(\mathbf a'\mathbf X) = \mathbf a'\boldsymbol\Sigma\mathbf a,$ the unit vector $\mathbf a$ with the largest variance is $\mathbf e_1,$ and the variance is $\lambda_1.$ That is the first principal component. Because $\operatorname{Cov}(\mathbf a'\mathbf X, \mathbf e_1'\mathbf X) = \mathbf a'\boldsymbol\Sigma\mathbf e_1 = \lambda_1\,\mathbf a'\mathbf e_1,$ being uncorrelated with the first component means ${\mathbf a \perp \mathbf e_1}$ when ${\lambda_1 > 0},$ so the second component is $\mathbf e_2'\mathbf X$ with variance $\lambda_2,$ and so on.

A ratio of two forms. For symmetric $\mathbf B$ and positive definite $\mathbf W,$

is the largest eigenvalue of $\mathbf W^{-1}\mathbf B,$ reached at a matching eigenvector. Substituting ${\mathbf y = \mathbf W^{1/2}\mathbf x}$ turns it into the previous problem with the symmetric matrix $\mathbf W^{-1/2}\mathbf B\mathbf W^{-1/2},$ which has the same eigenvalues as $\mathbf W^{-1}\mathbf B.$ This is the problem behind Fisher’s discriminant functions for several groups.

Extended Cauchy-Schwarz inequality

Let $\mathbf B$ be positive definite. For any vectors $\mathbf b$ and $\mathbf d,$

with equality exactly when ${\mathbf b = c\,\mathbf B^{-1}\mathbf d}$ or ${\mathbf d = c\,\mathbf B\mathbf b}$ for some number $c.$ It is the ordinary inequality applied to $\mathbf B^{1/2}\mathbf b$ and $\mathbf B^{-1/2}\mathbf d,$ and with ${\mathbf B = \mathbf I}$ it is the ordinary inequality.

The version used most is a maximization. For positive definite $\mathbf B$ and a fixed ${\mathbf d \ne \mathbf 0},$

reached at ${\mathbf x = c\,\mathbf B^{-1}\mathbf d}$ for any ${c \ne 0}.$

• With ${\mathbf B = \mathbf S}$ and ${\mathbf d = \bar{\mathbf x} - \boldsymbol\mu_0},$ it says Hotelling’s $T^2$ is the largest squared $t$ statistic over all combinations $\mathbf a'\bar{\mathbf x},$ which is where the simultaneous intervals come from.
• With ${\mathbf B = \mathbf S_{\text{pooled}}}$ and ${\mathbf d = \bar{\mathbf x}_1 - \bar{\mathbf x}_2},$ it gives Fisher’s linear discriminant $\mathbf a = \mathbf S_{\text{pooled}}^{-1}(\bar{\mathbf x}_1 - \bar{\mathbf x}_2).$

Partitioned matrices

Split a ${p \times p}$ matrix after its first $q$ rows and columns, so $\mathbf A_{11}$ is ${q \times q}$ and $\mathbf A_{22}$ is ${(p - q) \times (p - q)}.$

Blocks multiply like entries, as long as the sizes fit and you keep the order of each product. The transpose $\mathbf A'$ has $\mathbf A_{11}'$ and $\mathbf A_{22}'$ on the diagonal, $\mathbf A_{21}'$ top right and $\mathbf A_{12}'$ bottom left, so a symmetric $\mathbf A$ has ${\mathbf A_{21} = \mathbf A_{12}'}.$

Block determinant

The first row comes from the factorization

in which the first factor has determinant 1 and the second is block triangular.

Block inverse

When $\mathbf A_{11}$ is invertible, the matrix

is the Schur complement of $\mathbf A_{11}.$ By the determinant rule, $\mathbf A$ is invertible exactly when $\mathbf A_{22\cdot 1}$ is, and then

Swapping the roles of the blocks gives the version for invertible $\mathbf A_{22},$ with

in place of $\mathbf A_{22\cdot 1}.$ So when $\mathbf A,$ $\mathbf A_{11}$ and $\mathbf A_{22}$ are all invertible, the bottom right block of $\mathbf A^{-1}$ is $\mathbf A_{22\cdot 1}^{-1}$ and the top left block is $\mathbf A_{11\cdot 2}^{-1},$ inverses of Schur complements rather than of the matching blocks of $\mathbf A.$ If ${\mathbf A_{12} = \mathbf 0}$ and ${\mathbf A_{21} = \mathbf 0},$ the inverse is block diagonal with blocks $\mathbf A_{11}^{-1}$ and $\mathbf A_{22}^{-1}.$

In statistics. With ${\mathbf A = \boldsymbol\Sigma}$ positive definite, the Schur complement

is the covariance of the first block $\mathbf X^{(1)}$ given the second block $\mathbf X^{(2)}$ for a multivariate normal, as in conditional distributions. The determinant rule gives $\lvert \boldsymbol\Sigma \rvert = \lvert \boldsymbol\Sigma_{22} \rvert\,\lvert \boldsymbol\Sigma_{11\cdot 2} \rvert,$ which is how the joint normal density factors into a marginal density times a conditional one.

Matrix calculus

Gradients

For a function $f$ of ${\mathbf x \in \mathbb R^p},$ the gradient $\partial f/\partial\mathbf x$ is the column of partial derivatives $\partial f/\partial x_1, \dots, \partial f/\partial x_p.$ Setting it to $\mathbf 0$ finds the stationary points. Every local maximum and minimum of a differentiable $f$ on $\mathbb R^p$ is one, but a stationary point can also be a saddle, so deciding which it is needs a separate argument, such as the concavity used for the normal likelihood.