PaperlayerView exact revision on arXiv

Source-linked AI summary

BRiDCT: Fast Two-Dimensional DCTs Using SIMD: SIMD Organization, Register Blocking, and Numerical Verification

Antoine Moevus, Max Mignotte

arXiv:2609.28519v1eess.IVeess.SP

TL;DR

Repeated two-dimensional DCTs need fast implementations across sizes, but arithmetic counts alone do not predict CPU time. BRiDCT combines an established factorization with SIMD and storage strategies, and on the tested Apple M3 Max outperforms the evaluated general-purpose routes across their available cases.

  • Problem

    Fast float32 DCTs are needed across multiple array sizes, but fewer arithmetic operations alone do not guarantee shorter CPU time.

    Source: §1.1 · pp. 1–2

  • Method

    BRiDCT combines the Shao–Johnson factorization with SIMD lane organization, banded intermediate storage, padding, and register blocking.

    Source: §1.3 · p. 3

  • Results

    On the tested Apple M3 Max, BRiDCT is approximately 1.3–16.0× faster than tested general-purpose routes across available native cases from 8 × 8 to 256 × 256, and its compiled NumPy interface beats OpenCV in every tested contract.

    Source: §9 · p. 12

  • Takeaways & Limitations

    The measured gains come from organizing SIMD execution and intermediate storage rather than introducing a new arithmetic factorization.

    Source: §9 · p. 12

    Boundary diagnostics include overflow and underflow or rounding errors, and weak-component loss can pass the global norm criterion, so weak spectral detail may remain unreliable.

    Source: §6 · p. 10

Abstract

from arXiv · show

Repeated two-dimensional discrete cosine transforms (DCTs) require implementations that remain fast across several array sizes while preserving the mathematical transform. Fewer arithmetic operations do not guarantee a faster DCT: execution time also depends on how the computation uses the CPU's single-instruction, multiple-data (SIMD) units, registers and memory hierarchy. We implement established DCT algorithms in SIMD and compare execution strategies for forward, inverse, and normalized round-trip transforms on float32 arrays from $8\times8$ to $256\times256$. The resulting library, BRiDCT, combines the Shao--Johnson factorization with SIMD lane organization, banded intermediate storage, padding, and register blocking in C11/ARM NEON. On one Apple M3 Max, the final native library has lower execution times than the tested general-purpose routes---Apple's Accelerate framework (vDSP), FFTW and Ooura---across their available cases from $8\times8$ to $256\times256$. It also outperforms the tested libjxl and floating-point libjpeg-turbo routes. The compiled NumPy interface outperforms OpenCV in all 130 tested contracts per session, including squares up to $1024\times1024$ and two rectangles, with call and allocation costs included. Complementary comparisons evaluate SciPy, DUCC and cached/planned pyFFTW. The measured gains come from organizing SIMD execution and intermediate storage, rather than introducing a new arithmetic factorization. A reproducible 474-input verification corpus probes cancellation, neighboring floats, dynamic range, and extreme amplitudes: all 1 260 core checks of the final configuration pass, while boundary diagnostics reveal overflow and weak-component limitations. Larger native transforms show size-dependent limits.

1 Introduction

The paper targets float32 two-dimensional DCTs across six square sizes, where arithmetic counts alone do not predict CPU time. It combines established transforms with measured SIMD and memory organization strategies.

1.1 The computational problem

  • The target is fast forward, inverse, and normalized round-trip DCTs on row-major float32 arrays from 8 × 8 to 256 × 256.A fast 8 × 8 routine alone does not cover the stated six-size problem.Source: §1.1 · p. 1
  • SIMD use, register pressure, dependencies, and data access can offset arithmetic savings, so the preferred execution organization may vary with size and architecture.The approach compares factorizations together with SIMD and memory organization.Source: §1.1 · p. 2

1.2 Existing algorithms and implementations

  • Shao–Johnson lowers scalar arithmetic counts, while SPIRAL establishes prior work on joint algorithm and SIMD optimization; BRiDCT uses its own generators and policies.The attempted SPIRAL integration did not yield a validated float32 NEON DCT for this protocol, so SPIRAL is absent from the timing comparison.Source: §1.2 · p. 2

