# Training Robust Deep Models for Time-Series Domain:
Novel Algorithms and Theoretical Analysis

Taha Belkhouja, Yan Yan, Janardhan Rao Doppa

###### Abstract

Despite the success of deep neural networks (DNNs) for real-world applications
over time-series data such as mobile health, little is known about how to
train robust DNNs for time-series domain due to its unique characteristics
compared to images and text data. In this paper, we fill this gap by proposing
a novel algorithmic framework referred as RObust Training for Time-Series (RO-
TS) to create robust deep models for time-series classification tasks.
Specifically, we formulate a min-max optimization problem over the model
parameters by explicitly reasoning about the robustness criteria in terms of
additive perturbations to time-series inputs measured by the global alignment
kernel (GAK) based distance. We also show the generality and advantages of our
formulation using the summation structure over time-series alignments by
relating both GAK and dynamic time warping (DTW). This problem is an instance
of a family of compositional min-max optimization problems, which are
challenging and open with unclear theoretical guarantee. We propose a
principled stochastic compositional alternating gradient descent ascent
(SCAGDA) algorithm for this family of optimization problems. Unlike
traditional methods for time-series that require approximate computation of
distance measures, SCAGDA approximates the GAK based distance on-the-fly using
a moving average approach. We theoretically analyze the convergence rate of
SCAGDA and provide strong theoretical support for the estimation of GAK based
distance. Our experiments on real-world benchmarks demonstrate that RO-TS
creates more robust deep models when compared to adversarial training using
prior methods that rely on data augmentation or new definitions of loss
functions. We also demonstrate the importance of GAK for time-series data over
the Euclidean distance. The source code of RO-TS algorithms is available at
https://github.com/tahabelkhouja/Robust-Training-for-Time-Series

## Introduction

Predictive analytics over time-series data enables many important real-world
applications including mobile health, smart grid management, smart home
automation, and finance. In spite of the success of deep neural networks
(DNNs) (Wang, Yan, and Oates 2017), very little is known about the robustness
of DNNs for time-series data. Recent work on adversarial examples for image
(Kolter and Madry 2018) and text data (Wang, Singh, and Li 2019) exposed the
brittleness of DNNs and motivated methods to improve their robustness.
Therefore, training robust deep models is a necessary requirement when we
deploy DNNs over time-series data in critical and high-stakes applications
(e.g., mobile health (Belkhouja and Doppa 2020)). Real-world deployment of
deep models for time-series data pose several robustness challenges. First,
labeled training set for time-series domain tend to be smaller when compared
to image and text domains. Consequently, the learned DNNs may not perform well
on unseen data drawn from the same distribution. Second, disturbances due to
noisy sensor observations or different sampling frequencies can potentially
result in poor performance of DNN models. Third, adversarial attacks for
malicious purposes to break DNN models can pose security threats.

Prior work on robustness for image data considers training DNNs to be robust
to adversarial attacks in a $L_{p}$ ball and input perturbations. Algorithms
to improve robustness of DNNs fall into two broad categories. First,
adversarial training using data augmentation (e.g., adversarial examples or
input perturbations). Second, optimizing an explicit loss function for
robustness criterion (e.g., similar input images should produce similar DNN
outputs). Almost all prior methods are designed for images and are based on
the $L_{p}$ norm distance. Since time-series data has unique characteristics
(e.g., sparse peaks, fast oscillations), $L_{p}$ distance can rarely capture
the true similarity between time-series pairs and prior methods are likely to
fail on time-series data, as demonstrated by results in Fig. 1. The main
question of this paper is: what are good methodologies to train robust DNNs
for time-series domain by handling its unique challenges?

To answer this question, we propose a novel and principled framework referred
as RObust Training for Time-Series (RO-TS) to create robust DNNs for time-
series data. We employ additive noise variables to simulate perturbations
within a small neighborhood around each training example. We incorporate these
additive noise variables to formulate a min-max optimization problem to reason
about the robustness criteria in terms of disturbances to time-series inputs
by minimizing the worst-case risk. To capture the special characteristics of
time-series signals, we employ the global alighnment kernel (GAK) based
distance (Cuturi et al. 2007) to define neighborhood regions for training
examples. We show the generality and advantage of our formulation using the
summation structure over time-series alignments by relating both GAK and
dynamic time warping (DTW) (Berndt and Clifford 1994).

Unfortunately, the min-max optimization problem with GAK based distance is
challenging and may fall in the family of compositional min-max problems due
to lack of theoretical guarantees. To efficiently solve this general family of
optimization problems, we develop a principled stochastic compositional
alternating gradient descent ascent (SCAGDA) algorithm by carefully leveraging
the underlying structure of this problem. Another key computational challenge
is that time-series distance measures including the GAK based distance involve
going through all possible alignments between pairs of time-series inputs,
which is expensive, e.g., $O(T^{2})$ for GAK where $T$ is the length of time-
series signal. As a consequence, the computational cost grows significantly
for iterative optimization algorithms where we need to repeatedly compute the
distance between time-series signals. SCAGDA randomly samples a constant
number of alignments at each iteration to approximate the GAK based distance
on-the-fly using a moving average approach and reduces the computational-
complexity to $O(T)$.

Our theoretical analysis shows that SCAGDA achieves an $\epsilon$-primal gap
in $O(1/\epsilon^{2})$ iterations for a family of nonconvex-nonconcave
compositional min-max optimization problems. To the best of our knowledge,
this is the first convergence rate established for nonconvex-nonconcave min-
max problems with a compositional structure. We also prove that SCAGDA
approximates GAK based distance during the execution of algorithm: the
approximation error converges as the primal gap converges. This result
provides strong theoretical support for our algorithm design that leverages
the structure of the RO-TS formulation. Our experiments on real-world datasets
show that DNNs learned using RO-TS are more robust than prior methods; and GAK
based distance has a more appropriate bias for time-series than $L_{2}$.

Contributions. The key contribution of this paper is the development and
evaluation of the RO-TS algorithmic framework to train robust deep models for
time-series domain.

  * •

Novel SCAGDA algorithm for a family of nonconvex-nonconcave compositional min-
max problems that covers RO-TS as a special case with GAK based distance.
SCAGDA approximates the GAK distance by randomly sampling a constant number of
time-series alignments.

  * •

Theoretical analysis of SCAGDA shows that the iteration complexity to achieve
$\epsilon$-primal gap and the $\epsilon$-approximation error of GAK distance
is $O(1/\epsilon^{2})$.

  * •

Comprehensive experimental evaluation of RO-TS on diverse real-world benchmark
datasets and comparison with state-of-the-art baselines. The source code of
RO-TS algorithms is available at https://github.com/tahabelkhouja/Robust-
Training-for-Time-Series

## Problem Setup and Formulation

We consider the problem of learning robust DNN classifiers over time-series
data. We are given a training set of $n$ input-output pairs
$\left\\{\left(x_{i},y_{i}\right)\right\\}_{i=1}^{n}$. Each input
$x_{i}\in\mathcal{X}$ is a time-series signal, where
$\mathcal{X}\subseteq\mathbb{R}^{C\times T}$ with $C$ denoting the number of
channels and $T$ being the window-size of the signal; and
$y_{i}\in\mathcal{Y}$ is the associated ground-truth label, where
$\mathcal{Y}\in\\{1,\cdots,\mathcal{C}\\}$ is a set of $\mathcal{C}$ discrete
class labels. For example, in a health monitoring application using
physiological sensors for patients diagnosed with cardiac arrhythmia, we use
the measurements from wearable devices to predict the likelihood of a cardiac
failure. Traditional empirical risk minimization learns a DNN classifier
$f:\mathcal{X}\times\Theta\rightarrow\mathcal{Y}$ with weights $w\in\Theta$
that maps time-series inputs to classification labels for a hypothesis space
$\Theta$ and a loss function $\ell$:

$\displaystyle\min_{w\in\Theta}\frac{1}{n}\sum_{i=1}^{n}\ell(f(x_{i},w),y_{i}).$

Training for robustness. We would like the learned classifier $f(x,w)$ to be
robust to disturbances in time-series inputs due to noisy observations or
adversarial attacks. For example, a failure in the prediction task for the
above health monitoring application due to such disturbances can cause injury
to the patient without the system notifying the needed assistance. Therefore,
we want the trained DNN classifier to be invariant to such disturbances.
Mathematically, for an appropriate distance function $d(x,x^{\prime})$ over
time-series inputs $x$ and $x^{\prime}$, we want the classifier $f$ to predict
the same classification label as $x$ for all inputs $x^{\prime}$ such that
$d(x,x^{\prime})<\varepsilon$, where $\varepsilon$ stands for the bound on
allowed disturbance to input $x$. This goal can be achieved by reasoning about
the worst-case empirical risk over possible perturbations
$a_{i}\in\mathbb{R}^{C\times T}$ of $x_{i}$ such that
$d(x_{i},x_{i}+a_{i})\leq\varepsilon$. The resulting min-max optimization
problem is given below.

$\displaystyle\min_{w\in\Theta}~{}$
$\displaystyle~{}\frac{1}{n}\sum_{i=1}^{n}\max_{a_{i}}\ell(f(x_{i}+a_{i},w),y_{i})$
s.t. $\displaystyle~{}d(x_{i},x_{i}+a_{i})\leq\varepsilon$ (1)

In practice, instead of solving the above hard constrained problem, one can
solve an equivalent soft constrained problem using regularization as follows

$\displaystyle\min_{w\in\Theta}\max_{a_{i}}\frac{1}{n}\sum_{i=1}^{n}\ell(f(x_{i}+a_{i},w),y_{i})-\lambda
d(x_{i},x_{i}+a_{i})$ (2)

There is a natural interpretation of this optimization problem. The inner
maximization problem serves the role of an attacker whose goal is to find
adversarial examples that achieves the largest loss. The outer minimization
problem serves the role of a defender whose goal is to find the parameters of
the deep model $w$ by minimizing the adversarial loss from the inner attack
problem. This formulation is applicable to all types of data by selecting an
appropriate distance function $d$. For example, $L_{p}$-norm distance is
usually used in the image domain (Kolter and Madry 2018).

Typical stochastic approaches to solving the above adversarial training
problem include alternating stochastic optimization (Junchi Y. 2020) and
stochastic gradient descent ascent (GDA) (Lin, Jin, and Jordan 2020; Yan et
al. 2020). The alternating method first fixes $w$ and solves the inner
maximization approximately to get each $a_{i}$ (e.g., using stochastic
gradient descent). Next, $a_{i}$ is fixed and the outer minimization is solved
over $w$. These two steps are performed alternatively until convergence. The
GDA method computes the gradient of $w$ and $a_{i}$ simultaneously at each
iteration, and then use these gradients to update $w$ and $a_{i}$. Both
methods require the ability to compute the unbiased estimation of the
gradients w.r.t. $w$ and $a_{i}$. When $d$ is decomposable, e.g.,
$L_{p}$-norm, then its stochastic gradient can be easily computed. However, in
time series domain, commonly used distance measures may not be decomposable,
so its stochastic gradients are not accessible. Consequently, one has to
calculate the exact gradient of $d(x_{i},x_{i}+a_{i})$ w.r.t, $a_{i}$. We will
investigate this key challenge in Section SCAGDA Optimization Algorithm.

We summarize in Table 1 the main mathematical notations used in this paper.

Variable | Definition
---|---
$\mathcal{X}$ | Input space of time-series data
$\mathcal{Y}$ | Discrete set of class labels
$f$ | Deep neural network classifier
$w\in\Theta$ | Weights of classifier $f$
$l(\cdot,\cdot,\cdot)$ | Surrogate loss function
$k_{\mathrm{GAK}}(\cdot,\cdot)$ | Global alignment kernel
$\pi$ | Alignment between two time-series
| defined as a pair $(\pi_{1},\pi_{2})$
$\mathcal{A}$ | Set of all possible alignments
$d_{\pi}(\cdot,\cdot)$ | Distance function according
| to an alignment $\pi$
Table 1: Mathematical notations used in this paper.

