:
1 2022 Machine Learning for Health (ML4H) 2022

# A Temporal Fusion Transformer for Long-term Explainable Prediction of
Emergency Department Overcrowding ††thanks: This work was partially supported
by the strategic project NOVA LINCS (UIDB/04516/2020), the FCT project
DSAIPA/AI/0087/2018 and the Carnegie Mellon University - Portugal FCT project
CMU/TIC/0016/2021

Francisco M. Caldas<EMAIL_ADDRESS>Cláudia Soares
<EMAIL_ADDRESS>
NOVA School of Science and Technology  Caparica  Portugal

###### Abstract

Emergency Departments (EDs) are a fundamental element of the Portuguese
National Health Service, serving as an entry point for users with diverse and
very serious medical problems. Due to the inherent characteristics of the ED,
forecasting the number of patients using the services is particularly
challenging. And a mismatch between affluence and the number of medical
professionals can lead to a decrease in the quality of the services provided
and create problems that have repercussions for the entire hospital, with the
requisition of healthcare workers from other departments and the postponement
of surgeries. ED overcrowding is driven, in part, by non-urgent patients that
resort to emergency services despite not having a medical emergency,
representing almost half of the total number of daily patients. This paper
describes a novel deep learning architecture, the Temporal Fusion Transformer,
that uses calendar and time-series covariates to forecast prediction intervals
and point predictions for a 4-week period. We have concluded that patient
volume can be forecast with a Mean Absolute Percentage Error (MAPE) of 9.87%
for Portugal’s Health Regional Areas (HRA) and a Root Mean Squared Error
(RMSE) of 178 people/day. The paper shows empirical evidence supporting the
use of a multivariate approach with static and time-series covariates while
surpassing other models commonly found in the literature.

###### keywords:

Time Series Emergency Department Machine Learning Temporal Fusion Transformer
Forecasting Manchester Triage System Neural Network Explainable ML National
Health Service

##  1 Introduction

The forecast of the number of patients who use emergency services daily is
essential to determine in advance the human resources needed at hospital
Emergency Departments (ED). Multi-step ahead predictions allow hospital
managers to organise rotation schedules and diminish waiting times in urgent
care facilities (Wargon et al., 2009; Jones et al., 2008). When not accounted
for, overcrowding can lead to a decrease in the quality of patient care and
worse clinical outcomes (Bernstein et al., 2009; Hurwitz et al., 2014). From a
macro point of view, the influx in the emergency department combines an
expected number of people who are taken to the emergency room with a very
serious illness, for example, heart attack, with people that use the emergency
hospital to deal with non urgent problems, such as common cold, strained
muscles, or to deal with problems associated with chronic illness (Zachariasse
et al., 2019; Hurwitz et al., 2014). The most serious cases are reasonably
constant over time, and, predominantly, people in life threatening conditions
have no choice but to go to emergency care, thus the indicators of a rise in
patients with serious illnesses might not be the same for non urgent users. A
large number of patients that resort to urgent care are not, however, urgent,
according to the Manchester Triage system, used in the Portuguese National
Healthcare System. Roughly 40% of the patients are classified during triage at
the green/blue level, which means not urgent. Unlike more urgent patients, the
influx of green/blue patients has several factors that follow well-defined
cycles. For example, it is easy to identify that the day with the most influx
of non-urgent patients is Monday, with a smaller number of patients pursuing
emergency care during the weekend (Batal et al., 2001; Holleman et al., 1996;
Rathlev et al., 2007). To combine the predictive power of Deep Neural Networks
with the explainability usually reserved for simpler algorithms, we will use a
recently developed machine learning model to predict the influx of non-urgent
patients: the Temporal Fusion Transformer (TFT) (Lim et al., 2021); and study
which variables, time-series or not, had the most impact on the model, and
thus which are most relevant to predict daily patient volume.

Previous studies have examined the multi-step forecasting of daily patient
volumes (Jones et al., 2008; Diehl et al., 1981). Most focus is on the use of
classical statistical tools for temporal linear regression such as moving
averages (Milner, 1988), and their many extensions, namely ARIMA, SARIMA or
VARIMA (Afilal et al., 2016; Schweigler et al., 2009; Whitt and Zhang, 2019;
Carvalho-Silva et al., 2018). In recent years, with the advent of machine
learning, newer studies have been conducted that use neural networks (Jones et
al., 2008; Zhou et al., 2018), or otherwise other machine learning techniques
to tackle the same problem (Rocha and Rodrigues, 2021; Navares et al., 2018;
Tuominen et al., 2022). From the use of Feed-forward Neural Networks (Jones et
al., 2008; Navares et al., 2018), to Recurrent Neural Networks (Kadri and
Abdennbi, 2020; Harrou et al., 2020), 1-D Convolution Neural Networks
(Sharafat and Bayati, 2021), and later to Long Short-Term Unit (LSTM)
(Sudarshan et al., 2021; Harrou et al., 2020), there has been a constant
advance in the field, from linear models to deep neural network models.

In most studies using ARIMA and its variants, it was found that calendar
variables (day, day of the week, holidays) have a significant contribution to
model results (Boyle et al., 2012; Hertzum, 2017; Jones et al., 2008; Wargon
et al., 2009). A more detailed Literature Review can be found in Annex A.

##  2 Data Analysis

fig:blue_white

Figure 1: Time series from January, 2019 to February 20, 2020. We can observe
the weekly cycle, as well as annual trends and volume variation according to
RHA.

In this section, we will present the dataset used in this work. The data was
obtained from the public database “Transparência
SNS”111https://transparencia.sns.gov.pt/explore/dataset/atividade-sindrome-
gripal-csh
https://transparencia.sns.gov.pt/explore/dataset/atendimentos-nos-csp-gripe
and refers to daily data of care in primary health centres together with daily
data of consultations and waiting times in hospitals’ emergency departments
(ED) across Portugal, divided by Regional Health Area (RHA). The dataset
covers the time period from November 1st, 2016 to July 31th, 2022; with 16
variables per observation that define, among other things, the Regional Health
Area (RHA), Área Regional de Saúde (ARS) in Portuguese, of the observation. In
total, the dataset contains information regarding the daily volume of patients
in emergency care, the number of scheduled and unscheduled consultations in
primary health facilities, the daily number of patients arriving at the
Emergency Department (ED) with respiratory issues, the waiting times between
triage and the first medical evaluation and categorical variables pertaining
to calendar information, such as weekend, day of week or national holidays.

In fig:blue_white we can observe the weekly variation in the number of non-
urgent patients, as well as the volume shift according to RHA. It is visible
that despite having different levels of affluence, the different RHA follow
the same trend, with peaks of affluence occurring on Monday, and reduced
volume on weekends and during the Summer months, usually associated with
vacations. This aspect of the data served as motivation for the application of
a non-linear model over multiple time series, unlike well-established models
such as ARIMA.

##  3 Study Setting

Now that we have presented the data used in this paper, let us define, and
expose, the reasoning behind the rules by which we will create and evaluate
the model.

  * •

Multivariate forecasting: we want a model that leverages data and forecasts
across different Regional Health Areas. In most research in the area, models
are usually restricted to a single Hospital, and a more general model, capable
of forecasting across different regions, might be able to uncover new
interactions in the data and increase robustness.

  * •

Long-time forecast: In order to add value at the hospital management level,
the forecast of the number of patients should not be limited to the following
day or week. In this paper, we have chosen a 4-week (28-day) forecast,
considering that it allows breathing room for management and personnel
decisions. To the best of our knowledge, few papers have worked on such an
extended forecast horizon Boyle et al. (2012); Carvalho-Silva et al. (2018),
with only partial success.

  * •

Probability prediction: besides obtaining an estimate of the most likely value
in the future, a model that presents a probability density function on the
prediction conveys much more information. Of special value is, for example,
the definition of confidence intervals, which can transmit to those who use
the model an idea of the confidence, or precision, of the model in its
estimation.

  * •

Explanatory variables: Importantly, we want to evaluate the predictive
capacity of different variables, determining up to which passed time-step the
model finds predictive value or which covariates, categorical or numerical,
have a significant impact on the prediction.

##  4 Models

The first and simpler method used for forecasting is the replication of the
last $k$ time-steps. This technique, which is used as Baseline in this paper,
is also referred to as the naïve algorithm. By evaluating this model on the
validation, the optimal value for $k$ was estimated to be $7$, thus
representing the weekly periodicity in the data.

For comparison, other models commonly used in this area were also applied,
namely AutoRegressive Integrated Moving Average (ARIMA) with a seasonal
component (Jones et al., 2008; Schweigler et al., 2009; Eyles et al., 2022),
and its multivariate variant Vector AutoRegressive Integrated Moving Average
(VARIMA) (Kadri et al., 2014).

Also used was the exponential Smoothing algorithm, a simple method that has
also shown good results in the literature (Champion et al., 2007). Finally, to
gauge the performance of common machine learning models, the XGBoost model was
used. Out of these models, the XGBoost (Ke et al., 2017) (a Decision Tree
Boosting algorithm), is the only model capable of using past and future
covariates, with the disadvantage of not being specifically tailored for time-
series data. All these baselines models have statistical variants, or
extensions, that can estimate a prediction interval, or otherwise a
probability distribution that can in turn be used to define prediction
intervals.

Initially introduced by Lim et al. (2021), the Temporal Fusion Transformer
(TFT) model instantiated a novel architecture, combining a few mechanisms
previously only used separately, in a single model. The key features of the
TFT are:

  * •

Variable Selection Network: three independent Selection Networks, one for each
variable set (static, past and future), to select only relevant variables at
each time-step. This module removes noisy variables that do not add predictive
value, while giving some level of insight into the variables that are more
significant to the prediction;

  * •

A Gating Mechanism to skip any other element of the architecture. For specific
cases where exogenous variables are not useful or there is no need for non-
linear processing (e.g. in very simple forecasts), the Gating Mechanism, also
referred to as Gated Residual Network (He et al., 2016), allows the model to
only use non-linear processing when needed;

  * •

Static variables encoding to combine static information with time-series data;

  * •

Temporal Dependency Processing to capture short-term dependency, with an LSTM
encoder-decoder (Fan et al., 2019), and long-term dependency using a Multi-
Head Attention mechanism (Vaswani et al., 2017). By an additive aggregation of
the different heads, this mechanism gains explainability, as the weights in
the aggregated Multi-head represent time-step importance;

  * •

Confidence Intervals: the output of the models are quantiles, that define
prediction intervals, at each forecast time-step.

The hyperparameters of the model, alongside a more in-depth explanation of its
features and a diagram of its architecture, are detailed in Annex A.4.

##  5 Results

In this section, we present the results of the TFT and five other models for a
4 week forecast window over a 6 months period of time, from January, 24 to
July, 31 2022. tab1 illustrates how the TFT outperforms other common models in
the literature for long time prediction, with a Root Mean Squared Error (RMSE)
of 178.19, or approximately 178 people per day. This metric, however, might be
deceptive, as it is scale dependent, meaning that RHAs with a larger daily
volume will necessarily yield a higher RMSE, and skew the results. The Mean
Absolute Percentage Error (MAPE) on the other hand, is scale invariant, and it
better depicts the overall predictive power of the models, with the TFT
obtaining a 9.87% percentage error. To better compare the models, we utilised
an empirical CDF for each model, as seen in fig:ecdf. In this Figure,
depicting Absolute Error, the TFT shows overall slightly better performance.
We also acknowledge that the VARIMA algorithm obtains somewhat good results
and that the Gradient Boosting Tree obtains a better score with its 95%
confidence interval, however this is in part due to the fact that no
observation falls outside of the confidence, which might indicate that the
interval is to wide. By analyzing the results, the authors also found that the
statistical methods, such as ARIMA or Exponential Smoothing, had degraded
results in the last 2 weeks of the forecast window, unlike the Gradient
Boosted Tree method and the TFT. In Table LABEL:tab2 this effect can be
observed by evaluating the MAPE error per week.

This illustrates the superior capability of machine learning to perform long
term prediction, as the complexity of the model helps to identify long term
patterns, and identify break points where the trend shifts.

tab1  Models MAE RMSE MAPE MSE MIS Baseline 160.85 192.01 11.14 56615 - Exp.
S. 189.32 222.71 12.40 83522 1114 ARIMA 160.48 189.79 10.60 61234 2848 VARIMA
159.01 187.43 10.32 56776 - XGBoost 192.48 233.34 12.17 103569 730 TFT 151.60
178.19 9.87 55559 768

Table 1: Forecast error in the 24/01/2022 – 31/07/2022 period. To evaluate the
models, five metrics are used: Mean Absolute Error (MAE), Root Mean Squared
Error (RMSE), Mean Absolute Percentage Error (MAPE) ,Mean Squared Error (MSE)
and Mean Interval Score (MIS). Bold indicates the best result; TFT is
consistently more accurate than the baselines.

fig:ecdf

