# Nonlinear Time-Series Embedding by Monotone Variational Inequality

Jonathan Y. Zhou  H. Milton Stewart School of Industrial and Systems
Engineering, Georgia Institute of Technology. Yao Xie  H. Milton Stewart
School of Industrial and Systems Engineering, Georgia Institute of Technology.

###### Abstract

In the wild, we often encounter collections of sequential data such as
electrocardiograms, motion capture, genomes, and natural language, and
sequences may be multichannel or symbolic with nonlinear dynamics. We
introduce a new method to learn low-dimensional representations of nonlinear
time series without supervision and can have provable recovery guarantees. The
learned representation can be used for downstream machine-learning tasks such
as clustering and classification. The method is based on the assumption that
the observed sequences arise from a common domain, but each sequence obeys its
own autoregressive models that are related to each other through low-rank
regularization. We cast the problem as a computationally efficient convex
matrix parameter recovery problem using monotone Variational Inequalities and
encode the common domain assumption via low-rank constraint across the learned
representations, which can learn the geometry for the entire domain as well as
faithful representations for the dynamics of each individual sequence using
the domain information in totality. We show the competitive performance of our
method on real-world time-series data with the baselines and demonstrate its
effectiveness for symbolic text modeling and RNA sequence clustering.

##  1 Introduction

Collections of time-series data, where each sequence is represented by a
series of points indexed over time, are ubiquitous and increasingly prevalent.
Notable examples include physiological signals [1, 2], power systems [3],
financial data [4, 5], computer networks [6], and electronic health records
[7, 8]. In addition to traditional time series data, other sequential data
like genes and proteins [9, 10] as well as natural language have garnered
significant recent attention, particularly with the advent of large language
models [11, 12, 13, 14].

Learning high-quality representations of sequences and time series [15] is an
essential building block [16] and important for understanding the dynamics
underlying observed sequences, enabling informed decisions and downstream
machine learning tasks [17]. A key paradigm underlying the unsupervised
representation learning has been that of self-supervised learning [18], where
we first solve some auxiliary task (e.g., autoregression or reconstruction),
which leads implicitly to a compressed representation of the input data [19,
20]. The development of self-supervised methods for natural language has been,
in turn, paralleled by embedding methods for other types of sequential data,
and there is now a burgeoning literature on time series and sequence
representation [17, 21, 22].

A lasting challenge in bringing representation learning to time series is how
to learn the information common to the entire domain in conjunction with
faithful individual representations of each sequence. In contrast to the
natural language setting where there exists a shared “universal” embedding
space for all languages [23] backed up by a large amount of available training
data, time series data are often highly domain-specific in the sense that each
domain is distinct from another (e.g., electrocardiogram (ECG) vs power
generation). Additionally, the individual processes which we receive
observations for may, in fact, be substantially different among themselves
(e.g., differences between patients and power plants). From the above, there
is a difficulty in balancing learning the common dynamics of a set of observed
sequences in addition to faithfully representing the dynamics of each
sequence. This is especially the case when observations of all sequences
individually are limited. To this end, we take inspiration from an area where
the above challenges are common and well known — low-rank matrix recovery [24,
25, 26, 27, 28] — and develop it towards the general sequential and time
series representation learning setting, enabling us to bring provable recovery
guarantees to the modeling of a broad class of sequences with autoregressive
character.

To this end, we introduce an approach to the unsupervised learning of low-
dimensional representations for nonlinear time series based on the assumption
that each sequence behaves according to its own autoregressive model, but that
the sequences are related to each other through low-rank regularization. We
cast the problem as a computationally efficient convex matrix parameter
recovery problem using monotone VIs. By enforcing the low-rank assumption
across the learned representations, we capture the geometry of the entire
domain and accurately represent the dynamics of each individual sequence
utilizing the entire domain information. We apply our method to real-world
time-series data and demonstrate its effectiveness in symbolic text modeling
and RNA sequence clustering.

###  1.1 Related work

We first review related work in time-series representation learning,
clustering, and classification. A number of methods exist to address the above
problems, the simplest of which is to extract features [29] or establish a
distance metric between time-series [30, 31, 32]. Another approach is to
construct a model for each series — our work falls within this purview of
model-based time-series representation [33, 34]. Most recent time-series
representation learning approaches adopt a contrastive learning framework to
disambiguate between sequences (e.g., by taking sub-samples of the same
sequence as positive samples and from differing sequences as negative ones)
using deep networks [35, 36, 37, 38, 39]. Recent work in this line has focused
particularly on the choice of neural architecture, the data augmentation to
ensure robustness, and the selection scheme for contrastive learning. Other
auxiliary tasks for self-supervised representation learning of time series
include reconstruction and context prediction [40, 41], also popular for
language modeling [42].

In our work, we adopt auto-regression as the auxiliary task. Unlike recent
methods that use contrastive learning to learn an encoder for the latent space
indirectly, we do not assume inherent similarities or differences across the
sequences. Instead, we directly constrain the geometry of the representations
themselves to be low-rank. We motivate our work from the perspective of low-
rank matrix recovery [24], common in other areas of machine learning and which
serves as the foundation for principal component analysis [43], classical
methods in natural language processing (topic modeling) [44, 45] and
collaborative filtering [46]. Problems from this area often admit convex
formulations and come with provable recovery guarantees [28]. Most recently, a
line of work has loosened the structural assumptions needed for signal (time-
series) recovery in an autoregressive context while still maintaining problem
convexity, using VIs with monotone operators the main tool [47, 48, 49].
However, to our knowledge, we are the first to apply these ideas in
conjunction with low-rank matrix recovery for modeling and representing
collections of general time-series observations.

##  2 Problem Setup

We aim to representing observations into $N$ vector-valued time series of
length $T$ each of the form $\\{\mathbf{x}_{i,t}\\}$, where
$\mathbf{x}\in\mathbb{R}^{C}$, $t\in[T]$, and $i\in[N]$. The sequences are
sampled from a common domain independently of each other across $i$, but have
temporal dependence across $t$. We refer to the history of events for sequence
$i$ up to time-point $t$ as $\mathcal{H}_{i,t}:=\\{\mathbf{x}_{i,s}\mid
s<t\\}$. We expect the behavior at event $\mathbf{x}_{i,t}$ to be a function
of past observations. In particular, for time-point $\mathbf{x}_{i,t}$, we
suppose the dependence on $\mathcal{H}_{i,t}$ is sufficiently captured an
order $d$ nonlinear vector autoregressive model with $C$ channels on
$\bm{\xi}_{i,t}=\operatorname{vec}(1,\\{\mathbf{x}_{i,t-s}\\}_{s=1}^{d})\in\mathbb{R}^{Cd+1}$,
the values from the $d$ preceding observations and a bias term so that

$\mathbb{E}[\mathbf{x}_{i,t}\mid\mathcal{H}_{t}]=\eta(\mathbf{R}_{i}\bm{\xi}_{i,t}).$
(1)

The matrices $\mathbf{R}_{i}\in\mathbb{R}^{C\times(Cd+1)}$ each serve as
weights for the prediction of the focal observation $\mathbf{x}_{i,t}$ that we
aim to learn. We allow
$\mathbf{b}_{i}=\operatorname{vec}({\mathbf{R}_{i}})\in\mathbb{R}^{C^{2}d+C}$
to be the parameters corresponding to the $i^{\text{th}}$ sequence arranged as
a vector and sufficient to capture the dynamics of the $i^{\text{th}}$ time-
series. The function $\eta:\mathbb{R}^{C}\to\mathbb{R}^{C}$ is a monotone link
function. The choice of link function $\eta$ naturally corresponds to the
character of the recovered sequence dynamics, which we illustrate via the
following examples:

Vector auto-regression

$\eta(\mathbf{x})=\mathbf{x};\mathbf{x}\in\mathbb{R}^{n}$, e.g. motion
capture, electrocardiogram (ECG)signals.

Symbolic sequences

$\eta(\mathbf{x})=\exp({\mathbf{x}})/\sum_{i}\exp(x_{i});\mathbf{x}\in[\Sigma]^{n}$,
e.g. natural language, genes.

Count processes

$\eta(\mathbf{x})=\exp(\mathbf{x});\mathbf{x}\in\mathbb{Z}_{\geq 0}^{n}$, e.g.
traffic intensity, call center arrival rates.

Bernoulli Processes

$\eta(\mathbf{x})=\exp(\mathbf{x})/(1+\exp(\mathbf{x}));\mathbf{x}\in\mathbb{B}^{n}$,
e.g. wildfire presence, neuron firing.

We do not restrict the mechanics of the link function $\eta$ beyond the
monotone property. We also remark that each vector $\mathbf{b}_{i}$
corresponding to each sequence may itself be high dimensional. The key aspect
of our method for low dimensional representation learning lies in the common
domain assumption, which should limit how the sequences are similar
(different). We leverage this information by a low rank assumption on the
space of parameters by which each sequence is described. In this way, we
constrain the individual $\mathbf{b}_{i}$ to lie approximately on a low
dimensional subspace of the possible full parameter space of
$\mathbf{b}_{i}$s. The representation of each sequence within this subspace
may be taken as an embedding and used for downstream tasks such as clustering,
classification, and anomaly detection.

In particular, we consider the autoregressive sequence model introduced in
(1), allowing $\mathbf{b}_{i}$ to be those parameters unique to the $i$th
sequence. Allow the matrix
$\mathbf{B}=\begin{bmatrix}\mathbf{b}_{1}&\ldots&\mathbf{b}_{N}\end{bmatrix}\in\mathbb{R}^{m\times
N}$ denote the parameters across all the sequences. We aim to recover a good
choice of the matrix $\mathbf{B}$ without supervision. We should balance two
goals: (1) we desire each $\mathbf{b}_{i}$ to be as faithful to the generating
dynamics of their respective observed data as possible; (2) we hope to
leverage the common domain assumption about the sequences and use information
from the other sequences to inform the prediction of the focal sequence. To
express the corresponding low-rank constraint, consider the rank $r$ Singular
Value Decomposition (SVD) of $\mathbf{B}$

$\mathbf{B}=\mathbf{U}\mathbf{\Sigma}\mathbf{V^{*}}=\sum_{k=1}^{r}\sigma_{k}\mathbf{u}_{k}\mathbf{v}_{k}^{*}$
(2)

where $\mathbf{\Sigma}=\operatorname{diag}\\{\sigma_{k}\\}_{k=1}^{r}$
corresponds to the singular values, columns of
$\mathbf{U}=[\mathbf{u}_{k}]_{k=1}^{r}\in\mathbb{R}^{m\times r}$ form an
orthobasis in $\mathbb{R}^{d}$, and columns of
$\mathbf{V}^{*}=[\mathbf{v}_{k}]_{k=1}^{r}\in\mathbb{R}^{r\times N}$ form an
orthobasis in $\mathbb{R}^{N}$. The recovered columns
$\mathbf{C}:=\mathbf{\Sigma}\mathbf{V}^{*}=[\mathbf{c}_{i}]_{i=1}^{N}\in\mathbb{R}^{r\times
N}$ give an $r$-dimensional representation for each of the $N$ sequences.
Likewise, the geometry of the subspace
$\operatorname{sp}\\{\mathbf{u}_{k}\\}_{k=1}^{r}$ describes the common domain
from which the generating processes of the sequences arise. We consider the
low dimensional representation for the $i$th sequence
$\mathbf{b}_{i}=\mathbf{U}\mathbf{c}_{i}$ as the embedding of the dynamics for
the $i^{\text{th}}$ sequence.

Because rank-constrained optimization is, in general, an NP-hard problem [50],
to enforce the low-rank requirement on $\mathbf{B}$, we instead constrain our
setup to a nuclear norm ball. The nuclear norm is given by
$\|\mathbf{X}\|_{*}=\sum_{j=1}^{r}\sigma_{i}(\mathbf{X})$ where $\sigma_{i}$
is the $i^{\text{th}}$ singular value of the matrix $\mathbf{X}$. The nuclear
norm is also the tightest convex relaxation to matrix rank [51] leading to
tractable parameter recovery and allows us to leverage a long line of work
from convex optimization and matrix recovery [52, 24, 53].

In terms of the modeling of sequences themselves, we treat each sequence as
arising from an individual stochastic source from which we acquire
observations and aim to learn a sequence-level autoregressive model. We cast
this parameter recovery problem as a monotone variational inequality, the most
general form of convex problem for which there exist efficient solution
methods [28, 54, 49, 48]. The objective is designed so that the dynamics
should be faithful in the sense that, for a focal sequence, the learned
dynamics can be used to generate sequences similar to the focal time-series
observations while still maintaining problem convexity.

##  3 Method

In the following, we present our method, first for linear auto-regressive
models and then for general non-linear auto-regressive models, which can be
applied to categorical sequences.

###  3.1 Low rank time-series embedding for linear auto-regressive models

First, suppose events $\mathbf{x}_{i,t}\in\mathbb{R}^{n}$ obey a linear
autoregressive model. Recall that we allow
$\bm{\xi}_{i,t}=\operatorname{vec}(1,\\{\mathbf{x}_{i,t-s}\\}_{s=1}^{d})\in\mathbb{R}^{Cd+1}$,
the values from the $d$ preceding observations and bias term so that

$\mathbb{E}[\mathbf{x}_{i,t}|\mathcal{H}_{t}]=\mathbf{R}_{i}\bm{\xi}_{i,t}.$
(3)

