Gradient Verification for Matrix-Valued Likelihoods

Why Verify a Matrix Gradient?

An optimizer can only use the gradient it is given. If an analytic gradient has a sign error, a missing factor, or an incorrect treatment of matrix symmetry, the resulting optimization may be inaccurate even when the likelihood itself is evaluated stably. A numerical check compares the proposed gradient with changes in the objective under small, controlled perturbations.

For a matrix parameter, checking every entry separately is often less informative than checking directional derivatives. A direction tests a weighted combination of entries, and several well-chosen directions can reveal errors that a single coordinate check misses. The challenge is to choose perturbations large enough to be distinguishable from floating-point noise, but small enough that curvature does not dominate.

Definition: For a differentiable scalar function $\ell(A)$ of a real matrix $A$, its gradient $G$ under the Frobenius inner product is defined by $$ d\ell(A)[H]=\langle G,H\rangle_F=\operatorname{tr}(G^T H) $$ for every admissible perturbation direction $H$. Thus the directional derivative in direction $H$ is $\operatorname{tr}(G^T H)$. The admissible directions are determined by the parameter space; for a symmetric matrix parameter, they are symmetric matrices.

The gradient is therefore tied to a convention: the displayed matrix $G$ represents derivatives under the Frobenius inner product. For a symmetric parameter, only its symmetric part affects derivatives along symmetric directions. Indeed, if $H=H^T$, then $\operatorname{tr}(G^T H)=\operatorname{tr}(((G+G^T)/2)^T H)$. A candidate gradient can be checked in any admissible direction by comparing its inner product with a finite-difference estimate.

Directional Finite Differences and Their Error

For a direction $H$, the central directional difference at step $h>0$ is $$ D_h(H)=\frac{\ell(A+hH)-\ell(A-hH)}{2h}. $$ The analytic gradient predicts the value $\langle G,H\rangle_F$. Central differences are generally preferable to forward differences because their leading truncation error is quadratic in $h$, rather than linear.

Result: Suppose $t\mapsto\ell(A+tH)$ is three times continuously differentiable for $-h\leq t\leq h$, and the absolute value of its third derivative is at most $M$ on this interval. Then $$ \left|D_h(H)-\langle G,H\rangle_F\right|\leq \frac{Mh^2}{6}. $$

To prove the bound, write $g(t)=\ell(A+tH)$. Taylor’s theorem about zero gives $g(h)=g(0)+hg'(0)+h^2g''(0)/2+r_+$ and $g(-h)=g(0)-hg'(0)+h^2g''(0)/2+r_-$, where $|r_+|\leq Mh^3/6$ and $|r_-|\leq Mh^3/6$. Subtract the expansions and divide by $2h$. The constant and quadratic terms cancel, while the two remainder terms contribute at most $(|r_+|+|r_-|)/(2h)\leq Mh^2/6$. Since $g'(0)=\langle G,H\rangle_F$, the result follows. $\square$

This bound describes truncation error, not all numerical error. If each computed likelihood value has absolute error at most $\varepsilon C$, where $C$ describes the scale of its evaluation error, differencing two computed values adds an error of at most $\varepsilon C/h$. A useful qualitative model is therefore $$ \text{total error}\ \lesssim\ \frac{Mh^2}{6}+\frac{\varepsilon C}{h}. $$ A step that is too large emphasizes curvature; one that is too small amplifies evaluation noise. Balancing these terms gives $h$ of order $(\varepsilon C/M)^{1/3}$. In practice, $M$ and $C$ are rarely known, so a step-size sweep is more informative than treating one step as definitive.

Worked Examples

Worked Example: Checking a Gaussian Variance Derivative

Let $X_1,X_2,X_3$ be independent $N(0,v)$ observations with values $-1,0,2$. Ignoring terms constant in $v$, the log-likelihood is $\ell(v)=-\tfrac{3}{2}\log v-\tfrac{5}{2v}$, since the sum of squares is $1+0+4=5$. At $v=2$, its derivative is $$ \ell'(2)=-\frac{3}{2(2)}+\frac{5}{2(2^2)}=-\frac34+\frac58=-\frac18=-0.125. $$