Figure 2: Empirical Cumulative Distribution function for the absolute
difference between the true value and the predicted value, for all RHA.

But the strength of the Temporal Fusion Transformer used goes beyond the
precision of the model. The ability to give some insight into how the forecast
was obtained is just as important. Two key findings can be summarized from the
results: first, we can observe the attention given to each time-step and
determine that the most recent time-steps are given in average a higher
weight, which is intuitively expected and shows that old information has less
value to the model. Secondly, the two most relevant covariates are the
percentage of patients in the emergency room with respiratory problems and
holidays. It is possible that an increase in the percentage of patients with
respiratory problems might be indicative of a spike in influenza/COVID-19
transmission, thus indicating future positive trends in the number of non-
urgent cases. During holidays the number of patients using the emergency
department drops significantly, so we would expect the model to give
particular importance to this one feature. A more detailed analysis of the
covariates and the accompanying Figures is presented in Annex B.

tab2  Models week 1 week 2 week 3 week 4 Total Baseline 7.16 10.01 12.84 14.55
11.14 Exp. S. 5.72 10.30 14.95 18.65 12.40 ARIMA 8.57 9.40 11.72 12.70 10.60
VARIMA 8.09 9.47 11.22 12.50 10.32 XGBoost 10.64 11.94 12.44 13.66 12.17 TFT
8.81 9.91 9.88 10.82 9.87

Table 2: MAPE error by week; Exponential Smoothing showed the best results in
the first week, followed by quick degradation for longer forecast periods.

##  6 Conclusion

This paper presented a novel application of the Temporal Fusion Transformer
(TFT) model to predict non-urgent patient volume in Portuguese public
hospitals by Health Regional Areas (HRA). The results were encouraging,
surpassing other models commonly found in the literature (Jones et al., 2008;
Kadri et al., 2014). The forecasting of an entire month is seldom done in the
literature (Carvalho-Silva et al., 2018; Aboagye-Sarfo et al., 2015), and it
could be further enhanced by a thorough analysis of model deterioration over
time. The introduction of a multivariate model with good results across groups
is a positive prospect, since one limitation of univariate time-series is the
natural low-data regimen, while multivariate models can merge information from
multiple sources, thus increasing the total amount of data fed to the models.
In the future, this model can increase in granularity, forecasting at the
hospital level instead of aggregated values by HRAs. Although a greater
challenge, due to the increased noise and randomness that comes from the
decrease in the study population, we expect that the combination of a large
number of time-series could improve the robustness and global quality of the
model, specifically if we add more relevant static variables. T For this
paper, only HRA and time-series statistics were used as static covariates, but
as noted in Farmer and Emami (1990), across different regions there is
distinct demand for emergency care, thus impacting the scale and variance of
the time-series. In future work, we plan to introduce other factors that might
contribute to encode region-specific information as static covariates, such as
demographics, modes of transport available, socio-economic characterisation of
the patient population and number and capacity of private healthcare providers
in the region. All these elements might help to represent each class, and
ultimately be used for a generalisation of the model to unseen hospitals,
where these variables might help to represent how similar a new unseen
hospital/RHA is to hospitals/RHAs in the training data.

## References

  * Aboagye-Sarfo et al. (2015) Patrick Aboagye-Sarfo, Qun Mai, Frank M. Sanfilippo, David B. Preen, Louise M. Stewart, and Daniel M. Fatovich.  A comparison of multivariate and univariate time series approaches to modelling and forecasting emergency department demand in western australia.  _Journal of Biomedical Informatics_ , 57:62–73, 2015.  ISSN 1532-0464.  https://doi.org/10.1016/j.jbi.2015.06.022.
  * Afilal et al. (2016) Mohamed Afilal, Farouk Yalaoui, Frédéric Dugardin, Lionel Amodeo, David Laplanche, and Philippe Blua.  Forecasting the emergency department patients flow.  _Journal of Medical Systems_ , 40(7):175, 2016\.  ISSN 1573-689X.  10.1007/s10916-016-0527-0.
  * Batal et al. (2001) Holly Batal, Jeff Tench, Sean McMillan, Jill Adams, and Phillip S Mehler.  Predicting patient visits to an urgent care clinic using calendar variables.  _Academic Emergency Medicine_ , 8(1):48–53, 2001\.  10.1111/j.1553-2712.2001.tb00550.x.
  * Bernstein et al. (2009) Steven L. Bernstein, Dominik Aronsky, Reena Duseja, Stephen Epstein, Dan Handel, Ula Hwang, Melissa McCarthy, K. John McConnell, Jesse M. Pines, Niels Rathlev, Robert Schafermeyer, Frank Zwemer, Michael Schull, Brent R. Asplin, and Emergency Department Crowding Task Force Society for Academic Emergency Medicine.  The effect of emergency department crowding on clinically oriented outcomes.  _Academic Emergency Medicine_ , 16(1):1–10, 2009\.  https://doi.org/10.1111/j.1553-2712.2008.00295.x.
  * Boyle et al. (2012) Justin Boyle, Melanie Jessup, Julia Crilly, David Green, James Lind, Marianne Wallis, Peter Miller, and Gerard Fitzgerald.  Predicting emergency department admissions.  _Emergency Medicine Journal_ , 29(5):358–365, 2012.  ISSN 1472-0205.  10.1136/emj.2010.103531.
  * Carvalho-Silva et al. (2018) Miguel Carvalho-Silva, M. Teresa T. Monteiro, Filipe de Sá-Soares, and Sónia Dória-Nóbrega.  Assessment of forecasting models for patients arrival at emergency department.  _Operations Research for Health Care_ , 18:112–118, 2018\.  ISSN 2211-6923.  https://doi.org/10.1016/j.orhc.2017.05.001.
  * Champion et al. (2007) Robert Champion, Leigh D. Kinsman, Geraldine A. Lee, Kevin A. Masman, Elizabeth A. May, Terence M. Mills, Michael D. Taylor, Paulett R. Thomas, and Ruth J. Williams.  Forecasting emergency department presentations.  _Australian Health Review_ , 31(1):83–90, 2007\.  10.1071/AH070083.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota, 2019. Association for Computational Linguistics.  10.18653/v1/N19-1423.
  * Diehl et al. (1981) Andrew K Diehl, Max D Morris, and Stephen A Mannis.  Use of calendar and weather data to predict walk-in attendance.  _Southern medical journal_ , 74(6):709–712, 1981\.  10.1097/00007611-198106000-00020.
  * Ekström et al. (2015) Andreas Ekström, Lisa Kurland, Nasim Farrokhnia, Maaret Castrén, and Martin Nordberg.  Forecasting emergency department visits using internet data.  _Annals of Emergency Medicine_ , 65(4):436–442.e1, 2015.  ISSN 0196-0644.  https://doi.org/10.1016/j.annemergmed.2014.10.008.
  * Eyles et al. (2022) Emily Eyles, Maria Theresa Redaniel, Tim Jones, Marion Prat, and Tim Keen.  Can we accurately forecast non-elective bed occupancy and admissions in the nhs? a time-series msarima analysis of longitudinal data from an nhs trust.  _BMJ Open_ , 12(4), 2022.  ISSN 2044-6055.  10.1136/bmjopen-2021-056523.
  * Fan et al. (2019) Chenyou Fan, Yuze Zhang, Yi Pan, Xiaoyue Li, Chi Zhang, Rong Yuan, Di Wu, Wensheng Wang, Jian Pei, and Heng Huang.  Multi-horizon time series forecasting with temporal attention learning.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’19, page 2527–2535, New York, NY, USA, 2019. Association for Computing Machinery.  ISBN 9781450362016.  10.1145/3292500.3330662.
  * Farmer and Emami (1990) R D Farmer and J Emami.  Models for forecasting hospital bed requirements in the acute sector.  _Journal of Epidemiology & Community Health_, 44(4):307–312, 1990.  ISSN 0143-005X.  10.1136/jech.44.4.307.
  * Gneiting and Raftery (2007) Tilmann Gneiting and Adrian E. Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American Statistical Association_ , 102(477):359–378, 2007.  10.1198/016214506000001437.
  * Harrou et al. (2020) Fouzi Harrou, Abdelkader Dairi, Farid Kadri, and Ying Sun.  Forecasting emergency department overcrowding: A deep learning framework.  _Chaos, Solitons & Fractals_, 139:110247, 2020.  ISSN 0960-0779.  10.1016/j.chaos.2020.110247.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, 2016.  10.1109/CVPR.2016.90.
  * Hertzum (2017) Morten Hertzum.  Forecasting hourly patient visits in the emergency department to counteract crowding.  _The Ergonomics Open Journal_ , 10(1), 2017.  10.2174/1875934301710010001.
  * Holleman et al. (1996) Donald R Holleman, Renee L Bowling, and Charlane Gathy.  Predicting daily visits to a waik-in clinic and emergency department using calendar and weather data.  _Journal of General Internal Medicine_ , 11(4):237–239, 1996.
  * Hurwitz et al. (2014) Joshua E. Hurwitz, Jo Ann Lee, Kenneth K. Lopiano, Scott A. McKinley, James Keesling, and Joseph A. Tyndall.  A flexible simulation platform to quantify and manage emergency department crowding.  _BMC Medical Informatics and Decision Making_ , 14(1):50, 2014.  ISSN 1472-6947.  10.1186/1472-6947-14-50.
  * Jones et al. (2008) Spencer S. Jones, Alun Thomas, R. Scott Evans, Shari J. Welch, Peter J. Haug, and Gregory L. Snow.  Forecasting daily patient volumes in the emergency department.  _Academic Emergency Medicine_ , 15(2):159–170, 2008.  10.1111/j.1553-2712.2007.00032.x.
  * Kadri and Abdennbi (2020) Farid Kadri and Kahina Abdennbi.  Rnn-based deep-learning approach to forecasting hospital system demands: application to an emergency department.  _International Journal of Data Science_ , 5:1–25, 2020\.  10.1504/IJDS.2020.10031621.
  * Kadri et al. (2014) Farid Kadri, Fouzi Harrou, Sondès Chaabane, and Christian Tahon.  Time series modelling and forecasting of emergency department overcrowding.  _Journal of Medical Systems_ , 38(9):107, 2014\.  ISSN 1573-689X.  10.1007/s10916-014-0107-0.
  * Ke et al. (2017) Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu.  Lightgbm: A highly efficient gradient boosting decision tree.  In _Advances in Neural Information Processing Systems_ , volume 30 of _NIPS’17_ , pages 3149––3157, Red Hook, NY, USA, 2017. Curran Associates, Inc.  URL https://proceedings.neurips.cc/paper/2017/file/6449f44a102fde848669bdd9eb6b76fa-Paper.pdf.
  * Koenker and Hallock (2001) Roger Koenker and Kevin F Hallock.  Quantile regression.  _Journal of Economic Perspectives_ , 15(4):143–156, December 2001.  10.1257/jep.15.4.143.
  * Lim et al. (2021) Bryan Lim, Sercan O. Arık, Nicolas Loeff, and Tomas Pfister.  Temporal fusion transformers for interpretable multi-horizon time series forecasting.  _International Journal of Forecasting_ , 37(4):1748–1764, 2021.  ISSN 0169-2070.  https://doi.org/10.1016/j.ijforecast.2021.03.012.
  * Makridakis (1993) Spyros Makridakis.  Accuracy measures: theoretical and practical concerns.  _International Journal of Forecasting_ , 9(4):527–529, 1993.  ISSN 0169-2070.  https://doi.org/10.1016/0169-2070(93)90079-3.
  * Milner (1988) PC Milner.  Forecasting the demand on accident and emergency departments in health districts in the trent region.  _Statistics in medicine_ , 7(10):1061–1072, 1988\.  10.1002/sim.4780071007.
  * Navares et al. (2018) Ricardo Navares, Julio Díaz, Cristina Linares, and José L. Aznarte.  Comparing arima and computational intelligence methods to forecast daily hospital admissions due to circulatory and respiratory causes in madrid.  _Stochastic Environmental Research and Risk Assessment_ , 32(10):2849–2859, 2018.  ISSN 1436-3259.  10.1007/s00477-018-1519-z.
  * Pines and Hollander (2008) Jesse M Pines and Judd E Hollander.  Emergency department crowding is associated with poor care for patients with severe pain.  _Annals of emergency medicine_ , 51(1):1–5, 2008\.  10.1016/j.annemergmed.2007.07.008.
  * Pulkkinen (2020) Eetu Pulkkinen.  _forecasting emergency department arrivals with neural networks_.  Bachelor’s thesis, Tampere University, Tampere, Finland, 2020.
  * Rathlev et al. (2007) Niels K Rathlev, John Chessare, Jonathan Olshaker, Dan Obendorfer, Supriya D Mehta, Todd Rothenhaus, Steven Crespo, Brendan Magauran, Kathy Davidson, Richard Shemin, et al.  Time series analysis of variables associated with daily mean emergency department length of stay.  _Annals of emergency medicine_ , 49(3):265–271, 2007.  doi:10.1016/j.annemergmed.2006.11.007.
  * Rocha and Rodrigues (2021) Carlos Narciso Rocha and Fátima Rodrigues.  Forecasting emergency department admissions.  _Journal of Intelligent Information Systems_ , 56(3):509–528, 2021.  ISSN 1573-7675.  10.1007/s10844-021-00638-9.
  * Schweigler et al. (2009) Lisa M. Schweigler, Jeffrey S. Desmond, Melissa L. McCarthy, Kyle J. Bukowski, Edward L. Ionides, and John G. Younger.  Forecasting models of emergency department crowding.  _Academic Emergency Medicine_ , 16(4):301–308, 2009.  https://doi.org/10.1111/j.1553-2712.2009.00356.x.
  * Sharafat and Bayati (2021) Ali R. Sharafat and Mohsen Bayati.  Patientflownet: A deep learning approach to patient flow prediction in emergency departments.  _IEEE Access_ , 9:45552–45561, 2021.  10.1109/ACCESS.2021.3066164.
  * Sudarshan et al. (2021) Vidya K. Sudarshan, Mikkel Brabrand, Troels Martin Range, and Uffe Kock Wiil.  Performance evaluation of emergency department patient arrivals forecasting models by including meteorological and calendar information: A comparative study.  _Computers in Biology and Medicine_ , 135:104541, 2021.  ISSN 0010-4825.  https://doi.org/10.1016/j.compbiomed.2021.104541.
  * Tuominen et al. (2022) Jalmari Tuominen, Francesco Lomio, Niku Oksala, Ari Palomäki, Jaakko Peltonen, Heikki Huttunen, and Antti Roine.  Forecasting daily emergency department arrivals using high-dimensional multivariate data: a feature selection approach.  _BMC Medical Informatics and Decision Making_ , 22:134, 2022\.  10.1186/s12911-022-01878-7.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.
  * Wargon et al. (2009) M Wargon, B Guidet, T D Hoang, and G Hejblum.  A systematic review of models for forecasting the number of emergency department visits.  _Emergency Medicine Journal_ , 26(6):395–399, 2009.  ISSN 1472-0205.  10.1136/emj.2008.062380.
  * Whitt and Zhang (2019) Ward Whitt and Xiaopei Zhang.  Forecasting arrivals and occupancy levels in an emergency department.  _Operations Research for Health Care_ , 21:1–18, 2019.  ISSN 2211-6923.  https://doi.org/10.1016/j.orhc.2019.01.002.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander Rush.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45. Association for Computational Linguistics, 2020.  10.18653/v1/2020.emnlp-demos.6.
  * Zachariasse et al. (2019) Joany M Zachariasse, Vera van der Hagen, Nienke Seiger, Kevin Mackway-Jones, Mirjam van Veen, and Henriette A Moll.  Performance of triage systems in emergency care: a systematic review and meta-analysis.  _British Medical Journal Open_ , 9(5), 2019.  ISSN 2044-6055.  10.1136/bmjopen-2018-026471.
  * Zhou et al. (2018) Lingling Zhou, Ping Zhao, Dongdong Wu, Cheng Cheng, and Hao Huang.  Time series model for forecasting the number of new admission inpatients.  _BMC medical informatics and decision making_ , 18(1):39, 2018.  ISSN 1472-6947.  10.1186/s12911-018-0616-8.

