11institutetext: ADVANCED APPLIED AI RESEARCH
PALO ALTO NETWORKS
https://www.paloaltonetworks.com/

# Boosted Embeddings for Time-Series Forecasting

Sankeerth Rao Karingula  Nandini Ramanan  Rasool Tahmasbi  Mehrnaz Amjadi
Deokwoo Jung  Ricky Si  Charanraj Thimmisetty  Luisa Polania Cabrera  Marjorie
Sayer  Claudionor Nunes Coelho Jr

(March 2021)

###### Abstract

Time-series forecasting is a fundamental task emerging from diverse data-
driven applications. Many advanced autoregressive methods such as ARIMA were
used to develop forecasting models. Recently, deep learning based methods such
as DeepAR, NeuralProphet, Seq2Seq have been explored for the time-series
forecasting problem. In this paper, we propose a novel time-series forecast
model, DeepGB. We formulate and implement a variant of gradient boosting
wherein the weak learners are deep neural networks whose weights are
incrementally found in a greedy manner over iterations. In particular, we
develop a new embedding architecture that improves the performance of many
deep learning models on time-series data using a gradient boosting variant. We
demonstrate that our model outperforms existing comparable state-of-the-art
methods using real-world sensor data and public data sets.

###### Keywords:

Time-series Forecasting Deep Learning Gradient Boosting Embedding

##  1 Introduction

Time-series forecasting plays a key role in many business decision-making
scenarios and is one of the central problems in engineering disciplines. In
particular, many prediction problems arising in financial data[32], weather
data[17], econometrics[10] and medical data[18] can be modeled as a time-
series forecasting problem.

Time-series forecasting models can be developed using various autoregressive
(AR) methods. Classical linear models such as autoregressive integrated moving
average (ARIMA) [6] are used to explain the past behavior of a given time-
series data and then used to make predictions of the time-series. ARIMA is one
of the most widely used forecasting methods for univariate time-series data
forecasting. To account for seasonality in time-series data, ARIMA models can
be further extended to seasonal autoregressive integrated moving average
(SARIMA) [6]. In turn, SARIMA models can be extended with covariates or other
regression variables to Seasonal AutoRegressive Integrated Moving Averages
with eXogenous regressors, referred to as SARIMAX model, where the X added to
the end stands for “exogenous”. The exogenous variable can be a time-varying
measurement like the inflation rate, or a categorical variable separating the
different days of the week, or a Boolean representing special festive periods.
One limitation of AR models is that they become impractically slow when
attempting to model long-range dependencies and do not scale well for large
volumes of training data due to the strong assumptions they impose on the
time-series [33].

Recurrent neural networks (RNNs), and in particular Long Short Term Memory
(LSTM) networks, have achieved success in time-series forecasting due to their
ability to capture long-range dependencies and to model nonlinear functions.
Hewamalage et. al. [15] ran an extensive empirical study of the existing RNN
models for forecasting. They concluded that RNNs are capable of modeling
seasonality directly if the series contains homogeneous seasonal patterns;
otherwise, they recommended a deseasonalization step. They demonstrated that
RNN models generally outperform ARIMA models. However, RNN models require more
training data than ARIMA models as they make fewer assumptions about the
structure of the time series and they lack interpretability.

A popular RNN-based architecture is the Sequence to Sequence (Seq2Seq)
architecture [16], which consists of the encoder and the decoder, where both
act as two RNN networks on their own. The encoder uses the encoder state
vectors as an initial state, which is how the decoder gets the information to
generate the output. The decoder learns how to generate target $y[t+1,\ldots]$
by matching the given target $y[\ldots,t]$ to the input sequence. The DeepAR
model for probabilistic forecasting, recently proposed by Salinas et. al.
[27], uses a Seq2Seq architecture for prediction. DeepAR can learn seasonal
behavior and dependencies on given covariates and makes probabilistic
forecasts in the form of Monte Carlo samples with little or no history at all.
Also, it does not assume Gaussian noise and the noise distribution can be
selected by users.

In an effort to combine the best of traditional statistical models and neural
networks, the AR-Net was proposed [33]. It is a network that is as
interpretable as Classic-AR but also scales to long-range dependencies. It
also eliminates the need to know the true order of the AR process since it
automatically selects the important coefficients of the AR process. In terms
of computational complexity with respect to the order of the AR process, it is
only linear for AR-Net, as compared to quadratic for Classic-AR.

Facebook Prophet [30] is another forecasting method, which uses a decomposable
time-series model with three main model components: trend, seasonality, and
holidays. Using time as a regressor, Prophet attempts to fit several linear
and non-linear functions of time as components. Prophet frames the forecasting
problem as a curve-fitting exercise rather than explicitly looking at the
time-based dependence of each observation within a time series. NeuralProphet
[1], which is inspired by Facebook Prophet [30] and AR-Net [33], is a neural
network-based time-series model. It uses Fourier analysis to identify the
seasons in a particular time series and can model trends (autocorrelation
modeling through AR-Net), seasonality (yearly, weekly, and daily), and special
events.

Inspired by the success of learning using Gradient Boosting (GB) [9] and Deep
Neural Networks [11, 35], we present a novel technique for time-series
forecasting called Deep Neural Networks with Gradient Boosting or DeepGB.
Neural networks can represent non-linear regions very well by approximating
them with a combination of linear segments. It is proven by the universal
approximation theorem [7] that neural networks can approximate any non-linear
region but it could potentially need an exponential number of nodes or
variables in order to do this approximation. On the other hand, methods such
as gradient boosting that build decision trees are very good at representing
non-linear regions but they cannot handle complex arithmetic relations, such
as multiplication of elements. The main idea of the proposed approach is to
combine the strengths of gradient boosting and neural networks.

The proposed approach consists of building a series of regressors or
classifiers and solving the problem for the residual at each time. This
approach enables the generation of a number of small parallel models for a
single task, instead of creating a large deep learning model that attempts to
learn very complex boundary regions. Because each subsequent model attempts to
solve only the gradient of the loss function, the task becomes simpler than
attempting to perform a full regression or classification on the output range.
Eventually, the error becomes small enough and indistinguishable from the data
noise.

Our main contributions are as follows.

  * •

We propose DeepGB, an algorithm for learning temporal and non-linear patterns
for time-series forecasting by efficiently combining neural network embedding
and gradient boosting.

  * •

We propose boosted embedding, a computationally efficient embedding method
that learns residuals of time-series data by incrementally freezing embedding
weights over categorical variables.

  * •

In our empirical evaluation, we demonstrate how the proposed approach DeepGB
scales well when applied to standard domains and outperforms state-of-the-art
methods including Seq2Seq and SARIMA in terms of both efficiency and
effectiveness.

The rest of the paper is organized as follows: First, the background on
gradient boosting is introduced in Section 2. Section 3 outlines the need for
neural network embedding for time-series modelling, followed by the detailed
description of the proposed approach, DeepGB. In Section 4, we present the
experimental evaluation of the proposed method. Finally, conclusions and
directions for future work are presented in Section 5.

##  2 Gradient boosting

In Gradient Boosting, the solution is comprised of simpler parallel models
that are trained sequentially but added together at the end. One of the main
ideas of gradient boosting is that each subsequent model (which may be as
small as a small tree or as large as a full deep neural network) needs to only
solve for the residue between the output and the previous regressors built,
thus making it a much easier problem to solve. Gradient Boosting [13] is
motivated by the intuition that finding multiple weak hypotheses to estimate
local probabilistic predictions can be easier than finding a highly accurate
model. Friedman [13] in the seminal work, proposed a gradient boosting
framework to train decision trees in sequence, such that each tree is modeled
by fitting to the gradient/error from the tree in the previous iteration.
Consider for $Data$ = $\\{\langle\mathbf{x},y\rangle\\}$, in the gradient
boosting process, a series of approximate functions $F_{m}$ are learned to
minimize the expected loss $Loss:=E_{x}[L(y,F_{m}(x))]$ in a greedy manner:

$\displaystyle F_{m}:=F_{m-1}+\rho_{m}*\psi_{m}$

where $\rho$ is the step size and $\psi_{m}$ is a tree selected from a series
of $\mathbf{\Psi}$ functions to minimize $L$ in the $m$-th iteration:

$\displaystyle\psi_{m}=\operatorname*{arg\,min}_{\psi\in\mathbf{\Psi}}E[L(y,F^{m-1}(x))+\mathbf{\psi}(x)]$

where the loss function is usually least-squares in most works and a negative
gradient step is used to solve for the minimization function.

In recent years three highly efficient and more successful gradient-based
ensemble methods became popular, namely, XGBoost, LightGBM, and CatBoost [20,
22]. However, both XGBoost and LightGBM suffer from overfitting due to biased
point-wise gradient estimates. That is, gradients at each iteration are
estimated using the same instances that were used by the current model for
learning, leading to a bias. Unbiased Boosting with Categorical Features,
known as CatBoost, is a machine learning algorithm that uses gradient boosting
on decision trees [8]. CatBoost gains significant efficiency in parameter
tuning due to the use of trees that are balanced to predict labels. The
algorithm replaces the gradient step of the traditional gradient boosting with
ordered boosting, leading to reduced bias in the gradient estimation step. It
additionally transforms the categorical features as numerical characteristics
by quantization, i.e., by computing statistics on random permutations of the
data set and clustering the labels into new classes with lower cardinality. In
this paper, we employ CatBoost, which is known to demonstrate effectiveness
from literature, as the gradient boosting algorithm of choice.

In the field of Neural Networks (NNs), He et al.([14]) introduced a deep
Residual Network (ResNet) learning architecture where trained ResNet layers
are fit to the residuals. Although ResNet and gradient boosting are methods
designed for different spaces of problems, there has been significant research
that has gone into formalizing the gradient boosting perspective of ResNet
[24, 34]. Inspired by the complementary success of deep models and gradient
boosting, we in this work propose gradient boosting of NNs, DeepGB, where we
fit a sequence of NNs using gradient boosting. NNs have previously been used
as weak estimators although the methods mostly focused on majority voting for
classification tasks [12] and uniform or weighted averaging for regression
tasks [25, 26]. However, we are the first to investigate gradient boosting of
NNs in the context of time-series forecast in this work.

1:function GradientBoosting($Data$,$M$,$Models$)

2: where $Data$ = $\\{\langle\mathbf{x}_{i},y_{i}\rangle\\}$, $M$: data
samples

3: $models$ = []

4: $F_{0}$ := $y$

5: for $1\leq m\leq|Models|$ do $\triangleright$ Iterate over the Models

