matchrangestartf#1

Satellite Quantum Modelling & Analysis Software
Version 1.1: Documentation

J. S. Sidhu, T. Brougham, D. McArthur, R. G. Pousa and D. K. L. Oi

University of Strathclyde
Department of Physics
John Anderson Building, 107 Rottenrow East, Glasgow, G4 0NG
December 2021

##  Chapter 1 Introduction

###  1.1 Mission statement

We provide a numerical key rate analysis, which determines the amount of
expected key generation in satellite-based quantum key distribution protocols.
This key length analysis will help develop an intuition on the effects of
different operational scenarios on the key rate and inform the development of
source and receiver systems. This numerical toolkit will provide a guide to
future satellite missions.

###  1.2 Scope of current version

The current release of SatQuMA, version 1.1, is an update to the first working
version that could calculate finite key lengths for a limited set of systems
and circumstances. The new version, primarily, provides greatly improved
optimisation routines and introduces asymptotic expressions. We implement an
optimised, asymmetric two-decoy state BB84 protocol with weak coherent pulses.
Elevation, and time, dependent system losses are defined externally and read-
in from a data file; an example file is supplied with the software. We
consider multiple satellite overpasses, but only for identical orbits.
Protocol parameters can be either optimised or specified, collectively.

###  1.3 List of updates

The latest release contains many small improvements over the original v1.0
release, listed here in no particular order:

  * •

The system parameters Pec ($p_{ec}$), the probability of extraneous counts,
and QBERI ($\text{QBER}_{\text{I}}$), the intrinsic quantum bit error rate
(previously called Pdc and PolError respectively) are now iterable by default.

  * •

Due to the introduction of Pec and QBERI as iterable parameters the naming
convention for output files has changed such that each file is appended with
the relevant values of these parameters. Users are now asked to provide the
base name (or prefix) for all output files, if requested.

  * •

We have introduced a parameter, shift_elev0, which allows the user to shift
the centre of the transmission window from zenith by a specified number of
degrees.

  * •

We have removed the boolean flag tChernoff, which switched between the
Chernoff and Hoeffding tail bounds, and introduced a list of bounds that can
be selected from - including ‘Chernoff’, ‘Hoeffding’, and ‘Asymptotic’ (i.e.
none). Selecting the ‘Asymptotic’ option for the tail bounds also results in
some other default parameters being automatically selected, regardless of the
users choices, so as to ensure that the asymptotic key rate is calculated
properly.

  * •

We have incorporated the other two available constrained optimisation methods
from the SciPy package, namely ‘COBYLA’ (Constrained Optimisation BY Linear
Approximation) and ‘SLSQP’ (Sequential Least SQuares Programming). The former
provides a much more robust optimisation than the previous algorithm (‘trust-
constr’) and has been set as the default.

  * •

Optimisations for a particular set of parameters are now performed in a while
loop, primarily allowing for a brute force approach to finding non-zero keys
in high loss/QBER systems. This can also greatly improve the performance of
the optimiser but at the cost of additional computations. The user can now
specify the minimum number of optimisations to be performed via the parameter
NoptMin. Additionally, there are two boolean flags which can be used to try
and minimise the number of optimisations performed. First, the flag tStopZero
will break out of the while loop if the first NoptMin optimisations all return
a value of zero – this stops the software from wasting time looking for key in
parameter regimes where there is none. Second, the flag tStopBetter simply
breaks out of the while loop if a larger amount of key has been found and the
minimum number of optimisations has been satisfied.

  * •

The way in which the initial optimisation parameters are set has also been
changed. If a user has requested that these parameters are initialised from a
set of specified values, then these values are used for each calculation –
preciously they were only used for the first calculation. If a user has
instead requested that the parameters are to be initialised randomly then, if
it is not the first calculation, and the previous calculation did not produce
zero key, the parameters from the previous calculation are tried first before
being randomised. If it is the first calculation, or the previous calculation
did return zero key, then a new set of random parameters are generated.

  * •

We have removed the boolean flag tSortData, which enabled sorting of the data
to be written out, as it is not necessary without parallel computations. We
have also introduced the boolean flag tMetrics which allows the user to
explicitly specify whether the optimisation metrics for each (final)
calculation should be written to file.

  * •

We have introduced two calculation timers; one for the total calculation
runtime and one for each block of losses and transmission time windows.

  * •

We have corrected a minor bug in the method of determining the size of the
data arrays. This bug was only an issue when using an even-valued step size
for either the losses (ls) or transmission time half-windows (dt) arrays with
an odd-valued end point.

###  1.4 Installation

The latest version of the SatQuMA software can be found at
https://github.com/cnqo-qcomms/SatQuMA.

####  1.4.1 Required packages

The current version of SatQuMA requires an installation of Python 3.* and the
following standard packages:

  * •

scipy (SciPy)

  * •

numpy (NumPy)

  * •

sys

  * •

time

##  Chapter 2 Theoretical summary

###  2.1 Background

Here we present a high-level summary of the equations required to calculate
the secret key length (SKL) which appear in the current SatQuMA release. We do
not, however, provide any form of derivation for these relations here.

####  2.1.1 Protocol and statistics

In our protocol, Alice randomly prepares a state in the basis $X$ or $Z$,
where we usually assume $X=\\{D,A\\}$ and $Z=\\{H,V\\}$, with one of three
intensities $\mu=\\{\mu_{1},\mu_{2},\mu_{3}\\}$, each with a probability of
being selected $P_{\mu}=\\{P_{\mu_{1}},P_{\mu_{2}},P_{\mu_{3}}\\}$ .

Once Alice has sent her signals to Bob and the reconciliation process, error
correction, and post-processing has been completed we can define some
measurement statistics from the sifted key. We define the number of events,
for each basis, for each intensity Alice could prepare

$\displaystyle n_{\text{X},\mu}$
$\displaystyle:=\\{n_{\text{X},\mu_{1}},n_{\text{X},\mu_{2}},n_{\text{X},\mu_{3}}\\},$
(2.1) $\displaystyle n_{\text{Z},\mu}$
$\displaystyle:=\\{n_{\text{Z},\mu_{1}},n_{\text{Z},\mu_{2}},n_{\text{Z},\mu_{3}}\\},$
(2.2)

and similarly we define the number of bit errors, for each basis, for each
intensity

$\displaystyle m_{\text{X},\mu}$
$\displaystyle:=\\{m_{\text{X},\mu_{1}},m_{\text{X},\mu_{2}},m_{\text{X},\mu_{3}}\\},$
(2.3) $\displaystyle m_{\text{Z},\mu}$
$\displaystyle:=\\{m_{\text{Z},\mu_{1}},m_{\text{Z},\mu_{2}},m_{\text{Z},\mu_{3}}\\}.$
(2.4)

####  2.1.2 Secure key length

The length of the secure key is given by [1]

$\ell=\left\lfloor
s_{\text{X},0}+s_{\text{X},1}\left[1-h(\phi_{\text{X}})\right]-\lambda_{\text{EC}}-6\log_{2}\left(\frac{21}{\epsilon_{\text{s}}}\right)-\log_{2}\left(\frac{2}{\epsilon_{\text{c}}}\right)\right\rfloor,$
(2.5)

where the outer brackets indicate that we should take the floor of this
expression. Here, $s_{\text{X},0}$ is the number of vacuum events,
$s_{\text{X},1}$ is the number of single-photon events, and $\phi_{\text{X}}$
is the phase error rate in the sifted $X$ basis. The parameter
$\lambda_{\text{EC}}$ provides an estimate, or a bound, on the number of bits
required for error correction although this should be replaced with the
_actual_ number of bits used when this is known. The security parameters
$\epsilon_{\text{c}}$ and $\epsilon_{\text{s}}$ are the prescribed security
parameters which define the _correctness_ and _secrecy_ of the resulting key
respectively. The binary entropy function used above is defined as

$h\left(x\right)=-x\log_{2}x-\left(1-x\right)\log_{2}\left(1-x\right).$ (2.6)

####  2.1.3 Number of vacuum events

The number of vacuum events in a particular basis is evaluated as, for
example,

$s_{\text{X},0}\geq\tau_{0}\frac{\mu_{2}n_{\text{X},\mu_{3}}^{-}-\mu_{3}n_{\text{X},\mu_{2}}^{+}}{\mu_{2}-\mu_{3}},$
(2.7)

where the probability that Alice sends an $n$-photon state is given by the
Poisson distribution,

$\tau_{n}=\sum\limits_{j=1}^{3}\frac{\mathrm{e}^{-\mu_{j}}\mu_{j}^{n}p_{j}}{n!}.$
(2.8)

In order to account for statistical fluctuations in the expected number of
$n$-photon events, we apply the Chernoff bound and define the functions [2]

$\displaystyle n_{\text{X},\mu_{j}}^{+}$
$\displaystyle=\frac{\mathrm{e}^{\mu_{j}}}{p_{\mu_{j}}}\left[n_{\text{X},\mu_{j}}+\log_{\mathrm{e}}\left(\frac{21}{\epsilon_{\text{s}}}\right)+\sqrt{2n_{\text{X},\mu_{j}}\log_{\mathrm{e}}\left(\frac{21}{\epsilon_{\text{s}}}\right)+\log_{\mathrm{e}}\left(\frac{21}{\epsilon_{\text{s}}}\right)^{2}}\right],$
(2.9a) $\displaystyle n_{\text{X},\mu_{j}}^{-}$
$\displaystyle=\frac{\mathrm{e}^{\mu_{j}}}{p_{\mu_{j}}}\left[n_{\text{X},\mu_{j}}-\frac{1}{2}\log_{\mathrm{e}}\left(\frac{21}{\epsilon_{\text{s}}}\right)-\sqrt{2n_{\text{X},\mu_{j}}\log_{\mathrm{e}}\left(\frac{21}{\epsilon_{\text{s}}}\right)+\frac{1}{4}\log_{\mathrm{e}}\left(\frac{21}{\epsilon_{\text{s}}}\right)^{2}}\right].$
(2.9b)

Note, we could also use these expressions to determine the number of vacuum
events in the $Z$ basis by exchanging the $n_{\text{X},\mu_{j}}$ terms for the
corresponding $n_{\text{Z},\mu_{j}}$ terms.

####  2.1.4 Number of single-photon events

The number of single-photon events in a particular basis is similarly
evaluated as,

$s_{\text{X},1}\geq\tau_{1}\frac{\mu_{1}\left[n_{\text{X},\mu_{2}}^{-}-n_{\text{X},\mu_{3}}^{+}-\frac{\mu_{2}^{2}-\mu_{3}^{2}}{\mu_{1}^{2}}\left(n_{\text{X},\mu_{1}}^{+}-\frac{s_{\text{X},0}}{\tau_{0}}\right)\right]}{\mu_{1}\left(\mu_{2}-\mu_{3}\right)-\mu_{2}^{2}+\mu_{3}^{2}},$
(2.10)

where $\tau_{1}$ is given by (2.8). As with the expressions for the number of
vacuum events, we can use the same expression to determine the number of
single-photon events in the $Z$ basis by again exchanging the
$n_{\text{X},\mu_{j}}$ terms for the corresponding $n_{\text{Z},\mu_{j}}$
terms.

####  2.1.5 The phase error rate

We evaluate the phase error rate in the $X$ basis according to

$\phi_{\text{X}}\leq\frac{v_{\text{Z},1}}{s_{\text{Z},1}}+\gamma\left(\epsilon_{\text{s}},\frac{v_{\text{Z},1}}{s_{\text{Z},1}},s_{\text{Z},1},s_{\text{X},1}\right),$
(2.11)

where we use the function

$\gamma\left(a,b,c,d\right)=\sqrt{\frac{\left(c+d\right)\left(1-b\right)b}{cd\log_{\mathrm{e}}2}\log_{2}\left[\frac{c+d}{bcd\left(1-b\right)}\frac{21^{2}}{a^{2}}\right]},$
(2.12)

and the single-photon events are defined as above in Sec. 2.1.4. We have also
introduced the number of bit errors associated with single-photon events in
$Z$

$v_{\text{Z},1}\leq\tau_{1}\frac{m_{\text{Z},\mu_{2}}^{+}-m_{\text{Z},\mu_{3}}^{-}}{\mu_{2}-\mu_{3}},$
(2.13)

