Transformed Gaussian Processes (s) are stochastic processes specified by transforming samples from the joint distribution from a prior process (typically a ) using an invertible transformation; increasing the flexibility of the base process.

Furthermore, they achieve competitive results compared with Deep Gaussian Processes (s), which are another generalization constructed by a hierarchical concatenation of s. In this work, we propose a generalization of s named Deep Transformed Gaussian Processes (s), which follows the trend of concatenating layers of stochastic processes. More precisely, we obtain a multi-layer model in which each layer is a . This generalization implies an increment of flexibility with respect to both s and s. Exact inference in such a model is intractable. However, we show that one can use variational inference to approximate the required computations yielding a straightforward extension of the popular inference algorithm [Salimbeni and Deisenroth, 2017]. The experiments conducted evaluate the proposed novel s in multiple regression datasets, achieving good scalability and performance.

§ INTRODUCTION

Although neural networks present highly accurate results on classification and regression tasks, they do not offer uncertainty estimations associated with the predictions made, which are mandatory in some fields. For example, in medicine, a typical problem is cancer detection using histopathological images, where given an image (or a set of images), the models try to determine if the patient has cancerous tissue or not [Schmidt et al., 2023]. In this case, the model must offer an interpretable (well-calibrated) output for the pathologists, who may use the output of the model for diagnosis. Bayesian learning offers a solution for this problem as it automatically outputs an estimation of the associated prediction uncertainty. However, it comes at a cost: computational tractability. The posterior distribution is most of the time intractable, so it has to be approximated using techniques such as variational inference.

Gaussian Processes (s) are a very powerful, non-parametric model that allows inference in the function space by placing a prior distribution on the target latent function [Rasmussen and Williams, 2005]. These models have been studied extensively in the literature [Titsias, 2009, Hensman et al., 2013], leading to different generalizations. The most popular are Deep Gaussian Processes (s) [Damianou and Lawrence, 2013], which use the output of a  as the input to another , increasing the expressiveness of the resulting model. Also, the usage of transformations on the prior and likelihood of the  has been explored. Using Normalizing Flows [Rezende and Mohamed, 2015], the Transformed Gaussian Process () [Maroñas et al., 2021] extends standard s by transforming  the prior distribution of the , which is no longer Gaussian, using an invertible transformation.

In this work, Deep Transformed Gaussian Processes (s) are introduced as a novel form of extending Transformed Gaussian Processes by concatenating the output of a to another . This model aims to improve the performance and the uncertainty estimation of previous models by being more flexible than them and thus being able to model more complex data. Inference in s is intractable, so the Evidence Lower Bound  is approximated using Monte Carlo samples, using a straightforward extension of the algorithm  proposed for inference in s by [Salimbeni and Deisenroth, 2017]. The usage of normalizing flows between each layer, leads to an increment in flexibility, adding a small term to the total computational cost.

To validate the proposed model, we conduct extensive experimentation with s. With this goal, we employ a specifically designed toy dataset and eight real datasets from the repository, which have been used to test the performance of s in comparison with other state-of-the-art models. Our results show that s obtain better or comparable results in almost all datasets.

§ PRELIMINARIES

The problem we aim to solve is to infer an unknown function \(f: \R^D \to \R\) given noisy observations \(\Ysamples = (y_1\cdots,y_N)^T\) at locations \(\Xsamples = (\X{1},\cdots, \X{N})\), where in our problem $\X{} \in\Xspace \subseteq \mathbb{R}^D$. Gaussian Processes place a prior on \(f\) such that the distribution of all function values is jointly Gaussian, with a mean function \(\mu: \mathcal X \to \R\) and covariance function \(K: \mathcal X \times \mathcal X \to \R\) [Rasmussen and Williams, 2005]. Since the computational cost of exact inference in s scales cubically with \(N\), a set of inducing locations \(\Zsamples{} = (\Z{1},\cdots,\Z{M})\) with \(M << N\) is considered [Titsias, 2009], aiming to reduce this cost. We denote \(\f{} = f(\Xsamples)\) and \(\uu{} = f(\Zsamples{})\). The joint distribution of \(\f{}, \uu{}\) and \(\Ysamples\) is

\begin{equation}
p(\Ysamples, \f{}, \uu{}) = \underbrace{\prod_{i=1}^N p(y_i \mid f_i)}_{\text{likelihood}} \underbrace{p(\f{} \mid \uu{}; \ \Xsamples, \Zsamples{}) p( \uu{} ; \ \Zsamples{})}_{\text{\GP prior}},
\end{equation}

where \(\uu{}\) is supposed to follow the same prior as \(\f{}\): \(p(\uu{}) = \mathcal N (\uu{} \mid m(\Zsamples{}), K(\Zsamples{}, \Zsamples{}))\).

Some authors have explored the possibility of transforming the prior or the likelihood distributions [Maroñas et al., 2021].
In this sense, using two mappings \(\T,\G\), the whole modelling process can be generalized as

\begin{equation}
\begin{rcases}        \fnull{} \sim \mathcal{GP} (\mu(\Xsamples), K(\Xsamples, \Xsamples));  &\ \fK{} = \G(\fnull{}) \\        \T(\Ysamples) = \fK{} + \epsilon; &\epsilon \sim \mathcal N(0, \Sigma)
\end{rcases}
\end{equation}

In this work, we are interested in transforming the prior using \(\G\) and maintaining \(\T\) as the identity. For this we consider  $\G$ to be given by an invertible transformation (Normalizing Flow), as in [Maroñas et al., 2021], followed by a hierarchical concatenation  as in [Damianou and Lawrence, 2013]. However, from now on $\G$  will just denote the invertible transformation. Since the composition of invertible transformation remains invertible and differentiable, we define $\G$ to be the composition of $K$ invertible transformations \(\Gcomp = \Gk{0} \circ \Gk{1} \circ \cdots \circ \Gk{K-1}\). This composition helps to increase the flexibility of the flow as much as we want. Also, in each of the steps of the flow the parameters \(\pmb{\theta}_k\) may depend on the input via a transformation such as a Neural Network (), giving rise to Input-Dependent () normalizing flows, which yield a non-stationary process with well-inductive biases [Maroñas and Hernández-Lobato, 2023, Maroñas et al., 2021]. In this case, parameters are given by functions $\theta:\Wspace\times\Xspace \to \mathbb{R}$ which we denote as \(\idparams{}{\Xsamples}\).

The Transformed Gaussian Process () is defined by the generative process given by composing a sample \(\f{0}\) of a with a normalizing flow:

\begin{equation}
\label{eq:tgp:definition}
\fnull{} \sim \mathcal{GP}(\mu( \Xsamples), K( \Xsamples, \Xsamples)), \quad \fK{} = \Gcomp(\fnull{}).
\end{equation}

Following [Maroñas et al., 2021] we consider element-wise mappings, which produce diagonal Jacobians. This, together with the application of the inverse function theorem and change of variable formula, leads to the distribution:

\begin{equation}\label{eq:tgp:transformed:distr}
p(\fK{} \mid \Gcomp , \Xsamples) = p(\fnull{} \mid \Xsamples) \Jacobian{\f{k}}.
\end{equation}

We denote \(\J{\pmb{a}}= \Jacobian{\pmb{a}}\). Unlike in standard s, inference in s is intractable. The posterior distribution is efficiently approximated using a variational distribution \(q(\fK{},\uK{}) = p(\fK{} \mid \uK{}) q(\uK{})\), which is chosen to contain a factor equal to the exact conditional prior $p(\fK{} \mid \uK{})$ and a marginal variational distribution \(q(\uK{}) = \mathcal N(\unull{} \mid \varm{}, \varS{})\J{\uK{}}\). A cancellation of both the conditional prior and the Jacobian \(\J{\uK{}}\) is achieved in the variational giving:

\begin{equation}\label{eq:tgp:elbo:simple}
\mathcal L_{\text{\TGP}} =  \Expectation{q(\fnull{})}{\log p(\Ysamples \mid \Gcomp (\fnull{})} -  \KL{q(\unull{})}{p(\unull{})}.
\end{equation}

This provides a bound that can be evaluated in \(\complexity{NM^2 + M^3 + NLK}\), where \(L\) is the computational cost of the normalizing flow. When flows are used, a Bayesian treatment can be considered by assigning a prior distribution \(p(\W{})\) to the weights of the , see [Maroñas et al., 2021].

§ DEEP TRANSFORMED GAUSSIAN PROCESSES

In this section, we present Deep Transformed Gaussian Processes (s) which generalize s through their hierarchical composition. Following the fashion of s, the output of the is used as the input of another , recursively defining a stochastic process:

A Deep Transformed Gaussian Process () is a collection of random variables \(\{\fpos{K}{h,l}\}_{h=1,l=1}^{H^l,L}\) with a hierarchical dependency such that
\(\fpos{K}{h,l}(\{\fpos{K}{h,l-1}\}_{h=1}^{H^{l-1}})\). Each of those functions follows the generative process:

\begin{equation}
\begin{rcases}
\fnull{h,l} \sim \GP (\mu(\cdot), K(\cdot, \cdot)); &  \W{h,l}\sim p^l_\lambda(\W{h,l}) \\
\btheta = c^l(\Xsamples, \W{h,l}); &\fpos{K}{h,l} \mid \btheta,\fpos{0}{h,l} = \Gcompc{h,l}(\fpos{0}{h,l})
\end{rcases}
\end{equation}

where \(\fnull{0} = \Xsamples\), \(c\) is a transformation that generates coefficients for the normalizing flow $\Gcompcinpdep{h,l}{\Xsamples}$, and \(h=1,\dots,H^l\) is the depth of each layer. By this construction, the joint distribution is given by:

\begin{equation}
\begin{split}
p(\{\fK{h,l}\}_{h=1,l=1}^{H^l,L}) &= \myprod{l=1}{L} \myprod{h=1}{H^l}\Ngauss{\fnull{h,l}\mid \mu\left(\{\fK{h,l-1}\}_{h=1}^{H^{l-1}}\right),K\left(\{\fK{h,l-1}\}_{h=1}^{H^{l-1}},\{\fK{h,l-1}\}_{h=1}^{H^{l-1}}\right)} \times\\
\end{split}
\end{equation}

mynode = [draw, circle, minimum size=1.2cm]

[mynode] (gp) at (0, 0) ;

[mynode] (dgp) at (2.2, 0) ;
[->] (gp) –node[midway, above] Concatenation(dgp);

[mynode] (tgp) at (0, -2.2) ;
[->] (gp) –node[midway, rotate=90, above] Normalizing Flow (tgp);

[mynode] (dtgp) at (2.2, -2.2) ;
[->] (gp) –node[midway, rotate=-45, above]  Normalizing Flow + Concatenation   (dtgp);
[->, dashed] (tgp) – (dtgp);
[->, dashed] (dgp) – (dtgp);

(left): intuitive explanation of how our model was built by combining two generalizations of Gaussian Processes.  (right): Example of with 3 layers. In each box, a flow is compounded with the output of a  (colored circle).

<ref> shows the structure of a . Clearly, it can be shown that the proposed architecture is very flexible and contains
as special cases a , a , and a , which can all be recovered by selecting an appropriate number of layers and normalizing flows. The usage of normalizing flows between the layers of the is expected to allow the model to encode prior expert knowledge about the problem, as it already occurs in s [Maroñas et al., 2021]. From now on,
we will consider \(H^l = 1\) for all \(l\) to simplify the notation. The subindex \(\btheta^l_h\) will also be omitted as it can be understood from the context. With these assumptions, the joint prior distribution of a  is given by:

\begin{equation}\label{eq:dtgp:prior}
p(\Ysamples,\deeptransformedset) = \underbrace{\myprod{n=1}{N} p(y_n \mid\Gcompc{L} (\fpos{0,n}{L}))}_{\text{Likelihood}}\underbrace{\myprod{l=1}{L} p(\fnull{l} \mid  \unull{l}) \J{\fK{l} \mid \uK{l}} p(\unull{l}) \J{\uK{l}}}_{\text{ \DTGP Prior}},
\end{equation}

where the diagonal Jacobian appears in each of the layers, representing the transformation of the prior as the novelty from 's prior definition. Since the posterior distribution is intractable, a variational distribution that maintains the exact conditional prior [Titsias, 2009] is chosen following Maroñas et al., 2021:

\begin{equation}\label{eq:dtgp:variational:distr}
q(\deeptransformedset)  = \myprod{l=1}{L} p(\fnull{l} \mid \unull{l}) \J{\fK{l} \mid \uK{l}} q(\unull{l})\J{\uK{l}}.
\end{equation}

This choice of variational distribution allows for term cancellations in the which implies a gain in computational cost. In virtue of the Law of the Unconscious Statistician (), the expression of the is the following:

\begin{equation}\label{eq:elbo:dtgp}
\mathcal L_{\text{\DTGP}} = \underbrace{\mysum{n=1}{N} \Expectation{q(\fnull{L} \mid \fK{L-1})\myprod{l=1}{L-1}q(\fK{l} \mid \fK{l-1})
}{ \log p \left(y_n \mid  \Gcomp^L (\fnull{L})\right)}}_{\text{ELL}} + \underbrace{\mysum{l = 1}{L} \KL{ q(\unull{l})}{ q(\unull{l})}}_{\text{\KLD}}.
\end{equation}

Since \(q(\unull{l})\) and \(p(\unull{l})\) are both Gaussian, the divergence can be computed analytically. The Expected Log Likelihood () term must be approximated since the expectation under the variational distribution is intractable. We approximate this term using Monte Carlo samples from \(\myprod{l=1}{L} q(\fK{l} \mid \fK{l-1}, )\). Thus, we achieve computational tractability using two sources of stochasticity as it was done in [Salimbeni and Deisenroth, 2017]. Firstly, the ELL is a sum across data points, so it can be evaluated using mini-batches. Also, Monte Carlo samples are used to approximate this term. This amounts a computational cost of \(\complexity{(NM^2 + M^3 + NLK)(H^1 + \cdots + H^L)}\) which, compared to the cost of s, adds the cost of computing the normalizing flow \(NLK\). As it will be shown later, the inference algorithm (Alg. <ref>) is a straightforward extension of the one presented in [Salimbeni and Deisenroth, 2017].

§.§ Using Bayesian Priors on Flows

Input Dependent flows compute the parameters of the normalizing flow using a neural network that outputs a set of parameters for each of the flows \(\btheta^l_h = \text{\NN}(\Xsamples, \W{l})\). A Bayesian treatment can be given to these parameters \(\W{l}\). Firstly, it is assumed that the distribution of the parameters of the flows are independent between the layers, that is:

\begin{equation}\label{eq:dtgp:prior:weights:flows}
p(\setWs) = \myprod{l=1}{L} p(\W{l})
\end{equation}

Now, following the observations in [Maroñas et al., 2021], the new joint prior model and the variational posterior factorizes as:

\begin{align}
\begin{split}
p(\Ysamples,\deeptransformedset, \setWs) & = p(\Ysamples,\deeptransformedset \mid \setWs) p(\setWs),\\
q(\deeptransformedset, \setWs) & = q(\deeptransformedset \mid \setWs) q(\setWs).
\end{split}
\end{align}

Using these expressions, the in Equation (<ref>) slightly changes. A new term appears, acting as a regularizer for the distribution of the flow parameters. The complete notation \(\G^L_{\btheta^L(\Xsamples,\W{l})}\) is recovered to remark the dependence of the coefficients obtained by \(\btheta\). Then, expression takes the form:

\begin{align}
\label{eq:dtgp:elbo:flows:prior}
\begin{split}
\mathcal L_{\text{\DTGP}} =  \underbrace{\Expectation{q\left( \deeptransformedset, \setWs \right)}{p(\Ysamples \mid  \Gcompcinpdep{L}{\Xsamples} (\fnull{L}))}}_{\ELL}  & - \underbrace{\mysum{l = 1}{L} \KL{ q(\unull{l})}{ p(\unull{l})}}_{\KLD{}_1} \\
& - \underbrace{\mysum{l = 1}{L} \KL{q(\W{l})}{p(\W{l})}}_{\KLD{}_2}\,.
\end{split}
\end{align}

§.§ Predictions

In this work, we considered a that always uses an identity flow in the last layer \(\Gcompc{L}(\fnull{L}) = \fnull{L}\). First, this allows us to compare directly the modeling advantage that the flows provide in the inner layers compared to identity flows (). Note that since both inference algorithms rely on the same assumptions, we can attribute the performance difference to additional expressiveness provided by the flows, and not to an improved approximation through a better inference algorithm (using [Havasi et al., 2018]).

The aforementioned simplification has another important advantage: the latent function values at the last layer \(\fnull{L}\) remain Gaussian as a consequence of the linear transformation. This simplifies computing expectations w.r.t. the likelihood. More precisely, it leads to a model that has more flexibility than the s due to the normalizing flows between the layers, but that allows closed-form marginalization of the latent function values \(\fnull{L}\) at the last layer (given \(\fK{L-1}\)), as in s. Importantly, note that the only difference between the proposed inference algorithm for s and the algorithm for s relies on the fact that the samples of each of the hidden layers are passed through a non-linearity. Algorithm <ref> shows the evaluation, where the remarked part is the added difference with respect to the .

Having an input \(\X{}\), this input is propagated through the layers and \(S\) Monte Carlo samples are used to approximate \(q(\fpos{K}{L})\), which has the form:

\begin{equation}\label{eq:predict:dtgp:simple:variational:posterior}
q(\fpos{K}{L}) \approx \frac{1}{S} \mysum{s=1}{S} \mathcal N \left(\fpos{K}{L} \mid \varm{}_{qf}\left(\fpos{K,(s)}{L-1}\right), \varS{}_{qf}\left(\fpos{K,(s)}{L-1},\fpos{K,(s)}{L-1}\right)\right),
\end{equation}

where, if we name \(\alpha(\X{}) = K^L(\Zsamples{L}, \Zsamples{L})^{-1}K^L(\Zsamples{}, \X{})\),

\begin{align}
\varm{}_{qf}(\X{}) & = m^L(\X{}) + \alpha(\X{})^T(\varm{L} - m(\Zsamples{L}))\\
\varS{}_{qf}(\X{i},\X{j}) & = K^L(\X{i},\X{j}) -  \alpha(\X{i})^T\left( K^L(\Zsamples{L},\Zsamples{L}) - \varS{L}\right) \alpha(\X{j}),
\end{align}

where \(m^L, K^L, \Zsamples{L}, \varm{L}, \varS{L}\) denote the mean function, kernel, inducing locations, variational mean and variational variance of the layer \(L\), respectively.

The same happens with the predictive distribution for the labels \(y\). We can now assume a Gaussian likelihood to obtain:

\begin{align*}
\begin{split}
p(y \mid \X{}) & = \int q(\fpos{K}{L} \mid \fpos{K}{L-1}) p(y \mid \fpos{K}{L}) \ d \fpos{K}{L}, \\
& = \frac{1}{S}\sum_{s=1}^S \mathcal N(\fpos{K}{L} \mid\varm{}_{qf}\left(\fpos{K,(s)}{L-1}\right), \varS{}_{qf}\left(\fpos{K,(s)}{L-1},\fpos{K,(s)}{L-1}\right) + \pmb{\sigma}^2 \pmb{I}),
\end{split}
\end{align*}

where \(\pmb{\sigma}^2\) is the variance of the likelihood and \(\pmb{I}\) is the identity matrix.

Training Step using a mini batch of size \(B\). We highlight the differences with  algorithm.

Receive a batch \(\mathcal D_B = \{(\X{1}, y_1), \cdots, (\X{B}, y_B)\}\)

\(\fK{0} \gets \mathcal \{\X{1},\dots, \X{B}\}\)

Iterate through layers
\(l = 1,\dots, L\)
Compute marginal predictive distribution
\(\pmb{\mu}_0^l, \pmb{\Sigma}_0^l,  \gets q(\fnull{l} \mid \fK{l-1})\)

Obtain \(S\) samples from the marginal predictive
\(\fnull{l} \gets Sample(\pmb{\mu}_0^l, \pmb{\Sigma}_0^l, S)\)

Input Dependent Flow
\(\btheta^l(\Xsamples, \W{l}) \gets Transformation_l(\Xsamples)\)

Warp the samples using the Normalizing Flow
\(\fK{l} = \G^l_{\btheta^l(\Xsamples, \W{l}}(\fnull{l})\)

Compute ELL
\(\ELL \gets \frac{N}{B} \mysum{b=1}{B} \Expectation{q(\fK{L})}{\log(p(y_b\mid \fK{L})}\)

Compute , analytical form
\(\KLD \gets \mysum{l=1}{L} \KL{q(\unull{l})}{p(\unull{l})} \)

\(\ELBO \gets \ELL - \KLD\)

Optimize parameters using automatic differentiation
\(\btheta \gets Optimizer(\ELBO, \btheta)\)

§ RELATED WORK

The sparse approaches to Gaussian Processes have allowed these models to be computationally tractable when the number of training data points becomes large [Hensman et al., 2013]. Some works have also studied the usage of Harmonic Features in s and their relation to deep models has been studied by Dutordoir et al., 2020.  Deep Gaussian Processes [Damianou and Lawrence, 2013] have been extensively studied in recent literature as a generalization of the s that increase their flexibility.

Also, the usage of Normalizing Flows [Kobyzev et al., 2021] to transform probability density functions has become an active research field [Bånkestad et al., 2023]. Some works warp the  likelihood to increase its flexibility [Snelson et al., 2003, Rios and Tobar, 2019]. In our case, we are more interested in the models that warp the prior distribution of the , leading to an improvement in  the performance of transformed model [Wilson and Ghahramani, 2010, Wauthier and Jordan, 2010, Maroñas et al., 2021, Maroñas and Hernández-Lobato, 2023]. Efficient Transformed Gaussian Processes (s) also increase the speed of standard sparse s methods in multi-class classification tasks [Maroñas and Hernández-Lobato, 2023] so they could be also considered as the base model for concatenation, which we will do in future work.

Our work builds on [Maroñas et al., 2021] and extends it by applying a concatenation of the Transformed Gaussian Process following the fashion of the previously presented s. This same idea has been used in different works that use different generalizations of s to define more complex models. In [Ortega et al., 2023], the concatenated model is the Implicit Processes (), and in [Rudner et al., 2020] the way the covariance is computed is changed, leading to a more efficient model than s. Our work differs from the previous ones in the form of extending the s by using Normalizing flows and also having the possibility of using input-dependent transformations. Since inference remains intractable in our model, we follow a similar approach to [Salimbeni and Deisenroth, 2017], taking samples from the posterior that in our case are warped using the normalizing flow.

§ EXPERIMENTS

s have been carefully examined in an environment where s do not perform at their best. To this end, [Rudner et al., 2020] presents a toy dataset with a step-wise function where the s have difficulties due to the smoothness that the RBF kernel imposes in the prior. <ref> shows the visual comparison of s versus s in this toy dataset. s use the Steptanh flow, given by a linear combination of the hyperbolic tangent function [Maroñas et al., 2021, Maroñas and Hernández-Lobato, 2023]. We observe that the proposed model shows not only a better mean prediction of the data but also a good uncertainty estimation where the function value changes. Also, the functional form of the normalizing flow  in <ref> confirms that using expert information about the problem in the normalizing flows (in this example, a flow that induces a step) can be very helpful for s.

On the left, 2 layer , which has difficulties in modeling the jumps between the steps the data makes. On the right, the proposed  (2 layers) which is capable of fitting the data  better thanks to the nonlinearity in the form of a Steptanh Normalizing flow.

The last experiments using s have been conducted using 8 real datasets from the UCI repository. In particular, the datasets used in [Salimbeni and Deisenroth, 2017] have been chosen, exchanging the Naval dataset by Yacht. The goal of these experiments is to compare the s with the s, again with the objective of testing if adding normalizing flows in between the layers results in better performance.
As per usual, a \(10\%\) test size has been chosen and each experiment has been performed using \(20\) different random seeds, averaging the results. The initializations of the variational and kernel parameters have followed the ones in [Salimbeni and Deisenroth, 2017], to replicate their experimental environment.  The number of layers used is 2, 3, 4 and 5 for each of the models, using \(M = 100\) inducing points in each . All the models are trained for \(80.000\) iterations using a fixed batch size of \(200\) and a learning rate of \(10^{-2}\). Regarding normalizing flows, it must be remarked that in this work all the layers share the same functional form with different trainable parameters in each layer. Also, we have fixed noninput dependent arcsinh flow as the used normalizing flow, since it performed well in some initial experiments. Lastly, \(H_l = 1\) for every \(l=1,\dots,L\).

Using the configuration mentioned above, the results obtained in terms of Negative Log Likelihood () for each method are shown in Figure <ref>. The dot indicates the mean of the \(20\) splits, and the bars indicate the standard deviation divided by \(\sqrt{20}\). It can be observed that  s achieve better scores in two of the datasets (Kin8nm and Power), slightly better results in mean in other two (Boston, Protein), comparable results in three datasets (Energy, Redwine and Yacht) and slightly worse results (in some cases comparable) in Concrete. Numeric values are shown in Table <ref>.
These results also can give some more insight on s. While most of the time increasing the number of layers improves the performance in s, in s sometimes the performance becomes worse. This points out again the difficulty of training the s since the performance should normally be maintained as the number of layers is increased. Also, the increase in expressiveness may lead to data overfitting, which was already observed in s [Maroñas et al., 2021]. Further work includes Bayesian flows as a way to prevent this overfitting.

Negative Log Likelihood results (lower values, to the left, are better) comparing s with noninput dependent Arcsinh Flow with s, both using 2,3,4, and 5 Layers.

NLL             2cboston           2cconcrete             2cenergy              2ckin8nm              2cpower            2cprotein            2credwine              2cyacth

2-     2.442     0.049      3.150     0.016      0.718     0.02     -0.968     0.004     2.821     0.007     2.095     0.011      0.950     0.014     0.544     0.112

3-      2.453     0.050     3.143     0.016      0.725     0.020     -0.964     0.004      2.820     0.007      2.080     0.006      0.950     0.014     0.495     0.103

4-     2.521     0.071     3.148     0.016       0.720     0.020     -0.966     0.004     2.819     0.007     2.079     0.006     0.949     0.014      0.516     0.110

5-     2.412     0.045     3.143     0.016      0.714     0.020     -0.967     0.005     2.816     0.007     2.078     0.006     0.948     0.0140     0.464     0.071

2-      2.370     0.045     3.165     0.021     0.716     0.027     -0.998     0.004     2.795     0.007     2.069     0.007     0.957     0.015     0.491     0.141

3-     2.349     0.045      3.166     0.020     0.727     0.024     -0.997     0.004     2.798     0.007     2.068     0.007     0.957     0.014     0.434     0.107

4-      2.359     0.05     3.175     0.022      0.730     0.023     -0.994     0.004     2.799     0.007      2.07     0.007     0.954     0.014     0.523     0.127

5-      2.350     0.043     3.194     0.023     0.739     0.024     -0.994     0.005     2.799     0.007      2.070     0.007     0.955     0.014     0.643     0.127

Negative Log Likelihood exact results.

Another aspect that must be carefully examined is the time used by each method. To measure this, we measure the elapsed time in the first 1000 iterations.
The results are shown in Figure <ref>. We observe that the training time in seconds becomes more significant as the number of layers grows.
Specifically, a  with 4 layers needs almost the same computational time as a 5 layer .

CPU time elapsed by each model in 1000 iterations, averaged in 20 splits.

§ CONCLUSIONS

We have presented Deep Transformed Gaussian Processes (s), a model based on the concatenation of Transformed Gaussian Processes (s). s increase the expressiveness of both s and s, by transforming the predictive distribution in each layer. Also, we have derived a further extension of the model where Bayesian priors are placed on the transformation that computes the flow's parameters.

s inherit the intractability of their base models. Due to this variational inference is used to find an approximation to its true posterior distribution. The derivation of the Evidence Lower Bound for s has been presented. Furthermore, we have shown how to evaluate this lower bound using Monte Carlo samples. In the performed experiments in the toy data, the improvements that the nonlinearities of the s offer have been shown using different types of normalizing flows. This improvement is reflected in the real datasets, where the proposed implementation of the model achieves better or at worst comparable results to the model.

§ ACKNOWLEDGMENTS

The authors gratefully acknowledge the use of the facilities of Centro de Computacion Cientifica (CCC) at Universidad Autónoma de Madrid.  The authors also acknowledge financial support from the Spanish Plan Nacional I+D+i, PID2019-106827GB-I00 and PID2022-139856NB-I00, and from the Autonomous Community of Madrid
(ELLIS Unit Madrid).

§ DTGP DERIVATIONS

This is a complementary section in which the derivations of the Deep Transformed Gaussian Process model are presented.

§ DTGP PRIOR

The DTGP prior in Equation (<ref>) is obtained as follows:

\begin{align*}
\setlength{\jot}{10pt}
\begin{WithArrows}
p(\Ysamples,\deeptransformedset) & = p(\Ysamples \mid \deeptransformedset) p( \deeptransformedset) \Arrow {Conditional independence \\ + layer factorization}\\
& = p(\Ysamples{} \mid \fK{L}) \myprod{l=1}{L} p(\fK{l}, \uK{l}) \Arrow {Change of Variable} \\
& = p(\Ysamples \mid \Gcompc{L} (\fnull{L}))\myprod{l=1}{L} p(\fnull{l}, \unull{l}) \J{\fK{l}, \uK{l}} \\
& = \myprod{n=1}{N} p(y_n \mid \Gcompc{L} (\fpos{0,n}{L}))\myprod{l=1}{L} p(\fnull{l} \mid  \unull{l}) p(\unull{l}) \J{\fK{l}, \uK{l}}\\
& = \underbrace{\myprod{n=1}{N} p(y_n \mid\Gcompc{L} (\fpos{0,n}{L}))}_{\text{Likelihood}}\underbrace{\myprod{l=1}{L} p(\fnull{l} \mid  \unull{l}) p(\unull{l}) \J{\fK{l}, \uK{l}}}_{\text{ DTGP Prior}}
\end{WithArrows}
\end{align*}

§ EVIDENCE LOWER BOUND IN DTGPS

Recalling that when using marginal flows \(\J{\fK{}}\J{\uK{}} = \J{\fK{},\uK{}}\), the ELBO in Equation (<ref>) is derived as follows:

\[
\setlength{\jot}{20pt}
\begin{WithArrows}
\mathcal L(\Ysamples) & = \Expectation{q\left( \deeptransformedset\right)} { \log \frac{p\left(\Ysamples, \deeptransformedset\right)}{q\left(\deeptransformedset\right)}}\\
& =  \Expectation{q\left( \deeptransformedset\right)} { \log
\frac{p(\Ysamples \mid \Gcompc{L} (\fnull{L})) \myprod{l=1}{L} \cancel{p(\fnull{l} \mid  \unull{l})} p(\unull{l})\cancel{ \J{\fK{l}, \uK{l}}}}
{\myprod{l=1}{L} \cancel{p(\fnull{l} \mid \unull{l} ) } q(\unull{l}) \cancel{\J{\fK{l},\uK{l}}}  } } \\
& = \Expectation{q\left( \deeptransformedset\right)} { \log p(\Ysamples \mid \Gcompc{L} (\fnull{L}) )} + \mysum{l=1}{L}  \Expectation{q\left(\fK{l}, \uK{l}\right)} {\log \frac{p(\unull{l}) }{ q(\unull{l}) }} \Arrow{LOTUS} \\
& = \Expectation{q\left( \deeptransformedset\right)}{ \log p(\Ysamples \mid \Gcompc{L} (\fnull{L}))} + \mysum{l=1}{L}  \Expectation{q\left(\fnull{l}, \unull{l}\right)}{\log \frac{p(\unull{l}) }{ q(\unull{l}) }} \\
& = \underbrace{\Expectation{q\left( \deeptransformedset\right)}{ \log p(\Ysamples \mid \Gcompc{L} (\fnull{L}) )}}_{\text{Expected Log Likelihood}} - \underbrace{\mysum{l = 1}{L} \KL{ q(\unull{l})}{ p(\fnull{l})}}_{\KLD}
\end{WithArrows}
\]

The Expected Log likelihood term (ELL) requires further development. We operate with the expectation as an integral. Recalling that \( p(\Ysamples \mid \Gc{L}(\fnull{L}))\) factorizes over the data, we can observe that, naming \(\ELL =  \Expectation{Q}{\log p(\Ysamples \mid \Gcompc{L} (\fnull{L}) )}\)

\begin{align*}
\ELL & = \int q(\deeptransformedset) \log \myprod{n=1}{N}  p(\Ysamples_n \mid  \Gcompc{L} (\fnull{L})) \dd \fK{1,\dots,L} d \uK{1,\dots,L} \\
& = \int  \log \left( \myprod{n=1}{N} p(\Ysamples_n \mid  \Gcompc{L} (\fnull{L})) \right)\left\{\int \myprod{l=1}{L} p(\fK{l} \mid \uK{l}; \fK{l-1}, \Zsamples{l}) q(\uK{l}) d \uK{1,\dots,L} \right\}\dd \fK{1,\dots,L} \\
& = \int  \log  \left(\myprod{n=1}{N} p(\Ysamples_n \mid  \Gcompc{L} (\fnull{L})) \right)\left\{\int \myprod{l=1}{L} p(\fnull{l} \mid \unull{l}; \fK{l-1}, \Zsamples{l}) q(\unull{l}) \J{\fK{}} \dd \unull{1,\dots,L} \right\}d \fK{1,\dots,L} \\
& = \int \log \left(\myprod{n=1}{N} p(\Ysamples_n \mid  \Gcompc{L} (\fnull{L})) \right)\left\{ \myprod{l=1}{L}\int  p(\fnull{l} \mid \unull{l}; \fK{l-1}, \Zsamples{l}) q(\unull{l}) \J{\fK{}} \dd \unull{l} \right\}d \fK{1,\dots,L} \\
& = \int \log \left(\myprod{n=1}{N} p(\Ysamples_n \mid  \Gcompc{L} (\fnull{L})) \right)\left\{ \myprod{l=1}{L}\int  q(\fnull{l},\unull{l} \mid \fK{l-1}, \Zsamples{l}) \J{\fK{}} d \unull{1,\dots,L} \right\}\dd \fK{1,\dots,L} \\
& = \int \log \left(\myprod{n=1}{N} p(\Ysamples_n \mid  \Gcompc{L} (\fnull{L}))  \right)\myprod{l=1}{L}  q(\fK{l} \mid \fK{l-1}, \Zsamples{l}) d \fK{1,\dots,L} \\
& = \mysum{n=1}{N} \int \log p\left(\Ysamples_n \mid  \Gcompc{L} (\fnull{L}) \right) \myprod{l=1}{L}  q(\fK{l} \mid \fK{l-1}, \Zsamples{l}) d \fK{1,\dots,L} \\
& = \mysum{n=1}{N} \Expectation{\myprod{l=1}{L} q(\fK{l} \mid \fK{l-1}, \Zsamples{l})}{ \log p\left(\Ysamples_n \mid  \Gcompc{L} (\fnull{L})\right)}
\end{align*}

Where, to get from the second to the third equations we have applied the LOTUS rule and, from the third to the fourth equation it is used that each integral only depends on one of the factors.

§ EXPECTED LOG LIKELIHOOD IN DTGP USING BAYESIAN PRIORS

We can expand the ELL term in Equation (<ref>) following the same procedure done in the previous case. In this case, since the flow depends on the sampled weights, we will denote it as \(\Gcompcinpdep{l}{\Xsamples}\) with \(l = 1,\dots, L\). The derivations proceed as follows:

\begin{align*}
\ELL & = \int q(\deeptransformedset, \setWs) \log \myprod{n=1}{N}  p(\Ysamples_n \mid  \Gcompcinpdep{L}{\Xsamples}  (\fnull{L})) \dd  \fK{1,\dots,L} d \uK{1,\dots,L} d \W{1,\dots,L}\\
& = \int q(\deeptransformedset \mid \setWs) q(\setWs) \log \myprod{n=1}{N}  p(\Ysamples_n \mid  \Gcompcinpdep{L}{\Xsamples} (\fnull{L})) \dd  \fK{1,\dots,L} d \uK{1,\dots,L} d\W{1,\dots,L}\\
& = \int q(\setWs)\log \myprod{n=1}{N} p(\Ysamples_n \mid  \Gcompcinpdep{L}{\Xsamples} (\fnull{L}))  \left(\int  q(\deeptransformedset \mid \setWs)   d \uK{1,\dots,L}\right)  \dd \fK{1,\dots,L}  d\W{1,\dots,L}\\
& = \int q(\setWs)\log \myprod{n=1}{N} p(\Ysamples_n \mid \Gcompcinpdep{L}{\Xsamples} (\fnull{L}))  \left( \myprod{l=1}{L}q(\fK{l} \mid \fK{l-1}, \W{l})\right)  d \fK{1,\dots,L}  d\W{1,\dots,L}\\
& = \int \log \left(\myprod{n=1}{N}  p(\Ysamples_n \mid  \Gcompcinpdep{L}{\Xsamples} (\fnull{L}))\right) \ \myprod{l=1}{L} q(\W{l}) q(\fK{l} \mid \fK{l-1}, \W{l}) \dd \fK{1,\dots,L}  d\W{1,\dots,L}\\
& = \mysum{n=1}{N} \int \log \left( p(\Ysamples_n \mid \Gcompcinpdep{L}{\Xsamples} (\fnull{L}))\right) \ \myprod{l=1}{L} q(\W{l}) q(\fK{l} \mid \fK{l-1}, \W{l}) \dd \fK{1,\dots,L}  d\W{1,\dots,L}\\
& = \mysum{n=1}{N}\Expectation{\myprod{l=1}{L} q(\W{l}) q(\fK{l} \mid \fK{l-1}, \W{l})}{\log p\left( \Ysamples_n \mid  \Gcompcinpdep{L}{\Xsamples} (\fnull{L})\right)}\\
& \stackrel{(1)}{\approx} \frac{1}{S} \mysum{s=1}{S}\mysum{n=1}{N} \Expectation{\myprod{l=1}{L}  q(\fK{l} \mid \fK{l-1}, \W{l}_s)} {\log  p\left( \Ysamples_n \mid  \Gk{\btheta(\W{L}_{s},\Xsamples)}^{L}(\fnull{L})\right)}
\end{align*}

where the approximation in \((1)\) is done using \(S\) Monte-Carlo Samples \(\W{}_s \sim q(\W{})\).

[Bånkestad et al., 2023]
M. M. Bånkestad, J. Sjölund, J. Taghia, and T. B. Schön.
Variational elliptical processes.
Transactions on Machine Learning Research, 2023.

[Damianou and Lawrence, 2013]
A. Damianou and N. D. Lawrence.
Deep Gaussian processes.
In Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics, 2013.

[Dutordoir et al., 2020]
V. Dutordoir, N. Durrande, and J. Hensman.
Sparse gaussian processes with spherical harmonic features.
In Proceedings of the 37th International Conference on Machine Learning, 2020.

[Havasi et al., 2018]
M. Havasi, J. M. Hernández-Lobato, and J. J. Murillo-Fuentes.
Inference in deep gaussian processes using stochastic gradient hamiltonian monte carlo.
Advances in neural information processing systems, 2018.

[Hensman et al., 2013]
J. Hensman, N. Fusi, and N. D. Lawrence.
Gaussian processes for big data.
In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, 2013.

[Kobyzev et al., 2021]
I. Kobyzev, S. D. Prince, and M. A. Brubaker.
Normalizing flows: An introduction and review of current methods.
IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021.

[Maroñas and Hernández-Lobato, 2023]
J. Maroñas and D. Hernández-Lobato.
Efficient transformed gaussian processes for non-stationary dependent multi-class classification.
In International Conference on Machine Learning, 2023.

[Maroñas et al., 2021]
J. Maroñas, O. Hamelijnck, J. Knoblauch, and T. Damoulas.
Transforming gaussian processes with normalizing flows.
In The 24th International Conference on Artificial Intelligence and Statistics, AISTATS, 2021.

[Ortega et al., 2023]
L. A. Ortega, S. R. Santana, and D. Hernández-Lobato.
Deep variational implicit processes.
In The Eleventh International Conference on Learning Representations, 2023.

[Rasmussen and Williams, 2005]
C. E. Rasmussen and C. K. I. Williams.
Gaussian Processes for Machine Learning.
The MIT Press, 2005.

[Rezende and Mohamed, 2015]
D. Rezende and S. Mohamed.
Variational inference with normalizing flows.
In Proceedings of the 32nd International Conference on Machine Learning, 2015.

[Rios and Tobar, 2019]
G. Rios and F. Tobar.
Compositionally-warped gaussian processes.
Neural Networks, 2019.

[Rudner et al., 2020]
T. G. Rudner, D. Sejdinovic, and Y. Gal.
Inter-domain deep gaussian processes.
In International Conference on Machine Learning, 2020.

[Salimbeni and Deisenroth, 2017]
H. Salimbeni and M. Deisenroth.
Doubly stochastic variational inference for deep gaussian processes.
In Advances in Neural Information Processing Systems, 2017.

[Schmidt et al., 2023]
A. Schmidt, P. Morales-Álvarez, and R. Molina.
Probabilistic attention based on gaussian processes for deep multiple instance learning.
IEEE Transactions on Neural Networks and Learning Systems, 2023.

[Snelson et al., 2003]
E. Snelson, Z. Ghahramani, and C. Rasmussen.
Warped gaussian processes.
In Advances in Neural Information Processing Systems, 2003.

[Titsias, 2009]
M. Titsias.
Variational learning of inducing variables in sparse gaussian processes.
In Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, 2009.

[Wauthier and Jordan, 2010]
F. L. Wauthier and M. Jordan.
Heavy-tailed process priors for selective shrinkage.
In Advances in Neural Information Processing Systems, 2010.

[Wilson and Ghahramani, 2010]
A. G. Wilson and Z. Ghahramani.
Copula processes.
In Advances in Neural Information Processing Systems, 2010.


# A walk through of time series analysis on quantum computers

Ammar Daskin
Dept. of Computer Engineering
Istanbul Medeniyet University, Istanbul, Turkey<EMAIL_ADDRESS>
(Received: date / Accepted: date)

###### Abstract

Because of the rotational components on quantum circuits, some quantum neural
networks based on variational circuits can be considered equivalent to the
classical Fourier networks. In addition, they can be used to predict the
Fourier coefficients of continuous functions. Time series data indicates a
state of a variable in time. Since some time series data can be also
considered as continuous functions, we can expect quantum machine learning
models to do many data analysis tasks successfully on time series data.
Therefore, it is important to investigate new quantum logics for temporal data
processing and analyze intrinsic relationships of data on quantum computers.

In this paper, we go through the quantum analogues of classical data
preprocessing and forecasting with ARIMA models by using simple quantum
operators requiring a few number of quantum gates. Then we discuss future
directions and some of the tools/algorithms that can be used for temporal data
analysis on quantum computers.

###### Index Terms:

quantum machine learning, quantum algorithms, quantum optimization, quantum
time series, quantum temporal data,

##  I Introduction

It is shown that nature itself is explained better in quantum physics than
classical one. The idea of quantum computers was put forward by Feynman [1]
and other pioneers [2] so as to simulate and understand intrinsic physical
properties of nature in a better way: Simulating quantum systems with many
particles is difficult on classical computers because the computational
complexity in general grows exponentially with the number of involved
particles. On the other hand, quantum computers are based on qubits
implemented by quantum particles which in turn indicates that in principle, we
can increase the computational power of quantum computers exponentially by
adding controllable new qubits. Therefore, the difficulty of simulating nature
for classical computers would become the main advantage of quantum computers.
Unfortunately since the input/output of quantum computers are still in
classical ones and zeros, it necessitates a state preparation that involves
classical processing of the data and a measurement that requires many
repetitions of the same computation. Therefore, the complexity boundaries in
classical computation mostly remain the same for quantum computing.
Nonetheless, it is still possible to design more efficient algorithms by
employing quantum entanglement and superposition (quantum parallelism): Some
known examples are Shor’s integer factoring algorithm [3], Grover’s search
algorithm in an unordered list [4], and the quantum principal component
analysis (PCA) [5].

###  I-A Classical neural networks

The representation theorem indicates any multivariate continuous function can
be represented as a superposition of functions of one variable and addition
[6, 7]. An artificial neural network can be considered as a composition of
simple nonlinear activation function[8]: Let $\displaystyle I^{n}$ be the
n-dimensional unit-cube $\displaystyle[0,1]^{n}$ and $\displaystyle C(I^{n})$
be the space of continuous functions on $\displaystyle I^{n}$. In addition,
let $\displaystyle\sigma$ be a sigmoidal function. Given any $\displaystyle
f\in C(I^{n})$ and $\displaystyle\epsilon>0$, there is a $\displaystyle G$ of
the following form:

$G(x)=\sum_{j=1}^{n}\alpha_{j}\sigma(w_{j}^{T}x+\theta_{j})$ (1)

for which $\displaystyle|G(x)-f(x)|<\epsilon$.

Gallant et al. [9] showed a single feed-forward neural network where the
output is squashed by using a cosine activation function: i.e. by considering
$\displaystyle\sigma$ as cosine, sine, or exponential (i.e. in the form
$\displaystyle e^{iw_{j}^{T}x+\theta_{j}}$) functions, one can approximate the
Fourier transform of a continuous function. These networks are often called
Fourier neural networks (FNN). On the other hand, Fourier analysis of data
also sheds lights on the generalization performances of deep neural
networks[10].

###  I-B Time series problems

Temporal data indicates a state of a variable in time, which could be
collected through a wide variety of monitoring devices. Examples include
financial market data; sensor data which may indicate temperature, humidity
and so on; and medical data collected by different monitoring medical devices.
Depending on the collection method, the data can be continuous or discrete.
Continuous temporal data sets are generally referred to as times series data.
[11] The data mining tasks on time series data include clustering,
classification, and forecasting: In general it is required to specialize the
known algorithms in data analysis in order to use in the analysis of time
series data: e.g. see the classification with deep learning for time series
data [12] and the distance based time series classification [13]. It is shown
that time series data can be also modeled by using Fourier deep neural
networks [14], which is beriefly explained above.

###  I-C Quantum machine learning

Since the quantum machine learning [15, 16] is a very new field, there are
only a small number of partial studies which can be found for the analysis of
time series data on quantum computers: the time series data analysis is done
by using a quantum approach based on the Schrödinger’s equation [17]. In a
very recent paper, the quantum machine learning approaches are investigated
for the problems in finance [18]. It is shown that the variational quantum
machine learning models can be used to approximate Fourier transforms of the
continuous functions [19]. It is also shown that the continuous variable (CV)
quantum neural networks can be used for curve fitting and function
approximations[20]. It is still important to study and build generic separate
models and packages for temporal data analysis.

Quantum circuits are in general based on quantum gates described by some form
of rotations in the Hilbert space. For instance, the rotation gates for one-
qubit can be defined in the following forms;

$R_{x}(2\theta)=\left(\begin{matrix}\cos(\theta)&-i\sin(\theta)\\\
-i\sin(\theta)&\cos(\theta)\end{matrix}\right),$ (2)

$R_{y}(2\theta)=\left(\begin{matrix}\cos(\theta)&-\sin(\theta)\\\
\sin(\theta)&\cos(\theta)\end{matrix}\right),$ (3)

$\text{and }R_{z}(2\theta)=\left(\begin{matrix}e^{-i\theta}&0\\\
0&e^{i\theta}\end{matrix}\right).$ (4)

Because of these rotational components, the quantum circuits can in theory
provide a natural equivalent to the classical Fourier networks (e.g., see our
previous paper[21] which gives an example quantum neural net with periodic
activation function.). In addition, they can be used to predict Fourier
coefficients of continuous functions as done in Ref.[20]

###  I-D Quantum temporal data analysis

Since some time series data can also be considered as continuous functions,
quantum machine learning models can successfully do many data analysis tasks
on the time series data. Therefore, it is important to investigate new quantum
logics for temporal data processing and analyze intrinsic relationships of
data on quantum computers. The processing data on quantum computers in general
requires the followings:

  * •

classical preprocessing and quantum state preparation,

  * •

quantum preprocessing,

  * •

quantum processing of data,

  * •

quantum measurement,

  * •

and analyses of the measurement results on classical computers

In general, all or some of these steps are repeated until a desired
result/model is obtained in data analysis. For continuous temporal data, these
steps are visualized in Fig. 1.

Figure 1:  General picture for processing temporal data on quantum computers
(the measurement step is not shown.)

In the following sections, we go through the quantum analogues of classical
data preprocessing and forecasting. Then we discuss some of the
tools/algorithms that can be used for temporal data analyses on quantum
computers and discuss future directions.

##  II Preparing data for quantum algorithms

The time series data is generally given in contextual representations where
attributes describe the context and behavioral values at any reference points
[11]. If more than one behavioral attribute is associated with the series,
then the data is referred as the multivariate time series: i.e. A time series
data with dimension $\displaystyle d$ and length $\displaystyle n$. The data
is formed by a set of the values received at the time stamp $\displaystyle
t_{i}$. In the $\displaystyle i$th time stamp, the received data is a vector
of dimension $\displaystyle d$:
$\displaystyle\mathbf{y_{i}}=(y_{i}^{1},\dots,y_{i}^{d})$.

In classical data analyses, the time series data is taken under preprocessing
operations: e.g., the missing values may be handled,the data may be
normalized, smoothed, or scaled. Then, a variety of methods are used to
transform and reduce the dimension of the data. After these steps, the
standard data analysis tasks such as classification and clustering are also
applicable to time series data. The classification and clustering tasks
involve determining the correlation across the series and classifying the
series based on their shapes. These can be offline or real time. And class
labels may need to be assigned to the individual time stamps (point labels) or
the whole series.

For quantum computers, the algorithms either use the data after some classical
preprocessing, or as done in PCA [5], it assumes the data is processed and
stored in quantum RAM [22]. In the former case, the data is prepared as a
quantum state by a state preparation circuit that converts a standard initial
state:e.g. $\displaystyle\left|0\right\rangle$, into the state that represents
the data $\displaystyle\left|x\right\rangle$. If the length of
$\displaystyle\left|x\right\rangle$ is $\displaystyle n$, then this kind of
state preparation requires a circuit with $\displaystyle O(n)$ number of
controlled-NOT quantum gates in general [23]. In some cases, the data is also
used as parameters that determine the rotation angle of the quantum gates(e.g.
in [24], the data is used with parameters to define quantum rotation gates.
The learning task is done by reusing the data and adjusting the value of the
parameters for quantum gates.). After the mapping,
$\displaystyle\left|x\right\rangle$ is a vector defined in the Hilbert space.
It is shown in Ref.[25] that the Hilbert space can be considered as a kernel
space used in classical machine learning to ease the training process.
Therefore, the training of a learning model can be achieved also in the
Hilbert space. Although this solves the representability of the data in the
Hilbert space, it is still necessary to design efficient data mapping methods
without disrupting the representability of the data (e.g. for MNIST data, it
is reduced to 8 qubits by using classical PCA [26] which reduces the training
results to as low as 70%. ).

Consequently, the representation is highly dependent on the quantum algorithm
to be used in the data processing step. However, in most cases, since the time
series data is multivariate in the form
$\displaystyle\mathbf{y_{i}}=(y_{i}^{1},\dots,y_{i}^{d})$, the quantum
superposition approach can be used to efficiently represent the data on
quantum computers. This in general may reduce the $\displaystyle d\times n$
data into an $\displaystyle n$-dimensional quantum state:

$\left|\mathbf{y}\right\rangle=\frac{1}{\sqrt{n}}\sum_{i}^{n}\left|\mathbf{y_{i}}\right\rangle$
(5)

Alternatively, one can build the following quantum state from the input data:

$\left|\mathbf{y}\right\rangle=\left(\begin{matrix}\mathbf{y_{1}}\\\
\mathbf{y_{2}}\\\ \vdots\\\ \mathbf{y_{n}}\end{matrix}\right)$ (6)

Here, the dimension of the quantum state is $\displaystyle\log(nd)$. Another
alternative can be the following:

$\left|\mathbf{y}\right\rangle=\left|\mathbf{y}_{1}\right\rangle\left|\mathbf{y}_{2}\right\rangle\dots\left|\mathbf{y}_{n}\right\rangle.$
(7)

Although the required number of qubit which is $\displaystyle n\log(d)$ is
much more than the first two methods, this representation may allow us to
apply quantum operations to individual set of values more easily. However,
since the required number of qubits depend on $\displaystyle n$, it may not be
possible on current quantum computers to use this representation for the
series with large $\displaystyle n$.

##  III Quantum preprocessing of temporal data

Below we show how to implement quantum analogues of some of the famous
classical preprocessing techniques: the exponential smoothing, binning and
data transformation. Note that in most cases, we go through the steps for
univariate data.

###  III-A Exponential smoothing

In the exponential smoothing, the smoothed value $\displaystyle
y^{\prime}_{i}$ is obtained through the combination of $\displaystyle y_{i}$
and the pervious values with the ratios determined by the decay factor
$\displaystyle\alpha$:

$y^{\prime}_{i}=\alpha y_{i}+(1-\alpha)y^{\prime}_{i-1}.$ (8)

The recursive substitution in the above equation yields:

$y^{\prime}_{i}=(1-\alpha)^{i}y^{\prime}_{0}+\alpha\sum_{j=1}^{i}y_{j}(1-\alpha)^{i-j}.$
(9)

In the analyses, the choice of $\displaystyle\alpha$ determines the impact of
the previous data points.

In quantum analyses, we can construct a similar exponential smoothing
technique:

  1. 1.

Different from the classical analysis,
$\displaystyle\left|\mathbf{y_{i}}\right\rangle$, the data observed at the
$\displaystyle i$th time stamp, can be a single data point or multiple data
points observed on $\displaystyle d$ dimension. Therefore, as a physical
qubit, this state can be implemented as a $\displaystyle d$-dimensional
quantum system with $\displaystyle O(\log(d))$ qubits or a $\displaystyle
d$-level qubit.

  2. 2.

The new observed data at time stamp $\displaystyle(i+1)$ is prepared as a
quantum state $\displaystyle\left|\mathbf{y_{i+1}}\right\rangle$ on different
qubits: i.e.,
$\displaystyle\left|\mathbf{y_{i+1}}\right\rangle$$\displaystyle\left|\mathbf{y_{i}}\right\rangle$).
Or we can create the following quantum state by using a control qubit and data
qubit.

$\left(\begin{matrix}\left|\mathbf{y_{i+1}}\right\rangle\\\
\left|\mathbf{y_{i}}\right\rangle)\end{matrix}\right)$ (10)

Note that this state can be formed with the help of control qubit by applying
the operators related to $\displaystyle\left|\mathbf{y_{i+1}}\right\rangle$
and $\displaystyle\left|\mathbf{y_{i}}\right\rangle$) to the data qubit in
$\displaystyle\left|\mathbf{0}\right\rangle$ state.

  3. 3.

