2.1 Simulated 2D throwing examples

To illustrate how our decomposition allows to identify the key features of the action distribution contributing to performance, we consider a toy model of a throwing task in which a projectile is released from a fixed position with given horizontal and vertical velocity components and must hit a target at a fixed height on a vertical plane at a fixed distance. Additionally, we assume that the horizontal component of the release velocity lies in the vertical plane containing the release position and the target, so that the action a is a two-dimensional (2D) vector, the projectile trajectory lies on a vertical plane, and the outcome (x, i.e., the arrival height) is a scalar (Fig 3A). The score π(x;x_T) is the squared distance of the arrival position from the target.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Examples of different throwing strategies in a simulated 2D throwing task.

(A) Toy model of a throwing task in which the action (a) is characterized by only two parameters, i.e. the horizontal (v_y0) and the vertical (v_z0) release velocity components, the outcome (x) is the arrival position of the projectile on a vertical plane at a distance d from the release position, and the score (π) is the squared distance of the arrival position from a target (x_T). (B) Action score (gray-shaded background), solution manifold (yellow line), and Hessian (red ellipses) in a region of release parameters. The score sensitivity for different throwing strategies is captured by the Hessian, which is represented with red ellipses whose major axis indicates, locally, the direction of maximum sensitivity (smaller tolerance to noise). (C) Five (simulated) individual strategies. For each strategy (S1-S5, different panels and colors) the individual throwing actions are shown together with mean (large black circle), their covariance (black ellipse, 95% C.L.) and the Hessian at the mean action (red ellipse). (D) Hessian-based decomposition of the mean score of the five strategies. For each strategy, the mean release parameters (v_y0 and v_z0), the mean score (), and the indicators from the decomposition (score of the mean action α, tolerance-variability index β, tolerance (τ), noise η, and alignment θ) are presented in the table. The four scatter plots illustrate how these indicators allow to differentiate strategies with different performance as well as strategies with same performance. See text for more details.

https://doi.org/10.1371/journal.pone.0253626.g003

As in Fig 2B, since the actions are 2D, in Fig 3B we can represent the action-to-score function as contours in the action plane (i.e., the plane v_y0-v_z0 of the horizontal and vertical release velocities) with gray shading (lighter shades for lower error) and the solution manifold, i.e. the set of optimal actions which results in 0 penalty score, with a yellow curve. Fig 3B also shows the local sensitivity of the score, i.e. the key element of our decomposition, which is captured by the Hessian of the action-to-score function. For different values of the mean action (i.e. and ), the sensitivity of the mean score is represented by a red ellipse. The longer axis of the ellipse, whose length is proportional to the largest eigenvalues of the Hessian matrix (see Methods), corresponds to the direction in which a given deviation of an action from the mean affects the mean score the most, i.e. the direction of maximal sensitivity or minimal tolerance to noise. The length of each axes of the ellipse indicates how much a deviation in that direction affects the mean score. Consequently, action distributions with the same amount of noise will have different performances depending on the alignment of their covariance with the direction of maximal sensitivity and on the lengths of the axes of the local sensitivity ellipses. Note that the orientation of the ellipses, their size, and the ratio of their axes depend on the position on the action plane. On the solution manifold, the ellipses have only one axis with non-zero length, indicating that there is only one direction in action space which is relevant for the score. This direction however, is not constant, but changes along the solution manifold. Also note that the ellipses are smaller for larger horizontal velocities , indicating that the score is more tolerant for faster throws.

2.1.2 Comparison with TNC, UCM, and GEM analyses.

The 2D throwing toy model also allows to compare the Hessian-based decomposition with existing computational methods that have been introduced for analyzing the relationship between an action distribution and its mean score (TNC), for decomposing action variability into task-relevant and task-irrelevant components (UCM), and for relating action variability and outcome variability (GEM).

The TNC-Cost analysis [20] characterizes an action distribution with three indicators (or costs) computed by comparing the mean score of the distribution with the score of a transformed distribution, optimal according to different criteria. The tolerance cost (T-cost) represents the difference between the mean score of the distribution and the mean score of the same distribution translated to an optimal location, i.e. a location with the lowest possible score. For the noise cost (N-cost) the optimal distribution is obtained by shrinking the distribution (i.e., reducing the noise η) without affecting the mean. The covariation cost (C-cost) is computed using an optimal distribution with the same mean and the same individual action components (e.g. and in our 2D throwing example) but with an optimal re-association of those components in different actions, i.e. with an optimal covariation. These indicators are related to the indicators of our decomposition (see Section A.2 in S1 Appendix) but they represent a re-combinations of those indicators and, more importantly, do not provide a direct decomposition of the mean score, thus they do not allow to describe the performance trade-offs between indicators presented above. Also, critically, as the computation of both the N-Cost and the C-cost relies on exhaustive searches, they require a number of operations that scales exponentially with the dimensions of the action space and they are not feasible for 6D actions, as in unconstrained throwing.