The matrices $\mathbf{R}_{i}\in\mathbb{R}^{C\times(Cd+1)}$ serve as weights
for the prediction of the focal observation $\mathbf{x}_{i,t}$. Recall also
that we allow
$\mathbf{b}_{i}=\operatorname{vec}({\mathbf{R}_{i}})\in\mathbb{R}^{C^{2}d+C}$.
To recover the parameter matrix
$\mathbf{B}=[\mathbf{b}_{i}]_{i=1}^{N}\in\mathbb{R}^{(C^{2}d+C)\times N}$, a
natural choice is to take least squares loss and write

$\displaystyle\min_{\widehat{\mathbf{B}}\in\mathbb{R}^{d\times
N}}\frac{1}{N}\sum_{i=1}^{N}\left(\frac{1}{T-d}\sum_{t=d+1}^{T}\|\mathbf{x}_{i,t}-\widehat{\mathbf{R}}_{i}\bm{\xi}_{i,t}\|_{2}^{2}\right)\qquad\text{
s.t. }\qquad\|{\widehat{\mathbf{B}}}\|_{*}\leq\lambda.$ (4)

##### Low-rank recovery and the nuclear norm regularization.

We now discuss Program (4) in the context of low-rank matrix recovery [24]. We
aim to recover matrix $\mathbf{B}$, but instead of observing it directly, we
receive indirect samples $\mathbf{y}\approx\mathcal{A}(\mathbf{B})$ through
random linear measurement operator $\mathcal{A}:\mathbb{R}^{(C^{2}d+C)\times
N}\to\mathbb{R}^{CN}$ such that
$\mathbb{E}[\mathbf{y}|\mathcal{A}]=\mathcal{A}(\mathbf{y})$. Namely, we
realize samples of form $(\mathcal{A}_{t},\mathbf{y}_{t})$ by selecting a
time-point $t$ across the observation horizon. The samples
$\mathbf{y}_{t}:=\operatorname{vec}([\mathbf{x}_{i,t}]_{i=1}^{N})\in\mathbb{R}^{CN}$
represent the values from the present time-point across all $N$ sequences. The
sampled operator $\mathcal{A}_{t}$ packages together the preceding length $d$
window of regressors across the $N$ sequences and adds also the corresponding
random noise $\bm{\epsilon}$ so that

$\mathcal{A}_{t}(\widehat{\mathbf{B}}):=\frac{1}{N}\operatorname{vec}([\mathbf{R}_{i}\bm{\xi}_{i,t}+\bm{\epsilon}]_{i=1}^{N}):\mathbb{R}^{(C^{2}d+C)\times
N}\to\mathbb{R}^{CN}.$ (5)

We take noisy indirect observations $\mathbf{y}\approx\mathcal{A}(\mathbf{B})$
to be drawn from $(\mathcal{A},\mathbf{y})\sim P$ and consider the expected
least squares loss via stochastic program

$\min_{\widehat{\mathbf{B}}}\ell(\widehat{\mathbf{B}}):=\mathbb{E}_{(\mathcal{A},\mathbf{y})\sim
P}\|\mathcal{A}(\widehat{\mathbf{B}})-\mathbf{y}\|_{2}^{2}\qquad\text{ s.t.
}\qquad\|\widehat{\mathbf{B}}\|_{*}\leq\lambda.$ (6)

Since we have only access to the observed temporal slices of size $d+1$
running up to time $T$, we now build the empirical analog to Program (6),

$\min_{\widehat{\mathbf{B}}}\hat{\ell}(\widehat{\mathbf{B}}):=\frac{1}{T-d}\sum_{t=d+1}^{T}\|\mathcal{A}_{t}(\widehat{\mathbf{B}})-\mathbf{y}_{t}\|_{2}^{2}\qquad\text{
s.t. }\qquad\|\widehat{\mathbf{B}}\|_{*}\leq\lambda,$ (7)

which is a Lipschitz smooth convex program on the nuclear ball of radius
$\lambda$ [55]. Program (7) is exactly the same as Program (4) except placed
in a matrix recovery context. We aim to recover the optimal $\mathbf{B}$ from
the samples while accounting for the global structure. When $\lambda$ is
arbitrarily large, there is no constraint on $\widehat{\mathbf{B}}$ and
Program (7) corresponds to fitting each sequence individually with no global
information. On the other extreme, forcing $\widehat{\mathbf{B}}$ to be rank
one constrains the model for each sequence to be multiples of each other. The
intermediate values of $\lambda$ correspond to various trade-offs between
learning the common global structure and attention to the individual
sequences.

Program (7) can be readily cast as a Semidefinite Program (SDP) solvable via
interior point methods [56]. However, as the size of $\mathbf{B}$ may reach
into the hundreds of thousands of decision variables, we turn our discussion
instead to solutions via efficient first-order proximal algorithms [57, 58]
analogous to those for linear inverse problems [59]. We will now describe the
proximal setup for nuclear norm minimization in the context of time series
embedding [53]. To illustrate, consider the following Mirror Descent (MD)
procedure

$\mathbf{B}_{k+1}=\operatorname{Prox}_{\mathbf{B}_{k}}(\gamma_{k}\nabla_{\mathbf{B}_{k}}[\ell(\mathbf{B}_{k})]),\qquad\mathbf{B}_{0}\in\\{\mathbf{X}\mid\|\mathbf{X}\|_{*}\leq\lambda\\},$
(8)

for an appropriately chosen sequence of steps $\\{\gamma_{k}\\}$. The solution
at step $k$ is given by the aggregate
$\tilde{\mathbf{B}}_{k}=(\sum_{\tau=1}^{k}\gamma_{\tau})^{-1}\sum_{\tau=1}^{k}\gamma_{\tau}\mathbf{B}_{\tau}$.
We take prox-mapping
$\operatorname{Prox}_{\mathbf{Z}}(\mathbf{X})=\operatorname*{arg\,min}_{\|\mathbf{Y}\|_{*}\leq\lambda}\omega(\mathbf{Y})+\langle\mathbf{X}-\nabla_{\mathbf{Z}}[\omega(\mathbf{Z})],\mathbf{Y}\rangle$
and Bregman divergence $\omega$ to be ones associated with the nuclear ball.
In particular, we can compute the nuclear-norm prox mapping by Singular Value
Thresholding (SVT), eliminating the small singular values at every step [52,
53]. Namely, with $m:=\min(d,N)$; $q:=\frac{1}{2\ln(2m)}$;
$\alpha:=\frac{4\sqrt{e}\log(2m)}{2^{q}(1+q)};\mathbf{M}=\mathbf{U}\operatorname{diag}\\{\sigma_{i}\\}_{i=1}^{m}\mathbf{V}^{*}$
we have the Bregman divergence and its subgradient as

$\omega(\mathbf{M})=\alpha\sum_{i=1}^{m}\sigma_{i}^{1+q}\implies\partial_{\mathbf{M}}[\omega(\mathbf{M})]=\sum_{i=1}^{m}[\alpha(1+q)\sigma_{i}^{q}]\mathbf{u}_{i}\mathbf{v}_{i}^{*}.$
(9)

To compute the prox-mapping, consider the SVD of
$\mathbf{X}-\partial_{\mathbf{Z}}[\omega(\mathbf{Z})]=\mathbf{P}\operatorname{diag}\\{\sigma_{k}\\}_{k=1}^{r}\mathbf{Q^{*}}$.
The optimal value of the linear program

$\mathbf{t}=\min_{\mathbf{t}\in\mathbb{R}^{n}}\\{\sum_{i=1}^{m}\frac{1}{2}t_{j}^{2}-\sigma_{i}t_{i}\mid\mathbf{t}\geq
0,\sum_{j=1}^{m}t_{j}\leq\lambda\\}$ (10)

gives
$\operatorname{Prox}_{\mathbf{X}}(\mathbf{B})=\mathbf{U}\operatorname{diag}\\{-\mathbf{t}\\}_{t=1}^{m}\mathbf{V}^{*}$
as the prox-mapping associated with the nuclear ball [53]. Notice that Linear
Program (10) can be solved in time $\mathcal{O}(m)$ in the worst case. The
time for each iteration of Mirror Descent (8) is dominated by the cost of the
SVD and the cost to compute the gradient of the loss function.

###  3.2 Nonlinear time-series embedding by monotone VI

We will now extend our discussion to the nonlinear case. Consider again the
events $\mathbf{x}_{i,t}\in\mathbb{R}^{C}$, and allow
$\bm{\xi}_{i,t}=\operatorname{vec}(1,\\{\mathbf{x}_{i,t-s}\\}_{s=1}^{d})\in\mathbb{R}^{Cd+1}$
to encode the past $d$ observations with bias term. Allow
$\eta:\mathbb{R}^{C}\to\mathbb{R}^{C}$ to be a fixed monotone link function.
Then consider the observation model

$\mathbb{E}[\mathbf{x}_{i,t}|\mathcal{H}_{t}]=\eta(\mathbf{R}_{i}\bm{\xi}_{i,t}).$
(11)

Our goal is to form a rank-constrained stochastic estimate to
$\mathbf{B}=[\operatorname{vec}(\mathbf{R}_{i})]_{i=1}^{N}\in\mathbb{R}^{(C^{2}d+C)\times
N}$. However, with arbitrary monotone link function, the Least Squares (LS)
approach outlined in (4) and (7) loses convexity and computational
tractability in general. Likewise, Maximum Likelihood Estimation (MLE) based
parameter estimation also becomes computationally difficult [28]. By contrast,
we shall cast the parameter recovery problem into a monotone VI formulation,
the most general type of convex program where there are known methods to
efficiently find high accuracy solutions [47, 49, 48].

##### Preliminaries on monotone VI

A monotone vector field on $\mathbb{R}^{m}$ with modulus of convexity $\beta$
is a vector field $G:\mathbb{R}^{m}\to\mathbb{R}^{m}$ such that

$\langle
G(\mathbf{x})-G(\mathbf{x}^{\prime}),\mathbf{x}-\mathbf{x}^{\prime}\rangle\geq\beta\|\mathbf{x}-\mathbf{x}^{\prime}\|\qquad\forall\
\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}$

when $\beta>0$, $G$ is strongly monotone. For some convex compact set
$\mathcal{X}\subseteq\mathbb{R}^{m}$, a point $\mathbf{x}^{*}$ is a weak
solution to the VI associated with $(G,\mathcal{X})$ if for all
$\mathbf{x}\in\mathcal{X}$ we have $\langle
G(\mathbf{x}),\mathbf{x}-\mathbf{x}^{*}\rangle\geq 0$. If $G$ is strongly
monotone and a weak solution exists, then the solution is unique. When
$\langle G(\mathbf{x}^{*}),\mathbf{x}-\mathbf{x}^{*}\rangle\geq 0$ for all
$\mathbf{x}\in\mathcal{X}$, we term $\mathbf{x}^{*}$ a strong solution to the
VI. When $G$ is continuous on $\mathcal{X}$, all strong solutions are weak
solutions and vice versa.

##### Monotone VI for nonlinear parameter recovery

We now turn our attention to the construction of a Monotone VI, which has as
its root optimal parameters corresponding to Model (11). We will use the same
operator $\mathcal{A}$ from the linear case from (5) together with its
associated adjoint $\mathcal{A}^{*}$. Recall that $\mathcal{A}$ is a random
operator that draws upon the concrete time-dependent samples $\mathcal{A}_{t}$
(which is coupled with observation $\mathbf{y}_{t}$ of values for the focal
time-point across all observed sequences) for some random choice of $t$. The
corresponding adjoint
$\mathcal{A}_{t}^{*}:\mathbb{R}^{CN}\to\mathbb{R}^{(C^{2}d+C)\times N}$ takes
the pre-image of the multichannel predictions (observations)
$\mathbf{y}=\operatorname{vec}([\mathbf{x}_{i}]_{i=1}^{N})$ and maps them back
to the parameter space, and may be computed using the below formula

$\mathcal{A}_{t}^{*}(\mathbf{y})=\frac{1}{N}[\operatorname{vec}([x_{i,c}\mathbf{R}_{i}^{T}\mathbf{e}_{c}])_{c=1}^{C})]_{i=1}^{N}:\mathbb{R}^{CN}\to\mathbb{R}^{(C^{2}d+C)\times
N}$ (12)

where $\mathbf{e}_{c}$ is $c^{\text{th}}$ standard basis. The adjoint, for
each entry and channel, multiplies the parameters by the value of the
observation corresponding to the channel.

We consider again the accompanying noisy observations $\mathbf{y}$ such that
in expectation
$\mathbb{E}[\mathbf{y}|\mathcal{A}]=\eta(\mathcal{A}(\mathbf{B}))$ as
discussed in Section 3.1 and where we extend the link function $\eta$ acts
sample wise. Consider now the vector field on the space of matrices

$\Psi(\widehat{\mathbf{B}})=\mathbb{E}_{(\mathcal{A},\mathbf{y})\sim
P}[\mathcal{A}^{*}(\eta(\mathcal{A}(\widehat{\mathbf{B}}))-\mathbf{y})]:\mathbb{R}^{(C^{2}d+C)\times
N}\to\mathbb{R}^{(C^{2}d+C)\times N}$ (13)

and notice that the matrix $\mathbf{B}$ of true generating parameters is a
zero of $\Psi$,