From the above state, a quantum entangling gate (e.g. the Hadamard gate) can
be used to generate the following recurrence:

$\begin{split}\left|\mathbf{y^{\prime}_{i+1}}\right\rangle=&\sqrt{\alpha}\left(\left|\mathbf{y_{i+1}}\right\rangle+\left|\mathbf{y^{\prime}_{i}}\right\rangle\right)+\sqrt{1-\alpha}\left(\left|\mathbf{y_{i+1}}\right\rangle-\left|\mathbf{y^{\prime}_{i}}\right\rangle\right)\\\
\\\
=&a\left|\mathbf{y_{i+1}}\right\rangle+b\left|\mathbf{y^{\prime}_{i}}\right\rangle\end{split}$
(11)

where $\displaystyle a$ $\displaystyle=$
$\displaystyle\left(\sqrt{\alpha}+\sqrt{1-\alpha}\right)$ and $\displaystyle
b$ $\displaystyle=$ $\displaystyle\left(\sqrt{\alpha}-\sqrt{1-\alpha}\right).$
In the case of the Hadamard gates, $\displaystyle\alpha$ is some power of
$\displaystyle 1/2$. Using the above recursion, if we substitute
$\displaystyle\left|\mathbf{y^{\prime}_{i}}\right\rangle$ into the above,
$\displaystyle\left|\mathbf{y^{\prime}_{i+1}}\right\rangle$ becomes:

$\begin{split}\left|\mathbf{y^{\prime}_{i+1}}\right\rangle=&a\left|\mathbf{y_{i+1}}\right\rangle+b\left(a\left|\mathbf{y_{i}}\right\rangle+b\left|\mathbf{y^{\prime}_{i-1}}\right\rangle\right)\\\
=&a\left|\mathbf{y_{i+1}}\right\rangle+ab\left|\mathbf{y_{i}}\right\rangle+ab^{2}\left|\mathbf{y_{i-1}}\right\rangle\dots\\\
=&ab^{i+1}\left|\mathbf{y^{\prime}_{0}}\right\rangle+\sum_{k=1}^{i}ab^{i-k}\left|\mathbf{y_{k+1}}\right\rangle,\end{split}$
(12)

If we choose an $\displaystyle\alpha$ which makes $\displaystyle b<1$, we can
have the quantum equivalent of the classical exponential smoothing. In this
case, the choice of $\displaystyle\alpha$ regulates the decay of the terms
through $\displaystyle b$; therefore, we can call $\displaystyle b$ as the
decay factor that determines the importance of the recent data points in any
analysis using $\displaystyle\left|\mathbf{y^{\prime}_{i}}\right\rangle$.

###  III-B Binning and moving average smoothing

The classical moving average is a special case of the binning method where the
date is divided into equally spaced time intervals of size $\displaystyle k$:
$\displaystyle[t_{1},t_{k}],[t_{k+1},t_{2k}]$. Then, the data points are
assigned into the corresponding bins and the mean values of the bins are used
for the analyses:

$y^{\prime}_{i+1}=\frac{\sum_{r=1}^{k}y_{i.k+r}}{k}$ (13)

Here, the number of the data points are reduced by a factor of $\displaystyle
k$.

Utilizing the Hadamard gates on certain qubits, one can obtain the average of
the neighboring states on certain quantum states. After finding the average,
we can also disregard the rest of the state by collapsing quantum state onto
the interested part of the system: One example could be as follows:

  * •

Let
$\displaystyle\left|\mathbf{y_{i}}\right\rangle=\left(\begin{matrix}y_{i1}\\\
\vdots\\\ y_{id}\end{matrix}\right)$. If we multiply this quantum state
$\displaystyle U_{bin}=I^{\otimes\log_{2}(d/k)}\otimes H^{\otimes\log_{2}k}$,
where $\displaystyle H$ is the Hadamard gate:

$H=1/\sqrt{2}\left(\begin{matrix}1&1\\\ 1&-1\end{matrix}\right)$ (14)

If we apply $\displaystyle U_{bin}$ to the quantum state
$\displaystyle\left|\mathbf{y_{i}}\right\rangle$, then when the first
$\displaystyle\log_{2}k$ qubits in $\displaystyle\left|0\dots 0\right\rangle$
state, we would have an equivalent $\displaystyle\mathbf{y^{\prime}_{i+1}}$
which could be considered reduced data points.

