# Assignment 4 - Inverse Kinematics

Due Wednesday, 10/21 Friday, 10/23 at 23:59:59. You must work individually.

## Goal

In this assignment, you will be implementing inverse kinematics with your own nonlinear solver:

• GDLS: Gradient Descent with (Backtracking) Line Search
• NM: Newton’s Method

Grades will be given out depending on how far you get along in the assignment. For example, if your program works with $$4$$ links, then you will get full points for the $$1$$- and $$2$$-link tasks.

Here is an example with $$10$$ links:

## Associated Labs

• Lab 0 (required): Setting Up Your Development Environment
• Lab 5 (required): Forward Kinematics / Hierarchical Modeling
• Lab 6 (required): Gradient Descent and Newton’s Method

The best starting point is Lab 5. Implement forward kinematics fully before attempting inverse kinematics.

## List of Parameters

For quick reference, this is a listing of all the parameters to be used in the optimizers:

• Objective Function
• $$w_\text{tar}$$: 1e3
• $$w_\text{reg}$$: 1e0
• Gradient Descent with [Backtracking] Line Search
• iterMax: 10*n, where n is the number of links
• alphaInit: 1e-4
• tol: 1e-6
• Newton’s Method
• iterMax: 10*n, where n is the number of links
• tol: 1e-6

We will start slowly. First, we will solve inverse kinematics with only $$1$$ joint using Gradient Descent with [Backtracking] Line Search (GDLS).

Each link is of length $$1.0$$, so place an end-effector at xlocal = [1 0]'. Keep in mind that there should be a homogeneous coordinate of $$1$$ for this point. In your code, it is a good idea to use Vector3d instead of Vector2d and explicitly store the homogeneous coordinate.

The world-space position, $$\vec{x}$$, of the end-effector and its derivatives are:

$\displaylines{ \vec{x}(\theta) = T R(\theta) \vec{x}_\text{local}\\ \vec{x}'(\theta) = T R'(\theta) \vec{x}_\text{local}\\ \vec{x}''(\theta) = T R''(\theta) \vec{x}_\text{local}, }$

where $$R(\theta)$$ is the usual rotation matrix about the Z-axis:

$\displaylines{ R(\theta) = \begin{pmatrix} \cos(\theta) & -\sin(\theta) & 0\\ \sin(\theta) & \cos(\theta) & 0\\ 0 & 0 & 1 \end{pmatrix}\\ R'(\theta) = \begin{pmatrix} -\sin(\theta) & -\cos(\theta) & 0\\ \cos(\theta) & -\sin(\theta) & 0\\ 0 & 0 & 0 \end{pmatrix}\\ R''(\theta) = \begin{pmatrix} -\cos(\theta) & \sin(\theta) & 0\\ -\sin(\theta) & -\cos(\theta) & 0\\ 0 & 0 & 0 \end{pmatrix}. }$

The second derivative, $$\vec{x}''(\theta)$$ and $$R''(\theta)$$, will not be needed until we use Newton’s Method.

In the original rest state (figure above), the numerical values for the $$1$$-link system are (homogeneous coordinate not shown):

$\displaylines{ \vec{x}(\theta) = \begin{pmatrix} 1 \\ 0 \end{pmatrix}\\ \vec{x}'(\theta) = \begin{pmatrix} 0 \\ 1 \end{pmatrix}\\ \vec{x}''(\theta) = \begin{pmatrix} -1 \\ 0 \end{pmatrix}. }$

For Gradient Descent, with or without line search, we need the objective function value $$f(\theta)$$ and its gradient, $$\vec{g}(\theta)$$. (We don’t need the Hessian, $$H(\theta),$$ until we get to Newton’s Method.) The objective value is

$\displaylines{ f(\theta) = \frac{1}{2} w_\text{tar} \Delta \vec{x}^\top \Delta \vec{x} + \frac{1}{2} w_\text{reg} \Delta \theta^2\\ \Delta \vec{x} = \vec{x}(\theta) - \vec{x}_\text{target} \qquad\qquad\quad\;\;\,\\ \Delta \theta = \theta - \theta_\text{rest}. \qquad\qquad\qquad\quad\;\; }$

The two weights, $$w_\text{tar}$$ and $$w_\text{reg}$$, can be set to arbitrary positive values, depending on the desired strength of the “targeting” energy and the “regularizer” energy.

• The targeting energy causes the optimizer to minimize the squared distance between the end-effector and the target position.
• The regularizer energy causes the optimizer to prefer specific angles. For this assignment, you should set $$\theta_\text{rest} = 0$$.

You can use $$w_\text{tar} = 1e3$$ and $$w_\text{reg} = 1e0$$ for all tasks of this assignment.

The gradient is a scalar for this $$1$$-link task:

$g(\theta) = w_\text{tar} \Delta \vec{x}^\top \vec{x}' + w_\text{reg} \Delta \theta.$

The Hessian (not needed yet) is also a scalar for this $$1$$-link task:

$H(\theta) = w_\text{tar} (\vec{x}'^\top \vec{x}' + \Delta \vec{x}^\top \vec{x}'') + w_\text{reg}.$

With the initial angle set to $$\theta = 0$$, and with the goal position set to $$\vec{x}_\text{goal} = (0 \; 1)^\top$$, you should get the following values:

$f(\theta) = 1000, \quad g(\theta) = -1000, \quad H(\theta) = 1.$

Run GDLS for $$10$$ iterations (i.e., 10*n iterations). The returned answer should be close to, but not exactly $$\theta = \pi/2$$. The exact answer will differ somewhat, based on the exact implementation of GDLS.

User Interface

Now let the end-effector goal be driven by the mouse. Set up your code so that pressing the space bar toggles between forward kinematics and inverse kinematics.

• While in forward kinematics mode, you should be able to increase/decrease the joint angle(s) by pressing the . and , keys.
• While in inverse kinematics mode, you should be able to move the mouse and drive the joint angles.
• In either mode, pressing the r key should reset all the angles to 0.0.

The answer returned by GDLS will be visibly wrong because GDLS converges slowly. Therefore, we’re going to run a few iterations of NM to improve the answer. However, we must be careful since Newton’s Method can often get stuck in a local minimum. We’re going to take the following approach:

1. Run GDLS to get $$\theta_\text{GDLS}$$ and $$f(\theta_\text{GDLS})$$
2. Run NM, starting the search with $$\theta_\text{GDLS}$$, to get $$\theta_\text{NM}$$ and $$f(\theta_\text{NM})$$
3. Check if NM improved the answer: $$f(\theta_\text{NM})$$ < $$f(\theta_\text{GDLS})$$
• If true, set $$\theta = \theta_\text{NM}$$
• If false, set $$\theta = \theta_\text{GDLS}$$

This ensures that we use the result of NM only if it improved upon the result of GDLS. With the combined approach above, you should get a visibly better result within a few Newton iterations.

Please output to the console the number of iterations taken by GDLS and by NM.

Next, we’re going to write a version that supports $$2$$ links. If you’re feeling ambitious, you can skip to the $$4$$-link or the $$n$$-link tasks.

With two links, the minimization variable, $$\vec{\theta}$$, is now a vector. The world-space position, $$\vec{x}$$, of the end-effector is:

$\vec{x}(\vec\theta) = T_0 R_0(\theta_0) T_1 R_1(\theta_1) \vec{x}_\text{local}.$

The derivatives are no longer scalars. The first derivative is a matrix of size $$2 \times 2$$ (or $$1 \times 2$$ block vector of $$2 \times 1$$ vectors):

$\displaylines{ X'(\vec\theta) = \begin{pmatrix} \vec{x}'_0 & \vec{x}'_1 \end{pmatrix}\\ \vec{x}'_0 = T_0 R'_0(\theta_0) T_1 R_1(\theta_1) \vec{x}_\text{local}\\ \vec{x}'_1 = T_0 R_0(\theta_0) T_1 R'_1(\theta_1) \vec{x}_\text{local}. }$

The second derivative is a matrix of size $$4 \times 2$$ (or $$2 \times 2$$ block matrix of $$2 \times 1$$ vectors):

$\displaylines{ X''(\vec\theta) = \begin{pmatrix} \vec{x}''_{00} & \vec{x}''_{01} \\ \vec{x}''_{10} & \vec{x}''_{11} \end{pmatrix}\\ \vec{x}''_{00} = T_0 R''_0(\theta_0) T_1 R_1(\theta_1) \vec{x}_\text{local}\\ \vec{x}''_{01} = T_0 R'_0(\theta_0) T_1 R'_1(\theta_1) \vec{x}_\text{local}\\ \vec{x}''_{10} = T_0 R'_0(\theta_0) T_1 R'_1(\theta_1) \vec{x}_\text{local}\\ \vec{x}''_{11} = T_0 R_0(\theta_0) T_1 R''_1(\theta_1) \vec{x}_\text{local}. }$

The objective function, $$f(\vec\theta)$$, is still a scalar. In fact, the objective function for an optimization problem always returns a scalar. Compared to the $$1$$-link case, the only difference is that $$\Delta \vec\theta$$, the deviation of the current joint angles from the rest joint angles, is now a vector. Therefore, instead of squaring this difference, we take the dot product by itself (i.e., transpose and multiply).

$f(\vec\theta) = \frac{1}{2} w_\text{tar} \Delta \vec{x}^\top \Delta \vec{x} + \frac{1}{2} w_\text{reg} \Delta \vec\theta^\top \Delta \vec\theta.$

Now the gradient is a vector:

$\displaylines{ \vec{g}(\vec\theta) = w_\text{tar} \Delta \vec{x}^\top X' + w_\text{reg} \Delta \vec\theta \qquad\qquad\quad\;\\ = w_\text{tar} \begin{pmatrix} \Delta \vec{x}^\top \vec{x}'_0\\ \Delta \vec{x}^\top \vec{x}'_1 \end{pmatrix} + w_\text{reg} \begin{pmatrix} \Delta \theta_0\\ \Delta \theta_1 \end{pmatrix}. }$

And the Hessian is a matrix:

$\displaylines{ H(\vec\theta) = w_\text{tar} (X'^\top X' + \Delta \vec{x}^\top \odot X'') + w_\text{reg} I \qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\;\,\\ = w_\text{tar} \left( \begin{pmatrix} \vec{x}_0^{'\top} \vec{x}'_0 & \vec{x}_0^{'\top} \vec{x}'_1\\ \vec{x}_1^{'\top} \vec{x}'_0 & \vec{x}_1^{'\top} \vec{x}'_1 \end{pmatrix} + \begin{pmatrix} \Delta \vec{x}^\top \vec{x}''_{00} & \Delta \vec{x}^\top \vec{x}''_{01}\\ \Delta \vec{x}^\top \vec{x}''_{10} & \Delta \vec{x}^\top \vec{x}''_{11} \end{pmatrix} \right) + w_\text{reg} \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}. }$

The $$\odot$$ symbol above is used to mean that $$\Delta \vec{x}$$ is dotted with each vector element of $$X''$$.

With the initial angle set to $$\vec{\theta} = (0 \; 0)^\top$$, and with the goal position set to $$\vec{x}_\text{goal} = (1 \; 1)^\top$$, you should get the following values:

$\vec{x}(\vec{\theta}) = \begin{pmatrix} 2\\ 0 \end{pmatrix}, \quad \vec{x}'(\vec{\theta}) = \begin{pmatrix} 0 & 0\\ 2 & 1\\ \end{pmatrix}, \quad \vec{x}''(\vec{\theta}) = \begin{pmatrix} -2 & -1\\ 0 & 0\\ -1 & -1\\ 0 & 0 \end{pmatrix}, \quad f(\vec{\theta}) = 1000, \quad g(\vec{\theta}) = \begin{pmatrix} -2000\\ -1000 \end{pmatrix}, \quad H(\vec{\theta}) = \begin{pmatrix} 2001 & 1000\\ 1000 & 1 \end{pmatrix}.$

Wrapping $$\theta$$

Since the minimization variables are joint angles, these variables may wrap around beyond (clockwise or counterclockwise) many times over. To prevent this, add the following pseudocode to map the angles to be between $$-\pi$$ and $$+\pi$$:

while angle > pi
angle -= 2*pi
while angle < -pi
angle += 2*pi

Debugging the Derivatives

One of the biggest sources of bugs is the derivatives, so it is helpful to have a routine that tests them against finite differences. To test the gradient vector, use the following code inside the main loop of GDLS:

double h = 1e-7;
VectorXd g_(n);
for(int i = 0; i < n; ++i) {
VectorXd x_ = x;
x_(i) += h;
double f_ = objective->evalObjective(x_);
g_(i) = (f_ - f)/h;
}

The finite differenced version of the gradient, g_, should closely match the analytical version of the gradient, g.

To test the Hessian matrix, use the following code inside the main loop of NM:

double h = 1e-7;
MatrixXd H_(n, n);
for(int i = 0; i < n; ++i) {
VectorXd x_ = x;
x_(i) += h;
VectorXd g_(n);
objective->evalObjective(x_, g_);
H_.col(i) = (g_ - g)/h;
}

The finite differenced version of the Hessian, H_, should closely match the analytical version of the Hessian, H.

In this task, you’re going to hard code a $$4$$-link system, which follows the same steps as in the $$2$$-link system. If you’re feeling ambitious, you can skip this task and go straight to the $$n$$-link system.

If you skipped $$2$$-link IK, make sure to read the Wrapping and Debugging sections.

The minimization variable is now $$\vec\theta = (\theta_0 \; \theta_1 \; \theta_2 \; \theta_3)^\top$$. The world-space position, $$\vec{x}$$, of the end-effector is:

$\vec{x}(\vec\theta) = T_0 R_0(\theta_0) T_1 R_1(\theta_1) T_2 R_2(\theta_2) T_3 R_3(\theta_3) \vec{x}_\text{local}.$

The first derivative is a matrix of size $$2 \times 4$$ (or $$1 \times 4$$ block vector of $$2 \times 1$$ vectors):

$\displaylines{ X'(\vec\theta) = \begin{pmatrix} \vec{x}'_0 & \vec{x}'_1 & \vec{x}'_2 & \vec{x}'_3 \end{pmatrix}\\ \vec{x}'_0 = T_0 R'_0(\theta_0) T_1 R_1(\theta_1) T_2 R_2(\theta_2) T_3 R_3(\theta_3) \vec{x}_\text{local}\\ \vec{x}'_1 = T_0 R_0(\theta_0) T_1 R'_1(\theta_1) T_2 R_2(\theta_2) T_3 R_3(\theta_3) \vec{x}_\text{local}\\ \vec{x}'_2 = T_0 R_0(\theta_0) T_1 R_1(\theta_1) T_2 R'_2(\theta_2) T_3 R_3(\theta_3) \vec{x}_\text{local}\\ \vec{x}'_3 = T_0 R_0(\theta_0) T_1 R_1(\theta_1) T_2 R_2(\theta_2) T_3 R'_3(\theta_3) \vec{x}_\text{local}. }$

The second derivative is a matrix of size $$8 \times 4$$ (or $$4 \times 4$$ block vector of $$2 \times 1$$ vectors):

$\displaylines{ X''(\vec\theta) = \begin{pmatrix} \vec{x}''_{00} & \vec{x}''_{01} & \vec{x}''_{02} & \vec{x}''_{03}\\ \vec{x}''_{10} & \vec{x}''_{11} & \vec{x}''_{12} & \vec{x}''_{13}\\ \vec{x}''_{20} & \vec{x}''_{21} & \vec{x}''_{22} & \vec{x}''_{23}\\ \vec{x}''_{30} & \vec{x}''_{31} & \vec{x}''_{32} & \vec{x}''_{33} \end{pmatrix}\\ \vec{x}''_{00} = T_0 R''_0(\theta_0) T_1 R_1(\theta_1) T_2 R_2(\theta_2) T_3 R_3(\theta_3) \vec{x}_\text{local}\\ \vec{x}''_{01} = T_0 R'_0(\theta_0) T_1 R'_1(\theta_1) T_2 R_2(\theta_2) T_3 R_3(\theta_3) \vec{x}_\text{local}\\ \vec{x}''_{02} = T_0 R'_0(\theta_0) T_1 R_1(\theta_1) T_2 R'_2(\theta_2) T_3 R_3(\theta_3) \vec{x}_\text{local}\\ \vdots\\ \vec{x}''_{32} = T_0 R_0(\theta_0) T_1 R_1(\theta_1) T_2 R'_2(\theta_2) T_3 R'_3(\theta_3) \vec{x}_\text{local}\\ \vec{x}''_{33} = T_0 R_0(\theta_0) T_1 R_1(\theta_1) T_2 R_2(\theta_2) T_3 R''_3(\theta_3) \vec{x}_\text{local}. }$

The scalar objective function is:

$f(\vec\theta) = \frac{1}{2} w_\text{tar} \Delta \vec{x}^\top \Delta \vec{x} + \frac{1}{2} w_\text{reg} \Delta \vec\theta^\top \Delta \vec\theta.$

The gradient, expressed as a column vector is:

$\displaylines{ \vec{g}(\vec\theta) = w_\text{tar} \Delta \vec{x}^\top X' + w_\text{reg} \Delta \vec\theta \qquad\qquad\quad\quad\\ = w_\text{tar} \begin{pmatrix} \Delta \vec{x}^\top \vec{x}'_0\\ \Delta \vec{x}^\top \vec{x}'_1\\ \Delta \vec{x}^\top \vec{x}'_2\\ \Delta \vec{x}^\top \vec{x}'_3 \end{pmatrix} + w_\text{reg} \begin{pmatrix} \Delta \theta_0\\ \Delta \theta_1\\ \Delta \theta_2\\ \Delta \theta_3 \end{pmatrix}. }$

And the Hessian is a matrix:

$\displaylines{ H(\vec\theta) = w_\text{tar} (X'^\top X' + \Delta \vec{x}^\top \odot X'') + w_\text{reg} I \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\\ = w_\text{tar} \left( \begin{pmatrix} \vec{x}_0^{'\top} \vec{x}'_0 & \vec{x}_0^{'\top} \vec{x}'_1 & \vec{x}_0^{'\top} \vec{x}'_2 & \vec{x}_0^{'\top} \vec{x}'_3\\ \vec{x}_1^{'\top} \vec{x}'_0 & \vec{x}_1^{'\top} \vec{x}'_1 & \vec{x}_1^{'\top} \vec{x}'_2 & \vec{x}_1^{'\top} \vec{x}'_3\\ \vec{x}_2^{'\top} \vec{x}'_0 & \vec{x}_2^{'\top} \vec{x}'_1 & \vec{x}_2^{'\top} \vec{x}'_2 & \vec{x}_2^{'\top} \vec{x}'_3\\ \vec{x}_3^{'\top} \vec{x}'_0 & \vec{x}_3^{'\top} \vec{x}'_1 & \vec{x}_3^{'\top} \vec{x}'_2 & \vec{x}_3^{'\top} \vec{x}'_3 \end{pmatrix} + \begin{pmatrix} \Delta \vec{x}^\top \vec{x}''_{00} & \Delta \vec{x}^\top \vec{x}''_{01} & \Delta \vec{x}^\top \vec{x}''_{02} & \Delta \vec{x}^\top \vec{x}''_{03}\\ \Delta \vec{x}^\top \vec{x}''_{10} & \Delta \vec{x}^\top \vec{x}''_{11} & \Delta \vec{x}^\top \vec{x}''_{12} & \Delta \vec{x}^\top \vec{x}''_{13}\\ \Delta \vec{x}^\top \vec{x}''_{20} & \Delta \vec{x}^\top \vec{x}''_{21} & \Delta \vec{x}^\top \vec{x}''_{22} & \Delta \vec{x}^\top \vec{x}''_{23}\\ \Delta \vec{x}^\top \vec{x}''_{30} & \Delta \vec{x}^\top \vec{x}''_{31} & \Delta \vec{x}^\top \vec{x}''_{32} & \Delta \vec{x}^\top \vec{x}''_{33} \end{pmatrix} \right) + w_\text{reg} \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{pmatrix}. }$

The $$\odot$$ symbol above is used to mean that $$\Delta \vec{x}$$ is dotted with each vector element of $$X''$$.

With the initial angle set to $$\vec{\theta} = (0 \; 0 \; 0 \; 0)^\top$$, and with the goal position set to $$\vec{x}_\text{goal} = (3 \; 1)^\top$$, you should get the following values:

$\vec{x}(\vec{\theta}) = \begin{pmatrix} 4\\ 0 \end{pmatrix}, \quad \vec{x}'(\vec{\theta}) = \begin{pmatrix} 0 & 0 & 0 & 0\\ 4 & 3 & 2 & 1 \end{pmatrix}, \quad \vec{x}''(\vec{\theta}) = \begin{pmatrix} -4 & -3 & -2 & -1\\ 0 & 0 & 0 & 0\\ -3 & -3 & -2 & -1\\ 0 & 0 & 0 & 0\\ -2 & -2 & -2 & -1\\ 0 & 0 & 0 & 0\\ -1 & -1 & -1 & -1\\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad f(\vec{\theta}) = 1000, \quad g(\vec{\theta}) = \begin{pmatrix} -4000\\ -3000\\ -2000\\ -1000 \end{pmatrix}, \quad H(\vec{\theta}) = \begin{pmatrix} 12001 & 9000 & 6000 & 3000\\ 9000 & 6001 & 4000 & 2000\\ 6000 & 4000 & 2001 & 1000\\ 3000 & 2000 & 1000 & 1 \end{pmatrix}.$

Convert your code to support an arbitrary number of links. Note that as you increase $$n$$, the code will get slower and slower. Make sure to run your code in release mode. It is recommended that you first work on GDLS and then NM.

If you skipped $$2$$-link IK, make sure to read the Wrapping and Debugging sections.

## Point breakdown

• [10] Forward kinematics
• [5] Hard-coded $$1$$-link IK using GDLS
• [5] Hard-coded $$1$$-link IK using GDLS+NM
• [10] Hard-coded $$2$$-link IK using GDLS
• [10] Hard-coded $$2$$-link IK using GDLS+NM
• [15] Hard-coded $$4$$-link IK using GDLS
• [15] Hard-coded $$4$$-link IK using GDLS+NM
• [20] $$n$$-link IK using GDLS
• [+20] $$n$$-link IK using GDLS+NM
• [10] Coding style and general execution

Total: 100 points + 20 bonus points

You can skip $$1$$-, $$2$$-, and $$4$$-link IK and go straight to $$n$$-link IK. However, if you only complete $$n$$-link IK using GDLS, you will not get the $$5+10+15=30$$ points for the GDLS+NM tasks.

## What to hand in

Please output to the console the number of iterations taken by GDLS and by NM.

If you’re using Mac/Linux, make sure that your code compiles and runs by typing:

> mkdir build
> cd build
> cmake ..
> make
> ./A4 <SHADER DIR>

If you’re using Windows, make sure your code builds using the steps described in Lab 0.

• Put your name on the program window.
• Include an ASCII README that includes
• Do not hand in the executable, old save files (*.~), or object files (*.o). You should hand in the minimum set of files you need to compile plus the README file. Your “resources” directory should only contain your GLSL files.
• Create a single zip file of all the required files. The filename of this zip file should be USERNAME.zip (e.g., sueda.zip). The zip file should extract everything into a folder named USERNAME/ (e.g. sueda/).
• When you unzip this file, it should extract src/, resources/, CMakeLists.txt, and your README file to the USERNAME/ directory.
• Use the standard .zip format (not .gz, .7z, .rar, etc.).