Let \( A \in \mathbb{R}^{n \times n} \) be a linear operator acting on vectors in \( \mathbb{R}^n \). A nonzero vector \( \mathbf{v} \in \mathbb{R}^n \) is called an eigenvector of \(A\) with associated eigenvalue \( \lambda \in \mathbb{C} \) if

\[ A \mathbf{v} = \lambda \mathbf{v}, \quad \mathbf{v} \neq \mathbf{0}. \]

The equation \(A\mathbf{v} = \lambda \mathbf{v}\) can be rewritten as \((A - \lambda I)\mathbf{v} = \mathbf{0}\). A nontrivial solution exists iff the matrix \(A - \lambda I\) is singular:

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

The polynomial \( p_A(\lambda) = \det(\lambda I - A) \) is the characteristic polynomial of \(A\). Its roots are the eigenvalues of \(A\) (counted with algebraic multiplicity). If \(A\) has \(n\) linearly independent eigenvectors, we can form the eigenvector matrix \( V = [\mathbf{v}_1 \ \cdots \ \mathbf{v}_n] \) and diagonal matrix \( \Lambda = \mathrm{diag}(\lambda_1,\dots,\lambda_n) \), and obtain the eigendecomposition

\[ A = V \Lambda V^{-1}, \quad AV = V\Lambda. \]

From the viewpoint of coordinate changes, the columns of \(V\) define a new basis in which the linear map \(A\) acts as simple coordinate-wise scaling by the eigenvalues.

flowchart TD
  A["Matrix A in R^n x R^n"] --> B["Solve (A - lambda I)v = 0"]
  B --> C["Eigenvalues lambda_i"]
  B --> D["Eigenvectors v_i"]
  C --> E["Scaling factors in modes"]
  D --> F["Invariant directions (modes)"]
  E --> G["Interpret behaviour of robot models"]
  F --> G
        

Consider a linear map \( A:\mathbb{R}^2 \rightarrow \mathbb{R}^2 \). It transforms lines through the origin to lines through the origin (because \(A\mathbf{0}=\mathbf{0}\)). An eigenvector direction \( \mathbf{v} \) is special because it is preserved by the transformation: all points on the line \( \{\alpha \mathbf{v} \mid \alpha \in \mathbb{R} \} \) are mapped back to the same line.

If \( \lambda > 1 \), the length of vectors along that direction is stretched; if \( 0 < \lambda < 1 \), vectors are contracted; and if \( \lambda < 0 \), the direction is flipped. In robotics, such privileged directions will correspond to principal axes of inertia, dominant directions of stiffness, or directions along which linearized motion grows or decays fastest.

For example, in \( \mathbb{R}^2 \), \[ A = \begin{bmatrix} 2 & 0 \\ 0 & 1 \end{bmatrix} \] has eigenvalues \( \lambda_1 = 2 \), \( \lambda_2 = 1 \) with eigenvectors \( \mathbf{e}_1 = [1,0]^\top \), \( \mathbf{e}_2 = [0,1]^\top \). The unit circle is mapped to an ellipse elongated by a factor 2 along the \(x\)-axis but unchanged along the \(y\)-axis.

In \( \mathbb{R}^3 \), three linearly independent eigenvectors of a diagonalizable matrix form a new orthogonal (or at least independent) coordinate system of principal directions. This picture is central for interpreting inertia matrices of links: eigenvectors give principal axes of rotation, eigenvalues give corresponding principal moments of inertia.

Many matrices in robot modeling are real and symmetric: inertia matrices, stiffness matrices, covariance matrices, and mass matrices of discretized structures. A fundamental result is the spectral theorem:

\[ \text{If } A = A^\top \in \mathbb{R}^{n \times n}, \text{ then} \quad \exists \text{ orthogonal } Q \in \mathbb{R}^{n \times n} \text{ and diagonal } \Lambda \text{ such that } A = Q\Lambda Q^\top. \]

Here \(Q^\top Q = QQ^\top = I\) and \(\Lambda = \mathrm{diag}(\lambda_1,\dots,\lambda_n)\) contains the real eigenvalues of \(A\). If, in addition, \(A\) is positive definite (all eigenvalues positive), it defines a quadratic form that can be interpreted as an energy:

\[ E(\mathbf{x}) = \tfrac{1}{2}\mathbf{x}^\top A \mathbf{x} = \tfrac{1}{2}\mathbf{z}^\top \Lambda \mathbf{z}, \quad \mathbf{z} = Q^\top \mathbf{x}. \]

The change of coordinates \( \mathbf{z} = Q^\top \mathbf{x} \) rotates the coordinate frame so that the energy splits into independent contributions: \( E(\mathbf{x}) = \tfrac{1}{2}\sum_i \lambda_i z_i^2 \). Each eigenvalue \(\lambda_i\) scales the contribution of coordinate \(z_i\).

• In a stiffness matrix \(K\), large \(\lambda_i\) correspond to directions that are stiff (large restoring forces), while small \(\lambda_i\) correspond to compliant directions.
• In an inertia matrix \(I\), eigenvectors are principal axes of rotation of a rigid body; eigenvalues are principal moments of inertia about those axes.
• In a covariance matrix \(\Sigma\), eigenvectors give principal directions of uncertainty; eigenvalues give variances along those directions.

flowchart TD
  X["State x (position, velocity, etc.)"] --> Q["Rotate: z = Q^T x via eigenvectors"]
  Q --> L["Quadratic form E = 0.5 * sum(lambda_i * z_i^2)"]
  L --> P1["Large lambda_i: stiff / high energy change"]
  L --> P2["Small lambda_i: compliant / low energy change"]
        

Eigenvalues and eigenvectors become especially important when we consider linear time-invariant (LTI) systems, which will appear systematically in later chapters. For now, consider the linear homogeneous system

\[ \dot{\mathbf{x}}(t) = A \mathbf{x}(t), \]

where \(A\) is diagonalizable with eigenpairs \((\lambda_i,\mathbf{v}_i)\). If we expand the initial condition as \( \mathbf{x}(0) = \sum_i c_i \mathbf{v}_i \), the solution is

\[ \mathbf{x}(t) = \sum_{i=1}^{n} c_i e^{\lambda_i t} \mathbf{v}_i. \]

Thus each eigenvalue \(\lambda_i\) determines the growth or decay rate of the corresponding mode \(\mathbf{v}_i\). If \( \mathrm{Re}(\lambda_i) < 0 \), the mode decays; if \( \mathrm{Re}(\lambda_i) > 0 \), it grows. In mechanical systems linearized around an equilibrium, these modes correspond to oscillatory or exponential behaviors along particular directions in state space.

In robotic manipulators, linearization of the dynamics around a stationary configuration yields a system where the eigenvalues/eigenvectors of the state matrix encode how perturbations in joint positions and velocities evolve. Although detailed derivations belong to later chapters, the key idea is that eigenvalues quantify rates and frequencies, while eigenvectors quantify shapes or directions of motion.

In principle, eigenvalues are roots of the characteristic polynomial \( p_A(\lambda) = \det(\lambda I - A) \). For \( n \leq 4 \), closed-form formulas exist; however, for higher dimensions direct polynomial root-finding is numerically unstable. Modern numerical linear algebra uses iterative schemes and matrix factorizations.

One simple algorithm illustrating the role of eigenvalues is the power iteration, which approximates the eigenvector associated with the eigenvalue of largest magnitude. Starting from a nonzero vector \(\mathbf{v}^{(0)}\):

\[ \mathbf{w}^{(k+1)} = A \mathbf{v}^{(k)}, \quad \mathbf{v}^{(k+1)} = \frac{\mathbf{w}^{(k+1)}}{\|\mathbf{w}^{(k+1)}\|_2}. \]

Under mild conditions (dominant eigenvalue magnitude strictly larger than the others and nonzero projection of \(\mathbf{v}^{(0)}\) onto the dominant eigenvector), the sequence \( \mathbf{v}^{(k)} \) converges to that eigenvector, and the Rayleigh quotient

\[ \rho(\mathbf{v}^{(k)}) = \frac{(\mathbf{v}^{(k)})^\top A \mathbf{v}^{(k)}}{(\mathbf{v}^{(k)})^\top \mathbf{v}^{(k)}} \]

converges to the corresponding eigenvalue. In robotics, such iterations appear in model reduction, stiffness analysis, and approximating dominant modes of motion.

