Inverse Kinematics for Robot Arms: Jacobians, Damped Least Squares, Singularities and Analytical Solvers (2026)

Ask a robot arm to put its gripper at a pose and you have asked a question with no guaranteed answer. There may be eight joint configurations that reach the pose, there may be one, there may be none, and there may be a continuum of them if the arm has a spare joint. The computation that turns “where the hand should be” into “what the joints should do” is inverse kinematics, and it sits underneath every pick-and-place cell, every teleoperation rig and every motion planner you will ever run on a manipulator.

It matters more in 2026 than it did a few years ago because the loop above it has changed. Learned policies and vision-language-action models emit end-effector targets at tens of hertz, and something deterministic has to convert them into joint commands at one kilohertz without jerking the arm through a singularity. Most production failures I have seen traced back to the IK layer, not the planner: a branch flip, a quaternion sign, an undamped pseudo-inverse.

You will leave with the math, working code and a decision rule. The math covers forward kinematics, the manipulator Jacobian, and the pseudo-inverse and damped least squares updates. The code is a runnable NumPy solver for a planar arm. The decision rule tells you when to choose closed-form solvers and when to choose numerical ones.

What this covers: forward kinematics with DH parameters and product of exponentials, the Jacobian, closed-form solutions via spherical-wrist decoupling and IKFast, numerical methods from Newton-Raphson to Levenberg-Marquardt, singularities and manipulability, null-space redundancy resolution, joint limits, orientation error, the solver ecosystem, 1 kHz considerations, and failure modes.

Context and Background

A serial manipulator is a chain of rigid links connected by one-degree-of-freedom joints, almost always revolute in industrial arms. Its configuration is a vector of joint angles q with n entries. Its end-effector pose lives in SE(3), the six-dimensional space of positions and orientations. Forward kinematics is the map f from q to pose; it is a straightforward product of transforms and always has exactly one answer. Inverse kinematics inverts f, and inverting a nonlinear, many-to-one, partly non-surjective map is where all the difficulty comes from.

For most of robotics history the field split into two camps. Industrial vendors shipped closed-form solvers for the common 6-axis geometry, because a controller must solve IK thousands of times per second and cannot afford an iteration that might not converge. Research and service robotics, with 7-DOF arms, mobile manipulators and humanoids, moved to numerical solvers that treat IK as an optimisation problem. The tooling reflects the split: MoveIt’s default plugin has long been the KDL numerical solver, with IKFast generating closed-form C++ for supported arms. If you are planning your stack around a specific MoveIt generation, our MoveIt Pro 10.0 vs 9.4 migration guide covers how the kinematics plugin interface fits into the newer releases.

Three shifts define the 2026 picture. First, GPU motion generation has made batched IK cheap: planners such as the one we examine in cuRobo versus OMPL solve thousands of seeds in parallel rather than one seed carefully. Second, differential IK formulated as a quadratic program, as in the open-source Pink library built on Pinocchio, has become the default for teleoperation and humanoid whole-body control. Third, the robot model itself is increasingly a shared artifact between simulator, planner and digital twin, which makes the description format a real dependency; see our comparison of URDF, SDF, MJCF and OpenUSD.

The classic reference for the numerical side is Samuel Buss’s survey, “Introduction to Inverse Kinematics with Jacobian Transpose, Pseudoinverse and Damped Least Squares methods” (University of California San Diego), which I lean on for the damped least squares formulas and step clamping below. For the geometry, Lynch and Park’s Modern Robotics is the standard modern text for the product of exponentials formulation.

Before the machinery, fix the vocabulary. The joint space is the space of q. The task space is the space of end-effector poses (or just positions, if orientation does not matter). The workspace is the set of poses the arm can reach at all. The dexterous workspace is the subset reachable in every orientation. A pose outside the workspace has no IK solution, and a robust solver must say so rather than return the least-bad garbage.

The Reference Model: From Joints to Pose and Back

Direct answer: Forward kinematics multiplies one homogeneous transform per joint to get the end-effector pose from joint angles, always uniquely. Inverse kinematics reverses this and may return zero, one, several or infinitely many solutions. Practical solvers either exploit special geometry for closed-form answers or iterate with the Jacobian, the matrix that maps joint velocities to end-effector velocity.

Figure 1: The inverse kinematics pipeline. Forward kinematics and the Jacobian come from the robot model; the solver is either closed-form or iterative, followed by branch selection and joint-limit checks.

Figure 1 shows the data flow. The robot description supplies the kinematic parameters. From them you build the forward map and its Jacobian. A solver consumes the target pose and a seed configuration, and returns candidate joint vectors. A final stage picks a branch, checks joint limits and rejects solutions that collide. Every serious IK bug lives in one of those boxes, and the rest of this article walks them in order.

Forward kinematics: DH parameters and product of exponentials

