PaperlayerView exact revision on arXiv

Source-linked AI summary

BF16 Component-Product Emulation of FP32 and FP64 GEMM on Intel AMX

Bing Cui, Yu Liu

arXiv:2609.04663v1cs.MScs.DC

TL;DR

Scientific computing still needs FP32 and FP64 GEMM, whereas modern CPU matrix engines emphasize low-precision AI arithmetic. This paper uses Intel AMX BF16 component products to emulate higher-precision GEMM, reporting higher FP32 throughput than oneMKL SGEMM and an FP64 accuracy–performance tradeoff across 6, 10, 15, and 21 products.

  • Problem

    Scientific applications require FP32 and FP64 GEMM for numerical accuracy, but AI-oriented hardware primarily provides higher-throughput low-precision matrix operations.

  • Method

    The method decomposes operands into BF16 components, evaluates selected AMX products, and reconstructs results in wider precision for FP32 and FP64 GEMM emulation.

  • Results

    AMX-FP32 exceeds oneMKL SGEMM throughput, while AMX-FP64 exposes accuracy–performance tradeoffs through 6, 10, 15, and 21 retained products.

  • Takeaways & Limitations

    The approach provides a CPU-oriented route for using AI-style low-precision hardware to accelerate higher-precision GEMM within the tested configurations and supported domains.

  • Takeaways & Limitations

    The FP64 path is restricted to inputs whose extracted components remain representable in BF16 and does not implement scaling or fallback processing for excluded cases.

Abstract

from arXiv · show

Modern CPUs increasingly integrate high-throughput matrix engines optimized for low-precision AI workloads, while many scientific computing applications still rely on FP32 and FP64 GEMM to meet their numerical accuracy requirements. This mismatch motivates an algorithmic bridge that uses low-precision matrix products to emulate higher-precision GEMM. This paper presents a CPU-oriented method based on Intel Advanced Matrix Extensions (AMX) and BF16 matrix products. For FP32, each operand is decomposed into three BF16 components and six selected component products are evaluated, targeting FP32-level accuracy relative to oneMKL SGEMM without claiming elementwise or bitwise identity. For FP64 inputs within the supported BF16 exponent range, the method uses a simplified fixed six-slice Ozaki decomposition. Each retained BF16 product is first produced in FP32, then widened and accumulated in FP64. Four product-count settings retain 6, 10, 15, or 21 component products, exposing the accuracy--performance tradeoff relative to oneMKL DGEMM. The implementation combines precomputed packed component buffers, VNNI-packed $B$ panels, and an FP32 tile-resident operand-reuse schedule. On the tested square matrices, AMX-FP32 exceeds oneMKL SGEMM throughput. For AMX-FP64, low-product-count variants can exceed DGEMM at sufficiently large orders, while retaining more products improves accuracy at additional cost.

1 Introduction

The paper asks whether AI-oriented low-precision matrix engines can accelerate FP32 and FP64 GEMM needed by scientific computing. It proposes decomposing high-precision operands into BF16 components, selecting component products, and reconstructing results in wider precision on Intel AMX.

  • Scientific applications continue to rely on FP32 and FP64 arithmetic for stability, convergence, and reproducibility, while low-precision engines offer much higher throughput.
  • GEMM is a natural target because it is central to machine learning and appears across dense and sparse scientific-computing workloads.
  • Intel AMX is used to evaluate whether AI-style low-precision CPU hardware can accelerate high-precision GEMM through component decomposition and wider reconstruction.
  • High-precision emulation requires numerical product selection and wider accumulation alongside conversion, packing, blocking, scheduling, and reconstruction.
  • The FP32 path uses three BF16 components and six selected products, while the FP64 path uses six BF16 slices with a tunable retained-product count and compares against oneMKL baselines.

2 Related Work

Prior work connects low-precision matrix hardware with mixed-precision algorithms, floating-point expansions, and component-product GEMM emulation. This paper builds on those directions while focusing on CPU AMX execution and controllable product counts.

  • Mixed-precision methods reserve wider arithmetic for sensitive operations and have reduced cost or data motion in several scientific applications.
  • Floating-point expansions represent values with leading and residual components, while accurate summation addresses cancellation and rounding during reconstruction.
  • Expansion-based matrix multiplication decomposes operands into low-precision components and reconstructs results from selected pairwise products.
  • Henry, Tang, and Heinecke showed that three BF16 components and six selected products can recover FP32-oriented matrix-multiplication accuracy.
  • Ozaki-style methods extend component-product emulation to higher precision, with scaling policies and retained product counts governing accuracy and cost.
  • Subsequent work explicitly controls product counts, analyzes accuracy requirements, and develops reduced-precision routes to DGEMM accuracy.
  • The paper positions its AMX study within this broader connection between low-precision hardware, floating-point expansions, and accurate GEMM emulation.

