Skip to content

Commit 111ee8e

Browse files
petercorkeclaude
andauthored
perf(base): speed up qqmul and qvmul (#215)
## Summary - `qqmul` (Hamilton product) and `qvmul` (vector rotation) 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` in #213/#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 through `qqmul`/`qpure`/`qconj`, each with its own `getvector` re-validation, and each Hamilton product wasting work on a zero scalar part) 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 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 (with `v` scaled to magnitude ~10), well within floating-point noise. - Net effect: `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` drops 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), a separate, un-addressed candidate for a possible future fix in the same spirit, not attempted here. ## Test plan - [x] `pytest tests/` - 343 passed - [x] `qqmul` verified against the prior implementation over 2000 random quaternion pairs (max diff ~1 ULP, i.e. bit-identical modulo floating-point operation ordering) - [x] `qvmul` verified against the prior implementation over 2000 random (quaternion, vector) pairs, vector magnitude ~10 (max abs diff ~1.4e-14) 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
1 parent abe79e4 commit 111ee8e

1 file changed

Lines changed: 30 additions & 7 deletions

File tree

‎spatialmath/base/quaternions.py‎

Lines changed: 30 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -326,12 +326,19 @@ def qqmul(q1: ArrayLike4, q2: ArrayLike4) -> QuaternionArray:
326326
"""
327327
q1 = smb.getvector(q1, 4)
328328
q2 = smb.getvector(q2, 4)
329-
s1 = q1[0]
330-
v1 = q1[1:4]
331-
s2 = q2[0]
332-
v2 = q2[1:4]
329+
s1, x1, y1, z1 = q1
330+
s2, x2, y2, z2 = q2
333331

334-
return np.r_[s1 * s2 - np.dot(v1, v2), s1 * v2 + s2 * v1 + np.cross(v1, v2)]
332+
# explicit scalar arithmetic avoids the generic-dispatch overhead of
333+
# np.dot/np.cross on a 3-vector, which dominates cost at this size
334+
return np.array(
335+
[
336+
s1 * s2 - x1 * x2 - y1 * y2 - z1 * z2,
337+
s1 * x2 + x1 * s2 + y1 * z2 - z1 * y2,
338+
s1 * y2 - x1 * z2 + y1 * s2 + z1 * x2,
339+
s1 * z2 + x1 * y2 - y1 * x2 + z1 * s2,
340+
]
341+
)
335342

336343

337344
def qinner(q1: ArrayLike4, q2: ArrayLike4) -> float:
@@ -398,8 +405,24 @@ def qvmul(q: ArrayLike4, v: ArrayLike3) -> R3:
398405
"""
399406
q = smb.getvector(q, 4)
400407
v = smb.getvector(v, 3)
401-
qv = qqmul(q, qqmul(qpure(v), qconj(q)))
402-
return qv[1:4]
408+
s, x, y, z = q
409+
vx, vy, vz = v
410+
411+
# closed-form v' = v + 2s(w x v) + 2 w x (w x v), for q = (s, w) unit;
412+
# mathematically equivalent to q * pure(v) * conj(q) but avoids two full
413+
# Hamilton products (each wasting work on a zero scalar part) and the
414+
# np.cross/np.dot dispatch overhead within them
415+
tx = 2 * (y * vz - z * vy)
416+
ty = 2 * (z * vx - x * vz)
417+
tz = 2 * (x * vy - y * vx)
418+
419+
return np.array(
420+
[
421+
vx + s * tx + (y * tz - z * ty),
422+
vy + s * ty + (z * tx - x * tz),
423+
vz + s * tz + (x * ty - y * tx),
424+
]
425+
)
403426

404427

405428
def vvmul(qa: ArrayLike3, qb: ArrayLike3) -> R3:

0 commit comments

Comments
 (0)