## RO-TS Algorithmic Framework

In this section, we describe the technical details of our proposed RO-TS
framework to train robust DNN classifiers for time-series domain. First, we
instantiate the min-max formulation with GAK based distance as it
appropriately captures the similarity between time-series signals. Second, we
provide an efficient algorithm to solve the GAK-based formulation to learn
parameters of DNN classifers.

### Distance Measure for Time-Series

Unlike images and text, time-series data exhibits unique characteristics such
as sparse peaks, fast oscillations, and frequency/time shifting which are
important for pattern matching and data classification. Hence, measures such
as Euclidean distance that do not account for these characteristics usually
fail in recognizing the similarities between time-series signals. To address
this challenge, elastic measures have been introduced for pattern-matching
tasks for time-series domain (Cuturi et al. 2007), where one time-step of a
signal can be associated with many time-steps of another signal to compute the
similarity measure.

Time-series alignment. Given two time series $x$=$(x_{1},\cdots,x_{T_{1}})$
and $x^{\prime}$=$(x^{\prime}_{1},\cdots,x^{\prime}_{T_{2}})$ for
$T_{1},T_{2}\in\mathbb{N}_{+}$, the alignment $\pi$ = $(\pi_{1},\pi_{2})$ is
defined as a pair of increasing integral vectors of length $r\leq
T_{1}+T_{2}-1$ such that $1$ = $\pi_{1}(1)\leq\cdots\leq\pi_{1}(r)$ = $T_{1}$
and $1$ = $\pi_{2}(1)\leq\cdots\leq\pi_{2}(r)$ = $T_{2}$ with unitary
increments and without simultaneous repetitions, which presents the
coordinates of $x$ and $x^{\prime}$. This alignment defines the one-to-many
alignment between $x$ and $x^{\prime}$ to measure their similarity. Using a
candidate alignment $\pi$, we can compute their similarity as follows:

$\displaystyle
d_{\pi}(x,x^{\prime})=\sum_{i=1}^{|\pi|}\mathrm{dist}(x_{\pi_{1}(i)},x^{\prime}_{\pi_{2}(i)})$
(3)

where $|\pi|$=$r$ denotes the length of alignment and
$\mathrm{dist}(\cdot,\cdot)$ in the above equation is the Minkowski distance:

$\mathrm{dist}(x_{\pi_{1}(i)},x^{\prime}_{\pi_{2}(i)})=\|x_{\pi_{1}(i)}-x^{\prime}_{\pi_{2}(i)}\|_{p},~{}p\in\\{1,\cdots,\infty\\}$

Global alignment kernel (GAK) based distance. The concept of alignment allows
us to take into consideration the intrinsic properties of time-series signals,
such as frequency shifts, to compute their similarity. There are some well-
known approaches to define distance metrics using time-series alignment. For
example, dynamic time warping (DTW) (Berndt and Clifford 1994) selects the
alignment with the minimum distance:

$\displaystyle
D_{\mathrm{DTW}}(x,x^{\prime})=\min_{\pi\in\mathcal{A}}d_{\pi}(x,x^{\prime}),$

where $\mathcal{A}$ denotes the set of all possible alignments.

While DTW only takes into account one candidate alignment, global alignment
kernel (GAK) (Cuturi et al. 2007) takes all possible alignments into
consideration:

$\displaystyle
k_{\mathrm{GAK}}(x,x^{\prime})=\sum_{\pi\in\mathcal{A}}\exp\Big{(}-\frac{d_{\pi}(x,x^{\prime})}{\nu}\Big{)}$
(4)

where $\nu$ is a hyper-parameter and $d_{\pi}(\cdot,\cdot)$ is defined in
Equation (3). In practice, to handle the diagonally dominance issue (Cuturi et
al. 2007; Wu et al. 2018; Cuturi 2011),
$D_{\mathrm{GAK}}:=-\nu\log(k_{\mathrm{GAK}})$ is typically used as a distance
measure for a pair of time-series signals. GAK enjoys several advantages over
DTW (Cuturi 2011): (i) differentiable, (ii) positive definite, (iii) coherent
measure over all possible alignments. Therefore, $k_{\mathrm{GAK}}$ (or
$D_{\mathrm{GAK}}$) is a better fit to train robust DNNs for the time-series
domain.

On the other hand, GAK can also be a more general measure than DTW due to its
summation structure, as $\lim_{\nu\rightarrow
0}D_{\mathrm{GAK}}(x,x^{\prime})=D_{\mathrm{DTW}}(x,x^{\prime})$, i.e.,
arbitrarily close to DTW by changing $\nu$. The following proposition shows
the tight approximation of the soft minimum of GAK to the hard minimum of DTW
(proof and details in Appendix).

######  Proposition 1.

For a time-series pair $(x,x^{\prime})$, we have:

$\displaystyle 0\leq
D_{DTW}(x,x^{\prime})-D_{\mathrm{GAK}}(x,x^{\prime})\leq\nu\log(|\mathcal{A}|).$

As shown, $D_{\mathrm{GAK}}$ converges to $D_{\mathrm{DTW}}$ in
$\nu\log(|\mathcal{A}|)$ as $\nu$ decreases. Due to the above advantages and
the approximation ability of $D_{\mathrm{GAK}}$ to $D_{\mathrm{DTW}}$, we
consider the more general $k_{\mathrm{GAK}}$ and $D_{\mathrm{GAK}}$ in our RO-
TS method.

### SCAGDA Optimization Algorithm

By plugging $k_{\mathrm{GAK}}$ from Equation (4) to replace $d$ into the min-
max formulation in Equation (2), we reach the following objective function of
our RO-TS framework:

$\displaystyle\min_{w\in\Theta}\max_{a_{i}}$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\ell(f(x_{i}+a_{i},w),y_{i})$
$\displaystyle+\lambda\overbrace{\log\big{(}k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i})\big{)}}^{=d(x_{i},x_{i}+a_{i})}$
(5)

where $\nu$ outside $\log$ in $D_{\mathrm{GAK}}$ can be merged into $\lambda$.
The above problem is decomposable over individual training examples (i.e.,
index $i$), so we can compute stochastic gradients by randomly sampling a
batch of data and employ stochastic gradient descent ascent (SGDA) (Lin, Jin,
and Jordan 2020; Yan et al. 2020), a family of stochastic algorithms for
solving min-max problems.

Key challenge. The second term $\log(k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i}))$ has
a compositional structure due to the outer $\log$ function. By chain rule, its
gradient w.r.t. the dual variable $a_{i}$ is

$\displaystyle\nabla_{a_{i}}\log(k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i}))=\frac{\nabla_{a_{i}}k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i})}{k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i})}$

where one has to go through all possible alignments to compute
$k_{\mathrm{GAK}}$ and $\nabla_{a_{i}}k_{\mathrm{GAK}}$ (see Equation (4)) and
there is no unbiased estimation (i.e., stochastic gradients) for it.

Consequently, at each iteration, SGDA has to compute the exact value of
$k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i})$ and
$\nabla_{a_{i}}k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i})$ according to the chain
rule, which leads to an additional time-complexity of $O(CT^{2})$ per SGDA
iteration, where $C$ and $T$ denote the number of channels and window-size
respectively. This computational bottleneck will lead to extremely slow
training algorithm when $C$ and/or $T$ is large, which is the case in many
real-world applications.

One candidate approach to alleviate the computational challenge due to
$\log(k_{\mathrm{GAK}})$ part of the objective is to make use of the inner
summation structure of $k_{\mathrm{GAK}}$. Since $k_{\mathrm{GAK}}$ involves a
summation over all alignments, as shown in (4), we can use only a subset of
alignments for estimating the full summation. This procedure will give an
unbiased estimation of $k_{\mathrm{GAK}}$, but the outer logarithmic function
makes it a biased estimation for
$\nabla_{a_{i}}\log(k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i}))$. However, such
biased estimation violates the assumption in SGDA studies, so their
theoretical analysis cannot hold.

There is another line of research investigating stochastic compositional
gradient methods for minimization problems with compositional structure (Wang,
Fang, and Liu 2017; Chen, Sun, and Yin 2020). However, min-max optimization
with compositional structure, including our case shown in Equation (SCAGDA
Optimization Algorithm), is not studied yet. It is unclear whether these
techniques and analysis hold for min-max problems.

Algorithm 1 SCAGDA (Stochastic Compositional Alternating Gradient Descent
Ascent)

1: Initialize $w_{0},a_{i}^{0}$ for $i=1,\cdots,n$ and $\omega_{i}^{0}=0$ for
$i=1,\cdots,n$, step sizes $\\{\eta_{k}\\}_{k=1}^{K}$ and
$\\{\gamma_{k}\\}_{k=1}^{K}$.

2: for $k=0,\cdots,K-1$ do

3: Randomly sample an index $i_{1}$ to compute stochastic gradient
$\nabla_{w}f_{i_{1}}(w_{k},a_{i_{1}}^{k})$

4: Set: $w_{k+1}=w_{k}-\eta_{k}\nabla_{w}f_{i}(w_{k},a_{i_{1}}^{k})$

5: Randomly sample an index $i_{2}$ to compute stochastic gradient
$\nabla_{a}f_{i_{2}}(w_{k+1},a_{i_{2}}^{k})$

6: Randomly sample two independent indices $j_{1},j_{2}$ of $h_{i_{2}}$ to
compute $h_{i_{2},j_{1}}(a_{i_{2}}^{k})$ and $\nabla
h_{i_{2},j_{2}}(a_{i_{2}}^{k})$