Fig 4A presents an example of three simulated strategies with the same performance and an α-β trade-off. Increasing the mean horizontal release velocity , a larger bias α is compensated by a lower tolerance-variability index β achieved by a lower alignment θ and higher tolerance τ with constant noise η. The T-Cost is constant, as the three distributions have the same mean score and the same covariance matrix, i.e. they differ only for the position of the mean action, and therefore they have the same difference with the mean score of the optimally translated distribution. Thus, the T-Cost does not distinguish iso-performing strategies with different bias when this is compensated by better alignment and tolerance. Moreover, the strategies have different N-cost, even if their action distribution have the same covariance matrix. This is due to the fact that when the action-to-score function is smooth and the dispersion of the actions not too high, as in this case, the optimal distribution used to compute the N-cost has zero noise, i.e. is concentrated on the mean action, and the N-cost corresponds to .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Comparison of the Hessian-based decomposition with existing computational methods for the analysis of throwing strategies in a 2D throwing task.

(A) Comparison with TNC-Cost analysis. Three strategies (S1 blue, S2 green, S3 magenta) are simulated by generating a distribution of 50 release parameters (v_y0, v_z0) by randomly sampling from a bi-dimensional Gaussian distribution with a given mean and covariance (see parameters in the table on the left). The strategies have the same performance because they trade-off bias (α) with tolerance-variability (β), as can be clearly observed in the α-β scatter plot. Moreover, the strategies have equal noise (η). However, the strategies have equal T-cost but different N-Cost, which then represent a recombination of the indicators extracted from the Hessian-based decomposition. Moreover, as the TNC costs do not sum up to the performance, the α-β tradeoff cannot be observed. (B) Comparison with UCM analysis. Three simulated strategies with no bias (i.e. with a mean action on the solution manifold) differ in performance because of differences in β, due to differences in the alignment of the covariance (black ellipses) with the Hessian (red degenerate ellipses), since the noise is equal. The decomposition of the action variability onto an uncontrolled manifold (UCM) subspace and an orthogonal (ORT) subspace does not discriminate the three strategies. (C) Comparison with GEM analysis. Three simulated strategies with equal noise but different performance due to differences in α are not discriminated by the goal-relevant sensitivity Σ_GR and the goal-relevant variability fraction Φ_GR.

https://doi.org/10.1371/journal.pone.0253626.g004

In contrast to both Hessian-based decomposition and TNC analysis, UCM analysis [15, 16] aims at decomposing action variability rather than characterizing the dependence of mean score on strategy and variability. It decomposes variability into a component within an uncontrolled manifold (UCM) subspace, which does not affect the task because such subspace is the null space of the Jacobian of the action-to-outcome function, and a component within the orthogonal (ORT) subspace. The ratio , where σ_UCM is the square root of the variance in the UCM and σ_ORT in the ORT, then indicates whether a controller is successful in pushing variability onto task-irrelevant dimensions, thus improving performance. However, as shown by our decomposition, the relationship between action variability and performance, quantified by mean score, depends on the alignment of the action covariance with the Hessian and on the size of the Hessian rather than on the alignment with the Jacobian. Thus, a decomposition of variability onto the UCM and ORT subspaces does not discriminate strategies with different performance but equal σ_UCM and σ_ORT (and therefore equal ratio).

Fig 4B shows an example of three simulated strategies with different mean horizontal release velocity and different performances, due to different β, but equal σ_UCM and σ_ORT. The three action distributions have means on the solution manifold, i.e. α = 0 and equal noise η. They also have the same orientation (with major and minor axes parallel to the and axes respectively) of the covariance ellipses but different ratios between the axes of the ellipses. As the Jacobian as well as the associated UCM and ORT directions change with the mean release velocity, the different amounts of variability along the horizontal and vertical velocities have the same amount of variability projected along the UCM and ORT directions. In contrast, the changes in the alignment of the covariance with the Hessian, which also changes with the mean release velocity, as well the changes in the size of the Hessian, do capture the features of the action distribution that affect performance. In particular, the differences in mean score, are due to different ratios of θ and τ as shown by the iso-lines in the β − η plane.

Finally, the GEM analysis [17] extends the UCM analysis by considering the noise sensitivity of the action-to-outcome function in addition to the directions of the Jacobian and its null space. The total body-goal sensitivity at the mean action, i.e., the ratio , where σ_e is the square root of the total variation of the outcome ( in our notation) and σ_u is the square root of the total variation of the action (), gives the amplification of variability between the action and outcome levels and can be decomposed as (9) where the goal-relevant sensitivity , with σ_R = σ_ORT, gives the amplification of variability in the ORT subspace and the goal-relevant variability fraction is the fraction of action variability available for mapping into the outcome and hence that can have an effect on performance. However, as for the UCM analysis, the GEM analysis characterizes the effect of the first order approximation of the action-to-outcome function, i.e., the Jacobian, on the variability in outcome space rather than the effect of the interplay between the geometry of the action-to-score function and the action distribution (both mean and covariance) on the mean score, i.e. the performance. Indeed, the GEM does not capture the effect on performance of a non-optimal mean in the distribution, which gives rise to the α-β trade-off presented above.

Fig 4C shows an example of three simulated strategies with the same horizontal release velocity and the same noise but different mean vertical release velocities and, thus, different α and different performances, as clearly shown in the α-β scatter plot. In contrast, the GEM decomposition does not discriminate between the three strategies, as they all have nearly equal Φ_GR and Σ_GR. This is due to the fact that covariances are also nearly equal and the Jacobian does not change much when the mean vertical release velocity deviates from the optimal value on the solution manifold. However, such deviation of the mean action has a substantial effect on performance.

2.2 Unconstrained 6D throwing experimental data