6: $model$ = $Models.select()$

7: $model.fit(x,F_{m-1})$ $\triangleright$ Fit the selected model

8: $F_{m}$ = $F_{m-1}-model.predict(x)$ $\triangleright$ Compute the residual

9: if $abs(F_{m}-F_{m-1})<\epsilon$ then $\triangleright$ Check termination
condition

10: break

11: end if

12: $models.append(model)$

13: end for

14: return $models$

15:end function

Algorithm 1 Gradient Boosting Algorithm

Algorithm 1 describes a generic algorithm of gradient boosting. It shows how a
combination of simple models (but not trees) can be composed to create a
better model by just adding them together, effectively assuming each
subsequent model can solve the difference. Interestingly enough, nonlinear
autoregressive moving average model with exogenous inputs - NARMAX [28] has
proposed this approach by clustering increasingly more complex models
together. Note that Algorithm 1 does not make any distinction between how to
select the order of the models. In NARMAX [28], for example, it is suggested
that we should seek simpler models first, and start using more complex models
later on.

##  3 DeepGB Algorithm

In this section, we present our main algorithm DeepGB that is a combination of
gradient boosting with embeddings.

It is often not straightforward to select the best model $F_{1}\cdots F_{m}$
in Algorithm 1 as the models need to account for complex temporal dependency
and non-linear patterns in the time-series. Our approach is to use a neural
network with an embedding layer to learn the underlying residual features from
data with time categorical variables extracted from variable time stamps. The
embedding layers are introduced to learn distinct time-series at once by
embedding their IDs, or to encode categorical features in a meaningful way
(e.g., time-related variables like a month, week of the year, day of the week,
day of the month, holidays, etc.) in a lower dimension and extract valuable
information.

The first modification we will perform in the basic Gradient Boosting
algorithm using general models is to consider

$Models=\bf{list}(Embedding_{1},\ldots,Embedding_{N},Residual),$

where $Embedding_{i}$ is an embedding model that models categorical data, such
as dates, holidays, or user-defined categories, and $Residual$ is a
traditional machine learning model, such as deep neural networks, Gradient
Boosting or Support Vector Machines.

