Book contents
Contents
raw Math

A matrix \(\mathbf{A}\) transforms a vector \(\mathbf{x}\) into another vector \(\mathbf{Ax}\). In general, \(\mathbf{Ax}\) points in a different direction than \(\mathbf{x}\), since some sort of rotation may be part of the transformation.

However, special vectors called eigenvectors keep their direction under a transformation with \(\mathbf{A}\) (i.e. they remain invariant in direction). In this case, \(\mathbf{Ax}\) and \(\mathbf{x}\) are parallel, differing only in a constant scale factor (only stretched, compressed, or flipped), which we call the associated eigenvalue.

Illustration for the matrix \(\mathbf{A}=\bigl(\begin{smallmatrix}2&1\\1&2\end{smallmatrix}\bigr)\). A generic vector \(\mathbf{x}\) maps to \(\mathbf{Ax}\) pointing in a visibly different direction, marked by a small arc showing the angle change. An eigenvector \(\mathbf{v}\), by contrast, maps to \(\mathbf{Av}=3\mathbf{v}\) lying on the exact same ray from the origin, only three times as long, illustrated with a dashed guide line through both.

Definition

For a square matrix \(\mathbf{A}\), a to be found nonzero vector \(\mathbf{x}\) and a complex or real scalar \(\lambda\) are eigenvectors and associated eigenvalues, iff they satisfy

\[\mathbf{Ax}=\lambda\mathbf{x}\]

There are infinetely many solutions, since \(\frac{\lambda}{c}\) with \(c\mathbf{x}\) for any \(c\neq 0\) also satisfies this equation. Consequently, eigenvectors are assumed to be normalized, i.e., satisfy the constraint \(\mathbf{x}^T\mathbf{x}=1\).

The original equation can be re-arranged a bit:

\[\begin{array}{rrl} &\mathbf{Ax}&=\lambda\mathbf{x}\\ \Leftrightarrow&\mathbf{Ax} - \lambda\mathbf{x}&=\mathbf{0}\\ \Leftrightarrow&(\mathbf{A} - \lambda\mathbf{I})\mathbf{x}&=\mathbf{0} \end{array}\]

This is a homogeneous linear system in \(\mathbf{x}\). It always has the trivial solution \(\mathbf{x}=\mathbf{0}\), but we are explicitly interested in nonzero solutions. A nonzero solution exists if and only if the matrix \(\mathbf{A}-\lambda\mathbf{I}\) is singular, i.e. not invertible. For square matrices, singularity is equivalent to a zero determinant. Therefore, eigenvalues \(\lambda\) are precisely those scalars for which

\[\det(\mathbf{A}-\lambda\mathbf{I})=0.\]

Once an eigenvalue \(\lambda\) is known, eigenvectors are obtained by solving the corresponding linear system

\[(\mathbf{A}-\lambda\mathbf{I})\mathbf{x}=\mathbf{0}\]

and then choosing any nonzero solution and (optionally) normalizing it.

The Characteristic Polynomial

For \(\mathbf{A}\in\mathbb{C}^{n\times n}\), the expression \(\det(\mathbf{A}-\lambda\mathbf{I})\) is a polynomial in \(\lambda\) of degree \(n\). This polynomial is called the characteristic polynomial of \(\mathbf{A}\):

\[p_{\mathbf{A}}(\lambda):=\det(\mathbf{A}-\lambda\mathbf{I}).\]

The equation \(p_{\mathbf{A}}(\lambda)=0\) is called the characteristic equation. Its roots \(\lambda_1,\ldots,\lambda_n\) (counted with algebraic multiplicity) are the eigenvalues of \(\mathbf{A}\).

Some texts instead define the characteristic polynomial as \(\det(\lambda\mathbf{I}-\mathbf{A})\). Since \(\det(\lambda\mathbf{I}-\mathbf{A})=\det\bigl(-(\mathbf{A}-\lambda\mathbf{I})\bigr)=(-1)^n\det(\mathbf{A}-\lambda\mathbf{I})\), the two conventions differ only by the constant factor \((-1)^n\) and therefore have exactly the same roots, i.e. the same eigenvalues.

To see directly why a polynomial appears, it helps to look at an explicit example and observe what changes when we subtract \(\lambda\mathbf{I}\): only the diagonal entries are shifted by \(-\lambda\). When you compute the determinant, you multiply and add these entries in a structured way; as soon as \(\lambda\) appears on the diagonal, the determinant becomes an expression containing powers of \(\lambda\), i.e. a polynomial.

A concrete \(2\times 2\) example

Let

\[\mathbf{A}=\begin{pmatrix}a&b\\c&d\end{pmatrix}.\]

Then

\[\mathbf{A}-\lambda\mathbf{I}= \begin{pmatrix}a-\lambda&b\\c&d-\lambda\end{pmatrix}\]

and the determinant expands to

\[\det(\mathbf{A}-\lambda\mathbf{I}) =\det\begin{pmatrix}a-\lambda&b\\c&d-\lambda\end{pmatrix} =(a-\lambda)(d-\lambda)-bc.\]

Multiplying out gives

