Let \( A \in \mathbb{R}^{n\times n} \) represent a linear transformation (Lesson 2): \( T(\mathbf{x}) = A\mathbf{x} \). An eigenvector is a nonzero vector whose direction is preserved by \( A \), with only a scaling by an eigenvalue.

Eigenpair definition. A scalar \( \lambda \in \mathbb{C} \) is an eigenvalue of \( A \) if there exists \( \mathbf{v}\neq \mathbf{0} \) such that \( A\mathbf{v} = \lambda \mathbf{v} \). Then \( \mathbf{v} \) is an eigenvector associated with \( \lambda \).

In control engineering, eigenvalues will later appear as intrinsic “modes” of linear dynamics; diagonalization is the algebraic mechanism that can decouple those modes in a suitable coordinate basis (without yet invoking any state-space solution formulas).

flowchart TD
  A["Input: square matrix A"] --> B["Solve det(A - lambda I) = 0"]
  B --> C["Eigenvalues: lambda_1,...,lambda_m (with multiplicities)"]
  C --> D["For each lambda: solve (A - lambda I) v = 0"]
  D --> E["Eigenvectors / eigenspaces"]
  E --> F{"Do we have n linearly \nindependent eigenvectors?"}
  F -->|yes| G["Diagonalize: A = V D V^{-1}"]
  F -->|no| H["Not diagonalizable (needs more structure later)"]
  G --> I["Use diagonal form for powers, \nfunctions, and decoupled action"]
        

Throughout, we use the standard basis \( \{\mathbf{e}_1,\dots,\mathbf{e}_n\} \) and the identity matrix \( I \). We assume familiarity with determinants, rank/null space, and invertibility from Lesson 2.

The eigenvalue condition \( A\mathbf{v}=\lambda\mathbf{v} \) is equivalent to \( (A-\lambda I)\mathbf{v}=\mathbf{0} \) with \( \mathbf{v}\neq \mathbf{0} \). Thus, \( A-\lambda I \) must be singular.

Definition (characteristic polynomial). Define

\[ p_A(\lambda) \;=\; \det(A-\lambda I). \]

Theorem 2.1 (eigenvalues as roots). A scalar \( \lambda \) is an eigenvalue of \( A \) if and only if \( p_A(\lambda)=0 \).

Proof.

(⇒) If \( A\mathbf{v}=\lambda\mathbf{v} \) for some \( \mathbf{v}\neq \mathbf{0} \), then \( (A-\lambda I)\mathbf{v}=\mathbf{0} \), so \( A-\lambda I \) is singular and hence \( \det(A-\lambda I)=0 \).

(⇐) If \( \det(A-\lambda I)=0 \), then \( A-\lambda I \) is singular, so its null space is nontrivial: there exists \( \mathbf{v}\neq \mathbf{0} \) with \( (A-\lambda I)\mathbf{v}=\mathbf{0} \), i.e. \( A\mathbf{v}=\lambda\mathbf{v} \). □

Since \( p_A(\lambda) \) is a degree-\( n \) polynomial in \( \lambda \), the Fundamental Theorem of Algebra implies that \( A \) has \( n \) complex eigenvalues counting algebraic multiplicities.

Trace and determinant identities. If the eigenvalues (with multiplicity) are \( \lambda_1,\dots,\lambda_n \), then

\[ \operatorname{tr}(A) = \sum_{i=1}^n \lambda_i, \qquad \det(A) = \prod_{i=1}^n \lambda_i. \]

These follow by comparing coefficients of \( p_A(\lambda)=\det(A-\lambda I) \) to its factored form \( \prod_{i=1}^n (\lambda_i-\lambda) \) (up to a sign convention), a standard result from determinant expansions.

For a fixed eigenvalue \( \lambda \), the set of all eigenvectors (plus the zero vector) forms a subspace:

Definition (eigenspace). The eigenspace associated with \( \lambda \) is

\[ \mathcal{E}_{\lambda}(A) \;=\; \ker(A-\lambda I). \]

The geometric multiplicity of \( \lambda \) is \( \dim(\mathcal{E}_\lambda(A)) \). The algebraic multiplicity is the multiplicity of \( \lambda \) as a root of \( p_A(\lambda) \).

Theorem 3.1 (basic inequality). For any eigenvalue \( \lambda \),

\[ 1 \;\le\; \dim(\mathcal{E}_{\lambda}(A)) \;\le\; \text{algmult}(\lambda). \]

