Defines eigenvalues, eigenvectors, eigenspaces and the characteristic polynomial, distinguishes algebraic from geometric multiplicity, and gives the diagonalizability criterion together with trace and determinant as the sum and product of the eigenvalues.
Definition(Eigenvalue, eigenvector, eigenspace)AdaptedAdapted from Hassani (2013), §6.2, Defs 6.2.1 and 6.2.3, Prop 6.2.2, pp. 172-173; Riley, Hobson & Bence (2006), §8.13, Eqs (8.67)-(8.68), pp. 272-273, and §8.13.1 (the word 'degenerate'), p. 275
Let A be an n×n matrix with complex entries, n≥1, equivalently a linear operator on an n-dimensional complex vector space with a basis fixed. A number λ∈C is an eigenvalue of A if
Ax=λxfor some x∈Cn,x=0,
and any such x is an eigenvector of Aforλ. The zero vector is excluded; without that exclusion every λ would qualify. Any non-zero scalar multiple of an eigenvector for λ is again an eigenvector for λ, so a single eigenvector never comes with a preferred length; for that reason eigenvectors are usually normalised to x†x=1. (When λ has more than one independent eigenvector, as happens below, normalisation leaves far more freedom than a scalar: the choice is then a choice of basis in the eigenspace.)
The set
Eλ=ker(A−λI)={x∈Cn:Ax=λx}
is a subspace of Cn, called the eigenspace of λ; its non-zero vectors are exactly the eigenvectors for λ. Its dimension g(λ)=dimEλ≥1 is the geometric multiplicity of λ, and λ is simple when g(λ)=1 and degenerate when g(λ)≥2. The set of all eigenvalues is the spectrumσ(A).
Ax=λx has a non-zero solution exactly when A−λI fails to be invertible. Since that happens precisely when its determinant vanishes, σ(A) is the zero set of the characteristic polynomial
χA(λ)=det(A−λI),
a polynomial of degree exactly n in λ with leading term (−1)nλn. The equation χA(λ)=0 is the characteristic (or secular) equation. Over C, if λ1,…,λp are the distinct roots then
χA(λ)=j=1∏p(λj−λ)mj,j=1∑pmj=n,
and the exponent mj=m(λj)≥1 is the algebraic multiplicity of λj.
A is diagonalizable if Cn has a basis consisting of eigenvectors of A; collecting such a basis v1,…,vn as the columns of an invertible S gives S−1AS=diag(μ1,…,μn), where μk is the eigenvalue of vk. The μk are the eigenvalues listed with repetitions, each λj occurring g(λj) times; they are not the distinct roots λ1,…,λp above unless p=n. A matrix that is not diagonalizable is called defective.
Let A be an n×n complex matrix, n≥1, with distinct eigenvalues λ1,…,λp of algebraic multiplicities m1,…,mp. Then:
A has at least one eigenvalue, and ∑jmj=n.
Eigenvectors belonging to distinct eigenvalues are linearly independent, and Eλi∩Eλj={0} for i=j.
detA=j=1∏pλjmj and trA=j=1∑pmjλj: the determinant is the product and the trace the sum of the n eigenvalues counted with algebraic multiplicity.
A is invertible if and only if 0∈/σ(A).
Statement 1 uses that a non-constant polynomial over C has a root; over R it can fail. Part 3 is asserted for every complex A, which is the generality of Hassani's Eq. (6.9) for the determinant and of Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Eq. (8.87), p. 280) for the trace; Arfken's §6.4 gives the pair only for a Hermitian matrix, as invariants of a unitary transformation.
Idea of the proofOriginal
The characteristic polynomial carries three of the four claims. Its degree is exactly n because, in the sum over permutations that defines the determinant, only the identity permutation keeps all n diagonal entries and so reaches that degree; the fundamental theorem of algebra then hands over a root, and a root is an eigenvalue because a vanishing determinant means a non-trivial kernel. The determinant and the trace are two of the polynomial's coefficients, each computable in two ways — from the matrix entries, and from the factorization into roots — and equating the two readings is exactly the assertion that they are the product and the sum of the eigenvalues. Independence of eigenvectors for distinct eigenvalues is a separate and shorter argument: applying A−λpI to a relation annihilates its last term and leaves a shorter relation.
ProofAdaptedAdapted from Hassani (2013), §6.2, Thm 6.2.5, Prop 6.2.4, Prop 6.2.6 and eqs. (6.7)-(6.9), pp. 173-174; Riley, Hobson & Bence (2006), §8.14, the coefficient of lambda^(N-1) and Eq. (8.87), p. 280
Part 1. By the permutation formula for the determinant (Determinants),
χA(λ)=det(A−λI)=σ∑sgn(σ)i=1∏n(A−λI)iσ(i).
Each factor is a polynomial in λ of degree at most 1, and a factor has degree 1 only when σ(i)=i. A permutation other than the identity moves at least two indices, so its term has degree at most n−2; the identity contributes ∏i(Aii−λ), of degree n with leading term (−1)nλn. Hence χA has degree exactly n with leading coefficient (−1)n, and is non-constant because n≥1. The fundamental theorem of algebra gives a root λ1∈C; then det(A−λ1I)=0, so A−λ1I is not invertible and its kernel contains some x=0 with Ax=λ1x. So A has an eigenvalue. Its eigenvalues are exactly the roots of χA, by the same equivalence read in both directions. Factoring χA completely over C and collecting equal roots gives χA(λ)=c∏j(λj−λ)mj; comparing degrees gives ∑jmj=n, and comparing leading coefficients gives c(−1)n=(−1)n, so c=1.
Part 2. Let x1,…,xp be eigenvectors for the distinct eigenvalues λ1,…,λp and suppose ∑j=1pcjxj=0. We argue by induction on p. For p=1, c1x1=0 with x1=0 gives c1=0. For p≥2, apply A−λpI to the relation; since (A−λpI)xj=(λj−λp)xj and the xp term dies,
j=1∑p−1cj(λj−λp)xj=0.
By the inductive hypothesis cj(λj−λp)=0 for j<p, and λj=λp, so cj=0 for j<p. The original relation then reads cpxp=0, giving cp=0 as well. For the intersection: if v∈Eλi∩Eλj with i=j then λiv=Av=λjv, so (λi−λj)v=0 and v=0.
Part 3. For the determinant, evaluate both expressions for χA at λ=0: the definition gives χA(0)=detA, and the factorization gives χA(0)=∏jλjmj.
For the trace, compare the coefficients of λn−1. On the matrix side, only the identity permutation contributes in that degree, as observed in part 1, and
i=1∏n(Aii−λ)=(−1)nλn+(−1)n−1(i=1∑nAii)λn−1+⋯.
On the factorized side, χA(λ)=∏j(λj−λ)mj=(−1)n∏j(λ−λj)mj, a product of n linear factors λ−λj in which λj appears mj times; expanding it, the coefficient of λn−1 is minus the sum of those n roots, so the coefficient of λn−1 in χA is (−1)n(−∑jmjλj)=(−1)n−1∑jmjλj. Equating the two gives trA=∑jmjλj.
Part 4. By part 3, detA=∏jλjmj, which is non-zero exactly when no λj equals 0; and A is invertible exactly when detA=0 (Determinants). □
Definition(The two multiplicities and diagonalizability)Original
With the notation above:
1≤g(λj)≤m(λj) for every j.
A is diagonalizable if and only if g(λj)=m(λj) for every j, equivalently if and only if ∑jg(λj)=n.
If A has n distinct eigenvalues, then every mj=1, so A is diagonalizable. This is sufficient but not necessary: the identity matrix has one eigenvalue and is diagonal.
Similar matrices have the same characteristic polynomial, hence the same eigenvalues with the same algebraic multiplicities, the same trace and the same determinant; corresponding eigenspaces have the same dimension. If B=S−1AS and Ax=λx, then B(S−1x)=λ(S−1x).
Idea of the proofOriginal
Take part 4 first, since the rest leans on it: a similarity leaves det(A−λI) alone, because the two copies of S cancel inside the determinant. Then the bound g≤m becomes visible by choosing coordinates well: put a basis of the eigenspace first, and in those coordinates the matrix begins with g columns that are λ times the first g standard basis vectors, which forces g factors of (λ−μ) out of the characteristic polynomial. The diagonalizability criterion is then a counting argument: the eigenspaces are independent of one another, so the only way to accumulate n independent eigenvectors is for each eigenspace to be as large as its algebraic multiplicity allows.
ProofOriginal
Part 4 comes first, because the other parts use it. Let B=S−1AS with S invertible. Then B−λI=S−1(A−λI)S, so similarity invariance of the determinant (Determinants) gives
χB(λ)=det(S−1(A−λI)S)=det(A−λI)=χA(λ)
for every λ. Equal characteristic polynomials have the same roots with the same exponents, which is the same spectrum with the same algebraic multiplicities; the trace and determinant then agree by part 3 of the previous statement. For the eigenspaces, the map x↦S−1x is linear and injective, and it carries Eλ(A) into Eλ(B): if Ax=λx then
B(S−1x)=S−1ASS−1x=S−1Ax=λ(S−1x).
The same computation with A and B exchanged, and S−1 in place of S, shows y↦Sy carries Eλ(B) into Eλ(A). The two maps are mutually inverse, so Eλ(A) and Eλ(B) are isomorphic and have the same dimension.
Part 1.g(λj)≥1 because an eigenspace contains an eigenvector, which is non-zero by definition.
For the upper bound, write λ=λj and g=g(λ), choose a basis v1,…,vg of Eλ, and extend it to a basis v1,…,vn of Cn (Vector spaces and bases). Let S be the invertible matrix with these columns, so Sei=vi, and put B=S−1AS. For i≤g the ith column of B is
Bei=S−1Avi=S−1(λvi)=λei,
so the first g columns of B−μI are (λ−μ)e1,…,(λ−μ)eg. Expand det(B−μI) along its first column by cofactors (Determinants): only the (1,1) entry is non-zero there, so det(B−μI)=(λ−μ)detM, where M is B−μI with its first row and column deleted. That deletion turns the old columns 2,…,g into (λ−μ)e1,…,(λ−μ)eg−1, so M has the same shape with g−1 in place of g. Iterating the expansion g times,
χB(μ)=(λ−μ)gq(μ)
for some polynomial q. By part 4, χA=χB, so (λ−μ)g divides χA(μ), and therefore the exponent m(λ) in the factorization of χA is at least g.
Part 2. Suppose A is diagonalizable and fix a basis of Cn of eigenvectors. Sort its members by eigenvalue, say kj of them lie in Eλj; then ∑jkj=n, and those kj members are linearly independent elements of Eλj, so kj≤g(λj). With part 1 and ∑jmj=n from the previous statement,
n=j∑kj≤j∑g(λj)≤j∑m(λj)=n,
so both inequalities are equalities. Since g(λj)≤m(λj) term by term, equality of the sums forces g(λj)=m(λj) for every j.
Conversely, suppose ∑jg(λj)=n; by part 1 and ∑jmj=n this is the same condition as g(λj)=m(λj) for all j. Choose a basis of each Eλj and concatenate them into a list of n vectors. The list is linearly independent: group a vanishing linear combination by eigenvalue as ∑juj=0 with uj∈Eλj. If some uj were non-zero, the non-zero ones would be eigenvectors for distinct eigenvalues obeying a non-trivial relation, which part 2 of the previous statement forbids. So every uj=0, and then every coefficient vanishes because each chosen basis is independent. A list of n independent vectors in Cn is a basis, and this one consists of eigenvectors, so A is diagonalizable.
Part 3. If A has n distinct eigenvalues then p=n, and mj≥1 with ∑j=1nmj=n forces mj=1 for each j. By part 1, 1≤g(λj)≤mj=1, so the two multiplicities agree and part 2 makes A diagonalizable. The converse fails: In has the single eigenvalue 1 with m=g=n, and is already diagonal. □
NoteOriginal
An eigenvector is a direction that A only stretches, by the factor λ. If there are enough such directions to make a basis, A acts as n independent one-dimensional stretches and every question about A — powers, inverses, exponentials, long-time behaviour of x˙=Ax — becomes n scalar questions. The characteristic polynomial is how you find the stretch factors; the gap between algebraic and geometric multiplicity is what can stop the basis from existing.
Remark(Where the hypotheses matter)AdaptedAdapted from Hassani (2013), §6.2, the rotation R_z(theta), p. 172, and Example 6.2.8, p. 174; Arfken, Weber & Harris (2013), §6.5, 'Defective Matrices', p. 324
The field matters. A rotation of the plane through an angle that is not a multiple of π moves every non-zero vector off its own line, as Hassani (Hassani (2013), §6.2, p. 172) observes of Rz(θ); concretely, (0−110) has χ(λ)=λ2+1 and so no real eigenvalue and no real eigenvector, although as a complex matrix it has the eigenvalues ±i. All statements above are for C; a real matrix is handled by regarding it as complex, and its non-real eigenvalues then come in conjugate pairs.
Algebraic multiplicity does not bound geometric multiplicity from below. For (1011) the only eigenvalue is λ=1 with m=2, but (A−I)x=0 forces the second component of x to vanish, so g=1. The matrix is defective: it has no eigenbasis, and the shortfall m−g is what the Jordan form of a defective matrix records.
Conventions. Many texts use det(λI−A), which is monic, instead of det(A−λI); the two differ by (−1)n and have the same roots. "Degenerate" (physics) and "repeated" or "multiple" (mathematics) both signal that an eigenvalue occurs more than once, but they count different things: "repeated" refers to algebraic multiplicity, m(λ)≥2, and "degenerate" to geometric multiplicity, g(λ)≥2. On a defective matrix the two conditions come apart, so the words are not interchangeable there; the second example below is repeated but not degenerate. The spectrum is sometimes defined as a set and sometimes as a list with multiplicities; we take it as a set.
Example(A diagonalizable 3 by 3 matrix, worked through)Original
Let
A=421151−1−22.
The characteristic polynomial. Expanding det(A−λI) along the first row,
χA(λ)=(4−λ)[(5−λ)(2−λ)+2]−[2(2−λ)+2]−[2−(5−λ)].
Since (5−λ)(2−λ)+2=λ2−7λ+12, the three brackets give −λ3+11λ2−40λ+48, then 2λ−6, then 3−λ, so
χA(λ)=−λ3+11λ2−39λ+45=−(λ−5)(λ−3)2.
The factorization is checked by expanding: (λ−5)(λ2−6λ+9)=λ3−11λ2+39λ−45. So σ(A)={5,3} with m(5)=1 and m(3)=2.
First check.trA=4+5+2=11 and 5+3+3=11; detA=χA(0)=45 and 5⋅3⋅3=45. Computing detA from the matrix independently, 4(10+2)−1(4+2)−1(2−5)=48−6+3=45.
The eigenspaces. For λ=5,
A−5I=−121101−1−2−3,
whose second row gives x1=x3 and whose first row then gives x2=x1+x3=2x1; the third row is satisfied automatically, since x1+2x1−3x1=0. So E5=span{(1,2,1)} and g(5)=1=m(5).
For λ=3,
A−3I=121121−1−2−1,
all of whose rows are multiples of (1,1,−1), so the single condition is x1+x2−x3=0. Hence E3 is a plane, g(3)=2=m(3), with basis (−1,1,0) and (1,0,1).
Diagonalization. Since ∑jg(λj)=1+2=3, part 2 of the second statement says A is diagonalizable, and the three eigenvectors found are a basis. Collecting them as columns,
S=121−110101,detS=1(1)−(−1)(2)+1(−1)=2=0,
and S−1AS=diag(5,3,3). Rather than invert S, verify the equivalent form AS=Sdiag(5,3,3) column by column:
A121=5105,A−110=−330,A101=303,
which are 5, 3 and 3 times the respective columns.
Example(One entry away: a defective matrix)Original
Change the (2,2) entry of the previous matrix from 5 to 6:
B=421161−1−22.
The same expansion along the first row changes in two places: the first bracket becomes (6−λ)(2−λ)+2=λ2−8λ+14, and the third cofactor det(216−λ1)=λ−4 contributes 4−λ where it contributed 3−λ before. The middle bracket, 2(2−λ)+2, is unchanged. Hence
χB(λ)=−λ3+12λ2−45λ+54=−(λ−6)(λ−3)2,
so m(6)=1 and m(3)=2, and the two checks hold again: trB=12=6+3+3 and detB=54=6⋅9. The spectrum looks much like before. The eigenspaces do not.
For λ=6, the rows of B−6I=−221101−1−2−4 give x1=x3 from the second and then x2=2x1+x3=3x1 from the first, so E6=span{(1,3,1)} and g(6)=1.
For λ=3,
B−3I=121131−1−2−1.
Rows one and three are identical and give x1+x2−x3=0; subtracting twice that from the second row leaves x2=0, and then x1=x3. So E3=span{(1,0,1)} is only a line: g(3)=1<2=m(3).
Hence ∑jg(λj)=1+1=2<3, and by part 2 of the second statement B has no eigenbasis: it is defective. The two multiplicities differ by exactly one, which is the amount of Jordan structure B carries (Jordan normal form). Note how little the spectrum revealed: the trace and determinant tests passed, the eigenvalues were real and simple to find, and nothing but computing dimE3 exposed the defect.
ApplicationOriginal
Decoupling a linear system. The reason to want an eigenbasis is that it reduces one n-dimensional problem to n one-dimensional ones. For x˙=Ax with A diagonalizable, substituting x=Sy turns the system into y˙=diag(λk)y, whose solutions are yk(t)=yk(0)eλkt; the general solution is then a superposition x(t)=∑kckeλktvk over the eigenvectors. Every qualitative question becomes a question about where the λk sit: the equilibrium is stable when all Reλk<0, the slowest decay rate is maxkReλk, and a complex conjugate pair α±iβ is an oscillation of frequency β with envelope eαt. The same substitution in xt+1=Axt turns a discrete model — a population projection matrix, a Markov chain, a numerical iteration — into independent geometric sequences, and the eigenvalue of largest modulus governs the long-run behaviour.
Where the defect shows up physically. The gap m−g of the second example is not a curiosity. A damped oscillator u¨+2γu˙+ω02u=0 written as a first-order system has a defective coefficient matrix exactly at critical damping γ=ω0, where the two roots coincide but only one eigenvector survives; the solution there is (c1+c2t)e−γt, and the factor of t is what stands in for the eigenvector that is missing. The same thing produces secular growth in perturbation theory at resonance and in orbital mechanics: a repeated root with a shortfall in eigenvectors always contributes polynomial-times-exponential behaviour, which is why engineers avoid designing near such a coincidence.
Observables and degeneracy. In quantum mechanics the measurable values of a quantity are the eigenvalues of its operator and the states with a definite value are the eigenvectors, so the spectrum is literally the set of possible outcomes of a measurement. The geometric multiplicity is what physics calls the degree of degeneracy — the 2ℓ+1 states of a given angular momentum, the two spin states of an electron — and it is almost always forced by a symmetry of the Hamiltonian. Splitting a degeneracy by an applied field (the Zeeman and Stark effects) is the statement that a perturbation has reduced some g(λ) by breaking the symmetry that protected it. For such operators the eigenbasis is guaranteed and can even be chosen orthonormal; that is the subject of the spectral theorem for Hermitian matrices.
A Hermitian matrix has real eigenvalues, orthogonal eigenvectors for distinct eigenvalues, and an orthonormal eigenbasis, so some unitary matrix of eigenvector columns diagonalizes it; a real symmetric matrix can be diagonalized by a real orthogonal matrix.
For a complex m×n matrix A, the Hermitian conjugate (adjoint) A† is the conjugate transpose, (A†)ij=Aji. On Cn we use the inner product
(x,y)=i=1∑nxiyi=y†x,
which is linear in the first slot and conjugate-linear in the second, with ∥x∥2=(x,x).
A square matrix A is Hermitian if A=A†, equivalently if (Ax,y)=(x,Ay) for all x,y∈Cn. For a matrix with real entries this is the condition A=AT, i.e. that A is symmetric. A matrix U is unitary if U†U=I, and a real matrix Q is orthogonal if QTQ=I; see unitary and orthogonal matrices.
Orthogonality. If Ax=λx, Ay=μy and λ=μ, then (x,y)=0.
Orthonormal eigenbasis. There is an orthonormal basis u1,…,un of Cn with Auk=λkuk, the λk real and listed with repetitions. Within a degenerate eigenspace the basis is not forced on you: any linearly independent set spanning the eigenspace can be orthonormalised by the Gram-Schmidt process, and the result is still an eigenbasis of that eigenspace.
Unitary diagonalization. Equivalently, the matrix U=(u1⋯un) whose columns are such a basis is unitary and
U†AU=Λ=diag(λ1,…,λn),soA=UΛU†.
Spectral resolution. Writing this out columnwise, A=k=1∑nλkukuk†, and trA=∑kλk, detA=∏kλk.
Real symmetric case. If in addition A is real, then A=AT, the eigenvalues are real, and the uk may be chosen real, so that S=(u1⋯un) is a real orthogonal matrix with STAS=Λ.
Idea of the proofOriginal
Parts 1 and 2 are one manoeuvre done twice: move A from one side of the inner product to the other, which is free because A=A†, and compare the two readings of the same number. A conjugate appears on one side and not on the other, and that mismatch is what makes the eigenvalue real and makes eigenvectors belonging to different eigenvalues perpendicular. Part 3 is where the work is. It is an induction that peels off one direction at a time: over C there is always at least one eigenvalue, and the hyperplane perpendicular to its eigenvector is again carried into itself by A — this is the third use of A=A† — so the same argument runs in one dimension lower. Parts 4 and 5 then only restate part 3: an orthonormal basis written as the columns of a matrix is a unitary matrix, and A=UΛU† read entry by entry is the spectral resolution.
ProofAdaptedAdapted from Arfken, Weber & Harris (2013), §6.3, Eqs. (6.17)-(6.21), pp. 310-311, and §6.4, Eqs. (6.22)-(6.24), pp. 311-312; Riley, Hobson & Bence (2006), §8.13.2, Eqs. (8.81)-(8.84), pp. 276-277; Hassani (2013), §6.4, Thm 6.4.3, p. 177, Cor 6.4.9, p. 181, and Cor 6.4.13, p. 182
Throughout we use hermiticity in the form (Ax,y)=(x,Ay) recorded in the definition, and conjugate-linearity of (⋅,⋅) in its second slot.
Part 1. Let Ax=λx with x=0. Then
λ(x,x)=(λx,x)=(Ax,x)=(x,Ax)=(x,λx)=λ(x,x),
the second-to-last step using A=A† and the last conjugate-linearity in the second slot. Since x=0 gives (x,x)=∥x∥2>0, we may cancel it and conclude λ=λ, that is λ∈R.
Part 2. Let Ax=λx, Ay=μy with λ=μ; both are real by part 1. Then
λ(x,y)=(Ax,y)=(x,Ay)=(x,μy)=μ(x,y)=μ(x,y),
where the last equality uses that μ is real. Hence (λ−μ)(x,y)=0, and λ=μ forces (x,y)=0.
Part 3. We prove, by induction on d, the following statement: if W is a complex inner product space with dimCW=d≥1 and T:W→W is linear with (Tx,y)=(x,Ty) for all x,y∈W, then W has an orthonormal basis of eigenvectors of T, and the corresponding eigenvalues are real. Taking W=Cn and T=A gives part 3.
For d=1, any unit vector u of W spans it, so Tu is a multiple of u and (u) is such a basis; the eigenvalue is real by the computation of part 1, which used nothing but (Tx,y)=(x,Ty).
Let d≥2 and assume the statement in dimension d−1. Fixing any basis of W represents T by a d×d complex matrix, which has an eigenvalue by part 1 of Eigenvalues and eigenvectors; so there are λ1∈C and u1∈W with u1=0 and Tu1=λ1u1, and we may scale u1 to ∥u1∥=1. By part 1, λ1 is real. Put
W1={w∈W:(w,u1)=0}.
W1 has dimension d−1. Extend (u1) to an orthonormal basis (u1,f2,…,fd) of W, which is possible by part 3 of the Gram-Schmidt statement in Orthonormal bases and Gram-Schmidt. By the expansion statement there, w=(w,u1)u1+∑i≥2(w,fi)fi for every w∈W, so (w,u1)=0 holds exactly when w∈span{f2,…,fd}. Hence W1=span{f2,…,fd}, an inner product space of dimension d−1 under the restricted inner product.
T maps W1 into itself. For w∈W1,
(Tw,u1)=(w,Tu1)=(w,λ1u1)=λ1(w,u1)=0,
the first equality being the hypothesis on T. So Tw∈W1.
The restriction of T to W1 is linear and still satisfies (Tx,y)=(x,Ty) for x,y∈W1, so the induction hypothesis supplies an orthonormal basis u2,…,ud of W1 consisting of eigenvectors of T with real eigenvalues. Each ui with i≥2 lies in W1 and is therefore orthogonal to u1, so (u1,…,ud) is orthonormal; it has d members in a space of dimension d, hence is a basis by part 1 of the expansion statement. This completes the induction.
The freedom inside a degenerate eigenspace. Let λ be an eigenvalue of A and let (w1,…,wd) be any linearly independent list spanning Eλ. Gram-Schmidt applied inside the inner product space Eλ returns an orthonormal basis of Eλ, and each output is a linear combination of w1,…,wd, hence still an element of the subspace Eλ, hence still an eigenvector for λ. Doing this in each eigenspace and taking the union gives an orthonormal set, because vectors drawn from eigenspaces of different eigenvalues are orthogonal by part 2 and those drawn from one eigenspace are orthonormal by construction.
It remains to see that this set is a basis of Cn, that is that ∑λdimEλ=n, the sum running over the distinct eigenvalues. One inequality comes from the eigenbasis u1,…,un just constructed: sorting its members by eigenvalue, the kλ of them lying in Eλ are linearly independent elements of Eλ, so kλ≤dimEλ, and summing over λ gives n=∑λkλ≤∑λdimEλ. The other comes from the union itself: being orthonormal it is linearly independent, by part 1 of the expansion statement in Orthonormal bases and Gram-Schmidt, and it has ∑λdimEλ members in Cn, so ∑λdimEλ≤n. Hence the two are equal to n, and an orthonormal list of n vectors in Cn is an orthonormal basis.
Part 4. Let U=(u1⋯un) for such a basis. Then
(U†U)ij=l=1∑nUliUlj=(uj,ui)=δij,
so U†U=I: U is unitary and U−1=U†. Column j of AU is Auj=λjuj, which is also column j of UΛ; so AU=UΛ, and multiplying on the left by U† gives U†AU=Λ, on the right by U† gives A=UΛU†.
For the last two identities, Λ=U−1AU is similar to A, so it has the same trace and the same determinant (Eigenvalues and eigenvectors, second statement, part 4); and for a diagonal matrix those are ∑kλk and ∏kλk read off the diagonal.
Part 6. Let A be real with A=AT; then A=A†, so every eigenvalue of A is real by part 1. Give Rn the real inner product x⋅y=∑ixiyi, which is symmetric, so that A=AT reads (Aw)⋅z=w⋅(Az).
We prove by induction on d: if W⊆Rn is a subspace with dimRW=d≥1 and AW⊆W, then W has an orthonormal basis of real eigenvectors of A. Part 6 is the case W=Rn, which satisfies ARn⊆Rn because A is real.
One real eigenvector in W. The point of the induction is that the restriction of A to W must be made to look symmetric again, and for that the basis of W has to be orthonormal. So take an orthonormal basis (f1,…,fd) of W, which exists by part 3 of the Gram-Schmidt statement in Orthonormal bases and Gram-Schmidt applied to the real inner product space W (the statement is proved there over R as well as over C). Let B∈Md×d(R) be the matrix of A∣W in that basis, so Afj=∑iBijfi; pairing with fi and using orthonormality gives Bij=(Afj)⋅fi. Then
Bji=(Afi)⋅fj=fi⋅(Afj)=(Afj)⋅fi=Bij,
the middle equality being A=AT and the last the symmetry of the real inner product. So B is real symmetric, hence Hermitian as a complex matrix. Read as a complex matrix it has an eigenvalue (Eigenvalues and eigenvectors, first statement, part 1, which applies because d≥1), and that eigenvalue λ is real by part 1 above. Therefore B−λId is a real matrix with det(B−λId)=0, so its columns are linearly dependent over R (Determinants): there is a real c=0 in Rd with (B−λId)c=0. Setting u=∑jcjfj, a non-zero real vector of W,
Au=j∑cjAfj=i∑(j∑Bijcj)fi=λi∑cifi=λu.
Scale u to ∥u∥=1. If d=1 then (u) is the required basis and we are done.
Peeling it off. Let d≥2 and put W′={w∈W:w⋅u=0}. Expanding any w∈W in an orthonormal basis of W whose first member is u — which Gram-Schmidt supplies, by part 3 of the Gram-Schmidt statement applied inside W — shows dimRW′=d−1, exactly as in part 3. And AW′⊆W′: for w∈W′ we have Aw∈W because AW⊆W, and
(Aw)⋅u=w⋅(Au)=λ(w⋅u)=0.
By the induction hypothesis W′ has an orthonormal basis u2,…,ud of real eigenvectors of A; each lies in W′ and so is orthogonal to u, so (u,u2,…,ud) is an orthonormal list of d vectors in W, hence a basis of W, and all of its members are real eigenvectors of A. This closes the induction.
Applying it to W=Rn gives a real orthonormal basis u1,…,un of Rn of eigenvectors of A, with real eigenvalues λ1,…,λn. Then S=(u1⋯un) satisfies STS=I and STAS=Λ by the computation of part 4 with every conjugation removed. □
NoteOriginal
A Hermitian matrix is, after a rotation of the coordinate axes, just a list of real numbers: it stretches n mutually perpendicular directions by n real factors and does nothing else. This is why Hermitian operators model observables (real measured values, orthogonal states) and why a real symmetric matrix always has principal axes. Reach for it whenever a matrix built as a Gram matrix, a Hessian, an energy or an adjacency matrix turns out to equal its own adjoint.
Remark(Where the hypotheses are used, and conventions)Original
Real eigenvalues are not enough.(1011) has the single real eigenvalue 1 but a one-dimensional eigenspace, so no eigenbasis exists; it is not symmetric. Part 3 is a statement about A=A†, not about where the spectrum lies.
Diagonalizability is not enough; orthogonality is exactly enough.(1002) and (1012) have the same eigenvalues, and the second is diagonalizable, but only the first is symmetric and only the first has orthogonal eigenvectors. The converse of parts 3-4 holds in the sharp form: if U†AU is diagonal with real entries for some unitary U, then A=UΛU† satisfies A†=A. Dropping "real entries" gives exactly the normal matrices; see normal matrices.
No Hermitian matrix is defective. Part 3 supplies n independent eigenvectors, so geometric and algebraic multiplicity agree for every eigenvalue of a Hermitian matrix.
Conventions. Some books state part 4 as "UAU−1 is diagonal", with the columns of U−1, not of U, being the eigenvectors; that is the same statement with U replaced by U†. Mathematical texts write A∗ for A† and say "self-adjoint" where physics texts say "Hermitian". Our inner product is linear in the first argument; texts that take it linear in the second obtain (x,Ay)=(Ax,y) in the same form, because A is Hermitian, so no statement above changes.
Example(A real symmetric matrix with a degenerate eigenvalue)Original
Let
H=62−423−2−4−26=HT.
Eigenvalues. Expanding det(H−λI) along the first row, with (3−λ)(6−λ)−4=λ2−9λ+14,
which factors as −(λ−11)(λ−2)2. The eigenvalues 11,2,2 are real, as part 1 requires, and the two checks hold: trH=6+3+6=15=11+2+2 and detH=44=11⋅2⋅2.
The simple eigenvalue. The rows of H−11I=−52−42−8−2−4−2−5 reduce to x1=−x3 and 2x2=−x3, so E11=span{(2,1,−2)}. Since ∥(2,1,−2)∥=3, take
u1=31(2,1,−2),Hu1=31(22,11,−22)=11u1.
The degenerate eigenvalue. Every row of H−2I=42−421−2−4−24 is a multiple of (2,1,−2), so E2 is the plane 2x1+x2−2x3=0, of dimension 2: here g=m=2, as part 3 requires, since an eigenbasis of R3 has to take two of its three members from E2. A convenient spanning pair is
v1=(1,0,1),v2=(0,2,1),
and these are not orthogonal: v1⋅v2=1. Part 3 says Gram-Schmidt inside E2 fixes that. First ∥v1∥=2, so q1=21(1,0,1). Then v2⋅q1=1/2 and
v2−21q1=(0,2,1)−21(1,0,1)=21(−1,4,1),
of norm 2118=23, so q2=321(−1,4,1). Both are still in E2 — check 2(−1)+4−2(1)=0 — hence still eigenvectors for λ=2.
Orthogonality across eigenspaces is automatic by part 2, and visibly so: u1⋅q1=322−2=0 and u1⋅q2=92−2+4−2=0.
The diagonalizing matrix. Put
S=2/31/3−2/31/201/2−1/(32)4/(32)1/(32).
Its columns are orthonormal: each was normalised to length 1 above, the two pairings u1⋅q1 and u1⋅q2 were just computed, and q1⋅q2=0 holds by the Gram-Schmidt construction of q2 — which one can also read off directly, q1⋅q2=6−1+0+1=0. So STS=I3 and STHS=diag(11,2,2), which is part 6.
The spectral resolution, checked. Part 5 reads H=11u1u1T+2q1q1T+2q2q2T. Because the three projectors sum to I, this is the same as H=2I+9u1u1T, and
9u1u1T=9⋅9142−421−2−4−24=H−2I,
which is exactly the matrix H−2I written out above.
Example(A complex Hermitian matrix of size two)Original
Let
H=(11−i1+i0).
It is Hermitian: the diagonal is real and the off-diagonal entries are conjugates. Since (1+i)(1−i)=2,
det(H−λI)=(1−λ)(−λ)−2=λ2−λ−2=(λ−2)(λ+1),
so the eigenvalues are 2 and −1: real, and summing to trH=1 with product detH=−2.
For λ=2, the first row of H−2I gives −x1+(1+i)x2=0, so x=(1+i,1) spans E2; its squared norm is 2+1=3. For λ=−1, the first row of H+I gives 2x1+(1+i)x2=0, hence x2=−1+i2x1=−(1−i)x1, so x=(1,−1+i) spans E−1, again of squared norm 3. Normalising,
u1=31(1+i,1),u2=31(1,−1+i).
Part 2 is confirmed directly: (u1,u2)=31[(1+i)1+1⋅(−1+i)]=31[(1+i)+(−1−i)]=0. Note where the conjugate sits — without it the two vectors would not come out orthogonal.
Inside the bracket each slot is 2(u1u1†)ij−(u2u2†)ij with the common factor 31 pulled out; in the (2,1) slot that is 2(1−i)−(−1+i)=2−2i+1−i=3−3i, where the minus sign in front of u2u2† and the minus signs inside it combine, and in the (1,2) slot it is the conjugate 3+3i.
Equivalently U=31(1+i11−1+i) satisfies U†U=I2 and U†HU=diag(2,−1).
ApplicationOriginal
Observables and measurement. Quantum mechanics represents a measurable quantity by a Hermitian operator, and the theorem is what makes that representation coherent: part 1 says the possible measured values are real, part 2 says states with different values are distinguishable (orthogonal), and part 3 says every state can be expanded in states of definite value, with ∣(ψ,uk)∣2 the probability of the outcome λk. Part 5 is the resolution of the identity that textbook derivations insert between operators, and the geometric multiplicity of an eigenvalue is the degeneracy of that energy level. For 2×2 problems the matrix of the second example is the generic case: any Hermitian 2×2 matrix is aI+b⋅σ with a and b real, so its eigenvalues are a±∣b∣ and its eigenvectors are the spin states along b.
Principal axes. A real symmetric matrix met as a physical tensor acquires, through part 6, a frame in which it is diagonal and which is still a genuine rotation or reflection of space. For the inertia tensor these are the principal axes and the eigenvalues are the principal moments, which is why a rigid body has three special axes about which it can spin without wobble; for the stress tensor they are the principal stresses and the directions of zero shear; for the strain tensor they give the directions of pure extension; for the Hessian of a potential they give the directions of pure curvature, with the signs of the eigenvalues classifying the critical point. In each case the content of part 6 is that the diagonalizing change of coordinates is orthogonal, so lengths and angles — and therefore the physics — are unchanged.
Small oscillations and data. A system near a stable equilibrium has a symmetric stiffness matrix and a symmetric positive-definite mass matrix, and the simultaneous diagonalization of the two is what produces normal modes with definite frequencies, each oscillating independently (Generalized eigenproblems and normal modes). The same theorem underlies principal component analysis: a covariance matrix is symmetric by construction, its orthonormal eigenvectors are the principal directions of a cloud of data, and the eigenvalues are the variances along them, so part 5 decomposes the covariance into uncorrelated contributions ranked by size. Because the eigenvalues are real they can be ordered, which is what makes the Rayleigh quotient a usable bound on them (Quadratic forms and the Rayleigh quotient).
Defines unitary and orthogonal matrices, shows they are exactly the matrices that preserve the inner product and have orthonormal columns, and records that their eigenvalues have modulus one and that the real ones reduce to rotations and reflections.
Definition(Unitary and orthogonal matrices)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.12.4 and §8.12.6, pp. 270-272
Let U be an n×n complex matrix. U is unitary if it is invertible with
U†=U−1,equivalentlyU†U=I,
and for a square matrix U†U=I already implies UU†=I. A real matrix Q with QT=Q−1 is orthogonal; an orthogonal matrix is exactly a unitary matrix with real entries. The real orthogonal n×n matrices form the group O(n), those with detQ=+1 the subgroup SO(n), and the unitary n×n matrices the group U(n): in each case the identity belongs to the set, and products and inverses stay in it.
For an n×n complex matrix U the following are equivalent:
U†U=I;
the columns of U are orthonormal in Cn, i.e. (Uej,Uek)=δjk — and likewise for the rows;
(Ux,Uy)=(x,y) for all x,y∈Cn;
∥Ux∥=∥x∥ for all x∈Cn.
The same four statements are equivalent for a real matrix Q with QTQ=I, real vectors and the real inner product xTy. That 4 implies 3 is the polarization identity. Over R one identity suffices,
2xTy=∥x+y∥2−∥x∥2−∥y∥2,
and over C the real and imaginary parts are recovered separately,
the second being the first applied with y replaced by iy; each right-hand side is built only from norms, so hypothesis 4 — used at x, at y and at x+y and x+iy — reproduces (Ux,Uy)=(x,y) in full.
Idea of the proofOriginal
The single identity (Ux,Uy)=y†U†Ux turns the matrix condition U†U=I into a statement about inner products, and reading it at the basis vectors x=ej, y=ek turns it back again. So conditions 1, 2 and 3 are the same equation tested on all vectors, on basis vectors, and as a matrix. Condition 4 looks weaker, because it tests only x=y; the polarization identities say it is not, because an inner product on Cn is already determined by the norms of x, y, x+y and x+iy.
ProofOriginal
Throughout, (x,y)=y†x, which is linear in the first slot and conjugate-linear in the second, since (x,cy)=(cy)†x=c(x,y).
1 implies 3. For any x,y, the associativity of matrix multiplication gives
(Ux,Uy)=(Uy)†(Ux)=y†U†Ux=y†Ix=(x,y),
where hypothesis 1 was used at the third equality.
3 implies 2. Take x=ej and y=ek: then (Uej,Uek)=(ej,ek)=ek†ej=δjk, which is condition 2 for the columns, since Uej is the j-th column of U.
2 implies 1. The (k,j) entry of U†U is ek†U†Uej=(Uej,Uek), which is δkj by hypothesis 2. Every entry of U†U therefore agrees with the corresponding entry of I.
The rows. Condition 1 gives det(U†U)=detI=1, so detU=0 by the multiplicativity of the determinant (Determinants) and U is invertible with U−1=U†; substituting that identification into UU−1=I gives UU†=I. Its (j,k) entry is ∑mUjmUkm, which is exactly the inner product of the j-th row with the k-th row, so the rows are orthonormal as well. Conversely, orthonormal rows give UU†=I and hence, by the same argument with U and U† exchanged, U†U=I. This is the step that uses squareness: for a non-square matrix U†U=I and UU†=I are not equivalent, and only one of them can hold.
3 implies 4. Put y=x: then ∥Ux∥2=(Ux,Ux)=(x,x)=∥x∥2, and both norms are non-negative, so ∥Ux∥=∥x∥.
4 implies 3. Fix x,y. Hypothesis 4 applied at the four vectors x, y, x+y and x+iy gives
and U is linear, so U(x+y)=Ux+Uy and U(x+iy)=Ux+iUy. Substituting these four equalities into the two polarization identities displayed in the statement — once for the pair (Ux,Uy) and once for the pair (x,y) — makes their right-hand sides identical, so Re(Ux,Uy)=Re(x,y) and Im(Ux,Uy)=Im(x,y). A complex number is fixed by its real and imaginary parts, so (Ux,Uy)=(x,y).
The four conditions are therefore equivalent: 1, 2 and 3 are equivalent by the first three paragraphs, and 3 and 4 by the last two.
The real case. Replace † by T and Cn by Rn throughout. Every step above is unchanged except the last, where only the first polarization identity is available and is all that is needed, because a real inner product has no imaginary part to recover. □
Definition(Determinant and spectrum)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.12.4, p. 270, §8.12.6, p. 271 and §8.13.3, p. 278
Let U be unitary and Q real orthogonal, both n×n. Then:
∣detU∣=1, and detQ=±1.
Every eigenvalue λ of U satisfies ∣λ∣=1, so λ=eiθ for some real θ. In particular the only possible real eigenvalues are +1 and −1.
U−1=U† is unitary, and Q−1=QT is orthogonal.
Part 1 does not characterise unitarity: (1011) has determinant 1 and is not orthogonal.
Idea of the proofOriginal
Take determinants of U†U=I and the modulus of detU falls out, because conjugating a matrix conjugates its determinant. The eigenvalues are controlled by the characterisation rather than by any new computation: an eigenvector is stretched by ∣λ∣ and not stretched at all, so ∣λ∣=1. The third part is the observation that the defining relation is symmetric in U and U†.
Part 1. Taking determinants in U†U=I and using multiplicativity (Determinants),
1=detI=det(U†)detU.
Now U†=UT, and the determinant is unchanged by transposition and conjugated by entrywise conjugation, so det(U†)=detU. Hence 1=detUdetU=∣detU∣2, that is ∣detU∣=1. For a real orthogonal Q the same computation reads 1=det(QT)detQ=(detQ)2, and detQ is real, so detQ=±1.
Part 2. Let Ux=λx with x=0. Part 4 of the characterisation, applied to this x, gives ∥Ux∥=∥x∥; on the other hand ∥Ux∥=∥λx∥=∣λ∣∥x∥ by homogeneity of the norm. Since x=0 we may divide by ∥x∥>0 to get ∣λ∣=1, and every complex number of modulus one is eiθ with θ real. If moreover λ is real then λ=±1, those being the only real numbers of modulus one. (Riley, Hobson and Bence reach the same conclusion by sandwiching: x†x=x†U†Ux=λλx†x, so ∣λ∣2=1.)
Part 3. By part 1, detU=0, so U is invertible, and U†U=I identifies U−1=U†; substituting that into UU−1=I gives UU†=I. Writing V=U†, this says V†V=UU†=I, so V is unitary. The real case is the same computation with transposes. □
Definition(Rotations and reflections)AdaptedAdapted from Arfken, Weber & Harris (2013), §3.3, Eqs (3.23)-(3.28) and 'Reflections', pp. 134-136; Hassani (2013), Thm 6.6.12 and Eq. (6.26), p. 201
Plane rotations. The rotation of R2 about the origin by an angle θ has matrix
R(θ)=(cosθsinθ−sinθcosθ),detR(θ)=+1.
Every 2×2 real orthogonal matrix is either R(θ), with determinant +1, or
F(θ)=(cosθsinθsinθ−cosθ),detF(θ)=−1,
for some θ∈[0,2π). F(θ) fixes the line through the origin at angle θ/2 and negates the perpendicular direction, so it is the reflection in that line.
3. Real canonical form. For any real orthogonal Q there is a real orthogonal P with
a block-diagonal matrix with N++N−+2m=n. A real orthogonal matrix therefore cannot in general be diagonalized over R; this canonical form is as close as one can get.
Idea of the proofOriginal
In two dimensions orthogonality is almost no freedom at all. The first column is a unit vector of the plane, so it is (cosθ,sinθ) for one angle θ; the second must be a unit vector on the line perpendicular to it, and such a line contains exactly two unit vectors. One choice gives R(θ) and the other F(θ), which is the whole of part 2. That F(θ) is a reflection is then read off from the angle-addition formulas, which show at once that it fixes the direction at angle θ/2 and reverses the direction perpendicular to it.
Part 1. Write a point of R2 in polar form as (rcosφ,rsinφ). By definition the rotation by θ about the origin leaves r alone and sends φ to φ+θ, so it sends this point to (rcos(φ+θ),rsin(φ+θ)). Expanding with the angle-addition formulas,
rcos(φ+θ)=xcosθ−ysinθ,rsin(φ+θ)=xsinθ+ycosθ,
where x=rcosφ, y=rsinφ. Both are linear in (x,y), and reading off the coefficients gives exactly R(θ)(x,y)T. Its columns (cosθ,sinθ) and (−sinθ,cosθ) each have norm cos2θ+sin2θ=1 and inner product −cosθsinθ+sinθcosθ=0, so R(θ) is orthogonal by part 2 of the characterisation, and detR(θ)=cos2θ+sin2θ=+1.
Part 2, the two forms. Let Q be a 2×2 real orthogonal matrix with columns c1,c2. By the characterisation these are orthonormal. Since ∥c1∥=1, there is exactly one θ∈[0,2π) with c1=(cosθ,sinθ)T: this is the polar form of a unit vector of the plane. The vectors orthogonal to c1 in R2 form the solution set of the single non-trivial linear equation c1Tv=0, hence a line, and (−sinθ,cosθ)T is a unit vector on it; a line through the origin contains exactly two unit vectors, so
c2=±(−sinθ,cosθ)T.
With the upper sign Q=R(θ). With the lower sign c2=(sinθ,−cosθ)T and Q=F(θ), whose determinant is cosθ(−cosθ)−sinθsinθ=−1. These are the only possibilities, so every 2×2 real orthogonal matrix is one or the other, and the determinant says which.
Part 2, why F(θ) is a reflection. Set v=(cos(θ/2),sin(θ/2))T, the unit vector at angle θ/2, and w=(−sin(θ/2),cos(θ/2))T, the unit vector perpendicular to it. Multiplying out and using the angle-subtraction formulas with θ−θ/2=θ/2,
So F(θ) fixes every vector of the line Rv, which makes angle θ/2 with the first axis, and negates every vector of the perpendicular line Rw. Since {v,w} is a basis of R2, these two actions determine F(θ), and a linear map that fixes a line through the origin and negates its orthogonal complement is by definition the reflection in that line. □
Note(Part 3 is not proved here)AdaptedAdapted from Hassani (2013), §6.6, Thm 6.6.11, pp. 200-201, Thm 6.6.12 and Eq. (6.26), pp. 201-202
The real canonical form, part 3 of the statement above, is quoted without proof: it rests on the real spectral decomposition of a normal operator, which decomposes a real inner product space into one- and two-dimensional subspaces invariant under the operator, and that machinery is outside this unit. (Hassani (2013), Thm 6.6.11, pp. 200-201) proves the decomposition and then deduces the canonical form, with the orthogonality of the transforming matrix, in (Hassani (2013), Thm 6.6.12, pp. 201-202). Hassani also gives, in Example 6.6.13 on p. 202, a short non-rigorous route through the complexification: regard Q as unitary, write it as eA with A real antisymmetric, note that the eigenvalues of A are 0 and pure imaginary pairs ±iθj, and exponentiate. Our worked example below carries out the construction by hand in a single case, which is a check on the form but not a proof of it.
NoteOriginal
A unitary matrix is a rigid motion of Cn fixing the origin: it moves vectors around but changes no length and no angle. That is why unitary matrices are the allowed changes of orthonormal basis, the matrices that diagonalize Hermitian and normal matrices, the time evolutions of quantum mechanics, and — in the real case — the rotations and reflections of rigid-body kinematics.
Remark(Conventions and traps)Original
Active and passive. Expressing the same vector in axes rotated by +θ multiplies its component column by R(θ)T=R(−θ), whereas rotating the vector by +θ in fixed axes multiplies it by R(θ). Books that set up rotations by rotating the coordinate system therefore display the transpose of the matrix above. Always say which of the two you mean.
O(n) is not the same as SO(n). Orthogonality forces only detQ=±1; proper rotations are the detQ=+1 half. The two halves are not connected: the determinant is continuous and takes no value between −1 and +1 on O(n), so no continuous family of rotations starting at the identity ever reaches a reflection.
"Orthogonal matrix" is not "matrix with orthogonal columns". The columns must also be unit vectors. A matrix Q whose columns are pairwise orthogonal and all non-zero satisfies QTQ=D with D diagonal and every diagonal entry positive; such a Q stretches the j-th coordinate direction by Djj, so unless all the Djj are equal it changes lengths by different factors in different directions, and unless they all equal 1 it changes lengths at all. Dropping "non-zero" loses even this: the zero vector is orthogonal to everything, so a matrix with a zero column has pairwise orthogonal columns while being singular.
Preserving the norm requires all x.(0010) sends e2 to e1, so it preserves the length of e2, yet it kills e1; part 4 of the characterisation is a statement about every vector, not about one.
The real block form is not unique as displayed. Part 3 of "Rotations and reflections" leaves each θj free in [0,2π), and R(0)=I2 while R(π)=diag(−1,−1): a block R(0) is two +1 entries and a block R(π) two −1 entries, so the same Q admits displays with different N± and m. Uniqueness, up to the order of the blocks, needs θj∈(0,π), which forces sinθj=0; then the 2×2 blocks carry exactly the non-real eigenvalue pairs e±iθj and N± are the multiplicities of ±1. (Hassani (2013), Thm 6.6.12, p. 201) states the form without that restriction, and part 3 follows the source.
Real versus complex diagonalization. A unitary matrix is normal, hence unitarily diagonalizable over C. R(θ) with sinθ=0 illustrates the difference: its eigenvalues e±iθ are not real, so it is diagonalizable over C but not over R, which is what part 3 of the previous statement expresses.
Example(Classifying two orthogonal matrices)Original
Part A: a 2×2 matrix. Let
Q1=51(344−3).
Its columns are (3,4)/5 and (4,−3)/5, each of norm 9+16/5=1, with inner product (12−12)/25=0; so Q1 is orthogonal by part 2 of the characterisation. Its determinant is (−9−16)/25=−1, so by part 2 of "Rotations and reflections" it is F(θ), not a rotation, with cosθ=3/5 and sinθ=4/5. The mirror line is the one at angle θ/2, and
tan2θ=1+cosθsinθ=1+3/54/5=8/54/5=21,
so the mirror is the line y=x/2, spanned by (2,1). Checking directly: Q1(2,1)T=51(6+4,8−3)T=(2,1)T, fixed; and Q1(−1,2)T=51(−3+8,−4−6)T=(1,−2)T=−(−1,2)T, negated. The eigenvalues are +1 and −1, as part 2 of "Determinant and spectrum" allows for a real matrix.
Part B: a 3×3 matrix, reduced to the canonical form. Let
Q=3121−2−22−1122.
Step 1: it is orthogonal. The columns are c1=(2,1,−2)/3, c2=(−2,2,−1)/3, c3=(1,2,2)/3. Each has squared norm (4+1+4)/9=(4+4+1)/9=(1+4+4)/9=1, and
Step 2: norms are preserved. As a check on part 4 of the characterisation, take x=(1,2,3)T. Then Qx=31(2−4+3,1+4+6,−2−2+6)T=31(1,11,2)T, and
∥Qx∥2=91+121+4=9126=14=1+4+9=∥x∥2.
Step 3: the determinant. Expanding along the first row of the integer matrix 3Q,
det(3Q)=2(4+2)+2(2+4)+1(−1+4)=12+12+3=27,
so detQ=27/33=+1 and Q∈SO(3): a rotation, not a reflection.
Step 4: the spectrum. The trace is (2+2+2)/3=2, the sum of the three 2×2 principal minors of Q is 2, and detQ=1, so the characteristic polynomial is
λ3−2λ2+2λ−1=(λ−1)(λ2−λ+1),
with roots 1 and 21±23i=e±iπ/3. All three have modulus one, as part 2 of "Determinant and spectrum" requires. Comparing the trace with the canonical form diag(1,R(θ)), whose trace is 1+2cosθ, gives cosθ=21: the rotation angle is θ=π/3, or 60∘.
Step 5: the axis. Solving (Q−I)x=0, i.e. 31(−x1−2x2+x3,x1−x2+2x3,−2x1−x2−x3)=0, gives the one-dimensional solution space spanned by n=(−1,1,1)T; indeed Qn=31(−2−2+1,−1+2+2,2−1+2)T=31(−3,3,3)T=n. Normalising, n^=(−1,1,1)T/3.
Step 6: the canonical form. Pick any unit vector perpendicular to n^, say e^=(1,1,0)T/2 (note n^Te^=(−1+1+0)/6=0), and complete to a right-handed frame with
With P=(n^e^f^), which is orthogonal because its columns are orthonormal by construction and has detP=+1,
PTQP=10001/23/20−3/21/2=diag(1,R(π/3)),
the canonical form of part 3 with N+=1, N−=0, m=1 and θ1=π/3. Since sin(π/3)=0 this display is the unique one up to the order of the blocks, as the remark above requires. The upper-left 1 is the axis and the 2×2 block the rotation by 60∘ in the plane perpendicular to it.
ApplicationOriginal
Unitarity is conservation of probability. In quantum mechanics a state is a unit vector ψ of a complex inner product space, and ∣(ψ,ϕ)∣2 is the probability of finding the state ψ in the state ϕ; the total probability is ∥ψ∥2=1. A time evolution is therefore allowed only if it preserves the norm of every state — part 4 of the characterisation — which by that characterisation means exactly that it is unitary. For a finite-dimensional system with Hermitian Hamiltonian H, the evolution U(t)=e−iHt/ℏ satisfies U(t)†=e+iH†t/ℏ=eiHt/ℏ=U(t)−1, so it is unitary, and part 2 above says its eigenvalues are the phases e−iEkt/ℏ attached to the energy levels: time evolution rotates phases and changes no amplitude. The same argument constrains a scattering matrix S: unitarity of S is the statement that no probability is lost between the incoming and outgoing channels, and it is the hypothesis behind the optical theorem. Changing between two orthonormal bases is likewise a unitary matrix, which is why the predictions of the theory do not depend on the basis chosen.
Orthogonality turns a rigid rotation into an antisymmetric angular velocity. Let Q(t)∈SO(3) be the matrix carrying body-fixed axes to laboratory axes at time t, a differentiable function of t. Differentiating the constraint QTQ=I gives
Q˙TQ+QTQ˙=0,soΩ=QTQ˙ satisfies ΩT=−Ω.
A 3×3 antisymmetric matrix has only three independent entries, and those three numbers are the components of the angular velocity vector: a rigid body in three dimensions has three rotational degrees of freedom, not nine, and this is why. The determinant is what forbids the body from passing continuously into a mirror image of itself, as the remark on O(n) versus SO(n) explains: detQ(t) is continuous in t and confined to {±1}, so it cannot change during a motion.
Orthogonal transformations are the numerically safe ones. In numerical linear algebra the standard algorithms — QR factorisation, the reduction to Hessenberg or tridiagonal form, the QR eigenvalue iteration, the singular value decomposition — act by multiplying on the left by orthogonal matrices, built either as Givens plane rotations R(θ) embedded in a pair of coordinates or as Householder reflections, which are the n-dimensional analogues of F(θ). The reason is part 4 of the characterisation: because ∥Qx∥=∥x∥ exactly, an orthogonal step neither magnifies the vector nor magnifies the error already in it, so the rounding errors of the previous steps are carried along rather than amplified. A general invertible step, Gaussian elimination for instance, can multiply an error by a large factor. In the language of conditioning, an orthogonal matrix has condition number ∥Q∥∥Q−1∥=1 in the Euclidean norm, the smallest value possible, which is the quantitative form of the same statement.
A complex matrix has an orthonormal eigenbasis exactly when it is normal, and two normal matrices can be given a common orthonormal eigenbasis exactly when they commute.
An n×n complex matrix A is normal if it commutes with its Hermitian conjugate:
AA†=A†A.
Hermitian matrices (A=A†), anti-Hermitian matrices (A†=−A) and unitary matrices (A†=A−1) are normal, since in each case A† is a function of A that commutes with it; in the real case this covers symmetric, antisymmetric and orthogonal matrices. If A is normal and invertible then A−1 is normal.
Theorem(Spectral theorem for normal matrices)AdaptedAdapted from Hassani (2013), Thm 6.4.8 and Cor 6.4.9, pp. 180-181; Riley, Hobson & Bence (2006), §8.13.1, pp. 274-276
Let A be a normal n×n complex matrix, with distinct eigenvalues λ1,…,λp and eigenspaces Eλj. Then:
Cn=Eλ1⊕⋯⊕Eλp, the summands being mutually orthogonal, and Cn has an orthonormal basis u1,…,un of eigenvectors of A. Equivalently, the matrix U with these columns is unitary and U†AU=Λ=diag(λ1′,…,λn′), where the λk′ are the eigenvalues listed with multiplicity.
Spectral resolution.A=j=1∑pλjPj=k=1∑nλk′ukuk†, where Pj is the orthogonal projection onto Eλj, so that Pj=Pj†=Pj2, PiPj=0 for i=j, and ∑jPj=I.
The adjoint shares the eigenvectors. If Ax=λx then A†x=λx; the same x serves both matrices, and σ(A†)={λ:λ∈σ(A)}.
Eigenvectors belonging to distinct eigenvalues are orthogonal, and within one eigenspace any spanning independent set may be replaced by an orthonormal one by the Gram-Schmidt process.
Idea of the proofOriginal
A normal matrix is a complex combination of two commuting Hermitian ones: put H=21(A+A†) and K=2i1(A−A†), so that A=H+iK, and then normality of A is exactly the statement HK=KH. So the whole theorem reduces to one question about Hermitian matrices, which the spectral theorem already answers one matrix at a time: can two commuting Hermitian matrices be diagonalized in the same orthonormal basis? They can, because each eigenspace of H is carried into itself by K, so the Hermitian spectral theorem can be applied a second time inside each eigenspace, where H is already a scalar. A common orthonormal eigenbasis of H and K is then an eigenbasis of A=H+iK, with eigenvalues αk+iβk. The remaining parts — that A† has the same eigenvectors with conjugated eigenvalues, and that distinct eigenvalues give orthogonal eigenvectors — come from the single observation that ∥Bx∥=∥B†x∥ whenever B is normal, applied to B=A−λI.
A restriction lemma. Let W⊆Cn be a subspace of dimension d≥1 that is carried into itself by a matrix T and by T†. Choose an orthonormal basis w1,…,wd of W, which exists by Gram-Schmidt, and set Sjk=(Twk,wj). Since Twk lies in W, expanding it in that orthonormal basis gives Twk=∑jSjkwj: the matrix S represents T on W. The matrix representing T† on W in the same basis has entries
(T†wk,wj)=(wk,Twj)=(Twj,wk)=Skj=(S†)jk,
where the first equality is the definition of the adjoint and the second the conjugate symmetry of the inner product (Inner products and norms). So T† is represented by S†. Because W is invariant under both, the restriction of a product is the product of the restrictions, so TT†=T†T on Cn forces SS†=S†S, and T=T† forces S=S†: a Hermitian matrix restricts to a Hermitian matrix, and a normal one to a normal one.
Eigenspaces from an eigenbasis. If T has an orthonormal eigenbasis u1,…,un with Tuk=τkuk, then for each distinct eigenvalue τ the span Fτ of those uk with τk=τ is the whole eigenspace Eτ(T). Indeed Fτ⊆Eτ(T); and if Tx=τx, then writing x=∑k(x,uk)uk and applying T gives ∑kτk(x,uk)uk=∑kτ(x,uk)uk, so (τk−τ)(x,uk)=0 for every k and the coefficients with τk=τ vanish, putting x in Fτ. The distinct Fτ are mutually orthogonal and together span Cn.
Lemma: two commuting Hermitian matrices have a common orthonormal eigenbasis. Let H and K be Hermitian with HK=KH. By the Hermitian spectral theorem (Spectral theorem for Hermitian matrices) and the previous paragraph, Cn is the orthogonal direct sum of the eigenspaces Eμ(H) over the distinct eigenvalues μ of H. Each is invariant under K: if Hx=μx then
H(Kx)=(HK)x=(KH)x=K(Hx)=μ(Kx),
which is where the commutation hypothesis enters, and it is invariant under K†=K for the same reason. By the restriction lemma, K is represented on Eμ(H) by a Hermitian matrix, so the Hermitian spectral theorem applied a second time — now inside Eμ(H) — produces an orthonormal basis of Eμ(H) made of eigenvectors of K. Every vector of Eμ(H) is an eigenvector of H, so these are eigenvectors of both. Collecting them over all μ gives n vectors, orthonormal because they are orthonormal inside each eigenspace and the eigenspaces are mutually orthogonal, hence an orthonormal basis.
Part 1. Let A be normal and put
H=21(A+A†),K=2i1(A−A†),soA=H+iK,A†=H−iK.
Both are Hermitian: H†=21(A†+A)=H, and K†=−2i1(A†−A)=K. Multiplying out,
AA†=H2+K2+i(KH−HK),A†A=H2+K2+i(HK−KH),
so AA†−A†A=2i(KH−HK) and the normality hypothesis says precisely that HK=KH. The lemma then gives an orthonormal basis u1,…,un with Huk=αkuk and Kuk=βkuk, the αk,βk real, whence
Auk=(H+iK)uk=(αk+iβk)uk=:λk′uk.
So Cn has an orthonormal basis of eigenvectors of A. The matrix U=(u1⋯un) has orthonormal columns, so it is unitary (Unitary and orthogonal matrices, part 2 of the characterisation), and
(U†AU)jk=uj†Auk=λk′uj†uk=λk′δjk,
that is U†AU=Λ. The eigenspace decomposition Cn=Eλ1⊕⋯⊕Eλp with mutually orthogonal summands is then the "eigenspaces from an eigenbasis" paragraph applied to A.
Part 3. First, for any normal B and any x,
∥Bx∥2=x†B†Bx=x†BB†x=∥B†x∥2.
Next, B=A−λI is normal whenever A is: expanding,
B†B=A†A−λA−λA†+∣λ∣2I,BB†=AA†−λA−λA†+∣λ∣2I,
and these agree because A†A=AA†. Combining the two displays with B=A−λI and B†=A†−λI,
Ax=λx⟺∥(A−λI)x∥=0⟺∥(A†−λI)x∥=0⟺A†x=λx,
using that a vector of norm zero is zero. The same non-zero x witnesses both, so σ(A†)={λ:λ∈σ(A)}. (This is the route of Riley, Hobson and Bence, who note that A−λI inherits normality and then read off B†x=0 from x†B†Bx=0.)
Part 4. Let Ax=λx and Ay=μy with λ=μ. Using part 3 at y, so that A†y=μy, and the conjugate-linearity of (⋅,⋅) in its second slot,
λ(x,y)=(Ax,y)=y†Ax=(x,A†y)=(x,μy)=μ(x,y).
Hence (λ−μ)(x,y)=0, and λ=μ forces (x,y)=0. Within a single eigenspace there is no such constraint, but Eλ is a subspace, so any linearly independent spanning set of it can be replaced by an orthonormal basis of it by Gram-Schmidt, and every vector of Eλ is still an eigenvector for λ.
Part 2. From U†AU=Λ we get A=UΛU†, and reading the product columnwise, A=∑k=1nλk′ukuk†. For each distinct eigenvalue λj let Ij={k:λk′=λj} and Pj=∑k∈Ijukuk†. Then Pj†=Pj, and PiPj=∑k∈Ii∑l∈Ijuk(uk†ul)ul†=δijPj by orthonormality, so each Pj is idempotent and PiPj=0 for i=j. Also ∑jPj=∑k=1nukuk†=UU†=I. The image of Pj is exactly Eλj: it is contained in it because each uk with k∈Ij lies there, and it contains it because x∈Eλj has (x,uk)=0 for k∈/Ij by the "eigenspaces from an eigenbasis" paragraph, so x=∑k(x,uk)uk=Pjx. A Hermitian idempotent with image Eλj is the orthogonal projection onto Eλj. Finally, grouping the terms of A=∑kλk′ukuk† by eigenvalue gives A=∑jλjPj. □
Theorem(Normality is exactly unitary diagonalizability)Original
Conversely, if U†AU=Λ is diagonal for some unitary U, then A=UΛU† and
AA†=UΛΛU†=UΛΛU†=A†A,
because diagonal matrices commute. So for an n×n complex matrix A:
A is normal⟺Cn has an orthonormal basis of eigenvectors of A⟺U†AU is diagonal for some unitary U.
Reading off the diagonal gives the sharper classification: a normal matrix is Hermitian exactly when σ(A)⊂R, anti-Hermitian exactly when σ(A)⊂iR, unitary exactly when σ(A) lies on the unit circle, and an orthogonal projection exactly when σ(A)⊆{0,1}.
Idea of the proofOriginal
Once A=UΛU†, every identity between A and A† becomes the same identity between Λ and Λ, because conjugating by a fixed unitary matrix is injective and respects products. So each of the four special classes is read off the diagonal entry by entry: being Hermitian is λ=λ, being unitary is ∣λ∣=1, and so on.
ProofOriginal
The equivalence. If A is normal, part 1 of the previous statement gives an orthonormal eigenbasis and the unitary U with U†AU=Λ. If Cn has an orthonormal eigenbasis of A, the matrix U of those columns is unitary and the computation (U†AU)jk=uj†Auk=λk′δjk of that proof makes U†AU diagonal. If U†AU=Λ is diagonal for some unitary U, then the display above shows AA†=A†A, since A†=UΛU† and diagonal matrices commute. The three conditions therefore imply one another in a cycle.
The classification. Let A=UΛU† with U unitary, so A†=UΛU†, and note that UMU†=UNU† implies M=N, because U is invertible. Each line below is read entry by entry on the diagonal, and σ(A) is the set of diagonal entries of Λ.
A=A†⟺Λ=Λ⟺ every λk′ is real ⟺σ(A)⊂R.
A†=−A⟺Λ=−Λ⟺Reλk′=0 for every k⟺σ(A)⊂iR.
A†A=I⟺ΛΛ=I⟺∣λk′∣2=1 for every k⟺σ(A) lies on the unit circle.
A is an orthogonal projection, i.e. A=A†=A2, ⟺Λ=Λ and Λ2=Λ⟺ every λk′ is real and satisfies λ2=λ⟺σ(A)⊆{0,1}, the only solutions of λ2=λ being 0 and 1.
In each case the hypothesis that A is normal is what supplies the U and Λ in the first place; without it none of these spectral conditions says anything, as (0010) shows — its spectrum is {0}⊆{0,1} and it is not a projection. □
Let A and B be normal n×n complex matrices. Then some single orthonormal basis of Cn is made of vectors that are eigenvectors of Aand eigenvectors of B — equivalently, one unitary U makes both U†AU and U†BU diagonal — if and only if
[A,B]=AB−BA=0.
In particular this holds for two commuting Hermitian matrices, and then the common basis may be chosen real whenever A and B are real symmetric. The criterion extends to any commuting family of normal matrices.
Idea of the proofOriginal
One direction is immediate: two matrices that are diagonal in the same basis commute, because diagonal matrices do. The other direction turns commutation with A into commutation with each of A's spectral projections Pj — the reason being that B maps each eigenspace of A into itself, so it cannot mix the projections up. Once B's projections Qα and A's projections Pj all commute with each other, the products PjQα are again orthogonal projections, onto the joint eigenspaces Eλj∩Nμα, and they still sum to the identity. So Cn splits orthogonally into subspaces on each of which A and B are both scalars, and any orthonormal basis adapted to that splitting is a common eigenbasis.
Write A=∑jλjPj and B=∑αμαQα for the spectral resolutions of part 2 of the spectral theorem, with Pj the orthogonal projection onto Eλj=ker(A−λjI) and Qα that onto Nμα=ker(B−μαI); these exist because A and B are normal.
Necessity. Suppose u1,…,un is a basis with Auk=αkuk and Buk=βkuk. Then
ABuk=βkAuk=αkβkuk=αkBuk=BAuk
for every k, and two matrices that agree on a basis are equal, so AB=BA. (Only the existence of a common eigenbasis is used here, not its orthonormality.)
A lemma. If a matrix T commutes with the normal matrix A, then T commutes with every spectral projection Pj of A. For if Ax=λjx then A(Tx)=T(Ax)=λj(Tx), so each Eλj is carried into itself by T. Now take any x and write x=∑iPix with Pix∈Eλi. For i=j we have PjTPix=0, because TPix∈Eλi is fixed by Pi and PjPi=0; and PjTPjx=TPjx, because TPjx lies in Eλj, on which Pj is the identity. Hence
PjTx=i∑PjTPix=PjTPjx=TPjx,
which is the claim. This is the lemma of Hassani, in the form we need.
All four families of projections commute. Applying the lemma to the normal matrix A and T=B gives BPj=PjB for every j. Applying it now to the normal matrix B and T=Pj — legitimate because Pj commutes with B — gives PjQα=QαPj for every j and α.
the middle step of the second line using the commutation just proved. So each Rjα is a Hermitian idempotent, hence the orthogonal projection onto its image, and that image is
Wjα=Eλj∩Nμα:
every Rjαx=Pj(Qαx) lies in the image of Pj, which is Eλj, and equals Qα(Pjx), which lies in Nμα; conversely a vector in both is fixed by Pj and by Qα, hence by Rjα. Moreover
j,α∑Rjα=(j∑Pj)(α∑Qα)=I⋅I=I,
and RjαRj′α′=PjPj′QαQα′=0 unless j=j′ and α=α′. Therefore every x decomposes as x=∑j,αRjαx into vectors of the Wjα, and these subspaces are mutually orthogonal, since for (j,α)=(j′,α′) and x∈Wjα, y∈Wj′α′ we get (x,y)=(Rjαx,Rj′α′y)=y†Rj′α′Rjαx=0. Hence Cn is the orthogonal direct sum of the non-zero Wjα.
The common basis. Pick an orthonormal basis of each non-zero Wjα (Gram-Schmidt). The union is orthonormal, because the Wjα are mutually orthogonal, and spans Cn, so it is an orthonormal basis; each of its members lies in some Eλj∩Nμα and is therefore an eigenvector of A with eigenvalue λjand of B with eigenvalue μα. Taking U to be the matrix of these columns, U is unitary and both U†AU and U†BU are diagonal.
The real symmetric case. If A and B are real symmetric and commute, then each Wjα is the solution space of the real homogeneous linear system (A−λjI)x=0, (B−μαI)x=0 — real because the eigenvalues of a real symmetric matrix are real — so it is the complexification of the real solution space Wjα∩Rn, which has the same dimension. Running Gram-Schmidt over R inside each real solution space gives a common real orthonormal eigenbasis, and the matrix of those columns is real orthogonal.
A commuting family. Let F be a set of pairwise commuting normal matrices. Refine as above: start from the eigenspace decomposition of one member and intersect repeatedly with the eigenspaces of further members, choosing at each step a member that splits at least one summand in two. Each such step strictly increases the number of non-zero summands, and that number never exceeds n, so after at most n−1 steps no member of F splits any summand. Call the resulting orthogonal decomposition Cn=W1⊕⋯⊕Wr. Each Wm is an intersection of eigenspaces of finitely many members, so by the lemma — applied once for each of those members — every A∈F commutes with the projection onto Wm and hence carries Wm into itself; and since A splits no summand, Wm lies inside a single eigenspace of A, i.e. A acts on Wm as a scalar. An orthonormal basis of each Wm is then a common eigenbasis of the whole family. □
NoteOriginal
Normal is the exact condition for "diagonal in some orthonormal frame": the matrix is a list of complex stretch factors along n perpendicular axes. The commutation criterion is the algebraic reason two observables can be measured at once in quantum mechanics — a common eigenbasis exists precisely when the operators commute — and the practical reason that a symmetry of a system lets you block-diagonalize its dynamical matrix.
Remark(Where the hypotheses are used)Original
"Orthonormal" cannot be dropped.(1012) has distinct eigenvalues 1,2 and so is diagonalizable, but AAT=(2224)=(1115)=ATA, so it is not normal and no orthonormal eigenbasis exists. Diagonalizable by a general similarity and diagonalizable by a unitary are different conditions, and only the second is normality.
Normality cannot be dropped from the commutation criterion.A=B=(1011) commute, yet neither is diagonalizable at all. Commuting buys a common invariant decomposition, not a common eigenbasis, unless each matrix separately has an eigenbasis.
Degeneracy leaves a choice, and the choice matters. If λ is a degenerate eigenvalue of A, not every orthonormal eigenbasis of A is an eigenbasis of a commuting B: one must diagonalize B restricted to each eigenspace Eλ of A. The theorem asserts that some such choice works, not that any does.
Real normal is weaker than real symmetric.R(θ)=(cosθsinθ−sinθcosθ) is real and normal, but for sinθ=0 it has no real eigenvector, so it cannot be diagonalized by a real orthogonal matrix. Among real matrices, orthogonal diagonalizability over R characterises the symmetric ones, not the normal ones; see unitary and orthogonal matrices for the real canonical form.
Conventions. Physics texts usually write the condition as [A,A†]=0 and call A† the adjoint; mathematical texts write A∗. Some books state part 1 of the spectral theorem as "a normal operator is diagonalizable"; read in such a statement that the diagonalizing basis is orthonormal, which is what the proof produces and what the converse above needs.
Example(A normal matrix that is neither Hermitian nor unitary, and a commuting partner)Original
Let C be the cyclic shift of C3 and A=I+C:
C=001100010,A=101110011.
Step 1: A is normal, and why.C is a permutation matrix, so its columns are the standard basis vectors in some order and are orthonormal: C is unitary, with C†=C−1=C2 (shifting three times returns to the start, C3=I). Hence
A†=I+C2,AA†=(I+C)(I+C2)=I+C+C2+C3=2I+C+C2,
and A†A=(I+C2)(I+C) gives the same four terms in the other order, so AA†=A†A. Written out, both equal 211121112. The structural reason is worth keeping: A and A† are both polynomials in the one matrix C, and any two polynomials in the same matrix commute.
A is not Hermitian, since AT=A, and not unitary, as step 2 will show by exhibiting an eigenvalue of modulus 2.
Step 2: the eigenvalues and an orthonormal eigenbasis. Let ω=e2πi/3=−21+23i, so ω3=1 and 1+ω+ω2=0. For k=0,1,2 put
vk=(1,ωk,ω2k)T.
Then Cvk=(ωk,ω2k,1)T=ωkvk, using ω3k=1 for the last entry, so
Two checks: trA=3=2+2cos(π/3)=2+1, and detA=λ0λ1λ2=2∣eiπ/3∣2=2, which matches the direct expansion detA=1(1−0)−1(0−1)=2. Since ∣λ0∣=2=1, A is not unitary, as promised.
The vk are orthogonal: by the geometric-sum formula,
because for j=k the sum is 1+ζ+ζ2=0 with ζ=ωj−k=1 a cube root of unity. So uk=vk/3 is an orthonormal eigenbasis — exactly what part 1 of the spectral theorem promises — and with U=(u0u1u2), the discrete Fourier matrix of order 3,
U†AU=diag(2,eiπ/3,e−iπ/3).
Step 3: the adjoint and the spectral resolution. Since C2uk=ω2kuk=ωkuk (because ∣ω∣=1, so ω2k=ω−k=ωk), we get A†uk=(1+ωk)uk=λkuk: the same eigenvectors, the conjugated eigenvalues, as part 3 requires. The eigenvalues are distinct, so each eigenspace is one-dimensional, Pk=ukuk†, and part 2 reads
A=2u0u0†+eiπ/3u1u1†+e−iπ/3u2u2†,
which multiplies out to A (one may check the (1,1) entry: 31(2+eiπ/3+e−iπ/3)=31(2+1)=1, and the (1,2) entry: 31(2+eiπ/3ω+e−iπ/3ω2)=31(2+1)=1, since eiπ/3ω=eiπ/3e−2πi/3=e−iπ/3 and the two conjugate terms sum to 2cos(π/3)=1).
Step 4: a commuting partner and the role of degeneracy. Let
B=A+A†=2I+C+C2=211121112,
real symmetric, and commuting with A because both are polynomials in C. By the simultaneous-diagonalization theorem the two share an orthonormal eigenbasis, and indeed Buk=(λk+λk)uk=2(Reλk)uk, giving the eigenvalues
4,1,1.
So B has a degenerate eigenvalue 1, with the two-dimensional eigenspace N1=span{v1,v2}, and this is where the third trap above becomes concrete. Solving (B−I)x=0 by hand gives the obvious real basis (−1,1,0)T, (−1,0,1)T of N1, and
A(−1,1,0)T=(0,1,−1)T,
which is not a multiple of (−1,1,0)T: that perfectly good eigenvector of B is not an eigenvector of A. Only a particular basis of N1 — up to scalars, {v1,v2} — diagonalizes A as well. The theorem says such a choice exists, not that every choice works.
The same example illustrates the fourth trap. A is real and normal, but its eigenvalues e±iπ/3 are not real, so it has no real eigenvector for them and cannot be diagonalized by a real orthogonal matrix, although the real symmetric B can be: for instance by the real orthonormal basis (1,1,1)T/3, (1,−1,0)T/2, (1,1,−2)T/6, which diagonalizes B but not A.
ApplicationOriginal
Commuting observables and good quantum numbers. In quantum mechanics an observable is a Hermitian operator and a measurement leaves the system in an eigenvector of it. Two observables can therefore be assigned definite values at the same time exactly when they have a common eigenvector, and a complete set of such states exists exactly when they commute — which is the simultaneous-diagonalization theorem, and the reason [⋅,⋅] is the central object of the theory. Position and momentum do not commute, so no state has both sharp; the classification of the hydrogen atom, by contrast, rests on the fact that H, L2 and Lz do commute pairwise. The proof above explains why the labels come in a tuple: the joint eigenspaces are the subspaces Wjα=Eλj∩Nμα, and the quantum numbers (n,ℓ,m) are just the names of the eigenvalues whose intersection the state lies in. Degeneracy is the practical point. An energy level En of hydrogen is highly degenerate, so "an eigenbasis of H" is far from determined; the extra commuting operators are precisely what selects a basis inside each level, and the fourth trap above warns that an arbitrary basis of a degenerate level will not do.
Symmetry block-diagonalizes a dynamical matrix. The worked example is, read physically, a ring of three identical masses with identical nearest-neighbour springs: the matrix B is (up to a constant) the dynamical matrix of that ring, and C is the operator that rotates the ring by one site. Rotating the ring is a symmetry of it, so C commutes with B, and the theorem says the normal modes can be chosen to be eigenvectors of C as well. Those eigenvectors are the discrete plane waves vk, with Cvk=e2πik/3vk — the finite version of Bloch's theorem, the phase 2πk/3 playing the role of the wavenumber. Note what has been bought: the modes were found by diagonalizing the symmetry, which was trivial, rather than the dynamical matrix. For a ring of N masses this replaces an N×N eigenproblem by inspection, and the same argument in the general form — decompose into joint eigenspaces of the symmetry group — is how the vibrational modes of a molecule are sorted by irreducible representation before any force constant is computed. The degeneracy of the ring (λ1=λ2 in the example) is itself a consequence of the symmetry, since v1 and v2 are exchanged by reflecting the ring.
Why numerical analysts care whether a matrix is normal. The eigenvalues of a normal matrix are perfectly conditioned: because the diagonalizing matrix U is unitary, a perturbation A+ΔA moves each eigenvalue by at most ∥ΔA∥ in the Euclidean operator norm. For a non-normal matrix the eigenvector matrix can be badly conditioned and the eigenvalues correspondingly violent: the n×n matrix with ones on the superdiagonal and ε in the bottom-left corner has eigenvalues of modulus ε1/n, which is 10−1 for ε=10−10 and n=10, although the unperturbed matrix is nilpotent with every eigenvalue 0. So "is this matrix normal?" decides whether computed eigenvalues can be trusted at face value, and it is the reason library routines dispatch Hermitian and unitary input to separate, more accurate algorithms.
Every square matrix whose eigenvalues lie in the ground field is similar to a block diagonal matrix of Jordan blocks, unique up to the order of the blocks; the block sizes are read off the nullities of the powers of A−λI, and the form turns powers and power series of a defective matrix into polynomial-times-exponential formulas.
Definition(Jordan block)AdaptedAdapted from Dummit & Foote (2004), §12.3, the definition on p. 492
For λ in a field F and an integer m≥1, the Jordan blockJm(λ) is the m×m matrix carrying λ in every diagonal position, 1 in every position of the first superdiagonal, and 0 everywhere else:
Dummit and Foote (Dummit & Foote (2004), p. 492) call this the elementary Jordan matrix with eigenvalue λ. A matrix is in Jordan canonical form when it is block diagonal with Jordan blocks as its diagonal blocks. The nilpotent part Nm shifts the standard basis down one step, Nme1=0 and Nmet=et−1 for t≥2, so Nmm=0 while Nmm−1=0.
The block Jm(λ) is upper triangular, so its characteristic polynomial is (x−λ)m: the eigenvalue λ has algebraic multiplicity m. Its eigenspace, however, is the line spanned by e1, because (Jm(λ)−λI)x=0 reads xt=0 for t≥2. So a single block of size m≥2 is defective in the sense of the eigenvalue problem, with geometric multiplicity 1 against algebraic multiplicity m.
Definition(Generalized eigenvectors, generalized eigenspaces and Jordan chains)Original
Let A be n×n over F and λ∈F. A vector x=0 is a generalized eigenvector of A for λ if
(A−λI)jx=0for some j≥1,
and the smallest such j is its rank. Rank 1 means an ordinary eigenvector. The generalized eigenspace is
Gλ=ker(A−λI)n,
which contains every generalized eigenvector for λ, because the kernels of the powers stop growing after at most n steps (step 1 of the proof below). Write
dj(λ)=dimker(A−λI)j=n−rank(A−λI)j,d0(λ)=0,
so that d1(λ) is the geometric multiplicity of λ and dn(λ)=dimGλ.
A Jordan chain of length s for λ is a list x1,…,xs of non-zero vectors with
(A−λI)x1=0,(A−λI)xt=xt−1(2≤t≤s),
that is, xt=(A−λI)s−txs with xt of rank t. Written as the columns of a matrix in this order, a chain is exactly a basis on which A acts as Js(λ): Ax1=λx1 and Axt=λxt+xt−1.
Finally, some non-zero polynomial does annihilate A: the n2+1 matrices I,A,A2,…,An2 cannot be independent in the n2-dimensional space of n×n matrices, so a non-trivial linear relation among them exhibits a non-zero polynomial q of degree at most n2 with q(A)=0, and dividing by its leading coefficient makes it monic. The minimal polynomialμA of A is therefore well defined as the monic polynomial of least degree with μA(A)=0. It divides every polynomial that annihilates A, since dividing such a polynomial by μA leaves a remainder of lower degree that also annihilates A, which forces the remainder to be 0.
Theorem(Jordan normal form)AdaptedAdapted from Dummit & Foote (2004), §12.3, Thm 23 and Cor 24, p. 493
Let A be an n×n matrix over a field F, and suppose every eigenvalue of A already lies in F: the characteristic polynomial is a product of linear factors there. (Over C this costs nothing.) Then:
Existence. Some invertible matrix P with entries in F puts A into Jordan canonical form by similarity:
each νi an eigenvalue of A, the same eigenvalue permitted in several blocks. The columns of P are b Jordan chains laid end to end.
2. Uniqueness.J is determined by A up to the order of its blocks.
3. Should a diagonal matrix D happen to be similar to A, then D is itself the Jordan form, since diagonal entries are 1×1 blocks and part 2 leaves no second answer (Dummit & Foote (2004), Cor 24, p. 493). So nothing is lost by always asking for the Jordan form rather than for a diagonalization.
Theorem(The block structure is read off the nullities)Original
With dj(λ) as above, in the Jordan form of A:
the number of blocks carrying the eigenvalue λ is d1(λ), the geometric multiplicity of λ;
the number of those blocks of size at leastj is dj(λ)−dj−1(λ), and the number of size exactly j is 2dj(λ)−dj−1(λ)−dj+1(λ);
the total of their sizes is dimGλ, which equals the algebraic multiplicity of λ;
the largest of their sizes is the multiplicity of λ as a root of the minimal polynomial μA; consequently μA(x)=∏λ(x−λ)s(λ) with s(λ) that largest size.
In particular A is diagonalizable exactly when every block has size 1, which by 4 happens exactly when μA has no repeated root (Dummit & Foote (2004), Cor 25, p. 494), and by 1 and 3 exactly when d1(λ) equals the algebraic multiplicity for every λ.
Idea of the proofOriginal
Two independent mechanisms are at work, and the proof separates them.
The first is that different eigenvalues do not interfere. If T=A−λI, then the kernels of T,T2,T3,… grow and must stall, and at the moment they stall the space splits as kerTn⊕imTn into a part where T is nilpotent and a part where T is invertible. Peeling off one eigenvalue at a time this way decomposes Fn into the generalized eigenspaces Gλ, on each of which A is λI plus a nilpotent matrix.
The second is a statement about one nilpotent matrix N: its domain has a basis made of chains Ns−1u,…,Nu,u. This is proved by induction on the dimension, applying the result to the smaller space N(V) and then lengthening every chain found there by one step — possible precisely because the chains in N(V) start from vectors that are themselves images.
Uniqueness needs no new idea: the numbers dj(λ) are unchanged by similarity, and a direct count shows a block of size s contributes min(j,s) to dj(λ). Inverting that count recovers the block sizes from the dj(λ), so no two different block patterns can describe the same matrix.
ProofOriginal
Throughout, V=Fn and matrices act on column vectors. Dummit and Foote reach this theorem as a corollary of the elementary divisor decomposition of V as a module over the polynomial ring F[x](Dummit & Foote (2004), §12.3, pp. 490-493); the argument below is self-contained and uses nothing beyond rank and nullity.
Step 1: the kernels of the powers stall, and the space splits. Let T:V→V be linear, d=dimV. Since Tjx=0 implies Tj+1x=0, the kernels increase: kerT⊆kerT2⊆⋯. Suppose kerTj=kerTj+1 for some j. If Tj+2x=0 then Tj+1(Tx)=0, so Tx∈kerTj+1=kerTj, so Tj+1x=0: the chain is constant from j on. As each strict increase costs at least one dimension and dimV=d, equality is reached at the latest at j=d. Hence
kerTd=kerTd+1=⋯=:K,R:=imTd.
Then V=K⊕R. Indeed dimK+dimR=d by rank-nullity applied to Td, and the sum is direct: if x=Tdy lies in K, then T2dy=0, so y∈kerT2d=kerTd and therefore x=Tdy=0. Both summands are carried into themselves by T (and by any matrix commuting with T), T is nilpotent on K because Td=0 there, and T is injective on R — a vector of R killed by T lies in K∩R={0} — hence invertible on R.
Step 2: V is the direct sum of the generalized eigenspaces. We show by induction on dimV that if the characteristic polynomial of A splits over F, with distinct roots λ1,…,λp, then V=Gλ1⊕⋯⊕Gλp.
The case dimV=0 is empty. Otherwise apply step 1 to T=A−λ1I, giving V=K⊕R with K=ker(A−λ1I)n=Gλ1 and both summands invariant under A (because A commutes with T). Here K={0}, since λ1 is a root of the characteristic polynomial and so has an eigenvector, whence dimR<dimV. Choosing a basis adapted to the splitting makes A block diagonal, so its characteristic polynomial is the product of those of A∣K and A∣R; in particular the latter splits over F as well. On K we have A=λ1I+T with T nilpotent, hence strictly upper triangular in a suitable basis of K — choose a basis adapted to the flag {0}⊆kerT⊆kerT2⊆⋯⊆K, which exhausts K by step 1, and note that T carries kerTj into kerTj−1, so each basis vector is sent into the span of the strictly earlier ones (step 3 below gives the sharper chain basis, and does not depend on the present step). Hence the characteristic polynomial of A∣K is (x−λ1)dimK and the remaining roots λ2,…,λp are exactly the eigenvalues of A∣R. Also, for λ=λ1 the matrix A−λI=(λ1−λ)I+T is invertible on K — a non-zero scalar plus a nilpotent matrix is invertible, its inverse being the finite sum ∑i≥0(−1)i(λ1−λ)−i−1Ti. Consequently, for j≥2,
Gλj⊆R:
write x∈Gλj as x=y+z with y∈K, z∈R. Because K and R are invariant, (A−λjI)ny lies in K and (A−λjI)nz lies in R, and their sum is 0; the decomposition of 0 being unique, each vanishes separately. Then y=0 by the invertibility just noted, so x=z∈R. The inductive hypothesis applied to A∣R decomposes R into the generalized eigenspaces of A∣R, and the one for λj is ker(A−λjI)dimR∩R, which equals Gλj: it is contained in Gλj because the kernels only grow with the exponent, and contains it because Gλj⊆R and the kernels have stalled by exponent dimR inside R, by step 1 applied there. Adding K=Gλ1 gives the claim.
Two consequences. First, the computation just made of the characteristic polynomial of A∣K applies to every Gλ, so dimGλ is the algebraic multiplicity of λ. Second, it suffices to produce a Jordan basis of each Gλ separately.
Step 3: a nilpotent matrix has a basis of chains. Let N:V→V be nilpotent, dimV=d. We claim V has a basis that is a disjoint union of chains Ns−1u,…,Nu,u with Nsu=0, and that the ends Ns−1u of those chains form a basis of kerN.
Induct on d. For d=0 there is nothing to prove. If N=0, any basis works, each vector a chain of length one, and the ends span kerN=V. Otherwise set W=N(V) and k=dimkerN, so dimW=d−k by rank-nullity, and k≥1 because a nilpotent matrix on a non-zero space has non-trivial kernel (take the largest i with Ni=0 and any x with Nix=0; then Nix∈kerN). Hence dimW<d, and N restricts to a nilpotent map of W. By the inductive hypothesis W has a chain basis, say chains
Ci:Nsi−1ui,…,Nui,ui(i=1,…,r),Nsiui=0,
whose ends Nsi−1ui form a basis of kerN∩W=ker(N∣W); in particular r=dim(kerN∩W).
Now lengthen. Each ui lies in W=N(V), so choose vi∈V with Nvi=ui; the list
Ci′:Nsivi(=Nsi−1ui),…,Nvi(=ui),vi
is a chain of length si+1, since Nsi+1vi=Nsiui=0. Extend the r vectors Nsi−1ui to a basis of kerN by adjoining z1,…,zk−r, each of which is a chain of length one.
The resulting list has ∑i(si+1)+(k−r)=(d−k)+r+(k−r)=d vectors, so it is enough to show it is independent. Suppose
i=1∑rt=0∑siai,tNtvi+j=1∑k−rbjzj=0.
Apply N. The terms t=si die because Nsi+1vi=0, and the zj die because they lie in kerN, leaving ∑i∑t=0si−1ai,tNtui=0. These vectors are precisely the chain basis of W, so every ai,t with t≤si−1 vanishes. What is left of the relation is ∑iai,siNsi−1ui+∑jbjzj=0, a relation among a basis of kerN, so the remaining coefficients vanish too. Finally, the ends of the new chains are the Nsi−1ui together with the zj, which is the chosen basis of kerN, so the second half of the claim is inherited.
Step 4: existence. Apply step 3 to N=(A−λI)∣Gλ for each eigenvalue λ, which is legitimate by step 2. A chain x1,…,xs of N satisfies Ax1=λx1 and Axt=λxt+xt−1, so on the span of that chain, in that order, A is the matrix Js(λ). Concatenating the chains over all λ gives a basis of V by step 2; let P be the matrix of its columns. Then AP=PJ with J block diagonal of Jordan blocks, that is P−1AP=J.
Step 5: uniqueness. Let J be any Jordan matrix similar to A, say P−1AP=J. Since P−1(A−λI)jP=(J−λI)j, similar matrices have the same nullities, so dj(λ) may be computed from J. A kernel of a block diagonal matrix is the direct sum of the kernels of its blocks, so dj(λ) is the sum over blocks Js(ν) of dimker(Js(ν)−λI)j. For ν=λ that block contributes nothing: Js(ν)−λI is triangular with the non-zero number ν−λ on its diagonal, hence invertible, and so is its j-th power. For ν=λ the block contributes dimkerNsj, and Nsj maps et↦et−j (zero when t≤j), so its kernel is spanned by e1,…,emin(j,s). Therefore
dj(λ)=blocks Js(λ) of J∑min(j,s),hencedj(λ)−dj−1(λ)=#{blocks Js(λ):s≥j},
because min(j,s)−min(j−1,s) is 1 when s≥j and 0 otherwise. The right-hand side determines, for every j, how many blocks for λ have size at least j, hence the whole multiset of block sizes for every eigenvalue, hence J up to the order of the blocks. This proves both the uniqueness statement and parts 1, 2 and 3 of the structure statement — part 1 being the case j=1, part 3 following from dn(λ)=∑ss=dimGλ and step 2, and the "exactly j" count in part 2 being the difference of the counts for ≥j and ≥j+1.
For part 4, note that for a block, (Js(λ)−λI)i=Nsi vanishes exactly when i≥s, while Js(λ)−νI is invertible for ν=λ. So a polynomial q annihilates J if and only if it annihilates every block, and q(Js(λ))=0 holds if and only if (x−λ)s divides q: writing q(x)=∑ici(x−λ)i, we get q(Js(λ))=∑iciNsi, whose superdiagonals are the ci for i<s, so it vanishes exactly when c0=⋯=cs−1=0. The monic polynomial of least degree divisible by (x−λ)s(λ) for every λ, with s(λ) the largest block size, is ∏λ(x−λ)s(λ), and μA=μJ because similarity preserves polynomial identities. □
Theorem(Powers and power series of a Jordan block)Original
Let Jm(λ)=λI+Nm as above. For every integer k≥0,
the entry being 0 at every other position — that is, for q<p and also for q−p>k, so that no negative power of λ is ever called for. (Writing the formula with the convention (ik)=0 for i>k would make it read 0⋅λnegative when λ=0 and q−p>k, which is why those positions are excluded outright rather than by a convention. At the far edge of that range, q=p+k with λ=0, the entry is (kk)λ0=1, using 00=1; this is the statement that Jm(0)k=Nmk.) More generally, if f(z)=∑i≥0ci(z−λ)i converges in a disc about λ, then
f(Jm(λ))=i=0∑m−1i!f(i)(λ)Nmi,
a finite sum: the block sees only the first m Taylor coefficients of f at λ. Applying this blockwise through A=PJP−1 gives, for an n×n matrix A with largest block size s(λ) at each eigenvalue,
(Ak)pq=λ∈σ(A),λ=0∑φpqλ(k)λk,
with each φpqλ a polynomial in k of degree at most s(λ)−1. The sum omits a zero eigenvalue, and the formula then holds for all k≥s(0), where s(0) is the largest block size at 0 (and for all k≥1 when A is invertible, the omission being vacuous): the blocks at 0 are nilpotent of index at most s(0), so they contribute to Ak only while k<s(0), and a term φ(k)0k could not record that contribution. Hence Ak→0 if and only if every eigenvalue satisfies ∣λ∣<1; and Ak stays bounded as k→∞ if and only if every eigenvalue satisfies ∣λ∣≤1 and every eigenvalue with ∣λ∣=1 has s(λ)=1. The same formula with f(z)=etz gives the polynomial-times-exponential shape of exp(tA); the computational route is in the matrix exponential.
ProofOriginal
λI commutes with Nm, so the binomial theorem holds for (λI+Nm)k exactly as for numbers — the usual proof by induction only ever needs that the two terms commute. Terms with i≥m drop out because Nmm=0, and the sum runs only to i=k, so no term carries a negative power of λ. The entry form follows from (Nmi)pq=1 when q=p+i and 0 otherwise, which is the statement that Nmi shifts the basis i places: the position (p,q) receives a contribution from the single term i=q−p when 0≤q−p≤min(k,m−1), and from no term at all otherwise, which is why the entry vanishes for q<p and for q−p>k whatever λ may be.
For the power series, substitute z=Jm(λ) into f(z)=∑ici(z−λ)i, so that (z−λ)i becomes Nmi, which vanishes for i≥m; only i≤m−1 survives, the sum is finite, and ci=f(i)(λ)/i! is Taylor's formula. (Convergence of the series is therefore not at issue for a single block, though it is for the eigenvalue to lie in the disc of convergence of f at all.)
For A=PJP−1 we get Ak=PJkP−1, and Jk is block diagonal with the blocks Jsi(νi)k just computed. Every entry of Ak is therefore a fixed linear combination of the entries of the blocks, each of which is (ik)λk−i for some i≤s(λ)−1; and (ik)=k(k−1)⋯(k−i+1)/i! is a polynomial in k of degree i. Grouping by eigenvalue and writing (ik)λk−i=((ik)λ−i)λk gives the stated form with degφpqλ≤s(λ)−1. That rewriting needs λ=0, which is why the eigenvalue 0 is excluded from the sum; once k≥s(0) every block at 0 has already died, since Nsk=0 for s≤s(0)≤k, so the omitted blocks contribute nothing and the formula is exact from that k on.
For the limit: if every ∣λ∣<1 then each ki∣λ∣k→0, since the geometric factor beats any polynomial, so Ak→0. Conversely, if Ax=λx with x=0 then Akx=λkx, which tends to 0 only if ∣λ∣<1. For boundedness both conditions are necessary, and here Jk=P−1AkP lets us argue on J: the two are bounded or unbounded together. Suppose Ak is bounded. If some eigenvalue had ∣λ∣>1, an eigenvector x=0 would give ∥Akx∥=∣λ∣k∥x∥→∞; so every eigenvalue has ∣λ∣≤1. If some eigenvalue with ∣λ∣=1 had s(λ)=m≥2, the entry (m−1k)λk−m+1 in the corner of a largest block for that λ would have modulus (m−1k)→∞, making Jk and hence Ak unbounded; so every eigenvalue of modulus 1 has s(λ)=1. Conversely, if every eigenvalue of modulus 1 has s(λ)=1 and all others have ∣λ∣<1, each entry of Jk is either λk with ∣λ∣=1 or one of the terms (ik)λk−i with ∣λ∣<1, which tend to 0; so Jk, and with it Ak, is bounded. □
The field must contain the eigenvalues. This hypothesis of Dummit and Foote (Dummit & Foote (2004), Thm 23, p. 493) enters at step 2, where an eigenvalue λ1∈F is needed to start the peeling, and again in the induction. It is not a formality: over F=R the rotation (01−10) has characteristic polynomial x2+1, no real eigenvalue and hence no real Jordan form, although over C it is diag(i,−i). What a real matrix does have is the real canonical form of Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), §5.5, Proposition, p. 94): a block diagonal real matrix whose blocks are either real Jordan blocks or, for a conjugate pair α±iβ with β=0, block-upper-bidiagonal matrices built from C2=(α−ββα) on the diagonal and 2×2 identity matrices one place above. One obtains it from the complex Jordan form by replacing each conjugate pair of chains with the real and imaginary parts of one of them.
Uniqueness is only up to the order of the blocks, and P is far from unique. Even for A=λI, every invertible P works. In general, scaling a whole chain by a non-zero constant, or swapping two chains of the same length for the same eigenvalue, changes P and leaves J alone; the form is an invariant of A, the basis is not.
The eigenvalues and the geometric multiplicities do not determine J. Blocks of sizes (3,1) and (2,2) for the same λ both give algebraic multiplicity 4 and geometric multiplicity 2; they are told apart only by d2(λ), which is 3 in the first case and 4 in the second. This is why the structure statement needs the whole sequence d1,d2,… and not just d1.
The form is discontinuous in A, so it cannot be computed in floating point. For ε=0 the matrix (0ε10) has distinct eigenvalues ±ε and is diagonalizable; at ε=0 it is a single 2×2 Jordan block. Arbitrarily small perturbations therefore change the answer completely, and the ranks rank(A−λI)j that the method below uses are exactly the quantities that rounding destroys. Jordan form is a tool for exact arithmetic and for theory; numerical work uses the Schur form or the singular value decomposition instead, and the relevant sensitivity is measured by condition numbers. Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), §5.6, p. 100) make the complementary point that a matrix met at random has distinct eigenvalues, so the non-trivial Jordan structure that matters in practice arises from symmetry or from a tuned parameter, not by accident.
Normal matrices never need it. A normal matrix is unitarily diagonalizable (Normal matrices and simultaneous diagonalization), so all its Jordan blocks have size 1; non-trivial Jordan structure is a strictly non-normal phenomenon.
Conventions. Many texts, including Hassani (Hassani (2013), §17.7, p. 539), write the 1s on the superdiagonal as here; others put them on the subdiagonal, which is the transpose and equally valid. Some books, among them Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §6.6, p. 270), present the form with one block per distinct eigenvalue, of size equal to that eigenvalue's multiplicity; that is the special case in which each eigenvalue has geometric multiplicity 1, and part 1 of the structure statement is what allows several blocks in general. We follow the several-block statement of Dummit and Foote (Dummit & Foote (2004), Thm 23, p. 493) throughout.
Method(Computing the Jordan form and a Jordan basis)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §5.5, the proof on pp. 95-96 and the worked matrices on pp. 96-99
Given an n×n matrix A over C with exact entries:
Eigenvalues. Factor the characteristic polynomial χA and record each distinct eigenvalue λ with its algebraic multiplicity a(λ).
Nullities. For each λ, compute rj=rank(A−λI)j for j=1,2,… until rj stops decreasing, and set dj=n−rj. Stop at the first j with dj=a(λ); this j is the largest block size for λ.
Block sizes. The number of blocks for λ of size at least j is dj−dj−1. Equivalently, the list d1−d0,d2−d1,… (automatically non-increasing, since a block of size ≥j+1 has size ≥j) is the conjugate partition of the list of block sizes. Check that the sizes sum to a(λ).
Chains, longest first. For the largest size s, pick v∈ker(A−λI)s not in ker(A−λI)s−1 and generate the chain downwards by repeated multiplication:
v,(A−λI)v,(A−λI)2v,…,(A−λI)s−1v.
Every later chain for the same λ, whether another of length s or one of a smaller length s′, is started from a vector chosen by one and the same rule. Abbreviate T=A−λI and let U be the span of all chain vectors produced so far; then take
v∈kerTs′withv∈/kerTs′−1+U,
and run the chain v,Tv,…,Ts′−1v downwards as before.
The second condition is about the sum of the two subspaces, not about each of them. Requiring only v∈/kerTs′−1 and v∈/U is genuinely too weak. Take λ=0 and A the direct sum of J3(0) and J2(0), with the chains a1,a2,a3 and b1,b2 as its basis, so Aa3=a2, Aa2=a1, Aa1=0, Ab2=b1, Ab1=0. Steps 1 to 3 give d1=2, d2=4, d3=5 and block sizes 3,2. Choosing v=a3 produces U=span{a1,a2,a3}, and then v′=a2+b1 satisfies all three of the weaker conditions — A2v′=0, Av′=a1=0, v′∈/U — yet its chain is a1,a2+b1, which repeats a1: the assembled list has rank 4, not 5, and P is singular. The rule above rejects v′, since a2+b1∈U+kerA, and accepts b2, which is the correct choice.
Why the rule is the right one. For v∈kerTs′, the condition v∈/kerTs′−1+U says precisely that the bottom vector Ts′−1v of the new chain is not in the span of the bottom vectors Tsi−1vi of the chains already built. (All chains built so far have length si≥s′, since the sizes are taken in decreasing order. Suppose v=y+u with Ts′−1y=0 and u∈U. From Ts′v=0 and Ts′y=0 we get Ts′u=0; expanding u=∑i∑t<siai,tTtvi in the chain vectors, the only terms of Ts′u that survive are those with t<si−s′, and they are distinct members of the chain list, so ai,t=0 for every t<si−s′. Applying Ts′−1 to what remains kills every term with t>si−s′ and gives Ts′−1v=Ts′−1u=∑iai,si−s′Tsi−1vi, a combination of bottom vectors. Conversely, if Ts′−1v=∑iciTsi−1vi, then v−∑iciTsi−s′vi lies in kerTs′−1, so v∈kerTs′−1+U.) And chains whose bottom vectors are independent are jointly independent: in a vanishing combination ∑i∑t<siai,tTtvi=0 with some coefficient non-zero, let ci be the least t with ai,t=0 and put q=maxi(si−ci) over the chains that occur; applying Tq−1 annihilates every term except those with si−ci=q and t=ci, and leaves a vanishing combination of distinct bottom vectors with a non-zero coefficient — impossible. Finally, a legal choice is available exactly as often as step 3 predicts: the vectors Tsi−s′vi taken from the chains already built have independent images in kerTs′/kerTs′−1 (apply Ts′−1 and use the independence of the bottoms), that quotient has dimension ds′−ds′−1, the number of blocks of size at least s′, and by the equivalence just proved kerTs′−1+U meets kerTs′ in exactly the preimage of the span of those images. So the supply of choices runs out only when the chains of length ≥s′ already number ds′−ds′−1.
Equivalently, and this is the device Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), §5.5, pp. 95-96) use in their worked 3×3 and 4×4 examples, start from an eigenvector and back-solve(A−λI)xt=xt−1, going up the chain instead of down.
5. Assemble. Write each chain as consecutive columns ending with the top vector v, that is in the order (A−λI)s−1v,…,(A−λI)v,v; concatenate the chains over all eigenvalues to form P. Then J=P−1AP, with the blocks in the order the chains were laid down.
6. Check without inverting P: verify AP=PJ column by column, and verify trA=trJ and detA=detJ.
If only the form is wanted — enough for powers, growth rates and the minimal polynomial — stop after step 3; the eigenvectors are needed only for P. Never run this on inexact data: see the remark above.
Example(A defective 4 x 4 matrix and its powers)Original
Let
A=2−124−1224111−4−1−1−23.
Step 1: eigenvalues. Expanding det(xI−A) gives
χA(x)=x4−8x3+18x2−27=(x−3)3(x+1),
the factorization being confirmed by multiplying out: (x3−9x2+27x−27)(x+1)=x4−8x3+18x2−27, the linear terms −27x and +27x cancelling. The two standard tests agree: trA=2+2+1+3=8=3+3+3−1, and detA=χA(0)=−27=3⋅3⋅3⋅(−1). So a(3)=3 and a(−1)=1.
Call the rows of B in order R1,…,R4. Then R2=R1, and subtracting multiples of R1 from the last two rows leaves the same vector twice, R3+2R1=R4+4R1=(0,0,0,−4), so R4=R3−2R1. The row space is therefore spanned by R1 and R3, which are independent because R3 has a non-zero fourth entry after R1 has been cleared from it: rankB=2 and d1(3)=2. For B2 the two non-zero rows are equal, so rankB2=1 and d2(3)=3=a(3); the chain of kernels has stalled.
Step 3: block sizes. The number of blocks for 3 is d1=2; the number of size at least 2 is d2−d1=1; the number of size at least 3 is d3−d2=0. So the blocks for 3 have sizes 2 and 1, summing to 3=a(3) as required. The eigenvalue −1 is simple, giving one block of size 1. Hence
J=300013000030000−1,
and the minimal polynomial is μA(x)=(x−3)2(x+1), the largest block size at 3 being 2. Indeed B2(A+I)=0 while B(A+I)=0: that product is −4−4−80−4−4−804480−4−4−80.
Step 4: the chains. Take the top of the long chain to be v=(0,0,0,1)T, which lies in kerB2 (the fourth column of B2 is zero) but not in kerB (the fourth column of B is (−1,−1,−2,0)T=0). The chain is
x1=Bv=(−1,−1,−2,0)T,x2=v=(0,0,0,1)T,
and Bx1=0, as the four rows confirm: 1+1−2+0=0 from R1 and again from R2, then −2−2+4+0=0 from R3, then −4−4+8+0=0 from R4.
For the second block at 3 we need an eigenvector independent of x1. Write a general vector as ξ=(ξ1,ξ2,ξ3,ξ4)T and solve Bξ=0: rows 1 and 2 give −ξ1−ξ2+ξ3−ξ4=0 and row 4 gives ξ1+ξ2−ξ3=0, whence ξ4=0 and ξ3=ξ1+ξ2; row 3 is then automatic. A basis of kerB is (−1,1,0,0)T and (1,0,1,0)T, and we may take
w=(−1,1,0,0)T,
which is independent of x1=−(−1,1,0,0)T−2(1,0,1,0)T.
For λ=−1: (A+I)ξ=0 has the one-dimensional solution space spanned by z=(0,0,1,1)T, as one checks on the first row, 2⋅0−1⋅0+1⋅1−1⋅1=0, and similarly on the others.
Step 5: assemble and check. Laying the chain down first,
P=−1−1−200001−11000011,detP=−2=0.
The identity AP=PJ holds column by column: Ax1=3x1, Av=3v+x1 (which is Bv=x1, verified above), Aw=3w and Az=−z. Hence P−1AP=J.
Powers. By the block-power formula, Jk=diag(J2(3)k,3k,(−1)k) with J2(3)k=(3k0k3k−13k), so Ak=PJkP−1 splits into three fixed matrices:
Here Π and Π′ are the projections onto G3 and G−1 along each other (Π2=Π, Π′2=Π′, Π+Π′=I) and M=BΠ is the nilpotent part, M2=0. Note that Π is not an orthogonal projection: A is not normal, and nothing forces the two generalized eigenspaces to be perpendicular.
Checking the formula. At k=1 it returns 3Π+M−Π′=A; for instance the (3,1) entry is 3⋅1+1⋅(−2)−(−1)=2, which is A31. At k=2, squaring A directly gives
A2=3−6−48−63−486613−8−6−6−129,
and the formula gives for the (3,1) entry 9⋅1+2⋅3⋅(−2)+(−1)=9−12−1=−4, which agrees.
which is what repeated squaring of A produces. The structure is visible in the arithmetic: the term −393660, carrying the factor k, already dominates 310=59049 at k=10, and the whole entry grows like −2k3k−1 rather than like 3k. That extra factor of k is the single 2×2 Jordan block, and nothing in the spectrum {3,−1} alone would have predicted it.
ApplicationOriginal
How fast a finite Markov chain forgets its start. For a finite irreducible aperiodic chain with transition matrix S and stationary distribution π, the n-step probabilities are the entries of Sn. When the eigenvalues are distinct one diagonalizes, and Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §6.6, p. 270) read off the geometric rate: with ϱ the largest modulus among the eigenvalues other than 1 (they write μ for it), there is a constant c with ∣pij(n)−πj∣≤cϱn. They then note that repeated eigenvalues block this route, and that the Jordan form takes its place: S=B−1JB in our notation (they write M for the Jordan matrix), with Sn=B−1JnB. Their Perron-Frobenius theorem (Grimmett & Stirzaker (2020), §6.6, Thm (1), p. 269) lists the eigenvalue 1 once for an irreducible aperiodic chain and puts every other eigenvalue strictly inside the unit circle, so the block at 1 is 1×1 and only the other blocks can be large. The block-power formula of this unit then says exactly what changes: an eigenvalue λ=1 whose largest block has size s contributes terms (in)λn−i with i≤s−1, so with s the largest block size over all eigenvalues other than 1 the clean bound acquires a polynomial factor and becomes ∣pij(n)−πj∣≤cns−1ϱn. Since ns−1ϱn≤c′τn for any τ strictly between ϱ and 1, geometric convergence survives, but the pre-asymptotic behaviour does not: for n up to about (s−1)/log(1/ϱ) the polynomial is still winning, and a chain with a large Jordan block mixes visibly more slowly than its second eigenvalue suggests. This is the discrete-time counterpart of the warning that the spectrum alone does not control ∥Ak∥.
Defective linear systems and the factor of t. For x˙=Ax with a repeated eigenvalue and too few eigenvectors, the missing solutions are supplied by the chain. Lee (Lee (2013), Appendix D, pp. 673-674) gives the 2×2 case explicitly: when A has the single eigenvalue λ with eigenvector Z0 and A−λI=0, there is a generalized eigenvector Z1 with (A−λI)Z1=Z0, and the second independent solution is eλt(tZ0+Z1). In the language above, Z0,Z1 is a Jordan chain of length two, and the solution is the second column of Pexp(tJ). The same bookkeeping in general: a chain x1,…,xs for λ carries the s solutions ur(t)=eλt∑i=0r−1(ti/i!)xr−i for r=1,…,s, polynomial in t of degree up to s−1. Differentiating ur and using Axr−i=λxr−i+xr−i−1 (with x0=0) reproduces u˙r=Aur term by term. Physically this is where secular terms come from — critically damped oscillators, resonant forcing, marginal stability in a linearized flow — and it is why a stability analysis must look at the Jordan structure of any eigenvalue sitting on the imaginary axis, not merely at its real part. See linear systems for the solution theory and the matrix exponential for computing exp(tA) in practice.
Perturbation theory near a degeneracy. In quantum mechanics the Hamiltonian is Hermitian and therefore never defective, but the matrices that govern open and driven systems are not. An effective non-Hermitian Hamiltonian describing a resonance, a gain-loss optical cavity or a damped coupled oscillator typically has isolated parameter values — exceptional points — where two eigenvalues and their eigenvectors coalesce, which is precisely a 2×2 Jordan block. The splitting of the two eigenvalues near such a point goes like ε in the detuning ε, not like ε, exactly as in the matrix (0ε10) of the remark above, and that square-root response is what makes exceptional points attractive as sensors. The Jordan form is the right language here because the whole phenomenon is the failure of diagonalizability at one point of parameter space.
A reason to stop short of step 4. Because steps 1 to 3 use only ranks and the characteristic polynomial, the Jordan form is often all one computes: it is what decides whether Ak or exp(tA) carries polynomial factors, what the minimal polynomial is (hence how many powers of A are genuinely independent, which bounds the length of a Krylov recurrence), and whether a projection or involution identity can hold. The eigenvectors, which cost far more, are needed only when the solution itself, and not its shape, is wanted.
How to reduce a quadratic or Hermitian form to its principal axes, decide its definiteness from the eigenvalue signs or from Sylvester's criterion on the leading principal minors, and bound the eigenvalues by the Rayleigh quotient.
Definition(Quadratic and Hermitian forms, definiteness)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.17, Eqs (8.105)-(8.112), pp. 288-290
A quadratic form on Rn is a function Q(x)=xTAx with A a real n×n matrix. Only the symmetric part of A contributes: writing A=21(A+AT)+21(A−AT), the antisymmetric part B satisfies xTBx=0 for every x, because xTBx is a scalar and equals its own transpose −xTBx. We therefore always take A=AT, at the cost of no generality; the symmetric matrix representing a given form is then unique. The complex analogue is the Hermitian formH(x)=x†Ax with A=A†, which is real-valued for every x∈Cn since H(x)=x†A†x=H(x).
A Hermitian (or real symmetric) A and its form are called
positive definite if x†Ax>0 for all x=0;
positive semidefinite if x†Ax≥0 for all x;
negative definite / negative semidefinite if −A is positive definite / semidefinite;
indefinite if x†Ax takes both signs.
For 1≤k≤n, the k-th leading principal minor of A is Dk=detA(k), where A(k) is the top-left k×k submatrix of A. A principal minor is the determinant of any submatrix obtained by keeping the same set of rows and columns.
Method(Principal axes)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.17, Eqs (8.110)-(8.111), p. 289 and §8.17.2, p. 292
Let A be real symmetric n×n with eigenvalues λ1,…,λn (with multiplicity) and a real orthonormal eigenbasis u1,…,un, as supplied by the spectral theorem, and let S=(u1⋯un), so STAS=Λ. Substituting x=Sx′ gives
Q(x)=xTAx=(x′)TΛx′=λ1x1′2+⋯+λnxn′2,
a sum of squares with no cross terms; the new axes, along the uk, are the principal axes of the form. The same holds for a Hermitian form with S unitary and xk′2 replaced by ∣xk′∣2.
The level set L={x:xTAx=1} is non-empty if and only if at least one λk is positive. If some λk>0, then L meets the line Ruk, on which the form takes the values λkt2, at t=±λk−1/2. If instead every λk≤0, then Q(x)=∑kλkxk′2≤0<1 for everyx, not merely on the axis lines, so L is empty.
Assume L=∅ and write p≥1, q, z for how many of the eigenvalues are respectively positive, negative and zero, so p+q+z=n. Then, with the axes taken along the uk:
z=0 and q=0 (so p=n): L is an ellipsoid, compact, with semi-axis λk−1/2 in the direction uk.
z=0 and q≥1: L is a hyperboloid, unbounded — for n=2 a hyperbola.
z≥1: L is a cylinder. It is invariant under translation by any multiple of a uk with λk=0, because such a translation does not change xk′'s contribution to Q; its cross-section in the span of the uk with λk=0 is the ellipsoid (q=0) or hyperboloid (q≥1) of those p+q eigenvalues.
The three cases are disjoint, and the indices k with λk≤0 are exactly those for which L misses the axis line Ruk.
Both names are meant in their degenerate senses as well. An "ellipsoid" in n=1 is the two points ±λ1−1/2 and in n=2 an ellipse; a "cylinder" is the product of a cross-section with the z-dimensional space spanned by the null directions, so for n=2 and λ=(1,0) it is the pair of parallel lines x1′=±1, and for z=n−1, p=1 a pair of parallel hyperplanes.
Idea of the proofOriginal
The spectral theorem supplies an orthogonal change of coordinates in which A becomes diagonal, and a diagonal matrix has no cross terms: in those coordinates the form is λ1x1′2+⋯+λnxn′2, one independent contribution per axis. Because the change of coordinates is orthogonal it changes no length and no angle, so every geometric statement about the level set may be read off in the new coordinates and transported back unchanged. The classification is then a one-variable question asked n times: a positive λk confines the level set in the uk direction, a negative one lets it run off to infinity, and a zero one lets it slide freely.
ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §8.17, Eqs (8.109)-(8.111), pp. 288-289 and §8.17.2, p. 292
Step 1: the change of coordinates. Since A is real symmetric, part 6 of the spectral theorem gives a real orthonormal eigenbasis u1,…,un and makes S=(u1⋯un) real orthogonal with STAS=Λ. Then S−1=ST, so x=Sx′ and x′=STx say the same thing, and x↦x′ is a bijection of Rn that preserves the norm, by part 4 of the characterisation in Unitary and orthogonal matrices. Substituting,
with xk′=ukTx the component of x along uk. The hypothesis that A is symmetric enters only through the spectral theorem — it is what makes the eigenbasis real and orthonormal, hence S orthogonal rather than merely invertible, and that is what keeps the geometry undistorted. For a Hermitian form the same lines with S unitary, † for T, give x†Ax=∑kλk∣xk′∣2, which is real because the λk are.
Step 2: the axis lines and non-emptiness. On the line Ruk a point is x=tuk, whose coordinates are x′=tek because ujTuk=δjk; so Q(tuk)=λkt2. The equation λkt2=1 has a real solution exactly when λk>0, and then exactly the two solutions t=±λk−1/2. This proves both that L meets Ruk precisely for λk>0, at those two parameter values, and that L=∅ as soon as some λk>0. Conversely if every λk≤0 then Q(x)=∑kλk(xk′)2≤0<1 for every x, so L=∅.
Step 3: adapted coordinates. Assume L=∅ and renumber the eigenvalues, carrying the columns of S with them, so that λ1,…,λp>0, then q negative ones, then z zeros, with p≥1 by step 2. In the x′ coordinates
L={x′:k≤p∑λk(xk′)2−p<k≤p+q∑∣λk∣(xk′)2=1}.
Because x↦x′ is a norm-preserving bijection, boundedness, compactness, distances from the origin and the lines contained in L are the same computed in either set of coordinates.
Step 4: z=0 and q=0, so p=n. Put λmin=minkλk>0. For x∈L,
λmin∥x′∥2≤k∑λk(xk′)2=1,so∥x∥=∥x′∥≤λmin−1/2,
and L is bounded. It is also closed, being the preimage of the closed set {1} under the continuous map Q, hence compact. It is the ellipsoid with semi-axes λk−1/2: the diagonal map y↦x′ with xk′=λk−1/2yk is a bijection of Rn carrying the unit sphere ∑kyk2=1 exactly onto L, since ∑kλk(xk′)2=∑kyk2; and by step 2 its intersection with the k-th axis is at distance λk−1/2 from the origin.
Step 5: z=0 and q≥1. Here λ1>0 and some λm<0. For each real t let x′(t) have xm′=t, x1′=((1+∣λm∣t2)/λ1)1/2 and all other coordinates zero; the square root is of a positive number, so this is defined for every t, and ∑kλk(xk′)2=(1+∣λm∣t2)−∣λm∣t2=1. So x′(t)∈L for all t, while ∥x′(t)∥≥∣t∣→∞: L is unbounded. For n=2 the defining equation is λ1(x1′)2−∣λ2∣(x2′)2=1, the standard equation of a hyperbola with respect to its own axes.
Step 6: z≥1. Let Z be the span of the uk with λk=0 and Y the span of the others, so Rn=Y⊕Z orthogonally. The expression ∑kλk(xk′)2 contains no term in a coordinate with λk=0, so Q(y+cuk)=Q(y) whenever λk=0: translating by any vector of Z moves no point of L off L, which is the stated invariance. Writing a general x as y+z with y∈Y, z∈Z gives Q(x)=Q(y), so
x∈L⟺y∈L∩Y,that isL=(L∩Y)+Z.
Now L∩Y is the level set at 1 of the restriction of the form to the (p+q)-dimensional space Y, on which the eigenvalues are exactly the non-zero λk with the same eigenvectors; steps 4 and 5 applied inside Y make it an ellipsoid when q=0 and a hyperboloid when q≥1. So L is that cross-section times Z.
Step 7: disjointness and the degenerate readings. The triple (p,q,z) counts the positive, negative and zero eigenvalues of A, so it is determined by A and exactly one of the three cases applies. For the names: in step 4 with n=1, L={±λ1−1/2}, a pair of points, and with n=2 it is an ellipse, both being the image of the unit sphere of the relevant dimension under a diagonal stretch. In step 6 with p+q=1 and that one eigenvalue positive, L∩Y is a pair of points at distance λ−1/2 from the origin, so L is a pair of parallel hyperplanes; for n=2 with λ=(1,0) the equation is (x1′)2=1, the two lines x1′=±1. □
Method(Eigenvalue test for definiteness)AdaptedAdapted from Hassani (2013), §6.4, Cor 6.4.10, p. 181; Riley, Hobson & Bence (2006), §8.17, the paragraph following Eq. (8.112), p. 290
Let A be Hermitian n×n (in particular, real symmetric). Then A is positive definite if and only if every eigenvalue of A is positive, and positive semidefinite if and only if every eigenvalue is ≥0.
The strict half is (Riley, Hobson & Bence (2006), §8.17, p. 290), which states that Q>0 for every x forces every eigenvalue to be positive. The corresponding result in (Hassani (2013), Cor 6.4.10, p. 181) is the semidefinite one, because "positive" there is Hassani's non-strict notion (Def 4.3.9, p. 117); his strict notion is written A>0. The semidefinite half is the same reading of the sum of squares Q=∑kλkxk′2 with ≥ in place of >, which is how (Riley, Hobson & Bence (2006), §8.17, p. 290) arrives at both terms. The consequence that trA>0 and detA>0 for a positive definite A is ours — it is immediate, since the trace is the sum and the determinant the product of the n eigenvalues — and is on neither cited page, so do not attribute it there. It is only a one-way test: see the traps below.
Idea of the proofOriginal
In the principal-axis coordinates the form is the weighted sum ∑kλk∣xk′∣2 with non-negative weights ∣xk′∣2, and such a sum is positive for every non-zero x exactly when every coefficient is positive: if one λj fails, the single eigenvector uj exhibits the failure, and if none fails, the smallest eigenvalue bounds the whole sum from below.
By the spectral theorem choose an orthonormal eigenbasis u1,…,un of A with Auk=λkuk and the λk real, let U be the unitary matrix of those columns and put x′=U†x, so that by step 1 of the previous proof
x†Ax=k=1∑nλk∣xk′∣2,∥x′∥=∥x∥.
The map x↦x′ is a bijection, so x=0 if and only if x′=0.
All eigenvalues positive implies positive definite. Let λmin=minkλk>0 and let x=0. Then
x†Ax=k∑λk∣xk′∣2≥λmink∑∣xk′∣2=λmin∥x∥2>0,
the first inequality because every ∣xk′∣2≥0 and the last because ∥x∥>0.
Positive definite implies all eigenvalues positive. Suppose some λj≤0. The eigenvector uj is non-zero and uj†Auj=λjuj†uj=λj≤0, so the defining inequality fails at x=uj and A is not positive definite.
The semidefinite half. If every λk≥0 then x†Ax=∑kλk∣xk′∣2≥0 for every x, each term being a product of two non-negative numbers. If some λj<0 then uj†Auj=λj<0. This is the same reading of the sum of squares with ≥ in place of >, exactly as (Riley, Hobson & Bence (2006), §8.17, p. 290) obtains both terms.
The trace and determinant. If A is positive definite then all λk>0, and part 5 of the spectral theorem gives trA=∑kλk>0 as a sum of positive numbers and detA=∏kλk>0 as a product of them. Neither converse holds: diag(−1,−1) has trace −2 but determinant +1, and diag(3,−1) has trace +2 with determinant −3, so each of the two conditions can hold while A fails to be positive definite. □
Let A be real symmetric n×n, with leading principal minors Dk=detA(k) for k=1,…,n. Then positive definiteness of the form xTAx is equivalent to positivity of the whole descending chain
Each of the n inequalities is needed, and no condition beyond them is.
Applying this to −A and using Dk(−A)=(−1)kDk(A): A is negative definite if and only if (−1)kDk>0 for k=1,…,n, i.e. the leading principal minors alternate in sign beginning with D1<0.
Idea of the proofOriginal
One direction is easy: restricting a positive definite form to the first k coordinates leaves it positive definite, so each A(k) has only positive eigenvalues and its determinant Dk is positive. The other direction is an induction in which A(n−1) is already known to be positive definite. That means the form is positive on an (n−1)-dimensional subspace, and two subspaces of Rn whose dimensions add to more than n must share a non-zero vector; so A cannot have two non-positive eigenvalues, because the plane they span would meet that subspace. At most one eigenvalue is therefore in doubt, and Dn=detA>0 settles its sign.
ProofOriginal
Throughout, A(k) is the top-left k×k block and Vk=span{e1,…,ek}. For x∈Vk, only the first k entries of x are non-zero, so
xTAx=i,j≤k∑aijxixj=x~TA(k)x~,x~=(x1,…,xk)T:
the form of A(k) is the form of A restricted to Vk. Each A(k) is symmetric, being a block of a symmetric matrix cut out by the same rows and columns.
Lemma (overlapping subspaces). If V,W⊆Rn are subspaces with dimV+dimW>n, then V∩W contains a non-zero vector. Indeed, let v1,…,va be a basis of V and w1,…,wb one of W, with a+b>n. More than n vectors in an n-dimensional space are linearly dependent (Vector spaces and bases), so ∑icivi+∑jdjwj=0 with the coefficients not all zero. Put x=∑icivi=−∑jdjwj, which lies in V and in W. If x were 0, independence of the vi would force every ci=0 and independence of the wj every dj=0, contradicting the choice; so x=0.
Necessity. Suppose A is positive definite. For 1≤k≤n and x~∈Rk non-zero, the vector x∈Vk with those first k entries is non-zero, so x~TA(k)x~=xTAx>0: each A(k) is positive definite. By the eigenvalue test, its eigenvalues are all positive, so Dk=detA(k), their product, is positive.
Sufficiency, by induction on n. Suppose D1,…,Dn>0.
Base case n=1. Here A=(a11) with a11=D1>0, and xTAx=a11x12>0 for x1=0.
Induction step. Let n≥2 and assume the criterion for n−1. The leading principal minors of A(n−1) are D1,…,Dn−1 — a leading block of a leading block is a leading block — so the induction hypothesis makes A(n−1) positive definite, and therefore
Suppose, for contradiction, that λ2≤0, hence also λ1≤0. Let W=span{u1,u2}, of dimension 2. For x=c1u1+c2u2∈W, orthonormality gives
xTAx=λ1c12+λ2c22≤0.
Since dimW+dimVn−1=2+(n−1)=n+1>n, the lemma supplies a non-zero x∈W∩Vn−1. That x has xTAx≤0 by the display and xTAx>0 by (∗) — a contradiction. Hence λ2>0, and so λ2,…,λn are all positive.
Finally Dn=detA=λ1λ2⋯λn by part 5 of the spectral theorem. The hypothesis Dn>0 and the positivity of λ2⋯λn force λ1>0. All eigenvalues are then positive and the eigenvalue test makes A positive definite, completing the induction.
Every one of the n inequalities is needed. Dropping the k-th leaves a false criterion, as these diagonal matrices show. For k≤n−1 take A diagonal with akk=ak+1,k+1=−1 and every other diagonal entry +1: then Dj=1 for j<k, Dk=−1, and Dj=(−1)(−1)=1 for j≥k+1, so only the k-th inequality fails, while A is not positive definite because ekTAek=−1. For k=n take A=diag(1,…,1,−1), which has Dj=1 for j<n and Dn=−1. And no further condition is needed, because the induction above derives positive definiteness from these n alone.
The negative definite form.A is negative definite exactly when −A is positive definite, which by the above is exactly Dk(−A)>0 for every k. Now (−A)(k)=−A(k) and scaling a k×k matrix by −1 scales its determinant by (−1)k (the determinant is homogeneous of degree k in the rows, Determinants), so Dk(−A)=(−1)kDk(A) and the condition reads (−1)kDk>0: D1<0, D2>0, D3<0, and so on. □
Let A be Hermitian n×n. For x∈Cn, x=0, the Rayleigh quotient of A at x is
R(x)=x†xx†Ax,
which is real, because the numerator and denominator are, and unchanged by rescaling x by any non-zero scalar.
Method(Rayleigh bounds)Original
Let A be Hermitian n×n with eigenvalues ordered λ1≤λ2≤⋯≤λn and Eλ the eigenspace of λ. Then λ1≤R(x)≤λn for every x=0; the lower bound is an equality exactly for x∈Eλ1∖{0} and the upper bound exactly for x∈Eλn∖{0}. Hence
λ1=x=0minR(x),λn=x=0maxR(x),
both attained. (Expand x in an orthonormal eigenbasis: R(x) is the weighted average of the λk with weights ∣xk′∣2/∥x∥2 summing to 1.)
Idea of the proofOriginal
In an orthonormal eigenbasis the Rayleigh quotient is literally a weighted average of the eigenvalues, the weight on λk being the fraction of ∥x∥2 that x puts along uk. An average of numbers lies between the smallest and the largest of them, and equals the smallest only if all the weight sits on the smallest — which is what "x lies in the eigenspace of λ1" says.
each summand being a product of two non-negative numbers, by the ordering of the eigenvalues. So λ1≤R(x)≤λn.
The equality cases. In the first display, equality forces every summand to vanish, so wk=0 — that is ck=0 — for every k with λk>λ1. Then x=∑k:λk=λ1ckuk, a combination of eigenvectors for λ1, so x∈Eλ1. Conversely if x∈Eλ1 and x=0 then Ax=λ1x and R(x)=λ1x†x/x†x=λ1. The same argument with λn gives the upper equality case. Note that the whole eigenspace appears, not one vector: if λ1 is degenerate, every non-zero vector of Eλ1 attains the minimum.
Attainment.u1∈Eλ1 and un∈Eλn are non-zero, so the bounds are attained and the infimum and supremum are a minimum and a maximum. □
Method(Stationary values of the quotient)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.17.1, Eqs (8.114)-(8.116) and the worked result for lambda(x), p. 291
Let A be Hermitian n×n. Then R is stationary at x=0 — its first variation vanishes for every variation Δx — exactly when Ax=R(x)x, that is, exactly when x is an eigenvector of A; the stationary value is then the eigenvalue belonging to x. Equivalently, the critical points of x†Ax on the sphere x†x=1 are the normalised eigenvectors of A, and the critical value at each is its eigenvalue.
Idea of the proofOriginal
Differentiate the quotient once. Because A is Hermitian, the two cross terms produced by varying x in x†Ax are complex conjugates of each other, so the first variation collapses to 2Re of a single inner product against the vector Ax−R(x)x. Demanding that this vanish for every variation is demanding that the vector itself vanish, which is the eigenvalue equation with R(x) as the eigenvalue. Nothing in the argument needs the eigenvalues to be known in advance; the eigenvalue equation is the output.
ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §8.17.1, Eqs (8.113)-(8.116) and the variational computation for lambda(x), pp. 290-291
Write N(x)=x†Ax and D(x)=x†x, so R=N/D on Cn∖{0}, and fix x=0 and a variation h∈Cn.
Step 1: the first variations of N and D. Expanding,
N(x+h)−N(x)=h†Ax+x†Ah+h†Ah.
The first two terms are complex conjugates of each other, because x†Ah=h†A†x=h†Ax using A=A† — this is the one place the Hermitian hypothesis is used, and without it the two terms do not combine. So
Now x†Ah=(Ax)†h, again because A is Hermitian, so with
r=Ax−R(x)x
the variation is δR=x†x2Re(r†h).
Step 3: stationarity is the eigenvalue equation. If r=0 then δR=0 for every h. Conversely, if δR=0 for every h, take h=r: then Re(r†r)=∥r∥2=0, so r=0. Hence R is stationary at x exactly when Ax=R(x)x, that is exactly when x is an eigenvector, and then the eigenvalue is R(x), which is the stationary value. Conversely any eigenvector Ax=λx has R(x)=λ and r=0, so it is stationary; the quantifier "for every variation" is essential, since Re(r†h) vanishes for the h orthogonal to r however large r is.
Step 4: the version on the sphere.R is unchanged by rescaling, so R(x)=N(x) whenever ∥x∥=1, and δR vanishes on the radial variation h=x for every x: indeed r†x=x†Ax−R(x)x†x=0 by the definition of R. A variation tangent to the sphere at x is one with Re(x†h)=0, obtained by differentiating ∥x∥2=1; for such h, step 2 with D(x)=1 gives
So the critical points of N on the sphere are exactly the unit vectors at which Re(r†h) vanishes for all tangent h. That is the same as vanishing for all h∈Cn, because any h splits as h=ht+cx with c=Re(x†h) real and ht=h−cx tangent — check Re(x†ht)=Re(x†h)−c=0 — while Re(r†x)=0 always. Hence it is the same as r=0. The critical points on the sphere are therefore the normalised eigenvectors, with critical values N(x)=R(x)=λ. Riley, Hobson and Bence reach the same conclusion with a Lagrange multiplier imposing x†x=1, the multiplier turning out to be the eigenvalue. □
Method(Courant-Fischer min-max)Original
With A and the ordering λ1≤⋯≤λn as above, for every k=1,…,n
the inner extrema being attained and the outer ones attained at V=span{u1,…,uk} and W=span{uk,…,un} respectively. The cases k=1 and k=n are the Rayleigh bounds above. Restricting the minimum to a smaller family of subspaces can only raise it, which is what makes this a practical source of upper bounds for λk.
Idea of the proofOriginal
Two inequalities, in opposite directions, pin λk down. One is a calculation: on the particular subspace spanned by the first k eigenvectors the quotient never exceeds λk and reaches it, so the minimum over all k-dimensional subspaces is at most λk. The other is a counting argument: any k-dimensional subspace must meet the (n−k+1)-dimensional span of the lastn−k+1 eigenvectors, because their dimensions add to n+1; on that span the quotient is at least λk, so the maximum over any k-dimensional subspace is at least λk. The max-min formula is the min-max formula read for −A.
ProofOriginal
Keep the orthonormal eigenbasis u1,…,un with λ1≤⋯≤λn and the weighted-average form R(x)=∑kwkλk of the Rayleigh proof. Two consequences of it will be used repeatedly: if x=0 lies in Vk=span{u1,…,uk} then the weights wj vanish for j>k, so
R(x)=j≤k∑wjλj≤λk;(A)
and if x=0 lies in Wk=span{uk,…,un} then the weights vanish for j<k, so
R(x)=j≥k∑wjλj≥λk.(B)
The inner extrema are attained. Let V be a subspace with dimV=d≥1. R is continuous on Cn∖{0}, the unit sphere ΣV={x∈V:∥x∥=1} is closed and bounded in the finite-dimensional space V, hence compact, and R is unchanged by rescaling, so the extrema of R over the non-zero vectors of V are the extrema over ΣV and are attained there.
Upper bound for the min-max. Take V=Vk, which has dimension k because the uj are independent. By (A), maxx∈Vk∖{0}R(x)≤λk, and uk∈Vk gives R(uk)=λk, so the maximum over Vk equals λk. Hence the minimum over all k-dimensional subspaces is at most λk.
Lower bound for the min-max. Let V be any subspace of dimension k. Then dimV+dimWk=k+(n−k+1)=n+1>n, so by the overlapping-subspaces lemma of the Sylvester proof — whose argument is unchanged over C — there is a non-zero y∈V∩Wk. By (B), R(y)≥λk, and therefore
x∈V∖{0}maxR(x)≥R(y)≥λk.
As V was arbitrary, the minimum over all k-dimensional subspaces is at least λk. With the previous paragraph it equals λk and is attained at V=Vk.
The max-min form. Apply what is proved to −A, which is Hermitian with eigenvalues −λn≤⋯≤−λ1, so that the j-th smallest eigenvalue of −A is −λn−j+1, and with R−A=−RA. For j=n−k+1 the min-max formula for −A reads
using max(−f)=−minf twice. Multiplying by −1 gives the stated max-min formula over subspaces W of dimension n−k+1. The optimal subspace transports too: for −A the optimum was the span of its first n−k+1 eigenvectors, which are un,un−1,…,uk, that is W=span{uk,…,un}.
The cases k=1 and k=n. For k=1 the only 1-dimensional subspaces are the lines Cx, on each of which R is constant, so the min-max reduces to minx=0R(x)=λ1; for k=n the only n-dimensional subspace is Cn itself and the formula reads maxx=0R(x)=λn. These are the Rayleigh bounds.
Why restricting the family gives upper bounds. A minimum over a subfamily of the k-dimensional subspaces is at least the minimum over all of them, which is λk. So evaluating maxx∈VR(x) for any single conveniently chosen V of dimension k — a span of k trial vectors, say — delivers a number ≥λk, and equality holds exactly when that V is optimal. □
Method(Analysing a quadratic form)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.17 and §8.17.1, pp. 288-291
Given a form written as xTMx with M any real square matrix, or as an explicit polynomial such as x12+4x1x2−x22:
Symmetrise. Replace M by A=21(M+MT). From a polynomial, set aii to the coefficient of xi2 and aij=aji to half the coefficient of xixj. Skip this and you will diagonalize the wrong matrix: the form determines only the symmetric part.
If you only need definiteness, apply Sylvester's criterion: compute D1,…,Dn. All positive means positive definite; alternating from D1<0 means negative definite. Any other pattern with all Dk=0 means indefinite. If some Dk=0 the criterion is silent — go to step 3.
Otherwise diagonalize. Compute the eigenvalues of A from det(A−λI)=0 for small n, or with a symmetric eigensolver (numpy.linalg.eigh) for larger n. The signs of the eigenvalues give the definiteness class, with no exceptional cases; the triple (number positive, number negative, number zero) is the inertia of the form.
For the principal axes, compute an orthonormal eigenbasis, assemble S with those columns, and read off Q=∑kλkxk′2 where x′=STx. Order the eigenvalues as you like, but keep the columns of S in the matching order.
For the geometry of a level setxTAx=c with c>0: first test non-emptiness — the set is non-empty if and only if some λk>0 (for A=diag(−1,−1) and c=1 it is empty, and there is no geometry to describe). Given that, the axis directions are the uk; for each λk>0 the set meets the line Ruk at distance (c/λk)1/2 from the origin, and for each λk≤0 it does not meet that line at all. The set is unbounded precisely when some λk≤0; it contains, through each of its points, the whole line in the direction uk — equivalently, it is invariant under translation along uk — precisely when λk=0. These are two different conditions: a hyperbola (A=diag(1,−1), c=1) is unbounded along both of its axis directions while containing no line at all, since no eigenvalue vanishes. When every λk>0 it is an ellipsoid with those semi-axes; otherwise classify it by the counts (p,q,z) as in "Principal axes" — the set for level c is the set for level 1 scaled by c1/2, so every semi-axis is multiplied by c1/2 and the classification is unchanged.
For quick eigenvalue bounds, evaluate R at any convenient trial vectors: each value is an upper bound for λ1 and a lower bound for λn, and the bound is good when the trial vector is close to the eigenvector, because R is stationary there.
The recipe needs A Hermitian throughout; steps 2-6 are false without it. For a Hermitian complex A use x†Ax, A† and a unitary S in place of the real transposes.
NoteOriginal
A quadratic form is the second-order part of any smooth function near a critical point, so its definiteness decides whether a potential has a minimum, an equilibrium is stable or an energy is bounded below. Diagonalizing it chooses coordinates in which those questions become n questions about single numbers.
Remark(Traps)Original
Sylvester's criterion does not test semidefiniteness.A=diag(0,−1) has D1=0 and D2=0, both ≥0, yet A is not positive semidefinite. Positive semidefiniteness is equivalent to all principal minors (not just the leading ones) being ≥0, or simply to all eigenvalues being ≥0. For definiteness the leading minors suffice; for semidefiniteness they do not.
A positive determinant is not positive definiteness.diag(−1,−1) has determinant 1. Only the full chain D1,…,Dn settles the question.
The equality case needs the whole eigenspace. If λ1 is degenerate, R(x)=λ1 for every non-zero x∈Eλ1, not only for one eigenvector.
Min-max is over subspaces, not vectors. In the Courant-Fischer formula the outer minimum runs over k-dimensional subspaces; replacing it by a minimum over vectors would give λ1 for every k.
Sylvester's criterion for Hermitian A. Stated above for real symmetric A, as (Gelfand & Fomin (1963), §30, p. 127) gives it. It holds word for word for Hermitian A: each leading submatrix A(k) is itself Hermitian, so Dk=detA(k) is real, and the criterion reads the same. That extension is not on the cited page, so do not attribute it there.
Conventions: the factor 21 doubles the eigenvalues, it does not halve them. Fix a form Q. Writing it our way, Q(x)=xTAx, makes Ahalf the Hessian of Q. Texts that write Q(x)=21xTA~x — common in mechanics and optimisation, so that A~ is the Hessian — have A~=2A, whose eigenvalues are twice ours. For the polynomial Q=x12+4x1x2−x22 used to illustrate step 1 of the method, our A=(122−1) has eigenvalues ±5, while the Hessian A~=(244−2) has ±25. Going the other way — holding the matrix fixed and reading it in the half convention — halves the form and hence the eigenvalues; what the factor does depends on which of the two you hold fixed, so say which. Either way the definiteness class, the inertia and the principal axes are unchanged, since the two matrices differ by a positive scalar; only the numbers λk and therefore the semi-axis lengths λk−1/2 change.
Ordering. Order the eigenvalues increasingly, as here, or decreasingly, as is common in numerical linear algebra — but then the min-max formula swaps its two halves.
Example(A form in three variables, start to finish)Original
Suppose a calculation has produced the form Q(x)=xTMx with
M=3−2−24310−35,
which is not symmetric. Follow the method.
Step 1: symmetrise.A=21(M+MT) has entries aii=Mii and aij=21(Mij+Mji):
Reading the polynomial back off A is the check: the diagonal gives the squares and each off-diagonal entry is half the coefficient of the corresponding cross term, so a12=1 goes with 2x1x2.
Step 2: definiteness by Sylvester. The leading principal minors are
each verified directly: A(1,−1,0)T=(2,−2,0)T, A(1,1,1)T=(3,3,3)T, A(1,1,−2)T=(6,6,−12)T. The three are mutually orthogonal, as the spectral theorem requires for distinct eigenvalues, and all positive, confirming step 2 independently.
Step 4: the principal axes. With S=(u1u2u3) and x′=STx,
Step 5: the level set Q=1. The inertia is (p,q,z)=(3,0,0), so this is an ellipsoid, with semi-axes
λ1−1/2=21≈0.7071 along u1,31≈0.5774 along u2,61≈0.4082 along u3.
A spot check on the first: the point u1/2=(1,−1,0)T/2 has Q=41(3+3−2)=1, and its distance from the origin is 2/2=1/2. For the level Q=4 every semi-axis is multiplied by 4=2.
Step 6: Rayleigh bounds from trial vectors. Every trial vector gives 2≤R(x)≤6. With x=(2,−1,0)T, a rough guess at u1,
R(x)=4+13⋅4+3⋅1+2⋅2⋅(−1)=512+3−4=511=2.2,
an upper bound for λ1=2, high by 10% although the trial vector points 18.4∘ away from u1. Stationarity sharpens this: perturbing the true eigenvector to x=(1,−1,δ)T gives
R(x)=2+δ24+5δ2=2+23δ2+O(δ4),
so at δ=0.1 — a 7% tilt of the vector — the quotient is 135/67=2.01493, high by only 0.75%. The error in the eigenvalue is quadratic in the error of the eigenvector, which is exactly the content of the stationarity statement and the reason trial-vector bounds are useful.
Step 7: Courant–Fischer for the middle eigenvalue. Take the 2-dimensional trial subspace V=span{e1,e2}. For x∈V the quotient is the Rayleigh quotient of the block A(2)=(3113), whose eigenvalues are 2 and 4, so
x∈V∖{0}maxR(x)=4≥λ2=3,
an upper bound for λ2 from a subspace chosen without any eigenvector information. The optimal choice V=span{u1,u2} does give 3: in the orthonormal basis {u1,u2} of that plane, A is represented by diag(2,3), whose largest eigenvalue is 3.
Step 8: a one-entry change makes it indefinite. Replace the entry a33=5 by 0, i.e. drop the 5x32 term. Now D1=3, D2=8 and
D3=3(0−1)−1(0−1)+(−1)(−1+3)=−3+1−2=−4,
so the pattern is (+,+,−) with no zero: by step 2 of the method the form is indefinite, and Sylvester tells us so without any eigenvalue computation. For confirmation, the characteristic polynomial is −(λ−2)(λ2−4λ−2), with roots 2 and 2±6≈4.4495,−0.4495 — inertia (2,1,0). The level set is then a hyperboloid, unbounded along u for the negative eigenvalue, and Q takes both signs: Q(e3)=0 for the modified form while Q(0,0,1) was 5 before, and Q(1,1,1)=3+3+0+2−2−2=4>0 while Q(1,1,3)=3+3+0+2−6−6=−4<0.
ApplicationOriginal
The second-derivative test, and why stability is an eigenvalue question. Let f be twice continuously differentiable near a critical point a, so ∇f(a)=0. Taylor's theorem gives
f(a+h)=f(a)+21hTHh+o(∥h∥2),Hij=∂xi∂xj∂2f(a),
and H is symmetric, so the leading behaviour near a is governed by a quadratic form and this unit settles it: H positive definite makes a a strict local minimum, H negative definite a strict local maximum, and H indefinite a saddle, with directions of increase and decrease given by eigenvectors for eigenvalues of either sign. Sylvester's criterion is what one actually computes, since it needs only n determinants rather than the full spectrum; for n=2 it is the familiar test fxx>0 and fxxfyy−fxy2>0. In mechanics the same statement is the stability of an equilibrium: with f the potential energy, positive definiteness of H is the condition for a restoring force in every direction, and the eigenvalues of H (divided by the masses) are the squared frequencies of the small oscillations about it — see normal modes. A semidefinite H leaves the question open at second order, which is where the trap about Sylvester's criterion and semidefiniteness bites: f(x,y)=x2+y3 and f(x,y)=x2+y4 have the same Hessian at the origin and only the second has a minimum there.
Rayleigh–Ritz: a variational upper bound on a ground-state energy. In quantum mechanics the energy levels of a system are the eigenvalues λ1≤λ2≤⋯ of a Hermitian Hamiltonian, and the Rayleigh bound says the ground-state energy is
λ1=ψ=0minψ†ψψ†Hψ.
So any trial state gives an upper bound on the ground-state energy — never an underestimate, which is what makes the method trustworthy — and because the quotient is stationary at the true eigenvector, a trial state wrong by a relative error ε gives an energy wrong only by order ε2. That is the whole of the variational method: pick a family of trial states depending on parameters, minimise the quotient over the parameters, and report the result as an upper bound. Courant–Fischer extends it to excited states: taking a k-dimensional trial space and computing the largest eigenvalue of the k×k matrix of H restricted to it gives an upper bound for λk, which is exactly what step 7 of the example does with k=2. Enlarging the trial space can only lower each bound, and this is why the method converges as the basis is refined. The same mathematics appears as Rayleigh's method for the fundamental frequency of a vibrating structure and as the finite element method for eigenvalue problems.
Principal axes of a physical tensor. Many physical quantities are symmetric matrices whose quadratic form has a direct meaning, and diagonalizing them is the first step of every calculation. The inertia tensor I of a rigid body has ωTIω=2T twice the rotational kinetic energy; it is positive definite for a genuine three-dimensional body, so its eigenvectors are the three mutually perpendicular principal axes about which rotation produces no wobble, its eigenvalues the principal moments, and the level set ωTIω=1 is the ellipsoid of inertia with semi-axes Ik−1/2. The Cauchy stress tensor is symmetric, and its eigenvalues are the principal stresses with eigenvectors the directions of pure tension or compression — a negative eigenvalue is compression, so the inertia of the form is the statement of what kind of loading the material is under. The strain tensor, the electric susceptibility of an anisotropic crystal and the diffusion tensor of an anisotropic medium all work the same way: the eigenvectors give the material's own axes, and the numbers one quotes are the eigenvalues.
States the spectral theory of the pencil Ax = lambda Bx for Hermitian A and positive definite B, and applies it to the small oscillations of a conservative system, whose normal frequencies and normal modes are the eigenvalues and B-orthonormal eigenvectors of the pencil formed by the potential and kinetic energy matrices.
Definition(The generalized eigenvalue problem)AdaptedAdapted from Riley, Hobson & Bence (2006), §9.1, Eqs (9.9)-(9.10), p. 319
Let A and B be Hermitian n×n matrices with Bpositive definite. A number λ∈C is an eigenvalue of the pencil(A,B) if
Ax=λBxfor some x=0,
and such an x is a corresponding eigenvector. The equation Ax=λBx has a non-zero solution exactly when A−λB is singular, that is exactly when
det(A−λB)=0,
the characteristic equation of the pencil. That equivalence is just the singularity criterion and needs nothing of B; what invertibility of B contributes is that the left-hand side is a polynomial in λ of degree exactly n, so that the pencil has n eigenvalues and not fewer (part 1 of the next statement). Two vectors x,y are B-orthogonal if y†Bx=0; because B is Hermitian positive definite, (x,y)B=y†Bx is an inner product on Cn, and B-orthogonality is orthogonality in it. Taking B=I recovers the ordinary eigenproblem.
Application(Spectral theorem for a definite pencil)Original
Let A be Hermitian and B Hermitian positive definite, both n×n. Then:
det(A−λB) is a polynomial in λ of degree exactly n, with leading coefficient (−1)ndetB=0; so the pencil has exactly n eigenvalues counted with multiplicity, and all of them are real.
There is a basis x1,…,xn of Cn with
Axk=λkBxk,(xj)†Bxk=δjk,(xj)†Axk=λkδjk.
The eigenvectors are thus simultaneously B-orthonormal and A-orthogonal. If A and B are real symmetric, the xk may be taken real.
3. The problem is equivalent to an ordinary Hermitian eigenproblem: with B1/2 the Hermitian positive definite square root of B, the matrix C=B−1/2AB−1/2 is Hermitian, and Ax=λBx holds if and only if Cy=λy with y=B1/2x. In particular the pencil and C have the same eigenvalues with the same multiplicities, and x↦B1/2x carries a B-orthonormal eigenbasis of the pencil to an orthonormal eigenbasis of C. (A Cholesky factor B=LL† serves equally well, with C=L−1AL−†.)
4. Quotient bounds. With λ1≤⋯≤λn,
λ1≤x†Bxx†Ax≤λnfor every x=0,
both bounds attained, and the quotient is stationary at x exactly when x is an eigenvector of the pencil, with stationary value the corresponding eigenvalue. A is positive definite if and only if every λk>0, and positive semidefinite if and only if every λk≥0.
Idea of the proofOriginal
A positive definite B has a positive definite square root, so the substitution y=B1/2x is available and it is invertible. Under it the pencil equation Ax=λBx becomes the ordinary eigenvalue equation Cy=λy for the Hermitian matrix C=B−1/2AB−1/2, and the B-inner product (u,v)B=v†Bu becomes the ordinary one, since (u,v)B=(B1/2v)†(B1/2u). So the whole theorem — reality of the spectrum, the count, the double orthogonality, the quotient bounds — is the Hermitian spectral theorem and the Rayleigh bounds transported through one change of variable. Definiteness of B is what makes the square root, and hence the whole argument, exist; mere invertibility is not enough.
It is Hermitian, its eigenvalues βk are positive so it is positive definite, (B1/2)2=VΔV†=B, and it is invertible with B−1/2=Vdiag(βk−1/2)V†, also Hermitian. It is the only Hermitian positive definite square root: if X is one, then X has an orthonormal eigenbasis with positive eigenvalues γ, and X2=B makes each such eigenvector an eigenvector of B with eigenvalue γ2; since t↦t2 is injective on the positive reals, X must act as β on the eigenspace Eβ(B) for every β, which determines X completely.
Step 2: part 3, the reduction. Put C=B−1/2AB−1/2. It is Hermitian, since B−1/2 and A are: C†=B−1/2A†B−1/2=C. Because B1/2CB1/2=A and B1/2(λI)B1/2=λB,
A−λB=B1/2(C−λI)B1/2.(∗)
Hence, with y=B1/2x — a bijection of Cn carrying 0 to 0 and nothing else there —
Ax=λBx⟺B1/2(C−λI)B1/2x=0⟺(C−λI)y=0⟺Cy=λy,
the middle step because B1/2 is invertible. For the Cholesky variant, if B=LL† with L invertible then A−λB=L(L−1AL−†−λI)L† by the same algebra, and C′=L−1AL−† is Hermitian; everything below goes through with L† in place of B1/2.
Step 3: part 1. Taking determinants in (∗) and using multiplicativity (Determinants),
det(A−λB)=(detB1/2)2det(C−λI)=detB⋅det(C−λI),
the last equality because detB1/2detB1/2=det((B1/2)2)=detB. Now det(C−λI) is a polynomial in λ of degree n with leading coefficient (−1)n, so det(A−λB) is a polynomial of degree exactly n with leading coefficient (−1)ndetB, which is non-zero because B is positive definite and hence detB>0. A polynomial of degree n over C has exactly n roots with multiplicity, and by step 2 those roots are the eigenvalues of the pencil; they are the eigenvalues of the Hermitian matrix C, with the same multiplicities, hence all real.
Step 4: part 2. By the spectral theorem C has an orthonormal eigenbasis y1,…,yn, Cyk=λkyk. Put xk=B−1/2yk. These are a basis of Cn, being the images of a basis under an invertible matrix, and by step 2 they satisfy Axk=λkBxk. Moreover, using that B−1/2 is Hermitian,
So the xk are B-orthonormal and A-orthogonal at once. If A and B are real symmetric, then V may be taken real orthogonal (part 6 of the spectral theorem), hence B±1/2 and C are real symmetric, the yk may be taken real, and so may the xk.
Step 5: part 4, the bounds. For x=0 put y=B1/2x=0. Then
x†Ax=y†B−1/2AB−1/2y=y†Cy,x†Bx=y†B−1/2BB−1/2y=y†y,
so the generalized quotient at x equals the ordinary Rayleigh quotient RC(y). The denominator x†Bx=∥y∥2 is positive for x=0, which is where positive definiteness of B is needed to make the quotient defined at all. The Rayleigh bounds of Quadratic forms and the Rayleigh quotient now give λ1≤RC(y)≤λn with both bounds attained, at y=y1 and y=yn, that is at x=x1 and x=xn.
Step 6: part 4, stationarity and definiteness. The map x↦y=B1/2x is a linear bijection, so for any variation h at x the corresponding variation at y is B1/2h, and h↦B1/2h runs over all of Cn as h does. Hence the first variation of the generalized quotient vanishes for every h exactly when the first variation of RC vanishes for every variation at y, which by the stationary-value statement of Quadratic forms and the Rayleigh quotient happens exactly when Cy=RC(y)y, that is exactly when x is an eigenvector of the pencil, with the stationary value its eigenvalue. Finally, the first display of step 5 shows that x†Ax is positive at every non-zero x exactly when y†Cy is positive at every non-zero y, since x↦y is a bijection fixing only 0; and by the eigenvalue test that holds exactly when every eigenvalue of C — that is, every λk — is positive. The same sentence with ≥ gives the semidefinite case. □
Definition(Small oscillations)AdaptedAdapted from Riley, Hobson & Bence (2006), Ch. 9 introduction, Eqs (9.1)-(9.2), pp. 316-317, and §9.1, Eq. (9.7), p. 318
Let a conservative system be described by generalized coordinates q=(q1,…,qn)T, chosen to vanish at an equilibrium configuration, with a potential V(q) that does not depend on the velocities and with V(0)=0. For small q and q˙ the kinetic and potential energies are quadratic forms
T=21q˙TMq˙,V=21qTKq+O(∥q∥3),
where M and K are constant real symmetric n×n matrices, Kij=∂2V(0)/∂qi∂qj the Hessian of V at the equilibrium, and the linear term is absent because ∂V(0)/∂qi=0 at an equilibrium. M is the mass or kinetic-energy matrix and is positive definite, since T>0 whenever q˙=0; K is the stiffness or potential-energy matrix. Keeping only these quadratic terms, Lagrange's equations dtd∂q˙i∂T−∂qi∂T+∂qi∂V=0 give the linear system
Mq¨+Kq=0.
This is Eq. (9.7) of Riley, Hobson and Bence, who reach it instead from conservation of energy, dtd(T+V)=0, and defer the Lagrangian derivation to their chapter 22; the route through Lagrange's equations used here is ours.
The middle term disappears because M has been frozen at its equilibrium value: in general T=21q˙TM(q)q˙, so that
∂qi∂T=21q˙T(∂iM(q))q˙,
which carries two factors of q˙ and, since ∂iM(0) need not vanish, no factor of q: it is of second order in the small quantities q,q˙, one order beyond the first-order terms Mq¨ and Kq that are kept, and is therefore dropped in the linearization, so that only dtd(Mq˙)=Mq¨ survives from the kinetic part. Counted in the energies rather than in the equations of motion, the same statement is that the neglected part of T itself, 21q˙T(M(q)−M(0))q˙=O(∥q∥∥q˙∥2), is third order against the second-order forms kept. (The example below has exactly such a q-dependent T, through a factor cos(q1−q2), which step 1 there evaluates at q=0; because that factor is even in q1−q2 it is the special case ∂iM(0)=0, in which ∂T/∂qi is third order and not merely second.) The displayed system is the equation of motion of the system linearized about the equilibrium; everything below is a statement about it, not about the exact nonlinear motion.
Application(Normal frequencies and normal modes)AdaptedAdapted from Riley, Hobson & Bence (2006), §9.1, Eqs (9.8)-(9.13), pp. 318-321 and §9.3, Eqs (9.16), (9.19), pp. 327-328
With M,K as above, substituting q(t)=xe−iωt (or xcosωt) into Mq¨+Kq=0 gives
(K−ω2M)x=0,
the pencil (K,M) with λ=ω2. By the previous statement there are n real values λ1≤⋯≤λn and real M-orthonormal vectors x1,…,xn with Kxk=λkMxk. Then:
The normal coordinatesηk(t)=(xk)TMq(t) satisfy the uncoupled normal equations
η¨k+λkηk=0,k=1,…,n,
and q(t)=∑kηk(t)xk. In these coordinates T=21∑kη˙k2 and V=21∑kλkηk2.
2. Hence every solution is
q(t)=k∑(akcosωkt+bksinωkt)xkover the k with λk=ωk2>0,
plus a term (ak+bkt)xk for each k with λk=0 and a term (akeγkt+bke−γkt)xk with γk=−λk for each k with λk<0. The positive ωk=λk are the normal frequencies (eigenfrequencies) and the motions q(t)=ηk(t)xk are the normal modes.
3. Every solution is bounded for all t, for every choice of initial data, if and only if every λk>0, i.e. if and only if K is positive definite; this is the case when the equilibrium is a strict local minimum of V to second order. A negative λk gives exponential growth along xk (unstable equilibrium), and a zero λk gives a drift linear in t along xk.
4. Orthogonality and bounds.(xj)TMxk=δjk and (xj)TKxk=λkδjk: modes of distinct frequencies are orthogonal with respect to both energies, and within a degenerate eigenspace the basis can be chosen so that the two displayed relations hold for every j=k as well, degenerate or not — that is what part 2 of the pencil theorem delivers and what step 5 of the method below constructs. In neither case are the modes in general orthogonal with respect to the plain dot product (xj)Txk. For any trial vector x=0,
λ1≤xTMxxTKx≤λn,
with λ1=ω12 and λn=ωn2 whenever those eigenvalues are positive — parts 2 and 3 allow λk≤0, and no frequency ωk is defined for such a k. So at a stable equilibrium a guessed mode shape gives an upper bound for the lowest frequency and a lower bound for the highest.
Idea of the proofOriginal
Everything follows from the one change of variables ηk=(xk)TMq, which reads off the amount of the k-th mode present in the configuration q. Because the modes are M-orthonormal, this change of variables inverts itself as q=∑kηkxk; and because they are also K-orthogonal, applying it to Mq¨+Kq=0 kills every cross term and leaves n separate equations η¨k+λkηk=0. A single scalar equation of that form is the harmonic oscillator, whose solution everyone knows and whose character depends on the sign of λk alone: oscillation, drift, or exponential growth. So the behaviour of the whole coupled system is read off from n signs.
The eigenvalue problem. Substituting q(t)=xe−iωt into Mq¨+Kq=0 and using q¨=−ω2xe−iωt gives (−ω2M+K)xe−iωt=0; the scalar factor never vanishes, so (K−ω2M)x=0. With λ=ω2 this is the pencil (K,M), and M is positive definite by hypothesis while both are real symmetric, so the previous statement applies: there are n real λ1≤⋯≤λn and real M-orthonormal x1,…,xn with Kxk=λkMxk, (xj)TMxk=δjk and (xj)TKxk=λkδjk. The trial form xcosωt leads to the same equation, since it too has q¨=−ω2q.
Part 1, the change of variables inverts itself. The xk are a basis, so any q can be written q=∑kckxk; then
(xj)TMq=k∑ck(xj)TMxk=k∑ckδjk=cj,
using M-orthonormality. So cj=ηj and q(t)=∑kηk(t)xk, as claimed. This is the only place the normalisation (xk)TMxk=1 is used; with a different normalisation the formula would carry a factor.
Part 1, the normal equations. Differentiating ηk=(xk)TMq twice in t — the xk and M are constant — and then using the equation of motion in the form Mq¨=−Kq,
η¨k=(xk)TMq¨=−(xk)TKq=−j∑ηj(xk)TKxj=−λkηk,
the last step by K-orthogonality, which is what removes every term but j=k. So η¨k+λkηk=0, and the n equations are uncoupled.
Part 1, the energies. Substituting q=∑kηkxk and q˙=∑kη˙kxk,
where in the second computation the double sum collapses by K-orthogonality to the terms j=k, each contributing λkηk2.
Part 2. The scalar equation η¨+λη=0 is a second-order linear equation with constant coefficients, so its solution space is two-dimensional — a solution is determined by the two numbers η(0),η˙(0), and the general statement for linear equations of order n belongs to the theory of linear differential equations, which comes later in the handbook and is not needed as a prerequisite here. A two-parameter family of solutions is therefore all of them. For λ=ω2>0, cosωt and sinωt are solutions and are independent, giving η=acosωt+bsinωt. For λ=0 the equation is η¨=0, so η=a+bt. For λ<0, put γ=−λ>0; then e±γt are independent solutions and η=aeγt+be−γt. Assembling q=∑kηkxk gives the stated general solution, and the motion with ηj≡0 for j=k is q(t)=ηk(t)xk: the whole system moves in the fixed pattern xk, with every coordinate oscillating at the single frequency ωk=λk and passing through zero together. That is the normal mode.
Part 3. Suppose every λk>0. Then ∣ηk(t)∣=∣akcosωkt+bksinωkt∣≤(ak2+bk2)1/2 for all t, by the Cauchy-Schwarz inequality in R2 applied to (ak,bk) and (cosωkt,sinωkt), so
∥q(t)∥≤k∑∣ηk(t)∣∥xk∥≤k∑(ak2+bk2)1/2∥xk∥,
a bound independent of t: every solution is bounded. Conversely suppose λm≤0 for some m, and choose the initial data ηm(0)=0, η˙m(0)=1 and ηk(0)=η˙k(0)=0 for k=m, which is legitimate because the map from (q(0),q˙(0)) to (ηk(0),η˙k(0)) is a bijection. Then q(t)=ηm(t)xm with ηm(t)=t if λm=0 and ηm(t)=γm−1sinh(γmt) if λm<0, both unbounded, and xm=0. So boundedness of every solution is equivalent to all λk>0, which by part 4 of the pencil theorem applied to (K,M) is equivalent to K being positive definite. Since V=21qTKq to second order, K positive definite says exactly that the quadratic part of V is positive in every direction, i.e. the equilibrium is a strict local minimum of V to second order. The two remaining cases are the computations just made: λm<0 gives eγmt growth along xm, and λm=0 gives drift linear in t along xm.
Part 4. The orthogonality relations are part 2 of the pencil theorem, restated for (K,M). They do not imply (xj)Txk=0: M-orthogonality is orthogonality in the inner product (x,y)M=yTMx, and unless M is a multiple of the identity that is a different relation — the example below exhibits two modes whose unnormalised shapes u=(1,3)T and v=(1,−3)T satisfy uTMv=0 while uTv=−2=0; rescaling to the M-orthonormal x1,x2 changes only the size of that dot product, to −2/(ml224)=−0.408248/(ml2), not its vanishing. The bound is part 4 of the pencil theorem with A=K, B=M, and λk=ωk2 for those λk that are positive. □
Method(Finding the normal modes of a system)AdaptedAdapted from Riley, Hobson & Bence (2006), Ch. 9 introduction and §9.1, pp. 316-321
Choose coordinatesqi that vanish at the equilibrium, as many as the system has degrees of freedom and no more. A redundant coordinate makes M singular and the pencil ill-posed.
Expand the energies to second order. Write T in terms of the q˙i, keeping terms quadratic in q˙, to get M; expand V about q=0 to second order to get K. Both come out symmetric, or can be symmetrised as in quadratic forms.
Check that M is positive definite (Sylvester's criterion, or its eigenvalues). It will be, if the coordinates are independent.
Solvedet(K−ω2M)=0 for the n roots λk=ωk2. For n≤3 expand the determinant by hand; otherwise use a symmetric generalized eigensolver (scipy.linalg.eigh(K, M)), not eig(inv(M) @ K).
For each root solve (K−λkM)x=0 for the mode shape and normalise it by (xk)TMxk=1. For a degenerate λk the solution space has dimension >1; pick any basis of it and make it M-orthonormal by Gram-Schmidt in the M-inner product, which is the ordinary process with (x,y)M=yTMx in place of the dot product.
Check∑kλk=tr(M−1K) and (xj)TMxk=δjk before interpreting the answer.
Read off the physics.λk=0 means a free motion with no restoring force (a rigid translation or rotation, or a soft direction of V); λk<0 means an unstable equilibrium; the signs within xk say which parts move in phase.
The recipe requires M and K symmetric with M positive definite. If M is only semidefinite, or if the forces are not derived from a potential (damping, magnetic forces, gyroscopic terms), the pencil is no longer definite, its eigenvalues need not be real, and this method does not apply.
NoteOriginal
A coupled oscillating system looks complicated only in the wrong coordinates. The pencil (K,M) finds the coordinates in which it falls apart into n independent harmonic oscillators: the eigenvalues are the squared frequencies you can measure, and each eigenvector is the pattern of motion that one frequency belongs to.
The factor of one half. Riley, Hobson and Bence write T=q˙TAq˙ and V=qTBq with no 21, so their A and B are half of our M and K. Since det(B−ω2A)=2−ndet(K−ω2M), the normal frequencies and mode shapes are unaffected; only the normalisation constant of the modes changes. Use either convention, but not both in one calculation.
Positive definite K is about strict minima. A minimum of V gives K positive semidefinite; definiteness is stronger and fails exactly when V is flat to second order in some direction. Such directions are common: the linear triatomic molecule of Riley, Hobson and Bence §9.1 (p. 321) has ω2=0 for the mode in which all three masses translate together, since V depends only on the bond extensions. Part 3 above therefore separates "bounded motion" (K>0) from "no exponential growth" (K≥0). Note that Riley, Hobson and Bence restrict Chapter 9 from the start to a potential with a local minimum at the equilibrium and state there (p. 317) that this makes the potential matrix positive definite, like the kinetic one; the cases λk≤0 that parts 2 and 3 above admit therefore lie outside their setting, and are covered here by the pencil theorem, which assumes definiteness of M only.
Do not throw away the symmetry. Multiplying by M−1 gives the ordinary eigenproblem M−1Kx=ω2x, with the same eigenvalues and eigenvectors; but M−1K is in general neither symmetric nor normal, as Arfken, Weber and Harris point out for their three-mass chain (§6.5, pp. 322-324). Its eigenvectors are then not orthogonal in the plain dot product, and the left and right eigenvectors differ. M-orthogonality is the correct replacement, and solvers that take M as a second argument preserve it.
B must be definite, not merely invertible. For A=(0110) and B=(100−1) the equation det(A−λB)=−λ2−1=0 gives λ=±i: an indefinite B destroys the reality of the spectrum, which is why part 1 assumes B>0 and not just detB=0.
Example(A double pendulum with unequal masses)Original
Two point masses hang in a vertical plane on light rigid rods of equal length l: the upper mass is 2m, at the end of the first rod, and the lower mass is m, at the end of the second rod hinged to the first. Let q1 and q2 be the angles the two rods make with the downward vertical, so that q=0 is the hanging equilibrium. Write μ=g/l throughout.
Step 1: the energies to second order. The upper mass has speed lq˙1; the lower mass has velocity lq˙1 along the first rod's perpendicular plus lq˙2 along the second's, with angle q1−q2 between them, so its squared speed is l2(q˙12+q˙22+2q˙1q˙2cos(q1−q2)). Hence
using cosq=1−q2/2+O(q4). The constant is V(0)=−4mgl, the energy of the hanging configuration measured from the pivot; dropping it changes no force and no mode, and puts V in the normalisation V(0)=0 that the definition of small oscillations above asks for. Hence
V=21qTKq,K=mgl(3001).
Step 2: M is positive definite. Its leading principal minors are 3ml2>0 and (3−1)(ml2)2=2m2l4>0, so Sylvester's criterion (Quadratic forms and the Rayleigh quotient) applies and the pencil is definite. K is diagonal with positive entries, so it is positive definite too, and the equilibrium will turn out stable.
Step 3: the characteristic equation. Factor out ml2 from both matrices, which does not change the roots:
Expanding, 3(μ−λ)2−λ2=2λ2−6μλ+3μ2, a quadratic of degree exactly n=2 with leading coefficient 2=(−1)2det(3111), as part 1 of the pencil theorem predicts. Its roots are
Also λ1λ2=(9−3)μ2/4=23μ2, which matches detK/detM=3μ2/2. (Note that M−1K is not symmetric, exactly as the third remark above warns.)
Step 5: the mode shapes. The first row of (K−λM)x=0 reads 3(μ−λ)x1=λx2, so x2/x1=3(μ−λ)/λ. For λ1, μ−λ1=23−1μ and
x1x2=(3−3)/23(3−1)/2=3(3−1)3(3−1)=3,
while for λ2, μ−λ2=−21+3μ and the same computation gives x2/x1=−3. So, before normalisation,
x1∝(1,3)T,x2∝(1,−3)T.
In the slow mode the two rods swing the same way, the lower one through 3≈1.732 times the angle; in the fast mode they swing in opposition. Both are visible in the signs, which is the reading asked for in step 7 of the method.
Step 6: M-orthonormalisation and the orthogonality checks. With M~=(3111) so that M=ml2M~,
They are not orthogonal in the plain dot product: (1,3)⋅(1,−3)=1−3=−2=0. And K-orthogonality holds: (1,3)diag(3,1)(1,−3)T=3−3=0, while (1,±3)diag(3,1)(1,±3)T=6 for both, and indeed 6μ/(6+23)=λ1 and 6μ/(6−23)=λ2, as (xk)TKxk=λk demands.
Step 7: a trial-vector bound. Guess the slow mode as x=(1,1)T — both rods swinging together through equal angles, which is in the right ballpark but has the wrong amplitude ratio. Then
an upper bound, high by 5.2% in λ1 and only 2.5% in ω1 because of the square root. That one line of arithmetic gets the fundamental frequency to better than three percent with no eigenvalue calculation at all.
Step 8: the motion. Since both λk>0, K is positive definite and part 3 says every motion is bounded; explicitly
using (3+3)/(3−3)=(3+3)2/6=2+3 and ((2+6)/2)2=(8+43)/4=2+3. This is irrational, so a general motion never exactly repeats: the linearized double pendulum is quasi-periodic, not periodic, unless one of the two amplitudes vanishes.
ApplicationOriginal
Vibrational spectroscopy. A molecule of N atoms is a conservative system with 3N coordinates, and the apparatus of this unit is how its infrared and Raman spectrum is computed. The kinetic-energy matrix M in Cartesian displacements is diagonal, each atom contributing its mass three times; the stiffness matrix K is the Hessian of the electronic energy at the equilibrium geometry, which a quantum chemistry program computes by differentiating the energy twice. Diagonalizing the pencil (K,M) gives 3N eigenvalues, of which several vanish: the rigid motions change no bond and so leave V flat, precisely the λk=0 case of part 3. For a nonlinear molecule there are six such directions — three rigid translations and three rigid rotations — leaving 3N−6 vibrational frequencies ωk=λk to compare with the measured absorption lines. For a linear molecule only five vanish, because every atom sits on the axis and a rotation about that axis displaces no atom at all, so it is not a direction in configuration space; such a molecule has 3N−5 vibrations. Thus a linear triatomic, CO2 for instance, has 3⋅3−5=4 vibrational frequencies rather than the 3 that 3N−6 would give, and a diatomic has the one frequency of its stretching mode, where 3N−6 would predict none. The three-mass chain of the second remark above is the axial part of this calculation for such a molecule: only the three displacements along the axis are retained, so it yields three modes — the rigid translation with ω2=0 and two axial vibrations — while the remaining two of the four vibrations are bending motions transverse to the axis, which that reduced model does not carry. The eigenvectors are the pictures that accompany each line in a spectroscopy paper: which atoms move, how far, and in or out of phase. A negative eigenvalue in such a calculation is the standard diagnostic that the geometry found is not a minimum but a saddle point of the energy surface — a transition state — and the corresponding eigenvector is the reaction coordinate along which the molecule falls apart. Part 3 is exactly that statement: a negative λk means exponential departure along xk.
Modal analysis of structures, and the practical use of the quotient bound. A building, a bridge or a turbine blade is modelled by a finite element mesh whose nodal displacements q carry a mass matrix M and a stiffness matrix K, both large, sparse and symmetric with M positive definite. Engineers compute the lowest few eigenvalues of (K,M) because resonance with those frequencies is what destroys structures: a design is checked by confirming that the expected excitation — traffic, wind, rotor imbalance, the dominant frequencies of regional earthquakes — is far from every ωk. The method's warning in step 4 matters here: a solver given M as a second argument preserves the symmetry of the problem and the M-orthogonality of the computed modes, whereas forming M−1K destroys both, and on a large mesh it also destroys the sparsity that made the computation feasible. Rayleigh's method is the classical hand version of step 7 of the example: assume a plausible deflection shape and evaluate xTKx/xTMx. What the quotient bounds is λ1=ω12, from above; the square root of the quotient is then an upper bound on the fundamental frequency itself. The stationarity statement is what makes the estimate sharp: since the quotient is stationary at the true mode, its error is quadratic in the error of the shape, so a shape wrong by ten percent gives λ1 wrong by about one percent and ω1 wrong by only about half of one percent, the square root halving the relative error. Step 7 of the example is that arithmetic in miniature — a crude shape overestimates λ1 by 5.2% and ω1 by 2.5%.
Instability as a negative eigenvalue. The same pencil analysed about an unstable equilibrium delivers the growth rate rather than a frequency. An inverted pendulum, a column under a compressive load past its buckling threshold, or a plasma column in an adverse field all have a stiffness matrix K with at least one negative eigenvalue, and part 3 then gives the fastest-growing disturbance as the eigenvector xm for the most negative λm, with growth rate γm=−λm. Linear stability analysis in this form — assemble K, test its definiteness, and if it fails read off the shape of the unstable mode — is the standard first calculation for any equilibrium whose stability is in question, and the buckling load of a column is found by increasing the load until the smallest eigenvalue of K first crosses zero.
Every complex m-by-n matrix factors as A = USV* with U, V unitary and S non-negative diagonal; the singular values are the square roots of the min(m,n) largest eigenvalues of A*A, truncating the factorisation gives the best approximation of each rank, and inverting the non-zero part gives the pseudoinverse that solves least squares with minimum norm.
Definition(Singular values and the singular value decomposition)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.18.3, Eqs (8.131)-(8.133), pp. 302-303
Let A be a complex m×n matrix and p=min(m,n). A singular value decomposition (SVD) of A is a factorisation
A=USV†
in which U is m×m unitary, V is n×n unitary, and S is m×n with Sij=0 for i=j and diagonal entries
s1≥s2≥⋯≥sp≥0,
the singular values of A. The columns u1,…,um of U are the left singular vectors and the columns v1,…,vn of V the right singular vectors. For real A one may take U and V real orthogonal.
Every complex m×n matrix A has a singular value decomposition, and:
Where the singular values come from.A†A (n×n) and AA† (m×m) are Hermitian and positive semidefinite, so their eigenvalues are real and ≥0. The numbers s12,…,sp2 are the p largest eigenvalues of each, in decreasing order; the remaining eigenvalues of the larger of the two matrices are zero. Equivalently si=λi with λi the i-th largest eigenvalue of A†A. One may take for V any unitary whose columns are an orthonormal eigenbasis of A†A ordered by decreasing eigenvalue; every column of U belonging to a non-zero singular value is then forced, ui=Avi/si for each i with si>0, and the remaining ones are any orthonormal completion. The columns of U so obtained are an orthonormal eigenbasis of AA†. The two families cannot be chosen independently of each other.
Coupling of the two families. For i=1,…,p,
Avi=siui,A†ui=sivi,
while Avi=0 for i>p and A†ui=0 for i>p.
3. Rank, range, null space, dyadic form. Let r be the number of non-zero singular values. Then r=rankA; u1,…,ur are an orthonormal basis of the range of A; vr+1,…,vn are an orthonormal basis of its null space; v1,…,vr are an orthonormal basis of the range of A†; and
A=i=1∑rsiuivi†.
Uniqueness. The singular values are uniquely determined by A. The singular vectors are not, and the non-zero and the zero parts are free in different ways. Let r=rankA.
Non-zero singular values. Let W⊆{1,…,r} be a maximal block of indices on which si takes one common value s>0. Writing UW and VW for the matrices of the corresponding columns, the pair may be replaced by UWQ and VWQ for any ∣W∣×∣W∣ unitary Q — the sameQ on both sides, as the coupling Avi=sui of part 2 forces. The common phase ui↦eiαui, vi↦eiαvi is the case ∣W∣=1.
Zero singular values. The trailing columns ur+1,…,um and vr+1,…,vn are constrained only by being orthonormal completions, of a basis of kerA† and of kerA respectively. They contribute nothing to A=∑i≤rsiuivi†, so here no matching is required: the two completions may be changed by unrelated unitaries, of sizes m−r and n−r, which need not even be equal.
Every SVD of A arises from any one of them by such changes.
Relation to the eigenvalues. For a general square A the singular values are the eigenvalues of ∣A∣=(A†A)1/2, and bear no simple relation to the eigenvalues of A itself. If A is normal — that is, AA†=A†A, which holds in particular for Hermitian and for unitary A — with eigenvalues λ1,…,λn listed with multiplicity, then the singular values are the numbers ∣λk∣ sorted decreasingly.
Idea of the proofOriginal
There is only one Hermitian matrix in sight, A†A, and the spectral theorem hands us an orthonormal eigenbasis for it. Its eigenvalues measure ∥Ax∥22 along that basis, so their square roots are the stretching factors we are after, and dividing each image Avi by its own length produces an orthonormal family on the target side. Between those two bases A is diagonal, and that is the factorisation. The rest of the statement — the second eigenvalue list, the coupling, rank, range, null space, the dyadic form, and exactly how much freedom is left — is read off from that one diagonalisation.
A†A is Hermitian and positive semidefinite.(A†A)†=A†A, since (XY)†=Y†X† and (A†)†=A. For every x∈Cn,
x†A†Ax=(Ax)†(Ax)=∥Ax∥22≥0.
Each eigenvalue λ of A†A is real by the spectral theorem, and taking x to be an eigenvector for it the display reads λ∥x∥22≥0, so λ≥0. The same two lines with A† in place of A settle AA†.
An orthonormal eigenbasis on the right. By the spectral theorem pick an orthonormal basis v1,…,vn of Cn with A†Avi=λivi, numbered so that λ1≥λ2≥⋯≥λn≥0. Put si=λi≥0, and let r be the number of indices with λi>0; by the ordering these are i=1,…,r.
At most p of the si are non-zero. If A†Ax=0 then ∥Ax∥22=x†A†Ax=0 by step 1, so ker(A†A)=kerA and hence rank(A†A)=rankA (least squares, part 4). Since A†A is diagonalisable with exactly r non-zero eigenvalues, r=rankA≤min(m,n)=p. In particular λi=0 for every i>p, so the numbers s1≥⋯≥sp carry all of the information and fit on the diagonal of an m×n matrix.
Normalising the images. For i≤r put ui=Avi/si, which is legitimate because si>0. These vectors are orthonormal:
by the definition of ui, then the definition of the adjoint, then step 2, then the orthonormality of the vj. For i>r we have λi=0, so ∥Avi∥22=vi†A†Avi=0 and therefore Avi=0.
Completing the basis gives the factorisation. Since r≤m by step 3, the orthonormal set u1,…,ur extends to an orthonormal basis u1,…,um of Cm: take any basis of (span{u1,…,ur})⊥ and orthonormalise it by Gram-Schmidt. Let U and V be the matrices with those columns; both are unitary, their columns being orthonormal bases. Column j of U†AV is U†Avj, which by step 4 is sjU†uj=sjej when j≤r and 0 when j>r. As sj=0 for r<j≤p, that is precisely the matrix S of the definition, so U†AV=S, that is, A=USV†. The construction also establishes the recipe asserted in part 1: V may be any unitary whose columns are an orthonormal eigenbasis of A†A in decreasing order of eigenvalue, and then ui=Avi/si is forced whenever si>0, the remaining columns being any orthonormal completion.
Part 1, the two eigenvalue lists. Conversely, from any SVD A=USV†, unitarity of U and V gives
A†A=V(S†S)V†,AA†=U(SS†)U†.
Since Sij=0 off the diagonal, (S†S)jj=∑i∣Sij∣2 is sj2 for j≤p and 0 for j>p, and likewise (SS†)ii is si2 for i≤p and 0 for i>p; both matrices are diagonal. A unitary similarity does not change eigenvalues, so the eigenvalues of A†A are s12,…,sp2 together with n−p zeros, those of AA† the same numbers with m−p zeros, and the columns of V and of U are orthonormal eigenbases of the two matrices.
Part 2, the coupling. Multiplying A=USV† on the right by V gives AV=US. Column i of the left-hand side is Avi; column i of the right-hand side is siui for i≤p and 0 for i>p, because column i of S then vanishes. That is the first identity. Taking adjoints, A†=VS†U†, so A†U=VS†, and reading columns the same way — S is real, so S†=ST — gives A†ui=sivi for i≤p and A†ui=0 for i>p.
Part 3, rank, range, null space and the dyadic form. By part 2 together with si=0 for r<i≤p, Avi=0 for every i>r. Expanding an arbitrary x as x=∑i=1n(vi†x)vi,
valid for every x, which is the dyadic form. It gives imA⊆span{u1,…,ur}, and the reverse inclusion holds because Avi=siui with si=0; the ui are orthonormal, hence independent, so rankA=r and u1,…,ur is an orthonormal basis of the range. For the null space, Ax=0 says ∑i≤rsi(vi†x)ui=0, which by independence of the ui and si=0 forces vi†x=0 for every i≤r, that is x∈span{vr+1,…,vn}; and conversely those vi are in kerA. Applying what has just been proved to A†=∑i≤rsiviui† shows that v1,…,vr is an orthonormal basis of imA†.
Part 4, uniqueness of the singular values. By step 6 the numbers s12≥⋯≥sp2 are the p largest eigenvalues of A†A, hence determined by A alone, and the constraint si≥0 fixes the square roots.
Part 4, the freedom in the singular vectors. Let A=USV†=U′SV′† be two SVDs, with the same S by step 9. By step 6 the columns of V and of V′ are both orthonormal eigenbases of A†A in decreasing order of eigenvalue. Two orthonormal bases adapted to the same eigenspaces in the same order differ exactly by a unitary change of basis inside each eigenspace. Those eigenspaces are span{vi:i∈W} for each maximal block W of indices carrying one common value s>0, and span{vr+1,…,vn}=kerA for the eigenvalue 0. So, columnwise, VW′=VWQW with QW unitary of size ∣W∣, while vr+1′,…,vn′ is an arbitrary orthonormal basis of kerA. For i≤r the left vectors are then forced, ui′=Avi′/si by step 5, so
UW′=AVW′/s=AVWQW/s=UWQW
with the same unitary QW. The trailing columns ur+1′,…,um′ complete u1′,…,ur′, a basis of imA, to an orthonormal basis of Cm, so they form an arbitrary orthonormal basis of (imA)⊥=kerA†, of dimension m−r and unrelated to the n−r dimensions available on the other side. Conversely each such change does produce an SVD, because inside a block UWQW(VWQW)†=UWQWQW†VW†=UWVW† leaves the dyadic sum of part 3 unchanged, and the trailing columns never enter it.
Part 5, the normal case. For square A (so m=n=p), step 6 gives A†A=Vdiag(s12,…,sn2)V†, and B=Vdiag(s1,…,sn)V† is positive semidefinite with B2=A†A. It is the only such matrix, so the notation ∣A∣=(A†A)1/2 is unambiguous: a positive semidefinite C with C2=A†A commutes with A†A, hence maps each eigenspace of A†A into itself, and on the eigenspace for μ it is a positive semidefinite matrix whose square is μI, so it is Hermitian with every eigenvalue equal to μ and therefore equals μI there. The eigenvalues of ∣A∣=B are thus s1,…,sn. If in addition A is normal, normal matrices supplies a unitary W with A=WΛW† and Λ=diag(λ1,…,λn); then A†A=WΛΛW†=Wdiag(∣λ1∣2,…,∣λn∣2)W†, so the eigenvalues of A†A are the ∣λk∣2 and the singular values, their square roots in decreasing order, are the ∣λk∣ sorted decreasingly. □
Theorem(Spectral norm and best low-rank approximation)Original
Let A be m×n with singular values s1≥⋯≥sp≥0 and rank r, let ∥A∥2=maxx=0∥Ax∥2/∥x∥2 and ∥A∥F=(∑i,j∣Aij∣2)1/2. Then ∥A∥2=s1 and ∥A∥F=(∑i=1psi2)1/2, and for each integer k with 1≤k<r the truncation
Ak=i=1∑ksiuivi†
(of rank exactly k) is a best rank-k approximation to A in both norms: over all m×n matrices B with rankB≤k,
The Frobenius case is the Eckart-Young theorem (1936), whose integral-operator analogue goes back to Schmidt (1907); Mirsky (1960) extended it to every unitarily invariant norm, which is where the spectral-norm case comes from. The Frobenius minimiser is unique exactly when sk>sk+1; the spectral-norm minimiser is in general not unique. For k≥r the minimum is 0, attained by A itself.
Idea of the proofOriginal
Both norms are blind to the unitary factors, so each is a function of the singular values alone, and the error made by chopping the dyadic sum is immediate. The two lower bounds are the real content. In the spectral norm, a matrix of rank at most k must annihilate some direction inside the (k+1)-dimensional span of v1,…,vk+1, and A stretches every direction there by at least sk+1. In the Frobenius norm, split A−B along the range of B: the part of A pointing out of that range survives untouched, and a subspace of dimension k can capture at most the k largest eigenvalues' worth of AA†.
ProofOriginal
Throughout, fix an SVD A=USV† as in the existence statement, so that A=∑i≤rsiuivi†, and write (X,Y)F=tr(Y†X), which is an inner product on m×n matrices with ∥X∥F2=(X,X)F and is linear in its first slot.
The dyads are orthonormal in the Frobenius inner product. For i,j≤p, using tr(xy†)=y†x for column vectors,
by the orthonormality of the u's and then of the v's. Hence ∑iciuivi†F2=∑i∣ci∣2 for any scalars ci.
The Frobenius norm and the Frobenius truncation error. Applying step 1 to A=∑i≤rsiuivi† gives ∥A∥F2=∑i≤rsi2=∑i≤psi2, the omitted terms being zero. Since A−Ak=∑k<i≤rsiuivi†, the same step gives ∥A−Ak∥F2=∑i>ksi2.
The spectral norm and the spectral truncation error. For any x put y=V†x, so ∥y∥2=∥x∥2 as V is unitary. Then ∥Ax∥2=∥USy∥2=∥Sy∥2, again by unitarity, and
∥Sy∥22=i≤p∑si2∣yi∣2≤s12∥y∥22,
with equality at x=v1 (then y=e1); so ∥A∥2=s1. Running the same computation for A−Ak, whose dyadic form has only the indices k<i≤r, gives ∥(A−Ak)x∥22=∑k<i≤rsi2∣vi†x∣2≤sk+12∥x∥22 with equality at x=vk+1, so ∥A−Ak∥2=sk+1. By the range argument of part 3 of the existence statement, applied to Ak=∑i≤ksiuivi† with s1,…,sk>0 (here k<r), imAk=span{u1,…,uk}, so rankAk=k exactly.
Spectral lower bound. Let rankB≤k. Then dimkerB≥n−k, while dimspan{v1,…,vk+1}=k+1, and (n−k)+(k+1)>n, so the two subspaces of Cn intersect in a non-zero vector. Normalise it: w=∑i≤k+1civi with ∑i≤k+1∣ci∣2=1 and Bw=0. Then, by part 2 of the existence statement and the orthonormality of the ui,
the last step because si≥sk+1 for every i≤k+1 and the ∣ci∣2 sum to 1. Since ∥w∥2=1, this gives ∥A−B∥2≥sk+1, which with step 3 proves the spectral-norm minimum.
A capture lemma.Let M be an m×m Hermitian positive semidefinite matrix with eigenvalues μ1≥⋯≥μm≥0 and an orthonormal eigenbasis w1,…,wm, and let q1,…,qj be orthonormal in Cm with j≤k. Then ∑l≤jql†Mql≤∑i≤kμi. To see it, expand M=∑iμiwiwi† and set ti=∑l≤j∣wi†ql∣2, so that ∑l≤jql†Mql=∑iμiti. Now ti is the squared norm of the orthogonal projection of wi onto span{q1,…,qj}, hence 0≤ti≤∥wi∥22=1 (orthogonal projection, parts 3 and 4), while summing over i instead gives ∑iti=∑l≤j∥ql∥22=j≤k, because w1,…,wm is an orthonormal basis. Therefore
where the first inequality uses ti−1≤0 with μi≥μk for i≤k and ti≥0 with μi≤μk for i>k, and the second uses μk≥0 with ∑iti≤k.
Frobenius lower bound. Let rankB≤k, let L=imB, of dimension j≤k, and let Q be the orthogonal projection of Cm onto L, so QB=B. Split
A−B=(I−Q)A+(QA−B).
Every column of the first term lies in L⊥ and every column of the second in L, so ((I−Q)A,QA−B)F=0, the trace being the sum of the columnwise inner products. Hence
by step 5 with M=AA†, whose eigenvalues are s12,…,sp2 and zeros by part 1 of the existence statement. Combining, ∥A−B∥F2≥∥A∥F2−∑i≤ksi2=∑i>ksi2, which with step 2 proves the Frobenius minimum.
The Frobenius minimiser is unique when sk>sk+1. Suppose rankB≤k and ∥A−B∥F2=∑i>ksi2. Equality in step 6 forces B=QA and ∑l≤jql†AA†ql=∑i≤ksi2, so both inequalities of step 5 are tight for M=AA†, μi=si2. Here μk=sk2>0 because k<r, so tightness of the second gives ∑iti=j=k. Tightness of the first says
i≤k∑(μi−μk)(1−ti)+i>k∑(μk−μi)ti=0,
a sum of non-negative terms. Choosing wi=ui as the eigenbasis of AA†, the hypothesis sk>sk+1 makes μk−μi>0 for every i>k, so ti=0 there; then ∑i≤kti=k with each ti≤1 forces ti=1 for i≤k, that is ui∈L for i≤k. As dimL=k, this gives L=span{u1,…,uk} and Q=∑i≤kuiui†. Finally ui†A=(A†ui)†=sivi† by part 2 of the existence statement, so
B=QA=i≤k∑ui(ui†A)=i≤k∑siuivi†=Ak.
This also shows Ak does not depend on which SVD was fixed at the start, as it must not: it is the unique minimiser.
It is not unique when sk=sk+1. Write s for the common value, which is positive because k<r. For real θ let
uk′=cosθuk+sinθuk+1,uk+1′=−sinθuk+cosθuk+1,
and define vk′,vk+1′ by the same formulas. By part 4 of the existence statement — the case of a block W⊇{k,k+1} of equal singular values, rotated by the same real orthogonal Q on both sides — this is again an SVD of A. Its rank-k truncation Ak′=∑i<ksiuivi†+suk′vk′† is a best rank-k approximation, by steps 2 and 6 applied to that SVD. But vk′†vk+1=sinθ, so Ak′vk+1=ssinθuk′, while Akvk+1=0. For sinθ=0 the two differ, so there are infinitely many minimisers.
The case k≥r. Then rankA=r≤k, so B=A is admissible and both minima are 0. □
Definition(Pseudoinverse)AdaptedAdapted from Stakgold & Holst (2011), Exercises 5.5.3-5.5.4, pp. 325-326
With A=USV† as above and r non-zero singular values, the Moore-Penrose pseudoinverse of A is the n×m matrix
so that A+=∑i=1rsi−1viui†. It does not depend on which SVD of A is used. For any b∈Cm, the vector x+=A+b is the minimum-norm least-squares solution: it minimises the residual ∥Ax−b∥2 over all x∈Cn, and among the minimisers it is the unique one of smallest ∥x∥2.
That last property is what (Stakgold & Holst (2011), Exercises 5.5.3-5.5.4, pp. 325-326) takes as the definition of the pseudoinverse, for an operator with closed range and with no factorisation in sight. The formula A+=VS+U†, and the identification of the two, are ours; the cited pages do not contain them. Consequently
AA+ is the orthogonal projection of Cm onto the range of A, and A+A the orthogonal projection of Cn onto the range of A†, the orthogonal complement of kerA;
if A is square and invertible then A+=A−1;
if A has full column rank n then A+=(A†A)−1A†, recovering the normal equations of least squares.
Idea of the proofOriginal
In the singular bases A is a diagonal matrix, and for a diagonal matrix the two demands are transparent: the residual splits into coordinates the unknown can set exactly and coordinates it cannot touch at all, and the norm of the unknown splits into coordinates the data forces and coordinates that are free. Setting the free ones to zero is the whole content of the pseudoinverse. Since this description of A+b never mentions a basis, it cannot depend on which SVD was used.
Proof(Properties of the pseudoinverse)Original
Fix an SVD A=USV† with r non-zero singular values, let b∈Cm, and for x∈Cn put y=V†x and c=U†b. Unitarity gives ∥x∥2=∥y∥2, and the correspondence x↔y is a bijection of Cn.
The residual in the singular bases. Since U is unitary, ∥Ax−b∥2=∥U(Sy−c)∥2=∥Sy−c∥2. The matrix S has (Sy)i=siyi for i≤p and (Sy)i=0 for p<i≤m, and si=0 for r<i≤p, so
∥Ax−b∥22=i=1∑r∣siyi−ci∣2+i=r+1∑m∣ci∣2.
The second sum does not involve x at all. Hence x minimises the residual if and only if yi=ci/si for every i≤r, the coordinates yr+1,…,yn being unconstrained, and the minimal residual is (∑i>r∣ci∣2)1/2.
Minimum norm. Among those minimisers, ∥x∥22=∥y∥22=∑i≤r∣ci/si∣2+∑i>r∣yi∣2, which is smallest exactly when yi=0 for every i>r, and then uniquely so. That choice of y is precisely S+c, so the minimum-norm least-squares solution is x=VS+U†b=A+b, as claimed. Taking b to run over Cm identifies A+ with the operator of the cited exercises.
Independence of the SVD. Steps 1 and 2 characterise A+b as the unique minimum-norm minimiser of ∥Ax−b∥2, a description that mentions only A and b. So A+ is the same matrix for every SVD of A.
The two projections.SS+ is the m×m diagonal matrix with r ones followed by zeros, since (SS+)ij=Sii(S+)ij=sisi−1δij for i≤r and vanishes otherwise. Hence
AA+=USV†VS+U†=U(SS+)U†=i=1∑ruiui†,
which is Hermitian and idempotent with range span{u1,…,ur}=imA, so it is the orthogonal projection onto imA (orthogonal projection, part 5). Identically, S+S is n×n diagonal with r ones, so A+A=∑i≤rvivi† is the orthogonal projection onto span{v1,…,vr}=imA†=(kerA)⊥.
Square and invertible. Then r=m=n, so S is square with all si>0 and S+=S−1; hence A+=VS−1U†=(USV†)−1=A−1.
Full column rank. If rankA=n then r=n≤m, and S†S=diag(s12,…,sn2) is invertible, so A†A=V(S†S)V† is invertible with inverse V(S†S)−1V†. Therefore
(A†A)−1A†=V(S†S)−1V†⋅VS†U†=V((S†S)−1S†)U†,
and (S†S)−1S† is the n×m matrix with (i,i) entry si/si2=1/si and all other entries zero, which is S+. So (A†A)−1A†=VS+U†=A+, and by step 1 the minimiser is unique because no coordinates yi with i>r=n exist. □
NoteOriginal
The SVD says that every linear map, of any shape and any rank, is the same thing in disguise: choose the right orthonormal basis in the source and in the target, and the map simply scales r perpendicular directions by s1,…,sr and annihilates the rest. Because the scaling factors come out sorted, truncating the list is the optimal way to simplify the map — which is why the SVD underlies principal component analysis, image and model compression, rank determination, and the regularisation of ill-posed least-squares problems.
Remark(Traps and conventions)Original
Singular values are not eigenvalues. For A=(0020) both eigenvalues are 0 while the singular values are 2 and 0. The two lists agree up to absolute values and ordering exactly for normalA, meaning AA†=A†A (see normal matrices) — not only for Hermitian A: a 2×2 rotation matrix and diag(i,1) are normal without being Hermitian, and both have singular values 1,1 and eigenvalues of modulus 1. (Normality is necessary as well as sufficient, since ∑isi2=∥A∥F2≥∑k∣λk∣2 for every square A, with equality precisely when A is normal.) The two lists coincide, with no absolute values, exactly for positive semidefinite A.
si≥0 and the ordering are conventions, and they are what makes the factorisation nearly unique. A sign can always be moved into ui or vi; reversing the order permutes the columns of U and V together. Do not compare two SVDs column by column without fixing both conventions.
Full versus reduced. Many texts and all numerical libraries offer a "thin" or "economy" form in which U is m×p and V is n×p with orthonormal columns and S is p×p square. The thin factors are isometries, so U†U=Ip and V†V=Ip hold in both forms, while the identities with the factors the other way round can fail. Which one fails depends on the shape: since p=min(m,n), the thin U is square — hence already unitary, with UU†=Im — exactly when m≤n, and the thin V is square exactly when n≤m. For m>n the thin U is m×n and UU† is the rank-n orthogonal projection onto span{u1,…,un}, not Im; for n>m the same happens to VV†. At most one of the two can fail, and when m=n exactly one does; when m=n the thin form is the full form, both factors are unitary and neither fails. So "the thin factors are not unitary" is a statement about one named factor, decided by the shape, never about both at once.
The pseudoinverse is not continuous in A. A singular value that passes through zero makes A+ blow up: for Aϵ=diag(1,ϵ), Aϵ+=diag(1,1/ϵ) for ϵ=0 but A0+=diag(1,0). This is why practical computations zero out all singular values below a tolerance rather than inverting them, and it is the quantitative content of conditioning.
A+=(A†A)−1A† needs full column rank. With dependent columns A†A is singular and the normal equations have many solutions; the pseudoinverse selects the one of least norm, which the normal equations by themselves do not.
Notation.A+, A† and A# all appear in the literature for the pseudoinverse; we reserve A† for the Hermitian conjugate throughout, as in the spectral theorem. Singular values are written si here and σi in much of the numerical literature.
Example(A rank-deficient 3 x 3 matrix, worked through completely)Original
Take the real matrix
A=6069−666−120,
so m=n=p=3 and A†=AT.
1. Form ATA. Its (i,j) entry is the dot product of columns i and j of A:
ATA=7290369015312636126180.
2. Its eigenvalues. The characteristic polynomial of a 3×3 matrix is λ3−e1λ2+e2λ−e3 with e1 the trace, e2 the sum of the three 2×2 principal minors and e3 the determinant. Here e1=72+153+180=405,
Each has norm 311+4+4=1, and the three are mutually orthogonal (1⋅2+2⋅1−2⋅2=0, and similarly for the other two pairs). Verify one eigenvector equation: ATA(1,2,2)T=(72+180+72,90+306+252,36+252+360)T=(324,648,648)T=324(1,2,2)T, as required.
4. Left singular vectors. By the construction, ui=Avi/si for i≤r=2:
The third column is any unit vector orthogonal to these two; u3=31(2,1,−2)T serves. Checks demanded by the statement: Av3=0, so v3 spans kerA; ATu3=0, so u3 spans kerAT=(imA)⊥; ATu1=(6,12,12)T=18v1 and ATu2=(6,3,−6)T=9v2, which is part 2.
5. The decomposition. With U=(u1u2u3), V=(v1v2v3) and S=diag(18,9,0),
and the two matrices do add up to A. As a check on the singular values, ∥A∥F2=36+81+36+36+144+36+36=405 (the three zero entries of A contribute nothing) and s12+s22=324+81=405.
6. Best rank-one approximation.A1=18u1v1T is the first matrix displayed above, and A−A1 is the second. Its spectral and Frobenius norms are both s2=9: indeed A−A1=9u2v2T is a single dyad, and ∥A−A1∥F2=4+1+4+16+4+16+16+4+16=81. No matrix of rank 1 comes closer in either norm, and because s1>s2 this A1 is the only Frobenius minimiser.
7. Pseudoinverse.A+=181v1u1T+91v2u2T, which works out to
A+=54122020−432−2.
Check the projection property: AA+=I−u3u3T=915−24−282425, which is symmetric, squares to itself, and kills u3 — the orthogonal projection onto imA.
The residual is 36u3, so it is orthogonal to imA as it must be, and its norm is 36. Every least-squares solution is x++tv3 with t∈R, and since x+⋅(2,−2,1)=4−4+0=0 the vector x+ is orthogonal to kerA, giving
∥x++tv3∥22=∥x+∥22+t2=8+t2,
minimal only at t=0. So ∥x+∥2=22≈2.8284 is the smallest norm among all least-squares solutions, and ∥b−Ax+∥2=36 is the common residual of all of them.
ApplicationOriginal
Empirical orthogonal functions in climate science. Arrange a gridded field — monthly sea-surface temperature anomalies, say — as a matrix with one row per time and one column per grid point, and compute its SVD. The right singular vectors are the spatial patterns (EOFs), the left singular vectors scaled by the singular values are their time series (principal components), and si2/∑jsj2 is the fraction of variance that pattern carries. The best low-rank approximation statement is what justifies keeping the first few: no other combination of that many patterns reproduces the data field better, in the least-squares sense. The same computation under the name principal component analysis reduces the dimension of any data matrix.
Quantum entanglement. A pure state of two subsystems, ∣ψ⟩=∑jkCjk∣j⟩⊗∣k⟩, has a coefficient matrix C; its SVD is the Schmidt decomposition ∣ψ⟩=∑isi∣ui⟩⊗∣vi⟩. The singular values are the Schmidt coefficients, rankC is the Schmidt rank, the state is a product state exactly when that rank is 1, and the reduced density matrices CC† and CTC have the common eigenvalues si2, from which the entanglement entropy −∑isi2logsi2 (for ∑isi2=1) is computed. Truncating to the largest few si is the compression step at the heart of matrix-product-state and DMRG methods.
Regularising an ill-posed inversion. In tomography, deconvolution and geophysical inversion one must solve Ax=b for a discretised operator whose singular values decay towards zero. The pseudoinverse multiplies the data component along ui by 1/si, so noise in the directions with tiny si is amplified without limit. Truncated SVD regularisation keeps only the terms with si above a threshold, which replaces A+ by the pseudoinverse of the best low-rank approximation Ak and trades resolution for stability; the quantitative version of the trade-off is the subject of condition numbers.
Rank, and how many digits of it you can believe. The theoretical rank of a measured matrix is almost never attained numerically, because rounding and measurement error perturb every zero singular value away from zero. The spectral-norm statement says that A is within sk+1 of a matrix of rank k, so counting the singular values above the noise level is the defensible way to determine rank — the standard test for the order of a linear model, the number of independent reaction pathways in a kinetics matrix, or the dimension of an attractor reconstructed by delay embedding.
Defines a function of a matrix by its power series, states when the series converges, collects the properties of exp(A) including the exact condition — that A and B commute — for exp(sA+tB) = exp(sA)exp(tB) to hold for all s and t near 0, which is more than the single identity exp(A+B) = exp(A)exp(B) requires, and gives the recipes for computing exp(tA) by diagonalization and for defective matrices.
Powers of a square matrix are defined by A0=I and Ak=AAk−1, so for scalars c0,c1,… one can form
f(A)=k=0∑∞ckAk,
a matrix of the same size as A, provided the series converges. Convergence means entrywise convergence, which is the same as convergence in any norm on the space of matrices. The matrix exponential is the case ck=1/k!:
expA=eA=k=0∑∞k!Ak.
The same construction defines sinA and cosA, which converge for every A, and the geometric (Neumann) series ∑kAk, which converges to (I−A)−1 when ρ(A)<1 and diverges when ρ(A)>1. Here the spectral radius of A is ρ(A)=max{∣λ∣:λ∈σ(A)}, and the convergence statement below is what supplies these conditions.
Method(Convergence)Original
Let ∑kckzk have radius of convergence R∈(0,∞] and let A be an n×n complex matrix.
If ρ(A)<R, then ∑kckAk converges absolutely, entry by entry. If ρ(A)>R, it diverges.
The boundary case ρ(A)=R may go either way, and is not decided by ρ(A) alone. At A=I, where ρ(A)=1, both the coefficient sequence ck=1 and the coefficient sequence c0=0, ck=1/k2 for k≥1, have R=1; the first matrix series diverges and the second converges.
The exponential has R=∞, so part 1 covers it for every A. The next statement covers it again with an explicit bound and no reference to the spectrum.
Idea of the proofOriginal
A matrix power series can only do to A what the scalar series does to each of its eigenvalues, so what matters is how fast the entries of Ak can grow. The Jordan form answers that: they grow like a polynomial in k times ρ(A)k, and a polynomial factor cannot defeat the geometric margin available strictly inside the disc of convergence. Outside the disc a single eigenvector is enough to wreck everything, because on that vector the matrix series is the scalar series evaluated at z=λ. On the circle the geometric margin is gone, and the polynomial factor that it used to absorb becomes decisive instead: at ρ(A)=R the answer depends on the Jordan structure of A and not on the spectrum or the coefficients alone, so no condition on ρ(A) can settle it. The degenerate case A=I, where that structure drops out entirely and the matrix series is just (∑kck)I, is therefore the place to look when one wants two coefficient sequences of the same radius to behave differently at the same A.
ProofOriginal
Entrywise growth of the powers.Claim: for every r>ρ(A) there is a constant C, depending on A and r only, with ∣(Ak)il∣≤Crk for all k≥0 and all i,l. Write A=TJT−1 with J the Jordan normal form of A, so Ak=TJkT−1; since T and T−1 are fixed, each entry of Ak is a fixed finite combination of entries of Jk, and it is enough to bound those. J is a direct sum of blocks λIq+Nq with q≤n, Nq the nilpotent shift (Nqq=0, entries 0 or 1) and λ an eigenvalue, so ∣λ∣≤ρ(A)<r. As λIq commutes with Nq, the binomial theorem gives
(λIq+Nq)k=j=0∑min(k,q−1)(jk)λk−jNqj,
so every entry of the block of Jk is (jk)λk−j for some j≤q−1, or 0. If λ=0 only j=k survives, which needs k≤q−1: finitely many k, absorbed into C. If λ=0 then ∣λk−j∣≤max(1,∣λ∣−(n−1))∣λ∣k and (jk)≤kn−1 for k≥1, so the entry is at most a constant times kn−1∣λ∣k≤kn−1ρ(A)k. Finally kn−1ρ(A)k/rk→0 as k→∞, because ρ(A)<r makes the geometric factor win, so that sequence is bounded and the claim follows.
Part 1, convergence. Suppose ρ(A)<R and choose r with ρ(A)<r<R. A power series converges absolutely at every point strictly inside its disc of convergence (power series), so ∑k∣ck∣rk<∞. By step 1,
k=0∑∞∣ck∣∣(Ak)il∣≤Ck=0∑∞∣ck∣rk<∞
for every i,l: the matrix series converges absolutely, entry by entry.
Part 1, divergence. Suppose ρ(A)>R and pick an eigenvalue λ with ∣λ∣=ρ(A) and an eigenvector v=0. Then Akv=λkv by induction, so the partial sums PK=∑k≤KckAk satisfy PKv=(∑k≤Kckλk)v. If the matrix series converged entrywise then each (PKv)i, a fixed linear combination of entries of PK, would converge; choosing i with vi=0 makes ∑k≤Kckλk=(PKv)i/vi a convergent scalar sequence. But ∣λ∣>R, and a power series diverges at every point outside its disc of convergence. So the matrix series diverges.
Part 2, the boundary. For A=I every power is I, so ∑kckAk=(∑kck)I converges exactly when the scalar series ∑kck does, while ρ(I)=1. Taking ck=1 gives R=1 (the coefficients have ∣ck∣1/k=1) and a divergent ∑k1. Taking c0=0 and ck=1/k2 gives R=1 as well, since (1/k2)1/k→1, and ∑k≥11/k2 converges. So ρ(A)=R decides nothing.
The Neumann series. For ∑kzk we have R=1, so by part 1 the series ∑kAk converges when ρ(A)<1 and diverges when ρ(A)>1. In the convergent case, telescoping gives (I−A)∑k≤KAk=I−AK+1, and by step 1 with any r∈(ρ(A),1) the entries of AK+1 are at most CrK+1→0. Letting K→∞ gives (I−A)∑kAk=I, and the same computation on the other side gives ∑kAk(I−A)=I, so the sum is (I−A)−1. □
Let A be an n×n complex matrix and put a=maxi,j∣Aij∣. Then
(Ak)ij≤nk−1ak(k≥1),
and consequently every entry of ∑kAk/k! is an absolutely convergent scalar series, with
k=0∑∞k!∣(Ak)ij∣≤1+nena−1.
So expA is defined for every A, under no hypothesis at all. The same bound makes ∑kckAk converge absolutely whenever na<R: a crude but spectrum-free test, which needs no eigenvalue of A to be known. A sharp test of that kind uses a submultiplicative norm — one with ∥XY∥≤∥X∥∥Y∥, so that ∥Ak∥≤∥A∥k — together with the fact that ρ(A)≤∥A∥; such norms are the subject of a later unit, matrix norms and conditioning, and nothing below depends on them.
Idea of the proofOriginal
Each entry of Ak+1 is a sum of n products of an entry of Ak with an entry of A, so passing from one power to the next costs at most one factor n and one factor a; iterating that is the bound. Dividing by k! then leaves, up to the factor 1/n, a term of the scalar series for ena, and the comparison test finishes the job. Nothing about the spectrum of A is used anywhere, which is what makes the resulting test crude but usable with the entries of A alone.
The bound on the powers. Induction on k. For k=1 the claim reads ∣Ail∣≤a, which is the definition of a. Assume it for k. Each entry of Ak+1 is ∑q=1n(Ak)iqAql, a sum of n terms, so by the triangle inequality, the inductive hypothesis and ∣Aql∣≤a,
∣(Ak+1)il∣≤q=1∑nnk−1ak⋅a=nkak+1.
Comparison with the scalar exponential series. Split off the term k=0, where (A0)il=δil has modulus at most 1. For k≥1 step 1 gives ∣(Ak)il∣/k!≤(na)k/(nk!), so
the middle equality because ∑k≥1(na)k/k!=ena−1. Each entry series is therefore dominated by a convergent series of non-negative numbers, hence absolutely convergent by the comparison test (convergence tests). This is the route of (Hirsch, Smale & Devaney (2013), §6.4, p. 125), which bounds the entry by ena without separating the k=0 term.
The spectrum-free test. Run the same chain with the coefficients in place: ∑k≥1∣ck∣∣(Ak)il∣≤n1∑k≥1∣ck∣(na)k by step 1. If na<R then na lies strictly inside the disc of convergence of ∑kckzk, where a power series converges absolutely (power series), so the right-hand side is finite and ∑kckAk converges absolutely entry by entry.
Relation to the previous statement. The conclusion for expA is also the case R=∞ of part 1 above. The bound just proved is nevertheless the one used in the rest of this unit, because it is explicit in the entries of A and survives being applied to sA+tB with unknown spectrum. □
Let A,B be n×n complex matrices and T invertible. Then:
exp(T−1AT)=T−1(expA)T.
If AB=BA, then exp(A+B)=(expA)(expB).
exp(−A)=(expA)−1; in particular expA is invertible for every A.
If Av=λv with v=0, then (expA)v=eλv.
dtdexp(tA)=Aexp(tA)=exp(tA)A, and x(t)=exp(tA)x0 is a solution of x˙=Ax with x(0)=x0 and is the only one ((Hirsch, Smale & Devaney (2013), §6.4, Theorem, p. 129)).
det(expA)=etrA.
The commuting hypothesis is what makes the proof of property 2 work: that proof goes through the binomial theorem for (A+B)k, which is valid only when AB=BA. It is not necessary for the identity itself — see part 3 of the statement "When the product rule holds" below, where a non-commuting pair satisfying exp(A+B)=exp(A)exp(B) is exhibited.
Idea of the proofOriginal
Each property is a scalar identity about ex transplanted to matrices, and the transplant is always the same two moves: do the algebra on the partial sums, where everything is a finite sum, then pass to the limit, which is legitimate because multiplying by a fixed matrix, or applying a fixed linear functional, is continuous entrywise. Two of them need more than that. The product rule has two infinite series to multiply, which is the Cauchy product, available because the exponential series converges absolutely; it is the binomial theorem afterwards that insists on AB=BA. And the derivative in property 5 needs an estimate rather than a limit: a difference quotient has to be controlled as h→0, and what controls it is the quantitative entrywise tail bound of the previous statement applied to (exp(hA)−I)/h−A. The determinant identity is of a different kind: put A into triangular form and both sides can be read off the diagonal.
Property 1.(T−1AT)k=T−1AkT, by induction: the inductive step cancels the adjacent TT−1 in (T−1AkT)(T−1AT). Hence the partial sums satisfy ∑k≤K(T−1AT)k/k!=T−1(∑k≤KAk/k!)T. Every entry of the right-hand side is one fixed linear combination of the entries of ∑k≤KAk/k!, with coefficients built from T and T−1, so it converges to the same combination of the entries of expA. Letting K→∞ gives exp(T−1AT)=T−1(expA)T.
Property 2. Fix entries i,l. By the statement "The exponential converges for every matrix" each scalar series ∑j(Aj)iq/j! and ∑k(Bk)ql/k! converges absolutely, so the Cauchy product of two absolutely convergent series (absolute convergence) applies to each of the n pairs indexed by q; summing the n resulting identities,
with no hypothesis on A and B so far. Now the hypothesis enters: since AB=BA, also BAj=AjB for every j, and the binomial theorem holds in the form (A+B)N=∑j+k=Nj!k!N!AjBk (induction on N, the step using BAj=AjB to move every B past every A). So the inner sum is (A+B)N/N! and the right-hand side is exp(A+B).
Property 3.A and −A commute, so property 2 gives (expA)(exp(−A))=exp(A−A)=exp0=I, and with the roles exchanged (exp(−A))(expA)=I. Hence expA is invertible with inverse exp(−A).
Property 4. From Av=λv we get Akv=λkv by induction, so (∑k≤KAk/k!)v=(∑k≤Kλk/k!)v. The left side converges to (expA)v because expA converges entrywise, the right to eλv because the scalar exponential series converges. Hence (expA)v=eλv.
Property 5. For fixed t and h the matrices tA and hA commute, so property 2 gives exp((t+h)A)=exp(tA)exp(hA) and therefore
hexp((t+h)A)−exp(tA)=exp(tA)⋅hexp(hA)−I.
Now hexp(hA)−I−A=∑k≥2hk−1Ak/k!, whose entries are bounded, for 0<∣h∣≤1, using the entrywise bound ∣(Ak)il∣≤nk−1ak with a=maxi,l∣Ail∣, by
which tends to 0 with h. So hexp(hA)−I→A entrywise, and the difference quotient converges to exp(tA)A: that is the derivative. Moreover A commutes with every Ak, hence with every partial sum of the series, hence with exp(tA) in the limit, so the derivative equals Aexp(tA) too. Consequently x(t)=exp(tA)x0 satisfies x˙=Aexp(tA)x0=Ax(t) and x(0)=x0. For uniqueness, let y be any solution with y(0)=x0 and set z(t)=exp(−tA)y(t). By the product rule and the derivative just computed, applied with −A,
z˙=−Aexp(−tA)y+exp(−tA)y˙=exp(−tA)(−Ay+Ay)=0,
where the commuting of A with exp(−tA) was used to move −A across. So z is constant, z(t)=z(0)=x0, and multiplying by exp(tA) — invertible by property 3, with inverse exp(−tA) — gives y(t)=exp(tA)x0. This is the existence-and-uniqueness theorem cited in the statement.
Property 6. Let A=TJT−1 with J the Jordan form, which is upper triangular with the eigenvalues λ1,…,λn of A on its diagonal (listed with algebraic multiplicity). Multiplying two upper triangular matrices X,Y leaves the result upper triangular and multiplies the diagonals entrywise, as (XY)ii=∑qXiqYqi=XiiYii shows — only q=i survives, since Xiq=0 for q<i and Yqi=0 for q>i. So Jk is upper triangular with (Jk)ii=λik. Each partial sum of expJ is therefore upper triangular with (i,i) entry ∑k≤Kλik/k!, and these properties survive the entrywise limit: expJ is upper triangular with diagonal entries eλi. The determinant of a triangular matrix is the product of its diagonal entries, so
detexpJ=i=1∏neλi=e∑iλi=etrA,
the last step because trA=trJ is the sum of the eigenvalues. By property 1, expA=T(expJ)T−1, whence detexpA=detT⋅detexpJ⋅(detT)−1=detexpJ by the product rule for determinants. This is the route of (Riley, Hobson & Bence (2006), §8.16, Eq. (8.104), pp. 287-288), which runs it with a diagonalising similarity; the triangular form is what covers defective A as well. □
Method(Adjoints, unitarity and positivity)Original
For every n×n complex A, exp(A†)=(expA)†. Consequently:
if A is anti-Hermitian (A†=−A) then expA is unitary;
if A is Hermitian then expA is Hermitian, and it is positive definite: x†(expA)x>0 for every x=0.
Neither cited source states these; they are ours.
Idea of the proofOriginal
Conjugate transposition does two convenient things at once: it is continuous entrywise, being conjugation composed with a relabelling of indices, and it turns Ak into (A†)k. So it slides through the partial sums and hence through the limit, which is the whole of the first claim. The two consequences are then bookkeeping with the properties already proved: an anti-Hermitian exponent makes the adjoint equal to the inverse, and a Hermitian exponent lets the exponential be split as B†B with B invertible, which is what positive definiteness asks for.
ProofOriginal
Powers and the adjoint.(Ak)†=(A†)k for every k≥0, by induction: the case k=0 is I†=I, and (Ak+1)†=(AAk)†=(Ak)†A†=(A†)kA†=(A†)k+1, using (XY)†=Y†X† and then the inductive hypothesis.
The adjoint passes to the limit. Fix i,l. By definition (X†)il=Xli, so
by step 1. The outer series on the left converges to (expA)li=((expA)†)il, because complex conjugation is continuous, and the right-hand side converges to (exp(A†))il. Hence exp(A†)=(expA)†.
Anti-Hermitian exponent. If A†=−A, then by step 2 and property 3,
(expA)†(expA)=exp(A†)exp(A)=exp(−A)exp(A)=I,
and exchanging the two factors gives (expA)(expA)†=I in the same way. So expA is unitary.
Hermitian exponent. If A†=A then (expA)†=exp(A†)=expA by step 2, so expA is Hermitian. Also 21A is Hermitian, so B=exp(21A) satisfies B†=B by the same step, and B is invertible by property 3. Since 21A commutes with itself, property 2 gives
expA=exp(21A+21A)=B2=B†B.
Therefore x†(expA)x=x†B†Bx=∥Bx∥22 for every x, and for x=0 invertibility of B gives Bx=0, so the value is strictly positive: expA is positive definite. □
Method(When the product rule holds)Original
For n×n complex A,B:
AB=BA implies exp(sA+tB)=exp(sA)exp(tB) for all s,t∈C.
Conversely, if exp(sA+tB)=exp(sA)exp(tB) holds for all s,t in a neighbourhood of 0 in R2, then AB=BA.
Commuting is not necessary for the single identity exp(A+B)=exp(A)exp(B). For
A=(0002πi),B=(0012πi)
one has exp(A+B)=I=exp(A)exp(B) while [A,B]=(00−2πi0)=0. Consistently with part 2, the identity fails for this pair at other values of s,t, for instance at s=1, t=21.
4. For A and B small enough, exp(A)exp(B)=exp(A+B+21[A,B]+⋯), the Baker-Campbell-Hausdorff series, whose further terms are iterated commutators. The series converges only under such a smallness condition — and the matrix logarithm is not single-valued, as the remark below records — so for general A,B it is to be read as a formal expansion. What does hold for all A,B is its leading content: the discrepancy between the two sides is of second order in the size of A and B, with 21[A,B] as the second-order term, which is the error committed by operator-splitting schemes.
Idea of the proofOriginal
One direction is property 2 applied to the pair sA, tB, which commute as soon as A and B do. The converse is a matter of reading off coefficients: expand both sides to second order in (s,t) and the terms of order st differ by exactly 21[A,B], so an identity holding on a whole neighbourhood of the origin forces the commutator to vanish. The counterexample shows why "for all s,t" is not decoration: at one isolated pair of values the identity can hold by accident, because the scalar exponential is periodic along the imaginary axis and can hide a non-zero commutator.
ProofOriginal
Part 1. If AB=BA then (sA)(tB)=stAB=stBA=(tB)(sA) for all scalars s,t, so property 2 of the previous statement applies to the commuting pair sA, tB and gives exp(sA+tB)=exp(sA)exp(tB).
A second-order expansion. Let c=maxi,lmax(∣Ail∣,∣Bil∣) and σ=∣s∣+∣t∣, and put X=sA+tB, so maxi,l∣Xil∣≤σc. By the entrywise bound, applied to X, the entries of the tail expX−I−X−21X2=∑k≥3Xk/k! are bounded in modulus by
where O(σ3) means entries bounded by a constant times σ3 for σ≤1. Applying the same bound to exp(sA) and to exp(tB) separately and multiplying the two expansions — all discarded terms have total degree at least 3 in (s,t), and each factor has entries bounded by a constant for σ≤1, so they too are O(σ3) —
Part 2. Assume the identity holds for all (s,t) in a neighbourhood of the origin in R2. Subtracting the two displays of step 2, for all small (s,t)
0=st(21(AB+BA)−AB)+O(σ3)=−21st[A,B]+O(σ3).
Take s=t=h with h real and small, divide by h2 (so σ3/h2=8∣h∣→0), and let h→0: every entry of [A,B] vanishes, that is AB=BA.
Part 3, the counterexample. For M=(00βd) with d=0, induction gives Mk=(00βdk−1dk) for k≥1, so summing the series termwise (legitimate because the exponential series converges absolutely entry by entry),
expM=(10β(ed−1)/ded).
With d=2πi we have ed=1, so this is I: that covers B (take β=1, d=2πi) and A+B=(0014πi) (take β=1, d=4πi, again ed=1), while A is diagonal with expA=diag(1,e2πi)=I. Hence exp(A+B)=I=(expA)(expB). The commutator does not vanish:
At s=1, t=21 the formula with β=21 and d=3πi (so ed=−1) gives exp(A+21B)=(10i/(3π)−1), using −1/(3πi)=i/(3π), while exp(A)exp(21B)=exp(21B)=(10i/π−1) by the same formula with d=πi. The two differ, so for this pair exp(sA+tB)=exp(sA)exp(tB) is not an identity in s and t: the coincidence at (s,t)=(1,1) is isolated. Part 2 says more than this one computation does — it forces failure at points arbitrarily close to the origin, which (1,21) is not — but a single explicit pair of values is the cheapest way to see that the identity does not persist.
Part 4 is not proved here. The expansion of step 2 does give its leading content: replacing A by ϵA and B by ϵB and setting s=t=1,
exp(ϵA)exp(ϵB)−exp(ϵA+ϵB)=21ϵ2[A,B]+O(ϵ3),
which is the quantitative statement that the error of operator splitting is second order in the step size. The full Baker-Campbell-Hausdorff series, and the proof that all of its terms are iterated commutators, is beyond the scope of this unit; the standard reference is the theory of Lie algebras and the Lie group-Lie algebra correspondence. □
Method(Computing exp(tA))Original
If A is diagonalizable, A=SΛS−1 with Λ=diag(λ1,…,λn) and the columns of S an eigenbasis. Then by property 1,
exp(tA)=Sexp(tΛ)S−1=Sdiag(eλ1t,…,eλnt)S−1.
The same formula computes any f(A) with f(A)=Sdiag(f(λk))S−1, provided every λk lies inside the disc of convergence.
2. If A is normal, S may be taken unitary (see normal matrices), so S−1=S† and no inverse has to be computed; this is also the numerically stable case.
3. If A is defective, split it as A=D+N with DN=ND, D diagonalizable and N nilpotent (Nq=0 for some q≤n). Those three properties of D and N are all the step uses, however the splitting was obtained. One source of it is the Jordan normal form: if S−1AS=Λ+N with Λ its diagonal part and N its superdiagonal nilpotent part, then in the original basis D=SΛS−1 and N=SNS−1, so D is diagonalizable but not in general diagonal. Then by property 2
exp(tA)=exp(tD)exp(tN)=exp(tD)k=0∑q−1k!tkNk,
a finite sum, with exp(tD) computed by step 1. The entries of exp(tA) are then sums of terms (polynomial in t) ×eλt over the eigenvalues λ, the degree of the polynomial multiplying eλt being at most one less than the size of the largest Jordan block for λ — at most, because cancellation in the particular entry can lower it, and the degree is exactly that only for the entries where the top power of N survives.
4. For n small, or to avoid eigenvectors altogether, use that An is a combination of I,A,…,An−1 by the Cayley-Hamilton theorem — a matrix satisfies its own characteristic equation, that is, the characteristic polynomial p(λ)=det(λI−A) of the eigenvalue problem satisfies p(A)=0 ((Arfken, Weber & Harris (2013), §6.4, Exercise 6.4.1, Note, p. 317)) — so exp(tA)=∑j=0n−1αj(t)Aj for scalar functions αj. Determine the αj by requiring eλt=∑jαj(t)λj at each distinct eigenvalue, together with its first m−1 derivatives with respect to λ there when that eigenvalue has algebraic multiplicity m. That is n linear conditions on the n unknowns α0(t),…,αn−1(t), and what makes them determine the αj is that the coefficient matrix — a confluent Vandermonde matrix in the distinct eigenvalues — is nonsingular: if ∑jαjλj and its first m−1 derivatives vanish at an eigenvalue of multiplicity m, then the characteristic polynomial divides this polynomial of degree <n, which forces it to be 0, so the homogeneous system has only the zero solution.
5. Numerically, call a library routine (scipy.linalg.expm), which uses scaling and squaring with a rational approximation. Verify with property 6: detexp(tA) should equal ettrA.
How each route fails. Summing the defining series directly loses accuracy badly when ∥tA∥ is large, because the terms grow before they decay (and catastrophically for A with eigenvalues of both signs). Route 1 is useless for a defective A and unreliable for a nearly defective one, since S is then ill-conditioned and S−1 amplifies the error; see the later unit on condition numbers. Route 3 requires the Jordan structure, which is itself not computable stably. Never use exp(A)exp(B) in place of exp(A+B) unless [A,B]=0.
NoteOriginal
exp(tA) is the propagator of the linear system x˙=Ax: it carries the state at time 0 to the state at time t. When A is diagonalizable it does so by multiplying each eigendirection by eλt, so the eigenvalues alone read off the dynamics — decay for Reλ<0, growth for Reλ>0, oscillation at frequency Imλ. When A is defective that picture is incomplete: each Jordan block of size m contributes eλt times a polynomial in t of degree m−1, which is why A=(0010) has exp(tA)=(10t1) and grows without bound although its only eigenvalue is 0. The Traps below say what survives of the eigenvalue reading in that case. The same object appears as the time evolution e−iHt/ℏ of quantum mechanics and as the transition matrix of a continuous-time Markov chain.
Remark(Traps)Original
expA is not the matrix of entrywise exponentials. For A=(0010), expA=(1011), not (11e1).
exp is not surjective onto invertible matrices over R.detexpA=etrA>0 for real A, so no real matrix of negative determinant, such as diag(1,−1), is a real exponential. Over C every invertible matrix is an exponential.
Eigenvalues, not the spectral radius, govern growth.∥exp(tA)∥ is bounded as t→+∞ only if Reλ≤0 for every eigenvalue and every eigenvalue with Reλ=0 has equal algebraic and geometric multiplicity: A=(0010) has ρ(A)=0 yet exp(tA)=(10t1) grows.
The eigenvalues of expA are the eλ, but the map on eigenvalues is not injective.λ and λ+2πi give the same eλ, which is why there is no single-valued matrix logarithm and why the counterexample in part 3 of the product-rule statement exists at all.
Conventions. Physics writes e−iHt/ℏ and mathematics exp(tA); both are this object. A few texts define f(A) for non-diagonalizable A by the Jordan form or by a Cauchy integral rather than by a power series. The three definitions agree wherever more than one applies, but only the series version is used above.
Example(exp(tA) for a defective 3 x 3 matrix)Original
Take
A=00110−3013.
1. Eigenvalues.A is the companion matrix of the polynomial whose coefficients stand in its last row, so
det(λI−A)=λ3−3λ2+3λ−1=(λ−1)3,
one eigenvalue λ=1 of algebraic multiplicity 3.
2. A is defective, so route 1 is unavailable. Put
N=A−I=−1011−1−3012.
Its third row is −1 times the first plus 2 times the second, and its first two rows are independent, so rankN=2 and dimkerN=1. The eigenvalue 1 therefore has geometric multiplicity 1 against algebraic multiplicity 3: there is no eigenbasis, and A has a single Jordan block of size 3.
3. The splitting of step 3 of the method. Take D=I and N=A−I. Then D is diagonalizable, DN=ND because D=I, and N is nilpotent:
N2=111−2−2−2111,N3=0,
all three rows of N2 being equal, every column of N2 is a multiple of (1,1,1)T, and N(1,1,1)T=(−1+1+0,0−1+1,1−3+2)T=0, so N3=N⋅N2=0.
4. Assemble. Since tD=tI commutes with tN, property 2 gives exp(tA)=et(I+tN+21t2N2), a finite sum:
Every entry is a quadratic polynomial times et, the degree being one less than the size of the Jordan block, exactly as step 3 of the method predicts.
5. Two checks. At t=0 the matrix is I. Differentiating, the derivative of et(I+tN+21t2N2) is et(I+tN+21t2N2)+et(N+tN2), which at t=0 is I+N=A, as property 5 requires.
Pulling the factor e/2 out of all three rows multiplies the determinant by (e/2)3, and expanding the integer matrix along its first row,
det1130−2−8137=1⋅(−14+24)−0+1⋅(−8+6)=10−2=8,
so detexpA=(e/2)3⋅8=e3=20.085537, which is etrA with trA=0+0+3=3: property 6 confirmed. Note in passing that (expA)12=0 although A12=1, and that the entrywise exponential of A would have no zero entry at all: the exponential redistributes the entries completely.
7. What it solves.x˙=Ax for this companion matrix is the scalar equation y′′′−3y′′+3y′−y=0 written for x=(y,y′,y′′)T, whose general solution is (c0+c1t+c2t2)et. The first column of exp(tA) must therefore be the solution with y(0)=1, y′(0)=y′′(0)=0, namely y(t)=et(1−t+21t2); differentiating gives y′(t)=21t2et and y′′(t)=(t+21t2)et, which are precisely the second and third entries of that column. The matrix exponential has solved the third-order equation without any eigenvector being computed.
ApplicationOriginal
Unitarity of quantum time evolution. For a time-independent Hamiltonian, the state obeys iℏψ˙=Hψ, so ψ(t)=exp(−iHt/ℏ)ψ(0) by property 5. If H is Hermitian then (−iHt/ℏ)†=+iHt/ℏ=−(−iHt/ℏ) for real t, so the exponent is anti-Hermitian and the adjoints statement makes the propagator unitary. That is exactly the conservation of total probability, ∥ψ(t)∥=∥ψ(0)∥, and it is a consequence of the algebra of the exponential rather than of anything about the particular H. The eigenvalue property 4 then says that an energy eigenstate only picks up the phase e−iEt/ℏ, which is why stationary states are stationary.
Stability of a linear system, read off the spectrum. For x˙=Ax the solution is x(t)=exp(tA)x0, and property 1 together with step 3 of the method says that every entry of exp(tA) is a sum of terms, each a polynomial in t times eλt, one group of terms for each eigenvalue λ. So all solutions decay when every eigenvalue has Reλ<0, and the trap noted above — a purely imaginary eigenvalue whose Jordan block is larger than 1×1 produces secular growth like teiωt — is the mathematics behind resonance in a driven oscillator and behind the marginal cases in the stability analysis of a steady state.
Operator splitting in numerical simulation. Split the generator of an evolution equation into parts that can each be exponentiated exactly: a Laplacian handled in Fourier space and a nonlinear or reaction term handled pointwise, or the kinetic and potential parts of a Hamiltonian. The scheme advances by exp(hA)exp(hB) in place of exp(h(A+B)), and part 5 of the product-rule proof shows the error per step is 21h2[A,B]+O(h3). This is the error estimate behind the split-step Fourier method for the nonlinear Schrodinger equation, Strang splitting in atmosphere and ocean models, and the leapfrog and higher-order symplectic integrators of molecular dynamics, where symmetrising the composition as exp(21hA)exp(hB)exp(21hA) cancels the commutator term and leaves an O(h3) error per step.
Rotations from angular velocity. A real antisymmetric 3×3 matrix Ω with Ωx=ω×x satisfies ΩT=−Ω, so exp(tΩ) is real orthogonal by the adjoints statement, and detexp(tΩ)=ettrΩ=1 by property 6: it is a rotation, through angle ∣ω∣t about ω. Rigid-body kinematics is the statement R˙=ΩR, solved by this exponential, and the failure of exp(A)exp(B)=exp(A+B) for non-commuting A,B is the familiar fact that rotations about different axes do not commute.
Continuous-time Markov chains. A generator Q has non-negative off-diagonal entries and rows summing to zero, and the transition matrix over time t is P(t)=exp(tQ). Both defining properties of a stochastic matrix come straight out of this unit: Q1=0 gives Qk1=0 for k≥1 and hence P(t)1=1, so the rows sum to one; and writing exp(tQ)=e−ctexp(t(Q+cI)) with c=maxi∣Qii∣, permitted by property 2 because cI commutes with everything, exhibits P(t) as a positive multiple of the exponential of a matrix with non-negative entries. For t≥0 — the only range in which a transition matrix is asked for — every term tk(Q+cI)k/k! of that exponential series is then entrywise non-negative, and so is P(t). The Chapman-Kolmogorov relation P(s+t)=P(s)P(t) is property 2 again.
Defines induced matrix norms and the spectral norm, defines the condition number and expresses it through the singular values, and states how far the solution of a linear system can move when the data are perturbed.
Fix one norm ∥⋅∥, given on Ck for every k and so available on the domain and on the codomain alike; as in (Stakgold & Holst (2011), §10.6, p. 772), where the operator acts on a single space, it is the same norm on both sides. The induced (or operator) norm of an m×n complex matrix A is
∥A∥=x=0max∥x∥∥Ax∥=∥x∥=1max∥Ax∥,
the maximum being attained because the unit sphere is compact and x↦∥Ax∥ continuous (compactness). It is the smallest number c with ∥Ax∥≤c∥x∥ for all x. Because one and the same norm is used on both spaces, induced norms satisfy ∥I∥=1 and are submultiplicative: ∥AB∥≤∥A∥∥B∥ whenever the product is defined. Both properties need that agreement: if different norms are allowed on the two spaces, ∥⋅∥∞ on the domain and ∥⋅∥1 on the codomain give ∥In∥=n, and ∥⋅∥1 on the domain with ∥⋅∥∞ on the codomain give ∥M∥=maxi,j∣Mij∣, for which the 2×2 all-ones A has ∥A2∥=2>1=∥A∥2.
Taking the vector norms ∥x∥q=(∑i∣xi∣q)1/q for 1≤q<∞, and ∥x∥∞=maxi∣xi∣, gives
and the spectral norm∥A∥2, induced by the Euclidean norm. The Frobenius norm∥A∥F=(∑i,j∣Aij∣2)1/2 is submultiplicative but is not an induced norm, since ∥In∥F=n=1 for n≥2. It is not induced even by a pair of norms, one on each space — but that value n does not show it, being itself ∥In∥ for the pair (∥⋅∥∞ on the domain, ∥⋅∥2 on the codomain). What shows it is the behaviour on rank-one matrices; see step 0 of the proof below.
Definition(The spectral norm and the spectrum)AdaptedAdapted from Stakgold & Holst (2011), §10.6, pp. 771-772 (including Lemma 10.6.1)
Let A be a complex matrix and ρ(A)=max{∣λ∣:λ∈σ(A)} its spectral radius (defined for square A).
If A is square then ρ(A)≤∥A∥ for every induced norm.
If A is Hermitian then ∥A∥2=ρ(A)=maxk∣λk∣, and λmin(A)≤x†Ax/x†x≤λmax(A) for x=0.
Inequality 2 can be strict: A=(0010) has ρ(A)=0 and ∥A∥2=1.
Idea of the proofOriginal
Everything here comes from one change of coordinates. In the singular bases A is a diagonal matrix, so its effect on the unit sphere is a list of independent scalings, and the norm is the largest of them. The spectral radius, by contrast, only tests the single direction of an eigenvector, so it can never exceed a quantity defined as a maximum over all directions; it misses the stretching that happens off the eigenvectors, which is why a nilpotent matrix can have radius zero and norm one. When A is Hermitian there is no off-eigenvector stretching left to miss, and the two agree.
ProofAdaptedAdapted from Stakgold & Holst (2011), §10.6, pp. 771-772, including Lemma 10.6.1
The claims of the definition.∥I∥=1 because ∥Ix∥=∥x∥ for every x — which is a statement about one norm, and is where the convention that the same norm is used on both spaces is needed. Submultiplicativity: ∥ABx∥≤∥A∥∥Bx∥≤∥A∥∥B∥∥x∥, by applying the defining inequality ∥My∥≤∥M∥∥y∥ first to M=A with y=Bx and then to M=B; here too the intermediate vector Bx is measured in the one norm shared by the two spaces. For the 1-norm, ∥Ax∥1=∑i∑jAijxj≤∑j∣xj∣∑i∣Aij∣≤(maxj∑i∣Aij∣)∥x∥1 by the triangle inequality, with equality at x=ej∗ for a column j∗ of largest sum. For the ∞-norm, ∥Ax∥∞=maxi∑jAijxj≤(maxi∑j∣Aij∣)∥x∥∞, with equality at the x having ∥x∥∞=1 and xj=Ai∗j/∣Ai∗j∣ (or 1 where that entry vanishes) for a row i∗ of largest sum, since then (Ax)i∗=∑j∣Ai∗j∣. Finally ∥In∥F=n=1 for n≥2, and every induced norm has ∥I∥=1, so the Frobenius norm is not an induced norm.
Not even a pair of norms induces ∥⋅∥F, for n≥2. This last argument is ours; the cited page does not discuss the Frobenius norm. Suppose some norm a on the domain Cn and some norm b on the codomain Cn had max∥x∥a=1∥Ax∥b=∥A∥F for every A. Write ∥v∥a∗=max∥x∥a=1∣v†x∣, a norm on Cn. Test the hypothesis on the rank-one matrices A=uv†: since (uv†)x=(v†x)u,
the second because ∑i,j∣uivj∣2=∥u∥22∥v∥22. So ∥u∥b∥v∥a∗=∥u∥2∥v∥2 for all u,v. Fixing a v0 with ∥v0∥2=1 and putting c=∥v0∥a∗>0 gives ∥u∥b=∥u∥2/c for every u; feeding that back with any fixed u=0 gives ∥v∥a∗=c∥v∥2 for every v. Now take any x with ∥x∥a=1 and use v=x in the definition of ∥⋅∥a∗: c∥x∥2=∥x∥a∗≥∣x†x∣=∥x∥22, so ∥x∥2≤c on the whole unit sphere of a. Hence
∥x∥a=1max∥Inx∥b=∥x∥a=1maxc∥x∥2≤1,
whereas the hypothesis demands the value ∥In∥F=n>1. No such pair exists. (This is the argument the definition block promises; the value ∥In∥F=n alone does not settle the question, because n is the pair-induced norm of In for ∥⋅∥∞ on the domain and ∥⋅∥2 on the codomain.)
the first equality again by unitarity. Equality holds at x=v1, where y=e1. Hence ∥A∥2=s1, and s12 is the largest eigenvalue of A†A by part 1 of the existence statement for the SVD, so ∥A∥2=λmax(A†A).
Part 2. Let λ be an eigenvalue of A with ∣λ∣=ρ(A) and x=0 an eigenvector. For any induced norm, ∥A∥≥∥Ax∥/∥x∥=∥λx∥/∥x∥=∣λ∣ by absolute homogeneity of the vector norm, so ∥A∥≥ρ(A).
Part 3. A Hermitian matrix is normal, so by part 5 of the SVD existence statement its singular values are the numbers ∣λk∣ sorted decreasingly. With part 1 this gives ∥A∥2=maxk∣λk∣=ρ(A), which is the matrix case of Lemma 10.6.1 of the cited section. The two-sided Rayleigh estimate λmin(A)≤x†Ax/x†x≤λmax(A) for x=0 is the Rayleigh bounds.
Part 4. For A=(0010) the characteristic polynomial is λ2, so σ(A)={0} and ρ(A)=0. But A†A=diag(0,1), whose largest eigenvalue is 1, so ∥A∥2=1 by part 1. The inequality of part 2 is therefore strict here, and it cannot be reversed in general.
A lemma for the blocks that follow: the inverse. Nothing in parts 1-4 needs this step; it is recorded here, where the singular bases are at hand, because the condition number defined in the next block and the sensitivity bounds after it both use it. Let A be square and invertible, so p=n and sn>0. Inverting A=USV† gives A−1=VS−1U†, and reversing the order of the columns of V and of U turns this into an SVD of A−1 with singular values 1/sn≥⋯≥1/s1. So ∥A−1∥2=1/sn by part 1, and
κ2(A)=∥A∥2∥A−1∥2=sns1,
which is the formula asserted in the definition below. If A is Hermitian positive definite its singular values are its eigenvalues, all positive — by part 5 of the SVD existence statement (the singular value decomposition), as in step 3 above, with the absolute values dropped because the eigenvalues are positive — so κ2(A)=λmax/λmin. A unitary Q has Q†Q=I, so every eigenvalue of Q†Q is 1 and every singular value of Q is 1, giving κ2(Q)=1. Two further consequences used below: applying A−1 to Avi=siui gives A−1ui=vi/si, and the computation of part 1 read from the other end gives ∥Ax∥22=∑isi2∣yi∣2≥sn2∥x∥22, that is ∥Ax∥2≥sn∥x∥2 for every x. □
For an invertible square matrix A and an induced norm ∥⋅∥, the condition number is
κ(A)=∥A∥∥A−1∥≥∥AA−1∥=∥I∥=1.
For a singular A one sets κ(A)=∞. In the spectral norm, with singular values s1≥⋯≥sn>0,
κ2(A)=sns1,
and if A is Hermitian positive definite this is λmax/λmin. A matrix with κ close to 1 is called well-conditioned, and one with κ large ill-conditioned. Unitary and orthogonal matrices have κ2=1.
What the cited page displays is the more general A-condition numberκA(M)=∥M∥A∥M−1∥A, built from the norm of an A-inner product and introduced there as an extension of "the standard notion"; the κ above is that standard notion, the case in which the inner product is the Euclidean one.
Definition(Sensitivity of a linear system)Original
Let A be an invertible n×n matrix, b=0, and x the solution of Ax=b. Fix any induced norm and write κ=κ(A).
Perturbed right-hand side. If A(x+δx)=b+δb, then
κ1∥b∥∥δb∥≤∥x∥∥δx∥≤κ∥b∥∥δb∥.
In the spectral norm the upper bound is attained. So κ is exactly the worst-case amplification of relative error, and no sharper constant exists.
2. Perturbed matrix. If A+δA is invertible and (A+δA)x^=b, then
∥x^∥∥x^−x∥≤κ∥A∥∥δA∥.
Distance to singularity. In the spectral norm,
min{∥δA∥2:A+δA is singular}=sn=κ2(A)∥A∥2,
the minimum being attained at δA=−snunvn†. An ill-conditioned matrix is therefore literally a matrix that is relatively close to a singular one.
4. Rule of thumb. Solving Ax=b in arithmetic with unit roundoff ε gives a relative error of order κ(A)ε: about log10κ(A) of the available decimal digits are lost, whatever algorithm is used.
Idea of the proofOriginal
Both bounds are the defining inequality of an induced norm used twice, once for A and once for A−1 — and the condition number is precisely the product of those two uses. Which way round you use them decides whether you get the upper or the lower bound, which is why the true amplification is trapped between κ and 1/κ and not merely below κ. Sharpness is then a matter of aiming: put the data along the direction A stretches most and the error along the one it stretches least. The distance to singularity is the same arithmetic read geometrically: A moves every vector by at least sn, so no smaller perturbation can manufacture a null vector, and one of size exactly sn can.
ProofOriginal
Write ∥⋅∥ for the fixed induced norm and κ=∥A∥∥A−1∥. Note x=0, since b=Ax=0.
Part 1, upper bound. Subtracting Ax=b from A(x+δx)=b+δb gives Aδx=δb, so δx=A−1δb and ∥δx∥≤∥A−1∥∥δb∥. Separately ∥b∥=∥Ax∥≤∥A∥∥x∥, which rearranges to 1/∥x∥≤∥A∥/∥b∥. Multiplying the two inequalities,
∥x∥∥δx∥≤∥A−1∥∥δb∥⋅∥b∥∥A∥=κ∥b∥∥δb∥.
Part 1, lower bound. Run the same two steps the other way round: δb=Aδx gives ∥δb∥≤∥A∥∥δx∥, and x=A−1b gives ∥x∥≤∥A−1∥∥b∥, that is 1/∥b∥≤∥A−1∥/∥x∥. Multiplying,
∥b∥∥δb∥≤∥A∥∥δx∥⋅∥x∥∥A−1∥=κ∥x∥∥δx∥,
which is the lower bound after dividing by κ>0.
Part 1, the upper bound is attained in the spectral norm. Let A=USV† be an SVD. Take b=s1u1; then x=A−1b=s1⋅v1/s1=v1 by step 5 of the previous proof, so ∥x∥2=1 and ∥b∥2=s1. Take δb=ϵun with ϵ=0; then δx=A−1δb=ϵvn/sn, so ∥δx∥2=∣ϵ∣/sn. Hence
so no constant smaller than κ2 can replace it in part 1. (For n=1 the two directions coincide, but then κ2=1 and the bound is an equality for every δb anyway.)
Part 2. Subtracting Ax=b from (A+δA)x^=b gives Ax^+δAx^=Ax, hence A(x^−x)=−δAx^ and x^−x=−A−1δAx^. Taking norms, ∥x^−x∥≤∥A−1∥∥δA∥∥x^∥, and x^=0 because (A+δA)x^=b=0. Dividing by ∥x^∥ and inserting ∥A∥/∥A∥,
∥x^∥∥x^−x∥≤∥A−1∥∥A∥⋅∥A∥∥δA∥=κ∥A∥∥δA∥.
The bound is in terms of ∥x^∥, the computed vector, not ∥x∥; that is what makes it an identity-free consequence of one line of algebra.
Part 3, no singular matrix is nearer than sn. Suppose ∥δA∥2<sn. For every x=0, the triangle inequality and ∥Ax∥2≥sn∥x∥2 (step 5 of the previous proof) give
∥(A+δA)x∥2≥∥Ax∥2−∥δAx∥2≥(sn−∥δA∥2)∥x∥2>0,
so A+δA has trivial kernel and is invertible. Hence every δA making A+δA singular has ∥δA∥2≥sn.
Part 3, and one at exactly that distance. Put δA=−snunvn†. For any x, ∥unvn†x∥2=∣vn†x∣∥un∥2=∣vn†x∣≤∥x∥2 by Cauchy-Schwarz, with equality at x=vn, so ∥δA∥2=sn. And (A+δA)vn=snun−snun=0 with vn=0, so A+δA is singular. The minimum is therefore sn, attained there, and sn=s1/κ2(A)=∥A∥2/κ2(A) by step 5 of the previous proof.
Part 4 is a rule of thumb, and the only part of this statement not proved. Its basis is part 2 together with backward stability: a backward-stable solver returns an x^ that is the exact solution of a nearby system, (A+δA)x^=b with ∥δA∥/∥A∥ of order ε (times a modest factor depending on n and on the algorithm). Part 2 then bounds the relative error of x^ by κ(A) times that, and log10 of the bound is the count of lost decimal digits. The clause "whatever algorithm is used" rests on step 3: merely rounding b into the machine perturbs it by a relative ε, and for the worst b and δb that perturbation alone is amplified by exactly κ, so no algorithm can be uniformly better. What the rule of thumb adds to this is the empirical claim that the orders of magnitude are typically attained, which is not a theorem. □
NoteOriginal
A matrix norm measures how much a matrix can stretch a vector; the condition number measures how much worse the stretching can be in one direction than in another. That ratio, not the size of the determinant, is what decides whether a linear system can be solved from imperfect data: κ(A) is the factor by which relative errors in b or in A may be magnified in the answer, and the quantity to inspect before trusting a numerical solution or a fitted parameter.
A small determinant does not mean ill-conditioned, and a large one does not mean well-conditioned.εIn has determinant εn, as small as you like, yet κ2=1: it is perfectly conditioned, merely scaled. Conversely diag(106,10−6) has determinant 1 and κ2=1012. Conditioning is about the ratio of the extreme singular values, which is invariant under scaling A↦cA, while the determinant is not.
Conditioning is a property of the problem, not of the algorithm. Strauss (§1.5, pp. 26-27) places it as the stability requirement in well-posedness: a nonsingular matrix with one very small eigenvalue gives a unique solution that nevertheless moves a great deal when b is perturbed slightly, and no choice of solver can repair that. The sharp condition is a small smallest singular value, which is the weaker requirement: sn≤∣λ∣ for every eigenvalue λ, since ∥Ax∥2≥sn∥x∥2 applied to an eigenvector gives ∣λ∣≥sn, so a small eigenvalue forces a small sn but not conversely. A standard example is the Hilbert matrix Aij=1/(i+j+1) with i,j=0,…,n−1, which arises when a function is fitted by a power series in the mean-square sense and whose determinant and smallest singular value shrink extremely fast with n (Arfken, Weber and Harris note detA=3.7×10−12 already at n=5, Exercise 5.1.3, p. 267).
κ depends on the norm. Different induced norms give different values, though they are equivalent up to factors depending only on n (for example n−1κ∞≤κ2≤nκ∞), so "ill-conditioned" is norm-independent as a qualitative statement. Always say which norm you mean; κ with no subscript usually means κ2.
ρ(A) is not a norm. It fails the triangle inequality and vanishes on non-zero nilpotent matrices, so it cannot be used in the bounds above; inequality 2 of the spectral-norm statement runs one way only.
Non-square matrices.κ2(A)=s1/sr with r=rankA is also used for rectangular A and governs least-squares sensitivity, but then ∥A∥2∥A+∥2 replaces ∥A∥∥A−1∥ and the bounds of parts 1-2 need modification.
Example(A 2 x 2 matrix with condition number exactly 100)Original
Take
A=(1019999101).
1. Spectrum and singular values.A is real symmetric, and A(1,1)T=(200,200)T, A(1,−1)T=(2,−2)T, so the eigenvalues are 200 and 2 with orthonormal eigenvectors
w1=21(1,1)T,w2=21(1,−1)T.
Both eigenvalues are positive, so A is symmetric positive definite and its singular values are its eigenvalues: s1=200, s2=2, with u1=v1=w1 and u2=v2=w2.
2. The norms. By part 1 of the spectral-norm statement, ∥A∥2=200. The largest column sum and the largest row sum are both 101+99=200, so ∥A∥1=∥A∥∞=200 as well, and ∥A∥F=2(1012+992)=40004=210001≈200.0100.
3. The inverse and the condition numbers.detA=1012−992=(101−99)(101+99)=400, so
A−1=4001(101−99−99101).
Then ∥A−1∥2=1/s2=21 and ∥A−1∥1=(101+99)/400=21, giving
κ2(A)=s2s1=2200=100,κ1(A)=200⋅21=100.
Note that detA=400 is not small: conditioning and the determinant are unrelated, as the remark above says.
4. The worst case is attained. Take b=(200,200)T, which is a multiple of u1; then x=(1,1)T, since Ax=(200,200)T. Perturb along u2: δb=(0.1,−0.1)T. Then
a ratio of exactly 100=κ2(A): a perturbation of the data in the fourth decimal place has moved the answer in the second.
5. The best case, for contrast. Keep the same b and the same ∥δb∥2, but point the perturbation along u1 instead: δb=(0.1,0.1)T. Then δx=4001(0.2,0.2)T=(5×10−4,5×10−4)T and ∥δx∥2/∥x∥2=5×10−4, an amplification of 1. The two experiments differ only in the direction of δb, and part 1 confines every intermediate case between the factors 1/100 and 100.
6. Distance to singularity. Part 3 says the nearest singular matrix is at spectral distance s2=2 and is reached at δA=−s2u2u2T. Here
which is indeed singular, and ∥δA∥2=2. Relative to ∥A∥2=200, this matrix is only 1%=1/κ2(A) away from being singular — the geometric meaning of κ2=100.
7. Digits.log10κ2(A)=2, so part 4 predicts the loss of about two decimal digits. In IEEE double precision, where the unit roundoff is ε≈1.1×10−16, a solution of Ax=b should be trusted to a relative accuracy of about 10−14 and no better.
ApplicationOriginal
Never form the normal equations if you can avoid it. The singular values of A†A are the squares of those of A, so κ2(A†A)=κ2(A)2. Solving a least-squares problem through A†Ax=A†b therefore squares the condition number and doubles the number of decimal digits lost — from log10κ2(A) to 2log10κ2(A), so that double precision and κ2(A)=104 leave about eight good digits instead of twelve — while a QR or SVD factorisation of A itself works with κ2(A). This single inequality is why numerical libraries fit models by orthogonal factorisation rather than by the textbook normal equations.
Fitting a polynomial, and why high degree fails. Least-squares fitting of a function on [0,1] by a power series leads to the Hilbert matrix Aij=1/(i+j+1) with i,j=0,…,n−1, whose singular values decay geometrically: κ2≈4.8×105 already at n=5 and ≈1.6×1013 at n=10. At n=10 double precision leaves about three correct digits in the fitted coefficients, whatever the solver. The cure is not a better solver but a better basis — orthogonal polynomials, for which the corresponding Gram matrix is diagonal, and the identity once they are normalised, so that κ2=1. Without normalisation κ2 is the ratio of the extreme squared norms, which stays modest: for the Legendre polynomials, with ∫−11Pk2=2/(2k+1), degrees up to 10 give κ2=21.
Discretised differential operators get worse as the mesh is refined. The standard second-difference approximation of −d2/dx2 on [0,1] with Dirichlet conditions and N interior points, spacing h=1/(N+1), is symmetric positive definite with eigenvalues (4/h2)sin2(kπh/2), so
κ2=cot2(2πh)≈π2h24,
which is 4.13×103 at N=100 and grows like N2. Accuracy and conditioning therefore pull in opposite directions, the iteration count of conjugate gradients grows like κ2∼N, and this is the reason multigrid and other preconditioners exist: they replace the system by an equivalent one whose condition number does not blow up with h.
Inverse problems in the geosciences and in imaging. Recovering a source distribution from remote measurements — seismic velocity from travel times, atmospheric composition from radiances, an image from projections — discretises to an Ax=b whose singular values decay towards zero, so κ2 is enormous and the relation of part 1 says the inversion will amplify measurement noise by that factor. The practical response is to accept a different problem: truncate the singular value expansion or add a Tikhonov penalty, replacing A by a nearby operator of moderate condition number and trading resolution for stability. The condition number is the number that quantifies the trade.
Trusting an eigenvector computation. Route 1 of the matrix exponential computes f(A)=Sf(Λ)S−1 from an eigenbasis. The relevant conditioning is not κ(A) but κ(S), which blows up as A approaches a defective matrix, because two eigenvectors then become nearly parallel. The same quantity controls the sensitivity of computed eigenvalues, which is why nearly degenerate normal modes of a structure, or nearly crossing energy levels, demand care that well-separated ones do not.