∇ ANNA ZHOU / RESEARCH NOTES
A visual research tutorial · 14 notebook demonstrations

Matrix Functions
as Dynamics

From \(e^A\) and \(\log A\) to matrix sine, cosine, Fréchet derivatives, and stable numerical computation—explained through concise experiments.

Higham · Chapters 10–12Numerical linear algebraNumPy / SciPyReproducible figures
OverviewMaster ideaMatrix exponentialFréchet derivativesStable algorithmsMatrix logarithmSine & cosineSO(2) demoReproduce
Overview

One framework, four matrix functions.

The exponential propagates dynamics, the logarithm recovers generators on a chosen branch, and sine/cosine describe oscillation. Fréchet derivatives explain sensitivity; Padé approximation and scaling provide practical algorithms.

3.01×Peak state norm in the stable nonnormal example (from initial norm 1)
5.52 × 10⁻¹⁴Relative error of an illustrative Padé exponential computation
≈10⁻¹³Cosine/sine double-angle errors at 8 halvings

Numbers are measured outputs of the supplied educational notebook, not performance guarantees.

Mathematical foundation

The master equation

\[e^A=\sum_{k\ge0}\frac{A^k}{k!},\qquad e^{\log A}=A,\qquad e^{iA}=\cos A+i\sin A.\]

The last two identities use the principal logarithm where it exists; \(\log(e^A)=A\) additionally needs a compatible branch.

\[f(A+E)=f(A)+L_f(A,E)+o(\|E\|).\]
01

Define \(f(A)\)

02

Analyze sensitivity

03

Approximate locally

04

Reconstruct and verify

Experiments 01–03 · Chapter 10

Eigenvalues do not tell the whole story.

A Jordan block adds polynomial factors to the exponential. A stable but nonnormal matrix can show large transient growth before decaying. Noncommuting operators also invalidate scalar addition rules.

\[e^{t(\lambda I+N)}=e^{\lambda t}\sum_{k=0}^{m-1}\frac{t^kN^k}{k!},\qquad x(t)=e^{tA}x_0.\]
Figure 1 · Jordan off-diagonal polynomial factor and transient norm growth.
Figure 1 · Jordan off-diagonal polynomial factor and transient norm growth. ↗ expand
Figure 2 · Stable normal versus nonnormal state evolution.
Figure 2 · Stable normal versus nonnormal state evolution. ↗ expand

Composition requires commutators

\[\log(e^Ae^B)=A+B+\tfrac12[A,B]+\cdots,\qquad[A,B]=AB-BA.\]
Figure 3 · BCH commutator correction reduces local composition error.
Figure 3 · BCH commutator correction reduces local composition error. ↗ expand
Takeaway. Modal decay, transient amplification, and operator composition are distinct phenomena; nonnormality and commutators matter.
Experiments 04–05 · Chapter 10

Matrix derivatives are linear operators.

\[L_{\exp}(A,E)=\int_0^1 e^{(1-s)A}Ee^{sA}\,ds,\qquad \kappa_f(A)=\frac{\|L_f(A)\|\,\|A\|}{\|f(A)\|}.\]

In general \(L_{\exp}(A,E)\ne e^AE\). The notebook compares integral quadrature, exact Fréchet routines, finite differences, and the worst perturbation direction.

Figure 4 · Finite-difference derivative error versus step size.
Figure 4 · Finite-difference derivative error versus step size. ↗ expand
Figure 5 · Most sensitive perturbation direction.
Figure 5 · Most sensitive perturbation direction. ↗ expand
Takeaway. Finite differences exhibit a truncation-versus-rounding tradeoff; the most sensitive direction need not be an obvious perturbation.
Experiment 06 · Chapter 10

Scale → approximate → reconstruct.

Compute \(e^{A/2^s}\) near the origin using a Padé approximant, then square repeatedly.

\[e^A=\big(e^{A/2^s}\big)^{2^s}\approx\big(r_m(A/2^s)\big)^{2^s}.\]
Figure 6 · Padé scaling-and-squaring versus unscaled Taylor.
Figure 6 · Padé scaling-and-squaring versus unscaled Taylor. ↗ expand