The classical moving average smoothing uses overlapping bins which can be
easily done on quantum case by utilizing a quantum operator of the form:
$\displaystyle H^{\otimes\log_{2}k}\otimes I^{\otimes\log_{2}(d/k)}$.

###  III-C Data transformation

In classical data analyses, often the temporal data is made smaller or
discreticezed through data transformation methods such as discrete wavelet and
Fourier transforms (DWT and DFT respectively).

Similar to the Fourier transform, DWT represents a data (function) in terms of
an orthonormal basis. One of the simplest wavelet transform is the Haar
transform which can be described a orthogonal matrix: e.g., $\displaystyle
4\times 4$ Haar transform is

$H_{4}={\frac{1}{2}}{\begin{bmatrix}1&1&1&1\\\ 1&1&-1&-1\\\
{\sqrt{2}}&-{\sqrt{2}}&0&0\\\ 0&0&{\sqrt{2}}&-{\sqrt{2}}\end{bmatrix}}$ (15)

Similarly to Walsh-Hadamard transform, the unnormalized Haar transform matrix
can be obtained by the following recurrence:

$H_{{2N}}={\begin{bmatrix}H_{{N}}\otimes[1,1]\\\
I_{{N}}\otimes[1,-1]\end{bmatrix}}.$ (16)

The transformation on an input vector $\displaystyle\mathbf{y^{\prime}}$ is
described by $\displaystyle\mathbf{y^{\prime}}=H_{N}\mathbf{x}$. From this
product, the rows of the Haar matrix allow us to draw different properties of
the temporal data: For instance, for $\displaystyle H_{4}$, the first row
takes the average of the input vector, while the second row measures the low
frequencies. The rest of the rows measures different parts of the input
vector.

In the wavelet transform, $\displaystyle\mathbf{y}$ can be represented with
the basis vectors that form $\displaystyle H_{n}$:

$\mathbf{y^{\prime}}=\sum_{i=1}^{q}a_{i}\mathbf{w_{i}},$ (17)

where $\displaystyle\mathbf{w_{i}}$ is the $\displaystyle i$th basis vector
and $\displaystyle a_{i}$ is the normalized coefficient. The dimension
reduction is done by disregarding the lowest $\displaystyle a_{i}$s.

In quantum computing, since the Haar matrix is orthogonal and can be built
from the recurrence that is described by the Kronecker product, it can be
implemented as a quantum gate in a similar fashion to the Hadamard gate.
Therefore, given $\displaystyle\left|\mathbf{y}\right\rangle$, one can easily
obtain $\displaystyle\left|\mathbf{y^{\prime}}\right\rangle$. The elimination
of lowest parts of $\displaystyle\left|\mathbf{y^{\prime}}\right\rangle$ can
be done through either sampling or measuring a few qubits. This leads to a
problem of finding the best qubits that gives the maximum reduction with the
minimum information loss.

####  III-C1 Quantum Fourier transform (QFT)

The discrete Fourier transform (DFT) can be used to transform a data vector
into a linear combination of sinusoidal series. The similarity of two vectors
can be measured by taking the Euclidean distance of the Fourier coefficients.
In quantum computing, DFT-which has the complexity bound $\displaystyle
O(nlogn)$ for $\displaystyle n$-dimensional vector-can be implemented with
$\displaystyle O(logn)$ computational steps. It provides the main speedup of
many quantum algorithms such as integer factoring. The quantum Fourier
transform for the data $\displaystyle\left|\mathbf{y}\right\rangle$ can be
obtained by applying the operation $\displaystyle QFT$. If we have two data
vectors $\displaystyle\mathbf{y_{1}}$ and $\displaystyle\mathbf{y_{2}}$ of the
same dimension $\displaystyle n$ (we assume this for convenience and brevity
in notations, it can be easily generalized to any dimension.), we first
construct the following quantum state:

$\left|\mathbf{y}\right\rangle=\left(\begin{matrix}\mathbf{y_{1}}\\\
\mathbf{y_{2}}\end{matrix}\right)$ (18)

Then we apply the quantum Fourier transform to each part of the vectors:

$\left(\begin{matrix}QFT&\\\
&QFT\end{matrix}\right)\left|\mathbf{y}\right\rangle=\left(\begin{matrix}QFT\mathbf{y_{1}}\\\
QFT\mathbf{y_{2}}\end{matrix}\right)$ (19)

