Numerical Computing Fundamentals

Precision can be lost during a computation when predictors are highly correlated, when fitting a high-degree polynomial, or when the data has a large mean and small variance. This page explains how that happens and what MIDAS does about it. It also covers the concepts behind the Numerical Accuracy page.

Floating-Point Numbers and Significant Digits

MIDAS approximates real numbers with IEEE 754 double-precision floating-point numbers. This format represents a real number as a binary significand paired with an exponent, storing 52 bits of the significand. A nonzero value is stored with its exponent adjusted so that the leading digit of the significand is 1. Since that digit is always 1, it does not need to be stored and is implied. The 52 stored bits plus the implicit 1 bit give 53 significant bits, corresponding to about 15.9 decimal significant digits:

\log_{10}(2^{53}) \approx 15.95

This is the upper bound on the significant digits a single double-precision number can represent.

The accuracy checks on the Numerical Accuracy page use data with known correct values and report how many significant digits MIDAS's computed value shares with the correct value, as an index called LRE (Log Relative Error). The closer the LRE is to this upper bound, the fewer significant digits were lost during the computation. That page covers how to judge LRE values.

When converting a real number to a floating-point number, any fraction that cannot be represented in a finite number of bits is rounded. This rounding error is tiny for a single operation, but can accumulate over successive computations. How much it accumulates depends on the algorithm and the data. The following sections explain catastrophic cancellation and the condition number, the two main mechanisms that amplify the accumulation.

Catastrophic Cancellation

Catastrophic cancellation occurs when subtracting two nearly equal floating-point numbers, causing a large loss of significant digits.

Consider $a = 1.234567890123\mathbf{45}$ and $b = 1.234567890123\mathbf{00}$ . Both have 15 significant digits, but their difference $a - b = 0.00000000000045$ has only 2. The leading zeros merely indicate position and are not significant, so only the trailing 45 carries meaningful information. Once these decimals are stored in double precision their trailing digits already carry rounding error, so the 13 digits shared between $a$ and $b$ cancel out, leaving only those rounding-error-laden low digits.

In statistics, variance computation is a classic example. The algebraic form $(\sum x_i^2 - n\bar{x}^2) / (n-1)$ concentrates the loss into a single subtraction: when the mean is large and the variance small, $\sum x_i^2$ and $n\bar{x}^2$ are nearly equal, and most of the significant digits vanish at once. The definitional form $\sum(x_i - \bar{x})^2 / (n-1)$ also suffers cancellation in each deviation $x_i - \bar{x}$ , but the squared deviations are summed without any subtraction, so the loss stays limited to the digits lost in each individual deviation, and most digits are never lost at once as in the algebraic form.

MIDAS computes variance with Welford's online algorithm (Welford, 1962). Updating the mean and the sum of squared deviations one data point at a time avoids the subtraction of large intermediate values seen in the algebraic form. In addition, the values are shifted by the first data point before the incremental update, which curbs the precision loss that a large offset would otherwise cause in the running mean. This achieves precision comparable to the definitional form in a single pass over the data.

What Welford's algorithm prevents is the single subtraction of the algebraic form and the degradation of the running mean. The cancellation in the individual deviations remains, just as in the definitional form, so when the mean is extremely large and the variance extremely small, the standard deviation loses precision. NumAcc3 and NumAcc4 on the Numerical Accuracy page show how much.

Design Matrix

In the linear regression model $Y = X\beta + \varepsilon$ , $X$ is the design matrix. $Y$ is the response vector, $\beta$ is the coefficient vector, and $\varepsilon$ is the error term. The design matrix has $n$ rows and $p$ columns, with each row corresponding to one observation. When an intercept is included, one of the $p$ columns is a constant column of ones, and the remaining columns correspond to the predictors ( $p$ is the total column count including the intercept). See OLS Fundamentals for the formulation and the derivation of the OLS (ordinary least squares) estimator.

The coefficients $\hat\beta$ are computed from the design matrix, so their numerical accuracy depends on its properties. The more strongly the predictors are correlated, the larger the condition number explained in the next section becomes, and the more the rounding errors are amplified. Before fitting, you can check the correlations among predictors in the Relationships section of the Statistics tab; after fitting, the VIF in the Linear Regression and GLM coefficient tables measures multicollinearity.

Condition Number

The condition number measures how much small perturbations in the input are amplified in the output. For a nonsingular matrix $A$ , the condition number is defined as:

\kappa(A) = \|A\| \cdot \|A^{-1}\|

$\|A\|$ is the matrix 2-norm (largest singular value). Since $\|A\| \cdot \|A^{-1}\| \ge \|AA^{-1}\| = \|I\| = 1$ , the condition number has a lower bound of 1. A matrix with a condition number close to 1 is called well-conditioned; one with a large condition number is called ill-conditioned. There is no sharp threshold; the distinction depends on the precision required.