1.3 Approach and contribution

  • BRiDCT combines Shao–Johnson with four-row bands, padded storage, and register blocking; small kernels retain simpler organizations when faster.The contribution is an implementation and evaluation connecting complete source changes to measured gains, with numerical errors checked independently.Source: §1.3 · p. 3

2 Transform Contract and Algebraic Basis

The paper fixes the target as the orthonormal two-dimensional DCT and its adjoint inverse, with the normalized round trip returning the input. It also derives scaled-core identities to place normalization efficiently.

2.1 Forward, inverse, and round trip

  • The forward and inverse operators define the exact transform target, while the timed normalized round trip is C−1(C(X)) and has exact result X.Exactness excludes approximate matrices, not floating-point error; the inverse is the adjoint under the Frobenius inner product, not a memory transpose.Source: §2.1 · p. 3
  • Adapters absorb unnormalized DCT scale factors so every comparison delivers the same output normalization; impulse responses check those factors before timing.For FFTW’s REDFT10 and REDFT01 pair, the unnormalized one-dimensional factor is 2n.Source: §2.1 · p. 3

2.2 Scaled cores and adjoints

  • The scaled-core formulation folds forward scaling into final stores and round-trip scaling into one correction table; incorporating a solver’s diagonal operator is algebraically possible but not benchmarked.The adjoint of a generated linear graph is obtained by reversing dependencies, accumulating through additions, and propagating constant multiplications with the same constants.Source: §2.2 · p. 4

3 Implementation Organization

BRiDCT applies separable one-dimensional transforms with SIMD lanes arranged for the vertical and horizontal passes on ordinary row-major arrays. Padding and register blocking address execution costs that vary across sizes.

3.1 Two SIMD passes on a natural-layout array

  • Interleaving loads and stores pack and unpack rows, removing a full-array transpose but not all memory rearrangement; smaller policies omit padding when its cost outweighs its savings.Packing only permutes values, so the two passes preserve the target transform; the inverse uses the transposed core and corresponding input scales.Source: §3.1 · p. 4
  • In the banded route, four neighboring columns occupy vertical-pass lanes, then four intermediate rows are packed into a band and transformed horizontally.The band occupies scratch memory; this is not a register-only transpose.Source: §3.1 · p. 4, Figure 1 · p. 5

3.2 Blocking, padding, and specialization

  • Blocked kernels limit generated leaves to at most 16 points and use loops and coefficient tables for larger combinations, reducing live values that can spill from registers.This changes execution organization, not the underlying Shao–Johnson factorization.Source: §3.2 · p. 5
  • Padded variants use an n + 4 float stride to perturb cache-set access patterns, but timing does not reveal physical cache indexing or isolate cache conflicts.Padding adds scratch capacity and copies rather than reducing memory footprint.Source: §3.2 · p. 5

4 Experimental Protocol

The evaluation compares a fixed final BRiDCT release against documented native and Python routes, with separate protocols for each. Timing includes specified adapter or call costs while excluding planning and preparation.

4.1 Native comparisons and provenance

  • Native comparisons use one Apple M3 Max, float32 orthonormal outputs, one worker, and a rotating pool of four inputs; adapter input copies count, while preparation and planning are excluded.The primary comparisons evaluate the final distributed BRiDCT version 0.2.0-dev.7 without selecting policies from those measurements.Source: §4.1 · p. 6
  • Two fresh-process native sessions use 21 randomized paired blocks per arm and a 30 ms calibration target, retaining all observations; exclusive desktop use is not guaranteed.FFTW PATIENT planning is capped at one second per plan.Source: §4.1 · p. 6

4.2 Python calls and native execution

  • Python timing covers complete calls through the compiled NumPy extension, including dispatch, validation, output allocation, execution, normalization, copies, and output destruction.Imports, input generation, explicit planning, and warmup precede timing; BRiDCT reuses its plan and workspace.Source: §4.2 · pp. 6–7
  • Python comparisons include SciPy, DUCC, cached and planned pyFFTW, and 130 OpenCV contracts; their interfaces and absolute timings are evaluated separately.The OpenCV comparison is a practical interface comparison rather than an identical-wrapper or pure native-kernel test.Source: §4.2 · p. 7

