Hybrid Attention Model Using Feature Decomposition and Knowledge Distillation for Glucose Forecasting

The task of blood glucose forecasting with multi-modal input data can be formalized as a time series prediction problem. Let $X=\{x_{0},x_{1},\ldots,x_{n}\}$ represent the set of $n$ feature/sensor observations in the Sensing module. The $j^{th}$ feature observation $x_{j,0:t}=\{x_{j,0},x_{j,1},\\ \ldots,x_{j,t}\}\in\mathbb{R}^{n}$ represent the feature vector of $j^{th}$ sensor’s output up to time $t$. In (Shuvo and Islam, 2023) identified that BGL from CGM devices, insulin, and Carbohydrate (Carbs) intake features have prominent effects on forecasting blood glucose levels. However, combinations of other features, such as skin temperature and heart rate, either have no improvement or are sometimes overfitting for forecasting BGLs. Therefore, Our proposed method (GlucoNet) aims to utilize combinations of these measurements, i.e. BGL data, carbohydrate (Carbs) intake, and insulin dosage, to develop a forecasting model for blood glucose prediction. The corresponding output time series for multi-output forecasting is denoted as $x_{g,t+1:t+h}^{\prime}=\{x_{g,t+1}^{\prime},x_{g,t+2}^{\prime},\ldots,x_{g,t+h}^{\prime}\}\in\mathbb{R}$, assuming $g^{th}$ sensor is CGM, which represents multiple future blood glucose levels values across a given prediction horizon (PH). The prediction horizon refers to the length of time into the future for which predictions are made based on past data. In our evaluation, we test our approach using a prediction horizon of $5$, $30$, and $60$ minutes. Given that CGM data is usually recorded at $5$-minute intervals, a PH of $5$ minutes corresponds to $h=1$ samples, a PH of $30$ minutes corresponds to $h=6$ samples, while a PH of $60$ minutes corresponds to $h=12$ samples.

To ensure consistent evaluation across both output settings, we calculate the error metrics, such as root mean square error (RMSE), by comparing the actual future glucose level with the predicted values in the estimated multi-output sequence. Mathematically, we can express the forecasting task as $x_{g,t+1:t+h}^{\prime}=f(X)$, where $f(\cdot)$ and $\hat{x}_{g,t+1:t+h}^{\prime}$ represent the forecasting model and predicted future glucose values, respectively.

To improve the accuracy of our proposed forecasting of blood glucose, we combined raw CGM values, which are time-series features, with event-based features (e.g., carbohydrate and insulin intake). We proposed transforming event-based features into continuous time-series data. This transformation aims to better represent the physiological effects of carbohydrate intake and insulin on blood glucose levels over time and decrease the effect of sampling irregularities of these data.

We transformed the sparse meal events into a continuous ‘operative carbohydrate’ feature that represents the ongoing effect of the intake of carbohydrates on blood glucose (Kraegen et al., 1981). The carbohydrate intake is transformed from an event-based feature to continuous time series data by Eq. 1. The transformation assumes that carbohydrates begin affecting blood glucose $15$ minutes after the meal, peak at $60$ minutes, and gradually decrease over $3$ hours.

Here, $t_{s}$, $t_{meal}$, and $t_{peak}$ present the sampling time, time when the meal is taken, and time when $C_{op}$ reaches its maximum value, respectively. $C_{op}$ is the effective carbohydrates at any given time, and $C_{meal}$ is the total amount of carbohydrates in the meal. Based on the (Kraegen et al., 1981), $\alpha_{inc}$ is the increasing rate, which is set to $0.11$, and $\alpha_{dec}$ is the decreasing rate set to $0.028$.

The other variable significantly impacting blood glucose levels is insulin dosage. Similar to meal carbohydrates, we transform insulin dosage into a continuous active insulin feature. This transformation extended from (Boiroux et al., 2010), based on the insulin activity curve. This transformation accounts for insulin activity’s gradual onset, peak, and decay, providing a more accurate representation of insulin’s effect on blood glucose over time. The active insulin at any given time is modeled using Eq. 2.

Where $t_{s}$ is the sampling time, $t_{d}$ is the total duration of insulin activity, $\tau$ is the time constant of exponential decay, $a$ is the rise time factor, $S$ is the auxiliary scale factor

