Numerical Analysis

Language revisions

C++20

std::bit_cast gives defined access to the representation; std::midpoint and std::lerp avoid classical overflow and cancellation traps.

C++23

Extended floating-point types (float16_t, bfloat16_t, float128_t) and constexpr <cmath> in part.P1467R9

C++26

std::simd exposes the same rounding semantics across lanes; contracts can carry numerical preconditions.P1928

Definition

Numerical analysis is the study of the difference between the problem posed and the problem solved.

The floating-point number system

DefinitionBinary floating-point format

A binary floating-point format with precision $t$ and exponent range $[e_{\min}, e_{\max}]$ represents numbers

$$x = \pm\, m \times 2^{e - t + 1}, \qquad 2^{t-1} \le m < 2^{t}, \quad e_{\min} \le e \le e_{\max},$$

together with signed zeros, subnormal numbers with $e = e_{\min}$ and $m < 2^{t-1}$, two infinities, and quiet and signalling NaNs.

The formats in use

Format	$t$	$u = 2^{-t}$	Decimal digits	Range
binary16	11	$4.9 \times 10^{-4}$	about 3	$10^{\pm 5}$
bfloat16	8	$3.9 \times 10^{-3}$	about 2	$10^{\pm 38}$
binary32	24	$6.0 \times 10^{-8}$	about 7	$10^{\pm 38}$
binary64	53	$1.1 \times 10^{-16}$	about 16	$10^{\pm 308}$
binary128	113	$9.6 \times 10^{-35}$	about 34	$10^{\pm 4932}$

bfloat16 is the instructive entry: it has the range of binary32 and the precision of nothing much, which is exactly the trade that machine-learning workloads want and that numerical linear algebra usually does not.

Spacing and the meaning of the unit roundoff

The spacing between consecutive representable numbers near $x$ is about $2u|x|$, so the relative spacing is constant across the range and the absolute spacing is not. This is why relative tolerances are the meaningful ones, and why an absolute tolerance appropriate near $1$ is either useless or unattainable near $10^{6}$ or $10^{-6}$.

Subnormals, infinities, and NaN

Subnormals fill the gap between zero and the smallest normal number, which preserves the property $x - y = 0 \iff x = y$. Without them, two distinct numbers could subtract to zero, and every algorithm that tests a difference against zero would be wrong.

Infinities and NaN make the arithmetic total, so that a computation produces a value rather than a trap. The cost is that NaN is unordered with everything including itself, so equality is not reflexive on the full carrier, and every ordering-based algorithm has undefined behaviour when NaN is present.

The standard model

DefinitionThe standard model

IEEE-754 arithmetic with round-to-nearest satisfies, for $\circ \in \{+, -, \times, \div\}$ and representable $x, y$ with no overflow or underflow,

$$\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta), \qquad |\delta| \le u,$$

where $u$ is the unit roundoff. Square root satisfies the same bound, and the fused multiply-add satisfies $\mathrm{fl}(ab + c) = (ab + c)(1 + \delta)$ with a single such $\delta$.

What the model does and does not say

It is a relative error bound, so it says nothing about quantities that are near zero by cancellation: an absolute error that is negligible relative to the operands can be enormous relative to their difference.

It holds for each operation independently, so an error analysis is an accumulation of $n$ such bounds, conventionally collected into

$$\prod_{i=1}^{n}(1 + \delta_i)^{\pm 1} = 1 + \theta_n, \qquad |\theta_n| \le \gamma_n = \frac{nu}{1 - nu}.$$

It excludes overflow and underflow, so an analysis based on it is silent about the circumstances in which a computation fails catastrophically rather than inaccurately.

Rounding modes

Round-to-nearest with ties to even is the default and the one every bound in this reference assumes. The directed modes exist for interval arithmetic, where the lower and upper bounds must be rounded outward, and changing the mode globally invalidates the error-free transformations below, which is a good reason to change it only in a scoped and documented way.

Conditioning and stability

DefinitionConditioning and stability

For a problem $f$ and input $x$, the relative condition number is

