Sketching the Error, Not the Product: Post Hoc Fault Recovery for Half Precision GPU Matrix Multiplication Pranav Napolean
arXiv:2609.19758v1 [cs.DC] 17 Sep 2026
Department of Computer Science and Engineering National Institute of Technology Warangal Warangal, India [email protected]
Vikas Srivastava
Napolean Periathambi
Department of Mathematics Executive Director National Institute of Technology Warangal AthenaHealth Warangal, India [email protected] [email protected]
Abstract. Silent data corruption (SDC) from defective accelerators now interrupts large scale training, yet deployed mitigations act on whole nodes. Algorithm based fault tolerance (ABFT) for a single GEMM has to be fused into the kernel or encode the operands, and it localizes at most one error per checksum. We present FP-Sketch, a verifier that runs after an unmodified tensor core GEMM whose half precision operands are accumulated and delivered at FP32. A sum sketch detects corruption on every call. Hashed first moment sketches, confirmed by independent recomputation, then localize several corrupted entries with no false positives by construction, and each fault yields a coordinate and a magnitude for fleet diagnosis. In floating point, sketch noise rather than bucket collisions limits localization. We measure that noise and find that its constant depends on the BLAS and the operand format and that the bucket count must grow as n2.57 for a square product. Sizing the bucket count by measured noise rather than by a fitted power of n raises recovery on eight transformer shapes from 0.402 to 1.000, and measuring the noise at run time adapts the bucket count to the kernel and the model. Instruction level injection with NVBit shows that upsets in a live accumulator are often only 2 to 9% of a typical entry, a population that output side injection cannot produce. Output side injection recovers every fault, while under NVBit the same engine sized for faults of typical magnitude recovers 0.550, and sizing for the measured magnitudes restores 1.000. On Llama-2-7B, guarding the MLP down projections removes 99.4% (BF16) and 99.9% (FP16) of the perplexity damage caused by 2048 bit flips, and the clean path probe costs 0.78 to 3.06 ms against GEMMs of 0.35 to 12.47 ms. Index Terms. Silent data corruption, algorithm based fault tolerance, matrix multiplication, mixed precision, fault injection, sketching
I. I NTRODUCTION Silent data corruption is a compute fault that returns a valid but wrong value with no exception, NaN or crash. It has become a first order reliability problem for large scale AI [1], [2]. Training multiplies the exposure. A single run occupies tens of thousands of accelerators for weeks of synchronous work, and since every step builds on the state the previous one produced, one wrong value is carried into every later step [3], [4]. Reported incidents understate the damage. Instruction level injection into the matrix multiply instructions of BF16 LLaMA training shows exponent bit upsets raising evaluation loss [4],
and permanent GPU faults can distort parameters with no NaN, no Inf, and no loss anomaly [5]. The most damaging corruption is thus the kind that is never counted. Production mitigations work at the level of a node: they detect a misbehaving node, evict it, and roll back [3]. They do not say which product was wrong, or where. Below the node, ABFT embeds checksums in the matrix product [6]–[8]. GPU instances fuse the checksum into a custom GEMM kernel [9], [10] or encode the operands [11], [12], and the position weighted form localizes one error per checksum [7], [12]. Mixed precision training makes all of this hard to deploy. The GEMMs that dominate training run in closed vendor libraries that a fault tolerance layer cannot modify, and a checksum over output stored in BF16 is swamped by storage rounding. ABFT forced to BF16 raised false positives on healthy nodes [3]. V-ABFT shows that verifying inside a fused kernel, before quantization, restores FP32 level sensitivity [12]. Our question is how to obtain that sensitivity, together with entry level localization of several simultaneous errors, without owning the kernel. FP-Sketch verifies the product after the fact. The unmodified GEMM delivers its FP32 accumulator instead of a BF16 store, the verifier reads A, B and C, and the caller narrows C afterwards. A sum sketch of the error E = AB − C, computed without forming AB, decides on every call whether C is wrong. When it fires, first moment sketches weighted by row and column index recover the coordinates of each isolated error as a ratio. This is the weighted checksum device [7], applied inside hash buckets so that several errors separate. A candidate is accepted only if an independent recomputation of that entry disagrees with C. Two effects that do not exist over exact rings decide whether this works on real hardware. First, the ratio must resolve an index among thousands of rows, so the rounding noise of the sketch, not bucket collisions, sets the bucket count, and that noise depends on the summation order of the library. Second, a fault in a live accumulator is buried by the additions that follow it, and hardware faults are often far smaller than the flipped bit suggests. The standard evaluation method, writing
a flipped value into the finished product, never produces them. V-ABFT makes the same observation for a checksum fused into the kernel [12]. Here the same point is reached without Contributions. modifying the kernel. • A design for ex post facto verification of half precision tensor core GEMMs (§III). A cheap S only probe B. Threat Model runs on every call, while localization and repair run In scope are compute faults in a processing element, FMA, or only on the rare dirty call. Localization combines FP32 accumulator that leave the stored operands intact and produce moment sketches with scaled index weights, peeling across wrong entries in C, with ∥E∥0 = s ≪ n1 n3 . The evaluation hash rounds, a bucket count and neighbourhood search covers single and multiple bit upsets, stuck at lines, whole sized from the noise measured on each product, and an word corruption, and nonfinite results. exact scan for nonfinite entries (Algorithms 1 and 2). Out of scope is corruption already present in A or B, since False positives are impossible under a trusted recompute a recomputation reads the same wrong operand and ECC and assumption (§IV). checkpoint validation apply there. Faults outside the matrix • A measured noise law for hashed moment localization product are also out of scope, as are products that are never (§III.F). On cuBLAS it bounds all 65 H100 measurements, materialized. Fused attention kernels never store QK ⊤ , so it makes the bucket count scale as n2.57 , and assigns FP16 cannot be verified after the fact. Within a training step, an operands 1.76 times the BF16 constant. Its n2 exponent undetected fault in the output of layer k becomes the input differs between libraries. Sizing by measured noise rather of layer k+1. Checking every guarded GEMM stops the fault than by a fitted power of n raises recovery on eight before it is laundered into legitimate data. transformer shapes from 0.402 to 1.000. Trusted recomputation. The guarantee of zero false • A fault magnitude result (§VI.C). Measured against NVBit positives (Lemma 1) requires that the recomputation validating accumulator injection, output side injection overstates a candidate is not corrupted by the same fault. Our localization recovery. Sized for the measured magnitudes, implementation recomputes at FP32 on the same device, so the engine recovers 1.000 (held out campaigns: 0.975 a permanent fault that corrupts both the product and its with BF16 operands, 1.000 with FP16), and the flipped recomputation identically is not covered. bit position cannot be recovered (0/80). C. Related Work • An H100 evaluation (§VI) with instruction level injection across single bit, multiple bit, adjacent burst, stuck at and ABFT. Huang and Abraham introduced checksum based whole word upsets, real models (including accumulator ABFT for matrix operations [6], and the position weighted faults on a real Llama-2-7B layer, recovered at 0.975 checksum of Jou and Abraham recovers the index of a to 1.000), per shape deployable costs, and baselines on single error from a ratio of syndromes [7]. Gunnels et al. identical faults. brought detection and rollback to high performance matrix multiplication [8], and Ding et al. located and corrected errors online on GPUs [9]. FT-GEMM fuses ABFT into a custom GPU A. Mixed Precision GEMM and Where to Verify SGEMM kernel at 8.89% average overhead over cuBLAS [10]. For A ∈ Rn1 ×n2 and B ∈ Rn2 ×n3 , a tensor core GEMM A-ABFT and V-ABFT derive tighter detection thresholds [12], loads FP16 or BF16 operands, accumulates each partial [13], ApproxABFT loosens them to cut recomputation [14], sum in an FP32 register, and by default stores C back at ATTNChecker specializes to attention and extreme values [15], the operand precision. The instrumented instruction stream and grid like codes correct errors confined to at most two rows confirms this: FP16 and BF16 operands both emit an HMMA and columns by encoding the operands [11]. Table I contrasts whose destination is an FP32 partial sum (HMMA.16816.F32, these approaches. Verifying and correcting products. Freivalds’ test decides HMMA.*.F32.BF16). The store discards sixteen mantissa bits, and reading them whether C = AB in O(n2 ) [16], and coding theoretic tests ˛ et al. correct s wrong back returns zeros, not data. A verifier shown the stored detect sparse errors [17]. Gasieniec value sees a number from which the evidence has already entries after the fact over a ring by sketching the error with been removed, whatever its threshold. We therefore fix the Pagh’s compressed multiplication [18], [19]. Roche and Wu deployment point: the product is verified at accumulation and Wang give exact algorithms over fields and integers [20], precision, before the store narrows it. The GEMM is asked [21]. Hashed sketches [22], [23] and invertible Bloom lookup for an FP32 output (cuBLASLt compute type 32F with an tables [24] supply the bucket structure. All of these assume FP32 result, or a kernel that stores the accumulator it already exact arithmetic, where the sizing problem driven by noise holds), and the caller narrows C after verification. What this (§III.F) does not arise. costs is memory traffic. Delivering C at FP32 and narrowing SDC characterization and injection. Field and injection it afterwards adds 6n1 n3 bytes to a product of 2n1 n2 n3 flops, studies characterize SDC in LLM training and GPU a fraction 3π/(βn2 ) of the GEMM for achieved tensor core kernels [3]–[5], [25]–[27]. Error injection at a high level is throughput π and bandwidth β. On an H100 that is about known to misestimate resilience [28], and single bit flips are a 443/n2 : 11% at n2 = 4096, 3.1% at 14336, and 58% at 768. minority of gate level GPU faults [25]. Section VI.C measures II. BACKGROUND , T HREAT M ODEL , AND R ELATED W ORK
Unmodified GEMM half prec. in, FP32 out
C
Probe sum sketch S
clean
Narrow and store C
dirty Fleet diagnosis
faults
Localize S, R, T , recompute
Apply repair, record
Fig. 1. FP-Sketch runs after the GEMM and before the store. The probe runs on every call, and localization and repair run only when it fires.
one specific consequence of this for localizers. Floating point reproducibility. Summation order, and with it rounding, depends on the parallel decomposition of a library [29], and the rounding noise of GPU matmul is structured [30]. Section III.F quantifies what that noise does to localization. III. D ESIGN A. Overview Figure 1 shows the three stages. The probe answers whether C is wrong and runs on every call. Localize answers where and by how much, and it does not modify C. Apply writes the repairs and emits one record per fault. We keep the stages as separate calls because their costs differ by orders of magnitude (§VI.E). A deployment that wants only a correct product can recompute on detection and still localize for telemetry. B. Moment Sketches Draw hashes h1 : [n1 ] → [m] and h2 : [n3 ] → [m] together with Rademacher signs v1 , v2 ∈ {±1} [23]. Let H1 ∈ Rm×n1 hold v1 (i) at (h1 (i), i), and define H2 in the same way. The sum sketch S = H1 EH2⊤ = (H1 A)(BH2⊤ ) − H1 CH2⊤
(1)
aggregates the error in each bucket and is formed without AB. Its m × n2 and n2 × m factors come from scatter adds over A and B, and their product costs m2 n2 . The first moment sketches R and T repeat (1) with the rows of A and C weighted by wi = i/2⌈log2 n1 ⌉ and the columns of B and C weighted by wj = j/2⌈log2 n3 ⌉ , with i and j counted from 1. If bucket (a, b) holds exactly one error at (i, j), then i = 2⌈log2 n1 ⌉ Rab /Sab ,
j = 2⌈log2 n3 ⌉ Tab /Sab .
(2)
The weights are scaled by a power of two so that no term of R is larger than the corresponding term of S. Without the scaling, an error with |Eij | > fp32max/n1 would overflow R while S stayed finite. Dividing by a power of two only shifts an exponent, so every rounding in the accumulation is unchanged. C. The Cheap S Only Probe Production GEMMs are clean on the overwhelming majority of calls, so the stage that decides deployed cost is the one that runs on a clean product. Deciding whether C is wrong needs only S. If no bucket of S exceeds the probe threshold, no error large enough to localize is present, and R and T need not be formed. The probe therefore builds S alone from a scatter over
the rows of A, a scatter over the rows of C folded over its columns, and one m × n2 × m product with the factor BH2⊤ of the B side, which is cached because weights stay the same across calls. It forms three of the seven accumulated outputs that a localization round needs, reads A once and C once, and recomputes no entry. Its output is a single bit, thresholded as in §III.E together with an explicit finiteness test (§III.I). That bit can stay on the device until the caller reads it once per step. The probe costs O(n2 ) (Theorem 3) and never walks the neighbourhood search. Its measured cost appears in §VI.E. D. Validation A bucket holding several errors yields a meaningless ratio, so every candidate (i, j) is validated. The inner product γ = Ai,∗ B∗,j is recomputed with FP32 accumulation, and the candidate P is accepted only if |γ − Cij | > τij , where τij = c ε k |Aik ||Bkj | bounds the rounding of that inner product [31], c = 100, and ε = 2−23 is the FP32 machine epsilon. The accepted correction is δij = γ − Cij . E. Thresholds The analytic threshold τ = 100 n2 ε max |Aik | max |Bkj | ik
kj
bounds the worst case accumulation noise in a bucket. Probe and localizer threshold the same kind of sketch but run different tests. The output of the probe is a single bit with no filter after it, so the threshold must√clear the maximum over m2 clean buckets: thr = min(τ, µ 2 ln m2 σ̂) with σ̂ = 1.2533 · mean(min(|S|, 5 mean|S|)) and margin µ = 4. This clipped mean needs no sort, so the estimate stays on the device. The localizer selects candidates, and a false candidate costs one recomputation, so its threshold may be looser. It uses τc = min(τ, k σ̂MAD ) with k = 2, where σ̂MAD = 1.4826 · median |S| − median|S| estimates the clean bucket noise from the median absolute deviation of |S|. About 4.6% of clean buckets clear τc by chance. Validation discards all of them, but each one would still walk the whole neighbourhood search. A lone fault of size δ decodes its row to within about e max(n1 , n3 ) σ/δ, where e is the normalized index error (median 0.58, 95th percentile 2.5). A bucket can therefore localize within radius r only if its signal is at least of order max(n1 , n3 ) σ/r. The localizer adds this as a floor, τd = max(n1 , n3 ) σ̂MAD / max(r, 12 ), keeps the buckets with |Sab | > max(τc , τd ), and retains at most K = max(64, 8s) of them, ranked by |Sab |. F. Sketch Noise Sets the Bucket Count Over a ring, √ collisions alone set m. The condition m ≥ max(3c∆, 3cs) isolates an error with probability at least 1 − 1/c, where ∆ bounds the errors per row or column (Lemma 3). In floating point, (1) is the difference of two independently rounded accumulations, and its noise σ grows with the number of entries per bucket. On an H100 (cuBLAS, BF16 operands, FP32 accumulation), normalizing each measured σ √ by 2−24 n1 n3 n2 rms(C)/m gives 0.1726 to 0.1749 across
TABLE I FAULT TOLERANCE FOR A SINGLE GEMM. Vendor GEMM ASKS WHETHER THE METHOD CAN PROTECT AN UNMODIFIED LIBRARY CALL . N / S : NOT STATED IN THE CITED SOURCE . Method
Redundancy
Checksum ABFT [6] Weighted checksum [7] FT-GEMM [10] V-ABFT [12] Grid like ECC [11] Gasieniec ˛ et al. [18] This work
row and column checksums any 1 post hoc form position weighted checksums any 1 per checksum post hoc form fused into custom kernel FP32 n/s no encoded, fused for FP32 threshold BF16 to FP64 1 per row no encoded, enlarged product real ≤ 2 rows, columns no hashed sketch exact ring s yes hashed moment sketches BF16/FP16 s (sized) yes (FP32 output)
Operands
Errors localized
Vendor GEMM
Reported cost O(n2 ) O(n2 ) 8.89% over cuBLAS n/s 1.24 to 1.37× latency Õ(n2 + sn) probe 0.78 to 3.06 ms
n3 ∈ [1024, 16384] and 0.1717 to 0.1783 across m ∈ [32, 512]. G. Fault Magnitude, Calibration, and the Neighbourhood Across n2 ∈ [1024, 65536], however, the same normalized Search value climbs from 0.162 to 0.219. The n2 exponent is 0.57 Equation (4) is written for a fault the size of a typical entry, rather than 1/2, and with it the normalized value stays flat to ρ ≡ δ/rms(C) = 1. Faults in a live accumulator can be much within 1.10× over that sweep. FP16 operands are noisier than smaller (§VI.C), which multiplies B by 1/ρ. A larger m and a BF16 at every one of 19 identical shapes, by 1.20 to 1.76× wider validation window both absorb that factor. Validating a with a median of 1.47. Fitting all 65 H100 measurements (43 (2r+1)2 neighbourhood of the rounded index costs that many BF16, 22 FP16) as an upper envelope gives inner products per candidate, and on the GPU this dominates: 0.07 −24 √ n1 n3 n2 (n2 /8192) rms(C)/m (3) at the Llama-2-7B down_proj shape (4096×11008×4096), σ ≲ cf 2 localization with radius 16 took 1153 to 1196 ms at every m with cf = 0.19 for BF16 and cf = 0.33 (1.76 × 0.19) from 112 to 960, while radius 4 took 96 to 97 ms. FP-Sketch for FP16. No measurement exceeds (3). The earlier fit therefore declares the smallest magnitude it must localize, ρmin √ 0.17 · 2−24 n1 n3 n2 rms(C)/m, with no format term, was (default 0.02), sizes the bucket count for a target radius r⋆ = 4, exceeded by 17 of the 43 BF16 measurements, by up to 1.29×. and lets the radius grow toward rmax = 16 only when m Rounding in (2) perturbs the recovered row by (νR −i νS )/δ, cannot grow further. Calibration. Equation (3) is fitted on one GPU, one library where νR and νS are the bucket noises and δ = |Eij |. The scatter is Θ(max(n1 , n3 ) σ/δ) with zero bias, and we call and Gaussian operands. Other kernels and real activations B = max(n1 , n3 ) σ/δ the index budget. Keeping B below a change the noise. On a Llama-2-7B down_proj layer, products from the tiled kernel used for injection carried about target B ⋆ requires 8× the noise of cuBLAS products on the same operands. The max(n1 , n3 ) −24 √ plan for a dirty call is therefore computed from the noise mnum ≳ c 2 n n n f 1 3 2 (4) B⋆ δ measured on that call. Let σ̂ be the probe estimate of §III.E if 0.07 · (n2 /8192) rms(C) the probe has just run on this C, and otherwise the estimate ms(C) be the rms of the finite which is Θ(n2.57 ) for a square product. FP-Sketch sizes for σ̂MAD of a freshly built S. Let rd δ = rms(C) with B ⋆ √ = 0.3. With the collision requirement entries of C taken at a stride of max(1, ⌊n1 n3 /65536⌋), with magnitudes clipped at 64 times their median. Then mcomb = ⌈max(3c∆, 3cs)⌉ at c = 2, it starts from max(n1 , n3 ) σ̂ m = min max(mcomb , mnum , 16), B= , (5) ρmin rd ms(C) max(mcomb , mmax ) ( 0 B ≤ 0.05, p rB = The cap mmax = max(16, ⌊ n1 n3 /16⌋) keeps an average max(2, ⌈2.5 B⌉) otherwise. of 16 entries in each bucket. Below that population the noise estimates collapse and clean products read as dirty, which we The factor 2.5 covers the 95th percentile of measured index error. If rB ≤ r⋆ , localization keeps mloc = m and uses observed at m = 966 on a 1024 × 1024 output. The constant belongs to the BLAS, not to the algorithm. r = rB . Otherwise it builds its sketches at σ measures the disagreement between two summation orders, mloc = max m, min(16⌈⌈m rB /r⋆ ⌉/16⌉, mmax , mmem ) , and the blocking of a library sets how far it departs from sequential summation [29]. The n2 exponent is +0.587 for where mmem is the largest m with 4(3m2 +6mn2 ) bytes within an explicitly sequential sum, +0.200 through a blocked CPU 2 GiB. The radius at that size follows from B ′ = B m/mloc matmul (oneDNN) and +0.133 for the sketch built on it, against as r = 0 if B ′ ≤ 0.05 and r = min(rmax , max(2, ⌈2.5 B ′ ⌉)) +0.577 for the sketch on cuBLAS. A rule fitted on one library otherwise. If σ̂ or rd ms(C) is not a positive finite number, the and deployed on the other underprovisions m by 4.4× at plan falls back to the same rule with σ taken from (3). n2 = 8192. For this reason FP-Sketch uses (3) only as a Recalibrating the probe. The probe keeps its own m, which starting point and measures σ at run time (§III.G). sets its detection floor. It averages σ̂ over the calls it reports
Algorithm 1 One round of floating point localization Require: A, B, C (FP32), round (h1 , v1 , h2 , v2 ) at mloc , radius r, confirmed set F with absolute entries Fabs ⊆ F 1: S, R, T ← sketches of E = AB − C with scaled weights (§III.B) ′ 2: for all (i, j, δ) ∈ F \ Fabs do ▷ peel across rounds When m lies outside [m/1.25, 1.25 m], the engine rebuilds ′ 3: subtract v (i)v (j) δ · (1, w , w ) from (S, R, T ) at 1 2 i j its hashes, buffers, caches and graphs at m . The formula (h (i), h (j)) 1 2 therefore matters only until the first measurement. After repair, 4: end for FP-Sketch probes the product once more with a fresh hash 5: σ̂ ← MAD estimate of S round and recomputes it if the probe still fires, so a missed 6: τc ← min(τ, 2σ̂), τd ← max(n1 , n3 ) σ̂/ max(r, 12 ) fault that the probe can still see does not reach the output. 7: Q ← the K largest |Sab | with |Sab | > max(τc , τd ) None of these steps runs on a clean product. 8: U ← {(a, b) ∈ Q : Sab ̸= 0 and Sab , Rab , Tab finite} H. Peeling Across Hash Rounds 9: for all (a, b) ∈ U do (ı̂ab , ȷ̂ab ) ← (2), rounded to integers A fault occupies exactly one bucket per hash round, so 10: subtracting a recovered fault inside the round in which it was 11: end for found unlocks nothing. Across rounds it does. FP-Sketch draws 12: for all offsets (di , dj ) with |di |, |dj | ≤ r, by Chebyshev then Manhattan distance do k = 3 independent hash rounds and, before selecting candidates for all (a, b) ∈ U , in a fixed order do in a round, subtracts every finite fault confirmed so far from 13: (i, j) ← (ı̂ab + di , ȷ̂ab + dj ) the S, R and T of that round. A bucket that would still hold 14: if (i, j) lies outside C or (i, j) ∈ F then continue a mixture then becomes a singleton. Only validated faults are 15: end if subtracted, and every new candidate is validated again, so 16: 17: γ ← Ai,∗ B∗,j ▷ FP32 accumulation peeling cannot create a false correction. 18: if |γ − Cij | > τij then I. Nonfinite and Out of Range Entries 19: F ← F ∪ {(i, j, γ − Cij )}, U ← U \ {(a, b)} An infinite or NaN entry is contagious in a linear sketch. 20: end if One inf in the row scatter poisons a bucket row, and every 21: end for fault sharing that row becomes unrecoverable. FP-Sketch finds 22: end for such entries, and any entry with |Cij | > fp32max/16, with an 23: return F exact elementwise scan and keeps at most K of them. Each is recomputed as an absolute value rather than a delta, because NaN + δ is NaN, and C holds the recomputed value while Theorem 1 (Floating point recovery). Let error (i, j) be the sketches are built. These entries are not subtracted when isolated in bucket (a, b) with |Sab | > max(τc , τd ) and among peeling, since their error is already absent from the sketch. the K largest such buckets, and let the rounded index of (2) lie The probe tests the finiteness of S explicitly: (|S| > thr ) is within r of (i, j) in each coordinate. Then Algorithm 1 returns false for a NaN bucket, so a probe that only compared values (i, j, δ) with |δ − Eij | ≤ τij and returns no uncorrupted entry. would report such a product clean. Since the index scatter is Θ(B) (§III.F), the index condition holds with high probability once r ≥ 2.5B. IV. A NALYSIS Lemma 1 (Validation soundness). If the recomputation of γ Proof. The search reaches (i, j), and validation accepts it is not affected by the fault and τij bounds its FP32 rounding, because |Eij | exceeds the rounding bound. Lemma 1 excludes every clean position visited. every accepted (i, j) is a corrupted entry. clean, without a host synchronization, and after 16 calls and then every 1024 it computes B1 = max(n1 , n3 ) σ̄/d rms(C) and m′ = max mcomb , 16, min(⌈m B1 /B ⋆ ⌉, mmax ) .
Proof. By assumption |γ − (AB)ij | ≤ τij , so |γ − Cij | > τij implies Cij ̸= (AB)ij . Lemma 2 (Collision). For a fixed error (i, j) under independent uniform hashes, Pr[(i, j) not isolated] ≤ 2(∆ − 1)/m + s/m2 . Proof. An error in the same row or column collides with probability 1/m. Any other error collides only if both hashes agree, which happens with probability 1/m2 . A union bound completes the argument. Lemma 3 (Isolation per round). If m ≥ 3c∆ and m2 ≥ 3cs with c > 1, each error is isolated in a round with probability at least 1 − 1/c.
Theorem 2 (Recovery over rounds). Under Lemma 3 and the index condition of Theorem 1, k independent rounds recover all s errors with probability at least 1 − s c−k . Proof. An error is missed in all rounds with probability at most c−k , and a union bound over the s errors gives the result. Remark 1 (Peeling). Peeling across rounds subtracts only validated faults, so a bucket isolated for an error in the independent process stays isolated, up to the rounding residual of each subtracted δ and the candidate cap K. A proof that covers both effects is open. We use Theorem 2 as the bound, and the evaluation measures the implemented decoder.
Algorithm 2 Verify, localize, and apply Require: A, B, C, probe size m, rounds k, ρmin , r⋆ , rmax , K S only probe (every GEMM, §III.C) 1: reject C unless it is FP32 ▷ §II.A 2: S ← sum sketch (1) at m ▷ B side cached 3: σ̂ ← 1.2533 · mean(min(|S|, 5 mean|S|)) √ 4: thr ← min(τ, µ 2 ln m2 σ̂) 5: if max |S| ≤ thr and S finite then 6: add σ̂ to the recalibration average (§III.G) 7: return CLEAN 8: end if Localize (dirty GEMM only) 9: Fabs ← {(i, j, Ai,∗ B∗,j ) : Cij nonfinite or |Cij | > fp32max/16}, at most K 10: F ← Fabs , and set those Cij to their recomputed values 11: (mloc , r) ← plan from σ̂ and rd ms(C) ▷ §III.G 12: for t = 1, . . . , k do 13: F ← Algorithm 1 with a fresh round at mloc , r, F 14: end for 15: restore the entries of Fabs in C Apply 16: for all (i, j, δ) ∈ F do 17: Cij ← δ if (i, j, δ) ∈ Fabs , else Cij ← Cij + δ 18: emit record (i, j, δ) 19: end for 20: if the probe fires on C with a fresh round then C ← AB at FP32 21: end if
TABLE II R ECOVERY ON TRANSFORMER LAYER SHAPES , ONE OUTPUT SIDE FAULT PER TRIAL , 60 TRIALS PER SHAPE AND FORMAT. Fitted: BF16, BUCKET COUNT FROM A RULE FITTED AS n3/2 , LOCALIZATION AT THAT m WITH EXACT ROUNDING AND NO MEASURED PLAN . Law: THE CALIBRATED ENGINE , WITH IDENTICAL RECOVERY IN BF16 AND FP16. m IS THE PROBE ’ S BUCKET COUNT. FPR IS 0 IN EVERY CELL . Fitted n1 ×n2 ×n3
m
4096×4096×4096 48 8192×4096×4096 48 8192×4096×14336 105 8192×14336×4096 48 16384×8192×8192 128 32768×4096×4096 363 4096×4096×32768 363 4096×32768×4096 48 Pooled
Law Rec. m BF16 m FP16
Rec.
0.750 0.367 0.183 0.383 0.383 0.333 0.400 0.417
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
48 110 358 224 649 874 874 127
0.402 [0.359, 0.447]
68 193 630 393 1142 1538 1538 223
1.000 [0.992, 1.000]
the measured sketch noise floor 419×. Code and raw logs will be released. VI. E VALUATION A. Methodology
All measurements use an NVIDIA H100 80 GB (driver 580.95.05, CUDA 13.0, PyTorch 2.13) with TF32 disabled. Operands are BF16 or FP16 as stated, and products are accumulated and delivered at FP32. Output side faults flip IEEE 754 bit 26 (256×) of an entry of the finished C. Instruction level faults use NVBit [32]. A tool instruments every HMMA in a rectangular BF16 WMMA GEMM with 163 tiles and Theorem 3 (Cost). With the sketch of the B side cached, corrupts a live FP32 accumulator register partway through the probe performs O(n1 n2 + n1 n3 + m2 n2 ) operations, the multiply. A flip of bit 26 on an accumulator holding 16 2 produced exactly 4096. We use a custom NVBit tool rather and localization performs O k(n1 n2 + n1 n3 + mloc n2 ) + 2 than NVBitFI [33], which was released for NVBit 1.5.5 and kK(2r+1) n2 , against O(n1 n2 n3 ) for recomputing AB. CUDA 11.2 before the Hopper architecture, so that we can The repaired entry is a fresh FP32 inner product. It is accurate target the FP32 accumulator of HMMA instructions with up to the rounding of that recomputation rather than equal to to four simultaneous sites, arbitrary bit masks, and stuck (AB)ij . at modes. Ground truth comes from differencing against an uninstrumented run, not from the report of the injector. A trial V. I MPLEMENTATION whose difference is nonfinite has no defined ground truth and FP-Sketch is a PyTorch engine. The row and column scatters is not scored. Recovery is the fraction of injected faults located run as segmented kernels without atomic operations, compiled and repaired, reported with 95% Wilson intervals. Every result uses the calibrated engine of §III.G with k = 3, at run time through NVRTC with native BF16 and FP16 ⋆ 2 ρ variants, so the O(n ) passes never materialize a weighted min = 0.02, r = 4, rmax = 16, and τc = min(τ, 2σ̂MAD ), unless stated otherwise. The only exception is the fitted baseline temporary. The sketch of the B side is cached and keyed of Table II, which is described with it. Figure 3 and the timings on tensor identity and version, so an optimizer update in of the dirty path in §VI.E were measured with the GPU to place invalidates it. Activations change on every call and themselves. The timings reported with Tables IV and VII were are not cached. The probe can be replayed from a CUDA measured while NVBit campaigns shared the same GPU. graph. The graph writes its noise estimate to a static buffer that recalibration reads, and it can accumulate the verdict into a flag on the device, so a caller synchronizes once per B. RQ1: Accuracy on Deployable Shapes step rather than once per GEMM. probe, localize and Table II shows that sizing by measured noise is necessary apply_corrections are separate calls, and a narrowed C at deployable shapes. The fitted rule recovers 193/480, and is rejected rather than tolerated. TF32 must be disabled for the none of its 287 misses is a detection failure. Every one was contractions of the verifier itself, because enabling it raised detected and then mislocalized, which is the Θ(n) failure
Recovery
1 0.9 0.8 0.7 0.6 0.5 0.4
TABLE III R ECOVERY UNDER NVB IT INJECTION BY FAULT PATTERN , 4096×2048×4096, TWO SITES PER TRIAL , ρmin = 0.02. Above: BF16 RECOVERY AMONG FAULTS WHOSE ERROR EXCEEDS THE ROUNDING BOUND τij OF THE ENTRY. B ELOW THAT BOUND NO VERIFIER CAN SEPARATE A FAULT FROM ROUNDING . 38 TO 40 SCORED TRIALS PER CELL . Stuck at 1: ONE BF16 TRIAL OF 40 AND NO FP16 TRIAL PRODUCED AN OBSERVABLE CORRUPTION .
BF16 FP16 1.00
0.20
0.05
0.02
ρmin (smallest fault localized, ×rms(C)) Fig. 2. Recovery under NVBit accumulator injection with the calibrated engine, 4096×2048×4096, 80 faults per point (40 trials, two sites each), Wilson 95% intervals. FP16 points are offset slightly along the axis for legibility. There were no detection failures and no false positives. Every miss was a localization failure except one FP16 fault at ρmin = 1.00, which lay below its rounding bound. The same engine recovers every output side fault (Table II).
of the index budget described in §III.F. The worst shape is 8192×4096×14336 at 0.183. The calibrated engine recovers all 480 faults in each format, so 960 of 960 in total, with no false positives. The law also exposes a scaling limit. At 655363 it asks for m = 48010, which would need 25.8 GB for S, R and T , and the population cap stops m at 16384 (§VII). C. RQ2: Hardware Faults and Fault Magnitude
Fault pattern
BF16 Above FP16 FPR
1 bit, bits 26 and 27 1 bit, anywhere 2 random bits 3 random bits 4 random bits 2 adjacent bits [33] 3 adjacent bits 2 exponent bits Stuck at 0 Stuck at 1 Whole word replacement
0.975 0.423 0.671 0.775 0.936 0.440 0.380 0.987 1.000 1/1 0.988
0.975 0.767 0.841 0.849 0.948 0.786 0.682 0.987 1.000 1/1 0.988
0.975 0.382 0.662 0.787 0.962 0.372 0.449 0.987 1.000 silent 1.000
0 0 0 0 0 0 0 0 0 0 0
TABLE IV P ERPLEXITY ON WIKITEXT-2 (4096 TOKENS ), 8 OUTPUT SIDE FAULTS ON BIT 26 PER GUARDED GEMM PER CALL . G UARDED : ALL 32 L LAMA -2-7B D O W N _ P R O J AND ALL 36 GPT-2 LARGE C _ P R O J LAYERS , FP32 DELIVERY. E ACH MODEL IS HOSTED ENTIRELY IN THE STATED OPERAND FORMAT. Model
Op.
Clean Unguarded Guarded Recovered FP
Llama-2-7B BF16 7.4629 7.6965 7.4614 2047/2048 0 Instruction level injection produces a fault population that Llama-2-7B FP16 7.4650 7.6132 7.4651 2048/2048 0 output side injection cannot. NVBit corrupts the accumulator GPT-2 large BF16 26.2301 26.7666 26.2301 2303/2304 0 at an arbitrary point of the n2 loop, so an early hit perturbs GPT-2 large FP16 26.2223 26.6885 26.2223 2304/2304 0 a partial sum that is still small. The resulting |Eij | span two orders of magnitude, and the smallest are 2 to 9% of a typical entry. A value written into a finished C is always O(1) relative bound separately so that these are not counted as losses. Flips to that entry. that reach high exponent bits produce infinite, NaN or out of Figure 2 shows what follows. Declaring ρmin = 1, which range values, which the exact scan of §III.I recovers without sizes the plan for faults of typical magnitude, gives a recovery a sketch. Everything in between is localized by Algorithm 1, of 0.550 with BF16 operands and 0.487 with FP16, although subject to the index budget of §III.G. the same engine recovers every output side fault in Table II. Table III tests one to four random bits anywhere in the word, The plan sets the bucket count, the radius and the floor τd adjacent bursts including the double bit model of NVBitFI [33], for that magnitude, and the smaller accumulator faults do not two exponent bits, stuck at lines and whole word replacement fit it. Declaring the measured magnitude, ρmin = 0.02, raises under instruction level injection. Faults above the rounding recovery to 1.000 and 0.963. Recovery figures obtained by bound are recovered at 0.682 to 1.000 with BF16 operands writing into a completed product are therefore upper bounds and 0.690 to 1.000 with FP16, and no pattern produced a false for a localizer. Because ρmin was chosen with this campaign positive. Recovery rises with the number of random bits flipped, in mind, the held out campaigns on fresh fault sites are the from 0.423 for one bit anywhere to 0.936 for four, because ones that validate it. They recover 78 of 80 faults with BF16 more flipped bits yield larger errors. Adjacent bursts are the operands (0.975, [0.913, 0.993]) and 80 of 80 with FP16 (1.000, hardest pattern, at 0.440 and 0.380, since a burst confined [0.954, 1.000]), with no false positives, and every one of the to low mantissa bits often stays below the bound. Two bit 80 products was delivered correct after repair and the second exponent upsets are recovered at 0.987 in both formats. Stuck probe. at 1 lines are almost always silent. One BF16 trial in 40 Upsets of several bits. Repair does not depend on how produced an observable, finite error, which was recovered, and many bits flipped. A localized entry is replaced by a fresh no FP16 trial produced one. For values of normal magnitude, recomputation, so a two bit, stuck at or whole word corruption the exponent bits such a line pins are already set, which is the is repaired exactly as a single flip is. The number of flipped bits data dependence that lets stuck at faults survive in the field. matters only through the magnitude of the error it produces, and each regime is handled by a different mechanism. Flips D. RQ3: Real Models confined to low mantissa bits can leave an error below the Guarding every MLP down projection of Llama-2-7B rounding bound of the entry, where no verifier can separate the (Table IV) recovers 2047 of 2048 faults with BF16 operands fault from rounding, and Table III reports recovery above the and all 2048 with FP16, and every guarded call is flagged. The
Operands, fault
BF16 FP16
Llama, 1 bit 0.988 0.988 Llama, 2 exponent bits 1.000 0.975 Synthetic, 1 bit 0.988 0.988
TABLE VI T RAINING DAMAGE VERSUS FAULT COUNT, GPT-2 H .5. M L P . C _ P R O J , BF16 OPERANDS , FP32 DELIVERY, OUTPUT SIDE FAULTS ON BIT 26, CALIBRATED ENGINE WITH A DECLARED BUDGET OF s = 16. DAMAGE IS max |LOSS − LOSS CLEAN | AFTER INJECTION .
Milliseconds (log)
TABLE V R ECOVERY UNDER NVB IT ACCUMULATOR INJECTION ON A REAL L LAMA -2-7B D O W N _ P R O J LAYER (4096×11008×4096, ACTIVATIONS FROM A WIKITEXT-2 FORWARD PASS ) AND ON SYNTHETIC OPERANDS OF THE SAME SHAPE , CALIBRATED ENGINE AT ρmin = 0.02, TWO SITES PER TRIAL , 40 TRIALS PER CELL , FPR 0 IN EVERY CELL . P RODUCTS DELIVERED CORRECT AFTER REPAIR AND THE SECOND PROBE : 157 OF 160 ON THE REAL OPERANDS AND 78 OF 80 ON SYNTHETIC ONES .
101
BF16 GEMM
probe
100
10−1 qkv
ffn up
attn 70B
square
Fig. 3. Deployable probe cost against the guarded GEMM, BF16 operands, FP32 delivery, with activations rotating so the sketch of the activation side cannot be reused. Shapes: qkv 8192×4096×4096, ffn up 8192×4096×14336, attn 70B 16384×8192×8192, square 163843 . FP16: probe 0.780 to 6.675 ms against GEMMs of 0.360 to 12.736 ms.
E. RQ4: Cost
Clean path. The S only probe (§III.C) is the only stage paid on every call. Figure 3 measures it with activations rotating, Faults Repaired Unguarded Guarded which is the only setting that corresponds to a forward pass. −6 16 16/16 1.93 × 10 5.45 × 10−7 The probe takes 0.777 to 3.062 ms against GEMMs of 0.346 to 64 64/64 1.24 × 10−5 1.09 × 10−6 12.474 ms, that is, 25 to 225% of the GEMM and 8.1 to 21.0× 256 256/256 1.45 × 10−4 8.94 × 10−7 the bandwidth floor (π/β)(1/n3 + 2/n2 ) for reading A once 1024 379/1024 7.62 × 10−2 7.44 × 10−5 and C once. The term 1/n3 +2/n2 governs the ratio, so shallow 4096 381/4096 3.50 × 10−1 6.69 × 10−4 projections are expensive and a decision per layer follows directly. Segmented scatters and replay from a CUDA graph with a deferred verdict account for most of the reduction from guard removes 99.4% and 99.9% of the perplexity damage. the unoptimized path. The remaining gap to the floor is launch GPT-2 large recovers 2303 of 2304 and 2304 of 2304, and and synchronization overhead rather than bandwidth. With FP16 its guarded perplexity matches the clean one to four digits in operands, recalibration moved the probe at 16384×8192×8192 both formats. No run produced a false positive. With every one from the law’s m = 1142 to m = 656, where it costs 1.259 ms. of the 256 Llama calls and 288 GPT-2 large calls corrupted The FP16 square shape ran at the law’s m = 2397. Dirty path. Localization time follows the search radius, as and localized, a guarded pass took 48.2 s (BF16) and 41.5 s §III.G describes. At the Llama-2-7B down_proj shape, the (FP16) for Llama-2-7B against 3.2 to 3.5 s clean, and 75.9 s calibrated engine localizes in 8.5 to 9.0 ms on cuBLAS products and 25.8 s for GPT-2 large against 1.2 to 1.8 s. (mloc = 448 to 464, r = 4) and in 13.6 ms on products of the Table VI shows a dose threshold in training. Sixteen faults tiled kernel (mloc = 1024, r = 15), against a GEMM of 0.46 leave the loss essentially unchanged, and the damage becomes to 0.47 ms. The floor τd accounts for the second case. At almost large by 1024. Repair is complete up to 256 faults, which the same radius, localizing those products took 1162 to 1176 ms is 16× the declared budget. Beyond that, localization keeps before it was added, because clean buckets that cleared τc at most K = 128 candidates in each of its three rounds and walked the whole search. At 16384×8192×8192, localization repairs 379 and 381 faults. Even so, the guard lowers the takes 42 to 43 ms (mloc = 2896, r = 8 to 9) against a GEMM damage about 1000× at 1024 faults and 520× at 4096. Both of 2.86 to 2.96 ms, down from 1179 to 1191 ms at r = 16. experiments use output side injection and, by §VI.C, bound A deployment that wants only a correct product should still recovery from above. recompute on detection. Localization earns its place through Table V repeats the measurement under NVBit on the telemetry (§VI.G), not through repair economics. operands of a real Llama-2-7B down_proj layer, whose inner dimension of 11008 raises sketch noise. At the default F. RQ5: Comparison with Baselines ρmin = 0.02, with no manual setting, the calibrated engine Table VII compares FP-Sketch with checksum ABFT, recovers 0.975 to 1.000 on the real operands and 0.988 on weighted checksum ABFT, a Freivalds probe, and recompute synthetic operands of the same shape. For these products of the and compare on the same half precision operands and the same tiled kernel, the measured noise sets mloc = 1024. The probe corrupted products. All five detect every corrupted product, does not slow down as m adapts. Measured with the GPU to and all except Freivalds localize a single fault. Two faults per itself, it costs 0.62 to 0.69 ms on cuBLAS products and 0.93 product separate them. Checksum ABFT localizes 2 of 80, to 0.94 ms at the recalibrated m = 726 on products of the tiled weighted checksum ABFT localizes none and reports 2 false kernel, beside a GEMM of 0.46 to 0.47 ms, because reading positives, recompute and compare localizes 78, and FP-Sketch A and C dominates the probe at this shape. localizes all 80. On a corrupted call at qkv, FP-Sketch took
TABLE VII BASELINES ON IDENTICAL BF16 OPERANDS AND IDENTICAL FAULTS , WITH THE SAME ANALYTIC THRESHOLD FORM . D ET.: FRACTION OF TRIALS DETECTED . s=1: RECOVERY WITH ONE OUTPUT SIDE FAULT AT QKV, 20 TRIALS . s=2: RECOVERY WITH TWO NVB IT ACCUMULATOR FAULTS AT 4096×2048×4096, 40 TRIALS . FP: FALSE POSITIVES OVER BOTH WORKLOADS . MS : MEDIAN COST PER CALL AT QKV TO DETECT AND LOCALIZE , MEASURED WHILE NVB IT CAMPAIGNS SHARED THE GPU. Method
Det.
s=1
s=2
FP-Sketch Checksum ABFT Weighted checksum Freivalds Recompute
1.000 1.000 1.000 1.000 1.000
1.000 1.000 1.000 n/a 1.000
1.000 0 9.884 0.025 0 0.462 0.000 2 0.697 n/a n/a 0.253 0.975 0 5.824
FP
ms
9.9 ms to detect and localize, against 5.8 ms for a recompute measured under the same shared conditions. Row and column checksums flag two rows and two columns for two errors, which leaves four candidate positions, and a weighted syndrome averages two indices. Both therefore localize only single errors by construction. FT-GEMM is not measured. It is an FP32×FP32 kernel of its own rather than a wrapper around the vendor GEMM, so it cannot run this workload, and its overhead would be relative to a GEMM roughly an order of magnitude slower than the tensor core path. Table I compares it qualitatively. G. RQ6: What Localization Adds A recompute yields a correct product and a checksum yields one bit. Neither tells an operator which device to pull. Each localized fault carries its coordinate, its magnitude, and its direction, and aggregated records support three diagnoses that need progressively more data. The first is a repeated coordinate across independent products. Under a transient model, a coincidence among F faults on n1 n3 entries has probability ≈ 1−e−F (F −1)/2n1 n3 , about 1.5×10−4 for F = 18 on an output of 220 entries. The second is concentration within a tile of (i mod t1 , j mod t3 ), which needs roughly 640 faults for a χ2 test over a 16 × 8 tile. The third is a device whose rate is an outlier against the fleet. The flipped bit position is not recoverable, and the reason is structural. Inferring it from Cij and γ identifies the bit in 400/400 faults written into a finished product and in 0/80 under NVBit. An upset in the middle of accumulation is followed by the remaining partial products, so the delivered entry is not the clean value with one bit inverted. We therefore report coordinate, magnitude and direction only. Because every candidate is validated, a degraded localizer yields incomplete telemetry rather than wrong telemetry, which is another reason to verify at accumulation precision. VII. D ISCUSSION AND L IMITATIONS Trust. The guarantee against false positives assumes a recomputation that the fault cannot reach (§II.B), so a permanent fault on the verifying device is not covered. Scaling. Equation (4) makes m grow as n2.57 for square products. At 655363 the law asks for 25.8 GB of sketches, and the population
cap stops m at 16384, which bounds the method before it bounds accuracy. Transformer projections are rectangular and far below that size (Table II). Portability. The constants in (3) are specific to the BLAS, the device and the operand format, and FP-Sketch measures the noise at run time rather than relying on them (§III.G). Economics. The probe is the steady state cost and is expensive on shallow projections. Localization costs 15 to 30 GEMMs at the shapes measured, so repair still costs more than a recompute and localization is justified by the diagnostic record. Coverage. Operand corruption, faults outside matrix products, and fused kernels that never materialize their products are out of reach of any verifier that runs after the fact. Evaluation. Real model perplexity uses output side injection. Instruction level injection on the operands of a real Llama-2-7B layer recovers 0.988 of single bit faults with the calibrated engine, the same as on synthetic operands of that shape. VIII. C ONCLUSION FP-Sketch verifies an unmodified mixed precision GEMM after the fact, at the precision of its accumulator, and localizes several corrupted entries with no false positives by construction. Making hashed moment localization work in floating point required a measured noise law whose constant depends on the library, and deploying it required measuring that noise at run time. Evaluating it honestly required instruction level injection, because output side injection cannot produce the small faults that a live accumulator yields and so overstates what a localizer recovers. Localization is not the cheapest way to obtain a correct product, but of the three responses it is the only one that says where the corruption was. R EFERENCES [1] N. George et al., “Silent data corruption in AI,” Open Compute Project (OCP), White Paper, 2025. [Online]. Available: https: //www.opencompute.org/documents/sdc-in-ai-ocp-whitepaper-final-pdf [2] N. George, S. Gurumurthi, V. Sridharan, H. D. Dixit, E. Goksu, B. Parthasarathy, A. Huffman, T. Macieira, A. Sinha, D. Liberty, L. Minwell, and R. S. Chappell, “Silent data corruption in artificial intelligence: A growing challenge for large-scale machine learning,” IEEE Micro, vol. 46, no. 1, pp. 66–72, Jan.–Feb. 2026. [3] J. J. Ma, H. Pei, L. Lausen, and G. Karypis, “Understanding silent data corruption in LLM training,” arXiv preprint arXiv:2502.12340, 2025. [4] A. Altenbernd, P. Wiesner, and O. Kao, “Exploring silent data corruption as a reliability challenge in LLM training,” in Proceedings of the IEEE/ACM International Symposium on Cluster, Cloud and Internet Computing (CCGrid), 2026. [5] A. Tyagi, S. Hukerikar, N. Saxena, Y. Huang, P. Shirvani, C.-H. Tung, and Y. Zhu, “LLM-PRISM: Characterizing silent data corruption from permanent GPU faults in LLM training,” arXiv preprint arXiv:2604.10390, 2026. [6] K.-H. Huang and J. A. Abraham, “Algorithm-based fault tolerance for matrix operations,” IEEE Transactions on Computers, vol. C-33, no. 6, pp. 518–528, June 1984. [7] J.-Y. Jou and J. A. Abraham, “Fault-tolerant matrix operations on multiple processor systems using weighted checksums,” in Real-Time Signal Processing VII, Proceedings of SPIE, vol. 495, 1984, pp. 94–101. [8] J. A. Gunnels, D. S. Katz, E. S. Quintana-Ortí, and R. A. van de Geijn, “Fault-tolerant high-performance matrix multiplication: Theory and practice,” in Proceedings of the International Conference on Dependable Systems and Networks (DSN), 2001, pp. 47–56. [9] C. Ding, C. Karlsson, H. Liu, T. Davies, and Z. Chen, “Matrix multiplication on GPUs with on-line fault tolerance,” in Proceedings of the IEEE International Symposium on Parallel and Distributed Processing with Applications (ISPA), 2011.
[10] S. Wu, Y. Zhai, J. Liu, J. Huang, Z. Jian, B. M. Wong, and Z. Chen, “Anatomy of high-performance GEMM with online fault tolerance on GPUs,” in Proceedings of the 37th ACM International Conference on Supercomputing (ICS), 2023. [11] H. Shi, Z. Jiang, Z. Huang, B. Bai, G. Zhang, and H. Hou, “Grid-like error-correcting codes for matrix multiplication with better correcting capability,” arXiv preprint arXiv:2508.04355, 2025. [12] Y. Gao, Q. Hua, and Z. Chen, “V-ABFT: Variance-based adaptive threshold for fault-tolerant matrix multiplication in mixed-precision deep learning,” arXiv preprint arXiv:2602.08043, 2026. [13] C. Braun, S. Halder, and H.-J. Wunderlich, “A-ABFT: Autonomous algorithm-based fault tolerance for matrix multiplications on graphics processing units,” in Proceedings of the 44th Annual IEEE/IFIP International Conference on Dependable Systems and Networks (DSN), 2014. [14] X. Xue, C. Liu, H. Huang, B. Liu, Y. Wang, B. Yang, T. Luo, L. Zhang, H. Li, and X. Li, “ApproxABFT: Approximate algorithm-based fault tolerance for vision transformers,” arXiv preprint arXiv:2302.10469, 2023. [15] Y. Liang, X. Li, J. Ren, A. Li, B. Fang, and J. Chen, “ATTNChecker: Highly-optimized fault tolerant attention for large language model training,” in Proceedings of the 30th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming (PPoPP), 2025. [16] R. Freivalds, “Probabilistic machines can use less running time,” in Proceedings of the IFIP Congress, 1977, pp. 839–842. [17] H. Bennett, K. Gajulapalli, A. Golovnev, and E. Warton, “Matrix multiplication verification using coding theory,” in Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (APPROX/RANDOM 2024), ser. Leibniz International Proceedings in Informatics (LIPIcs), vol. 317. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2024, pp. 42:1–42:20. [18] L. Gasieniec, ˛ C. Levcopoulos, A. Lingas, R. Pagh, and T. Tokuyama, “Efficiently correcting matrix products,” Algorithmica, vol. 79, no. 2, pp. 428–443, 2017. [19] R. Pagh, “Compressed matrix multiplication,” ACM Transactions on Computation Theory, vol. 5, no. 3, pp. 9:1–9:17, 2013. [20] D. S. Roche, “Error correction in fast matrix multiplication and inverse,” in Proceedings of the ACM International Symposium on Symbolic and Algebraic Computation (ISSAC). ACM, 2018. [21] Y.-L. Wu and H.-L. Wang, “Correcting matrix products over the ring of integers,” arXiv preprint arXiv:2307.12513, 2024. [22] G. Cormode and S. Muthukrishnan, “An improved data stream summary: The count-min sketch and its applications,” Journal of Algorithms, vol. 55, no. 1, pp. 58–75, 2005. [23] D. Achlioptas, “Database-friendly random projections: Johnson–Lindenstrauss with binary coins,” Journal of Computer and System Sciences, vol. 66, no. 4, pp. 671–687, 2003. [24] M. T. Goodrich and M. Mitzenmacher, “Invertible bloom lookup tables,” arXiv preprint arXiv:1101.2245, 2011. [25] C.-H. Tung, Y. Huang, N. Saxena, P. Shirvani, S. Hukerikar, T. Jain, A. Tyagi, and S. Gongalore, “The anatomy of silent data corruption: GPU error pattern study and modeling guidance,” arXiv preprint arXiv:2605.04213, 2026. [26] B. Fang, X. Li, H. Dam, C. Tan, S. K. S. Hari, T. Tsai, I. Laguna, D. Tao, G. Gopalakrishnan, P. Nair, K. Barker, and A. Li, “MPGemmFI: A fault injection technique for mixed precision GEMM in ML applications,” in Proceedings of the IEEE International Conference on Cluster Computing (CLUSTER), 2024. [27] D. Chai, Z. Liu, S. Wang, S. Pei, C. Liu, H. Li, and S. Wang, “Analysis of LLM vulnerability to GPU soft errors: An instruction-level fault injection study,” arXiv preprint arXiv:2601.19912, 2026. [28] H. Cho, S. Mirkhani, C.-Y. Cher, J. A. Abraham, and S. Mitra, “Quantitative evaluation of soft error injection techniques for robust system design,” in Proceedings of the 50th Annual Design Automation Conference (DAC), 2013. [29] W. Ahrens, H. D. Nguyen, and J. Demmel, “Algorithms for efficient reproducible floating point summation,” ACM Transactions on Mathematical Software, vol. 46, no. 3, 2020. [30] T. S. Yashwanth, “On the structure of floating-point noise in batch-invariant GPU matrix multiplication,” arXiv preprint arXiv:2511.00025, 2025. [31] N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed. Society for Industrial and Applied Mathematics (SIAM), 2002.
[32] O. Villa, M. Stephenson, D. Nellans, and S. W. Keckler, “NVBit: A dynamic binary instrumentation framework for NVIDIA GPUs,” in Proceedings of the 52nd Annual IEEE/ACM International Symposium on Microarchitecture (MICRO-52). ACM, 2019, pp. 372–383. [33] T. Tsai, S. K. S. Hari, M. Sullivan, O. Villa, and S. W. Keckler, “NVBitFI: Dynamic fault injection for GPUs,” in Proceedings of the 51st Annual IEEE/IFIP International Conference on Dependable Systems and Networks (DSN), 2021, pp. 284–291.