the other terms of Eq. 2 are calculated as follows, $\tau=t_{p}\times\left(1-\frac{t_{p}}{t_{d}}\right)/\left(1-2\times\frac{t_{p}}{t_{d}}\right)$, $a=2\times\tau/t_{d}$, and
$S=1/\left(1-a+(1+a)\times e^{-\frac{t_{d}}{\tau}}\right)$.

Fig. 3 shows an example of transforming the event-based features (e.g., Carb intake and insulin dosage) to the continuous active features using the SSR module. By converting meal and insulin events into continuous values, we aim to capture the dynamic relationships between these factors and blood glucose levels to improve the accuracy of our prediction model.

Variational Mode Decomposition (VMD) is an adaptive, non-wavelet multi-resolution decomposition technique that decomposes a complex signal into a set of band-limited intrinsic mode functions. VMD is an improved version of Empirical Mode Decomposition (EMD) and does not decompose the time series data with fixed functions. However, unlike EMD, VMD is a non-recursive algorithm that concurrently extracts multiple modes, with the number of modes ($m$) determined by the user based on the complexity of the problem (Liu et al., 2022). These $m$ modes are identified by signal peaks in the frequency domain. Each mode is centered around a specific pulsation frequency ($\omega_{m}$), which is essential for accurate analysis (Dragomiretskiy and Zosso, 2014).

For a time series $X(t)$ where $X$ is a vector of size $T$, VMD decomposes the time series into $m$ modes, resulting in a new matrix $X_{new}$ of size $m\times T$, where $m$ presents the number of modes. we extend the VMD decomposition technique proposed in (Dragomiretskiy and Zosso, 2014) to decompose the blood glucose time-series data to decrease the complexity of this signal. The VMD formulation for a single time series $X(t)$ can be expressed as:

subject to $\sum_{k=1}^{K}u_{k}=X(t)$, where $u_{k}(t)$ represents the decomposed modes, $\omega_{k}$ are the center frequencies for each mode, $\delta(t)$ is the Dirac distribution, and $*$ denotes convolution.

To solve this optimization problem, it is converted into an unconstrained Lagrangian form:

where $\mathscr{L}$ is the Lagrangian function, $\lambda$ represents the Lagrange multipliers, and $\alpha$ is the balancing parameter for the data-fidelity constraint.

The Alternate Direction Method of Multipliers (ADMM) is employed to solve this problem iteratively. The updates for each mode, center frequency, and Lagrange multiplier are performed in the frequency domain as follows (Dragomiretskiy and Zosso, 2014):

Detailed explanations of VMD are presented in (Dragomiretskiy and Zosso, 2014). Note that, in our proposed method, we apply the VMD to the training dataset and testing dataset separately to avoid information leakage from the test dataset to the training dataset. After applying VMD on data, the decomposed blood glucose signals are divided into two groups. The first group has less fluctuation and includes the trend of the signal, while the second group has higher fluctuation. The first group’s modes are summed up to generate a new signal with a trend and low fluctuation of blood glucose values called low-frequency BGL ($x_{g}^{low}$). The second group’s modes are also summed up to generate a new signal with high fluctuation of blood glucose features called high-frequency BGL ($x_{g}^{high}$). Both $x_{g}^{low}$ and $x_{g}^{high}$ are then combined with effective carbs intake ($X_{c}$) and insulin dosage ($X_{ins}$), which are computed in our SSR module (Section 3.2), to create two new datasets. These two datasets are called low-frequency features dataset ($LFFD$) and high-frequency features dataset ($HFFD$), respectively.