##  Appendix A Literature Review, Data and Models

###  A.1 Literature Review

As mentioned before in the main section of the paper, much work as been done
using ARIMA and other statistical models to forecast ED overcrowding, with
more recent work using Deep Learning techniques suited to time series such as
RNNs or LSTM (Kadri and Abdennbi, 2020; Harrou et al., 2020), and many works
pointing to the impact of calendar variables (Boyle et al., 2012; Jones et
al., 2008; Wargon et al., 2009) in improving forecasts. Weather data, such as
temperature and rain, has also shown predictive power for ED arrivals with
respiratory problems (Navares et al., 2018), but in others studies that
analysed the whole spectrum of ED visitors, it is either not a significant
variable, or it could be replaced by calendar variables, e.g. month of the
year (Hertzum, 2017; Wargon et al., 2009). Covariate interpretability is one
of the frequent drawbacks of Neural Networks, alongside the failure to
recognize long-term dependencies in time-series. One specific device that
addresses both problems is the Attention Mechanism (Vaswani et al., 2017):
simply put, it evaluates long-term dependencies and also represents how each
time-step impacts the model’s prediction. Attention has been used as part of a
specific Neural Network family of architectures called Transformers, that has
shown impressive results in the Natural Language Processing field (Devlin et
al., 2019; Wolf et al., 2020). In the literature, we found only one example
that used a Temporal Fusion Transformer model to predict Emergency Department
(ED) volume in one hospital for one day ahead (Pulkkinen, 2020). While not
being the only work that performed only daily predictions (Rocha and
Rodrigues, 2021; Tuominen et al., 2022), we find that a longer forecasting
window produces increased value for hospital management and poses a different
challenge from a machine learning perspective, as seasonal fluctuation needs
to be fully represented, and common forecasting models tend to decrease in
predictive quality as the forecast period becomes wider.

###  A.2 Data Analysis

Due to the time period and the type of data, it is important to consider the
Covid Lockdown period. Portugal’s response to the Covid-19 Pandemic resulted
in a series of containment measures that reduced travel and in person work: in
this period (10/03/2020 – 1/08/2021) the percentage of non urgent visits in
the RHA of Lisboa e Vale do Tejo (Lisbon and Tagus Valley) dropped from the
normal value of 48% of the total to 40%, with more dramatic drops for example
in the Algarve RHA from 45% to 30% at the beginning of the pandemic. This
dramatic period influenced the way people used emergency services, and it can
demonstrate how external factors influence people going to the emergency room.
This shift, associated with the general decline in the number of people in
urgent care, represents a distribution change in the time series, therefore
making it exceedingly difficult to predict the COVID period using only pre-
COVID information. In the same way, we can conclude that this COVID period
does not have useful information about the post-COVID future, and, in fact, we
have experimentally verified that the quality of the models decreased with the
introduction of the COVID period, thus leading to the decision to exclude this
temporal section from the training set.

It is, in a certain way, clear that the prediction of the influx in emergency
rooms can be useful for a more efficient management of hospital services, but
there is visible value added at user level, in the sense that they will get
better and faster care (Pines and Hollander, 2008). To sustain this claim, we
can observe the impact that the number of non-urgent patients has on the
waiting time before being treated in the Emergency Department. In
fig:time_blue, for the RHAs with the highest daily affluence, we observe a
positive, moderate to strong correlation between the number of non-urgent
patients and the waiting time. This is an indicator, not entirely unexpected,
that ED overcrowding of non urgent patients can lead to a substantial increase
in the average waiting time for all patients, urgent or non urgent.

fig:time_blue

Figure 3: Correlation between waiting times and non urgent patient volume. A
strong to moderate correlation exist between these two variables, therefore
implying that overcrowding increases waiting times.

###  A.3 Metrics

In this paper, we use five metrics to evaluate the models. The Mean Absolute
Error (MAE):

$\text{MAE}=\frac{1}{5}\sum_{k=1}^{5}\frac{1}{n}\sum_{i=1}^{n}|y_{i}^{k}-\hat{y}_{i}^{k}|$
(1)

the Root Mean Squared Error (RMSE),

$\text{RMSE}=\frac{1}{5}\sum_{k=1}^{5}\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_{i}^{k}-\hat{y}_{i}^{k})^{2}}$
(2)

the Mean Absolute Percentage Error (MAPE),

$\text{MAPE}=\frac{1}{5}\sum_{k=1}^{5}\frac{1}{n}\sum_{i=1}^{n}\frac{|y_{i}^{k}-\hat{y}_{i}^{k}|}{|y_{i}^{k}|}$
(3)

the Mean Squared Error (MSE),

$\text{MSE}=\frac{1}{5}\sum_{k=1}^{5}\frac{1}{n}\sum_{i=1}^{n}(y_{i}^{k}-\hat{y}_{i}^{k})^{2}.$
(4)

and finally the Mean Interval Score (MIS) (Gneiting and Raftery, 2007) to
evaluate the confidence interval with $\alpha=0.95$,

$\displaystyle\text{MIS}_{\alpha}=$
$\displaystyle\frac{1}{5}\sum_{k=1}^{5}\frac{1}{n}\sum_{i=1}\big{(}(\hat{U}_{i}^{k}-\hat{L}_{i}^{k})$
$\displaystyle+\frac{2}{\alpha}(\hat{L}_{i}^{k}-\hat{y}_{i}^{k})\mathbbm{1}(y_{i}^{k}<\hat{L}_{i}^{k})$
$\displaystyle+\frac{2}{\alpha}(\hat{y}_{i}^{k}-\hat{U}_{i}^{k})\mathbbm{1}(y_{i}^{k}>\hat{U}_{i}^{k})\big{)},$
(5)

where $\hat{L}_{i}^{k}$ and $\hat{U}_{i}^{k}$ are, respectively, the predicted
lower and upper bound and $\mathbbm{1}$ is the characteristic function. This
score rewards a small confidence interval, while applying a penalty whenever
the predicted result falls outside of the interval. Since we are evaluating
the predictions over several groups (RHAs), the total error will be the
average across RHAs. The most common metric across the literature for ED
forecasting is the Mean Absolute Percentage Error (MAPE) (Carvalho-Silva et
al., 2018; Ekström et al., 2015), however, when the true value is close to
zero, this metric becomes unreliable. It also places a heavier penalty on
negative errors (when the predicted value is higher than the true value)
(Makridakis, 1993). To overcome that, outliers values very close to zero are
removed for this particular metric. To correctly evaluate the out-of-sample
predictive capacity of the model, the dataset is divided into three subsets:
train, validation, and test. The training and validation sets combined
represent roughly 4 years, while the test set is roughly 6 months. The
validation set is used to optimise hyperparameters and to identify overfitting
during training, while the test set is unseen until the end and is only used
to produce the final results.

###  A.4 Temporal Fusion Transformer

fig:tft

Figure 4: TFT architecture. The inputs are static metadata, time-varying past
inputs (including past target values) and known future information. The
Variable selection unit selects the most relevant features, while the Gated
Residual Network allows to skip over unused sections of the architecture. The
interpretable multi-head attention is used to evaluate the most relevant time-
steps. Image adapted from Lim et al. (2021).

The novel model used in this paper is the Temporal Fusion Transformer. In
Section 4 of the main paper the key features of the model are defined,
starting with the Variable Selection Network. To understand this mechanism, we
first need to define the model input, and separate the input variables into
static, target and time-dependent. Static covariates, such as time-series
variance or mean, are specific to each group, i.e. RHA, and are defined as
$s_{i}$ with $i={0,...,4}$. $y_{i,t}$ is the target for group $i$ at time-step
$t$ and $x_{i,t}=[p_{i,t}^{T},f_{i,t}^{T}]^{T}$ the time-dependent covariates,
with $p$ representing past covariates, meaning covariates that are only known
until the present, and $f$ future covariates, that can be assumed to be known
in the past and the future, in our case, holidays and weekends. The Variable
Selection Network is independent for each type of input, and it produces
feature weights that will be used to interpret the results, with higher
weights ascribed to the most important covariates. The Masked Interpretable
Multi-head Attention is used to learn long-term relationships between
different time-steps and to determine which time-steps are more relevant to
the model prediction. Multi-head means that multiple attention mechanisms are
applied over the sequence, and it was presented by Vaswani et al. (2017) to
improve the learning capacity of the standard attention mechanism, however,
the use of multiple heads over different sections of the input sequence
hinders interpretability. A solution presented by the authors of the TFT model
(Lim et al., 2021) is to use multiple heads over the same sequence, and then
ensembled them into a single interpretable head, thus combining the learning
capacity of multi-head attention with the explainability of a single attention
mechanism.

The prediction function is defined as (Lim et al., 2021):

$\hat{y}_{i}(q,t,\tau)=f_{q}(\tau,y_{i,t-k:t},x_{i,t-k:t},s_{i})$ (6)

where $\hat{y}_{i}(q,t,\tau)$ is the predicted $q$th quantile for the
$\tau\in\\{1,..,\tau_{max}\\}$ value in group $i$, at time $t$. For the
specific case of this work, $\tau_{max}=28$, as we want to forecast
simultaneously $28$ days ahead. By predicting quantiles, we obtain a quasi-
distribution of the expected value, and gain the capacity to define confidence
intervals.

Table 3: Hyperparameters for TFT model after tuning. hyperparameter | value
---|---
encoder length | 42
batch size | 40
prediction length | 28
gradient clipping | 0.022730
learning rate | 0.0011149
hidden size | 33
number of attention heads | 8
dropout | 0.19230
hidden continuous size | 19

To obtain the quantile predictions, a specific loss, the Quantile Loss, is
defined as (Koenker and Hallock, 2001):

