How to decompose 2D trajectory data into submovements using matlab

There is significant evidence that movements are made up of elementary building blocks (see Flash & Hochner, 2005 for a review). When looking at arm reaching movements, or drawing movements, these building blocks are often called "submovements". If we assume that these submovements have a particular form, then we can decompose an arbitrary movement trajectory into its constituent submovements.

Here, I will assume that the submovements are minimum jerk submovements, that is, they minimize mean squared jerk (other options are available). These submovements are described as 5th order polynomials with 3 free parameters (which is the minimum possible) - amplitude, duration, and starting time. The equation for the position (relative to the starting point) and the velocity are given by:

where D is the duration, A is the amplitude and t0 the starting time, and t0 < t < t0 + D. The Matlab function MJxy.m computes this velocity profile.

We can have more than one submovement being executed at a time, in which case we just add the change of positions of both submovements. Our goal is to take a trajectory and decompose it one or more submovements. To do this, we need to find the 3 parameters (or 4 if we are working in 2D, where we have an extra amplitude parameter) that describe each submovement (and decide how many submovements we need for a particular movement). To decide how good a reconstruction is, we need to define a cost, which is based on the difference between the observed trajectory and the reconstructed trajectory. Here will we used the squared difference (MSE). Minimising the MSE is equivalent to minimising the root mean squared error (RMSE), we use the MSE for computational reasons. Once we have the cost, we need to vary the parameters to find the minimum cost (i.e., best fit).

Rather than just minimising the difference between the observed and predicted velocity profiles, we also include a term for the tangential velocity, which we calculate for each submovement and then sum. This is to prevent the algorithm finding two submovements which start at approximately the same time, but have approximately equal and opposite amplitudes. While the resultant velocity will be approximately zero (and so won't affect the MSE), from a theoretical point of view, these seem highly unlikely. So our cost function E is:

where Fxi are Fyi are the predicted ith submovement velocities in the x and y directions, and Gx and Gy are the observed x and y velocities. The cost function, in Matlab, can be found in the file calculateerrorMJxy.m. Due to the potentially large number of parameters (4 per submovement), the optimisation can be very slow. In order to help Matlab along, I calculated the gradient (partial derivatives of the error function) and the Hessian (second-order partial derivatives), which otherwise Matlab has to estimate numerically. These expressions are long and complicated, so I used the program maxima to find them [the maxima code, and a perl script to turn the maxima output into legal matlab code, and some code to check that there are no mistakes in calculating the gradient and Hessian].

To minimise this function in Matlab, we can use fmincon (part of the Optimization toolbox). We use this rather than fminsearch so that we can set constraints on some of the variables, in particular:

• 0 <= t0 < tf - 0.167
• 0.167 < D < 1

These enforce that submovements need to have a duration of at least 167 ms. We also constrain the amplitudes to reasonable values. After running the optimization and getting back the parameters and the error, we need to decide how many submovements we want. I usually just use a threshold (0.02), and select the minimum number of submovements that give an error of less than 0.02. Before running the optimization, we need to select a starting "guess". Here, I randomly select a starting point within the legal ranges for each parameter. However, we may end up in a local minimum. To minimize the chance of this occurring, we start from 10 different starting points, as suggested in Rohrer & Hogan (2006).

The function decompose.m performs the decomposition (calls fmincon, etc). I have provided an example (example.m) to demonstrate how to use the code (which needs a sample data file: sampledata.mat). The output of the example program is shown below:

In this case, two submovements were needed to get the reconstruction error less than 0.02. In this case, this consisted of a small submovement to the right (in green, starting at 378 ms with an amplitude of [0.088,0.125] and duration of 345 ms), followed by another submovement at 551 ms (in red, with an amplitude of [-0.432,0.346], and a duration of 689ms). Note that these numbers which describe the submovements are taken from the best fit parameters.

All files needed to run the demo can be downloaded here in a zip file.