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, producing formulas such as \[ d(X^3) \;=\; (dX)X^2 + X(dX)X + X^2(dX), \qquad d(X^{-1}) \;=\; -X^{-1}(dX)X^{-1}. \] Symbolic derivations like these are short but error-prone — miss a term, 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 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\) is the linear part of the increment — equivalently, the quantity that satisfies \[ f(X + dX) - f(X) \;=\; df + o(\|dX\|). \] Two standard finite-difference approximations of \(df\) are \[ \text{Forward difference:} \quad f(X + dX) - f(X), \qquad \text{Backward difference:} \quad f(X) - f(X - dX). \] Each equals \(df\) up to an \(o(\|dX\|)\) term, so in the limit \(dX \to 0\) both agree with the analytical formula.

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 — 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. Random non-integer values exercise the full logic.

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}\) — 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 indicate either a bug in the analytical formula or an inappropriate perturbation scale, and the validation would have caught it.