$QL(y,\hat{y},q)=\max\\{q(y-\hat{y}),(q-1)(y-\hat{y})\\}$ (7)

for each quantile $q$. The final Loss is the average QL across quantiles and
for the entire prediction horizon [0,$\tau_{max}$]. In this work, the
quantiles used were [0.2,0.1,0.25,0.5, 0.75,0.9,0.98]. When $q=0.5$ the Loss
is equal to MAE divided by 2, and $q=0.5$ (the median) is the value used for
the point-wise prediction of the model. The complete hyperparamenters of the
model are defined in tab:hyp and the overall architecture of the TFT can be
seen in fig:tft.

##  Appendix B Results

fig:pred

Figure 5: Test sample forecast. Over the input vector, we can see the grey
line representing attention. In orange is the median predictive value (q=0.5),
with different quantiles shown as shaded area.

The two elements of explainability that are given by the model are time-step
attention and variable importance. In Figure LABEL:fig:pred we have a sample
forecast from the test set. In grey are the attention weights given to each
time-step. These values change according to the input, so if we average the
attention for all predictions we can get a sense of the most relevant time-
steps. In Figure LABEL:fig:att, it can be distinguished how the model values
the most recent time-steps with a higher weight, which is intuitively expected
and shows that old information has less value to the model. This validates a
common assumption in linear models, that ascribe more weight to more recent
observations, as is the case of the Exponential Smoothing model. In addition,
we observe another more intriguing feature, already observable in the sample
forecast, which are spikes in attention during the weekend; these spikes may
occur because particular attention is given to one or two previous weekends to
define patient volume in future weekends.

fig:att

Figure 6: Average attention attributed over the input vector. More recent
time-steps have higher values.

The model categorises covariates into three types: static, past, and future.
And, as a consequence of using the Variable Selection Network, features can be
listed by importance. In fig:previsoes,fig:previsoes2 we have, respectively,
the past and future covariates importance. The most relevant past feature is
the percentage of patients in the emergency room with respiratory problems.
For this period, excess affluence in hospital emergency rooms could be
attributed to peaks in influenza/COVID-19 transmission, it therefore makes
sense that this variable can be a predictive indicator of future positive
trends in the number of non-urgent cases. The second most important variable
is patient waiting time, which is in line with the positive relationship
presented in chap:dataa between the increase in waiting time and the increase
in non-urgent patients. However, we should not focus our attention solely on
the variables relevant to the model. There is interest in observing the
variables that did not add value to the model; here we can observe that the
information regarding healthcare centres (n_cons_total, prog) did not add
value to the model, meaning that there is no clear interaction between patient
volume in healthcare centres, mostly used for primary healthcare and minor
health issues, and non-urgent patients in Emergency Departments.

The most important future covariate feature is the categorical variable
indicating public holidays in Portugal. The model has attributed such an
importance to holidays because they have a severe impact on patient volume,
not only on the day, but also on the next day, when close to the weekend.
Furthermore, the other future covariates have a non-negligible importance both
as past and future covariates, thus supporting the claim found in the
literature that calendar variables have a significant impact on the
prediction. Although these elements of explainability certainly add value to
the model, there are still some limitations when interpreting these results.
In the case of the variable importance, it is not possible to ascertain if the
impact of each variable is positive or negative. Furthermore, because these
three networks are separate, it is not possible to compare past with future
variables, nor the total impact of a given input type (static, past or
future).

fig:previsoes

Figure 7: Past Covariates Importance. The most relevant past feature to the
model is the percentage of patients in ED with indication of respiratory
problems.

fig:previsoes2

Figure 8: Future Covariates Importance. For future covariates, variables that
are known in the future, the most relevant one is a feature that indicates
public holidays in Portugal.


# Differentiable Linear Bandit Algorithm

Kaige Yang
University College London
<EMAIL_ADDRESS>
&Laura Toni
University College London
<EMAIL_ADDRESS>

###### Abstract

Upper Confidence Bound (UCB) is among the most commonly used methods for
linear multi-arm bandit problems. While conceptually and computationally
simple, this method highly relies on the confidence bounds, failing to strike
the optimal exploration-exploitation if these bounds are not properly set. In
the literature, confidence bounds are typically derived from concentration
inequalities based on assumptions on the reward distribution, e.g., sub-
Gaussianity. The validity of these assumptions however is unknown in practice.
In this work, we aim at learning the confidence bound in a data-driven
fashion, making it adaptive to the actual problem structure. Noting that
existing UCB-typed algorithms are not differentiable with respect to
confidence bound, we first propose a novel differentiable linear bandit
algorithm. Then, we introduce a gradient estimator, which allows permit to
learn the confidence bound via iterative gradient ascent. Theoretically, we
show that the proposed algorithm achieves a
$\tilde{\mathcal{O}}(\hat{\beta}\sqrt{dT})$ upper bound of $T$-round regret,
where $d$ is the dimension of arm features and $\hat{\beta}$ is the learned
size of the confidence bound. Empirical results show that $\hat{\beta}$ is
significantly smaller than its theoretical upper bound and proposed algorithms
outperform baseline ones on both synthetic and real-world datasets.

##  1 Introduction

Multi-Arm Bandit (MAB) [5] is an online decision making problem, in which an
agent selects arms sequentially and observes stochastic rewards as feedback.
The goal of the agent is to maximize the expected cumulative reward over a
number of trials. The expected reward of each arm is unknown a priori and it
is learned from experience by the agent. As a consequence, the agent needs to
balance the selection of arms to improve its knowledge (exploration) and the
selection of the highest rewarding arm given the knowledge acquired till thus
far (exploitation). This is formalized as the so-called exploration-
exploitation trade-off. Bandit algorithms are designed to strike this trade-
off. One class of MAB problems is the linear MAB [8], in which each arm is
described by a feature vector and the expected reward follows a linear model
over its feature vector and an unknown parameter vector. Each arm’s feature
vector is known a priori by the agent and it is considered as a hint on the
arm reward. The learning problem boils down to the agent inferring the unknown
parameter vector, based on the history (selected arms and received rewards)
and selecting arms accordingly.

One popular algorithm to solve linear MAB is the Upper Confidence Bound (UCB)
[3] [8] [1]. Its popularity is motivated by its conceptual simplicity and
strong theoretical guarantees. UCB-typed algorithms rely on the construction
of an upper confidence bound, which is the estimated reward inflated based on
the level of uncertainty of the estimate. At each decision opportunity, the
agent selects the arm with the highest upper confidence bound. This reflects
the Optimism in Face of Uncertainty principle. In such way, either the arm
with high estimated reward (exploitation) or high uncertainty (exploration) is
selected. However, to properly balance between exploration and exploitation,
it is fundamental to establish a tight confidence bound [16].

In most existing works, confidence bounds are derived from concentration
inequalities [1] [4] [18] given a priori assumptions on the reward
distribution (e.g., sub-Gaussinaity). These bounds achieve strong minimax
theoretical guarantees, outperforming competitor algorithms such as LinTS [2].
While these bounds are essential for a theoretical analysis, they do not
necessarily translate into practice. In fact, these constructed confidence
bounds are typically conservative in practice, as noted in [19] [14]. This is
because concentration inequalities are usually built based on given reward
distributions instead of the actual data (or problem structure). This results
in non-adaptive and potentially wide confidence bounds which in turn lead to
suboptimal performance in practice.

Alternatively, in this work we aim to learn the confidence bound in a data-
driven fashion making it adaptive to the actual problem structure. Inspired by
[6], we aim at having a parametrized and differentiable cumulative reward
function with respect to the confidence bound, which can then be optimized.
The key challenge is that existing UCB-typed algorithms are non-differentiable
with respect to the confidence bound, mainly due to maximization of the UCB
index (i.e., due to the presence of the $\arg\max$ operator in the OFUL [1],
LinUCB [8]). To address this, we propose a novel differentiable UCB-typed
linear bandit algorithm and introduce a gradient estimator which enables the
confidence bound to be learned via gradient ascent.

Our proposed algorithm contains two core components. First, we consider a more
informative UCB-based index than the classical UCB index used in OFUL [1],
LinUCB [8], which not only summarizes the history of each arm but also
differentiates arms to be suboptimal arms and non-suboptimal arms. Second, we
consider a softmax function, which transforms each index into a probability
distribution, where the probability for each suboptimal arm to be selected is
arbitrary small. Conversely, the probability for a non-suboptimal arm to be
selected is greater for arms with larger index. The key idea is that the
exploration is conducted by selecting arms with large index more often than
others. The exploitation is achieved by soft-eliminating suboptimal arms
(arbitrary small probability to be selected). The softmax function ensures the
differentiabilty of the reward function, paving the way to learn confidence
bound via gradient ascent. Based on this, we provide two linear bandit
algorithms for learning confidence bound in both offline and online settings.
Theoretically, we provide a regret upper bound for the offline learning
setting.

In summary, we make the following contributions:

  * •

We propose a novel UCB-typed linear bandit algorithm where the expected
cumulative reward is a differentiable function of the confidence bound.

  * •

We introduce a gradient estimator and show how the confidence bound can be
learned via gradient ascent both in offline/online settings.

  * •

Theoretically, we prove a $\tilde{\mathcal{O}}(\hat{\beta}\sqrt{dT})$ upper
bound of $T$-rounds regret where $\hat{\beta}$ is the learned size of the
confidence bound in the offline setting.

  * •

Empirically, we show $\hat{\beta}$ is significantly smaller than its
theoretical upper bound, leading to substantially lower cumulative regrets
with respect to state-of-the-art baselines on synthetic and real-world
datasets.

Notation: $[K]$ mean the set $\\{1,2,...,K\\}$. Arm is indexed by
$i,j\in\mathcal{A}$. We use boldface lower letter, e.g., $\mathbf{x}$, to
denote vector and boldface upper letter. e.g., $\mathbf{M}$, to denote matrix.
For a positive definite matrix $\mathbf{M}\in\mathbb{R}^{d\times d}$ and a
vector $\mathbf{x}\in\mathbb{R}^{d}$, we denote the weighted 2-norm by
$||\mathbf{x}||_{\mathbf{M}}=\sqrt{\mathbf{x}^{T}\mathbf{M}\mathbf{x}}$. Each
arm $k$ is represented by the feature vector
$\mathbf{x}_{k}\in\mathbb{R}^{d}$. We denote by $\mathbb{P}$ and $\mathbb{E}$
the probability distribution and the expectation operator, respectively.

##  2 Related work

Our work is inspired by [6], which was the first attempt in addressing policy-
gradient optimization of bandit policies via differentiable bandit algorithm.
However, there are fundamental differences between [6] and our work. First,
authors proposed a differential bandit framework for Bayesian MAB problem,
which is not directly applicable to linear MAB problems. Conversely, we
propose a differentiable UCB-typed linear bandit algorithm. Second, the main
goal of [6] is to learn the learning rate (coldness-parameter) of the softmax
function, while our algorithm aims at learning the size of the confidence
bound. Third, we propose algorithms for both offline and online settings,
while [6] covered the offline setting only. Moreover, in [6] a regret analysis
was provided for MAB with two arms. In contrast, we provide a regret analysis
for linear MAB with arbitrary finite number of arms in offline setting.

Another work focused on data-dependent UCB is [12]. Authors proposed an
algorithm called bootstrapedUCB. In [12], the stochastic reward is assumed to
be sub-Weibull random variable. Multiplier bootstrap was employed to
approximate the reward distribution. The boostrapped quantile acted as UCB to
facilitie exploration. Their algorithm was deployed on both MAB and linear MAB
problems, while regret analysis covered MAB only. Similar to this work, other
bootstrap techniques were employed [7] [9] [22]. Although aiming to the same
goal (data-dependent UCB), these works are fundamentally differnt from our
approach. Our algorithm is a differentiable bandit algorithm where we rely on
gradient estimator to learn UCB. Their algorithm is non-differentiable,
relying on the boostrapped quantile of the assumed reward distribution to
construct UCB.

Bootstrap techniques were used also for Thompson Sampling exploration in [19],
in which author proposed the BoostrapThompson algoritihm for MAB. Bootstrap
techniques were used to sample observations from historical and pseudo
observations to approximate the posterior distribution which was then used to
encourage exploration. As an extension, [23] generalized this technique to
Gaussian reward MAB, while [13] and [14] proposed an extension to contextual
linear bandit, achieving the same regret bound of LinTS [2]. The problem they
aimed to address was the computational infeasibility of inferring posterior
distribution when reward follows nonlinear models. This departs from our goal,
which is rather learning the confidence bound from data.

