Uncertainty-aware Self-supervised Learning for Cross-domain Technical Skill Assessment in Robot-assisted Surgery

Formally, we denote the source domain as a group of labeled samples, $\mathcal{D}_{L}:=\{(x_{1},y_{1}),\cdots,(x_{i},y_{i})\}=\{(X,Y)\}$, that were collected from basic skill training. Each $i$-th sample pair is composed by the input kinematics data $x_{i}\in\mathbb{R}^{d}$, associated with a corresponding skill class label $y_{i}\in\{1,\cdots,K\}$. Here, $d$ is the dimension of input data and $K$ is the number of classes. Similarly, we denote the target domain as a set of unlabeled samples $\mathcal{D}_{U}:=\{x_{1},\cdots,x_{j}\}=\{X^{\prime}\}$, $x_{j}\in\mathbb{R}^{d}$. Target domain was collected from a previously unseen exercise involving robotic surgical skills, where no ground truth labels of target domain are available.

In general, domain adaption refers to a set of approaches predicting the labels of samples from target domain $\mathcal{D}_{U}$, using both labeled samples from source domain and unlabeled samples from the target domain itself. As follows, the goal of transferable skill assessment is to train a model $f_{\theta}(\cdot)$ overs $\mathcal{D}_{L}$ and $\mathcal{D}_{U}$ that enables learning of both transferable features and domain-invariant classifier. Here, the output of $f_{\theta}(\cdot)$ is the predictive class probability of skills and $\theta$ denotes the learned network parameters.

Self-supervised learning, also known as self-learning, refers to a protocol that is specifically designed to learn from both labeled and unlabeled data [39]. Typically, self-learning approaches first construct an initial classifier based on the labeled data. Then, the classifier is applied to some of unlabeled data and make predictions using the model. These predictions, namely $pseudo$ labels, are included to re-train the existing classifier in an iterative fashion until a certain stopping criterion is satisfied.

We set up a similar procedure and devoted to domain adaptation for transferable skill assessment: train a classifier on the labeled source domain, produce approximate pseudo labels of surgical exercises in target domain. Next, these pseudo-labeled target samples are treated as labeled ones and used for re-training a new classifier in target domain. Then, the classifier explores the remaining set of target domain until model learning converges or the maximum iteration number is reached. In essence, self-learning considers the information from both domains and enables to propagate knowledge from the labeled source to the unlabeled target. Therefore, it could not only exploit the hidden structures in motion kinematic data that regularize learning, but help to find a shared feature space that matches data distributions of both the domains.

We train the whole network (with parameters $\theta$) by minimizing the loss between the model output $f_{\theta}(x)$ and the labels (true labels from source and $pseudo$ labels from target), while penalizing large activation. Here, $y$ denotes the true labels in source domain $x\in\mathcal{D}_{L}$, and $y^{\prime}$ denotes the pseudo labels of unlabeled samples in target domain $x\in\mathcal{D}_{U}$, Thus, the overall loss $\mathcal{L(\theta)}$ is a weighted mixture of the cross-entropy loss of labeled source domain, $\mathcal{L}_{labeled}$, the cross-entropy loss of pseudo-labeled ones from target domain, $\mathcal{L}_{pseudo}$, and a $L2$ weight regularizer term:

where $\alpha$ and $\lambda$ are two non-negative hyper-parameters to balance contributions of three loss components. In particular, the parameter $\alpha$ controls the importance of $pseudo$-labeled target samples relative to the true labeled source samples. A larger $\alpha$ could encourage selecting more $pseudo$-labels for training, whereas a smaller $\alpha$ reduces the impact of target samples for training.

Jointly learning classifiers and optimizing pseudo labels on unlabeled target domain is naturally difficult since it is not possible to completely guarantee the correctness of generated pseudo labels given an existing model. The uncertainty could come from the variability in the relationship between observed tool kinematics data and the corresponding skill class outputs. Conventionally, $pseudo$ labels were selected for self-learning by taking the class which has the maximum softmax class probability or the one with corresponding probability that is larger than a predefined threshold (chosen experimentally). The softmax function approximates relative class probabilities, but it did not capture any uncertainty or confidence information regarding the prediction outputs [40]. To cope with this issue, our strategy is first to estimate the model uncertainty of skill predictions and then seek for pseudo labels from the most confident ones in unlabeled target domain. These high-confidence pseudo labels presumably approximate the underlying true target labels and could thus better help self-learning adaptation.