7: Set:  $\omega_{i}^{k+1}=\left\\{\begin{array}[]{ll}\omega_{i}^{k}&\text{
for }i\neq i_{2}\\\ (1-\beta)\omega^{k}_{i_{2}}+\beta
h_{i_{2},j_{1}}(a_{i_{2}}^{k})&\text{ for }i=i_{2}\end{array}\right.$

8: Set:
$a_{i_{2}}^{k+1}=a_{i_{2}}^{k}+\gamma_{k}(\nabla_{a}f_{i_{2}}(w_{k+1},a_{i_{2}}^{k})-\nabla
h_{j2}(a_{i_{2}}^{k})^{\top}\nabla g(\omega^{k+1}_{i_{2}}))$

9: end for

10: return final solution $w_{K}$

SCAGDA algorithm. We propose a novel stochastic compositional alternating
gradient descent ascent (SCAGDA) algorithm to solve a family of nonconvex-
nonconcave min-max compositional problems, which include RO-TS (Equation
(SCAGDA Optimization Algorithm)) as a special case. We summarize SCAGDA in
Algorithm 1. Specifically, we consider solving the following family of
problems:

$\displaystyle\min_{w}\max_{a_{i}}\frac{1}{n}\sum_{i=1}^{n}\phi_{i}(w,a_{i})$
(6)

where
$\phi_{i}(w,a_{i}):=f_{i}(w,a_{i})-g(\frac{1}{m}\sum_{j=1}^{m}h_{i,j}(a_{i}))$.

Mapping Problem (6) to RO-TS (SCAGDA Optimization Algorithm). As mentioned
above, RO-TS for time-series in Equation (SCAGDA Optimization Algorithm) is a
special case of Problem (6) as shown below. The variables
$\phi_{i}(w,a_{i}),f_{i},g,h_{i,j}$ in Problem (6) can be instantiated by the
following mappings:

  * •

$f_{i}$ in $\phi_{i}$ of (6) $\Rightarrow$ the loss $\ell$ on the $i$-th data
in (SCAGDA Optimization Algorithm)

  * •

$-g(\cdot)$ in $\phi_{i}$ of (6) $\Rightarrow$ $\lambda\log(\cdot)$ in (SCAGDA
Optimization Algorithm)

  * •

$\frac{1}{m}\sum_{j=1}^{m}h_{i,j}(a_{i})$ in $\phi_{i}$ of (6) $\Rightarrow$
$k_{\mathrm{GAK}}(x_{i},x_{i}+a_{i})=\sum_{\pi\in\mathcal{A}}\exp(-d_{\pi}(x_{i},x_{i}+a_{i})/\nu)$
in (SCAGDA Optimization Algorithm), where $m$ and $j$ corresponds to the total
number of alignment paths $|\mathcal{A}|$ and the index of alignment path,
respectively. Note that the summation form of $k_{\mathrm{GAK}}$ can be easily
converted to an average form due to $\log(x)=\log(x/m)+\log(m)$.

Algorithmic analysis of SCAGDA. To introduce and analyze Algorithm 1 for
solving Problem (6), we first introduce some notations. Denote
$P(w):=\max_{a_{i}}\frac{1}{n}\sum_{i=1}^{n}\phi_{i}(w,a_{i})$ as the primal
function of the above min-max optimization problem, where we are interested in
analyzing the convergence of the primal gap after the $K$-th iteration:

$P(w_{K})-\min_{w}P(w).$

Let $a:=(a_{1},a_{2},...,a_{n})\in\mathbb{R}^{C\times T\times n}$ be the
concatenation of $a_{i}$ for $i=1,\cdots,n$. We also use the following
notations to improve the technical exposition and ease of readability.

$\displaystyle\phi(w,a):=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\phi_{i}(w,a_{i}),$ $\displaystyle
h_{i}(a_{i}):=$ $\displaystyle\frac{1}{m}\sum_{j=1}^{m}h_{i,j}(a_{i}),$
$\displaystyle h(a):=$
$\displaystyle\frac{1}{n}(h_{1}(a_{1}),\cdots,h_{n}(a_{n})),$

where the last term $h(a)$ is the concatenation of all $h_{i}$ for
$i=1,\cdots,n$. In Appendix , we provide the details of how (6) is
specifically viewed as a stochastic problem.

As mentioned above while discussing the key challenge of the compositional
structure in RO-TS (SCAGDA Optimization Algorithm), conventional SGDA methods
for Problem (6) require us to compute the full gradient of the compositional
part $g(h_{i}(a_{i}))$, i.e., $\nabla h_{i}(a_{i})^{\top}\nabla
g(h_{i}(a_{i}))$, which involves all alignments in the case of RO-TS.

In contrast, SCAGDA only samples a constant number of $h_{i,j}(a_{i})$ over
$j$ (i.e., over a subset of alignments for RO-TS) and $\nabla h_{i,j}(a_{i})$
(Line 6). Subsequently, SCAGDA employs a simple iterative moving average (MA)
approach to accumulate $h_{i,j}(a_{i})$ into $\omega^{k+1}_{i}$ at iteration
$k$ for estimating $h_{i}(a_{i})$ (Line 7). The key idea behind moving average
method is to control the variance of the estimation for $h_{i}(a_{i}^{k+1})$
using a weighted average from the previous estimate $h_{i}(a_{i}^{k})$. Even
though $\nabla g(\omega^{k+1}_{i})$ is a biased estimation of $\nabla
g(h_{i}(a_{i}^{k}))$, we can still use smoothness condition (introduced in
Section Theoretical Analysis later) and bound the approximation error
$\mathbb{E}[\|\omega_{i}^{k}-h_{i}(a_{i}^{k-1})\|^{2}]$ where
$\omega_{k}:=(\omega_{1}^{k},...,\omega_{n}^{k})$ is the concatenation of all
$\\{\omega_{i}^{k}\\}_{i=1}^{n}$ at iteration $k$, as shown in Theorem 2.

Therefore, instantiation of SCAGDA for RO-TS does not require us to perform
computation over all alignments contained in $h_{i}(a_{i})$ for each time-
series training sample, which leads to a more efficient algorithm with high
scalability on large datasets. As shown in Line 3 and 5, SCAGDA updates the
primal variable $w$ and dual variable $a$ in an alternating scheme, which
means that $w_{k+1}$ is updated based on $a_{k}$, while $a_{k+1}$ is updated
based on $w_{k+1}$. This is different from SGDA, which updates $a_{k}$ based
on $w_{k}$ instead. We instantiate SCAGDA for the proposed RO-TS framework as
shown in Algorithm 2. The primal variable update is provided in Line 6, and
the dual variable update is provided in Line 11. In particular, Line 8 and 8
correspond to the moving average step for estimating $k_{\mathrm{GAK}}$ using
randomly sampled alignment subset $\widehat{\mathcal{A}}_{i}^{k}$ for the
$i$-th time-series training example.

Algorithm 2 RO-TS Instantiation of SCAGDA

Input: A training set
$\\{(x,y)\in\mathcal{X}\times\mathcal{Y}\\}^{n_{\text{train}}}$; mini-batch
size $s$, deep neural network $f(w,x,y)$; learning rates $\eta_{k}$ and
$\gamma_{k}$, loss function $l(\cdot)$, distance function $D(\cdot,\cdot)$.
Output: Classifier weights $w\in\Theta$

1: Randomly initialize weights of the DNN classifier: $w_{0}\in\Theta$ //
vector of worst-case perturbations, one for each time-series

2: Initialize $a_{0}=0$ and $\omega_{0}=0$ // Multiple iterations of SCAGDA

3: for $k=0,\cdots,K-1$ do

4: Randomly sample a mini-batch of data samples indexed by $\mathcal{I}_{k}$
s.t. $|\mathcal{I}_{k}|=s$

5: Compute the stochastic gradient w.r.t. $w_{k}$
$\mathcal{G}_{w,k}=\frac{1}{s}\sum_{i\in\mathcal{I}_{k}}\nabla_{w}l(f(w_{k},x_{i}+a_{i}^{k},y_{i})$

6: Perform stochastic gradient descent on $w_{k}$
$w_{k+1}=w_{k}-\eta_{k}\mathcal{G}_{w,k}$

7: Randomly sample a mini-batch of alignments indexed by
$\widehat{\mathcal{A}}_{i}^{k}$ for each data index $i\in\mathcal{I}_{k}$

8: Moving average for
$i\in\mathcal{I}_{k}$$\omega_{i}^{k+1}=(1-\beta)\omega_{i}^{k}+\beta\sum_{\pi\in\widehat{\mathcal{A}}_{i}^{k}}\exp(-d_{\pi}(x_{i},x_{i}+a_{i}^{k})/\nu)$

9: $\omega_{i}^{k+1}=\omega_{i}^{k}$ for $i\notin\mathcal{I}_{k}$.

10: Compute
$\mathcal{G}_{a,k,i}=\nabla_{a}\ell(f(w_{k+1},x_{i}+a_{i}^{k},y_{i})-\sum_{\pi\in\widehat{\mathcal{A}}_{i}^{k}}\exp(-d_{\pi}(x_{i},x_{i}+a_{i}^{k})/\nu)\cdot\nabla_{a}d_{\pi}(x_{i},x_{i}+a_{i}^{k})\cdot\frac{\lambda}{\omega^{k+1}_{i}\nu}$
for $i\in\mathcal{I}_{k}$

11: Perform stochastic gradient ascent over perturbations
$a_{i}^{k+1}=a_{i}^{k}+\gamma_{k}\mathcal{G}_{a,k,i}$

12: end for

13: return weights of the learned DNN classifier, $w_{K}$

In the next section, we show that our algorithm can converge to primal gap
$P(w_{K})-\min_{w}P(w)\leq\epsilon$ with iteration complexity
$O(1/\epsilon^{2})$, where $\epsilon$ is a pre-defined threshold. To the best
of our knowledge, this is the first optimization algorithm and convergence
analysis for the famaily of compositional min-max optimization problems shown
in (6).

## Theoretical Analysis

In this section, we present novel theoretical convergence analysis for SCAGDA
algorithm. As mentioned in the previous section, for the problem (6), existing
theoretical analysis of SGDA (Lin, Jin, and Jordan 2020; Yan et al. 2020),
stochastic alternating gradient descent ascent (SAGDA) (Junchi Y. 2020)
require us to compute exact gradient of $g(h_{i}(a_{i}))$ at each iteration.
On the other hand, it is unclear if stochastic compositional alternating
gradient algorithms for minimization problems (Wang, Fang, and Liu 2017; Chen,
Sun, and Yin 2020) can handle the complex min-max case.

Summary of results. We answer the following question: can we establish
convergence guarantee of our SCAGDA algorithm for nonconvex-nonconcave
compositional min-max optimization problems?

Theorem 1 proves that SCAGDA shown in Algorithm 1 converges to an
$\epsilon$-primal gap in $O(\frac{1}{\epsilon^{2}})$ iterations. Theorem 2
demonstrates the efficacy of the moving average strategy to approximate GAK
based distance: the approximation error
$\|\omega^{K}_{i}-h_{i}(a_{i}^{K-1}))\|^{2}$ is also bounded by $\epsilon$ in
expectation when the $\epsilon$-primal gap is achieved.

### Main Results

The following commonly used assumptions are used in our analysis. Due to the
space limit, definitions, proofs, and detailed constant dependencies can be
found in Appendix .

######  Assumption 1.

Suppose $\mu,L,C_{g},C_{h},L_{g},L_{h}\geq 0$.
(i) $\phi(w,a)$ satisfies two side $\mu$-PL (Polyak-Lojasiewicz) condition:

$\displaystyle\|\nabla_{w}\phi(w,a)\|^{2}\geq
2\mu(\phi(w,a)-\min_{w^{\prime}}f(w,a)),$
$\displaystyle\|\nabla_{a}\phi(w,a)\|^{2}\geq
2\mu(\max_{a^{\prime}}f(w,a)-\phi(w,a)).$

(ii) $\phi(w,a)$ is $L$-smooth in $w$ for fixed $a$.

(iii) $\phi(w,a)$ is $L$-smooth in $a_{i}$ for fixed $w$.

(iv) $g$ (resp. $h$) is $C_{g}$ (resp. $C_{h}$)-Lipschitz continuous.

(v) $g$ (resp. $h$) is $L_{g}$ (resp. $L_{h}$)-smooth.

(vi) $\exists$ $\sigma>0$ s.t.
$\mathbb{E}[\|\nabla_{w}\phi_{i}(w,a_{i})-\nabla_{w}\phi(w,a)\|^{2}]\leq\sigma^{2}$,

$\mathbb{E}[\|\nabla_{a}f_{i}(w,a_{i})-\nabla_{a}f_{i}(w,a_{i})\|^{2}]\leq\sigma^{2}$,
$\mathbb{E}[\|h_{i,j}(a_{i})-h(a)\|^{2}]\leq\sigma^{2}$, and
$\mathbb{E}[\|\nabla h_{i,j}(a_{i})-\nabla h(a)\|^{2}]\leq\sigma^{2}$

We present our main results for SCAGDA below.

######  Theorem 1.

Suppose Assumption 1 holds. In Algorithm 1, set
$\eta_{k}=\eta=O(1/\epsilon^{2})$, $\gamma_{k}=\gamma=O(1/L^{2})$ and
$\beta=\sqrt{18\mu\eta}$. After running Algorithm 1 for $K$ iterations where
$K=\tilde{O}(1/\epsilon^{2})$ ($\tilde{O}$ hides logarithmic factor), we have

$\displaystyle\mathbb{E}[P(w_{K})-P^{*}]+\frac{1}{8}\mathbb{E}[P(w_{K})-\phi(w_{K},a_{K})]$
$\displaystyle+\Big{(}\frac{4C_{h}^{4}L_{g}^{4}\eta_{K}}{\mu^{5}}\Big{)}^{1/2}\mathbb{E}[\|\omega_{K}-h(a_{K-1})\|^{2}]\leq\epsilon$

Remark 1. The above theorem gives us two critical observations of the behavior
of SCAGDA. (1) After running $K$ iterations of SCAGDA, the primal gap
$P(w_{K+1})-P^{*}$ converges to $\epsilon$ in expectation, since all terms in
the left hand side of the inequality are non-negative. This result shows that
SCAGDA is able to effectively solve the compositional min-max optimization
problem shown in Equation (6). (2) The required iteration complexity of SCAGDA
is $O(1/\epsilon^{2})$. To put this result in perspective, we compare it with
related theoretical results. The rate for nonconvex-nonconcave min-max problem
without compositional structure is shown to be $O(1/\epsilon)$ (Junchi Y.
2020). However, this improvement requires unbiased estimation (or exact value)
of the gradient and computing the exact $g(h_{i}(a_{i}))$ at each iteration.
Our iteration complexity is in the same order of that for (Chen, Sun, and Yin
2020), whose convergence result is $O(1/\epsilon^{2})$ for nonconvex
compositional minimization problems instead of min-max ones. The difference is
that their convergence metric is the average squared norm of gradients, while
ours is for the primal gap. Importantly, this is the first result on
convergence rate for stochastic compositional min-max problems.

######  Theorem 2.

After $K=\tilde{O}(1/\epsilon^{2})$ iterations of Algorithm 1, we have:
$\mathbb{E}[\|\omega^{K}_{i}-h_{i}(a^{K-1}_{i})\|^{2}]\leq O(\epsilon).$

Remark 2. The above result shows that as SCAGDA algorithm is executed, the
approximation error of $\|\omega^{K}_{i}-h_{i}(a^{K-1}_{i})\|^{2}$ converges
to $\epsilon$ in the expectation as it is achieving the $\epsilon$-primal gap.
For the condition numbers, we always have $L\geq\mu$. In practice, we usually
set the accuracy level $\epsilon$ to a very small value, so the condition
$\epsilon\leq O(L^{3}/\mu^{2})$ will generally hold. This result provides
strong theoretical support that if we apply SCAGDA to optimize our RO-TS
problem in (SCAGDA Optimization Algorithm), it is able to approximate
$k_{\mathrm{GAK}}$ on-the-fly, where we only need a constant number of
alignments, rather than all possible alignments for computing
$k_{\mathrm{GAK}}$ in each iteration of SCAGDA. When we have $\epsilon$-primal
gap, we also achieve $\epsilon$-accurate estimation of $k_{\mathrm{GAK}}$ at
the same time.

## Related Work

Prior work on robustness of DNNs is mostly focused on image/text domains; and
can be classified into two categories.

Adversarial training employs augmented data such as adversarial examples
(Kolter and Madry 2018; Wang, Singh, and Li 2019) and input perturbations.
Methods to create adversarial examples include general attacks such as Carlini
& Wagner attack (Carlini and Wagner 2017), boundary attack (Brendel, Rauber,
and Bethge 2018), and universal attacks (Moosavi-Dezfooli et al. 2017). Recent
work regularizes adversarial example generation methods to obey intrinsic
properties of images (Laidlaw and Feizi 2019; Xiao et al. 2018; Hosseini et
al. 2017). There are also specific adversarial methods for NLP domain (Samanta
and Mehta 2017; Gao et al. 2018). There is little to no prior work on
adversarial techniques for time-series domain. Fawaz et al. (Fawaz et al.
2019) employed the standard Fast Gradient Sign method (Kurakin, Goodfellow,
and Bengio 2017) to create adversarial noise for time-series. Network
distillation was also employed to train a student model for creating
adversarial attacks (Fazle, Somshubra, and Houshang 2020). This method is
severely limited: it can generate adversarial examples for only a small number
of target labels and cannot guarantee generation of adversarial example for
every input.

Training via explicit loss function employ an explicit loss function to
capture the robustness criteria and optimize it. Stability training (Zheng et
al. 2016a; Li et al. 2019) for images is based on the criteria that similar
inputs should produce similar DNN outputs. Adversarial training can be
interpreted as min-max optimization, where a hand-designed optimizer such as
projected gradient descent is employed to (approximately) solve inner
maximization. (Xiong and Hsieh 2020) train a neural network to guide the
optimizer. Since characteristics of time-series (e.g., fast-pace oscillations,
sharp peaks) are different from images/text, $L_{p}$ distance based methods
are not suitable for time-series domain.

In summary, there is no prior work111In a concurrent work, (Belkhouja and
Doppa 2022) proposed an adversarial framework for time-series domain using
statistical features and also provided robustness certificates. to train
robust DNNs for time-series domain in a principled manner. This paper
precisely fills this important gap in our scientific knowledge.

## Experiments and Results

We present experimental evaluation of RO-TS on real-world time-series
benchmarks and compare with prior methods.

### Experimental Setup

Datasets. We employed diverse univariate and multi-variate time-series
benchmark datasets from the UCR repository (Bagnall et al. 2020). Table 2
describes the details of representative datasets for which we show the results
(due to space limits) noting that our overall findings were similar on other
datasets from the UCR repo. We employ the standard training/validation/testing
splits for these datasets.

Name | Classes |  Input Size ($C\times T$)
---|---|---
ECG200 | 2 | 1$\times$97
BME | 3 | 1$\times$129
ECG5000 | 5 | 1$\times$141
MoteStrain | 2 | 1$\times$85
SyntheticControl | 6 | 1$\times$61
RacketSports | 4 | 6$\times$30
ArticularyWR | 25 | 9$\times$144
ERing | 6 | 4$\times$65
FingerMovements | 2 | 28$\times$50
Table 2: Description of different datasets.

Algorithmic setup and baselines. We employ a 1D-CNN architecture (Bai, Kolter,
and Koltun 2018) as the deep model for our evaluation. The details of the
neural architecture are provided in the Appendix. We ran RO-TS algorithm shown
in Appendix for a maximum of 500 iterations to train robust models. To
estimate GAK distance within RO-TS, we employed 15 percent of the total
alignments noting that larger sample sizes didn’t improve the optimization
accuracy and increased the training time. We also employ adversarial training
to create models using baseline attacks that are not specific to image domain
for comparison: Fast Gradient Sign method (FGS) (Kurakin, Goodfellow, and
Bengio 2017) that was used by Fawaz et al. (2019) and Projected Gradient
Descent (PGD)(Madry et al. 2017). We also compare RO-TS against stability
training (STN) (Zheng et al. 2016b).

Evaluation metrics. We evaluate the robustness of created models using
different attack strategies on the testing data. The prediction accuracy of
each model (via ground-truth labels of time-series) is used as the metric. To
ensure robustness, DNN models should be least sensitive to different types of
perturbations over original time-series signals. We measure the accuracy of
each DNN model against: 1) Adversarial noise is introduced by FGS and PGD
baseline attacks; and 2) Gaussian noise $\sim~{}\mathcal{N}(0,\Sigma)$ that
may naturally occur to perturb time-series. The covariance matrix $\Sigma$
diagonal elements (i.e., variances) are all equal to $\sigma$. DNNs are
considered robust if they are successful in maintaining their accuracy
performance against such noises.