$$\kappa(f, x) = \lim_{\varepsilon \to 0} \sup_{\lVert \delta x \rVert \le \varepsilon \lVert x \rVert} \frac{\lVert f(x + \delta x) - f(x) \rVert / \lVert f(x) \rVert}{\lVert \delta x \rVert / \lVert x \rVert}.$$

An algorithm $\hat{f}$ is backward stable when for every $x$ there is a perturbed input with $\hat{f}(x) = f(x + \delta x)$ and $\lVert \delta x \rVert = O(u)\lVert x \rVert$.

The governing inequality

$$\frac{\lVert \hat{f}(x) - f(x) \rVert}{\lVert f(x) \rVert} \lesssim \kappa(f,x)\, u.$$

Accuracy is the product of a property of the problem and a property of the algorithm. A backward stable algorithm on an ill-conditioned problem returns an inaccurate answer, and no implementation effort changes this.

Why backward error is the right target

Forward error is what a user wants and is usually not attainable, because it is bounded below by $\kappa u$ whatever the algorithm. Backward error is attainable, is checkable by computing a residual, and composes: a sequence of backward stable steps is backward stable with the constants added. Designing for backward stability and reporting conditioning separately is therefore the only arrangement in which both quantities are available.

Mixed and forward stability

Some algorithms are neither backward stable nor useless. An algorithm is mixed stable when it computes nearly the right answer to nearly the right problem, and forward stable when its forward error is comparable to that of a backward stable method without itself being backward stable. Fast matrix multiplication is the standard example of the latter: Strassen is norm-wise forward stable and not componentwise backward stable.

Condition numbers of the standard problems

Problem	Condition number
$x + y$	$(\lvert x\rvert + \lvert y\rvert)/\lvert x + y\rvert$
$xy$, $x/y$	$1$
$\sqrt{x}$	$1/2$
$\exp(x)$	$\lvert x \rvert$
$\log(x)$	$1/\lvert \log x \rvert$
Root of $f$	$1/\lvert f'(x^{*}) \rvert$
$Ax = b$	$\kappa(A)$
Eigenvalue of symmetric $A$	$1$
Eigenvalue of general $A$	$1/\lvert y^{*}x \rvert$ for unit left and right vectors

Two rows carry most of the practical weight. Addition is the only elementary operation that can be ill conditioned, which is why cancellation dominates error analysis. And symmetric eigenvalues are perfectly conditioned while general ones are not, which is why the symmetric case admits far stronger guarantees.

Motivation

Every performance argument in this reference is conditional on producing the right answer, and floating-point arithmetic makes "the right answer" a quantitative rather than a binary notion.

Which reformulations are legitimate

Reassociating a sum, using a fused multiply-add, replacing division by reciprocal multiplication, and vectorizing a reduction all change the computed result. Each is acceptable when its effect on the error bound is understood, and unacceptable when it is not. The essay Writing the Reassociation develops the point at length.

Which precision suffices

Mixed-precision methods, computing a factorization in binary32 and refining in binary64, halve memory traffic and double vector throughput. Whether they deliver the accuracy of the higher precision depends on the conditioning, and iterative refinement supplies the criterion: refinement converges while $\kappa(A)\, u_{\text{low}} < 1$.

When a result is meaningless

A computed residual smaller than $u \lVert A \rVert \lVert x \rVert$ carries no information, and an iteration driven past that point is fitting noise. Recognising the floor is part of using a method, not an optional refinement.

Why the worst case is still worth computing

Worst-case bounds are pessimistic by roughly $\sqrt{n}$, since rounding errors behave statistically. They are computed anyway because they establish the scaling: a bound of $\gamma_n$ against $\gamma_{\log n}$ is the difference between an algorithm that degrades with size and one that does not, and that distinction survives the pessimism.

Formal formulation

Summation

Summation is the canonical example, and the one whose bounds are quoted throughout this reference. For the sequential sum of $n$ terms, the computed result satisfies

$$\hat{s}_n = \sum_{i=1}^{n} x_i (1 + \theta_i), \qquad |\theta_i| \le \gamma_{n},$$

giving the forward bound