In Python, the de facto standard libraries for eigenvalue computations are NumPy and SciPy, which wrap highly optimized LAPACK routines. In robotics-oriented toolboxes (e.g. roboticstoolbox-python), eigenvalue decompositions appear when analyzing inertia matrices, stiffness matrices, and covariance.

import numpy as np

# Example: 3x3 symmetric inertia-like matrix (kg*m^2)
I = np.array([[0.20, 0.01, 0.00],
              [0.01, 0.15, 0.00],
              [0.00, 0.00, 0.10]])

# Use eigh for symmetric (Hermitian) matrices: more stable than eig
eigvals, eigvecs = np.linalg.eigh(I)

print("Eigenvalues (principal moments):")
print(eigvals)
print("Eigenvectors (principal axes as columns):")
print(eigvecs)

# Verify eigen-equation A v = lambda v numerically
for i in range(len(eigvals)):
    lam = eigvals[i]
    v = eigvecs[:, i]
    lhs = I @ v
    rhs = lam * v
    print(f"Mode {i}: ||A v - lambda v|| = {np.linalg.norm(lhs - rhs):.3e}")
      

We can also implement a simple power iteration to approximate the dominant eigenpair of a (not necessarily symmetric) matrix:

import numpy as np

def power_iteration(A, max_iters=1000, tol=1e-8):
    n = A.shape[0]
    v = np.random.randn(n)
    v = v / np.linalg.norm(v)
    lam_old = 0.0

    for k in range(max_iters):
        w = A @ v
        v = w / np.linalg.norm(w)
        lam = float(v.T @ (A @ v))  # Rayleigh quotient

        if abs(lam - lam_old) < tol:
            break
        lam_old = lam

    return lam, v

# Example matrix: could be a local linearization of a 2D state update
A = np.array([[0.9, 0.1],
              [0.0, 0.8]])

lam_dom, v_dom = power_iteration(A)
print("Dominant eigenvalue approximation:", lam_dom)
print("Dominant eigenvector approximation:", v_dom)
      

In a robotics context, applying power_iteration to a stiffness matrix can reveal the stiffest mode of a joint or end-effector compliance model.

In C++ robotics software (e.g. ROS, Orocos KDL), the Eigen library is ubiquitous for linear algebra. It provides high-level wrappers for eigenvalue problems.

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

int main() {
    // Symmetric "inertia" matrix (for a link in its body frame)
    Eigen::Matrix3d I;
    I << 0.20, 0.01, 0.00,
          0.01, 0.15, 0.00,
          0.00, 0.00, 0.10;

    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(I);
    if (solver.info() != Eigen::Success) {
        std::cerr << "Eigen decomposition failed\n";
        return 1;
    }

    std::cout << "Eigenvalues (principal moments):\n"
              << solver.eigenvalues() << "\n\n";
    std::cout << "Eigenvectors (principal axes as columns):\n"
              << solver.eigenvectors() << "\n";

    // Check A v = lambda v for the first eigenpair
    Eigen::Vector3d v0 = solver.eigenvectors().col(0);
    double lambda0 = solver.eigenvalues()(0);
    Eigen::Vector3d lhs = I * v0;
    Eigen::Vector3d rhs = lambda0 * v0;
    std::cout << "Residual norm for first eigenpair: "
              << (lhs - rhs).norm() << "\n";
    return 0;
}
      

This pattern is typical in robot dynamics: once the inertia matrix is formed (often from CAD data or identification), one uses an eigen solver to obtain principal axes and moments, which are then used for simplification or analysis of rotational dynamics.

For Java-based robotics or simulation frameworks, typical choices include EJML (Efficient Java Matrix Library) or Apache Commons Math. The following example uses EJML's SimpleMatrix interface.

import org.ejml.simple.SimpleMatrix;
import org.ejml.data.Complex_F64;
import org.ejml.dense.row.decomposition.eig.EigenOps_DDRM;