### Results and Discussion

Random Noise

Adversarial Noise

Figure 1: Comparison of RO-TS algorithm vs. adversarial training algorithm
using baseline FGS and PGD attacks.

Random Noise

Adversarial Noise

Figure 2: Comparison of RO-TS algorithm using GAK distance
($k_{\mathrm{GAK}}$) vs. RO-TS using Euclidean distance ($L_{2}$).

Random Noise

Adversarial Noise

Figure 3: Comparison of RO-TS vs. stability training (STN).

RO-TS vs. adversarial training. One of our key hypothesis is that Euclidean
distance-based perturbations do not capture the appropriate notion of
invariance for time-series domain to improve the robustness of the learned
model. We show that using baseline attacks to create augmented data for
adversarial training does not create robust models. From Figure 1, we can
observe that models from our RO-TS algorithm achieve significantly higher
accuracy than the baselines. For example, on MoteStrain dataset, RO-TS has a
steady performance against both types of noises, unlike the baselines. On the
other datasets, we can clearly observe that in most cases, RO-TS outperforms
the baselines. We conclude that adversarial training using prior methods and
attack strategies is not as effective as our RO-TS method, where we perform
explicit primal-dual optimization to create robust models.

RO-TS vs. RO-TS with L2 distance. We want to demonstrate that choosing the
right distance metric to compute similarity between time-series signals is
critical to create robust models. Therefore, we compare models created by RO-
TS by using two different distance metrics: 1) The standard Euclidean distance
$\|\cdot\|_{2}$ used in image domains and prior work; and 2) Using the GAK
distance $D_{\mathrm{GAK}}$. From Figure 2, we can clearly see that the GAK
distance is able to explore the time-seris input space better to improve
robustness. The Euclidean distance either performs significantly worse than
GAK (e.g., on ECG5000, ERing, and RacketSports datasets) or performs
comparably to GAK (e.g., on SyntheticControl or ArticularyWR datasets). This
experiment concludes that GAK is a suitable distance metric for time-series
domain.

RO-TS vs. stability training. Unlike adversarial training, stability training
(STN) employs the below loss function (Zheng et al. 2016b) to introduce
stability to the deep model.

$Loss_{STN}=L_{0}(x,\theta)+0.01\times L_{stability}(x,x^{\prime},\theta)$

where $x$ is the original input, $x^{\prime}$ is a perturbed version of $x$
using additive Gaussian noise $\sim~{}\mathcal{N}(0,0.04^{2})$, $L_{0}$ is the
cross-entropy loss, and $L_{stability}$ relies on KL-divergence. We
experimentally demonstrate that RO-TS formulation is more suitable than STN
for creating robust DNNs for time-serirs domain. Figure 3 shows a comparison
between DNNs trained using STN and RO-TS. We observe that for most datasets,
RO-TS creates significantly more robust DNNs when compared to STN for both
types of perturbations. RO-TS algorithm is specifically designed for time-
series domain by making appropriate design choices, whereas STN is designed
for image domain. Hence, RO-TS allows us to create more robust DNNs for time-
series domain.

Empirical convergence. We demonstrate the efficiency of RO-TS algorithm by
observing the empirical rate of convergence. Figure 4 shows the optimization
objective over iterations on some representative datasets noting that we
observe similar patterns on other datasets. We can observe that RO-TS
converges roughly before 150 iterations for most datasets.

Figure 4: Empirical convergence of RO-TS algorithm.

Figure 5 shows the accuracy gap results and the computational runtime when
comparing RO-TS with sampled alignments and original GAK (i.e., all alignment
paths). The results clearly match with our theoretical analysis that accuracy
gap decreases over training iterations leading to convergence. We conclude
from these results that RO-TS converges quickly in practice and supports our
theoretical analysis.

Figure 5: The accuracy gap in the gradients over weights $\Delta G_{W}$ and
over perturbations $\Delta G_{a}$ using 5% of alignments and GAK using all
alignments for RO-TS training on ERing (Left) and the comparison of the
computational runtime between both settings of RO-TS (Right).

## Conclusions

We introduced the RO-TS algorithm to train robust deep neural networks (DNNs)
for time-series domain. The training problem was formulated as a min-max
optimization problem to reason about the worst-case risk in a small
neighborhood defined by the global alignment kernel (GAK) based distance. Our
proposed stochastic compositional alternating gradient descent and ascent
(SCAGDA) algorithm carefully leverages the structure of the optimization
problem to solve it efficiently. Our theoretical and empirical analysis showed
that RO-TS and SCAGDA are effective in creating more robust DNNs over prior
methods and GAK based distance is better suited for time-series over the
Euclidean distance.

## Acknowledgments

This research is supported in part by the AgAID AI Institute for Agriculture
Decision Support, supported by the National Science Foundation and United
States Department of Agriculture - National Institute of Food and Agriculture
award #2021-67021-35344.