$$\left| \hat{s}_n - \sum_i x_i \right| \le \gamma_{n-1} \sum_i |x_i|.$$

The condition number of summation

The ratio $\sum |x_i| / |\sum x_i|$ is the condition number of summation. For terms of like sign it equals one and the sum is computed to full accuracy; for cancelling terms it is unbounded, and no summation algorithm can produce an accurate result from inaccurate inputs.

This is the cleanest instance of the general principle: the accuracy attainable is set by the problem, and the algorithm determines only whether that accuracy is achieved.

The methods compared

Method	Error bound	Cost
Sequential	$\gamma_{n-1}\sum \lvert x_i \rvert$	$n$ additions
Increasing order	$\gamma_{n-1}\sum \lvert x_i \rvert$, smaller in practice	$n \log n$ for the sort
Pairwise	$\gamma_{\lceil \log_2 n\rceil}\sum \lvert x_i \rvert$	$n$ additions
Kahan or Neumaier	$(2u + O(nu^2))\sum \lvert x_i \rvert$	about $4n$
Exact accumulator	correctly rounded	$n$ plus a wide accumulator

Pairwise summation is the default in this reference: its accuracy is within a factor of $\log n / n$ of compensated summation at the cost of the naive form, and its recursive structure vectorises and parallelises without further change.

Error-free transformations

Certain rounding errors are themselves representable, which allows them to be recovered exactly.

PropositionTwo-sum

For representable $a, b$ with $s = \mathrm{fl}(a + b)$, the error $e = (a + b) - s$ is representable, and Knuth's six-operation sequence computes it exactly with no branch.

PropositionTwo-product

For representable $a, b$ with $p = \mathrm{fl}(ab)$, the error $e = ab - p$ is representable and equals $\mathrm{fma}(a, b, -p)$ exactly.

What they are used for

These are the primitives on which compensated summation, compensated dot products and double-double arithmetic are built. A double-double number is an unevaluated sum of two floating-point numbers whose exponents do not overlap, giving about 32 decimal digits at roughly twenty times the cost of binary64 arithmetic and without hardware support.

The fragility

They are correct only under round-to-nearest and are destroyed by compiler options permitting reassociation, which is the strongest single argument against enabling such options globally. A build that turns on fast-math does not make these routines slower; it makes them return zero.

Cancellation

Cancellation is a property of the problem, not of the arithmetic: subtracting nearly equal quantities is exact, and what it reveals is error already present in the operands.

PropositionSterbenz

If $y/2 \le x \le 2y$ for representable $x, y > 0$, then $x - y$ is computed exactly.

Standard reformulations

The remedies are reformulations that avoid forming the cancelling difference:

Unstable form	Stable form
$x^2 - y^2$	$(x-y)(x+y)$
$\sqrt{x+1} - \sqrt{x}$	$1/(\sqrt{x+1}+\sqrt{x})$
$1 - \cos x$	$2\sin^2(x/2)$
$\log(1+x)$ for small $x$	std::log1p(x)
$e^{x} - 1$ for small $x$	std::expm1(x)
Smaller quadratic root	$x_2 = c/(a x_1)$ from the larger root
Variance by $\sum x^2 - n\bar{x}^2$	Welford's recurrence

Each replaces an unstable expression with an algebraically equal, numerically stable one. The presence of log1p and expm1 in the standard library is the language acknowledging this list.

Linear systems

For Gaussian elimination with partial pivoting the computed solution $\hat{x}$ satisfies $(A + \Delta A)\hat{x} = b$ with

$$\lVert \Delta A \rVert_\infty \le \gamma_{3n} \rho_n \lVert A \rVert_\infty,$$

where $\rho_n$ is the growth factor.

The growth factor

The bound on $\rho_n$ is $2^{n-1}$, attained only by contrived matrices; in practice $\rho_n$ grows like $n^{1/2}$, which is why the method is used despite its exponential worst case. The gap between the worst case and the observed behaviour is one of the oldest open questions in the subject, and it is the reason partial pivoting is universal while complete pivoting, whose bound is $\Theta(n^{1/2}\log n)$, is not.

Iterative refinement

