The Jacobian: velocities, forces, and singularities

Last lesson you built forward kinematics: a pure function $f(q)$ from six joint angles to a gripper pose. Today you take its derivative — and in your world, the derivative of a forward pass is never just one thing. Forward-mode gives you Jacobian-vector products, reverse-mode gives you vector-Jacobian products, and the conditioning of that Jacobian decides whether optimization converges or explodes. Manipulation is the same picture with a body attached: the JVP is the gripper's velocity, the VJP is the torque your motors must produce, and ill-conditioning is the arm physically unable to move in some direction no matter how fast the motors spin. One matrix, $J(q) \in \mathbb{R}^{6\times 6}$ for your WidowX AI, refit at every control tick, carries all three roles.

The Jacobian, built column by column

Fix a base frame and let $p_e(q)$ be the gripper position and $R_e(q)$ its orientation. Differentiate the pose with respect to time and you get a twist: linear velocity $v_e$ stacked on angular velocity $\omega_e$, six numbers. Because differentiation is linear in $\dot q$, the twist is a linear function of joint velocities at any fixed configuration — that linear map is the manipulator Jacobian.

The columns have a clean physical derivation. Freeze every joint except joint $i$. Everything distal to joint $i$ — links, wrist, gripper, the mug it holds — now rotates as one rigid body about joint $i$'s axis: a unit vector $z_i$ passing through a point $o_i$, both expressed in the base frame. Rigid-body rotation at rate $\dot q_i$ moves any point $p$ with velocity $\dot q_i \, z_i \times (p - o_i)$ and contributes angular velocity $\dot q_i \, z_i$. Set $p = p_e$, divide out $\dot q_i$, and you have column $i$. Because velocities superpose linearly, running all joints at once just sums the columns — which is exactly what the matrix-vector product $J\dot q$ computes.

Everything on the right-hand side is a byproduct of the forward kinematics you wrote last lesson. As you walk the chain multiplying transforms, the world-frame axis $z_i$ and joint origin $o_i$ fall out of the intermediate frames for free — building $J$ costs one FK pass plus six cross products. Units matter and are worth internalizing: the top three rows carry meters (m/s of gripper motion per rad/s of joint motion), the bottom three are dimensionless. Each linear entry is bounded by the distance from that joint's axis to the gripper, so on the 0.769 m-reach WidowX AI no entry exceeds 0.77. Concretely: with the gripper 0.3 m from the base yaw axis, the joint-0 column's linear part has magnitude 0.3 — spin joint 0 at 1 rad/s and the gripper sweeps sideways at 0.3 m/s. Connect that to the latency lesson: a command that is 100 ms stale during that sweep is a 30 mm error, an entire gripper-width.

Geometric vs analytic, and driving the arm with J

One subtlety will bite you the first time you check your Jacobian against autodiff, so meet it now. The geometric Jacobian above outputs angular velocity $\omega$. But $\omega$ is not the time derivative of any three-parameter orientation representation — there is no function $\phi(t)$ of Euler angles with $\dot\phi = \omega$ globally. If instead you differentiate an FK that returns roll-pitch-yaw (or a quaternion), you get the analytic Jacobian $J_a$, whose bottom rows are rates of your chosen coordinates. The two are related by the representation's rate-mapping matrix:

Two practical consequences. First, if you autodiff your FK (JAX or PyTorch over the pose output) and compare against the geometric construction, the position rows will match and the orientation rows will "disagree" — you are comparing $J_a$ to $J_g$ without the $B(\phi)$ conversion. Second, never confuse a representation singularity (RPY pitch near $\pm 90^{\circ}$, fixed by switching to quaternions) with a kinematic singularity (a property of the mechanism, fixed by nothing). Only the geometric Jacobian sees the second kind honestly, which is why everything after this section uses $J_g$.

Now use the map in the forward direction. Resolved-rate motion control (Whitney, 1969) is the oldest trick in manipulation and still the workhorse under your classical pick-and-place baseline: specify the twist you want — say, descend at 5 cm/s while holding orientation — and solve for the joint rates that produce it, once per control tick.