The Denavit-Hartenberg (DH) convention describes each joint with four numbers: link length a, link twist alpha, link offset d and joint angle theta. The transform from frame i-1 to frame i is the product of a rotation about z by theta, a translation along z by d, a translation along x by a and a rotation about x by alpha. Chain n of them and you have the base-to-flange transform. DH is compact and ubiquitous in vendor datasheets, but it has two well-known annoyances. The “standard” and “modified” (Craig) variants place frames differently and are routinely confused, and DH is singular for parallel adjacent joint axes, where the common normal is not unique.

The product of exponentials (POE) formulation avoids both. You pick a single fixed base frame, describe each joint by a screw axis S_i (a unit rotation axis plus a point on it), and write the end-effector pose as the product of matrix exponentials of the screw axes times the home configuration M:

T(q) = exp([S_1] q_1) exp([S_2] q_2) … exp([S_n] q_n) M

There is no per-link frame bookkeeping and no singular parameterisation, which is why URDF-style descriptions, which already give each joint an origin and an axis, map onto POE almost directly. For a revolute joint the exponential is just Rodrigues’ rotation formula plus the translation induced by the axis not passing through the origin.

A worked two-link example

Take a planar 2-link arm with link lengths L1 = 0.5 m and L2 = 0.4 m (illustrative values). Forward kinematics is
x = L1 cos(q1) + L2 cos(q1 + q2)
y = L1 sin(q1) + L2 sin(q1 + q2)

Now ask for the point (0.6, 0.3). The squared distance is 0.45. The law of cosines gives cos(q2) = (0.45 – 0.25 – 0.16) / (2 x 0.5 x 0.4) = 0.04 / 0.4 = 0.1, so q2 = plus or minus 84.26 degrees. For each sign q1 follows from the angle to the target minus the angle the second link makes within the first link’s frame. For q2 = +84.26 degrees, q1 = -9.83 degrees; for q2 = -84.26 degrees, q1 = +62.96 degrees. Both reach exactly the same point. This is the elbow-up and elbow-down pair, the simplest case of IK having multiple solutions.

If the target were at distance greater than L1 + L2 = 0.9 m, cos(q2) would exceed 1 and there would be no real solution; if exactly 0.9 m, q2 = 0 and the arm is fully stretched, the first singularity we will meet. If the target is closer than |L1 – L2| = 0.1 m, cos(q2) falls below -1 and the point is inside the unreachable hole. Even this trivial arm has a workspace that is an annulus, and every edge of the annulus is a place where solvers misbehave.

The manipulator Jacobian

Differentiate forward kinematics with respect to time and you get v = J(q) q_dot, where v stacks linear and angular velocity of the end effector and J is 6 x n for a spatial arm. Each column is the contribution of one joint: for a revolute joint i, the angular part is its axis z_i expressed in the base frame and the linear part is z_i cross (p_e – p_i), where p_e is the end-effector position and p_i the joint origin. This geometric Jacobian is cheap to build while you are already computing forward kinematics.

For the two-link arm,
J = [[-L1 sin q1 – L2 sin(q1+q2), -L2 sin(q1+q2)], [L1 cos q1 + L2 cos(q1+q2), L2 cos(q1+q2)]]
and its determinant simplifies to L1 L2 sin(q2). At the worked solution above that is 0.5 x 0.4 x sin(84.26 degrees) = 0.199 (illustrative units: square metres), close to the maximum 0.2. At q2 = 0 or 180 degrees it is zero, and the matrix cannot be inverted. The same fact is visible in the singular values. At the elbow-up solution I computed them numerically as roughly 0.73 and 0.27, so the arm moves comfortably in every direction; as q2 approaches zero the smaller singular value goes to zero and one direction of motion disappears.

The Jacobian is also the bridge between numerical IK and dynamics: the same matrix maps joint torques to end-effector wrench through its transpose, which is the basis of operational-space control. You build it once and use it everywhere.

Why six degrees of freedom is the magic number

A rigid body in space has six degrees of freedom, so an arm needs at least six independent joints to reach arbitrary poses in its workspace. With exactly six, the Jacobian is square, and in general the pose equations have a finite set of isolated solutions. With seven the system is under-determined: a one-parameter family of joint vectors reaches the same pose, which is redundancy, and which we handle with null-space projection later. With five or fewer, the arm cannot reach arbitrary orientations and the solver must decide which pose components to sacrifice.

The finite-solution property for six joints is a classical result: a general 6R serial chain has at most sixteen inverse kinematic solutions (Raghavan and Roth, 1990), but you will never meet that number on an industrial arm, because manufacturers choose geometry that makes the equations factor.

Closed-Form Solutions: Wrist Decoupling and IKFast

Closed-form (analytical) IK gives every solution in a fixed number of arithmetic operations, with no seed, no iteration count and no convergence tolerance. For the right geometry it is the best answer by every metric that matters in production: determinism, speed and completeness.

Spherical-wrist decoupling (Pieper)

Pieper’s observation, from his 1968 thesis, is that if three consecutive joint axes of a 6-DOF arm intersect at a point, or are parallel, the position and orientation problems separate. Almost every industrial six-axis arm has a spherical wrist, with joints 4, 5 and 6 axes meeting at the wrist centre. The procedure runs in three steps.