Bayesian probability theory has offered us a statistically principled way to infer the uncertainty within deep learning tools [41]. Research has shown that any deep learning model can be arbitrarily interpreted into Bayesian ones, which are often referred as Bayesian neural networks. In Bayesian settings, the model uncertainty, also referred to as epistemic uncertainty, could be readily derived from the variances of prediction with respect to the distribution of network parameters $\theta$, i.e., posterior $p(\theta|X,Y)$. The distribution characterizes the uncertainty or confidence of model that conforms to the observed data $(X,Y)\in\mathcal{D}_{train}$. Although the posterior $p(\theta|X,Y)$ is analytically intractable, a popular approach, variational inference, can approximate it with a simpler approximating distribution that is much easier to sample with, while minimizing the Kullback-Leibler (KL) divergence between the approximating distribution and the true posterior [42]. Kendal et al. [43] show that the variational inference using stochastic Monte-Carlo dropout (MC dropout) could allow to approximate Bayesian posterior and minimize the cross-entropy loss of a network with dropout is equivalent to minimizing the KL divergence. Differently from the standard dropout that is used for regularization in training, MC approach stochastically drops out parameters of existing model at test time, which is equivalent to impose a Bernoulli distribution on parameters, to obtain the posterior approximate.

To allow uncertainty estimate for self-learning, we extend and cast our skill model into the Bayesian setting and implement MC dropout to obtain the uncertainty estimates over the current network. In practice, estimating prediction uncertainty was accomplished as follows: (1) Given an existing classifier with network parameters $\theta$, we run the network with a chosen number of stochastic forward passes $T$, where each unit has a probability $p$ of being set to 0. (2) We collect the softmax probability outputs for each skill class in the $t$-th forward pass ($t=1,\cdots,T$). (3) From the outputs over $T$ forward passes, the final prediction of class probability and uncertainty are taken to be the average and variance of predictions of each class, as defined in Eq. 2 and Eq. 3, respectively. Further, we derive the entropy of per-class uncertainty as the overall uncertainty measure given the input data.

where $f_{\hat{\theta}}(\cdot)$ denote $f_{\theta}(\cdot)$ when including MC dropout, $\hat{\theta}$ is the current network parameter $\theta$ masked by the dropout during inference.

Taken together, Algorithm 1 depicts the overall procedure in details. Initially, the inputs are the labeled samples from source domain $\mathcal{D}_{L}$ and unlabeled samples of target domain $\mathcal{D}_{U}$. After initialization with pre-trained network parameters $\theta^{\prime}$, our approach performs the following major steps: (1) Given current model parameters, MC-dropout variational inference was conducted to query the uncertainty of predictions over the unlabeled samples $x\in\mathcal{D}_{U}$. (2) Each unlabeled sample are assigned with $pseudo$ labels by taking the class which has the maximum average class probability. (3) Subset of target domain $\mathcal{D}_{U}^{\prime}\subset\mathcal{D}_{U}$ in which the target samples with the first $k$ lowest amounts of uncertainty are added to the current training set $\mathcal{D}_{train}\leftarrow\mathcal{D}_{train}\cup\mathcal{D}_{U}^{\prime}$ and simultaneously removed from unlabeled target domain $\mathcal{D}_{U}\leftarrow\mathcal{D}_{U}-\mathcal{D}_{U}^{\prime}$. (3) Re-train a new classifier using the updated training set $D_{train}$, and update network parameters $\theta$. The classifier explores the remaining unlabeled data $\mathcal{D}_{U}$ with smaller confidence until learning converges or maximum iteration is reached.

In our experiment, we fixed MC-dropout probability as $p=0.5$ for each dropout layer, the total number of forward passes for inference was set as $T=50$ as we sample $T=50$ prediction outputs for each input, and the number of selected target samples ($pseudo$-labeled) in each training iteration was $k=50$. The hyper-parameters were chosen empirically as they have shown a faster training convergence for modeling. As our model may no longer achieve the best classification for source domain after adapting to the target domain, we adopt the term ”converge” to denote the best model for which it is confident enough to make a prediction for the target domain.

Fig 2 shows the overview of our proposed architecture. The network input is the motion kinematics measured from robot end-effectors and the output is an estimated softmax probability of skill classes that the input data is drawn from, along with an estimate of prediction uncertainty.

The network consists of two main components where the input data flows through in parallel: a convolutional component and a recurrent component. The convolution component comprises two stacked convolutional layers that extract local patterns along the length of time-series data. Each convolutional layer is followed by a ReLU activation and a standard dropout layer with 0.2 dropout rate to reduce the risk of overfitting. A global average pooling layer is then applied to sum out the spatial information. It enables the network to handle input sequences of varying lengths. Also, the global average pooling could largely reduce the number of network parameters and thus help to alleviate overfitting. In parallel, the recurrent component consisted of two bi-directional LSTMs followed by a standard dropout layer. The two sets of features from both components are concatenated and then fed through a dense layer. Finally, the network is forked at the last layer to have two outputs: the skill probability output through a layer with softmax nonlinearity, and an uncertainty estimate obtained from the aforementioned Bayesian variational inference.