\[(a-\lambda)(d-\lambda)-bc =\lambda^2-(a+d)\lambda+(ad-bc).\]

So the characteristic polynomial is

\[p_{\mathbf{A}}(\lambda)=\lambda^2-\mathrm{Tr}(\mathbf{A})\,\lambda+\det(\mathbf{A}).\]

Its roots are the eigenvalues. For higher-dimensional matrices the determinant expansion is more involved, but the same mechanism applies: subtracting \(\lambda\) on the diagonal forces the determinant to become a polynomial in \(\lambda\) whose roots are the eigenvalues.

How to compute eigenvectors once \(\lambda\) is known

For each eigenvalue \(\lambda\), form \(\mathbf{A}-\lambda\mathbf{I}\) and solve

\[(\mathbf{A}-\lambda\mathbf{I})\mathbf{x}=\mathbf{0}.\]

The solution set is the null space (kernel) \(\mathcal{N}(\mathbf{A}-\lambda\mathbf{I})\), called the eigenspace associated with \(\lambda\). Any nonzero vector in this eigenspace is an eigenvector. If the eigenspace has dimension greater than \(1\), then \(\lambda\) has multiple linearly independent eigenvectors; choosing an orthonormal basis of the eigenspace is often convenient.

A Numerical Example

The procedure above becomes concrete for the symmetric matrix

\[ \mathbf{A}=\begin{pmatrix}4&1&1\\1&4&1\\1&1&4\end{pmatrix}. \]

Expanding \(\det(\mathbf{A}-\lambda\mathbf{I})\) along the first row,

\[ p_{\mathbf{A}}(\lambda)=(4-\lambda)\bigl[(4-\lambda)^2-1\bigr]-\bigl[(4-\lambda)-1\bigr]+\bigl[1-(4-\lambda)\bigr]. \]

The last two bracketed terms are negatives of each other, so together they simplify to \(-2\bigl[(4-\lambda)-1\bigr]=2(\lambda-3)\). Using \((4-\lambda)^2-1=(3-\lambda)(5-\lambda)\), this leaves

\[ p_{\mathbf{A}}(\lambda)=(4-\lambda)(3-\lambda)(5-\lambda)-2(3-\lambda)=(3-\lambda)\bigl[(4-\lambda)(5-\lambda)-2\bigr], \]

and since \((4-\lambda)(5-\lambda)-2=\lambda^2-9\lambda+18=(\lambda-3)(\lambda-6)\), the characteristic polynomial factors completely:

\[ p_{\mathbf{A}}(\lambda)=(3-\lambda)(\lambda-3)(\lambda-6)=-(\lambda-3)^2(\lambda-6). \]

So \(\mathbf{A}\) has the eigenvalue \(\lambda_1=6\) and a repeated eigenvalue \(\lambda_2=\lambda_3=3\) with algebraic multiplicity \(2\). For \(\lambda_1=6\), row reduction gives

\[ \mathbf{A}-6\mathbf{I}=\begin{pmatrix}-2&1&1\\1&-2&1\\1&1&-2\end{pmatrix} \;\longrightarrow\; \begin{pmatrix}1&0&-1\\0&1&-1\\0&0&0\end{pmatrix}, \]

so \(x=z\) and \(y=z\); choosing \(z=1\) gives the eigenvector \(\mathbf{q}_1=(1,1,1)^T\). For \(\lambda_2=3\),

\[ \mathbf{A}-3\mathbf{I}=\begin{pmatrix}1&1&1\\1&1&1\\1&1&1\end{pmatrix} \;\longrightarrow\; \begin{pmatrix}1&1&1\\0&0&0\\0&0&0\end{pmatrix}, \]

which imposes only the single constraint \(x+y+z=0\), so the eigenspace \(\mathcal{N}(\mathbf{A}-3\mathbf{I})\) is two-dimensional: the geometric multiplicity of \(\lambda=3\) equals its algebraic multiplicity, so \(\mathbf{A}\) is diagonalisable despite the repeated eigenvalue. Two convenient eigenvectors spanning this plane are \(\mathbf{q}_2=(1,-1,0)^T\) and \(\mathbf{q}_3=(1,1,-2)^T\), chosen so that \(\mathbf{q}_2\cdot\mathbf{q}_3=0\).

Because \(\mathbf{A}\) is symmetric, \(\mathbf{q}_1\), \(\mathbf{q}_2\), and \(\mathbf{q}_3\) are pairwise orthogonal, exactly as the orthogonality property below predicts, and normalizing them gives the orthogonal matrix \(\mathbf{Q}\) of the spectral decomposition \(\mathbf{A}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T\) with \(\mathbf{\Lambda}=\mathrm{diag}(6,3,3)\).

Properties

Orthogonal Matrices

An orthogonal matrix is a real square matrix \(\mathbf{A}\) such that

\[\mathbf{A}^T=\mathbf{A}^{-1}.\]

Equivalently,

\[\mathbf{A}\mathbf{A}^T=\mathbf{I}\quad\text{and}\quad \mathbf{A}^T\mathbf{A}=\mathbf{I}.\]

This means the rows (and columns) of \(\mathbf{A}\) form an orthonormal set:

Geometrically, orthogonal matrices represent length- and angle-preserving transformations (rotations and reflections). In particular, \(\|\mathbf{Ax}\|_2=\|\mathbf{x}\|_2\) for all \(\mathbf{x}\).

This length-preservation forces every eigenvalue of an orthogonal matrix to have absolute value \(1\), even when it is complex. Writing \(\mathbf{v}^*\) for the conjugate transpose of an eigenvector \(\mathbf{Av}=\lambda\mathbf{v}\), and using that \(\mathbf{A}\) is real so \(\mathbf{A}^*=\mathbf{A}\),

\[|\lambda|^2\,\mathbf{v}^*\mathbf{v}=(\lambda\mathbf{v})^*(\lambda\mathbf{v})=(\mathbf{Av})^*(\mathbf{Av}) =\mathbf{v}^*\mathbf{A}^T\mathbf{A}\mathbf{v}=\mathbf{v}^*\mathbf{v},\]

and since \(\mathbf{v}\neq\mathbf{0}\) gives \(\mathbf{v}^*\mathbf{v}>0\), this forces \(|\lambda|^2=1\): every eigenvalue of an orthogonal matrix, real or complex, lies exactly on the unit circle.

Rotation Matrices and Complex Eigenvalues

A particularly beautiful example is the two-dimensional rotation matrix

\[ \mathbf{R}(\theta)= \begin{pmatrix} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{pmatrix}. \]

This matrix rotates every vector in \(\mathbb{R}^2\) by the angle \(\theta\) about the origin. Since rotations preserve lengths and angles, \(\mathbf{R}(\theta)\) is orthogonal and satisfies

\[ \mathbf{R}(\theta)^T\mathbf{R}(\theta)=\mathbf{I}. \]

This single matrix equation already reveals why a \(2\times 2\) rotation has only one free parameter: the right-hand side \(\mathbf{I}\) is symmetric, so it imposes only \(3\) independent scalar constraints (two diagonal entries and one off-diagonal entry) on the \(4\) entries of \(\mathbf{R}(\theta)\), leaving exactly \(4-3=1\) degree of freedom, matched by the single angle \(\theta\). The same counting argument applied to an \(n\times n\) orthogonal matrix gives \(n^2-\tfrac{n(n+1)}{2}=\tfrac{n(n-1)}{2}\) degrees of freedom; for \(n=3\) this is \(3\), which is why every three-dimensional rotation can be built by composing three elementary rotations about fixed coordinate axes. Requiring \(\det(\mathbf{R})=+1\) (rather than \(-1\)) is what singles out genuine rotations among all length-preserving orthogonal transformations, excluding reflections.

To find its eigenvalues, compute the characteristic polynomial:

\[ \det(\mathbf{R}(\theta)-\lambda\mathbf{I}) = \det\begin{pmatrix} \cos\theta-\lambda & -\sin\theta\\ \sin\theta & \cos\theta-\lambda \end{pmatrix}. \]

Expanding the determinant gives

\[ (\cos\theta-\lambda)^2+\sin^2\theta = 0, \]

hence

\[ \lambda^2 - 2\cos\theta\,\lambda + 1 = 0. \]

The two eigenvalues are therefore

\[ \lambda_{1,2}=\cos\theta \pm i\sin\theta = e^{\pm i\theta}. \]

This shows the eigenvalues of a planar rotation are themselves points on the unit circle in the complex plane. In other words, the rotation reappears inside its own spectrum.

To see why this is so natural, identify a vector \(\begin{pmatrix}x\\y\end{pmatrix}\in\mathbb{R}^2\) with the complex number \(z=x+iy\). Then applying the matrix \(\mathbf{R}(\theta)\) is exactly the same as multiplying by \(e^{i\theta}\):

\[ \mathbf{R}(\theta) \begin{pmatrix} x\\y \end{pmatrix} = \begin{pmatrix} x\cos\theta-y\sin\theta\\ x\sin\theta+y\cos\theta \end{pmatrix} \quad\longleftrightarrow\quad e^{i\theta}(x+iy). \]

So a planar rotation can be understood in two equivalent ways:

From this viewpoint, the complex eigenvalues are exciting because they reveal a more natural coordinate system for the transformation. What looks like a coupled two-dimensional motion in real coordinates becomes a simple multiplication in complex coordinates.

Over \(\mathbb{C}\), corresponding eigenvectors do exist. For example, for \(\lambda_1=\cos\theta+i\sin\theta\), one may choose

\[ \mathbf{v}_1= \begin{pmatrix} 1\\ -i \end{pmatrix}, \qquad \text{and for } \lambda_2=\cos\theta-i\sin\theta, \quad \mathbf{v}_2= \begin{pmatrix} 1\\ i \end{pmatrix}. \]

Thus the eigendecomposition does not merely produce numbers: it uncovers a change in representation. A real rotation matrix becomes diagonal once the plane is viewed through the complex-number model, where rotation is simply multiplication by \(e^{i\theta}\).

Reflections: The Case \(\det(\mathbf{R})=-1\)

Dropping the determinant condition, the general \(2\times 2\) orthogonal matrix with \(\det=-1\) is