5 Timing Results

The final BRiDCT library is faster than tested native comparators across the available primary cases, and its compiled NumPy interface has lower medians than the tested Python alternatives. Development ablations show that staged storage and recursion changes reduce forward times at larger sizes.

5.1 Native performance and implementation choices

  • BRiDCT is faster than FFTW PATIENT and general Ooura in all 18 size/task cases and vDSP in all 15 available cases, with ratios of 1.31–16.04.It also beats the study-specific Ooura SIMD adaptation and its argument-check-free control in all 6/6 small-block cases in each session.Source: Figure 2 · p. 7, §5.1 · p. 7, Table 2 · p. 8

5.2 Python interfaces and their overhead

  • Python gains combine kernel execution with interface and allocation behavior, so they do not rank pure native kernels.The comparisons use complete-call contracts rather than an isolated kernel boundary.Source: §5.2 · p. 8

5.1 Native performance and implementation choices

  • At 128, forward time falls 41.69→22.59→18.61→15.15 µs; at 256, 165.34→104.37→83.36→63.13 µs across banding, padding and blocking.Each effect is conditional on the preceding variant, so the differences are not independent cache, permutation or spill costs.Source: §5.1 · p. 8, Figure 3 · p. 9

5.2 Python interfaces and their overhead

  • BRiDCT has lower medians in all 60 tested cases per interface against SciPy, DUCC, cached pyFFTW and planned pyFFTW, and wins all 130 OpenCV contracts per session.The minimum OpenCV/BRiDCT ratio across contracts in both sessions is 3.90; complete-call measurements include call, allocation, copies and normalization.Source: §5.2 · p. 8, Figure 4 · p. 9

6 Numerical Verification

The verification uses a deterministic float32 corpus and a float64 cosine-matrix oracle to check forward, inverse, and round-trip transforms independently. The final library passes all core checks, while boundary diagnostics expose numerical-range and weak-component limits.

Verification design

  • The 474-array corpus tests all three transform tasks independently against a float64 cosine-matrix oracle, covering cancellation, neighboring floats, weak perturbations, and extreme amplitudes.It contains 420 core cases and 54 boundary diagnostics across six primary sizes.Source: §6 · pp. 9–10

Core and procedural checks

  • All 1 260 core checks pass, with maximum discrepancy 2.251 × 10−7; 2 412 required C checks pass with maximum discrepancy 2.325 × 10−7.The final-library audit reports maximum core error by side length and boundary failures separately; the compiled Python interface also passes 1 152 cases.Source: §6 · p. 10, Table 3 · p. 11

Criterion sensitivity

  • All 78 injected scaling, sign, orientation, and component-loss fault combinations are detected on at least one input, testing criterion sensitivity rather than universal defect detection.Neither intact buffer guards nor this finite corpus proves universal numerical or memory safety.Source: §6 · pp. 10–11

Boundary diagnostics

  • 90 of 162 boundary checks fail, including overflow in unnormalized intermediate stages, subnormal errors, and weak-component loss masked by the global norm criterion.A weak-perturbation case has about 8.95% relative discrepancy in its nonconstant coefficients despite passing the global norm criterion.Source: §6 · p. 10

7 Discussion and Limits

Larger native comparisons show that relative performance varies with size, task, and orientation, while non-power-of-two axes remain unsupported. The measured scope is limited to the tested hardware, interfaces, and warmed workloads.

7.1 Larger squares, rectangles, and unsupported dimensions

  • At 512 × 512, BRiDCT is faster than vDSP in all three tasks in session 2; at 1024 × 1024, it is faster only for the inverse.At 512 × 1024, vDSP leads for the forward transform, and 8 × 64 versus 64 × 8 round trips take 0.597 versus 0.834 µs.Source: §7.1 · p. 11
  • These observations do not establish a ranking for every rectangle.The tested cases include larger squares, reversed rectangles, and selected rectangular shapes.Source: §7.1 · p. 11
  • BRiDCT rejects 32 × 26 and 26 × 32, which FFTW computes; zero-padding is not a substitute because it changes the DCT.This is a compatibility limit, not a numerical failure.Source: §7.1 · p. 11