An example of an embedding model can be seen in Figure 1 using Keras. In this
example, $\tt{embedding\\_size}$ can be computed experimentally as
$(\tt{number\\_of\\_categories}+1)^{0.25}$ [31] or as
$\min(50,(\tt{number\\_of\\_categories}+1)//2)$ [3], where $1$ was added to
represent one more feature that is required by the Keras Embedding layer.

categorical_input_shape = (window_size, 1) embedding_size = min(50,
(number_of_categories+1)/ 2) x_inp = Input(categorical_input_shape) x_emb =
Embedding( number_of_categories + 1, embedding_size,
input_length=window_size)(x_inp[…, -1:, 0]) x_out = Dense(1)(x_emb)

Figure 1: Keras model for single output dense layer

A standard gradient boosting approach to train these simple embedding models,
$Embedding_{1},\ldots,Embedding_{N}$, is to select one of those models at each
iteration, train it and then compute its residual error. In the next
iteration, the next selected embedding model is trained on such residual
error.

The key insight of this paper is that we can greatly reduce computational
complexity by freezing embedding weights as we iterate through embedding
training. The motivation behind this iterative freezing is that we preserve
the gradient boosting notion of training weak learners on residuals. Each time
an embedding is learned, we remove the effect of that embedding from the
target data and train the next iteration on the residual. This idea can be
best understood from Figure 2.

Figure 2: Training sequence by freezing of the embedding layers

In Figure 2, we are showing the sequential process in which the different
embedding layers are trained. In particular, at any step $i$, we are
essentially freezing the weights of the embedding layers
$E_{1},E_{2},\ldots,E_{i-1}$ and then training the full model that
concatenates all the $i$ embedding layers and applies the dense layers, D, on
the concatenated output.

####  DeepGB Algorithm description:

Pseudocode for the DeepGB algorithm is presented in Algorithm 2. The steps of
the algorithm are described as follows

  * •

Line 6: We first initiate a sequence of embedding models,
$E_{1},\ldots,E_{N}$.

  * •

Line 7: We specify the input shape of the categorical variable that goes into
the current embedding model.

  * •

Line 8: We concatenate the current embedding with the previously concatenated
embedding models.

  * •

Line 10 - 13: We add dense layers that take the concatenated output of the
embedding models as input.

  * •

Line 15: In the current iteration, the model learns the weights of the current
embedding layer while the previous embedding models are frozen.

  * •

Line 18: We don’t need to learn for all N iteration and stop once the error
converges below a certain threshold $\epsilon$.

  * •

Line 21: The current frozen weights of the embedding model are appended to the
embedding list.

  * •

Line 24: This is where the external model gets trained. In this paper,
CatBoost is used as the external model.

1:function
DeepGB($Data$,$M$,$Models=\left[E_{1},\ldots,E_{N},Residual\right]$)

2: where $Data$ = $\\{\langle\mathbf{x}_{i},y_{i}\rangle\\}$, $M$: data
samples

3: models = [], embeddings = []

4: $F_{0}$ = $y$

5: for $1\leq m\leq N$ do $\triangleright$ Iterate over N embedding layers

6: $c,embedding$ = $Models.select(E_{1},\ldots,E_{N})$ $\triangleright$ Select
an embedding layer at position $c$

7: $\tt{x\\_inp}=\tt{Input(categorical\\_input\\_shape)}$

8:
$\tt{x\\_emb}=\tt{Concatenate}(embeddings+\left[embedding(x\\_inp[\ldots,-1:,c])\right])$

9:$\triangleright$ Adding the next embedding layer

10: for $1\leq d\leq num\\_dense\\_layers$ do

11: $\tt{x\\_emb}=\tt{Dense}(size,activation)(\tt{x\\_emb})$$\triangleright$
Adding dense layers

12: end for

13: $\tt{x\\_out}=\tt{Dense}(1)(\tt{x\\_emb})$

14: $model$ = $\tt{Model}(x\\_inp,x\\_out)$

15: $model.fit(x,y)$

16: Freeze weights of $embedding$

17: $F_{m}$ = $F_{m-1}-model.predict(x)$

18: if $abs(F_{m}-F_{m-1})<\epsilon$ then$\triangleright$ Check termination
condition

19: break

20: end if

21: $embeddings.append(embedding)$

22: end for

23: $models$ = $\left[model\right]$

24: $Residual.fit(x,F_{m})$ $\triangleright$ Fit Residual (Catboost) model on
the error

25: $models.append(Residual)$

26: return $models$

27:end function

Algorithm 2 DeepGB Algorithm

(a) Data set 1

(b) Data set 2

Figure 3: Embedding comparison over data sets

Figure 3 corresponds to plots of the outputs of the algorithm. In particular,
Figure 3 shows embedding layers’ prediction outputs ($model.predict(x)$ in
Algorithm 2) for two different data sets of server CPU utilization. In the
figures, embedding layers are incrementally trained with data input from two
time categorical variables, namely $dayofweek$ (days of a week) and $hour$
(hours of a day). In the first iteration, the embedding layer learns weekly
patterns using $dayofweek$ (i.e., $F_{1}$) as a categorical variable. Then, it
uses the next categorical variable $hour$ (i.e., $F_{2}$) to further improve
the model by learning hourly patterns. The bottom plots shows the final
residual (or irregular) time-series feature after embedding $F_{m}$ in
Algorithm 2.

Figure 3 (a) shows that CPU utilization follows a highly correlated weekly
pattern and that the embedding model accurately captures the weekly and hourly
regular temporal pattern through 1st and 2nd iterations. Note that the final
residual features capture irregular daily pattern in 2020-09-08 (Monday). In
the bottom plot of Figure 3 (b), three irregular hourly patterns in
2020-09-10, 09-11, and 09-23 are accurately captured in the residual feature
and shown as spikes.

##  4 Experimental Evaluation

###  4.1 Datasets

We conducted experiments on the Wikipedia[2] data available to the public via
the Kaggle platform and internal networking device data. The internal data
measures connections per seconds to a device, a rate-based variable which in
general is hard to forecast. Figure 4(a) depicts the time-series for
connections per seconds over a month, from Dec 2020 to Jan 2021. For the sake
of simplicity, we use $I1,\dots,I9$ for representing internal data sets.

The time series in the Wikipedia dataset represents the number of daily views
of a specific Wikipedia article during a certain timeframe. Figure 4(b)
depicts a few time-series plots of the number of accesses of a particular wiki
page starting from July 2015 to September 2019. The Wikipedia dataset also
represents a rate variable, the number of page access per day. Similar to the
internal dataset, we use $P1,\dots,P9$ for representing public data sets.

(a)

(b)

Figure 4: A few sample time-series plots from the internal and public data
sets

###  4.2 Model setup

For benchmarking, DeepGB is compared with SARIMA and Seq2Seq. To keep
comparisons as fair as possible, the following protocol is used: while
employing SARIMA, the implementation in statsmodels [29] is used. For SARIMA,
the best parameters are selected via grid search and SARIMA(p=3, d=1,
q=0)(P=1, D=0, Q=1) is the final model employed for our analysis. For Seq2Seq,
we started with the architecture in Hwang et al. [16] and further tuned to the
setting that gives the best empirical performance. The implementation of the
Seq2Seq model includes a symmetric encoder-decoder architecture with an LSTM
layer of dimension 128, a RepeatVector layer that acts as a bridge between the
encoder and decoder modules, a global soft attention layer [5] for many-to-one
sequence tasks, and a self-attention layer [21]. Also, Time2Vec [19] and
categorical embeddings [35], that enable the use of categorical features based
on the properties of the time-series signature, were added to Seq2Seq. These
enhancements improved the performance of Seq2Seq. The proposed DeepGB model
consists of two parts - layers in the first half of the model implement
embeddings and the later half uses a Boosting model. In particular, CatBoost
with $800$ trees of depth $3$ was used as the gradient boosting model for the
experiments. The proposed DeepGB model consists of 4 embedding layers followed
by a concatenation layer and 4 dense layers with ReLu activation function
after each layer. A dropout layer is also added to prevent overfitting.
Time2Vec [19] and categorical embedding [35] were also used in DeepGB.
Additionally, Root Mean Square Propagation with step size=0.0002 was used as
the optimization algorithm for the experiments. For CatBoost, the
implementation on scikit-learn 0.24.1 with default parameter setting was used.

For all the methods, 30 days and 3 days of data were used for training and
testing, respectively. Only SARIMA used 14 days of data for training, as its
time for convergence increases significantly with the length of data. The
widely used metric for time-series error analysis, Symmetric mean absolute
percentage error (SMAPE), is employed to evaluate forecast accuracy in the
experiments [4]. SMAPE is the main metric used in the M3 forecasting
competition [23], as it computes the relative error measure that enables
comparison between different time series, making it suitable for our setting.

###  4.3 Results

The results of our experiments for the task of time-series forecasting are
summarized in Tables 1 and 2 for the Wikipedia and internal data sets,
respectively. Table 1 shows SMAPE results for DeepGB, Seq2Seq and SARIMA. Bold
font indicates the lowest SMAPE for the given dataset. It can be observed that
DeepGB is significantly better than SARIMA and Seq2Seq in the majority of the
cases. These results help in stating that DeepGB is on par or significantly
better than state-of-the-art approaches, in real domains, at scale. The SMAPE
scores presented in Table 2 indicate that DeepGB outperforms the other two
methods in 6 out of 8 data sets. A closer inspection of the results suggests
that the efficiency of DeepGB is significantly higher (about $3$ times faster)
for dataset P1.

In general, it can be observed that deep models take more time than
statistical methods. However, DeepGB is significantly faster than the other
deep model, Seq2Seq, which is reported to have the second-best performance
after DeepGB. The results indicate that the proposed DeepGB approach does not
sacrifice efficiency (time) for effectiveness (error percentage).

Table 1: SMAPE(lower $\Rightarrow$ better) and training time(lower $\Rightarrow$ better) for DeepGB, Seq2Seq and SARIMA methods on the internal dataset. Dataset | Error (SMAPE) | Training time(secs)
---|---|---
DeepGB | Seq2seq | SARIMA | DeepGB | Seq2Seq | SARIMA
I1 | 1.04 | 6.12 | 1.44 | 10.38 | 85.80 | 3.19
I2 | 1.75 | 2.00 | 1.10 | 9.40 | 86.48 | 6.32
I3 | 5.44 | 21.95 | 20.77 | 8.66 | 86.63 | 5.99
I4 | 2.71 | 6.92 | 29.59 | 8.28 | 87.63 | 2.49
I5 | 6.33 | 6.11 | 6.58 | 8.74 | 81.76 | 5.64
I6 | 4.59 | 9.39 | 11.01 | 8.16 | 82.28 | 3.10
I7 | 6.98 | 19.03 | 17.61 | 11.72 | 193.54 | 5.00
I8 | 6.81 | 4.45 | 17.61 | 13.35 | 195.19 | 4.95
I9 | 61.08 | 62.97 | 59.14 | 12.95 | 196.26 | 3.37
Table 2: SMAPE(lower $\Rightarrow$ better) and training time(lower $\Rightarrow$ better) for DeepGB, Seq2Seq and SARIMA methods on the public dataset (wikipedia.org_all-access_spider). Dataset | Original Series | Error (SMAPE) | Training time(secs)
---|---|---|---
DeepGB | Seq2Seq | SARIMA | DeepGB | Seq2Seq | SARIMA
P1 | 2NE1_zh | 7.98 | 7.93 | 16.12 | 5.86 | 17.04 | 0.22
P2 | 3C_zh | 16.85 | 6.11 | 15.43 | 6.40 | 20.71 | 0.36
P3 | 4minute_zh | 2.34 | 4.38 | 5.90 | 6.22 | 20.36 | 0.36
P4 | 5566_zh | 4.39 | 8.01 | 12.44 | 6.05 | 20.73 | 0.30
P5 | AND_zh | 5.36 | 13.16 | 25.86 | 8.29 | 24.51 | 0.29
P6 | AKB48_zh | 2.27 | 5.81 | 6.78 | 6.30 | 19.48 | 0.34
P7 | ASCII_zh | 2.64 | 8.90 | 7.93 | 7.46 | 20.53 | 0.23
P8 | Ahq_e-Sports_Club_zh | 3.44 | 5.02 | 12.54 | 6.26 | 21.07 | 0.27

##  5 Conclusion

Time-series forecasting is a central problem in machine learning that can be
applied to several real-world scenarios, such as financial forecasts and
weather forecasts. To the best of our knowledge, this is the first work
employing gradient boosting of deep models, infused with embeddings, in the
context of time-series forecasting. To validate the performance of the
proposed approach, we evaluated DeepGB on a public Wikipedia data set and an
internal networking device data set. The experimental results showed that
DeepGB outperforms SARIMA and Seq2Seq using SMAPE as performance metric. It
was also shown in the empirical evaluations that the proposed method scales
well when applied to standard domains as it offers a faster convergence rate
when compared to other deep learning techniques.

Finally, this paper opens up several new directions for further research.
Extending the proposed approach to multivariate time-series forecasting and
deriving theoretical bounds and convergence properties for DeepGB, remain open
problems and are interesting to us from a practical application standpoint.

## References

  * [1] Neural prophet, https://github.com/ourownstory/neural_prophet
  * [2] Wikipedia web traffic time series forecasting, https://www.kaggle.com/c/web-traffic-time-series-forecasting/
  * [3] Arat, M.M.: How to use embedding layer and other feature columns together in a network using keras? (2019), https://mmuratarat.github.io/2019-06-12/embeddings-with-numeric-variables-Keras
  * [4] Armstrong, J.S.: Long-range forecasting. Wiley (1985)
  * [5] Bahdanau, D., Chorowski, J., Serdyuk, D., Brakel, P., Bengio, Y.: End-to-end attention-based large vocabulary speech recognition. In: ICASSP. pp. 4945–4949 (2016)
  * [6] Box, G., Jenkins, G.M.: Time Series Analysis: Forecasting and Control. Holden-Day (1976)
  * [7] Csáji, B.C.: Approximation with artificial neural networks. Faculty of Sciences; Eötvös Loránd University, Hungary (2001)
  * [8] Dorogush, A.V., Ershov, V., Gulin, A.: Catboost: gradient boosting with categorical features support. arXiv:1810.11363 (2018)
  * [9] Friedman, J.: Stochastic gradient boosting. Computational Statistics & Data Analysis 38, 367–378 (2002)
  * [10] Fuleky, P.: Macroeconomic forecasting in the era of big data: Theory and practice, vol. 52. Springer (2019)
  * [11] Goodfellow, I., Bengio, Y., Courville, A.: Deep Learning. The MIT Press (2016)
  * [12] Hansen, L.K., Salamon, P.: Neural network ensembles. IEEE transactions on pattern analysis and machine intelligence (1990)
  * [13] Hastie, T., Tibshirani, R., Friedman, J.: Boosting and additive trees. In: The elements of statistical learning, pp. 337–387. Springer (2009)
  * [14] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: CVPR (2016)
  * [15] Hewamalage, H., Bergmeir, C., Bandara, K.: Recurrent neural networks for time series forecasting: Current status and future directions. International Journal of Forecasting 37(1), 388–427 (2021)
  * [16] Hwang, S., Jeon, G., Jeong, J., Lee, J.: A novel time series based seq2seq model for temperature prediction in firing furnace process. Procedia Computer Science 155, 19–26 (2019)
  * [17] Karevan, Z., Suykens, J.A.: Transductive LSTM for time-series prediction: An application to weather forecasting. Neural Networks 125, 1–9 (2020)
  * [18] Kaushik, S., Choudhury, A., Sheron, P.K., Dasgupta, N., Natarajan, S., Pickett, L.A., Dutt, V.: AI in Healthcare: Time-Series Forecasting Using Statistical, Neural, and Ensemble Architectures. Frontiers in Big Data 3, 4 (2020)
  * [19] Kazemi, S.M., Goel, R., Eghbali, S., Ramanan, J., Sahota, J., Thakur, S., Wu, S., Smyth, C., Poupart, P., Brubaker, M.: Time2vec: Learning a vector representation of time (2019)
  * [20] Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., Liu, T.: Lightgbm: A highly efficient gradient boosting decision tree. In: Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA. pp. 3146–3154 (2017)
  * [21] Lin, Z., Feng, M., Santos, C.N.d., Yu, M., Xiang, B., Zhou, B., Bengio, Y.: A structured self-attentive sentence embedding. arXiv:1703.03130 (2017)
  * [22] Makridakis, S., Spiliotis, E., Assimakopoulos, V.: The m5 accuracy competition: Results, findings and conclusions. Int J Forecast (2020)
  * [23] Makridakis, S., Hibon, M.: The m3-competition: results, conclusions and implications. International Journal of Forecasting 16(4), 451–476 (2000). https://doi.org/https://doi.org/10.1016/S0169-2070(00)00057-1, https://www.sciencedirect.com/science/article/pii/S0169207000000571, the M3- Competition
  * [24] Nitanda, A., Suzuki, T.: Functional gradient boosting based on residual network perception. In: International Conference on Machine Learning. pp. 3819–3828. PMLR (2018)
  * [25] Opitz, D.W., Shavlik, J.W.: Actively searching for an effective neural network ensemble. Connection Science (1996)
  * [26] Perrone, M.P., Cooper, L.N.: When networks disagree: Ensemble methods for hybrid neural networks. Tech. rep., Brown University Institute for Brain and Neural Systems (1992)
  * [27] Salinas, D., Flunkert, V., Gasthaus, J., Januschowski, T.: DeepAR: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting 36(3), 1181–1191 (2020)
  * [28] S.Chen, S.A.Billings: Representations of non-linear systems: the NARMAX model. International Journal of Control 49(3), 1013–1032 (1989)
  * [29] Seabold, S., Perktold, J.: Statsmodels: Econometric and statistical modeling with Python. In: 9th Python in Science Conference (2010)
  * [30] Taylor, S.J., Letham, B.: Forecasting at scale. The American Statistician 72(1), 37–45 (2018)
  * [31] Team, T.: Introducing tensorflow feature columns (2017), https://developers.googleblog.com/2017/11/introducing-tensorflow-feature-columns.html
  * [32] Timmermann, A.: Forecasting methods in finance. Annual Review of Financial Economics 10, 449–479 (2018)
  * [33] Triebe, O., Laptev, N., Rajagopal, R.: AR-Net: A simple auto-regressive neural network for time-series. arXiv:1911.12436 (2019)
  * [34] Veit, A., Wilber, M.J., Belongie, S.: Residual networks behave like ensembles of relatively shallow networks. NIPS 29, 550–558 (2016)
  * [35] Wen, Y., Wang, J., Chen, T., Zhang, W.: Cat2vec: Learning distributed representation of multi-field categorical data (2016), http://openreview.net/pdf?id=HyNxRZ9xg


Study of the CP property of the Higgs boson to electroweak boson coupling in
the VBF $H\to\gamma\gamma$ channel with the ATLAS detector A test of CP
invariance in Higgs boson production via vector-boson fusion has been
performed in the $H\rightarrow\gamma\gamma$ channel using 139 fb-1 of
proton–proton collision data at $\sqrt{s}=13\,\mathrm{TeV}$ collected by the
ATLAS detector at the LHC. The Optimal Observable method is used to probe the
CP structure of interactions between the Higgs boson and electroweak gauge
bosons, as described by an effective field theory. No sign of CP violation is
observed in data. Constraints are set on the parameters describing the
strength of the CP-odd component in the coupling between the Higgs boson and
the electroweak gauge bosons in two effective field theory bases: $\tilde{d}$
in the HISZ basis and $c_{H\tilde{W}}$ in the Warsaw basis. The results
presented are the most stringent constraints on CP violation in the coupling
between Higgs and weak bosons. The constraints on $\tilde{d}$ are further
tightened through combination with results from the $H\rightarrow\tau\tau$
channel. HIGG-2020-08 CERN-EP-2022-134 Phys. Rev. Lett.

The violation of the charge-conjugation and parity (CP) symmetry is one of the
three Sakharov conditions [1] needed to explain the observed baryon asymmetry
of the universe. The only established CP violation source is the complex phase
in the quark mixing matrix [2], from which the derived magnitude of CP
violation in the early universe is insufficient to explain the observed value
of the baryon asymmetry [3, 4, 5]. The discovery of the Higgs boson by the
ATLAS and the CMS experiments [6, 7] at the Large Hadron Collider (LHC) [8]
opened a new direction to search for sources of CP violation: the interactions
of the Higgs boson. The Standard Model (SM) Higgs boson ($H$) is even under
simultaneous charge-conjugation and parity inversion. Any deviation from a
pure CP-even interaction of the Higgs boson with other SM particles could be a
new source of CP violation and also a direct indication of physics beyond the
SM (BSM). The CP structure of Higgs boson couplings to electroweak gauge
bosons and fermions has been studied extensively by the ATLAS and the CMS
experiments [9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. The results are
consistent with the SM prediction, and no sign of CP violation has been found
yet.

A CP-odd component in the Higgs boson coupling to electroweak bosons
($HVV,V=W/Z$) can be described by adding dimension-6 operators to the SM
Lagrangian, using an effective field theory (EFT) approach. The total matrix
element ($\mathscr{M}$) can be written as

$\displaystyle|\mathscr{M}|^{2}=$
$\displaystyle|\mathscr{M}_{\mathrm{SM}}|^{2}+2\cdot
c_{i}\cdot\mathrm{Re}(\mathscr{M}^{\ast}_{\mathrm{SM}}\mathscr{M}_{\textrm{CP-
odd}})$ (1) $\displaystyle+c^{2}_{i}\cdot|\mathscr{M}_{\textrm{CP-odd}}|^{2}.$

The first term describes the SM contribution. The second term (interference
term) is CP-odd, representing a new source of CP violation in Higgs boson
couplings, and is parameterized by the Wilson coefficient $c_{i}$. The third
term (quadratic term) describes a CP-even BSM contribution parameterized by
$c_{i}^{2}$. The interference term only affects CP-odd observables and does
not contribute to CP-even observables, e.g. the inclusive cross-section [19].

Several methods were developed to construct CP-odd observables that can
distinguish CP violation contributions, e.g. in Refs. [12, 17]. This study
adopts the Optimal Observable [20, 21, 22, 23, 24] defined as
$\mathit{OO}=2\cdot\mathrm{Re(\mathscr{M}^{*}_{\mathrm{SM}}\cdot\mathscr{M}_{\mathrm{CP-
odd}})}/{|\mathscr{M}_{\mathrm{SM}}|^{2}}$ to test the CP structure of the
Higgs boson coupling to electroweak bosons in vector-boson-fusion (VBF)
production. This method combines event-based information from a
multidimensional phase space into a single CP-sensitive observable. In the SM,
the $\mathit{OO}$ distribution is expected to be symmetric with a mean value
of zero, and any asymmetrical effects would indicate contributions from the CP
violation term, in the absence of rescattering by new light particles in loops
[25]. For a given event, the matrix elements in the $\mathit{OO}$ definition
are calculated using the four-momenta of the Higgs boson and the two forward
VBF jets, and have no dependence on the decay mode of the Higgs boson. This
method was first introduced in the $H\rightarrow\tau\tau$ analysis [9] by
ATLAS and can be used in all Higgs boson decay channels.

This Letter reports an analysis to test the CP invariance of the $HVV$
coupling by using the Optimal Observable method in the VBF $H\to\gamma\gamma$
channel, using the 139 fb-1 of proton–proton ($pp$) collision data at
$\sqrt{s}=13\leavevmode\nobreak\ \text{TeV}$ recorded during 2015–2018 with
the ATLAS detector. The VBF signal yield in $\mathit{OO}$ bins is extracted
from a simultaneous fit to the diphoton invariant mass spectra split into the
$\mathit{OO}$ bins, which is then used to determine the CP violation
contributions to the $HVV$ coupling.

Results are interpreted in two EFT bases: the HISZ [26] and Warsaw [27, 28,
29] bases. The HISZ basis is used in order to combine the results with the
previous measurement from the $H\to\tau\tau$ channel [9], whereas the Warsaw
basis is used to provide measurements for future combinations with other Higgs
boson measurements. In both bases, three Wilson coefficients multiplying CP-
odd operators describe possible CP-odd couplings between the Higgs boson and
electroweak gauge bosons. In the HISZ basis, $\tilde{d}$ is constrained by
assuming $\tilde{d}=\tilde{d}_{B}$ and setting the third coefficient to zero,
as in Ref. [9]. In the Warsaw basis, $c_{H\tilde{W}}$ is constrained by
setting $c_{H\tilde{B}}$ and $c_{H\tilde{W}B}$ to zero. In both bases, all CP-
even operators coefficients are set to zero. Constraints on all three
coefficients in the Warsaw basis were obtained previously in the $H\rightarrow
ZZ$ channel [13, 16] and $H\rightarrow\gamma\gamma$ channel using differential
cross-sections [12]. The measurements have significant correlations since
these channels cannot distinguish between the three operators. The VBF
topology in this analysis is mainly sensitive to $c_{H\tilde{W}}$ and could
help to reduce this correlation.

The ATLAS detector [30, 31, 32] is a multipurpose particle detector with a
forward-backward symmetric cylindrical geometry and near 4$\pi$ coverage in
solid angle.111ATLAS uses a right-handed coordinate system with its origin at
the nominal interaction point (IP) in the center of the detector and the
$z$-axis along the beam pipe. The $x$-axis points from the IP to the center of
the LHC ring, and the $y$-axis points upwards. Cylindrical coordinates
($r,\phi$) are used in the transverse plane, $\phi$ being the azimuthal angle
around the z-axis. The pseudorapidity is defined in terms of the polar angle
$\theta$ as $\eta$ = $\mathrm{-ln\leavevmode\nobreak\ tan(\theta/2)}$. Angular
distance is measured in units of $\Delta
R\equiv\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$. The trigger system consists
of a hardware-based first-level trigger and a software-based high-level
trigger [33]. Events used in this analysis were accepted by a diphoton trigger
requiring the leading and subleading photons to have transverse energies
($E_{\mathrm{T}}$) greater than 35 GeV and 25 GeV, respectively, during the
whole data-taking period. This trigger had a Loose photon identification
requirement in 2015–2016 [34], but due to the increasing instantaneous
luminosity the identification requirement was tightened for data-taking in
2017–2018 [34]. In addition, a single-photon trigger with Loose identification
criteria and an $E_{\mathrm{T}}$ threshold of 120 (140) GeV in 2015–2016
(2017–2018) was used to recover events with collimated diphoton pairs with
very high transverse momentum ($p_{\text{T}}$) [34]. The average trigger
efficiency is over 98% for events passing the full diphoton event selection
for this analysis [34]. An extensive software suite [35] is used in the
reconstruction and analysis of real and simulated data, in detector
operations, and in the trigger and data acquisition systems of the experiment.

Higgs boson production via VBF was simulated with Powheg Box v2 [36] using the
PDF4LHC15nlo [37] parton distribution function (PDF) set. The generation is
accurate to next-to-leading-order (NLO) in QCD, and the total cross-section is
normalized to a calculation including QCD corrections at full NLO and
approximate next-to-next-to-leading-order (NNLO) accuracy as well as
electroweak (EW) corrections at full NLO accuracy [38, 39, 40]. Higgs boson
production via gluon-gluon fusion (ggF) was modeled at NNLO accuracy in QCD
using Powheg Box v2 [41, 42] and the NNLO family of PDF4LHC15 PDFs. The
simulation achieves NNLO accuracy for arbitrary inclusive $gg\to H$
observables by reweighting the Higgs boson rapidity spectrum in Hj-MiNLO [43,
44, 45] to that in HNNLO [46], and the total cross-section is normalized to a
prediction calculated at next-to-next-to-next-to-leading-order (N3LO) accuracy
in QCD and including NLO EW corrections [47, 48, 49, 50, 51, 52, 53, 54, 55,
56, 57, 57]. Other Higgs boson production processes, e.g. in association with
a vector boson ($VH$) or top quark(s) ($ttH$, $tH$), were also modeled using
Powheg Box v2. Prompt diphoton production ($\gamma\gamma$) was simulated with
the Sherpa 2.2.4 [58] generator. More details can be found in Ref. [12].

To simulate the effects of nonzero values of $\tilde{d}$ and $c_{H\tilde{W}}$
in the $HVV$ vertex, a reweighting method is implemented for the HISZ basis
and Warsaw basis, respectively, and applied to the aforementioned SM VBF
signal sample. For the $\tilde{d}$ coefficient in the HISZ basis, as detailed
in Ref. [9], two weights are calculated by the HAWK program [38, 39, 59] for
each event using generator-level information:
$w_{1}=2\cdot\mathrm{Re(\mathscr{M}^{*}_{\mathrm{SM}}\cdot\mathscr{M}_{\textrm{CP-
odd}})}/{|\mathscr{M}_{\mathrm{SM}}|^{2}}$ and
$w_{2}=\mathrm{|\mathscr{M}_{\textrm{CP-
odd}}|^{2}}/{|\mathscr{M}_{\mathrm{SM}}|^{2}}$ with a specific amount of CP
mixing (given in terms of $\tilde{d}$), to model the contribution from the
interference term and the quadratic term, respectively, as shown in Eq. (1).
For the interpretation in the Warsaw basis, a reweighting of the reconstructed
$\mathit{OO}$ distribution at different values of $c_{H\tilde{W}}$ is obtained
by

$\frac{{d}\sigma}{{d}\mathit{OO}}=\left(\frac{{d}\sigma}{{d}\mathit{OO}}\right)^{\mathrm{NLO}}\cdot\left(\frac{{d}\sigma}{{d}\mathit{OO}}\right)_{c_{H\tilde{W}}}^{\mathrm{MG5}}/\left(\frac{{d}\sigma}{{d}\mathit{OO}}\right)_{c_{H\tilde{W}}=0}^{\mathrm{MG5}},$

where ‘MG5’ labels the prediction from MadGraph [60, 61] using SMEFTSim [27,
28], and ‘NLO’ labels the aforementioned SM VBF signal sample. MadGraph events
for nonzero values of $c_{H\tilde{W}}$ were generated setting the scale of new
physics $\Lambda=1\leavevmode\nobreak\ \text{TeV}$ and fixing all other Wilson
coefficients to zero. For both interpretations, higher-order QCD and
electroweak corrections are assumed to factorize from the new-physics effects.
Limits in the two bases are extracted from the effect of the interference-only
term and also from the effect of the interference-plus-quadratic terms. The
$\mathit{OO}$ value is calculated using HAWK because the $HWW$ operators in
the two EFT bases are similar. HAWK uses the HISZ basis assuming
$\tilde{d}=\tilde{d}_{B}$, which corresponds to
$c_{H\tilde{W}}=c_{H\tilde{B}}$ for the Warsaw basis. However, since
$c_{H\tilde{B}}$ has negligible impact on VBF, only $c_{H\tilde{W}}$ is varied
(setting $c_{H\tilde{B}}=0$) and the computed $\mathit{OO}$ is assumed to be
equally optimal for $c_{H\tilde{W}}$ only.

All generated events were passed through a full simulation of the ATLAS
detector response [62] using Geant4 [63], except the Sherpa $\gamma\gamma$
sample, which was passed through a fast parametric simulation of the detector
response [62]. The effects of multiple $pp$ interactions in the same or
neighboring bunch crossings (pileup) are included by overlaying events
generated with Pythia 8 [64]. Events are weighted such that the distribution
of the average number of interactions per bunch crossing matches that observed
in data.

Photons are reconstructed from variable-size topological clusters formed from
electromagnetic calorimeter cells with significant energy deposits and from
tracks, initiated by converted photons, measured in the inner detector (ID)
[65]. Events must have at least two photon candidates outside the
calorimeter’s barrel/endcap transition region, $1.37<|\eta|<1.52$, and within
$|\eta|<2.37$, where the two leading (highest-$E_{\text{T}}$) photons are used
to reconstruct the Higgs boson candidate and the primary vertex of the event
[66]. The diphoton invariant mass $m_{\gamma\gamma}$ is required to be in the
range 105–160 GeV. The leading and subleading photons are further required to
have $E_{\mathrm{T}}/m_{\gamma\gamma}$ greater than 0.35 and 0.25,
respectively, and fulfill the Tight identification selection and Tight
calorimetric and track-based isolation requirement [65]. Jets are
reconstructed using the anti-$k_{t}$ algorithm [67, 68] with a radius
parameter $R=0.4$ from inputs formed with a particle-flow algorithm [69],
which uses information from both the calorimeter and the ID. Jet candidates
are required to have $p_{\mathrm{T}}>30\,\mathrm{GeV}$ and $|\eta|<4.4$. To
suppress jets from pileup collisions, jet candidates with $|\eta|<2.4$ and
$p_{\mathrm{T}}<60\,\mathrm{GeV}$ are required to pass the Tight jet vertex
tagger (JVT) selection [70]. For jets with $|\eta|\geq 2.4$, the Loose forward
JVT selection [71] is applied to remove pileup jet contamination. To construct
the region enriched with VBF signal events, two loose criteria are applied:
events must have at least two jets with pseudorapidity separation
$|\Delta\eta_{jj}|>2$ and Zeppenfeld variable [72]
$\eta^{\mathrm{Zepp}}=|\eta_{\gamma\gamma}-(\eta_{j1}+\eta_{j2})/2|<5$.

To increase the VBF signal purity, two boosted decision trees (BDT) [73] are
trained. $\mathrm{BDT_{VBF/ggF}}$ is used to separate VBF signal from ggF
events, which are the major background from Higgs boson production.
$\mathrm{BDT_{VBF/Continuum}}$ is used to distinguish VBF
$H\rightarrow\gamma\gamma$ events from continuum background events, which
consist of the prompt diphoton events ($\gamma\gamma$) and events where one or
two of the photon candidates originate from jets misidentified as photons
($\gamma j$ or $jj$). The $\gamma\gamma$ events, which are the dominant
component of the continuum background, are obtained from simulation, while
$\gamma j$ and $jj$ events are obtained from dedicated data control regions,
as described later. The two BDTs use the same input variables: invariant mass
of the dijet system formed by the two leading jets ($m_{jj}$), pseudorapidity
separation of the dijet system ($\Delta\eta_{jj}$), perpendicular projection
of the diphoton $p_{\text{T}}$ onto the diphoton thrust axis
($p^{\gamma\gamma}_{\mathrm{Tt}}$) [74], $p_{\text{T}}$ of the Higgs boson and
the leading two jets ($p_{\mathrm{T}}^{Hjj}$), azimuthal angle between the
diphoton and dijet systems ($\Delta\phi(\gamma\gamma,jj)$), minimum angular
separation between the photons and the two leading jets ($\Delta
R^{\mathrm{min}}_{\gamma j}$), and $\eta^{\mathrm{Zepp}}$. These input
variables are designed to be insensitive to the CP property of the VBF signal
and to have negligible correlation with $m_{\gamma\gamma}$. Figure 1 shows the
BDT output distributions of the VBF signal, the ggF background, the continuum
background, and the data in the $m_{\gamma\gamma}$ sideband
($m_{\gamma\gamma}\in[105,118]\leavevmode\nobreak\ \text{GeV}$ or
$[132,160]\leavevmode\nobreak\ \text{GeV}$). The comparison between the
continuum background and the sideband data shows the continuum background used
in the BDT training is well modeled. Events are categorized as follows:
firstly, a requirement is placed on $\mathrm{BDT_{VBF/ggF}}$ to separate
events into ‘tight’ (T) and ‘loose’ (L) regions. The ratio of VBF signal to
ggF background is improved by a factor of ten in the ‘tight’ region. Then, two
independent requirements on $\mathrm{BDT_{VBF/Continuum}}$ are applied to the
‘tight’ and ‘loose’ regions to maximize the combined significance of the VBF
signal. Three signal regions are defined: TT, TL, and LT, where the first
(second) letter corresponds to the $\mathrm{BDT_{VBF/ggF}}$
($\mathrm{BDT_{VBF/Continuum}}$) separation type. More details on the BDT
requirements can be found in the Appendix. In the TT and TL categories, the
dominant Higgs boson backgrounds are from the ggF process, and the
contributions from non-ggF Higgs processes, e.g. $VH$, $ttH$ and $tH$, are
found to be negligible. In the LT category, Higgs boson backgrounds are still
mostly from the ggF process, while those from non-ggF Higgs processes increase
to about 1%–3% of the VBF event yield.

Figure 1: Distribution of the output of $\mathrm{BDT_{VBF/ggF}}$ (left) and
$\mathrm{BDT_{VBF/continuum}}$ (right). The comparison between the continuum
background and the sideband data indicates the continuum background used in
the BDT training is well modeled.

The $\mathit{OO}$ value is computed from the reconstructed four-momenta of the
Higgs boson and the two VBF jets, from which the momentum fractions of the two
initial-state partons are derived in the beam directions by exploiting energy-
momentum conservation, as detailed in Ref. [9]. The signal yield is extracted
via a combined unbinned maximum-likelihood estimation applied to the
$m_{\gamma\gamma}$ distribution of observed data in each $\mathit{OO}$ bin, as
shown in Figure 2. Both the signal and background shapes are modeled with
analytic functions. The $H\to\gamma\gamma$ signal shape is described by a
double-sided Crystal Ball (DSCB) function [12], consisting of a Gaussian
distribution in the region around the peak, continued by power-law tails at
lower and higher $m_{\gamma\gamma}$ values. The parameters of the DSCB
function in each category are obtained by a fit to the simulated VBF sample,
as well as other Higgs boson production modes in proportion to their SM cross-
sections.

The modeling of the continuum background relies on both simulation and data-
driven methods. The $m_{\gamma\gamma}$ shape of the $\gamma\gamma$ component
is estimated using the Sherpa sample, while the $m_{\gamma\gamma}$ shapes of
the $\gamma j$ and $jj$ components are obtained using data control regions
formed by inverting the Tight photon identification and isolation
requirements. The template is then built by summing the $\gamma\gamma$,
$\gamma j$ and $jj$ components, where their fractions are measured in data
using a two-dimensional double-sideband method [75]. The composition of the
continuum background is found to be approximately 85% $\gamma\gamma$ events
and 15% $\gamma j+jj$ events. The background templates are smoothed using
Gaussian process regression (GPR) [76] with the Gibbs kernel to reduce
fluctuations due to the limited sample size. The $m_{\gamma\gamma}$
distribution of the continuum background is found to have a smoothly falling
shape. The analytic function chosen to model the continuum background is
either a power-law function, a Bernstein polynomial [12], or an exponential
function of a polynomial, and it is selected for each $\mathit{OO}$ bin
independently. The selected function should have the smallest spurious signal,
defined as the systematic bias in the fitted signal yield due to differences
between the fit function choice and the background template. The coefficients
of these functions are considered to be independent across categories, and in
all cases are treated as free parameters in the fits to data. More details can
be found in Ref. [12].

Figure 2:  Distribution of the optimal observable $\mathit{OO}$ for events
with $m_{\gamma\gamma}\in[118,132]\leavevmode\nobreak\ \text{GeV}$.
Contributions in three signal regions are summed together with a weight of
$\ln(1+S/B)$ for each signal region, where $S$ and $B$ are the expected yields
of signal and background events with
$m_{\gamma\gamma}\in[118,132]\leavevmode\nobreak\ \text{GeV}$. The overflow
and underflow are included in the highest and lowest bin, respectively. The
uncertainty band shown includes all systematic uncertainties. The weighted
summed $m_{\gamma\gamma}$ distribution of data events is shown in the inner
panel along with the signal and background contributions. The lower panel is
the $\mathit{OO}$ distribution in data after subtraction of all backgrounds,
in comparison with the SM VBF process, and VBF processes with $\tilde{d}=0.06$
and $\tilde{d}=-0.06$. The sensitivity to $\tilde{d}$ is dominated by the
tails of the $\mathit{OO}$ distribution.

An unbinned likelihood is constructed with the $m_{\gamma\gamma}$ spectra of
each $\mathit{OO}$ bin in signal regions TT, TL and LT. The negative log-
likelihood ($\mathrm{NLL}$) is evaluated for various $\tilde{d}$ and
$c_{H\tilde{W}}$ hypotheses. Confidence intervals are obtained by reading
values off the $\mathrm{NLL}$ curve, which is constructed by interpolating
between the points with spline functions. The normalization of the signal is
allowed to float in the fit. The analysis therefore exploits only the shape of
the distribution of the Optimal Observable, and ignores the potential
dependence of the inclusive cross-section on CP-mixing scenarios. If present,
any BSM CP-even effects would mainly change the normalization, and produce
very small symmetric changes in the $\mathit{OO}$ distribution, which are
found to not bias the parameter of interest for the CP-odd effect. All other
Higgs boson production modes are considered as backgrounds and are normalized
to their SM predicted yields. The expected $\Delta\mathrm{NLL}$ curve is
obtained using a pseudo-dataset where the event yields and distributions in
the signal regions are set to the SM expectations for both the signal and
background processes.

Both the theoretical and experimental systematic uncertainties are
incorporated into the likelihood model of the measurement as nuisance
parameters. Theoretical uncertainties arise from the modeling of VBF and ggF
processes because of the missing higher-order terms in the perturbative QCD
calculations, the modeling of the underlying event and parton shower, the
parton distribution functions, and the value of $\alpha_{\text{s}}$. These
uncertainties are estimated by following the procedure described in Ref. [12].
The experimental uncertainties include the uncertainties in the photon energy
scale and resolution [65], the jet energy scale and resolution [77], the
luminosity measurement, and the modeling of pileup events and the photon
identification and isolation criteria [78]. The spurious signal that could
arise from mismodeling of continuum background is estimated in each
$\mathit{OO}$ bin.

Figure 3 shows the $\Delta\mathrm{NLL}$ curves as functions of $\tilde{d}$ or
$c_{H\tilde{W}}$. Here, the $\tilde{d}$ results use the interference-plus-
quadratic terms in Eq. (1), while the $c_{H\tilde{W}}$ results use only the
interference term. The confidence intervals for the two scenarios,
interference-only and interference-plus-quadratic, are shown in Table 1. The
difference between the results in the two scenarios is found to be small. The
results are compatible with the SM and the precision is limited by the
statistical uncertainty of the data. For example, the total impact on the 95%
confidence intervals of $\tilde{d}$ from the systematic uncertainty is less
than 2%. The measurement is sensitive enough to determine an observed 95%
confidence interval for $\tilde{d}$, which was not achieved in previous
analyses. The expected 68% confidence interval shown for the
$H\rightarrow\tau\tau$ channel in Table 1 differs slightly from that presented
in Ref. [9], where the expected $H\rightarrow\tau\tau$ results were obtained
with the nuisance parameters constrained only by the control regions. In the
present analysis, the expected results are obtained with the nuisance
parameters constrained by both the control regions and signal regions.

The 95% confidence interval for $c_{H\tilde{W}}$ obtained using the
interference-only term is a factor of five more restrictive than in the
$H\rightarrow\gamma\gamma$ differential measurement reported in Ref. [12]
because of the dedicated BDTs for the VBF signal selection and the use of the
Optimal Observable. The 68% confidence interval for $c_{H\tilde{W}}$ is about
twice as restrictive as that from either the ATLAS or CMS $H\rightarrow ZZ$
four-lepton analysis [13, 16]. Compared to the previous results for
$H\rightarrow\tau\tau$ [9], based on the interference-plus-quadratic scenario,
the expected 68% confidence interval for $\tilde{d}$ has been tightened by
about 20%. After combining the present results with the $H\rightarrow\tau\tau$
results, the expected confidence intervals are further improved by about 15%.
The luminosity uncertainty of the data in 2015–2016, the uncertainties of the
electron/photon energy scale and resolution, and the theoretical uncertainty
of VBF and ggF processes are correlated in the combination. The jet-related
uncertainties are not correlated since a different jet reconstruction
technique was used in the $H\rightarrow\tau\tau$ analysis.

Figure 3:  $\Delta\mathrm{NLL}$ curves as a function of (a) $\tilde{d}$ and (b) $c_{H\tilde{W}}$. In figure (a), the $\Delta\mathrm{NLL}$ of $\tilde{d}$ considers the interference-plus-quadratic terms, whereas in figure (b) the $\Delta\mathrm{NLL}$ of $c_{H\tilde{W}}$ considers the interference-only term. The solid lines are the observed results, while the dashed lines are the expected results. In figure (a), the blue lines represent the results of this analysis, while the red lines represent the results from the $H\rightarrow\tau\tau$ analysis [9]. The black lines show the combination of these two analyses. For all figures, the dashed horizontal lines show the values of $\Delta\mathrm{NLL}$ used to define the 68% and 95% confidence intervals.  Table 1:  Observed and expected 68% and 95% confidence intervals for $\tilde{d}$ and $c_{H\tilde{W}}$. Results for scenarios with the interference-only (noted as ‘inter. only’) term and interference-plus-quadratic terms (noted as ‘inter.+quad.’) are both presented. Combined results for $\tilde{d}$ including the $H\rightarrow\tau\tau$ analysis [9] are shown. The expected results of $H\rightarrow\tau\tau$ are slightly different from Ref. [9] due to the different correlation scheme between their signal region and control region.  | 68% (exp.) | 95% (exp.) | 68% (obs.) | 95% (obs.)
---|---|---|---|---
$\tilde{d}$ (inter. only) | [$-$0.027, 0.027] | [$-$0.055, 0.055] | [$-$0.011, 0.036] | [$-$0.032, 0.059]
$\tilde{d}$ (inter.+quad.) | [$-$0.028, 0.028] | [$-$0.061, 0.060] | [$-$0.010, 0.040] | [$-$0.034, 0.071]
$\tilde{d}$ from $H\rightarrow\tau\tau$ | [$-$0.038, 0.036] | $-$ | [$-$0.090, 0.035] | -
Combined $\tilde{d}$ | [$-$0.022, 0.021] | [$-$0.046, 0.045] | [$-$0.012, 0.030] | [$-$0.034, 0.057]
$c_{H\tilde{W}}$ (inter. only) | [$-$0.48, 0.48] | [$-$0.94, 0.94] | [$-$0.16, 0.64] | [$-$0.53, 1.02]
$c_{H\tilde{W}}$ (inter.+quad.) | [$-$0.48, 0.48] | [$-$0.95, 0.95] | [$-$0.15, 0.67] | [$-$0.55, 1.07]

In conclusion, a test of CP invariance in Higgs boson production via vector-
boson fusion is performed in the $H\rightarrow\gamma\gamma$ channel using 139
fb-1 of $\sqrt{s}=13\,\text{TeV}$ proton–proton collision data collected by
the ATLAS detector at the LHC. The Optimal Observable method is used to probe
CP-violating interactions between the Higgs boson and electroweak gauge bosons
described by an effective field theory. The results are compatible with the
SM. No sign of CP violation is observed in the Optimal Observable
distributions. The constraints on CP-violating effects in the $HVV$ coupling
are the most stringent to date. They allow 68% and 95% confidence intervals to
be set for parameters describing the strength of the CP-odd component in the
$HVV$ coupling in two effective field theory bases: $\tilde{d}$ in the HISZ
basis and $c_{H\tilde{W}}$ in the Warsaw basis. The sensitivity is sufficient
to set a 95% confidence interval for $\tilde{d}$ for the first time, and the
constraints on $\tilde{d}$ are tightened further by combining them with
previous results from the $H\rightarrow\tau\tau$ channel. The constraints on
$c_{H\tilde{W}}$ are about twice as restrictive as those from either the ATLAS
or CMS four-lepton analysis.

## Appendix

Appendix: BDT based event categorization. —This analysis uses two BDTs to
enrich the VBF signal purity: $\mathrm{BDT_{VBF/ggF}}$ and
$\mathrm{BDT_{VBF/Continuum}}$, which are trained to separate the VBF signal
from ggF background and continuum background, respectively. Figure 4 shows the
distributions of the BDT output, as well as the requirements for the event
categorization. The left panel shows the two-dimensional distributions of
$\mathrm{BDT_{VBF/continuum}}$ v.s. $\mathrm{BDT_{VBF/ggF}}$ for both data in
the $m_{\gamma\gamma}$ sideband and VBF signal which are normalized to unity.
The vertical line at 0.15 indicates the requirement on
$\mathrm{BDT_{VBF/ggF}}$ to separate the tight and loose regions. The ratio of
VBF signal over ggF background is improved by a factor of ten in the tight
region. The horizontal line at 0.32 (0.09) indicates the requirement on
$\mathrm{BDT_{VBF/Continuum}}$ in the $\mathrm{BDT_{VBF/ggF}}$ tight (loose)
region, to maximize the combined significance for the VBF signal. The events
are finally categorized into four regions: TT, TL, LT, and LL. The right panel
shows the $\mathrm{BDT_{VBF/continuum}}$ distributions in the
$\mathrm{BDT_{VBF/ggF}}$ tight and loose regions, while the inclusive
$\mathrm{BDT_{VBF/ggF}}$ distribution is shown in Figure 1.

Figure 4:  Left: Normalized two-dimensional distributions of
$\mathrm{BDT_{VBF/continuum}}$ v.s. $\mathrm{BDT_{VBF/ggF}}$ for both data in
the $m_{\gamma\gamma}$ sideband and VBF signal. The vertical line indicates
the requirement on $\mathrm{BDT_{VBF/ggF}}$ to separate the tight and loose
regions, while the horizontal lines indicate the requirements on
$\mathrm{BDT_{VBF/continuum}}$. The events are categorized into four regions:
TT, TL, LT, and LL. Right: $\mathrm{BDT_{VBF/continuum}}$ distributions in the
$\mathrm{BDT_{VBF/ggF}}$ tight and loose regions. The vertical lines indicate
the requirements on $\mathrm{BDT_{VBF/continuum}}$.

### References

  * [1] A.. Sakharov  “Violation of CP Invariance, C asymmetry, and baryon asymmetry of the universe”  In _Pisma Zh. Eksp. Teor. Fiz._ 5, 1967, pp. 32–35  DOI: 10.1070/PU1991v034n05ABEH002497
  * [2] Makoto Kobayashi and Toshihide Maskawa  “CP Violation in the Renormalizable Theory of Weak Interaction”  In _Prog. Theor. Phys._ 49, 1973, pp. 652–657  DOI: 10.1143/PTP.49.652
  * [3] M.. Gavela, P. Hernandez, J. Orloff and O. Pene  “Standard model CP violation and baryon asymmetry”  In _Mod. Phys. Lett. A_ 9, 1994, pp. 795–810  DOI: 10.1142/S0217732394000629
  * [4] Andrew G. Cohen, D.. Kaplan and A.. Nelson  “Progress in Electroweak Baryogenesis”  In _Ann. Rev. Nucl. Part. Sci._ 43, 1993, pp. 27–70  DOI: 10.1146/annurev.ns.43.120193.000331
  * [5] Patrick Huet and Eric Sather  “Electroweak baryogenesis and standard model CP violation”  In _Phys. Rev. D_ 51, 1995, pp. 379–394  DOI: 10.1103/PhysRevD.51.379
  * [6] ATLAS Collaboration  “Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC”  In _Phys. Lett. B_ 716, 2012, pp. 1  DOI: 10.1016/j.physletb.2012.08.020
  * [7] CMS Collaboration  “Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC”  In _Phys. Lett. B_ 716, 2012, pp. 30  DOI: 10.1016/j.physletb.2012.08.021
  * [8] Lyndon Evans and Philip Bryant  “LHC Machine”  In _JINST_ 3, 2008, pp. S08001  DOI: 10.1088/1748-0221/3/08/S08001
  * [9] ATLAS Collaboration  “Test of CP invariance in vector-boson fusion production of the Higgs boson in the $H\to\tau\tau$ channel in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector”  In _Phys. Lett. B_ 805, 2020, pp. 135426  DOI: 10.1016/j.physletb.2020.135426
  * [10] ATLAS Collaboration  “$CP$ Properties of Higgs Boson Interactions with Top Quarks in the $t\bar{t}H$and $tH$ Processes Using $H\to\gamma\gamma$ with the ATLAS Detector”  In _Phys. Rev. Lett._ 125, 2020, pp. 061802  DOI: 10.1103/PhysRevLett.125.061802
  * [11] ATLAS Collaboration  “Constraints on Higgs boson properties using $WW^{*}(\rightarrow e\nu\mu\nu)jj$ production in $36.1\,\text{fb}^{-1}$ of $\sqrt{s}=13\,\text{TeV}$ $pp$ collisions with the ATLAS detector”, 2021  arXiv:2109.13808 [hep-ex]
  * [12] ATLAS Collaboration  “Measurements of the Higgs boson inclusive and differential fiducial cross-sections in the diphoton decay channel with $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector”, 2022  arXiv:2202.00487 [hep-ex]
  * [13] ATLAS Collaboration  “Higgs boson production cross-section measurements and their EFT interpretation in the $4\ell$ decay channel at $\sqrt{s}=$13 TeV with the ATLAS detector” [Erratum: Eur. Phys. J. C 81 (2021) 29, Erratum: Eur. Phys. J. C 81 (2021) 398]  In _Eur. Phys. J. C_ 80.10, 2020, pp. 957  DOI: 10.1140/epjc/s10052-020-8227-9
  * [14] CMS Collaboration  “Constraints on anomalous $HVV$ couplings from the production of Higgs bosons decaying to $\tau$ lepton pairs”  In _Phys. Rev. D_ 100, 2019, pp. 112002  DOI: 10.1103/PhysRevD.100.112002
  * [15] CMS Collaboration  “Measurements of $t\bar{t}H$ Production and the $CP$ Structure of the Yukawa Interaction between the Higgs Boson and Top Quark in the Diphoton Decay Channel”  In _Phys. Rev. Lett._ 125, 2020, pp. 061801  DOI: 10.1103/PhysRevLett.125.061801
  * [16] CMS Collaboration  “Constraints on anomalous Higgs boson couplings to vector bosons and fermions in its production and decay using the four-lepton final state”  In _Phys. Rev. D_ 104, 2021, pp. 052004  DOI: 10.1103/PhysRevD.104.052004
  * [17] CMS Collaboration  “Analysis of the CP structure of the Yukawa coupling between the Higgs boson and $\tau$ leptons in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$”, 2021  arXiv:2110.04836 [hep-ex]
  * [18] CMS Collaboration  “Constraints on anomalous Higgs boson couplings to vector bosons and fermions from the production of Higgs bosons using the $\tau\tau$ final state”, 2022  arXiv:2205.05120 [hep-ex]
  * [19] Aleksandr Azatov, Roberto Contino, Camila S. Machado and Francesco Riva  “Helicity selection rules and noninterference for BSM amplitudes”  In _Phys. Rev. D_ 95.6, 2017, pp. 065014  DOI: 10.1103/PhysRevD.95.065014
  * [20] ALEPH Collaboration  “Improved measurement of the triple gauge-boson couplings $\gamma$WW and ZWW in $e^{+}e^{-}$ collisions”  In _Phys. Lett. B_ 614, 2005, pp. 7–26  DOI: 10.1016/j.physletb.2005.03.058
  * [21] OPAL Collaboration  “Measurement of $W$ boson polarizations and CP violating triple gauge couplings from $W^{+}W^{-}$ production at LEP”  In _Eur. Phys. J. C_ 19, 2001, pp. 229–240  DOI: 10.1007/s100520100602
  * [22] M. Davier, L. Duflot, F. Le Diberder and A. Rouge  “The Optimal method for the measurement of tau polarization”  In _Phys. Lett. B_ 306, 1993, pp. 411–417  DOI: 10.1016/0370-2693(93)90101-M
  * [23] D. Atwood and A. Soni  “Analysis for magnetic moment and electric dipole moment form factors of the top quark via ${\mathit{e}}^{+}$${\mathit{e}}^{\mathrm{-}}$$\rightarrow$tt¯”  In _Phys. Rev. D_ 45 American Physical Society, 1992, pp. 2405–2413  DOI: 10.1103/PhysRevD.45.2405
  * [24] M. Diehl and O. Nachtmann  “Optimal observables for the measurement of three gauge boson couplings in $e^{+}e^{-}\rightarrow W^{+}W^{-}$”  In _Z. Phys. C_ 62, 1994, pp. 397–412  DOI: 10.1007/BF01555899
  * [25] Johann Brehmer, Felix Kling, Tilman Plehn and Tim M.. Tait  “Better Higgs-CP tests through information geometry”  In _Phys. Rev. D_ 97.9, 2018, pp. 095017  DOI: 10.1103/PhysRevD.97.095017
  * [26] Kaoru Hagiwara, S. Ishihara, R. Szalapski and D. Zeppenfeld  “Low energy effects of new interactions in the electroweak boson sector”  In _Phys. Rev. D_ 48, 1993, pp. 2182–2203  DOI: 10.1103/PhysRevD.48.2182
  * [27] Ilaria Brivio, Yun Jiang and Michael Trott  “The SMEFTsim package, theory and tools”  In _JHEP_ 12, 2017, pp. 070  DOI: 10.1007/JHEP12(2017)070
  * [28] Ilaria Brivio  “SMEFTsim 3.0 — a practical guide”  In _JHEP_ 04, 2021, pp. 073  DOI: 10.1007/JHEP04(2021)073
  * [29] B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek  “Dimension-six terms in the Standard Model Lagrangian”  In _JHEP_ 10, 2010, pp. 085  DOI: 10.1007/JHEP10(2010)085
  * [30] ATLAS Collaboration  “The ATLAS Experiment at the CERN Large Hadron Collider”  In _JINST_ 3, 2008, pp. S08003  DOI: 10.1088/1748-0221/3/08/S08003
  * [31] B. Abbott  “Production and integration of the ATLAS Insertable B-Layer”  In _JINST_ 13, 2018, pp. T05008  DOI: 10.1088/1748-0221/13/05/T05008
  * [32] ATLAS Collaboration  “ATLAS Insertable B-Layer: Technical Design Report”, 2010  URL: https://cds.cern.ch/record/1291633
  * [33] ATLAS Collaboration  “Performance of the ATLAS trigger system in 2015”  In _Eur. Phys. J. C_ 77, 2017, pp. 317  DOI: 10.1140/epjc/s10052-017-4852-3
  * [34] ATLAS Collaboration  “Performance of electron and photon triggers in ATLAS during LHC Run 2”  In _Eur. Phys. J. C_ 80, 2020, pp. 47  DOI: 10.1140/epjc/s10052-019-7500-2
  * [35] ATLAS Collaboration  “The ATLAS Collaboration Software and Firmware”, ATL-SOFT-PUB-2021-001, 2021  URL: https://cds.cern.ch/record/2767187
  * [36] Paolo Nason and Carlo Oleari  “NLO Higgs boson production via vector-boson fusion matched with shower in POWHEG”  In _JHEP_ 02, 2010, pp. 037  DOI: 10.1007/JHEP02(2010)037
  * [37] Jon Butterworth  “PDF4LHC recommendations for LHC Run II”  In _J. Phys. G_ 43, 2016, pp. 023001  DOI: 10.1088/0954-3899/43/2/023001
  * [38] M. Ciccolini, Ansgar Denner and S. Dittmaier  “Strong and Electroweak Corrections to the Production of a Higgs Boson + 2 Jets via Weak Interactions at the Large Hadron Collider”  In _Phys. Rev. Lett._ 99, 2007, pp. 161803  DOI: 10.1103/PhysRevLett.99.161803
  * [39] Mariano Ciccolini, Ansgar Denner and Stefan Dittmaier  “Electroweak and QCD corrections to Higgs production via vector-boson fusion at the CERN LHC”  In _Phys. Rev. D_ 77, 2008, pp. 013002  DOI: 10.1103/PhysRevD.77.013002
  * [40] Paolo Bolzoni, Fabio Maltoni, Sven-Olaf Moch and Marco Zaro  “Higgs Boson Production via Vector-Boson Fusion at Next-to-Next-to-Leading Order in QCD”  In _Phys. Rev. Lett._ 105, 2010, pp. 011801  DOI: 10.1103/PhysRevLett.105.011801
  * [41] Keith Hamilton, Paolo Nason, Emanuele Re and Giulia Zanderighi  “NNLOPS simulation of Higgs boson production”  In _JHEP_ 10, 2013, pp. 222  DOI: 10.1007/JHEP10(2013)222
  * [42] Keith Hamilton, Paolo Nason and Giulia Zanderighi  “Finite quark-mass effects in the NNLOPS POWHEG+MiNLO Higgs generator”  In _JHEP_ 05, 2015, pp. 140  DOI: 10.1007/JHEP05(2015)140
  * [43] Keith Hamilton, Paolo Nason and Giulia Zanderighi  “MINLO: multi-scale improved NLO”  In _JHEP_ 10, 2012, pp. 155  DOI: 10.1007/JHEP10(2012)155
  * [44] John M. Campbell et al.  “NLO Higgs boson production plus one and two jets using the POWHEG BOX, MadGraph4 and MCFM”  In _JHEP_ 07, 2012, pp. 092  DOI: 10.1007/JHEP07(2012)092
  * [45] Keith Hamilton, Paolo Nason, Carlo Oleari and Giulia Zanderighi  “Merging H/W/Z + 0 and 1 jet at NLO with no merging scale: a path to parton shower + NNLO matching”  In _JHEP_ 05, 2013, pp. 082  DOI: 10.1007/JHEP05(2013)082
  * [46] Stefano Catani and Massimiliano Grazzini  “Next-to-Next-to-Leading-Order Subtraction Formalism in Hadron Collisions and its Application to Higgs-boson Production at the Large Hadron Collider”  In _Phys. Rev. Lett._ 98, 2007, pp. 222002  DOI: 10.1103/PhysRevLett.98.222002
  * [47] D. Florian  “Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector”, 2016  DOI: 10.23731/CYRM-2017-002
  * [48] Charalampos Anastasiou et al.  “High precision determination of the gluon fusion Higgs boson cross-section at the LHC”  In _JHEP_ 05, 2016, pp. 058  DOI: 10.1007/JHEP05(2016)058
  * [49] Charalampos Anastasiou et al.  “Higgs Boson Gluon-Fusion Production in QCD at Three Loops”  In _Phys. Rev. Lett._ 114, 2015, pp. 212001  DOI: 10.1103/PhysRevLett.114.212001
  * [50] Falko Dulat, Achilleas Lazopoulos and Bernhard Mistlberger  “iHixs 2 – Inclusive Higgs cross sections”  In _Comput. Phys. Commun._ 233, 2018, pp. 243–260  DOI: 10.1016/j.cpc.2018.06.025
  * [51] Robert V. Harlander and Kemal J. Ozeren  “Finite top mass effects for hadronic Higgs production at next-to-next-to-leading order”  In _JHEP_ 11, 2009, pp. 088  DOI: 10.1088/1126-6708/2009/11/088
  * [52] Robert V. Harlander and Kemal J. Ozeren  “Top mass effects in Higgs production at next-to-next-to-leading order QCD: Virtual corrections”  In _Phys. Lett. B_ 679, 2009, pp. 467–472  DOI: 10.1016/j.physletb.2009.08.012
  * [53] Robert V. Harlander, Hendrik Mantler, Simone Marzani and Kemal J. Ozeren  “Higgs production in gluon fusion at next-to-next-to-leading order QCD for finite top mass”  In _Eur. Phys. J. C_ 66, 2010, pp. 359–372  DOI: 10.1140/epjc/s10052-010-1258-x
  * [54] Alexey Pak, Mikhail Rogal and Matthias Steinhauser  “Finite top quark mass effects in NNLO Higgs boson production at LHC”  In _JHEP_ 02, 2010, pp. 025  DOI: 10.1007/JHEP02(2010)025
  * [55] Stefano Actis, Giampiero Passarino, Christian Sturm and Sandro Uccirati  “NLO electroweak corrections to Higgs boson production at hadron colliders”  In _Phys. Lett. B_ 670, 2008, pp. 12–17  DOI: 10.1016/j.physletb.2008.10.018
  * [56] Stefano Actis, Giampiero Passarino, Christian Sturm and Sandro Uccirati  “NNLO computational techniques: The cases $H\to\gamma\gamma$ and $H\to gg$”  In _Nucl. Phys. B_ 811, 2009, pp. 182–273  DOI: 10.1016/j.nuclphysb.2008.11.024
  * [57] Marco Bonetti, Kirill Melnikov and Lorenzo Tancredi  “Higher order corrections to mixed QCD-EW contributions to Higgs boson production in gluon fusion”  In _Phys. Rev. D_ 97.5, 2018, pp. 056017  DOI: 10.1103/PhysRevD.97.056017
  * [58] Enrico Bothmann  “Event Generation with Sherpa 2.2”  In _SciPost Phys._ 7.3, 2019, pp. 034  DOI: 10.21468/SciPostPhys.7.3.034
  * [59] Ansgar Denner, Stefan Dittmaier, Stefan Kallweit and Alexander Mück  “HAWK 2.0: A Monte Carlo program for Higgs production in vector-boson fusion and Higgs strahlung at hadron colliders”  In _Comput. Phys. Commun._ 195, 2015, pp. 161–171  DOI: 10.1016/j.cpc.2015.04.021
  * [60] J. Alwall et al.  “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”  In _JHEP_ 07, 2014, pp. 079  DOI: 10.1007/JHEP07(2014)079
  * [61] Rikkert Frederix and Stefano Frixione  “Merging meets matching in MC@NLO”  In _JHEP_ 12, 2012, pp. 061  DOI: 10.1007/JHEP12(2012)061
  * [62] ATLAS Collaboration  “The ATLAS Simulation Infrastructure”  In _Eur. Phys. J. C_ 70, 2010, pp. 823  DOI: 10.1140/epjc/s10052-010-1429-9
  * [63] GEANT4 Collaboration  “Geant4–a simulation toolkit”  In _Nucl. Instrum. Meth. A_ 506, 2003, pp. 250–303  DOI: 10.1016/S0168-9002(03)01368-8
  * [64] Torbjörn Sjöstrand et al.  “An introduction to PYTHIA 8.2”  In _Comput. Phys. Commun._ 191, 2015, pp. 159  DOI: 10.1016/j.cpc.2015.01.024
  * [65] ATLAS Collaboration  “Electron and photon performance measurements with the ATLAS detector using the 2015–2017 LHC proton–proton collision data”  In _JINST_ 14, 2019, pp. P12006  DOI: 10.1088/1748-0221/14/12/P12006
  * [66] ATLAS Collaboration  “Measurement of Higgs boson production in the diphoton decay channel in $pp$ collisions at center-of-mass energies of $7$ and $8\,\text{TeV}$ with the ATLAS detector”  In _Phys. Rev. D_ 90, 2014, pp. 112015  DOI: 10.1103/PhysRevD.90.112015
  * [67] Matteo Cacciari, Gavin P. Salam and Gregory Soyez  “The anti-$k_{t}$ jet clustering algorithm”  In _JHEP_ 04, 2008, pp. 063  DOI: 10.1088/1126-6708/2008/04/063
  * [68] Matteo Cacciari, Gavin P. Salam and Gregory Soyez  “FastJet user manual”  In _Eur. Phys. J. C_ 72, 2012, pp. 1896  DOI: 10.1140/epjc/s10052-012-1896-2
  * [69] ATLAS Collaboration  “Optimisation of large-radius jet reconstruction for the ATLAS detector in $13\,\text{TeV}$ proton–proton collisions”  In _Eur. Phys. J. C_ 81, 2020, pp. 334  DOI: 10.1140/epjc/s10052-021-09054-3
  * [70] ATLAS Collaboration  “Performance of pile-up mitigation techniques for jets in $pp$ collisions at $\sqrt{s}=8\,\text{TeV}$ using the ATLAS detector”  In _Eur. Phys. J. C_ 76, 2016, pp. 581  DOI: 10.1140/epjc/s10052-016-4395-z
  * [71] ATLAS Collaboration  “Identification and rejection of pile-up jets at high pseudorapidity with the ATLAS detector”  In _Eur. Phys. J. C_ 77, 2017, pp. 580  DOI: 10.1140/epjc/s10052-017-5081-5
  * [72] Jeppe R. Andersen, Ken Arnold and Dieter Zeppenfeld  “Azimuthal angle correlations for Higgs boson plus multi-jet events”  In _JHEP_ 06, 2010, pp. 091  DOI: 10.1007/JHEP06(2010)091
  * [73] Andreas Hoecker  “TMVA - Toolkit for Multivariate Data Analysis”, 2007  arXiv:physics/0703039 [physics.data-an]
  * [74] ATLAS Collaboration  “Measurements of Higgs boson production and couplings in diboson final states with the ATLAS detector at the LHC”  In _Phys. Lett. B_ 726, 2013, pp. 88  DOI: 10.1016/j.physletb.2013.08.010
  * [75] ATLAS Collaboration  “Measurement of the isolated diphoton cross section in $pp$ collisions at $\sqrt{s}=7\,\text{TeV}$ with the ATLAS detector”  In _Phys. Rev. D_ 85, 2012, pp. 012003  DOI: 10.1103/PhysRevD.85.012003
  * [76] Meghan Frate et al.  “Modeling Smooth Backgrounds and Generic Localized Signals with Gaussian Processes”, 2017  arXiv:1709.05681 [physics.data-an]
  * [77] ATLAS Collaboration  “Jet energy scale and resolution measured in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector”  In _Eur. Phys. J. C_ 81, 2020, pp. 689  DOI: 10.1140/epjc/s10052-021-09402-3
  * [78] ATLAS Collaboration  “Measurement of the photon identification efficiencies with the ATLAS detector using LHC Run 2 data collected in 2015 and 2016”  In _Eur. Phys. J. C_ 79, 2019, pp. 205  DOI: 10.1140/epjc/s10052-019-6650-6
  * [79] ATLAS Collaboration, 2012  URL: https://cds.cern.ch/record/1451888
  * [80] In _Phys. Rev. D_ 97 American Physical Society, 2018, pp. 099906  DOI: 10.1103/PhysRevD.97.099906
  * [81] ATLAS Collaboration  In _Eur. Phys. J. C_ 77, 2017, pp. 712  DOI: 10.1140/epjc/s10052-017-5245-3
  * [82] ATLAS Collaboration  In _Phys. Lett. B_ 734, 2014, pp. 406  DOI: 10.1016/j.physletb.2014.05.011