Compute the wrist centre: subtract the tool offset along the flange approach axis from the target position. This point depends only on joints 1 to 3.
Solve joints 1 to 3 from the wrist-centre position. Joint 1 is an atan2 of the wrist centre’s horizontal coordinates (with a second solution 180 degrees away for the “backward” posture). Joints 2 and 3 reduce to a planar two-link problem like the one above, with its own elbow-up and elbow-down pair.
Compute the rotation produced by joints 1 to 3, call it R03, and solve for the residual wrist rotation R36 = R03 transposed times R_target. Read joints 4, 5 and 6 from R36 as Euler angles (a ZYZ decomposition for the common wrist layout), again with two solutions, “flip” and “no flip”.

Two choices at joint 1, two at the elbow, two at the wrist yield 2 x 2 x 2 = 8 solutions, the number you will see quoted for typical industrial 6R arms. Many are discarded by joint limits, and in practice a controller usually has two or three feasible ones.

The pitfalls of a hand-derived solver

Doing this by hand is a rite of passage and a standing source of bugs. The ones I see most often are:

Using acos instead of atan2 for the elbow, which loses sign information and is numerically poor near full extension. Use atan2(sqrt(1 – c^2), c) with c clamped to [-1, 1].
Failing to clamp c at all, so a target at the exact workspace boundary produces NaN due to floating-point error making c = 1.0000000002.
Not handling the wrist singularity, where joint 5 is zero and joints 4 and 6 become collinear. Only their sum is defined, so the solver must fix one of them (often to its current value) and compute the other.
Mixing up modified and standard DH when copying parameters from a datasheet.
Forgetting that offsets such as a shoulder offset or an elbow offset in the arm geometry break the textbook formulae, and need an extra atan2 term.

IKFast

IKFast, part of the OpenRAVE project, automates the derivation. According to the MoveIt documentation, it “automatically analyses any complex kinematic chain for common patterns that allow for an analytic solution” and generates C++ code that MoveIt wraps as a kinematics plugin to replace the default KDL solver. The same documentation describes the generated solutions as stable and found “in a few microseconds on recent processors”. Those are the project’s own claims; measure on your hardware before depending on them.

The documentation also lists the catches. It supports 6-DOF and 7-DOF arms and not more than 7. Generation is per planning group, with explicit base and end-effector links. Installing OpenRAVE is described as tricky, its public documentation is no longer maintained, and the recommended workflow is a Docker image used only to generate the solver source. Floating-point trouble in the URDF-to-Collada conversion can stop generation from finding solutions. For 7-DOF arms IKFast does not give a pure closed form: it solves for the pose given one joint fixed, so you sweep a “free parameter”, typically by discretising it, which is why it is fast but not exhaustive. Verify details against the current MoveIt version you use, since the tutorial’s own tested environment is dated.

Closed-form is the right default when the arm has a spherical wrist or equivalent decoupling, when you need all solutions to choose by criteria such as continuity, and when you must meet hard real-time budgets. It is the wrong choice for redundant arms, mobile manipulators or tasks that mix constraints (keep the tool level, avoid an obstacle, bias toward a posture). For those, you iterate.

Numerical Inverse Kinematics: From Newton-Raphson to Damped Least Squares

Direct answer: Numerical IK treats the target as an error e = x_target – f(q) and updates joints by delta_q = J-sharp e, repeating until the error is small. The choice of J-sharp (transpose, pseudo-inverse or damped inverse) determines the behaviour near singularities. Damped least squares is the standard choice because it stays bounded when the Jacobian loses rank.

Figure 2: One iteration of damped least squares inverse kinematics: compute error, clamp it, build the Jacobian, solve with damping, limit the step and project onto joint limits.

Newton-Raphson and the pseudo-inverse

Linearise forward kinematics around the current configuration: f(q + dq) is approximately f(q) + J dq. Setting this equal to the target gives J dq = e. For a square, full-rank J, dq = J inverse times e, which is Newton-Raphson on the root-finding problem f(q) = target, quadratically convergent near the answer. For a non-square J, the Moore-Penrose pseudo-inverse J-dagger gives the minimum-norm least-squares dq, which is exactly what you want for a redundant arm: of all joint motions that produce the requested end-effector motion, the smallest one.

The trouble is in the singular value decomposition. Write J = U Sigma V transposed with singular values sigma_i. The pseudo-inverse scales each direction by 1/sigma_i. As an arm approaches a singularity, one sigma_i heads to zero, and the pseudo-inverse response to any error component in that direction heads to infinity. Buss puts it plainly: near singularities the pseudo-inverse is unstable, with large joint changes for small target movements, and since exact singularities are almost never hit in floating point, you cannot just test for zero.

Jacobian transpose