To gain insights on individual strategies and their relationship to performance in a complex real-world motor skill we applied our Hessian-based decomposition to unconstrained throwing. We instructed 20 participants to perform unconstrained overarm throws to 4 circular targets arranged on a vertical target board placed at 6 m from the initial position (see Methods). We considered the release position and velocities (Fig 1A) and the squared distance of the arrival position from the center of the target as the score (loss) of each trial. As we found that performance varied significantly across participants (Fig 1B, p < 0.001 linear mixed-effect model with participant and target identity as categorical predictors), our aim was to characterize the key features of the individual action distributions that relate to performance and to identify features that differentiate participants without affecting performance.

As actions are six-dimensional (6D), it is impossible to visualize both the score and the individual strategies in a single plot, as in the 2D throwing case. To illustrate some of the characteristics of the individual action distributions, Fig 5 shows the distribution of the three pairs of release velocity components in five exemplar participants (columns) aiming at target T1. As there are 15 different pairs of 6 action variables (3 position and 3 velocity components) these plots do not fully illustrate the individual strategies (see Section B in S1 Appendix for a display of 9 of the 15 possible pairs). However, these examples show that the action distributions and the local geometry of the score function differ significantly across participants even when we consider a subset of action variables.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Examples of distribution of release velocity components for five representative participants (target T1).

Each row illustrates a pairs of release velocity components (, , ). Colored circles represent release parameters of individual throws. Black circles and ellipses represent mean and covariance (95% c.l.) of each parameter distribution. The gray-level shading of the contours indicates the local score, as a function of the release parameters. Note that in each row, the range of horizontal and vertical axes are the same for all participants, and the different shapes of the score across participant, reflect individual differences in the average action (mean position and velocity vectors at release). The wider the white area around the mean action, the more tolerant the score is to noise. Participants have been sorted from left to right according to their average release speed.

https://doi.org/10.1371/journal.pone.0253626.g005

As for the simulated 2D strategies illustrated in Fig 3C, participants have been sorted from left to right according to their mean longitudinal release velocity (), hence P10 is the slowest thrower and P11 the fastest. The domain of each velocity variable corresponds to the population mean ± 3 standard deviations. The score associated to each pair of velocity components, with the remaining velocity component and the three position components fixed at the subject-specific mean value, is indicated by the gray-scale shading of the contours. Note that the regions corresponding to actions with the same score (i.e. same gray shading) have different shapes (or geometry) across planes and participants, as they depend on the individual mean release action. The wider the white areas around the mean release parameter (black large circle), the more tolerant is the action-to-score function to action variability. In the plane (third row), as in the 2D case, tolerance increases with mean . However, the mean longitudinal release velocity also affects the tolerance along the horizontal transverse release velocity (), as can be observed in the (first row) and (second row) planes.

The distributions of the release parameters (colored circles), summarized in each plot of Fig 5 in terms of mean (black large circle) and covariance (two-standard deviations, black ellipse), also differs across planes and participants. In addition to differences in mean longitudinal velocity (), there are differences in mean vertical velocity (). In particular, P15 (second column), on average, releases the ball with a faster vertical velocity than P1, who has a similar longitudinal velocity. Remarkably, individual covariances differ both in size and orientation. In the (third row) plane, P10 shows a small covariance ellipse with the major axis oriented along the axis. In contrast, P11 shows a much larger covariance ellipse with the major axis oriented diagonally, i.e. with a negative correlation between the two velocity components.

In sum, in all three pairs of velocity variables (as well as in the other pairs of variable shown in Section B in S1 Appendix), both the local geometry of the action-to-score function and the structure of action variability suggest that specific features of the individual action distributions affect performance, as they determine the amount of overlap of the distribution with the regions with lowest penalty. However, as the distribution is 6-dimensional and there are 15 different pairs of action variables, it is not possible to identify by visual inspection a unique source of the inter-individual differences in throwing strategies and to systematically explain the relationship between action distribution and throwing performance. These limitations can be overcome by the Hessian-based decomposition.

2.2.1 Validity of the Hessian-based decomposition.

The Hessian-based performance analysis is based on the assumption that the individual action distribution is sufficiently localized around the mean action such that higher order terms of the Taylor expansion can be neglected in the action score approximation. To validate this assumption with our experimental data, we tested whether Eq (8) can be considered as an acceptable model of the mean action score. Fig 6 shows the relationship between the mean squared error and the sum α + β across participants and targets. We found that the sum α + β could explain 99% of the variance for targets T2, T3, and T4 and 97% for target T1, due to the larger error observed in P9. P9 had in fact the worst performance (see Fig 1B) and the largest bias (α) and variability (η), which violates the assumption of the action distribution being sufficiently localized. In sum, α + β can explain quite accurately the individual performance across participants and conditions, except for P9, where the sum α + β tends to overestimate the average score for T1, T2 and T3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Validity of the quadratic approximation.

Mean score vs. quadratic approximation (8) across targets (T1-T4, different panels) and participants (P1-P20) as in Fig 1. Note that P9 is the only participant deviating substantially from the identity line (dashed line).

https://doi.org/10.1371/journal.pone.0253626.g006

2.2.2 Hessian-based decomposition of the mean score reveals key performance-related features of individual throwing strategies.

The Hessian-based decomposition of the mean score relies on the straightforward computation of the action mean, action covariance and Hessian of the action-to-score function at the mean action. While the analysis of the structure of the Hessian and the action covariance matrices provides useful information on the task geometry and on the features of the action variability shared across subjects (see Section C in S1 Appendix), the parameters extracted from the Hessian-based decomposition of the mean score compactly describe the key performance-related features of the individual action distribution and allow to fully characterize how performance depends on such features.