The linearization is only valid locally, but locality is cheap at control rates. At 50 Hz, a 5 cm/s descent moves 1 mm per tick; the error of treating $f$ as linear over a step is second order in the step, so at 1 mm steps on a 0.77 m arm the deviation is micrometers — far below the arm's ~1 mm repeatability. Stretch the tick to 200 ms (a 10 mm step) and the commanded straight line visibly bows near the workspace edge, because $J$ changed along the step and you kept using the stale one. Same failure mode as taking too large a step with a stale gradient — and the fix is the same: smaller steps or re-evaluate more often.

Force duality: τ = JᵀF from virtual work

The same matrix runs backward and becomes a force map. Setup: the gripper presses on a table with wrench $F$ (three forces, three torques, in the base frame) while the arm holds still. What torques $\tau$ must the six joint actuators exert? Use the principle of virtual work. Imagine a virtual displacement $\delta q$ — an infinitesimal joint motion consistent with the kinematics but not actually executed. In static equilibrium with a lossless chain, the net virtual work of all forces vanishes:

This one line explains half of what your hardware does. Gravity intuition: gravity on the full 1.5 kg rated payload is a downward 14.7 N force at the gripper; $J^{\top}$ maps it to joint torques, and the shoulder entry is exactly the horizontal moment arm. At full 0.769 m horizontal extension that is $14.7 \times 0.769 \approx 11.3$ N·m from the payload alone — add the links' own weight and you are pressing toward the 27 N·m effort limit of the shoulder-class joints. That is why payload ratings collapse at extension, and why the arm needs active gravity compensation the moment position control lets go. Feeling contact: the driver streams each joint's measured effort at roughly 500 Hz, so $\hat F = J^{-\top}\tau_{\text{measured}}$ is a free, uncalibrated force sensor. Gear-train friction makes it crude — expect ±20–30% — but it reliably detects a few newtons of table contact within a couple of control cycles, long before anything bends. Your guarded-descend primitive in the pick-and-place baseline is exactly this: resolved-rate downward until the estimated vertical contact force crosses ~2 N.

Foundations: Frames, Twists, and the Wrench

Before deriving the Jacobian, we must rigorously define the coordinate frames and vector quantities involved. The Jacobian is a linear map between two specific vector spaces: joint velocity space and task-space twist. A common source of error is conflating the frame in which the angular velocity is expressed. We distinguish between the spatial twist and the body twist. The spatial twist, denoted $\mathcal{V}_s$, expresses the end-effector's linear velocity $v_s$ and angular velocity $\omega_s$ in the fixed base frame. The body twist, denoted $\mathcal{V}_b$, expresses the same physical motion in the end-effector's own moving frame. The geometric Jacobian $J_g$ maps joint velocities to the spatial twist: $\mathcal{V}_s = J_g \dot{q}$. The analytic Jacobian $J_a$ maps joint velocities to the body twist (or a specific parameterization of orientation): $\mathcal{V}_b = J_a \dot{q}$. This distinction is critical because the angular velocity vector $\omega$ is not the time derivative of Euler angles; it is a physical vector that must be transformed between frames using rotation matrices, not simple differentiation.

The dual quantity to the twist is the wrench, denoted $\mathcal{F}$. A wrench consists of a force vector $f$ and a torque vector $\tau$. In the context of the geometric Jacobian, we define the spatial wrench $\mathcal{F}_s = [f_s; \tau_s]$ where $f_s$ is the force applied at the end-effector origin, expressed in the base frame, and $\tau_s$ is the torque about the end-effector origin, also expressed in the base frame. The principle of virtual work states that for a static system, the work done by joint torques $\tau$ equals the work done by the external wrench $\mathcal{F}_s$ against a virtual displacement. Mathematically, $\tau^T \delta q = \mathcal{F}_s^T \delta x$. Since $\delta x = J_g \delta q$, we substitute to get $\tau^T \delta q = \mathcal{F}_s^T J_g \delta q$. Because this holds for any $\delta q$, we derive the force duality: $\tau = J_g^T \mathcal{F}_s$. Note the sign convention: $\tau$ represents the torque applied by the motors to resist the external load. If the external force is a load, the motor torque opposes it, hence the balance equation.

