Floating-Point Reduction Trees (July 2026), 1–20
RESEARCH PREPRINT
A Second-Moment Theory for Floating-Point Reduction Trees Piyush Sao , Narasinga Miniskar , Pedro Valero-Lara , Keita Teranishi , and Sudip Seal
arXiv:2607.18758v1 [math.NA] 21 Jul 2026
Oak Ridge National Laboratory, Oak Ridge, TN, USA Emails (author order): [email protected]; [email protected]; [email protected]; [email protected]; [email protected].
Abstract Summation error depends on partial-sum order, which standard worst-case bounds omit. To capture this dependence, we derive an exact mean-square error (MSE) recurrence for a binary tree T under conditionally unbiased rounding. With unit roundoff u, the constant-ν model (Section 2.2) sets the local variance at the pre-rounding value x to νu2 x2 . Its leading tree-dependent cost for the input vector p is p⊤ KT p, where the common-ancestor kernel KT (Section 2.3) counts internal ancestors shared by leaves i and j. For i.i.d. inputs of mean µ and variance τ2 , this expected cost is τ2 Λ1 ( T ) + µ2 Λ2 ( T ) ; the scalar tree statistics Λ1 and Λ2 are defined in Section 2.4. Here Λ1 is total leaf depth, while Λ2 sums squared internal-subtree sizes. Thus Λ1 governs centered inputs, while Λ2 captures nonzero means. We use these statistics to characterize optimal tree topologies and schedules. Balanced and sequential trees attain the centered extrema. For k inputs, optimal two-stage sequential blocking yields root-mean-square (RMS) error that scales as k3/4 . For fixed-stage hierarchies, geometric schedules are optimal for centered inputs, whereas the optimal noncentered stage exponents halve successively. For independent centered inputs with unequal variances, Huffman coding minimizes variance-weighted depth over free leaf assignments and unconstrained trees. We extend the kernel to matrix multiplication through operand Gram matrices. To evaluate these predictions, we test the approximation under round-to-nearest using exact residuals. Across binary64, binary32, and software-emulated binary16 and bfloat16, the model recovers the ordering among tree topologies; KT tracks firstorder autoregressive (AR(1)) partial-sum costs. For general matrix–matrix multiplication (GEMM), independently calibrated predictions differ from measurements by at most 3% on the tested grid. A reduction tree extracted from an array library predicts the measured RMS scaling. However, stagnation and bias in positive low-precision sums limit the model’s applicability.
Keywords: floating-point arithmetic; reduction tree; rounding error; variance analysis; numerical stability
1.
Introduction
Every parallel sum has a reduction tree: a binary tree whose leaves are the inputs and whose internal nodes represent the pairwise additions that combine them. Consider Í inputs p1 , . . . , pk with exact sum s = ki=1 pi , and let bs be the computed result in floating-point arithmetic with unit roundoff u. An implementation may choose a sequential, pairwise, blocked, or hybrid tree for locality, synchronization, or bandwidth. We call this tree structure the reduction geometry. Although often chosen for performance, the ge-
This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).
ometry determines which intermediate values are rounded and how their errors accumulate. Two trees with the same leaves and operation count can therefore produce different errors. By contrast, the standard tree-independent forward-error bound ignores this geometry: |bs − s | ≤ γk−1
k ∑︁ i=1
| pi | ,
γn :=
nu , 1 − nu
for ( k − 1) u < 1. This bound provides a worst-case guarantee but no tree-specific mean-square estimate for random inputs. It therefore cannot rank a mixed library tree against a pairwise tree on random data; a second-moment model can. Our analysis keeps the geometry that the bound discards, as Figure 1 illustrates. Internal subtree sizes determine two scalar statistics, while shared ancestry determines a matrix-valued kernel. On a computational tree, the sum of squared exact partial sums, rather than operation count alone, controls typical error growth [1, 2]. To make this quantity precise, write p = ( p1 , . . . , pk ) ⊤ . For each internal node v of T ,
2
et al.
Piyush Sao Sequential
Two-stage blocked (b = 2)
Balanced pairwise
node label = subtree size |(v)|
8 7 6 5
8 6
3
4
2
leaf i
4
2
1
8
2
3
4
5
6
7
8
1
2
2
2
2
3
4
5
6
7
4
2
8
1
2
2
2
2
3
4
5
6
7
8
1
7
7
6
5
4
3
2
1
1
4
4
3
3
2
2
1
1
1
3
3
2
2
1
1
1
1
7
2
7
7
6
5
4
3
2
1
2
4
4
3
3
2
2
1
1
2
3
3
2
2
1
1
1
1
6
3
6
6
6
5
4
3
2
1
3
3
3
4
4
2
2
1
1
3
2
2
3
3
1
1
1
1
4
5
5
5
5
4
3
2
1
4
3
3
4
4
2
2
1
1
4
2
2
3
3
1
1
1
1
5
4
4
4
4
4
3
2
1
5
2
2
2
2
3
3
1
1
5
1
1
1
1
3
3
2
2
6
3
3
3
3
3
3
2
1
6
2
2
2
2
3
3
1
1
6
1
1
1
1
3
3
2
2
7
2
2
2
2
2
2
2
1
7
1
1
1
1
1
1
2
2
7
1
1
1
1
2
2
3
3
2
8
1
1
1
1
1
1
1
1
8
1
1
1
1
1
1
2
2
8
1
1
1
1
2
2
3
3
1
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
5 4
leaf j
leaf j
leaf j
tr KT = Λ1 = 35
tr KT = Λ1 = 26
tr KT = Λ1 = 24
1 ⊤KT1 = Λ2 = 203
1 ⊤KT1 = Λ2 = 132
1 ⊤KT1 = Λ2 = 112
3
shared internal ancestors (KT)ij = #{v : i, j ∈ (v)}
4
Figure 1. Shared ancestry determines the second-moment cost for three reduction trees on the same k = 8 inputs. Top: each internal-node label is its subtree size |L ( v)|; the blocked tree accumulates sequentially within blocks and across block totals. Bottom: ( KT )ij counts the internal ancestors shared by leaves i and j. The displayed diagonal and all-entry sums give Λ1 = tr KT and Λ2 = 1⊤ KT 1. Thus the sequential and blocked √︁ √︁ centered-RMS scales are 35/24 = 1.21 and 26/24 = 1.04 times the balanced-pairwise scale. Section 2.3 defines KT , and Section 2.4 derives the two statistics.
let qv be the exact sum of the leaves below v. The treeÍ dependent cost is v q2v , where the sum ranges over all internal nodes. We represent this cost as a quadratic form using the common-ancestor kernel KT , the k × k matrix whose entry ( KT )ij counts the internal nodes that are ancestors of both leaves i and j. Then ∑︁
q2v = p⊤ KT p.
v
For random inputs with second-moment matrix M := Ep [pp⊤ ] , the expectation over p yields the scalar contraction " Ep
# ∑︁ v
q2v
= Ep [ p⊤ KT p] = tr ( KT M ) = ⟨KT , M ⟩ ,
where ⟨X, Y ⟩ := tr ( X ⊤ Y ) is the Frobenius inner product. This contraction captures how the reduction geometry interacts with the input statistics. For i.i.d. inputs, it reduces to the two scalar statistics Λ1 and Λ2 ; the same construction extends to matrix multiplication. To analyze the tree√︁error eT := bs − s, define MSE ( eT ) := MSE ( eT ) . We call inputs centered when their means are zero and noncentered otherwise. The analysis separates three layers of increasing specialization. E[ eT2 ] and RMS ( eT ) :=
Layer 1 (exact recurrence). Standard stochastic rounding randomly selects one of the adjacent floating-point values that bracket an operand, with probabilities chosen to preserve that operand in expectation. Because stochastic
Floating-Point Reduction Trees
3
Table 1. Prediction workflow for a fixed binary reduction tree. Step Quantity 1
Extract T and its arithmetic at every node.
2
Accumulate Λ1 , Λ2 in one tree traversal.
3
Estimate µ, τ2 (or M for ⟨KT , M ⟩; §2.3).
4
Calibrate ν independently for the arithmetic regime.
5
Evaluate (1) and check model conditions.
rounding is conditionally unbiased, its accumulated error satisfies an exact recurrence with operand-dependent local variances (Theorem 2.2). Layer 2 (constant-ν model). The exact recurrence depends on the pre-rounding value at every node. We replace each operand-dependent local variance by νu2 x2 , where ν is a dimensionless local second-moment coefficient. This substitution preserves the tree structure while reducing the variance model to one coefficient (Corollary 2.3). Layer 3 (round-to-nearest calibration). Under deterministic round-to-nearest, local errors do not generally have conditional mean zero. We therefore treat the same leading expression as an empirical MSE approximation over random inputs and calibrate ν from independent nodelocal errors. Section 6 quantifies both its accuracy and its failure regimes. A fitted round-to-nearest value of ν serves as an empirical calibration parameter rather than an arithmetic constant.
1.1 Prediction Workflow For pure summation, the constant-ν model gives the following leading prediction on an extracted tree with i.i.d. summands of mean µ and variance τ2 : √︁ RMS ( eT ) ≈ u ν{τ2 Λ1 ( T ) + µ2 Λ2 ( T )}.
(1)
Evaluating this prediction requires three ingredients: the tree topology, the input statistics, and an independent calibration. We organize them into the validation workflow in Table 1. The same two-term decomposition identifies which tree statistic controls the prediction. Comparing its terms gives |µ| < τ
√︄
Λ1 ( T ) , Λ2 ( T )
|µ| Λ2 dominates if > τ
√︄
Λ1 ( T ) . Λ2 ( T )
Λ1 dominates if
For k = 4096, the closed forms in Section 2.4 give thresholds of approximately 0.019 for a sequential tree and 0.038 for a pairwise tree. Positive reductions such as softmax
denominators, squared norms, and batch statistics can therefore enter the Λ2 -dominated regime even when their mean is modest relative to their standard deviation. The criterion supplies the tree-dependent rounding component of a precision policy that also accounts for arithmetic, conditioning, and performance. We develop this framework in stages: Section 2 derives the exact recurrence and the common-ancestor kernel; Section 3 optimizes the resulting tree statistics; Section 4 extends the analysis to matrix multiplication; Section 5 maps implementations to trees; and Section 6 validates the predictions experimentally. Table 3 in Appendix 2 summarizes the core notation used throughout the scalar theory. We define symbols specific to hierarchy design, correlated inputs, GEMM, and calibration when they first appear.
2.
Conditional Second-Moment Model
Let T be a full binary reduction tree with leaf values p = ( p1 , . . . , pk ) ⊤ . Write int ( T ) for its internal nodes and h ( T ) for its height. For a node v, let L ( v) be its leaf set, Í qv = i ∈ L ( v ) pi its exact partial sum, b qv its computed value, and ev = b qv − qv its accumulated error. Expectations over rounding variables for fixed p are denoted by Eξ [· | p] . Averaging over random inputs then uses Ep . For deterministic round-to-nearest experiments, only Ep remains.
2.1 An Exact Conditional Recurrence To derive the exact recurrence for conditionally unbiased rounding errors and control its cross terms, we condition on the information available just before node v is rounded. We denote this information by the σ-algebra Fv− : it includes the fixed input p and all rounding draws in the two child subtrees of v, but excludes the fresh draw at v itself, hence the superscript “−”. Conditioning on Fv− fixes both the pre-rounding operand xv and the inherited error xv − qv . The cross term between that inherited error and the local error ξv therefore vanishes in expectation. This conditioning underlies the exact recurrence; the remaining cross term between the two child errors vanishes by subtree independence. Definition 2.1 (Conditionally unbiased rounding errors). Let flv denote the possibly random rounding map at node v. Let Fv− be the σ-algebra generated by the fixed input p and all rounding draws in the two child subtrees. With children v1 , v2 , the pre-rounding value xv = b qv1 + b qv2 = qv + ev1 + ev2
is Fv− -measurable. The local error ξv = flv ( xv ) − xv
4
Piyush Sao satisfies Eξ [ξv | Fv− ] = 0,
Eξ [ξ2v | Fv− ] = ψv ( xv ) ≥ 0.
Thus ψv ( x) is the conditional local variance for operand x. Rounding draws in disjoint subtrees are independent conditional on p. Rounded product leaves, when present, obey analogous conditions and are independent of later accumulation draws. Keeping the operand- and node-dependent functions
ψv explicit preserves the exact stochastic-rounding recurrence. The constant-ν model below replaces them with
homogeneous approximations.
Theorem 2.2 (Exact conditional second-moment recurrence). Define Vv = Eξ [ev2 | p] . Pure summation leaves have Vi = 0; a rounded product leaf has Vi = (p) Eξ [(ξi ) 2 | pi ] . Every internal node satisfies Vv = Vv1 + Vv2 + Eξ [ψv ( xv ) | p] .
et al.
For any node w, let d ( w) be its depth, measured as the number of edges from the root. Substituting the homogeneous model into (2) introduces the one-level propagation factor ω = 1 + νa u2a .
This factor measures how child-level variance is amplified as it passes through a rounded parent. Unrolling the recurrence from the root assigns each local term at w the geometric weight ωd ( w ) . Under h ( T ) u2a ≪ 1, these weights differ from unity by only O ( h ( T ) u2a ) . For binary64, ω − 1 = νa u2a ≈ 1.23 × 10−32 νa , so the weights are approximately one under the stated condition. Corollary 2.3 (Constant-ν model recurrence and solution). With the notation above, Vv = ω( Vv1 + Vv2 ) + νa u2a q2v , Vroot = νp u2p
(2)
k ∑︁
(4) ∑︁
ωd ( i ) p2i + νa u2a
i=1
ωd ( v ) q2v .
(5)
v ∈ int ( T )
If h ( T ) u2a → 0 and νa remains uniformly bounded, then Proof. Because ev = ev1 + ev2 + ξv , conditioning first on Fv− makes the cross term with ξv vanish. The child errors have zero conditional means and are independent on disjoint subtrees, so Eξ [ev1 ev2 | p] = 0. The remaining three second moments give (2). □ Standard stochastic rounding supplies a concrete ψv . If x ∈ [ x− , x+ ] lies between adjacent floating-point numbers, h = x+ − x− and θ = ( x − x− )/h, then Eξ [ξ | x] = 0,
ψ( x) = θ( 1 − θ) h2 .
(3)
Thus Theorem 2.2 is exact for standard stochastic rounding, including the operand-dependent interpolation coordinate θv .
2.2
Vroot = νp u2p
k ∑︁
! p2i + νa u2a
i=1
∑︁
q2v
h
i 1 + O h ( T ) u2a .
v ∈ int ( T )
(6)
Proof. Equation (4) follows from Eξ [xv2 | p] = q2v + Vv1 + Vv2 . Expanding from the root assigns each local term the weight ωd ( w ) . Since 0 ≤ d ( w) ≤ h ( T ) , this weight is 1+O ( h ( T ) u2a ) uniformly. Nonnegativity then gives (6). □ Theorem 2.2 gives the exact stochastic-rounding recurrence. Corollary 2.3 supplies the constant-ν model that Section 6 evaluates as a round-to-nearest MSE approximation over random inputs.
Constant-ν Model
We now specialize the exact recurrence by replacing operand-dependent local variances with a homogeneous model. Pure sums round only additions, whereas inner products may also round products. Subscripts p and a distinguish these operations: up , ua are their unit roundoffs and νp , νa are their nonnegative, dimensionless secondmoment coefficients. Because relative error scales with the operand, the homogeneous model uses variance proportional to u2 x2 : ψv ( x) = νa u2a x2 ,
Vi = νp u2p p2i ,
with separately calibrated product and accumulation coefficients.
2.3
Common-Ancestor Kernel
The quadratic cost v q2v can be expressed as a quadratic form in the input vector. This representation defines the common-ancestor kernel. Specifically, the coefficient of pi pj counts the internal nodes containing both leaves: their common internal ancestors. Figure 1 illustrates the resulting structure. Í
Proposition 2.4 (Common-ancestor kernel). Let bv ∈ {0, 1}k indicate the leaves below v, and define KT =
∑︁ v ∈ int ( T )
bv b⊤ v,
( KT )ij = #{v : i, j ∈ L ( v)}.
Floating-Point Reduction Trees
5 Where the mean takes over
Then
∑︁
q2v = p⊤ KT p.
dominance boundary |μ|/τ = √ Λ1/Λ2 ; exact family values joined by lines
(7)
v ∈ int ( T )
101
k = 4096
⊤
For the input second-moment matrix M = Ep [pp ] , " Ep
Λ2-dominated (mean-driven) for every family: softmax denominators, squared norms, batch statistics 100
# ∑︁
q2v = tr ( KT M ) = ⟨KT , M ⟩ .
(8)
v
Proof. Since qv = b⊤v p, summing q2v = p⊤ bv b⊤v p proves the first identity. Taking Ep and using the trace inner product proves the second. □
|μ|/τ
10−1
blocked b≍ √ k
0.038
sequential
0.019
10−2
The kernel KT is a leaf-indexed specialization of graphtheoretic ancestral matrices [3]: it counts internal common Í ancestors. Because a single traversal accumulates v q2v and the scalar statistics below, KT need not be formed densely.
pairwise
10−3 Λ1-dominated (variance-driven) for every family: centered gradients, zero-mean noise 10−4 102
104
103
105
106
reduction size k
2.4 I.i.d. Compression For i.i.d. inputs, the kernel contraction simplifies to two scalar tree statistics because M has only an identity and a rank-one component. Corollary 2.5 (I.i.d. tree statistics). If Ep [pi ] = µ and
Varp ( pi ) = τ2 , define ∑︁ Λ1 ( T ) = |L ( v)| , v ∈ int ( T )
Then
∑︁
Λ2 ( T ) =
|L ( v)| 2 .
v ∈ int ( T )
⟨KT , M ⟩ = τ2 Λ1 ( T ) + µ2 Λ2 ( T ) .
Figure 2. Dominance boundary for √︁ i.i.d. pure summation. Each curve connects exact values of |µ|/τ = Λ1 /Λ2 : integer sizes for sequential trees, powers of two for pairwise trees, and square sizes k = b2 for two-stage blocked trees. Above a curve, the µ2 Λ2 term exceeds τ2 Λ1 for that family; below it, the reverse holds. The marked k = 4096 thresholds are 0.019 (sequential) and 0.038 (pairwise). The workload labels illustrate typical regimes without imposing a distribution. Section 2.4 derives the two terms, and Section 1.1 states the criterion.
Their optimization therefore produces different schedules for centered and noncentered inputs.
(9)
Including rounded product leaves, the leading second moment under the constant-ν model is νp u2p k (τ2 + µ2 ) + νa u2a {τ2 Λ1 + µ2 Λ2 }.
Equivalently, Λ1 = tr KT and Λ2 = 1⊤ KT 1 are the diagonal and all-entry sums displayed in Figure 1. Figure 2 shows where the variance and squared-mean terms exchange dominance for three tree families. For a two-stage blocked tree with sequential accumulation within blocks and across block totals, the formulas below apply to the divisible case nb = k/b. For nondivisible cases, we evaluate Λ1 and Λ2 on the exact tree, including its shorter final block.
3.
Tree Statistics and Optimization
We optimize Λ1 and Λ2 over several restricted families of reduction trees to identify structures that minimize the predicted RMS error in centered and noncentered regimes and connect these choices with classical blocking and scheduling results. For centered√︁i.i.d. inputs, the leading RMS scale is proportional to Λ1 ( T ) . To formalize this scale asymptotically, let K ⊆ N be an unbounded set of input sizes, and let T• = ( Tk )k ∈ K be a family with one k-leaf tree for each k ∈ K . When the limit exists, we define the exponent α( T• ) by α( T• ) =
log Λ1 ( Tk ) 1 . lim 2 k→∞ log k
(10)
k∈ K
Tree
Λ1 ( T )
Λ2 ( T )
Sequential Pairwise, k = 2j Two-stage blocked
k ( k + 1 )/2 − 1 k log2 k k2 /( 2b ) + kb/2 + O ( k )
k ( k + 1 ) ( 2k + 1 )/6 − 1 2k2 − 2k k3 /( 3b ) + kb2 /3 + O ( k2 + kb )
The statistics have different meanings: Λ1 is total leaf depth, while Λ2 weights large subtrees quadratically.
The factor 1/2 converts second-moment growth into an RMS exponent. This asymptotic exponent α differs from b α, a log–log slope fitted to a tree-statistic scale or measured RMS over a specified finite grid. Logarithmic and finitesize drift can separate them; pairwise RMS, for example, √︁ has exponent 1/2 but includes a log k factor. Interpreting the combinatorial prediction arithmetically requires k →
6
Piyush Sao
∞, u → 0, and h ( Tk ) u2 → 0, with the local coefficient
bounded above and below. In the subsections below, we derive exponents 1/2, 2/3, 3/4, and 1 for pairwise, threestage, two-stage, and sequential families, respectively.
3.1 Extremal Trees We first identify the trees that minimize and maximize Λ1 . A counting argument gives an alternative formula. Each internal node v contributes its leaf count |L ( v)| to Λ1 . Equivalently, each leaf i is counted once for each of its di internal ancestors. Summing over leaves gives Λ1 ( T ) =
k ∑︁
di .
(11)
i=1
The right-hand side is the tree’s external path length, the sum of all root-to-leaf path lengths. This formulation reveals the extremal structure: nearly balanced trees minimize the external path length, whereas the sequential tree maximizes it. Because the centered RMS scale is propor√ tional to Λ1 , the same ordering applies to the predicted error. Proposition 3.1 (Λ1 extrema). Among full binary trees with k leaves,
leading-order constant-ν regime, where ωdi = 1 + o ( 1) uniformly in i, Huffman coding with weights τ2i miniÍ Í mizes i τ2i di [4, 5]. If W = i τ2i , wi = τ2i /W , and H2 ( w) is binary entropy, then WH2 ( w) ≤ min T
k ( k + 1) − 1. 2
Nearly complete balanced trees attain the minimum; the sequential tree attains the maximum. Consequently, every family with a well-defined exponent has α ∈ [1/2, 1] . Every intermediate exponent is also attainable. To attain a target α ∈ ( 1/2, 1] , combine pairwise blocks of size g = k2−2α sequentially. The outer stage contributes Θ( k2α ) to Λ1 , while the inner stages contribute O ( k log g) . The balanced family attains the lower extremum.
3.2 Heterogeneous Variances When independent centered inputs have heterogeneous variances, the optimal tree minimizes a variance-weighted external path length. Specifically, Equation (8) becomes Í 2 i τi di . This objective is mathematically identical to the expected-codeword-length objective for a binary prefix code in which symbol i has probability proportional to τ2i and codeword length di . The equivalence lets us import the classical Huffman optimality result directly.
3.3
τ2i di < W {H2 ( w) + 1}.
i
Two-Stage Sequential Blocking
In two-stage sequential blocking, each block and the block totals are summed sequentially. Here and below, f ( k) ≍ g ( k) means that cg ( k) ≤ f ( k) ≤ Cg ( k) for positive constants c, C independent of k. With b = kβ , Λ1 ∼
k2 kb + ≍ k2−β + k1+β . 2 2b
The first power decreases with β and the second increases, so their maximum is V-shaped. Balancing the two terms yields the centered optimum. Proposition 3.3 (Centered blocked exponent). For 0 ≤ β ≤ 1, ατ (β) =
1 max ( 2 − β, 1 + β) . 2
The two branches meet at β = 1/2, where b ≍ ατ = 3/4.
√
k and
Under the expected-RMS objective, Proposition 3.3 √ recovers the balance b ≍ k found in classical worst-case analyses of blocked summation [8, 9] and in the two-stage member of the Superblock family [10]. This balance yields the exponent ατ = 3/4 for sequential accumulation within blocks and across their totals. For noncentered inputs, the mean contribution instead uses Λ2 ∼
Proposition 3.2 (Huffman optimum and entropy bound). Assume independent centered inputs, an unconstrained full binary tree, and free leaf assignment. In the
∑︁
Kao and Wang optimize a deterministic magnitudeweighted objective. For same-sign data, it reduces to weighted external path length with weights | pi | ; for mixed signs, it depends on signed partial sums [6]. Proposition 3.2 instead optimizes the variance-weighted expected Í squared-partial-sum objective i τ2i di for independent centered inputs. With a maximum-depth constraint, the package-merge algorithm—a length-limited Huffman construction—gives the optimal tree [7].
Λmin = k ⌊ log2 k⌋ + 2{k − 2 ⌊ log2 k ⌋ }, 1 Λmax = 1
et al.
k3 kb2 + , 3 3b
αµ (β) =
1 max ( 3 − β, 1 + 2β) . 2
The exponent αµ (β) is minimized at β = 2/3, where it equals 7/6. The subscripts τ and µ denote the inputvariance and input-mean contributions, respectively.
Floating-Point Reduction Trees
7
Thus a single blocked construction requires different block growth rates in the Λ1 - and Λ2 -dominated regimes.
3.4 Fixed-Stage Hierarchies We extend the optimization to fixed-stage hierarchies with L sequential reduction levels and L − 1 nested blocking Î levels. Let rℓ be the stage ratios, Lℓ=1 rℓ = k, and let sℓ = Î j ≤ℓ rj . Summing the stage contributions yields L
Λ1 ∼
L
k ∑︁ rℓ , 2 ℓ=1
Λ2 ∼
k ∑︁ sℓ rℓ . 3 ℓ=1
(12)
L 1+1/L k , 2
ατ ( L) =
1 1 1+ . 2 L
(13)
Optimizing Λ1 independently recovers the centered geometric schedule selected by the Superblock optimizer [10]. Superblock controls worst-case exposure, whereas the schedule defined by (13) minimizes the expected partialsum-square statistic. The fixed-L limit is distinct from the pairwise regime L ≍ log k. For noncentered inputs, however, the cumulative factors sℓ in (12) make equal stage ratios suboptimal. We seek exponents a1 , . . . , aL ≥ 0 such that L ∑︁
rℓ = kaℓ ,
Set A0 = 0 and Aℓ = scales as kEℓ ,
j = 1 aj .
The ℓth term ksℓ rℓ in (12)
Eℓ = 1 + Aℓ−1 + 2aℓ .
Thus the polynomial exponent of Λ2 is maxℓ Eℓ . We therefore choose the stage ratios to minimize maxℓ Eℓ . Proposition 3.4 (Noncentered fixed-stage optimum). The smallest polynomial exponent of Λ2 over these feasible stage exponents is 1 + ( 1 − 2−L ) −1 , attained by aℓ =
2−ℓ . 1 − 2−L
The corresponding RMS exponent is αµ ( L) =
2L+1 − 1 . 2L+1 − 2
(14)
Proof. Let t bound every stage-dependent part, so that Aℓ−1 + 2aℓ ≤ t. Since Aℓ = Aℓ−1 + aℓ , this inequality is equivalent to Aℓ ≤
Aℓ−1 + t . 2
rℓ = k2
/( 2L − 1 )
.
Equivalently, the cumulative-exponent increments Aℓ − Aℓ−1 = aℓ double when read from the last stage back toward the first. This factor-of-two pattern motivates the name doubling-gap schedule. Within this fixed-stage family, the centered geometric schedule optimizes Λ1 , whereas the doubling-gap schedule optimizes the polynomial order of Λ2 . The sequential, optimally blocked two-stage, three-stage geometric, and pairwise families have centered exponents 1, 3/4, 2/3, and 1/2 (with a logarithmic factor), respectively. Figure 3 collects the exact centered growth curves and the corresponding centered-to-noncentered exponent changes.
4. General Matrix–Matrix Multiplication (GEMM)
aℓ = 1.
ℓ=1
Íℓ
Because later stages act on larger accumulated partial sums, the Λ2 objective forces their ratios to shrink. At the optimum, the stage exponents satisfy aℓ+1 = aℓ /2. Thus they halve from one level to the next, yielding L −ℓ
For centered inputs, the arithmetic–geometric mean inequality shows that Λ1 is minimized when the stage ratios are equal, rℓ = k1/L . This choice yields Λ★1 ∼
Iterating from A0 = 0 yields Aℓ ≤ t ( 1 − 2−ℓ ) . The feasibility condition AL = 1 therefore requires t ≥ ( 1 − 2−L ) −1 . We attain this lower bound by requiring equality at every stage, which yields the stated aℓ . The full Λ2 exponent is 1 + t; halving it to pass from the second moment to RMS gives (14). □
We extend the scalar theory to general matrix–matrix multiplication (GEMM). Let C = AB with A ∈ Rm×k and B ∈ Rk×n . Each output entry is an inner product whose k products are accumulated on the same tree T . The scalar recurrence applies entrywise, while a matrix formulation exposes how the operand data interact with the tree geometry. In particular, an algebraic identity expresses the exact-partial-matrix scale as a contraction of KT with a matrix G formed from the operand Gram matrices. Combined with ensemble variance comparability, this Hadamard-product structure controls the expected squared Frobenius error.
4.1 Exact Identity and Scaling Conditions By summing the scalar recurrence over output entries, we obtain an exact second-moment identity for GEMM without a scaling model.
8
Piyush Sao (a) Centered RMS growth, exact tree statistics sequential: k 1
3/2 1.4
b = 32: k 3/4 →k 1
105
b = √ k : k 3/4 3-stage: k
104
2/3
pairwise: k 1/2√ log2 k
103
(b) Exponent ladder, both regimes 1.6
RMS exponent RMS ∼ k α
√ Λ1
(centered RMS ÷ uτ√ ν )
10
6
et al.
1.2
1.0
0.8
7/6 15/14 1
sequential 1
two-stage 3/4 three-stage 2/3
0.6
102
pairwise 1/2
0.4
k = b 2: local slope = 34 exactly
102
104
103
105
centered ατ (Λ1)
106
reduction size k
noncentered αμ (Λ2)
√ Figure 3. Tree family and stage schedule jointly determine the RMS power law. Left: exact Λ1 for sequential, fixed-b = 32 blocked, optimized two-stage blocked, three-stage geometric, and pairwise families. The fixed-b curve has local slope 3/4 at k = b2 and approaches slope one, √ whereas b ≍ k retains exponent 3/4. Right: centered exponents ατ from Λ1 and noncentered exponents αµ from Λ2 ; the two- and three-stage noncentered values use their doubling-gap optima. Equation (10), Proposition 3.3, and Proposition 3.4 yield the displayed values.
Theorem 4.1 (Exact GEMM second moment). Under Definition 2.1, fix A, B. Let Σp ( A, B) sum the expected conditional variances of all rounded product leaves and let Σa ( A, B) sum those of all accumulation nodes. Then h
i
b − C ∥ 2 | A, B = Σp ( A, B) + Σa ( A, B) . Eξ ∥ C F
(15)
Proof. Apply Theorem 2.2 to every output entry and sum the conditional second moments. □
exact partial matrix accumulated at node v. Define P ( A, B) =
k ∑︁
ℓ=1
AT ( A, B) =
∥ A:,ℓ ∥ 22 ∥ Bℓ,: ∥ 22 ,
∑︁ v ∈ int ( T )
∥ Qv ∥ 2F .
Here P is the total squared product-leaf scale and AT the squared exact-partial-matrix scale. There exist positive comparison constants, independent of m, n, k, and T , such that EA,B [Σp ( A, B)] ≍ u2p EA,B [P ( A, B)] ,
We need an ensemble-level condition to convert this identity into a two-sided u2 x2 comparison. Such a comparison can fail when many nonzero operations are exact. For example, take m = n = k and A = B = Ik in binary arithmetic. Every product and partial sum is exactly representable, so all local variances vanish even though the squared product-leaf and exact-partial-matrix scales are positive. We exclude this aggregate degeneracy by requiring ensemble variance comparability. Assumption 1 (Ensemble variance comparability). Let ( A, B) be drawn from an operand ensemble independently of the rounding draws. For Sv ⊆ {1, . . . , k}, the set of inner indices below v, let Qv = A:,Sv BSv ,: be the
EA,B [Σa ( A, B)] ≍ u2a EA,B [AT ( A, B)] .
When we do not separately round products, we set both Σp and u2p EA,B [P ( A, B)] to zero. Here X ≍ Y means cY ≤ X ≤ CY , with positive constants that may depend on the ensemble and arithmetic but not on m, n, k, T . The upper comparisons impose the usual u2 x2 scale, while the lower comparisons exclude the aggregate degeneracy described above. Thus Assumption 1 makes the ensemble-and-arithmetic dependence explicit. In Section 6, we evaluate the resulting comparison on each stated operand ensemble and size grid using independent calibration.
Floating-Point Reduction Trees
9
Corollary 4.2 (GEMM scaling form). Under Assumption 1, h
b − C ∥2 EA,B,ξ ∥ C F
i
=Θ u2p EA,B [P ( A, B)]
+ u2a EA,B [AT ( A, B)] .
Proof. Take EA,B in Theorem 4.1 and apply Assumption 1. □
4.2
h i b − C ∥ 2 = Θ mnτ4 {u2p k + u2a Λ1 ( T )} . EA,B,ξ ∥ C F
Proposition 4.3 (GEMM common-ancestor contraction). Define G = ( A⊤ A)◦( BB⊤ ) ,
Then
∑︁
Gℓr = ( A:,ℓ ) ⊤ A:,r Bℓ,: ( Br,: ) ⊤ . ∥ Qv ∥ 2F = ⟨KT , G⟩ ,
(16)
v ∈ int ( T )
and the product-leaf scale is tr ( G) . Proof. Expanding ∥ Qv ∥ 2F yields ℓ,r ∈ Sv Gℓr . Summing over nodes counts each pair by its number of common ancestors, ( KT )ℓr . □ Í
Both Gram matrices are positive semidefinite. The Schur product theorem states that their Hadamard product Í is positive semidefinite, so G ⪰ 0. Because KT = v bv b⊤v is also positive semidefinite, ⟨KT , G⟩ ≥ 0. Forming G densely requires O ( k2 ( m + n)) work and O ( k2 ) storage, though structured representations can reduce this cost. For fixed or correlated operands, evaluating the tree-dependent scale requires the full contraction ⟨KT , G⟩ whenever the off-diagonal structure of G is material. For independent centered entries with variance τ2 , however, EA,B [G] = mnτ4 Ik ; the expected accumulation scale collapses to mnτ4 Λ1 ( T ) and can be evaluated without forming either dense matrix.
Centered Random Matrices
For centered random matrices with independent entries, the GEMM contraction simplifies and yields explicit RMS forward-error scales. Under Assumption 1, entries with
(17)
For sequential fused multiply–add (FMA) accumulation, Equation (17) gives an RMS error relative to ∥ C ∥ F of order √ u k. Standard forward-error bounds use the cancellationfree scale D = | A || B | , with entrywise absolute values, so Í that Dij = ℓ | Aiℓ Bℓj | . This cancellation-free scale measures the magnitude exposed to summation even when terms in Cij cancel. An experiment produces a separate forwarderror ratio on each trial, so we define
Kernel Contraction
We now express the exact-partial-matrix scale as a contraction of KT with a quantity derived from two Gram matrices. The matrix A⊤ A is the column Gram matrix of A, and BB⊤ is the row Gram matrix of B. Their Hadamard product, denoted ◦, multiplies corresponding entries.
4.3
variance τ2 and finite fourth moment give
ηF =
b − C ∥F ∥C . ∥D∥F
The second-moment theory instead compares a ratio of RMS quantities, b − C ∥ 2 1/2 EA,B,ξ ∥ C F RF := 1/2 . EA,B ∥ D ∥ 2F
Thus RF is a ratio of ensemble RMS values, with one population value for each ( m, n, k) configuration. By contrast, ηF is a per-trial random variable that changes with each draw of ( A, B) . The second-moment theory predicts RF , not E[ηF ] , because the latter is the expectation of a ratio of correlated random variables. In the experiments, we report RF to validate the model and the distribution of ηF to assess practical error. Proposition 4.4 (RMS forward-error scale). Under Assumption 1, let the entries of A and B be mutually independent copies of a centered random variable X with variance τ2 > 0 and E| X | > 0. For sequential FMA accumulation with u = ua , RF = Θ( u) .
Proof. For the matrix D defined above, each entry satisfies EA,B [ Dij2 ] = k (E| X | 2 ) 2 + k ( k − 1) (E| X |) 4 = Θ( k2 τ4 ) .
Hence EA,B ∥ D ∥ 2F = Θ( mnk2 τ4 ) , which has the same k2 scale as the sequential accumulation term in (17) after multiplication by u2 . □ These normalizations yield the following forwarderror scales for common inner-product accumulators. Tree Sequential √ Blocked, b ≍ k Pairwise
relative to ∥ C ∥ F √ Θ( u k) Θ( uk1/4 ) √︁ Θ( u log k)
relative to ∥| A || B |∥ F
Θ( u) Θ( uk−1/4 ) √︁ Θ( u log k/k)
10
Piyush Sao When off-diagonal entries contribute negligibly to
⟨KT , G⟩ , the tree-design objective reduces to the diagonal
surrogate, and the Huffman construction from Section 3 applies to diag ( G) . Otherwise, the objective remains the full contraction ⟨KT , G⟩ .
5.
Mapping Implementations to Trees
To predict errors, we map an implementation to its arithmetic tree. Library kernels, vector lanes, accelerator blocks, and communication collectives may introduce different trees at different scales. FPRev, a tool that infers floating-point accumulation order from compiled binaries, can automate parts of this reconstruction [11]. In the NumPy study, we reconstruct the tree manually and validate it bitwise. Once a tree is extracted, evaluating its statistics is inexpensive. Computing Λ1 and Λ2 requires one postorder traversal, while data-dependent analysis can accumulate Í 2 v qv during the reduction. Dense representations of KT , M , and G, however, cost O ( k2 ) , making them practical only for moderate-size or structured analyses.
5.1 A Library Reduction To test the model on a production library, we follow the workflow in Table 1: we reconstruct the reduction tree used by NumPy’s np.sum on a specific platform, compute its statistics, calibrate independently, and predict its error after bitwise validation. We reconstructed NumPy 1.26.4 on an Apple M3 Max arm64 system running macOS 26.5.1 and Python 3.12.4. On this build, the inner kernel accumulated segments of fewer than eight terms sequentially and used eight unrolled accumulators within a 128-element block. For longer segments, recursive halving split the segment near its midpoint, aligned the split to eight elements, reduced both halves, and combined their totals. An outer sequence combined 8192-element working buffers from NumPy’s universal-function (ufunc) machinery; np.setbufsize controls the default buffer size [12]. The resulting four-scale hierarchy is <8 ≤ 128 ≤ 8192 > 8192 −→ −→ −→ sequential eight lanes recursive halving buffer chain
To validate this hierarchy, we compare reconstructed and observed bit patterns on fixed-seed vectors, vary the ufunc buffer size with np.setbufsize, and repeat versioned builds. Section 6 reports the results of these validation steps and evaluates the error prediction on the extracted 1.26.4 tree. A version or platform change requires repeating both source reconstruction and bitwise validation.
et al.
5.2 Communication and Accelerator Hierarchies For a scalar collective across P ranks, a flat or ring-like arithmetic chain has sequential statistics, whereas a binomial or recursive-doubling tree has pairwise statistics. Message chunking and algorithm switching create nested trees; we extract their per-chunk arithmetic order to obtain the tree to which the formulas apply. We use the same extraction procedure for tensor-core blocks, single-instruction multiple-data (SIMD) lanes, GPU warpshuffle instructions that exchange values among threads within a warp, and epilogue accumulation. The rounding law and internal accumulation order determine whether Definition 2.1 represents the resulting arithmetic. 5.3 Workload Statistics Tree extraction and input modeling remain separate: after the tree is extracted, the input distribution determines the appropriate statistic. We use Λ1 for approximately centered inner-product terms, both Λ1 and Λ2 for noncentered reductions, and ⟨KT , M ⟩ for correlated gradients or spatial data. Softmax denominators, sums of squares used in norms, and batch-normalization statistics are common positive reductions whose mean-to-standard-deviation ratio may exceed the threshold in Section 1. Mixedprecision systems often retain selected reductions in binary32 [13]; the threshold supplies the tree-dependent rounding component of that precision-policy choice.
6.
Numerical Experiments
The exact recurrence (2) applies to conditionally unbiased stochastic rounding. The constant-ν model replaces operand-dependent local variances with constant coefficients and empirically approximates the MSE of deterministic round-to-nearest. Our experiments test structure and calibration separately: whether a tree statistic predicts how MSE changes, and whether an independently estimated coefficient predicts its level.
6.1 Design and Reproducibility Every reported summation error is the computed result minus its exact value. We use precision-specific strategies to evaluate these errors accurately. For binary64, Knuth’s error-free TWOSUM transformation records the exact residual at each node, and math.fsum combines those residuals at the roundoff-error scale. Python 3.12.4’s standard-library Decimal, at 80-digit precision, independently checks selected cases. For binary32, the exact sum of every node pair on the reported grids is representable in binary64. We therefore form its residual in binary64 before rounding to binary32.
Floating-Point Reduction Trees
11
GEMM references also use 80-digit Decimal, and we check the system libm fused multiply–add (FMA) bitwise against one-round Decimal emulation. Bootstrap intervals resample trials independently at each size. The three-stage (L = 3) ratios use a paired bootstrap because all schedules receive the same input. Table 2 summarizes the grids. We call an addition node stagnant when its rounded output equals one operand while the other is nonzero, so that the addition does not change the retained operand. To isolate arithmetic effects without emulating a hardware kernel, we implement a software binary16/bfloat16 emulator with round-to-nearest-even and gradual underflow. We validate random binary16 values against NumPy and bfloat16 values against an independent bit-rounding implementation. Both emulator modes pass exact-value, midpoint-tie, subnormal, overflow, and Decimal residualidentity tests. For pure-summation experiments, we define the effective coefficient νeff :=
MSE ( eT ) Í . u2 Ep [ v ∈ int ( T ) q2v ]
For i.i.d. inputs, the factor multiplying u2 in the denominator is τ2 Λ1 + µ2 Λ2 ; for correlated centered inputs it is ⟨KT , M ⟩ ; and for heterogeneous independent centered Í inputs it is i τ2i di . We define the squared-bias fraction as fbias :=
{Ep [ eT ]}2 Ep [ eT2 ]
,
and estimate it by dividing the squared sample mean by the sample mean square. Values near one indicate that systematic bias dominates MSE; values near zero indicate that trial-to-trial variation dominates.
6.2 Centered Scaling Across Precisions Figure 4 combines four precisions. For centered normal inputs, every precision preserves the ordering sequential > blocked > pairwise. The binary64 pairwise finite-grid theory slope is 0.579 and the measured slopes are 0.581 (normal) and 0.592 (uniform); sequential measurements √ are 1.000 and 0.980. The blocked- k measurements lie near the predicted 0.750. Binary32, software binary16, and software bfloat16 produce the same structural pattern. On this grid, the binary32 tests cover sampled uniform [−1, 1] and standard-normal operands. The positive-input panel identifies the onset of lowprecision stagnation. At k = 8192, a sequential binary16 sum returns one operand unchanged at 66.6% of internal nodes; bfloat16 reaches 95.8%. The corresponding squared-bias fractions exceed 0.9998, and the fitted slopes rise to 1.97 and 2.21. By contrast, pairwise accumulation has negligible stagnation and remains near its noncentered
exponent 1. Thus, on this grid, low precision preserves the centered topology law. Monotone positive sums leave the leading-model regime when low precision causes stagnation and bias. Across all centered configurations, the signed mean errors stored with the summaries contribute little to MSE. Large squared-bias fractions occur in the positive lowprecision configurations.
6.3 Blocked Families Fixed-b and growing-b predictions differ. For b = 32, the local slope passes through 3/4 near k = b2 and moves toward one as larger sizes enter the fit window. When b ≍ kβ , the fitted points follow the V-shaped curve of Proposition 3.3 (Figure 5). These results identify 3/4 as the optimum within the two-stage sequential blocking family, not as a universal summation exponent. A separate four-stage experiment uses nine sizes k = r 4 , r = 4, . . . , 12, 1000 normal trials per size, and 2000 bootstrap replicates. Its fitted slope is b α = 0.602 with interval [0.591, 0.614] , compared with the exact finite-grid √ Λ1 slope 0.619. Appendix 1.2 analyzes how νeff varies. 6.4 Correlated and Heterogeneous Inputs Correlated inputs require the full kernel contraction, whereas heterogeneous independent variances lead to the Huffman objective. For centered AR(1) inputs at k = 64, the empirical partial-sum-square cost remains within 3.6% of ⟨KT , M ⟩ for both sequential and pairwise trees (Figure 6). The i.i.d. statistic instead underpredicts the sequential cost by 14.2× at ρ = 0.9 and overpredicts it at ρ = −0.5. With the independent centered binary64 calibration, we predict the absolute RMS within 9.1% across all eight configurations. We test Proposition 3.2 using heterogeneous independent variances at k = 256. Figure 7 normalizes every weighted leaf-depth cost by the entropy lower bound WH2 ( w) . For τ2i ∝ i −2 , the Huffman, pairwise, and naturalorder sequential ratios are 1.02, 3.45, and 109; for τ2i ∝ e −i/8 , they are 1.01, 1.80, and 55.9. Relative to Huffman, pairwise trees cost 1.79–3.37× as much and sequential trees cost 55.5–106× as much. Exact-reference round-to-nearest measurements give νeff ∈ [0.160, 0.232] across the six configurations, and the deterministic cost ordering persists across this coefficient range. 6.5 Noncentered Inputs and Stage Schedules Noncentered normal, uniform, and chi-square inputs shift the fitted slopes toward the Λ2 predictions. Figure 8 compares measured and finite-grid slopes for normal ( 1, 1) and then tests the three-stage schedule optimizer. At k = 214 ,
12
Piyush Sao
et al.
Table 2. Experiment design. AR(1) denotes a first-order autoregressive input model. We fix trial counts in advance; all scripts, seeds, summaries, and compressed per-trial errors accompany the supplementary material. Study
Arithmetic and inputs
Size grid
Trials
Centered scaling
binary64/binary32; uniform, normal
64–8192
800/size
Primary comparison √ b α vs. Λ1
Software low precision
binary16/bfloat16; centered, positive
64–8192
800/400
MSE, bias, stagnation
Correlation
binary64 AR(1), ρ ∈ {−0.5, 0, 0.5, 0.9}
k = 64
4000
KT vs. i.i.d. model
Heterogeneous variance
binary64, two variance profiles
k = 256
400
Huffman/entropy bound
14
400
paired L=3 ratios
Noncentered schedules
binary64 normal( 1, 1)
k=2
Stochastic rounding
sums and 8 × k × 8 GEMM
stated grids
24–600
exact recurrence
GEMM
binary64, 32 × k × 32
32–256
50
independently calibrated MSE
Noncentered scaling
binary64; normal, uniform, chi-square
64–2048
300/size
finite-grid Λ1 –Λ2 slopes
Coefficient calibration
binary64; centered uniform, normal
{64, 256, 1024}
600/size
νeff
Local diagnostics
binary64; selected inputs and trees
selected 256–2048
binary64; centered uniform
300–500 1000 through 215 ; 300 thereafter
cross terms; significands
7
Extracted NumPy tree
the doubling-gap schedule ( 256, 16, 4) has the smallest deterministic Λ2 . The geometric schedule has predicted RMS ratio 1.192 and measured ratio 1.066 with paired interval [0.973, 1.168] . A nearby perturbation has predicted ratio 1.032 and measured ratio 0.954 with interval [ 0.877, 1.038] . Both paired intervals contain unity, so the measurements support the predicted cost scale without statistically distinguishing the schedule ranking under round-to-nearest. To determine whether deterministic correlations explain why the schedule rankings do not separate, we decompose the MSE as in Appendix 1.3 on related diagnostic grids. Let v < w enumerate each unordered node √ pair once. Across sequential, pairwise, and blocked- k Í trees at k = 256, 1024, 2048, the cross terms 2 v<w E[ξv ξw ] range from −10.2% to +5.2% of MSE. Their limited size does not by itself explain the unresolved ranking and instead leaves the configuration-dependent local coefficient as the larger observed source of variation. That appendix also reports the centered decomposition and supporting residual-identity tests.
6.6
Exact Stochastic-Rounding Check
In Figure 9, we compare measured second moments with the exact nodewise recurrence (2) using the local variance (3). Four of the six initial 95% intervals contain unity. For the two initial misses, we conducted post hoc diagnostic runs at larger fixed sample sizes; these runs also yielded intervals containing unity. Because we selected the follow-ups after inspecting the initial results, we report the initial and diagnostic intervals separately; they do not form a combined confirmatory confidence statement.
2 –2
20
extracted vs. pairwise tree
6.7 Calibrating νeff For centered i.i.d. round-to-nearest summation, the definition in Section 6.1 specializes to νeff =
MSE ( eT ) . u2 τ2 Λ1 ( T )
Across the binary64 calibration grid, the pooled mean is 0.1895, while individual configurations span approximately 0.15–0.25 (Figure 10). We use 0.19 as the calibration for centered binary64 configurations. All coefficients use unit roundoff u; normalization by machine epsilon ϵ = 2u would divide the reported coefficient by four. To explain the spread in νeff , we examine representable spacing. Write x = z2e with significand z ∈ [1, 2) ; its binade is the interval [2e , 2e+1 ) . Under the unit-roundoff convention used here, the spacing within this binade is h = 2u2e . Approximating roundoff as uniform on [−h/2, h/2] gives variance h2 /12 and hence relative local MSE u2 /( 3z2 ) . Averaging this z−2 coefficient under uniform and logarithmic significand laws gives 1/6 ≈ 0.167 and 1/( 8 ln 2) ≈ 0.180, respectively. Figure 10 (right) tests the cell-level prediction by normalizing out z−2 . The measured local second moments show the predicted dependence, with a mild downward trend across four bins. This mechanism explains the scale and configuration dependence, so the pooled coefficient remains configuration-dependent.
6.8 GEMM The FMA experiment uses m = n = 32, matrices with independent standard-normal entries (τ = 1), and k = 32, 64, 128, 256. Independent local-error runs calibrate accumulation coefficients νa ∈ [0.176, 0.184] . A separate product-and-addition rounding comparison calibrates
Floating-Point Reduction Trees
13
Centered normal inputs finite-grid theory binary64 binary32
1.0
Positive inputs expose low-precision stagnation 67% / 96% stagnant nodes at k = 8192
2.2
binary16 bfloat16
2.0
0.9
fitted α ̂
fitted α ̂
1.8 0.8
binary16 bfloat16
1.6
finite-grid √ Λ2
0.7
1.4
0.6
1.2 1.0
0.5
tial
uen
Seq
wise
Pair
ked Bloc
32
√ ked Bloc
k
Seq
tial
uen
wise
Pair
ked
Bloc
32
√ ked Bloc
k
√ Figure 4. Left: fitted centered-normal slopes with bootstrap 95% intervals over k = 64, . . . , 8192; black ticks are finite-grid fits √ of Λ1 . Right: fitted slopes for uniform [ 0, 1] inputs in software binary16 and bfloat16; black ticks are the noncentered finite-grid slopes of Λ2 . Sequential sums that stagnate at the largest size develop strong bias and depart from the constant-ν model. All blocked results use sequential accumulation within blocks and across block totals. Section 6.2 discusses the topology ordering and low-precision stagnation boundary.
νp ∈ [0.176, 0.184] and νa ∈ [0.180, 0.186] . All coefficients
are calibrated from node-local errors independently of the root error used for validation. For this τ = 1 ensemble, Figure 11 shows that Dmodel = mn u2 {νp k + νa Λ1 ( T )}
predicts mean squared Frobenius error within 3% for all eight configurations. The RMS ratio RF remains nearly constant in k, while the separate product rounding adds the predicted leaf contribution. Per-trial ηF central 95% ranges lie within [4.70, 5.75] × 10−17 . For the corresponding non-FMA dot product, the product/accumulation crossover is the value of k at which product rounding and accumulation contribute equally to the leading MSE. To calculate this crossover, we use independently calibrated constants rather than setting them equal. Appendix 1.1 reports this auxiliary finite-size effect.
6.9
Extracted NumPy Tree
Before testing MSE, we use intervention runs to verify the reconstructed tree. Changing np.setbufsize moves the outer boundary as predicted, and the NumPy 1.26.4 tree matches np.sum bit for bit on more than 1500 random vectors, including non-powers of two. Under NumPy 2.5.1, the observed sum matches the pure-pairwise reconstruction in all 20 trials at every tested size from 4096 to 262144, consistent with the contiguous-reduction fast path introduced during that development cycle [14]. Figure 12 evaluates Equation (1) on the reconstructed NumPy 1.26.4 node set. With ν = 0.19 fixed from the independent summation calibration, the model predicts
the error within 7% across k = 27 , . . . , 220 . The extracted tree captures the bend caused by sequentially combined 8192-element buffers; a pure-pairwise model does not. The result applies to the specified NumPy 1.26.4 build and platform.
7.
Related Work
7.1
Stability and Accurate Summation
Classical analyses distinguish recursive, pairwise, and compensated summation and provide worst-case forward or backward guarantees [8, 15]. Compensated algorithms, expansions, and superaccumulators change the arithmetic to obtain higher accuracy or reproducibility [16, 17, 18, 19]. FABsum combines blocked summation with an accurate final stage to obtain an operation-count-independent backward bound [9]. By contrast, we study the typical second moment of uncompensated accumulation on a fixed tree. The Superblock family provides a structural parallel for hierarchical dot-product schedules, but under a worstcase exposure objective [10]. Its centered geometric optimizer coincides with our Λ1 -optimal fixed-stage hierarchy, whereas our second-moment objective additionally distinguishes the Λ2 -optimal doubling-gap schedule. Earlier ordering analyses for nonnegative inputs anticipate different quadratic and cubic mean-square contributions [20, 21].
14
Piyush Sao V-shaped b ≍ k β continuum
Fixed block size b = 32 0.85
0.95 0.90
0.75
α̂
windowed α ̂
1.00
finite-grid theory uniform normal
0.80
et al.
0.70
0.85
0.65
0.80
0.60
0.75
theory measured
0.70
k 64--1
k 28--2
k 56--4
1
0.0
k 12--8
2
5
0.2
0.4
0.6
0.8
1.0
block exponent β
Figure 5. Binary64 round-to-nearest, centered inputs. Left: four overlapping fit windows for fixed b = 32. Right: block growth b ≍ kβ . Error bars are bootstrap 95% intervals; dashed curves are finite-grid or asymptotic tree-statistic predictions. Section 6.3 discusses the fixed- and growing-block fits.
i.i.d. Λ1 prediction
⟨KT, M⟩ prediction sequential pairwise
101
empirical / i.i.d.
empirical / kernel
1.03 1.02 1.01 1.00
100
0.99 0.98 −0.4
−0.2
0.0
0.2
0.4
0.6
0.8
−0.4
AR(1) correlation ρ
−0.2
0.0
0.2
0.4
0.6
0.8
AR(1) correlation ρ
Figure 6. First-order autoregressive (AR(1)) inputs. The full common-ancestor kernel tracks the empirical partial-sum cost (left), while the i.i.d. Λ1 compression departs from the empirical cost as correlation increases (right, logarithmic scale). Section 6.4 discusses the correlation test.
7.2
Probabilistic Roundoff
Probabilistic models of roundoff date to von Neumann and Goldstine [22]. Henrici, Hull and Swenson, and Barlow and Bareiss studied how roundoff errors propagate and how they are distributed [23, 24, 25]. Chatelin and Brunet applied a probabilistic propagation model to the eigenvalue problem [26]. Higham and Mary impose a centering assumption on roundoff errors, whereas we distinguish centered from noncentered input summands [1]. At the computational-tree level, Hallman and Ipsen √︁Í derive the first-order functional u v q2v [2]. We take its expectation under an input second-moment model to obtain ⟨KT , M ⟩ . For i.i.d. inputs, this contraction reduces to Λ1 and Λ2 . We then optimize those statistics and extend the contraction to GEMM. For stochastic rounding,
local errors are conditionally unbiased by construction [27, 28]. El Arar et al. give a general stochastic-rounding variance framework [29]. For reduction trees, we derive the common-ancestor kernel, its Λ1 , Λ2 compression, and the resulting optimization and GEMM formulas.
7.3
Ordering, Extraction, and Hardware
Kao and Wang optimize deterministic magnitudeweighted objectives for same-sign summation and study the mixed-sign problem [6, 30]. Those objectives differ from the variance-weighted depth and squared-partialsum objectives here. ReproBLAS and ExBLAS target bitwise reproducibility or exact accumulation rather than a typical-error prediction [31, 32]. Statistical MPI studies show that implementation topology changes observed errors [33]. FPRev demonstrates that the topology can be
Floating-Point Reduction Trees
15
2 i
i 2
eff = 0.16
e i/8
entropy bound
103
eff = 0.17
i
2 i di
102
i
2 i di
entropy bound
2 i
eff = 0.16
101
102
eff = 0.23 eff = 0.20
eff = 0.23
Huffman Pairwise Sequential
Huffman Pairwise Sequential
Figure 7. Variance-weighted leaf-depth costs for two heterogeneous profiles (k = 256). The shaded interval is the entropy bracket [WH2 ( w) , W {H2 ( w) + 1}]; labels above the bars give exact-reference round-to-nearest MSE coefficients. Proposition 3.2 gives the bracket, and Section 6.4 discusses the heterogeneous-variance comparison.
reconstructed [11]. Tensor-core and block-FMA analyses model arithmetic inside particular accelerator units [34]. Mixedprecision training and stochastic-rounding studies motivate low-precision reductions [13, 35]. Recent work studies Horner evaluation and pairwise summation under limited-precision stochastic rounding [36]. We analyze conditionally unbiased stochastic rounding exactly and evaluate the constant-ν model under deterministic round-to-nearest-even in the binary16/bfloat16 experiments. These emulation experiments isolate precision effects; applying the framework to a deployed accelerator begins by extracting its arithmetic tree.
8.
Summary and Limitations
The common-ancestor kernel combines a reduction tree with an input second-moment matrix to produce the expected partial-sum-square cost. For i.i.d. inputs, the resulting cost depends on the tree only through Λ1 and Λ2 , which separate the centered and noncentered regimes. This representation yields extremal-tree, blocking, fixedstage, and variance-weighted Huffman results and extends to GEMM through a Hadamard product of Gram matrices. Stochastic-rounding diagnostics test the exact recurrence. Under round-to-nearest, an independently calibrated constant-ν model predicts topology ordering and absolute MSE within the tested configurations. The low-precision and four-stage results identify regimes in which coefficient drift, bias, or stagnation produces departures from that approximation.
8.1
Data and Code Availability
An archived supplementary package contains the pinned environments, scripts, seeds, summary data, compressed per-trial errors, validation tests, figure-generation code, and reproduction instructions for all reported results.
8.2 Limitations • The leading second-moment formulas characterize random inputs on fixed binary trees and complement worstcase methods for adversarial or data-dependent orderings. • The constant-ν model represents operand-dependent local variance through a calibrated coefficient. Under round-to-nearest, the experiments quantify its dependence on distribution, tree, correlation, and significand stream. • The implementation evidence covers CPU software arithmetic, a reconstructed NumPy kernel, and validated binary16/bfloat16 emulation with gradual underflow. Applying the framework to GPU, tensor-core, and MPI implementations begins with extracting and validating their arithmetic trees. • The leading expansion applies when h ( T ) u2 ≪ 1; the positive low-precision experiments identify stagnation and bias as the dominant corrections outside this regime. The framework opens extensions to concentration bounds, compensated trees, combined Λ1 –Λ2 optimization, structured GEMM contractions, and extracted multi-level hardware hierarchies.
Appendix 1.
Auxiliary Results
16
et al.
Piyush Sao Nonzero-mean Λ2 regime
Three-stage schedule test Seq.
ratio to doubling gap
1.5
measured α ̂
1.20
1.4 1.3 Blk. √ k
1.2
Blk. k 2/3
1.1 Pair.
1.0
1.15 1.10 1.05 1.00 0.95
predicted RMS ratio measured RMS ratio
0.90
0.9 0.9
1.0
1.1
1.2
1.3
1.4
1.5
finite-grid theory
etric
ap
ling g
Doub
bed
r Pertu
Geom
Figure 8. Left: noncentered normal fitted slopes versus exact finite-grid tree-statistic slopes. Right: deterministic and measured three-stage √︃ schedule ratios. The predicted RMS ratio is Λ2 /Λ2,dg ; measured error bars are paired bootstrap 95% intervals. Section 6.5 discusses the slope and schedule comparisons.
Stochastic-rounding recurrence
obtain
Sum seq., μ = 0
αloc ( k) =
Sum pair., μ = 0
k3/4 =
2cp ca
.
(18)
Dividing the pooled MSE coefficients by u2 yields cp = 0.1801 and ca = 0.1831, giving k3/4 = 1.97. Figure 13 plots this curve alongside the three measured local-slope intervals. This finite-size crossover is distinct from the asymptotic 3/4 optimum of two-stage sequential blocking.
Sum seq., μ = 1 Sum pair., μ = 1 GEMM, μ = 0 initial selected follow-up
GEMM, μ = 1 0.9
1.0
1.1
1.2
measured [e ] / recurrence
1.3
2
Figure 9. Measured stochastic-rounding second moment divided by the exact nodewise recurrence. Circles show all six initial bootstrap 95% intervals; open squares show post-hoc larger-sample diagnostics for the two initial misses. Section 6.6 discusses the initial and diagnostic intervals.
Appendix 1.1
1 cp + ca k , 2 cp + ca k/2
Product and Accumulation Crossover
For a centered sequential non-FMA dot product, let cp and ca denote the independently calibrated product and accumulation coefficients with their common input-variance and unit-roundoff factors divided out. The leading MSE is proportional to cp k +
ca 2 k . 2
Define the local RMS exponent as the instantaneous log– d log RMS ( k ) log slope αloc ( k) := d log k . Taking the square root of the leading MSE and differentiating logarithmically, we
Appendix 1.2 Four-Stage Scaling Experiment The fixed-design four-stage experiment uses geometric trees with k = r 4 for r = 4, . . . , 12, centered normal inputs, 1000 trials per size, 2000 bootstrap replicates, fixed seeds, and exact residual references. It records the exact Λ1 , the signed mean, the squared-bias fraction, νeff , and compressed per-trial errors. The fitted slope is 0.602 with 95% interval [0.591, 0.614] , below the exact finite-grid √ Λ1 slope 0.619 and the asymptotic value 0.625 (Figure 14). Despite intermediate fluctuations, the effective coefficient νeff decreases from 0.217 at k = 256 to 0.190 at k = 20736. This decrease shifts the finite-grid round-to-nearest slope. Because Equation (13) depends only on the tree’s combinatorics, the observed drift in νeff on this grid explains the discrepancy. Appendix 1.3 Local-Error Covariances and Unit Tests For deterministic round-to-nearest, the MSE expands to MSE ( eT ) =
∑︁ v
Ep [ξ2v ] + 2
∑︁
Ep [ξv ξw ] ,
v<w
where v<w ranges over each unordered pair of distinct internal nodes exactly once. For centered k = 512, Í
Floating-Point Reduction Trees
17 Consistent with z −2 dependence
Configuration-dependent calibration uniform, sequential uniform, pairwise normal, sequential normal, pairwise
νeff
0.22
1.05
uniform-significand 1/6 Benford 0.180 pooled 0.1895
([ξ 2]/x 2)z 2/(u 2/3)
0.24
0.20 0.18 0.16
1.00 0.95 0.90 0.85 0.80 0.75
102
1.1
103
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
significand z ∈ [1, 2)
k
Figure 10. Left: effective centered binary64 MSE coefficient across size, tree, and input distribution, with three reference levels. Right: normalized local second moment versus significand; a horizontal value of one is the uniform-in-cell prediction. Section 6.7 discusses the calibration spread and significand mechanism.
Two-term model
Forward-error scale FMA separate product/add
5.5
1.02
5.4
1017RF
measured / model MSE
1.04
1.00 0.98 0.96
FMA separate product/add
0.94 32
64
128
256
5.3
5.2
32
64
k
128
256
k
Figure 11. Fused multiply–add (FMA) and separate product/addition rounding. Left: measured MSE divided by the independently calibrated two-term model. Right: RF , the ratio of RMS quantities. Error bars are bootstrap 95% intervals. Section 6.8 discusses the independent calibration and validation.
the covariance contribution ranges from −4% to +13% across two distributions and two trees. For noncentered normal ( 1, 1) at k = 256, 1024, 2048, it ranges from√−10.2% to +5.2% across sequential, pairwise, and blocked- k trees. These bounded contributions support the local-variance approximation, while the remaining variation reflects the operand-dependent effective coefficient. The unit tests cover both tree statistics and arithmetic. They reproduce every closed-form Λ1 value exactly on the tested trees. The low-precision software tests compare binary16 with NumPy on 12,000 random values and bfloat16 with an independent bit-level reference on 10,000 binary32 values. Separate tests exercise midpoint, subnormal, and overflow boundaries and verify the complete reduction residual against Decimal arithmetic.
Appendix 2.
Core Notation
Table 3 summarizes the core notation reused across the scalar theory, including the conditioning and asymptotic conventions. Symbols specific to hierarchy design, correlated inputs, significand calibration, communication costs, or GEMM are defined locally.
Acknowledgement This research was supported by the 2026 Laboratory Directed Research and Development (LDRD) Program at Oak Ridge National Laboratory through the initiative “CCSD Core: Foundational Research for Smart Extremescale Ecosystems.” Oak Ridge National Laboratory is managed by UT-Battelle, LLC, for the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.
18 Table 3. Core notation used across Sections 2–3.
Product/accumulation crossover k = 1.97
Symbol
Meaning
T, L ( v) , int ( T ); h ( T ) , d ( v)
tree, leaf set, internal nodes, height, and depth from the root
pi , qv , b qv
exact input, exact partial sum, and computed partial sum
flv , Fv− , ψv
rounding map, pre-rounding information, and local variance
ev , ξv , Vv
accumulated error, local roundoff, and conditional second moment
up , ua ; νp , νa
product and accumulation unit roundoffs and the corresponding local secondmoment coefficients
ω
one-level variance propagation factor 1 + νa u2a
KT , M ; ⟨·, ·⟩
common-ancestor kernel, input second moment, and trace inner product
Λ1 , Λ2
the sum of internal subtree sizes and the sum of their squares
µ, τ2
i.i.d. input mean and variance
α, b α, αloc ( k)
asymptotic exponent, finite-grid fit, and local log–log RMS slope
f ≍g
two-sided comparison up to positive sizeindependent constants
local RMS slope αloc
1.0
Nicholas J. Higham and Theo Mary. “A New Approach to Probabilistic Rounding Error Analysis”. In: SIAM Journal on Scientific Computing 41.5 (2019), A2815–A2835. DOI: 10.1137/ 18M1226312.
RMS error
10
measured extracted tree pure pairwise
10−14
103
104
105
Figure 12. Measured exact-reference RMS error of contiguous binary64 np.sum, the prediction from the extracted tree using the independently calibrated ν = 0.19, and the idealized pure-pairwise prediction. Section 5.1 sketches the extracted hierarchy. Shading is the bootstrap 95% interval. Section 6.9 discusses the validation and error prediction.
0.6
independent calibration measured
101
102
103
Figure 13. Local non-FMA dot-product slopes and the curve from independently calibrated product and accumulation coefficients. The derivation appears in Appendix 1.1.
[2]
Eric Hallman and Ilse C. F. Ipsen. “Precision-aware deterministic and probabilistic error bounds for floating point summation”. In: Numerische Mathematik 155.1–2 (2023), pp. 83–119. DOI: 10.1007/s00211-023-01370-y.
[3]
Eric O. D. Andriantiana, Kenneth Dadedzi, and Stephan Wagner. “The Ancestral Matrix of a Rooted Tree”. In: Linear Algebra and its Applications 575 (2019), pp. 35–65. DOI: 10.1016/j.laa. 2019.04.004.
[4]
David A. Huffman. “A Method for the Construction of Minimum-Redundancy Codes”. In: Proceedings of the IRE 40.9 (1952), pp. 1098–1101. DOI: 10.1109/JRPROC.1952.273898.
[5]
Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. 2nd ed. Wiley-Interscience, 2006. ISBN: 978-0-47124195-9. DOI: 10.1002/047174882X.
[6]
Ming-Yang Kao and Jie Wang. “Linear-Time Approximation Algorithms for Computing Numerical Summation with Provably Small Errors”. In: SIAM Journal on Computing 29.5 (2000), pp. 1568–1576. DOI: 10.1137/S0097539798341594.
[7]
Lawrence L. Larmore and Daniel S. Hirschberg. “A Fast Algorithm for Optimal Length-Limited Huffman Codes”. In: Journal of the ACM 37.3 (1990), pp. 464–473. DOI: 10.1145/ 79147.79150.
[8]
Nicholas J. Higham. “The Accuracy of Floating Point Summation”. In: SIAM Journal on Scientific Computing 14.4 (1993), pp. 783–799. DOI: 10.1137/0914050.
[9]
Pierre Blanchard, Nicholas J. Higham, and Theo Mary. “A Class of Fast and Accurate Summation Algorithms”. In: SIAM Journal on Scientific Computing 42.3 (2020), A1541–A1557. DOI: 10.1137/19M1257780.
[10]
Anthony M. Castaldo, R. Clint Whaley, and Anthony T. Chronopoulos. “Reducing Floating Point Error in Dot Product Using the Superblock Family of Algorithms”. In: SIAM Journal on Scientific Computing 31.2 (2009), pp. 1156–1174. DOI: 10.1137/070679946.
106
k
0.7
k
10−15 102
0.8
100
Predicting np.sum from its tree
−13
0.9
0.5
References [1]
et al.
Piyush Sao
Floating-Point Reduction Trees
19
RMS slopes: 0.602 measured, 0.619 theory
[16]
Douglas M. Priest. “Algorithms for arbitrary precision floating point arithmetic”. In: Proceedings of the 10th IEEE Symposium on Computer Arithmetic. IEEE Computer Society Press, 1991, pp. 132–143. DOI: 10.1109/ARITH.1991.145549.
[17]
James Demmel and Yozo Hida. “Fast and Accurate Floating Point Summation with Application to Computational Geometry”. In: Numerical Algorithms 37.1–4 (2004), pp. 101–112. DOI: 10.1023/B:NUMA.0000049458.99541.38.
[18]
Takeshi Ogita, Siegfried M. Rump, and Shin’ichi Oishi. “Accurate Sum and Dot Product”. In: SIAM Journal on Scientific Computing 26.6 (2005), pp. 1955–1988. DOI: 10.1137/030601818.
[19]
Caroline Collange et al. “Numerical reproducibility for the parallel reduction on multi- and many-core architectures”. In: Parallel Computing 49 (2015), pp. 83–97. DOI: 10.1016/j.parco. 2015.09.001.
[20]
Thomas G. Robertazzi and Stephen C. Schwartz. “Best “Ordering” for Floating-Point Addition”. In: ACM Transactions on Mathematical Software 14.1 (1988), pp. 101–110. DOI: 10.1145/ 42288.42343.
[21]
Terje O. Espelid. “On Floating-Point Summation”. In: SIAM Review 37.4 (1995), pp. 603–607. DOI: 10.1137/1037130.
[22]
John von Neumann and Herman H. Goldstine. “Numerical inverting of matrices of high order”. In: Bulletin of the American Mathematical Society 53.11 (1947), pp. 1021–1099. DOI: 10.1090/ S0002-9904-1947-08909-6.
[23]
Peter Henrici. “Test of probabilistic models for the propagation of roundoff errors”. In: Communications of the ACM 9.6 (1966), pp. 409–410. DOI: 10.1145/365696.365698.
[24]
T. E. Hull and J. R. Swenson. “Tests of probabilistic models for propagation of roundoff errors”. In: Communications of the ACM 9.2 (1966), pp. 108–113. DOI: 10.1145/365170.365212.
[25]
Jesse L. Barlow and Erwin H. Bareiss. “On roundoff error distributions in floating point and logarithmic arithmetic”. In: Computing 34.4 (1985), pp. 325–347. DOI: 10.1007/BF02251833.
[26]
Françoise Chatelin and Marie-Christine Brunet. “A probabilistic round-off error propagation model. Application to the eigenvalue problem”. In: Reliable Numerical Computation. Ed. by M. G. Cox and S. Hammarling. Oxford University Press, 1990, pp. 139–160. DOI: 10.1093/oso/9780198535645.003.0009.
[27]
Michael P. Connolly, Nicholas J. Higham, and Theo Mary. “Stochastic Rounding and Its Probabilistic Backward Error Analysis”. In: SIAM Journal on Scientific Computing 43.1 (2021), A566–A585. DOI: 10.1137/20M1334796.
[28]
Matteo Croci et al. “Stochastic rounding: implementation, error analysis and applications”. In: Royal Society Open Science 9.3 (2022), p. 211631. DOI: 10.1098/rsos.211631.
measured
RMS error
scaled √ Λ1
10−14
Effective coefficient drift 0.24
0.217 → 0.190
νeff
0.22 0.20 0.18
103
104
k
Figure 14. Four-stage √round-to-nearest scaling experiment. Left: measured RMS and the Λ1 curve scaled at the first point; their finitegrid slopes are 0.602 and 0.619. Right: νeff = MSE/( u2 Λ1 ) with bootstrap 95% intervals. Its decrease from 0.217 to 0.190 across the grid accounts for the observed slope difference. Appendix 1.2 describes the fixed-design experiment.
[11]
[12]
Peichen Xie et al. “Revealing Floating-Point Accumulation Orders in Software/Hardware Implementations”. In: 2025 USENIX Annual Technical Conference (USENIX ATC 25). Boston, MA: USENIX Association, July 2025, pp. 1425–1440. ISBN: 978-1-939133-48-9. URL: https : / / www . usenix . org / conference/atc25/presentation/xie. NumPy contributors. NumPy 1.26.4 source: pairwise_sum in numpy/core/src/umath/loops_utils.h.src. Version 1.26.4. The 8192-element chunking arises from the ufunc reduction buffer, whose size is controlled by numpy.setbufsize, outside the pairwise-summation routine. NumPy. Feb. 6, 2024. URL: https://github.com/numpy/numpy/ blob/v1.26.4/numpy/core/src/umath/loops_utils.h.src (visited on 07/20/2026).
[13]
Paulius Micikevicius et al. “Mixed Precision Training”. In: 6th International Conference on Learning Representations (ICLR 2018). 2018. arXiv: 1710.03740 [cs.AI]. URL: https://openreview. net/forum?id=r1gs9JgRZ.
[29]
El-Mehdi El Arar et al. “Stochastic Rounding Variance and Probabilistic Bounds: A New Approach”. In: SIAM Journal on Scientific Computing 45.5 (2023), pp. C255–C275. DOI: 10.1137/ 22M1510819.
[14]
Pieter Eendebak. ENH: fast path for full contiguous reductions. NumPy pull request 31274; merged 29 April 2026; released in NumPy 2.5.0 on 21 June 2026. NumPy. Apr. 19, 2026. URL: https://github.com/numpy/numpy/pull/31274 (visited on 07/20/2026).
[30]
Ming-Yang Kao and Jie Wang. “Minimizing Roundoff Errors of Prefix Sums via Dynamic Construction of Huffman Trees”. In: Theoretical Computer Science 262.1–2 (2001), pp. 101–115. DOI: 10.1016/S0304-3975(00)00183-3.
[15]
Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms. 2nd ed. SIAM, 2002. DOI: 10.1137/1.9780898718027.
[31]
James Demmel and Hong Diep Nguyen. “Fast Reproducible Floating-Point Summation”. In: 2013 IEEE 21st Symposium on Computer Arithmetic. IEEE, 2013, pp. 163–172. DOI: 10.1109/ ARITH.2013.9.
20
Piyush Sao
[32]
Roman Iakymchuk et al. “Reproducible and Accurate Matrix Multiplication”. In: Scientific Computing, Computer Arithmetic, and Validated Numerics. Vol. 9553. Lecture Notes in Computer Science. Cham: Springer International Publishing, 2016, pp. 126–137. DOI: 10.1007/978-3-319-31769-4_11.
[33]
Samuel D. Pollard and Boyana Norris. “A Statistical Analysis of Error in MPI Reduction Operations”. In: 2020 IEEE/ACM 4th International Workshop on Software Correctness for HPC Applications (Correctness). 2020, pp. 49–57. DOI: 10 . 1109 / CORRECTNESS51934.2020.00011.
[34]
Pierre Blanchard et al. “Mixed Precision Block Fused MultiplyAdd: Error Analysis and Application to GPU Tensor Cores”. In: SIAM Journal on Scientific Computing 42.3 (2020), pp. C124– C141. DOI: 10.1137/19M1289546.
[35]
Suyog Gupta et al. “Deep Learning with Limited Numerical Precision”. In: Proceedings of the 32nd International Conference on Machine Learning. Vol. 37. Proceedings of Machine Learning Research. 2015, pp. 1737–1746. URL: https://proceedings.mlr. press/v37/gupta15.html.
[36]
El-Mehdi El Arar et al. Probabilistic Error Analysis of LimitedPrecision Stochastic Rounding: Horner’s Algorithm and Pairwise Summation. Mar. 25, 2026. DOI: 10.48550/arXiv.2603.24161. arXiv: 2603.24161 [math.NA]. URL: https://arxiv.org/abs/2603. 24161 (visited on 07/20/2026).
et al.