The simplest alternative is dq = alpha J transpose e. It needs no inversion, it cannot blow up, and for a small enough scalar alpha it reduces the error, because the inner product of e with J J transpose e is never negative. It is also slow: convergence is linear at best, and the right alpha changes as the arm moves. It earns its place in very cheap controllers and as the gradient direction inside larger optimisers, but I rarely recommend it as a primary solver.

Damped least squares

Damped least squares (DLS), also called the Levenberg-Marquardt method in the nonlinear least-squares literature and the singularity-robust inverse in control papers from the 1980s, fixes the pseudo-inverse by penalising large joint motions. It minimises

|| J dq – e ||^2 + lambda^2 || dq ||^2

for a damping constant lambda. The solution, in Buss’s equations, is
dq = J transpose (J J transpose + lambda^2 I)^-1 e
or equivalently dq = (J transpose J + lambda^2 I)^-1 J transpose e.

The first form inverts an m x m matrix where m is the task dimension, 3 for position-only IK and 6 for full pose, so it is preferable when m is smaller than the number of joints. In singular-value terms, each direction is scaled by sigma_i / (sigma_i^2 + lambda^2) instead of 1/sigma_i. For sigma_i much larger than lambda this is almost the pseudo-inverse, and as sigma_i goes to zero the factor goes to zero rather than infinity. The arm simply stops trying to move in directions it cannot move in. The maximum of that scaling is 1/(2 lambda), reached at sigma_i = lambda, so lambda directly caps how violently the solver reacts.

An illustrative number: with sigma = 0.01 and lambda = 0.05, the pseudo-inverse gain is 100, while DLS gives 0.01 / (0.0001 + 0.0025) which is about 3.85. A target jitter of 1 mm would demand 100 mm of joint-space-equivalent motion from the pseudo-inverse and under 4 mm from DLS. That factor of roughly twenty-six is why DLS is the workhorse.

Choosing lambda, and why a constant is not enough

Damping trades accuracy for stability. A large lambda is safe but slow, and it leaves a persistent error because the penalty biases dq toward zero. A small lambda converges quickly but oscillates near singularities. Buss notes that the right value depends on the multibody and the targets, and that in his experiments minimising average error produced heavy oscillation, so he raised lambda until oscillation became rare, accepting some loss in tracking accuracy.

Two refinements are standard. Adaptive damping sets lambda from the current manipulability, for example lambda^2 = lambda_max^2 times (1 – w/w0)^2 for w below a threshold w0 and zero otherwise (the singularity-robust inverse of Nakamura and Hanafusa is a close relative in spirit), so you pay the accuracy penalty only near singularities. Selectively damped least squares, Buss and Kim’s SDLS, damps each singular direction according to how far the end effector would move if that direction were followed. Error clamping, below, lets you keep lambda small in the first place.

Clamping: the part tutorials skip

Newton-style steps assume the linearisation is valid, which it is only for small moves. Buss recommends clamping the error vector so a single update never asks the end effector to move farther than a bound D_max, suggesting about half the length of a typical link, and also clamping the largest joint angle change per update. Both are one-liners, both prevent the solver from leaping across the workspace into a different branch, and both are absent from most blog code. The listing below includes them.

A runnable NumPy implementation

The code solves a three-link planar arm (redundant for a 2D position target, so it also shows null-space behaviour). Link lengths and tolerances are illustrative.

import numpy as np

L = np.array([0.5, 0.4, 0.3])  # three-link planar arm, metres

def fk(q):
    a = np.cumsum(q)
    return np.array([np.sum(L*np.cos(a)), np.sum(L*np.sin(a))]), a[-1]

def jac(q):
    a = np.cumsum(q)
    J = np.zeros((2, len(q)))
    for i in range(len(q)):
        J[0, i] = -np.sum(L[i:]*np.sin(a[i:]))
        J[1, i] =  np.sum(L[i:]*np.cos(a[i:]))
    return J

def clamp(v, d):
    n = np.linalg.norm(v)
    return v if n <= d else v*(d/n)

def dls_ik(q0, target, lam=0.05, dmax=0.2, max_step=0.3,
           tol=1e-5, iters=200, qmin=None, qmax=None):
    q = q0.copy()
    for k in range(iters):
        p, _ = fk(q)
        if np.linalg.norm(target - p) < tol:
            return q, k, True
        e = clamp(target - p, dmax)
        J = jac(q)
        dq = J.T @ np.linalg.solve(J@J.T + lam**2*np.eye(2), e)
        m = np.max(np.abs(dq))
        if m > max_step:
            dq *= max_step/m
        q = q + dq
        if qmin is not None:
            q = np.clip(q, qmin, qmax)
    return q, iters, False

q, k, ok = dls_ik(np.array([0.3, 0.3, 0.3]), np.array([0.7, 0.4]))
print(q, k, ok, fk(q)[0])