$\displaystyle\Psi({\mathbf{B}})$
$\displaystyle=\mathbb{E}_{(\mathcal{A},\mathbf{y})\sim
P}[\mathcal{A}^{*}(\eta(\mathcal{A}({\mathbf{B}}))-\mathbf{y})]=\mathbb{E}_{(\mathcal{A},\mathbf{y})\sim
P}[\mathcal{A}^{*}(\eta(\mathcal{A}({\mathbf{B}})))-\mathcal{A}^{*}(\mathbf{y})]$
$\displaystyle=\mathbb{E}_{(\mathcal{A},\mathbf{y})\sim
P}[\mathcal{A}^{*}(\eta(\mathcal{A}({\mathbf{B}})))-\mathcal{A}^{*}(\mathbb{E}[{\mathbf{y}\mid\mathcal{A}}])]$
$\displaystyle=\mathbb{E}_{(\mathcal{A},\mathbf{y})\sim
P}[\mathcal{A}^{*}(\eta(\mathcal{A}({\mathbf{B}})))-\mathcal{A}^{*}(\eta(\mathcal{A}({\mathbf{B}})))]=0.$

Since we have only access to the given observations, we take solutions to the
empirical version of the VI that takes slices from the time-series

$\widehat{\Psi}(\widehat{\mathbf{B}})=\frac{1}{T-d}\sum_{t=d}^{T}[\mathcal{A}_{t}^{*}(\eta(\mathcal{A}_{t}(\widehat{\mathbf{B}})))-\mathcal{A}_{k}^{*}(y_{k})]=\frac{1}{T-d}\sum_{t=d+1}^{T}\mathcal{A}_{t}^{*}[\eta(\mathcal{A}_{t}(\widehat{\mathbf{B}}))-\mathbf{y}_{t}].$
(14)

At each time window $t$, $\mathbf{y}_{t}$ represents the sequence
observations. $\mathcal{A}_{t}$ then takes in as input an estimate to matrix
$\mathbf{B}$ and uses the data from the previous $d$ time-points to output a
prediction for $\mathbf{y}_{t}$. We note that both $\mathcal{A}_{t}$ and
$\mathcal{A}_{t}^{*}$ may be computed in time $\mathcal{O}(NC^{2}d)$. In most
of our computations, we form an approximation to $\mathcal{A}$ by averaging
across the entire time horizon, giving a cost of $\mathcal{O}(TNC^{2}d)$. We
also illustrate averaging using smaller random sub-windows of the data in
Section 4.3. Analogous to Program (7), the VI associated with (14) likewise
admits solutions by MD. To illustrate, consider the recurrence:

${\mathbf{B}}_{k+1}=\operatorname{Prox}_{{\mathbf{B}}_{k}}(\gamma_{k}\Psi(\widehat{\mathbf{B}}_{k}))\qquad\mathbf{B}_{0}\in\\{{\mathbf{X}}\mid\|{\mathbf{X}}\|_{*}\leq\lambda\\}$
(15)

with step sizes $\\{\gamma_{k}\\}$. The aggregate solution at step $k$ is
given by
$\tilde{\mathbf{B}}_{k}=(\sum_{\tau=1}^{k}\gamma_{\tau})^{-1}\sum_{\tau=1}^{t}\gamma_{\tau}{\mathbf{B}}_{\tau}$
[60]. Note if $\eta:=\mathbf{Id}$, the identity function, then the VI
corresponds exactly to the gradient field of Program (7). In this case,
$\Psi({\mathbf{X}})=\nabla_{\mathbf{X}}[\ell(\mathbf{X})]$ and the MD
procedure for VI and LS are the same.

##### First order methods for monotone VI

To concretely solve the monotone VIs outlined in (13) and (14), we detail an
accelerated mirror-prox scheme with backtracking for nuclear norm constrained
VI in Algorithm 1 of Appendix A, which addresses the following general problem

$\langle\Psi(\mathbf{B}),\mathbf{B}-\mathbf{B}^{*}\rangle\geq
0\qquad\forall\mathbf{B}\in\mathcal{X}:=\\{\mathbf{B}\mid\|\mathbf{B}\|_{*}\leq\lambda\\}$
(16)

where $\Psi$ is an (unbiased estimator of) $\kappa$-lipschitz continuous
monotone vector field [61, 53], which addresses the difficulty that $\kappa$
is in most cases is unknown to us beforehand. The convergence results for this
class of algorithm are typical and presented in [53, 61, 60, 62]. Namely, for
$\epsilon$ error as measured by
$\epsilon_{\Psi,\mathcal{X}}(\tilde{\mathbf{B}}_{t})=\sup_{\mathbf{Z}\in\mathcal{X}}\langle\Psi(\mathbf{Z},\tilde{\mathbf{B}}_{t}-\mathbf{Z})\rangle$
solution requires $\mathcal{O}({\kappa}/{\epsilon})$ iterations for
deterministic VIs and
$\mathcal{O}(\kappa/{\epsilon}+\sigma^{2}/{\epsilon}^{2})$ for
$\mathbb{E}[\|\widehat{\Psi}(\mathbf{B})-\Psi(\mathbf{B})\|^{2}]\leq\sigma^{2}$
where $\widehat{\Psi}(\mathbf{B})$ is a stochastic approximation to the true
field $\Psi$ with bounded variance. Each iteration requires $\mathcal{O}(1)$
evaluations of the monotone operator $\Psi$ and $\operatorname{Prox}$
operator, which in turn is dominated by the cost of the SVD at each step
$\mathcal{O}(NC^{2}d\min(N,C^{2}d))$, with the trade-off being in the variance
of the approximation to the true VI field. The convergence of the algorithm as
applied to the vector fields given in (13) and (14) in particular may be
established similarly to [49], and when the data follow the true model (11),
parameter recovery guarantee can be established similarly to [48].

##### Parameter recovery for symbolic sequences

As a special case of (13), consider now that each channel represents the
probability of emitting a token taken from syllabary $\\{s_{c}\\}_{c=1}^{C}$
of size $C$. Then each $\mathbf{x}_{i,t}$ represents a probability vector
$\sum_{c}^{C}x_{i,t,c}=1$ and we have
$\mathbb{E}[x_{i,t,c}|\mathcal{H}_{t}]=\mathbb{P}[x_{i,t,c}=s_{c}]$. We take
the softmax activation function
$\sigma(\mathbf{y})=\operatorname{vec}([\|\exp(\mathbf{y}^{(i)})\|_{1}^{-1}\exp(\mathbf{y}^{(i)})]_{i=1}^{N})$,
where $\mathbf{y}^{(i)}$ corresponds to values from the $i^{\mathrm{th}}$
sequence. This problem corresponds to learning representations for different
sequences. We illustrate in Section 4.3 the above as applied to learning
dynamics for genomics data and natural language, for which autoregressive
models have become increasingly popular.

##  4 Experiments

We first illustrate parameter recovery using synthetic univariate time-series
in Section 4.1. We investigate the choice of nuclear penalty $\lambda$ and the
geometry of the recovered parameter matrix. Section 4.2 describes benchmarks
using real-world time-series data from the UCR Time Series Classification
Archive [63]. We report classification and runtime performance against a
number of baselines. Section 4.3 provides two illustrations on the embedding
of real-world sequence data. In the first example, we consider a language
representation task where we embed without supervision a series of excerpts
taken either from the works of Lewis Carroll or abstracts scraped from arXiv
[64, 65, 66]. In the second illustration, we apply our method to the
clustering of gene sequences for strains of the Influenza A and Dengue viruses
[67].

We evaluated all experiments and illustrations using a cluster with 24 core
Intel Xeon Gold 6226 CPU (2.7 GHZ) processors, and NVIDIA Tesla V100 Graphics
coprocessors (16 GB VRAM), and 384 GB of RAM. However, we could reproduce the
results on a standard personal computer. We implement the routines described
in Sections 3.1 and 3.2 and Appendix A in the Julia programming language. We
describe the experimental setup and results in detail in Appendix B.

###  4.1 Parameter recovery with synthetic autoregressive sequences

To illustrate parameter recovery across autoregressive sequences, we
synthetically generated a set of ten baseline parameters for linear
autoregressive sequences of order $d=15$. Within each class, we then created
$N=300$ sequences of each type, perturbing the baseline coefficients by adding
a small amount of noise according to a fixed rule for the set of parameters.
We then generated $T=250$ observations for each sequence according to the
autoregressive model in (3). We formed all 120 combinations of $k=3$ type of
sequences from the ten classes and recovered the underlying parameter matrix
by solving the program given in (7) to optimality. We report the data-
generating procedure and experimental details in Appendix B.1.

Table 1 reports averages and standard deviations for the relative
reconstruction error
$\|\mathbf{B}-\widehat{\mathbf{B}}\|_{F}/\|\mathbf{B}\|_{F}$, the least
squares error of the objective function given in (4), the Adjusted Rand Index
(ARI) using k-means clustering with $k=3$ clusters in the embedding space
across the runs, and the number of large singular values (singular values
within $10^{-2}$ of the principal singular value). We first learn
representations for sequences without nuclear norm constraint.

To illustrate the low-rank matrix recovery, we search across the values of
nuclear regularization $\lambda$ via Brent search [68] and report the
performance for a close to the optimal value of $\lambda$ with respect to the
reconstruction error. When the underlying dynamics share a common low-rank
structure, the nuclear constraint effectively leverages the common information
shared across the sequences to recover the true parameters more faithfully
using the common domain information with comparable error in the least squares
sense. Furthermore, we observe the nuclear regularization procedure driving
the singular value spectrum to sparsity, with the number of large singular
values being much smaller than in the unconstrained case.

To further illustrate, Figure 1 depicts parameter recovery for a collection of
sequences with autoregressive order $d=15$. In the leftmost pane, we report
the relative reconstruction error and the number of large singular values
across different values of nuclear constraint $\lambda$. In the central pane,
we depict the singular value spectra of the true parameter matrix (which has
approximately low-rank structure plus noise) and recovered matrices with
differing numbers of large singular values. In the third pane, we show the
first two principal components (total explained variance $=0.868$) of the
sequence embedding with the smallest reconstruction error [69] along with the
original three generating classes of the data. With the introduction of
sufficient nuclear regularization, we drastically reduced the reconstruction
error of our recovered solution and observed that the solutions with low
reconstruction error were of approximately low rank. In the plots of the
singular value spectra and the projection of the learned sequence embeddings,
we observe that those rank-constrained recoveries effectively recover those
large singular values in the spectrum of the true parameter matrix. By
contrast, the parameter recoveries performed without or insufficient nuclear
constraint fit the noise component of the data, as evident in the distribution
of singular values.

Table 1: Time series parameter recovery for synthetic autoregressive time-series. | Relative Err. | LS Err. |  Cluster ARI [70] |  $\approx$ Rank
---|---|---|---|---
$\lambda$ Selection | Avg. | Std. | Avg. | Std. | Avg. | Std. | Avg. | Std.
Unconstrained | 0.341 | (0.016) | 79.333 | (0.169) | 0.967 | (0.049) | 14.442 | (1.203)
$\operatorname*{arg\,min}_{\lambda}\|\widehat{\mathbf{B}}_{\lambda}-\mathbf{B}\|_{F}$ | 0.158 | (0.020) | 80.709 | (0.233) | 0.997 | (0.008) | 7.392 | (4.255)
Figure 1: Simulation results: Parameter recovery for a collection of
univariate time series drawn from $k=3$ classes. Left: Relative reconstruction
error and approximate rank of recovered parameter matrices across levels of
nuclear constraint $\lambda$. Center: Singular values of the true parameter
matrix $\mathbf{B}$ and the singular values of the recovered solutions of
varying dimensions. Right: First two principal components of the recovered
matrix with the smallest reconstruction error and original class labels.

###  4.2 Real time-series classification

Following in the line of [71, 40, 36, 72], we conduct real data experiments
using 36 of the UCR time series datasets [63]. We report statistics for the
time series datasets in Appendix B.2.1. Each dataset has a default test/train
split. We first re-encoded each of the time series as a multichannel signal
consisting of the signal itself and its first finite-difference. We then embed
all the data without supervision by solving (14) using the adaptive mirror
prox scheme described in Algorithm 1 of Appendix A with a look-back length of
$d=20$ and running the algorithm for $256$ steps, taking the link function to
be linear. To pick the value of $\lambda$, we find a value when the solution
becomes rank one via bisection and perform a grid search across values of
$\lambda$, which results in rank-constrained parameters. We then report our
results using the choice of $\lambda$ with the best performance on the
training set. We report the ARI [70], Normalized Mutual Information (NMI)
[73], macro-F1 score [74], and accuracy on the test set. We also report the
average runtime.

We compare our method with five representative time series embedding and
classification methods: K-nearest neighbors (KNN) using Euclidean ($\ell_{2}$)
distance [75], KNN with Dynamic Time Warping (DTW) as the distance metric
[31], shapeDTW (another method based on DTW but with additional features) [72,
29], a dictionary-based method MultiROCKET+Hydra [76], and one deep-learning
representation method based on contrastive learning (TS2Vec) [36]. In line
with [36, 77], to evaluate the classification performance on test set for
methods which produce embeddings (TS2Vec and our method), we perform cross-
validated grid search across KNNs with $k=\\{2^{i}\mid i\in[0,4]\\}$ neighbors
or SVMs with RBF kernels with penalty values $c\in\\{2^{i}\mid
i\in[-10,15]\\}\cup{\infty}$. We defer all further details of our experimental
setup to Appendix B. Table 2 displays the mean and standard deviation across
the metrics across the datasets. We provide the detailed results per dataset
in Tables 5 and 6 of Appendix B.2. We observe superior performance to baseline
methods based on distance metrics, such as Euclidean distance or DTW, and
observe comparable performance to TS2Vec. We note that for this class of
univariate sequence, the heuristic dictionary-based ensemble (MR-HYDRA)
outperforms both our approach and the deep-learning-based approaches. However,
this method has been tuned specifically for this type of classification
problem. By contrast, similar to [36], we consider classification only as one
potential downstream task and do not rely on handcrafted features.

