Numerical Ill-Conditioning Versus Statistical Uncertainty

Two Different Reasons an Estimate Can Look Imprecise

A large standard error can indicate that the data genuinely provide little information about a parameter. It can also arise because the numerical calculation used to estimate that standard error is unstable. These are different problems: collecting more informative data may reduce genuine statistical uncertainty, while changing numerical scaling or improving the Hessian calculation may resolve numerical instability. A poorly conditioned calculation can also obscure genuine uncertainty, so the two explanations should not be treated as mutually exclusive.

In population pharmacokinetics, precision summaries commonly come from an estimated covariance matrix of the parameter estimates. That matrix is often calculated by inverting a Hessian, the matrix of second derivatives of an objective function at the fitted estimates. Inverting a matrix that is nearly singular can amplify small numerical errors. Conversely, a stable calculation can reveal a genuinely flat likelihood direction: the available data do not sharply distinguish certain parameter combinations.

The earlier material on condition numbers and weakly estimated parameter combinations provides useful vocabulary for describing these patterns. Here the emphasis is on deciding what a problematic precision estimate means, and on how convergence, scaling, and covariance conventions affect that decision.

Definition: Numerical ill-conditioning occurs when small perturbations in a matrix or its entries can cause large changes in a computed solution, such as an inverse or a covariance estimate. Statistical uncertainty is variability in an estimate due to limited information in the data under the assumed model. Ill-conditioning concerns the stability of a calculation; statistical uncertainty concerns what the data and model can estimate.

What the Inverse Hessian Represents

Let \(\hat{\theta}\) denote a fitted parameter vector, and let \(\ell(\theta)\) be the log likelihood. Under a regular likelihood model, the observed information matrix is the negative Hessian of the log likelihood at the estimate:

H_{\ell}(\hat{\theta})=-\nabla^2\ell(\hat{\theta})

When the likelihood is sufficiently smooth, the estimate is an interior optimum, the parameters are locally identifiable, and the sample contains adequate information, a quadratic approximation to the log likelihood gives the approximate covariance matrix

\widehat{\operatorname{Cov}}(\hat{\theta})\approx H_{\ell}(\hat{\theta})^{-1}.

This is an approximation, not a direct measurement of uncertainty. Its reliability depends both on the statistical approximation and on computing the Hessian and inverse accurately. The approximation can be poor near a parameter boundary, in a strongly non-quadratic likelihood, or when the model is not identifiable. A matrix can also be mathematically invertible but so ill-conditioned that its numerical inverse is unreliable.

Derivation: At an interior maximum, the score is zero, and a second-order Taylor expansion gives \(\ell(\theta)\approx\ell(\hat{\theta})-\frac{1}{2}(\theta-\hat{\theta})^T H_{\ell}(\hat{\theta})(\theta-\hat{\theta})\). Exponentiating this quadratic approximation gives a normal-shaped likelihood with covariance \(H_{\ell}(\hat{\theta})^{-1}\), provided the curvature matrix is positive definite. Thus small curvature corresponds to large variance along that direction. This argument justifies the local approximation; it does not guarantee that the fitted Hessian or its numerical inverse is accurate.

Always identify which objective was differentiated. If \(Q(\theta)=-2\ell(\theta)\), then its Hessian is twice the negative-log-likelihood Hessian, and the corresponding approximation is \(2[\nabla^2 Q(\hat{\theta})]^{-1}\). Applying the inverse-Hessian formula without checking the objective's scale can therefore introduce a factor-of-two error in the covariance.

Scaling Can Improve Computation Without Adding Information

A condition number depends on the coordinates and units used to represent a parameter vector. Changing from parameters measured in litres and litres per hour to dimensionless scaled parameters can make a matrix much easier to invert. If the transformation is exact and the covariance is transformed back correctly, it does not change the uncertainty represented in the original parameters. Scaling improves numerical behavior; it does not create information in the data.

Worked Example: A Change of Scale, Not a Change in Precision

Suppose a hypothetical fit estimates clearance \(CL=5\) L/h and volume \(V=50\) L. The Hessian of the negative log likelihood, in these original coordinates, is \(H_{\ell}=\operatorname{diag}(100,1)\). For this calculation, assume the Hessian is accurate, positive definite, and the local quadratic approximation is adequate.