PropositionRefinement in mixed precision

Let $\hat{x}$ solve $Ax = b$ using a factorization computed at unit roundoff $u_{\text{low}}$, and refine by computing $r = b - A\hat{x}$ at precision $u_r$, solving $Ad = r$ with the same factorization, and setting $\hat{x} \leftarrow \hat{x} + d$. The iteration converges while $\kappa(A)\,u_{\text{low}} < 1$, and the limiting forward error is $O(\kappa(A) u_r)$.

The proposition is what makes low-precision factorization respectable: the factorization is done once at low cost, and the accuracy of the working precision is recovered by $\Theta(n^2)$ work per refinement step against $\Theta(n^3)$ for the factorization. Computing the residual in a precision higher than the working one extends the range of $\kappa$ for which the method applies.

C++ representation

Numerical properties in the interface

#include <algorithm>
#include <cmath>
#include <concepts>
#include <limits>

namespace hpc {

template <std::floating_point T>
inline constexpr T unit_roundoff = std::numeric_limits<T>::epsilon() / T{2};

// Relative-plus-absolute comparison: an absolute tolerance alone is
// meaningless without a scale, a relative one alone fails near zero.
template <std::floating_point T>
constexpr bool approximately_equal(T x, T y, T relative, T absolute) noexcept {
    const T diff = std::abs(x - y);
    return diff <= absolute || diff <= relative * std::max(std::abs(x), std::abs(y));
}

// Sterbenz: exact when the operands are within a factor of two.
template <std::floating_point T>
constexpr bool subtraction_is_exact(T x, T y) noexcept {
    return y / T{2} <= x && x <= T{2} * y;
}

}  // namespace hpc

Comparing by units in the last place

For tests that must be tight, a tolerance in units in the last place is more informative than a relative one, because it counts representable numbers rather than a fraction:

#include <bit>
#include <cmath>
#include <cstdint>

namespace hpc {

// Distance in representable numbers. Requires both operands finite and of the
// same sign; the mapping from bit pattern to ordinal is monotone there.
inline std::uint64_t ulp_distance(double a, double b) noexcept {
    const auto ia = std::bit_cast<std::int64_t>(a);
    const auto ib = std::bit_cast<std::int64_t>(b);
    return static_cast<std::uint64_t>(ia > ib ? ia - ib : ib - ia);
}

}  // namespace hpc

A routine documented as accurate to 4 ulp makes a stronger and more checkable claim than one documented as accurate to $10^{-15}$, because the second depends on the magnitude of the result.

The error-free transformations

#include <cmath>

namespace hpc {

struct two_sum_result { double sum; double error; };

// Knuth's two-sum: sum + error is the exact sum of a and b.
// Correct only under round-to-nearest with no reassociation.
constexpr two_sum_result two_sum(double a, double b) noexcept {
    const double s = a + b;
    const double bb = s - a;
    const double err = (a - (s - bb)) + (b - bb);
    return {s, err};
}

struct two_product_result { double product; double error; };

// Exact with a fused multiply-add: product + error == a * b exactly.
inline two_product_result two_product(double a, double b) noexcept {
    const double p = a * b;
    return {p, std::fma(a, b, -p)};
}

}  // namespace hpc

Building these requires that the compiler be prohibited from simplifying (a - (s - bb)) to zero. The requirement belongs in the build configuration and should be asserted there, since a silent change to the flags removes the correctness of every routine built on them.

Stating what a routine guarantees

#include <cstddef>
#include <span>

namespace hpc {

// Computes the sum of xs.
//
// Guarantees: the computed s satisfies
//   |s - sum(xs)| <= gamma_{ceil(log2 n)} * sum(|xs|),
// with gamma_k = k u / (1 - k u), u the unit roundoff.
//
// The bound is independent of the order of the input, unlike the sequential
// sum, whose bound has n in place of log2 n.
double pairwise_sum(std::span<const double> xs) noexcept;

}  // namespace hpc

A numerical routine whose documentation states what it returns but not the accuracy with which it returns it has specified nothing testable. Every numerical entry in this reference carries a bound of this shape.

Implementation

The three summations