Here note that

$\left(\begin{matrix}QFT&\\\ &QFT\end{matrix}\right)=I\otimes QFT,$ (20)

where $\displaystyle I$ is a $\displaystyle 2\times 2$ identity matrix. We
then apply the Hadamard gate to the first qubit:

$\left(H\otimes I^{\otimes
log(n)}\right)\left(\begin{matrix}QFT\mathbf{y_{1}}\\\
QFT\mathbf{y_{2}}\end{matrix}\right)=\frac{1}{\sqrt{2}}\left(\begin{matrix}QFT\mathbf{y_{1}}+QFT\mathbf{y_{2}}\\\
QFT\mathbf{y_{1}}-QFT\mathbf{y_{2}}\end{matrix}\right)$ (21)

We can also start with a superposition of the input state:

$\left|\mathbf{y}\right\rangle=\frac{1}{\sqrt{2}}\left(\left|\mathbf{y_{1}}\right\rangle-\left|\mathbf{y_{2}}\right\rangle\right).$
(22)

Then we apply $\displaystyle QFT$:

$QFT\left|\mathbf{y}\right\rangle=\frac{1}{\sqrt{2}}\left(QFT\left|\mathbf{y_{1}}\right\rangle-
QFT\left|\mathbf{y_{2}}\right\rangle\right).$ (23)

This also generates a kind of the desired distance.

As mentioned before, we can also use a separate register for each data point
at different time: e.g., the quantum state
$\displaystyle\left|\mathbf{y_{1}}\right\rangle$$\displaystyle\left|\mathbf{y_{2}}\right\rangle$
is for the data measured at time t=1 and t=2. In this case we apply
$\displaystyle QFT$ to both registers and then use the SWAP test to create the
distance measure between data points.

$\left(QFT\otimes
QFT\right)\left|\mathbf{y_{1}}\right\rangle\left|\mathbf{y_{2}}\right\rangle=QFT\left|\mathbf{y_{2}}\right\rangle
QFT\left|\mathbf{y_{2}}\right\rangle$ (24)

In all cases, the Euclidean distance can be obtained by the measurement
results. However, if $\displaystyle d$ number of $\displaystyle\mathbf{y_{i}}$
are given as classical input, the first approach requires $\displaystyle
O(log(dn))$ number of qubits and the third approach requires $\displaystyle
O(dlog(n))$ since a separate register is used for each
$\displaystyle\mathbf{y_{i}}$. Since the second approach uses the
superpositioned state, the number of qubits is $\displaystyle O(\log(n))$
which is less than the first and third approaches. On the other hand, the
classical processing time is almost the same for all approaches and is
$\displaystyle O(nd)$ since it requires linear scanning of the every data
points. However, if $\displaystyle\left|\mathbf{y_{i}}\right\rangle$ are
stored on a random accessed quantum memory (qRAM), then all approaches
requires only $\displaystyle O(d)$ number of queries to qRAM to load the data
into quantum registers.

##  IV Time series forecasting

If the statistical parameters of time series such as the mean and variance do
not change with time, it is called stationary: i.e. the probabilistic
distribution of the parameters in any time interval is the same as or very
close to the time interval found by shifting from this interval. And non-
stationary if the distributions change with time, which is mostly the case
when we deal with real world data: e.g., if we look at the number of daily
COVID cases, we see that the daily cases changes based on the season and
different factors for the past few years. A nonstationary series can be made
stationary by different operations on the times series data.

###  IV-A Differencing

Time series can be converted by considering differences between the data
$\displaystyle y_{i}$ at time $\displaystyle t_{i}$ and the data
$\displaystyle y_{i-1}$ at the previous time $\displaystyle t_{i-1}$. This
operation is called differencing and defined as follows:

$y^{\prime}_{i}=y_{i}-y_{i-1}.$ (25)

If the data is still not stationary, then a second order differencing also can
be applied:

$y^{\prime\prime}_{i}=y^{\prime}_{i}+y^{\prime}_{i-1}=y_{i}-2y_{i-1}+y_{i-2}$
(26)

To implement this as a quantum operation, first we observe that the closer
$\displaystyle\left|y^{\prime}\right\rangle$ to the equal superposition state,
$\displaystyle\left|H\right\rangle=H^{\otimes n}\left|0\right\rangle$, the
more stationary it is. We can measure the closeness of the two quantum state
by the swap test[27] or by some statistical sampling from the quantum
state(i.e., measuring some qubits to obtain probabilities of 0 and 1. The same
or close values indicate an equal superposition state).

To construct $\displaystyle\left|y^{\prime}\right\rangle$ first we use an
auxiliary qubit in $\displaystyle\left|1\right\rangle$, then apply a Hadamard
gate to this qubit:

$\left(H\left|1\right\rangle\right)\left|y\right\rangle=\frac{1}{\sqrt{2}}\left(\begin{matrix}\left|y\right\rangle\\\
-\left|y\right\rangle\end{matrix}\right)$ (27)

Then we apply laddered Toffoli operations (similar to a bit-adder: See quantum
adders in Ref.[23]) to the second part of the state to convert into the
following:

$\frac{1}{\sqrt{2}}\left(\begin{matrix}\left|y\right\rangle\\\
-U_{adder}\left|y\right\rangle\end{matrix}\right)=\left(\begin{matrix}y_{0}\\\
y_{1}\\\ \vdots\\\ y_{n-1}\\\ -y_{n-1}\\\ -y_{0}\\\ \vdots\\\
-y_{n-2}\end{matrix}\right)$ (28)

Then, we apply a Hadamard gate again to the first qubit to obtain the
following

$\left(H\otimes I^{\otimes
n-1}\right)\frac{1}{\sqrt{2}}\left(\begin{matrix}y_{0}\\\ y_{1}\\\ \vdots\\\
y_{n-1}\\\ -y_{n-1}\\\ -y_{0}\\\ \vdots\\\
-y_{n-2}\end{matrix}\right)=\frac{1}{2}\left(\begin{matrix}y_{0}-y_{n-1}\\\
y_{1}-y_{0}\\\ \vdots\\\ y_{n-1}-y_{n-2}\\\ y_{0}+y_{n-1}\\\ y_{1}+y_{0}\\\
\vdots\\\ y_{n-1}+y_{n-2}\\\ \end{matrix}\right)$ (29)

Note that by applying Hadamard gates on different qubits and some swap
operations we can easily built
$\displaystyle\left|y^{\prime\prime}\right\rangle$ or any high order
differencing from the above state. For instance, one can obtain some seasonal
differencing as follows: For a natural number $\displaystyle s<i$, let the
seasonal difference be defined as

$y^{\prime}_{i}=y_{i}-y_{i-s}.$ (30)

In Eq.(29); instead of applying the Hadamard gate to the first qubit, if we
apply to the second or any other qubit based on the value of $\displaystyle
s$, we basically construct the equivalent quantum state for seasonal
difference. Also note that the desired quantum state is the half part of the
above state where the first qubit is in $\displaystyle\left|0\right\rangle$
state. In the further computations, one can either collapse the state to the
desired part by measuring the first qubit in
$\displaystyle\left|0\right\rangle$ state111For positive
$\displaystyle\mathbf{y}$, we can also apply an oblivious-ampltidue
amplification [28] to the first qubit or consider a different order of
differencing.

###  IV-B Forecasting Models

Differencing mean the time series is drifting in time; therefore, it can be
modeled by simply the following:

$y_{i+1}=y_{i}+c+e_{i+1},$ (31)

where $\displaystyle y_{i+1}$ is the predicted value, $\displaystyle e_{i+1}$
represents the white noise with zero mean, and $\displaystyle c$ is the time
drift.

Forecasting on time series data can be made through autoregressive models. The
forecasting with the univariate time series data is generally based on
autoregression models where the output value at time $\displaystyle t$,
$\displaystyle y_{t}$, is determined as a linear combination of the values
preceding window of some length $\displaystyle p$ (AR($\displaystyle p$)
model)[11]:

$y_{t}=\sum_{i=1}^{p}a_{i}\cdot y_{t-i}+c+\epsilon_{t},$ (32)

where $\displaystyle c$ and the coefficients $\displaystyle a_{1}\dots a_{p}$
are learned through a learning process.

In the moving average model (MA($\displaystyle q$)), the behavioral attribute
value at any timestamp is determined from the unexpected historical variations
in the time series data (shocks):

$y_{t}=\sum_{i=1}^{q}b_{i}\cdot\epsilon_{t-i}+c+\epsilon_{t},$ (33)

where $\displaystyle c$ is the mean value and the coefficients $\displaystyle
b_{i}$s are learned from the data. A more powerful and general model can be
obtained by combining the two aforementioned models (ARMA($\displaystyle
p,q$)):

$y_{t}=\sum_{i=1}^{p}a_{i}\cdot
y_{t-i}+c+\epsilon_{t}+\sum_{i=1}^{q}b_{i}\cdot\epsilon_{t-i}+c+\epsilon_{t}.$
(34)

Autoregressive moving average models (ARMA) are used best with the stationary
data (when the mean, variance, and autocorrelation does not change in time).
The non-stationary data are handled by integrating a series of differences
into the model (ARIMA($\displaystyle p,q,q$)). For instance with the first-
order difference ($\displaystyle d=1$) the model can be written as:

$y^{\prime}_{t}=\sum_{i=1}^{p}a_{i}\cdot
y^{\prime}_{t-i}+c+\epsilon_{t}+\sum_{i=1}^{q}b_{i}\cdot\epsilon_{t-i}+c+\epsilon_{t}.$
(35)

We can easily generate the above by using a variational quantum circuit
including two and single qubit gates (as done in many of quantum machine
learning models) and apply to the quantum state prepared in one of the
following forms:

$\left|\mathbf{y}\right\rangle\left|\mathbf{\epsilon}\right\rangle,\left(\theta_{y}\left|\mathbf{y}\right\rangle+\theta_{\epsilon}\left|\mathbf{\epsilon}\right\rangle\right),\text{
or }\left(\begin{matrix}\left|\mathbf{y}\right\rangle\\\
\left|\mathbf{\epsilon}\right\rangle\end{matrix}\right).$ (36)

Note that $\displaystyle\left|\mathbf{y}\right\rangle$ may also represent the
state after differencing. Also note that the angle values of the gates in the
circuit are related to the parameters $\displaystyle\mathbf{a}$ and
$\displaystyle\mathbf{b}$ so that we can obtain the parameterized model in
Eq.(35) or any of the previous models. Furthermore, note that following the
Ref.[29] where it is shown how to map the Schmidt decomposition of a vector
into a quantum circuit, one can also prepare a quantum operator whose first
row is the parameters $\displaystyle\mathbf{a}$ and $\displaystyle\mathbf{b}$.