Both $LFFD$ and $HFFD$ are then fed to the forecasting models separately to predict the blood glucose levels. The forecasting employs a one-step model for various ranges of prediction horizons (PHs) such as Short-range ($5minutes$), medium-range ($30minutes$), and long-range ($60minutes$). A CNN-LSTM model is used for forecasting $LFFD$. The reason behind selecting this model is that the $LFFD$ has lower fluctuation compared to $HFFD$. Therefore, forecasting these data is possible using a model such as CNN-LSTM. This model is called the low-frequency forecasting model. The low-frequency forecasting model architecture consists of a series of layers designed to process and analyze time-series data for blood glucose prediction. The input to the model is a windowed sample of past data with a length of $180$ minutes, which is determined to be the optimal input length for all prediction horizons. The input data first passes through a stack of one-dimensional convolutional and pooling layers. This initial stage is designed to extract relevant features from the time-series data. We used two groups of two back-to-back convolutional layers to extract important features based on their relevance to the prediction task. Following the convolutional layers, the extracted features are flattened into a single vector. This vector is then fed into an LSTM (Long Short-Term Memory) block, which can made up of different LSTM layers with different numbers of memory cells or units. These settings offer multiple implementation configurations for the low-frequency forecasting model. Therefore, considering a set of numbers of memory cells can provide flexibility in adapting GlucoNet to various accuracy requirements and computational constraints. In our low-frequency forecasting model of GlucoNet, each configuration of the model has $k$ LSTM layers. Each layer is defined using an input number of memory cells and an output number of memory cells configuration vector, $L_{vec}=[(L_{1,1},L_{1,2}),(L_{2,1},L_{2,2}),...,(L_{k,1},L_{k,2})]$. Here, $(L_{i,1},L_{i,2})$ represents the input number of memory cells and the output number of memory cells in the $i^{th}$ LSTM layer. Hence, the generic low-frequency forecasting model of GlucoNet representation, GlucoNet $\{[(L_{1,1},L_{1,2}),(L_{2,1},L_{2,2}),...,(L_{k,1},L_{k,2})]\}$, defines some possible low-frequency forecasting model of GlucoNet configuration. Note that we considered two layers for our convolutional layers and two fully connected layers in the low-frequency model to decrease the size of design space and make it feasible to explore our model’s configuration. For instance, one of our low-frequency forecasting models of GlucoNet configurations has an LSTM layer with a $128$ input number of memory cells and $64$ output number of memory cells represented by GlucoNet $\{[(128,64)]\}$. Then, the LSTM block output is fed into two fully connected layers. The number of neurons in these layers varies depending on the prediction horizon, e.g., for $5-minutePH:[64\times 32,32\times 1],30-minutePH:[64\times 32,32\times 6],$ and $60-minutePH:[64\times 32,32\times 12]$ neurons. The low-frequency model architecture is able to progressively refine the extracted features and forecast BGL accurately for different time horizons.

Due to high fluctuations in $HFFD$, we need to propose a complex model to forecast high-frequency blood glucose levels. The Transformer, developed by Google, introduced a novel approach to sequence modeling by replacing the traditional RNN structure with self-attention (Vaswani, 2017). By considering self-attention, Transformers addresses the potential issues in RNN regarding long-sequence dependencies and inefficiencies in backpropagation through time. In order to maintain the order of elements in a sequence, the Transformer uses positional encoding, especially for time-dependent tasks. Therefore, this adaptation makes it suitable for time series applications such as blood glucose monitoring, where tracking important temporal features can enhance the detection of irregular glucose levels  (Lee et al., 2023). The Transformer is made up of the encoder and the decoder. The encoder transforms input data into a set of contextual representations using self-attention followed by a multi-layered feed-forward network to capture complex relationships across the sequence. Afterward, the output from the encoder is fed to the decoder, which generates an output sequence. This process involves self-attention to identify relevant aspects of the generated output. Taking advantage of these capabilities, in this study, we used the Transformer to forecast high-frequency fluctuations in blood glucose, which is hard to predict due to the rapid and unpredictable changes in glucose levels. However, the complexity of transformers makes these models unsuitable for real-time applications such as forecasting blood glucose.

Knowledge Distillation (KD) is a technique that facilitates the transfer of knowledge from one model, typically a larger teacher model, to another, usually a smaller student model. We used the benefit of using knowledge distillation approaches for training deep learning models, particularly in forecasting tasks. Knowledge distillation allows a student model to learn from a more complex teacher network. This process allows the student to learn from more informative sources, such as the predictive probabilities provided by the teacher. Hence, the student can achieve performance levels comparable to the teacher despite being a more compact model. In some cases, when the student and teacher have similar capacities, the student may even outperform the teacher (Hinton, 2015). Therefore, we extended Knowledge Distillation (KD) to compact the large Transformer (Teacher model) to the small Transformer (Student model) to achieve more accurate forecasting of blood glucose levels for high-frequency features.