where the bounds on the number of bit errors due to statistical fluctuations,
$m_{\text{Z},\mu_{j}}^{\pm}$ are given by (2.9a) and (2.9b) where the
number(s) of events in the $X$ basis, $n_{\text{X},\mu_{j}}$, should be
substituted with the number(s) of bit errors in the $Z$ basis,
$m_{\text{Z},\mu_{j}}$.

####  2.1.6 Estimating the amount of error correction

We estimate the number of bits that need to be sacrificed to perform the error
correction in two ways: the first is more accurate but also much more complex;
the second is simple to implement but provides only a lower bound.

##### Method 1

We can evaluate the number of bits that we need to sacrifice for error
correction as [3]

$\displaystyle\lambda_{\text{EC}}$ $\displaystyle\approx
n_{X}h\left(\text{QBER}_{\text{X}}\right)+\left[n_{X}\left(1-\text{QBER}_{\text{X}}\right)-F^{-1}\left(\epsilon_{\text{c}};\left\lfloor
n_{\text{X}}\right\rfloor,1-\text{QBER}_{\text{X}}\right)-1\right]$
$\displaystyle\times\log_{\mathrm{e}}\left[\frac{\left(1-\text{QBER}_{\text{X}}\right)}{\text{QBER}_{\text{X}}}\right]-\frac{1}{2}\log_{\mathrm{e}}n_{\text{X}}-\log_{\mathrm{e}}\left(\frac{1}{\epsilon_{\text{c}}}\right),$
(2.14)

where we define the quantum bit error rate in the $X$ basis as

$\text{QBER}_{\text{X}}=\frac{\sum_{j}m_{\text{X},\mu_{j}}}{\sum_{j}n_{\text{X},\mu_{j}}},\;\text{for}\;j\in\\{1,2,3\\},$
(2.15)

and $F^{-1}(k;n,p)$ is the inverse (or quantile function) of the binomial
cumulative distribution function

$F(k;n,p)=\sum\limits_{i=0}^{\left\lfloor
k\right\rfloor}\left(\begin{matrix}n\\\
i\end{matrix}\right)p^{i}\left(1-p\right)^{n-i}.$

##### Method 2

Another method, based upon the block size, estimates the lower bound on the
error correction as

$\lambda_{\text{EC}}\geq
1.16\sum\limits_{j=1}^{3}n_{\text{X},\mu_{j}}h\left(\text{QBER}_{\text{X}}\right).$
(2.16)

##### Method 3

We can simply estimate the lower bound on the error correction based upon the
total number of bit errors in the $X$ basis

$\lambda_{EC}\geq 1.16\sum\limits_{j=1}^{3}m_{X,\mu_{j}}.$ (2.17)

##  Chapter 3 Satellite overpass geometry

In this chapter we will discuss the way in which we define a satellite
overpass, and how SatQuMA expects a transmission window to be specified for
the purposes of a key length calculation.

In SatQuMA, the transmission window is assumed to be symmetric in time (and
elevation) about the local zenith by default and discretized into time-slots.
An illustration of an ideal, zenith satellite overpass is shown in Fig. 3.1.
To define the overpass geometry we set a minimum elevation angle,
$\theta_{\text{min}}$, for which signals can be transmitted between the
satellite and OGS/receiver and choose a transmission window half-width,
$\Delta t$. The software will attempt to generate secret key over each time-
slot while $-\Delta t\leq t\leq\Delta t$.

We can also consider non-ideal satellite overpasses, where the orbital
geometry can be defined in terms of either the maximum elevation of that
orbit, relative to the local horizon of the receiver, or the orbit rotation
angle $\xi$, relative to the centre of the Earth, required to transform an
ideal zenith overpass into the non-ideal overpass specified. The relation
between these two angles is illustrated in Fig. 3.2. The maximum elevation is
related to the zenith orbit rotation angle $\xi$ as

$\theta_{\text{max}}=\cos^{-1}\left\\{\frac{\left(R_{\text{E}}+h_{\text{sat}}\right)\sin\xi}{\left[\begin{array}[]{c}\left(R_{\text{E}}+h_{\text{sat}}\right)^{2}+\left(R_{\text{E}}+h_{\text{OGS}}\right)^{2}-\\\
2\left(R_{\text{E}}+h_{\text{sat}}\right)\left(R_{\text{E}}+h_{\text{OGS}}\right)\cos\xi\end{array}\right]}\right\\},$
(3.1)

where $R_{\text{E}}$ is the radius of the Earth, $h_{\text{sat}}$ is the
orbital altitude of the satellite, and $h_{\text{OGS}}$ is the altitude of the
receiver/ground station.

Figure 3.1: Satellite overpass geometry. A satellite passes over an OGS, where
the satellite elevation angle $\theta\in\left[0^{\circ},90^{\circ}\right]$ is
measured from the local horizon of the OGS. The satellite and OGS can only
close a link when the satellite is above the minimum elevation angle
$\theta_{\text{min}}$. We assume that the actual transmission window is
symmetric about the local zenith, which we label as $t=0$, and that signals
are sent while $-\Delta t\leq t\leq\Delta t$.
OGS photo: ESA
Globe: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO, Landsat/Copernicus IBCAO
U.S. Geological Survey  Figure 3.2: Non-zenith satellite orbits. The satellite
orbit can be defined in terms of the maximum elevation that the satellite will
reach, $\theta_{\text{max}}\in\left[0^{\circ},90^{\circ}\right]$, above the
local horizon and also the angle with respect to the centre of the Earth,
$\xi$, that an ideal zenith orbit needs to be rotated to have equivalent
geometry.
OGS photo: ESA
Globe: Google, Data SIO, NOAA, U.S. Navy, NGA, GEBCO, Landsat/Copernicus IBCAO
U.S. Geological Survey

##  Chapter 4 Example of use

Here we go through the process of setting up an optimisation calculation using
SatQuMA, with script excerpts taken from SatQuMA_1.1.py as indicated by the
specified line numbers.

Many of the following parameter flags are set by switching between indices of
a simple boolean array.

⬇

54F_or_T = [False, True] # List used to switch between False or True values

###  4.1 Optimisation parameters

SatQuMA allows the main protocol parameters to be either optimised or
specified, collectively. To enable the optimisation of the protocol parameters
we must set the relevant boolean flag to True.

⬇

72tOptimise = F_or_T[1] # False (0) or True (1)

The optimiser requires these parameters to be given a range, with an upper and
lower bound, which we define using a numpy array.

⬇

83xb = np.array([[0.3,1.0],[0.6,0.9999],[0.0,0.4],[0.3,1.0],[0.1,0.5]])

Each pair of numbers (columns) of the numpy array are the (non-inclusive)
lower and upper bounds for the parameters in order
$P_{x},P_{\mu_{1}},P_{\mu_{2}},\mu_{1},\mu_{2}$.

We must now set initial values for the optimised parameters, keeping within
the pre-defined, respective parameter ranges. These initial values can either
be directly specified or they can be selected for us randomly. We note that
initialising the parameters can be an important step as the returned secret
key length (SKL) can be zero over relatively broad parameter regions, or
indeed in smaller localised regions, which may cause the optimiser to report
spurious issues with the supplied function (that it’s derivatives appear to be
zero). In order to specify the initial values, we first set the relevant
boolean flag to be True.

⬇

93tInit = F_or_T[1] # False (0) or True (1)

Next, we set the values for each parameter.

⬇

99Px_i = 0.5 # Asymmetric polarisation probability

100pk1_i = 0.7 # Probability Alice prepares intensity 1

101pk2_i = 0.1 # Probability Alice prepares intensity 2

102mu1_i = 0.8 # Intensity 1

103mu2_i = 0.3 # Intensity 2

If we instead want the initial parameters to be set randomly we set the
following boolean flag to False.

⬇

93tInit = F_or_T[0] # False (0) or True (1)

###  4.2 Calculation parameters

####  4.2.1 Input file options

SatQuMA requires a time/elevation vs link efficiency data file (in CSV format)
to be specified, hereafter referred to as the ‘loss file’. An example file has
been supplied with the software, the name of which we must specify.

⬇

132loss_file = 'FS_loss_XI0.csv'

We can also specify the directory containing the loss file, if it is not
located in the current working directory. If the loss file is in the current
work directory, then we can set this path to be an empty string.

⬇

130loss_path = ''

Additionally, we can specify the column within the loss file that contains the
link efficiency (losses), however the default is column 3.

⬇

133lc = 3 # Column containing loss data in file (counting from 1)

####  4.2.2 System parameters

We now specify the relevant parameters which characterise the performance of
the system we wish to model. First, we specify the orbit offset angle $\xi$
(rad.) which define the smallest angle between the orbit plane of the
satellite and the zenith plane of the optical ground station (OGS).

⬇

139xi = 0.0 # Angle between OGS zenith and satellite (from Earth's centre)
[rad.]

This value should be the same as that used when producing the input loss file,
as such we typically include the value of $\xi$ in the name of this file. At
present, however, this value is simply included in the output data and doesn’t
factor into the calculations directly.

Next, we specify the remaining protocol parameters: the intensity of the third
weak coherent pulse (second decoy state) $\mu_{3}$ and the prescribed errors
in both correctness and secrecy, $\epsilon_{\text{c}}$ and
$\epsilon_{\text{s}}$ respectively.

⬇

141mu3 = 0 # Intensity of pulse 3 (fixed)

⬇

143eps_c = 10**(-15) # Correctness parameter

144eps_s = 10**(-9) # Secrecy parameter

Here, by setting $\mu_{3}=0$ we have chosen the third pulse in our protocol to
be the vacuum state.

The intrinsic Quantum Bit Error Rate ($\text{QBER}_{\text{I}}$) is an iterable
parameter, and the values to calculate should be passed as an iterable object:
list, array, tuple/singleton, or generator

⬇

146QBERI_list = [0.001,0.003,0.005] # list, array, tuple or singleton

The extraneous count probability ($P_{\text{ec}}$) is also iterable.

⬇

148Pec_list = [1e-8,1e-7,1e-6] # list, array, tuple or singleton

The after-pulse probability of each detector must also be specified.

⬇

150Pap = 0.001 # After-pulse probability

Finally, we specify the number of (identical) satellite overpasses to include
and the repetition rate of the transmission source (Hz).

⬇

152NoPass = 1 # Number of satellite passes

⬇

154Rrate = 1*10**(9) # Source rate (Hz)

####  4.2.3 Time window and system loss

The first version of SatQuMA has been designed to calculate (loop) over the
duration of the overpass time half-window (s) and the relative system loss
(dB). The time window loop is principally controlled by a numpy array
specifying the start, stop and step indices relating to the time slots (as
specified in the input loss file, see 4.2.1).

⬇

162dt_range = np.array([200, 350, 10]) # Start, stop, step

Note, these time slots are labelled relative to the overpass zenith for which
we arbitrarily set $t=0$. Here we have requested SKL calculations with
overpass transmit time half-windows of, initially, up to 200 s then rising to
half-windows of a duration of 350 s in 10 s increments. We can also set a
minimum elevation for transmission (in degrees) which will override the values
we have just specified where necessary.

⬇

165min_elev = 10.0 # Minimum elevation transmission angle (degs)

Here, we have specified that no transmission is possible for elevations below
$10^{\circ}$. We can also shift the elevation angle taken as the centre of the
pass ($t=0$) in degrees.

⬇

166shift_elev = 0.0 # Shift the elevation angle taken as t = 0 (degs)

Next, we define the start, stop and step values for a loop over the excess
system losses. That is, we can add additional loss to those specified in the
loss file (converting from system efficiency).

⬇

169ls_range = np.array([0, 12, 2]) # Start, stop, step value

Here, we have considered systems which have 0 to 12 dB of excess loss above
the losses specified in the input loss file.

####  4.2.4 Output file options

During a calculation, SatQuMA can write to two different output streams: a
local file and the standard output (IDE/terminal/screen/etc). We need to
specify each output file that we wish to generate. To request that the full
calculation data for a given set of system loss metrics and time windows is
written to file (in CSV format) we first set the relevant boolean flag to be
True.

⬇

186tFullData = F_or_T[1] # False (0) or True (1)

If we want only the optimal time window data, we can set the following flag.

⬇