3 Background and Problem Setting

The paper defines FP32, FP64, and BF16-range input domains, then models GEMM emulation as decomposition, packing, multiple AMX products, reconstruction, and parallel overhead. AMX supplies BF16 tile products whose results are widened or reconstructed for the target format.

  • 3.1 Floating-Point Formats: BF16 has 8 significand bits and the FP32 exponent range, compared with 24 bits for FP32 and 53 bits for FP64.
  • 3.1 Floating-Point Formats: The FP64 path is restricted to finite nonzero values within the normal BF16 exponent range, while retaining FP64 significand precision.
  • 3.1 Floating-Point Formats: The unscaled six-slice FP64 prototype further excludes inputs whose extracted residual components become subnormal or zero in BF16.
  • 3.2 Problem Definition and Performance Model: The modeled emulation time combines decomposition, packing, Np component GEMMs, reconstruction, and parallel overhead.
  • 3.2 Problem Definition and Performance Model: AMX-FP32 targets FP32-level accuracy relative to oneMKL SGEMM rather than bitwise reproduction, while AMX-FP64 compares accuracy and performance against oneMKL DGEMM as products increase from 6 to 21.
  • 3.3 CPU Matrix Engines and Intel AMX: BF16 is selected over FP16 because its exponent range supports direct FP32 decomposition and the supported BF16-range FP64 path without per-panel scaling.
  • 3.3 CPU Matrix Engines and Intel AMX: AMX uses explicit tile loads, BF16 dot products into FP32 destination tiles, and tile stores; the implementation adds VNNI packing, blocking, operand reuse, and reconstruction.

4 BF16-Based High-Precision GEMM

The method emulates FP32 and FP64 GEMM with selected BF16 component products on Intel AMX, using decomposition and wider-format reconstruction. It trades component-product count and implementation cost against accuracy, with explicit limits for difficult FP32 cases and unsupported FP64 residuals.

  • AMX-FP32: The FP32 path retains six triangular component products, following the Henry-style schedule, and accumulates the AMX results in FP32.The schedule keeps the first three diagonals and omits later products and residual-containing terms.
  • Accuracy boundary: The six-product FP32 scheme targets FP32-level accuracy but does not guarantee elementwise or bitwise identity with oneMKL SGEMM.Cancellation, underflow, boundary cases, and FP32 accumulation order can produce larger errors for difficult matrices.
  • Decomposition: FP32 operands use three BF16 residual components, while FP64 operands use six fixed BF16 slices in a simplified Ozaki-style decomposition.The FP32 split uses direct residual rounding; the FP64 path is fixed-depth and unscaled.
  • FP64 scope: Direct BF16-slice storage removes separate scale metadata and rescaling, but requires extracted FP64 slices to remain normal BF16 values.Subnormal or zero lower residual slices are outside the current unscaled prototype and require a scaled or fallback extension.
  • AMX-FP64: FP64 uses four triangular schedules with 6, 10, 15, or 21 retained products, forming a controlled progression from lower cost and accuracy to higher cost and accuracy.The selected pairs are ordered by increasing slice-index sum; more products increase reconstruction work.
  • AMX-FP64: Each FP64 component GEMM produces an FP32 tile result that is materialized, widened to FP64, and added before the next product.This avoids FP32 accumulation across distinct component products but adds tile stores, conversion, and FP64 additions.

5 Numerical Accuracy Considerations