The left inequality holds because an eigenvalue implies a nontrivial null space. The right inequality is deeper (and will be revisited when discussing canonical forms later); it is consistent with all standard constructions in finite-dimensional linear algebra.

Theorem 3.2 (independence of eigenvectors for distinct eigenvalues). If \( \lambda_1,\dots,\lambda_k \) are distinct eigenvalues of \( A \) with corresponding eigenvectors \( \mathbf{v}_1,\dots,\mathbf{v}_k \) (each nonzero), then \( \{\mathbf{v}_1,\dots,\mathbf{v}_k\} \) is linearly independent.

Proof. Suppose a linear combination satisfies

\[ c_1\mathbf{v}_1 + \cdots + c_k\mathbf{v}_k = \mathbf{0}. \]

Apply \( A \) to both sides and use \( A\mathbf{v}_i=\lambda_i\mathbf{v}_i \):

\[ c_1\lambda_1\mathbf{v}_1 + \cdots + c_k\lambda_k\mathbf{v}_k = \mathbf{0}. \]

Subtract \( \lambda_k \) times the original equation:

\[ c_1(\lambda_1-\lambda_k)\mathbf{v}_1 + \cdots + c_{k-1}(\lambda_{k-1}-\lambda_k)\mathbf{v}_{k-1} = \mathbf{0}. \]

Because the eigenvalues are distinct, each factor \( \lambda_i-\lambda_k \neq 0 \) for \( i=1,\dots,k-1 \). Thus we obtain a dependence relation among \( \mathbf{v}_1,\dots,\mathbf{v}_{k-1} \). By repeating this elimination step (or using induction on \( k \)), we conclude \( c_1=\cdots=c_k=0 \). Therefore the vectors are independent. □

A key takeaway is that “many distinct eigenvalues” immediately implies “many independent eigenvectors,” which is precisely what makes diagonalization likely.

Diagonalization is a similarity representation that expresses a linear transformation in a basis of eigenvectors. It is the algebraic prototype of “decoupling” a linear action into independent scalar scalings.

Definition (diagonalizable matrix). A matrix \( A\in\mathbb{C}^{n\times n} \) is diagonalizable if there exist an invertible \( V \) and a diagonal \( D \) such that

\[ A = V D V^{-1}, \qquad D = \operatorname{diag}(\lambda_1,\dots,\lambda_n). \]

If \( V=[\mathbf{v}_1\,\cdots\,\mathbf{v}_n] \) collects eigenvectors and \( D \) collects the associated eigenvalues, then the identity \( AV=VD \) is equivalent to \( A\mathbf{v}_i = \lambda_i\mathbf{v}_i \) for each column.

Theorem 4.1 (diagonalization criterion). A matrix \( A \) is diagonalizable if and only if \( A \) has \( n \) linearly independent eigenvectors.

Proof.

(⇒) If \( A=VDV^{-1} \) with \( V \) invertible and \( D \) diagonal, write \( V=[\mathbf{v}_1\cdots\mathbf{v}_n] \). Then \( AV=VD \), so the \( i \)-th column gives \( A\mathbf{v}_i=\lambda_i\mathbf{v}_i \). Since \( V \) is invertible, its columns \( \mathbf{v}_i \) are linearly independent, hence we have \( n \) independent eigenvectors.

(⇐) If \( \mathbf{v}_1,\dots,\mathbf{v}_n \) are linearly independent eigenvectors with \( A\mathbf{v}_i=\lambda_i\mathbf{v}_i \), define \( V=[\mathbf{v}_1\cdots\mathbf{v}_n] \) and \( D=\operatorname{diag}(\lambda_1,\dots,\lambda_n) \). Then \( AV=VD \). Multiply by \( V^{-1} \) on the right to obtain \( A=VDV^{-1} \). □

Corollary 4.2 (distinct eigenvalues imply diagonalizable). If \( A \) has \( n \) distinct eigenvalues, then \( A \) is diagonalizable.

Proof. With distinct eigenvalues, Theorem 3.2 gives \( n \) linearly independent eigenvectors, and Theorem 4.1 applies. □

Algebraic payoff. If \( A=VDV^{-1} \), then for any integer \( k\ge 0 \),

\[ A^k = V D^k V^{-1}, \qquad D^k = \operatorname{diag}(\lambda_1^k,\dots,\lambda_n^k). \]