\[ \mathbf{R}(\theta)= \begin{pmatrix} \cos\theta & \sin\theta\\ \sin\theta & -\cos\theta \end{pmatrix}, \]

which represents a reflection rather than a rotation. Its characteristic equation

\[(\cos\theta-\lambda)(-\cos\theta-\lambda)-\sin^2\theta=\lambda^2-1=0\]

gives the real eigenvalues \(\lambda=\pm1\) directly — consistent with the general \(|\lambda|=1\) fact above, but now attained by real numbers rather than a complex-conjugate pair. The eigenvector for \(\lambda=1\) points along the line of reflection itself (vectors on that line are left unchanged), while the eigenvector for \(\lambda=-1\) is perpendicular to it (such vectors are exactly flipped to their negative). Solving \(\mathbf{R}(\theta)\mathbf{v}=\mathbf{v}\) with the half-angle identities \(\cos\theta-1=-2\sin^2\tfrac\theta2\) and \(\sin\theta=2\sin\tfrac\theta2\cos\tfrac\theta2\) shows this line has angle \(\theta/2\) to the \(x\)-axis, so \(\mathbf{v}=(\cos\tfrac\theta2,\sin\tfrac\theta2)\): every \(2\times 2\) orthogonal matrix with negative determinant is a pure reflection across a line through the origin, tilted at half the parameter angle \(\theta\). Since two reflections cancel, \(\mathbf{R}(\theta)^2=\mathbf{I}\), matching \(\lambda^2=1\) for both eigenvalues.

mirror line (θ/2) p Rp v (λ=−1) −v Rv=−v

Three-Dimensional Rotations and Euler's Rotation Theorem

The pattern extends to three dimensions, where a proper rotation (orthogonal, \(\det=+1\)) is described by an axis of rotation \(\hat{\mathbf{n}}\) and an angle \(\theta\). Choosing coordinates so the \(z\)-axis points along \(\hat{\mathbf{n}}\), the rotation takes the familiar block form

\[ \mathbf{R}= \begin{pmatrix} \cos\theta & -\sin\theta & 0\\ \sin\theta & \cos\theta & 0\\ 0 & 0 & 1 \end{pmatrix}, \]

whose trace is \(1+2\cos\theta\). A change of basis (a similarity transform \(\mathbf{PRP}^{-1}\)) never changes the trace, so \(\mathrm{Tr}(\mathbf{R})=1+2\cos\theta\) holds for any \(3\times 3\) rotation matrix, however it happens to be written in the standard basis — a direct way to read off the rotation angle from the matrix entries alone.

The block form also makes the eigenvalues transparent: the third row/column contributes the eigenvalue \(1\) with eigenvector \(\hat{\mathbf{n}}\) itself (a vector along the axis is left fixed by a rotation about that axis), while the leading \(2\times 2\) block contributes the same pair \(e^{\pm i\theta}\) already found above. So every \(3\times 3\) rotation matrix has eigenvalues \(1,\,e^{i\theta},\,e^{-i\theta}\).

That the eigenvalue \(1\) — and hence a fixed rotation axis — must always exist, for any \(3\times 3\) proper rotation matrix, whatever its origin, is known as Euler's rotation theorem. It follows from two facts already established here: every eigenvalue of an orthogonal matrix has \(|\lambda|=1\), and the characteristic polynomial of a real \(3\times3\) matrix has real coefficients of odd degree \(3\), so its non-real roots must occur in a single complex-conjugate pair, leaving exactly one real root \(\lambda_1\). Being real and of modulus \(1\), \(\lambda_1=\pm1\); since \(\det(\mathbf{R})=\lambda_1\cdot e^{i\theta}\cdot e^{-i\theta}=\lambda_1=+1\) for a proper rotation, only \(\lambda_1=1\) survives. Every rotation of \(\mathbb{R}^3\), no matter how it is generated, therefore fixes some line through the origin pointwise — its rotation axis.

Eigenvalue Decomposition (EVD)

If \(\mathbf{A}\in\mathbb{R}^{n\times n}\) has \(n\) linearly independent eigenvectors \(\mathbf{q}_i\) (for \(i=1,\ldots,n\)) and is therefore diagonalisable, then \(\mathbf{A}\) can be factorized as

\[\mathbf{A} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^{-1}\]

where \(\mathbf{Q}\) is a \(n\times n\) matrix whose ith column is the eigenvector \(\mathbf{q}_i\) of \(\mathbf{A}\) and \(\mathbf{\Lambda}\) is the diagonal matrix whose diagonal elements are the corresponding eigenvalues \({\Lambda}_{ii}=\lambda_i\). The decomposition can directly be derived from the initial statement about eigenvalues and eigenvectors:

\[\begin{array}{rrl} &\mathbf{A} \mathbf{q}_i &= \lambda_i \mathbf{q}_i \quad (i=1,\ldots,n)\\ \Leftrightarrow&\mathbf{A} \mathbf{Q} &= \mathbf{Q} \mathbf{\Lambda} \\ \Leftrightarrow&\mathbf{A} &= \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^{-1} . \end{array}\]