A critical misconception is that a singular Jacobian implies the arm is 'uncontrollable.' This is false. A singularity means the Jacobian loses rank, so there exists at least one direction in task space that cannot be achieved by any combination of joint velocities. However, the arm remains fully controllable in the subspace spanned by the remaining columns. Furthermore, the arm can often hold significant force in the singular direction with minimal joint torque, because the load is transmitted through the link structure rather than the actuators. This is the duality of mobility and force: directions where mobility is lost are directions where passive force resistance is maximized.

Worked Example: 2-DOF Planar Arm

Consider a 2-DOF planar arm with link lengths $L_1 = 0.5$ m and $L_2 = 0.3$ m. The joint angles are $q_1$ and $q_2$. The end-effector position is $p_e = [L_1 \cos q_1 + L_2 \cos(q_1+q_2), L_1 \sin q_1 + L_2 \sin(q_1+q_2)]^T$. The geometric Jacobian for the position is the derivative of $p_e$ with respect to $q$. At the configuration $q = [0, \pi/2]^T$, the arm is fully extended vertically. We compute the Jacobian, its singular values, and the joint velocities required to achieve a specific Cartesian velocity.

Substituting $q_1 = 0$ and $q_2 = \pi/2$: $\sin(0) = 0, \cos(0) = 1, \sin(\pi/2) = 1, \cos(\pi/2) = 0$. The Jacobian becomes $J = \begin{bmatrix} -0.3 & -0.3 \\ 0.8 & 0 \end{bmatrix}$. The singular values are found via SVD. $J J^T = \begin{bmatrix} 0.18 & -0.24 \\ -0.24 & 0.64 \end{bmatrix}$. The eigenvalues of $J J^T$ are $\lambda_1 \approx 0.7225$ and $\lambda_2 \approx 0.1025$. Thus, $\sigma_1 = \sqrt{0.7225} \approx 0.85$ and $\sigma_2 = \sqrt{0.1025} \approx 0.32$. The condition number $\kappa = \sigma_1 / \sigma_2 \approx 2.66$, indicating a well-conditioned configuration.

Suppose we command a Cartesian velocity $v = [0, 0.1]^T$ m/s (upward). The required joint velocities are $\dot{q} = J^{-1} v$. $J^{-1} = \frac{1}{0.24} \begin{bmatrix} 0 & 0.3 \\ -0.8 & -0.3 \end{bmatrix}$. Thus, $\dot{q} = \frac{1}{0.24} \begin{bmatrix} 0.03 \\ -0.08 \end{bmatrix} = [0.125, -0.333]^T$ rad/s. If we apply damping with $\lambda = 0.1$ m/rad, the damped solution is $\dot{q}^* = J^T (J J^T + \lambda^2 I)^{-1} v$. This reduces the joint velocities slightly, trading tracking accuracy for robustness against noise and singularity proximity.

Singularities: rank loss, manipulability, and the damped fix

$J(q)$ is a different matrix at every configuration, and at some configurations it loses rank. At those points, some twist direction is unreachable: no combination of joint velocities, however large, produces motion along it. The clean lens is the SVD — the same lens you use on attention matrices and low-rank adapters, now with a mechanical meaning for every factor.

Each singular value $\sigma_i$ is a gain: joint effort in the direction $v_i$ (right singular vector, joint space) produces gripper speed $\sigma_i$ along $u_i$ (left singular vector, task space). As the arm approaches a singularity, $\sigma_{\min}\to 0$, the ellipsoid flattens into a pancake, and achieving even modest speed along $u_{\min}$ demands $\dot q = (u_{\min}^{\top}v / \sigma_{\min})\, v_{\min}$ — which diverges. Numbers on your arm: at a comfortable mid-workspace pose, the position block of $J$ has singular values around 0.5, 0.3, and 0.08 m/rad. Near full extension $\sigma_{\min}$ falls below 0.005; requesting 5 cm/s along the weak direction then demands over 10 rad/s of joint speed, past the 6.3 rad/s (360°/s) velocity limit of the WidowX AI's proximal joints. The command gets clipped at the configured joint limits, the direction of motion distorts, and a protective fault can stop the trajectory outright. One units caveat before you compute anything: the full 6×6 $J$ mixes meters and radians, so its raw condition number depends on your unit choice. Either monitor the 3×$n$ position and orientation blocks separately, or scale the angular rows by a characteristic length (0.3 m is reasonable for this arm).

