Linear Systems & Matrix Exponential

A linear system of first-order ordinary differential equations with constant coefficients is an equation

$\mathbf{y}'(t) = A\mathbf{y}(t),$

where $A \in \mathbb{R}^{n \times n}$ is a constant matrix and $\mathbf{y}: I \to \mathbb{R}^n$ is the unknown vector-valued function. The initial value problem is $\mathbf{y}'(t) = A\mathbf{y}(t)$, $\mathbf{y}(0) = \mathbf{y}_0$. The system is homogeneous if there is no additional forcing term; the non-homogeneous version is $\mathbf{y}'(t) = A\mathbf{y}(t) + \mathbf{g}(t)$.

Two masses connected by springs: $y_1' = y_2$, $y_2' = -\omega^2 y_1$. In matrix form:

$\mathbf{y}' = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix} \mathbf{y}.$

The eigenvalues are $\pm i\omega$, producing oscillatory solutions $y_1(t) = A\cos(\omega t + \phi)$ — the system rotates in the phase plane. This is the harmonic oscillator, the most fundamental dynamical system in physics and engineering.

The second-order equation $y'' + 3y' + 2y = 0$ converts to a first-order system via $y_1 = y$, $y_2 = y'$:

$\mathbf{y}' = \begin{pmatrix} 0 & 1 \\ -2 & -3 \end{pmatrix} \mathbf{y}.$

Every $n$-th order linear ODE becomes an $n \times n$ first-order system. This reduction is not just a notational convenience — it places all linear ODEs (regardless of order) into a single framework where eigenvalue methods apply uniformly.

When $A$ is constant, the system is autonomous: the velocity field $A\mathbf{y}$ depends only on the state $\mathbf{y}$, not on time. This is the typical case in ML — the loss landscape does not change during training (for fixed data). Time-dependent systems $\mathbf{y}' = A(t)\mathbf{y}$ are harder and appear in §7 via the non-homogeneous formulation.

For any $n \times n$ matrix $A$ and scalar $t \in \mathbb{R}$, the matrix exponential is

$e^{At} = \sum_{k=0}^{\infty} \frac{(At)^k}{k!},$

where $A^0 = I$ (the identity matrix) and $A^k$ denotes the $k$-th matrix power. The series converges absolutely in any matrix norm.

For any $A \in \mathbb{R}^{n \times n}$ and $t \in \mathbb{R}$, the series $\sum_{k=0}^{\infty} \frac{(At)^k}{k!}$ converges absolutely, and $\|e^{At}\| \leq e^{\|A\||t|}$.

For any submultiplicative norm $\|\cdot\|$, we have

$\left\|\frac{(At)^k}{k!}\right\| \leq \frac{\|A\|^k |t|^k}{k!}.$

The scalar series $\sum_{k=0}^{\infty} \frac{(\|A\||t|)^k}{k!} = e^{\|A\||t|}$ converges (it is the exponential function evaluated at $\|A\||t|$). Therefore the matrix series converges absolutely by comparison.

Absolute convergence in a finite-dimensional normed space implies convergence (because $\mathbb{R}^{n \times n}$ with any matrix norm is a complete normed space), so $e^{At}$ is well-defined.

The bound follows from the triangle inequality applied to the partial sums:

$\|e^{At}\| = \left\|\sum_{k=0}^{\infty} \frac{(At)^k}{k!}\right\| \leq \sum_{k=0}^{\infty} \frac{\|A\|^k |t|^k}{k!} = e^{\|A\||t|}. \qquad \blacksquare$

Let $A, B \in \mathbb{R}^{n \times n}$ and $t, s \in \mathbb{R}$. Then:

(a) $e^{A \cdot 0} = I$.

(b) $\displaystyle\frac{d}{dt}e^{At} = Ae^{At} = e^{At}A$.

(c) $e^{A(t+s)} = e^{At}e^{As}$ for all $t, s$.

(d) $e^{At}$ is invertible with $(e^{At})^{-1} = e^{-At}$.

(e) If $AB = BA$, then $e^{(A+B)t} = e^{At}e^{Bt}$.

(f) If $AB \neq BA$, then $e^{(A+B)t} \neq e^{At}e^{Bt}$ in general.

Property (e) fails without commutativity. The Baker-Campbell-Hausdorff formula gives the correction: $e^A e^B = e^{A + B + \frac{1}{2}[A,B] + \cdots}$ where $[A, B] = AB - BA$ is the commutator. This non-commutativity is why matrix exponentials are harder than scalar exponentials, and why Lie theory is the natural language for matrix groups. In ML, this matters when composing transformations that don’t commute — rotation followed by scaling is not the same as scaling followed by rotation.

