Stochastic rounding: implementation, error analysis and applications

Abstract

Stochastic rounding (SR) randomly maps a real number x to one of the two nearest values in a finite precision number system. The probability of choosing either of these two numbers is 1 minus their relative distance to x. This rounding mode was first proposed for use in computer arithmetic in the 1950s and it is currently experiencing a resurgence of interest. If used to compute the inner product of two vectors of length n in floating-point arithmetic, it yields an error bound with constant $\sqrt{n}u$ with high probability, where u is the unit round-off. This is not necessarily the case for round to nearest (RN), for which the worst-case error bound has constant nu. A particular attraction of SR is that, unlike RN, it is immune to the phenomenon of stagnation, whereby a sequence of tiny updates to a relatively large quantity is lost. We survey SR by discussing its mathematical properties and probabilistic error analysis, its implementation, and its use in applications, with a focus on machine learning and the numerical solution of differential equations.

Keywords: IEEE 754; bfloat16; binary16; floating-point arithmetic; machine learning; rounding error analysis.

Figure 1.

Stochastic rounding rounds the real number x to the next smaller or the next larger value in F, which we denote by $\lfloor x\rfloor$ and $\lceil x\rceil$, respectively. In the example on the left, x is one quarter of the way between $\lfloor x\rfloor$ and $\lceil x\rceil$, thus RN will round x to $\lfloor x\rfloor$, while mode 1 SR will round it to $\lceil x\rceil$ with probability q(x) = 0.25 and to $\lfloor x\rfloor$ with probability 1 − q(x) = 0.75. In the example on the right, w is three-quarters of the way between $\lfloor w\rfloor$ and $\lceil w\rceil$, thus RN will round w to $\lceil w\rceil$, while mode 1 SR will round it to $\lceil w\rceil$ with probability q(w) = 0.75 and to $\lfloor w\rfloor$ with probability 1 − q(w) = 0.25.

Figure 2.

Alignment of bits in algorithms for stochastic rounding based on sums. The random bits are added to the significand m_t followed by the truncation of it. How the bits are generated and added depends on the implementation—we may only add the k bits to the top k bits of the bottom part of the significand and then use a carry-out bit to control the rounding of m_r after the truncation, or we may pack the k random bits in a word as long as m_t and add it to m_t using integer arithmetic: the propagating carry will cause rounding in the top p bits.

Figure 3.

Example: rounding a binary32 number y to a format with p = 11 significant digits (including the implicit bit) using the SR algorithm implemented in QPyTorch. The number n is the integer sum of the two bit strings representing y and m; fl(y) is obtained by zeroing out the trailing 24 − p = 13 trailing bits of n. For floating-point numbers, the three binary strings represent the sign (s), the unbiased exponent (e) and the integer significand without the implicit bit (m); the last group is further divided into a group of p − 1 bits to keep and 24 − p bits to zero out. We also report the hexadecimal string representing the numbers, and for floating-point numbers the corresponding exact decimal representations.

Figure 4.

Relative errors for computing $\sum _{i=1}^{n}1/i$ with RN and SR. The densely dashed and dash-dotted lines are the worst-case error bound for RN and the probabilistic error bound for SR (with λ = 1), respectively. Stochastic rounding experiments are repeated 10 times; the solid line represents the average error, the edges of the shaded area the minimum and maximum error. (a) binary16 arithmetic, (b) bfloat16 arithmetic.

Figure 5.

Backward error for computing y = Ax with RN and SR, where $A\in {\mathbb{R}}^{100\times n}$ has entries drawn from the uniform distribution over [0, 10⁻³] and $x\in {\mathbb{R}}^n$ has entries sampled from the uniform distribution over [0, 1]. The densely dashed and dash-dotted lines are the worst-case error bound for RN and the probabilistic error bound for SR (with λ = 1), respectively. Stochastic rounding experiments are repeated 10 times; the solid line represents the average error, and the edges of the shaded area the minimum and maximum error. (a) binary16 arithmetic, (b) bfloat16 arithmetic.

Figure 6.

Absolute errors in the forward Euler method for an ODE with exponentially decaying solutions with different floating-point arithmetics and rounding modes. Stochastic rounding experiments are repeated 10 times; the solid line represents the average error, the edges of the shaded area the minimum and maximum error. The step size is the interval length divided by n. The experiment is adapted from [37]. (a) y′ = −y, y(0) = 2⁻⁶, over [0, 1], (b) $y^{\prime }=-y/20,y(0)=1\hspace{0.17em}\hspace{0.17em}\mathrm{over}\hspace{0.17em}[0,2^{-6}].$

Figure 7.

Computed solutions from the forward Euler method with SR and RN for the ODE system (8.4). The exact solution is the unit circle. The x- and y-axis represent u and v, respectively. Note that in (d) and (h) only a small part of the solution computed with round-to-nearest is visible (marked with an arrow) since the ODE solver failed because of stagnation. The experiment is adapted from [37]. Stochastic rounding experiments are repeated 10 times; the solid line represents the average trajectory, the edges of the shaded area the points that are farthest from the exact solution in the Euclidean distance. (a) n = 2⁵, (b) n = 2⁹, (c) n = 2¹¹, (d) n = 2¹³, (e) n = 2⁵, (f) n = 2⁹, (g) n = 2¹⁴, (h) n = 2¹⁶.

Figure 8.

Comparison between the numerical steady-state solutions obtained with RN and SR with forward Euler and the bfloat16 format for different initial conditions. All SR solutions essentially converge to the same steady state. On the other hand, when RN is used different initial conditions lead to different steady-state solutions. The noise term in the initial condition has been obtained by sampling independent standard Gaussian random variables at each mesh node.

Figure 9.

Plot of the relative global rounding error in the L² norm for the solution of (8.5) in 1D (left) and 2D (right) with forward and backward Euler (FE and BE, respectively) and the bfloat16 format. We circled the RN data points for which the solution stagnates at the initial condition. The error behaviour matches the theoretical predictions from [43]. The SR errors are average errors computed within 2 digits of accuracy as in [43]. The worst-case SR errors were only a small constant factor larger than the average and are thus omitted.