Fig 7 shows the distributions of all the parameters of the Hessian-based decomposition of the mean score (α, β, η, τ and θ defined in Eq 8) across participants and relative to target T1. To visualize the inter-individual differences in terms of these five parameters, as for 2D throwing, we consider four scatter plots. The five exemplar participants of Fig 5 are identified by the color of the edges of the circular markers and show 6D throwing strategies that are similar to those illustrated in the 2D throwing examples of Fig 3, with matching colors.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Decomposition of the mean score across participants for target T1.

Circular markers are filled with areas that are gray-shaded according to the score (lighter gray for lower mean score). The edge of the markers corresponding to the five participants of Fig 5 are colored with matching colors. To increase the readability of the figure, P9 is not shown. (A) The α − β plane and the iso-performance lines. The dashed line represent the direction of maximal change in performance (orthogonal to the iso-performance lines). (B) The β − η plane. The lines have slopes corresponding to different values of . Participants along the same line, such as P1, P2, and P12, have the similar alignment-to-tolerance ratio and the difference in their β are only due to the noise η. (C) The θ − τ (tolerance-alignment) plane. The lines have slopes corresponding to different values of . The negative correlation between θ and τ indicates that participants with mean release action in a region with higher tolerance also tend to have lower alignment with the most noise sensitive directions. (D) The log(η)-log(θ/τ) plane and iso-β lines. As β is the product of θ/τ and η and log(β) = log(θ/τ) + log(η), in the log(η)-log(θ/τ) plane participants with the same β lie on a line with a −1 slope and maximal change of the tolerance-variability index occurs along the orthogonal direction.

https://doi.org/10.1371/journal.pone.0253626.g007

The α − β plane (Fig 7A) characterizes the performance of each participant and the trade-off between bias and tolerance-variability. The diagonal continuous lines with slope −1, corresponding to equal sum α + β, indicate iso-performing strategies. Participants closer to the origin of the plane are the best performers (e.g., P18 and P4, compare with Fig 1B). Moving away from the origin in the direction of the dashed line performance decreases. Participants with the same performance, i.e. on the same diagonal line, differ on their position on the line, thus showing a different trade-off between bias (α) and tolerance-variability (β). For example, P10, P1, and P11 have low α and P15 high α but they are all close to the same iso-performance line and thus have similar performance. In sum, the α − β plane shows that our participants differ in their performance because of their specific values of the α and β indicators and that participants with the same performance trade-off α with β.

The η-β plane (Fig 7B) characterizes the role of the amount of action variability, i.e. the noise η, in determining the tolerance-variability index β. The same contribution to performance due to tolerance-variability may be obtained with very different amounts of variability. For example, P10, P1, and P11 have similar β values but very different η values. This is because larger variability is compensated by lower alignment and higher tolerance. In fact, the slope of the lines in the η-β plane corresponds to and participants on lines with lower slope have a lower alignment or higher tolerance.

The interplay between tolerance and alignment is illustrated in the τ-θ plane (Fig in Fig 7C). Here the slope of the lines corresponds to a constant , i.e. each line corresponds to a line in the η-β plane (as ). While the relationship between τ and θ is not fixed, participants with a higher tolerance also have a lower alignment. This indicates that, when the mean action is located in a more tolerant region of the action space (e.g. higher horizontal velocity, as for P11) it is also possible to further optimize performance by reducing the alignment.

Finally, the log(η)-log(θ/τ) plane (Fig 7D) characterizes the trade-off between the amount of action variability and the other features of the action distribution determining the tolerance-variability index. In the log-log plane the multiplicative contribution of (determined by the mean action and the alignment of the action covariance) and of the amount of action variability η to β are represented by a sum () and, thus, strategies with equal β are located on lines with slope −1. Then, this plot clearly shows how participants differs in their tolerance-variability index β (participants with lower β such as P18 and P15 are on lines in the bottom-left) and how participants with the same β trade-off noise with tolerance and alignment (e.g. P10 vs. P1 vs. P11).

2.2.3 Consistency of the decomposition of individual action strategies across targets.

To test whether the individual performance-related features of action distributions described by the Hessian-based decomposition are consistent, we compared the α, β, η, τ and θ parameters across both participants and targets. If these parameters are robust indicators of the inter-individual differences in throwing strategy, we expect them to vary significantly across participants but not across targets.

Fig 8 shows the distributions across participants and targets of the parameters of the Hessian-based decomposition in the same combinations of parameters as in Fig 7, which included only target T1. The quadrilaterals (dashed edge lines and color-shaded areas), with the parameters for the four different targets of each participant as vertexes, are colored according to the mean score rank of each subject. In most cases the parameters do not show large variations across targets, i.e. they represent consistent features that characterize individual strategies. In the α-β plane (Fig 8A), a few participants show large across-target variations (e.g. P14, who has a larger α for T2 than for the other targets; P7, who has a larger β for T4; P20, who has both α and β larger for T2). However, most participants occupy small and distinct regions of the η-β (Fig 8B) and log(η)-log(θ/τ) (Fig 8D) planes. Finally, in the τ-θ plane (Fig 8C) about half of the participants occupy the region with intermediate values of τ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Decomposition of the mean score across participants for all targets.