Two datasets in the field of robotic surgical training were used: the first dataset is denoted as the source domain, which contains labeled data with the skill annotations from basic surgical training tasks; the other is the target domain which contains unlabeled data collected from a VR simulation-based surgical training task that was characterized by the application of adaptive haptic guidance during its execution. Both datasets contain tool kinematics that were collected from daVinci robotic surgical systems. The kinematics of two master tool manipulators (MTMs) and two corresponding patient-side manipulators (PSMs) were monitored and captured as multi-dimensional time-series. All kinematic variables were collected with 30 Hz sampling frequency. Specifically, the source domain data contains the joint kinematic variables from each robot manipulator, including Cartesian position (3), rotation matrix (9), linear velocity (3), angular velocity (3), and gripper angle (1), resulting 76-dimensional time-series per trial, while the target domain data contains all variables except for the measures of joint linear velocity, resulting 48-dimensional time-series per trial. To align with the data measurements in two datasets, we only considered the 48 common kinematic variables. The two datasets are described below and summarized in Table I.

Before feeding the motion kinematics into the network, we performed pre-processing for the raw data in two steps. First, we down-sampled data samples from both source and target domains by a factor of 30 (from the original sampling rate 30 Hz to 1 Hz). The purpose of down-sampling is to reduce unnecessary computational load. Then, normalization was also applied to each channel of all kinematics variables by subtracting the mean and dividing by the corresponding standard deviation.

We implemented our methodology in Keras with Tensorflow 1.9 API as backend based on Python 3.6. All experiments in this study were remotely run on the UTSouthwestern BioHPC server and utilized a computing cluster equipped with a NVIDIA Tesla P100 GPU with 16 GB memory and 72 CPUs with 256 GB memory running a Linux kernel. To initialize the network, we first pre-trained the network on the source domain data given the known skill labels in JIGSAWS as ground-truth. We selected the best model as the one with the minimized loss on the validation set. For training, we utilized the following hyper-parameters: 100 training epochs, stochastic gradient decent as chosen optimizer with initial learning rate as 0.001, first and second momentum of 0.9 and 0.999, and weight decay of $10^{-8}$. Before fitting the model in each epoch, the training set was randomly shuffled to achieve a robust process of model learning.

Since no ground-truth skill labeling is currently available in the target domain dataset, we perform a statistical analysis in order to check for the effectiveness of pseudo labeling for extracting skill patterns. We hypothesize that the generated pseudo labels, given the skill knowledge has been generalized from the source to the target domain in the model, can provide meaningful significant results to reveal the motor skill learning pattern during a surgical training, and also be in align with prior studies that robotic assistance would help improve the proficiency of surgical training in general. We performed a two-way ANOVA analysis to determine the significance of skill probability between the training groups of subjects, the training sessions, and the interaction between the two. Similarly, another two-way ANOVA analysis was performed to check the significance differences between the training groups (Assisted and Non-assisted), the training stage (Pre-training and Post-training), and the interaction between the two. A post-hoc Tukey comparison was used to compare different levels and highlight the significant effects within each group.

Fig. 6 (A) and (B) present the expert probability of robot-assisted and non-assisted subjects per repetition (left, averaged across corresponding group) and per session (right, averaged across all repeated trials), respectively. Table II summarizes all the statistics.

In general, skill learning of both assisted and non-assisted subjects was characterized by an increase of expert class probability (or, equivalently, a decrease of beginner class probability) as the surgical training proceeded. Significant difference was found between the assisted and non-assisted subjects on the skill probability, $p<0.05$. Specifically, the subjects with robotic assistance were associated with considerably higher probability of expert after the initial three sessions. On the other hand, non-assisted subjects who did not receive any haptic assistance showed a continuous increase of technical skills within the first four sessions. However, non-assisted subjects had the larger variances of skill probability after the $4^{th}$ session. This result can be explained by the fact that the non-assisted group underwent a management for self-improvement as the surgical training repetitions proceeded. Nevertheless, it might be more difficult for trainees to pick up the most representative, expert-relevant motion patterns without any assistance, and thus resulted in a less efficient learning outcome at the final stage of training. Importantly, the assisted subjects had a trend with continuous decrease of expert probability in the first three sessions (maximum haptic assistance). This result might be due to the fact that the advent of haptics may change the trainees’ intended motion and consequently interfere their initial learning momentum, and thus trainees would need to adapt to the changes for an improvement. Nevertheless, the assisted subjects demonstrated a considerable improvement in their skills in and after the $4^{th}$ session. The differences of skill probability in the first three session and the $4^{th}$ session can be explained by the increased familiarity of subjects such that significant improvements are shown.