Running this on my machine from the seed (0.3, 0.3, 0.3) rad to the target (0.7, 0.4) converged in 5 iterations to q of about (-0.32, 1.05, 1.09) rad with position error under 1e-5. The iteration count is a property of this seed and these gains, not a benchmark. Try seeds near a stretched pose, and raise lambda from 0.0001 to 0.3: on the same machine, a start at (0, 0.001, 0) with target (0.9, 0.2) converged in 6 iterations at the two smaller damping values and 18 at lambda = 0.3, which illustrates the speed cost of heavy damping. Note that the step clamps are already doing a lot of the singularity work here; remove them and the undamped version is where the trouble shows up.

Singularities, Manipulability and Redundancy

A singularity is a configuration where J loses rank, so some direction of end-effector motion is instantaneously impossible and some joint motion produces no end-effector motion at all. The arm is not broken there; it simply cannot do what you ask along one axis, and an inverse that pretends otherwise commands enormous joint velocities.

Figure 3: Singularity handling. The Jacobian singular values feed a manipulability test; near a singularity the solver raises damping, while redundant arms also project a secondary objective into the null space.

The three families of singularity for a 6R arm

Industrial six-axis arms share three recurring singularities.

Boundary (elbow) singularity: the arm is fully stretched or fully folded, so the end effector cannot move radially. This is the sin(q2) = 0 case from our planar example, now with the shoulder and elbow aligned.
Wrist singularity: joints 4 and 6 become collinear when joint 5 is zero (or 180 degrees). Rotation about their common axis is then achievable by either joint, and rotation about the perpendicular axis by neither. This is the one that bites in practice, because it can sit in the middle of the useful workspace.
Shoulder (overhead) singularity: the wrist centre lies on the axis of joint 1, so joint 1 can swing the wrist centre through no change in position and the solution for joint 1 is undefined.

Straight-line Cartesian moves that pass near the wrist or shoulder singularity require joint velocities that grow without bound as the path approaches. Controllers respond with either a hard speed limit, which makes the tool crawl, or a protective stop. Neither is a solver bug; both reflect geometry.

Measuring the distance to a singularity

The simplest scalar is Yoshikawa’s manipulability measure, w = sqrt(det(J J transpose)), which equals the product of the singular values. It is zero at a singularity and, for a square J, equals |det J|. For the planar two-link arm w = L1 L2 |sin q2|, a smooth, interpretable quantity. Its geometric meaning is the volume of the velocity ellipsoid: the set of end-effector velocities reachable with unit-norm joint velocity.

Manipulability has one weakness: it multiplies singular values with different units for a mixed linear-angular Jacobian, so its absolute value is arbitrary. Two better-behaved alternatives are the smallest singular value sigma_min, which is what actually limits the worst-case direction, and the condition number sigma_max / sigma_min, which measures how elongated the ellipsoid is. In the worked example earlier, the singular values of roughly 0.73 and 0.27 give a condition number near 2.7, comfortable. A common rule of thumb is to start worrying when the condition number goes beyond some tens, but that threshold depends on the arm and units; calibrate it on your robot rather than adopting a number from a blog.

Singularity avoidance strategies

The options, from least to most invasive:

Damped inverse everywhere. The default, as above. It is cheap and keeps commands bounded, but it silently degrades tracking.
Adaptive damping tied to sigma_min. Raise lambda smoothly as sigma_min falls below a threshold. This preserves accuracy in the healthy workspace.
Path-level avoidance. Plan the Cartesian path to avoid singular regions, or choose a configuration branch that keeps them away. This is the best answer when you know the task in advance.
Gradient projection on manipulability. For redundant arms, use the null space to climb the manipulability gradient.
Accept a pose error. Some controllers deliberately let orientation drift near a wrist singularity and recover afterward.

Redundancy resolution with null-space projection

A seven-joint arm tracking a 6D pose has a one-dimensional family of solutions, the self-motion manifold, along which the elbow swivels while the hand holds still. The pseudo-inverse solution plus a null-space term exploits it:

dq = J-dagger e + (I – J-dagger J) z

The matrix (I – J-dagger J) projects an arbitrary joint-velocity vector z onto the null space of J, so that term changes joint angles without moving the end effector at first order. Choosing z as the gradient of a secondary cost lets you pursue a second objective for free. Typical choices are staying near the middle of the joint ranges, increasing manipulability, keeping a preferred posture, or moving the elbow away from an obstacle.

The cost is mild for the primary task and real for the secondary one: the projection is exact only at first order, and drift accumulates in the null space on long motions. Another subtlety is that the strict pseudo-inverse form is not damped; if you damp the primary task, the null-space projector derived from the damped inverse is no longer an exact projector, and the secondary task leaks into the primary one. Careful implementations compute the projector from the undamped SVD with a singular-value cutoff, or use a task-priority formulation with a small QP. A 7-DOF arm also has an additional self-motion singularity (internal singularities where the elbow aligns with the wrist), so redundancy reduces but does not remove the problem.

Constraints in Practice: Joint Limits, Orientation Error and Branch Selection

Joint limits

Real joints have limits, and unconstrained IK happily returns 400 degrees on a joint that only moves plus or minus 170. The three standard treatments are:

Clamp after each step. The simplest. It is what the listing does with np.clip. The weakness is that clamping a joint breaks the linear relationship between dq and e, so the solver can stall against a limit while still reporting an error. TRAC-IK’s documentation describes using random jumps in its Newton-based solver precisely to escape local minima caused by joint limits.
Weighted or active-set Jacobian. Freeze joints at their limits when the step would push further, drop their columns from J and re-solve for the rest. This converges to the true constrained optimum for a convex local model.
Quadratic programming. Formulate the step as a QP with bounds on dq, which is what differential-IK libraries such as Pink do: per its project description, tasks and limits form one quadratic program, minimising weighted squared task errors subject to velocity bounds and configuration limits read from the model. This handles limits, multiple prioritised tasks and posture regularisation in one place, at the price of a QP solve per tick.

For continuous-rotation joints, remember the wrap-around: a joint that goes from +179 to -179 degrees is a 2-degree move, not 358. Normalise differences into (-pi, pi] when computing errors.

Orientation error: quaternions versus axis-angle

Position error is a subtraction. Orientation error is not, and getting it wrong is the most common reason a 6D solver oscillates or spins.

Given a current rotation R and a target R_d, the error rotation is R_e = R_d R transposed (expressed in the world frame). To feed it to the solver you need a three-vector. The axis-angle vector, from the matrix logarithm, is the cleanest choice: its magnitude is the angle of the shortest rotation and it is smooth except at a 180-degree error. Euler-angle differences are unsafe because of gimbal lock and discontinuities. Quaternions work well but carry their own traps:

Double cover. q and -q represent the same rotation. If the dot product of the current and target quaternions is negative, you must negate one before computing the error, or the solver rotates the long way round (up to 360 degrees instead of up to 180).
Vector part as error. The vector part of the error quaternion is sin(theta/2) times the axis, which is roughly half the angle for small errors. That is fine for a proportional gain as long as you account for the factor of two in your Jacobian weighting.
Frame consistency. The orientation error must be expressed in the same frame as the angular rows of J. Mixing a world-frame Jacobian with a body-frame error is a silent bug that works for small motions and fails for large ones.
Normalisation drift. Integrating quaternions without renormalising causes slow drift; renormalise after each update.

Weighting also matters. Position is in metres and rotation in radians; a 1 cm position error and a 1 degree orientation error are comparable in practical effect but numerically 0.01 versus 0.017. Scale the rows of the task-space error (and Jacobian) so the solver’s tolerance means the same thing for both, or allow a looser tolerance on the rotational axes that do not matter for the task, such as the spin about a symmetric tool’s axis.

Branches and continuity

Closed-form IK returns a set; numerical IK returns whichever solution is nearest the seed. Both approaches can flip branches when run in a streaming loop. Suppose you track a moving target: at each tick you seed with the previous solution, and the numerical solver stays on the branch it began on, which is exactly the continuity you want. A closed-form solver run from scratch each tick can legitimately pick a different branch when two solutions cost the same, producing a sudden multi-radian joint move. Fix it by selecting among the returned solutions the one with the smallest weighted joint-space distance from the current configuration, with the largest-moving joint weighted most.

Even then, some Cartesian paths cannot be followed on one branch: moving through a singularity may force a branch change. Conversely, a few unusual arm geometries (called cuspidal in the kinematics literature) can change solution branch without ever crossing a singularity. The planner above the IK layer has to know this, which is part of why planners treat joint space, not task space, as primary.

The Solver Ecosystem in 2026

You rarely write IK from scratch. The practical question is which library fits your robot and your loop.

Solver	Method	Strengths	Watch out for
KDL (Orocos)	Newton-style Jacobian iteration (pseudo-inverse family) with joint-limit clamping	Ubiquitous, generic, ships as the default MoveIt plugin	Can fail near joint limits and local minima; seed-dependent
TRAC-IK	Runs a KDL-style Newton solver with random restarts concurrently with an SQP optimiser; returns whichever converges first	Higher success on constrained chains than plain KDL, per its documentation; selectable distance and manipulability criteria	Two threads per call; non-deterministic restarts; verify gains on your own arm
IKFast	Generates closed-form C++ per robot	Microsecond-class solves reported in the MoveIt tutorial, all branches returned	Per-robot generation workflow; 6 and 7 DOF only; legacy toolchain
BioIK	Memetic (evolutionary plus gradient) search with multi-objective goals	Handles multiple goals and awkward kinematics	Stochastic; heavier per call
Pink with Pinocchio	Differential IK as a quadratic program with task weighting and limits	Excellent for teleoperation, humanoids, streaming targets; Apache-2.0 licence	Local method that converges to a nearby optimum and does not solve global IK
cuRobo	GPU-batched optimisation for IK and trajectory optimisation	Many seeds in parallel; couples IK with collision checking	Needs a CUDA GPU; see our cuRobo article for details