In each panel, the quadrilateral connects the points corresponding to the values of a pair of indicators for the four targets of each participant (excluding P9) and they are colored according to the mean score of that participant (color scale on the top right). The targets are indicated by different colors of the circular markers and the participant number is indicated close to the marker for T1 (except for P7 and P14). (A) The α − β plane and the iso-performance lines. (B) The β − η plane. (C) The θ − τ (tolerance-alignment) plane. (D) The log(θ/τ)-log(η) plane.

https://doi.org/10.1371/journal.pone.0253626.g008

A linear mixed-effect model with with participant and target identities as categorical predictor (Table 1) of all five parameters supports these qualitative observations, showing that the effect of participant is significant for all parameters while the effect of target is significant only for α (p = 0.02) and τ. Thus, α and τ are the only parameters that are not consistent across targets. This is not surprising for τ since it depends on the individual strategy only through the mean release position and velocity, which, in turn, depend on the target position, and, thus, τ varies with the target. Concerning α, the contrast between pairs of targets (fixed effect coefficients) shows that its value is significantly higher for T2 than T1 (Δα = −0.02, p < 10⁻²) and T3 (Δα = −0.01, p = 0.01), indicating that the bias of a throwing strategies also depends slightly on the target.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Robustness of decomposition parameters across targets.

https://doi.org/10.1371/journal.pone.0253626.t001

3.1 Comparison with related approaches

In redundant, high-dimensional, and noisy tasks it is not enough to characterize the mean and the covariance of the action distribution to fully capture the relationship between action variability and performance. In agreement with earlier computational approaches addressing variability in multivariate actions [7, 15–17], our method highlights the key role of the geometry of the action-to-score function (captured by the Hessian) to assess how action variability affects performance. Differently from first-order methods such as UCM, GEM, or the more recent approach in [24], which characterize the local geometry with a linear approximation (expressed through the Jacobian matrix or the gradient vector), our method relies on a second-order approximation (based on the Hessian matrix). The main reason for which our method does not depend on the first-order term of the Taylor expansion of the action-to-score function and requires a second-order approximation is the fact that we are considering the mean score rather than the variability in action or outcome space as a measure of performance. As indicated in Eqs (14) and (15), the mean score does not depend on the gradient of the action-to-score function computed at the mean action, the reason being that the first-order term of the expansion is multiplied by the mean deviation from the mean action, which is null by definition. In other terms, changes in score (with respect to the mean) associated to actions that deviate from the mean action sum up to zero in the linear approximation of the action-to-score function. Indeed, for a linear action-to-score mapping the mean score is given simply by the score of the mean action, as all higher order derivatives in the expansion are null. Thus, in a quadratic approximation, it is only the local curvature of the action-to-score function, captured by the Hessian matrix, that affects the mean score.

How action distribution affects the mean score in a goal-directed behavior has been addressed by the TNC methods [7, 11, 13, 20]. The methods have been developed for, and applied to, a two-dimensional throwing task inspired by the skittle game, in which participants have to hit a target by releasing through a rotating joint (i.e., the action parameters are the release angle and tangential velocity) a virtual ball that could rotate in a plane around a pole. The TNC methods take into account the geometry of the action-to-score mapping implicitly, by evaluating the effects on performance of different action distributions through surrogate data (when comparing pairs of distributions [7]) or by exhaustive search (when comparing with optimal distributions [20]). In contrast, our method explicitly decomposes the contribution of different features of the action distribution through a Taylor expansion. Such analytical approach overcomes the disadvantage of the TNC methods concerning the use of numerical procedures necessary to generate surrogate data or to search for optimal distributions, which limits its applicability to high-dimensional actions. Moreover, our method allows to determine the contribution of the local geometry and of the action variability independently on each other. As illustrated in the 2D throwing examples of Fig 4A and detailed in the Appendix (see Section A in S1 Appendix), under the assumption of a smooth action-to-score function, for which the Hessian matrix is well defined, the tolerance, noise, and covariation terms of the TNC decomposition correspond to specific combinations of the terms in our decomposition. However, importantly, all the three terms of the TNC decomposition depend on both the Hessian and the action covariance, thus they do not indepentently characterize the contributions of the local action-to-score geometry, of the total action variability, and of the alignment between action variability and score tolerance τ, η, and θ. Moreover, when assessing a single distribution, the TNC-Cost analysis does not provide a unique decomposition of performance, as the three costs do not sum up to the mean score. Then, it is not possible to directly characterize the trade-off between the different components of the variability and local task geometry as with the Hessian-based decomposition. As an advantage, however, the TNC method does not rely on any assumption on the action-to-score function, such as smoothness and adequateness of a second-order approximation.

A fundamental difference between GEM and both the TNC and our approach, is the choice of performance measure. In the study of goal-directed motor skills performance is typically associated with a measure of error with respect to a predefined goal. Reduction of task errors and their variability is broadly recognized as an indicator of skilled performance. However, more recent views of human motor control, based on decision-making theory, propose an alternative definition of performance which is based on the concept of a score/loss function that essentially assigns a number (or a cost) to a given task error. In this perspective, the mean or expected loss is taken as measure of performance over repeated trials of a given motor task. Central to the GEM analysis is a goal function e(a) = x_T − f(a) that expresses the error between a desired goal/target x_T and the outcome f(a) of a given action a. By linearizing this error/goal function around the mean action of a strategy, GEM quantifies the overall contribution of tolerance, noise and alignment with the task/goal-level variability, i.e var(e). The Hessian-based decomposition proposed in this work extends the GEM analysis to decomposition of the mean of a performance indicator. In particular, we have shown that is the Hessian and not the Jacobian that affects the expected score/loss around the mean action of a strategy (see Eq (30)).

