diff --git a/src/rqm_optimize/fusion.py b/src/rqm_optimize/fusion.py
index d76bc64..4a23c6a 100644
--- a/src/rqm_optimize/fusion.py
+++ b/src/rqm_optimize/fusion.py
@@ -11,7 +11,13 @@
 
 import numpy as np
 
-from .geometry import remove_global_phase
+from .geometry import (
+    quaternion_canonicalize,
+    quaternion_multiply,
+    quaternion_to_su2,
+    remove_global_phase,
+    su2_to_quaternion,
+)
 from .qiskit_adapter import extract_matrix, is_fuseable_single_qubit
 
 if TYPE_CHECKING:
@@ -107,24 +113,35 @@ def flush_all() -> None:
 
 
 def _fuse_matrices(run: list["CircuitInstruction"]) -> np.ndarray:
-    """Multiply gate matrices left-to-right (circuit order) and normalise.
+    """Accumulate a run of single-qubit gates as quaternion products on S³.
 
-    In quantum circuit convention the first gate applied is the leftmost
-    matrix.  We accumulate U_n @ ... @ U_1 by applying right-to-left
-    composition so that the overall effect is the same as applying gates
-    sequentially.
+    Each ``CircuitInstruction`` in *run* is inspected for its 2×2 unitary
+    matrix, which is then mapped to a unit quaternion via the exact SU(2)
+    isomorphism.  The quaternions are multiplied in circuit order (gate[0]
+    applied first), and the canonical accumulated result is converted back to
+    a 2×2 SU(2) matrix.
+
+    Using quaternion multiplication instead of raw matrix multiplication keeps
+    intermediate results on the unit 3-sphere S³, avoids accumulated complex
+    phase drift, and produces a canonical output via
+    :func:`~rqm_optimize.geometry.quaternion_canonicalize`.
+
+    Circuit order: gate[0] is applied first.  Quaternion product ``q_n * … *
+    q_1`` is accumulated with the last gate on the left (outer position),
+    matching the matrix convention ``U_n @ … @ U_1``.
 
     Args:
         run: List of ``CircuitInstruction`` objects in circuit order.
 
     Returns:
-        2×2 SU(2)-normalised combined unitary matrix.
+        2×2 SU(2)-normalized combined unitary matrix.
     """
-    # Start with identity and compose in circuit order.
-    # circuit order: gate[0] first, gate[1] second, ...
-    # combined unitary = gate[-1] @ ... @ gate[1] @ gate[0]
-    combined = np.eye(2, dtype=complex)
+    # Identity quaternion: q = (1, 0, 0, 0).
+    q_accum = np.array([1.0, 0.0, 0.0, 0.0])
     for instr in run:
         mat = extract_matrix(instr.operation)
-        combined = mat @ combined
-    return remove_global_phase(combined)
+        q_gate = su2_to_quaternion(remove_global_phase(mat))
+        # Apply q_gate after the current accumulation: q_gate * q_accum.
+        q_accum = quaternion_multiply(q_gate, q_accum)
+    q_accum = quaternion_canonicalize(q_accum)
+    return quaternion_to_su2(q_accum)
diff --git a/src/rqm_optimize/geometry.py b/src/rqm_optimize/geometry.py
index 58a5b83..50bfea7 100644
--- a/src/rqm_optimize/geometry.py
+++ b/src/rqm_optimize/geometry.py
@@ -2,6 +2,28 @@
 
 Keeps math helpers small and rigorous.  Does not duplicate rqm-core; this
 module only contains the minimum needed to support circuit-level optimization.
+
+Quaternion convention
+---------------------
+A unit quaternion is represented as a length-4 NumPy array ``q = [w, x, y, z]``
+where *w* is the scalar part and *x, y, z* are the pure-imaginary components.
+
+The exact isomorphism between unit quaternions and SU(2) matrices used here is::
+
+    q = (w, x, y, z)  <->  U = [[ w - iz,  -(y + ix)],
+                                  [ y - ix,    w + iz ]]
+
+This convention is consistent with the standard quantum rotation gates:
+
+* ``Rx(theta)``  ->  ``(cos(theta/2), sin(theta/2), 0,            0          )``
+* ``Ry(theta)``  ->  ``(cos(theta/2), 0,            sin(theta/2), 0          )``
+* ``Rz(theta)``  ->  ``(cos(theta/2), 0,            0,            sin(theta/2))``
+
+Quaternion multiplication ``quaternion_multiply(q2, q1)`` represents applying
+*q1* first, then *q2*, matching matrix composition ``U2 @ U1``.
+
+For canonicalization the representative with non-negative scalar part (``w >= 0``)
+is preferred, keeping the quaternion on the shortest geodesic arc on S³.
 """
 
 from __future__ import annotations