Proof is by associativity: \( A^2=(VDV^{-1})(VDV^{-1})=VD^2V^{-1} \), and induction completes the argument. This identity later becomes essential whenever matrix powers (or matrix functions) appear.

flowchart TD
  X["Vector x (standard coordinates)"] -->|"Change basis: z = inv(V) * x"| Z["Coordinates z (eigenvector basis)"]
  Z -->|"Diagonal action: w = D * z"| W["Scaled coordinates w"]
  W -->|"Return: y = V * w"| Y["Output y = A * x"]
  N["Interpretation: in eigenvector basis, \nA acts as independent scalings by lambda_i"] --- Z
  

A particularly important (and common) case is \( A=A^\top \). Symmetry yields stronger structure: real eigenvalues and orthogonal eigenvectors, enabling numerically stable diagonalization via orthonormal bases.

Lemma 5.1 (orthogonality for distinct eigenvalues). If \( A=A^\top \) and \( A\mathbf{v}_i=\lambda_i\mathbf{v}_i \), \( A\mathbf{v}_j=\lambda_j\mathbf{v}_j \) with \( \lambda_i\neq \lambda_j \), then \( \mathbf{v}_i^\top \mathbf{v}_j = 0 \).

Proof. Compute two ways:

\[ \mathbf{v}_i^\top A \mathbf{v}_j = \mathbf{v}_i^\top (\lambda_j \mathbf{v}_j) = \lambda_j \mathbf{v}_i^\top \mathbf{v}_j, \]

and using symmetry \( A^\top=A \):

\[ \mathbf{v}_i^\top A \mathbf{v}_j = (A^\top \mathbf{v}_i)^\top \mathbf{v}_j = (A\mathbf{v}_i)^\top \mathbf{v}_j = (\lambda_i \mathbf{v}_i)^\top \mathbf{v}_j = \lambda_i \mathbf{v}_i^\top \mathbf{v}_j. \]

Hence \( \lambda_j \mathbf{v}_i^\top \mathbf{v}_j = \lambda_i \mathbf{v}_i^\top \mathbf{v}_j \), so \( (\lambda_j-\lambda_i)\mathbf{v}_i^\top \mathbf{v}_j=0 \). Since \( \lambda_i\neq \lambda_j \), we conclude \( \mathbf{v}_i^\top \mathbf{v}_j=0 \). □

Theorem 5.2 (Spectral Theorem, real symmetric form). If \( A\in\mathbb{R}^{n\times n} \) and \( A=A^\top \), then there exists an orthogonal matrix \( Q \) (so \( Q^\top Q=I \)) and a real diagonal matrix \( \Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_n) \) such that

\[ A = Q \Lambda Q^\top. \]

Proof sketch (university-level, constructive idea).

Consider the Rayleigh quotient \( \rho(\mathbf{x})=\frac{\mathbf{x}^\top A \mathbf{x}}{\mathbf{x}^\top \mathbf{x}} \) for \( \mathbf{x}\neq \mathbf{0} \). On the unit sphere \( \{\mathbf{x}:\|\mathbf{x}\|_2=1\} \), \( \rho(\mathbf{x})=\mathbf{x}^\top A\mathbf{x} \) is continuous and attains a maximum at some \( \mathbf{q}_1 \) (compactness). Using Lagrange multipliers for the constrained maximization, one obtains \( A\mathbf{q}_1=\lambda_1 \mathbf{q}_1 \) with \( \lambda_1=\rho(\mathbf{q}_1)\in\mathbb{R} \).

Next, restrict \( A \) to the orthogonal complement of \( \mathbf{q}_1 \). By Lemma 5.1, this subspace is invariant under \( A \) in the sense needed to continue the construction. Apply the same argument inductively to obtain an orthonormal set \( \mathbf{q}_1,\dots,\mathbf{q}_n \) of eigenvectors. Let \( Q=[\mathbf{q}_1\cdots\mathbf{q}_n] \) and \( \Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_n) \). Then \( AQ=Q\Lambda \) and, since \( Q \) is orthogonal, \( A=Q\Lambda Q^\top \). □

The symmetric case is central in optimization and energy methods (later), because quadratic forms \( \mathbf{x}^\top A\mathbf{x} \) are diagonalized by orthogonal eigenbases without ill-conditioned inverses.