188tOptiData = F_or_T[1] # False (0) or True (1)

Likewise, if we want the optimal time window data for each calculation in a
single file

⬇

190tMultiOpt = F_or_T[1] # False (0) or True (1)

We can also have SatQuMA write out the optimiser metrics for each calculation

⬇

192tMetrics = F_or_T[1] # False (0) or True (1)

then choose the path and base filename for the output.

⬇

195outpath = '' # Path for output file (empty = current directory)

196outbase = 'out' # Name for output file (minus .csv)

Finally, we request that SatQuMA print out data to the standard output stream
as they are calculated.

⬇

200tPrint = F_or_T[1] # False (0) or True (1)

####  4.2.5 Advanced parameters

We can also specify some of the advanced parameters; these parameters have
been set to default values which should only be changed by more advanced users
looking for additional control over the calculations as they may strongly
affect the software performance and resulting SKL. We will specify that the
protocol uses the Chernoff bounds, for all tail bounds, by selecting the
default option from the list of bound functions.

⬇

214boundOpts = ['Chernoff','Hoeffding','Asymptotic']

215boundFunc = boundOpts[0] # Select an option from the list above.

Next we can select the method used to approximate the number of bits required
for error correction (EC). Again, we select the default option from the list
of EC functions, which are presented in decreasing order of complexity (i.e.
difficulty to optimise) and accuracy.

⬇

224errcorrOpts = ['logM','block','mXtot','None']

225errcorrFunc = errcorrOpts[0] # Select a method from the list above.

In order to evaluate the effect of an EC function on the optimisation process,
we can use a flag to compare the results obtained when performing EC during
and after optimisation. By default this option is turned off.

⬇

228tCompareEC = F_or_T[0] # False (0) or True (1)

We note that when comparing the effect of including the estimate of the number
of bits required for error correction both during and after optimisation only
the former case is included in the output data.

We can choose the SKL optimisation method from the three constrained
optimisation algorithms provided with scipy however we strongly recommend that
users select the default option ‘COBYLA’ to ensure robust performance.

⬇

233opt_methods = ['COBYLA','SLSQP','trust-constr']

234method = opt_methods[0] # Select a optimisation method

Due to the strong influence of the initial parameters on the final optimised
output we include an option to set the minimum number of optimisations
performed for a given set of parameters. The default value is 10, which tends
to be enough for most parameter ranges, and we would recommend a minimum of 5
for users looking to speed up their calculations.

⬇

235NoptMin = 10

The actual number of optimisations performed is also partially determined by
the number of function evaluations reported by the optimisation algorithm.
There are some additional flags which will stop the repeat optimisation loops
for a given parameter set which are best employed when the user expects that
there will be regions of zero SKL or where the optimised SKL varies little
with the initial parameters.

⬇

237tStopZero = F_or_T[1] # Stop optimizing if the first NoptMin return SKL =
0?

238tStopBetter = F_or_T[1] # Stop after NoptMin optimizations if SKL improved?

The first flag above allows the user to specify if SatQuMA should exit the
optimisation loop if the optimised SKL is consistently returned as zero. The
second flag controls whether the optimisation loop should be exited if the
optimised SKL has improved during the loop.

###  4.3 Visualising data

Once the software has successfully completed the calculation we can visualise
the output. Here we will focus on the optimised, sorted output data file
out_opt.csv. For example, we can plot the secret key length as a function of
the total system loss, both including and excluding the error correction
estimation, using the following (minimum working) python code.

⬇

1import numpy as np

2data = np.loadtxt('out_opt.csv',skiprows=1,delimiter=',')

3x = data[:,0] # Total system loss

4y1 = data[:,2] # SKL

5y1_ = y1 + data[:,7] # SKL + lambda_EC

6import matplotlib.pyplot as plt

7fig1, ax1 = plt.subplots(1,1,figsize=(5,5))

8ax1.semilogy(x,y1,'-',x,y1_,'--')

The resulting graph is shown in Fig. 4.1(a). We may also extend our plotting
code to visualise the error rates associated with these finite keys, as shown
in Fig. 4.1(b).

⬇

9fig2, ax2 = plt.subplots(1,1,figsize=(5,5))

10y2 = data[:,3:5] # QBERx and phi_x

11ax2.plot(x,y2[:,0],'-',x,y2[:,1],'--')

Finally, we may also plot the change in the optimised protocol parameters, as
shown in Fig. 4.2, using the following code.

⬇

12fig3, ax3 = plt.subplots(1,1,figsize=(5,5))

13y3 = data[:,20:27] # Protocol parameters

14ls = [':','-','-','-','--','--','--'] # Linestyles

15for ii in range(0,7,1):

16 ax3.plot(x,y3[:,ii],ls[ii])

The various other parameters and output data from this file are listed in
Table 4.1.

Figure 4.1: Total system loss in decibels against (a) secret key length, with and without error correction estimation, and (b) the phase and quantum bit error rates for the X basis.  Figure 4.2: Optimised protocol parameters (for asymmetric BB84) as a function of the total system loss in decibels. Note, $\mu_{3}=0$.  Index | Variable | Symbol |  Description
---|---|---|---
0 | ls+sysLoss | - |  Total system loss (dB).
1 | dt | $\Delta t$ |  Transmission half-window duration (s).
2 | SKL | $\ell$ |  Secret key length (bits).
3 | QBERx | $\text{QBER}_{\text{X}}$ |  Quantum bit error rate for X basis (bits/s).
4 | phi_x | $\phi_{\text{X}}$ |  Phase error rate for X basis (bits/s).
5 | nX | $n_{\text{X}}$ |  Number of events in the X basis.
6 | nZ | $n_{\text{Z}}$ |  Number of events in the Z basis.
7 | lambdaEC | $\lambda_{\text{EC}}$ |  Estimate of number of bits used for error correction.
8 | sX0 | $s_{\text{X},0}$ |  No. of vacuum events for X basis.
9 | sX1 | $s_{\text{X},1}$ |  No. of single photon events for X basis.
10 | vz1 | $v_{\text{Z},1}$ |  No. of bit errors associated with single-photon events in Z basis.
11 | sZ1 | $s_{\text{Z},1}$ |  No. of single photon events for Z basis.
12 | mpn | $\sum_{j}\mu_{j}/3$ |  Mean transmitted photon number.
13 | QBERI | $\text{QBER}_{\text{I}}$ |  Intrinsic quantum bit error rate (%).
14 | Pec | $P_{\text{ec}}$ |  Extraneous count probability.
15 | Pap | $P_{\text{ap}}$ |  Probability of after-pulse event.
16 | NoPass | $M$ |  Number of satellite overpasses.
17 | Rrate | $f_{\text{s}}$ |  Source repetition rate (Hz).
18 | eps_c | $\epsilon_{\text{c}}$ |  Correctness parameter.
19 | eps_s | $\epsilon_{\text{s}}$ |  Secrecy parameter.
20 | Px | $P_{\text{X}}$ |  Polarisation bias for X basis.
21 | P1 | $P_{\mu_{1}}$ |  Probability of sending pulse 1.
22 | P2 | $P_{\mu_{2}}$ |  Probability of sending pulse 2.
23 | P3 | $P_{\mu_{3}}$ |  Probability of sending pulse 3.
24 | mu1 | $\mu_{1}$ |  Intensity of pulse 1.
25 | mu2 | $\mu_{2}$ |  Intensity of pulse 2.
26 | mu3 | $\mu_{3}$ |  Intensity of pulse 3.
27 | xi | $\xi$ |  Offset angle of satellite orbital plane from OGS zenith (deg).
28 | min_elev | $\theta_{\text{min}}$ |  Minimum elevation of satellite for transmission (deg).
29 | max_elev | $\theta_{\text{max}}$ |  Maximum elevation of satellite overpass (deg).
30 | shift_elev | $\theta_{\text{shift}}$ |  Angle to shift transmission window central axis (deg).
Table 4.1: Table of data written to main output file(s) giving the index (data
column), python variable name, associated mathematical symbol and a brief
description.

##  Chapter 5 Bugs and future releases

SatQuMA is still very much under development with improved versions planned
for release in the near future. The current version employs an out-of-the-box
optimisation algorithm, trust-constr from scipy.

###  5.1 Known bugs

Unfortunately, as the function being minimized does not have derivatives that
can be continuously defined for all regions of the relevant parameter space
the scipy optimiser can report errors when it encounters locally flat regions
of the function space, typically where the SKL drops to zero.

###  5.2 Reporting bugs

If you encounter any behaviour you consider unexpected, or anytime python
raises an exception (excluding for user error), please send a copy of the main
SatQuMA python file, and include any relevant output files, to the development
team at the University of Strathclyde. It is always good practice to add
python files to an archive (zip/tar/etc) before sending via email to avoid
mail being rejected by mail scanning algorithms.

###  5.3 Suggested content

If there are any features that you would like to see added to this software
then please feel free to contact the development team at the University of
Strathclyde. Otherwise, if you are interested in collaborating on a project to
either develop or utilise SatQuMA then please contact the development team!

###  5.4 Contact

Please contact the development team through via Dr D. K. L. Oi at
<EMAIL_ADDRESS>
## Bibliography

  * [1] C. C. W. Lim, M. Curty, N. Walenta, F. Xu, and H. Zbinden, “Concise security bounds for practical decoy-state quantum key distribution,” Phys. Rev. A, vol. 89, p. 022307, February 2014.
  * [2] H.-L. Yin, M.-G. Zhou, J. Gu, Y.-M. Xie, Y.-S. Lu, and Z.-B. Chen, “Tight security bounds for decoy-state quantum key distribution,” Sci. Rep., vol. 10, p. 14312, August 2020.
  * [3] M. Tomamichel, J. Martinez-Mateo, C. Pacher, and D. Elkouss, “Fundamental finite key limits for one-way information reconciliation in quantum key distribution,” Quant. Inf. Proc., vol. 16, p. 280, October 2017.


#  : a Transformation Invariant Loss Function for Time-Series Forecasting

Hyunwook Lee1, Chunggi Lee2, Hongkyu Lim1, Sungahn Ko1
1Ulsan National Institute of Science and Technology, 2NAVER Webtoon
<EMAIL_ADDRESS><EMAIL_ADDRESS>Corresponding Author

###### Abstract

Time-series forecasting has caught increasing attention in the AI research
field due to its importance in solving real world problems across different
domains, such as energy, weather, traffic, and economy. As shown in various
types of data, it has been a must-see issue to deal with drastic changes,
temporal patterns, and shapes in sequential data that previous models are weak
in prediction. This is because most cases in time-series forecasting aims to
minimize $L_{p}$ norm distances as loss functions, such as mean absolute error
(MAE) or mean square error (MSE). These loss functions are vulnerable to not
only consider temporal dynamics modeling, but also capture the shape of
signals. In addition, these functions often make models misbehave and return
uncorrelated results to the original time-series. To become an effective loss
function, it has to be invariant to the set of distortions between two time-
series data instead of just comparing exact values. In this paper, we propose
a novel loss function, called (Transformation Invariant Loss function with
Distance EQuilibrium), that not only considers the distortions in amplitude
and phase but also allows models to capture the shape of time-series
sequences. In addition, supports modeling periodic and non-periodic temporal
dynamics at the same time. We evaluate the effectiveness of by conducting
extensive experiments with respect to periodic and non-periodic conditions of
data, from naive models to state-of-the-art models. The experiment results
indicate that the models trained with outperforms those trained with other
training metrics (e.g., MSE, dynamic time warping (DTW), temporal distortion
index (TDI), and longest common subsequence (LCSS)).

##  1 Introduction