#include <cstddef>
#include <span>

namespace hpc {

// Error <= gamma_{n-1} * sum |x_i|.
constexpr double naive_sum(std::span<const double> xs) noexcept {
    double s = 0.0;
    for (const double x : xs) { s += x; }
    return s;
}

// Error <= gamma_{ceil(log2 n)} * sum |x_i|. The cutoff is where the call
// overhead exceeds the accuracy benefit; it is measured, not guessed.
constexpr double pairwise_sum(std::span<const double> xs) noexcept {
    constexpr std::size_t cutoff = 128;
    if (xs.size() <= cutoff) { return naive_sum(xs); }
    const std::size_t half = xs.size() / 2;
    return pairwise_sum(xs.first(half)) + pairwise_sum(xs.subspan(half));
}

// Error about 2u * sum |x_i|, at roughly four times the arithmetic.
constexpr double kahan_sum(std::span<const double> xs) noexcept {
    double s = 0.0;
    double c = 0.0;
    for (const double x : xs) {
        const double y = x - c;
        const double t = s + y;
        c = (t - s) - y;
        s = t;
    }
    return s;
}

}  // namespace hpc

A compensated dot product

#include <cmath>
#include <cstddef>
#include <span>

namespace hpc {

// Ogita, Rump and Oishi: the result is as accurate as a dot product computed
// in twice the working precision and rounded once, at about five times the
// arithmetic of the naive form.
inline double compensated_dot(std::span<const double> x,
                              std::span<const double> y) noexcept
{
    double s = 0.0;
    double c = 0.0;
    for (std::size_t i = 0; i < x.size(); ++i) {
        const double p = x[i] * y[i];
        const double e = std::fma(x[i], y[i], -p);   // exact product error
        const double t = s + p;
        const double b = t - s;
        c += ((s - (t - b)) + (p - b)) + e;          // exact sum error
        s = t;
    }
    return s + c;
}

}  // namespace hpc

On a bandwidth-bound kernel the extra arithmetic hides under the loads, so the accuracy is often obtained at almost no cost. This is the case worth checking for before accepting a poor bound.

Iterative refinement

#include <cstddef>
#include <span>

namespace hpc {

// Given a factorization computed in low precision, refine the solution in
// working precision:
//   r <- b - A x        (in working or higher precision)
//   solve M d = r       (using the low-precision factorization)
//   x <- x + d
//
// Converges while kappa(A) * u_low < 1, at Theta(n^2) per step against
// Theta(n^3) for the factorization.
struct refinement_result {
    std::size_t steps;
    double residual_norm;
    bool converged;
};

}  // namespace hpc

Testing numerical code

Exact comparison is not available, and a tolerance chosen to make a test pass is a defect in the test. Three defensible practices:

Compare against a higher-precision reference on small inputs, with a tolerance derived from the stated bound rather than from observation.

Check the backward error on large inputs: form the residual, compare it against the bound the specification claims, and fail if it exceeds it. This is universally applicable and requires no reference solution.

Test on inputs designed to be ill conditioned, with the expected loss of accuracy stated in advance. A routine that passes only on well-conditioned data has been tested on the cases where every implementation succeeds.

Complexity

Complexity; n terms, binary64

Operation	Time	Space	Remarks
Naive summation	n flops	Θ(1)	Error ≤ γ_{n−1} Σ|xᵢ|
Pairwise summation	n flops	Θ(log n) stack	Error ≤ γ_{⌈log₂ n⌉} Σ|xᵢ|
Kahan or Neumaier summation	about 4n flops	Θ(1)	Error about 2u Σ|xᵢ|
Exact accumulator summation	Θ(n) with a wide accumulator	Θ(1) fixed	Correctly rounded result
two_sum	6 flops	Θ(1)	Exact; no branch
two_product with fma	2 flops	Θ(1)	Exact
Compensated dot product	about 10n flops	Θ(1)	As if computed in twice the precision
Double-double arithmetic	10 to 25× binary64	2× storage	About 32 decimal digits
Iterative refinement, one step	Θ(n²)	Θ(n)	Recovers accuracy while κu < 1
Scaled norm against naive	about 2× slower	Θ(1)	Avoids spurious overflow