Beyond closed-form solutions for small matrices, eigenvalues are typically computed numerically. A conceptually simple starting point is power iteration, which approximates the dominant eigenpair under mild conditions.

Assumption (dominant eigenvalue). Suppose \( A \) has an eigenvalue \( \lambda_1 \) such that \( |\lambda_1| > |\lambda_2| \ge \cdots \ge |\lambda_n| \) (strict spectral gap), and the initial vector \( \mathbf{x}_0 \) has a nonzero component along the eigenvector associated with \( \lambda_1 \).

Algorithm (power iteration).

\[ \mathbf{x}_{k+1} \;=\; \frac{A\mathbf{x}_k}{\|A\mathbf{x}_k\|_2}, \qquad \hat{\lambda}_{k} \;=\; \mathbf{x}_k^\top A \mathbf{x}_k \quad (\text{Rayleigh estimate, symmetric case}). \]

Convergence statement (informal but precise). Under the spectral gap assumption,

\[ \mathbf{x}_k \;\rightarrow\; \pm \mathbf{v}_1 \quad\text{and}\quad \hat{\lambda}_k \;\rightarrow\; \lambda_1, \quad\text{with an error rate proportional to}\quad \left|\frac{\lambda_2}{\lambda_1}\right|^k. \]

The proof follows by expressing \( \mathbf{x}_0 \) in an eigenvector basis (when diagonalizable), observing \( A^k\mathbf{x}_0 \) scales components by \( \lambda_i^k \), and normalizing. We avoid deeper canonical-form details here; the diagonalizable case already makes the mechanism transparent via \( A^k = V D^k V^{-1} \) from Section 4.

Libraries commonly used in control workflows include numpy and scipy.linalg. Here we: (i) compute eigenvalues/eigenvectors, (ii) attempt diagonalization, and (iii) implement power iteration from scratch.

import numpy as np
from numpy.linalg import norm

# Example matrix (choose one with real, distinct eigenvalues for clean diagonalization)
A = np.array([[4.0, 1.0, 0.0],
              [1.0, 3.0, 1.0],
              [0.0, 1.0, 2.0]])

# 1) Eigen-decomposition
w, V = np.linalg.eig(A)  # w: eigenvalues, V: eigenvectors (columns)
print("Eigenvalues:", w)
print("Eigenvectors (columns):\n", V)

# 2) Diagonalization check: A ?= V D V^{-1}
D = np.diag(w)
A_recon = V @ D @ np.linalg.inv(V)
print("Reconstruction error ||A - VDV^{-1}||_F =", norm(A - A_recon, ord='fro'))

# 3) Symmetric case: orthogonal diagonalization using eigh (more stable for symmetric A)
ws, Q = np.linalg.eigh(A)  # for symmetric matrices: real eigenvalues, orthonormal eigenvectors
Lambda = np.diag(ws)
A_recon2 = Q @ Lambda @ Q.T
print("Symmetric recon error ||A - QΛQ^T||_F =", norm(A - A_recon2, ord='fro'))

# 4) Power iteration (from scratch) for dominant eigenpair
def power_iteration(A, x0=None, max_iter=2000, tol=1e-10):
    n = A.shape[0]
    if x0 is None:
        x = np.random.randn(n)
    else:
        x = np.array(x0, dtype=float).reshape(n)
    x = x / norm(x)

    lam_old = 0.0
    for k in range(max_iter):
        y = A @ x
        x = y / norm(y)

        # Rayleigh quotient estimate (works best when A is symmetric)
        lam = float(x.T @ (A @ x))

        if abs(lam - lam_old) < tol:
            return lam, x, k+1
        lam_old = lam
    return lam, x, max_iter

lam_hat, v_hat, iters = power_iteration(A)
print("Power iteration: lambda_hat =", lam_hat, "iters =", iters)
print("v_hat =", v_hat)

Practical note: for real symmetric matrices, prefer eigh (or scipy.linalg.eigh) because it exploits symmetry for improved numerical stability and guarantees real outputs when \( A=A^\top \).

A widely used C++ linear algebra library is Eigen. The snippet below computes eigenvalues/eigenvectors for a symmetric matrix, then implements power iteration from scratch.

#include <iostream>
#include <Eigen/Dense>

