Intro to Numerical Computation

Example: numerical validation of \(d(X^2) = X\,(dX) + (dX)\,X\)

In the previous page we derived the differential of matrix-valued matrix functions symbolically. That derivation produced formulas such as \[ \begin{align*} d(X^3) &= (dX)X^2 + X(dX)X + X^2(dX), \\\\ d(X^{-1}) &= -X^{-1}(dX)X^{-1}. \end{align*} \]

Symbolic derivations like these are short but error-prone. Miss a term or flip a factor, and the whole downstream optimization is silently wrong. A basic hygiene step of computational mathematics is therefore to validate each derivation numerically before trusting it.

In this page we work through the simpler case \(f(X) = X^2\), for which, expanding \((X + dX)^2 - X^2\) as in the cubic case, the differential is \[ df = f'(X)[dX] = X(dX) + (dX)X, \] and we verify it against finite differences.

Finite differences

For a differentiable function \(f\), the differential \(df\) at \(X\) (matrices being identified with the vectors of their entries) is the linear part of the increment, that is, \(df = f'(X)[dX]\) for the unique linear map \(f'(X)\) that satisfies \[ f(X + dX) - f(X) = df + o(\|dX\|). \]

Two standard finite-difference approximations of \(df\) are \[ \begin{align*} \text{Forward difference:} &\quad f(X + dX) - f(X), \\\\ \text{Backward difference:} &\quad f(X) - f(X - dX). \end{align*} \] Each equals \(df\) up to an \(o(\|dX\|)\) term. So if \(dX = tV\) for a fixed direction \(V\) with \(f'(X)[V] \neq 0\), the difference from \(df\), divided by \(\|df\|\), tends to \(0\) as \(t \to 0\).

When comparing a computed approximation against a reference value, we measure quality by relative error in a norm (we use the Frobenius norm throughout): \[ \text{relative error} = \frac{\|\,\text{true} - \text{approx}\,\|_F}{\|\,\text{true}\,\|_F}. \] A small absolute difference does not necessarily mean a small error, since scale matters. Relative error normalizes by the size of the quantity being approximated.

Good practice for test inputs

Draw test matrices from a continuous distribution (here, np.random.randn) rather than hand-picking "nice" integer-valued matrices. Integer matrices can accidentally satisfy identities (commuting \(X\) and \(dX\), exact zero entries, singular matrices) that mask bugs in the general-position code path. For Gaussian entries, as drawn here, each of these coincidences occurs with probability zero.

Implementation and run

import numpy as np

# Random test matrix X in R^{2x2}
X = np.random.randn(2, 2)

# Perturbation dX at the sqrt(eps_machine) scale
dX = np.random.randn(2, 2) * 1e-8

# f(X) = X^2
f_X = X @ X

# Forward difference:  f(X + dX) - f(X)
f_approx = (X + dX) @ (X + dX) - f_X

# Backward difference: f(X) - f(X - dX)
b_approx = f_X - (X - dX) @ (X - dX)

# Analytical differential:  df = X dX + dX X
exact = X @ dX + dX @ X

# Relative errors in Frobenius norm
f_rel_err = np.linalg.norm(f_approx - exact, 'fro') / np.linalg.norm(exact, 'fro')
b_rel_err = np.linalg.norm(b_approx - exact, 'fro') / np.linalg.norm(exact, 'fro')

print(f"Input matrix X:\n{X}\n")
print(f"Perturbation dX:\n{dX}\n")
print(f"Forward difference  f(X + dX) - f(X):\n{f_approx}\n")
print(f"Backward difference f(X) - f(X - dX):\n{b_approx}\n")
print(f"Analytical df = X dX + dX X:\n{exact}\n")
print(f"Relative error (forward) : {f_rel_err:.2e}")
print(f"Relative error (backward): {b_rel_err:.2e}")

Typical output shows relative errors around \(10^{-8}\). That size is consistent with the optimal accuracy predicted by the balance of truncation and rounding errors discussed above. If one of these errors suddenly spiked to \(10^{-2}\) or larger, that would most likely point to one of three causes. The analytical formula may contain a bug, the perturbation scale may be inappropriate, or the direction \(dX\) may be unlucky, with \(df\) nearly vanishing along it. Rerunning with a fresh random \(dX\) separates the unlucky direction from the other two causes. If the spike persists at a sound scale, the validation has most likely caught a bug in the formula.