Using $h=0.01$, the central difference is $$ D_h=\!-\frac{3}{2}\frac{\log(2.01)-\log(1.99)}{0.02} -\frac{5}{2}\frac{1/2.01-1/1.99}{0.02}. $$ The first quotient is approximately $0.500004167$. The second term simplifies to $2.5/(4-h^2)=2.5/3.9999$, approximately $0.625015625$. Thus $D_h\approx-0.750006250+0.625015625=-0.124990625$, close to $-0.125$. The discrepancy is consistent with a small central-difference truncation error; it is not evidence of a wrong derivative.

Worked Example: A Directional Check for a Covariance Matrix

For $n$ centered multivariate normal observations with covariance parameter $\Sigma$ and scatter matrix $S$, the part of the log-likelihood depending on $\Sigma$ is $$ \ell(\Sigma)=-\frac{n}{2}\log\det(\Sigma)-\frac12\operatorname{tr}(\Sigma^{-1}S). $$ Assume $\Sigma$ is positive definite and $S$ is symmetric. The differential identities $d\log\det(\Sigma)=\operatorname{tr}(\Sigma^{-1}d\Sigma)$ and $d\Sigma^{-1}=-\Sigma^{-1}(d\Sigma)\Sigma^{-1}$ give $$ d\ell=-\frac{n}{2}\operatorname{tr}(\Sigma^{-1}d\Sigma) +\frac12\operatorname{tr}(\Sigma^{-1}(d\Sigma)\Sigma^{-1}S). $$ Cyclically rearranging the second trace and using symmetry identifies the Frobenius gradient as $G=-\tfrac n2\Sigma^{-1}+\tfrac12\Sigma^{-1}S\Sigma^{-1}$.

Take $n=4$, $\Sigma=\begin{pmatrix}2&0\\0&1\end{pmatrix}$, and $S=\begin{pmatrix}4&1\\1&2\end{pmatrix}$. Then $\Sigma^{-1}=\begin{pmatrix}1/2&0\\0&1\end{pmatrix}$ and $\Sigma^{-1}S\Sigma^{-1}=\begin{pmatrix}1&1/2\\1/2&2\end{pmatrix}$, so $$ G=\begin{pmatrix}-1/2&1/4\\1/4&-1\end{pmatrix}. $$ Check the off-diagonal direction $H=\begin{pmatrix}0&1\\1&0\end{pmatrix}$. The predicted derivative is $\langle G,H\rangle_F=1/4+1/4=1/2$; both off-diagonal entries contribute.

For $\Sigma+tH=\begin{pmatrix}2&t\\t&1\end{pmatrix}$, the determinant is $2-t^2$, which is positive for $|t|<\sqrt{2}$, so these perturbations remain positive definite. Direct inversion gives $\operatorname{tr}((\Sigma+tH)^{-1}S)=(8-2t)/(2-t^2)$. The determinant contribution to $\ell(t)$ is even in $t$, so it cancels in the central difference. The remaining term yields $$ \frac{\ell(h)-\ell(-h)}{2h}=\frac{1}{2-h^2}. $$ At $h=0.01$, this is approximately $0.500025001$, close to the predicted $0.5$. The calculation also shows why an off-diagonal check must account for both symmetric entries.

Worked Example: Choosing Steps for Parameters on Different Scales

Consider one bivariate normal observation $x=(1000,0.001)^T$ with known covariance $\operatorname{diag}(10^6,10^{-6})$ and mean $\mu=(a,b)^T$. At $a=1000.2$, $b=0.0015$, the log-likelihood, apart from a constant, is $$ \ell(a,b)=-\frac12\left\{\left(\frac{a-1000}{1000}\right)^2+ \left(\frac{b-0.001}{0.001}\right)^2\right\}. $$ Differentiation gives $\partial\ell/\partial a=-0.2/10^6=-2\times10^{-7}$ and $\partial\ell/\partial b=-0.0005/10^{-6}=-500$. These derivatives differ greatly because the parameters have different units and scales.

A scale-aware starting choice is $h_j=\varepsilon^{1/3}s_j$, where $\varepsilon$ is machine precision and $s_j$ is a meaningful scale for coordinate $j$. With $\varepsilon$ approximately $2.22\times10^{-16}$, $\varepsilon^{1/3}$ is approximately $6.06\times10^{-6}$. Taking $s_a=1000$ and $s_b=0.001$ gives $h_a\approx0.00606$ and $h_b\approx6.06\times10^{-9}$. The expected difference between the two objective evaluations is approximately $2h_j|\partial\ell/\partial j|$: about $2.42\times10^{-9}$ for $a$ and $6.06\times10^{-6}$ for $b$. Both perturbations are small relative to their parameter scales, while their signals are not needlessly suppressed.

