# Simulating an Inverse Kinematics Chain

Inverse kinematics is a fun application of linear algebra and calculus that is extremely relevant in computer animation. Given a rigged skeleton, inverse kinematics can solve the joint angles (or similar parameters for other types of joints) necessary to minimize the distance between the end effector (i.e. joint at the end of the chain) and target. Solving this system essentially boils down to the following system of equations:

$\vec{e}=J\vec{\Delta\theta}$

where $\vec{e}$ is called the “error vector,” and describes the velocity of the end effector necessary to converge on the target; $\Delta\vec{\theta}$ is the vector of joint angle velocities that will result in the correct end effector velocity (this is the unknown part), and the Jacobian matrix $J$, that encodes the change of basis from joint angle space to Cartesian coordinate space for the partial derivatives. Alternatively, you can view this as a transformation matrix that takes the rate of change of the joint angles and relates it to the rate of change of the Cartesian coordinates of the end effector joint.

It is common to have more columns than rows in a system like this (especially in a chain), because the there will be few end effectors and many joints. Because the matrix is not invertible, our next best thing is the Moore-Penrose psuedoinverse which gives a least-squares solution to the problem. This is desirable because there are many cases where the chain cannot converge on the target, so we want the best possible attempt. The results in the following equation:

$\Delta\vec{\theta}=J^T(JJ^T)^{-1}\vec{e}$

This form of the psuedoinverse assumes that you have full row rank in your matrix $J$, otherwise $JJ^T$ will be singular an non-invertible. This presents a problem because there do exist singularities in certain joint configurations. A common solution to this problem is to damp this matrix by a small factor to ensure invertibility. This method is called the Damped Least Squares method, and the requires only a simple change to the previous equation:

$\Delta\vec{\theta}=J^T(JJ^T + \lambda I)^{-1}\vec{e}$

Here, $\lambda$ is just an arbitrarily chosen constant that keeps the matrix from every going singular. The bigger the constant, the less realistic the motion near singular configurations. On the other hand, a higher constant results in more numerical stability in the equations. The demo I wrote allows you to set an arbitrary $\lambda$ value.

I also implemented two other solvers, a Jacobian Transpose solver and a Cyclic-Coordinate Decent solver. The former cheats by using the transpose of the Jacobian to solve the system rather than the psuedoinverse. It’s a simple solution but results in a lot of jittery motion and isn’t very realistic. CCD essentially loops through each joint multiple times and orients the child chain in the direction of the target. With enough iterations, this converges on a solution, but it only works in situations where joints are 1DOF. It is quite simple to understand and there are several articles online that explain it well.

Finally, the demo is written in Python using the pyglet library. User interactions is almost exclusively through a console interface (activated with ~). Included in the zip is a text file with the syntax for these commands.

Source Archive

For a more in-depth introduction to this topic, I found this paper extremely helpful.