Learned Motion Matching: Pitfalls and Implementation Notes

While implementing Learned Motion Matching, we ran into a number of problems that have to be solved before the result is usable. Below we walk through the methods we used to address them.

In the previous article, I reviewed the current methods for character animation based on machine learning. There is, however, a set of techniques that go unmentioned in that overview but without which the final result falls apart. This article shares what we learned implementing Learned Motion Matching.

Preface

All formulas assume a left-handed coordinate system:

Quaternions are stored in the $(x, y, z, w)$ format (scalar last).

Quaternion unroll

Because the same rotation can be represented by two quaternions $q$ and $-q$, the data can end up with sharp discontinuities (e.g. -179° and +181° encode the same rotation), and the same pose — or very similar poses — can be encoded in completely different ways. This makes the data meaningless for neural network training, and any attempt to interpolate across such opposing values produces garbage. To fix this, we have to unroll the quaternions:

def unroll(x):
    y = x.clone()
    for i in range(1, len(x)):
        d0 = torch.sum( y[i] * y[i-1], dim=-1)
        d1 = torch.sum(-y[i] * y[i-1], dim=-1)
        y[i][d0 < d1] = -y[i][d0 < d1]

    return y

The idea is simple: if the dot product with the previous frame's quaternion is negative, we flip the sign. This keeps the invariant $q_i \cdot q_{i-1} \ge 0$ across the entire sequence.

A visual comparison of quaternion component values before and after unrolling:

Once quaternions are unrolled, the next critical task is to create a new artificial bone that describes the motion trajectory — the root transformation.

Root transformation

This bone can be authored manually (e.g. in Blender) or generated procedurally, but it must exist — without it, the entire algorithm does not work. When building it, keep in mind that it has to capture the character of the motion as accurately as possible. After it is added, it becomes the root of the skeleton. Its position is computed by projecting the upper spine bone onto the ground plane and applying a Savitzky-Golay smoothing filter:

from scipy import ndimage, signal

g_pos, g_rot = quat_utils.fk(y_pos, y_rot, parents) # forward kinematics

sim_bone_pos = torch.from_numpy(
    signal.savgol_filter(g_pos[:, spine_node_id] * torch.tensor((1.0, 1.0, 0.0)),
                         31, 3, axis=0, mode='interp')
)

On a calm-walk dataset, this filter can produce a difference of about 0.6 meters between the original and the smoothed position — which is significant.

Next, we take the forward direction of the hip bone, project it onto the ground plane, and smooth it with the same filter. We then convert it into a rotation around the vertical axis — that gives us the bone's rotation. Again, we apply a Savitzky-Golay filter:

# direction taken from projected hip forward direction
# 0.0, 0.0, 1.0 is local up axis
sim_bone_dir = (torch.tensor((1.0, 1.0, 0.0)) *
                    quat_utils.mul_vec(g_rot[:, hips_node_id:hips_node_id+1],
                                       torch.tensor((0.0, 0.0, 1.0)).unsqueeze(0).unsqueeze(0))
                    )
sim_bone_dir = quat_utils.normalize(sim_bone_dir)
sim_bone_dir = torch.from_numpy(signal.savgol_filter(sim_bone_dir, 61, 3, axis=0, mode='interp'))
sim_bone_dir = quat_utils.normalize(sim_bone_dir)

sim_bone_rot = quat_utils.normalize(
        quat_utils.between(torch.tensor([1.0, 0.0, 0.0]).unsqueeze(0).unsqueeze(0), sim_bone_dir)
    ).squeeze(1)

Always remember to normalize quaternions after every operation.

After that, we make the new bone the root of the skeleton and recompute the positions of all other bones under the new hierarchy (and update the hierarchy itself):

# Transform first joint to be local to sim and append sim as root bone
y_pos[:, 0] = quat_utils.inv_mul_vec(sim_bone_rot, y_pos[:, 0] - sim_bone_pos)
y_rot[:, 0] = quat_utils.inv_mul(sim_bone_rot, y_rot[:, 0])

y_pos = torch.concat([sim_bone_pos.unsqueeze(1), y_pos], dim=1)
y_rot = torch.concat([sim_bone_rot.unsqueeze(1), y_rot], dim=1)

parents = torch.concat([torch.tensor((-1,)), parents + 1])

Once the new bone is in place, the next step is computing the velocities of every bone in the skeleton.

Bone velocities

Large spikes in the values cause animation jitter downstream, so even though velocity can be computed with the standard formula, it's better to use central difference interpolation to suppress the spikes: $v_i = (x_{i+1} - x_{i-1}) / 2h$

vel[1:-1] = (pos[2:] - pos[1:-1] +
            pos[1:-1] - pos[:-2]) * 0.5 * fps