Our work can be viewed as a subtle combination of EXP3 [5] and Phased
Elimination111Algorithm: Phased elimination with G-optimal exploration page.
258 [17] [17]. EXP3 was designed for MAB, where arms with higher empirical
averaged reward are signed with larger probability by softmax function. The
coldness-parameter of softmax function is a tunable hype-parameter chosen by
the user. In our work, we propose a novel scheme to set this parameter
automatically in a data-driven fashion. Moreover, although Exp3 is a
differentiable bandit algorithm, it is not an UCB-typed algorithm. Phased
Elimination eliminates suboptimal arms based on the same index as ours and
selects non-suboptimal arms uniformly (pure exploration). There are several
fundamental differences between this approach and our work: i) the confidence
bound in our work is learned from data and not from concentration inequalities
– leading to a less conservative bound; ii) Phased Elimination is a non-
differentiable algorithm; iii) Phased Elimination achieves optimality in a
worst case scenario (minmax regret) while our algorithm get an empirical gain
being data dependent.

In summary, to the best of our knowledge, our work is the first differentiable
UCB-typed linear bandit algorithm which enables confidence bound to be learned
purely from data without relying on concentration inequalities and assumptions
on the form of reward distribution.

##  3 Problem setting

We consider the stochastic linear bandit with an arm set $\mathcal{A}$ and a
time horizon of $T$-rounds. The arm set contains $K$ arms, i.e.,
$|\mathcal{A}|=K$, where $K$ could be large. Each arm $i\in\mathcal{A}$ is
associated with a known feature vector $\mathbf{x}_{i}\in\mathbb{R}^{d}$. The
expected reward of each arm $\mu_{i}=\mathbf{x}_{i}^{T}\boldsymbol{\theta}$
follows a linear relationship over $\mathbf{x}_{i}$ and an unknown parameter
vector $\boldsymbol{\theta}$. Similarly to other works in the bandit
literature, we assume that arm feature and parameter vector are bounded
$||\mathbf{x}||_{2}\leq L$ and $||\boldsymbol{\theta}||_{2}\leq C$, where
$L>0$ and $C>0$. At the beginning of each decision opportunity $t\in[T]$, the
learning agent selects one arm $i\in\mathcal{A}$ within the arm set
$\mathcal{A}$. Upon this selection, the agent observes the instantaneous
reward $y_{t}\in[0,1]$, which is drawn independently from a distribution with
unknown mean $\mu_{i}=\mathbf{x}_{i}^{T}\boldsymbol{\theta}$. The agent aims
to maximize the expected cumulative reward over the time horizon $T$. Namely,

$Y_{T}=\sum_{t=1}^{T}\mathbb{E}[y_{t}]$ (1)

This is equivalent to minimize the expected cumulative regret which measures
the difference between the expected cumulative reward if the optimal arm were
always selected and the agent’s expected cumulative reward. Denoting by
$\mu_{*}=\max_{i\in\mathcal{A}}\mathbf{x}_{i}^{T}\boldsymbol{\theta}$ the
expected reward of the optimal arm, we get

$R_{T}=T\mu_{*}-\sum_{t=1}^{T}\mathbb{E}[y_{t}]\,.$ (2)

Upper Confidence Bound (UCB). The upper confidence bound algorithm, e.g., OFUL
[1], is designed based on the Optimism in Face of Uncertainty principle. The
key aspect is to construct a confidence bound of the estimated reward of each
arm. Formally, at each round $t$, the confidence bound is defined as

$|\hat{\mu}_{i,t}-\mu_{i}|\leq\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}},\
\ \forall i\in\mathcal{A}$ (3)

where $\hat{\mu}_{i,t}$ is the estimate of the reward of arm $i$ at round $t$
and $\mathbf{V}_{t}=\sum_{t=1}^{T}\mathbf{x}_{t}\mathbf{x}_{t}^{T}$ is the
Gram matrix up to round $t$. Then, the agent selects the arm with the highest
upper confidence bound as follows

$i_{t}=\arg\max_{i\in\mathcal{A}}\hat{\mu}_{i,t}+\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
(4)

It is well known that the tighter the bound in Eq. 3, the better the balance
between exploration and exploitation [16]. Most existing confidence bounds are
established based on concentration inequalities. e.g., Hoeffding inequality
[4], self-normalized [1], Azuma Inequality [17], Bernstein inequality [18]. As
a specific example, under the assumption of the stochastic reward to be a
$R$-sub-Gaussian variable, one of the state-of-the-art high probability upper
bound of $\beta$, derived based on properties of self-normalized martingale,
was given by [1]:

$\beta\leq
R\sqrt{2\log\left(\frac{1}{\delta}\right)+d\log\left(1+\frac{T}{d}\right)}+\sqrt{\alpha}C$
(5)

where $\alpha$ is a regularizer parameter of least-square estimator,
$1-\delta$ is the probability of which Eq. 3 holds and
$||\boldsymbol{\theta}||_{2}\leq C$. The tightness of this (and other bounds)
relies on the validity of assumptions on the reward distribution, which is
unfortunately unknown in practice. Alternatively, we aim at learning the
confidence bound, i.e., $\beta$, in a data-driven fashion without any a priori
assumption on the unknown reward distribution except the linearity function of
the mean reward, i.e., is
$\mu_{i}=\mathbf{x}_{i}^{T}\boldsymbol{\theta},\forall i\in\mathcal{A}$.

##  4 Algorithms

In this section, we first present our proposed algorithm whose expected
cumulative reward is a differentiable function of the confidence bound. Then,
we provide a gradient estimator which enables confidence bound to be learned
via gradient ascent. Next, we propose two algorithms to learn the confidence
bound in offline and online settings, respectively. Finally, we prove a regret
upper bound for offline learning setting.

###  4.1 Differentiable Algorithm

Input : $\beta$, $\mathcal{A}$, $K$, $T$, $\alpha$.

Initialization : $\mathbf{V}_{0}=\alpha\mathbf{I}\in\mathbb{R}^{d\times d}$,
$\mathbf{b}_{0}=\mathbf{0}\in\mathbb{R}^{d}$,
$\hat{\boldsymbol{\theta}}_{0}=\mathbf{0}\in\mathbb{R}^{d}$, $\gamma_{0}=0$.

for  _$t\in[1,T]$_ do

  1. 1.

Find $S_{i,t},\forall i\in\mathcal{A}$ via Eq. 7 with $\beta$.

  2. 2.

Find $\boldsymbol{\pi}_{t}$ via Eq. 8 with $\gamma_{t-1}$.

  3. 3.

Select arm $i_{t}\in\mathcal{A}$ randomly following $\boldsymbol{\pi}_{t}$ and
receive payoff $y_{t}$.

  4. 4.

Update
$\mathbf{V}_{t}\leftarrow\mathbf{V}_{t}+\mathbf{x}_{t}\mathbf{x}_{t}^{T}$,
$\mathbf{b}_{t}\leftarrow\mathbf{b}_{t-1}+\mathbf{x}_{t}y_{t}$ and
$\hat{\boldsymbol{\theta}}_{t}=\mathbf{V}_{t}^{-1}\mathbf{b}_{t}$.

  5. 5.

Update $\gamma_{t}$ via Eq. 9 .

end for

Algorithm 1 SoftUCB

Our proposed algorithm named SoftUCB is shown in Algorithm 2. SoftUCB contains
two core components: an UCB-based index $S_{i,t}$ and an arm selection policy
$\boldsymbol{\pi}_{t}$. Formally, for $i\in\mathcal{A}$,
$\hat{\mu}_{i,t}=\mathbf{x}_{i}^{T}\hat{\boldsymbol{\theta}}_{t}$ where
$\hat{\boldsymbol{\theta}}_{t}=\mathbf{V}^{-1}_{t}\sum_{s=1}^{t}\mathbf{x}_{s}y_{s}$
is the least-square estimator and
$\mathbf{V}^{-1}_{t}=\sum_{s=1}^{t}\mathbf{x}_{s}\mathbf{x}_{s}^{T}$ is the
Gram matrix up to round $t$. Let denote by
$i_{*}=\arg\max_{i\in\mathcal{A}}\hat{\mu}_{i,t}-\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
the arm with the largest lower confidence bound at round $t$. Let us also
define

$\displaystyle\phi_{i,t}=||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i_{*}}||_{\mathbf{V}^{-1}_{t}}\text{\
\ and \ \ }\hat{\Delta}_{i,t}=\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t}$ (6)

where $\beta$ is the confidence bound defined in Eq. 3 and
$\hat{\Delta}_{i,t}$ is the estimated reward gap between $i_{*}$ and $i$.
Equipped with the above notations, we are now ready to introduce the UCB-based
index $S_{i,t}$ defined as

$S_{i,t}=\beta\phi_{i,t}-\hat{\Delta}_{i,t}\,.$ (7)

It is worth noting that $S_{i,t}$ is more informative than classical UCB index
provided Eq. 4, because of the following two key properties: i), $S_{i,t}$
differentiates arms into suboptimal arms and non-suboptimal arms.
Specifically, $S_{i,t}<0$ identifies arms which are suboptimal,
$\Delta_{i}=\mu_{*}-\mu_{i}>0$, and therefore could be eliminated (i.e., not
selected by the agent); ii), $S_{i,t}\geq S_{j,t}\geq 0$ implies that the
upper confidence bound
$\hat{\mu}_{i,t}+\beta|\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}\geq\hat{\mu}_{j,t}+\beta||\mathbf{x}_{j}||_{\mathbf{V}^{-1}_{t}}$
and therefore arm $i$ is more likely to be selected, in line with the Optimism
in Face of Uncertainty principle. These two properties are stated formally in
Lemma 1.

######  Lemma 1.

If $S_{i,t}<0$, arm $i$ is a suboptimal arm, i.e., $\mu_{*}-\mu_{i}>0$. If
$S_{i,t}\geq S_{j,t}\geq 0$, then the upper confidence bound
$\hat{\mu}_{i,t}+\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}\geq\hat{\mu}_{j,t}+\beta||\mathbf{x}_{j}||_{\mathbf{V}^{-1}_{t}}$.
The proof is provided in Appendix A.

We now describe the arm selection strategy. At each round $t\in[T]$, the
probability for arm $i$ to be selected is defined as

$p_{i,t}=\frac{\exp(\gamma_{t}S_{i,t})}{\sum_{j=1}^{K}\exp(\gamma_{t}S_{j,t})}$
(8)

where $\gamma_{t}>0$ is the coldness-parameter controlling the concentration
of the distribution (policy)
$\boldsymbol{\pi}_{t}=[p_{1,t},p_{2,t},...,p_{K,t}]$, and it is set as

$\gamma_{t}=\frac{\log\left(\frac{\delta|\mathcal{L}_{t}|}{1-\delta}\right)}{\tilde{S}_{\text{max},t}}$
(9)

where at each round $t$, the arm set $\mathcal{A}$ is divided into two subsets
$\mathcal{U}_{t}$ and $\mathcal{L}_{t}$ with
$\mathcal{U}_{t}\cup\mathcal{L}_{t}=\mathcal{A}$ and
$\mathcal{U}_{t}\cap\mathcal{L}_{t}=\emptyset$. Namely, $\mathcal{L}_{t}$ is
the set of suboptimal arms (i.e., $i\in\mathcal{L}_{t}$ if $S_{i,t}<0$) and
$\mathcal{U}_{t}$ is the set of non-suboptimal arms (i.e.,
$i\in\mathcal{U}_{t}$ if $S_{i,t}\geq 0$).
$\tilde{S}_{\text{max},t}=\max_{i\in\mathcal{U}_{t}}S_{i,t}$,
$|\mathcal{L}_{t}|$ is the cardinality of $\mathcal{L}_{t}$ and $\delta$ is a
probability hyper-parameter explained in the following Lemma.

######  Lemma 2.

At any round $t\in[T]$, for any $\delta\in(0,1)$, setting
$\gamma_{t}\geq\log(\frac{\delta|\mathcal{L}_{t}|}{1-\delta})/\tilde{S}_{\text{max},t}$
guarantees that
$p_{\mathcal{U}_{t}}=\sum_{i\in\mathcal{U}_{t}}p_{i,t}\geq\delta$ and
$p_{\mathcal{L}_{t}}=\sum_{i\in\mathcal{L}_{t}}p_{i,t}<1-\delta$. The proof is
provided in Appendix B.

According to Lemma 2, Eq. 9 guarantees that suboptimal arms
($i\in\mathcal{L}_{t}$) are selected with an arbitrary small probability
(i.e., $p_{\mathcal{L}_{t}}<1-\delta\approx 0$ when $\delta\approx 1$). This
leads to a soft-elimination of suboptimal arms. Furthermore, a positive
$\gamma_{t}$ guarantees $p_{i,t}\geq p_{j,t}$ if $S_{i,t}\geq S_{j,t}\geq 0$,
$\forall i,j\in\mathcal{U}_{t}$ which obeys the Optimism in Face of
Uncertainty principle.

Overall, SoftUCB (soft-) eliminates suboptimal arms and selects non-suboptimal
arms according to the index in Eq. 7 which favors the selection of arms with
either high estimated reward or high uncertainty.

###  4.2 Gradient Estimator of $\beta$

