Companion page: Chapter 5 — Eigenvalue problems — the notes that this page goes with.
Commands for the material of Chapter 5 so far: eigenvalues and eigenvectors, the characteristic polynomial, trace and determinant, the two multiplicities, diagonalization, powers, commuting matrices, the Cayley—Hamilton theorem, difference equations, Hermitian and symmetric matrices, unitary matrices, and similar matrices. Everything below also runs in GNU Octave, and the places where the two differ are called out. The running example is the \(2\times 2\) matrix of the notes,
A = [4 -5; 2 -3];Eigenvalues and eigenvectors
eig(A) % the eigenvalues: 2 and -1
[V, D] = eig(A) % eigenvectors in the columns of V, eigenvalues on the diagonal of D
norm(A*V - V*D) % of order 1e-16: A*v = lambda*v, column by columneig scales each eigenvector to length 1, so the columns of V are \((5,2)^T/\sqrt{29}\) and \((1,1)^T/\sqrt 2\) rather than the \((5,2)^T\) and \((1,1)^T\) of the notes. An eigenvector is fixed only up to a scalar multiple, so both are right. To see the notes' vectors, rescale every column so that its first entry is 1:
V ./ V(1,:) % one column is (1, 0.4) = (5,2)/5, the other (1, 1)
diag(D).' % the eigenvalue of each column, in the same orderThe order in which eig returns the eigenvalues is not specified, and MATLAB and Octave do not always agree on it (the projection below is an example). Match eigenvalues and eigenvectors through D and V, never by position alone.
The characteristic polynomial
c = poly(A) % [1 -1 -2], the coefficients of lambda^2 - lambda - 2
roots(c) % 2 and -1poly returns the coefficients of \(\det(\lambda I - A)\), whose leading coefficient is \(1\). The notes use \(p_A(\lambda) = \det(A - \lambda I)\), which is \((-1)^n\) times it: the same polynomial when \(n\) is even, the negative when \(n\) is odd. The roots are the same either way.
Triangular matrices
U = [1 4 5; 0 2 6; 0 0 3];
eig(U) % 1, 2, 3: the diagonal entries
[L, U2] = lu(A);
diag(U2) % the pivots of A: 4 and -0.5, which are not its eigenvaluesProjections
a = [1; 1];
P = a*a'/(a'*a); % projection onto the line through a
eig(P) % 0 and 1
[VP, DP] = eig(P) % (-1,1)/sqrt(2) belongs to 0, (1,1)/sqrt(2) to 1; see diag(DP) for the order
P4 = diag([1 0 0 1]);
eig(P4) % 0, 0, 1, 1
M = [1 0; 1 1; 1 2];
PM = M*((M'*M)\M'); % projection onto the column space of M
eig(PM) % 1, 1 and 0, up to round-offIn the last line MATLAB lists \(1, 0, 1\) and Octave lists \(0, 1, 1\), and the "0" comes out as a number of order \(10^{-17}\) rather than exactly zero.
A rotation, and complex numbers
R = [0 -1; 1 0];
eig(R) % 0 + 1i and 0 - 1i
[VR, DR] = eig(R);
VR ./ VR(2,:) % columns (1i, 1) and (-1i, 1): the notes' (j, 1), and (j, -1) times -1
diag(DR).' % the matching eigenvalues, j and -jMATLAB and Octave accept 1j as well as 1i on input (z = 3 + 4j), but both print the imaginary unit as i.
For complex vectors, x' is the conjugate transpose \(x^H\) and x.' is the plain transpose \(x^T\). The difference is the one Chapter 3 makes about lengths in \(\mathbb{C}^n\):
x = [1; 1j];
x'*x % 2, the squared length
x.'*x % 0: without the conjugate a nonzero vector would get zero lengthTrace and determinant
[trace(A) sum(eig(A))] % 1 and 1
[det(A) prod(eig(A))] % -2 and -2
[trace(R) sum(eig(R))] % 0 and 0
[det(R) prod(eig(R))] % 1 and 1Algebraic and geometric multiplicity
eig lists each eigenvalue as often as its algebraic multiplicity. The geometric multiplicity is the dimension of the null space of \(A - \lambda I\), which null returns as a basis, one vector per column:
N = [0 1; 0 0];
eig(N) % 0, 0: algebraic multiplicity 2
null(N) % one column, (-1, 0): geometric multiplicity 1
size(null(N), 2) % 1
2 - rank(N) % 1, the same count as n - rank(A - lambda*I)eig still returns two columns in V here, but they are parallel to within a negligible error:
[VN, DN] = eig(N) % second column is (-1, about 2e-292)
rank(VN) % 1: rank uses a tolerance, and the columns agree to within itThe other two examples of the notes:
n = 4;
size(null(eye(n) - 1*eye(n)), 2) % 4: every nonzero vector is an eigenvector of I
B = diag([1 0 0 0]);
size(null(B), 2) % 3, the geometric multiplicity of 0
size(null(B - eye(4)), 2) % 1, the geometric multiplicity of 1Diagonalization
A = [4 -5; 2 -3];
[S, L] = eig(A); % S: eigenvectors in its columns, L: the diagonal matrix Lambda
norm(S*L/S - A) % of order 1e-16: A = S*Lambda*inv(S)
S\A*S % Lambda again, up to round-off, in the order eig choseS*L/S is MATLAB’s way of writing \(S\Lambda S^{-1}\) without forming the inverse, and S\A*S is \(S^{-1}AS\). The rotation of the notes is diagonalizable over the complex numbers:
R = [0 -1; 1 0];
[SR, LR] = eig(R); % complex S and Lambda
norm(SR*LR/SR - R) % 0, or of order 1e-16Powers and the inverse
norm(A^5 - S*L^5/S) % of order 1e-15: A^k = S*Lambda^k*inv(S)
[sort(eig(inv(A))) sort(1./eig(A))] % both columns -1 and 0.5: reciprocals of 2 and -1Products of matrices
A1 = [0 1; 0 0]; B1 = [0 0; 1 0];
eig(A1).' % 0 0
eig(B1).' % 0 0
sort(eig(A1*B1)).' % 0 and 1: not products of the eigenvalues of A1 and B1
B = A^2 + 3*A; % a polynomial in A commutes with A
norm(A*B - B*A) % 0, or of order 1e-15
S\B*S % diagonal: the S that diagonalizes A also diagonalizes BThe Cayley—Hamilton theorem
polyvalm evaluates a polynomial at a matrix, in the sense of the notes: each power of \(\lambda\) becomes the same power of the matrix and the constant term becomes a multiple of eye.
polyvalm(poly(A), A) % the zero matrix: A^2 - A - 2I = 0
N = [0 1; 0 0];
polyvalm(poly(N), N) % zero as well, although N is not diagonalizable(poly returns the coefficients of \(\det(\lambda I - A)\), which is \((-1)^n p_A(\lambda)\); the sign makes no difference to a polynomial that evaluates to zero.)
Fibonacci numbers and stability
F = [1 1; 1 0];
u = [1; 0]; % u_0 = (f_1, f_0)
for k = 1:10
u = F*u; % u_k = (f_(k+1), f_k)
end
u(2) % 55 = f_10
eig(F).' % 1.6180 and -0.6180, in some order
phi = (1 + sqrt(5))/2;
round(phi^10/sqrt(5)) % 55 again, from the formula of the notesFloating point stores whole numbers exactly only up to \(2^{53} \approx 9\times 10^{15}\). The matrix route stays exact until \(f_{79}\) passes that limit. The formula route goes wrong earlier, at \(f_{71}\), because \(\phi\) and its powers are themselves rounded:
format long g
u50 = F^50*[1; 0]; u50(2) % 12586269025, exact
round(phi^50/sqrt(5)) % 12586269025
u80 = F^80*[1; 0]; u80(2) % close to, but not exactly, f_80 = 23416728348467685
format shortFor \(f_{100}\) exact integer arithmetic is needed. In MATLAB the Symbolic Math Toolbox provides it (this block is MATLAB only; Octave would need its separate symbolic package):
u = sym([1 1; 1 0])^100*[1; 0];
u(2) % 354224848179261915075 = f_100, exactStability is read off the sizes of the eigenvalues:
abs(eig(F)).' % 1.6180 and 0.6180: one is larger than 1, so unstable
M = [0.5 0.2; 0.1 0.3];
abs(eig(M)).' % 0.5732 and 0.2268: both smaller than 1, so stable
norm(M^50) % about 9e-13: M^k u_0 goes to zero for every u_0Hermitian and symmetric matrices
H = [2 1-1j; 1+1j 3];
ishermitian(H) % 1 (true): H' equals H, where ' is the conjugate transpose
diag(H).' % 2 3: the diagonal of a Hermitian matrix is real
x = [1+2j; 3-1j];
x'*H*x % 28, a real number, as Theorem 5.16 says
[V, D] = eig(H);
diag(D).' % 1 and 4: real eigenvalues
norm(V'*V - eye(2)) % of order 1e-16: the unit eigenvectors are orthonormalA real symmetric matrix with a repeated eigenvalue still has orthonormal eigenvectors, which eig returns:
A = [2 1 1; 1 2 1; 1 1 2];
[Q, L] = eig(A);
diag(L).' % 1, 1 and 4: the eigenvalue 1 is repeated
norm(Q'*Q - eye(3)) % of order 1e-16: Q is an orthogonal matrix
norm(Q*L*Q' - A) % of order 1e-15: A = Q*Lambda*Q'
A3 = L(1,1)*Q(:,1)*Q(:,1)' + L(2,2)*Q(:,2)*Q(:,2)' + L(3,3)*Q(:,3)*Q(:,3)';
norm(A3 - A) % of order 1e-15: the sum of rank-one matrices, (17) of the notesThe Hilbert matrix of Chapter 3 is symmetric, and its eigenvalues fall off quickly, so two terms of (17) of the notes already approximate it well. sort puts the eigenvalues in order of size, whatever order eig returned them in:
G = hilb(5);
[Q, L] = eig(G);
lam = diag(L);
[~, idx] = sort(abs(lam), 'descend');
lam(idx).' % 1.5671 0.2085 0.0114 0.0003 0.0000
q1 = Q(:, idx(1)); q2 = Q(:, idx(2));
G2 = lam(idx(1))*(q1*q1') + lam(idx(2))*(q2*q2');
max(max(abs(G - G2))) % 0.0060: largest error in any entry
sum(abs(lam(idx(3:5)))) % 0.0117: the bound from the notesOctave prints the same five eigenvalues in e-notation. The last one is \(3.3\times 10^{-6}\), not zero: the Hilbert matrix is invertible, only badly conditioned (Chapter 3).
Unitary and similar matrices
The Fourier matrix (19) of the notes, built entry by entry. (fft(eye(n)) gives a Fourier matrix with \(e^{-2\pi j/n}\) in place of \(w\), the conjugate of this one.)
n = 4;
w = exp(2*pi*1j/n); % w = j for n = 4
F = w.^((0:n-1)' * (0:n-1)) % F(k+1,l+1) = w^(k*l), as in (19) of the notes
U = F/sqrt(n);
norm(U'*U - eye(n)) % of order 1e-16: U is unitary (' is the conjugate transpose)
abs(eig(U)).' % all 1: the eigenvalues lie on the unit circleA rotation through \(30^\circ\), and the conjugate in the inner product on \(\mathbb{C}^2\):
th = pi/6;
R = [cos(th) -sin(th); sin(th) cos(th)];
eig(R).' % 0.8660 + 0.5000i and 0.8660 - 0.5000i, i.e. exp(+-j*pi/6)
x = [1; -1j]; y = [1; 1j];
x'*y % 0: orthogonal in C^2
x.'*y % 2: the plain transpose misses itA matrix similar to the running example has the same eigenvalues. M\A*M computes \(M^{-1}AM\) without forming the inverse, and sort removes any dependence on the order in which eig returns them:
A = [4 -5; 2 -3]; M = [1 2; 3 4];
B = M\A*M % 15 16; -13 -14
sort(eig(A)).', sort(eig(B)).' % -1 2 for both