Repository navigation
perf(base): speed up qqmul and qvmul - #215
Merged
taughz merged 1 commit intoOct 13, 2026
Merged
Conversation
Both used np.dot/np.cross on 3-vectors, whose generic dispatch overhead - built for arbitrary shapes/broadcasting - dominates cost at this size, same pattern as isR/trnorm/tr2adjoint (rai-opensource#213, rai-opensource#214). qqmul: replaced with explicit scalar Hamilton-product arithmetic, bit-identical to the prior result (verified to ~1 ULP over 2000 random trials). qvmul: replaced the q * pure(v) * conj(q) sandwich (two full Hamilton products, each wasting work on a zero scalar part, via qqmul/qpure/ qconj with their own getvector re-validation) with the closed-form rotation identity v' = v + 2s(w x v) + 2 w x (w x v) for q = (s, w). This is a different, well-known equivalent formula rather than a re-expression of the same one, so it is not bit-identical - verified to ~1e-14 absolute over 2000 random trials (v scaled to magnitude ~10), well within floating-point noise. qqmul ~10x faster (15.0us -> 1.4us isolated), qvmul ~24x faster (33.4us -> 1.4us isolated). End to end: Q1 * v drops from ~33us to ~4us (~9x), Q1 * Q2 from ~19us to ~10us (~2x - qqmul itself is no longer the bottleneck there; the remainder is UnitQuaternion construction overhead, in particular qunit()'s np.linalg.norm/np.r_ calls, which is a separate, un-addressed candidate for a future fix in the same spirit). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
|
Codecov Report✅ All modified and coverable lines are covered by tests. 📢 Thoughts on this report? Let us know! |
3 tasks done
17 tasks done
taughz
approved these changes
Oct 13, 2026
taughz
pushed a commit
that referenced
this pull request
Oct 13, 2026
…#216) ## Summary `trexp`/`trlog` are much slower than `rodrigues` (~20us and ~17us vs ~7us), but unlike the previous three perf PRs (#213, #214, #215) they don't have one dominant wasteful generic-NumPy call to fix - they're deep call chains (`trexp` -> `isskewa` -> `vexa` -> `iszerovec` -> `unittwist_norm` -> `rodrigues` -> `skew` -> `rt2tr` -> `ishom` -> `isR`) where cost is spread thin across many small layers, each independently re-validating/re-extracting. This PR pulls out the two safe, verified wins found in that chain; it does **not** close the full gap to `rodrigues` - see "further speedup opportunities" below. - **`isskewa`** (the validity check `trexp` runs on its se(3) input) had the exact same bug `ishom` had before #213: `np.linalg.norm()` on a small fixed-size matrix, plus an `all(S[-1,:] == 0)` array-allocation-and-compare for checking the bottom row is zero. Fixed the same way (explicit sum-of-squares residual + scalar comparisons). ~1.9x faster standalone. - **`np.eye(3)`** was rebuilt from scratch on every call in `rodrigues` (1x), `trexp`'s `V`-matrix construction (1x), and `trlog` (2x) - despite only ever being used as a read-only operand in an addition/multiplication that produces a new array. Replaced with a module-level constant `_EYE3` at those call sites only. **Not** used at the four zero-motion/zero-rotation early-return `np.eye(N)` calls in `rodrigues`/`trexp` (e.g. `if iszerovec(w): return np.eye(3)`) - those hand the array object directly to the caller, and aliasing them to a shared constant would risk a real bug if a caller ever mutated the returned matrix in place. Left as fresh arrays. - **`isskew`** (the plain, non-augmented so(n) check) was also tried with the same explicit-arithmetic treatment as `isskewa`, but did **not** show a reliable win under careful (min-of-repeats) benchmarking - the only cost it avoids is `np.linalg.norm` on an already-cheap `S + S.T`, not enough margin to reliably beat by hand-unrolling the indexing. Reverted to the original implementation rather than keep an unproven "fix" - same judgment call as an earlier fully-scalar `isR` variant that was tried and discarded for the same reason. Net effect, measured old-vs-new in the same process, min of 9 repeats each (single-run `timeit` numbers were misleading here - see commit message): **isskewa ~1.9x, rodrigues ~1.13x, trexp ~1.05x, trlog ~1.10x**. Real, but modest compared to #213-#215. ## Further speedup opportunities (not attempted here) A bigger win on `trexp`/`trlog` would mean collapsing the redundant `getvector`/re-validation layers in the hot path into a more direct, fused implementation - e.g. `trexp(se3_matrix)` currently re-extracts and re-validates the same data through `isskewa` -> `vexa` -> `iszerovec` -> `unittwist_norm` independently, each with its own `getvector` call. That's a larger, riskier change: this code has clearly-deliberate handling for numerically delicate edge cases (near-zero rotation, near-pi singularity, caller-supplied `theta`), so a fused rewrite needs much more careful edge-case test coverage than the drop-in fixes here. Not scoped or attempted in this PR. Two smaller, unrelated candidates surfaced but out of scope for this PR: - `np.linalg.norm` is used directly (instead of the project's own `norm()`, which is faster for small vectors per its own docstring) in several other places - `iszerovec`, `isunitvec`, `qnorm`/`qunit` in `quaternions.py`, and a few spots in `geom2d.py`/`geom3d.py`/`pose3d.py`. Worth noting: `iszerovec` was checked directly as part of this investigation, and swapping it to the project's `norm()` did **not** show a reliable win either (consistent with the project's own timing script, which shows `base.norm(R6)` measuring *slower* than `np.linalg.norm(R6)` for a 6-vector) - so this isn't a guaranteed win either, just an inconsistency worth someone measuring properly before touching. - `np.eye(N)` is rebuilt from scratch in ~20 other places across the codebase (`pose2d.py`, `pose3d.py`, `transforms2d.py`, `transformsNd.py`, and other unrelated functions in `transforms3d.py`). None are in the `trexp`/`trlog`/`rodrigues` path touched here; a broader sweep would need the same operand-vs-return-value aliasing care taken in this PR, function by function. ## Test plan - [x] `pytest tests/` - 343 passed - [x] `isskewa` verified against the prior implementation over 3000 randomized trials (skew/non-skew, near-boundary perturbations) - [x] `rodrigues`, `trexp`, `trlog` (both matrix and twist-vector return forms) verified bit-identical (max diff 0.0) against reference copies of the prior implementations over 1000-2000 randomized trials each 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
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
3 participants
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.
Summary
qqmul(Hamilton product) andqvmul(vector rotation) both usednp.dot/np.crosson 3-vectors, whose generic dispatch overhead - built for arbitrary shapes/broadcasting - dominates cost at this size. Same pattern asisR/trnorm/tr2adjointin perf(base): speed up isR/ishom/ishom2 orthogonality checks #213/perf(base): speed up trnorm and tr2adjoint #214.qqmul: replaced with explicit scalar Hamilton-product arithmetic. Bit-identical to the prior result (verified to ~1 ULP over 2000 random trials).qvmul: replaced theq * pure(v) * conj(q)sandwich (two full Hamilton products throughqqmul/qpure/qconj, each with its owngetvectorre-validation, and each Hamilton product wasting work on a zero scalar part) with the closed-form rotation identityv' = v + 2s(w x v) + 2 w x (w x v)forq = (s, w). This is a genuinely different (though mathematically equivalent) formula, not just a re-expression of the same one - so it is not bit-identical, verified instead to ~1e-14 absolute over 2000 random trials (withvscaled to magnitude ~10), well within floating-point noise.qqmul~10x faster (15.0us -> 1.4us isolated),qvmul~24x faster (33.4us -> 1.4us isolated). End to end:Q1 * vdrops from ~33us to ~4us (~9x).Q1 * Q2drops from ~19us to ~10us (~2x) -qqmulitself is no longer the bottleneck there; the remainder isUnitQuaternionconstruction overhead (in particularqunit()'snp.linalg.norm/np.r_calls), a separate, un-addressed candidate for a possible future fix in the same spirit, not attempted here.Test plan
pytest tests/- 343 passedqqmulverified against the prior implementation over 2000 random quaternion pairs (max diff ~1 ULP, i.e. bit-identical modulo floating-point operation ordering)qvmulverified against the prior implementation over 2000 random (quaternion, vector) pairs, vector magnitude ~10 (max abs diff ~1.4e-14)🤖 Generated with Claude Code