A common absolute step of $10^{-9}$ would produce an expected difference of only about $4\times10^{-16}$ for $a$, comparable to rounding changes in objective values of ordinary magnitude. The same step has a very different relative size for $b$. The example illustrates why a single absolute step for every matrix entry or parameter is often a poor default. It does not guarantee that the scale-aware steps are optimal: repeating the check at nearby step sizes is still necessary.

A Practical Verification Procedure

A useful check separates the mathematical comparison from the numerical conditions under which it is made. First choose admissible directions. For a symmetric covariance matrix, use symmetric directions; for a positive-definite parameter, ensure that both perturbed matrices remain in the positive-definite region. A sufficient condition is $h\|H\|_2<\lambda_{\min}(\Sigma)$, where $\|H\|_2$ is the spectral norm and $\lambda_{\min}(\Sigma)$ is the smallest eigenvalue. Alternatively, perturb unconstrained parameters in a factorization or other valid parameterization, and compare derivatives in those coordinates.

Set the gradient convention.
Confirm that the candidate $G$ satisfies $d\ell[H]=\operatorname{tr}(G^TH)$, and restrict $H$ to directions allowed by the parameter space.

Choose scaled directions.
Normalize directions or set their entrywise magnitudes relative to parameter scales. Use more than one direction, including directions that involve off-diagonal entries when relevant.

Compare over a step-size sweep.
For each direction, compare the finite difference with $\operatorname{tr}(G^TH)$ at several decreasing and increasing steps. Look for a range where the estimates agree before very small steps become unstable.

Interpret the pattern.
A discrepancy that remains similar across a stable range suggests a derivative or convention error. A discrepancy that grows mainly as $h$ shrinks suggests floating-point noise; one that grows as $h$ increases suggests truncation error or a nearby domain boundary.

For a formal summary across directions, an absolute discrepancy is $e_{\mathrm{abs}}=|D_h(H)-\langle G,H\rangle_F|$. A relative discrepancy can be useful when the derivative is not near zero, for example $e_{\mathrm{rel}}=e_{\mathrm{abs}}/\max(1,|D_h(H)|,|\langle G,H\rangle_F|)$. When the true directional derivative is close to zero, the absolute discrepancy is more interpretable: a relative error can become large simply because its denominator is small. In either case, the pattern across steps and directions matters more than one tolerance applied blindly.

Common Pitfalls and Why the Check Matters

Using only one step size. A small discrepancy at one step can occur by chance; a large discrepancy can be caused by either truncation or roundoff. Check a range of steps.
Forgetting the matrix inner product. For symmetric off-diagonal entries, $\operatorname{tr}(G^TH)$ includes both positions. Comparing only one entry can create an apparent factor-of-two error.
Perturbing outside the parameter space. A covariance perturbation that is not positive definite may make the likelihood undefined. Choose admissible directions and steps, or work in valid coordinates.
Treating every mismatch as a calculus error. If reducing $h$ makes the check less stable, floating-point noise is a plausible cause. If increasing $h$ worsens agreement, curvature may be dominating.
Confusing verification with proof of correctness. Finite differences test the implementation at selected points and directions; they do not establish correctness everywhere or validate the likelihood itself.

Key takeaway: A reliable matrix-gradient check compares the proposed Frobenius directional derivative with central finite differences in admissible, appropriately scaled directions. Use several directions and a step-size sweep to distinguish derivative mistakes from truncation and floating-point noise.

Check Your Understanding

Use the definitions and checks developed above to answer the following.

What does the Frobenius-gradient convention say about the directional derivative in direction $H$?
Why does the central difference have a leading truncation error of order $h^2$, rather than order $h$?
For a symmetric matrix parameter, why can a directional check expose a missing factor of two in off-diagonal derivatives?
What pattern across a step-size sweep would make floating-point noise more plausible than a persistent analytic-gradient error?
Give one condition that ensures a small covariance perturbation remains positive definite, and explain why this matters for likelihood evaluation.