3.2 Assumptions and limitations

Our decomposition method requires a smooth action-to-score function and it provides an accurate estimate of the mean score only if the non-linearities in such function are adequately approximated by the second-order term of the Taylor expansion over the domain spanned by the actions. The assumption of smoothness (or at least continuity of the function and all partial derivatives up to the second order) is valid for a broad class of score functions, such as most penalty or reward functions usually employed to quantify task performance. For goal-directed tasks, requiring to minimize the distance from a target, the squared distance is a good choice because it leads to an action-to-score function which is twice differentiable everywhere the action-to-outcome function is smooth. The squared distance is preferable over the Euclidean distance because the latter has a singularity in the second derivative at zero, i.e. on the solution manifold. However, if the subject is attempting to minimize (maximize) the score, as the distance and the squared distance have the same minimum (maximum), both functions capture the control strategy equally well.

Another key assumption in our approach is that the second-order Taylor expansion of the action-to-score function around the mean action provides an acceptable approximation. As shown in Fig 6, for almost all participants and targets, the estimation of the mean score based on such approximation (α + β) is close to the actual mean score (E[π]). The only exception is participant P9 who had a poor performance and a very large variability in the ball release parameters. Indeed, the validity of the quadratic approximation depends on the nature of the non-linearities of the action-to-score function and the range of the deviations from the mean, i.e. from the relative spatial scales characterizing the concentration of the action distribution and the Hessian. Thus, if behavior is very erratic, our decomposition may become inaccurate for the entire set of actions and may be restricted to a more concentrated subset. However, considering that participants in our sample were untrained throwers, it is noticeable that the quadratic approximation was good for all but one out of twenty participants. This suggests that our methods could be safely applied to more controlled tasks, e.g. in evaluating athletes performances (as athletes do not typically exhibit high variability in motor actions) or in assessing motor skill learning, where training tends to quickly reduce motor variability. However, in cases where the variability is high, such as during the initial exploration phase when learning a novel motor skill, one could apply the decomposition on local clusters of data and still provide a compact characterization through the parameters of Hessian-based decomposition of each cluster. We plan to develop such approach based on a mixture of Hessian-based decompositions and to apply it to the investigation of throwing skills learning in future work.

Our decomposition method relies on the computation of the action covariance and the Hessian of the action-to-score function. These matrices and some of the parameters of the decomposition depend on the choice of the coordinate system in action space. In particular, the noise η and the tolerance τ, being defined as traces of the covariance and Hessian matrices, respectively, change under coordinate transformations (unless a metric is chosen [25]). However, the α term is a scalar (i.e. is a single number corresponding to the score associated to the mean action) and it does not depend on coordinates. The β term is the trace of the product of the covariance and the Hessian matrices and it is invariant under affine coordinate transformations, given that Σ^a and H^a transform in opposite ways (see Section B in S1 Appendix). Thus, re-scaling of positional and velocity coordinates due to different choices of measurement units does not affect the decomposition of mean score as a sum of α and β. However, β is not invariant, in general, for non-linear coordinate transformations, such as the transformation from Cartesian to polar coordinates. Indeed, the dependence on action coordinates has raised concerns about the reliability of the TNC decomposition [26, 27] and of UCM and GEM methods [28]. While such dependence may provide an opportunity to evaluate the role of different coordinate systems for control [29], it has also been noticed that geometric properties of the action-to-score function such as the solution manifold do not depend on coordinates [28]. In our decomposition, if the mean of the action distribution is on the solution manifold (α = 0), β is invariant also under non-linear transformations, because the non-linear term in the transformation of H^a depends on the gradient of the action-to-score function, which is null on the solution manifold. Moreover, if the action distribution is not centered on the solution manifold but it is concentrated (i.e. η is small) the change in β due to non-linear coordinate transformations may be negligible.

A further limitation of our approach, which is shared with GEM and TNC, is that the decomposition cannot reveal the temporal structure of inter-trial fluctuations, as multiple trials are needed to compute the mean and the variance of an individual’s action strategy. Recent work addresses this issue with an inter-trial error correction model that predicts both the temporal and geometric structure of variability near the goal equivalent manifold (GEM) of a simplified shuffle board task [30]. Furthermore, the variability analysis is shown to be coordinate-independent as the characterization is performed in the eigenspace of the error-correcting controller matrix.

3.3 Applications to motor skill learning

In this work we have focused on characterizing steady-state performance and individual action distribution during short experimental sessions rather than on skill improvement over multiple sessions. Future work will include longitudinal studies to understand if and how the observed inter-individual differences are related to the time course and the magnitude of individual performance improvements and skill learning. Current theories of human sensorimotor control suggest the existence of two distinct mechanisms underlying motor skill learning: a model-based system that improves motor performance guided by an internal forward model of the body and the environment, which is updated based on prediction errors [31]; and a model-free system in which learning is driven by reinforcement and punishment of successful/erroneous actions [32, 33]. Motor adaptation studies, in which a systematic perturbation of the environment is introduced by means of force fields of visuomotor rotations, suggest that the model-based system is responsible for the quick adaptation/compensation of the mean error. The model-free system, driven by reinforcement and punishment, regulates instead motor variability, and is hence responsible for the slow reduction of the variable errors. However, the interplay between this two learning mechanisms, remains poorly understood.