Spectral theorem (important special case): If \(\mathbf{A}\) is real symmetric, it is always diagonalisable with an orthogonal eigenvector matrix. Hence

\[\mathbf{A}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T,\qquad \mathbf{Q}^T\mathbf{Q}=\mathbf{I}.\]

Note, however, that not every real matrix is diagonalisable over \(\mathbb{R}\). For example, a nontrivial two-dimensional rotation matrix has no real eigenvectors and therefore only becomes diagonalisable when the field is extended to \(\mathbb{C}\). Even a defective matrix with too few eigenvectors to diagonalise at all, or a rectangular matrix for which "eigenvector" is not even defined, always has a singular value decomposition, the natural generalisation of this decomposition that exists for every matrix.

Matrix Powers and the Square Root of a Symmetric Matrix

The "Powers" property above defines \(\mathbf{A}^t\) for integer \(t\) through the ordinary matrix product, but the spectral decomposition makes it possible to raise a symmetric matrix to any real power. For \(\mathbf{A}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T\) and a real exponent \(\alpha\), define

\[\mathbf{A}^{\alpha}:=\mathbf{Q}\mathbf{\Lambda}^{\alpha}\mathbf{Q}^T,\qquad \mathbf{\Lambda}^{\alpha}=\mathrm{diag}(\lambda_1^{\alpha},\ldots,\lambda_n^{\alpha}).\]

For this to be a real matrix, every \(\lambda_i^{\alpha}\) must itself be real, so a non-integer \(\alpha\) requires \(\mathbf{A}\) to be positive semidefinite (\(\lambda_i\ge 0\)), and a negative \(\alpha\) additionally requires \(\mathbf{A}\) to be invertible (\(\lambda_i\neq 0\)).

The most important case is \(\alpha=\frac12\), the symmetric square root of a positive semidefinite matrix. It is again symmetric, since \(\mathbf{Q}\) is orthogonal,

\[(\mathbf{A}^{1/2})^T=(\mathbf{Q}\mathbf{\Lambda}^{1/2}\mathbf{Q}^T)^T=\mathbf{Q}\mathbf{\Lambda}^{1/2}\mathbf{Q}^T=\mathbf{A}^{1/2},\]

and squaring it recovers \(\mathbf{A}\) exactly, using \(\mathbf{Q}^T\mathbf{Q}=\mathbf{I}\):

\[\mathbf{A}^{1/2}\mathbf{A}^{1/2} =\mathbf{Q}\mathbf{\Lambda}^{1/2}\underbrace{\mathbf{Q}^T\mathbf{Q}}_{\mathbf{I}}\mathbf{\Lambda}^{1/2}\mathbf{Q}^T =\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T=\mathbf{A}.\]

For the symmetric example above, \(\mathbf{\Lambda}=\mathrm{diag}(6,3,3)\), so \(\mathbf{\Lambda}^{1/2}=\mathrm{diag}(\sqrt6,\sqrt3,\sqrt3)\). Because the repeated eigenvalue \(3\) contributes the same scalar \(\sqrt3\) in both of its dimensions, the result does not depend on which orthonormal basis \(\mathbf{q}_2,\mathbf{q}_3\) of that eigenspace was chosen, and working out \(\mathbf{Q}\mathbf{\Lambda}^{1/2}\mathbf{Q}^T\) gives the closed form

\[ \mathbf{A}^{1/2}= \begin{pmatrix} d&o&o\\ o&d&o\\ o&o&d \end{pmatrix}, \qquad d=\frac{\sqrt6+2\sqrt3}{3}\approx1.9712, \qquad o=\frac{\sqrt6-\sqrt3}{3}\approx0.2391, \]

and squaring this matrix numerically reproduces \(\mathbf{A}=\begin{pmatrix}4&1&1\\1&4&1\\1&1&4\end{pmatrix}\) exactly.

Gram factors

Any matrix \(\mathbf{F}\) with \(\mathbf{F}\mathbf{F}^T=\mathbf{\Sigma}\) is called a Gram factor of the positive semidefinite matrix \(\mathbf{\Sigma}\). The symmetric square root itself is one such factor, since \(\mathbf{F}=\mathbf{\Sigma}^{1/2}\) satisfies \(\mathbf{F}\mathbf{F}^T=\mathbf{\Sigma}^{1/2}\mathbf{\Sigma}^{1/2}=\mathbf{\Sigma}\) by the identity above. Gram factors are never unique: if \(\mathbf{T}\) is any orthogonal matrix (\(\mathbf{T}\mathbf{T}^T=\mathbf{I}\)), then \(\mathbf{F}^*=\mathbf{F}\mathbf{T}\) is also a Gram factor, since

\[\mathbf{F}^*(\mathbf{F}^*)^T=\mathbf{F}\mathbf{T}\mathbf{T}^T\mathbf{F}^T=\mathbf{F}\mathbf{F}^T=\mathbf{\Sigma}.\]

One particular choice, obtained by Gaussian elimination instead of an eigendecomposition, restricts \(\mathbf{F}\) to be lower triangular with a positive diagonal; this special Gram factor is called the Cholesky decomposition of \(\mathbf{\Sigma}\).