I have kept the claims in that table to what each project’s own documentation says; benchmark on your arm and your targets before committing, because IK success rates depend heavily on joint limits, seeds and tolerance. The TRAC-IK page describes concurrent execution of a Newton-based solver with random jumps to escape joint-limit local minima and an SQP optimiser, returning when either converges, with optional secondary criteria (distance, manipulability, collision) for choosing a “best” solution. It does not give benchmark numbers on that page, and neither will I.

Choosing in practice

Industrial 6-axis with spherical wrist, fixed program: closed-form (vendor or IKFast). It is deterministic and gives you all branches.
7-DOF research arm on ROS 2: a numerical plugin (TRAC-IK or the framework’s default), with null-space posture control, or Pink for streaming targets. The ROS 2 distribution you run affects which plugin APIs exist, so check the notes in our ROS 2 Kilted to Lyrical migration guide before upgrading.
Teleoperation or VLA end-effector deltas at 30 to 100 Hz: differential IK (damped least squares or QP) seeded from the current joint state, running as the inner loop at 500 to 1000 Hz.
Planning with thousands of candidate goals: batched or GPU solvers, because you want many diverse solutions rather than one polished one.

Real-Time IK at 1 kHz

Direct answer: At 1 kHz you have one millisecond per tick for everything. A damped least squares step for a 7-DOF arm is a handful of small matrix operations and fits comfortably in a compiled language; the risks are allocation, iteration counts that vary with the input, and discontinuities. Run one bounded step per tick, warm-start from the previous solution, and never loop to convergence inside the control cycle.

Figure 4: A 1 kHz inverse kinematics control loop. The target filter and warm-started single-step solver run every tick; a slower thread handles branch selection and fallback.

What fits in the budget

Computing a geometric Jacobian for a 7-joint chain costs on the order of a few hundred floating-point operations. A 6 x 6 DLS solve via Cholesky is roughly n^3/3 operations for the factorisation, about seventy for n = 6 (illustrative arithmetic, not a benchmark). These are microsecond-scale operations in optimised C++ with Eigen or Pinocchio. The expensive parts are not the maths but the surrounding software: dynamic allocation, locks, logging and cache misses. In Python with NumPy, per-call overhead is typically tens of microseconds, so a loop of several NumPy calls is feasible at 1 kHz only on a non-real-time system; for hard real-time you want C++ or Rust.

Design rules

Single bounded step per tick. Differential IK is naturally incremental: the target moves a small amount, the joints follow. Convergence happens across ticks, not inside one.
Warm start. Seed every solve from the last commanded joint vector, not from measured joints, to avoid noise feeding back into the command.
Filter the target. A step change in the target is a step in error. Rate-limit position and use slerp for orientation, so the clamped error stays inside the linear regime.
Limit velocity and acceleration after IK. Treat the solver output as a desired velocity and pass it through the same limiters the drive expects. Jerk limits matter for real mechanisms.
Monitor sigma_min and residual. If sigma_min falls below a threshold or the residual grows for N ticks, raise damping, slow the target or stop. Make this explicit and logged.
Pre-allocate everything. No heap allocation in the loop; fixed-size matrices when the DOF is known at compile time.
Separate global from local. Global IK (all branches, collision filtering) belongs in a lower-rate thread. The 1 kHz loop only does local tracking.

Tolerances and convergence

Convergence tolerance should be set by the application, not by float precision. A typical position tolerance for assembly might be 0.1 mm, but a teleoperated pick may tolerate several millimetres. A tighter tolerance than your robot’s repeatability (which, for many industrial arms, is on the order of tens of micrometres, but check the datasheet) is meaningless. Use separate position and orientation tolerances, cap the iteration count and report failure explicitly. A solver that returns its best attempt without a success flag is a bug that waits for production.

For Newton-style methods, convergence near the solution is rapid because the linearisation becomes exact. For DLS the damping bounds the rate: with a constant lambda, convergence is linear with a rate that worsens as lambda grows relative to the singular values. The sample run above showed the practical effect, where 18 iterations replaced 6 when lambda increased by a factor of 6 over the middle setting.

Trade-offs, Gotchas, and What Goes Wrong

Local minima. Numerical IK is gradient-like and can stall. Joint limits, self-collision constraints and a poor seed all create configurations where the error cannot be reduced by any small step. Mitigations are random restarts (the TRAC-IK approach), multiple diverse seeds on a GPU, or a global planner that proposes the seed.

Branch flips. A solution continuous in the task space but jumping in the joint space shows up as a sudden multi-radian move. It typically follows from re-solving from scratch, from a target passing near a wrist singularity, or from a closed-form solver choosing among equal-cost branches. Seed from the previous state, and select by joint-space distance.

Silent failure. A solver that clamps the output into joint limits and returns the clamped joints will report success while the end effector is nowhere near the target. Always recompute forward kinematics on the returned solution and check the pose error against tolerance.