We now show that the expected cumulative reward of SoftUCB is a differentiable
function over $\beta$ and introduce a gradient estimator. Formally, given the
expected cumulative reward defined in Eq. 1 and SoftUCB described above, we
have the optimization objective defined as

$\begin{split}\max_{\beta}Y_{T}=\max_{\beta}\sum_{t=1}^{T}\mathbb{E}[y_{t}]=\max_{\beta}\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}\mu_{i},\
\ s.t.\
|\mu_{i}-\hat{\mu}_{i,t}|\leq\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}},\ \
\forall i\in\mathcal{A},t\in[T]\end{split}$ (10)

The imposed constraint ensures that
$\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$ is indeed an actual upper
confidence bound (UCB) at any round $t\in[T]$ for any arm $i\in\mathcal{A}$.
Applying the Lagrange multipliers gives the new objective:

$\begin{split}\max_{\beta}\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}\mu_{i}-\eta(|\mu_{i}-\hat{\mu}_{i,t}|-\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}),\
\ s.t.\ \ \eta>0\end{split}$ (11)

The gradient of $\beta$, denoted as $g(\beta)$, can be derived as (proof in
Appendix C):

$g(\beta)=\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}\mu_{i}\left(\gamma_{t}\phi_{i,t}-\frac{\sum_{j=1}^{K}\gamma_{t}\phi_{j,t}\exp(\gamma_{t}S_{j,t})}{\sum_{j=1}^{K}\exp(\gamma_{t}S_{j,t})}\right)+\eta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
(12)

Note that $\mu_{i}$ is unknown in practice and it is therefore replaced by its
empirical estimate $\hat{\mu}_{i,t}$, leading to the following gradient
estimator

$\hat{g}(\beta)=\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}\hat{\mu}_{i,t}\left(\gamma_{t}\phi_{i,t}-\frac{\sum_{j=1}^{K}\gamma_{t}\phi_{j,t}\exp(\gamma_{t}S_{j,t})}{\sum_{j=1}^{K}\exp(\gamma_{t}S_{j,t})}\right)+\eta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
(13)

The gradient estimator $\hat{g}(\beta)$ in Eq. 13 enables $\beta$ to be
learned via gradient ascent. As a stochastic gradient method, under standard
condition of learning rate, e.g., RM [20], it is expected that $\hat{\beta}$
converges to local optimum.

###  4.3 Training Settings

Equipped with the gradient estimator $\hat{g}(\beta)$(Eq. 13), we now show how
to learn $\beta$ in offline and online settings. The corresponding algorithms
named SoftUCB offline and SoftUCB online are presented in Appendix E.

Offline setting. In this setting, multiple $T$-rounds trajectories of the
bandit problem with the same arm set $\mathcal{A}$ are used to train $\beta$,
which is refined after each $T$-rounds trajectory. The key steps are to
initialize $\hat{\beta}_{0}$ and run SoftUCB on $\mathcal{A}$ for $N$ training
trajectories – each trajectory containing $T$-rounds. After each trajectory
$n\in[N]$, update
$\hat{\beta}_{n}\leftarrow\hat{\beta}_{n-1}+\lambda\hat{g}(\beta)$ via Eq. 13
where $\lambda$ is the learning step. At the end of the training, run SoftUCB
on $\mathcal{A}$ with $\hat{\beta}=\hat{\beta}_{N}$.

As a result of the training, the value of $\hat{\beta}$ is optimized in such a
way that it maximizes the expected cumulative reward of arm set $\mathcal{A}$.
Empirically, the $\hat{\beta}$ to which the algorithm converges is substantial
less than its theoretical upper bound Eq. 5. This translates into a
significant regret reduction. In the following subsection we provide a
theoretical regret upper bound of SoftUCB offline

While the above method is fully adaptive to the structure of $\mathcal{A}$, it
provides a burden on the computational complexity. Specifically, the
computational complexity SoftUCB offline is $\mathcal{O}(NKT)$, since we run
SoftUCB $N$ trajectories with $K$ arms and $T$ rounds in each trajectory. This
is much higher than other linear algorithms such as LinUCB [1] and LinTS [2].
To mitigate this issue, we propose SoftUCB online which learns $\beta$ within
one trajectory in an online fashion.

Online setting. In this setting, $\hat{\beta}$ is updated online during one
$T$-rounds trajectory. Specifically, $\hat{\beta}_{0}$ is initialized and
SoftUCB on $\mathcal{A}$ is run for $T$ rounds. At the end of each round
$t\in[T]$, update
$\hat{\beta}_{t}\leftarrow\hat{\beta}_{t-1}+\lambda\hat{g}_{t}(\beta)$ where
$\lambda$ is the learning step and $\hat{g}_{t}(\beta)$ is the gradient
estimator (Eq. 15 defined blow). This reduces the computationally complexity
to $\mathcal{O}(KT)$ since it does not require the $N$-training trajectories,
which is at the same level of OFUL [1], LinUCB [8] and LinTS [2].

In this setting, $Y_{T}=\sum_{t=1}^{T}\mathbb{E}[y_{t}]$, the objective
function we aim at maximizing, is not available before the end of the
trajectory. To obviate to this problem, similarly to policy gradient methods
for non-episodic reinforcement learning problems [21], we update $\hat{\beta}$
to maximizes the average reward per round $\hat{Y}_{t}$. Formally, at each
round $t$, $\hat{Y}_{t}$ consists of two parts: the observed cumulative reward
up to round $t$ and bootstrapped future reward under the current policy
$\boldsymbol{\pi}_{t}=[p_{1,t},p_{2,t},...,p_{K,t}]$. This translates in the
following problem formulation

$\begin{split}\max_{\beta}\hat{Y}_{t}&=\max_{\beta}\left(\sum_{s=1}^{t}\sum_{i=1}^{K}p_{i,s}\hat{\mu}_{i,s}+(T-t)\sum_{i=1}^{K}p_{i,t}\hat{\mu}_{i,t}\right)/T\\\
&s.t.\ \
|\hat{\mu}_{i,t}-\mu_{i,t}|\leq\beta||\mathbf{x}_{i}||_{\mathbf{V}_{t}^{-1}},\forall
i\in\mathcal{A}\end{split}$ (14)

The gradient estimator $\hat{g}_{t}(\beta)$ at round $t$ can be derived as

$\hat{g}_{t}(\beta)=\frac{1}{T}\bigg{(}\sum_{s=1}^{t}\sum_{i=1}^{K}\hat{\mu}_{i,s}\bigtriangledown_{\beta}p_{i,s}+(T-t)\sum_{i=1}^{K}\hat{\mu}_{i,t}\bigtriangledown_{\beta}p_{i,t}+\eta||\mathbf{x}_{i}||_{\mathbf{V}_{t}^{-1}}\bigg{)}$
(15)

It is worth noting that, at the end of trajectory $t=T$, the $\hat{Y}_{t}$
converges to. $Y_{T}$ in the offline setting.

###  4.4 Theoretical Analysis

######  Theorem 1.

Define $\mathbb{E}[r_{t}]=\mathbb{E}[\mu_{*}-\sum_{i=1}^{K}p_{i,t}\mu_{i}]$ be
the expected regret at round $t\in[T]$. Let $\hat{\beta}=\beta_{N}$ be the
confidence bound learned from the offline training setting after $N$
$T$-rounds trajectories. Let assume that $\gamma_{t}$ follows Lemma 2 and
$\delta\approx 1$. The cumulative regret of SoftUCB is bounded as

$R_{T}=\sum_{t=1}^{T}\mathbb{E}[r_{t}]\leq
4\sqrt{2}\hat{\beta}\delta\sqrt{Td\log\left(\alpha+\frac{T}{d}\right)}=\tilde{\mathcal{O}}\left(\hat{\beta}\sqrt{dT\log\left(1+\frac{T}{d}\right)}\right)$
(16)

where $\tilde{\mathcal{O}}(\cdot)$ hides absolute constant. The proof is
contained in Appendix D.

Theorem 1 provides a regret upper bound of SoftUCB in the offline setting. To
compare the regret bound with that of other algorithms, we show $\hat{\beta}$
explicitly in the upper bound. Our regret bound scales with $d$ and $T$ as the
regret bound $\mathcal{O}(\beta\sqrt{dT})$ of existing UCB-typed algorithms,
e.g., OFUL [1], LinUCB [8], Giro [13]. Since we make no assumption on the
reward distribution, we can not derive a theoretical upper bound on
$\hat{\beta}$. However, it is worth to noting that empirical results (in next
section) show that $\hat{\beta}$ is significantly smaller than its theoretical
upper bound Eq. 5. The theoretical analysis for the online setting is left for
future works.

##  5 Experiments

Our experimental evaluation aims to answer the following questions: (1) Does
the learning curve of $\hat{\beta}$ converge in offline and online settings?
(2) Is $\hat{\beta}$ lower than its theoretical counterpart? (3) How do our
proposed algorithms perform compare to baseline ones?

Table 1: The comparison between $\hat{\beta}$ (offline) and theoretical bound
$\tilde{\beta}$

.  $d=5,T=2^{8}$ $d=5,T=2^{9}$ $d=5,T=2^{10}$ $d=10,T=2^{10}$ $d=15,T=2^{10}$
$\hat{\beta}=\textbf{0.5}$ $\hat{\beta}=\textbf{0.6}$
$\hat{\beta}=\textbf{0.9}$ $\hat{\beta}=\textbf{1.1}$
$\hat{\beta}=\textbf{1.2}$ $\tilde{\beta}=2.56$ $\tilde{\beta}=2.66$
$\tilde{\beta}=2.76$ $\tilde{\beta}=3.25$ $\tilde{\beta}=3.61$

(a) $\hat{\beta}$ (offline)

(b) $R_{T}$ (offline)

(c) $\hat{\beta}$ (online)

(d) $R_{T}$ (online)

Figure 1: Learning curves of SoftUCB offline and SoftUCB online

(a) $d=10$

(b) $d=20$

(c) MovieLens

(d) Jester

Figure 2: Performance of algorithms on synthetic and real-world datasets

In synthetic datasets, there are $K=50$ arms with feature vector drawn
uniformly from $[-1,1]$. The dimension of arm feature is set as $d=10,20$. Arm
feature vectors are normalized to be unit vectors. The parameter vector
$\boldsymbol{\theta}$ is generated as a random unit vector. The noise level is
set as $0.5$ and the regularizer parameter is $\alpha=1$. We use two real-
world datasets: Jester [11] and Movielens [15] (see Appendix G for more
details). We compare the proposed algorithms with baseline ones, namely LinUCB
[1], LinTS [2] and $\epsilon$-greedy [21]. The $\beta$ in LinUCB is set as Eq.
5, LinTS follows [2], and $\epsilon=0.05$ in $\epsilon$-greedy.

Fig. 1 depicts the learning curves of $\hat{\beta}$ and the corresponding
$R_{T}$ in both offline and online settings for the synthetic datasets. The
feature dimension $d=10$. In both settings, $\hat{\beta}$ and $R_{T}$ achieve
convergence. Note that in offline setting, $\hat{\beta}$ is optimized to
maxmize the expected cumulative reward Eq. 10, while in online setting,
$\hat{\beta}$ is optimized to maximize the average reward per round Eq. 14.

In Table 1, we compare $\hat{\beta}$ obtained from offline training and its
theoretical suggested $\tilde{\beta}$ given by Eq. 5. Clearly, $\hat{\beta}$
is significantly less than $\tilde{\beta}$ consistently in all cases. This is
because $\hat{\beta}$ is adaptive to the structure of $\mathcal{A}$, while
$\tilde{\beta}$ is derived based on worst-case (minimax analysis). Note that
the value of $\hat{\beta}$ is highly data-dependent. The value we report here
only valid for our experimental data. However, it is reasonable to expect
$\hat{\beta}$ less that $\tilde{\beta}$ in general. The corresponding learning
curves are shown in Appendix F.

In Fig. 2 , the proposed algorithms converge to lower cumulative regret
comparing with baselines. There are two reasons: First, the confidence bound
$\hat{\beta}$ is optimized. Second, the proposed algorithm eliminates (softly)
suboptimal arms which accelerates the rate of convergence. It is worth noting
that the regret of SoftUCB online is large at the initial phase. This is
because at the beginning, when $\gamma_{0}=0,|\mathcal{L}_{t}|=0$, SoftUCB
online selects arms uniformly which results in large regret. Later, when
suboptimal arms are identified, $|\mathcal{L}_{t}|>0$, $\gamma_{t}>0$
according to Eq. 9. Suboptimal arms are soft-eliminated and non-suboptimal
arms are selected following index Eq. 7 which controls the regret.