Gram factors give a direct way to manufacture a random vector with a prescribed covariance matrix. If \(\mathbf{x}\) has \(\mathrm{Cov}(\mathbf{x})=\mathbf{I}\) (uncorrelated, unit-variance components) and \(\mathbf{\Sigma}=\mathbf{F}\mathbf{F}^T\), then \(\mathbf{y}:=\mathbf{F}\mathbf{x}\) satisfies

\[\mathrm{Cov}(\mathbf{y})=\mathbf{F}\,\mathrm{Cov}(\mathbf{x})\,\mathbf{F}^T=\mathbf{F}\mathbf{I}\mathbf{F}^T=\mathbf{F}\mathbf{F}^T=\mathbf{\Sigma},\]

turning independent, unit-variance noise into correlated data with a prescribed covariance structure, such as the one behind a covariance error ellipse.

Running this construction in reverse gives a whitening transform: if \(\mathbf{x}\) has \(\mathrm{Cov}(\mathbf{x})=\mathbf{\Sigma}\) for an arbitrary positive definite \(\mathbf{\Sigma}\), then \(\mathbf{\Sigma}^{-1/2}\mathbf{x}\) has covariance

\[\mathrm{Cov}(\mathbf{\Sigma}^{-1/2}\mathbf{x}) =\mathbf{\Sigma}^{-1/2}\mathbf{\Sigma}\mathbf{\Sigma}^{-1/2} =\mathbf{\Sigma}^{-1/2}\mathbf{\Sigma}^{1/2}\mathbf{\Sigma}^{1/2}\mathbf{\Sigma}^{-1/2}=\mathbf{I},\]

since \(\mathbf{\Sigma}^{1/2}\) and \(\mathbf{\Sigma}^{-1/2}\) are built from the same eigenvectors \(\mathbf{Q}\) and therefore commute. More generally, to turn data with covariance \(\mathbf{A}\) into data with a different prescribed covariance \(\mathbf{B}\), first whiten with \(\mathbf{A}^{-1/2}\) and then recolor with \(\mathbf{B}^{1/2}\):

\[\mathbf{y}:=\mathbf{B}^{1/2}\mathbf{A}^{-1/2}\mathbf{x} \quad\Longrightarrow\quad \mathrm{Cov}(\mathbf{y})=\mathbf{B}^{1/2}\underbrace{\mathbf{A}^{-1/2}\mathbf{A}\mathbf{A}^{-1/2}}_{\mathbf{I}}\mathbf{B}^{1/2}=\mathbf{B}.\]

The Rayleigh Quotient

For a real symmetric matrix \(\mathbf{A}\), the decomposition above gives an orthonormal eigenbasis \(\mathbf{q}_1,\ldots,\mathbf{q}_n\) with eigenvalues \(\lambda_1\ge\cdots\ge\lambda_n\). For a nonzero vector \(\mathbf{x}\), define the Rayleigh quotient

\[R(\mathbf{x}):=\frac{\mathbf{x}^T\mathbf{A}\mathbf{x}}{\mathbf{x}^T\mathbf{x}}.\]

Writing \(\mathbf{x}=\mathbf{Q}\mathbf{y}\) in eigenbasis coordinates \(\mathbf{y}=\mathbf{Q}^T\mathbf{x}\), and using \(\mathbf{A}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T\) together with \(\mathbf{Q}^T\mathbf{Q}=\mathbf{I}\),

\[\mathbf{x}^T\mathbf{A}\mathbf{x}=\mathbf{y}^T\mathbf{\Lambda}\mathbf{y}=\sum_{i=1}^n\lambda_i y_i^2, \qquad \mathbf{x}^T\mathbf{x}=\mathbf{y}^T\mathbf{y}=\sum_{i=1}^n y_i^2,\]

so the Rayleigh quotient is a weighted average of the eigenvalues,

\[R(\mathbf{x})=\frac{\sum_{i=1}^n\lambda_i y_i^2}{\sum_{i=1}^n y_i^2}.\]

A weighted average never exceeds the largest weight nor falls below the smallest, which immediately gives Rayleigh's principle:

\[\lambda_n\le R(\mathbf{x})\le\lambda_1\quad\text{for all }\mathbf{x}\neq\mathbf{0},\]

with equality \(R(\mathbf{x})=\lambda_1\) exactly when \(\mathbf{x}\) is proportional to \(\mathbf{q}_1\) (all weight on \(y_1\)), and \(R(\mathbf{x})=\lambda_n\) exactly when \(\mathbf{x}\) is proportional to \(\mathbf{q}_n\). So the extreme eigenvalues are exactly the extreme values of the quotient,

\[\lambda_1=\max_{\mathbf{x}\neq\mathbf{0}}R(\mathbf{x}),\qquad \lambda_n=\min_{\mathbf{x}\neq\mathbf{0}}R(\mathbf{x}).\]

The symmetric example above illustrates this directly: for \(\mathbf{q}_1=\frac{1}{\sqrt3}(1,1,1)^T\), \(R(\mathbf{q}_1)=6\), the largest eigenvalue; for either \(\mathbf{q}_2\) or \(\mathbf{q}_3\) spanning the repeated eigenspace, \(R(\mathbf{q}_2)=R(\mathbf{q}_3)=3\), the smallest.