Time-series forecasting has been a core problem across various domains,
including traffic domain (Li et al., 2018; Lee et al., 2020), economy (Zhu &
Shasha, 2002), and disease propagation analysis (Matsubara et al., 2014). The
crucial part of the time-series forecasting is modeling of the complex
temporal dynamics (e.g., non-stationary signal, periodicity). Temporal
dynamics, intuitively, shape, has always been one of the most attention-
getting keywords in time-series domains, such as rush hour of traffic data or
abnormal usage of the electricity (Keogh et al., 2003; Bakshi &
Stephanopoulos, 1994; Weigend & Gershenfeld, 1994; Wu et al., 2021; Zhou et
al., 2022). Deep learning methods are one of the appealing solutions to model
complex non-linear temporal dependencies and non-stationary signals, but
recent work reveals that even deep learning is often insufficient to model
temporal dynamics. To properly model the temporal dynamics, Wu et al. (2021);
Zhou et al. (2022) have proposed a novel deep learning approaches with input
sequence decomposition. Le Guen & Thome (2019) try to model sudden changes
timely and accurately with dynamic time warping (DTW). Bica et al. (2020)
adopts domain adversarial training to learn balanced representations, which is
a treatment invariant representations over time. Wu et al. (2021); Zhou et al.
(2022) have less attention to the essence of the problem: a shape, in other
words, temporal dynamics. Le Guen & Thome (2019); Bica et al. (2020) try to
capture the shape but still have some limitations like Fig. 1 (c).

A shape is a part of patterns in time-series data with a given time interval
that could give valuable information, such as rise, drop, trough, peak, and
plateau. We call the prediction is informative when it could properly consider
the shape. Time-series forecasting models should aim to both accurately
forecast the value for each time-step and the prediction should have similar
shapes as those in original time-series, but existing models do not consider
learning shape (Wu et al., 2021; Zhou et al., 2022; Bica et al., 2020; Le Guen
& Thome, 2019), so the forecasting results are often inaccurate and
uninformative, because deep learning model tends to learn in easy way (Karras
et al., 2019). Fig. 1 shows three real forecasting results with same model,
different training metrics. When we utilize mean squared error (MSE) as an
objectives, the model only aims to reduce gap between prediction and ground
truth for each time-step. As a results, the model generates relatively easy
prediction regardless of temporal dynamics (Fig. 1 (b)). It rarely gives
information about original time-series. In contrast, if we consider both gap
and shape of prediction and ground truth, the model could achieve both
accuracy and temporal dynamics, as shown in Fig. 1 (a).

In this work, we aim to design a novel objective function that guides models
to improve forecasting performance by learning the shapes in time-series data.
To design such shape-aware loss function, we review existing literature
(Esling & Agon, 2012; Bakshi & Stephanopoulos, 1994; Keogh, 2003) and
investigate the notions of shapes and distortions that interrupt measurement
for recognizing similarity of two time-series data in terms of shapes (Sec.
3.1, Sec. 3.2, and Sec. 3.3). Based on the investigation, we newly propose
required conditions for constructing an objective function for shape-aware
time-series forecasting (Sec. 3.4). We then present a novel loss function,
(Transformation Invariant Loss function with Distance EQualibrium), that
enables shape-aware representation learning with three different loss terms,
which are invariant to the distortions (Sec. 4). For evaluation, we conduct
extensive experiments with state-of-the-art deep learning models for time-
series forecasting with . The results indicate that is model-agnostic and
could improve accuracy of existing models, compared to MSE and DILATE.

Figure 1: Ground-truth and forecasting results with three metrics (a) , (b)
MSE, and (c) DTW-based loss function. (b) MSE tends to generate non-
informative forecasting results, similar to an average value of data and (c)
DTW often produces misaligned results.

#### Contributions

We make the following contributions: (1) To understand shape-awareness and
distortion invariances in time-series forecasting, we investigate existing
distortions in amplitude and phase; (2) we implement that has invariances to
many existing distortions and achieves shape-awareness and informative
forecasting in a timely manner; and (3) we show that the proposed allows
models to have higher accuracy compared to those with existing metrics such as
DTW, TDI, and LCSS on average.

##  2 Related Work

###  2.1 Time-Series Forecasting

There are many methods for time-series forecasting from traditional ones, such
as ARIMA model (Box et al., 2015) and hidden markov model (Pesaran et al.,
2004) to recent deep learning models. In this section, we briefly describe the
recent deep learning models for time-series forecasting. Starting with the
huge success of the recurrent neural networks (RNNs) (Clevert et al., 2016; Li
et al., 2018; Yu et al., 2017), researchers have developed novel deep learning
architectures, improving forecasting performance. To effectively capture long-
term dependency, which is a weakness of RNNs, Stoller et al. (2020) have
proposed convolutional neural networks (CNNs). However, it is required to
stack lots of the same CNNs to capture long-term dependency (Zhou et al.,
2021). Attention-based approaches have been another popular research direction
in time-series forecasting, including Transformer (Vaswani et al., 2017) and
Informer (Zhou et al., 2021). Although the attention-based models effectively
capture temporal dependencies, they require high computational cost and often
struggle to find proper temporal information (Wu et al., 2021). To cope with
the problem, Wu et al. (2021); Zhou et al. (2022) utilize the input
decomposition method that helps models better encode appropriate information.
The other state-of-the-art models adopt neural memory networks (Kaiser et al.,
2017; Sukhbaatar et al., 2015; Madotto et al., 2018; Lee et al., 2022), which
refer to historical data stored in memory to generate meaningful
representation.

###  2.2 Training Metrics

Conventionally, mean squared error (MSE), $L_{p}$ norm and its variants are
the mainstream to optimize forecasting models. However, they are not the best
metric to train forecasting models (Esling & Agon, 2012) because time-series
is temporally continuous. Additionally, $L_{p}$ norm gives less information
about temporal correlation among time-series data. To better model temporal
dynamics in time-series data, researchers have used differentiable,
approximated dynamic time warping (DTW), as an alternative metric of MSE
(Cuturi & Blondel, 2017; Abid & Zou, 2018; Mensch & Blondel, 2018). However,
using DTW as a loss function results in ignoring temporal localization of
changes. Recently, Le Guen & Thome (2019) suggests DILATE, a training metric
to timely catch sudden changes of non-stationary signals with smooth
approximation of DTW and penalized temporal distortion index (TDI). To
guarantee to work in a timely manner, Le Guen & Thome (2019) introduce a loss
function that gives a harsh penalty when predictions show high temporal
distortion. However, TDI relies on the DTW path, and DTW often shows
misalignment because of its noise- and scale-sensitive. Thus, DILATE often
loses its advantage with complex data, showing disadvantages at the beginning
of the training. In this work, we discuss distortions and transformation
invariances and design a new loss function that allows models to learn shapes
in the data and produce noise-robust forecasting results.

##  3 Preliminary

In this section, we aim to investigate common distortions without losing the
goal of time-series forecasting (i.e., modeling temporal dynamics and accurate
forecasting). To help understand the concepts, we first define notations and
terms (Sec. 3.1). We then discuss common distortions in time-series in
transformation perspectives that need to be considered for building a shape-
aware loss function (Sec. 3.2) and describe how other loss functions (e.g.,
DTW and TDI) handle shapes during learning (Sec. 3.3). Last, we explain the
conditions for effective time-series forecasting (Sec. 3.4).

###  3.1 Notations and Definitions

Let $X_{t}$ denote a data point at a time step $t$. Then, we can define a
time-series forecasting problem as:

#### Definition 1.

Given $T$-length historical time-series
$\mathbf{X}=[X_{t-T+1},\dots,X_{t}],X_{i}\in\mathbb{R}^{F}$ at time $t$ and
corresponding $T^{\prime}$-length future time-series
$\mathbf{Y}=[Y_{t+1},\dots,Y_{t+T^{\prime}}],Y_{i}\in\mathbb{R}^{C}$, time-
series forecasting aims to learn mapping function $f:\mathbb{R}^{T\times
F}\rightarrow\mathbb{R}^{T^{\prime}\times C}$.

To distinguish the label (i.e., ground-truth) and prediction time-series data,
we note the label data as $\mathbf{Y}$ and prediction data as
$\mathbf{\hat{Y}}$. Next, we set up two goals for time-series forecasting,
which require not only precise, but also informative forecasting Wu et al.
(2021); Zhou et al. (2022); Le Guen & Thome (2019) as follow:

  * •

Mapping function $f$ should be learnt to point-wisely reduce distance between
$\mathbf{\hat{Y}}$ and $\mathbf{Y}$; and

  * •

The output $\mathbf{\hat{Y}}$ should have similar temporal dynamics with
$\mathbf{Y}$.

Temporal dynamics are informative patterns in time-series, such as rise, drop,
trough, peak, and plateau. We define the temporal dynamics as follows:

#### Definition 2.

Temporal dynamics (or shapes) are the informative periodic and non-periodic
patterns in time-series data.

In this work, we aim to design a shape-aware loss function that satisfies both
goals. To this end, we first discuss distortions that two time-series with
similar shapes can have.

#### Definition 3.

Given two time-series $F$ and $G$ in a similar shape, distortion is a
difference between $F$ and $G$.

Distortion generally occurs in different aspects. Distortions are defined as
temporal distortion (i.e., warping) and amplitude distortion (i.e.,scaling)
with respect to its relevance of dimension, time and amplitude. Existing
distortion in data leads to misbehavior of the model, as measurements are
interrupted by the distortion. For example, if we have two time-series
$\mathbf{F}$ and $\mathbf{G}=\mathbf{F}+k$, which have a similar shape but
different means, $\mathbf{G}$ could represent many temporal dynamics of
$\mathbf{F}$. However, measurements often evaluate $\mathbf{F}$ and
$\mathbf{G}$ are different (e.g., measuring with MSE) and causes misguidance
of the model in training. As such, it is important to have measurements that
consider similar shape invariant to distortion. We define a measurement for a
distortion as follow:

#### Definition 4.

Let transformation $\mathcal{H}$ represents a distortion $H$. Then, we call
measurement $\mathcal{D}$ invariant to $\mathcal{H}$, if
$\exists\delta>0:\mathcal{D}(T,\mathcal{H}(T))<\delta$ for any time-series
$\mathbf{T}$.

Figure 2: Example of the six distortions on the amplitude axis (top) and
temporal axis (bottom).

###  3.2 Time-Series Distortions in Transformation Perspectives

Distortion, a gap between two similar time-series, affects on capturing shapes
in time-series data. As such, it is important to investigate different
distortions and their impact on representation learning aspects. There are six
common time-series distortions that models encounter during learning (Esling &
Agon, 2012; Batista et al., 2014; Berkhin, 2006; Warren Liao, 2005; Kerr et
al., 2008)–Amplitude Shifting, Phase Shifting, Uniform Amplification, Uniform
Time Scaling, Dynamic Amplification, and Dynamic Time Scaling. Next, we
explain each common time-series distortion in terms of transformation with
$n$-length time-series $\mathbf{F}(t)=[f(t_{1}),f(t_{2}),\dots,f(t_{n})]$,
where t = $[t_{1},t_{2},...,t_{n}]$. Fig. 2 presents example distortions,
categorized by amplitude and time dimensions.

  * •

Amplitude Shifting describes how much a time-series shifts against another
time-series. This can be described with two time-series and the degree of
shifting (k): $\mathbf{G}(t)=\mathbf{F}(t)+k=[f(t_{1})+k,\dots,f(t_{n})+k]$,
where $k\in\mathbb{R}$ is constant.

  * •

Phase Shifting is the same type of transformation (i.e., translation) as
amplitude shifting, but it occurs along with the temporal dimension. This
distortion can be represented with two time-series functions with the degree
of shift (k): $\mathbf{G}(t)=\textbf{F}(t+k)=[f(t_{1}+k),\dots,f(t_{n}+k)]$,
where $k\in\mathbb{R}$ is constant. Cross-correlation (Paparrizos & Gravano,
2015; Vlachos et al., 2005) is the most popular measure method that is
invariant to this distortion.

  * •

Uniform Amplification is a transformation that changes the amplitude by
multiplication of $k\in\mathbb{R}$. This distortion can be described with two
functions and a multiplication factor (k):
$\mathbf{G}(t)=k\cdot\mathbf{F}(t)=[k\cdot f(t_{1}),\dots,k\cdot f(t_{n})]$.

  * •

Uniform Time Scaling means a uniformly shortened or lengthened $\mathbf{F}(t)$
on the temporal axis. This distortion can be represented as
$\mathbf{G}(t)=[g(t_{1}),\dots,g(t_{m})]$, where $g(t_{i})=f(t_{\lceil k\cdot
i\rceil})$ and $k\in\mathbb{R}^{+}$. Although Keogh et al. (2004) propose
uniform time warping methods to handle this distortion, it still remains one
of the difficult distortion types to measure, due to the difficulty in finding
the scaling factor $k$ without testing all possible cases (Keogh, 2003).

  * •

