The Coronavirus disease 2019 (COVID-19) has been impacting the human society since early 2020. At the time of this writing, it is an ongoing public health crisis in over 187 countries and territories around the world, with more than 30 million confirmed cases, and a growing death toll exceeding 1, 000, 000. During this crisis, reliable forecasting of COVID-19 cases becomes important as it will help (1) healthcare institutes to allocate sufficient supply and resources, (2) policy-makers to consider new and further administrative interventions, (3) general public to be aware of the situation and to follow rules against the epidemic. Therefore, the Center for Disease Control and Prevention (CDC) has been actively collecting and publishing data about confirmed cases, hospitalization and deaths related to COVID-19, and hosting forecasting results in the coming weeks.
The US has been suffering the most severe loss from the pandemic, in which more than 200, 000 lives were lost. To encourage and to bring together efforts 1. Compartmental models These models explicitly compartmentalize the population in groups based on their status of infection and recovery, and simulate the transmission process using differential equations. As of today, most of the CDC-featured forecasting methods fall into this category. Examples includes [2, 22, 34] that are built upon classic SIR or SEIR models [11] . Compartmental models describe disease spreading dynamics; however, it is quite hard to determine parameters in these models as they are influenced by many uncontrollable and dynamically changing factors.
Statistical models This type of methods fits the data to regression models directly, such as [1, 21, 31] . While they are more flexible in processing real data compared to compartmental models, they often assume a simplified model class such as generalized linear models [1] , or require sophisticated hand-crafted features from additional, and possibly proprietary, data sources [31] .
The forecasting of COVID-19 is even harder as various constantly changing factors, such as virus characteristics, social and cultural distinctions, public attitudes and behaviors, intervention policies and healthcare preparation, influence the contagious rate and death rate significantly. Will there be a better alternative that is solely data-driven without any assumptions about the underlying disease propagation mechanisms? In particular, we experimented a set leading neural forecasters [17, 18, 25] , but none of them gave the best result. Since the deep models are originally designed for sufficiently long time series with hundreds of points, the scarce historical data in this task might be the reason of their failures. A natural alternative is to exploit other time series in the dataset if they reveals similar dynamics. Fortunately, even if any two regions present different disease curves over long term, it is likely to find short periods in which different regions sharing similar patterns. Figure 1(a) shows surprisingly that the growth pattern of confirmed cases in Santa Barbara county, California, is highly similar to that in Mexico 11 days ago even though at different scales. Moreover, the further growth in Santa Barbara is also close to that within the corresponding time period in Mexico. In light of this key observation, it is intuitively possible to do better forecasting for Santa Barbara by referring to Mexico in this specific time window via proper transformations.
Based on this intuition, we propose to generalize the conventional auto-regressive forecasting to a novel paradigm: besides the local historical data, we also refer to the past reports in all other regions simultaneously in forecasting. Figures 1(b) and (c) illustrate the the fundamental difference between the two paradigms. With time series data of COVID-19 from various locations accumulating over time, we are able to deliver a model outperforming the existing methods by interseries modeling. Note that unlike other cross-location epidemic forecasters such as [8] , only certain time periods rather than the entire time series from other regions will be referred to.
In order to make the proposed paradigm work, it is critical to find small segments in reference time series that exhibit similarity with target time series. It turns out that the attention mechanism originated in natural language process [29] is a good choice for pattern matching. Moreover, it is found that solely applying attention does not work the best as the embedded small segments do not contain long-term trends that are not directly comparable. We filter out these trends and introduce a normalization step so that the small segments can be matched at a consistent scale. In the end, we put all of these components together and achieve global optimum by joint training. Our new model called ACTS (Attention Crossing multiple Time Series), is able to outperform leading forecasters hosted at CDC. Our main contributions are summarized as follows:
• We develop a new paradigm that leverages interseries similarity to improve COVID-19 forecasting. Our method makes no assumption about epidemiological dynamics.
• We extend the attention mechanism to capture inter-series similarity in time series data. Trend filtering is also introduced to complement the attention-based framework and it can be trained jointly to maximize the performance.
• In comparison with a wide range of existing forecasters, the outstanding performance of ACTS is demonstrated on COVID-19 data.
There has been a large body of work focusing on epidemic forecasting. To incorporate domain knowledge, mechanistic models [15, 16, 35] has been favored since they often consider various factors such as epidemiological and social properties, and they make forecasts based on simulation. Moreover, geographic information can also be incorporated into the mechanistic models to better illustrate the spreading process of an infectious disease [3, 4] . These models have excellent interpretability but often fail to fit real observed data due to their rigid and over-simplified assumptions without careful calibration. On the other hand, statistical methods explicitly fit historical data to a statistical model and use it to obtain predictions by extrapolation [5, 6] . For example, [24] relies on kernel density estimation, [19] uses seasonal ARIMA, [33] chooses particle filtering and [36] employs Gaussian process regression. These methods are either too simple or require laborious feature engineering. Hence, various deep learning techniques are also introduced to forecast disease spreading, such as symbol interpretation
The time series from s to t in location i W · Parameter matrices to be learned.
The concatenation of a and b. a, b
Inner product of a and b. s, t
Consecutive index set s, s + 1, · · · , t Table 1 : Used notations [32, 30, 8, 7, 13, 27, 28, 23, 9] . They use deep neural networks to extract complex temporal patterns from historical data and a selected set of additional features. [8, 9] are conceptually closer to our model, both of which employ attention mechanism to compare encoded temporal patterns across multiple locations. However, they require a fixed graph structure with geographic information and produce a similarity score between locations that is independent of time. Instead, in our model we generate embeddings of dynamical patterns for attention over both spatial and temporal dimensions so that the generated attention map are temporally dynamical and free from any predefined geographic structures.
In COVID-19 forecasting, there are three types of incidences, namely confirmed cases, hospitalizations and deaths, to be predicted. The historical data is reported on a daily basis, and we will predict them for the coming weeks. Table 1 summarizes the notations we use in the following sections. Note that throughout the paper, terms "location" and "region" will be used interchangeably. Problem definition is formulated as follows.
We denote by x i t the reported value of a certain type of incidence data at date t and location i, for t = 1, 2, · · · , T and i = 1, 2, · · · , N . Hence, the incidence time series of location i denoted by Problem Statement Given N time series X i 1:T (i ∈ [1, N ]) and additional features, we aim to predict future incidences in a target region i 0 ∈ [1, N ] over H consecutive days after T , i.e. x i0 T +1:T +H .
Traditionally, epidemic forecasts are made by analyzing only the growth pattern of incidences. [3, 4, 8, 9] take the incidences from neighboring regions into consideration as diseases spread through social interaction. Rather than explicitly modeling the disease spreading process, we take a bold step to directly compare the incidence curves across regions. Once the similarities between the current incidences in the target region with the past time segments in reference regions are identified, the following incidences in the reference regions can be used to forecast the future incidences in the target region. Hence, the critical challenge in implementing our idea is to (1) define a representation of a time segment; (2) identify similar segments in reference regions through the representations; and (3) aggregate their following incidences for forecasting. Formally, we introduce an embedding function φ(·) to encode a time series segment x t−l+1:t into a vector, and then use dot-product of vectors to measure similarity. The following incidences x t+1:t+h is also encoded by another embedding function ψ(·) for further aggregation. However, while there are comparable short-term patterns that can be extracted from time series segments, there are also non-stationary long-term trends that hinder reasonable comparison and aggregation of local patterns within segments.
We resolve the problem in two steps. First, we apply a trainable detrending module to the raw time series to remove long-term trends so that incidences across different regions are more comparable. Second, we take rolling windows from residual time series and transform them into a common feature space using normalized convolution as embedding functions φ(·) and ψ(·). The embedding of the recent window in the target region is then compared with windows from references to produce weights for combining the following incidences of each reference window. In such pairwise comparisons, differences in both time-dependent and time-independent features are taken into account so that the curves in corresponding windows can be better aligned. The combi- nations are then added to the extrapolation of filtered trends to generate the final prediction. We jointly train both modules in an end-to-end manner so that both the long-and short-term patterns can be decoupled in an adaptive way. Figure 2 gives an overview of the framework. In the following subsections, we introduce each component in details.
We adopt a learnable Holt smoothing model ( [26] ) to remove long-term trends from the raw time series. Specifically, we introduce a set of parameters θ i e = [a i 0 ; b i 0 ; α i ; β i ] per series, where a i 0 is the initial level, b i 0 is the initial trend, α i is the level smoothing coefficient and β i is the trend smoothing coefficient. Then Holt's equations ( [12] ) are launched to iteratively derive levels and piecewise linear slopes in x i 1:T ,
After detrending, the residual time seriesx 1:T will contain short-term patterns for further processing. Projection from the long-term trend is generated by simple linear extrapolation, (2) x i t+h = a i t + hb i t . A more sophisticated detrending process might further boost performance; we leave it for future study. The detrending process is applied to all the time series and the residual time series are fed into the following attention module.
As COVID-19 is a new disease, we do not have its historical data in the past seasons. Hence, it is critical to leverage limited data from the same season, but across different regions, i.e. model the correlations between regions that have been undergoing the pandemic. Without detailed information about spatial dynamics such as population movement, we instead employ attention mechanism to measure the relation of one region to other regions by directly comparing the incidence curves after trend filtering. Since there are many stages in a dynamical epidemiological process, it is necessary to learn a representation for each time period in a region for alignment in attention. In light of this idea, we apply a convolution layer to encode the residual time series segmentx t−l+1:t to a vector, based on which attention scores measuring similarity between regions are computed.
Even after detrending, the scales of reported numbers in residual time series are still quite different across regions. It is important to normalize residuals before embedding. We empirically find it better to apply min-max normalization to the cumulative sum of incidence time series, which can be regarded as a kind of smoothing. Specifically, for a rolling window of size l representing a period of time, i.e. x i t−l+1:t , t ∈ [l, T ], we compute its cumulative sums and apply the min-max normalization to the monotonically increasing series,
, for j ∈ [t − l + 1, t]. As a result, the first and last values of the normalized series will consistently be 0 and 1 respectively. We then instantiate the function φ(·) using a convolution layer with d feature maps to the scaled segment and time-dependent features. The kernel size is empirically selected, and when it is smaller than l, average pooling is applied in order to reduce a sequence to a vectorized embedding,
These segment embeddings are used to model similarity in different temporal periods across different regions. Likewise, we employ another convolution-pooling layer as ψ(·) to encode the following incidences over H days after each segment into so-called development embedding,
They represent the succeeding development after encoded segments and will be the references for the prediction of the given target region. In fact, we can pair the segments and references by aligning the time indices, i.e.
Given the embeddings, we use dot-product attention to compare segments and combine the values. Specifically, we linearly map the segment embeddings to query vectors q i t and key vectors k t i , from which the similarity score is computed. The development embeddings are projected to value vectors v i t . On the other hand, the additional time-independent features u i are also incorprated into queries and keys. 
For a target region i 0 , we take q i0 T for the last segment and compute its similarity with all the keys from other time segments across all the regions, which is then used to obtain a weighted sum of values.
where Ω = [1, N ] × [l, T − H]. In this way, the past observations in both the target region and reference regions are fully utilized. The weighted combination of valuesv i0 T is then linearly projected to an estimate ofc i0 T +1:T +H . We apply the inverse transformation of (3) to get an estimate ofx i0 T +1:T +H , denoted byŷ i0 T +1:H . In the end, the estimate from attention module is added to the extrapolations in the detrending module to produce the final forecast y i0 T +1:T +H , where
The model can be trained by minimizing the joint loss with respect to the parameters in all the modules. The joint loss is an aggregation of prediction error E(·, ·) computed in two steps. First, for a single target region, we compare our forecasts and ground truths for different T , i.e. lengths of history. Second, we take the aggregated loss in the first step for every region. Formally, the joint loss is defined as
where L is the total number of available historical reports, and l is the minimum required history length.
In our experiments, we choose Mean Absolute Error (MAE) to be the error metric E(·, ·), i.e.
In this section, we demonstrate the effectiveness of the proposed model on real COVID-19 datasets. We intend to answer the following questions:
• Can ACTS outperform the popular COVID-19 forecasters referred at CDC and other state-of-theart deep learning models?
• How much does each component of ACTS contribute to the model performance?
• What kind of similarity can inter-series attention capture?
The COVID-19 incidence data is publicly available at JHU-CSSE 2 and COVID tracking project 3 . Additional features are also publicly available 4 5 6 .
The features we used include total population, population density, ratios of age/gender/race, available hospital beds, and traffic mobility, which are proven to bring marginal accuracy gain in the hospitalization forecasting task in our experiments. The dataset covers the reports up to September 27, 2020 from 50 states and DC in the US.
Evaluation Protocol As required by CDC, we predict the incidence data over the next 4 weeks at a given date and compare the forecasts with the reported ground truths. Suppose we are predicting the new confirmed cases in the state of California starting from 08/16. As context, we are provided a daily time series consisting of incidences in all the states till 08/15. There are three forecasting tasks: daily forecasts for new hospitalizations, weekly forecasts for new confirmed cases and deaths.
The forecasting performance is evaluated in terms of Weighted Absolute Percentage Error (WAPE), defined by the ratio of Mean Absolute Error (MAE) and mean value of ground truths and frequently used in research [17, 18] . At each prediction date, we keep the data in the last 7 days for validation, and the remaining historical data for training. We use the validation data to tune the hyperparameters and to avoid overfitting by early stopping. Other implementation details can be found in Appendix too.
Baselines We compare the performance of the epidemic models featured at CDC, including The first four are compartmental models and the last two rely on statistical modelling. Other than these conventional models, we also evaluate three deep learning models for time series forecasting,
• DeepCOVID [25] An operational deep learning framework designed for real-time COVID-19 forecasting developed by Georgia Tech;
• ConvTrans [17] A self-attention based Transformer model that also employs convolutions for pattern representations;
• TFT [18] A self-attention based deep learning model with feature selection.
We implement the ConvTrans and TFT and tune the hyperparameters using the validation data. All of our implementations run on a server with an Intel i7-6700K CPU and a single GTX 1080Ti GPU. For other baselines, since their implementations are not open-sourced, we take their forecasts submitted to the challenge hosted by CDC 9 . Table 2 shows the forecasting performance on 6 different dates. Three types of incidence data, namely confirmed cases (C), hospitalizations (H) and deaths (D) are separately predicted. We have three key observations: (1) In 13 out of 18 cases, ACTS outperforms other algorithms by a considerable margin. On average, it improves 9%, 5% and 4% over the best of these algorithms for C, H and D, respectively.
(2) ACTS is more favorable on recent days when there are more abundant data available, showing that data-driven methods benefit from more data. (3) The two deep learning approaches ConvTrans and TFT do not exhibit strong performance. The main difference between ours and these approaches is the employment of attention across multiple time series, which dramatically boosts the performance. Note that our model can be trained in less than 5 minutes and inference takes only seconds.
For deeper understanding of our model, we disable each component of ACTS to examine its contribution:
• ACTS-d We remove the detrending module and obtain an attention-only forecaster;
• ACTS-n We remove the normalization in segment embedding; Table 2 : Forecasting performances across different time periods for different types of incidence data in terms of WAPE. A smaller value indicates better performance. We also include the variance of our model's performance by running 5 times with different random initalizations. "-" means the forecasting results of the corresponding baseline are not available. • ACTS-f We remove the additional features in the model and only rely on incidence data.
The hyperparameters of all variants are kept the same. We compare their performance against ACTS using training data up to August 30, 2020. Figure 3 depicts the results, based on which we have the following observations:
• Overall every component of ACTS has positive effects on forecasting accuracy, except that the introduction of additional features has mixed effect. We suspect that either better modelling could help or their effect has been absorbed by the incidence time series;
• Among all the components, inter-series attention has the most significant impact on the performance, which proves that our design of attention crossing multiple time series is valid. It can capture crossregion similarity in COVID-19 forecasting;
• The detrending module makes some contribution. We believe it has the potential for further improvement, e.g. employing advanced trend filtering or even epidemic models. 
A key feature of ACTS is that it can capture similarity between regions via attention from data. According to (7), the reference set Ω is common for any target regions i 0 , and the learned attention distribution is determined by q i0 T . Hence, we directly take those d-dim queries for every region and apply K-means clustering to group them. In this experiment, we use the death forecasting model as an example, where T is August 30, 2020, and K = 4 is selected based on the Elbow method [20] .
A colored map is shown in Figure 4 based on obtained clusters. We can see that California, Texas and Florida, the three states recently hit most seriously are grouped together. Furthermore, states like Arizona, Illinois, North Carolina and Georgia are recognized since they also suffer severe crisis. Interestingly, the states of Wyoming and Vermont are distinguished by our model, in which few deaths are observed for a long period. Overall, our method is able to identify similarities between regions to a certain degree.
In this paper, we present ACTS for COVID-19 forecasting, a purely data-driven framework for an urgent forecasting problem concerning the entire world. It extends the popular deep learning technique, namely attention mechanism, to learning inter-series similarity for time series forecasting. Above that, we also introduce a detrending component to model long-term trends that are difficult for attention model to capture. Both modules are learned jointly based solely on COVID-19 incidence data and a handful of simple features. Without any domain knowledge, our model can empirically outperform many strong forecasters that are featured by CDC. On the other hand, we find great potential for improvement for trend filtering and incorporating additional features, which is left to future work.
We implement our model and its variants using Py-Torch. The hyperparameters we used in all of our experiments are listed in the following 1 6 second. Since for each task we need to forecast 4 weeks, we separately predict each week using the same attention module to avoid long-term errors. To predict the k-th week, we replace (5) by
i.e. take the development (k − 1)H days after the corresponding segment. For case and death forecasting where forecasts are aggregated by weeks, we directly aggregate g i t within a week before applying the final transformation.
In generating final prediction, we avoid negative values by clipping partial predictions from both detrending module and attention module.
Our implementation is open-sourced at https:// github.com/Gandor26/covid-open.
We here show example forecasts of hospitalizations and deaths in three representative states, i.e. Florida, Maryland and Virginia, to give a qualitative demonstration. Our forecasts are accompanied by two best baselines in either task.
We can see that in most cases ACTS fits the ground truths better than baselines. An exception is the hospitalization forecast for Maryland, in which ACTS has systematical underestimation. This is because the downward trend captured by the detrending module significantly drags the final prediction down. It indicates that more advanced trend filtering method can further improve the performance of our model. 