So what does "damp before you divide" actually mean? Stop demanding the exact twist and ask instead for the best compromise between tracking and effort. Pose it as regularized least squares — you have written this exact objective before as weight decay — and solve in closed form:

The filter-factor form is the one to remember, because it hands you a hard guarantee: the gain from task velocity to joint velocity never exceeds $1/(2\lambda)$. With $\lambda = 0.05$ m/rad on the position block, a 5 cm/s request can never command more than 0.5 rad/s from any direction — comfortably inside the joint velocity limits — no matter how singular the pose. The price is bias: along the weak direction the arm moves slower than asked, and with fixed $\lambda$ you pay a small tracking tax even in healthy configurations. The standard refinement is adaptive damping: $\lambda = 0$ while $\sigma_{\min}$ is above a health threshold, ramping smoothly to $\lambda_{\max}$ as it falls below. Hold onto this operator — the next lesson iterates exactly this damped step on a position error to get numeric IK, where it goes by its optimization name, Levenberg–Marquardt.

Code: build it, verify it, monitor it

You would never ship a fused CUDA kernel without checking it against a reference implementation, and the Jacobian deserves the same discipline: one analytic construction, one dumb-but-trustworthy numeric construction, and a test that they agree to floating-point levels at random configurations. Three ways to get $J$, three roles:

import numpy as np

def skew(v):
    return np.array([
        [0.0, -v[2], v[1]],
        [v[2], 0.0, -v[0]],
        [-v[1], v[0], 0.0],
    ])

def rodrigues(axis, angle):
    k = skew(axis)
    return np.eye(3) + np.sin(angle) * k + (1.0 - np.cos(angle)) * (k @ k)

# Illustrative chain approximating a 0.77 m-reach 6-DOF arm (WidowX AI
# class; exact link geometry is not published). Each joint: a unit
# rotation axis and the offset from the previous joint origin, both
# expressed in the parent link frame.
AXES = [
    np.array([0.0, 0.0, 1.0]),   # joint 0: base yaw
    np.array([0.0, 1.0, 0.0]),   # joint 1: shoulder pitch
    np.array([0.0, 1.0, 0.0]),   # joint 2: elbow pitch
    np.array([1.0, 0.0, 0.0]),   # joint 3: wrist roll
    np.array([0.0, 1.0, 0.0]),   # joint 4: wrist pitch
    np.array([1.0, 0.0, 0.0]),   # joint 5: tool roll
]
OFFSETS = [
    np.array([0.0, 0.0, 0.110]),
    np.array([0.0, 0.0, 0.050]),
    np.array([0.06, 0.0, 0.30]),
    np.array([0.17, 0.0, 0.0]),
    np.array([0.115, 0.0, 0.0]),
    np.array([0.075, 0.0, 0.0]),
]
EE_OFFSET = np.array([0.11, 0.0, 0.0])

def fk(q):
    R = np.eye(3)
    p = np.zeros(3)
    axes_w, origins_w = [], []
    for i in range(6):
        p = p + R @ OFFSETS[i]
        axes_w.append(R @ AXES[i])     # joint axis in world, before its own rotation
        origins_w.append(p.copy())
        R = R @ rodrigues(AXES[i], q[i])
    p_ee = p + R @ EE_OFFSET
    return R, p_ee, axes_w, origins_w

def geometric_jacobian(q):
    _, p_ee, axes_w, origins_w = fk(q)
    J = np.zeros((6, 6))
    for i in range(6):
        z = axes_w[i]
        J[:3, i] = np.cross(z, p_ee - origins_w[i])
        J[3:, i] = z
    return J

def vee(m):
    m = 0.5 * (m - m.T)               # project onto skew-symmetric part
    return np.array([m[2, 1], m[0, 2], m[1, 0]])

