11. CURAC-Jahrestagung, 15.-16. November 2012, Düsseldorf Efficient SVR model update approaches for respiratory motion prediction for the 11. CURAC 2012 R. Dürichen, T. Wissel, A. Schweikard University of Luebeck, Institute for Robotics and Cognitive Systems, Luebeck, Germany Contact: duerichen@rob.uni-luebeck.de Abstract: In order to successfully ablate moving tumours in robotic radiosurgery, respiratory motion prediction is needed to compensate time delays. In this context, recent studies revealed a high potential of support vector regression (SVR). However, high computational cost is one major drawback, particularly caused by batch mode training. We evaluate two approaches to reduce the update rate as well as computation time, while keeping a low prediction error. The update rules are either based on information about the respiratory phase or based on the current prediction error. An evaluation on patient data sets revealed that the second approach on average decreases computation time by 88.53% compared to a batch mode implementation. The prediction error increased by 0.3%, hence indicating enhanced efficiency. Keywords: radiosurgery, respiratory motion prediction, support vector regression, model update 1 Problem In modern robot-based radiation surgery, precisely radiating moving tumours has become more and more feasible, while sparing surrounding critical structures. One state-of-the-art system is the CyberKnife® (Accuray Systems, Sunnyvale, CA). Apart from radiating the tumour from multiple locations, the system can compensate internal tumour motion that is caused by respiration. For compensation optical markers, attached to the patient's chest, are continuously tracked. Their position is correlated with internal landmarks by stereoscopic X-ray images [1]. However, due to kinematic limitations, data acquisition and processing, time latencies have to be taken into account. In case of CyberKnife® Synchrony this latency is 115ms. Thus, to reduce this systematic error, the position of the optical markers has to be predicted. Several studies have shown that support vector regression (SVR) can precisely predict respiratory motion and is capable of competing with state-of-the-art motion prediction algorithms [2–5]. Implemented in usual batch mode (BM) the SVR model is recomputed at every incoming sample point. This leads to high computational cost and is therefore unpractical for real-time applications. To resolve this, Ernst et al. [2] suggested a global update factor τglobal, that reduces the update rate, but also increases the prediction error. In this study, we investigate two SVR update methods to decrease the computation time while maintaining a low prediction error. First, we analyzed the occurrence of resulting Support Vectors (SV) over time to obtain an update factor depending on the respiratory phase. Second, the current prediction error at time t is used as an update criterion. If the prediction error exceeds a certain threshold, the SVR is recomputed. 2 Methods For simplicity, all further explanations are reduced to a one dimensional signal. The algorithms can easily be expanded to three dimensions. It is assumed that the respiratory signal is equidistantly sampled with a sampling rate fs and that the resulting sequence has length N. Let yt be the signal amplitude at time point t and yt+δ the amplitude of the predicted point. Here, δ represents the prediction horizon. The predicted value is denoted by ̂ and depends on . 2.1 Support Vector Regression (SVR) Given a training set 107 11. CURAC-Jahrestagung, 15.-16. November 2012, Düsseldorf {( ) ( )} with { } (1) where again ui represents the training data and yi+δ the training labels, the SVR fits a functional relationship y = f(u) to this data. Here, L is the amount of training samples and M the number of past observations used to predict a new label. The ɛ-SVR does not penalize a fitting error up to ɛ per sample and optimizes f(u) to be as smooth as possible. In the linear case with slope vector w and offset b, f(u) can then be used to predict a point ̂ ( ) 〈 〉 , where 〈 〉 denotes the dot product. As smoothness is equivalent to minimizing the l2-norm of w, a convex optimization problem, that is regularized by violations of the mentioned ɛ constraint, can be formulated and solved to finally derive f(u). A full account of the underlying principles and mathematics is given in [6]. Training samples exceeding an absolute error of ɛ after performing the optimization are called Support Vectors (SV). Essentially, only these samples will contribute to the final definition of f(u) and are hence most important for the training or updating process. The SV machine can be extended to nonlinear cases by using kernel functions used to implicitly transform the linear case above to higher dimensions. The SVR was implemented with a Gaussian radial basis kernel (width parameter ɣ = 2) using LIBSVM Toolbox [7]. Further, we set the regularization constant C for prediction errors on the training data to 30, while L = 1000 and M = 5. The tube size ɛ was set to the standard deviation of the first á trous wavelet scale [2]. 2.2 Updating of SVR model Phase factor update rule (PFUR) In contrast to general BM training, Ernst et al. [2] introduced a global update factor τglobal, updating the model if ( ) is true, which happens at frequent intervals. Thus, an update factor = 1 is equivalent to the BM implementation. We will refer to this method as constant factor update rule (CFUR). Here, this idea is extended to a phase specific update factor. Therefore, we assume that the occurrence of the SVs is not homogenously distributed over time and that the sets of SVs used by the SVR model at time t and t-1 are similar. Note again, that only SVs will contribute information to the SVR model after training. Consequently, the effective update rate of SVR adapts to the phase of the respiration period, e.g. maximum in- and expiration. We define two phases: maximum in- and expiration with τmax = 1 and in- as well as exhalation phase with τslope < 1. The length of the first phase is determined by a proportion parameter θ given in percent of the mean period duration per maximum. These phases are centered on the maxima. Error update rule (EUR) The second update algorithm depends on the current prediction error. If at time t a value yt is measured, the algorithm computes ̂ based on the current SVR model and input ut-δ. If | ̂ | exceeds a threshold eth, the SVR model is updated. This aims at updating the model once it becomes too inaccurate. In order to estimate a general threshold independent of the signal, eth is set to a multiple f of the chosen tube size ɛ. 3 Results To investigate the PFUR and the EUR, we carried out three experiments on real patient data. First, we experimentally verified the assumptions being the basis of PFUR. Second, by varying τslope, θ and f, the potential of these two update algorithms is evaluated and compared to CFUR for the same signal. Third, EUR is further explored on 33 patient files. According to the latencies of the CyberKnife system, the prediction horizon δ was set to 3 samples (≈115ms). All signals have been globally scaled to a range of [0, 1] and reduced to the first principle component for simplicity. Investigations were done on an i7-2600 CPU@3.40 GHz with 16GB RAM and all algorithms were implemented in Matlab. The performance was evaluated using computation time tcalc and number of SVR model updates nupdate. Whereas tcalc depends on the specific hardware of the computer and software implementation of the algorithm, nupdate only depends on the algorithm and is approximately governed by a linear relationship to tcalc. Prediction accuracy was measured in terms of relative root mean square error RMSrel [2]: √∑ ( ̂) ⁄ (2) √∑ ( ) ⁄ 108 11. CURAC-Jahrestagung, 15.-16. November 2012, Düsseldorf a) b) mean shape of respiration period 1 0.4 mean y [mm] y [mm] 0.2 0.5 0 0 0 50 100 150 200 0 0.5 1 1.5 2 time [sec] time [sec] c) d) Training Labels - number of SVs:182 SV distribution rel. number of SVs 0.5 1 y [mm] 0.5 Training labels 0 0 0 200 400 600 800 SVs 1000 5 10 15 20 Trainings samples section Figure 1: Analysis of the SV location; a) part of the analysed real signal; b) mean shape of the respiration period; c) number and location of SVs within a training batch for one time step t; d) number of SVs sorted into 20 bins relative to the maximum number of SVs among all bins 3.1 Support Vector distribution over time In the first experiment, 20000 data points (≈770s) of a respiratory signal were processed by ɛ-SVR using a batch mode implementation (part of the signal is shown in fig. 1.a). At each time point, the number and location of the SVs were extracted from the SVR model. Fig. 1.c shows training labels and SVs for a certain time point t. The SVs are mainly distributed around the maxima in- and expiration. Over time the SVR model used an average amount of 446 from 1000 possible training samples as SVs. Per model update 3.3 elements from the SV set changed on anverage. Finally, a mean respiratory period was calculated from the signal (fig. 1.b) and divided into 20 bins. Depending on their location within the whole signal, each SV was assigned to one of the bins. Fig. 1.d shows a histogram of the accumulated number of SVs per bin normalized by the maximum number of training samples among all bins. 3.2 Comparing performance and efficiency of PFUR and EUR We further analyzed the signal used in experiment one with PFUR and EUR. For PFUR, τslope was varied between 0 and 1 in steps of 0.1 and θ was incremented in steps of 5% between 5-20%. The impact of the multiplication factor f for EUR was tested by grid search using f = {0.05, 0.1, 0.25, 0.5, 0.667, 0.8, 1, 1.25, 1.5, 2, 4, 10, 20}. The results were compared with CFUR using a global update factor τglobal varied between 0 and 1 in steps of 0.1 and assessed using RMSrel, tcalc and nupdate (Fig.2.a-f). 3.3 Analysis of multiple patients with error update rule (EUR) Based on the promising results of EUR, 33 patient files were evaluated only with this method. These signals were ran domly selected from a pool of 304 signals (31 patients) using 20000 data points for each analysis. As in section 3.2, the multiplication factor f was varied to find the optimal fopt. To compare the results for different patients i, the average RMSrel, tcalc and nupdate were calculated relative to CFUR τglobal =1: ̅̅̅̅̅̅̅̅̅ ( ); ̅̅̅̅̅̅ ( ); ̅̅̅̅̅̅̅̅̅̅ ( ) (3) Figures 3.a - c show the results averaged across 33 patients. In addition the results of CFUR with τglobal = 0.1, 0.5 and 1 a) b) c) 4 x 10 1500 2  = 10% 0.6 [%]  = 20% [s] update 1000 rel 0.5 1  = 30% calc RMS 500 n  = 40% t 0.4 0 0 CFUR 0 0.5 1 0 0.5 1 0 0.5 1 slope / global slope / global slope / global 4 d) e) f) x 10 0.55 1500 2 EUR [%] 0.5 CFUR, global = 1 [s] update 1000 rel 0.45 1 calc CFUR, global = 0.5 RMS 0.4 500 n t 0.35 CFUR, global = 0.1 0 0 -1 0 1 -1 0 1 -1 0 1 10 10 10 10 10 10 10 10 10 f f f Figure 2: Evaluation of the relative RMS error RMSrel, computation time tcalc and number of model updates nupdate for one patient; (a-c) PFUR for different update factors τslope and size parameters θ; (d-f) EUR for different multiplication factor f compared to a CFUR with different τglobal 109 11. CURAC-Jahrestagung, 15.-16. November 2012, Düsseldorf a) 1.6 b) 1 c) 1  EUR    1.4 CFUR, global = 1     0.5 0.5 1.2 CFUR, global = 0.5 1 CFUR, global = 0.1 0 0 -1 0 1 -1 0 1 -1 0 1 10 10 10 10 10 10 10 10 10 f f f Figure 3: Evaluation of 33 patients (≈770s per patient) with EUR, a-c) shows the avg. relative RMS error, avg. computation time and avg. number of model updates with respect to the results for a batch mode implementation for different multiplication factor f are illustrated. On average, EUR (f = 2) reduces tcalc by 88.53 %, and increases RMSrel by 0.3%, compared to batch mode training. The mean computation time to predict approximately 770s of respiratory signal was 534.4s. 4 Discussion In this work, we presented two updating methods to decrease computation time of an ɛ-SVR for respiratory motion prediction – one based on the current respiratory phase (PFUR) and one on the prediction error (EUR). In an initial analysis, we showed that the location of the SVs strongly depends on the phase in the respiratory period. We found that training samples at the maximum in- or expiration are likely to comprise more SVs than other periods. This result supports the initial assumption of PFUR, i.e. that the update rate can be reduced during in- and exhalation. The concept of this approach was studied for various τslope and θ on one data set. For τslope = 0.3 and θ = 15%, tcalc was reduced by 317. 4s compared to BM (nupdate,PFUR = 15159 instead of 20000 model updates), while the relative RMS error only increased by 0.34%. However, one drawback of this method is that the respiration phases, especially maximum in- and expiration, have to be known in advance or have to be identified online. This restricts the use of the PFUR for practical applications and makes it dependent on the accuracy with which the location of such extrema can be estimated. In contrast, EUR is independent of prior information and reduced tcalc by 1278.7s compared to an equivalent BM implementation, while yielding the same RMSrel increase as for the PFUR (nupdate,EUR = 111 for f = 4). The result was confirmed in an evaluation of 33 patients. The average computation time was reduced by 88.53 % from tcalc,BM = 4292.2s (BM) to tcalc,EUR = 534.4s, while the average RMSrel was increased by 0.3% only. With this method, respiratory motion prediction with SVR is performed more efficiently compared to the BM implementation. Even though it is still not possible to perform online motion prediction with SVR on an ordinary office computer, because the computation time for updating the model is larger than the sample time. But knowing when to update the model can be useful for implementing efficient multi-thread, online SVR predictors in the future. An example could be a dual-thread predictor, where one thread performs online predictions, while the second thread updates the SVR model. To reduce the calculation time even further, the algorithms could be implemented in C/C++ for instance or on dedicated hardware. Finally, it should be pointed out that EUR is not restricted to respiratory motion signals and can be applied to other applications. 5 References [1] A. Schweikard, G. Glosser, M. Bodduluri, M. J. Murphy, J. R. Adler, „Robotic Motion Compensation for Respiratory Movement during Radiosurgery“, Computer Aided Surgery, Nr. 5, 2000. [2] F. Ernst, A. Schweikard, „Forecasting respiratory motion with accurate online support vector regression (SVRpred)“, Int J CARS, Nr. 5, 2009. [3] W. D. D’Souza, K. Malinowski, H. H. Zhang, „Machine Learning for Intra-Fraction Tumor Motion Modeling with Respiratory Surrogates“, International Conference on Machine Learning and Applications, 2009. [4] N. Riaz, P. Shanker, R. Wiersma, O. Gudmundsson, W. Mao, B. Widrow, L. Xing, „Predicting respiratory tumor motion with multi-dimensional adaptive filters and support vector regression“, Physics in Medicine and Biology, Nr. 19, 2009. [5] A. Krauss, S. Nill, U. Oelfke, „The comparative performance of four respiratory motion predictors for real-time tumour tracking“, Physics in Medicine and Biology, Nr. 16, 2011. [6] A. J. Smola, B. Schölkopf, „A tutorial on support vector regression“, Statistics and Computing, Nr. 3, 2004. [7] C.-C. Chang, C.-J. Lin, „LIBSVM: A library for support vector machines“, ACM Trans. Intell. Syst. Technol., Nr. 3, 2011. 110