Table 2: Time series classification performance on UCR time series data of our method vs a number of baselines (higher is better, except for runtime). We outperform simple approaches and perform close to classification using the embeddings from the neural network based TS2Vec but use only 37% of the runtime. The best performing method, MR-Hydra, is a ensemble based on handpicked features tuned specially for time series classification.  |  ARI [70] |  NMI [73] |  F1[74] | Accuracy | Runtime (Sec)
---|---|---|---|---|---
Method | Avg. | Std. | Avg. | Std. | Avg. | Std. | Avg. | Std. | Avg. | Std.
$\ell_{2}$+KNN [75] | 0.422 | (0.294) | 0.416 | (0.290) | 0.752 | (0.144) | 0.725 | (0.166) | 0.128 | (0.413)
DTW+KNN [31] | 0.447 | (0.303) | 0.435 | (0.300) | 0.766 | (0.150) | 0.738 | (0.170) | 42.988 | (134.320)
shapeDTW [78] | 0.470 | (0.307) | 0.460 | (0.299) | 0.773 | (0.152) | 0.746 | (0.178) | 21.871 | (71.655)
TS2Vec [36] | 0.606 | (0.282) | 0.580 | (0.287) | 0.840 | (0.138) | 0.814 | (0.178) | 1,085.092 | (1,408.458)
Ours | 0.602 | (0.282) | 0.562 | (0.293) | 0.817 | (0.180) | 0.788 | (0.193) | 400.031 | (677.486)
MR-Hydra [76] | 0.682 | (0.273) | 0.656 | (0.285) | 0.877 | (0.121) | 0.851 | (0.162) | 10.197 | (16.459)

###  4.3 Symbolic sequences: language and genomics

We now illustrate the capability of our method to learn meaningful
representations for sequences with nonlinear dynamics, namely autoregressive
sequence and language modeling.

##### Symbolic sequences and language: arXiv abstracts or “Alice in
Wonderland”?

We first consider an autoregressive language modeling task, drawing textual
sequences from three sources: two works by the same author Lewis Caroll —
Alice’s Adventures in Wonderland ($n=228$) [64] and Through the Looking Glass
($n=316$) [65] — and machine learning related abstracts scraped from ArXiv
($n=600$) [66] (details may be found in Appendix B.4). We embed the sequences
without supervision with a lookback of $d=75$, and in order to reduce the
number of symbols in our alphabet and avoid the blowup in the number of
channels, we converted each of the sequences into a $c=4$ symbol code via
Huffman coding, based on the overall frequencies of letters in the English
language [79]. We then solve Program (14) using the multichannel measurement
operator given in (5) to optimality and using the softmax activation discussed
in Section 3.2. We show in Figure 2(a) the space learned when $\lambda$ was
chosen to be sufficiently small as to give a rank three representation of the
data. We then project the learned representation via UMAP [80]. Two distinct
clusters form corresponding to the two different genres of writing, however,
whereas the paper abstracts are clearly separable from the works of Lewis
Caroll, the two books written by him are not as clearly disambiguable
considering they are from the same author.

##### Virus strain identification from genome sequences

For the final illustration, in line with [81], the problem of classifying
genetic sequences, which allows for the placing of species/variants in the
evolutionary context of others (phylogenetic classification). We consider gene
sequence data from segment six of the Influenza A virus genome ($n=949$,
$\text{average sequence length}=1409$) [82] and the complete genome the Dengue
virus ($n=1633$, $\text{average sequence length}=10559$) [83]. We consider
gene sequences from five strains of Influenza and the four strains of Dengue.
Likewise, we provide a detailed overview of the data and learning procedure in
Appendix B.5. We encode the genomes in a similar manner as for the natural
language illustration, assigning one channel to each nucleotide (A, C, T, G),
and encode the presence/absence of each nucleotide at each position via one-
hot encoding.

To recover the embedding, we adopt the same softmax activation scheme as
described in Section 3.2. Since the genomes are of variable length, we
consider a stochastic approximation to the monotone field $\Psi$ (13) by
taking the sample average of randomly selected length $G=800$ sub-windows from
each of sequences at each training step. We consider clustering the Influenza
and Dengue genome segments individually and report UMAP projections of the
learned representations in Figure 2(b) and 2(c), respectively. The dimensions
of the learned embeddings are 7 and 20, respectively. In these subspaces, we
note the clear grouping of viral strains obtained via solving the stochastic
approximation to the VI (13).

(a) Lewis Caroll or ArXiv abstracts? Genres form clusters.

(b) Clustering strains of Influenza A virus genome data (segment 6)

(c) Clustering using full genome for the four strains of Dengue Virus

Figure 2: Learned embeddings for symbolic sequences collections using our
method— visualised by UMAP projections shows clear groupings based on
sequences with similar underlying dynamics.

##  5 Discussion

We have introduced a method to learn low-dimensional representations for
collections of nonlinear time series with provable recovery guarantees. We
framed the problem as a low-rank matrix recovery problem that is representable
as a monotone variational inequality, the most general type of convex problem
for which efficient solution methods exist. We suppose the sequence dynamics
obey an autoregressive model via a monotone link function but make no further
assumptions about the sequences. This notably includes the important case of
the probabilistic modeling of symbolic sequences. Our method exhibits strong
performance when these assumptions are satisfied. However, the method is
limited when the data falls short of these assumptions, as is the case for
some UCR time-series data in Section 4.2. To address this, we note that we may
also design alternative objectives as long they are representable in the form
of monotone VI. The runtime performance of our method is dominated by (1) the
cost of the SVD to enforce the low-rank constraint and (2) computing a high-
quality estimate to the underlying VI by sample stochastic approximation.
However, as shown in the Section 4.3 genome representation example, we can
take random windows substantially shorter than the entire sequence to estimate
the VI field and still maintain competitive performance. Furthermore, if the
low-rank assumption is satisfied and the majority of the signal energies are
concentrated in the principal singular values, at each $\operatorname{Prox}$
(thresholding) step, we need only compute singular values up to $\lambda$. In
this setting, a good idea is to apply an iterative Lanczos method to find the
first few singular values/vectors in an on-line fashion [52, 84, 85]. However,
we discuss neither these acceleration methods nor loosening the class of the
model while still maintaining convexity in the scope of this work.