If $A = \text{diag}(\lambda_1, \ldots, \lambda_n)$, then $e^{At} = \text{diag}(e^{\lambda_1 t}, \ldots, e^{\lambda_n t})$ — each component evolves independently.

If $N$ is nilpotent ($N^m = 0$ for some $m$), then $e^{Nt} = I + Nt + \frac{(Nt)^2}{2!} + \cdots + \frac{(Nt)^{m-1}}{(m-1)!}$ — a matrix polynomial in $t$, since all higher powers vanish. This finite truncation makes nilpotent exponentials easy to compute.

Let $A \in \mathbb{R}^{n \times n}$ be diagonalizable with eigenvalues $\lambda_1, \ldots, \lambda_n$ and corresponding linearly independent eigenvectors $\mathbf{v}_1, \ldots, \mathbf{v}_n$. The general solution of $\mathbf{y}' = A\mathbf{y}$ is

$\mathbf{y}(t) = c_1 e^{\lambda_1 t}\mathbf{v}_1 + \cdots + c_n e^{\lambda_n t}\mathbf{v}_n,$

where $c_i \in \mathbb{R}$ (or $\mathbb{C}$) are constants determined by $\mathbf{y}(0) = \mathbf{y}_0$.

Each $e^{\lambda_i t}\mathbf{v}_i$ satisfies

$\frac{d}{dt}\bigl[e^{\lambda_i t}\mathbf{v}_i\bigr] = \lambda_i e^{\lambda_i t}\mathbf{v}_i = A\bigl(e^{\lambda_i t}\mathbf{v}_i\bigr)$

because $A\mathbf{v}_i = \lambda_i \mathbf{v}_i$. By linearity of differentiation, any linear combination $\sum c_i e^{\lambda_i t}\mathbf{v}_i$ is also a solution.

Since $\{\mathbf{v}_1, \ldots, \mathbf{v}_n\}$ is a basis for $\mathbb{R}^n$, every initial condition $\mathbf{y}_0 = \sum c_i \mathbf{v}_i$ can be uniquely represented, so the general solution captures all solutions. $\blacksquare$

Consider $A = \begin{pmatrix} -1 & 2 \\ 0 & -3 \end{pmatrix}$. The eigenvalues are $\lambda_1 = -1$, $\lambda_2 = -3$, with eigenvectors $\mathbf{v}_1 = (1, 0)^T$, $\mathbf{v}_2 = (1, -1)^T$. The general solution is

$\mathbf{y}(t) = c_1 e^{-t}\begin{pmatrix} 1 \\ 0 \end{pmatrix} + c_2 e^{-3t}\begin{pmatrix} 1 \\ -1 \end{pmatrix}.$

Both eigenvalues are negative, so all solutions decay to the origin — a stable node. The fast component ($e^{-3t}$) dies out first, and trajectories ultimately align with the slow eigendirection $\mathbf{v}_1$.

The rotation matrix $A = \begin{pmatrix} 0 & -\omega \\ \omega & 0 \end{pmatrix}$ has eigenvalues $\pm i\omega$. The solution is $\mathbf{y}(t) = c_1 \cos(\omega t) + c_2 \sin(\omega t)$ — pure rotation in the phase plane (a center). With damping, $A = \begin{pmatrix} -\alpha & -\omega \\ \omega & -\alpha \end{pmatrix}$ ($\alpha > 0$), the eigenvalues become $-\alpha \pm i\omega$, and solutions spiral inward — a stable spiral.

When $A$ is real and $\lambda = \alpha + i\beta$ is a complex eigenvalue with eigenvector $\mathbf{v} = \mathbf{a} + i\mathbf{b}$, the two real-valued solutions are

$e^{\alpha t}\bigl(\cos(\beta t)\,\mathbf{a} - \sin(\beta t)\,\mathbf{b}\bigr) \quad\text{and}\quad e^{\alpha t}\bigl(\sin(\beta t)\,\mathbf{a} + \cos(\beta t)\,\mathbf{b}\bigr).$

The factor $e^{\alpha t}$ controls growth or decay; $\cos(\beta t)$ and $\sin(\beta t)$ produce rotation. When $\alpha < 0$, we get spirals that decay inward; when $\alpha > 0$, spirals that grow outward; when $\alpha = 0$, closed ellipses.

The phase portrait of $\mathbf{y}' = A\mathbf{y}$ is the collection of all solution trajectories $\mathbf{y}(t)$ plotted in the phase plane $(y_1, y_2)$. The origin $\mathbf{0}$ is the only equilibrium point (assuming $A$ is invertible, i.e., $\det(A) \neq 0$).

