High-performance implementations of transformer inference (attention kernel) and options pricing on Apple Silicon ARM64, demonstrating when SIMD helps and when it doesn't.
Attention Kernel: 31.88 GFLOPS/s — 11.1× speedup over naive baseline
Black-Scholes Pricer: 103.8M options/sec — 6.3× faster than scipy
Two kernels, two different bottlenecks:
| Kernel | Bottleneck | SIMD Benefit | Result |
|---|---|---|---|
| Scaled dot-product attention | Memory-bound | 11.1× speedup | V-transpose + NEON |
| Black-Scholes option pricing | Transcendental-bound | 1.0× speedup | C++ beats scipy by 6.3× anyway |
The contrast between these two kernels is the core lesson: roofline analysis predicts vectorization benefit before you write a line of SIMD code.
Attention(Q, K, V) = softmax(QK^T / √d_k) V
Complexity: O(seq_len² × d_k) FLOPs — the bottleneck of transformer inference.
| seq_len | Naive (GFLOPS/s) | SIMD (GFLOPS/s) | Speedup | OMP (GFLOPS/s) | Total |
|---|---|---|---|---|---|
| 64 | 1.86 | 29.95 | 16.1× | 49.82 | 26.8× |
| 128 | 2.24 | 31.70 | 14.2× | 48.77 | 21.8× |
| 256 | 2.87 | 31.88 | 11.1× | 50.40 | 17.6× |
| 512 | 2.89 | 29.14 | 10.1× | 52.21 | 18.1× |
Thread scaling (seq_len=256):
| Threads | GFLOPS/s | Efficiency |
|---|---|---|
| 1 | 51.38 | 100% |
| 2 | 48.77 | 95% |
| 4 | 50.01 | 97% |
| 8 | 52.21 | 102% |
Minimal scaling across threads (1.01× from 1→8) confirms memory bandwidth is saturated, not compute.
1. NEON 8-wide dot product with dual accumulators (core speedup)
Process 8 floats per iteration using two independent NEON registers. Independent accumulators hide memory load latency. Fused multiply-add via vfmaq_f32().
float32x4_t acc1 = vdupq_n_f32(0.0f);
float32x4_t acc2 = vdupq_n_f32(0.0f);
for (size_t i = 0; i + 8 <= n; i += 8) {
acc1 = vfmaq_f32(acc1, vld1q_f32(&a[i]), vld1q_f32(&b[i]));
acc2 = vfmaq_f32(acc2, vld1q_f32(&a[i+4]), vld1q_f32(&b[i+4]));
}
float result = horizontal_sum_f32(vaddq_f32(acc1, acc2));2. V-transpose preprocessing (1.5–2× gain)
Column access to V requires stride=d_model, causing cache misses on every access. Transposing V to V^T converts strided access to sequential — trading O(seq_len × d_k) setup cost for O(seq_len²) speedup. This single optimization contributes ~50% of the total SIMD speedup.
3. Vectorized softmax (1.1× gain)
Fast exp via polynomial approximation (7 decimal places accuracy). Process 4 elements per NEON cycle with numerically stable implementation (subtract max before exp).
4. OpenMP outer loop parallelism (1.6× additional)
Each output row is independent. Thread creation overhead exceeds benefit at small sizes; add threshold check before enabling.
| Implementation | Achieved (GFLOPS/s) | Peak (GFLOPS/s) | Utilization | Bottleneck |
|---|---|---|---|---|
| Naive | 2.87 | 140 | 2.1% | Scalar overhead |
| SIMD | 31.88 | 140 | 22.8% | Memory-bound (35 GB/s realized) |
| OpenMP | 50.40 | 140 | 36.0% | Bandwidth saturated |
Hardware: Apple Silicon M-series — Peak: 140 GFLOPS/s, Bandwidth: 100 GB/s theoretical, 35 GB/s realized, Ridge point: 1.4 FLOP/byte
The kernel is memory-bound at all tested sizes. Adding threads does not increase total memory bandwidth — each additional thread contends for the same hardware bus. 36% peak utilization is excellent for a real-world memory-bound kernel.
Further optimization requires cache blocking (tile the computation to fit L3), quantization (fewer bytes per element), or a different algorithm with lower memory footprint.
To validate whether hand-coded SIMD is worth the effort, benchmarked against cblas_sgemm from Apple's Accelerate framework:
| seq_len | SIMD (GFLOPS/s) | Accelerate (GFLOPS/s) | Accelerate speedup |
|---|---|---|---|
| 64 | 17.0 | 97.2 | 5.7× |
| 128 | 24.7 | 79.0 | 3.2× |
| 256 | 32.2 | 109.3 | 3.4× |
| 512 | 29.1 | 115.7 | 4.0× |
Accelerate is 3–6× faster. It uses cache-blocking to tile computation into L2/L3, cycle-accurate instruction scheduling, and runtime dispatch across multiple optimization paths. Hand-coded NEON uses a single approach with no tiling.
Verdict: Use Accelerate in production. Hand-code SIMD only when you need custom kernel fusion, multi-ISA portability, or are measuring optimization potential (as this project does).
C = S·N(d1) - K·e^(-rT)·N(d2)
d1 = (ln(S/K) + (r + σ²/2)T) / (σ√T)
d2 = d1 - σ√T
| Num Options | Scalar (ms) | SIMD (ms) | Speedup | Options/sec |
|---|---|---|---|---|
| 1K | 0.010 | 0.010 | 1.00× | 102M |
| 10K | 0.136 | 0.215 | 0.64× | 47M |
| 100K | 1.723 | 1.442 | 1.19× | 69M |
| 1M | 9.632 | 9.374 | 1.03× | 107M |
Average SIMD speedup: 1.00× — SIMD provides no benefit.
scipy (vectorized Python): 0.061s → 16.4M options/sec
C++ scalar: 0.010s → 103.8M options/sec
Speedup: 6.3×
C++ beats scipy not through SIMD but through eliminated Python overhead and a faster polynomial approximation for N(x).
Runtime breakdown per option:
4 × normal CDF() 87% of runtime ← polynomial approx, scalar-friendly
4 × exp() 5% of runtime ← NEON available
2 × log() 5% of runtime ← no NEON instruction
2 × sqrt() 2% of runtime ← NEON available
arithmetic 1% of runtime ← fully vectorizable
89% of runtime is transcendental functions. SIMD can parallelize the remaining 11% at most, yielding a theoretical maximum of 1.13× — which is not worth the extraction/reconstruction overhead from AoS data layout.
Amdahl's Law applied: Speedup ≤ 1 / (0.89 + 0.11/8) = 1.12×
To get meaningful Black-Scholes speedup beyond C++ scalar:
- Lookup tables: Precompute N(d) for d ∈ [−5, 5] at resolution 0.001 → O(1) lookup, ~3× speedup
- Struct-of-arrays layout + vectorized math library: Intel MKL or ARM optimized math for batch transcendentals → 3–4×
- Accept scalar: 100M options/sec is sufficient for most non-HFT workloads
| Algorithm class | Vectorization benefit | Example |
|---|---|---|
| Compute-bound (arithmetic) | 3–8× | Matrix multiply, FFT, convolution |
| Memory-bound (cache-friendly) | 2–4× | Attention + V-transpose |
| Transcendental-bound | 0–1× | Black-Scholes, Monte Carlo |
| Bandwidth-saturated | 0.5–1× | Attention without V-transpose |
Roofline model analysis tells you which category you are in before benchmarking. Arithmetic intensity (FLOP/byte) determines whether compute or memory is your ceiling.
Current benchmark: seq_len=64–512, fits in L3 cache (~16 MB working set).
At seq_len=4096: Working set ~256 MB, exceeds L3 by 16×. Realized bandwidth drops to ~15 GB/s. Expected throughput ~8 GFLOPS/s.
At seq_len=16384: Working set ~4 GB, exceeds main memory. Single-machine approach requires streaming or distributed compute.
Distributed options:
- Data parallelism (independent batches): linear scaling, trivial communication. Suitable for portfolio risk across subsets, multi-request inference.
- Sequence parallelism (split seq_len across nodes): requires AllReduce after QK^T, communication-bound at high node counts.
- Head parallelism (split attention heads): no synchronization needed, limited to number of heads (8–64).
Data parallelism is the right choice whenever the workload fits the batch structure. Head parallelism is the right choice for single long-sequence inference.
# Prerequisites
brew install gcc
# Build
make all # all targets
make test # correctness tests (outputs nothing on pass)
make benchmark # attention benchmark table
make benchmark-csv # write results.csv
make roofline # roofline analysis to roofline.txt
make bench-bs # Black-Scholes benchmark + scipy comparison
make cleanExpected test output:
All tests passed
Expected attention benchmark (excerpt):
seq_len Version Threads Time (ms) GFLOPS/s Speedup
256 Naive 1 5.97 2.87 1.0x
256 SIMD 1 0.54 31.88 11.1x
256 OMP 8 0.33 52.21 18.2x
Expected Black-Scholes output:
1M options — Scalar: 9.63ms (103.8M/sec) | SIMD: 9.37ms (1.03x)
scipy comparison — C++: 0.010s | scipy: 0.061s | Speedup: 6.3x
cpp_simd_attention/
├── README.md
├── Makefile
├── include/
│ ├── matrix.h # matrix class with aligned allocation
│ ├── attention.h # function signatures
│ ├── simd_utils.h # NEON utility functions, fast exp approximation
│ └── black_scholes.h # BS pricer signatures
├── src/
│ ├── attention_naive.cpp # scalar baseline
│ ├── attention_simd.cpp # NEON vectorized
│ ├── attention_omp.cpp # NEON + OpenMP
│ ├── black_scholes_scalar.cpp # scalar BS (optimal for this bottleneck)
│ └── black_scholes_simd.cpp # attempted SIMD — reverts to scalar (see Part 2)
├── benchmarks/
│ ├── test.cpp # correctness tests
│ ├── bench_attention.cpp # timing harness
│ ├── bench_accelerate.cpp # vs Accelerate cblas_sgemm
│ └── bench_black_scholes.cpp # BS + scipy comparison
└── analysis/
├── roofline.py # generates roofline analysis
├── compare_scipy.py # scipy benchmark script
└── results.csv # generated by make benchmark-csv
| Version | Flags |
|---|---|
| Naive | -O2 -std=c++17 |
| SIMD | -O2 -march=armv8.2-a+simd -std=c++17 |
| OpenMP | -O2 -march=armv8.2-a+simd -fopenmp -std=c++17 |
-O2 not -O3 to prevent auto-vectorization from obscuring the naive vs SIMD comparison.
- ARM NEON Intrinsics Guide: https://developer.arm.com/architectures/instruction-sets/intrinsics/
- Roofline Model (Williams et al. 2009): https://crd.lbl.gov/assets/pubs_presos/parlab08-roofline.pdf
- Agner Fog CPU Optimization Manuals: https://www.agner.org/optimize/
- Flash Attention (Dao et al. 2022): https://arxiv.org/abs/2205.14135
Compiler: gcc-15 | Language: C++17 | Target: Apple Silicon ARM64