def numeric_jacobian(q, eps=1e-6):
    R0, _, _, _ = fk(q)
    J = np.zeros((6, 6))
    for i in range(6):
        dq = np.zeros(6)
        dq[i] = eps
        R_plus, p_plus, _, _ = fk(q + dq)
        R_minus, p_minus, _, _ = fk(q - dq)
        J[:3, i] = (p_plus - p_minus) / (2.0 * eps)
        dR = (R_plus - R_minus) / (2.0 * eps)
        J[3:, i] = vee(dR @ R0.T)     # omega-hat = dR R^T for the geometric Jacobian
    return J

rng = np.random.default_rng(0)
worst = 0.0
for _ in range(100):
    q = rng.uniform(-1.5, 1.5, size=6)
    err = np.max(np.abs(geometric_jacobian(q) - numeric_jacobian(q)))
    worst = max(worst, err)
print("worst |J_geo - J_num| over 100 poses:", worst)   # expect ~1e-9

Expected agreement is around $10^{-9}$: central differences at step $10^{-6}$ rad have $\varepsilon^2 \approx 10^{-12}$ truncation error and roughly $10^{-16}/\varepsilon = 10^{-10}$ roundoff, and the position entries are order-1. Two classic failure signatures. If position rows agree and angular rows disagree by a structured factor, you are comparing geometric against analytic — the $B(\phi)$ conversion from earlier. If one column is wrong everywhere, its axis or origin is being read after applying that joint's own rotation instead of before; off-by-one-frame is the FK bug that survives every visual inspection and dies instantly under this test.

The second piece of code is the one that lives inside the control loop: an SVD-based health monitor plus the damped step. A 6×6 SVD costs single-digit microseconds — you spend more on the numpy call overhead — so there is no excuse for not running it every tick and logging the results next to the latency probe from the instrumentation module.

import numpy as np

JOINT_VEL_LIMIT = 6.28   # rad/s: 360 deg/s, the WidowX AI proximal-joint spec
SIGMA_HEALTHY = 0.03     # m/rad threshold on the position block
LAMBDA_MAX = 0.05        # m/rad; caps velocity gain at 1/(2*lambda) = 10

def damped_step(J, v_des, dt):
    """One resolved-rate tick for a Cartesian velocity command.

    J: 6x6 geometric Jacobian at the current q.
    v_des: desired linear velocity (3,), m/s in the base frame.
    Returns the joint step dq (6,) and a health dict for logging.
    """
    Jp = J[:3, :]                          # position block: uniform units (m/rad)
    U, S, Vt = np.linalg.svd(Jp)
    sigma_min = S[-1]

    if sigma_min >= SIGMA_HEALTHY:
        lam = 0.0                          # healthy: exact least-squares solution
    else:
        ratio = sigma_min / SIGMA_HEALTHY  # ramp damping in smoothly (Nakamura-style)
        lam = LAMBDA_MAX * np.sqrt(1.0 - ratio * ratio)

    A = Jp @ Jp.T + (lam * lam) * np.eye(3)
    qdot = Jp.T @ np.linalg.solve(A, v_des)

    peak = np.max(np.abs(qdot))
    scale = 1.0 if peak <= JOINT_VEL_LIMIT else JOINT_VEL_LIMIT / peak
    qdot = qdot * scale                    # uniform scaling preserves direction

    health = {
        "sigma_min": float(sigma_min),
        "kappa": float(S[0] / S[-1]),
        "manipulability": float(np.prod(S)),
        "lambda": float(lam),
        "vel_scale": float(scale),
    }
    return qdot * dt, health

Where this goes next: Forward kinematics from scratch gave you the function; this lesson gave you its derivative and all three of its physical jobs — velocity map, force map, and conditioning oracle. Inverse kinematics: analytic, numeric, and constrained closes the loop: iterate the damped step you derived here on a pose error and you have Levenberg–Marquardt IK, the solver that turns every Cartesian target your planner or your π₀ policy emits into joint commands — with the singularity monitor from this lesson standing guard underneath.