public class EigenExample {
    public static void main(String[] args) {
        // 2x2 system matrix (e.g. linearization of a planar joint subsystem)
        double[][] dataA = {
                {0.9, 0.1},
                {0.0, 0.8}
        };

        SimpleMatrix A = new SimpleMatrix(dataA);

        // Compute eigen decomposition using EJML low-level API
        var eig = org.ejml.dense.row.decomposition.eig.SwitchingEigenDecomposition_DDRM
                .decompose(A.getDDRM(), true);

        int N = A.numRows();
        for (int i = 0; i < N; ++i) {
            Complex_F64 lambda = eig.getEigenvalue(i);
            System.out.println("Eigenvalue " + i + ": " + lambda.real + " + " + lambda.imaginary + "i");

            var v = eig.getEigenVector(i);
            if (v == null) {
                System.out.println("No eigenvector (defective eigenvalue)");
            } else {
                System.out.print("Eigenvector " + i + ": [");
                for (int j = 0; j < N; ++j) {
                    System.out.print(v.get(j) + (j + 1 < N ? ", " : ""));
                }
                System.out.println("]");
            }
        }
    }
}
      

In practice, Java-based robotic frameworks will wrap such matrix operations in higher-level abstractions, but the underlying mathematics remains the eigenvalue problem discussed here.

MATLAB is widely used in control engineering and robotics. The built-in function eig computes eigenvalues/eigenvectors, and the Robotics System Toolbox uses these operations internally in model analysis and linearization.

% Inertia-like symmetric matrix
I = [0.20 0.01 0.00;
     0.01 0.15 0.00;
     0.00 0.00 0.10];

% Eigen-decomposition
[V, D] = eig(I);   % I * V = V * D, columns of V are eigenvectors

lambda = diag(D);

disp('Eigenvalues (principal moments):');
disp(lambda);

disp('Eigenvectors (principal axes as columns):');
disp(V);

% Check orthogonality (spectral theorem)
disp('V''*V (should be identity for symmetric I):');
disp(V' * V);

% Example: transform angular velocity omega into principal coordinates
omega_body = [1; 2; 0.5];   % rad/s in body frame
omega_principal = V' * omega_body;
kinetic_energy = 0.5 * omega_principal' * D * omega_principal;
fprintf('Rotational kinetic energy: %f J\n', kinetic_energy);
      

In Simulink, linearization blocks can export an LTI model \( \dot{\mathbf{x}} = A\mathbf{x} + B\mathbf{u} \). Applying eig(A) in MATLAB reveals modal properties (eigenvalues/eigenvectors) of the linearized robot subsystem.

Wolfram Mathematica excels at symbolic and numeric eigenvalue computations. This can be particularly useful for deriving analytical expressions for simple manipulators or verifying symbolic derivations of inertia or stiffness matrices.

(* Define a 2x2 stiffness matrix symbolic in parameters k1, k2, alpha *)
K = { {k1 Cos[alpha]^2 + k2 Sin[alpha]^2,
      (k1 - k2) Sin[alpha] Cos[alpha]},
     {(k1 - k2) Sin[alpha] Cos[alpha],
      k1 Sin[alpha]^2 + k2 Cos[alpha]^2} };

(* Compute eigenvalues and eigenvectors symbolically *)
vals = Eigenvalues[K];
vecs = Eigenvectors[K];

vals
vecs

(* Substitute numeric values for a particular configuration *)
params = {k1 -> 1000, k2 -> 200, alpha -> Pi/6};
N[vals /. params]
N[vecs /. params]
      

This example models an effective planar stiffness rotated by an angle alpha. The eigenvalues are the principal stiffnesses, while the eigenvectors give principal directions. Such symbolic insight is valuable when designing compliant mechanisms or analyzing passive stiffness in robotic structures.

In this lesson we formalized eigenvalues and eigenvectors via the characteristic equation \( \det(\lambda I - A) = 0 \) and explored their geometric meaning as invariant directions and modal scaling factors of linear transformations. We emphasized the spectral theorem for real symmetric matrices, which guarantees orthogonal eigendecomposition and underlies the interpretation of inertia, stiffness, and covariance matrices in robotics. We previewed how eigenvalues of linearized dynamics govern stability and modal behavior, and we implemented eigenvalue computations in Python, C++, Java, MATLAB, and Mathematica, connecting them to physically meaningful quantities such as principal axes and energy ellipsoids. These ideas are foundational for later chapters on kinematics, dynamics, and model-based analysis.