## References

  * Bagnall et al. (2020) Bagnall, A.; Lines, J.; Vickers, W.; and Keogh, E. 2020.  The UEA & UCR Time Series Classification Rep.  www.timeseriesclassification.com.  Accessed: 2021-08-02.
  * Bai, Kolter, and Koltun (2018) Bai, S.; Kolter, J. Z.; and Koltun, V. 2018.  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.  _arXiv preprint arXiv:1803.01271_.
  * Belkhouja and Doppa (2020) Belkhouja, T.; and Doppa, J. R. 2020.  Analyzing Deep Learning for Time-Series Data Through Adversarial Lens in Mobile and IoT Applications.  _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_.
  * Belkhouja and Doppa (2022) Belkhouja, T.; and Doppa, J. R. 2022.  Adversarial Framework with Certified Robustness for Time-Series Domain via Statistical Features.  _Journal of Artificial Intelligence Research (JAIR)_.
  * Berndt and Clifford (1994) Berndt, D. J.; and Clifford, J. 1994.  Using dynamic time warping to find patterns in time series.  In _KDD workshop_.
  * Brendel, Rauber, and Bethge (2018) Brendel, W.; Rauber, J.; and Bethge, M. 2018.  Decision-Based Adversarial Attacks: Reliable Attacks Against Black-Box Machine Learning Models.  In _ICLR_.
  * Carlini and Wagner (2017) Carlini, N.; and Wagner, D. A. 2017.  Towards Evaluating the Robustness of Neural Networks.  In _IEEE Symposium on Security and Privacy_.
  * Chen, Sun, and Yin (2020) Chen, T.; Sun, Y.; and Yin, W. 2020.  Solving stochastic compositional optimization is nearly as easy as solving stochastic optimization.  _arXiv preprint arXiv:2008.10847_.
  * Cuturi (2011) Cuturi, M. 2011.  Fast global alignment kernels.  In _ICML_.
  * Cuturi et al. (2007) Cuturi, M.; Vert, J.; Birkenes, O.; and Matsui, T. 2007.  A kernel for time series based on global alignments.  In _ICASSP_.
  * Fawaz et al. (2019) Fawaz, H. I.; Forestier, G.; Weber, J.; Idoumghar, L.; and Muller, P. 2019.  Adversarial Attacks on Deep Neural Networks for Time Series Classification.  In _IJCNN_.
  * Fazle, Somshubra, and Houshang (2020) Fazle, K.; Somshubra, M.; and Houshang, D. 2020.  Adversarial attacks on time series.  _IEEE Transactions on pattern analysis and machine intelligence_.
  * Gao et al. (2018) Gao, J.; Lanchantin, J.; Soffa, M. L.; and Qi, Y. 2018.  Black-box generation of adversarial text sequences to evade deep learning classifiers.  In _IEEE Security and Privacy Workshops (SPW)_.
  * Hosseini et al. (2017) Hosseini, H.; Xiao, B.; Jaiswal, M.; and Poovendran, R. 2017.  On the limitation of convolutional neural networks in recognizing negative images.  In _ICMLA_.
  * Junchi Y. (2020) Junchi Y., N. H., Negar K. 2020.  Global Convergence and Variance Reduction for a Class of Nonconvex-Nonconcave Minimax Problems.  In _NeurIPS_.
  * Karimi, Nutini, and Schmidt (2016) Karimi, H.; Nutini, J.; and Schmidt, M. 2016.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In _ECML_.
  * Kolter and Madry (2018) Kolter, Z.; and Madry, A. 2018.  Tutorial adversarial robustness: Theory and practice.  _NeurIPS_.
  * Kurakin, Goodfellow, and Bengio (2017) Kurakin, A.; Goodfellow, I. J.; and Bengio, S. 2017.  Adversarial examples in the physical world.  In _ICLR, Workshop Track Proceedings_.
  * Laidlaw and Feizi (2019) Laidlaw, C.; and Feizi, S. 2019.  Functional adversarial attacks.  In _NeurIPS_.
  * Li et al. (2019) Li, B.; Chen, C.; Wang, W.; and Carin, L. 2019.  Certified Adversarial Robustness with Additive Noise.  In _NeurIPS_.
  * Lin, Jin, and Jordan (2020) Lin, T.; Jin, C.; and Jordan, M. I. 2020.  Near-optimal algorithms for minimax optimization.  In _Conference on Learning Theory_.
  * Madry et al. (2017) Madry, A.; Makelov, A.; Schmidt, L.; Tsipras, D.; and Vladu, A. 2017.  Towards deep learning models resistant to adversarial attacks.  _arXiv preprint arXiv:1706.06083_.
  * Moosavi-Dezfooli et al. (2017) Moosavi-Dezfooli, S.; Fawzi, A.; Fawzi, O.; and Frossard, P. 2017.  Universal Adversarial Perturbations.  In _CVPR_.
  * Papernot et al. (2018) Papernot, N.; Faghri, F.; Carlini, N.; Goodfellow, I.; Feinman, R.; Kurakin, A.; Xie, C.; Sharma, Y.; Brown, T.; Roy, A.; Matyasko, A.; Behzadan, V.; Hambardzumyan, K.; Zhang, Z.; Juang, Y.; Li, Z.; Sheatsley, R.; Garg, A.; Uesato, J.; Gierke, W.; Dong, Y.; Berthelot, D.; Hendricks, P.; Rauber, J.; and Long, R. 2018.  CleverHans v2.1.0 Adversarial Examples Library.  _arXiv preprint arXiv:1610.00768_.
  * Samanta and Mehta (2017) Samanta, D.; and Mehta, S. 2017.  Towards crafting text adversarial samples.  _arXiv preprint arXiv:1707.02812_.
  * Wang, Fang, and Liu (2017) Wang, M.; Fang, E. X.; and Liu, H. 2017.  Stochastic compositional gradient descent: algorithms for minimizing compositions of expected-value functions.  _Mathematical Programming, Springer_.
  * Wang, Singh, and Li (2019) Wang, W.; Singh, S.; and Li, J. 2019.  Deep Adversarial Learning for NLP.  In _NAACL-HLT, Tutorial Abstracts_.
  * Wang, Yan, and Oates (2017) Wang, Z.; Yan, W.; and Oates, T. 2017.  Time series classification from scratch with deep neural networks: A strong baseline.  In _IJCNN_.
  * Wu et al. (2018) Wu, L.; Yen, I. E.-H.; Yi, J.; Xu, F.; Lei, Q.; and Witbrock, M. 2018.  Random warping series: A random features method for time-series embedding.  In _International Conference on Artificial Intelligence and Statistics_. PMLR.
  * Xiao et al. (2018) Xiao, C.; Zhu, J.; Li, B.; He, W.; Liu, M.; and Song, D. 2018.  Spatially Transformed Adversarial Examples.  In _ICLR_.
  * Xiong and Hsieh (2020) Xiong, Y.; and Hsieh, C.-J. 2020.  Improved Adversarial Training via Learned Optimizer.  In _European Conference on Computer Vision_ , 85–100. Springer.
  * Yan et al. (2020) Yan, Y.; Xu, Y.; Lin, Q.; Liu, W.; and Yang, T. 2020.  Optimal Epoch Stochastic Gradient Descent Ascent Methods for Min-Max Optimization.  In _NeurIPS_.
  * Zheng et al. (2016a) Zheng, S.; Song, Y.; Leung, T.; and Goodfellow, I. J. 2016a.  Improving the Robustness of Deep Neural Networks via Stability Training.  In _CVPR_.
  * Zheng et al. (2016b) Zheng, S.; Song, Y.; Leung, T.; and Goodfellow, I. J. 2016b.  Improving the Robustness of Deep Neural Networks via Stability Training.  In _CVPR_.

##  Appendix A Experimental Details

We employed the standard benchmark training, validation, and test split on the
UCR datasets. Each experiment was repeated 10 times and the average has been
reported. The shaded region in the graphs depicts the minimum and maximum
value achieved within those experiments. We implemented our proposed RO-TS
algorithm using TensorFlow and the baselines using the CleverHans library
(Papernot et al. 2018).

DNN Architecture. We employed a 1D-CNN architecture for all our experiments
(both RO-TS and baselines). The architecture as follows: {C:100; K:5; C:50;
K:5; P:4; R:200; R:100} where C: Convolutional layers, K: kernel size, P: max-
pooling kernel size, and R: rectified linear unit layer.

RO-TS. We provide the pseudo-code of RO-TS to train robust DNNs in Algorithm
2. The $\nu$ hyper-parameter is estimated as the median distance of different
points observed in different time-series of the training set. It is then
scaled by the square root of the median length of time-series in the training
set (Cuturi 2011). For all our experiments reported in Section Experiments and
Results, we employed $\lambda=10\mathrm{e}^{-2}$ to prevent exploding gradient
problems and a fixed $\beta\in[0,1]$ value based on validation experiments to
compute the moving average for GAK distance. We employed a mini-batch size
$s\in[10,25]$ in all our experiments. The source code of RO-TS algorithms is
available at Robust-Training-for-Time-Series GitHub

Baselines. To create the PGD and FGS baselines, we use adversarial examples
created using perturbation factors $\epsilon<$ 1\. This bound is used for two
main reasons. First, larger perturbations significantly degrade the overall
performance of adversarial training. Second, we want to avoid the risk of
leaking label information (Madry et al. 2017). Additionally, the STN baseline
has been implemented as described in the original paper (Zheng et al. 2016b).
The adversarial training using the baseline attacks employ the same
architecture and use Adam optimizer.

GAK distance. For the computation of $k_{\mathrm{GAK}}(x,x^{\prime})$ for
time-series signals $x$ and $x^{\prime}$: 1) Defining the set of paths $\pi$:
Every path $\pi$ uses a set of two coordinates $\pi_{1}(i)$ and $\pi_{2}(i)$.
We randomly sample set of paths $\\{\pi\\}$ for GAK computation such that
$\|\pi_{1}(i)-\pi_{2}(i)\|<\frac{T}{2}$. This procedure avoids clearly bad
alignments, where few time-steps from the first signal are matched with many
time-steps from the second signal. 2) Fixing the value of $\nu$: As $\nu$ is
the divisor of $d_{\pi}(x,x^{\prime})$, we compute the pairwise distance
$d_{\pi_{diag}}(\cdot,\cdot)$ between all pairs of time-series signals from
the training data and set $\nu$ as the median of all the different distances.
$\pi_{diag}$ is the alignment path where we have a point-to-point match
($\|\pi_{1}(i)-\pi_{2}(i)\|=0$)

##  Appendix B Proofs of Technical Lemmas

Before we prove our results, we first present the detailed definitions for our
assumptions.

######  Definition 1.

If a function $f(x)$ is $C$-Lipschitz continuous, then we have

$\displaystyle|f(x_{1})-f(x_{2})|\leq C\|x_{1}-x_{2}\|\text{ and }\|\nabla
f(x)\|\leq C.$

######  Definition 2.

If a function $f(x)$ is $L$-smooth, then its gradient is $L$-Lipschitz
continuous:

$\displaystyle\|\nabla f(x_{1})-\nabla f(x_{2})\|\leq L\|x_{1}-x_{2}\|,$

and equivalently

$\displaystyle f(x_{1})-f(x_{2})\leq\langle\nabla
f(x_{2}),x_{1}-x_{2}\rangle+\frac{L}{2}\|x_{1}-x_{2}\|^{2}.$

######  Definition 3.

If a function $f(x)$ has a non-empty solution set and satisfies $\mu$-PL
condition, then we have

$\displaystyle\|\nabla f(w)\|^{2}\geq
2\mu(f(x)-\min_{x^{\prime}}f(x^{\prime})).$

We now present two important lemmas that can be directly derived from our
Assumption 1, whose proofs have been done by previous studies.

######  Lemma 1.

(Lemma A.2 and A.3 in (Junchi Y. 2020)) If $\phi$ is $L$-smooth and satisfies
two side $\mu$-PL condition, then $P(w)$ also satisfies $\mu$-PL conditions,
and $P(w)$ is $L_{P}$-smooth with $L_{P}=L+L^{2}/\mu\leq 2L^{2}/\mu$ when
$L\geq 1$ and $\mu\leq 1$.

######  Lemma 2.

