EWMA: hoist the filter axis to last before scipy's IIR loop - #216
Merged
Conversation
`_process` handed `scipy.signal.lfilter` whatever axis the message carried. When time is axis 0 of a C-contiguous (time, ch) array -- the layout every acquisition source emits -- each sample step strides `n_ch * itemsize`, and the cost grows superlinearly with the trailing dimension: at 300 samples, going from 256 to 1024 channels is 4x the data but 8.7x the time, versus linear when the axis is already last. Move the axis to the end, filter, and move back. The two copies pay for themselves at every size measured and never lose: 1.04x at 300x256, 1.78x at 30x1024, 4.41x at 300x1024, 1.72x at 3000x1024 -- all including the copies. `filter.py`'s SOS kernel already does this via util/sosfilt_direct, which is why ButterworthZeroPhase measures layout-neutral while this did not. `zi` is one sample slice, so it rides through the same hoist; the stored state stays in the caller's layout for anything that inspects it. AdaptiveStandardScaler runs two EWMAs per message and so paid this twice. At 300x1024 float32 it goes 7.676 -> 3.670 ms (2.09x), and its layout sensitivity drops from 2.64x to 1.28x. On a downstream intracortical feature chain that is 1.51x on total chain cost at 300-sample chunks; just as usefully, run-to-run spread falls from 15-22% to 2-5.7%, because the strided access was making the stage's cost unpredictable rather than merely slow. Output is bit-identical (max|diff| = 0.0) across both layouts, ragged chunk sequences, single-chunk, and float64, checked against a reference that reproduces the old un-hoisted lfilter with the same streaming zi. Deliberately materializes the result rather than returning a transposed view. The view saves a full-size pass and is 12-17% faster in isolation, but loses end to end at the sizes that matter (6% worse at 300 samples, 2% at 1000, 5% better only by 3000) because every downstream op then reads strided. Refs #215.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #215.
EWMATransformer._processhandedscipy.signal.lfilterwhatever axis the message carried. When time is axis 0 of a C-contiguous(time, ch)array -- the layout every acquisition source emits -- each sample step stridesn_ch * itemsize, and scipy's IIR recurrence degrades badly.AdaptiveStandardScalerruns two EWMAs per message (one onx, one onx**2), so it paid this twice.Change
Move the filter axis to the end, run
lfilter(axis=-1), move back.ziis one sample slice so it rides through the same hoist, and the stored state stays in the caller's layout.This is what
filter.py/util/sosfilt_directalready do -- and it is whyButterworthZeroPhasemeasures layout-neutral (341.6 vs 342.6 ms/s in a downstream chain) despite being a 4th-order SOS cascade doing strictly more work than a 1st-order EWMA.Correctness
Bit-identical,
max|diff| = 0.0, verified against a reference that reproduces the old un-hoistedlfilterwith the same streamingzi, across: both layouts, ragged chunk sequences ([7, 13, 300, 1, 1679, 1000]), single-chunk, and float64.New
test_ewma_result_is_independent_of_filter_axis_positionpins that axis-0 and axis-(-1) filtering agree exactly over ragged chunks, for 2-D and 3-D messages. 4037 unit tests pass.Measurements
AdaptiveStandardScalerat 300x1024 float32, time-major:The isolated
lfilterhoist, including both copies:Never loses; wins large where it matters.
Downstream, on an intracortical feature chain at 1024 ch, this is 1.51x on total chain cost at 300-sample chunks. Arguably more valuable for a real-time graph: run-to-run spread drops from 15-22% to 2-5.7% -- the strided access was making the stage's cost unpredictable, not just slow.
Why it materializes instead of returning a view
Returning
np.moveaxis(y, ...)as a view skips a full-size pass and is 12-17% faster in isolation -- the obvious optimization. It loses end to end, because every downstream op then reads strided:Only wins by 3000 samples. Recorded in a code comment so it does not get "optimized" back.
Not included
Cost per second of data still rises with chunk size in the hoisted path (355 -> 375 -> 419 ms/s at 300/1000/3000) where a natively-contiguous one falls (286 -> 285 -> 276): past ~1000 samples the hoist's own copies become the limiting term. Tiling bounds them and wins clearly by 9000 samples (41.8 vs 59.0 ms), losing below ~3000. Where that threshold belongs is a design call, so this PR is untiled -- see #215 for the data. The other two items from #215 (the ~10 full-size temporaries in the scaler's arithmetic tail, and merging the two
lfiltercalls onto a stacked[x, x**2]) are also untouched.