Benchmark: freqresp! optimization strategies (@turbo, threading) for SISO 48-state x 150 freqs

[baggepinnen]

Benchmarking freqresp! optimization strategies

Benchmark of strategies for maximizing performance of freqresp! on a representative use case: a 48-state SISO continuous state-space system evaluated over 150 frequencies. A self-contained benchmark environment was used; the package's own freqresp! serves as the correctness oracle and no package source was modified for the comparison.

Variants

ID	Strategy
V0	ControlSystemsBase.freqresp! (oracle, full call incl. factorization)
V1	serial hand-rolled loop, original ldiv2! (loop-only)
V2	Threads.@spawn over contiguous frequency chunks
V2b	Polyester.@batch over chunks
V2c	OhMyThreads.@tasks with TaskLocalValue buffers
V3	serial loop, @turbo ldiv2! (real/imag split)
V4	@spawn + @turbo ldiv2!
V4b	Polyester.@batch + @turbo ldiv2!

All variants operate on pre-factored matrices (A=F.H, C=complex.(sys.C*Q), B=Q\sys.B, D) so the "loop-only" timing isolates the per-frequency strategy from the one-time Hessenberg factorization. All pass isapprox(·, ·; rtol=1e-10) against the package baseline.

Results (24 logical cores)

Loop-only time, µs (lower is better):

Variant	-t 4	-t 8	-t 24
V1 serial	545	566	565
V3 turbo	360	372	372
V2 threaded	147	90	106
V2b polyester	203	108	95
V2c omt	154	96	127
V4 threaded+turbo	104	78	101
V4b polyester+turbo	142	80	93

Full-call baseline (V0, factorization included): ~625–655 µs.

Findings

1. The per-frequency loop, not the factorization, dominates. The serial loop is ~565 µs while the one-time Hessenberg factorization is only ~85 µs (full call ≈ loop + factorization ≈ 650 µs). Optimizing the loop is therefore worthwhile.

2. @turbo is a solid ~1.5× win (565 → ~370 µs serial), more than expected for such short loops, with no threads and negligible accuracy loss (~1e-14 vs baseline).

3. Threading gives ~6–7× but scales sub-linearly: at this size 8 threads beats 24 (150 freqs / 24 threads ≈ 6 freqs each — task overhead and oversubscription dominate). Sweet spot ~8 threads.

4. Best overall is threading + @turbo (~78 µs at 8 threads, a ~7.2× speedup over serial). @spawn+@turbo (V4) is consistently near-best across thread counts; Polyester+@turbo (V4b) edges ahead only at 24 threads.

5. Scheduler choice is secondary. @spawn, Polyester, and OhMyThreads are within ~1.4× of each other; @spawn is the most robust across thread counts and adds no dependency.

The @turbo / complex caveat

LoopVectorization's @turbo does not support Complex. A bare @turbo on the complex ldiv2! loops will not compile, and a naive reinterpret-to-real is silently wrong (complex multiply has re/im cross terms). The @turbo variant therefore carries the state as separate real/imag Float64 arrays and writes out the complex arithmetic by hand, so the two hot j = 1:k-2 loops (X update, u update) become valid real @turbo loops. Result matches the baseline to ~1e-14.

Recommendation

The lowest-risk improvement is dropping in the real/imag-split @turbo ldiv2! (≈1.5× for free, no new threading semantics). If the loop is a known hotspot, add an opt-in threaded path (@spawn over chunks with per-task buffers) — combined with @turbo it reaches ~7× — but keep the default serial, since threading helps only for many frequencies and benefits from a bounded thread count.

Caveats

Single-machine numbers (24 cores, JULIA_EXCLUSIVE=1); the threading sweet spot shifts with core count, and absolute times vary a few % with the random ssrand system.

[baggepinnen]

Follow-up: @tturbo (LoopVectorization built-in threading)

Tested @tturbo (the threaded variant of @turbo) on the two hot ldiv2! j-loops, with a serial outer frequency loop (variant V3t).

It is slower than plain @turbo:

Variant	-t 8	-t 24
V3 turbo	368	357
V3t tturbo	412	411

@tturbo threads within a single loop, but the inner loops here are only k-2 ≤ 46 elements long and sit inside the sequentially-dependent outer k loop, so per-loop thread-dispatch overhead dominates and there is no win to be had. The exploitable parallelism is the 150-frequency outer loop, which @spawn/Polyester exploit (~78–94 µs) and @tturbo structurally cannot reach. Correctness matches baseline (~1e-14).

Conclusion unchanged: use plain @turbo for vectorization, and thread the outer frequency loop separately when many frequencies are evaluated.

[baggepinnen]

Follow-up: Reactant.jl, and the modal/pole-residue reformulation

Reactant traces Julia into MLIR/XLA and wants vectorized array programs. The custom Givens-based ldiv2! (scalar indexing, data-dependent control flow, givensAlgorithm) does not trace into efficient XLA. To use Reactant at all, the computation must first be rewritten as batched array ops — and the natural rewrite is the modal / pole-residue form: eigendecompose A = V Λ V⁻¹ once, then

R(ω) = D + Σⱼ resⱼ / (iω − λⱼ),   resⱼ = (C·V)ⱼ · (V⁻¹·B)ⱼ

Every frequency is then a pure broadcast + reduction over an nx × nfreq grid — no per-frequency linear solve. Benchmarked in plain Julia:

Variant	µs (8 threads)	vs serial Hessenberg
serial Hessenberg (current)	546	1×
best Hessenberg micro-opt (threaded+@turbo)	73	7.5×
modal serial	42	13×
modal broadcast	59	9×
modal threaded	10	54×

Accuracy matched the baseline to ~1e-14 on this system (cond(V)=147).

Verdict on Reactant for this problem: not worth it. The data is tiny (48×150 ≈ 7200 complex numbers) and the modal form already runs in ~10 µs on CPU — below the per-call XLA dispatch / GPU kernel-launch latency (tens of µs). Reactant would add a heavy dependency and long compile times for no gain at this size. It would only pay off at much larger scale (thousands of states, very many frequencies, or batches of systems) or on GPU — and there the same modal array program is what you would feed it.

Caveat: the modal form trades robustness for speed — it relies on the eigendecomposition, so accuracy degrades with cond(V) for non-normal / near-defective A. That conditioning risk is exactly why the package defaults to the Hessenberg solve. Modal is an excellent fast path for well-conditioned systems, not a safe drop-in default.

Takeaway: the algorithm matters far more than the micro-optimizations here. The biggest win is not @turbo/threading the Hessenberg solve (~7×) but switching to the modal formulation for well-conditioned systems (~13–54×). Reactant is the wrong tool at this scale; its sweet spot (large batched array programs on GPU) is reached via the same modal reformulation, not by tracing the current solver.