int main() {
  Eigen::MatrixXd A(3,3);
  A << 4.0, 1.0, 0.0,
       1.0, 3.0, 1.0,
       0.0, 1.0, 2.0;

  // 1) Symmetric eigendecomposition (spectral theorem scenario)
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(A);
  if (es.info() != Eigen::Success) {
    std::cerr << "Eigen decomposition failed." << std::endl;
    return 1;
  }

  std::cout << "Eigenvalues:\n" << es.eigenvalues() << "\n";
  std::cout << "Eigenvectors (columns):\n" << es.eigenvectors() << "\n";

  // 2) Reconstruct A = Q Λ Q^T
  Eigen::MatrixXd Q = es.eigenvectors();
  Eigen::MatrixXd Lambda = es.eigenvalues().asDiagonal();
  Eigen::MatrixXd Arec = Q * Lambda * Q.transpose();
  std::cout << "Reconstruction error (Frobenius): "
            << (A - Arec).norm() << "\n";

  // 3) Power iteration from scratch
  Eigen::VectorXd x = Eigen::VectorXd::Random(3);
  x.normalize();

  double lam_old = 0.0;
  const int max_iter = 2000;
  const double tol = 1e-10;

  for (int k = 0; k < max_iter; ++k) {
    Eigen::VectorXd y = A * x;
    x = y.normalized();

    double lam = x.transpose() * A * x;  // Rayleigh estimate (best for symmetric A)
    if (std::abs(lam - lam_old) < tol) {
      std::cout << "Power iteration converged in " << (k+1) << " iterations.\n";
      std::cout << "lambda_hat = " << lam << "\n";
      std::cout << "v_hat^T = " << x.transpose() << "\n";
      break;
    }
    lam_old = lam;
  }

  return 0;
}

Control-related libraries in C++ ecosystems often build on such linear algebra kernels; Eigen is a typical dependency for state-space tooling, simulation, estimation, and optimization stacks.

In Java, a practical numerical linear algebra option is EJML. The following shows: (i) eigen-decomposition for a real matrix (general case), and (ii) a from-scratch power iteration.

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;
import org.ejml.dense.row.NormOps_DDRM;
import org.ejml.dense.row.factory.DecompositionFactory_DDRM;
import org.ejml.interfaces.decomposition.EigenDecomposition_F64;

public class EigenLesson {
  public static void main(String[] args) {
    // Example matrix
    DMatrixRMaj A = new DMatrixRMaj(new double[][]{
      {4.0, 1.0, 0.0},
      {1.0, 3.0, 1.0},
      {0.0, 1.0, 2.0}
    });

    // 1) Eigen-decomposition (general; eigenvectors may be complex in general problems)
    EigenDecomposition_F64<DMatrixRMaj> eig =
        DecompositionFactory_DDRM.eig(A.numRows, true);

    if (!eig.decompose(A)) {
      System.out.println("Eigendecomposition failed.");
      return;
    }

    System.out.println("Eigenvalues:");
    for (int i = 0; i < A.numRows; i++) {
      System.out.println("  " + eig.getEigenvalue(i));
    }

    // 2) Power iteration (from scratch) for dominant eigenpair (works well with spectral gap)
    DMatrixRMaj x = new DMatrixRMaj(A.numRows, 1);
    for (int i = 0; i < A.numRows; i++) x.set(i, 0, Math.random());
    normalize2(x);

    double lamOld = 0.0;
    int maxIter = 2000;
    double tol = 1e-10;

    DMatrixRMaj y = new DMatrixRMaj(A.numRows, 1);
    for (int k = 0; k < maxIter; k++) {
      CommonOps_DDRM.mult(A, x, y);
      normalize2(y);

      // Rayleigh quotient lambda = x^T A x (reuse y = A x if desired)
      CommonOps_DDRM.mult(A, y, x);  // x <- A*y (temporary)
      double lam = dot(y, x);        // y^T (A y)
      // restore x to current normalized eigenvector estimate:
      x.set(y);

      if (Math.abs(lam - lamOld) < tol) {
        System.out.println("Power iteration converged in " + (k+1) + " iterations.");
        System.out.println("lambda_hat = " + lam);
        System.out.println("v_hat (approx) = ");
        x.print();
        break;
      }
      lamOld = lam;
    }
  }

  static double dot(DMatrixRMaj a, DMatrixRMaj b) {
    double s = 0.0;
    for (int i = 0; i < a.getNumElements(); i++) s += a.get(i) * b.get(i);
    return s;
  }

  static void normalize2(DMatrixRMaj v) {
    double n2 = NormOps_DDRM.normF(v);
    if (n2 > 0) CommonOps_DDRM.scale(1.0 / n2, v);
  }
}

