Do you know your Quaternion ?

Quaternions are a tuple of a scalar & a vector \((w, v)\), \[\operatorname{scl} : (w, v) \mapsto w \in \mathbb{R},\quad \operatorname{vec} : (w, v) \mapsto v \in \mathbb{R}^3,\] with the group-product defined by, \[(w_1, v_1) (w_2, v_2) = (w_1 w_2 - v_1 \cdot v_2,\; w_1 v_2 + w_2 v_1 + v_1 \times v_2),\] and inverse by \(q^{-1} = q^* = (w, -v)\).

For the inverse relation to hold, we require the quaternions to be unitary, \[q\, q^* = w^2 + v \cdot v = 1.\]

Quaternions are used widely for parametrizing rotations using their action on \(\mathbb{R}^{3}\), \[\operatorname{rot}_q v \;\equiv\; q\, (0, v)\, q^{*} \;=\; \operatorname{Aut}_q(0, v) \;=\; (0, v'),\] where \(\operatorname{Aut}_g(x) = g x g^{-1}\).

eg. If \(q = \exp({\theta \over 2} \hat{k}) = (\cos({\theta \over 2}), \hat{k} \sin({\theta \over 2}))\), then \(v\) gets rotated by \(\theta\) along the axis defined by the unit vector \(\hat{k}\), \[\operatorname{rot}_q v = \exp({\theta \over 2} \hat{k})\; (0, v)\; \exp(-{\theta \over 2} \hat{k}) = (0, v \cos(\theta) + (\hat{k} \times v) \sin(\theta) + \hat{k} (\hat{k} \cdot v) (1 - \cos(\theta))).\]

The above is the so-called Rodrigues rotation-formula, \[\operatorname{rot}_q v = (I + \sin(\theta) K_{\times} + (1 - \cos(\theta)) K_{\times}^2) v,\] where, \[K_{\times} : x \mapsto \hat{k} \times x.\]

Note: we'll be abusing the notation in the following sections by considering \((0, v)\) and \(v\) to be the same object.

The next building block we need is the effect of small-variations on this action, \[ \delta (\operatorname{rot}_q x) = \delta q\; (0, x)\; q^* + q\; (0, x)\; \delta (q^*) + \operatorname{rot}_q \delta x.\]

Since \(q q^* = 1\), \[\begin{array}{l} \delta q \; q^* + q \; \delta (q^*) = 0, \\ \delta (q^*) = - q^* \; \delta q \; q^*.\end{array}\]

Hence, \[\begin{array}{r l} \delta (\operatorname{rot}_q x) &= \delta q\; (0, x)\; q^* + q\; (0, x)\; \delta (q^*) + \operatorname{rot}_q \delta x\,\\ &= \delta q\; (0, x)\; q^* - q\; (0, x)\; q^* \; \delta q\; q^* + \operatorname{rot}_q \delta x,\\ &= \delta q\; q^*\; q\; (0, x)\; q^* - q\; (0, x)\; q^* \; \delta q\; q^* + \operatorname{rot}_q \delta x,\\ &= \{\delta q\; q^* , \operatorname{rot}_q x\} + \operatorname{rot}_q \delta x\\ \\\end{array}\]

From the quaternion-product, \[\begin{array}{r l} \{a, b\} &= a b - b a ,\\ &= (0, 2 \operatorname{vec}[a] \times \operatorname{vec}[b]). \\\end{array}\]

Ergo, \[\begin{array}{r l} \delta (\operatorname{rot}_q x) &= \{\delta q\; q^* , \operatorname{rot}_q x\} + \operatorname{rot}_q \delta x,\\ &= 2 \operatorname{vec}[\delta q\; q^*] \times \operatorname{rot}_q x + \operatorname{rot}_q \delta x,\\ &= \boxed{\omega \times \operatorname{rot}_q x + \operatorname{rot}_q \delta x,\quad \omega = 2 \operatorname{vec}[\delta q\; q^*]},\\ &= q \{ q^*\; \delta q , x\} q^{*} + \operatorname{rot}_q \delta x,\\ &= \operatorname{rot}_q (2 \operatorname{vec}[q^{*}\; \delta q] \times x) + \operatorname{rot}_q \delta x,\\ &= \boxed{\operatorname{rot}_q (\Omega \times x + \delta x), \quad \Omega = 2 \operatorname{vec}[q^* \; \delta q]}.\\ \end{array}\]

Also, \[\delta q = (0, {\omega \over 2})\; q = q\; (0, {\Omega \over 2})\]

Typically, one parametrizes variations using either \(\omega\) or \(\Omega\) rather than \(\delta q\) in the embedding coordinates, since the latter must respect the embedding's unitary-quaternion constraint on the tangent space.

Suppose we have some function \(f: Q \rightarrow \mathbb{R}\). Variations of \(f\) are of the form,

\[\begin{array}{r l} \delta f &= df \cdot \delta q,\\ &= df \cdot (0, \omega / 2) \; (\epsilon, \xi), \quad \mbox{let } q = (\epsilon, \xi),\\ &= df \cdot {1 \over 2} (- \omega \cdot \xi, \epsilon \omega + \omega \times \xi). \end{array}\]

Suppose \(df = (df_w, df_v)\), \[\begin{array}{r l} \delta f &= (df_w, df_v) \cdot {1 \over 2} (- \omega \cdot \xi, \epsilon \omega + \omega \times \xi),\\ &= {1 \over 2} \left( - df_w \omega \cdot \xi + \epsilon df_v \cdot \omega + df_v \cdot (\omega \times \xi) \right),\\ &= {1 \over 2} \left( (\epsilon df_v - df_w \xi) \cdot \omega + \omega \cdot (\xi \times df_v) \right),\\ &= {1 \over 2} (\epsilon \; df_v - df_w\; \xi + \xi \times df_v) \cdot \omega,\\ &= {1 \over 2} \operatorname{vec}[(df_w, df_v) (\epsilon, - \xi)] \cdot \omega,\\ &= {1 \over 2} \operatorname{vec}[df \; q^*] \cdot \omega. \end{array}\]

For a concrete demonstration, let's look at how to deal with the derivatives of pose-quaternions in MuJoCo. The object-pose shows up in loss-functions ranging from IK to state-estimation to control, so their derivatives show up everywhere as well. There are multiple ways of measuring deviations from a target-pose: one is to minimize the vector part of the difference, \(\operatorname{vec}[q_s(\theta)\, q_t^{-1}]\); another is to minimize the (quaternion) log of the difference, \(\log(q_t\, q_s(\theta)^{-1})\). The latter is very popular in the MuJoCo ecosystem for IK, etc.

The loss-function is typically the $L²$-norm of these features, so that iteration can proceed via the Gauss-Newton method using the feature-Jacobian. The code-fragment below shows how to compute this Jacobian exactly.

The joint-space gradient chains through three layers, \[F \xrightarrow{\;dF\;} q_s \xrightarrow{\;\text{pull-back}\;} \omega \xrightarrow{\;J_\omega\;} \theta,\] with the middle arrow being exactly the identity derived above, \[\partial_\omega F_i = {1 \over 2} \operatorname{vec}[dF_i \; q^*]\] and \(J_\omega := \partial\omega_{\text{spatial}} / \partial\theta\) is MuJoCo's site Jacobian (mj_jacSite), and so \(\partial_\theta F_i = \partial_{\omega} F_i\, J_{\omega}\). We check this against the finite-difference ground truth.

First, let's begin with the group operations,

import jax, jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

def qmul(a, b):
    aw, av = a[0], a[1:]; bw, bv = b[0], b[1:]
    return jnp.concatenate([(aw*bw - av @ bv)[None],
                            aw*bv + bw*av + jnp.cross(av, bv)])

qconj = lambda q: q * jnp.array([1., -1., -1., -1.])
pure  = lambda v: jnp.concatenate([jnp.zeros(1), v])
rot   = lambda R, v: qmul(qmul(R, pure(v)), qconj(R))[1:]      # rot_R(v)

def qexp_pure(v):                                # exp((0, v)) = (cos|v|, v̂ sin|v|)
    n = jnp.linalg.norm(v)
    return jnp.concatenate([jnp.cos(n)[None], v * jnp.sinc(n / jnp.pi)])

def qlog(q):                                                 # ~ mujoco.mju_quat2Vel (dt=1)
    q    = jnp.where(q[0] < 0, -q, q)                        # canonical branch, angle ∈ [0, π]
    w, v = q[0], q[1:]
    n    = jnp.linalg.norm(v)
    return v * jnp.where(n > 1e-12, 2 * jnp.arctan2(n, w) / n, 2. / w)

We use the Unitree G1 humanoid (from Mujoco Menagerie) to build this example, and use the deviation of a site \(s\) on the body (left_foot) from an arbitrary target as our features.

import mujoco
import numpy as np
import jax, jax.numpy as jnp

model = mujoco.MjModel.from_xml_path('../mujoco_menagerie/unitree_g1/g1.xml')
data = mujoco.MjData(model)
site_id = model.site('left_foot').id

mujoco.mj_resetData(model, data)
mujoco.mj_forward(model, data)

def site_quat(site_id=site_id):
    q = np.empty(4);  mujoco.mju_mat2Quat(q, data.site_xmat[site_id])
    return q

rng = np.random.default_rng(0)
q_target = rng.standard_normal(4);  q_target /= np.linalg.norm(q_target)
q_target = jnp.array(q_target)

#######
# vec[q_s q_t^*]
def lstsq_features_0(q_s, q_t):
    return qmul(q_s, qconj(q_t))[1:]

# log(q_t q_s^*)     -- dm_control convention
def lstsq_features_1(q_s, q_t):
    return qlog(qmul(q_t, qconj(q_s)))

# concatenate
def lstsq_features(q_s, q_t):
    return jnp.concatenate([lstsq_features_0(q_s, q_t),
                            lstsq_features_1(q_s, q_t)])

J_feat_q = jax.jacobian(lstsq_features)

# wire features to theta
def lstsq_features_theta(qpos):
    orig = data.qpos.copy()
    data.qpos[:] = qpos;   mujoco.mj_forward(model, data)
    val = lstsq_features(jnp.array(site_quat()), q_target)
    data.qpos[:] = orig;   mujoco.mj_forward(model, data)
    return val

### FD Jacobian
theta0 = data.qpos.copy();  h = 1e-6
feat0 = lstsq_features_theta(theta0)

J_feat_theta_fd = np.zeros((feat0.shape[0], model.nv))
for j in range(model.nv):
    dtheta = np.zeros(model.nv);  dtheta[j] = 1.0
    theta_eps = theta0.copy();    mujoco.mj_integratePos(model, theta_eps, dtheta, h)
    J_feat_theta_fd[:, j] = (lstsq_features_theta(theta_eps) - feat0) / h

#### Analytic Jacobian
q_site  = site_quat()
df = np.asarray(J_feat_q(jnp.array(q_site), q_target))                    # 4-covector

# pull-back:  ∂_ω F = ½ vec(dF q^*)
J_feat_omega = 0.5 * jax.vmap(qmul, in_axes=(0, None))(df, qconj(q_site))[:, 1:]

# spatial angular-velocity Jacobian of the site (MuJoCo uses the right-invariant world-frame \omega)
jacp = np.zeros((3, model.nv))
J_omega_theta = jacr = np.zeros((3, model.nv))
mujoco.mj_jacSite(model, data, jacp, jacr, site_id)

J_feat_theta = J_feat_omega @ J_omega_theta                               # analytic ∂F/∂θ

### 
print(f"‖J_analytic − J_fd‖₂  = {np.linalg.norm(J_feat_theta - J_feat_theta_fd):.2e}")
print(f"‖J_analytic − J_fd‖_∞ = {np.max(np.abs(J_feat_theta - J_feat_theta_fd)):.2e}")

‖J_analytic − J_fd‖₂  = 8.03e-07
‖J_analytic − J_fd‖_∞ = 3.06e-07

Since \(\log(q_t\, q_s^*) \sim -\theta\, \hat{v}\) (where \(q_s = \exp({1 \over 2} \theta\, \hat{v})\)) and \(q_t\) is constant, many libraries in the MuJoCo ecosystem (e.g. dm_control) use \(-J_\omega\) directly as the Jacobian of the log-loss features. This approximation holds when the rotation axis is fixed (e.g. planar systems, where the dynamics lives in an abelian subgroup), but fails in general.

It fails for Unitree as well,

print(f"‖J_analytic − J_approx‖₂  = {np.linalg.norm(J_feat_theta[3:] - (-J_omega_theta)):.2e}\n")

print(f"J_logfeat_omega = \n{J_feat_omega[3:, ]}")

‖J_analytic − J_approx‖₂  = 3.54e+00

J_logfeat_omega = 
[[-0.29193585  0.36133125 -1.31425973]
 [-0.07684831 -0.95215787 -0.38885768]
 [ 1.36085755  0.16295916 -0.28109551]]

The pull-back of the log-loss features to $ω$-space isn't even diagonal!

It's surprising that this simple trick for computing analytical derivatives on quaternions isn't more widely known.

The notation above was very hand-wavy by choice. The usual route involves Lie groups and algebras, which requires some mathematical maturity to follow but ultimately doesn't add much for readers not dealing with general Lie groups.

For a more thorough, scenic treatment, please refer to,

• [1]. Olver, Peter J. Applications of Lie groups to Differential Equations. Vol. 107. Springer Science & Business Media, 1993.
• [2]. Bloch, Anthony M. Nonholonomic Mechanics and Control. Springer New York, NY, 2016.

Alternatively, one can take SICM's path of using rotation-matrices directly,

• [3]. Sussman, Gerald Jay, and Jack Wisdom. Structure and Interpretation of Classical Mechanics. MIT Press, 2015.
• [4]. Sussman, Gerald Jay, and Jack Wisdom. Functional differential geometry. MIT Press, 2013.

• <2026-07-17 शुक्र>: original version used \(\log(q_s\, q_t^*)\); replaced this with what dm_control uses \(\log(q_t\, q_s^*)\).