Let $A$ be a $2 \times 2$ real matrix with $\tau = \text{tr}(A)$ and $\Delta = \det(A)$, and define the discriminant $D = \tau^2 - 4\Delta$. The phase portrait of $\mathbf{y}' = A\mathbf{y}$ is classified as follows:

(i) $\Delta < 0$: saddle (eigenvalues of opposite sign).

(ii) $\Delta > 0$, $D > 0$, $\tau < 0$: stable node (distinct negative real eigenvalues).

(iii) $\Delta > 0$, $D > 0$, $\tau > 0$: unstable node (distinct positive real eigenvalues).

(iv) $\Delta > 0$, $D < 0$, $\tau < 0$: stable spiral (complex eigenvalues with negative real part).

(v) $\Delta > 0$, $D < 0$, $\tau > 0$: unstable spiral (complex eigenvalues with positive real part).

(vi) $\Delta > 0$, $\tau = 0$: center (pure imaginary eigenvalues).

(vii) $\Delta = 0$: degenerate (at least one zero eigenvalue; line of equilibria).

Place a dot at $(\tau, \Delta)$ for any $2 \times 2$ matrix, and the diagram immediately tells you the phase portrait type. The parabola $\Delta = \tau^2/4$ separates nodes from spirals. The $\tau$-axis ($\tau = 0$) separates stable from unstable. The $\Delta$-axis ($\Delta = 0$) separates saddles from non-saddles. Every $2 \times 2$ linear system lives at exactly one point on this map.

Compare with the second-derivative test: the Hessian’s eigenvalues classify critical points of a function as minima ($\lambda_1, \lambda_2 > 0$), maxima ($\lambda_1, \lambda_2 < 0$), or saddle points (mixed signs). The ODE system matrix’s eigenvalues classify equilibria as stable nodes, unstable nodes, or saddles — the same eigenvalue analysis, governing dynamics instead of statics. Gradient descent connects them: $\dot{\theta} = -\nabla^2 f \cdot \theta$ near a critical point.

A vector $\mathbf{w}$ is a generalized eigenvector of rank 2 for $A$ with eigenvalue $\lambda$ if $(A - \lambda I)\mathbf{w} = \mathbf{v}$, where $\mathbf{v}$ is an eigenvector ($(A - \lambda I)\mathbf{v} = \mathbf{0}$, $\mathbf{v} \neq \mathbf{0}$). More generally, a generalized eigenvector of rank $k$ satisfies $(A - \lambda I)^k \mathbf{w} = \mathbf{0}$ but $(A - \lambda I)^{k-1}\mathbf{w} \neq \mathbf{0}$.

If $A$ has a repeated eigenvalue $\lambda$ with a single eigenvector $\mathbf{v}$ and generalized eigenvector $\mathbf{w}$ (with $(A - \lambda I)\mathbf{w} = \mathbf{v}$), then the general solution is

$\mathbf{y}(t) = c_1 e^{\lambda t}\mathbf{v} + c_2 e^{\lambda t}(t\mathbf{v} + \mathbf{w}).$

$A = \begin{pmatrix} 2 & 1 \\ 0 & 2 \end{pmatrix}$. The repeated eigenvalue is $\lambda = 2$, with a single eigenvector $\mathbf{v} = (1, 0)^T$. The generalized eigenvector satisfies $(A - 2I)\mathbf{w} = \mathbf{v}$, giving $\mathbf{w} = (0, 1)^T$. The solution is

$\mathbf{y}(t) = c_1 e^{2t}\begin{pmatrix} 1 \\ 0 \end{pmatrix} + c_2 e^{2t}\left(t\begin{pmatrix} 1 \\ 0 \end{pmatrix} + \begin{pmatrix} 0 \\ 1 \end{pmatrix}\right).$

The phase portrait is a degenerate unstable node — all trajectories are tangent to the eigenvector direction at the origin.

If $A = PJP^{-1}$ where $J$ is in Jordan normal form with blocks $J_k = \lambda_k I + N_k$ ($N_k$ nilpotent), then $e^{At} = P e^{Jt} P^{-1}$ where

$e^{J_k t} = e^{\lambda_k t} e^{N_k t} = e^{\lambda_k t}\left(I + N_k t + \frac{N_k^2 t^2}{2!} + \cdots\right).$

This is a finite sum because $N_k$ is nilpotent ($N_k^m = 0$ for some $m$). Each Jordan block of size $m$ contributes terms $t^j e^{\lambda t}$ for $j = 0, 1, \ldots, m-1$.