((Karimi, Nutini, and Schmidt 2016)) If $f(x)$ is $L$-smooth and satisfies
$\mu$-PL condition, then it also satisfies $\mu$-error bound condition
$\|\nabla f(x)\|\geq\mu\|x_{p}-x\|$ where $x_{p}$ is the projection of $x$
onto the optimal set. Also it satisfies $\mu$-quadratic growth condition
$f(x)-\min_{x}f(x)\geq\frac{\mu}{2}\|x_{p}-x\|^{2}$. Conversely, if $f(x)$ is
$L$-smooth and satisfies $mu$-error bound, then it satisfies
$\frac{\mu}{L}$-PL condition.

From the next subsection, we give the detailed proofs for our Theorem 1 and 2
The rate of convergence rate is a quantity that represents how quickly the
sequence approaches its limit. In our convergence analysis for our proposed
optimization algorithm, this reveals the upper bound for the primal gap, i.e.,
$P(w)-P(w^{*})$ where $w^{*}\in\arg\min_{w\in\Theta}P(w)$, as a function of
$K$, the total number of iterations.

###  B.1  Proof of Theorem 1

Before we can start proof of the main theorem, we need the following three
technical lemmas.

######  Lemma 3.

Suppose Assumption 1 holds. Assume $P$ is $L_{P}$-smooth. For iteration $k$ of
Algorithm 1, we have

$\displaystyle\mathbb{E}[P(w_{k+1})-P^{*}]\leq$ $\displaystyle
P(w_{k})-P^{*}-\frac{\eta_{k}}{2}\|\nabla
P(w_{k})\|^{2}+\frac{L\eta_{k}^{2}\sigma^{2}}{2}$
$\displaystyle+\frac{\eta_{k}}{2}\|\nabla_{w}\phi(w_{k},a_{k})-\nabla
P(w_{k})\|^{2}$ (7)

######  Lemma 4.

Suppose Assumption 1 hold. For iteration $k$ of Algorithm 1, we have

$\displaystyle\mathbb{E}[P(w_{k+1})-\phi(w_{k+1},a_{k+1})]$
$\displaystyle\leq(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[P(w_{k})-\phi(w_{k},a_{k})]$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})(\frac{\eta_{k}}{2}+2\eta_{k}+L\eta_{k}^{2})\|\nabla
P(w_{k})-\nabla\phi_{w}^{k}\|^{2}$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})(-\frac{\eta_{k}}{2}+2\eta_{k}+L\eta_{k}^{2})\|\nabla
P(w_{k})\|^{2}$
$\displaystyle+2\gamma_{k}C_{h}^{2}L_{g}^{2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\sigma^{2}(\eta_{k}^{2}+\gamma^{2}(1+2C_{g}^{2}))$
(8)

######  Lemma 5.

(Lemma 2 of (Wang, Fang, and Liu 2017))

$\displaystyle\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]\leq$
$\displaystyle(1-\beta_{k})\|\omega_{k}-h(a_{k-1})\|^{2}$
$\displaystyle+\frac{C_{h}^{2}\|a_{k}-a_{k-1}\|^{2}}{\beta_{k}}+2\beta_{k}^{2}\sigma^{2}$
(9)

We now have the three important inequalities that are required for our proof
for Theorem 1. The key idea is to construct a Lyapunov function that combines
the above three inequalities together, and derive the convergence for this
Lyapunov function. Let

$\displaystyle L_{k+1}:=$
$\displaystyle\mathbb{E}[P(w_{k+1})-P^{*}]+\frac{1}{8}\mathbb{E}[P(w_{k+1})-\phi(w_{k+1},a_{k})]$
$\displaystyle+\Big{(}\frac{4C_{h}^{4}L_{g}^{4}L^{4}\eta_{k}}{\mu^{5}}\Big{)}^{1/2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]$

Set $\eta_{k}=\eta$, $\gamma_{k}=\gamma=\frac{36L^{2}\eta}{\mu^{2}}$,
$\beta=(18\mu\eta)^{1/2}$. Recall that $\|\nabla
P(w_{k})-\nabla\phi_{w}^{k}\|^{2}\leq\frac{2L^{2}}{\mu}(P(w_{k})-\phi(w_{k},a_{k}))$.
By Lemma 3 and 4, we have

$\displaystyle L_{k+1}\leq$ $\displaystyle
P(w_{k})-P^{*}-\frac{\eta}{2}\|\nabla P(w_{k})\|^{2}+\frac{\eta}{2}\|\nabla
P(w_{k})-\nabla\phi_{w}^{k}\|^{2}$
$\displaystyle+\frac{L_{P}\eta^{2}\sigma^{2}}{2}+\frac{1}{8}\Bigg{(}(1-\frac{\mu\gamma}{2})\mathbb{E}[P(w_{k})-\phi(w_{k},a_{k})]$
$\displaystyle+(1-\frac{\mu\gamma}{2})(\frac{\eta}{2}+2\eta+L\eta^{2})\|\nabla
P(w_{k})-\nabla\phi_{w}^{k}\|^{2}$
$\displaystyle+(1-\frac{\mu\gamma}{2})(-\frac{\eta}{2}+2\eta+L\eta^{2})\|\nabla
P(w_{k})\|^{2}$
$\displaystyle+L\sigma^{2}(\eta^{2}+\gamma^{2}(1+2C_{g}^{2}))\Bigg{)}$
$\displaystyle+\underbrace{\Big{(}2\gamma
C_{h}^{2}L_{g}^{2}+\Big{(}\frac{L^{4}\eta}{\mu^{5}}\Big{)}^{1/2}2C_{h}^{2}L_{g}^{2}\Big{)}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]}_{=A}$
(10)

Recall that $\eta\leq\frac{1}{4L}$. Before bounding term $A$, we first
simplify the coefficients of $\mathbb{E}[\|\nabla P(w_{k})\|^{2}]$ and
$\mathbb{E}[\|\nabla P(w_{k})-\nabla\phi_{w}(w_{k},a_{k})\|^{2}]$ as follows

$\displaystyle\text{for }[\|\nabla P(w_{k})\|^{2}]:$
$\displaystyle-\frac{\eta}{2}+\frac{1}{8}(1-\frac{\mu\gamma}{2})(-\frac{\eta}{2}+2\eta+L\eta^{2})\leq-\frac{\eta}{4},$
$\displaystyle\text{for }\mathbb{E}[\|\nabla
P(w_{k})-\nabla\phi_{w}(w_{k},a_{k})\|^{2}]:$
$\displaystyle\frac{\eta}{2}+\frac{1}{8}(1-\frac{\mu\gamma}{2})(\frac{\eta}{2}+2\eta+L\eta^{2})\leq\frac{7\eta}{8}.$

Due to $L$-smoothness of $\phi$ in $a$ and Lemma 2, we can merge
$\mathbb{E}[\|\nabla P(w_{k})-\nabla\phi_{w}(w_{k},a_{k})\|^{2}]$ with
$P(w_{k})-\phi(w_{k},a_{k})$. Then the coefficient of
$P(w_{k})-\phi(w_{k},a_{k})$ in (B.1) can be derived as follows

$\displaystyle\frac{1}{8}(1-\frac{\mu\gamma}{2})+\frac{2L^{2}}{\mu}\cdot\frac{7\eta}{8}=\frac{1}{8}(1-\frac{18L^{2}\eta}{\mu}+14\frac{L^{2}\eta}{\mu})\leq\frac{1}{8}(1-\frac{\mu\eta}{4}),$

where the inequality is due to $\mu\leq L$. Now we employ Lemma 5 to bound
term $A$ in (B.1) as follows

$\displaystyle\Big{(}2\gamma_{k}C_{h}^{2}L_{g}^{2}+\Big{(}\frac{L^{4}\eta}{\mu^{5}}\Big{)}^{1/2}2C_{h}^{2}L_{g}^{2}\Big{)}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]$
$\displaystyle\leq$ $\displaystyle
2C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(18\sqrt{\mu\eta}+1)(1-\beta)\|\omega_{k}-h(a_{k})\|^{2}$
$\displaystyle+2C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(18\sqrt{\mu\eta}+1)(\frac{C_{h}^{2}\gamma^{2}C_{h}^{2}C_{g}^{2}}{\beta}+2\beta^{2}\sigma^{2})$
$\displaystyle\leq$ $\displaystyle
2C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(1-\frac{\mu\eta}{4})\|\omega_{k}-h(a_{k})\|^{2}$
$\displaystyle+20C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(\frac{36^{2}C_{h}^{2}C_{g}^{2}L^{4}\eta^{2}}{18\mu^{4}\sqrt{\mu\eta}}+2\sigma^{2}18^{2}\mu\eta).$
(11)

By plugging the above simplified coefficients and merging (B.1) into (B.1), we
have the following

$\displaystyle L_{k+1}\leq$ $\displaystyle
P(w_{k})-P^{*}-\frac{\eta}{4}\|\nabla P(w_{k})\|^{2}$
$\displaystyle+\frac{1}{8}(1-\frac{\mu\eta}{4})(P(w_{k})-\phi(w_{k},a_{k}))$
$\displaystyle+\frac{L_{P}\eta^{2}\sigma^{2}}{2}+\frac{1}{8}L\sigma^{2}(\eta^{2}+\gamma^{2}(1+2C_{g}^{2}))$
$\displaystyle+2C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(1-\frac{\mu\eta}{4})\|\omega_{k}-h(a_{k})\|^{2}$
$\displaystyle+20C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(\frac{36^{2}C_{h}^{2}C_{g}^{2}L^{4}\eta^{2}}{18\mu^{4}\sqrt{\mu\eta}}+2\sigma^{2}18^{2}\mu\eta)$
$\displaystyle\leq$
$\displaystyle(1-\frac{\mu\eta}{4})(P(w_{k})-P^{*})+\frac{1}{8}(1-\frac{\mu\eta}{4})(P(w_{k})-\phi(w_{k},a_{k}))$
$\displaystyle+2C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(1-\frac{\mu\eta}{4})\|\omega_{k}-h(a_{k})\|^{2}+B$
$\displaystyle=$ $\displaystyle(1-\frac{\mu\eta}{4})L_{k}+B$

where the last inequality is due to $\mu$-PL condition of $P$ and we define

$\displaystyle
B:=\frac{L_{P}\eta^{2}\sigma^{2}}{2}+\frac{1}{8}L\sigma^{2}(\eta^{2}+\gamma^{2}(1+2C_{g}^{2}))$
$\displaystyle+20C_{h}^{2}L_{g}^{2}L^{2}/\mu^{2}\sqrt{\frac{\eta}{\mu}}(\frac{36^{2}C_{h}^{2}C_{g}^{2}L^{4}\eta^{2}}{18\mu^{4}\sqrt{\mu\eta}}+2\sigma^{2}18^{2}\mu\eta)$

Finally, we expand the above recursion from $K$ to $1$ and have

$\displaystyle
L_{K}\leq(1-\frac{\mu\eta}{4})^{K}L_{0}+B\sum_{k=0}^{K-1}(1-\frac{\mu\eta}{4})^{k}\leq\exp(-\frac{\mu\eta
K}{4})L_{0}+\frac{4B}{\mu\eta}$

We can verify that $\frac{4B}{\mu\eta}\leq 3\epsilon/4$ if $\eta$ satisfies
the following three inequalities:

$\displaystyle\eta\leq$
$\displaystyle\frac{\mu^{4}\epsilon}{16\sigma^{2}(163+324C_{g}^{2})L^{5}},~{}~{}\eta\leq\frac{\mu^{8}}{23040C_{h}^{6}C_{g}^{2}L_{g}^{2}L^{6}}$
$\displaystyle\eta\leq$ $\displaystyle\frac{\mu^{5}\epsilon^{2}}{(4\cdot
160\cdot 324)^{2}C_{h}^{2}L_{g}^{2}L^{4}}$

In practice, $\epsilon$ is usually set to a very small value, so we employ the
last inequality as the upper bound of $\eta$. Then, to make
$\exp(-\frac{\mu\eta K}{4})L_{0}\leq\epsilon/4$, $K$ has the following lower
bound.

