Matrix Functions and Special Matrices

The matrix exponential of $\mathbf{A} \in \mathbb{C}^{n \times n}$ is:

$$e^{\mathbf{A}} = \sum_{k=0}^{\infty} \frac{\mathbf{A}^k}{k!} = \mathbf{I} + \mathbf{A} + \frac{\mathbf{A}^2}{2!} + \frac{\mathbf{A}^3}{3!} + \cdots$$

from scipy.linalg import expm
E = expm(A)

SciPy uses the Pade approximation with scaling and squaring (the Al-Mohy & Higham algorithm), which is numerically stable.

Key property: $e^{\mathbf{A}} e^{\mathbf{B}} = e^{\mathbf{A}+\mathbf{B}}$ only when $\mathbf{A}\mathbf{B} = \mathbf{B}\mathbf{A}$ (i.e., $\mathbf{A}$ and $\mathbf{B}$ commute).

The matrix logarithm $\log(\mathbf{A})$ satisfies $e^{\log(\mathbf{A})} = \mathbf{A}$ (when $\mathbf{A}$ has no nonpositive real eigenvalues).

The matrix square root $\mathbf{A}^{1/2}$ satisfies $\mathbf{A}^{1/2}\mathbf{A}^{1/2} = \mathbf{A}$ (for positive definite $\mathbf{A}$).

from scipy.linalg import logm, sqrtm

L = logm(A)           # matrix logarithm
S = sqrtm(A)          # principal matrix square root

The matrix square root of a covariance matrix $\mathbf{R}$ is used to generate correlated random vectors: $\mathbf{x} = \mathbf{R}^{1/2}\mathbf{z}$ where $\mathbf{z} \sim \mathcal{CN}(\mathbf{0}, \mathbf{I})$.

A Toeplitz matrix has constant diagonals: $T_{ij} = t_{i-j}$

$$\mathbf{T} = \begin{bmatrix} t_0 & t_{-1} & t_{-2} & \cdots \\ t_1 & t_0 & t_{-1} & \cdots \\ t_2 & t_1 & t_0 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}$$

from scipy.linalg import toeplitz

c = [2, 1, 0, 0]     # first column
r = [2, -1, 0, 0]    # first row
T = toeplitz(c, r)

Toeplitz matrices represent linear time-invariant (LTI) filters: the output of a discrete convolution $y[n] = \sum_k h[k] x[n-k]$ can be written as $\mathbf{y} = \mathbf{T}_h \mathbf{x}$ where $\mathbf{T}_h$ is Toeplitz with the filter coefficients.

A circulant matrix is a Toeplitz matrix where each row is a cyclic shift of the previous row:

$$\mathbf{C} = \begin{bmatrix} c_0 & c_{n-1} & c_{n-2} & \cdots & c_1 \\ c_1 & c_0 & c_{n-1} & \cdots & c_2 \\ c_2 & c_1 & c_0 & \cdots & c_3 \\ \vdots & & & \ddots & \vdots \\ c_{n-1} & c_{n-2} & c_{n-3} & \cdots & c_0 \end{bmatrix}$$

Every circulant matrix is diagonalized by the DFT matrix: $$\mathbf{C} = \mathbf{F}^H \mathrm{diag}(\mathbf{F}\mathbf{c}) \mathbf{F}$$

This is why circular convolution can be computed via FFT in $O(n \log n)$ instead of $O(n^2)$.

from scipy.linalg import circulant
C = circulant([1, 2, 3, 4])

The $n \times n$ DFT matrix $\mathbf{F}$ has entries:

$$F_{kl} = \frac{1}{\sqrt{n}} e^{-j 2\pi kl/n}, \quad k, l = 0, \ldots, n-1$$

The DFT matrix is unitary: $\mathbf{F}\mathbf{F}^H = \mathbf{I}$.

from scipy.linalg import dft

F = dft(n, scale='sqrtn')   # unitary DFT matrix
# Or construct manually:
k = np.arange(n)
F_manual = np.exp(-1j * 2 * np.pi * np.outer(k, k) / n) / np.sqrt(n)

In practice, you never form $\mathbf{F}$ explicitly — the FFT computes $\mathbf{F}\mathbf{x}$ in $O(n \log n)$ without the matrix:

X = np.fft.fft(x) / np.sqrt(n)   # equivalent to F @ x