Inverting gives

\widehat{\operatorname{Cov}}(CL,V)=H_{\ell}^{-1} =\begin{pmatrix}0.01&0\\0&1\end{pmatrix}.

The standard errors are \(SE(CL)=\sqrt{0.01}=0.1\) L/h and \(SE(V)=\sqrt{1}=1\) L. The corresponding relative standard errors are \(0.1/5=0.02\), or 2%, for clearance and \(1/50=0.02\), or 2%, for volume.

Now define dimensionless coordinates \(q_1=CL/10\) and \(q_2=V/100\), so that \((CL,V)^T=Dq\), with \(D=\operatorname{diag}(10,100)\). The transformed Hessian is \(H_q=D^T H_{\ell}D=\operatorname{diag}(10000,10000)\). Its condition number is 1, compared with 100 for the original diagonal Hessian. In the new coordinates the covariance is \(H_q^{-1}=\operatorname{diag}(0.0001,0.0001)\). Transforming it back gives

D H_q^{-1}D^T =\begin{pmatrix}0.01&0\\0&1\end{pmatrix}.

The original standard errors and relative standard errors are unchanged. This example assumes an exact linear rescaling and a correctly transformed covariance. It shows why a better condition number after scaling is evidence of improved numerical representation, not evidence that the data became more informative.

Weak Curvature Can Be Statistical, Numerical, or Both

A small eigenvalue of the Hessian means that curvature is small along its corresponding parameter combination. If this persists in a reliable likelihood calculation, it signals weak information in that direction and hence genuine statistical uncertainty. If the eigenvalue changes greatly with finite-difference step size, numerical precision, or small perturbations to the Hessian, the estimated uncertainty may also be numerically unstable. The size of a standard error alone cannot distinguish these explanations.

Worked Example: A Weakly Estimated Parameter Combination

Suppose a hypothetical fit has a positive-definite Hessian with eigenvalues 25 and 0.002. Assume, for illustration, that the two unit-length eigenvectors are \((1,1)^T/\sqrt{2}\) and \((1,-1)^T/\sqrt{2}\), respectively, in two fitted log-parameter coordinates. The covariance eigenvalues are the reciprocals, \(1/25=0.04\) and \(1/0.002=500\). Thus the standard errors along these two directions are \(0.2\) and \(\sqrt{500}=22.36\).

Transforming back to the two coordinates gives the approximate covariance matrix

\widehat{\operatorname{Cov}}(\hat{\theta}) = \begin{pmatrix} 250.02&-249.98\\ -249.98&250.02 \end{pmatrix}.

Each coordinate has an approximate standard error of \(\sqrt{250.02}\approx15.81\), and the large negative covariance reflects their strong dependence. More importantly, the combination along \((1,-1)^T/\sqrt{2}\) is especially uncertain. If that weak curvature is stable under careful recalculation and a likelihood profile is broad in the same direction, the result supports genuine weak information. But if a small numerical change moves the eigenvalue from 0.002 to 0.001, the estimated variance in that direction doubles from 500 to 1000. That sensitivity is a warning that numerical instability is affecting the precision calculation as well. The Hessian alone does not settle which explanation dominates.

Convergence and Covariance Conventions

An optimizer's convergence report and a covariance calculation answer different questions. Convergence concerns whether the algorithm has reached a point meeting its stopping criteria. Covariance estimation concerns the curvature near that point. A small gradient is useful evidence that the fit is near a stationary point, but it does not establish that the point is a suitable maximum, that the Hessian is well estimated, or that its inverse is reliable. For a regular interior maximum, the negative log-likelihood Hessian should be positive definite. A negative eigenvalue can indicate a saddle point, a failed fit, or a numerical Hessian problem; it is not a standard error to interpret.

Worked Example: Checking the Objective-Function Convention

A hypothetical program reports that it minimized \(Q=-2\ell\), with Hessian \(\nabla^2 Q(\hat{\theta})=\operatorname{diag}(50,8)\). Assume an interior fit and an adequate local quadratic approximation. Since \(Q=-2\ell\), the covariance approximation is twice the inverse Hessian of \(Q\):