Finally, during our experiments, we noticed that the convergence of
$\hat{\beta}$ in both offline and online setting is sensitive to the Lagrange
multiplier $\eta$. With large $\eta$, the gradient ascent algorithm fails in
converging, this is because the gradient estimator Eq. 15 is dominated by
$\eta||\mathbf{x}_{i}||_{\mathbf{V}_{t}^{-1}}$. On the other hand, too small
$\eta$ does not ensure the key constraint
$|\hat{\mu}_{i,t}-\mu_{i}|\leq\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$.
This can lead to erroneously eliminating the optimal arm. Therefore, the
hyper-parameter $\eta$ needs to be tuned carefully during experiments.

##  6 Conclusion

We propose SoftUCB, a novel UCB-typed linear bandit algorithm based on an
adaptive confidence bound, resulting in a less conservative algorithm respect
to UCB-typed algorithms with constructed confidence bounds. The key novelty is
to propose an expected cumulative reward which is a differentiable function of
the confidence bound, and derive a gradient estimator, which enables
confidence bound to be learned via gradient ascent. The estimated confidence
bound $\hat{\beta}$ can be updated under offline/online training settings with
the proposed SoftUCB offline and SoftUCB online, respectively. Theoretically,
we provide a $\tilde{\mathcal{O}}(\hat{\beta}\sqrt{dT})$ regret upper bound of
SoftUCB in the offline setting. Empirically, we show that $\hat{\beta}$ is
significantly less that its theoretical counterpart leading to a reduction of
the cumulative regret compared to state-of-the-art baselines.

There are several directions for future work. First, our work can be combined
with meta-learning algorithms, e.g., MAML [10], to learn a confidence bound
which is adaptive to the common structure of a set of bandit tasks. Second, we
believe our work can be generalized to reinforcement learning (RL) tasks where
exploration and exploitation trade-off is a long standing challenge.

##  7 Broader Impact Discussion

Our work is an algorithm for multi-arm bandit (MAB) problem. On the novelty
side, our work automates the exploration in bandit problems. Such algorithm
could be used in recommendation system and clinic trials. On the positive
side, our work could balance the exploration and exploitation trade-off in a
problem dependently way, which might improve the customer satisfaction or
patient’s health care. On the negative side, depending to the deployed
application, the recommended contents might be unsuitable for some users. To
mitigate this issue, domain knowledge might be required to filter the
recommended contents before releasing to users. Regarding the health care
application, expert’s supervision is essential to avoid any potential hazard.

## References

  * [1] Yasin Abbasi-Yadkori, Dávid Pál, and Csaba Szepesvári.  Improved algorithms for linear stochastic bandits.  In Advances in Neural Information Processing Systems, pages 2312–2320, 2011.
  * [2] Shipra Agrawal and Navin Goyal.  Thompson sampling for contextual bandits with linear payoffs.  In International Conference on Machine Learning, pages 127–135, 2013.
  * [3] Peter Auer.  Using confidence bounds for exploitation-exploration trade-offs.  Journal of Machine Learning Research, 3(Nov):397–422, 2002.
  * [4] Peter Auer, Nicolo Cesa-Bianchi, and Paul Fischer.  Finite-time analysis of the multiarmed bandit problem.  Machine learning, 47(2-3):235–256, 2002.
  * [5] Peter Auer, Nicolo Cesa-Bianchi, Yoav Freund, and Robert E Schapire.  Gambling in a rigged casino: The adversarial multi-armed bandit problem.  In Proceedings of IEEE 36th Annual Foundations of Computer Science, pages 322–331. IEEE, 1995.
  * [6] Craig Boutilier, Chih-Wei Hsu, Branislav Kveton, Martin Mladenov, Csaba Szepesvari, and Manzil Zaheer.  Differentiable bandit exploration.  arXiv preprint arXiv:2002.06772, 2020.
  * [7] Richard Y Chen, Szymon Sidor, Pieter Abbeel, and John Schulman.  Ucb exploration via q-ensembles.  arXiv preprint arXiv:1706.01502, 2017.
  * [8] Wei Chu, Lihong Li, Lev Reyzin, and Robert Schapire.  Contextual bandits with linear payoff functions.  In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, pages 208–214, 2011.
  * [9] Adam N Elmachtoub, Ryan McNellis, Sechan Oh, and Marek Petrik.  A practical method for solving contextual bandit problems using decision trees.  arXiv preprint arXiv:1706.04687, 2017.
  * [10] Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 1126–1135. JMLR. org, 2017.
  * [11] Ken Goldberg, Theresa Roeder, Dhruv Gupta, and Chris Perkins.  Eigentaste: A constant time collaborative filtering algorithm.  information retrieval, 4(2):133–151, 2001.
  * [12] Botao Hao, Yasin Abbasi Yadkori, Zheng Wen, and Guang Cheng.  Bootstrapping upper confidence bound.  In Advances in Neural Information Processing Systems, pages 12123–12133, 2019.
  * [13] Branislav Kveton, Csaba Szepesvari, Mohammad Ghavamzadeh, and Craig Boutilier.  Perturbed-history exploration in stochastic multi-armed bandits.  arXiv preprint arXiv:1902.10089, 2019.
  * [14] Branislav Kveton, Csaba Szepesvari, Zheng Wen, Mohammad Ghavamzadeh, and Tor Lattimore.  Garbage in, reward out: Bootstrapping exploration in multi-armed bandits.  arXiv preprint arXiv:1811.05154, 2018.
  * [15] Shyong Lam and Jon Herlocker.  Movielens data sets.  Department of Computer Science and Engineering at the University of Minnesota, 2006.
  * [16] Tor Lattimore and Csaba Szepesvari.  The end of optimism? an asymptotic analysis of finite-armed linear bandits.  arXiv preprint arXiv:1610.04491, 2016.
  * [17] Tor Lattimore and Csaba Szepesvári.  Bandit algorithms.  preprint, 2018.
  * [18] Volodymyr Mnih, Csaba Szepesvári, and Jean-Yves Audibert.  Empirical bernstein stopping.  In Proceedings of the 25th international conference on Machine learning, pages 672–679, 2008.
  * [19] Ian Osband and Benjamin Van Roy.  Bootstrapped thompson sampling and deep exploration.  arXiv preprint arXiv:1507.00300, 2015.
  * [20] Herbert Robbins and Sutton Monro.  A stochastic approximation method.  The annals of mathematical statistics, pages 400–407, 1951.
  * [21] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [22] Liang Tang, Yexi Jiang, Lei Li, Chunqiu Zeng, and Tao Li.  Personalized recommendation via parameter-free contextual bandits.  In Proceedings of the 38th international ACM SIGIR conference on research and development in information retrieval, pages 323–332, 2015.
  * [23] Sharan Vaswani, Branislav Kveton, Zheng Wen, Anup Rao, Mark Schmidt, and Yasin Abbasi-Yadkori.  New insights into bootstrapping for bandits.  arXiv preprint arXiv:1805.09793, 2018.

## Appendidx A

This section contains the proof of Lemma 1.

###### Proof.

Suppose $S_{i,t}<0$, that is

$\beta(||\mathbf{x}_{i_{*}}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})<\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t}$
(17)

Rearrange terms gives

$\hat{\mu}_{i,t}+\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}\leq\hat{\mu}_{i_{*},t}-\beta||\mathbf{x}_{i_{*}}||_{\mathbf{V}^{-1}_{t}}$
(18)

Note that
$|\mu_{i}-\hat{\mu}_{i,t}|\leq\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}},\
\forall i\in\mathcal{A}$. Then,

$\hat{\mu}_{i_{*},t}-\beta||\mathbf{x}_{i_{*}}||_{\mathbf{V}^{-1}_{t}}\leq\mu_{i_{*}}$
(19)

and

$\mu_{i}\leq\hat{\mu}_{i,t}+\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
(20)

Combine together we have

$\mu_{i}\leq\mu_{i_{*}}\leq\mu_{*}$ (21)

Recall by definition
$i_{*}=\arg\max_{i\in\mathcal{A}}\hat{\mu}_{i,t}-\beta||\mathbf{x}_{i}||_{\mathbf{V}_{t}^{-1}}$
is the arm with largest lower upper bound at round $t$. Therefore,
$\Delta_{i}=\mu_{*}-\mu_{i}>0$. In words, arm $i$ is suboptimal.
Suppose $S_{j,t}\geq S_{i,t}\geq 0$

$\beta(||\mathbf{x}_{j_{*}}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{j}||_{\mathbf{V}^{-1}_{t}})-(\hat{\mu}_{j_{*},t}-\hat{\mu}_{j,t})\leq\beta(||\mathbf{x}_{i_{*}}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})-(\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t})$
(22)

Recall the definition of $i_{*}$,

$i_{*}=\arg\max_{j\in[K]}\hat{\mu}_{j,t}-\beta||\mathbf{x}_{j}||_{\mathbf{V}^{-1}_{t}}$
(23)

Thus, at each time $t$, $i_{*}=j_{*}$. Then,

$\beta||\mathbf{x}_{j}||_{\mathbf{V}^{-1}_{t}}+\hat{\mu}_{j,t}\leq\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}+\hat{\mu}_{i,t}$
(24)

∎

## Appendix B

This section contains the proof of Lemma 2.

###### Proof.

$p_{\mathcal{U}_{t}}=\frac{\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})}{\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})+\sum_{j\in\mathcal{L}_{t}}\exp(\gamma_{t}S_{j,t})}$
(25)

By definition, $S_{j,t}<0$, $\forall j\in\mathcal{L}$. Thus,

$\exp(\gamma S_{j,t})<1,\ \forall j\in\mathcal{L}$ (26)

Then,

$\sum_{j\in\mathcal{L}_{t}}\exp(\gamma S_{j,t})<|\mathcal{L}_{t}|$ (27)

Therefore,

$p_{\mathcal{U}_{t}}>\frac{\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})}{\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})+|\mathcal{L}_{t}|}$
(28)

For any probability $\delta\in(0,1)$, we can find a $\gamma_{t}$ such that
$p_{\mathcal{U}_{t}}\geq\delta$, namely

$\frac{\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})}{\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})+|\mathcal{L}_{t}|}\geq\delta$
(29)

Rearrange terms gives

$\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})\geq\frac{\delta|\mathcal{L}_{t}|}{1-\delta}$
(30)

Take logarithm on both sides,

$\log\left(\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})\right)\geq\log\left(\frac{\delta|\mathcal{L}_{t}|}{1-\delta}\right)$
(31)

The left side term is LogSumExp which can be approximated by

$\log\left(\sum_{i\in\mathcal{U}_{t}}\exp(\gamma_{t}S_{i,t})\right)\geq\max_{i\in\mathcal{U}_{t}}\gamma_{t}S_{i,t}=\gamma_{t}\max_{i\in\mathcal{U}_{t}}S_{i,t}$
(32)

Denote $\tilde{S}_{\text{max},t}=\max_{i\in\mathcal{U}_{t}}S_{i,t}$ and let

$\gamma_{t}\tilde{S}_{\text{max},t}\geq\log(\frac{\delta|\mathcal{L}_{t}|}{1-\delta})$
(33)

we have

$\gamma_{t}\geq\frac{\log(\frac{\delta\mathcal{L}_{t}}{1-\delta})}{\tilde{S}_{\text{max},t}}$
(34)

Therefore, if $\gamma_{t}$ satisfies Eq. 34,

$p_{\mathcal{U}_{t}}\geq\delta$ (35)

Clearly, $p_{\mathcal{L}_{t}}<1-\delta$ since
$p_{\mathcal{L_{t}}}+p_{\mathcal{U}_{t}}=1$. ∎

## Appendix C

This section contains the derive of gradients.

###### Proof.

$\begin{split}\max_{\beta}Y(T)&=\max_{\beta}\sum_{t=1}^{T}\mathbb{E}[y_{t}]=\max_{\beta,\gamma}\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}\mu_{i}\\\
&s.t.\ \
|\mu_{i}-\hat{\mu}_{i,t}|-\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}\leq
0,\ \ \forall i\in\mathcal{A},\ \ \forall t\in[T]\end{split}$ (36)

Apply the Lagrange multipliers, the optimization objective is

$\begin{split}\max_{\beta}\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}\mu_{i}-\eta(|\mu_{i}-\hat{\mu}_{i,t}|-\beta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})\
\ s.t.\ \ \eta>0\end{split}$ (37)

Apply the score function
$\bigtriangledown_{\theta}f(\theta)=f(\theta)\bigtriangledown_{\theta}\log
f(\theta)$ to $p_{i,t}$

$\log p_{i,t}=\gamma S_{i,t}-\log\sum_{j=1}^{K}\exp\gamma S_{j,t}$ (38)
$\bigtriangledown_{\beta}\log
p_{i,t}=\gamma_{t}\phi_{i,t}-\frac{\sum_{j=1}^{K}\gamma_{t}\phi_{j,t}\exp\gamma_{t}S_{j,t}}{\sum_{j=1}^{K}\exp\gamma_{t}S_{j,t}}$
(39)

Then, the gradient $g(\beta)$ is