## References

  * [1] Mike X Cohen  “Analyzing Neural Time Series Data: Theory and Practice”  The MIT Press, 2014  DOI: 10.7551/mitpress/9609.001.0001
  * [2] Erick A Perez Alday et al.  “Classification of 12-lead ECGs: the PhysioNet/Computing in Cardiology Challenge 2020”  In _Physiological Measurement_ 41.12 IOP Publishing, 2020, pp. 124003  DOI: 10.1088/1361-6579/abc960
  * [3] J.J. Van Wijk and E.R. Van Selow  “Cluster and calendar based visualization of time series data”  In _Proceedings 1999 IEEE Symposium on Information Visualization (InfoVis’99)_ , 1999, pp. 4–9  DOI: 10.1109/INFVIS.1999.801851
  * [4] Ruey S. Tsay  “Analysis of financial time series”, Wiley series in probability and statistics  Hoboken, NJ: Wiley-Interscience, 2005  URL: http://gso.gbv.de/DB=2.1/CMD?ACT=SRCHA&SRT=YOP&IKT=1016&TRM=ppn+483463442&sourceid=fbw_bibsonomy
  * [5] Wei Min et al.  “Explainable Deep Behavioral Sequence Clustering for Transaction Fraud Detection”  In _CoRR_ abs/2101.04285, 2021  arXiv: https://arxiv.org/abs/2101.04285
  * [6] S. Basu, A. Mukherjee and S. Klivansky  “Time series models for internet traffic”  In _Proceedings of IEEE INFOCOM ’96. Conference on Computer Communications_ 2, 1996, pp. 611–620 vol.2  DOI: 10.1109/INFCOM.1996.493355
  * [7] Matthew A Reyna et al.  “Early prediction of sepsis from clinical data: The PhysioNet/Computing in Cardiology Challenge 2019”  In _Crit. Care Med._ 48.2 Ovid Technologies (Wolters Kluwer Health), 2020, pp. 210–217
  * [8] Laila Rasmy et al.  “Med-BERT: pretrained contextualized embeddings on large-scale structured electronic health records for disease prediction”  In _npj Digital Medicine_ 4.1, 2021, pp. 86  DOI: 10.1038/s41746-021-00455-y
  * [9] R. Argelaguet et al.  “MOFA+: a statistical framework for comprehensive integration of multi-modal single-cell data”  In _Genome Biology_ 21.1, 2020, pp. 111
  * [10] John Jumper et al.  “Highly accurate protein structure prediction with AlphaFold”  In _Nature_ 596.7873, 2021, pp. 583–589  DOI: 10.1038/s41586-021-03819-2
  * [11] Peter F. Brown et al.  “Class-Based n-gram Models of Natural Language”  In _Computational Linguistics_ 18.4, 1992, pp. 467–480  URL: https://aclanthology.org/J92-4003
  * [12] Matthew E. Peters et al.  “Deep Contextualized Word Representations”  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ New Orleans, Louisiana: Association for Computational Linguistics, 2018, pp. 2227–2237  DOI: 10.18653/v1/N18-1202
  * [13] Nils Reimers and Iryna Gurevych  “Sentence-BERT: Sentence Embeddings using Siamese BERT-Networks”  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ Hong Kong, China: Association for Computational Linguistics, 2019, pp. 3982–3992  DOI: 10.18653/v1/D19-1410
  * [14] Daniel Cer et al.  “Universal Sentence Encoder for English”  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ Brussels, Belgium: Association for Computational Linguistics, 2018, pp. 169–174  DOI: 10.18653/v1/D18-2029
  * [15] Tomas Mikolov et al.  “Distributed Representations of Words and Phrases and their Compositionality”  In _Advances in Neural Information Processing Systems_ 26 Curran Associates, Inc., 2013  URL: https://proceedings.neurips.cc/paper_files/paper/2013/file/9aa42b31882ec039965f3c4923ce901b-Paper.pdf
  * [16] Yoshua Bengio, Aaron Courville and Pascal Vincent  “Representation Learning: A Review and New Perspectives”  In _IEEE Trans. Pattern Anal. Mach. Intell._ 35.8 USA: IEEE Computer Society, 2013, pp. 1798–1828  DOI: 10.1109/TPAMI.2013.50
  * [17] Patara Trirat et al.  “Universal Time-Series Representation Learning: A Survey”, 2024  arXiv:2401.03717 [cs.LG]
  * [18] Ravid Shwartz Ziv and Yann LeCun  “To Compress or Not to Compress—Self-Supervised Learning and Information Theory: A Review”  In _Entropy_ 26.3, 2024  DOI: 10.3390/e26030252
  * [19] Kevin P. Murphy  “Probabilistic Machine Learning: Advanced Topics”  MIT Press, 2023  URL: http://probml.github.io/book2
  * [20] Diederik P Kingma and Max Welling  “Auto-Encoding Variational Bayes”, 2022  arXiv:1312.6114 [stat.ML]
  * [21] Baptiste Lafabregue, Jonathan Weber, Pierre Gançarski and Germain Forestier  “End-to-end deep representation learning for time series clustering: a comparative study”  In _Data Mining and Knowledge Discovery_ 36.1, 2022, pp. 29–81  DOI: 10.1007/s10618-021-00796-y
  * [22] Rayan Krishnan, Pranav Rajpurkar and Eric J. Topol  “Self-supervised learning in medicine and healthcare”  In _Nature Biomedical Engineering_ 6.12, 2022, pp. 1346–1352  DOI: 10.1038/s41551-022-00914-1
  * [23] Yinfei Yang et al.  “Multilingual Universal Sentence Encoder for Semantic Retrieval”  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics: System Demonstrations_ Online: Association for Computational Linguistics, 2020, pp. 87–94  DOI: 10.18653/v1/2020.acl-demos.12
  * [24] Mark A. Davenport and Justin Romberg  “An Overview of Low-Rank Matrix Recovery From Incomplete Observations”  In _IEEE Journal of Selected Topics in Signal Processing_ 10.4 Institute of ElectricalElectronics Engineers (IEEE), 2016, pp. 608–622  DOI: 10.1109/jstsp.2016.2539100
  * [25] Emmanuel J. Candes and Terence Tao  “The Power of Convex Relaxation: Near-Optimal Matrix Completion”  In _IEEE Transactions on Information Theory_ 56.5, 2010, pp. 2053–2080  DOI: 10.1109/TIT.2010.2044061
  * [26] Emmanuel J. Candes and Yaniv Plan  “Matrix Completion With Noise”  In _Proceedings of the IEEE_ 98.6, 2010, pp. 925–936  DOI: 10.1109/JPROC.2009.2035722
  * [27] Ali Ahmed and Justin Romberg  “Compressive Multiplexing of Correlated Signals”  In _IEEE Transactions on Information Theory_ 61.1, 2015, pp. 479–498  DOI: 10.1109/TIT.2014.2366459
  * [28] Anatoli Juditsky and Arkadi Nemirovski  “Statistical Inference via Convex Optimization”  Princeton University Press, 2020  URL: http://www.jstor.org/stable/j.ctvqsdxqd
  * [29] Lexiang Ye and Eamonn Keogh  “Time series shapelets: a new primitive for data mining”  In _Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’09  Paris, France: Association for Computing Machinery, 2009, pp. 947–956  DOI: 10.1145/1557019.1557122
  * [30] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest and Clifford Stein  “Introduction to Algorithms”  The MIT Press, 2001
  * [31] “Dynamic Time Warping”  In _Information Retrieval for Music and Motion_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 69–84  DOI: 10.1007/978-3-540-74048-3_4
  * [32] Anthony Bagnall et al.  “The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances”  In _Data Mining and Knowledge Discovery_ 31.3, 2017, pp. 606–660  DOI: 10.1007/s10618-016-0483-9
  * [33] Padhraic Smyth  “Clustering Sequences with Hidden Markov Models”  In _Advances in Neural Information Processing Systems_ 9 MIT Press, 1996  URL: https://proceedings.neurips.cc/paper_files/paper/1996/file/6a61d423d02a1c56250dc23ae7ff12f3-Paper.pdf
  * [34] K. Kalpakis, D. Gada and V. Puttagunta  “Distance measures for effective clustering of ARIMA time-series”  In _Proceedings 2001 IEEE International Conference on Data Mining_ , 2001, pp. 273–280  DOI: 10.1109/ICDM.2001.989529
  * [35] Ling Yang and Shenda Hong  “Unsupervised Time-Series Representation Learning with Iterative Bilinear Temporal-Spectral Fusion”  In _Proceedings of the 39th International Conference on Machine Learning_ 162, Proceedings of Machine Learning Research  PMLR, 2022, pp. 25038–25054  URL: https://proceedings.mlr.press/v162/yang22e.html
  * [36] Zhihan Yue et al.  “TS2Vec: Towards Universal Representation of Time Series”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 36.8, 2022, pp. 8980–8987  DOI: 10.1609/aaai.v36i8.20881
  * [37] Zhiwen Xiao et al.  “Deep Contrastive Representation Learning With Self-Distillation”  In _IEEE Transactions on Emerging Topics in Computational Intelligence_ 8.1, 2024, pp. 3–15  DOI: 10.1109/TETCI.2023.3304948
  * [38] Archibald Felix Fraikin, Adrien Bennetot and Stephanie Allassonniere  “T-Rep: Representation Learning for Time Series using Time-Embeddings”  In _The Twelfth International Conference on Learning Representations_ , 2024  URL: https://openreview.net/forum?id=3y2TfP966N
  * [39] Yihe Wang, Yu Han, Haishuai Wang and Xiang Zhang  “Contrast Everything: A Hierarchical Contrastive Framework for Medical Time-Series”  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023  URL: https://openreview.net/forum?id=sOQBHlCmzp
  * [40] Qianli Ma, Jiawei Zheng, Sen Li and Gary W Cottrell  “Learning Representations for Time Series Clustering”  In _Advances in Neural Information Processing Systems_ 32 Curran Associates, Inc., 2019  URL: https://proceedings.neurips.cc/paper_files/paper/2019/file/1359aa933b48b754a2f54adb688bfa77-Paper.pdf
  * [41] Vincent Fortuin, Dmitry Baranchuk, Gunnar Rätsch and Stephan Mandt  “Gp-vae: Deep probabilistic time series imputation”  In _International conference on artificial intelligence and statistics_ , 2020, pp. 1651–1661  PMLR
  * [42] Jacob Devlin, Ming-Wei Chang, Kenton Lee and Kristina Toutanova  “BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding”  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ Minneapolis, Minnesota: Association for Computational Linguistics, 2019, pp. 4171–4186  DOI: 10.18653/v1/N19-1423
  * [43] Harold Hotelling  “Analysis of a complex of statistical variables into principal components.”  In _Journal of educational psychology_ 24.6 Warwick & York, 1933, pp. 417
  * [44] David M. Blei, Andrew Y. Ng and Michael I. Jordan  “Latent dirichlet allocation”  In _J. Mach. Learn. Res._ 3.null JMLR.org, 2003, pp. 993–1022
  * [45] David M. Blei  “Probabilistic topic models”  In _Commun. ACM_ 55.4 New York, NY, USA: Association for Computing Machinery, 2012, pp. 77–84  DOI: 10.1145/2133806.2133826
  * [46] Yehuda Koren, Robert Bell and Chris Volinsky  “Matrix Factorization Techniques for Recommender Systems”  In _Computer_ 42.8, 2009, pp. 30–37  DOI: 10.1109/MC.2009.263
  * [47] Anatoli Juditsky, Arkadi Nemirovski, Yao Xie and Chen Xu  “Generalized generalized linear models: Convex estimation and online bounds”, 2023  arXiv:2304.13793 [stat.ME]
  * [48] Anatoli Juditsky, Arkadi Nemirovski, Liyan Xie and Yao Xie  “Convex Parameter Recovery for Interacting Marked Processes”  In _IEEE Journal on Selected Areas in Information Theory_ 1.3, 2020, pp. 799–813  DOI: 10.1109/JSAIT.2020.3040999
  * [49] A.. Juditsky and A.. Nemirovski  “Signal Recovery by Stochastic Optimization”  In _Automation and Remote Control_ 80.10, 2019, pp. 1878–1893  DOI: 10.1134/S0005117919100084
  * [50] B.. Natarajan  “Sparse Approximate Solutions to Linear Systems”  In _SIAM Journal on Computing_ 24.2, 1995, pp. 227–234  DOI: 10.1137/S0097539792240406
  * [51] Benjamin Recht, Maryam Fazel and Pablo A. Parrilo  “Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization”  In _SIAM Review_ 52.3, 2010, pp. 471–501  DOI: 10.1137/070697835
  * [52] Jian-Feng Cai, Emmanuel J. Candès and Zuowei Shen  “A Singular Value Thresholding Algorithm for Matrix Completion”  In _SIAM Journal on Optimization_ 20.4, 2010, pp. 1956–1982  DOI: 10.1137/080738970
  * [53] Yurii Nesterov and Arkadi Nemirovski  “On first-order algorithms for l1/nuclear norm minimization”  In _Acta Numerica_ 22, 2013, pp. 509–575  DOI: 10.1017/S096249291300007X
  * [54] Francisco Facchinei and Jong-Shi Pang  “Finite-dimensional variational inequalities and complementarity problems” 1, Springer series in operations research  New York: Springer, 2003
  * [55] Alexander Shapiro, Darinka Dentcheva and Andrzej Ruszczynski  “Lectures on Stochastic Programming: Modeling and Theory, Third Edition”  Philadelphia, PA: Society for IndustrialApplied Mathematics, 2021  DOI: 10.1137/1.9781611976595
  * [56] Aharon Ben-Tal and Arkadi Nemirovski  “Lectures on Modern Convex Optimization”, MOS-SIAM Series on Optimization  Society for IndustrialApplied Mathematics, 2001
  * [57] Patrick L. Combettes and Jean-Christophe Pesquet  “Proximal Splitting Methods in Signal Processing”  In _Fixed-Point Algorithms for Inverse Problems in Science and Engineering_ New York, NY: Springer New York, 2011, pp. 185–212  DOI: 10.1007/978-1-4419-9569-8_10
  * [58] Neal Parikh and Stephen Boyd  “Proximal Algorithms”  In _Found. Trends Optim._ 1.3 Hanover, MA, USA: Now Publishers Inc., 2014, pp. 127–239  DOI: 10.1561/2400000003
  * [59] Amir Beck and Marc Teboulle  “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems”  In _SIAM Journal on Imaging Sciences_ 2.1, 2009, pp. 183–202  DOI: 10.1137/080716542
  * [60] Arkadi Nemirovski  “Prox-Method with Rate of Convergence O(1/t) for Variational Inequalities with Lipschitz Continuous Monotone Operators and Smooth Convex-Concave Saddle Point Problems”  In _SIAM Journal on Optimization_ 15.1, 2004, pp. 229–251  DOI: 10.1137/S1052623403425629
  * [61] Yunmei Chen, Guanghui Lan and Yuyuan Ouyang  “Accelerated schemes for a class of variational inequalities”  In _Mathematical Programming_ 165.1, 2017, pp. 113–149  DOI: 10.1007/s10107-017-1161-4
  * [62] Anatoli Juditsky, Arkadi Nemirovski and Claire Tauvel  “Solving Variational Inequalities with Stochastic Mirror-Prox Algorithm”  In _Stochastic Systems_ 1.1, 2011, pp. 17–58  DOI: 10.1287/10-SSY011
  * [63] Hoang Anh Dau et al.  “The UCR Time Series Classification Archive” https://www.cs.ucr.edu/~eamonn/time_series_data_2018/, 2018
  * [64] Lewis Carroll  “Alice’s Adventures in Wonderland”  Macmillan Publishers, 1865
  * [65] Lewis Carroll  “Through the Looking-Glass, and What Alice Found There”  Macmillan Publishers, 1871
  * [66] Kaggle Team  “Leveraging ML to Fuel New Discoveries with the arXiv Dataset”, 2020  URL: https://medium.com/@kaggleteam/leveraging-ml-to-fuel-new-discoveries-with-the-arxiv-dataset-981a95bfe365
  * [67] Eric W Sayers et al.  “Database resources of the national center for biotechnology information”  In _Nucleic Acids Res._ 50.D1 Oxford University Press (OUP), 2022, pp. D20–D26
  * [68] R.. Brent  “Algorithms for minimization without derivatives /”  Mineola, N.Y. :: Dover Publications, 2002  URL: http://catdir.loc.gov/catdir/description/dover031/2001047459.html
  * [69] Kevin P. Murphy  “Probabilistic Machine Learning: An introduction”  MIT Press, 2022  URL: probml.ai
  * [70] Lawrence Hubert and Phipps Arabie  “Comparing partitions”  In _Journal of Classification_ 2.1, 1985, pp. 193–218  DOI: 10.1007/BF01908075
  * [71] Matthew Middlehurst, Patrick Schäfer and Anthony Bagnall  “Bake off redux: a review and experimental evaluation of recent time series classification algorithms”  In _Data Mining and Knowledge Discovery_ , 2024  DOI: 10.1007/s10618-024-01022-1
  * [72] Jesin Zakaria, Abdullah Mueen and Eamonn Keogh  “Clustering time series using unsupervised shapelets”  In _IEEE 12th International Conference on Data Mining_ , 2012, pp. 785–794  IEEE
  * [73] Nguyen Xuan Vinh, Julien Epps and James Bailey  “Information Theoretic Measures for Clusterings Comparison: Is a Correction for Chance Necessary?”  In _Proceedings of the 26th Annual International Conference on Machine Learning_ , ICML ’09  Montreal, Quebec, Canada: Association for Computing Machinery, 2009, pp. 1073–1080  DOI: 10.1145/1553374.1553511
  * [74] Tom Fawcett  “An introduction to ROC analysis” ROC Analysis in Pattern Recognition  In _Pattern Recognition Letters_ 27.8, 2006, pp. 861–874  DOI: https://doi.org/10.1016/j.patrec.2005.10.010
  * [75] T. Cover and P. Hart  “Nearest neighbor pattern classification”  In _IEEE Transactions on Information Theory_ 13.1, 1967, pp. 21–27  DOI: 10.1109/TIT.1967.1053964
  * [76] Angus Dempster, Daniel F. Schmidt and Geoffrey I. Webb  “Hydra: competing convolutional kernels for fast and accurate time series classification”  In _Data Mining and Knowledge Discovery_ 37.5, 2023, pp. 1779–1805  DOI: 10.1007/s10618-023-00939-3
  * [77] Jean-Yves Franceschi, Aymeric Dieuleveut and Martin Jaggi  “Unsupervised Scalable Representation Learning for Multivariate Time Series”  In _Advances in Neural Information Processing Systems_ 32 Curran Associates, Inc., 2019  URL: https://proceedings.neurips.cc/paper_files/paper/2019/file/53c6de78244e9f528eb3e1cda69699bb-Paper.pdf
  * [78] Jiaping Zhao and Laurent Itti  “shapeDTW: Shape Dynamic Time Warping”  In _Pattern Recognition_ 74, 2018, pp. 171–184  DOI: https://doi.org/10.1016/j.patcog.2017.09.020
  * [79] Thomas M. Cover and Joy A. Thomas  “Data Compression”  In _Elements of Information Theory_ John Wiley & Sons, Ltd, 2005, pp. 103–158  DOI: https://doi.org/10.1002/047174882X.ch5
  * [80] Leland McInnes, John Healy, Nathaniel Saul and Lukas Grossberger  “UMAP: Uniform Manifold Approximation and Projection”  In _Journal of Open Source Software_ 3.29 The Open Journal, 2018, pp. 861  DOI: 10.21105/joss.00861
  * [81] Pablo Millan Arias, Fatemeh Alipour, Kathleen A. Hill and Lila Kari  “DeLUCS: Deep learning for unsupervised clustering of DNA sequences”  In _PLOS ONE_ 17.1 Public Library of Science, 2022, pp. 1–25  DOI: 10.1371/journal.pone.0261531
  * [82] Yiming Bao et al.  “The influenza virus resource at the National Center for Biotechnology Information”  In _J. Virol._ 82.2, 2008, pp. 596–601
  * [83] Eneida L Hatcher et al.  “Virus Variation Resource - improved response to emergent viral outbreaks”  In _Nucleic Acids Res._ 45.D1, 2017, pp. D482–D490
  * [84] Lloyd N. Trefethen and David Bau  “Numerical Linear Algebra, Twenty-fifth Anniversary Edition”  Philadelphia, PA: Society for IndustrialApplied Mathematics, 2022  DOI: 10.1137/1.9781611977165
  * [85] Gene H. Golub and Charles F. Van Loan  “Matrix Computations - 4th Edition”  Philadelphia, PA: Johns Hopkins University Press, 2013  DOI: 10.1137/1.9781421407944
  * [86] Chang Wei Tan, Angus Dempster, Christoph Bergmeir and Geoffrey I. Webb  “MultiRocket: multiple pooling operators and transformations for fast and effective time series classification”  In _Data Mining and Knowledge Discovery_ 36.5, 2022, pp. 1623–1646  DOI: 10.1007/s10618-022-00844-1
  * [87] “Nomenclature Committee of the International Union of Biochemistry (NC-IUB). Nomenclature for incompletely specified bases in nucleic acid sequences. Recommendations 1984”  In _Biochemical Journal_ 229.2, 1985, pp. 281–286  DOI: 10.1042/bj2290281