Performance considerations

Denormals

Denormal operands are handled by microcode on many implementations, at a penalty of tens to hundreds of cycles per operation. A kernel that decays into the denormal range, such as a damped oscillation or a decaying exponential, can slow by an order of magnitude with no change in instruction count.

Flush-to-zero removes the penalty and changes the arithmetic, including the property that $x - y = 0$ implies $x = y$. The choice is made deliberately and recorded, not left to ambient compiler flags.

Division and square root

Division and square root have latencies of ten to twenty cycles and are not pipelined at full rate. Replacing $n$ divisions by one division and $n$ multiplications is a standard transformation that changes the error bound by one additional rounding per element: usually acceptable, always to be stated.

Reciprocal-estimate instructions with Newton refinement are faster still and less accurate; they belong behind an explicit opt-in rather than behind a global flag.

Fused multiply-add

The fused operation improves both speed and accuracy, computing $ab + c$ with a single rounding. It also changes results relative to separate operations, so a routine whose correctness depends on the absence of fusion, such as two_sum, and one that depends on its presence, such as two_product, must be compiled with contraction controlled explicitly.

Reduced precision as a bandwidth optimisation

Binary32 halves memory traffic and doubles vector lane count relative to binary64. Whether the accuracy suffices is answered by $\kappa u < 1$ in the working precision, and where it does not, iterative refinement recovers the accuracy at $\Theta(n^2)$ cost on top of a $\Theta(n^3)$ factorization. For bandwidth-bound kernels the speed-up is close to the precision ratio; for compute-bound kernels it is the vector-width ratio.

Compensation is nearly free when memory-bound

Compensated summation performs about four times the arithmetic of the naive form. On a bandwidth-bound reduction the extra operations hide entirely under the loads, so the accuracy is obtained at almost no cost. This is the standard case where the better algorithm is also the cheaper one.

Reproducibility costs something

Bitwise-reproducible reductions require a fixed reduction tree independent of thread count and vector width, which forfeits some of the scheduling freedom that makes reductions fast. The cost is typically 10 to 30 percent, and it is worth paying only where reproducibility is a requirement rather than a preference.

Limitations

Backward stability is a statement about the algorithm and says nothing about accuracy for an ill-conditioned problem. Reporting a small residual as evidence of a good solution is the most common error in numerical software: the residual bounds the backward error, not the forward one.

Error bounds are worst-case and typically pessimistic by a factor of about $\sqrt{n}$, since rounding errors behave statistically. They establish scaling, not magnitudes.

The standard model excludes overflow and underflow, so an analysis based on it is silent about the circumstances in which a computation fails catastrophically. Scaling arguments must be made separately, and they are why production norm routines are written with an explicit scaling pass rather than as a plain sum of squares.

Comparison of floating-point results has no tolerance-free formulation. Every equality test in numerical code is a decision about scale, and the decision is part of the specification of the routine.

Condition numbers are themselves computed quantities, usually estimated rather than computed, and an estimate that under-reports by a factor of two changes the number of digits one is entitled to trust by less than one, which is why the practice is acceptable.

Related entries

Analysis for the truncation side; Algebra for the algebraic laws that rounding breaks; Linear Algebra for conditioning of linear problems; Numerical Linear Algebra for the factorization bounds; Vectorization for the reassociation questions.

References

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002.
• IEEE Std 754-2019, IEEE Standard for Floating-Point Arithmetic.
• D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys 23(1), 1991.
• J.-M. Muller et al., Handbook of Floating-Point Arithmetic, 2nd ed., Birkhäuser, 2018.
• T. Ogita, S. Rump and S. Oishi, Accurate Sum and Dot Product, SIAM J. Sci. Comput. 26(6), 2005.
• E. Carson and N. J. Higham, Accelerating the Solution of Linear Systems by Iterative Refinement in Three Precisions, SIAM J. Sci. Comput. 40(2), 2018.
• W. Kahan, Further Remarks on Reducing Truncation Errors, CACM 8(1), 1965.