Dynamic Amplification can be interpreted as any distortion occurred by non-
zero multiplication on the amplitude dimension. This distortion can be
described as follows:
$\mathbf{G}(t)=\mathbf{H}(t)\cdot\mathbf{F}(t)=[h(t_{1})\cdot
f(t_{1}),\dots,h(t_{n})\cdot f(t_{n})]$ with function $h(t)$ such that
$\forall_{t\in\mathbb{T}},h(t)\neq 0$. Local amplification is a representative
distortion of this type of distortions, which still remains challenging to
solve.

  * •

Dynamic Time Scaling means any transformation that dynamically lengthens or
shortens signals on the temporal dimension including local time scaling
(Batista et al., 2014) and occlusion (Batista et al., 2014; Vlachos et al.,
2003). It can be represented as follows:
$\mathbf{G}(t)=\mathbf{F}(h(t))=[f(h(t_{1})),\dots,f(h(t_{n}))]$, where $h(t)$
is a positive, strictly increasing function. Dynamic time warping (DTW)
(Bellman & Kalaba, 1959; Berndt & Clifford, 1994; Keogh & Ratanamahatana,
2005) is the most popular technique on this distortion. Das et al. (1997) also
introduce the longest common subsequence (LCSS) algorithm to tackle occlusion,
noise, and outliers in this distortion.

There are several studies on shape-aware clustering (Bellman & Kalaba, 1959;
Batista et al., 2014; Paparrizos & Gravano, 2015; Berkhin, 2006; Warren Liao,
2005; Kerr et al., 2008) and classification (Xi et al., 2006; Batista et al.,
2014; Srisai & Ratanamahatana, 2009) tasks with the consideration of shapes.
On the other hand, only a few studies exist for time-series forecasting tasks,
including Le Guen & Thome (2019) that utilizes dynamic time warping (DTW) and
temporal distortion index (TDI) for modeling temporal dynamics. Next we
describe mean square error (MSE) and DILATE, proposed by Le Guen & Thome
(2019), and discuss their invariance to the distortions.

###  3.3 Distortion Handling in Current Time-Series Forecasting Objectives

Many measurement metrics have been used in the time-series forecasting domain,
and those based on the $L_{p}$ distance, including Euclidean distance, are
widely used to handle time-series data. However, such metrics do not have
invariance to the aforementioned distortions (Ding et al., 2008; Le Guen &
Thome, 2019) due to its point-wise mapping. Specifically, since $L_{p}$
distance compares the values per time step, it cannot handle temporal
distortions appropriately and vulnerable to scaling of the data. Le Guen &
Thome (2019) propose a loss function, called DILATE, to overcome the
inadequate characteristic in the $L_{p}$ distance metrics by recognizing
temporal dynamics with DTW and TDI. In terms of transformation, DILATE handles
dynamic time scaling, especially, local time scaling with DTW, and phase
shifting with penalized TDI, defined as follows:

$\mathcal{L}_{{DILATE}}(\hat{y}_{i},y_{i}):=-\gamma\log\Big{(}\sum_{\mathbf{A}\in\mathcal{A}_{k,k}}\exp\Big{(}-\frac{\langle\mathbf{A},\alpha\Delta(\hat{y}_{i},y_{i})+(1-\alpha)\mathbf{\Omega}\rangle}{\gamma}\Big{)}\Big{)},$

where $\mathbf{A}$, $\delta(\hat{y}_{i},y_{i})$, $\mathbf{\Omega}$ are the
warping path, cost matrix, and squared penalization matrix, respectively.

While DILATE shows better performance than existing methods, there is a
missing point in invariance point of view. Basically, DTW computes the
Euclidean distance of two time-series after its temporal alignment in dynamic
programming and the alignment relies on the distance function. Consequently,
the dynamic alignment of the DTW can be properly achieved only when two time-
series have the same range (Esling & Agon, 2012; Bellman & Kalaba, 1959). That
means, it hardly achieves invariance on amplitude distortion without
appropriate pre-processing. Gong & Chen (2017) also show that DTW poorly
matches the prediction and target (i.e., ground-truth) time-series with
amplitude shifting. Even when the target time-series is aligned with
normalization, we cannot guarantee that the predicted and target time-series
are properly aligned due to DTW’s high sensitivity to noise. As a result,
DILATE can generate poor alignment results that can cause wrong optimization
of TDI, which produces instability during optimization steps and incorrect
results. To design an effective shape-aware loss function, we have to
understand measures and when the measures have transformation invariances. In
the next section, we discuss how we interpret transformations in time-series
forecasting point of view and which types of transformations should be
considered in objective function design.

###  3.4 Transformation Invariances in Time-Series Forecasting

In the time-series domain, data often have various distortions so measurements
are needed to satisfy a number of transformation invariances for meaningfully
modeling temporal dynamics. As discussed in Sec. 3.1, we set the goal of time-
series forecasting as (1) point-wisely reducing the gap between prediction and
target time-series and (2) preserving temporal dynamics of the target time-
series. To satisfy both of them, we have to consider (1) the method that
should not have a negative impact on the traditional goal of accurate time-
series forecasting and (2) the distortions that play a crucial role in
capturing the temporal dynamics of the target time-series. In this section, we
review all six distortions whether it is a feasible loss function or not,
discuss their benefits and trade-offs, and find appropriate distortions to be
considered in time-series forecasting.

#### Amplitude Shifting

In a wide range of situations, it is beneficial to capture the the trends of
time-series sequence in spite of shifts in terms of amplitude. Thus, being
invariant to amplitude shifting in a loss function takes many advantages in
time-series forecasting: (1) shape-awareness invariant to amplitude shifting,
(2) accurate deviation of values in modeling, and (3) effective on-time
prediction of the peak or sudden changes. To guarantee the amplitude shifting
invariant in the optimization stage, the loss function should induce an equal
gap $k$ between prediction and ground truth data in each step. Formally
speaking, the loss function with consideration of the amplitude shifting
should satisfy:

$\mathcal{L}(\mathbf{Y},\mathbf{\hat{Y}})=0\Leftrightarrow\forall_{i\in[1,\dots,n]},d(y_{i},\hat{y}_{i})=k,$
(1)

where $k\in\mathbb{R}$ is an arbitrary and equal gap, and
$d(y_{i},\hat{y}_{i})$ is a signed distance with a boundary
$y_{i}>\hat{y}_{i}$. By allowing tolerance between prediction and target time-
series, models can follow trends in time-series instead of tending to predict
exact values in point-wise. In short, unlike existing loss functions that
handle only point-wise distance (e.g., DTW), we should deal with both the
point-wise distance and its relational distance values to guarantee amplitude
shifting.

#### Phase Shifting

There are forecasting tasks, whose main objectives concern accurate
forecasting of peaks and periodicity in time-series (e.g., heart beat data and
stock price data). For such tasks, phase shifting invariance is one of the
best solutions for (1) modeling periodicity, regardless of translation on
temporal axis and (2) having precise statistics with shapes, such as peak and
plateau values. If a loss function is to be invariant to phase shifting, the
function should satisfy:

$\mathcal{L}(\mathbf{Y},\mathbf{\hat{Y}})=0\Leftrightarrow\mathbf{Y}\textrm{
and }\mathbf{\hat{Y}}\textrm{ have the same dominant frequency}.$ (2)

Note Eq. 2 allows a similar shape as target time-series in forecasting, not
exactly same shape (e.g., $\sin(x)\textrm{ and }2\sin(x+x_{0})$ with the same
dominant frequency).

#### Uniform Amplification

This proposition will be useful in case of sparse data that contains a
significant number of zeros. By adopting the uniform amplification invariance,
models are able to focus non-zero sequences, whereas this proposition allows
models to receive less penalty in zero sequences. Since it guarantees shape-
awareness with a multiplication factor in a timely manner as Fig. 2,
invariance for uniform amplification fits well. To have a model trained with
the uniform amplification invariance, the loss function should satisfy:

$\mathcal{L}(\mathbf{Y},\mathbf{\hat{Y}})=0\Leftrightarrow\forall_{i\in[1,\dots,n]},\frac{y_{i}}{\hat{y}_{i}}=k(\hat{y}_{i}\neq
0).$ (3)

#### Uniform Time Scaling, Dynamic Amplification, and Dynamic Time Scaling

After careful consideration, we conclude that uniform time scaling, dynamic
amplification, and dynamic time scaling are incompatible for optimization. We
describe the reason below.

To achieve invariance for the uniform time scaling, the loss function should
satisfy:

$\displaystyle\mathcal{L}(\mathbf{Y},\mathbf{\hat{Y}})=0$
$\displaystyle\Leftrightarrow\exists
c\in\mathbb{Z}^{+}:\\{c|y_{i}=\hat{y}_{ci}\\}\cup\\{c|y_{ci}=\hat{y}_{i}\\}\forall
i\in[0,1,\dots,T^{\prime}].$

This proposition will influence negatively original temporal dynamics,
considering that it gives the tolerance of mispredicting periodicity (e.g.,
daily periodic signals) and even cannot catch events (e.g., abrupt changing
values) in timely manner. In summary, it hinders models from capturing shape
and corrupts periodic information.

For both dynamic amplification and dynamic time scaling, loss functions always
are zero for all pairs when we do not set the limit of tolerance. For example,
if we do not limit tolerance, the proposition for dynamic amplification
invariance is as follow:

$\displaystyle\mathcal{L}(\mathbf{Y},\mathbf{\hat{Y}})=0\Leftrightarrow\forall
c_{i}\in\mathbb{R}:y_{i}=c_{i}\hat{y}_{i},$

If a loss function satisfies the proposition, it is always zero because there
always exists $c_{i}=y_{i}/\hat{y}_{i}$ except $\hat{y}_{i}=0$. Therefore, it
is not able to give any information because all random values could be an
optimal solution. The same situation happens with the dynamic time scaling if
we do not limit the window. Consequently, all of the uniform time scaling,
dynamic amplification, and dynamic time scaling are unsuitable to be
objectives in time-series forecasting.

##  4 Methods

In this section, we describe a novel loss function (a Transformation Invariant
Loss function with Distance EQuilibrium), which allows models to perform
shape-aware time-series forecasting based on the three distortion invariances.
To build a transformation invariant loss function, we have to design a loss
function that satisfies the proposition for amplitude shifting invariance (Eq.
1), phase shifting invariance (Eq. 2), and uniform amplification shifting
invariance (Eq. 3), as discussed in Sec. 3.4. We select them for our loss
function because they help models capture the shape and do not harm the goal
of the traditional time-series forecasting (i.e., minimize gap between
prediction and target time-series). Not only the loss function should satisfy
these propositions, but also it should consider correlations between the whole
sequence of outputs and ground truths rather than point-wisely optimizing the
model. It is not achieved by other loss functions, such as MSE or DILATE. To
handle all three distortions and the whole sequence of correlations, we build
three objective functions (a.shift, phase, and amp losses) that achieve one or
more invariance by utilizing softmax, Fourier coefficient, and auto-
correlation to design a loss function.

#### Amplitude Shifting Invariance with Softmax (Amplitude Shifting)

To strengthen amplitude shifting invariance, we design a loss function that
satisfies Eq. 1. This means, $d(y_{i},\hat{y}_{i})$ needs to be the same value
for all $i$. To satisfy the condition, we utilize the softmax function:

$\mathcal{L}_{a.shift}(\mathbf{Y},\mathbf{\hat{Y}})=T^{\prime}\sum_{i=1}^{T^{\prime}}|\frac{1}{T^{\prime}}-\textrm{Softmax}(d(y_{i},\hat{y}_{i})))|,\textrm{Softmax}(d(y_{i},\hat{y}_{i}))=\frac{e^{d(y_{i},\hat{y}_{i})}}{\sum_{j=1}^{T^{\prime}}e^{d(y_{j},\hat{y}_{j})}}$
(4)

where $T^{\prime}$, Softmax, and $d(\cdot,\cdot)$ are the length of sequence,
softmax function, and signed distance function, respectively. Because the
Softmax produces the proportion of each value, it only reaches to the optimal
solution when it satisfies Eq. 1. Also, if we utilize Softmax, there is no
need to know arbitrary equal gap $k$.