##  Appendix A Implementation Details

###  A.1 First order methods for monotone VI with nuclear ball setup

We present the concrete accelerated mirror-prox method with backtracking for
nuclear norm constrained VI based on [61, 53], in Algorithm 1.

Algorithm 1 Accelerated Mirror Prox Method with Backtracking for Nuclear Norm
constrained VI

1:procedure AcceleratedBacktrackingMirrorProx($\Psi$, $N$, $\lambda$,
$\kappa_{0}$) $\triangleright$ $\Psi$ monotone VI, $N>0$ number of steps,
$\lambda>0$ radius of nuclear ball, $\kappa_{0}>0$ initial guess of Lipschitz
constant

2:
$\mathbf{R}_{1}:=\mathbf{0},\mathbf{B}_{1}:=\mathbf{R}_{1},\tilde{\mathbf{B}}_{1}:=\mathbf{R}_{1}$

3: for $t:=1,N$ do

4: $\hat{\kappa}_{t}:=\kappa_{t-1}$

5: repeat

6: $\alpha_{t}:=\frac{2}{t+1},\hat{\gamma}_{t}:=\frac{t}{2\hat{\kappa}_{t}t}$

7:
$\mathbf{B}_{t+1}:=\text{\sc{ProxNuc}}(\mathbf{R}_{t},\hat{\gamma}_{t}\Psi(\mathbf{R}_{t}))$

8: if
$\|\Psi(\mathbf{B}_{t+1})-\Psi(\mathbf{R}_{t})\|_{F}>\hat{\kappa}_{t}\|\mathbf{B}_{t+1}-\mathbf{R}_{t}\|_{F}$
then

9: $\hat{\kappa}_{t}:=2\hat{\kappa}_{t}$

10: end if

11: until
$\|\Psi(\mathbf{B}_{t+1})-\Psi(\mathbf{R}_{t})\|_{F}\leq\hat{\kappa}_{t}\|\mathbf{B}_{t+1}-\mathbf{R}_{t}\|_{F}$

12: $\kappa_{t}:=\hat{\kappa}_{t}$

13: $\gamma_{t}=\hat{\gamma}_{t}$

14:
$\mathbf{R}_{t+1}:=\text{\sc{ProxNuc}}(\mathbf{R}_{t},{\gamma}_{t}\Psi(\mathbf{B}_{t+1}))$

15:
$\tilde{\mathbf{B}}_{t+1}:=(1-\alpha_{t})\tilde{\mathbf{B}}_{t}+\alpha_{t}\mathbf{B}_{t+1}$

16: return $\tilde{\mathbf{B}}_{N+1}$

17: end for

18:end procedure

19:procedure ProxNuc($\mathbf{Z}$,$\mathbf{X},\lambda$)$\triangleright$
Computes
$\operatorname{Prox}_{\mathbf{Z},\|\cdot\|_{*}\leq\lambda}(\mathbf{X})$,
$\mathbf{Z}\in\mathbb{R}^{m\times n}$

20: $r:=\min(m,n)$

21: $q:=(2\log 2r)^{-1}$

22: $\mathbf{Y}:=x-\partial\omega(\mathbf{Z})$

23: $\mathbf{U},\bm{\delta},\mathbf{V}^{*}:=\text{\sc{svd}}(\mathbf{Y})$

24:
$\mathbf{t}:=\min_{\mathbf{t}\in\mathbb{R}^{r}}\\{\frac{1}{2}\sum_{j=1}^{m}t_{j}^{2}-\delta_{j}t_{j}\mid
t\geq\mathbf{0},\sum_{j=1}^{m}t_{j}\leq\lambda\\}$

25: return $\mathbf{U}\operatorname{diag}(-\mathbf{t})\mathbf{V}^{*}$

26:end procedure

27:procedure $\partial\omega$($\mathbf{Z}$) $\triangleright$ Compute
subgradient of DGF $\omega$ from (9), $\mathbf{Z}\in\mathbb{R}^{m\times n}$

28: $r:=\min(m,n)$

29: $q:=(2\log 2r)^{-1}$

30: $c:=\frac{4\sqrt{e}\log(2r)}{2^{q}(1+q)}$

31: $\mathbf{U},\bm{\sigma},\mathbf{V}^{*}:=\text{\sc{svd}}(\mathbf{Y})$

32: return $c\sum_{i}^{r}(1+q)(\sigma_{i})^{q}\mathbf{u}_{i}\mathbf{v}^{*}$

33:end procedure

##  Appendix B Detailed Experimental Setup and Results

###  B.1 Synthetic sequences

####  B.1.1 Data Generation

For the synthetic sequence recovery experiment, we adopt the following data-
generating procedure: We take the order of the sequences to be $d=15$, and we
generate data according to the following procedure within each of the five
generated classes of observations

  1. 1.

Pick a baseline set of coefficients according to a given random distribution

  2. 2.

For each of the $N=300$ sequences to generate, perturb the coefficients
according to the pre-specified rule

  3. 3.

Generate the data matrix of size consisting of $T=250$ of the $N=300$
sequences according to the perturbed coefficients such that the data obeys
(3). To do so, we seat the first 15 observations using random noise such that
$x_{i,t}\sim\mathcal{N}(\mu=0,\sigma^{2}=1),\forall t\in[1,d],i\in[N]$. Then
each successive entry $x_{i,t},t\in[d+1,T]$ is then given by taking
$x_{i,t}=\sum_{s=1}^{d}{b}_{i,s}x_{i,t-s}+\epsilon_{i,t}$,
$\epsilon_{i,t}\sim\mathcal{N}(\mu=0,\sigma^{2}=0.02)$.

We draw the ten generated classes of data from the following five generation
procedures given in Table 3. We use each procedure twice to generate the ten
classes of data. We denote the coefficients common to the sequences (for some
class) as $\mathbf{b}_{\text{common}}$, and the coefficients for the
$i^{\text{th}}$ sequence in said class as $\mathbf{b}_{i}$.

The baseline coefficients generation methods are given as:

Exponentially Time Decaying:

$b_{\mathrm{common},s}=Z\gamma^{s}/(\sum_{j=1}^{d}\gamma^{j})\qquad
Z\sim\mathrm{Uniform}([0,1]),\forall s\in[d]$

Uniform:

$b_{\mathrm{common},s}=Z\qquad Z\sim\mathrm{Uniform}([0,1/2d]),\forall
s\in[d]$

and the perturbation methods are given as:

Gaussian:

$\mathbf{b}_{i}=\mathbf{b}_{\mathrm{common}}+\mathbf{Z}\qquad
Z_{j}\sim\mathcal{N}(\mu=0,\sigma^{2}=0.02)$

$d/3$ Most Recent:

$\mathbf{b}_{i}=\mathbf{b}_{\mathrm{common}}+\mathbf{Z}\qquad
Z_{j}\sim\begin{cases}\mathcal{N}(\mu=0,\sigma^{2}=0.02)&j<\lceil d/3\rceil\\\
0\end{cases},\forall i\in[N]$

Uniform $\times$ Fixed Vector:

$\mathbf{b}_{i}=\mathbf{b}_{\mathrm{common}}+\theta\mathbf{v}\qquad\theta\sim\mathrm{Uniform}([-1,1]),\|\mathbf{v}\|=1,\forall
i\in[N]$
($\mathbf{v}$ chosen uniformly on a unit hypersphere, and is the same for all
sequences generated in the class)

Table 3: Five classes of sequence generating procedure Baseline Coefficients | Perturbation Pattern
---|---
Exponentially Time Decaying | Gaussian
Exponentially Time Decaying |  $d/3$ Most Recent
Exponentially Time Decaying | Uniform $\times$ Fixed Vector
Uniform | Gaussian
Uniform | Uniform $\times$ Fixed Vector

####  B.1.2 Parameter Recovery

For the parameter recovery experiment, we take all $\binom{10}{3}=120$
combinations of $k=3$ sequences from the 10 classes and concatenate the
generated sequences to form a matrix of $900$ observations. We then recover
the baseline coefficient matrix $\mathbf{B}\in\mathbb{R}^{15\times 900}$. To
recover the parameters for each sequence, we solve the program given in (7) to
optimality for differing levels of $\lambda$ (using a standard convex solver
by contrast to Algorithm 1). To find which levels of $\lambda$ to solve for,
we first solve the unconstrained version of the problem. We then compute the
nuclear norm of recovered $\|\widehat{\mathbf{B}}\|_{*}$. We then successively
search for the optimal $\lambda^{*}$ on the interval
$[0,\|\widehat{\mathbf{B}}\|_{*}]$ using the relative reconstruction error
$\|\mathbf{B}-\widehat{\mathbf{B}}\|_{F}/\|\mathbf{B}\|_{F}$ as the objective
until we achieve an absolute tolerance of $10^{-3}$. We report the results in
Table 1 using the matrices $\widehat{\mathbf{B}}$ and
$\widehat{\mathbf{B}}_{\lambda^{*}}$ across the 120 runs. For Figure 1, we use
data drawn from the classes uniform baseline coefficients with Gaussian
perturbation, exponentially time-decaying baseline coefficients with Gaussian
perturbation; and exponentially time-decaying baseline coefficients with
uniform*fixed vector perturbation. We find the $\lambda^{*}$ via Brent search,
and we sweep 40 values of $\lambda\in[16.3,25.2]$ (the right bound
corresponding to the value of $\|\widehat{\mathbf{B}}\|_{*}$) to produce
Figure 1. When reporting the spectra of singular values, we pick solutions
that correspond to the largest value of $\lambda$ we have for some fixed
number of large singular values. Finally, we depict the principal components
according to the recovered matrix $\widehat{\mathbf{B}}_{\lambda^{*}}$, which
is optimal in the sense of the reconstruction error (found via Brent search).

###  B.2 UCR time series

####  B.2.1 Data Overview

We compare our method with the following representative time series clustering
methods. In line with [40] we selected 36 of the univariate UCR time series
classification datasets
111https://www.cs.ucr.edu/%7Eeamonn/time_series_data_2018/. We report basic
statistics (training samples, testing samples, length, and number of classes)
of the datasets in Table 4.

Table 4: Basic statistics for the UCR time-series classification benchmark data | Training Samples | Testing Samples | Length | Classes
---|---|---|---|---
ArrowHead | 36 | 175 | 251 | 3
Beef | 30 | 30 | 470 | 5
BeetleFly | 20 | 20 | 512 | 2
BirdChicken | 20 | 20 | 512 | 2
Car | 60 | 60 | 577 | 4
ChlorineConc. | 467 | 3840 | 166 | 3
Coffee | 28 | 28 | 286 | 2
DiatomsizeReduction | 16 | 306 | 345 | 4
Dist.Pha.Outln.AgeGrp. | 400 | 139 | 80 | 3
Dist.Pha.Outln.Correct | 600 | 276 | 80 | 2
ECG200 | 100 | 100 | 96 | 2
ECGFiveDays | 23 | 861 | 136 | 2
GunPoint | 50 | 150 | 150 | 2
Ham | 109 | 105 | 431 | 2
Herring | 64 | 64 | 512 | 2
Lightning2 | 60 | 61 | 637 | 2
Meat | 60 | 60 | 448 | 3
Mid.Pha.Outln.AgeGrp. | 400 | 154 | 80 | 3
Mid.Pha.Outln.Correct | 600 | 291 | 80 | 2
Mid.PhalanxTW | 399 | 154 | 80 | 6
MoteStrain | 20 | 1252 | 84 | 2
OSULeaf | 200 | 242 | 427 | 6
Plane | 105 | 105 | 144 | 7
Prox.Pha.Outln.AgeGrp. | 400 | 205 | 80 | 3
Prox.PhalanxTW | 400 | 205 | 80 | 6
SonyAIBORobotSurf.1 | 20 | 601 | 70 | 2
SonyAIBORobotSurf.2 | 27 | 953 | 65 | 2
SwedishLeaf | 500 | 625 | 128 | 15
Symbols | 25 | 995 | 398 | 6
ToeSegmentation1 | 40 | 228 | 277 | 2
ToeSegmentation2 | 36 | 130 | 343 | 2
TwoPatterns | 1000 | 4000 | 128 | 4
TwoLeadECG | 23 | 1139 | 82 | 2
Wafer | 1000 | 6164 | 152 | 2
Wine | 57 | 54 | 234 | 2
WordSynonyms | 267 | 638 | 270 | 25

####  B.2.2 Overview of Evaluation Methods

We evaluate the classification performance using the following methods:

ARI:

Similarity of learned and ground truth assignments [73]. For matched
clustering partitions

$\mathrm{RI}=(a+b)/\binom{N}{2},\qquad\mathrm{ARI}=\frac{\mathrm{RI}-\mathbb{E}[\mathrm{RI}]}{\max(\text{RI})-\mathbb{E}[RI]]}$

NMI:

The mutual information between the true class labels and the cluster
assignments, normalized by the entropy of the true labels and the cluster
assignments [73].

$\text{NMI}=\frac{2I(X;Y)}{H(X)+H(Y)}$

where $X,Y$ are the true and assigned labels, $H(X)$ is the entropy of $X$,
and $I(X;Y)$ is the mutual information between $X$ and $Y$.

Accuracy:

Proportion of correct predictions to a total number of predictions.

F1:

Harmonic mean of the precision and recall. In the multiclass case, we take the
macro average by calculating the metric for each label and computing their
unweighted mean.

$F_{1}=\frac{2\times\mathrm{TP}}{2\times\mathrm{TP}+\mathrm{FP}+\mathrm{FN}}$

Runtime:

Runtime of the algorithm in terms the user CPU time in the computational
setting described in 4. If the method is GPU accelerated we report the user
CPU/GPU time spent in the routine.

####  B.2.3 Overview of Methods

In addition to our method, we evaluate the performance of the following
baseline methods

$\ell_{2}$+KNN

K-Nearest Neighbors Classification with the distance metric as the Euclidean
distance between two time-series treating the entire observation sequence as a
high dimensional vector [75].

DTW+KNN

K-Nearest Neighbors Classification with the distance metric calculated
according to DTW [31], which aims to align the two given sequences by solving
the following program

$\mathrm{DTW}_{q}(\mathbf{x},\mathbf{x}^{\prime})=\min_{\pi\in\mathcal{A}(\mathbf{x},\mathbf{x}^{\prime})}\langle
A_{\pi},D_{q}(\mathbf{x},\mathbf{x}^{\prime})\rangle^{1/q}.$

The set $\mathcal{A}(\mathbf{x},\mathbf{x}^{\prime})$ is the set of all
admissible paths as represented by boolean matrices. Non-zero entries
correspond to matching time series elements in the path. A path is admissible
if the beginning and end of the time series are matched together, the sequence
is monotone in both $i$ and $j$, and all entries appear at least once. We take
$q=2$ as the Euclidean metric.

shapeDTW

Extension to DTW scheme by incorporating point-wise local structures into the
matching procedure [78]. Examples of such shape descriptors include data
itself, a rolling average of, a discrete wavelet transform, and a finite
difference/derivative. Finally, the encoded sequences are then aligned by DTW
and used for nearest neighbor classification.

MR-Hydra

Combination of dictionary-based Multirocket and Hydra algorithms for time
series classification, extracts and counts symbolic patterns using competing
convolutional kernels [76, 86].

TS2Vec

Construct an encoder network for time series embedding based on hierarchical
contrastive learning [36]. The discrimination is done both between sequences
and within the sequences themselves. The encoder network consists of an input
projection layer, a timestamp masking module, and a dilated convolutional
module, and is optimized jointly with temporal and cross-sequence contrastive
loss.

####  B.2.4 Detailed Experimental Procedure

We split our data into testing and training splits according to those given by
the UCR repository. For the methods that directly perform classification (KNN,
shapeDTW, Inception Time), we train on the test set and then report the the
performance on the training set. In line with [36, 77], to evaluate the
classification performance on test set for methods which produce embeddings
(TS2Vec and our method), we perform cross-validated grid search (based on
$k=5$ folds, ) across KNNs with $k=\\{2^{i}\mid i\in[0,4]\\}$ neighbors or
SVMs with RBF kernels with penalty values $c\in\\{2^{i}\mid
i\in[-10,15]\\}\cup{\infty}$. For the KNN-based methods, we do the same grid
search as outlined above. For our own method, we also grid search across
parameters of $\lambda$ and report the performance for the best choice under
rank constraint. To find the embedding, we run Algorithm 1 for 256 iterations.

###  B.3 Detailed results

In Tables 5 and 6, we present the classification performance for the discussed
metrics for the evaluated methods for each of the tested UCR datasets.

Table 5: Detailed results per UCR dataset (Part I) Dataset | Method | ARI | NMI | Acc. | F1 | RT | Method | ARI | NMI | Acc. | F1 | RT
---|---|---|---|---|---|---|---|---|---|---|---|---
ArrowHead |  $\ell_{2}$+ | 0.482 | 0.453 | 0.800 | 0.800 | 0.007 | MR | 0.629 | 0.585 | 0.863 | 0.863 | 3.183
Beef | KNN | 0.322 | 0.518 | 0.667 | 0.672 | 0.001 | Hydra | 0.454 | 0.627 | 0.767 | 0.768 | 1.967
BeetleFly |  | 0.219 | 0.344 | 0.750 | 0.733 | 0.001 |  | 0.621 | 0.619 | 0.900 | 0.899 | 1.532
BirdChicken |  | -0.044 | 0.007 | 0.550 | 0.549 | 0.001 |  | 0.621 | 0.619 | 0.900 | 0.899 | 1.536
Car |  | 0.403 | 0.477 | 0.733 | 0.737 | 0.003 |  | 0.830 | 0.860 | 0.933 | 0.933 | 4.023
ChlorineConc. |  | 0.231 | 0.157 | 0.650 | 0.610 | 0.633 |  | 0.472 | 0.373 | 0.789 | 0.753 | 57.36
Coffee |  | 1.000 | 1.000 | 1.000 | 1.000 | 0.001 |  | 1.000 | 1.000 | 1.000 | 1.000 | 1.411
DiatomsizeReduction |  | 0.872 | 0.830 | 0.935 | 0.883 | 0.004 |  | 0.921 | 0.896 | 0.964 | 0.947 | 6.537
Dist.Pha.Outln.AgeGrp. |  | 0.190 | 0.224 | 0.626 | 0.613 | 0.019 |  | 0.383 | 0.404 | 0.770 | 0.775 | 4.794
Dist.Pha.Outln.Correct |  | 0.181 | 0.137 | 0.717 | 0.684 | 0.054 |  | 0.366 | 0.286 | 0.804 | 0.790 | 8.056
ECG200 |  | 0.571 | 0.445 | 0.880 | 0.868 | 0.004 |  | 0.667 | 0.542 | 0.910 | 0.902 | 1.765
ECGFiveDays |  | 0.352 | 0.304 | 0.797 | 0.794 | 0.011 |  | 1.000 | 1.000 | 1.000 | 1.000 | 7.914
GunPoint |  | 0.681 | 0.578 | 0.913 | 0.913 | 0.004 |  | 1.000 | 1.000 | 1.000 | 1.000 | 2.064
Ham |  | 0.031 | 0.029 | 0.600 | 0.600 | 0.006 |  | 0.229 | 0.177 | 0.743 | 0.742 | 5.055
Herring |  | -0.015 | 0.003 | 0.516 | 0.516 | 0.003 |  | 0.207 | 0.155 | 0.734 | 0.726 | 3.825
Lightning2 |  | 0.246 | 0.193 | 0.754 | 0.750 | 0.003 |  | 0.104 | 0.084 | 0.672 | 0.665 | 5.208
Meat |  | 0.799 | 0.797 | 0.933 | 0.935 | 0.003 |  | 0.810 | 0.808 | 0.933 | 0.933 | 3.096
Mid.Pha.Outln.AgeGrp. |  | 0.055 | 0.026 | 0.519 | 0.443 | 0.021 |  | 0.092 | 0.071 | 0.591 | 0.491 | 4.599
Mid.Pha.Outln.Correct |  | 0.280 | 0.208 | 0.766 | 0.756 | 0.057 |  | 0.475 | 0.372 | 0.845 | 0.842 | 8.154
Mid.PhalanxTW |  | 0.379 | 0.367 | 0.513 | 0.382 | 0.020 |  | 0.383 | 0.433 | 0.513 | 0.339 | 5.382
MoteStrain |  | 0.573 | 0.467 | 0.879 | 0.877 | 0.015 |  | 0.794 | 0.699 | 0.946 | 0.945 | 7.062
OSULeaf |  | 0.298 | 0.383 | 0.521 | 0.525 | 0.023 |  | 0.921 | 0.919 | 0.963 | 0.956 | 10.34
Plane |  | 0.919 | 0.943 | 0.962 | 0.963 | 0.005 |  | 1.000 | 1.000 | 1.000 | 1.000 | 2.505
Prox.Pha.Outln.AgeGrp. |  | 0.492 | 0.422 | 0.785 | 0.693 | 0.027 |  | 0.662 | 0.564 | 0.868 | 0.797 | 5.489
Prox.PhalanxTW |  | 0.584 | 0.566 | 0.707 | 0.444 | 0.027 |  | 0.718 | 0.671 | 0.805 | 0.490 | 4.922
SonyAIBORobotSurf.1 |  | 0.148 | 0.280 | 0.696 | 0.688 | 0.007 |  | 0.598 | 0.570 | 0.887 | 0.887 | 3.396
SonyAIBORobotSurf.2 |  | 0.514 | 0.395 | 0.859 | 0.849 | 0.013 |  | 0.782 | 0.682 | 0.942 | 0.940 | 4.828
SwedishLeaf |  | 0.629 | 0.761 | 0.789 | 0.782 | 0.109 |  | 0.950 | 0.965 | 0.976 | 0.977 | 10.78
Symbols |  | 0.791 | 0.843 | 0.899 | 0.898 | 0.017 |  | 0.955 | 0.954 | 0.981 | 0.981 | 21.89
ToeSegmentation1 |  | 0.126 | 0.095 | 0.680 | 0.675 | 0.006 |  | 0.832 | 0.782 | 0.956 | 0.956 | 4.760
ToeSegmentation2 |  | 0.340 | 0.244 | 0.808 | 0.744 | 0.003 |  | 0.640 | 0.464 | 0.915 | 0.866 | 3.608
TwoPatterns |  | 0.770 | 0.726 | 0.907 | 0.906 | 1.328 |  | 1.000 | 1.000 | 1.000 | 1.000 | 50.14
TwoLeadECG |  | 0.244 | 0.217 | 0.747 | 0.741 | 0.015 |  | 0.993 | 0.983 | 0.998 | 0.998 | 6.182
Wafer |  | 0.971 | 0.923 | 0.995 | 0.988 | 2.088 |  | 0.998 | 0.993 | 1.000 | 0.999 | 76.27
Wine |  | 0.031 | 0.036 | 0.611 | 0.611 | 0.002 |  | 0.720 | 0.687 | 0.926 | 0.926 | 1.718
WordSynonyms |  | 0.537 | 0.571 | 0.618 | 0.465 | 0.069 |  | 0.725 | 0.753 | 0.777 | 0.658 | 15.76
ArrowHead | DTW+ | 0.312 | 0.282 | 0.703 | 0.700 | 2.139 | TS2 | 0.480 | 0.462 | 0.794 | 0.794 | 81.08
Beef | KNN | 0.276 | 0.490 | 0.633 | 0.629 | 1.158 | Vec | 0.284 | 0.494 | 0.667 | 0.670 | 109.0
BeetleFly |  | 0.131 | 0.275 | 0.700 | 0.670 | 0.618 |  | 0.800 | 0.761 | 0.950 | 0.950 | 79.59
BirdChicken |  | 0.212 | 0.221 | 0.750 | 0.744 | 0.606 |  | 0.621 | 0.619 | 0.900 | 0.899 | 80.60
Car |  | 0.446 | 0.501 | 0.733 | 0.728 | 7.565 |  | 0.709 | 0.787 | 0.867 | 0.867 | 298.9
ChlorineConc. |  | 0.231 | 0.154 | 0.648 | 0.607 | 247.313 |  | 0.432 | 0.333 | 0.764 | 0.730 | 2439
Coffee |  | 1.000 | 1.000 | 1.000 | 1.000 | 0.336 |  | 0.857 | 0.811 | 0.964 | 0.964 | 129.1
DiatomsizeReduction |  | 0.938 | 0.921 | 0.967 | 0.942 | 3.272 |  | 0.968 | 0.952 | 0.984 | 0.973 | 87.44
Dist.Pha.Outln.AgeGrp. |  | 0.389 | 0.368 | 0.770 | 0.763 | 1.804 |  | 0.272 | 0.277 | 0.705 | 0.699 | 2046
Dist.Pha.Outln.Correct |  | 0.183 | 0.132 | 0.717 | 0.690 | 5.382 |  | 0.246 | 0.176 | 0.750 | 0.737 | 2882
ECG200 |  | 0.280 | 0.192 | 0.770 | 0.749 | 0.468 |  | 0.540 | 0.417 | 0.870 | 0.858 | 463.5
ECGFiveDays |  | 0.286 | 0.252 | 0.768 | 0.763 | 1.848 |  | 0.991 | 0.979 | 0.998 | 0.998 | 80.84
GunPoint |  | 0.659 | 0.557 | 0.907 | 0.907 | 0.847 |  | 0.973 | 0.949 | 0.993 | 0.993 | 234.7
Ham |  | -0.005 | 0.003 | 0.467 | 0.467 | 12.116 |  | 0.210 | 0.168 | 0.733 | 0.733 | 526.1
Herring |  | -0.012 | 0.001 | 0.531 | 0.520 | 6.536 |  | 0.064 | 0.047 | 0.641 | 0.625 | 326.8
Lightning2 |  | 0.537 | 0.480 | 0.869 | 0.864 | 8.898 |  | 0.318 | 0.252 | 0.787 | 0.783 | 301.2
Meat |  | 0.799 | 0.797 | 0.933 | 0.935 | 4.497 |  | 0.687 | 0.714 | 0.883 | 0.883 | 300.1
Mid.Pha.Outln.AgeGrp. |  | 0.024 | 0.022 | 0.500 | 0.411 | 2.066 |  | 0.038 | 0.029 | 0.519 | 0.426 | 2042
Mid.Pha.Outln.Correct |  | 0.153 | 0.109 | 0.698 | 0.691 | 5.663 |  | 0.385 | 0.297 | 0.811 | 0.808 | 3055
Mid.PhalanxTW |  | 0.380 | 0.368 | 0.506 | 0.374 | 1.996 |  | 0.420 | 0.404 | 0.545 | 0.396 | 1997
MoteStrain |  | 0.448 | 0.351 | 0.835 | 0.834 | 0.900 |  | 0.528 | 0.424 | 0.863 | 0.862 | 86.43
OSULeaf |  | 0.309 | 0.392 | 0.591 | 0.588 | 50.818 |  | 0.644 | 0.671 | 0.822 | 0.799 | 1071
Plane |  | 1.000 | 1.000 | 1.000 | 1.000 | 1.172 |  | 1.000 | 1.000 | 1.000 | 1.000 | 541.7
Prox.Pha.Outln.AgeGrp. |  | 0.504 | 0.430 | 0.805 | 0.716 | 2.696 |  | 0.506 | 0.437 | 0.780 | 0.689 | 2051
Prox.PhalanxTW |  | 0.644 | 0.587 | 0.756 | 0.511 | 2.673 |  | 0.674 | 0.628 | 0.771 | 0.562 | 2052
SonyAIBORobotSurf.1 |  | 0.200 | 0.316 | 0.725 | 0.721 | 0.309 |  | 0.588 | 0.571 | 0.884 | 0.883 | 84.73
SonyAIBORobotSurf.2 |  | 0.435 | 0.324 | 0.831 | 0.817 | 0.571 |  | 0.671 | 0.584 | 0.910 | 0.907 | 125.5
SwedishLeaf |  | 0.639 | 0.770 | 0.792 | 0.787 | 25.957 |  | 0.875 | 0.916 | 0.936 | 0.937 | 2573
Symbols |  | 0.891 | 0.913 | 0.950 | 0.949 | 21.459 |  | 0.928 | 0.936 | 0.969 | 0.969 | 132.7
ToeSegmentation1 |  | 0.293 | 0.260 | 0.772 | 0.762 | 3.704 |  | 0.816 | 0.723 | 0.952 | 0.952 | 213.5
ToeSegmentation2 |  | 0.398 | 0.249 | 0.838 | 0.764 | 2.904 |  | 0.714 | 0.631 | 0.931 | 0.899 | 170.8
TwoPatterns |  | 1.000 | 1.000 | 1.000 | 1.000 | 330.318 |  | 0.973 | 0.964 | 0.990 | 0.990 | 5216
TwoLeadECG |  | 0.654 | 0.564 | 0.904 | 0.904 | 0.901 |  | 0.958 | 0.918 | 0.989 | 0.989 | 87.20
Wafer |  | 0.867 | 0.748 | 0.980 | 0.944 | 720.459 |  | 0.942 | 0.868 | 0.991 | 0.976 | 5386
Wine |  | 0.003 | 0.016 | 0.574 | 0.574 | 0.860 |  | 0.339 | 0.271 | 0.796 | 0.796 | 303.5
WordSynonyms |  | 0.575 | 0.600 | 0.649 | 0.533 | 66.754 |  | 0.349 | 0.421 | 0.522 | 0.309 | 1407
Table 6: Detailed results per UCR dataset (Part II) Dataset | Method | ARI | NMI | Acc. | F1 | RT | Method | ARI | NMI | Acc. | F1 | RT
---|---|---|---|---|---|---|---|---|---|---|---|---
ArrowHead | shape | 0.521 | 0.492 | 0.817 | 0.818 | 0.672 | Ours | 0.336 | 0.306 | 0.720 | 0.720 | 136.7
Beef | DTW | 0.322 | 0.518 | 0.667 | 0.672 | 0.214 |  | 0.453 | 0.654 | 0.733 | 0.736 | 70.59
BeetleFly |  | 0.219 | 0.344 | 0.750 | 0.733 | 0.107 |  | 1.000 | 1.000 | 1.000 | 1.000 | 51.50
BirdChicken |  | -0.044 | 0.007 | 0.550 | 0.549 | 0.106 |  | 1.000 | 1.000 | 1.000 | 1.000 | 51.92
Car |  | 0.560 | 0.585 | 0.817 | 0.815 | 1.077 |  | 0.372 | 0.457 | 0.667 | 0.648 | 173.1
ChlorineConc. |  | 0.199 | 0.133 | 0.628 | 0.587 | 120.800 |  | 0.537 | 0.414 | 0.811 | 0.782 | 2252
Coffee |  | 1.000 | 1.000 | 1.000 | 1.000 | 0.105 |  | 1.000 | 1.000 | 1.000 | 1.000 | 39.92
DiatomsizeReduction |  | 0.921 | 0.890 | 0.958 | 0.921 | 0.849 |  | 0.818 | 0.865 | 0.882 | 0.704 | 331.7
Dist.Pha.Outln.AgeGrp. |  | 0.209 | 0.251 | 0.633 | 0.615 | 1.584 |  | 0.323 | 0.402 | 0.741 | 0.748 | 103.9
Dist.Pha.Outln.Correct |  | 0.188 | 0.140 | 0.721 | 0.690 | 4.733 |  | 0.316 | 0.233 | 0.783 | 0.772 | 158.0
ECG200 |  | 0.541 | 0.420 | 0.870 | 0.860 | 0.397 |  | 0.668 | 0.554 | 0.910 | 0.904 | 43.80
ECGFiveDays |  | 0.705 | 0.605 | 0.920 | 0.920 | 1.279 |  | 0.765 | 0.666 | 0.937 | 0.937 | 306.1
GunPoint |  | 0.845 | 0.761 | 0.960 | 0.960 | 0.517 |  | 0.845 | 0.761 | 0.960 | 0.960 | 70.03
Ham |  | 0.031 | 0.029 | 0.600 | 0.600 | 2.714 |  | 0.160 | 0.124 | 0.705 | 0.704 | 265.6
Herring |  | -0.012 | 0.006 | 0.531 | 0.531 | 1.129 |  | 0.176 | 0.128 | 0.719 | 0.688 | 162.0
Lightning2 |  | 0.358 | 0.299 | 0.803 | 0.797 | 1.345 |  | 0.399 | 0.320 | 0.820 | 0.817 | 230.8
Meat |  | 0.799 | 0.797 | 0.933 | 0.935 | 0.860 |  | 0.856 | 0.841 | 0.950 | 0.950 | 133.3
Mid.Pha.Outln.AgeGrp. |  | 0.053 | 0.022 | 0.513 | 0.432 | 1.722 |  | 0.184 | 0.143 | 0.649 | 0.523 | 106.6
Mid.Pha.Outln.Correct |  | 0.281 | 0.207 | 0.766 | 0.759 | 5.110 |  | 0.410 | 0.324 | 0.821 | 0.813 | 177.9
Mid.PhalanxTW |  | 0.357 | 0.360 | 0.487 | 0.361 | 1.740 |  | 0.362 | 0.436 | 0.591 | 0.334 | 107.5
MoteStrain |  | 0.573 | 0.467 | 0.879 | 0.877 | 0.804 |  | 0.212 | 0.175 | 0.731 | 0.731 | 243.3
OSULeaf |  | 0.316 | 0.411 | 0.566 | 0.567 | 10.754 |  | 0.600 | 0.625 | 0.810 | 0.798 | 584.7
Plane |  | 0.937 | 0.961 | 0.971 | 0.972 | 0.644 |  | 1.000 | 1.000 | 1.000 | 1.000 | 74.50
Prox.Pha.Outln.AgeGrp. |  | 0.482 | 0.399 | 0.780 | 0.688 | 2.244 |  | 0.681 | 0.584 | 0.883 | 0.808 | 119.3
Prox.PhalanxTW |  | 0.585 | 0.565 | 0.702 | 0.426 | 2.240 |  | 0.752 | 0.728 | 0.834 | 0.436 | 118.7
SonyAIBORobotSurf.1 |  | 0.206 | 0.333 | 0.729 | 0.724 | 0.282 |  | 0.619 | 0.572 | 0.894 | 0.893 | 102.1
SonyAIBORobotSurf.2 |  | 0.589 | 0.468 | 0.885 | 0.876 | 0.560 |  | 0.767 | 0.655 | 0.938 | 0.934 | 139.3
SwedishLeaf |  | 0.697 | 0.806 | 0.830 | 0.827 | 16.027 |  | 0.805 | 0.864 | 0.899 | 0.897 | 441.6
Symbols |  | 0.823 | 0.864 | 0.918 | 0.917 | 4.945 |  | 0.903 | 0.917 | 0.956 | 0.956 | 1127
ToeSegmentation1 |  | 0.221 | 0.181 | 0.737 | 0.728 | 1.136 |  | 0.678 | 0.580 | 0.912 | 0.912 | 203.4
ToeSegmentation2 |  | 0.486 | 0.379 | 0.862 | 0.809 | 0.780 |  | 0.655 | 0.471 | 0.923 | 0.868 | 158.2
TwoPatterns |  | 0.908 | 0.870 | 0.965 | 0.964 | 204.000 |  | 0.679 | 0.145 | 0.143 | 0.514 | 1719
TwoLeadECG |  | 0.484 | 0.438 | 0.848 | 0.846 | 0.742 |  | 0.993 | 0.981 | 0.998 | 0.998 | 240.5
Wafer |  | 0.977 | 0.936 | 0.996 | 0.991 | 373.891 |  | 0.939 | 0.858 | 0.991 | 0.975 | 3329
Wine |  | 0.016 | 0.025 | 0.593 | 0.593 | 0.300 |  | 0.150 | 0.128 | 0.704 | 0.702 | 62.29
WordSynonyms |  | 0.578 | 0.600 | 0.639 | 0.487 | 20.963 |  | 0.246 | 0.325 | 0.395 | 0.220 | 764.7