The accuracy analysis attributes deviations to omitted component products, residual terms, accumulation precision, and reconstruction choices. It evaluates accuracy experimentally against oneMKL rather than claiming bitwise or complete FP64 agreement.

  • Accuracy scope: The method targets FP32-level accuracy relative to oneMKL SGEMM, not bitwise or elementwise identity.The claim is assessed experimentally, with omitted products and FP32 accumulation contributing to deviations.
  • FP32 omissions: Six selected FP32 products omit lower-order cross-products and terms involving residual matrices R_A and R_B.These omissions explain why the method targets accuracy relative to SGEMM rather than exact reproduction.
  • FP64 error sources: FP64 error arises from terminal residuals, ordinary FP64 residual-update rounding, and truncation of unretained component products.The fixed six-slice split leaves residuals that are not represented by stored BF16 components.
  • FP64 product schedules: Adding triangular diagonals reduces the omitted FP64 product tail but does not guarantee full FP64 accuracy for all matrices.Cancellation, input scaling, and residuals from the fixed six-slice split can still dominate the final error.
  • Reconstruction: AMX-FP32 combines the six products in FP32, retaining them in AMX accumulator tiles but fixing a particular summation order.This choice preserves the performance motivation while making accumulation order part of the numerical behavior.

6 AMX Implementation

The implementation uses blocked AMX microkernels with precomputed BF16 components, VNNI-packed B panels, and distinct FP32 and FP64 reconstruction paths. FP32 exploits tile-resident operand reuse, while FP64 stores and widens each component product before accumulation.

  • Blocked kernel organization: The prototype partitions GEMM into cache-resident BM × BN blocks, BK reduction panels, and TM × TN AMX microtiles.The selected microkernel uses TM = TN = 32 within an eight-tile register budget.
  • Tile mapping: A 32 × 32 microtile uses two A-input tiles, two VNNI-packed B tiles, and four FP32 accumulator quadrants.The B tiles represent distinct output-column groups, while the four accumulators cover 16 × 16 output quadrants.
  • Decomposition and packing: Both precision paths precompute decomposition and packing, then reuse the same AMX-facing component-buffer dataflow.FP32 uses a three-component residual split; FP64 uses an Ozaki-style split before common packing logic.
  • Product scheduling: The FP32 operand-reuse schedule reduces six-product traffic from 12 component-panel loads to 7, saving 10 explicit tile-load instructions.Adjacent products share one operand, while the FP32 accumulator tiles remain resident across the sequence.
  • Reconstruction and parallelization: FP32 keeps accumulator tiles live across six products, whereas FP64 stores each AMX-BF16 result, converts it to FP64, and accumulates outside AMX.The two paths share blocked parallelization but diverge after each component product.

7 Experimental Methodology

The study benchmarks AMX emulation against oneMKL SGEMM and DGEMM using controlled, repeated measurements on square GEMM problems. Accuracy uses complementary normwise and componentwise metrics, while performance uses wall-clock time and normalized throughput.

  • Hardware and configuration: The experiments run on a dual-socket Intel Xeon Platinum 8462Y+ server with Intel oneMKL baselines and fixed frequency policy.The platform uses Rocky Linux 9.3, oneAPI 2026.1, and compact OpenMP thread placement.
  • Measurement protocol: Each configuration uses 64 OpenMP threads, one warm-up invocation, and 10 timed repetitions reported as an arithmetic mean.The AMX warm-up includes operand decomposition.
  • Baselines and variants: The principal baselines are oneMKL SGEMM for AMX-FP32 and oneMKL DGEMM for AMX-FP64, with variants using common layouts, thread counts, and timing.Table 3 organizes the evaluated emulation variants under these corresponding baselines.
  • Input matrices: Validation considers random, scaled or log-uniform, and cancellation-sensitive matrices, but detailed supplementary results are not reported across tested orders.The reported Section 8.2 accuracy results use uniform random inputs.
  • Experimental scope: All reported experiments use square GEMM with m = n = k = N, leaving rectangular GEMM outside the study’s scope.The primary references are CS from oneMKL SGEMM and CD from oneMKL DGEMM.
  • Metrics: Performance uses wall-clock time, effective throughput 2N^3/T, and speedup over the corresponding oneMKL baseline.The normalization counts requested high-precision GEMM work rather than the larger number of BF16 component operations.

8 Results

The evaluation first measures raw AMX-BF16 headroom, then assesses accuracy and the performance implications of component expansion. Results show substantial but bounded throughput potential and an accuracy–cost tradeoff for retained FP64 products.

  • Design choices: Six-product FP32 emulation is motivated by the need to limit component-expansion overhead while retaining a plausible performance opportunity.The FP64 path instead evaluates 6, 10, 15, and 21 products as distinct accuracy/performance points.
  • Performance headroom: Roughly 4× SGEMM throughput and about 8× DGEMM throughput represent the best measured large-matrix AMX-BF16 headroom.This raw advantage is measured before decomposition, packing, and reconstruction costs are included.
  • Accuracy: FP32 Frobenius relative error remains on the order of 10^-7, with scaled componentwise error at the same or smaller order.The reported comparison uses uniform random inputs and oneMKL references.