7.2 Scope and reproducibility

  • The measured gains apply to specified hardware, interfaces, and warmed workloads; input layout, first-use costs, and transform share can change practical benefit.Energy consumption and dynamic memory traffic were not measured.Source: §7.2 · p. 11

8 Future Work

The conclusion reports faster BRiDCT timings across the tested native and compiled-Python comparisons, while future work targets portability, additional comparisons, and larger or unsupported transforms. Supplementary materials provide protocols, numerical diagnostics, code, and recorded measurements.

8 Future Work

  • Multiplatform priorities include a fresh final-release Ubuntu evaluation, large-format Windows performance, Linux ARM64, and additional compilers.CI can verify builds and accuracy, while performance requires controlled-machine measurements.Source: §8 · pp. 11–12
  • Future comparisons target Intel IPP, alternative Accelerate routes, FFTW/SPIRAL-generated kernels, and MATLAB, with precision and MEX call costs accounted for.Large and rectangular transforms, exact non-power-of-two support, numerical-range protection, and end-to-end application and GPU measurements remain targets.Source: §8 · p. 12

9 Conclusion

  • On the tested Apple M3 Max, BRiDCT is approximately 1.3–16.0× faster than tested general-purpose implementations across available forward, inverse, and round-trip cases from 8 × 8 to 256 × 256.Its compiled NumPy interface outperforms OpenCV in every tested call contract, including squares up to 1024 × 1024 and two rectangles, with gains of at least 3.90× including call and allocation costs.Source: §9 · p. 12

Supplementary Material, Code and Data Availability

S1 Native implementation and comparison contracts

The section defines native and Python call contracts and measurement protocols for comparing implementations. It also records an incomplete SPIRAL comparator and cautions that memory interventions are not isolated measurements of cache or spill costs.

S1 Native implementation and comparison contracts

  • The native campaign times one-array calls on a rotating four-array pool, including required copies and output normalization, while excluding planning and allocation.FFTW PATIENT planning was capped at one second per plan; the final policy was fixed before measurement.Source: S1 Native implementation and… · pp. 14–15
  • SPIRAL produced no compiled, independently validated float32 NEON DCT for this protocol, so its comparison was not completed.The attempted recipes failed during rule-tree construction or code generation; later float32 configuration recipes were not run.Source: S1 Native implementation and… · p. 14
  • The memory-intervention campaign changes band buffers, padding, and recursive blocking in sequence, but does not isolate cache misses, spills, or dynamic traffic.Padding adds scratch capacity and copies; the n + 4 stride motivates an intervention but does not establish physical cache indexing.Source: S1 Native implementation and… · p. 15

S2.1 Common configuration

  • The compiled Python interface times warmed repeated calls, including dispatch, checks, prescribed allocation, transform execution, normalization, copies, and allocated-output destruction.Plan creation, imports, input generation, and warmup are excluded; no empty-call cost is subtracted.Source: S2.1 Common configuration · p. 16

S2.2 SciPy, DUCC and pyFFTW

  • SciPy, DUCC, cached pyFFTW, and explicit FFTW plans are compared across 60 interface cases with fresh float32 outputs and separate forward and inverse calls.The round trips materialize the orthonormal spectrum; cached and planned pyFFTW are two interfaces to FFTW, not independent algorithm families.Source: S2.2 SciPy, DUCC and pyFFTW · p. 16

S2.3 OpenCV and output ownership

  • The OpenCV campaign covers 130 call contracts, timing BRiDCT batch processing in a C loop against OpenCV’s per-image calls, with allocation and batching included.OpenCV’s reported thread query does not establish actual worker activity.Source: S2.3 OpenCV and output ownership · p. 16

S3 Larger squares, rectangles, and unsupported dimensions

The larger-shape evaluation extends the final policies beyond the primary sizes and shows that relative performance varies with size and orientation. Unsupported dimensions are rejected rather than padded into a different transform.