$\displaystyle\displaystyle\left(\begin{matrix}a_{1}&\dots&a_{p}&b_{1}\dots&b_{q}\\\
\bullet&\dots&\bullet&\dots&\bullet\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
\end{matrix}\right)$ (37)

If the dimensions of $\displaystyle\mathbf{a}$ and $\displaystyle\mathbf{b}$
are $\displaystyle k$, this operator can be generated by using $\displaystyle
O(k)$ number of quantum operations. We map the parameters to the angles of the
rotation gates. Then, we adjust the rotation gate parameters by using a
classical optimization technique as done in the variational solvers.

Since the value of $\displaystyle p$ determines the seasonality (period) in
the time series. Therefore, ARIMA model indicates the auto-correlation between
a value and its neighbors. In some ARIMA like models, $\displaystyle p$ can be
considered as the period of the time series. On quantum computers the value of
$\displaystyle p$ can be found along with the other optimization parameters by
using a classical optimization. In addition, since it is the period of the
data, one can also adapt quantum phase estimation algorithm to find this
period as done in Shor’s factoring algorithm.

###  IV-C Forecasting Models As Circulant Matrix Operators

If the equation for $\displaystyle y_{t}$ in any of the models is stacked up
for the values $\displaystyle t=0,1,\dots,T-1$, the above generic models can
be rewritten as a system of equations constructed by circulant matrices [30].
The circulant matrices can be implemented very efficiently on quantum
computers since their eigenspace is formed by the Fourier transform (see
Ref.[31] for circuit implementation of these matrices). Therefore, one can use
those matrices in conjunction with the variational quantum circuits to
optimize the parameters of the forecasting models.

##  V Discussion and Future Directions

In this paper, we go through data processing and forecasting with ARIMA models
and show how to implement them on quantum computers. As future direction, more
example data analysis are needed to investigate algorithms and methods for
preprocessing and transformation of time series data into the Hilbert space in
which quantum computers work better than classical computers; and investigate
quantum machine learning methods and algorithms for univariate and/or
multivariate time series data to forecast, classify and cluster the data. In
addition, for time series data, there are many classical packages such as
Facebook’s- Prophet [32] and ARIMA (Autoregressive Integrated Moving Average)
models. There is also need to design a quantum packages-that can be used with
the current quantum libraries such as IBM-QISKIT[33]-for time series data to
forecast and do other tasks.

###  V-A Data representation and Measurement

Representation of data on quantum computers may not be without noise.
Therefore, the small changes in any feature of data may not be observed on
quantum state by using simple measurements. If one qubit is used for one
feature, then we expect any change in the input data to impact the measurement
statistics of a qubit. However, representing whole feature vector as a quantum
state with fewer qubits enforces an automatic dimension reduction in the
measurement and hence may cause a significant data loss that impedes the
prediction results.

###  V-B Quantum optimization for autoregressive models

Note that autoregressive models can be also combined with quantum neural
networks as done with the classical neural networks [34]. Also note that
autoregressive models can be formulated as combinatorial optimization problems
in the framework of graph learning [35]. In addition, it is shown that
quadratic unconstrained optimization formulation can be used for forecasting
in finance [36]. With the help of these and other similar studies,
autoregressive models can be formulated as combinatorial optimization in the
form of quadratic unconstrained binary optimization. Therefore, the quantum
optimization approaches such as adiabatic quantum computation [37], quantum
unconstrained bounded optimization algorithm [38], quantum power iteration
[39] can be used for these models.

## References

  * [1] R. P. Feynman, “Quantum mechanical computers,” _Foundations of physics_ , vol. 16, no. 6, pp. 507–531, 1986.
  * [2] P. Benioff, “The computer as a physical system: A microscopic quantum mechanical hamiltonian model of computers as represented by turing machines,” _Journal of statistical physics_ , vol. 22, no. 5, pp. 563–591, 1980.
  * [3] P. W. Shor, “Algorithms for quantum computation: discrete logarithms and factoring,” in _Proceedings 35th annual symposium on foundations of computer science_. Ieee, 1994, pp. 124–134.
  * [4] L. K. Grover, “Quantum mechanics helps in searching for a needle in a haystack,” _Physical review letters_ , vol. 79, no. 2, p. 325, 1997.
  * [5] S. Lloyd, M. Mohseni, and P. Rebentrost, “Quantum principal component analysis,” _Nature Physics_ , vol. 10, no. 9, pp. 631–633, 2014.
  * [6] A. N. Kolmogorov, “On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition,” in _Doklady Akademii Nauk_ , vol. 114, no. 5. Russian Academy of Sciences, 1957, pp. 953–956.
  * [7] R. Hecht-Nielsen, “Kolmogorov’s mapping neural network existence theorem,” in _Proceedings of the international conference on Neural Networks_ , vol. 3. IEEE Press New York, NY, USA, 1987, pp. 11–14.
  * [8] G. Cybenko, “Approximation by superpositions of a sigmoidal function,” _Mathematics of control, signals and systems_ , vol. 2, no. 4, pp. 303–314, 1989.
  * [9] A. R. Gallant and H. White, “There exists a neural network that does not make avoidable mistakes.” in _ICNN_ , 1988, pp. 657–664.
  * [10] Z.-Q. John Xu, Y. Zhang, T. Luo, Y. Xiao, and Z. Ma, “Frequency principle: Fourier analysis sheds light on deep neural networks,” _Communications in Computational Physics_ , vol. 28, no. 5, pp. 1746–1767, 2020. [Online]. Available: http://global-sci.org/intro/article_detail/cicp/18395.html
  * [11] C. C. Aggarwal _et al._ , _Data mining: the textbook_. Springer, 2015, vol. 1.
  * [12] H. Ismail Fawaz, G. Forestier, J. Weber, L. Idoumghar, and P.-A. Muller, “Deep learning for time series classification: a review,” _Data mining and knowledge discovery_ , vol. 33, no. 4, pp. 917–963, 2019.
  * [13] A. Abanda, U. Mori, and J. A. Lozano, “A review on distance based time series classification,” _Data Mining and Knowledge Discovery_ , vol. 33, no. 2, pp. 378–412, 2019.
  * [14] M. S. Gashler and S. C. Ashmore, “Modeling time series data with deep fourier neural networks,” _Neurocomputing_ , vol. 188, pp. 3–11, 2016.
  * [15] V. Dunjko and P. Wittek, “A non-review of quantum machine learning: trends and explorations,” _Quantum Views_ , vol. 4, p. 32, 2020.
  * [16] V. Dunjko and H. J. Briegel, “Machine learning & artificial intelligence in the quantum domain: a review of recent progress,” _Reports on Progress in Physics_ , vol. 81, no. 7, p. 074001, 2018.
  * [17] P. Singh, G. Dhiman, and A. Kaur, “A quantum approach for time series data based on graph and schrödinger equations methods,” _Modern Physics Letters A_ , vol. 33, no. 35, p. 1850208, 2018.
  * [18] D. Emmanoulopoulos and S. Dimoska, “Quantum machine learning in finance: Time series forecasting,” _arXiv preprint arXiv:2202.00599_ , 2022.
  * [19] M. Schuld, R. Sweke, and J. J. Meyer, “Effect of data encoding on the expressive power of variational quantum-machine-learning models,” _Physical Review A_ , vol. 103, no. 3, p. 032430, 2021.
  * [20] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Quesada, and S. Lloyd, “Continuous-variable quantum neural networks,” _Physical Review Research_ , vol. 1, no. 3, p. 033063, 2019.
  * [21] A. Daskin, “A simple quantum neural net with a periodic activation function,” in _2018 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_. IEEE, 2018, pp. 2887–2891.
  * [22] V. Giovannetti, S. Lloyd, and L. Maccone, “Quantum random access memory,” _Physical review letters_ , vol. 100, no. 16, p. 160501, 2008.
  * [23] M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” 2002\.
  * [24] A. Pérez-Salinas, A. Cervera-Lierta, E. Gil-Fuster, and J. I. Latorre, “Data re-uploading for a universal quantum classifier,” _Quantum_ , vol. 4, p. 226, 2020.
  * [25] M. Schuld and N. Killoran, “Quantum machine learning in feature hilbert spaces,” _Physical review letters_ , vol. 122, no. 4, p. 040504, 2019.
  * [26] E. Grant, M. Benedetti, S. Cao, A. Hallam, J. Lockhart, V. Stojevic, A. G. Green, and S. Severini, “Hierarchical quantum classifiers,” _npj Quantum Information_ , vol. 4, no. 1, pp. 1–8, 2018.
  * [27] H. Buhrman, R. Cleve, J. Watrous, and R. De Wolf, “Quantum fingerprinting,” _Physical Review Letters_ , vol. 87, no. 16, p. 167902, 2001.
  * [28] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, “Exponential improvement in precision for simulating sparse hamiltonians,” in _Forum of Mathematics, Sigma_ , vol. 5. Cambridge University Press, 2017.
  * [29] A. Daskin, A. Grama, G. Kollias, and S. Kais, “Universal programmable quantum circuit schemes to emulate an operator,” _The Journal of chemical physics_ , vol. 137, no. 23, p. 234112, 2012.
  * [30] D. S. G. Pollock, “Circulant matrices and time-series analysis,” _International Journal of Mathematical Education in Science and Technology_ , vol. 33, no. 2, pp. 213–230, 2002.
  * [31] A. Daskin, “Quantum implementation of circulant matrices and its use in quantum string processing,” _arXiv preprint arXiv:2206.09364_ , 2022.
  * [32] “Prophet,” accessed: Feb. 15, 2022. [Online]. Available: https://facebook.github.io/prophet/
  * [33] “Qiskit,” accessed: Feb. 15, 2022. [Online]. Available: https://qiskit.org/
  * [34] G. P. Zhang, “Time series forecasting using a hybrid arima and neural network model,” _Neurocomputing_ , vol. 50, pp. 159–175, 2003.
  * [35] Y. Peng, B. Choi, and J. Xu, “Graph learning for combinatorial optimization: A survey of state-of-the-art,” _Data Science and Engineering_ , vol. 6, no. 2, pp. 119–141, 2021.
  * [36] R. Orús, S. Mugel, and E. Lizaso, “Forecasting financial crashes with quantum computing,” _Physical Review A_ , vol. 99, no. 6, p. 060301, 2019\.
  * [37] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, “Quantum computation by adiabatic evolution,” _arXiv preprint quant-ph/0001106_ , 2000.
  * [38] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” _arXiv preprint arXiv:1411.4028_ , 2014.
  * [39] A. Daskin, “Combinatorial optimization through variational quantum power method,” _Quantum Information Processing_ , vol. 20, no. 10, pp. 1–14, 2021\.