In analyzing our free overam throwing data, we have highlighted the existence of iso-performing participants, such as P1, P11 and P15, which have the same mean score, but different contributions of α and β. Do inter-individual differences in terms of α and β translate into individual differences in terms of performance improvement? In future work we plan to use the proposed framework to study the acquisition of throwing skills in virtual reality environments in which we can alter both the dynamics of the ball, for instance by manipulating the (virtual) gravity field, as well as the task score geometry, in this work assumed quadratic and isotropic in both task directions. As adapting to an altered dynamics requires learning a new forward model while a new task geometry changes the reward function, the dissociation between these two contributions might allow us to dissociate between model-based and model-free learning and to understand how initial inter-individual differences in terms of performance, variability and score tolerance translate into individual performance improvement.

3.4 Summary and conclusions

We have introduced a novel method to characterize the key features of the distribution of high-dimensional and redundant actions that affect performance, defined as the mean of the score assigned at individual actions. We have applied the method to the investigation of inter-individual differences in unconstrained throwing. We found that the indicators derived from the Hessian-based decomposition allow to identify specific and consistent features relating individual throwing strategies to different performance level and to understand how different strategies achieve similar performances. Participants differed in their throwing performance because they consistently differed in either the score of their mean action (bias, α) or the level of their tolerance-variability index (β), a measure of the interplay between action variability and action-to-score geometry. The same performance could be achieved trading off α with β and the same β could be achieved trading off θ/τ with η. In sum, the compact characterization of the relation between high-dimensional, redundant, and noisy action distributions and performance provided by our Hessian-based decomposition may be applied to a variety of complex real-life motor skill, opening up new opportunities both for systematic investigations of inter-individual differences in real-life motor skills and for practical applications to training of complex motor tasks. As motor control investigators, we plan to address in future studies how the individual parameters of the Hessian-based decomposition relate to individual motor learning capabilities, a first step to unravel the fundamental mechanisms underlying individual learning differences. A sport trainer could also use the Hessian-based decomposition indicators of a trainee, based on the analysis of a few tenths of repetitions of a specific motor skill, to select an optimal individualized training strategy according to those indicators.

4.1 Performance as expected value of non-linear and high-dimensional action scores

Let us assume that, at every trial t, an individual generates an action a_t ∈ A ⊂ Rⁿ with some random noise such that the action distribution can be described by a probability density function (p.d.f.) p_A. Let us also assume that, at every trial, the action receives a score point π_t = s^a(a_t) through a score (or loss) function s^a: a ∈ Rⁿ → R⁺ which punish motor actions according to some optimality criteria such as task errors or metabolic cost of an action. For instance, in Fig 1, the score function assigns a penalty score to an action, that is the squared distance between the action outcome x(a) and the target position x_T. We define the performance of a motor strategy p_A as the expected score: (10) which is the sum over all possible actions of the probability of a taking a given action times its score (integrated in Rⁿ). In this framework, an individual motor strategy p_A is optimal if, over repeated attempts, minimizes the average score. Similarly, given two motor strategies we can establish if they are equivalent or if one is better than the other by simply comparing their expected score. In practice, the score function may be highly non-linear and the action space high dimensional (n > >1) making difficult to find an analytic solution to (10). Therefore, we do not usually solve (10) but instead approximate the expected score with its sample average as in Fig 1B, where best performers, such as P18 and P4 have the lowest average quadratic error. Nevertheless, the sample average by itself is not informative enough to understand why the individual action distributions of P18 and P4 are much more performing compared to the majority of participants.

4.2 Hessian-based decomposition of motor skills

4.3 Application to throwing tasks

Throwing skills, as many other motor skills, are usually assessed by means of score functions which essentially define the objective of the throwing task. For instance, for a javelin thrower, the score may be a function of the longitudinal distance travelled by the javelin. The further the javelin lands, the larger the score assigned to the throwing action. Conversely, for a dart thrower, the goal is not to throw the dart as far as possible, but as accurate as possible, and hence, as in Fig 1A, the score function could assign a penalty increasing with the distance between the landing position of the projectile and the center of the target.

It should be noted that in our experimental protocol [22] participants did not receive any explicit performance feedback (or score) at end of each throwing trial (but they could see the arrival position of the ball on the target board) and therefore, in this work, in line with computational and experimental evidences [1, 31], we assume that participants optimize an accuracy score which penalizes the squared error between the outcome x of a release action and the target position x_T [31]: (25)

Written in this form, the accuracy score, represents a task score s^x: R² → R⁺ which penalizes the two-dimensional action outcomes x with a scalar score π. To find the relationship between release actions and quadratic (task) error, i.e. the action-to-score function s^a: a ∈ R⁶ → R⁺, we need to express the task score (25) as a function of the release parameters.

Assuming a point-mass projectile and hence neglecting friction and Magnus’s forces, the projectile trajectory f(a, t) can be predicted from the release parameters a = [p₀, v₀] as: (26) with p₀ = [x₀, y₀, z₀]^T and v₀ = [v_x0, v_y0, v_z0]^T.

For a target board oriented as in Fig 1A, i.e. with the normal pointing in the longitudinal direction y, the time of impact of the ball with board can be estimated as , i.e. as the ratio between the longitudinal distance of the projectile with the board at release (y_b is the coordinate of the board with respect to the world frame) and the velocity of the ball along such direction (that is constant according to the model in (26)). Hence, at the time of impact, the projectile will hit the board at: (27)