#### Invariances with Fourier Coefficients (Phase Shifting)

As we discussed in Sec. 3.4, one candidate method to obtain phase shifting
invariance is to use Fourier coefficients. As described in prior studies (Ng &
Goldberger, 2007), we can reconstruct original time-series only with dominant
frequencies. In this way, we utilize the norm of dominant Fourier coefficient
of ground truth and prediction sequences as our additional objective function,
achieving phase shifting invariance. When it comes to the other frequencies,
we denote the norm of prediction sequence to reduce the value of Fourier
coefficient. Consequently, with the help of our loss function, this loss
function allows model to be noise robustness because the Fourier coefficients
of white noises in original time-series are relatively small. Simply, we
optimize the distance between Fourier coefficients of two time-series as:

$\mathcal{L}_{phase}(\mathbf{Y},\mathbf{\hat{Y}})=\begin{cases}||\mathcal{F}(\mathbf{Y})-\mathcal{F}(\mathbf{\hat{Y}})||_{p},&\text{if
dominant frequency}\\\
||\mathcal{F}(\mathbf{\hat{Y}})||_{p},&\text{otherwise}\end{cases}$ (5)

where $||\cdot||_{p}$ is the $L_{p}$ norm. This loss function obtains uniform
amplification invariance by utilizing a normalization technique to Fourier
coefficients. For example, $\sin x$ and $c\cdot\sin x$ have the same Fourier
coefficients if properly normalized. In summary, by Eq. 5, we could obtain (1)
invariance for phase shifting, (2) invariance for uniform amplification, and
(3) robustness to noise.

#### Invariances with auto-correlation (Uniform Amplification)

Although Fourier coefficients can be considered as a reasonable solution to
catch the periodicity of the target time-series, it is not fully invariant to
phase shifting for three reasons–(1) the statistics (e.g., mean and variance)
in data keep changing, (2) such changing statistics also cause the changes of
Fourier coefficients even in the same frequency, and (3) objectives only with
a norm of them cannot fully represent the original time-series. Thus, we
introduce an objective based on normalized cross-correlation, which satisfies
Eq. 2 for a periodic signal:

$\mathcal{L}_{amp}(\mathbf{Y},\mathbf{\hat{Y}})=||R(\mathbf{Y},\mathbf{Y})-R(\mathbf{Y},\mathbf{\hat{Y}})||_{p},$
(6)

where $R(\cdot,\cdot)$ is a normalized cross correlation function. This loss
function helps predicted sequences to mimic label sequences by calculating
difference between the auto-correlation of the label sequences and cross-
correlation between label and predicted sequences. Therefore, the label and
prediction have similar temporal dynamics regardless of phase shifting and
uniform amplification.

In summary, we introduce (Transformation Invariant Loss Function with Distance
Equilibrium), combining Eq. 4, Eq. 5, and Eq. 6 as follows:

$\mathcal{L}_{TILDEq}(\mathbf{Y},\mathbf{\hat{Y}})=\alpha\mathcal{L}_{a.shift}(\mathbf{Y},\mathbf{\hat{Y}})+(1-\alpha)\mathcal{L}_{phase}(\mathbf{Y},\mathbf{\hat{Y}})+\gamma\mathcal{L}_{amp}(\mathbf{Y},\mathbf{\hat{Y}}),$
(7)

where $\alpha\in[0,1]$ and $\gamma$ is hyperparameter.

##  5 Experiments

In this section, we present the results of our comprehensive experiments,
demonstrating the effectiveness of and importance of transformation
invariance.

Table 1: Experimental results of short-term time-series forecasting with
Seq2Seq GRU model.

Methods | GRU + MSE | GRU + DILATE | GRU +
---|---|---|---
Eval | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS
Synthetic | 0.0107 | 3.5080 | 1.0392 | 0.3523 | 0.0130 | 3.4005 | 1.1242 | 0.3825 | 0.0119 | 3.2873 | 1.1564 | 0.3811
ECG5000 | 0.2152 | 1.9718 | 0.8442 | 0.7743 | 0.8270 | 3.9579 | 2.0281 | 0.4356 | 0.2141 | 1.9575 | 0.7714 | 0.7773
Traffic | 0.0070 | 1.4628 | 0.2343 | 0.7209 | 0.0095 | 1.6929 | 0.2814 | 0.6806 | 0.0072 | 1.4600 | 0.2276 | 0.7220

#### Experimental Setup

We conduct the experiments with three state-of-the-art models, including
Informer (Zhou et al., 2021), N-Beats (Oreshkin et al., 2020), and Autoformer
(Wu et al., 2021) and one simple sequence-to-sequence gated recurrent unit
(GRU) model. We use five real-world datasets–ECG5000, Traffic, ETTh2, ETTm2,
and ECL, and one synthetic dataset–Synthetic for model training. We repeat
each experiment with a model and dataset 10 times in combinations with three
different objective functions. Appendix A provides detailed explanations on
the datasets, hyperparameter setting, and model architectures.

Table 2: Experimental results on three real-world datasets (four cases) with
Informer.

Methods | Informer + MSE | Informer + DILATE | Informer +
---|---|---|---
Metric | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS
ETTh2 | 96 | 0.2466 | 6.9254 | 3.6676 | 0.4633 | 0.3284 | 6.3109 | 3.5838 | 0.5037 | 0.1768 | 5.8437 | 1.6734 | 0.5379
192 | 0.2818 | 10.2654 | 11.1580 | 0.4254 | 0.4086 | 8.8262 | 7.1780 | 0.4893 | 0.2432 | 10.2134 | 9.9865 | 0.4317
336 | 0.3089 | 12.1822 | 18.7014 | 0.4434 | 0.4164 | 10.3779 | 13.2580 | 0.5062 | 0.2958 | 13.5586 | 20.2850 | 0.4165
720 | 0.2877 | 17.6369 | 38.4617 | 0.4425 | 0.4229 | 14.1196 | 23.9403 | 0.4815 | 0.3157 | 18.4617 | 43.3238 | 0.4262
ETTm2 | 96 | 0.0889 | 3.4007 | 1.5719 | 0.7386 | 0.1263 | 6.0144 | 2.7757 | 0.5129 | 0.0871 | 3.1354 | 1.3474 | 0.7817
192 | 0.1157 | 5.7964 | 2.8128 | 0.6705 | 0.2340 | 9.7004 | 7.8354 | 0.5266 | 0.1317 | 5.7093 | 2.9129 | 0.6983
336 | 0.1860 | 8.9971 | 6.7970 | 0.6365 | 0.2805 | 11.7889 | 13.3861 | 0.5025 | 0.1767 | 9.0866 | 7.4023 | 0.6555
720 | 0.2165 | 14.7685 | 24.6694 | 0.5768 | 0.3745 | 16.7734 | 29.2783 | 0.4747 | 0.2063 | 15.3057 | 24.1959 | 0.5860
ECL | 96 | 0.2709 | 2.8067 | 0.1720 | 0.7032 | 0.9856 | 3.6394 | 1.4794 | 0.6324 | 0.2800 | 2.9466 | 0.2473 | 0.7275
192 | 0.2793 | 4.1193 | 0.1508 | 0.7060 | 1.1209 | 5.2289 | 2.1749 | 0.6053 | 0.3077 | 4.2693 | 0.2978 | 0.7336
336 | 0.3203 | 5.9533 | 0.1642 | 0.7222 | 1.2331 | 7.8470 | 3.0415 | 0.5694 | 0.3271 | 5.8090 | 0.1984 | 0.7143
720 | 0.6414 | 15.8561 | 4.4284 | 0.4564 | 1.3706 | 12.5981 | 5.6720 | 0.5506 | 0.4676 | 11.4027 | 0.7107 | 0.6298

#### Evaluation Metrics

In the experiment, we evaluate with four evaluation metrics: mean squared
error (MSE), dynamic time warping (DTW), its corresponding temporal distortion
index (TDI), all of which are used in Le Guen & Thome (2019). As DTW is
sensitive to noise and generates incorrect paths when one of the time-series
data is noisy (as discussed in Sec. 3.3), we additionally use the longest
common subsequence (LCSS) for comparison, which is more robust to outliers and
noise (Esling & Agon, 2012). The longer the length of matched subsequences is
achieved, the better performance LCSS shows in modeling the shapes.

Table 3: Experimental results on three real-world datasets (four cases) with
N-Beats.

Methods | N-Beats + MSE | N-Beats + DILATE | N-Beats +
---|---|---|---
Metric | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS
ETTh2 | 96 | 0.1869 | 7.2379 | 2.3787 | 0.4688 | 0.3105 | 6.5849 | 3.6490 | 0.4879 | 0.1557 | 5.1011 | 1.3240 | 0.5862
192 | 0.2385 | 11.5667 | 4.9153 | 0.4505 | 0.6186 | 9.7254 | 7.0831 | 0.4637 | 0.1738 | 7.6334 | 2.4122 | 0.5819
336 | 0.2889 | 16.5255 | 11.5207 | 0.4544 | 1.1406 | 13.7328 | 14.6986 | 0.4584 | 0.2132 | 11.3351 | 5.3556 | 0.5373
720 | 0.3881 | 24.1570 | 18.8462 | 0.4381 | 1.6713 | 19.4392 | 23.7028 | 0.4575 | 0.3044 | 17.6006 | 9.6636 | 0.5287
ETTm2 | 96 | 0.0790 | 3.9685 | 2.0436 | 0.6721 | 0.1524 | 7.9302 | 5.5597 | 0.4379 | 0.0952 | 4.0110 | 2.1939 | 0.6902
192 | 0.1224 | 6.8695 | 3.2834 | 0.5762 | 0.2055 | 10.0393 | 8.5602 | 0.5107 | 0.1286 | 6.3556 | 4.9798 | 0.6160
336 | 0.1824 | 12.1438 | 8.5915 | 0.4587 | 0.2501 | 12.6342 | 16.1473 | 0.4819 | 0.1705 | 8.9377 | 8.3539 | 0.6195
720 | 0.2370 | 22.8676 | 17.8458 | 0.4929 | 0.4170 | 17.7764 | 24.6877 | 0.5836 | 0.2336 | 14.2715 | 19.0883 | 0.7070
ECL | 96 | 0.3666 | 3.5207 | 0.2989 | 0.6589 | 1.1156 | 5.1430 | 2.6613 | 0.5074 | 0.3183 | 2.9707 | 0.4844 | 0.7229
192 | 0.4307 | 5.7578 | 0.4253 | 0.6212 | 1.1859 | 7.3406 | 2.8488 | 0.4973 | 0.3383 | 4.1817 | 0.4229 | 0.7187
336 | 0.5199 | 8.5563 | 0.5384 | 0.5965 | 1.2460 | 9.5096 | 3.0517 | 0.5091 | 0.3831 | 5.6643 | 0.3024 | 0.7112
720 | 0.6240 | 13.9436 | 0.6510 | 0.5717 | 1.3061 | 13.1928 | 3.7279 | 0.5337 | 0.4540 | 8.9997 | 0.3251 | 0.6960

#### Results and Analysis

Table 1 shows the results of short-term forecasting performance of gated
recurrent unit (GRU) optimized with MSE, DILATE, and metrics. Synthetic,
ECG5000, and Traffic datasets are used for the experiment. With the Synthetic
dataset, every used metric shows its own benefits. This result indicates that
similarity of the shape and MSE measures have a clear advantage when a model
is trained and evaluated with themselves. Also, since the model is evaluated
with real-world datasets, it is revealed outperforms other objective functions
in most evaluation metrics. These results indicate our approach for learning
shapes in time-series data works better than existing methods for forecasting.
DILATE does not show impressive performance with ECG5000 due to its high
sensitiveness to noise, as discussed in Sec. 3.3.

Table 2 and Table 3 summarize the experiment results with two state-of-the-art
models, Informer and N-Beats. The models make predictions for both short-term
($L$=24) and long-term ($L$ up to 720), so that we can investigate their
performances with different forecasting difficulties. In most of datasets, the
models with outperform those with other training metrics. Especially for long-
term forecasting, we observe that For N-Beats and Informer with significantly
improve the performance with the other metrics. We provide some visual
examples in Appendix B and more detailed analysis, qualitative experiments
with example visualizations, ablation study results in Appendix B. This result
implies that improves performances of the models in learning temporal
dynamics, including LCSS of N-Beats (improved over 10%).