\widehat{\operatorname{Cov}}(\hat{\theta}) =2\begin{pmatrix}1/50&0\\0&1/8\end{pmatrix} =\begin{pmatrix}0.04&0\\0&0.25\end{pmatrix}.

The standard errors are \(\sqrt{0.04}=0.2\) and \(\sqrt{0.25}=0.5\), in the respective parameter units. If someone instead used the inverse Hessian without the factor of 2, the reported standard errors would be \(\sqrt{0.02}=0.1414\) and \(\sqrt{0.125}=0.3536\). Those values are too small for this stated objective convention. The arithmetic is stable; the error would be a mismatch between the covariance formula and the objective that was differentiated.

Practical Diagnostics and Alternatives

A useful diagnosis combines evidence about the fitted point, the matrix calculation, and the likelihood shape. No single diagnostic proves that an estimate is trustworthy. In particular, a successful optimizer status does not certify a good covariance estimate, and a large condition number does not by itself prove that a parameter is statistically unidentifiable.

Check the fitted point.
Review convergence status, gradient information, parameter bounds, and whether the reported estimate is an interior optimum. If the Hessian has negative curvature at a claimed maximum, investigate the fit before interpreting standard errors.

Check scaling and matrix stability.
Inspect the Hessian's eigenvalues and conditioning in appropriately scaled coordinates. Recalculate the Hessian using sensible finite-difference settings or higher numerical precision when available. Large changes in small eigenvalues or standard errors across reasonable settings suggest numerical sensitivity.

Check whether uncertainty persists in the likelihood.
Examine likelihood profiles for suspect parameters or combinations. A broad, shallow profile supports weak information; a profile that is well defined while the inverse Hessian varies sharply suggests a problem with the local curvature calculation or its inversion.

Compare appropriate uncertainty estimates.
Profile likelihood and bootstrap estimates can assess uncertainty with fewer reliance on a local quadratic approximation, although each has its own assumptions and computational limits. A sandwich covariance may address certain forms of model misspecification; it is not a repair for an inaccurate Hessian inversion.

A profile likelihood can be especially useful near boundaries or when the likelihood is asymmetric, cases in which a symmetric normal approximation may be misleading even if matrix inversion is numerically stable. Bootstrap estimates can reveal sampling variability and fitting failures across resamples, but they may themselves be unstable with sparse data or frequent failed fits. These alternatives complement, rather than automatically replace, Hessian-based estimates.

Common warning signs include a non-positive-definite Hessian, standard errors that change substantially with parameter scaling or numerical settings, extreme correlations, and failed or boundary-hitting bootstrap fits. A warning sign calls for investigation, not a reflexive decision to delete a parameter, report no standard error, or force the Hessian to be positive definite without explanation.

Key takeaway: First ask whether the curvature is being calculated and inverted reliably; then ask what that reliable curvature implies about information in the data. Rescaling can improve numerical conditioning without changing correctly transformed uncertainty, while a stable shallow likelihood direction remains a genuine limitation of the data and model.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

Why does rescaling parameters sometimes improve a Hessian's condition number without reducing statistical uncertainty?
A negative-log-likelihood Hessian has eigenvalues 10 and 0.01. What are the approximate standard errors along the corresponding unit eigenvectors, assuming the inverse-Hessian approximation applies?
A covariance estimate changes greatly when a finite-difference step is modestly adjusted. What does this suggest, and what should be checked next?
For an objective \(Q=-2\ell\), what factor relates the inverse Hessian of \(Q\) to the usual inverse-information covariance approximation?

Solutions. (1) Rescaling changes the numerical coordinates and therefore the matrix condition number. If the covariance is transformed back correctly, the uncertainty in the original parameters is unchanged. (2) The covariance eigenvalues are \(1/10=0.1\) and \(1/0.01=100\), so the standard errors are \(\sqrt{0.1}\approx0.3162\) and \(\sqrt{100}=10\), in the corresponding coordinate units. (3) The Hessian inversion may be numerically unstable or the curvature poorly resolved. Check convergence and gradients, scaled eigenvalues, and Hessian calculations under reasonable numerical settings; use likelihood profiles to assess whether a broad direction also reflects genuine weak information. (4) Multiply the inverse Hessian of \(Q\) by 2, because \(Q=-2\ell\).