$g(\beta)=\sum_{t=1}^{T}\sum_{i=1}^{K}\mu_{i}p_{i,t}\left(\gamma_{t}\phi_{i,t}-\frac{\sum_{j=1}^{K}\gamma_{t}\phi_{j,t}\exp\gamma_{t}S_{j,t}}{\sum_{j=1}^{K}\exp\gamma_{t}S_{j,t}}\right)+\eta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
(40)

The gradient estimator $\hat{g}(\beta)$ is obtained by repalcing $\mu_{i}$
with $\hat{\mu}_{i,t}=\mathbf{x}_{i}^{T}\hat{\boldsymbol{\theta}}_{t}$ where
$\hat{\boldsymbol{\theta}}_{t}=\mathbf{V^{-1}_{t}}\sum_{s=1}^{t}\mathbf{x}_{s}y_{s}$
is obtained via least-square estimator.

$\hat{g}(\beta)=\sum_{t=1}^{T}\sum_{i=1}^{K}\hat{\mu}_{i,t}p_{i,t}\left(\gamma_{t}\phi_{i,t}-\frac{\sum_{j=1}^{K}\gamma_{t}\phi_{j,t}\exp\gamma_{t}S_{j,t}}{\sum_{j=1}^{K}\exp\gamma_{t}S_{j,t}}\right)+\eta||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
(41)

∎

## Appendix D

This sextion contains the proof of Theorem 1.

###### Proof.

The probability of each arm is defined as

$p_{i,t}=\frac{\exp(\gamma_{t}S_{i,t})}{\sum_{j=1}^{K}\exp(\gamma_{t}S_{j,t})}$
(42)

$S_{i,t}$ is defined as

$S_{i,t}=\hat{\beta}\phi_{i,t}-\hat{\Delta}_{i,t}=\hat{\beta}(||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i_{*}}||_{\mathbf{V}^{-1}_{t}})-(\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t})$
(43)

The cumulative regret to be minimized is defined as

$\begin{split}R_{T}=\sum_{t=1}^{T}\mathbb{E}[r_{t}]=\sum_{t=1}^{T}\mu_{*}-\mathbb{E}[y_{t}]&=\sum_{t=1}^{T}(\mu_{*}-\sum_{i=1}^{K}p_{i,t}\mu_{i})\\\
&=\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}(\mu_{*}-\mu_{i})=\sum_{t=1}^{T}\sum_{i=1}^{K}p_{i,t}\Delta_{i}\end{split}$
(44)

where we use $\sum_{i=1}^{K}p_{i,t}=1$.
At each time $t$, trm set $\mathcal{A}$ is divided into two subsets
$\mathcal{U}_{t}$ and $\mathcal{L}_{t}$ with
$\mathcal{U}_{t}\cup\mathcal{L}_{t}=\mathcal{A}$. Arm $i\in\mathcal{U}_{t}$ if
$S_{i,t}\geq 0$ and arm $i\in\mathcal{L}_{t}$ if $S_{i,t}<0$.

$\mathbb{E}[r_{t}]=\sum_{i=1}^{K}p_{i,t}\Delta_{i}=\sum_{i\in\mathcal{U}_{t}}p_{i,t}\Delta_{i}+\sum_{i\in\mathcal{L}_{t}}p_{i,t}\Delta_{i}$
(45)

Suppose $\gamma_{t}$ follows Lemma 2, then
$\sum_{i\in\mathcal{L}_{t}}p_{i,t}<1-\delta$. Assume $\Delta_{i}\leq 1,\forall
i\in\mathcal{A}$. Then,

$\mathbb{E}[r_{t}]=\sum_{i\in\mathcal{U}_{t}}p_{i,t}\Delta_{i}+\sum_{i\in\mathcal{L}_{t}}p_{i,t}\leq\sum_{i\in\mathcal{U}_{t}}p_{i,t}\Delta_{i}+(1-\delta)$
(46)

By setting $\delta\approx 1$, we have $1-\delta\approx 0$. It means arms in
$\mathcal{L}_{t}$ are unlikely to be selected. So, the second term can be
dropped. Therefore,

$\mathbb{E}[r_{t}]\leq\sum_{i\in\mathcal{U}_{t}}p_{i,t}\Delta_{i}$ (47)

Thus,

$\begin{split}\mathbb{E}[r_{t}]\leq\sum_{i\in\mathcal{U}_{t}}p_{i,t}\Delta_{i}&=\sum_{i\in\mathcal{U}_{t}}p_{i,t}(\mu_{*}-\mu_{i})\end{split}$
(48)

Note that at each time $t$,
$|\hat{\mu}_{i,t}-\mu_{i}|\leq\hat{\beta}||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}},\
\forall i\in[K]$. Then

$\mu_{*}\leq\hat{\mu}_{*,t}+\hat{\beta}||\mathbf{x}_{*}||_{\mathbf{V}^{-1}_{t}}$
(49)

and

$\mu_{i}\geq\hat{\mu}_{i,t}-\hat{\beta}||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$
(50)

Thus,

$\mu_{*}-\mu_{i}\leq\hat{\beta}(||\mathbf{x}_{*}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})+(\hat{\mu}_{*,t}-\hat{\mu}_{i,t})$
(51)

Note that
$\hat{\mu}_{*,t}-\hat{\mu}_{i,t}\leq\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t}$ where
$i_{*}=\arg\max_{j\in[K]}\hat{\mu}_{j,t}-\hat{\mu}_{i,t}$. Therefore,

$\begin{split}\mu_{*}-\mu_{i}\leq\hat{\beta}(||\mathbf{x}_{*}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})+(\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t})\end{split}$
(52)

Since $i\in\mathcal{U}_{t}$, $S_{i,t}\geq 0$. That is
$\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t}\leq\beta(||\mathbf{x}_{*}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})$.
Then,

$\begin{split}\mu_{*}-\mu_{i}&\leq\hat{\beta}(||\mathbf{x}_{*}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})+(\hat{\mu}_{i_{*},t}-\hat{\mu}_{i,t})\\\
&\leq
2\hat{\beta}(||\mathbf{x}_{*}||_{\mathbf{V}^{-1}_{t}}+||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}})\end{split}$
(53)

Define $\psi_{t}=\max_{i\in[K]}||\mathbf{x}_{i}||_{\mathbf{V}^{-1}_{t}}$. We
have

$\mu_{*}-\mu_{i}\leq 4\hat{\beta}\psi_{t}$ (54)

Plugging this into Eq. 48 gives

$\mathbb{E}[r_{t}]\leq 4\hat{\beta}\sum_{i\in\mathcal{U}_{t}}p_{i,t}\psi_{t}$
(55)

Since we assume $\gamma_{t}$ follows Lemma 2, we have
$p_{\mathcal{U}_{t}}=\sum_{i\in\mathcal{U}_{t}}p_{i,t}=\delta$. Therefore,

$\mathbb{E}[r_{t}]\leq
4\hat{\beta}\sum_{i\in\mathcal{U}_{t}}p_{i,t}\psi_{t}=4\hat{\beta}\phi_{t}\sum_{i\in\mathcal{U}_{t}}p_{i,t}=4\hat{\beta}\psi_{t}p_{\mathcal{U}_{t}}\leq
4\hat{\beta}\delta\psi_{t}$ (56)

Thus, the cumulative regret

$\begin{split}R_{T}=\sum_{t=1}^{T}\mathbb{E}[r_{t}]\leq\sqrt{T\sum_{t=1}^{T}\mathbb{E}[r_{t}]^{2}}\leq
4\hat{\beta}\delta\sqrt{T\sum_{t=1}^{T}\psi_{t}^{2}}\end{split}$ (57)

From Lemma 3 (stated below), we have

$\sum_{t=1}^{T}\psi_{t}^{2}\leq 2d\log(\alpha+\frac{T}{d})$ (58)

Plugging in Eq. 57,

$R_{T}\leq
4\hat{\beta}\delta\sqrt{2Td\log(\alpha+\frac{T}{d})}=\tilde{\mathcal{O}}(\hat{\beta}\sqrt{Td\log(1+\frac{T}{d})})$
(59)

where $\delta$ is the probability parameter chosen by user.

######  Lemma 3.

(Lemma 11 in [1])

$\sum_{t=1}^{T}||\mathbf{x}||^{2}_{\mathbf{V}^{-1}_{t}}\leq\log
det(\mathbf{V}_{t})\leq 2d\log(\alpha+\frac{T}{d})$ (60)

∎

## Appendix E

This section contains the pseudo code of SoftUCB, SoftUCB offline and SoftUCB
online.

Input : $\beta$, $\mathcal{A}$, $K$, $T$, $\alpha$.

Initialization : $\mathbf{V}_{0}=\alpha\mathbf{I}\in\mathbb{R}^{d\times d}$,
$\mathbf{b}_{0}=\mathbf{0}\in\mathbb{R}^{d}$,
$\hat{\boldsymbol{\theta}}_{0}=\mathbf{0}\in\mathbb{R}^{d}$, $\gamma_{0}=0$.

for  _$t\in[1,T]$_ do

  1. 1.

Find $S_{i,t},\forall i\in\mathcal{A}$ via Eq. 7 with $\beta$.

  2. 2.

Find $\boldsymbol{\pi}_{t}$ via Eq. 8 with $\gamma_{t-1}$.

  3. 3.

Select arm $i_{t}\in\mathcal{A}$ randomly following $\boldsymbol{\pi}_{t}$ and
receive payoff $y_{t}$.

  4. 4.

Update
$\mathbf{V}_{t}\leftarrow\mathbf{V}_{t}+\mathbf{x}_{t}\mathbf{x}_{t}^{T}$,
$\mathbf{b}_{t}\leftarrow\mathbf{b}_{t-1}+\mathbf{x}_{t}y_{t}$ and
$\hat{\boldsymbol{\theta}}_{t}=\mathbf{V}_{t}^{-1}\mathbf{b}_{t}$.

  5. 5.

Update $\gamma_{t}$ via Eq. 9 .

end for

Algorithm 2 SoftUCB

Input : $\mathcal{A}$, $K$, $T$, $\lambda$, $\eta$

Initialization : $\beta_{0}=0$, $\hat{\beta}=0$.

for  _$n\in[1,N]$_ do

  1. 1.

Run SoftUCB on $\mathcal{A}$ rounds with $\beta=\beta_{n-1}$.

  2. 2.

Update $\beta_{n}\leftarrow\beta_{n-1}+\lambda\hat{g}(\beta)$ via Eq. 13

end for

Output : $\hat{\beta}\leftarrow\beta_{N}$

Run SoftUCB on $\mathcal{A}$ with $\beta=\hat{\beta}$.

Algorithm 3 SoftUCB offline

Input : $\mathcal{A}$, $K$, $T$, $\alpha$, $\lambda$, $\eta$

Initialization : $\beta_{0}=0$, $\mathbf{V}_{0}=\alpha I\in\mathbb{R}^{d\times
d}$, $\mathbf{b}_{0}=\mathbf{0}\in\mathbb{R}^{d}$,
$\hat{\boldsymbol{\theta}}_{0}=\mathbf{0}\in\mathbb{R}^{d}$, $\gamma_{0}=0$.

for  _$t\in[1,T]$_ do

  1. 1.

Select arm $i_{t}\in[K]$ randomly following $\boldsymbol{\pi}_{t}$ and receive
payoff $y_{t}$.

  2. 2.

Update
$\mathbf{V}_{t}\leftarrow\mathbf{V}_{t}+\mathbf{x}_{t}\mathbf{x}_{t}^{T}$,
$\mathbf{b}_{t}\leftarrow\mathbf{b}_{t-1}+\mathbf{x}_{t}\mathbf{y}_{t}$ and
$\hat{\boldsymbol{\theta}}_{t}=\mathbf{V}_{t}^{-1}\mathbf{b}_{t}$.

  3. 3.

Update $\beta_{t}\leftarrow\beta_{t-1}+\lambda\hat{g}_{t}(\beta)$ via Eq. 15.

end for

Algorithm 4 SoftUCB online

## Appendix F

This section contains the learning curves of SoftUCB offline.

(a) $d=5,T=2^{8}$

(b) $d=5,T=2^{9}$

(c) $d=5,T=2^{10}$

Figure 3: Learning curves of SoftUCB offline

(a) $d=10,T=2^{10}$

(b) $d=15,T=2^{10}$

Figure 4: Learning curves of SoftUCB offline

## Appendix G

The dataset Jester contains ratings of 40 jokes from 19891 users. We sample
$K=50$ users randomly as arms. Their rating to top 39 jokes are used as
feature vector. Then, to reduce the sparsity, we apply principle component
analysis algorithm to reduce the dimension $d=10$. Their rating on the 40th
jokes are used as rewards. At each round, the algorithm selects on user to
recommend the joke and the reward is the rating given by the user. MovieLens
contains 6k users and their ratings on 40k movies. Since not every user gives
ratings on all movies, there are a large mount of missing ratings. We
factorize the rating matrix to fill the missing values. The rest works the
same as in Jester.