Table 4: Experimental results on three real-world datasets (four cases) with
Autoformer.

Methods | Autoformer + MSE | Autoformer + DILATE | Autoformer +
---|---|---|---
Metric | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS | MSE | DTW | TDI | LCSS
ETTh2 | 96 | 0.1538 | 5.2227 | 2.1865 | 0.6187 | 0.2211 | 6.0453 | 2.5345 | 0.5315 | 0.1494 | 5.1060 | 1.9752 | 0.6317
192 | 0.1974 | 7.8730 | 3.3382 | 0.6019 | 0.2825 | 8.6696 | 5.6671 | 0.5335 | 0.2079 | 7.8917 | 3.7532 | 0.5984
336 | 0.2393 | 10.8002 | 7.3141 | 0.5954 | 0.3759 | 11.0335 | 13.1347 | 0.5257 | 0.2360 | 10.7212 | 7.0085 | 0.5971
720 | 0.2859 | 16.3502 | 15.9233 | 0.5772 | 0.4296 | 15.9819 | 22.2173 | 0.4924 | 0.2378 | 16.0002 | 13.7906 | 0.5795
ETTm2 | 96 | 0.0990 | 4.3498 | 2.5052 | 0.6756 | 0.1135 | 5.3097 | 2.2211 | 0.5936 | 0.0940 | 3.9078 | 2.2587 | 0.7075
192 | 0.1340 | 6.3207 | 3.3676 | 0.6512 | 0.1854 | 8.5209 | 3.7894 | 0.5506 | 0.1259 | 6.0979 | 2.9278 | 0.6810
336 | 0.1587 | 9.4374 | 6.9205 | 0.6036 | 0.2001 | 12.0265 | 8.8305 | 0.5370 | 0.1548 | 9.5223 | 7.2875 | 0.6169
720 | 0.1999 | 14.8332 | 11.9655 | 0.6064 | 0.2665 | 17.8025 | 17.4114 | 0.5001 | 0.1885 | 14.5844 | 9.9918 | 0.6277
ECL | 96 | 0.4209 | 3.5957 | 0.2461 | 0.6487 | 0.6813 | 3.6490 | 0.4780 | 0.6253 | 0.3515 | 3.2173 | 0.2298 | 0.6912
192 | 0.4206 | 4.9924 | 0.3416 | 0.6574 | 0.7319 | 5.5324 | 0.2775 | 0.6118 | 0.4032 | 4.8581 | 0.3301 | 0.6680
336 | 0.4621 | 6.6888 | 0.2795 | 0.6535 | 0.7895 | 7.5665 | 0.2503 | 0.6091 | 0.4637 | 6.7335 | 0.3923 | 0.6429
720 | 0.5005 | 10.8571 | 0.2383 | 0.6183 | 0.8630 | 12.1416 | 0.1877 | 0.6074 | 0.5049 | 9.8492 | 0.2525 | 0.6420

##  6 Conclusion and Future Work

We propose , a transformation invariant loss function with distance
equilibrium, which allows shape-aware time-series forecasting in a timely
manner. To design , we review existing transformations in time-series data and
discuss the conditions that ensure transformation invariances during
optimization tasks. The designed ensures a model to be invariant to the
amplitude shifting, phase shifting, and uniform amplification so that a model
better captures the shape in time-series data. To prove the effectiveness of ,
we conduct comprehensive experiments with state-of-the-art models and real-
world datasets. The results indicate that the model trained with generates
more timely, robust, accurate, and shape-aware forecasting in both short-term
to long-term forecasting tasks. We conjecture that this work can facilitate
future research on transformation invariances and shape-aware forecasting.

## References

  * Abid & Zou (2018) Abubakar Abid and James Y Zou.  Learning a warping distance from unlabeled time series using sequence autoencoders.  In _Advances in Neural Information Processing Systems_ , volume 31, pp. 10568–10578, 2018.
  * Bakshi & Stephanopoulos (1994) B.R. Bakshi and G. Stephanopoulos.  Representation of process trends—iv. induction of real-time patterns from operating data for diagnosis and supervisory control.  _Computers & Chemical Engineering_, 18(4):303–332, 1994.
  * Batista et al. (2014) Gustavo E. A. P. A. Batista, Eamonn J. Keogh, Oben Moses Tataw, and Vinícius M. A. de Souza.  CID: an efficient complexity-invariant distance for time series.  _Data Mining and Knowledge Discovery_ , 28(3):634–669, 2014.  doi: 10.1007/s10618-013-0312-3.
  * Bellman & Kalaba (1959) R. Bellman and R. Kalaba.  On adaptive control processes.  _IRE Transactions on Automatic Control_ , 4(2):1–9, 1959.
  * Berkhin (2006) Pavel Berkhin.  A survey of clustering data mining techniques.  In _Grouping Multidimensional Data - Recent Advances in Clustering_ , pp. 25–71. Springer, 2006.
  * Berndt & Clifford (1994) Donald J. Berndt and James Clifford.  Using dynamic time warping to find patterns in time series.  In _Proceedings of the International Conference on Knowledge Discovery and Data Mining_ , AAAIWS’94, pp. 359–370. AAAI Press, 1994.
  * Bica et al. (2020) Ioana Bica, Ahmed M Alaa, James Jordon, and Mihaela van der Schaar.  Estimating counterfactual treatment outcomes over time through adversarially balanced representations.  In _International Conference on Learning Representations_ , 2020.
  * Box et al. (2015) George E. P. Box, Gwilym M. Jenkins, Gregory C. Reinsel, and Greta M. Ljung.  _Time series analysis: forecasting and control_.  John Wiley, 2015.
  * Clevert et al. (2016) Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  In _Proceedings of the International Conference on Learning Representations_ , 2016.
  * Cuturi & Blondel (2017) Marco Cuturi and Mathieu Blondel.  Soft-dtw: A differentiable loss function for time-series.  In _Proceedings of the 34th International Conference on Machine Learning_ , ICML’17, pp. 894–903, 2017.
  * Das et al. (1997) Gautam Das, Dimitrios Gunopulos, and Heikki Mannila.  Finding similar time series.  In _Principles of Data Mining and Knowledge Discovery_ , pp. 88–100, 1997.
  * Dau et al. (2019) Hoang Anh Dau, Anthony Bagnall, Kaveh Kamgar, Chin-Chia Michael Yeh, Yan Zhu, Shaghayegh Gharghabi, Chotirat Ann Ratanamahatana, and Eamonn Keogh.  The ucr time series archive.  _IEEE/CAA Journal of Automatica Sinica_ , 6(6):1293–1305, 2019.
  * Ding et al. (2008) Hui Ding, Goce Trajcevski, Peter Scheuermann, Xiaoyue Wang, and Eamonn Keogh.  Querying and mining of time series data: Experimental comparison of representations and distance measures.  _Proceedings of the VLDB Endowment_ , 1(2):1542–1552, 2008.
  * Esling & Agon (2012) Philippe Esling and Carlos Agon.  Time-series data mining.  _ACM Computing Surveys_ , 45(1), 2012.
  * Gong & Chen (2017) Zhichen Gong and Huanhuan Chen.  Dynamic state warping.  _CoRR_ , abs/1703.01141, 2017.
  * Kaiser et al. (2017) Lukasz Kaiser, Ofir Nachum, Aurko Roy, and Samy Bengio.  Learning to remember rare events.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_ , 2017.
  * Karras et al. (2019) Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * Keogh (2003) Eamonn J. Keogh.  Efficiently finding arbitrarily scaled patterns in massive time series databases.  In _Knowledge Discovery in Databases: PKDD 2003_ , volume 2838 of _Lecture Notes in Computer Science_ , pp. 253–265, 2003.
  * Keogh & Ratanamahatana (2005) Eamonn J. Keogh and Chotirat (Ann) Ratanamahatana.  Exact indexing of dynamic time warping.  _Knowledge and Information Systems_ , 7(3):358–386, 2005.
  * Keogh et al. (2003) Eamonn J. Keogh, Jessica Lin, and Wagner Truppel.  Clustering of time series subsequences is meaningless: Implications for previous and future research.  In _Proceedings of the IEEE International Conference on Data Mining_ , pp. 115–122. IEEE Computer Society, 2003.
  * Keogh et al. (2004) Eamonn J. Keogh, Themis Palpanas, Victor B. Zordan, Dimitrios Gunopulos, and Marc Cardle.  Indexing large human-motion databases.  In _Proceedings of the International Conference on Very Large Data Bases_ , pp. 780–791, 2004.
  * Kerr et al. (2008) G. Kerr, H.J. Ruskin, M. Crane, and P. Doolan.  Techniques for clustering gene expression data.  _Computers in Biology and Medicine_ , 38(3):283–293, 2008.
  * Le Guen & Thome (2019) Vincent Le Guen and Nicolas Thome.  Shape and time distortion loss for training deep time series forecasting models.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Lee et al. (2020) Chunggi Lee, Yeonjun Kim, Seungmin Jin, Dongmin Kim, Ross Maciejewski, David Ebert, and Sungahn Ko.  A visual analytics system for exploring, monitoring, and forecasting road traffic congestion.  _IEEE Transactions on Visualization and Computer Graphics_ , 26(11):3133–3146, 2020.  doi: 10.1109/TVCG.2019.2922597.
  * Lee et al. (2022) Hyunwook Lee, Seungmin Jin, Hyeshin Chu, Hongkyu Lim, and Sungahn Ko.  Learning to remember patterns: Pattern matching memory networks for traffic forecasting.  In _International Conference on Learning Representations_ , 2022.
  * Li et al. (2018) Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  In _Proceedings of the International Conference on Learning Representations_. OpenReview.net, 2018.
  * Madotto et al. (2018) Andrea Madotto, Chien-Sheng Wu, and Pascale Fung.  Mem2seq: Effectively incorporating knowledge bases into end-to-end task-oriented dialog systems.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics, ACL 2018, Melbourne, Australia, July 15-20, 2018, Volume 1: Long Papers_ , pp. 1468–1478, 2018.
  * Matsubara et al. (2014) Yasuko Matsubara, Yasushi Sakurai, Willem G. van Panhuis, and Christos Faloutsos.  FUNNEL: automatic mining of spatially coevolving epidemics.  In _The ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pp. 105–114. ACM, 2014.
  * Mensch & Blondel (2018) Arthur Mensch and Mathieu Blondel.  Differentiable dynamic programming for structured prediction and attention.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 3462–3471, 2018.
  * Ng & Goldberger (2007) Jason Ng and Jeffrey J Goldberger.  Understanding and interpreting dominant frequency analysis of af electrograms.  _Journal of Cardiovascular Electrophysiology_ , 18(6):680–685, 2007.
  * Oreshkin et al. (2020) Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-beats: Neural basis expansion analysis for interpretable time series forecasting.  In _International Conference on Learning Representations_ , 2020.
  * Paparrizos & Gravano (2015) John Paparrizos and Luis Gravano.  K-shape: Efficient and accurate clustering of time series.  In _Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data_ , SIGMOD ’15, pp. 1855–1870, 2015.  doi: 10.1145/2723372.2737793.
  * Pesaran et al. (2004) M.H. Pesaran, D. Pettenuzzo, and A. Timmermann.  Forecasting time series subject to multiple structural breaks.  Cambridge Working Papers in Economics 0433, Faculty of Economics, University of Cambridge, 2004.
  * Srisai & Ratanamahatana (2009) Dararat Srisai and Chotirat Ann Ratanamahatana.  Efficient time series classification under template matching using time warping alignment.  In _Proceedings of the International Conference on Computer Sciences and Convergence Information Technology_ , pp. 685–690, 2009.
  * Stoller et al. (2020) Daniel Stoller, Mi Tian, Sebastian Ewert, and Simon Dixon.  Seq-u-net: A one-dimensional causal u-net for efficient sequence modelling.  In Christian Bessiere (ed.), _Proceedings of the International Joint Conference on Artificial Intelligence_ , pp. 2893–2900. ijcai.org, 2020\.
  * Sukhbaatar et al. (2015) Sainbayar Sukhbaatar, arthur szlam, Jason Weston, and Rob Fergus.  End-to-end memory networks.  In _Advances in Neural Information Processing Systems_ , volume 28, 2015.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Isabelle Guyon, Ulrike von Luxburg, Samy Bengio, Hanna M. Wallach, Rob Fergus, S. V. N. Vishwanathan, and Roman Garnett (eds.), _Advances in Neural Information Processing Systems_ , pp. 5998–6008, 2017.
  * Vlachos et al. (2003) Michail Vlachos, Marios Hadjieleftheriou, Dimitrios Gunopulos, and Eamonn Keogh.  Indexing multi-dimensional time-series with support for multiple distance measures.  In _Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’03, pp. 216–225, 2003.
  * Vlachos et al. (2005) Michail Vlachos, Philip S. Yu, and Vittorio Castelli.  On periodicity detection and structural periodic similarity.  In _Proceedings of the SIAM International Conference on Data Mining_ , pp. 449–460, 2005.
  * Warren Liao (2005) T. Warren Liao.  Clustering of time series data—a survey.  _Pattern Recognition_ , 38(11):1857–1874, 2005\.
  * Weigend & Gershenfeld (1994) Andreas S. Weigend and Neil A. Gershenfeld.  _Time Series Prediction: Forecasting the Future and Understanding the Past_.  Addison-Wesley, 1994.  ISBN 0-201-62601-2.
  * Wu et al. (2021) Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long.  Autoformer: Decomposition transformers with auto-correlation for long-term series forecasting.  In _Advances in Neural Information Processing Systems_ , volume 34, pp. 22419–22430, 2021.
  * Xi et al. (2006) Xiaopeng Xi, Eamonn Keogh, Christian Shelton, Li Wei, and Chotirat Ann Ratanamahatana.  Fast time series classification using numerosity reduction.  In _Proceedings of the International Conference on Machine Learning_ , ICML ’06, pp. 1033–1040. Association for Computing Machinery, 2006\.
  * Yu et al. (2017) Fisher Yu, Vladlen Koltun, and Thomas A. Funkhouser.  Dilated residual networks.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 636–644. IEEE Computer Society, 2017.
  * Zhou et al. (2021) Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang.  Informer: Beyond efficient transformer for long sequence time-series forecasting.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35(12):11106–11115, 2021.
  * Zhou et al. (2022) Tian Zhou, Ziqing Ma, Qingsong Wen, Xue Wang, Liang Sun, and Rong Jin.  FEDformer: Frequency enhanced decomposed transformer for long-term series forecasting.  In _Proceedings of the 39th International Conference on Machine Learning_ , volume 162 of _Proceedings of Machine Learning Research_ , pp. 27268–27286, 2022.
  * Zhu & Shasha (2002) Yunyue Zhu and Dennis Shasha.  Statstream: Statistical monitoring of thousands of data streams in real time.  In _Proceedings of the International Conference on Very Large Databases_ , pp. 358–369. Morgan Kaufmann, 2002.