###  B.4 Natural language embedding

To acquire the data, we retrieved the raw text of Alice’s Adventures in
Wonderland and Through the Looking Glass from Project Gutenberg
222https://www.gutenberg.org. For the paper abstracts, we used the training
portion of the ML-ArXiv-Papers dataset
333https://www.kaggle.com/datasets/Cornell-University/arxiv. For each dataset,
we stripped all non ASCII characters and uncommon punctuation (<,>, ‘, =, |,
?, &, [, ], *,  , !, #, @, and ").

After acquiring the data, we then encoded using a Huffman tree with $n=4$
symbols derived from the frequency of letters in our corpora. We treated each
abstract as a document and considered 500-character chunks of the two books.
We rejected abstracts containing less than 500 words. After encoding the
sequences using the Huffman code, we cut off each sequence at 1000 coded
symbols and rejected all sequences less than this length after coding. This
left us with $n=228$ samples from “Alice’s Adventures in Wonderland”, $n=316$
samples from “Through the Looking Glass”, and $n=600$ machine learning-related
ArXiv abstracts.

To learn the embedding, we use the method described in Appendix A.1 and grid
searched across values of $\lambda$ for $512$ steps using the softmax link
function described in Section 3.2.

###  B.5 Gene sequence embedding

##### Data acquisition and processing

In line with [81], we downloaded viral genome sequences for two different
kinds of human viruses: Influenza A virus and Dengue virus. We consider
different strains of each virus in addition to the species as a whole. We
provide a textual description below. In Table 7, we provide summary
statistics, including the number of sequences in the strain, the average and
standard deviation of the sequence lengths, and the length of the shortest and
longest sequences in the strains.

Influenza A

The Influenza A virus genome data ($n=949$) is acquired from the NCBI
Influenza Virus Resource [82]. We consider the genome of segment 6, which
encodes the neuraminidase protein, and include sequence samples belonging to
subtypes H1N1, H2N2, H5N1, H7N3, and H7N9.

Denuge

We consider $n=1562$ full Dengue virus genomes downloaded from the NCBI Virus
Variation Resource [83]. We consider all four subtypes of the virus DENV-1,
DENV-2, DENV-3, and DENV-4.

We encoded all the sequences as four-channel signals via one-hot encoding,
with each nucleotide (A,C,T,G) corresponding to one of the channels. In the
case we encounter incompletely specified bases in the nucleic acid sequences
[87], we give equiprobable weights to the possible corresponding nucleotides.

##### Learning procedure

Same as the natural language case, we represent the data as a four channel
signal and adopt the softmax activation scheme as described in Section 3.2.
Since the sequences are of considerable and variable length ($>1000$
nucleotides, see Table 7), we adopt a stochastic estimation to (13) by
randomly sampling length $G=800$ sub-windows from each of sequences. We take
the sample average for each of the sub-window observations similar to (14). We
run Algorithm 1 for $N=1024$ iterations, using the stochastic approximation
described above and grid searching across values of $\lambda$. To produce
Figures 2(b) and 2(c), we took the learned representations and projected them
into two dimensions via UMAP.

Table 7: Statistics for selected viral genomes. Virus | Strain | Count | Sequence Length
---|---|---|---
|  |  | Avg. | Std. | Min. | Max.
Influenza-A | H5N1 | 188 | 1368.521 | (21.682) | 1350 | 1457
| H1N1 | 191 | 1421.0 | (15.25) | 1350 | 1468
| H7N9 | 190 | 1403.521 | (12.048) | 1389 | 1444
| H2N2 | 187 | 1430.053 | (17.87) | 1376 | 1467
| H7N3 | 193 | 1423.15 | (21.537) | 1345 | 1468
| Total | 949 | 1409.326 | (28.512) | 1345 | 1468
Dengue | DENV-1 | 409 | 10577.812 | (194.4) | 10176 | 10821
| DENV-2 | 409 | 10592.504 | (196.308) | 10173 | 10991
| DENV-3 | 408 | 10614.137 | (132.911) | 10173 | 10810
| DENV-4 | 407 | 10452.469 | (205.208) | 10161 | 10772
| Total | 1633 | 10559.328 | (194.74) | 10161 | 10991