$\displaystyle K\geq\frac{4}{\mu\eta}\log(4L_{0}/\epsilon)=\frac{(8\cdot
160\cdot
324)^{2}C_{h}^{2}L_{g}^{2}L^{4}}{\mu^{6}\epsilon^{2}}\log(4L_{0}/\epsilon)$

Since $L_{k+1}$ is an upper bound of $\mathbb{E}[P(w_{k+1})-P^{*}]$, to
$\epsilon$-primal gap, we have the total number of iterations $K$ =
$\tilde{O}\Big{(}\frac{L^{4}}{\mu^{6}\epsilon^{2}}\Big{)}$.

###  B.2  Proof of Theorem 2

To prove Theorem 2, let
$\Delta_{K+1}=\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]$ and use Lemma 5:

$\displaystyle\Delta_{K}\leq$
$\displaystyle(1-\beta)\Delta_{K-1}+\frac{C_{h}^{2}\gamma^{2}C_{h}^{2}C_{g}^{2}}{\beta}+2\beta^{2}\sigma^{2}$
$\displaystyle\leq$
$\displaystyle(1-\beta)^{K-1}\Delta_{1}+\frac{C_{h}^{2}\gamma^{2}C_{h}^{2}C_{g}^{2}}{\beta^{2}}+2\beta\sigma^{2}\leq\epsilon,$

where the second inequality is due to
$\sum_{k=0}^{K-1}(1-\beta)^{k}\leq\frac{1}{\beta}$. The last inequality is due
to

$\displaystyle O((1-\beta)^{K-1})\leq O(\exp(-\beta(K-1)))$
$\displaystyle=O(\exp(-\sqrt{\mu\eta}\frac{L^{6}\log(1/\epsilon)}{\mu^{6}\epsilon^{2}}))\leq
O(\epsilon^{1/\epsilon})\leq O(\epsilon),$ $\displaystyle
O(\frac{\gamma^{2}}{\beta^{2}})\leq
O(\frac{L^{6}\eta^{2}}{\mu^{4}\mu\eta})\leq
O(\frac{L^{6}\eta}{\mu^{4}\mu})\leq O(\epsilon),$ $\displaystyle O(\beta)\leq
O(\sqrt{\mu\eta})\leq O(\epsilon).$

###  B.3 Proof of Lemma 3

###### Proof.

(of Lemma 3)

Conditioned on iteration $k$, by $L_{P}$-smoothness of $P$, we have

$\displaystyle\mathbb{E}[P(w_{k+1})-P(x_{k})]\leq\mathbb{E}[\langle\nabla
P(w_{k}),w_{k+1}-w_{k}\rangle]+\mathbb{E}[\frac{L_{P}}{2}\|w_{k+1}-w_{k}\|^{2}]$
$\displaystyle=$ $\displaystyle\mathbb{E}[\langle\nabla
P(w_{k}),-\eta_{k}\nabla_{w}\phi_{i_{1}}(w_{k},a_{i_{1}})\rangle]+\mathbb{E}[\frac{L_{P}}{2}\|-\eta_{k}\nabla_{w}\phi_{i_{1}}(w_{k},a_{i_{1}})\|^{2}]$
$\displaystyle\leq$ $\displaystyle-\eta_{k}\langle\nabla
P(w_{k}),\nabla_{w}\phi_{i_{1}}(w_{k},a_{i_{1}})\rangle+\frac{L_{P}\eta_{k}^{2}}{2}\mathbb{E}[\|\nabla_{w}\phi(w_{k},a)\|^{2}]+\frac{L_{P}\eta_{k}^{2}\sigma^{2}}{2}$
$\displaystyle\leq$ $\displaystyle-\frac{\eta_{k}}{2}\|\nabla
P(w_{k})\|^{2}+\frac{\eta_{k}}{2}\|\nabla
P(w_{k})-\nabla_{w}\phi(w_{k},a)\|^{2}+\frac{L_{P}\eta_{k}^{2}\sigma^{2}}{2},$

where the above first inequality is due to

$\displaystyle\mathbb{E}[\|\nabla_{w}\phi_{i_{1}}(w_{k},a_{i_{1}})-\nabla_{w}\phi(w_{k},a)+\nabla_{w}\phi(w_{k},a)\|^{2}]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\|\nabla_{w}\phi_{i_{1}}(w_{k},a_{i_{1}})-\nabla_{w}\phi(w_{k},a)\|^{2}]+\mathbb{E}[\|\nabla_{w}\phi(w_{k},a)\|^{2}]$
$\displaystyle+2\mathbb{E}[\langle\nabla_{w}\phi_{i_{1}}(w_{k},a_{i_{1}})-\nabla_{w}\phi(w_{k},a),\nabla_{w}\phi(w_{k},a)\rangle]$
$\displaystyle\leq$
$\displaystyle\sigma^{2}+\mathbb{E}[\|\nabla_{w}\phi(w_{k},a)\|^{2}]$

and the second inequality is due to $\eta_{k}\leq 1/(4L_{P})$. ∎

###  B.4 Proof of Lemma 4

###### Proof.

(of Lemma 4)

The following lemma is the standard analysis for one step update using
stochastic gradient descent except that we keep the biased estimation term.

######  Lemma 6.

Suppose Assumption 1 hold. For iteration $k$ of Algorithm 1, we have

$\displaystyle P(w_{k+1})-\phi(w_{k+1},a_{k+1})$ $\displaystyle\leq$
$\displaystyle(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[P(w_{k+1})-\phi(w_{k+1},a_{k})]+L\gamma_{k}^{2}\sigma^{2}(1+2C_{g}^{2})$
$\displaystyle+\gamma_{k}C_{h}^{2}L_{g}^{2}(\frac{1}{2}(1-\frac{L\gamma_{k}}{2})+2L\gamma_{k})\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}].$

Now we can prove Lemma 4. Denote
$\nabla\phi_{w}^{k}=\nabla_{w}\phi(w_{k},a_{k})$ and
$\hat{\nabla}\phi_{w}^{k}=\nabla_{w}\phi_{i_{1}}(w_{k},a_{i_{1}}^{k})$.

By Lemma 6, we have

