Identifiable Decomposition of Submovements in Human Hand Trajectories
Authors: Adrian Prados, James Hermus, Ramon Barber, Sylvain Calinon
Organizations: Universidad Carlos III de Madrid, Leganes, Spain · Idiap Research Institute, Martigny, Switzerland · Ecole Polytechnique F´ed´erale de Lausanne (EPFL), Lausanne, Switzerland
Voluntary movements have long been hypothesised to be comprised of discrete primitives called submovements, as a descriptive model of human motor behaviour. However, existing methods scale poorly, and no principled method exists to determine whether a decomposition is informative. We propose a spatiotemporal kernel correlation between primitive pairs as an identifiability criterion. Submovement-Identifiable Decomposition (Sub-ID) embeds this criterion in its adaptive-ridge regularisation, biasing the optimiser toward low-correlation solutions. Identifiability is lost when primitives become collinear and recovered when they diverge spatially. On synthetic data, Sub-ID recovers ground-truth parameter distributions where existing methods fail; furthermore, when primitives overlap too heavily to be distinguished, the method explicitly detects this ambiguity rather than outputting misleading results. Sub-ID extracts submovements from real three-dimensional, long-horizon movements, a regime no prior method addresses. This method has the potential to identify physiologically grounded primitives for motor control research and imitation learning.
Figures & tables
Figure 1: Sub-ID framework and identifiability criterion. (A) Submovement decomposition. The spatial trajectory is decoupled into scalar velocity to extract temporal velocities bases. These are spatially projected, integrated to position and scaled via Ridge Regression to prevent collinearity. A final global optimization minimizes positional error to yield the position reconstruction. (B) Spatiotemporal Kernel Correlation. Severe temporal overlap ( ρ→1 ) causes unresolvable aliasing. Sub-ID mitigates this by incorporating spatial alignment, ensuring highly overlapped primitives remain distinguishable as long as their spatial directions diverge.
Figure 2: Comparison of kinematic parameter recovery against the GT for ( K=3,8,15 ). Columns represent, from left to right: basis function durations, amplitudes, the number of extracted functions, and the refractory period. Results are shown for our method, Sub-ID (blue), SSSUMO (green), Scattershot (red) and Gowda (purple), against the original GT distribution (grey). The square symbol indicates that no valid solution was found by the algorithm for that specific scenario. Sub-ID consistently replicates the original distributions across all scenarios, whereas baseline methods exhibit significant parameter shifts, incorrect basis counts, or convergence failure as temporal density increases.
Figure 3: Evaluation of extracted kinematic parameters and kernel overlaps. (A) Overlap distributions for varying complexities ( K={3,8,15} ). Sub-ID consistently matches the GT low-overlap profiles. Gowda offers a resilient baseline but exhibits a moderate increase in overlapping bases at K=15 . SSSUMO exhibits increasing correlation, culminating in heavy overlap near 1.0 at K=15 . Scattershot suffers interference and fails to converge at K=15 . (B) Bimodal GT comparison ( K=15 ). Sub-ID robustly captures the underlying dual structure and maintains bounded overlaps, despite minor variance in the recovered basis count (14–16 vs. 15). (C) Unimodal GT extreme overlap scenario ( K=15 ). Severe temporal aliasing exposes identifiability limits: closely spaced submovements are merged, causing a significant underestimation of the basis count and a drastic increase in overlap. For (B) and (C) , blue presents our method and gray the Ground Truth (GT).
Figure 4: Qualitative evaluation of kinematic decompositions in 2/3D tasks. The figure presents spatial trajectory reconstructions alongside their tangential velocity profiles and extracted primitives. For the 2D tasks, the top row shows a simple stroke from the Letters dataset (’A’), and the bottom row a continuous sequence from the PushT dataset. For the 3D tasks, the top row displays a short trajectory from the Moving dataset, and the bottom row a continuous sequence from the PushT 3D dataset. Across both domains, Sub-ID (blue) fits the signal using a concise set of primitives. In contrast, SSSUMO (green) exhibits overfitting in all scenarios, fragmenting the movement into an excess of highly overlapped primitives, clearly observable in the PushT sequences. Scattershot (red) is included in the 2D evaluation—where it shows underfitting in long tasks by omitting secondary peaks and valleys—but is excluded from the 3D tasks due to its difficulty scaling to high-dimensional and complex sequences.
Overlaps >0.8 ( μ±σ )
Position RMSE (cm)
Velocity RMSE (m/s)
Co-lineality Φ
Spatiotemporal ρ
Method
K=3
K=8
K=15
K=3
K=8
K=15
K=3
K=8
K=15
K=3
K=8
K=15
K=3
K=8
K=15
GT
0.00 ± 0.00%
0.00 ± 0.00%
0.00 ± 0.00%
0.0
0.0
0.0
0.000
0.000
0.000
1.0
1.2
1.5
–
–
–
Sub-ID
0.84 ± 0.19%
2.08 ± 1.06%
3.00 ± 1.35%
0.2
0.4
1.2
0.005
0.007
0.009
1.5
2.8
8.2
0.32
0.35
0.45
Gowda
0.95 ± 0.25%
7.20 ± 2.15%
12.30 ± 3.20%
0.4
1.5
6.0
0.006
0.018
0.035
3.5
11.0
27.3
0.32
0.45
0.68
Scattershot
0.00 ± 0.00%
17.14 ± 6.98%
–
0.5
8.5
–
0.005
0.057
–
4.0
150.0
–
0.30
0.82
–
SSSUMO
4.67 ± 13.35%
15.46 ± 6.31%
51.83 ± 4.61%
1.5
4.5
12.0
0.012
0.073
1.673
12.0
85.0
450.0
0.72
0.88
0.96
Table 1: Performance comparison across synthetic tasks of varying complexity (top) and Comparative Evaluation on 2D and 3D Datasets (bottom).
Figure 5: Algorithmic pipeline of the Sub-ID framework. The method begins with Heuristic Peak-Based Detection , extracting the scalar tangential velocity from the human trajectory to identify the initial set of basis functions. The system then evaluates the residual error to trigger Greedy Residual Refinement , where new candidate primitives ( ϕvel(t;θtry) ) are iteratively generated, evaluated, and injected to minimize the remaining velocity error. Finally, the Adaptative-Ridge Optimization stage translates these temporal bases into the spatial domain via numerical integration ( Φpos,k(t) ). A Ridge Regression approach is applied to robustly compute the spatial amplitudes across all dimensions, ensuring numerical stability, preventing spatial amplitude errors, and finalizing the optimal spatio-temporal reconstruction.
Figure 6: (A) . Temporal initialization and greedy residual refinement. (Left) Heuristic Peak-Based Initialization algorithm successfully identifies the most dominant kinematic submovements (light blue), providing a coarse initial fit. Because these primary peaks are insufficient to fully capture the complex movement, a significant residual error remains. (Center) The greedy refinement process iteratively locates the maximum residual error and progressively introduces new basis functions (red) to explain the missing kinematic intensity. In each iteration, the amplitudes of all previously discovered bases (blue) are simultaneously re-optimized. (Right) The final converged solution accurately reconstructs the scalar tangential velocity profile. This refined temporal scaffolding serves as the robust starting point for the subsequent global optimization in the position domain. (B) Impact of Ridge Regression on basis cancellation. Comparison of the spatio-temporal global optimization of a 3D trajectory with L2 regularization (Left) and without (Right). Both approaches reconstruct the 3D spatial position and tangential velocity, but the unregularized solution exhibits severe multicollinearity: multiple nested primitives emerge simultaneously with extreme opposite signs (visible in x˙2 and x˙3 , marked in dashed red lines), nearly cancelling one another to fit small residual variations. Ridge Regression (Left) penalizes these extreme amplitudes, yielding a stable, sparse decomposition without altering the total number of basis functions.
Figure 7: Unimodal and bimodal synthetic data. The central panels depict the temporal organization of 15 overlapping basis functions generated. The right panels show the empirical distributions of duration, amplitude, and refractory period, parameterized from human demonstrations collected in extended 2D and 3D tasks. In the bimodal setting, these distributions exhibit a clear separation into two distinct modes, which is consistently reflected in the structure of the reconstructed tangential velocity profiles and the corresponding basis function activations.
Appendix figures & tables7 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 1: Example of data generated and results obtained by our method for the case of K=3 (top image), K=8 (middle image), and K=15 (bottom image) for unimodal distributions. It can be seen that the method is capable of generating solutions that produce basis functions practically identical to those used as ground truth to generate the synthetic data. An image of the 2D trajectory for the case K=3 is also included, as well as the decomposition into the x˙1 and x˙2 velocity components, where the solutions generated by our algorithm for the synthetic-data case can be observed.
Figure 2: Trajectory fitting pipeline via LGNB with asymmetry optimization. The shaded areas represent the contribution of each optimized LGNB base. The search space for the μ parameter enables the generation of asymmetric solutions. This geometric flexibility is crucial for accurately approximating real velocity peaks skewed toward temporal boundaries, achieving a robust fit against non-symmetric dynamics and drastically reducing overfitting thanks to the coupled Ridge regularization.
Figure 3: Comparative spatiotemporal decomposition using MinJerk and LGNB basis functions. The pipeline demonstrates its agnosticism to the underlying kinematic model by successfully decomposing the same 3D trajectory using Minimum Jerk (Left) and Bounded Lognormal (Right) profiles. Both models achieve a highly accurate 3D spatial fit (bottom left of each panel) with an approximate RMSE of 2×10−4 . However, the choice of basis introduces differences in optimization complexity and regularization needs. The symmetric MinJerk model (Left) resolves the trajectory using K=9 primitives and a Ridge penalty of α=0.1 in 2.7 seconds. The more complex LGNB model (Right), which optimizes additional parameters like skewness ( μ ), requires K=10 overlapping primitives and a highly controlled, smaller regularization penalty ( α=5×10−5 ) to manage multicollinearity. This increased dimensionality extends the computation time to 4.0 seconds but allows for the capture of asymmetric velocity profiles.
Figure 4: Experimental setup for the PushT dataset collection. (A) Overview of the complete experimental setup, including the Motion Capture (MOCAP) system used to record the tasks performed by the users. The image also shows the block to be pushed and the silhouette corresponding to its target final pose. (B) Participant performing the 2D and 3D PushT data collection tasks.
Metric
Reference
Sub-ID
Sync Error ( Esync )
151.4 ms
119.3 ms
Δtmin ( Tmed , r=0.95 )
—
0.0352 s
Δtmin ( Tmed , r=0.90 )
—
0.0492 s
Max Spatiotemporal ρ
—
0.6595
Mean Overlap ( μ±σ )
0.201 ± 0.021
0.171 ± 0.036
Overlaps >0.8 ( μ±σ )
0.00% ± 0.00%
2.90% ± 1.64%
Appendix
Table 1: Kinematic and Mathematical Identifiability Metrics for Bimodal Synthetic Data
Metric
Reference
Sub-ID
Sync Error ( Esync )
297.4 ms
212.3 ms
Δtmin ( Tmed , r=0.95 )
—
0.0449 s
Δtmin ( Tmed , r=0.90 )
—
0.0629 s
Max Spatiotemporal ρ
—
0.4982
Mean Overlap ( μ±σ )
0.251 ± 0.031
0.534 ± 0.029
Overlaps >0.8 ( μ±σ )
15.00% ± 5.76%
33.90% ± 8.08%
Appendix
Table 2: Kinematic and Mathematical Identifiability Metrics for Unimodal Failure Synthetic Data
Metric
Equation
Reference
Pos Error (RMSE)
RMSEpos=N1∑n=1NP(tn)−P^(tn)2
-
Vel R2
Rvel2=1−∑(v(tn)−vˉ)2∑(v(tn)−v^(tn))2
-
Vel RMSE
RMSEvel=N1∑n=1N(v(tn)−v^(tn))2
-
Primitive Rate
Rate=TtotalK
-
Temporal Overlap
Oi,j=min(Ti,Tj)intersect(ti,tj)
[ 50 ]
Spatiotemporal Kernel Correl.
ρi,j=cos(ψi,j)⋅ρi,jtemp
[ 50 , 12 ]
Appendix
Table 3: Summary of Statistical Analysis and Evaluation Metrics. Mathematical formulations used to assess kinematic fidelity and structural identifiability.
School of Future Technology, South China University of Technology · Shenzhen Key Laboratory of Ubiquitous Data Enabling, Tsinghua Shenzhen International Graduate School, Tsinghua University