For control-focused Java stacks, EJML (or Apache Commons Math) typically provides the eigenvalue routines used for stability checks, model transformations, and numerical conditioning diagnostics.

MATLAB provides robust eigenvalue solvers via eig (general) and specialized routines for symmetric matrices. Below we compute \( A=VDV^{-1} \) and \( A=Q\Lambda Q^\top \) (when symmetric), then implement power iteration from scratch.

% Example matrix
A = [4 1 0;
     1 3 1;
     0 1 2];

% 1) General eigen-decomposition
[V,D] = eig(A);
disp('Eigenvalues (diag(D)) = '); disp(diag(D));
disp('Check ||A - VDV^{-1}||_F = '); disp(norm(A - V*D/V, 'fro'));

% 2) Symmetric case: orthogonal diagonalization (A is symmetric here)
[Q,Lambda] = eig((A + A')/2);
disp('Check ||A - QΛQ^T||_F = '); disp(norm(A - Q*Lambda*Q', 'fro'));

% 3) Power iteration from scratch
x = randn(size(A,1),1);
x = x / norm(x);

lam_old = 0;
tol = 1e-10;
maxIter = 2000;

for k = 1:maxIter
    y = A*x;
    x = y / norm(y);
    lam = x' * A * x; % Rayleigh estimate

    if abs(lam - lam_old) < tol
        fprintf('Converged in %d iterations\\n', k);
        fprintf('lambda_hat = %.12f\\n', lam);
        disp('v_hat = '); disp(x);
        break;
    end
    lam_old = lam;
end

Simulink note (minimal workflow consistent with this chapter).

Without requiring state-space modeling yet, you can still integrate eigen-computation into a model: use a Constant block to output a matrix \( A \), feed it into a MATLAB Function block that calls eig, and route the eigenvalues to a Display or To Workspace block.

% MATLAB Function block code (input: A, output: lam)
function lam = fcn(A)
%#codegen
lam = eig(A);
end

Later chapters will use eigenvalues as diagnostic quantities for dynamic matrices (e.g., checking intrinsic stability properties), but the computational core is the same: eig or structured solvers when symmetry or other structure is present.

Mathematica provides symbolic and numeric eigen-analysis. The following demonstrates eigenpairs, diagonalization, and reconstruction checks.

A = { {4, 1, 0},
     {1, 3, 1},
     {0, 1, 2} };

(* 1) Eigenvalues and eigenvectors *)
{vals, vecs} = Eigensystem[A];
vals
vecs

(* 2) Construct V and D and verify A == V D V^{-1} *)
V = Transpose[vecs];              (* columns are eigenvectors *)
D = DiagonalMatrix[vals];
Arec = V . D . Inverse[V];
Norm[A - Arec, "Frobenius"]

(* 3) For symmetric A: orthonormal eigenvectors can be enforced *)
vecsOrtho = Orthogonalize[vecs];
Q = Transpose[vecsOrtho];
Lambda = DiagonalMatrix[Diagonal[Transpose[Q].A.Q]];
Norm[A - Q . Lambda . Transpose[Q], "Frobenius"]

(* 4) Power iteration from scratch *)
powerIteration[A_, x0_, maxIter_:2000, tol_:10^-10] := Module[
  {x = x0/Norm[x0], lamOld = 0, lam = 0, y, k},
  For[k = 1, k <= maxIter, k++,
    y = A.x;
    x = y/Norm[y];
    lam = x.(A.x); (* Rayleigh estimate *)
    If[Abs[lam - lamOld] < tol, Return[{lam, x, k}]];
    lamOld = lam;
  ];
  {lam, x, maxIter}
];

{lamHat, vHat, iters} = powerIteration[A, RandomReal[NormalDistribution[0,1], 3]];
lamHat
vHat
iters

Problem 1 (Eigenvalue condition via singularity). Prove that \( \lambda \) is an eigenvalue of \( A \) if and only if \( A-\lambda I \) is singular.

Solution. If \( A\mathbf{v}=\lambda\mathbf{v} \) with \( \mathbf{v}\neq \mathbf{0} \), then \( (A-\lambda I)\mathbf{v}=\mathbf{0} \), so the null space is nontrivial and the matrix is singular. Conversely, if \( A-\lambda I \) is singular, there exists \( \mathbf{v}\neq \mathbf{0} \) with \( (A-\lambda I)\mathbf{v}=\mathbf{0} \), hence \( A\mathbf{v}=\lambda\mathbf{v} \). □

