Quaternions are the tool we use to represent 3D orientation. Here, we build the algebra from scratch — definition, product, conjugate, norm, inverse, unit quaternions.
A quaternion $\mathbf q$ extends the complex numbers with three imaginary units $i,j,k$. Following the Cayley–Dickson construction, taking two complex numbers $A=a+bi$ and $C=c+di$ and forming $\mathbf q = A + Cj$ with $k\triangleq ij$ gives a number with one real and three imaginary parts:
$$\mathbf q = q_0 + q_1 i + q_2 j + q_3 k \in \mathbb H ,\qquad \{q_0,q_1,q_2,q_3\}\in\mathbb R,$$
where $\mathbb H$ is the space of quaternions. We split it into a scalar part $q_w$ and a vector part $\mathbf q_v=(q_x,q_y,q_z)$, and write it as a 4-vector in scalar-first order:
$$\mathbf q = q_w + \mathbf q_v = \begin{bmatrix} q_w \\ \mathbf q_v \end{bmatrix} = \begin{bmatrix} q_w \\ q_x \\ q_y \\ q_z \end{bmatrix}.$$
Convention. We use the Hamilton convention throughout: scalar part first, and $ij=+k$ (right-handed). This matches Eigen/ROS/Ceres and Joan Solà's reference.
The units $i,j,k$ are defined by Hamilton's multiplication rules:
$$i^2 = j^2 = k^2 = ijk = -1,$$
from which all pairwise products follow:
$$ij = k,\quad jk = i,\quad ki = j,\qquad ji = -k,\quad kj = -i,\quad ik = -j.$$
Because $ij=k$ but $ji=-k$, the units anti-commute, and quaternion multiplication is therefore non-commutative: in general $\mathbf q_a \otimes \mathbf q_b \neq \mathbf q_b \otimes \mathbf q_a$ (the symbol $\otimes$ is defined below). This is the algebraic root of the fact that 3D rotations do not commute.
Multiplying out $(q_{a,w}+\mathbf q_{a,v})(q_{b,w}+\mathbf q_{b,v})$ with the unit rules above and collecting scalar and vector parts gives the compact form
$$\mathbf q_a \otimes \mathbf q_b = \begin{bmatrix} q_{a,w}\,q_{b,w} - \mathbf q_{a,v}\!\cdot\!\mathbf q_{b,v} \\[2pt] q_{a,w}\,\mathbf q_{b,v} + q_{b,w}\,\mathbf q_{a,v} + \mathbf q_{a,v}\times\mathbf q_{b,v}\end{bmatrix}.$$
The cross-product term is exactly what makes the product non-commutative (swapping $a,b$ flips its sign). The product is also bilinear, so it can be written as a matrix–vector product in two ways — with the left and right product matrices:
$$\mathbf q_a \otimes \mathbf q_b = [\mathbf q_a]_L\,\mathbf q_b = [\mathbf q_b]_R\,\mathbf q_a ,\qquad [\mathbf q]_{L,R} = q_w\mathbf I + \begin{bmatrix} 0 & -\mathbf q_v^\top \\ \mathbf q_v & \pm[\mathbf q_v]_\times\end{bmatrix}$$
($+$ for $[\mathbf q]_L$, $-$ for $[\mathbf q]_R$; see the digest §B). These matrix forms are handy for the Kalman-filter Jacobians later.
A unit quaternion has $\lVert\mathbf q\rVert = 1$. For these the inverse is simply the conjugate, $\mathbf q^{-1}=\mathbf q^*$, and — because the norm is multiplicative — the product of two unit quaternions is again a unit quaternion. The unit quaternions thus form a group, and it is this group that represents 3D rotations:
$$\mathbf q = \begin{bmatrix}\cos(\phi/2)\\ \mathbf u\,\sin(\phi/2)\end{bmatrix}\quad\text{encodes a rotation of angle }\phi\text{ about axis }\mathbf u.$$
How a unit quaternion rotates a vector (the sandwich product) we be covered in the next chapter.
A few lines confirm the algebra: the product is non-commutative, the norm is multiplicative, the conjugate reverses product order, and $\mathbf q\otimes\mathbf q^{-1}$ is the identity.
import numpy as np
# Hamilton convention, scalar-first q = [q_w, q_x, q_y, q_z] = [q_w, q_v]
def qmul(a, b): # Hamilton product a ⊗ b
aw, av = a[0], a[1:]; bw, bv = b[0], b[1:]
return np.concatenate([[aw*bw - av @ bv],
aw*bv + bw*av + np.cross(av, bv)])
def qconj(q): return np.concatenate([[q[0]], -q[1:]]) # conjugate q*
def qnorm(q): return np.linalg.norm(q) # ||q||
def qinv(q): return qconj(q) / (q @ q) # q^{-1} = q* / ||q||^2
qa = np.array([0.5, 1.0, -2.0, 0.3])
qb = np.array([-1.0, 0.4, 0.7, 2.0])
identity = np.array([1.0, 0.0, 0.0, 0.0])
print("non-commutative : ||qa⊗qb - qb⊗qa|| =", round(np.linalg.norm(qmul(qa, qb) - qmul(qb, qa)), 4))
print("norm multiplies : | ||qa⊗qb|| - ||qa||·||qb|| | =",
abs(qnorm(qmul(qa, qb)) - qnorm(qa) * qnorm(qb)))
print("conjugate reverses : ||(qa⊗qb)* - qb*⊗qa*|| =",
np.linalg.norm(qconj(qmul(qa, qb)) - qmul(qconj(qb), qconj(qa))))
print("inverse : ||qa⊗qa^-1 - 1|| =", np.linalg.norm(qmul(qa, qinv(qa)) - identity))