The Cayley–Hamilton Theorem

A remarkable fact ties a matrix back to its own characteristic polynomial. Write \(p_{\mathbf{A}}(\lambda)=\lambda^n+c_{n-1}\lambda^{n-1}+\cdots+c_1\lambda+c_0\). Substituting the matrix \(\mathbf{A}\) itself for \(\lambda\) (and \(\mathbf{I}\) for the constant term) always yields the zero matrix:

\[\mathbf{A}^n+c_{n-1}\mathbf{A}^{n-1}+\cdots+c_1\mathbf{A}+c_0\mathbf{I}=\mathbf{O}.\]

A tempting but invalid shortcut is to simply substitute \(\lambda=\mathbf{A}\) directly into the defining determinant, \(\det(\mathbf{A}-\mathbf{A}\mathbf{I})=\det(\mathbf{O})=0\) — this is not a proof, since the \(0\) on the right is a scalar (a determinant), while \(p_{\mathbf{A}}(\mathbf{A})\) in the actual theorem is an \(n\times n\) matrix; the two sides are not even the same kind of object, let alone equal. A genuine proof must show that every entry of the matrix \(p_{\mathbf{A}}(\mathbf{A})\) vanishes, not merely that one particular determinant does.

This is the Cayley–Hamilton theorem. It is easiest to see first for a diagonalizable matrix \(\mathbf{A}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^{-1}\): since every power satisfies \(\mathbf{A}^k=\mathbf{Q}\mathbf{\Lambda}^k\mathbf{Q}^{-1}\) with the same \(\mathbf{Q}\), the whole polynomial collapses to \(p_{\mathbf{A}}(\mathbf{A})=\mathbf{Q}\,p_{\mathbf{A}}(\mathbf{\Lambda})\,\mathbf{Q}^{-1}\). But \(p_{\mathbf{A}}(\mathbf{\Lambda})\) is diagonal with entries \(p_{\mathbf{A}}(\lambda_i)\), which are all zero because \(\lambda_1,\ldots,\lambda_n\) are precisely the roots of \(p_{\mathbf{A}}\). Hence \(p_{\mathbf{A}}(\mathbf{A})=\mathbf{Q}\mathbf{O}\mathbf{Q}^{-1}=\mathbf{O}\).

To extend this from diagonalizable matrices to every matrix, note that each entry of \(p_{\mathbf{A}}(\mathbf{A})\) is a polynomial expression in the entries of \(\mathbf{A}\), hence a continuous function of them. Diagonalizable matrices are dense in \(\mathbb{C}^{n\times n}\): any matrix with repeated eigenvalues can be perturbed by an arbitrarily small amount to separate them and make it diagonalizable. Since the identity \(p_{\mathbf{A}}(\mathbf{A})=\mathbf{O}\) holds on this dense set and both sides depend continuously on \(\mathbf{A}\), it holds everywhere, including for matrices that are not diagonalizable.

One direct use of the theorem is expressing \(\mathbf{A}^{-1}\) purely in terms of lower powers of \(\mathbf{A}\). Since \(c_0=p_{\mathbf{A}}(0)=(-1)^n\det(\mathbf{A})\) is nonzero exactly when \(\mathbf{A}\) is invertible, multiplying the theorem by \(\mathbf{A}^{-1}\) and rearranging gives

\[\mathbf{A}^{-1}=-\frac{1}{c_0}\left(\mathbf{A}^{n-1}+c_{n-1}\mathbf{A}^{n-2}+\cdots+c_1\mathbf{I}\right).\]

More generally, the theorem lets any power \(\mathbf{A}^k\) with \(k\ge n\) be rewritten as a combination of only \(\mathbf{I},\mathbf{A},\ldots,\mathbf{A}^{n-1}\). This is precisely why the earlier fact that \(f(\mathbf{A})\) has eigenvalues \(f(\lambda_i)\) for a polynomial \(f\) reduces to a finite, explicitly computable matrix expression no matter how high the degree of \(f\) is.

Determinants from Traces of Powers

Taking the trace of the Cayley–Hamilton identity gives another useful shortcut. For \(n=2\), where \(p_{\mathbf{A}}(\lambda)=\lambda^2-\mathrm{Tr}(\mathbf{A})\lambda+\det(\mathbf{A})\) (the concrete example above), the theorem reads \(\mathbf{A}^2-\mathrm{Tr}(\mathbf{A})\mathbf{A}+\det(\mathbf{A})\mathbf{I}=\mathbf{O}\). Taking the trace of both sides and using \(\mathrm{Tr}(\mathbf{I})=2\),

\[\mathrm{Tr}(\mathbf{A}^2)-\mathrm{Tr}(\mathbf{A})^2+2\det(\mathbf{A})=0 \quad\Longrightarrow\quad \det(\mathbf{A})=\tfrac12\bigl(\mathrm{Tr}(\mathbf{A})^2-\mathrm{Tr}(\mathbf{A}^2)\bigr),\]