Problem 2 (Independence of distinct-eigenvalue eigenvectors). Let \( A\in\mathbb{C}^{n\times n} \) and suppose \( A\mathbf{v}_i=\lambda_i\mathbf{v}_i \) with \( \lambda_i\neq \lambda_j \) for \( i\neq j \). Prove \( \mathbf{v}_1,\dots,\mathbf{v}_k \) are linearly independent.

Solution. Use the elimination argument in Theorem 3.2: assume a dependence relation, apply \( A \), subtract \( \lambda_k \) times the original relation, and induct to conclude all coefficients are zero. □

Problem 3 (Diagonalization criterion). Show that \( A \) is diagonalizable if and only if \( \mathbb{C}^n \) has a basis of eigenvectors of \( A \).

Solution. If \( A=VDV^{-1} \), columns of \( V \) are eigenvectors forming a basis. Conversely, if eigenvectors form a basis, assemble them into \( V \), define \( D \) with corresponding eigenvalues, and obtain \( AV=VD \Rightarrow A=VDV^{-1} \). □

Problem 4 (Orthogonality in the symmetric case). Let \( A=A^\top \) and \( A\mathbf{v}_1=\lambda_1\mathbf{v}_1 \), \( A\mathbf{v}_2=\lambda_2\mathbf{v}_2 \) with \( \lambda_1\neq \lambda_2 \). Prove \( \mathbf{v}_1^\top\mathbf{v}_2=0 \).

Solution. Compare \( \mathbf{v}_1^\top A\mathbf{v}_2 \) computed as \( \lambda_2\mathbf{v}_1^\top\mathbf{v}_2 \) and via symmetry as \( (A\mathbf{v}_1)^\top\mathbf{v}_2=\lambda_1\mathbf{v}_1^\top\mathbf{v}_2 \). Subtract to get \( (\lambda_2-\lambda_1)\mathbf{v}_1^\top\mathbf{v}_2=0 \), hence orthogonality. □

Problem 5 (Power iteration rate in the diagonalizable case). Assume \( A=VDV^{-1} \) with \( D=\operatorname{diag}(\lambda_1,\dots,\lambda_n) \) and \( |\lambda_1| > |\lambda_2| \ge \cdots \). Let \( \mathbf{x}_0 \) satisfy \( \mathbf{x}_0 = V\mathbf{c} \) with \( c_1\neq 0 \). Show that the normalized iterates of power iteration converge in direction to the eigenvector of \( \lambda_1 \) at a rate proportional to \( \left|\frac{\lambda_2}{\lambda_1}\right|^k \).

Solution. Since \( \mathbf{x}_0=V\mathbf{c} \),

\[ A^k\mathbf{x}_0 = VD^kV^{-1}V\mathbf{c} = V D^k \mathbf{c} = c_1\lambda_1^k\mathbf{v}_1 + \sum_{i=2}^n c_i\lambda_i^k\mathbf{v}_i. \]

Factor out \( \lambda_1^k \):

\[ A^k\mathbf{x}_0 = \lambda_1^k\left(c_1\mathbf{v}_1 + \sum_{i=2}^n c_i\left(\frac{\lambda_i}{\lambda_1}\right)^k\mathbf{v}_i\right). \]

Because \( \left|\frac{\lambda_i}{\lambda_1}\right| \le \left|\frac{\lambda_2}{\lambda_1}\right| < 1 \) for \( i\ge 2 \), the non-dominant terms decay geometrically. Normalization removes the scalar \( \lambda_1^k \), leaving direction convergence to \( \mathbf{v}_1 \) with rate controlled by \( \left|\frac{\lambda_2}{\lambda_1}\right|^k \). □

We defined eigenvalues/eigenvectors via the singularity of \( A-\lambda I \), established core structural results (eigenspaces, multiplicities, independence for distinct eigenvalues), and proved the diagonalization criterion \( A=VDV^{-1} \) precisely when \( n \) independent eigenvectors exist. For symmetric matrices, we proved orthogonality and outlined the spectral theorem \( A=Q\Lambda Q^\top \). Finally, we implemented eigen-analysis across Python, C++, Java, MATLAB/Simulink, and Mathematica, including a from-scratch power iteration.