vel[0] = vel[1] - (vel[3] - vel[2])
vel[-1] = vel[-2] + (vel[-2] - vel[-3])

The method uses a symmetric window around each point. This has three effects:

• noise and spikes in a single frame are compensated by neighboring frames;
• the estimate becomes second-order accurate $O(h^2)$, rather than first-order as with forward/backward difference;
• the resulting velocity is smoother and less sensitive to individual outliers;

Note that this does not eliminate outliers completely — it only reduces sensitivity to noise.

The first and last points have to be handled separately via linear extrapolation, since they obviously don't have neighbors on both sides.

With velocities computed, we can move on to representing bone rotations and detecting foot contact.

2-column rotation matrix

As mentioned earlier (see quaternion unroll), having a neural network predict quaternions directly is not an optimal choice, so bone rotations are fed into and produced by the network as pairs of relative forward and up vectors.

def to_xform_xz(x):
    qx, qy, qz, qw = x[..., 0:1], x[..., 1:2], x[..., 2:3], x[..., 3:4]

    x2, y2, z2 = qx + qx, qy + qy, qz + qz

    # Precompute products for the rotation matrix
    xx, yy, wx = qx * x2, qy * y2, qw * x2
    xy, yz, wy = qx * y2, qy * z2, qw * y2
    xz, zz, wz = qx * z2, qz * z2, qw * z2

    return torch.concatenate([
        torch.concatenate([1.0 - (yy + zz), xz + wy], dim=-1)[..., torch.newaxis, :],
        torch.concatenate([xy + wz, yz - wx], dim=-1)[..., torch.newaxis, :],
        torch.concatenate([xz - wy, 1.0 - (xx + yy)], dim=-1)[..., torch.newaxis, :],
    ], dim=-2)

Foot contact detection

Foot contact can be derived from velocity with a threshold — if the velocity is below a threshold, we consider the foot planted. The L2 norm of the 3D velocity collapses the value to a scalar. We add a median filter to eliminate spurious short windows (false positives):

contact_velocity_threshold = 0.15 * 100  # to cm
# L2 norm of 3D velocity collapses all motion into single scalar speed
contact_velocity = torch.norm(g_vel[:, [left_toe_id, right_toe_id]], dim=-1)

extra = contact_velocity < contact_velocity_threshold

for ci in range(extra.shape[1]):
    extra[:, ci] = torch.from_numpy(ndimage.median_filter(extra[:, ci], size=6, mode='nearest'))

An example of the filter in action:

Finally, the last data augmentation step before training the networks is computing the motion trajectory.

Motion trajectory

The trajectory is computed at 20, 40, and 60 frames into the future, so the end of the clip has to be handled correctly:

def clamp_frame(frame, offset):
    for i in range(len(range_starts)):
        if range_starts[i] <= frame < range_stops[i]:
            return min(max(frame + offset, range_starts[i].item()), range_stops[i].item() - 1)


for i in range(start, stop):
    future_frame_20 = clamp_frame(i, 20)
    future_frame_40 = clamp_frame(i, 40)
    future_frame_60 = clamp_frame(i, 60)

This way, as we approach the end of the clip, the future frame is smoothly clamped to the last available frame.

Trajectory positions and rotations could be computed as a simple difference, but if we additionally multiply by the inverse of the root transformation quaternion, the character is rotated so that it always faces the same direction at every frame. This gives us invariance to the direction of motion.

traj_pos_20 = quat_utils.inv_mul_vec(y_rot[i, 0], y_pos[future_frame_20, 0] - y_pos[i, 0])
traj_pos_40 = quat_utils.inv_mul_vec(y_rot[i, 0], y_pos[future_frame_40, 0] - y_pos[i, 0])
traj_pos_60 = quat_utils.inv_mul_vec(y_rot[i, 0], y_pos[future_frame_60, 0] - y_pos[i, 0])

# Rotate the local forward vector (1.0, 0.0, 0.0) because X is forward to get world-space forward
traj_dir_20 = quat_utils.inv_mul_vec(y_rot[i, 0], quat_utils.mul_vec(y_rot[future_frame_20, 0], torch.tensor((1.0, 0.0, 0.0))))
traj_dir_40 = quat_utils.inv_mul_vec(y_rot[i, 0], quat_utils.mul_vec(y_rot[future_frame_40, 0], torch.tensor((1.0, 0.0, 0.0))))
traj_dir_60 = quat_utils.inv_mul_vec(y_rot[i, 0], quat_utils.mul_vec(y_rot[future_frame_60, 0], torch.tensor((1.0, 0.0, 0.0))))

The same applies to computing positions and rotations of the bones that go into the pose feature vector:

left_foot_pos = quat_utils.inv_mul_vec(rot[:, 0], g_pos[:, left_foot_id] - g_pos[:, 0])
right_foot_pos = quat_utils.inv_mul_vec(rot[:, 0], g_pos[:, right_foot_id] - g_pos[:, 0])
left_foot_vel = quat_utils.inv_mul_vec(rot[:, 0], g_vel[:, left_foot_id])
right_foot_vel = quat_utils.inv_mul_vec(rot[:, 0], g_vel[:, right_foot_id])
hips_vel = quat_utils.inv_mul_vec(rot[:, 0], g_vel[:, hips_node_id])

That concludes the data augmentation pipeline.

Projector nearest neighbor search

As a reminder, the Projector network emulates a nearest-neighbor search for $x$. The problem reduces to computing Euclidean distances and is quite resource-intensive — it takes up most of the runtime of the entire algorithm.

$d(p,q) = \sqrt{(p_1 - q_1)^2 + (p_2 - q_2)^2 + \ldots + (p_n - q_n)^2} = \sqrt{\sum_{k=1}^n(p_k - q_k)^2 }$

To speed up the search, we can use the sklearn.neighbors.BallTree structure:

from sklearn.neighbors import BallTree


tree = BallTree(x)
nearest = tree.query(x_hat, k=1, return_distance=False)[:, 0]

Other solutions exist — faiss, for instance, gives an even bigger performance boost, but it's harder to use and lacks official Windows support.

Saving and integrating trained models

After the LMM models are trained, the next step is integrating them into the target runtime. It's important to pick a format that is both inference-friendly and compatible with the engine.

There are two options:

1. Save the models in ONNX format
2. Write the weights of each layer and the normalization parameters (std and mean) directly into a binary format

ONNX is convenient because, if the target environment provides inference tooling for this format, there is nothing else to do. In the second case, you additionally have to implement data normalization, matrix multiplication, and the activation function (ReLU in LMM) — and do it efficiently.

The only issue we hit with ONNX was compatibility between the ONNX opset version supported by the target environment and the opset_version used when exporting the model.

For the LMM integration target, we picked Unreal Engine 5, which has built-in ONNX inference tooling — UNNERuntime. To use it, you just load the model and describe the inputs and outputs.

The ONNX export itself is trivial — just remember to switch the model to eval mode and disable gradients:

compressor.eval()

with torch.no_grad():
    input = (torch.cat([
            y_pos[:, 1:].reshape([num_frames, -1]),
            y_txz[:, 1:].reshape([num_frames, -1]),
            y_vel[:, 1:].reshape([num_frames, -1]),
            y_ang[:, 1:].reshape([num_frames, -1]),
            q_pos[:, 1:].reshape([num_frames, -1]),
            q_txz[:, 1:].reshape([num_frames, -1]),
            q_vel[:, 1:].reshape([num_frames, -1]),
            q_ang[:, 1:].reshape([num_frames, -1]),
            y_root_vel.reshape([num_frames, -1]),
            y_root_ang.reshape([num_frames, -1]),
            y_extra.reshape([num_frames, -1]),
        ], dim=-1) - compressor_mean_in) / compressor_std_in

    z = compressor(input.to(device))

torch.save(z.detach().cpu(), root_dir / 'z.pt')

torch.onnx.export(
    compressor,
    (compressor_mean_in,),
    f"{root_dir}/compressor.onnx",
    input_names=["y_q"],
    output_names=["z"],
    dynamo=True
)

Other models are exported the same way, except you don't need to save any intermediate values.

Once the model is saved, it can be loaded into UE5 and run.

Below is an example of running the Decompressor model:

bool NNModels::RunDecompressorModel(TArray<float>& features, TArray<float>& latent) {
   if (DecompressorModel.IsValid() == false || InputBindings.Num() == 0) {
      return false;
   }

   TArray<float> input;
   input.SetNumZeroed(features.Num() + latent.Num());

   for (int i = 0; i < features.Num(); i++) {
      input[i] = features[i];
   }

   for (int i = 0; i < latent.Num(); i++) {
      input[features.Num() + i] = latent[i];
   }

   InputBindings[0].Data = input.GetData();
   InputBindings[0].SizeInBytes = input.Num() * sizeof(float);

   OutputBindings[0].Data = DecompressorOutputData.GetData();
   OutputBindings[0].SizeInBytes = DecompressorOutputData.Num() * sizeof(float);

   return DecompressorModel->RunSync(InputBindings, OutputBindings) == UE::NNE::EResultStatus::Ok;
}

Conclusion

With these techniques in place, you can build a reliably working Learned Motion Matching system. That said, none of this should be treated as the only correct approach — there are alternatives for interpolation, for constructing the root transformation, for data representation, and so on.