a formula for the determinant of any \(2\times 2\) matrix purely from traces of its powers, without ever computing \(ad-bc\) directly. This is the smallest case of a general pattern: the coefficients of the characteristic polynomial can always be expressed in terms of the power sums \(t_k:=\mathrm{Tr}(\mathbf{A}^k)\) via the classical Newton's identities, which relate the elementary symmetric functions of the eigenvalues (the \(c_i\)) to their power sums (the \(t_k\)). For \(n=3\) this yields

\[\det(\mathbf{A})=\tfrac16\bigl(\mathrm{Tr}(\mathbf{A})^3-3\,\mathrm{Tr}(\mathbf{A})\,\mathrm{Tr}(\mathbf{A}^2) +2\,\mathrm{Tr}(\mathbf{A}^3)\bigr),\]

which the symmetric example \(\mathbf{A}=\left(\begin{smallmatrix}4&1&1\\1&4&1\\1&1&4\end{smallmatrix}\right)\) confirms directly: with eigenvalues \(6,3,3\), the power sums are \(\mathrm{Tr}(\mathbf{A})=12\), \(\mathrm{Tr}(\mathbf{A}^2)=6^2+3^2+3^2=54\), and \(\mathrm{Tr}(\mathbf{A}^3)=6^3+3^3+3^3=270\), giving \(\tfrac16(12^3-3\cdot12\cdot54+2\cdot270)=\tfrac16(1728-1944+540)=54=\det(\mathbf{A})\), matching \(6\cdot3\cdot3=54\) computed directly from the eigenvalues.

Power Iteration and Asymptotic Growth

Suppose \(\mathbf{A}\) has \(n\) linearly independent eigenvectors \(\mathbf{v}_1,\ldots,\mathbf{v}_n\) with eigenvalues ordered so that \(|\lambda_1|>|\lambda_2|\ge\cdots\ge|\lambda_n|\), and write an arbitrary starting vector in this eigenbasis,

\[\mathbf{x}=c_1\mathbf{v}_1+c_2\mathbf{v}_2+\cdots+c_n\mathbf{v}_n,\qquad c_1\neq0.\]

Applying \(\mathbf{A}\) repeatedly and using the powers property above,

\[\mathbf{A}^k\mathbf{x}=c_1\lambda_1^k\mathbf{v}_1+c_2\lambda_2^k\mathbf{v}_2+\cdots+c_n\lambda_n^k\mathbf{v}_n =\lambda_1^k\Bigl(c_1\mathbf{v}_1+c_2\Bigl(\frac{\lambda_2}{\lambda_1}\Bigr)^{\!k}\mathbf{v}_2+\cdots\Bigr).\]

Since \(|\lambda_2/\lambda_1|<1\), every term but the first vanishes as \(k\to\infty\), so

\[\frac{\mathbf{A}^k\mathbf{x}}{\lambda_1^k}\longrightarrow c_1\mathbf{v}_1.\]

The direction of \(\mathbf{A}^k\mathbf{x}\) converges to the eigenvector of the dominant eigenvalue, while its length grows (or shrinks) geometrically like \(|\lambda_1|^k\). This is the basis of the numerical power iteration method for finding the largest eigenvalue of a matrix.

Successive normalized directions \(\mathbf{A}^k\mathbf{x}/ \|\mathbf{A}^k\mathbf{x}\|\) for \(k=0,1,\ldots,6\), starting from an arbitrary \(\mathbf{x}\), spiral around the unit circle and settle onto the dominant eigenvector \(\mathbf{v}_1\) of the matrix used in the Fibonacci example below.

A classic illustration is the Fibonacci recursion \(F_{k+1}=F_k+F_{k-1}\), which can be written as the matrix iteration

\[ \begin{pmatrix}F_{k+1}\\F_k\end{pmatrix} = \underbrace{\begin{pmatrix}1&1\\1&0\end{pmatrix}}_{=:\mathbf{A}} \begin{pmatrix}F_k\\F_{k-1}\end{pmatrix}. \]

Its characteristic polynomial is \(\lambda^2-\lambda-1=0\), with roots

\[ \lambda_{1,2}=\frac{1\pm\sqrt5}{2}, \]

the golden ratio \(\varphi=\frac{1+\sqrt5}{2}\approx1.618\) and its conjugate \(\psi=\frac{1-\sqrt5}{2}\approx-0.618\), with \(|\psi|<1<\varphi\). Solving \((\mathbf{A}-\lambda\mathbf{I})\mathbf{v}=\mathbf{0}\) gives the eigenvectors \(\mathbf{v}_1=(\varphi,1)^T\) and \(\mathbf{v}_2=(\psi,1)^T\). Decomposing the initial values \(F_0=0\), \(F_1=1\) in this eigenbasis and reading off the first coordinate reproduces the closed-form expression known as Binet's formula,

\[F_k=\frac{\varphi^k-\psi^k}{\sqrt5}.\]

Because \(|\psi|<1\), the second term vanishes exponentially, so \(F_k\sim\varphi^k/\sqrt5\) and

\[\lim_{k\to\infty}\frac{F_{k+1}}{F_k}=\varphi,\]

exactly the dominant-eigenvalue growth rate predicted above.

Applications of Eigenvalues and Eigenvectors