For an input $\mathbf{x}$ and a $K$-dimensional target $\mathbf{y}$, a model generates a vector logits $\mathbf{z}(\mathbf{x})=[z_{1}(\mathbf{x}),\ldots,z_{K}(\mathbf{x})]$, which is then converted into predicted probabilities $P(\mathbf{x})=[p_{1}(\mathbf{x}),\ldots,p_{K}(\mathbf{x})]$ using a softmax function. Hinton et al. (Hinton, 2015) proposed temperature scaling to soften these probabilities for improved distillation using Eq. 6.

where $\tau$ is the temperature parameter. By adjusting the softmax output $P^{T}(\mathbf{x})$ of the teacher and $P^{S}(\mathbf{x})$ of the student, the student is trained using a combination of teacher probabilities and ground truth labels. Eq. 7 presents the loss function ($\mathcal{L}_{KD}$) used for training the student model.

where $H$ denotes Mean Squared Error (MSE) loss and $\alpha$ is a hyperparameter.

Knowledge distillation techniques address the issues related to using a complex model to achieve an acceptable performance in forecasting blood glucose model. In this regard, we proposed a large teacher model as a Transformer and a small student model to forecast blood glucose levels. The overall architecture of the teacher and student model is illustrated in Fig. 4.

We conducted a series of experiments to optimize the Transformer model architecture. For the teacher model, the best-performing configuration consisted of a single encoder layer, $64$ input dimensions, $4$ attention heads, and $128$ feed-forward units. Moreover, for the student model, the optimal configuration differed slightly. It utilized a single encoder layer, $32$ input dimensions, $2$ attention heads, and $64$ feed-forward units. These setups demonstrated superior performance compared to other configurations. Also, the high-frequency model of the proposed GlucoNet offers multiple implementation configurations. It can be deployed with or without knowledge distillation (KD) and implemented with the large Transformer (teacher Transformer) $(LT)$ or the small transformer (student Transformer) $(ST)$. These options provide flexibility in adapting GlucoNet to various accuracy requirements and computational constraints. Thus, the generic GlucoNet representation, GlucoNet $(+,)$ (KD,) $(ST,LT)$ $\{[(L_{1,1},L_{1,2}),(L_{2,1},L_{2,2}),...,(L_{k,1},L_{k,2})]\}$, defines some possible GlucoNet configuration. For example, GlucoNet $+$ KD $(ST)$ $\{[(32,16)]\}$ refers to a model that is made up of one LSTM layer with 32 input numbers of memory cells and 64 output numbers of memory cells, and implemented by KD with small Transformer as a student model and GlucoNet $(LT)$ $\{[(32,16)]\}$ refers to a model that is made up of one LSTM layer with 32 input numbers of memory cells and 64 output numbers of memory cells, followed by large Transformer (without KD).

Eventually, the outputs of the CNN-LSTM model and the Transformer are added to construct the forecasted blood glucose levels.

In summary, accurate forecasting of blood glucose levels for different ranges of PHs by GlucoNet is presented by Algorithm 1. First, CGM time series data are decomposed by using the VMD technique to low-frequency and high-frequency features (line $7$ of the Algorithm 1), and then the event-based features, such as carbohydrate intake and insulin dosage, are transformed into continuous time series data (line $9$ of the Algorithm 1). Then, the continuous time series data are combined by low frequency and high-frequency features and generate two sets of features. The $LFFD$ are fed to the LSTM model to predict the low frequency of blood glucose levels (line $10$ of the Algorithm 1), and $HFFD$ are fed to the knowledge distillation transformer-based model to forecast the high-frequency blood glucose levels (lines $11$ and $12$ of the Algorithm 1). Finally, the predicted blood glucose levels for a PH are obtained by combining the low-frequency and high-frequency blood glucose levels (line $17$ of the Algorithm 1).