@@ -69,3 +91,196 @@ def matrices_close(
     n = a.shape[0]
     inner = np.trace(a.conj().T @ b)
     return bool(abs(abs(inner) - n) < atol * n)
+
+
+# ---------------------------------------------------------------------------
+# Quaternion helpers
+# ---------------------------------------------------------------------------
+
+def su2_to_quaternion(matrix: NDArray[np.complexfloating]) -> NDArray[np.floating]:
+    """Return the unit quaternion ``[w, x, y, z]`` for an SU(2) matrix.
+
+    Uses the exact algebraic isomorphism::
+
+        U = [[ w - iz,  -(y + ix)],
+             [ y - ix,    w + iz ]]  <->  q = (w, x, y, z)
+
+    The result is normalized and has ``w >= 0`` (canonical shortest-path
+    representative on S³).
+
+    Args:
+        matrix: 2×2 complex SU(2) array (determinant +1, unitary).
+
+    Returns:
+        Length-4 float array ``[w, x, y, z]`` representing the unit quaternion.
+    """
+    # Extract components from matrix entries:
+    #   U[0,0] = w - iz  =>  w = Re(U[0,0]),  z = -Im(U[0,0])
+    #   U[0,1] = -(y+ix) =>  y = -Re(U[0,1]), x = -Im(U[0,1])
+    w = float(matrix[0, 0].real)
+    x = float(-matrix[0, 1].imag)
+    y = float(-matrix[0, 1].real)
+    z = float(-matrix[0, 0].imag)
+    q = np.array([w, x, y, z], dtype=float)
+    return quaternion_canonicalize(q)
+
+
+def quaternion_to_su2(q: NDArray[np.floating]) -> NDArray[np.complexfloating]:
+    """Return the SU(2) matrix for a unit quaternion ``[w, x, y, z]``.
+
+    Inverse of :func:`su2_to_quaternion`.  Uses the isomorphism::
+
+        q = (w, x, y, z)  <->  U = [[ w - iz,  -(y + ix)],
+                                      [ y - ix,    w + iz ]]
+
+    Args:
+        q: Length-4 float array ``[w, x, y, z]``.
+
+    Returns:
+        2×2 complex SU(2) array.
+    """
+    w, x, y, z = float(q[0]), float(q[1]), float(q[2]), float(q[3])
+    return np.array(
+        [
+            [complex(w, -z), complex(-y, -x)],
+            [complex(y, -x), complex(w,  z)],
+        ],
+        dtype=complex,
+    )
+
+
+def quaternion_multiply(
+    q2: NDArray[np.floating],
+    q1: NDArray[np.floating],
+) -> NDArray[np.floating]:
+    """Return the quaternion product ``q2 * q1``.
+
+    Represents applying gate *q1* first, then gate *q2*, matching matrix
+    composition ``U2 @ U1``.
+
+    Uses the standard Hamilton product formula::
+
+        (w2, x2, y2, z2) * (w1, x1, y1, z1) = (
+            w2*w1 - x2*x1 - y2*y1 - z2*z1,
+            w2*x1 + x2*w1 + y2*z1 - z2*y1,
+            w2*y1 - x2*z1 + y2*w1 + z2*x1,
+            w2*z1 + x2*y1 - y2*x1 + z2*w1,
+        )
+
+    Args:
+        q2: Unit quaternion for the second (outer) gate.
+        q1: Unit quaternion for the first (inner) gate.
+
+    Returns:
+        Product quaternion (not yet canonicalized or re-normalized).
+    """
+    w2, x2, y2, z2 = float(q2[0]), float(q2[1]), float(q2[2]), float(q2[3])
+    w1, x1, y1, z1 = float(q1[0]), float(q1[1]), float(q1[2]), float(q1[3])
+    return np.array(
+        [
+            w2 * w1 - x2 * x1 - y2 * y1 - z2 * z1,
+            w2 * x1 + x2 * w1 + y2 * z1 - z2 * y1,
+            w2 * y1 - x2 * z1 + y2 * w1 + z2 * x1,
+            w2 * z1 + x2 * y1 - y2 * x1 + z2 * w1,
+        ],
+        dtype=float,
+    )
+
+
+def quaternion_canonicalize(q: NDArray[np.floating]) -> NDArray[np.floating]:
+    """Return the canonical unit quaternion with non-negative scalar part.
+
+    Both *q* and *-q* represent the same physical SU(2) rotation.  This
+    function normalizes the quaternion and selects the representative with
+    ``w >= 0``, keeping it on the shortest geodesic arc (canonical
+    representative) on S³.
+
+    Args:
+        q: Length-4 float array ``[w, x, y, z]``.
+
+    Returns:
+        Normalized quaternion with ``w >= 0``.
+    """
+    norm = float(np.linalg.norm(q))
+    if norm < _ATOL:
+        return np.array([1.0, 0.0, 0.0, 0.0], dtype=float)
+    q = q / norm
+    if q[0] < 0.0:
+        q = -q
+    return q
+
+
+def quaternion_to_axis_angle(
+    q: NDArray[np.floating],
+) -> tuple[NDArray[np.floating], float]:
+    """Extract the rotation axis and angle from a unit quaternion.
+
+    The physical Bloch-sphere rotation angle is ``theta = 2 * arccos(w)``
+    (the factor of 2 arises from the half-angle spinor convention).
+
+    Args:
+        q: Length-4 float array ``[w, x, y, z]``, assumed to be a unit
+           quaternion.  Canonicalization is applied internally.
+
+    Returns:
+        A tuple ``(axis, theta)`` where:
+
+        * ``axis`` is a length-3 float array giving the unit rotation axis
+          ``[nx, ny, nz]``.  When the rotation is near-identity (``theta``
+          close to 0), ``axis`` defaults to ``[0, 0, 1]``.
+        * ``theta`` is the rotation angle in radians (in ``[0, pi]`` after
+          canonicalization).
+    """
+    q = quaternion_canonicalize(q)
+    w = float(np.clip(q[0], -1.0, 1.0))
+    theta = 2.0 * float(np.arccos(w))
+    sin_half = float(np.sqrt(max(0.0, 1.0 - w * w)))
+    if sin_half < _ATOL:
+        axis = np.array([0.0, 0.0, 1.0], dtype=float)
+    else:
+        axis = np.array(q[1:], dtype=float) / sin_half
+    return axis, theta
+
+
+# Reference unit vectors for the three Cartesian rotation axes.
+_CARTESIAN_AXES: tuple[tuple[str, NDArray[np.floating]], ...] = (
+    ("x", np.array([1.0, 0.0, 0.0])),
+    ("y", np.array([0.0, 1.0, 0.0])),
+    ("z", np.array([0.0, 0.0, 1.0])),
+)
+
+# Default tolerance for axis-alignment checks.
+_AXIS_ATOL = 1e-6
+
+
+def axis_aligned_rotation(
+    q: NDArray[np.floating],
+    atol: float = _AXIS_ATOL,
+) -> tuple[str, float] | None:
+    """Return ``(axis_name, theta)`` if *q* is a rotation about x, y, or z.
+
+    Checks whether the rotation axis of *q* lies within *atol* of the x, y, or
+    z unit vector.  This is the geometric test that drives the axis-aware
+    compression pass: when the result is non-``None``, the caller can emit a
+    single named rotation gate (``rx``, ``ry``, or ``rz``) instead of a
+    generic ``U`` gate.
+
+    Near-identity rotations (``|theta| < atol``) return ``None`` because the
+    axis is numerically undefined at that scale.
+
+    Args:
+        q: Length-4 float array ``[w, x, y, z]``, a unit quaternion.
+        atol: Absolute tolerance for axis alignment and near-identity detection.
+
+    Returns:
+        ``("x"/"y"/"z", theta)`` when the axis is Cartesian-aligned, or
+        ``None`` when the rotation is generic or near-identity.
+    """
+    q = quaternion_canonicalize(q)
+    axis, theta = quaternion_to_axis_angle(q)
+    if abs(theta) < atol:
+        return None
+    for name, ref in _CARTESIAN_AXES:
+        if np.allclose(axis, ref, atol=atol):
+            return name, theta
+    return None
diff --git a/src/rqm_optimize/qiskit_adapter.py b/src/rqm_optimize/qiskit_adapter.py
index 3c182c1..d60c6bc 100644
--- a/src/rqm_optimize/qiskit_adapter.py
+++ b/src/rqm_optimize/qiskit_adapter.py
@@ -178,6 +178,48 @@ def _try_extract_matrix(op: "Instruction") -> np.ndarray | None:
     return None
 
 
+def emit_axis_aligned_gate(
+    matrix: np.ndarray,
+    circuit: "QuantumCircuit",
+    qubit: "Qubit",
+) -> int:
+    """Emit a single Rx/Ry/Rz gate if *matrix* is a Cartesian-axis rotation.
+
+    Inspects the quaternion representation of *matrix* to determine whether the
+    fused rotation is aligned with the x, y, or z axis on S³.  When it is, a
+    single named rotation gate (``rx``, ``ry``, or ``rz``) is appended to
+    *circuit* and 1 is returned.  When the rotation is generic or near-identity,
+    the circuit is left unchanged and 0 is returned.
+
+    This is the axis-aware compression pass: it operates directly on the S³
+    quaternion representation rather than on matrices, and produces a semantically
+    named gate that is directly executable on hardware (in particular, ``rz`` is
+    virtual/"free" on many backends).
+
+    Args:
+        matrix: 2×2 SU(2)-normalized unitary matrix for the fused run.
+        circuit: Target ``QuantumCircuit`` to append the gate to.
+        qubit: The target qubit.
+
+    Returns:
+        1 if an axis-aligned gate was emitted, 0 otherwise.
+    """
+    from .geometry import axis_aligned_rotation, su2_to_quaternion
+
+    q = su2_to_quaternion(matrix)
+    result = axis_aligned_rotation(q)
+    if result is None:
+        return 0
+    axis_name, theta = result
+    if axis_name == "x":
+        circuit.rx(theta, qubit)
+    elif axis_name == "y":
+        circuit.ry(theta, qubit)
+    else:
+        circuit.rz(theta, qubit)
+    return 1
+
+
 def emit_euler_gate(
     unitary: np.ndarray,
     circuit: "QuantumCircuit",
@@ -190,7 +232,7 @@ def emit_euler_gate(
     produce a compact, exact decomposition.
 
     Args:
-        unitary: 2×2 complex SU(2)-normalised unitary matrix.
+        unitary: 2×2 complex SU(2)-normalized unitary matrix.
         circuit: The ``QuantumCircuit`` to append gates to.
         qubit: The target qubit.
         basis: Euler decomposer basis string (e.g. ``"U"``, ``"ZSX"``,
@@ -246,9 +288,19 @@ def build_optimized_circuit(
             matrix: np.ndarray = seg["matrix"]
             original_count: int = seg["original_count"]
 
-            # Build a temporary single-qubit sub-circuit to count decomposition.
             from qiskit import QuantumCircuit as QC
 
+            # First pass: axis-aware compression.
+            # If the fused rotation is aligned with the x, y, or z axis, emit
+            # a single rx/ry/rz gate directly from the S³ quaternion form.
+            # This always produces 1 gate (≤ original_count which is ≥ 2).
+            tmp = QC(1)
+            if emit_axis_aligned_gate(matrix, tmp, tmp.qubits[0]) > 0:
+                for instr in tmp.data:
+                    out.append(instr.operation, [qubit])
+                continue
+
+            # Second pass: Euler decomposer for generic (non-axis-aligned) rotations.
             tmp = QC(1)
             emitted = emit_euler_gate(matrix, tmp, tmp.qubits[0], basis=basis)
 
diff --git a/tests/test_axis_aware.py b/tests/test_axis_aware.py
new file mode 100644
index 0000000..dca92d6
--- /dev/null
+++ b/tests/test_axis_aware.py
@@ -0,0 +1,242 @@
+"""tests/test_axis_aware.py — Tests for axis-aware single-qubit compression.
+
+Validates that fused same-axis rotation runs emit named rx/ry/rz gates
+(rather than a generic U gate) and that the resulting circuits are equivalent
+to their inputs.
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+from qiskit import QuantumCircuit
+from qiskit.quantum_info import Operator
+
+from rqm_optimize import optimize
+from rqm_optimize.geometry import axis_aligned_rotation, quaternion_canonicalize
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def equal_up_to_phase(a: np.ndarray, b: np.ndarray, atol: float = 1e-6) -> bool:
+    n = a.shape[0]
+    inner = np.trace(a.conj().T @ b)
+    return bool(abs(abs(inner) - n) < atol * n)
+
+
+# ---------------------------------------------------------------------------
+# Unit tests for axis_aligned_rotation
+# ---------------------------------------------------------------------------
+
+class TestAxisAlignedRotation:
+    def _quat_rz(self, theta: float) -> np.ndarray:
+        return quaternion_canonicalize(
+            np.array([math.cos(theta / 2), 0.0, 0.0, math.sin(theta / 2)])
+        )
+
+    def _quat_rx(self, theta: float) -> np.ndarray:
+        return quaternion_canonicalize(
+            np.array([math.cos(theta / 2), math.sin(theta / 2), 0.0, 0.0])
+        )
+
+    def _quat_ry(self, theta: float) -> np.ndarray:
+        return quaternion_canonicalize(
+            np.array([math.cos(theta / 2), 0.0, math.sin(theta / 2), 0.0])
+        )
+
+    @pytest.mark.parametrize("theta", [0.3, 0.7, math.pi / 2, math.pi * 0.9])
+    def test_rz_detected(self, theta: float) -> None:
+        q = self._quat_rz(theta)
+        result = axis_aligned_rotation(q)
+        assert result is not None
+        name, angle = result
+        assert name == "z"
+        assert abs(angle - theta) < 1e-9
+
+    @pytest.mark.parametrize("theta", [0.4, 1.1, math.pi / 3])
+    def test_rx_detected(self, theta: float) -> None:
+        q = self._quat_rx(theta)
+        result = axis_aligned_rotation(q)
+        assert result is not None
+        name, angle = result
+        assert name == "x"
+        assert abs(angle - theta) < 1e-9
+
+    @pytest.mark.parametrize("theta", [0.5, 0.9, math.pi / 4])
+    def test_ry_detected(self, theta: float) -> None:
+        q = self._quat_ry(theta)
+        result = axis_aligned_rotation(q)
+        assert result is not None
+        name, angle = result
+        assert name == "y"
+        assert abs(angle - theta) < 1e-9
+
+    def test_generic_rotation_returns_none(self) -> None:
+        # A rotation about an off-axis direction must return None.
+        q = quaternion_canonicalize(
+            np.array([math.cos(0.3), math.sin(0.3) / math.sqrt(2), math.sin(0.3) / math.sqrt(2), 0.0])
+        )
+        assert axis_aligned_rotation(q) is None
+
+    def test_near_identity_returns_none(self) -> None:
+        q = quaternion_canonicalize(np.array([1.0 - 1e-9, 0.0, 0.0, 1e-9]))
+        result = axis_aligned_rotation(q)
+        # Near-identity: theta ≈ 0, should return None.
+        assert result is None
+
+    def test_custom_atol(self) -> None:
+        # A slightly off-z axis that passes with a looser tolerance.
+        # axis[x] ≈ epsilon / sin(0.5) ≈ epsilon / 0.479; choose epsilon so that
+        # axis[x] < 1e-3 (loose) but > 1e-6 (strict).
+        epsilon = 1e-4  # axis[x] ≈ 2e-4 — inside 1e-3, outside 1e-6
+        q = quaternion_canonicalize(
+            np.array([math.cos(0.5), epsilon, 0.0, math.sin(0.5)])
+        )
+        # Strict tolerance: not aligned.
+        assert axis_aligned_rotation(q, atol=1e-6) is None
+        # Loose tolerance: aligned.
+        result = axis_aligned_rotation(q, atol=1e-3)
+        assert result is not None
+        assert result[0] == "z"
+
+
+# ---------------------------------------------------------------------------
+# Integration: optimize() emits named gates for axis-aligned runs
+# ---------------------------------------------------------------------------
+
+class TestAxisAwareCompression:
+    @pytest.mark.parametrize("a,b", [(0.3, 0.5), (0.7, 1.1), (math.pi / 4, math.pi / 4)])
+    def test_rz_run_emits_rz(self, a: float, b: float) -> None:
+        """Two consecutive Rz gates fuse into a single rz gate."""
+        qc = QuantumCircuit(1)
+        qc.rz(a, 0)
+        qc.rz(b, 0)
+
+        opt = optimize(qc)
+
+        gate_names = [instr.operation.name for instr in opt.data]
+        assert gate_names == ["rz"], f"Expected ['rz'] but got {gate_names}"
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt).data)
+
+    @pytest.mark.parametrize("a,b", [(0.3, 0.5), (0.6, 0.9)])
+    def test_rx_run_emits_rx(self, a: float, b: float) -> None:
+        """Two consecutive Rx gates fuse into a single rx gate."""
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        qc.rx(b, 0)
+
+        opt = optimize(qc)
+
+        gate_names = [instr.operation.name for instr in opt.data]
+        assert gate_names == ["rx"], f"Expected ['rx'] but got {gate_names}"
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt).data)
+
+    @pytest.mark.parametrize("a,b", [(0.4, 0.6), (1.0, 0.5)])
+    def test_ry_run_emits_ry(self, a: float, b: float) -> None:
+        """Two consecutive Ry gates fuse into a single ry gate."""
+        qc = QuantumCircuit(1)
+        qc.ry(a, 0)
+        qc.ry(b, 0)
+
+        opt = optimize(qc)
+
+        gate_names = [instr.operation.name for instr in opt.data]
+        assert gate_names == ["ry"], f"Expected ['ry'] but got {gate_names}"
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt).data)
+
+    def test_rz_three_gates_emits_rz(self) -> None:
+        """Three consecutive Rz gates fuse into a single rz gate."""
+        qc = QuantumCircuit(1)
+        qc.rz(0.3, 0)
+        qc.rz(0.5, 0)
+        qc.rz(0.2, 0)
+
+        opt = optimize(qc)
+
+        gate_names = [instr.operation.name for instr in opt.data]
+        assert gate_names == ["rz"]
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt).data)
+
+    def test_rz_angle_matches_sum(self) -> None:
+        """The fused rz angle equals the sum of the individual angles."""
+        a, b = 0.4, 0.6
+        qc = QuantumCircuit(1)
+        qc.rz(a, 0)
+        qc.rz(b, 0)
+
+        opt = optimize(qc)
+
+        gate_names = [instr.operation.name for instr in opt.data]
+        assert gate_names == ["rz"]
+        emitted_angle = float(opt.data[0].operation.params[0])
+        # Angle should equal a + b (modulo 2π wrapping inside arccos, but close).
+        assert abs(emitted_angle - (a + b)) < 1e-6
+
+    def test_mixed_axes_still_emits_u(self) -> None:
+        """A run of mixed-axis gates (non-axis-aligned result) still emits U."""
+        qc = QuantumCircuit(1)
+        qc.rx(0.5, 0)
+        qc.ry(0.3, 0)
+        qc.rz(0.2, 0)
+
+        opt = optimize(qc)
+
+        gate_names = [instr.operation.name for instr in opt.data]
+        assert gate_names == ["u"]
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt).data)
+
+    def test_rz_run_with_cx_boundary(self) -> None:
+        """Rz gates on separate sides of a CX are each compressed independently."""
+        qc = QuantumCircuit(2)
+        qc.rz(0.3, 0)
+        qc.rz(0.5, 0)
+        qc.cx(0, 1)
+        qc.rz(0.4, 1)
+        qc.rz(0.6, 1)
+
+        opt = optimize(qc)
+
+        rz_ops = [instr for instr in opt.data if instr.operation.name == "rz"]
+        cx_ops = [instr for instr in opt.data if instr.operation.name == "cx"]
+        assert len(rz_ops) == 2
+        assert len(cx_ops) == 1
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt).data)
+
+    def test_axis_aware_with_metadata(self) -> None:
+        """OptimizationResult reflects the gate reduction from axis-aware compression."""
+        from rqm_optimize import OptimizationResult
+
+        qc = QuantumCircuit(1)
+        qc.rz(0.3, 0)
+        qc.rz(0.4, 0)
+        qc.rz(0.5, 0)
+
+        result = optimize(qc, return_metadata=True)
+        assert isinstance(result, OptimizationResult)
+        assert result.original_1q_gate_count == 3
+        assert result.optimized_1q_gate_count == 1
+        assert result.fused_runs == 1
+
+    def test_axis_aware_produces_fewer_gates_than_u_for_ibm_basis(self) -> None:
+        """For axis-aligned runs, axis-aware (1 gate) wins over ZSX (≥1 gate)."""
+        qc = QuantumCircuit(1)
+        qc.rz(0.4, 0)
+        qc.rz(0.5, 0)
+
+        # With default basis: should emit 1 rz via axis-aware.
+        opt_default = optimize(qc)
+        names_default = [i.operation.name for i in opt_default.data]
+        assert names_default == ["rz"]
+
+        # With IBM basis: axis-aware still fires first and emits 1 rz.
+        opt_ibm = optimize(qc, native_basis="ibm")
+        names_ibm = [i.operation.name for i in opt_ibm.data]
+        assert names_ibm == ["rz"]
+
+        # Both must be equivalent to the original.
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt_default).data)
+        assert equal_up_to_phase(Operator(qc).data, Operator(opt_ibm).data)
diff --git a/tests/test_geometry_quaternion.py b/tests/test_geometry_quaternion.py
new file mode 100644
index 0000000..492b123
--- /dev/null
+++ b/tests/test_geometry_quaternion.py
@@ -0,0 +1,337 @@
+"""tests/test_geometry_quaternion.py — Tests for quaternion math in geometry.py.
+
+Verifies the exact algebraic isomorphism between unit quaternions and SU(2)
+matrices, the quaternion product rule for gate composition, canonicalization,
+and axis-angle extraction.
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+
+from rqm_optimize.geometry import (
+    matrices_close,
+    quaternion_canonicalize,
+    quaternion_multiply,
+    quaternion_to_axis_angle,
+    quaternion_to_su2,
+    su2_to_quaternion,
+)
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def _rx(theta: float) -> np.ndarray:
+    """Standard Rx gate matrix."""
+    c, s = math.cos(theta / 2), math.sin(theta / 2)
+    return np.array([[c, -1j * s], [-1j * s, c]], dtype=complex)
+
+
+def _ry(theta: float) -> np.ndarray:
+    """Standard Ry gate matrix."""
+    c, s = math.cos(theta / 2), math.sin(theta / 2)
+    return np.array([[c, -s], [s, c]], dtype=complex)
+
+
+def _rz(theta: float) -> np.ndarray:
+    """Standard Rz gate matrix."""
+    c, s = math.cos(theta / 2), math.sin(theta / 2)
+    return np.array([[complex(c, -s), 0], [0, complex(c, s)]], dtype=complex)
+
+
+def _h() -> np.ndarray:
+    """Hadamard gate matrix."""
+    return np.array([[1, 1], [1, -1]], dtype=complex) / math.sqrt(2)
+
+
+def _qnorm(q: np.ndarray) -> float:
+    return float(np.linalg.norm(q))
+
+
+# ---------------------------------------------------------------------------
+# su2_to_quaternion / quaternion_to_su2 round-trip
+# ---------------------------------------------------------------------------
+
+def test_identity_roundtrip() -> None:
+    """Identity matrix maps to identity quaternion and back."""
+    eye = np.eye(2, dtype=complex)
+    q = su2_to_quaternion(eye)
+    assert np.allclose(q, [1.0, 0.0, 0.0, 0.0], atol=1e-10)
+    recovered = quaternion_to_su2(q)
+    assert np.allclose(recovered, eye, atol=1e-10)
+
+
+@pytest.mark.parametrize("theta", [0.0, math.pi / 4, math.pi / 2, math.pi, 2 * math.pi])
+def test_rx_quaternion_form(theta: float) -> None:
+    """Rx(theta) maps to quaternion (cos(theta/2), sin(theta/2), 0, 0)."""
+    mat = _rx(theta)
+    q = su2_to_quaternion(mat)
+    # Must be unit quaternion.
+    assert abs(_qnorm(q) - 1.0) < 1e-10
+    # Either q or -q should match the expected form.
+    expected = np.array([math.cos(theta / 2), math.sin(theta / 2), 0.0, 0.0])
+    # After canonicalization (w >= 0) both q and expected share the same sign.
+    expected = quaternion_canonicalize(expected)
+    assert np.allclose(q, expected, atol=1e-10), f"theta={theta}: q={q} expected={expected}"
+
+
+@pytest.mark.parametrize("theta", [0.0, math.pi / 3, math.pi / 2, math.pi])
+def test_ry_quaternion_form(theta: float) -> None:
+    """Ry(theta) maps to quaternion (cos(theta/2), 0, sin(theta/2), 0)."""
+    mat = _ry(theta)
+    q = su2_to_quaternion(mat)
+    expected = quaternion_canonicalize(
+        np.array([math.cos(theta / 2), 0.0, math.sin(theta / 2), 0.0])
+    )
+    assert np.allclose(q, expected, atol=1e-10)
+
+
+@pytest.mark.parametrize("theta", [0.0, math.pi / 6, math.pi / 2, math.pi])
+def test_rz_quaternion_form(theta: float) -> None:
+    """Rz(theta) maps to quaternion (cos(theta/2), 0, 0, sin(theta/2))."""
+    mat = _rz(theta)
+    q = su2_to_quaternion(mat)
+    expected = quaternion_canonicalize(
+        np.array([math.cos(theta / 2), 0.0, 0.0, math.sin(theta / 2)])
+    )
+    assert np.allclose(q, expected, atol=1e-10)
+
+
+def test_hadamard_quaternion_form() -> None:
+    """H is a pi-rotation about the (x+z)/sqrt(2) axis.
+
+    Expected quaternion: (0, 1/sqrt(2), 0, 1/sqrt(2)) up to sign/canonicalization.
+    Since H has det=-1 (not SU(2)), global phase must be removed first.
+    """
+    from rqm_optimize.geometry import remove_global_phase
+
+    mat = _h()
+    mat_su2 = remove_global_phase(mat)
+    q = su2_to_quaternion(mat_su2)
+    # Hadamard: theta = pi, axis = (1/sqrt2, 0, 1/sqrt2).
+    # q = (cos(pi/2), sin(pi/2)/sqrt2, 0, sin(pi/2)/sqrt2) = (0, 1/sqrt2, 0, 1/sqrt2)
+    # w=0 so sign ambiguity exists; verify via matrix reconstruction.
+    recovered = quaternion_to_su2(q)
+    assert matrices_close(mat_su2, recovered, atol=1e-10)
+
+
+def test_roundtrip_arbitrary_su2() -> None:
+    """A random SU(2) matrix round-trips through su2_to_quaternion/quaternion_to_su2."""
+    rng = np.random.default_rng(42)
+    for _ in range(20):
+        # Generate a random unitary with det=+1.
+        theta = rng.uniform(0, 2 * math.pi)
+        phi = rng.uniform(0, 2 * math.pi)
+        lam = rng.uniform(0, 2 * math.pi)
+        # Use standard U-gate form.
+        c, s = math.cos(theta / 2), math.sin(theta / 2)
+        mat = np.array(
+            [
+                [c, -np.exp(1j * lam) * s],
+                [np.exp(1j * phi) * s, np.exp(1j * (phi + lam)) * c],
+            ],
+            dtype=complex,
+        )
+        # Normalize to SU(2) (det → +1).
+        det = np.linalg.det(mat)
+        mat /= np.sqrt(det)
+
+        q = su2_to_quaternion(mat)
+        recovered = quaternion_to_su2(q)
+        # Must match up to global phase.
+        assert matrices_close(mat, recovered, atol=1e-9), "Round-trip failed for random SU(2)"
+
+
+# ---------------------------------------------------------------------------
+# quaternion_canonicalize
+# ---------------------------------------------------------------------------
+
+def test_canonicalize_positive_w() -> None:
+    """A quaternion with w>0 is returned unchanged (after normalization)."""
+    q = np.array([0.6, 0.2, 0.5, 0.3])
+    q = q / np.linalg.norm(q)
+    qc = quaternion_canonicalize(q)
+    assert qc[0] >= 0.0
+    assert np.allclose(qc, q)
+
+
+def test_canonicalize_negative_w_flipped() -> None:
+    """A quaternion with w<0 is negated so that w>=0."""
+    q = np.array([-0.6, 0.2, 0.5, 0.3])
+    q = q / np.linalg.norm(q)
+    qc = quaternion_canonicalize(q)
+    assert qc[0] >= 0.0
+    assert np.allclose(qc, -q / np.linalg.norm(q))
+
+
+def test_canonicalize_normalizes() -> None:
+    """quaternion_canonicalize normalizes an unnormalized quaternion."""
+    q = np.array([3.0, 0.0, 0.0, 0.0])
+    qc = quaternion_canonicalize(q)
+    assert abs(_qnorm(qc) - 1.0) < 1e-10
+    assert np.allclose(qc, [1.0, 0.0, 0.0, 0.0])
+
+
+# ---------------------------------------------------------------------------
+# quaternion_multiply — gate composition
+# ---------------------------------------------------------------------------
+
+def test_multiply_identity_left() -> None:
+    """Identity * q == q."""
+    q = su2_to_quaternion(_rx(0.5))
+    identity = np.array([1.0, 0.0, 0.0, 0.0])
+    result = quaternion_multiply(identity, q)
+    assert np.allclose(result, q, atol=1e-12)
+
+
+def test_multiply_identity_right() -> None:
+    """q * identity == q."""
+    q = su2_to_quaternion(_ry(0.7))
+    identity = np.array([1.0, 0.0, 0.0, 0.0])
+    result = quaternion_multiply(q, identity)
+    assert np.allclose(result, q, atol=1e-12)
+
+
+def test_multiply_matches_matrix_composition() -> None:
+    """quaternion_multiply(q2, q1) produces the same SU(2) as U2 @ U1."""
+    from rqm_optimize.geometry import remove_global_phase
+
+    angles = [(0.5, 0.3), (math.pi / 3, math.pi / 7), (1.2, 0.9)]
+    for theta1, theta2 in angles:
+        u1 = _rx(theta1)
+        u2 = _ry(theta2)
+        q1 = su2_to_quaternion(remove_global_phase(u1))
+        q2 = su2_to_quaternion(remove_global_phase(u2))
+        q_prod = quaternion_canonicalize(quaternion_multiply(q2, q1))
+        u_prod = quaternion_to_su2(q_prod)
+        u_direct = remove_global_phase(u2 @ u1)
+        assert matrices_close(u_prod, u_direct, atol=1e-9), (
+            f"Matrix mismatch for theta1={theta1}, theta2={theta2}"
+        )
+
+
+def test_multiply_three_gates_matches_matrix() -> None:
+    """Three-gate quaternion chain matches three-matrix product."""
+    from rqm_optimize.geometry import remove_global_phase
+
+    u1, u2, u3 = _rx(0.4), _ry(0.6), _rz(0.2)
+    q1 = su2_to_quaternion(remove_global_phase(u1))
+    q2 = su2_to_quaternion(remove_global_phase(u2))
+    q3 = su2_to_quaternion(remove_global_phase(u3))
+
+    q_total = quaternion_canonicalize(quaternion_multiply(q3, quaternion_multiply(q2, q1)))
+    u_quat = quaternion_to_su2(q_total)
+    u_matrix = remove_global_phase(u3 @ u2 @ u1)
+    assert matrices_close(u_quat, u_matrix, atol=1e-9)
+
+
+def test_multiply_noncommutativity() -> None:
+    """Rx then Ry is different from Ry then Rx (non-commutative)."""
+    from rqm_optimize.geometry import remove_global_phase
+
+    theta = math.pi / 4
+    q_x = su2_to_quaternion(_rx(theta))
+    q_y = su2_to_quaternion(_ry(theta))
+
+    q_xy = quaternion_multiply(q_y, q_x)  # Rx first, Ry second
+    q_yx = quaternion_multiply(q_x, q_y)  # Ry first, Rx second
+
+    # Products must differ.
+    assert not np.allclose(q_xy, q_yx, atol=1e-6)
+    # But their matrices must individually match the direct matrix products.
+    u_xy = quaternion_to_su2(quaternion_canonicalize(q_xy))
+    u_yx = quaternion_to_su2(quaternion_canonicalize(q_yx))
+    assert matrices_close(u_xy, remove_global_phase(_ry(theta) @ _rx(theta)), atol=1e-9)
+    assert matrices_close(u_yx, remove_global_phase(_rx(theta) @ _ry(theta)), atol=1e-9)
+
+
+def test_multiply_same_axis_adds_angles() -> None:
+    """Rx(a) followed by Rx(b) == Rx(a+b) up to global phase."""
+    a, b = 0.5, 0.7
+    q_a = su2_to_quaternion(_rx(a))
+    q_b = su2_to_quaternion(_rx(b))
+    q_sum = quaternion_canonicalize(quaternion_multiply(q_b, q_a))
+    q_direct = su2_to_quaternion(_rx(a + b))
+    assert np.allclose(q_sum, q_direct, atol=1e-9)
+
+
+# ---------------------------------------------------------------------------
+# quaternion_to_axis_angle
+# ---------------------------------------------------------------------------
+
+def test_axis_angle_identity() -> None:
+    """Identity quaternion gives theta=0 and a default axis."""
+    q = np.array([1.0, 0.0, 0.0, 0.0])
+    axis, theta = quaternion_to_axis_angle(q)
+    assert abs(theta) < 1e-10
+    assert abs(np.linalg.norm(axis) - 1.0) < 1e-10
+
+
+@pytest.mark.parametrize(
+    "theta, expected_axis",
+    [
+        (math.pi / 2, np.array([1.0, 0.0, 0.0])),
+        (math.pi / 3, np.array([1.0, 0.0, 0.0])),
+        (math.pi,     np.array([1.0, 0.0, 0.0])),
+    ],
+)
+def test_axis_angle_rx(theta: float, expected_axis: np.ndarray) -> None:
+    """Rx(theta) has rotation axis x̂ and rotation angle theta."""
+    mat = _rx(theta)
+    q = su2_to_quaternion(mat)
+    axis, angle = quaternion_to_axis_angle(q)
+    assert abs(angle - theta) < 1e-9, f"angle mismatch: {angle} vs {theta}"
+    assert np.allclose(axis, expected_axis, atol=1e-9), f"axis mismatch: {axis}"
+
+
+@pytest.mark.parametrize("theta", [math.pi / 4, math.pi / 2, math.pi])
+def test_axis_angle_ry(theta: float) -> None:
+    """Ry(theta) has rotation axis ŷ and rotation angle theta."""
+    mat = _ry(theta)
+    q = su2_to_quaternion(mat)
+    axis, angle = quaternion_to_axis_angle(q)
+    assert abs(angle - theta) < 1e-9
+    assert np.allclose(axis, [0.0, 1.0, 0.0], atol=1e-9)
+
+
+@pytest.mark.parametrize("theta", [math.pi / 6, math.pi / 2, math.pi])
+def test_axis_angle_rz(theta: float) -> None:
+    """Rz(theta) has rotation axis ẑ and rotation angle theta."""
+    mat = _rz(theta)
+    q = su2_to_quaternion(mat)
+    axis, angle = quaternion_to_axis_angle(q)
+    assert abs(angle - theta) < 1e-9
+    assert np.allclose(axis, [0.0, 0.0, 1.0], atol=1e-9)
+
+
+def test_axis_angle_reconstruction() -> None:
+    """axis_angle → quaternion round-trip via axis-angle-to-quaternion formula."""
+    for theta in [0.3, 0.7, math.pi / 2, math.pi]:
+        for axis_dir in [
+            np.array([1.0, 0.0, 0.0]),
+            np.array([0.0, 1.0, 0.0]),
+            np.array([0.0, 0.0, 1.0]),
+            np.array([1.0, 1.0, 0.0]) / math.sqrt(2),
+        ]:
+            # Build quaternion from axis-angle.
+            q_orig = np.array(
+                [math.cos(theta / 2), *(axis_dir * math.sin(theta / 2))],
+                dtype=float,
+            )
+            q_orig = quaternion_canonicalize(q_orig)
+            # Extract axis-angle back.
+            axis_out, theta_out = quaternion_to_axis_angle(q_orig)
+            # Rebuild quaternion.
+            q_rebuild = np.array(
+                [math.cos(theta_out / 2), *(axis_out * math.sin(theta_out / 2))],
+                dtype=float,
+            )
+            q_rebuild = quaternion_canonicalize(q_rebuild)
+            assert np.allclose(q_orig, q_rebuild, atol=1e-9), (
+                f"Round-trip failed: theta={theta}, axis={axis_dir}"
+            )