Inverse Jacobian Methods

You cannot invert the forward kinematics analytically for a general chain. But you can linearize it locally: near the current configuration, the Jacobian says how the hand moves per unit of joint motion. So take the error between where the hand is and where you want it, ask the Jacobian for the joint motion that would reduce that error, take a small step, and repeat. This is Newton-Raphson on a vector-valued function, and it is structurally identical to the gradient descent loop from track three.

e = target - forward_kinematics(q) qdot = J^+ e pseudoinverse of the Jacobian q <- q + alpha * qdot step, then recompute J at the new q

Iterative IK with damped least squares

import numpy as np

def inverse_kinematics(target, q0, lengths, damping=0.05,
                       alpha=1.0, tol=1e-4, max_iters=200):
    """Damped-least-squares IK. Returns (solution, error_history, converged)."""
    q = np.array(q0, dtype=float)
    history = []

    for _ in range(max_iters):
        error = target - forward_kinematics(q, lengths)
        norm = float(np.linalg.norm(error))
        history.append(norm)
        if norm < tol:
            return q, history, True

        J = jacobian_analytic(q, lengths)          # (2, n)

        # Damped least squares (Levenberg-Marquardt):
        #   J^T (J J^T + lambda^2 I)^{-1}
        # The damping term keeps the inverse finite when J loses rank. It trades
        # exactness for stability — near a singularity you accept a slightly wrong
        # direction instead of an unbounded joint velocity.
        JJt = J @ J.T
        damped = JJt + (damping ** 2) * np.eye(JJt.shape[0])
        qdot = J.T @ np.linalg.solve(damped, error)

        q = q + alpha * qdot

    return q, history, False

Three ways to invert, and when each is right

Jacobian transpose: `qdot = J^T e`. Cheap, always stable, no matrix inverse — but it is gradient descent on squared error, so convergence is slow. Use when compute is scarce or robustness matters more than speed.
Moore-Penrose pseudoinverse: `qdot = J^+ e`. Gives the minimum-norm exact solution, converges fast, and blows up near singularities. Use away from singular configurations.
Damped least squares: the code above. Adds `lambda^2 I` before inverting, which bounds the solution near singularities at the cost of a small steady-state error. This is the production default, and the damping term is exactly ridge regularization from track three.

One idea wearing three hats

The damping constant here, the ridge penalty in regression, and weight decay in a neural network are the same mathematical device: add a small multiple of the identity to make an ill-conditioned problem well-posed. Recognizing that is worth more than memorizing any of the three.

Convergence behavior is worth watching directly. The error history should fall roughly geometrically; a plateau means you are at a local structure the linearization cannot escape, and oscillation means alpha is too large — the same two readings you learned to make from a training loss curve.

Try each one before opening the solution. Getting it wrong first is most of where the learning happens.

import numpy as np def ik(target, q0, lengths, mode="damped", damping=0.05, alpha=1.0, tol=1e-4, max_iters=500): q = np.array(q0, dtype=float) for i in range(max_iters): error = target - forward_kinematics(q, lengths) if np.linalg.norm(error) < tol: return q, i J = jacobian_analytic(q, lengths) if mode == "transpose": qdot = J.T @ error elif mode == "pinv": qdot = np.linalg.pinv(J) @ error else: JJt = J @ J.T qdot = J.T @ np.linalg.solve(JJt + damping ** 2 * np.eye(2), error) q = q + alpha * qdot return q, max_iters for mode in ("transpose", "pinv", "damped"): q, iters = ik(np.array([1.2, 0.8]), [0.3, 0.4], np.array([1.0, 1.0]), mode=mode) print(f"{mode:10} {iters} iterations")

import numpy as np lengths = np.array([1.0, 1.0]) q, iters = ik(np.array([1.2, 0.8]), [0.3, 0.4], lengths, mode="damped") assert np.linalg.norm(forward_kinematics(q, lengths) - np.array([1.2, 0.8])) < 1e-3 print("ok")

import numpy as np for damping in (0.0, 0.1): q, peak, max_err = np.array([0.2, 0.3]), 0.0, 0.0 for s in np.linspace(0.5, 1.95, 300): target = np.array([s, 0.05]) error = target - forward_kinematics(q, lengths) J = jacobian_analytic(q, lengths) JJt = J @ J.T qdot = J.T @ np.linalg.solve(JJt + damping ** 2 * np.eye(2), error) peak = max(peak, float(np.abs(qdot).max())) q = q + 0.2 * qdot max_err = max(max_err, float(np.linalg.norm(error))) print(f"damping={damping}: peak qdot={peak:.1f} max tracking error={max_err:.4f}")

import matplotlib.pyplot as plt q, history, converged = inverse_kinematics(np.array([1.2, 0.8]), [0.3, 0.4], lengths) plt.semilogy(history); plt.xlabel("iteration"); plt.ylabel("|error|"); plt.show() print("converged:", converged, "in", len(history), "iterations")