A condition number of $\kappa$ causes, in the worst case, a loss on the order of $\log_{10}(\kappa)$ significant digits. Subtracting from double precision's roughly 15.9 digits gives an approximate estimate:

\text{significant digits} \approx 15.9 - \log_{10}(\kappa)

This is a worst-case order-of-magnitude estimate, not a strict guarantee. Actual precision depends on the data and algorithm details, and can substantially exceed the estimate. When computing $\hat\beta$ in linear regression, the relevant condition number depends on the solver: QR decomposition works directly with $X$ , so $\kappa(X)$ applies; the normal equations solve $(X'X)^{-1}X'Y$ , so $\kappa(X'X) = \kappa(X)^2$ applies. MIDAS uses QR decomposition for coefficient estimation, so the dominant factor in precision is $\kappa(X)$ .

When fitting a model via QR decomposition, MIDAS computes an estimate of the condition number and shows a warning with the results when the estimate exceeds $10^{10}$ . This estimate is the Frobenius condition number $\|R\|_F \cdot \|R^{-1}\|_F$ of the upper triangular factor $R$ , which satisfies $\kappa(X) \le \|R\|_F \cdot \|R^{-1}\|_F \le p\,\kappa(X)$ , where $p$ is the number of columns of the design matrix including the intercept. Because it never underestimates the condition number, an ill-conditioned design matrix is not missed. It overestimates by at most a factor of $p$ , a difference of $\log_{10} p$ digits in the significant-digit estimate.

For Linear Regression, two-way ANOVA, and DoE, the problem being solved is unweighted least squares, so this estimate corresponds to $\kappa(X)$ . GLMs are fit with IRLS (iteratively reweighted least squares), so the condition number is estimated from the weighted design matrix $\sqrt{W}\,X$ at each iteration, and the maximum across iterations drives the warning ( $W$ is the per-iteration weight matrix). With the Gaussian family and the identity link the weights are constant and the estimate coincides with $\kappa(X)$ , but for logistic or Poisson regression the weights vary substantially across observations, so the column correlations of $X$ alone may not explain the warning.

Even when the warning appears, the coefficients are still computed and shown; interpret the coefficients in light of the lost significant digits.

Ill-conditioning arises from two distinct factors: high correlation among predictors, and the magnitude of predictor values. Correlation can be checked with the Relationships section and the VIF described in the Design Matrix section.

The magnitude of predictor values raises the condition number in two ways. A predictor whose mean is far from zero becomes nearly collinear with the intercept column. And when the variances of the predictors differ greatly from each other, that difference alone drives up the condition number.

Which factor is the cause can be told apart by standardizing each predictor to mean 0 and variance 1 and refitting. Standardization removes only the effect of magnitude, so if the warning still appears, the cause is correlation; if the warning goes away, the cause was magnitude.

When magnitude is the cause, the loss of precision is far smaller than the estimate from the condition number, and the original model can usually be used as is¹. The reason and a worked example are given in Polynomial Regression.

Polynomial Regression

In polynomial regression, the columns of the design matrix are $1, x, x^2, \ldots, x^d$ . As the degree $d$ increases, the power columns become strongly correlated and their magnitudes spread far apart, so both factors from the previous section act together and the design matrix becomes ill-conditioned.

The NIST StRD datasets on the Numerical Accuracy page show this effect. The simple regression dataset Norris achieves a coefficient LRE of 12.3, but the 10th degree polynomial Filip drops to 7.3.

The significant-digit estimate from the condition number is a worst-case order of magnitude, and the actual precision often exceeds it. Filip's design matrix has $\kappa(X) \approx 2 \times 10^{15}$ . Plugging that into the estimation formula predicts less than one significant digit, yet the observed LRE of 7.3 far exceeds it.

The estimate is pessimistic because the rounding errors of the QR decomposition MIDAS uses for coefficient estimation enter relative to the scale of each column of the design matrix. Scale differences between columns drive up the raw condition number, but they do not contribute correspondingly to the actual loss of precision, so the raw condition number overestimates the difficulty by however much the columns' scales differ. Scaling each column to unit $\ell_2$ norm before computing the condition number gives about $5 \times 10^9$ for Filip, so the estimate becomes about 6 digits, close to the observed value. Correlation-driven ill-conditioning, where the columns already have comparable scales, does not get this overestimate.

To reduce the condition number, replace the monomial basis $1, x, x^2, \ldots$ with an orthogonal polynomial basis. When the design matrix consists only of the intercept and orthogonal polynomial columns, the columns are mutually orthogonal with equal norms, so the condition number is theoretically 1. Floating-point rounding makes the computed value deviate slightly from 1. The Orthogonal Polynomials tab in MIDAS applies this transformation.

References

Welford, B. P. (1962). Note on a method for calculating corrected sums of squares and products. Technometrics, 4(3), 419-420. https://www.jstor.org/stable/1266577