Investigate SIMD vectorization for pairwise distance computation

[stalep]

Summary

Investigate using Java's Vector API (JEP 508) to accelerate the pairwise distance computation in PairDistanceCalculator, which is the innermost hot loop of the eDivisive change-point detection algorithm. Initial analysis suggests 4-8x speedup on the core distance computation loops.

Motivation

The eDivisive algorithm spends the majority of its time computing pairwise distances between data points. The PairDistanceCalculator has three O(n^2) loops that are textbook SIMD candidates:

Loop 1: V (column sums) — computeVH() outer loop

for (int j = 1; j < n; j++) {
    double sum = 0;
    for (int i = 0; i < j; i++) {
        sum += distance(i, j);  // Math.abs(series[i] - series[j])
    }
    V[j - 1] = sum;
}

The inner loop computes sum += |series[i] - series[j]| for all i < j. For a fixed j, series[j] is a constant and the loop iterates contiguous series[i] values. This maps directly to the "common SIMD shape":

1. Broadcast: series[j] into all vector lanes
2. Loop: load series[i..i+lanes] one vector at a time
3. Operate: v.sub(broadcast_j).abs() — subtract and absolute value in parallel
4. Reduce: accumulate into a sum vector, then reduceLanes(ADD) at the end
5. Scalar tail: handle remaining elements

// Vectorized inner loop (conceptual)
DoubleVector broadcastJ = DoubleVector.broadcast(SPECIES_256, series[j]);
DoubleVector sumVec = DoubleVector.zero(SPECIES_256);
int i = 0;
for (; i + lanes <= j; i += lanes) {
    DoubleVector vi = DoubleVector.fromArray(SPECIES_256, series, i);
    sumVec = sumVec.add(vi.sub(broadcastJ).abs());
}
double sum = sumVec.reduceLanes(VectorOperators.ADD);
for (; i < j; i++) sum += Math.abs(series[i] - series[j]); // scalar tail

With AVX2 (4 doubles per vector): ~4x speedup on this inner loop.
With AVX-512 (8 doubles per vector): ~8x speedup.

Loop 2: H (row cumulative sums) — computeVH() second loop

for (int i = 0; i < n - 1; i++) {
    double cumsum = 0;
    for (int j = i; j < n - 1; j++) {
        cumsum += distance(i, j + 1);
        H[i][j] = cumsum;
    }
}

The inner loop is a running sum with stores — harder to vectorize because each H[i][j] depends on the previous cumulative sum. However, the distance(i, j+1) calls can be vectorized (same broadcast-subtract-abs pattern), and the prefix sum can use SIMD prefix-sum techniques.

Loop 3: Column prefix sums — computeColPrefix()

for (int j = 1; j < n; j++) {
    double sum = 0;
    for (int i = 0; i < j; i++) {
        sum += distance(i, j);
        colPrefix[j - 1][i] = sum;
    }
}

Same pattern as Loop 2 — distance computation is vectorizable, prefix sum requires more care.

Loop 4: Q matrix computation — getQVals()

for (int i = 0; i < size; i++) {
    for (int j = i; j < size; j++) {
        // arithmetic on cumsumV, cumsumH values
        Q[i][j] = A - B - C;
    }
}

This loop does arithmetic on precomputed arrays. The inner loop body is more complex but operates on contiguous memory, making partial vectorization possible.

Estimated Impact

For a time series of length N:

• V computation: O(N^2/2) distance calls. With 4x SIMD: ~O(N^2/8) effective operations.
• H computation: O(N^2/2) distance calls + prefix sums.
• Total: for N=1000, that's ~500K distance calls reduced to ~125K vector operations (AVX2).

For h5m's typical use case (change detection on 50-500 datapoints per fingerprint), the absolute time savings may be modest. But for larger series (1000+ points) or batch processing across many fingerprints, the improvement compounds.

The "Common SIMD Shape" (Hashimoto's Pattern)

All vectorizable loops follow the same five steps (reference):

1. Broadcast constants — DoubleVector.broadcast(species, series[j])
2. Loop one vector at a time — for (; i + lanes <= j; i += lanes)
3. Parallel operation — v.sub(broadcast).abs()
4. Reduce the result — sumVec.reduceLanes(ADD) or store to output array
5. Scalar tail — handle remaining 0 to lanes-1 elements with the original loop

The distance function |series[i] - series[j]| (when power=1.0, the common case) maps perfectly to steps 1-3: broadcast series[j], load series[i..i+lanes], subtract, absolute value.

When power != 1.0, Math.pow(diff, power) is harder to vectorize efficiently, but power=1.0 is the default and most common configuration.

Java Vector API Status

The Vector API (JEP 508) is in its 10th incubation round as of JDK 25. Key considerations:

• Requires --add-modules jdk.incubator.vector at runtime
• API may change between JDK versions (incubator status)
• Performance depends on Valhalla for value types (vectors are currently heap-allocated)
• jhunter has no zero-dependency constraint (unlike jjq-core), so using the incubator module is acceptable behind a feature flag

Proposed Approach

1. Create a VectorPairDistanceCalculator that extends or replaces PairDistanceCalculator when the Vector API is available
2. Feature flag: detect Vector API availability at runtime via Module API or class loading; fall back to scalar implementation
3. Benchmark: use the existing PairDistanceCalculatorBenchmark to compare scalar vs vectorized on series of length 100, 500, 1000, 5000
4. Start with Loop 1 (V column sums) — cleanest vectorization, highest impact per line of code changed

Acceptance Criteria

• Vectorized computeVH() for the power=1.0 (default) case
• Benchmark showing speedup vs scalar baseline at N=100, N=500, N=1000
• Runtime fallback to scalar when Vector API is not available
• All existing tests pass with both scalar and vectorized implementations
• No regression for the power != 1.0 case (stays scalar)

References

• Mitchell Hashimoto: "Everyone Should Know SIMD" — the common 5-step pattern
• JEP 508: Vector API (Tenth Incubator)
• simdjson — advanced SIMD techniques (not needed here; the common pattern suffices)
• jjq perf stat data (2026-07-02): confirmed branch misprediction is not a bottleneck; SIMD benefit comes from data parallelism, not branch elimination