The notebook uses an educational [6/6] Padé implementation, not a production-grade degree/scaling selector.

Experiments 07–10 · Chapter 11

The matrix logarithm has a branch.

The principal logarithm requires no eigenvalues on the closed negative real axis. In rotation examples, it wraps angles into \((-\pi,\pi)\).

\[\log A=\int_0^1(A-I)[I+t(A-I)]^{-1}\,dt.\]
Figure 7 · Principal rotation angle wraps by multiples of 2π.
Figure 7 · Principal rotation angle wraps by multiples of 2π. ↗ expand
Figure 8 · Mercator, Gregory and resolvent quadrature convergence.
Figure 8 · Mercator, Gregory and resolvent quadrature convergence. ↗ expand
Figure 9 · Square-root reduction and fixed-degree log approximation.
Figure 9 · Square-root reduction and fixed-degree log approximation. ↗ expand
\[\log A=2^s\log\big(A^{1/2^s}\big),\qquad L_{\exp}(\log A,L_{\log}(A,E))=E.\]
Takeaway. Repeated square roots move the argument toward \(I\), where local approximations of \(\log\) are accurate. The notebook verifies the derivative inverse-chain identity to about \(5.38\times10^{-16}\) relative residual.
Experiments 11–13 · Chapter 12

Oscillation is exponential dynamics in disguise.

\[\cos A=\frac{e^{iA}+e^{-iA}}2,\quad\sin A=\frac{e^{iA}-e^{-iA}}{2i},\quad \cos^2A+\sin^2A=I.\]
Figure 10 · Oscillatory eigenmodes of matrix sine and cosine.
Figure 10 · Oscillatory eigenmodes of matrix sine and cosine. ↗ expand
Figure 11 · Double-angle error versus number of halvings.
Figure 11 · Double-angle error versus number of halvings. ↗ expand

Double-angle reconstruction uses \(C_{k+1}=2C_k^2-I\) and \(S_{k+1}=2S_kC_k\). Too little scaling causes approximation error; too much can amplify roundoff.

\[\ddot x+Kx=0,\quad \Omega=K^{1/2},\quad x(t)=\cos(t\Omega)x_0+\Omega^{-1}\sin(t\Omega)v_0.\]
Figure 12 · Equivalent oscillator solutions and energy conservation.
Figure 12 · Equivalent oscillator solutions and energy conservation. ↗ expand
Takeaway. Sine/cosine and first-order state-transition exponentials give equivalent oscillator solutions; the notebook's maximum state discrepancy was \(1.77\times10^{-13}\).
Experiment 14 · Lie group application

Try the matrix exponential as a rotation.

\[J=\begin{bmatrix}0&-1\\1&0\end{bmatrix},\qquad e^{\theta J}=R(\theta)=\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}.\]

Rotate a unit vector

Adjust \(\theta\) to see \(e^{\theta J}\) in action. The principal logarithm recovers the wrapped angle, not the unbounded input.

Principal angle: 45.0°

R = [[0.707, −0.707], [0.707, 0.707]]

xy
Figure 13 · Notebook-generated SO(2) unit-vector rotations.
Figure 13 · Notebook-generated SO(2) unit-vector rotations. ↗ expand
Figure 14 · Eigenvector conditioning and spectral evaluation error.
Figure 14 · Eigenvector conditioning and spectral evaluation error. ↗ expand
Reproducibility

Run every experiment.

The article's plots come from the executed notebook supplied with this project. Figures illustrate particular 2×2 matrices and educational algorithms—not new theorems or general benchmarks.

Dependencies: NumPy, SciPy, Matplotlib · CPU only · no external dataset.

See all 14 experiment figures

Reference & attribution

Nicholas J. Higham, Functions of Matrices: Theory and Computation, SIAM, 2008. Chapters 10 (Matrix Exponential), 11 (Matrix Logarithm), and 12 (Matrix Cosine and Sine). Numerical experiments and plots are from the accompanying educational Colab notebook, not figures reproduced from the book.

Design direction follows the structure of using original page styling adapted for the previous topic.