Generically, matrices are diagonalizable — defective matrices form a set of measure zero in matrix space. But in applications, defective matrices arise naturally at bifurcation points — the boundary between qualitatively different behaviors. The transition from a stable spiral to a stable node passes through a defective node at $\tau^2 = 4\Delta$ on the trace-determinant diagram. Stability & Dynamical Systems will examine these transitions systematically.

A fundamental matrix for the system $\mathbf{y}' = A\mathbf{y}$ is a matrix-valued function $\Phi(t)$ satisfying $\Phi'(t) = A\Phi(t)$ whose columns $\boldsymbol{\phi}_1(t), \ldots, \boldsymbol{\phi}_n(t)$ are linearly independent solutions. The standard fundamental matrix is $\Phi(t) = e^{At}$, satisfying $\Phi(0) = I$.

If $\Phi(t)$ is a fundamental matrix for $\mathbf{y}' = A\mathbf{y}$, then

$\det \Phi(t) = \det \Phi(0) \cdot e^{\operatorname{tr}(A) \cdot t}.$

For the standard fundamental matrix, $\det e^{At} = e^{\operatorname{tr}(A) \cdot t}$.

We show that $\frac{d}{dt}\det \Phi(t) = \operatorname{tr}(A) \cdot \det \Phi(t)$ — a scalar linear ODE with solution $\det \Phi(t) = \det \Phi(0) \cdot e^{\operatorname{tr}(A)t}$.

For the derivative of the determinant: by the multilinearity of the determinant in its rows,

$\frac{d}{dt}\det \Phi = \sum_{i=1}^{n} \det \Phi_i(t),$