N AMX-FP32 AMX-FP64-6 AMX-FP64-10 AMX-FP64-15 AMX-FP64-21

Across the tested square matrices, AMX-FP32 remains faster than oneMKL SGEMM, while AMX-FP64 performance depends strongly on retained product count and matrix size. Runtime composition shows preprocessing becomes less dominant as order grows.

  • AMX-FP32 performance: AMX-FP32 is faster than oneMKL SGEMM for every tested square size, with speedup recovering to 1.16×–1.32× for N ≥8192.Intermediate sizes show a minimum speedup because decomposition, packing, six products, and reconstruction remain costly.
  • AMX-FP32 scaling: At N = 1024, both AMX-FP32 and SGEMM peak at 32 threads before throughput drops at 64 threads.For N = 4096 and N = 16384, throughput continues increasing, while AMX-FP32’s relative speedup becomes smaller at high thread counts.
  • AMX-FP64 performance: AMX-FP64-6 reaches about 1.7× for the two largest sizes after crossing over near N = 4096, whereas higher-product variants retain smaller margins.AMX-FP64-10 exceeds DGEMM only for the largest two matrices; the 15- and 21-product variants do not outperform DGEMM in the reported passage.
  • Runtime breakdown: Figure 10 normalizes each variant’s runtime into preprocessing and compute/reconstruction stages across matrix sizes.The same full-core configuration is used as in the corresponding performance experiments.
  • Runtime breakdown: AMX-FP32 preprocessing falls from one third of total runtime at N = 256 to less than 4% at N = 32768.AMX-FP64 preprocessing reaches 56%–75% at N = 256 because its six-slice decomposition generates and packs more components.

9 Discussion and Future Work

The study’s conclusions are bounded by square real GEMM evaluation, AMX tile-register resources, portability engineering requirements, and a narrower FP64 input domain than native DGEMM.

  • Scope: The evaluation covers real square C = AB problems and does not establish performance for rectangular, skinny, or more general GEMM interfaces.Those shapes are identified as targets for future blocking, packing, scheduling, and application-specific optimization.
  • Hardware boundary: The eight-tile AMX register file limits resident component panels, the 32 × 32 microkernel, and removable tile load/store traffic.More tile registers or larger tile capacity could increase the attainable speedup.
  • Portability: Adapting the component-expansion framework beyond x86 AMX requires platform-specific packing, blocking, and microkernel design.The passage identifies Arm SME/SME2 implementations as potential targets while retaining the decomposition and wider-accumulation principles.
  • Input domain: The unscaled AMX-FP64 method does not cover the full FP64 exponent range and requires extracted components to remain representable as normal BF16 values.Residuals near the lower BF16 boundary that become subnormal or zero are outside the supported domain because scaling and fallback processing are omitted.

10 Conclusion

The paper uses AMX-BF16 component products to emulate higher-precision GEMM, with product count controlling the FP64 accuracy–performance tradeoff. In square full-core experiments, FP32 emulation outperforms SGEMM, while FP64 gains narrow as more products are retained.

  • Method: AMX-FP32 combines three BF16 components, six selected products, and resident FP32 accumulator tiles to target FP32-level accuracy relative to oneMKL SGEMM.The target is not elementwise or bitwise identity.
  • Method: AMX-FP64 uses a fixed six-slice decomposition and retains 6, 10, 15, or 21 products, exposing an accuracy–performance tradeoff.Products are materialized in FP32, converted, and accumulated in FP64.
  • Results: 1.14×–2.56× AMX-FP32 speedup over oneMKL SGEMM is reported across tested square sizes while meeting FP32-level accuracy criteria.AMX-FP64-6 reaches up to 1.72× over DGEMM; AMX-FP64-10 gains only slightly at the two largest sizes, while 15 and 21 products do not outperform DGEMM.
  • Conclusion: Selective use is supported when input range, matrix shape, and application-level accuracy requirements match the method’s operating regime.Increasing product count consumes available AMX-BF16 headroom through additional stores, conversions, and FP64 accumulation.
Loading 2609.04663v1…