Fig. 7 compares the expert class probability between the initial training (first session) and post training (last session). The figure highlights the motor learning outcome as the result of surgical training. As shown, either assisted or non-assisted subjects have significantly higher expert probability after the training compared to the initial training session. This result confirms that repeated training practices in general could help trainees in achieving a potentially higher expertise. Note that the steeper learning curve, as revealed by the dotted lines, indicates that the robotic assistance potentially increased the rate of learning for performing the complex motor task. Even though the assistance intensity was adaptively diminished after $3^{th}$ session and no assistance was provided after $6^{th}$ session, trainees were able to improve their technical skills more efficiently and to retain this improvement even after the assistance removal.

Fig. 8 shows the averaged prediction uncertainties in each session of VR surgical training. The value of uncertainty is determined by the MC-dropout variational inference and reflects the level of predictive confidence when outputting a specific skill label given an input of motion kinematics. In general, surgical skills in the initial phase of surgical training are relatively more uncertain than the final stage. As the surgical training proceeds, the level of uncertainty for assessing skill decreases. The uncertainty could come from the variability between the observed motion data of trainees and the corresponding skill output. As the consequence of iterative practice, trainees could demonstrate distinguished motion patterns that provide a higher certainty of their skill levels, with or without robot assistance. We also note that non-assisted subjects are relatively more uncertain than the assisted subjects, especially at the end of surgical training. One interpretation is that between the assisted and non-assisted group, the motion profiles of assisted subjects could demonstrate more consistent signatures that are characteristic of high expertise; in contrast, subjects without robotic assistance might carry less distinguishable skill information in the motion that is hard to measure.

In the absence of skill labels in the target domain dataset, we evaluated the validity of the generated pseudo labels by comparing their probabilities with the task performance metrics. These metrics, which include completion time, translational root mean square errors, and rotational root mean square errors, were used as a proxy measure of technical skills to indicate motor learning in the surgical training task. The task performance metrics were collected in each training trial and the average performance for all subjects is presented in Fig. 9. As shown, a clear learning curve of the trainees is captured by the task performance metrics. Consistent with our skill predictions, the non-assisted subjects improved their performance after multiple training sessions. Additionally, the assisted subjects showed better performance improvement due to the contribution of haptic assistance.

However, the evolution of the users’ task performance metrics differs from the ones inferred from our skill predictions in terms of their slope and saturation. Specifically, the task performance metrics tend to be more sensitive to changes in assistance. We observed that, for the assisted subjects, the task performance measures rapidly improved when the maximum haptic assistance was activated, but showed relatively larger variances after the fourth session when the assistance intensity decreased. In contrast, as measured by the skill probability, the assisted subjects showed long-term persistence at the final stage of the training, indicating stable and consistent skill learning over time, even with a decrease in the intensity of haptic assistance. The discrepancy between the predicted skill and task performance metrics makes sense. The task performance metrics accurately reflect the direct impact of the physical assistance provided by the robotic system through its haptic guidance. The purpose of this guidance is to align the surgeon’s current pose with the ideal one, resulting in optimal task performance metrics, by minimizing the translation and orientation errors. On the other hand, the predicted skills reflect a more complex motor learning process and changes in these skills may not be immediately reflected in the surgeon’s performance. This could be considered an advantage of our approach for skill assessment, as it allows for the assessment of skills without any ”assistance bias” in training scenarios that include robotic feedback or guidance.

This study focuses on adapting skill models from a fully-labeled dataset to a VR simulated training exercise and aims to investigate the transferability of skills across physical and virtual environment. The source data consists of only the labeled data of standardized training exercises, but its limited size may impact the generalizability of the skill models. To address this, future work could involve collecting a larger dataset for pretraining. The selected surgical training exercises from the source and target domains in our study are different but limited in term of the complexity. For future work, we plan to extend the analysis to include more complex tasks in boarder fields of surgical training, such as physical box trainers and wet-lab exercises [46]. We believe that a comparison of the skill models presented in our study and those measured on physical box trainers could help to provide insights into the training proficiency and outcomes across different environments [47, 48]. Additionally, a more granular analysis of the skill models at surgical phase and task level might be helpful to reveal the corresponding surgeons’ skills and workflow patterns [7, 49]. Furthermore, it would be interesting to explore the potential of this approach to generalize from technical to non-technical skill assessment in surgical training [50, 51].