##  Appendix A Detailed Experiment Setup

#### Dataset

In our experiment, we utilize six datasets – Synthetic, ECG5000, and Traffic
dataset for the simple model (i.e., Sequence-to-Sequence Gated Recurrent Unit)
and ETTh2, ETTm2, and Electricity for the state-of-the-art model (i.e.,
Informer and N-Beats). For each dataset, we describe some metadata of them and
experimental setting, including the input length $n$ and prediction window
$L$.

Synthetic: As Le Guen & Thome (2019) describe, the Synthetic dataset is an
artificial dataset for measuring model performance on sudden changes (step
functions) with an input signal composed of two peaks. The amplitude and
temporal position of the two peaks are randomly selected. Then the selected
position and amplitude of the step are determined by a peak position and
amplitude. We use 500 time-series for training, 500 for validation and 500 for
testing. For the Synthetic dataset, we set input length as $n=20$ and
prediction window as $L=40$. The generation code is provided in DILATE
Github111https://github.com/vincent-leguen/DILATE.

ECG5000: This dataset is originally a 20-hour long ECG (Electrocardiogram),
downloaded from Physionet222https://physionet.org/ and archived in UCR Time
Series Classification Archive (Dau et al., 2019). The data is split by each
heartbeat and processed to be in equal lengths (140). In the training, we use
500 for training, 500 for validation, and 4000 for testing. We take first
$n=84$ steps as input and predict last $L=56$ steps.

Traffic: Traffic dataset is a collection of 48 months (2015-2016) hourly road
occupancy rate (between 0 to 1) data from the California Department of
Transportation333http://pems.dot.ca.gov. As Le Guen & Thome (2019) do, we
utilize univariate series of the first sensor, a total of 17544 data points.
We set our problem as forecasting $L=24$ future occupancy rates with $n=168$
historical data (past week). We use 60% of the data for training, 20% for
validation, and the rest for evaluation.

ETT: The ETT (Electricity Transformer Temperature) dataset, published by Zhou
et al. (2021), is 2-year data collected from two separated counties in China,
including ETTh2 and ETTm2 datasets. Each data point has a target value of “oil
temperature” and other 6 power load features. ETTh2 and ETTm2 datasets have
1-hour and 15-minute intervals, respectively. As Zhou et al. (2021) do, we
split them into 12/4/4 months for the training/validation/testing. Detailed
settings, such as the input and output length and hyperparameter setting, are
based on the information at Informer
Github444https://github.com/zhouhaoyi/Informer2020.

ECL: The ECL (Electricity Consuming Load) is a dataset recorded in kWh every
15-minutes from 2012 to 2014, for 321 clients. In our experiment, we split
them into 15/3/4 months for the train/validation/test, as Zhou et al. (2021)
do. Note that we use the same hyperparameter settings in the ETTh2 dataset.

#### Deep Learning Model Architectures

We perform experiments with three different model architectures, including
Sequence-to-Sequence GRU, Informer, and N-Beats. To induce models to predict
future time-series in a timely manner, we set $\alpha=0.99$ and $\gamma=0.5$
for . Other training metrics, including MSE and DILATE, are used as described
in their original papers. All models are trained with Early Stopping and ADAM
optimizer.

Sequence-to-Sequence GRU To evaluate in simple model, we utilize one layer
Sequence-to-Sequence GRU model. For the training of the GRU model, we set
learning rate of $1e-3$, hidden size of 128, trained by maximum 1000 epochs
with Early Stopping and ADAM optimizer.

Informer When we train Informer with ETTh2, ETTm2, and ECL dataset, we utilize
the official code and hyperparameter setting. In the case of ECL dataset, as
author answered in their official code${}^{\ref{footnote:informer}}$, we
utilize same hyperparameter and dataset splitting criteria as ETTh2 dataset.

N-Beats For N-Beats, we utilize two generic blocks with the hidden size of
128. Additionally, we set the learning rate as $1e-3$ for all three datasets.

Autoformer For Autoformer555https://github.com/thuml/Autoformer, we use the
official code and hyperparameter setting. For the ETTh2 dataset, we utilize
hyperparameter settings described in the official code of
FEDFormer666https://github.com/MAZiqing/FEDformer.

##  Appendix B Additional Evaluations

###  B.1 Detailed Experiment Analysis

At first, we observe that the model optimized with outperforms the same model
optimized with other objective functions in both short- and long-term
forecasting tasks. An interesting point in the results is the large increased
errors of TDI and DTW with long-term forecasting. For example, TDI of Informer
with DILATE shows dramatically increased error with ECL dataset, as the
forecasting window increases, while LCSS does not produce such large increased
error. We attribute this to the weakness of DTW-based loss functions, which
have a weakness due to high sensitiveness on noise. In contrast, does not show
such large performance drop and even achieves better performance in the long-
term forecasting (e.g., Table 2, ETTh2). Additionally, we can find that
Informer with on ECL data and N-Beats with on all three datasets show
significant improvements. It indicates that success to model shape, but other
metrics could not. We provide additional qualitative results below.

Next, we present qualitative analysis of the results. Fig. 3 shows how the
model with different training metrics forecast with different datasets. From
the figure, we have noticed that allows the model to generate more robust,
shape-aware forecasting, regardless of the amplitude shifting, phase shifting,
and uniform amplification. For example, in the case of N-Beats (Fig. 3 (b)
bottom), generate forecasting results, which are more robust, shape-aware
prediction compared to other metrics. We also see the strength in the Informer
case (Fig. 3 (b), top). Even when the model has not enough ability to capture
shape, tries to retrieve the shape. We provide additional qualitative results
with visualization below. When the model have enough ability to capture shape
(i.e., except ETTh2, Informer of $T^{\prime}\in[192,336,720]$), shown its
noise-robust, smooth forecasting with correctly modeled temporal dynamics. In
the most of N-Beats results and some of Informer results, reveals that these
models have enough ability to capture the temporal dynamics with proper loss
function. In summary, proves that it is model-agnostic, noise-robust, and able
to capture the shape.

Figure 3: Qualitative results with simple sequence-to-sequence GRU model (a)
and state-of-the-art model (b).

###  B.2 Additional Qualitative Examples

Figure 4: Qualitative results with simple sequence-to-sequence GRU model
Figure 5: Qualitative results with ETTh2 in short-term forecasting Figure 6:
Qualitative results with ETTh2 in long-term forecasting Figure 7: Qualitative
results with ETTm2 in short-term forecasting Figure 8: Qualitative results
with ETTm2 in long-term forecasting Figure 9: Qualitative results with ECL in
short-term forecasting Figure 10: Qualitative results with ECL in long-term
forecasting

###  B.3 Ablation Study

To evaluate the effect of the $alpha$ and measure the effect of each loss
function, we conduct a set of experiment with ETTh2 dataset and N-Beats on the
long-term forecasting problem. As we can see in the Fig. 11, the model tends
to predict phase- and amplification-free forecasting when $\alpha$ decreases.
This results indicate our motivation, “$\mathcal{L}_{a.shift}$ will return the
forecasting results with same standard deviation with timely manner but
without consideration of proper average value,” and “$\mathcal{L}_{phase}$
will catch the dominant frequency of the target time-series, but it hardly
matches its amplitude and phase.”

Furthermore, in the top of Fig. 12, we can observe three things: (1) if we
utilize $\mathcal{L}_{a.shift}$ only, as we intended, it have different
average (-1.19 vs. 0.11) but relatively similar standard deviation (0.408 vs.
0.299); (2) In the case of $\mathcal{L}_{phase}$ only, they can capture
dominant frequency and produce relatively less-noisy forecasting; (3)
$\mathcal{L}_{amp}$ have relatively similar average value (-1.195 vs. -0.319),
but it has far different standard deviation (0.408 vs. 8.592). In contrast,
forecasting results of the model trained with MSE is very noisy and hard to
interpret (Fig. 12, bottom). Note that we normalized the results in Fig. 12
because of the scale issue.

Table 5: Ablation study on with ETTh2, $L=672$, and N-Beats

Methods | N-Beats
---|---
Eval | MSE | DTW | TDI | LCSS
$\alpha=0.99$ | 0.2388 ± 0.0324 | 13.9859 ± 0.5772 | 18.5824 ± 3.1545 | 0.7063 ± 0.0083
$\alpha=0.9$ | 0.5767 ± 0.8339 | 17.4543 ± 0.4896 | 28.4352 ± 6.9420 | 0.6384 ± 0.0106
$\alpha=0.5$ | 1.9377 ± 0.0222 | 16.2181 ± 1.2729 | 19.8596 ± 2.6131 | 0.6678 ± 0.0197
$\alpha=0.2$ | 1.8361 ± 0.0652 | 16.7026 ± 0.8879 | 16.4266 ± 2.2085 | 0.6726 ± 0.0137
Amplitude only | 1.5123 ± 0.0324 | 16.8485 ± 2.8237 | 29.6971 ± 5.2713 | 0.5740 ± 0.0175
Phase only | 1.8453 ± 0.0470 | 13.7197 ± 0.1196 | 10.2519 ± 1.5417 | 0.6919 ± 0.0043
Amplification only | 674.9155 ± 505.1760 | 14.7763 ± 0.6338 | 30.8568 ± 5.1776 | 0.6120 ± 0.0182

Figure 11: Ablation study result visualization with different $\alpha$ on
ETTm2 dataset Figure 12: Ablation study result visualization of three proposed
loss function on ETTm2 dataset