The dataset that we used to demonstrate our proposed method is OhioT1DM (2018 and 2020 versions) (Marling and Bunescu, 2020). This is a clinical dataset with heterogeneous data, including diverse age, gender, and device variability, that reduces the risks of underfitting to deal with unseen data. Each version consists of comprehensive data on six individuals, which are a total of 12 participants (7 male, 5 female) with type 1 diabetes over an eight-week period. The cohort comprised two male and four female participants. Data collection was facilitated through Medtronic 530G insulin pumps and Medtronic Enlite CGM sensors, which were employed consistently throughout the study duration. The dataset encompasses various measurements such as Continuous Glucose Monitoring (CGM) readings at 5-minute intervals, insulin administration, including both bolus and basal doses, and self-reported data on meal times with estimated carbohydrate intake information. Moreover, other information includes exercise (E), sleep, work, stress, and illness. A summary of the dataset is provided in Table 2. Note that the missing values, especially for CGM values, are imputed with interpolation or extrapolation.

In this work, we used (1) Root Mean Square Error (RMSE) (Arefeen and Ghasemzadeh, 2023), (2) Mean Absolute Error (MAE) (Arefeen and Ghasemzadeh, 2023), and (3) $R^{2}$ Square (Syafrudin et al., 2022; Khadem et al., 2023) as the error metrics for comparing the accuracy of the different forecasting model. The definitions of these error metrics are presented below. Root Mean Square Error (RMSE) of forecasted values of blood glucose levels for specific Prediction Horizons (PHs) is defined as $\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{y_{i}})^{2}}$. Here, $y_{i}$ and $\hat{y_{i}}$ refer to the actual value and the predicted value, respectively. $n$ refers to the number of test samples.

Mean Absolute Error (MAE) of forecasted values of blood glucose levels for specific Prediction Horizons (PHs) is defined as $\text{MAE}=\frac{1}{n}\sum_{i=1}^{n}|y_{i}-\hat{y_{i}}|$.

$R^{2}$ of forecasted values of blood glucose levels for specific Prediction Horizons (PHs) is defined as $R^{2}=1-\frac{\sum_{i=1}^{N}(y_{i}-\hat{y_{i}})^{2}}{\sum_{i=1}^{N}(y_{i}-\overline{y_{i}})^{2}}$. Here, $\overline{y_{i}}$ presents the mean of the actual BGL.

We evaluated our proposed GlucoNet model compared with state-of-the-art models by implementing each model on the OhioT1DM dataset and comparing their performance using error metrics. The OhioT1DM dataset is divided into training and testing datasets for each participant. The models are trained on the training dataset. Prediction metrics are computed using the test dataset for each participant. The model’s input consists of a three-hour ($180$ minutes) sliding window of historical data, providing sufficient information for making good predictions for the next prediction horizon (PH). Various PHs are considered for comparison, which are $5$ minutes, $30$ minutes, and $60$ minutes. The LSTM model was trained for 300 epochs, repeated over 5 separate runs to ensure consistency and reliability of results. Moreover, our teacher model is trained for 500 epochs; then, the student model is trained by the knowledge distillation loss method for 500 epochs as well.

The training process began with an initial learning rate of $0.001$. The adaptive moment estimation (Adam) optimizer is chosen to minimize the linear combination of the MSE loss function in both LSTM and teacher models. The student model used the Adam optimizer to minimize the knowledge distillation loss function. Adam optimizer is selected due to its suitability for learning in non-stationary data, such as blood glucose. A minibatch training approach was adopted, using a batch size of 64. The model building and hyperparameter tuning were implemented using Python as the programming language and the Pytorch package.

Fig. 5 highlights a demonstration of forecasting the blood glucose model with $PH=60$ minutes for low-frequency, high-frequency, and blood glucose levels for a participant of the Ohio dataset. In Fig. 5, we can observe that the GlucoNet model forecasts the blood glucose levels accurately when PH is $60$ minutes.