$\displaystyle P(w_{k+1})-\phi(w_{k+1},a_{k+1})$ $\displaystyle\leq$
$\displaystyle(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[P(w_{k+1})-\phi(w_{k+1},a_{k})]$
$\displaystyle+\gamma_{k}C_{h}^{2}L_{g}^{2}(\frac{1}{2}(1-\frac{L\gamma_{k}}{2})+2L\gamma_{k})\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\gamma_{k}^{2}\sigma^{2}(1+2C_{g}^{2})$
$\displaystyle\leq$
$\displaystyle(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[P(w_{k})-\phi(w_{k},a_{k})+P(w_{k+1})-P(w_{k})+\phi(w_{k},a_{k})$
$\displaystyle-\phi(w_{k+1},a_{k})]+2\gamma_{k}C_{h}^{2}L_{g}^{2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\gamma_{k}^{2}\sigma^{2}(1+2C_{g}^{2})$
$\displaystyle\leq$
$\displaystyle(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[P(w_{k})-\phi(w_{k},a_{k})]$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})(-\frac{\eta_{k}}{2}\|\nabla
P(w_{k})\|^{2}+\frac{\eta_{k}}{2}\|\nabla
P(w_{k})-\nabla\phi_{w}^{k}\|^{2}+\frac{L\eta_{k}^{2}\sigma^{2}}{2})$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[(\langle-\nabla\phi_{w}^{k},w_{k+1}-w_{k}\rangle+\frac{L}{2}\|w_{k+1}-w_{k}\|^{2})]$
$\displaystyle+2\gamma_{k}C_{h}^{2}L_{g}^{2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\gamma_{k}^{2}\sigma^{2}(1+2C_{g}^{2})$
$\displaystyle=$
$\displaystyle(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[P(w_{k})-\phi(w_{k},a_{k})]$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})(-\frac{\eta_{k}}{2}\|\nabla
P(w_{k})\|^{2}+\frac{\eta_{k}}{2}\|\nabla
P(w_{k})-\nabla\phi_{w}^{k}\|^{2}+\frac{L\eta_{k}^{2}\sigma^{2}}{2})$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})(\mathbb{E}[\eta_{k}\langle\nabla\phi_{w}^{k},\hat{\nabla}\phi_{w}^{k}\rangle]+\frac{L\eta_{k}^{2}}{2}\|\hat{\nabla}\phi_{w}^{k}-\nabla\phi_{w}^{k}+\nabla\phi_{w}^{k}\|^{2})$
$\displaystyle+2\gamma_{k}C_{h}^{2}L_{g}^{2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\gamma_{k}^{2}\sigma^{2}(1+2C_{g}^{2})$
$\displaystyle+\frac{L\eta_{k}^{2}}{2}(\|\hat{\nabla}\phi_{w}^{k}-\nabla\phi_{w}^{k}\|^{2}+2\|\nabla\phi_{w}^{k}-\nabla
P(w_{k})\|^{2}+2\|\nabla P(w_{k})\|^{2}))$
$\displaystyle+2\gamma_{k}C_{h}^{2}L_{g}^{2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\gamma_{k}^{2}\sigma^{2}(1+2C_{g}^{2})$
$\displaystyle\leq$
$\displaystyle(1-\frac{\mu\gamma_{k}}{2})\mathbb{E}[P(w_{k})-\phi(w_{k},a_{k})]$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})(\frac{\eta_{k}}{2}+2\eta_{k}+L\eta_{k}^{2})\|\nabla
P(w_{k})-\nabla\phi_{w}^{k}\|^{2}$
$\displaystyle+(1-\frac{\mu\gamma_{k}}{2})(-\frac{\eta_{k}}{2}+2\eta_{k}+L\eta_{k}^{2})\nabla
P(w_{k})\|^{2}$
$\displaystyle+2\gamma_{k}C_{h}^{2}L_{g}^{2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\sigma^{2}(\eta_{k}^{2}+\gamma^{2}(1+2C_{g}^{2}))$

where the third inequality is due to Lemma 3 and $L$-smoothness of $\phi(w,a)$
in $a$. The fourth inequality is due to $(a+b)^{2}\leq 2a^{2}+2b^{2}$. ∎

###  B.5 Proof of Lemma 6

###### Proof.

(of Lemma 6)

Denote
$\nabla\phi_{a}^{k}=\nabla_{a}\phi(w_{k+1},a_{k})=\nabla_{a}f(w_{k+1},a_{k})-\nabla
h(a_{k})^{\top}\nabla g(\omega_{k+1})$ and
$\hat{\nabla}\phi_{a}^{k}=\nabla_{a}f_{i_{2}}(w_{k+1},a_{i_{2}}^{k})-\nabla
h_{j_{2}}(a_{i_{2}}^{k})^{\top}\nabla g(\omega_{i_{2}}^{k+1})$ for simplicity.
Note that
$\mathbb{E}[\nabla_{a}f_{i_{2}}(w_{k+1},a_{i_{2}}^{k})]=\nabla_{a}f(w_{k+1},a_{k})$,
$\mathbb{E}[\nabla h_{j_{2}}(a_{i_{2}}^{k})]=\nabla h(a_{k})$. Then we have

$\displaystyle\phi(w_{k+1},a_{k})-P(w_{k+1})+P(w_{k+1})-\phi(w_{k+1},a_{k+1})$
$\displaystyle=$ $\displaystyle-\phi(w_{k+1},a_{k+1})-(-\phi(w_{k+1},a_{k}))$
$\displaystyle\leq$
$\displaystyle\langle-\nabla\phi_{a}^{k},a_{k+1}-a_{k}\rangle+\frac{L}{2}\|a_{k+1}-a_{k}\|^{2}$
$\displaystyle=$
$\displaystyle\langle-\nabla\phi_{a}^{k},\gamma_{k}\hat{\nabla}\phi_{a}^{k}\rangle+\frac{L\gamma_{k}^{2}}{2}\|\hat{\nabla}\phi_{a}^{k}\|^{2}$
$\displaystyle=$
$\displaystyle-\gamma_{k}\langle\nabla\phi_{a}^{k},\nabla\phi_{a}^{k}-\nabla\phi_{a}^{k}+\gamma_{k}\hat{\nabla}\phi_{a}^{k}\rangle+\frac{L\gamma_{k}^{2}}{2}\|\nabla\phi_{a}^{k}-\nabla\phi_{a}^{k}+\hat{\nabla}\phi_{a}^{k}\|^{2}$
$\displaystyle=$
$\displaystyle-\gamma_{k}\|\nabla\phi_{a}^{k}\|^{2}-\gamma\langle\nabla\phi_{a}^{k},\hat{\nabla}\phi_{a}^{k}-\nabla\phi_{a}^{k}\rangle$
$\displaystyle+\frac{L\gamma_{k}^{2}}{2}(\|\hat{\nabla}\phi_{a}^{k}-\nabla\phi_{a}^{k}\|^{2}+\|\nabla\phi_{a}^{k}\|^{2}+2\langle\hat{\nabla}\phi_{a}^{k}-\nabla\phi_{a}^{k},\nabla\phi_{a}^{k}\rangle),$

where the first inequality is due to $L$-smoothness of $\phi(w,a)$ in $a$.

Taking expectation over randomness of $i_{2}$ and $j_{2}$, we have

$\displaystyle\mathbb{E}[\phi(w_{k+1},a_{k})-P(w_{k+1})+P(w_{k+1})-\phi(w_{k+1},a_{k+1})]$
$\displaystyle\leq$
$\displaystyle-\gamma_{k}(1-\frac{L\gamma_{k}}{2})\mathbb{E}[\|\nabla\phi_{a}^{k}\|^{2}]$
$\displaystyle-\gamma$
$\displaystyle(1-\frac{L\gamma_{k}}{2})\langle\nabla\phi_{a}^{k},\mathbb{E}[\nabla
h_{i_{2},j_{2}}(a_{i_{2}}^{k})]^{\top}\mathbb{E}[\nabla
g(\omega_{i_{2}}^{k+1})]-{\nabla h(a_{k})}^{\top}\nabla g(h(a_{k}))\rangle$
$\displaystyle+L\gamma_{k}^{2}\mathbb{E}\|\nabla_{a}f_{i_{2}}(w_{k+1},a_{i_{2}})-\nabla_{a}f(w_{k+1},a_{k})\|^{2}$
$\displaystyle+L\gamma_{k}^{2}\mathbb{E}\|\nabla
h_{j_{2}}(a_{i_{2}}^{k})^{\top}\nabla g(\omega^{k+1}_{i_{2}})-{\nabla
h(a_{k})}^{\top}\nabla g(h(a_{k}))\|^{2}$ $\displaystyle\leq$
$\displaystyle-\gamma_{k}(1-\frac{L\gamma_{k}}{2}-\frac{1}{2}(1-\frac{L\gamma_{k}}{2}))\mathbb{E}[\|\nabla\phi_{a}^{k}\|^{2}]$
$\displaystyle+\frac{\gamma_{k}}{2}(1-\frac{L\gamma_{k}}{2})C_{h}^{2}L_{g}^{2}\mathbb{E}[\|\omega_{k+1}-h(a_{k})\|^{2}]+L\gamma_{k}^{2}\sigma^{2}$
$\displaystyle+2L\gamma_{k}^{2}(C_{g}^{2}\sigma^{2}+C_{h}^{2}L_{g}^{2}\|\omega_{k+1}-h(a_{k})\|^{2}),$

where the first inequality is due to $(a+b)^{2}\leq 2a^{2}+2b^{2}$,the last
inequality is due to Young’s inequality $ab\leq a^{2}/2+b^{2}/2$, $E[\|\nabla
h_{i_{2},j_{2}}(a_{i_{2}}^{k})-\nabla h(a_{k})\|^{2}]\leq\sigma^{2}$ and
$\mathbb{E}[\|\nabla_{a}f_{i_{2}}(w_{k+1},a_{i_{2}})\\\
-\nabla_{a}f(w_{k+1},a_{k})\|^{2}]\leq\sigma^{2}$.

Finally, we finish the proof by setting $\gamma_{k}\leq 1/(4L)$ and the
$\mu$-PL condition
$2\mu(P(w_{k+1})-\phi(w_{k+1},a_{k}))\leq\|\nabla\phi_{a}^{k}\|^{2}$. ∎

###  B.6 Proof of Proposition 1

###### Proof.

(of Proposition 1)

Suppose $\\{y_{i}\\}_{i=1}^{n}$ where $y_{i}\geq 0$. We have

$\displaystyle\max_{i\in\\{1,...,n\\}}y_{i}\leq\sum_{i=1}^{n}y_{i}\leq
n\max_{i\in\\{1,...,n\\}}y_{i}.$

Let $y_{i}=\exp(x_{i}/\nu)$ ($x_{i}$ can be negative), we have

$\displaystyle\max_{i\in\\{1,...,n\\}}\exp(x_{i}/\nu)\leq\sum_{i=1}^{n}\exp(x_{i}/\nu)\leq
n\max_{i\in\\{1,...,n\\}}\exp(x_{i}/\nu).$

Taking logarithm, we have

$\displaystyle\log\Big{(}\max_{i\in\\{1,...,n\\}}\exp(x_{i}/\nu)\Big{)}\leq\log\Big{(}\sum_{i=1}^{n}\exp(x_{i}/\nu)\Big{)}\leq$
$\displaystyle\log\Big{(}n\max_{i\in\\{1,...,n\\}}\exp(x_{i}/\nu)\Big{)}$
$\displaystyle=$
$\displaystyle\log(n)+\log\Big{(}\max_{i\in\\{1,...,n\\}}\exp(x_{i}/\nu)\Big{)}.$

Since $\exp(x)$ is a monotone increasing function, we have
$\exp(x_{1})<\exp(x_{2})$ for $x_{1}<x_{2}$. Applying it in the above
inqualities, we have
$\log\Big{(}\max_{i\in\\{1,...,n\\}}\exp(x_{i}/\nu)\Big{)}=\log\Big{(}\exp(\max_{i\in\\{1,...,n\\}}x_{i}/\nu)\Big{)}$,
which leads to

$\displaystyle\log\Big{(}\exp(\max_{i\in\\{1,...,n\\}}x_{i}/\nu)\Big{)}\leq\log\Big{(}\sum_{i=1}^{n}\exp(x_{i}/\nu)\Big{)}\leq\log(n)+\log\Big{(}\exp(\max_{i\in\\{1,...,n\\}}x_{i}/\nu)\Big{)}.$

By cancelling $\log$ with $\exp$ and multiplying $\nu$ in both sides, we have

$\displaystyle\max_{i\in\\{1,...,n\\}}x_{i}\leq\nu\log\Big{(}\sum_{i=1}^{n}\exp(x_{i}/\nu)\Big{)}\leq\nu\log(n)+\max_{i\in\\{1,...,n\\}}x_{i}.$

By replacing $x_{i}$ with $-d_{\pi}(x,x^{\prime})$, we have

$\displaystyle\max_{i\in\\{1,...,n\\}}-d_{\pi}(x,x^{\prime})\leq\nu\log\Big{(}\sum_{i=1}^{n}\exp(-d_{\pi}(x,x^{\prime})/\nu)\Big{)}\leq\nu\log(n)+\max_{i\in\\{1,...,n\\}}-d_{\pi}(x,x^{\prime}).$

To better fit $D_{\mathrm{GAK}}$ and $D_{\mathrm{DTW}}$, we may also have the
following equivalent inequalities:

$\displaystyle-\max_{i\in\\{1,...,n\\}}-d_{\pi}(x,x^{\prime})\geq-\nu\log\Big{(}\sum_{i=1}^{n}\exp(-d_{\pi}(x,x^{\prime})/\nu)\Big{)}\geq-\nu\log(n)-\max_{i\in\\{1,...,n\\}}-d_{\pi}(x,x^{\prime})$
$\displaystyle\Rightarrow$
$\displaystyle\underbrace{\min_{i\in\\{1,...,n\\}}d_{\pi}(x,x^{\prime})}_{=D_{\mathrm{DTW}}}\geq\underbrace{-\nu\log\Big{(}\sum_{i=1}^{n}\exp(-d_{\pi}(x,x^{\prime})/\nu)\Big{)}}_{=D_{\mathrm{GAK}}}\geq-\nu\log(n)+\underbrace{\min_{i\in\\{1,...,n\\}}d_{\pi}(x,x^{\prime})}_{=D_{\mathrm{DTW}}}.$
(12)

∎

The above Equation (B.6) indicates the relationship between $D_{\mathrm{GAK}}$
and $D_{\mathrm{DTW}}$, as well as how fast $D_{\mathrm{GAK}}$ converges to
$D_{\mathrm{DTW}}$ when $\nu$ decreases and approaches to $0$.

First, for a pair of time-series signals $(x,x^{\prime})$, $D_{\mathrm{DTW}}$
(hard minimum) is always less than $D_{\mathrm{GAK}}$ (soft minimum). Second,
the gap between $D_{\mathrm{DTW}}$ and $D_{\mathrm{GAK}}$ can be upper bounded
by $\nu\log(n)$, where $n=|\mathcal{A}|$ in our case, i.e., the total number
of possible alignments between $x$ and $x^{\prime}$. Since the logarithm
function $\log(n)$ is dominated by $\nu$, we may set $\nu$ to a very small
value to approximate $D_{\mathrm{DTW}}$ using $D_{\mathrm{GAK}}$ by making
$\nu\log(n)$ very small. On the other hand, we can also tune $\nu$ in
$D_{\mathrm{GAK}}$ to provide more flexibility, which sometimes leads to
better practical performance.

###  B.7 Stochastic Approximation in Optimization Problem of Equation (6)

Denote the dual variable $a_{i}$ as a whole:
$a=[a_{1},\cdots,a_{n}]\in\mathbb{R}^{C\times T\times n}$. Then we can
concatenate all components $h_{i,j}(a_{i})$ for $i=1,\cdots,n$ and
$j=1,\cdots,m$ together:

$\displaystyle
h(a)=\frac{1}{n}\cdot\left[\begin{array}[]{cccc}\frac{1}{m}h_{1,1}(a_{1}),&\frac{1}{m}h_{1,2}(a_{1}),&...,&\frac{1}{m}h_{1,m}(a_{1})\\\
...\\\
\frac{1}{m}h_{i,1}(a_{i}),&\frac{1}{m}h_{i,2}(a_{i}),&...,&\frac{1}{m}h_{i,m}(a_{i})\\\
...\\\
\frac{1}{m}h_{n,1}(a_{n}),&\frac{1}{m}h_{n,2}(a_{n}),&...,&\frac{1}{m}h_{n,m}(a_{n})\end{array}\right]$

Note that if we sample an index $(i,j)$ to get $h(i,j)(a_{i})$, then it is a
stochastic approximation of $h(a)$. We can explicitly define the following
notation:

$\displaystyle\hat{h}_{i,j}(a_{i})=\left[\begin{array}[]{ccc}0,&...,&0\\\
...&...&...\\\ 0,&...,h_{i,j}(a_{i}),...,&0\\\ ...&...&...\\\
0,&...,&0\end{array}\right]$

Then we can compute the expectation of $\hat{h}_{i,j}(a_{i})$ as follows by
sampling each $(i,j)$ index with $\frac{1}{mn}$ probability:

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{1}{mn}\hat{h}_{i,j}(a_{i})=h(a),$

This means that $\hat{h}_{i,j}(a_{i})$ is an unbiased estimation of $h(a)$.
Similarly, we can find that $\hat{h}_{i,j}(a_{i})$ is an unbiased estimation
of $h_{i}(a_{i})$. Same development can be done for $\nabla h(a)$.