S3 Larger squares, rectangles, and unsupported dimensions

  • The final policies support all 64 power-of-two axis pairs from 8 through 1024, with larger-square and rectangle measurements limited to the tested shapes and orientations.The campaign samples eight square shapes and twelve rectangles, including reversed orientations.Source: S3 Larger squares, rectangles, and… · pp. 17–18
  • At 512 × 512, BRiDCT has lower vDSP medians in all three tasks; at 1024 × 1024, it leads in only one.For the 512 × 1024 forward transform, vDSP is faster: 729.50 µs versus 757.20 µs for BRiDCT.Source: S3 Larger squares, rectangles, and… · p. 18, Table S3 · p. 18
  • Padding a length-26 image to 32 would change its DCT, so unsupported dimensions are rejected and excluded from numerical-failure counts.Input layouts and copies required by real applications remain part of later application-level comparisons.Source: S3 Larger squares, rectangles, and… · p. 18

S4 Verification data and test scripts

The verification section combines a deterministic float32 corpus with procedural shape checks and records both required-check performance and boundary behavior. It also supplies scripts, data, and recorded results as public artifacts.

S4 Verification data and test scripts

  • The 474-array corpus has 420 core and 54 boundary inputs, while procedural tests cover all 64 supported shapes; the corpus itself does not include every rectangle.Inputs have deterministic generation, identifiers, and hashes.Source: S4 Verification data and test scripts · p. 18
  • All 1 260 core checks pass with maximum input-relative Frobenius error 2.251 × 10−7; 1 152 additional C shape checks pass, for 2 412 required checks.The compiled Python interface separately passes 1 152 numerical cases with allocated and supplied outputs.Source: S4 Verification data and test scripts · p. 18
  • Boundary diagnostics retain 90 failures among 162 checks, including extreme-range nonfinite results and subnormal errors exceeding the criterion.These boundary cases are excluded from the core accuracy claim, not removed from the dataset.Source: S4 Verification data and test scripts · pp. 18–19
  • The corpus detects each of 78 injected size/direction/fault combinations on at least one input, but this tests criterion sensitivity rather than universal defect detection.The public artifact includes source, build and test scripts, verification data, and recorded comparison results.Source: S4 Verification data and test scripts · p. 19

S5 Small-kernel ablation and portability in progress

The section reports a targeted inverse-kernel refinement and platform evaluations of the current package. The reported comparisons are bounded by specific machines, software routes, and incomplete portability validation.

S5.1 Inverse-16 × 16 refinement

  • The retained inverse-16 × 16 candidate gains more than 2% in both batch sizes and confirmation sessions, with intervals excluding equality.It becomes the final ARM small kernel; its targeted sanitizer check passes 23 472 assertions, but large entry points are stubbed.Source: S5.1 Inverse-16 × 16 refinement · p. 19, Table S4 · p. 19

S5.2 Platform support and validation scope

  • The package passes required checks on the Apple M3 Max and in the independent Windows evaluation, but native Intel/AMD performance and AVX2 execution are not established by the Mac checks.MSVC and Linux ARM64 remain unvalidated, and a fresh Ubuntu evaluation of this exact release is pending.Source: S5.2 Platform support and validation… · p. 20
  • Across two seven-block sessions per comparator family, BRiDCT had lower OpenCV 5.0.0 medians in 116/130 contracts per session; all 104 contracts for six square sizes 8–256 and two rectangles favored BRiDCT.The smallest ratio in session 1 was 1.0008, effectively parity; at 1024 × 1024, only 1/13 then 0/13 contracts favored BRiDCT, and the selected OpenCV DCT route was not profiled.Source: S5.2 Platform support and validation… · p. 20
  • Against SciPy, DUCC, and cached or planned pyFFTW, all 60 Python contracts per reference favor BRiDCT in each session.The separate native C comparison favors BRiDCT over public Ooura in 12/12 cases per session, while optimized study-specific Ooura controls remain faster in several small-block cases.Source: S5.2 Platform support and validation… · p. 20
Loading 2609.28519v1…