We evaluate GlucoNet with three types of forecasting horizons - short (5 minutes), medium (30 minutes), and long-term (60 minutes). The average value of error metrics for different configurations of our GlucoNet compared to the baseline model is shown in Table 4. The configurations of GlucoNet that we considered are GlucoNet with or without knowledge distillation (KD) and implemented with the large transformer (teacher transformer) or the small transformer (student transformer). These configurations are GlucoNet $(LT)$ $\{[(128,128),\\ (128,64)]\}$, GlucoNet $(ST)$ $\{[(128,128),(128,64)]\}$, and GlucoNet $+$ KD $(ST)$ $\{[(128,128),(128,64)]\}$ . Moreover, the baseline model architecture is 1D-CNN, followed by LSTM. The 1D-CNN consists of two 1D convolution layers, each with a kernel size of 3. The layers have configurations of $Conv1d(18,64)$ and $Conv1d(64,128)$, respectively. The LSTM has two layers, which include $LSTM(128,128)$ and $LSTM(128,64)$. Then, the output of LSTM is fed to three fully connected layers.

Table 4 shows that GlucoNet deployed with knowledge distillation achieves efficient models without compromising performance metrics compared to GlucoNet without KD and implemented with the large Transformer.

Moreover, The average values of RMSE, MAE, and $R^{2}$ Square for different PHs (5 minutes, 30 minutes, 60 minutes) and configurations of our GlucoNet are shown in Fig. 6. In this figure, different configurations are considered different LSTMs with different numbers of memory cells or units. The configurations that are considered for comparison are GlucoNet $+$ KD $(ST)$ $\{[(128,64)]\}$, }, GlucoNet $+$ KD $(ST)$ $\{[(64,32)]\}$,GlucoNet $+$ KD $(ST)$ $\{[(32,16)]\}$, and GlucoNet $+$ KD $(ST)$ $\{[(16,8)]\}$.

By noting Fig. 6, we can conclude that the accuracy of the model optimized for configurations of GlucoNet $+$ KD $(ST)$ $\{[(64,32)]\}$ is evident. It is noteworthy that increasing the size or number of layers in the LSTM may lead to overfitting, which results in a decrease in the model’s accuracy.

Furthermore, we compared Gluconet with state-of-the-art models. Table 3 presents a comparison of GlucoNet with other state-of-the-art papers in terms of RMSE, MAE, and $R^{2}$ Square for PHs of 30 minutes and 60 minutes. The configuration of GlucoNet for this comparison is defined as GlucoNet $+$ KD $(ST)$ $\{[(128,128),(128,64)]\}$. In this table, the error metrics values are reported as average for comparison. This table highlights the performance of our proposed forecasting model.

The state-of-the-art models that we compared in Table 3 included deep learning models such as long short-term memory (LSTM), 1D CNN-LSTM (GlySim), deep residual network, recurrent self-attention, and convolutional recurrent neural networks (CNN-RNN). we can note that when $PH=30$ minutes, GlucoNet with KD achieved the lowest RMSE of $8.59mg/dL$ and MAE of $4.97mg/dL$. This outperformed the next best model, MTL-LSTM, which had an RMSE of $16.06mg/dL$ and MAE of $10.64mg/dL$. Also, when $PH=60$ minutes, it achieved an RMSE of $11.77mg/dL$ and MAE of $7.06mg/dL$, which is better than the GlySim, which had RMSE of $24.2mg/dL$ and MAE of $16.5mg/dL$. Therefore, the GlucoNet showed significant improvements compared with another model, about $47\%$ and $53\%$ in terms of RMSE and MAE when $PH=30$ minutes. Furthermore, when $PH=60$ minutes GlucoNet improvements for RMSE is about $51\%$ and MAE is about $57\%$. These results highlight the effectiveness of the proposed GlucoNet in capturing complex blood glucose dynamics and forecasting for longer-term predictions. Furthermore, the computational efficiency of the GlucoNet model compared to other models will discussed in the following section.

We have also compared our GlucoNet with other State-of-the-Art forecasting models in terms of accuracy-efficiency. We computed accuracy metrics, e.g., RMSE and MAE. Furthermore, the efficiency metric is computed by the total number of model parameters. Fig. 7 shows the accuracy-efficiency trade-offs of GlucoNet compared with state-of-the-art models.

We can note that GlucoNet configurations offer the best accuracy-efficiency trade-off compared to the other models. Our GlucoNet with KD achieves high accuracy with only about $10,900$ total parameters, which makes it suitable for potential edge deployment (Shuvo et al., 2022). Moreover, this model can be integrated into continuous glucose monitoring devices for real-time decision support.