where $\Phi_i$ replaces the $i$-th row of $\Phi$ with its derivative (the $i$-th row of $\Phi' = A\Phi$).

The $i$-th row of $A\Phi$ is $\sum_j a_{ij} (\text{row } j \text{ of } \Phi)$. By the alternating property of the determinant, a matrix with two identical rows has determinant zero. Therefore, only the diagonal term $a_{ii}$ survives in each $\det \Phi_i(t)$:

$\frac{d}{dt}\det \Phi = \sum_{i=1}^{n} a_{ii} \det \Phi = \operatorname{tr}(A) \det \Phi.$

This is the scalar ODE $W' = \operatorname{tr}(A) \cdot W$ with $W(0) = \det \Phi(0)$, whose unique solution is $W(t) = \det \Phi(0) \cdot e^{\operatorname{tr}(A)t}$. $\blacksquare$

Abel’s formula says that the flow of $\mathbf{y}' = A\mathbf{y}$ preserves volume if and only if $\operatorname{tr}(A) = 0$. Hamiltonian systems (from classical mechanics) have $\operatorname{tr}(A) = 0$ and are volume-preserving — this is Liouville’s theorem. Dissipative systems ($\operatorname{tr}(A) < 0$) contract volumes, which is why stable systems have trajectories that converge to lower-dimensional attractors.

The solution of the IVP $\mathbf{y}' = A\mathbf{y} + \mathbf{g}(t)$, $\mathbf{y}(0) = \mathbf{y}_0$ is:

$\mathbf{y}(t) = e^{At}\mathbf{y}_0 + \int_0^t e^{A(t-s)}\mathbf{g}(s)\,ds.$

The first term is the homogeneous solution (free response); the integral is the convolution of the matrix exponential with the forcing (forced response).

Let $\mathbf{y}(t) = e^{At}\mathbf{c}(t)$ (variation of parameters ansatz). Differentiating using the product rule:

$\mathbf{y}' = Ae^{At}\mathbf{c} + e^{At}\mathbf{c}' = A\mathbf{y} + e^{At}\mathbf{c}'.$

Comparing with $\mathbf{y}' = A\mathbf{y} + \mathbf{g}(t)$ gives $e^{At}\mathbf{c}'(t) = \mathbf{g}(t)$, so $\mathbf{c}'(t) = e^{-At}\mathbf{g}(t)$.

Integrating with the initial condition $\mathbf{c}(0) = \mathbf{y}_0$:

$\mathbf{c}(t) = \mathbf{y}_0 + \int_0^t e^{-As}\mathbf{g}(s)\,ds.$

Substituting back: $\mathbf{y}(t) = e^{At}\mathbf{y}_0 + \int_0^t e^{A(t-s)}\mathbf{g}(s)\,ds$. $\blacksquare$

The system $\mathbf{y}' = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix}\mathbf{y} + \begin{pmatrix} 0 \\ \cos(\Omega t) \end{pmatrix}$ models a harmonic oscillator driven by an external force at frequency $\Omega$. The convolution integral produces resonance when $\Omega = \omega$ — the forced response grows without bound because the driving frequency matches the natural frequency.

Writing $\mathbf{y}' = A\mathbf{y} + B\mathbf{u}$, $\mathbf{z} = C\mathbf{y}$ (state-space form), the output is

$\mathbf{z}(t) = Ce^{At}\mathbf{y}_0 + C\int_0^t e^{A(t-s)}B\mathbf{u}(s)\,ds.$

The matrix $Ce^{At}B$ is the impulse response — the system’s output when the input is a delta function. This is exactly the structure of state-space models in ML (S4/Mamba): the matrices $A$, $B$, $C$ are learned, and the matrix exponential mediates between continuous and discrete representations.

The variation of parameters formula is a matrix convolution. The Laplace transform converts this to an algebraic equation: $\hat{\mathbf{y}}(s) = (sI - A)^{-1}\mathbf{y}_0 + (sI - A)^{-1}\hat{\mathbf{g}}(s)$. The transfer matrix $(sI - A)^{-1}$ encodes the system’s frequency response. The poles of the transfer function are the eigenvalues of $A$ — connecting back to the phase portrait classification.

Let $\mathbf{f}: D \to \mathbb{R}^n$ be continuous on an open set $D \subseteq \mathbb{R} \times \mathbb{R}^n$ and Lipschitz in $\mathbf{y}$:

$\|\mathbf{f}(t, \mathbf{y}_1) - \mathbf{f}(t, \mathbf{y}_2)\| \leq L\|\mathbf{y}_1 - \mathbf{y}_2\| \quad \text{for all } (t, \mathbf{y}_1), (t, \mathbf{y}_2) \in D.$

Then the IVP $\mathbf{y}' = \mathbf{f}(t, \mathbf{y})$, $\mathbf{y}(t_0) = \mathbf{y}_0$ has a unique local solution.

If $A \in \mathbb{R}^{n \times n}$ is constant and $\mathbf{g}: \mathbb{R} \to \mathbb{R}^n$ is continuous, then the IVP $\mathbf{y}' = A\mathbf{y} + \mathbf{g}(t)$, $\mathbf{y}(0) = \mathbf{y}_0$ has a unique solution on all of $\mathbb{R}$.

This follows because $\mathbf{f}(\mathbf{y}) = A\mathbf{y}$ is globally Lipschitz with constant $\|A\|$, and the Grönwall inequality prevents blow-up: $\|\mathbf{y}(t)\| \leq e^{\|A\|t}\|\mathbf{y}_0\| + \int_0^t e^{\|A\|(t-s)}\|\mathbf{g}(s)\|\,ds$.

Compare with Topic 21’s blow-up example $y' = y^2$, where the nonlinearity causes finite-time escape to infinity. For linear systems, the solution $\mathbf{y}(t) = e^{At}\mathbf{y}_0$ grows at most exponentially: $\|\mathbf{y}(t)\| \leq e^{\|A\|t}\|\mathbf{y}_0\|$. Exponential growth is fast, but it never reaches infinity in finite time. This is why linearization is always locally valid — the linear approximation never predicts blow-up.

In the phase plane, solution trajectories of $\mathbf{y}' = A\mathbf{y}$ never cross (except at the origin, which is an equilibrium). If two trajectories crossed at a point $\mathbf{y}^*$ at time $t^*$, that point would have two distinct solutions to the same IVP — contradicting uniqueness. This topological consequence constrains the geometry of phase portraits: trajectories can approach or depart from the origin, but they cannot intersect.

Consider $L(\theta) = \frac{1}{2}\theta^T H \theta$ with $H = \text{diag}(1, 10, 100)$. The gradient flow $\dot{\theta} = -H\theta$ has solution $\theta_i(t) = \theta_{i,0} e^{-\lambda_i t}$.

After time $t$, the residual error in each direction is $e^{-\lambda_i t}$: direction 3 converges 100× faster than direction 1. The loss $L(t) = \frac{1}{2}\sum \theta_{i,0}^2 e^{-2\lambda_i t}$ shows multi-scale decay — a fast initial drop (the large-eigenvalue components) followed by a long tail (the small-eigenvalue component).

Given continuous parameters $A, B$, the discrete parameters are $\bar{A} = e^{A\Delta t}$ and $\bar{B} = (e^{A\Delta t} - I)A^{-1}B$ (assuming $A$ is invertible). This exact discretization preserves the dynamics of the continuous system. The bilinear (Tustin) approximation $\bar{A} = (I + A\Delta t/2)(I - A\Delta t/2)^{-1}$ is a Padé approximant that preserves stability — mapping eigenvalues with negative real part to eigenvalues inside the unit circle.