Frame and convention bugs. Standard versus modified DH, a tool offset applied in the wrong frame, a world-frame Jacobian paired with a body-frame error, and quaternion order (xyzw in ROS and Eigen’s coefficient order, wxyz in some other libraries) are frequent causes of “it works except it’s mirrored” bugs. Write a test that compares your IK against forward kinematics on a thousand random reachable poses, and include a round trip through the URDF-loaded model rather than only hand-coded parameters.

Calibration error. Kinematic parameters in a URDF are nominal. Real arms differ by link-length errors and joint offsets of fractions of a millimetre or degree, and that error appears directly in the tool pose. IK cannot fix it; kinematic calibration can.

Overdamping. The failure that looks like success. A large lambda gives smooth, safe, wrong motion: the arm approaches the target and settles with a persistent error. Make the residual visible.

Over-trusting the Jacobian near workspace boundary. Targets slightly outside the workspace make an iterative solver hunt around the boundary. Detect unreachable targets cheaply by distance tests (for example, wrist-centre distance against the maximum reach) before iterating.

Anti-patterns. Calling a blocking multi-hundred-iteration IK from a controller callback. Using the pseudo-inverse with no damping in any loop that touches real hardware. Treating “IK returned true” as proof that the pose was reached. Re-seeding with a fixed home pose on every call, which makes the solver both slower and branch-unstable.

Practical Recommendations

Start from the robot, not the algorithm. If your arm has a spherical wrist and the task is fixed, use a closed-form solver and spend your effort on branch selection and joint-limit checking. If your arm is redundant, mobile or humanoid, use a differential IK solver with damping, task weights and limits, and treat posture regularisation as part of the problem statement. In both cases, wrap the solver in the same verification layer: forward kinematics on the result, joint-limit check, collision check, and a continuity check against the previous command.

Treat the digital twin as the test bench. The same URDF or OpenUSD model that drives simulation can generate the random reachable poses, singular configurations and workspace-boundary targets that exercise your solver; running IK regression tests in simulation before hardware is cheap and finds frame-convention bugs early.

Checklist:

[ ] Forward-kinematics round trip passes on 1,000+ random reachable poses.
[ ] Orientation error uses a rotation-matrix logarithm or sign-corrected quaternion, in the same frame as the Jacobian.
[ ] Damping is adaptive to sigma_min, or error and step clamping are in place.
[ ] Joint limits are enforced inside the solver and re-verified after.
[ ] The solver reports success, iteration count and residual, and the caller checks them.
[ ] Branch selection uses joint-space distance from the current state.
[ ] Unreachable targets are detected before iteration.
[ ] The 1 kHz loop does a single bounded step with no allocation.
[ ] Singularity proximity (sigma_min or condition number) is logged and alarmed.
[ ] Kinematic calibration data, not nominal parameters, is used for precision tasks.

Frequently Asked Questions

What is inverse kinematics in robotics?

Inverse kinematics computes the joint angles that place a robot’s end effector at a desired position and orientation. It is the reverse of forward kinematics, which maps known joint angles to a pose. Because the mapping can be many-to-one, an IK problem may have no solution, one solution, several discrete solutions or an infinite family, depending on the arm geometry, joint limits and target.

What is the difference between analytical and numerical inverse kinematics?

Analytical (closed-form) IK derives explicit formulas, typically by exploiting a spherical wrist, and returns all solutions in constant time with no seed. Numerical IK iterates, usually with the Jacobian, from an initial guess and returns one solution near that guess. Analytical is faster and complete but needs special geometry; numerical is general, handles redundancy and extra constraints, but can stall in local minima.

Why does the Jacobian pseudo-inverse fail near singularities?

The pseudo-inverse scales each direction by the reciprocal of its singular value. As an arm nears a singularity one singular value approaches zero, so the corresponding gain approaches infinity and small errors demand huge joint velocities. Damped least squares replaces 1/sigma with sigma/(sigma^2 + lambda^2), which tends to zero as sigma does, so commands stay bounded at the cost of some tracking accuracy.

How do I choose the damping factor in damped least squares?

Choose lambda large enough to avoid oscillation near singularities, and small enough that convergence stays acceptable. Buss reports that the right value depends on the robot and the targets, and that clamping the error vector allows a smaller lambda. A practical approach is adaptive damping that rises as the smallest singular value falls, tuned on recorded trajectories from your robot, with the residual logged to detect overdamping.

How many inverse kinematics solutions does a 6-axis robot arm have?

A general 6R chain has up to sixteen solutions, but common industrial arms with a spherical wrist usually have up to eight, from the combinations of shoulder left or right, elbow up or down, and wrist flipped or not. Joint limits typically remove several of these, so a controller often sees two to four feasible configurations. The exact count depends on offsets and limits in your model.

Should I use IKFast, TRAC-IK or a differential IK library?

Use closed-form (IKFast or the vendor solver) for a 6-DOF arm with fixed tasks where determinism and speed matter. Use TRAC-IK or a similar numerical plugin when you need a generic solver that copes with joint limits. Use a differential IK library such as Pink for streaming targets, redundant or humanoid robots, and prioritised tasks. Validate success rate and latency on your own arm.