Substituting the above system of equations into (25), let us writing the action score as: (28) i.e, as a scalar function a ∈ R⁶ → R⁺ of the throwing action a: (29)

Given an individual release strategy with mean action and covariance Σ^a and aiming at hitting a desired target x_T, the mean squared error (E[e²] = E[π]) can be approximated with (15): (30) where is just (29) evaluated at , i.e. the quadratic error of the outcome of the mean action, and β can be decomposed into the three components η, τ, and θ by using the covariance matrix of the action strategy Σ^a and the 6 × 6 Hessian of (29) evaluated at .

In this work, the 6 × 6 Hessian matrix of (29) is calculated with the Matlab Symbolic Toolbox for each target condition and for each individual strategy.

4.4 2D throwing simulation

To illustrate the Hessian-based decomposition in simple task with 2D actions, we considered a toy model of throwing. We assumed that a projectile is released from a fixed position and moves in a vertical plane to hit a target on a vertical axis at a given distance from the release position (Fig 3A). Thus, the action a corresponds to a 2D vector whose components are the horizontal release velocity () and the vertical release velocity (), the outcome x corresponds to the arrival height on the vertical axis, and the score the squared distance from the target (x_T).

We simulated distributions of release actions with a given mean and covariance Σ_a (Figs 3 and 4) by drawing 50 random samples from a multivariate Gaussian distribution. We simulated a target at a distance of 5 m from the release position and at height 0.7 m below it, reproducing the mean conditions for target T1 in our sample of 6D throwing participants (see below). In addition to decomposing the mean score according to Eq (30) we performed TNC-Cost, UCM, and GEM analyses. For the TNC-Cost analysis, we followed the procedure described in [20]. For the UCM analysis and GEM analysis, we computed the Jacobian of the action-to-outcome function at the mean action and the orthogonal direction.

4.5 6D throwing experimental protocol and data analysis

Individual release strategies were obtained from the experimental dataset acquired in our previous study [22], where twenty right-handed participants (10 females, 10 males; age: 28.2±6.8 years) performed a series of overarm throws, starting from a fixed initial position. All participants signed an informed consent form in accordance with the Declaration of Helsinki. The data collection was carried out in accordance with Italian laws and European Union regulations on experiments involving human participants. The protocol was approved by the Ethical Review Board of the Santa Lucia Foundation (Prot. CE/PROG.542). Participants were instructed to hit one of four circular targets arranged on a vertical target board placed at 6 m from the initial position (marked with a sign on the floor) and to start from a fixed posture (standing with the arm along the body). The four targets were custom made and consisted in white circles of 40 cm diameter, arranged on a rectangular layout on the target board. The distances between the centers were 70 cm vertically and 80 cm horizontally, similarly to Fig 1. Moreover, the targets midpoints in the horizontal direction were shifted with respect to the projected initial position of the participant: the left and right targets were centered respectively at 60 cm to the left and 20 cm to the right of the projected initial position of the thrower’s midline. An opto-electronic system (OptiTrack, NaturalPoint, Inc., Corvallis, OR, United States) operating at 120 Hz was used to capture whole-body kinematic information of the participants throwing actions and the corresponding ball trajectories.

For each trial and participant, the release action a, was obtained by fitting each of the three spatial components of the ball path with a 3rd-order polynomial function, and therefore the release position p₀ and velocity v₀ were obtained from the zero and first order coefficients, respectively. Then, this release action was used off-line in (27) and (25) to generate ‘ideal’ ball paths and scores, respectively, which were not influenced by friction and/or spinning effect of the ball. Trials in which the ball path did not intersect the target plane, or for which the ball was partially tracked by the optical system, were excluded from the analysis (13% of the total number of throws; we verified that differences in the number of trials across participants and conditions affect the values of some of the decomposition parameters only for sets with fewer trials than our smallest set). The error distribution, across trials, participants and target conditions, between experimental and ideal performance (mean squared error) is shown in Fig 9 (mean ± SD: 0.0012±0.0912m²). The dataset is available in S1 Dataset.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Error distribution between experimental score and score predicted with the no-drag model in Eq (29).

Black and red vertical dashed lines indicate mean and ±1SD of the distribution, respectively.

https://doi.org/10.1371/journal.pone.0253626.g009

We next assessed the validity of the assumption that the individual action distribution is sufficiently localized around the mean action such that the second-order approximation of the sample mean score given in Eq (29) as the sum α + β is adequate. To do so we computed the fraction of variance accounted for (VAF) by the approximation, defined as: (31)

VAF is therefore defined as the variance of the error between the individual performance (sample mean score) and the second order approximation normalized by the variance of the population performance.

In summary, for each participant and for each target x_T, we estimated the mean-quadratic-error E[π], the mean action , and the action covariance Σ^a, with their respective sample mean and covariance: , , and , where N is the total number of successful actions, or trials, executed for target x_T.

4.6 Statistical analysis

To assess the effect of target and participant identity on the score and on the indicators extracted from the Hessian-based decomposition, we fit a linear mixed-effects model (Matlab function fitlme) with target and participants identity (categorical variables) as fixed effects and participants intercept as random effect (32) and we performed an analysis of variance for linear mixed-effects model (Matlab function anova), thus taking into account the repeated measures design.