# Automatic Registration and Clustering of Time Series

###### Abstract

Clustering of time series data exhibits a number of challenges not present in
other settings, notably the problem of registration (alignment) of observed
signals. Typical approaches include pre-registration to a user-specified
template or time warping approaches which attempt to optimally align series
with a minimum of distortion. For many signals obtained from recording or
sensing devices, these methods may be unsuitable as a template signal is not
available for pre-registration, while the distortion of warping approaches may
obscure meaningful temporal information. We propose a new method for automatic
time series alignment within a _clustering problem_. Our approach, Temporal
Registration using Optimal Unitary Transformations (TROUT), is based on a
novel dissimilarity measure between time series that is easy to compute and
automatically identifies optimal alignment between pairs of time series. By
embedding our new measure in a optimization formulation, we retain well-known
advantages of computational and statistical performance. We provide an
efficient algorithm for TROUT-based clustering and demonstrate its superior
performance over a range of competitors.

Index Terms— time series clustering, convex clustering, time series
registration, majorization-minorization algorithms

University of Florida Informatics Institute
University of Florida
Gainesville, FL 32611.
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
##  1 Introduction

###  1.1 Convex Clustering

The previous two decades have witnessed a flurry of activity in statistical
machine learning and statistical signal processing, as classical methods have
been reinterpreted as convex penalized estimation problems and powerful new
theory has extended this methodology to the high-dimensional (under-sampled)
setting. While the most notable advances have been in supervised methods such
as regression and signal reconstruction, the convexity revolution has also
yielded powerful new techniques for unsupervised learning. Notably, a convex
formulation of clustering [2, 3] has received significant attention as it
allows for theoretical and computational guarantees not easily obtained for
classical clustering approaches [4, 5, 6]. Traditional convex clustering
combines a squared Frobenius norm data fidelity term, which keeps the
estimated cluster centroids near the original data, with a convex fusion
penalty, which shrinks the estimated centroids together, yielding the
following estimator:

$\hat{\bm{U}}=\operatorname*{arg\,min}_{\bm{U}\in\mathbb{R}^{n\times
p}}\frac{1}{2}\|\bm{X}-\bm{U}\|_{F}^{2}+\lambda\sum_{\begin{subarray}{c}i,j=1\\\
i<j\end{subarray}}^{n}w_{ij}\|\bm{U}_{i\cdot}-\bm{U}_{j\cdot}\|_{q}$

where $\bm{X}\in\mathbb{R}^{n\times p}$ is the data matrix,
$\hat{\bm{U}}\in\mathbb{R}^{n\times p}$ is the matrix of estimated centroids,
$\lambda>0$ is a regularization parameter controlling the overall degree of
fusion, and $\\{w_{ij}\\}$ are optional problem-specific non-negative fusion
weights. Observations $i$ and $j$ are said to belong to the same cluster if
$\hat{\bm{U}}_{i\cdot}=\hat{\bm{U}}_{j\cdot}$. In this paper, we only consider
the convex $\ell_{q}$ fusion penalty ($q=1,2,\infty$), though some research
suggests non-convex penalties also perform well [7].

Focusing our attention on the data fidelity term, we note that any distance
function can be used in lieu of the Frobenius norm without sacrificing the
advantages of the convex formulation. [8] replace the Frobenius loss with a
Huber loss to make the clustering solution robust to outliers, while [9]
propose a framework for incorporating arbitrary Bregman divergences and
exponential family-type distances in the convex clustering framework. [10] use
a Mahalanobis-type distance and embed convex clustering in a metric learning
scheme. Extending convex clustering to non-vector-valued data, [11] use an
earth-mover’s distance to cluster histogram-valued observations arising in
genomic sequencing. Building on this line of research, we propose a novel
automatic spectral-registration dissimilarity suitable for use with temporal
data and investigate its use in the clustering context.

Fig. 1: Simulation Scenarios Used in Section 4. The first scenario (top row)
consists of three signals which are clearly distinct in the time domain: a
sinusoid, a smoothed triangle wave, and a square wave. The second scenario
(middle row) consists sinusoids at three different frequencies. Third scenario
(bottom row) consists of combinations of pairs of sinusoids with different
relative phase.From left to right, the panels depict: the true cluster
centroid (true signal); a sample observation, here with $\text{SNR}=3$ and von
Mises($\kappa=10$) phase noise; centroids from DTW + hierarchical clustering
(HC), naïve time-domain HC, frequency-domain HC, and TROUT.

###  1.2 Time Series Clustering

Before introducing our proposal, we review some of the particular challenges
associated with clustering time series data. Specifically, we focus on what
[12] term “whole time-series clustering,” wherein multiple time series are
grouped into a small number of distinct clusters; they distinguish this from
“subsequence clustering,” identifying subsequences within a single long time
series, and “time point clustering,” the task of clustering individual points
in a way which respects their temporal ordering. Most commonly, whole series
clustering is performed by identifying a suitable distance measure between
series and using that distance measure within a classical clustering method
such as $k$-means or hierarchical clustering.

Perhaps the most frequently used time series distance measure is Dynamic Time
Warping (DTW), originally proposed by [13, 14]. DTW uses a dynamic programming
approach to optimally insert or remove additional time points to minimize the
Euclidean distance between the time-domain representation of two signals. DTW
is particularly well-suited for natural signals, such as those arising from
social, economic, or consumer-behavior processes, where the “shape” of the
signal is more meaningful than the exact temporal structure. The flexibility
of DTW is, for some applications, a double-edged sword so a maximum distortion
constraint of 10%, following [14], is imposed. [15] argue that this
flexibility is often quite detrimental for clustering performance and suggest
that a much lower deformation constraint be used in practice.

Alternatives to DTW include naïve (time-domain) metrics, where the temporal
structure is ignored, as well as clustering using the Fourier or wavelet
coefficients of the time series as input data. Spectral representations have
the advantage of limiting the temporal dependence and integrating well with
established signal processing techniques. See the review paper by [16] for a
thorough evaluation of time series metrics in the classification context.

##  2 Automatic Registration of Time Series

While the flexibility of DTW and related methods may be suitable for some
problems, we instead consider an alignment method tailored for stationary
signals. Specifically, we consider signals where the temporal information
_within_ a signal is meaningful, but the different signals are not aligned to
each other. Considering these signals in the spectral (Fourier) domain, only
the relative phase of the different components is meaningful, while the
absolute phase is not. We model this situation using “ _phase noise_ ,” where
a random shift is applied uniformly across the signal, corresponding to our
(random) observational window. To address phase noise, we propose a new
distance measure between which only depends on the relative phase of various
components, and not the absolute phase. As a generalization of Euclidean
distance, our method also handles observational (white) noise in the usual
manner.

Given two univariate time series, we can consider three types of distances
based on their spectral representations: i) _phase-sensitive_ , corresponding
to the standard Euclidean distance between their Fourier coefficients; ii)
_phase-oblivious_ , corresponding to the Euclidean distance between the
magnitudes of their Fourier coefficients (retaining power, but discarding
phase); and iii) _phase-adaptive_ , corresponding to the Euclidean distance
between their _optimally-aligned_ Fourier coefficients. Formally, the phase-
adaptive spectral distance can be defined as

$d_{\text{TROUT}}(\bm{u},\bm{x})=\min_{\theta\in[0,2\pi)}\|\bm{u}-e^{\mathfrak{i}\theta}\bm{x}\|_{2}.$
(1)

Compared to DTW, $d_{\text{TROUT}}$ preserves the temporal structure in the
signal, while identifying the optimal alignment between pairs of signals.

Given two matrices $\bm{X},\bm{U}\in\mathbb{C}^{n\times p}$, each row of which
represents the DFT of a separate series, the TROUT distance is defined as the
sum of row-wise distances:

$d_{\text{TROUT}}(\bm{U},\bm{X})^{2}=\sum_{i=1}^{n}d_{\text{TROUT}}(\bm{U}_{i\cdot},\bm{X}_{i\cdot})^{2}.$

Combining this matrix distance with a convex fusion penalty yields our
proposed method for simultaneous time series registration and clustering:

$\operatorname*{arg\,min}_{\bm{U}\in\mathbb{C}^{n\times
p}}\frac{1}{2}d_{\text{TROUT}}(\bm{U},\bm{X})^{2}+\lambda\sum_{\begin{subarray}{c}i,j=1\\\
i<j\end{subarray}}^{n}w_{ij}\|\bm{U}_{i\cdot}-\bm{U}_{j\cdot}\|_{q}.$ (2)

An important aspect of our formulation is that $d_{\text{TROUT}}$ implicitly
re-aligns the series at each step. Compared with standard pre-alignment
techniques, this can significantly improve performance and interpretability.
Specifically, $d_{\text{TROUT}}$ works by aligning each estimated centroid
(row of $\bm{U}$) with the observed data. This is preferable to pre-alignment
as it avoids the difficult question of selecting a registration target, and
aligns a centroid only to that data in its cluster. Similar advantages of
dynamic alignment have been observed in the shape analysis context by [17],
though our approach is simpler to implement than theirs. We note that
$d_{\text{TROUT}}$ fails to satisfy the triangle inequality and hence is not a
proper distance measure. As such, Problem (2) is not generally convex without
further assumptions on $\bm{X}$.

As a generalization of Euclidean convex clustering, TROUT clustering inherits
many of the attractive properties of convex clustering including: a continuous
solution path, indexed by $\lambda$, which smoothly interpolates from $n$ to
$1$ clusters [3]; the ability to quickly identify local optima [18]; the
ability to construct dendrogram representations of the solution path [19];
provable statistical consistency [4, 5, 6]; and a high-degree of robustness to
noise [9, 8]. While our focus is on the use of $d_{\text{TROUT}}$ in the
fusion clustering context, we also note that the TROUT distance, and
associated notion of TROUT Fréchet mean, can also be used within $K$-means or
hierarchical clustering.

Fig. 2: Comparitive performance of TROUT clustering against dynamic time
warping (DTW), time-domain, and frequency-domain based hierarchical clustering
(HC). The same scenario designs as shown in Figure 1 are used, with the degree
of phase noise–the concentration parameter $\kappa$ in a von Mises
distribution–and observational noise (SNR) varied. As can be clearly seen,
TROUT outperforms all methods in the high-phase noise (low $\kappa$) setting:
this outperformance is most pronounced in Scenario 3 (bottom row) where the
signals differ only in the relative phase of components. Time-domain
clustering with oracle alignment performs as well as frequency-domain
clustering, as expected.

##  3 Computational Considerations

The TROUT distance (1) is defined by a one-dimensional optimization problem.
This problem is non-convex in $\theta$, due to the sinusoidal behavior of the
complex exponential, but simple enough that bisection or grid-search methods
can be used to identify the minimum. Still, the non-convex problem may pose
difficulties for composite optimization schemes [20]. Instead, we recast
$d_{\text{TROUT}}$ as a complex-valued optimization problem

$d_{\text{TROUT}}(\bm{u},\bm{x})=\min_{z\in\mathbb{B}_{\mathbb{C}}(r=1)}\|\bm{u}-z\,\bm{x}\|_{2},$
(3)

where $\mathbb{B}_{\mathbb{C}}(r=1)$ is the complex unit sphere. This problem
is essentially complex least squares with a norm constraint and has closed-
form minimizer
$z=\angle[\bm{x}^{H}\bm{u}]=\bm{x}^{H}\bm{u}/\|\bm{x}\|\|\bm{u}\|$. In
multivariate settings, $z$ is optimized over unitary transforms, justifying
the TROUT acronym.

Now that we can calculate $d_{\text{TROUT}}$ efficiently, we can develop
algorithms for the TROUT clustering problem (2). Following [18] and [19], we
use an alternating direction method of multipliers approach, whose updates are
given by

$\displaystyle\bm{U}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{U}\in\mathbb{C}^{n\times
p}}\frac{1}{2}d_{\text{TROUT}}(\bm{U},\bm{X})^{2}+$
$\displaystyle\phantom{\operatorname*{arg\,min}}\quad\frac{\rho}{2}\left\|\bm{D}\bm{U}-\bm{V}^{(k)}+\bm{Z}^{(k)}\right\|_{F}^{2}$
$\displaystyle\bm{V}^{(k+1)}$
$\displaystyle=\operatorname{\textsf{prox}}_{\lambda/\rho
P(\cdot;\bm{w},q)}(\bm{D}\bm{U}^{(k+1)}+\bm{Z}^{(k)})$
$\displaystyle\bm{Z}^{(k+1)}$
$\displaystyle=\bm{Z}^{(k)}+\bm{D}\bm{U}^{(k+1)}-\bm{V}^{(k+1)}$

where $\bm{D}$ is a suitable difference matrix and $P(\cdot;\bm{w},q)$ is the
weighted convex fusion penalty function. For $q=1,2$, the proximal operator
yields a simple closed-form update for $\bm{V}$, while for $q=\infty$, the
efficient algorithm of [21] can be used.

Unlike standard convex clustering, the primal ($\bm{U}$) update does not have
a closed-form update and instead requires an iterative solver. Note that, for
two vectors $\bm{u}$ and $\bm{x}$, $d_{\text{TROUT}}(\bm{u},\bm{x})$ can be
equivalently expressed using the projection of $\bm{u}$ onto the set of all
“phase-adjusted” $\bm{x}$ vectors: that is,

$d_{\text{TROUT}}(\bm{u},\bm{x})=\|\bm{u}-P_{\tilde{\bm{x}}}(\bm{u})\|_{2}\text{
where }\tilde{\bm{x}}=\\{e^{\mathfrak{i}\theta}\bm{x}:\theta\in[0,2\pi)\\}$

and where $P_{\tilde{\bm{x}}}(\cdot)$ denotes the projection onto the set
$\tilde{\bm{x}}$. This sort of “distance-to-set” function was popularized in
machine learning by [22] and is well-behaved (single-valued and
differentiable) for almost all $\bm{u},\bm{x}$ [23, Sections 18.6 and 24.3].
With this formulation, we can use a majorization-minimization (MM) approach
[24] for the $\bm{U}$ subproblem. In particular, note that

$\frac{1}{2}d_{\text{TROUT}}(\bm{U},\bm{X})^{2}\leq\frac{1}{2}\|\bm{U}-P_{\mathcal{U}^{\\{t\\}}}(\bm{X})\|_{F}^{2}$

where $\mathcal{U}^{\\{t\\}}$ is the “phase-adjusted” equivalence class of the
previous iterate $\bm{U}^{\\{t\\}}$ and $P_{\mathcal{U}^{\\{t\\}}}(\cdot)$
denotes projections thereon. This majorizes the TROUT distance because
$P_{\mathcal{U}^{\\{t\\}}}(\bm{X})\in\mathcal{U}^{\\{t\\}}$ by construction
and the distance to a set is always less than or equal to the distance to a
fixed point within that set. This yields the inner iterates:

$\displaystyle\tilde{\bm{U}}^{\\{t+1\\}}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{U}\in\mathbb{C}^{n\times
p}}\frac{1}{2}\left\|\bm{U}-P_{\mathcal{U}^{\\{t\\}}}(\bm{X})\right\|_{F}^{2}+$
$\displaystyle\phantom{\operatorname*{arg\,min}}\quad\frac{\rho}{2}\left\|\bm{D}\bm{U}-\bm{V}^{(k)}+\bm{Z}^{(k)}\right\|_{F}^{2}$

which have a closed-form solution

$\tilde{\bm{U}}^{\\{t+1\\}}=(\bm{I}+\rho\bm{D}^{T}\bm{D})^{-1}(P_{\mathcal{U}^{\\{t\\}}}(\bm{X})+\rho\bm{D}^{T}(\bm{V}^{(k)}-\bm{Z}^{(k)})).$

Though non-convex, the $d_{\text{TROUT}}(\cdot,\bm{X})$ function is well-
behaved and global convergence to a stationary point follows from Scenario
1(a) in Theorem 1 of [Wang:2019].

Additional performance improvements can be obtained by taking only a single MM
step within each ADMM iteration [25], though this is typically unnecessary as
the inner iterates tend to converge quite quickly. If the $\mathcal{O}(n^{3})$
Cholesky factorization of $\bm{I}+\rho\bm{D}^{T}\bm{D}$ is cached, the per
iterate cost of this approach is $\mathcal{O}(n^{2}p)$, dominated by two
rounds of triangular back-substitution, which scales efficiently with high
sampling rates. For the simulations considered in the next section
($\bm{X}\in\mathbb{C}^{60\times 64}$), this algorithm solves Problem (2) to
numerical precision in an average of 0.0825 seconds per value of $\lambda$ in
a warm start scheme on a low-end desktop with 16 GB of RAM and an Intel
i7-9700 processor running at 3.0 GHz. When combined with the “one-step”
algorithmic regularization approach of [19], performance is increased by an
order of magnitude (0.000665 seconds per $\lambda$) without any noticeable
drop in clustering accuracy.

##  4 Simulation Study

In this section, we compare TROUT clustering with various time- and frequency-
domain clustering methods. We make repeated use of three simulation designs,
shown in the left column of Figure 1. These three families highlight different
challenges in time series clustering. The first scenario (top row) includes
both a square wave and a (smoothed) triangle wave, which have slowly decaying
spectra and exhibit Gibbs phenomena (rapid variation overshooting the target
function) at their discontinuities; because of the slow spectrum decay, these
signals have more relevant features than observations in each cluster, a
situation which is well-known to challenge many clustering methods, while the
discontinuities in the time-domain challenge warping or smoothing approaches.
The second scenario (middle row) consists of sinusoids of different
frequencies: these signals are easily separated in the frequency domain, but
time-domain approaches struggle to identify structure. The third scenario
(bottom row) is the most challenging: all three signals consist of pairs of
sinusoids, differing only in relative phase. Power-only approaches are unable
to distinguish these signals, and phase-noise poses a particular challenge for
both time- and frequency-domain methods.

From left-to-right, Figure 1 shows the true cluster centroids; typical
observations generated with a signal-to-noise ratio of 3 and with phase noise
drawn from a mean zero von Mises random variable with concentration
$\kappa=10$; centroids estimated by dynamic time warping followed by
hierarchical clustering (HC); centroids estimated by HC in the time-domain;
centroids estimated by HC in the frequency-domain; and centroids estimated by
TROUT. In each setting, signals are observed at 128 Hz, with 20 samples
generated from each cluster for a total of 60 observations, yielding a
spectral representation $\bm{X}_{C}\in\mathbb{C}^{60\times 64}$ or a time-
domain representation $\bm{X}_{R}\in\mathbb{R}^{60\times 128}$. For
legibility, a $5$-nearest neighbor smoother is super-imposed. While both TROUT
and frequency-domain HC accurately separate clusters, only TROUT is able to
correctly recover the shape of the true centroids (up to phase shift).

Figure 2 uses the same scenarios, but varies the signal to noise ratio and the
degree of phase noise. The sampling rate and number of observations is fixed
as before. The adjusted Rand index used to measure clustering performance, is
averaged over 50 replicates, with the oracle number of clusters made available
to all methods. As can be seen, TROUT and frequency-domain HC are consistently
among the best performing methods in all regimes, with the advantages of TROUT
increasing as the phase noise increases ($\kappa$ decreases). For
$\kappa=\infty$, corresponding to the case of no phase noise, the additional
flexibility of TROUT registration slightly decreases performance. Focusing on
the third scenario (bottom row) wherein signals vary only by their phase, we
note that frequency domain HC is unable to distinguish signals in the high-
phase noise setting, even with a high SNR. Time-domain methods including naïve
HC and DTW are unable to correctly discover true clusters in this scenario at
any SNR.

##  5 Discussion and Extensions

We have proposed a new approach to time series clustering which combines a
convex fusion penalty with automatic registration. We have derived a closed-
form solution for the optimal spectral-domain alignment and have used it to
develop an efficient MM algorithm which attains a local optimum provably and
efficiently. Simulation studies confirm the efficiency and stability of our
approach across several regimes. Our method is particularly well-suited for
situations of high “phase noise,” as would be expected when clustering
naturally occurring stationary signals, though it also performs competitively
in more standard domains. In the simulations shown in Section 4, we only
considered basic spectral analysis (whole-sample DFT without tapering or
smoothing), though our method is readily applied within more robust analytical
pipelines.

The heart of our proposed method is the automatic and optimal alignment
induced by the TROUT distance. While we have focused on clustering stationary
univariate signals, the TROUT distance can be used anywhere a standard
$\ell_{2}$ or Frobenius distance is used in machine learning. The basic
framework of TROUT-based clustering can also be extended in several useful
directions, including automatic amplitude adjustment, frequency isolation /
band identification, and multivariate time series clustering. We conclude with
a brief discussion of each of these extensions.

As defined above, the TROUT distance (1) only allows for phase adjustment. For
signals which may be observed on different scales, where only intra-signal
relative magnitude is meaningful, the complex form of the TROUT distance (3)
can be relaxed to allow arbitrary complex $z$. As with standard TROUT, the
relaxed problem has closed-form minimizer $z=\bm{x}^{H}\bm{u}/\|\bm{x}\|^{2}$
and the associated clustering problem can be analyzed with a MM algorithm.
Notably, this relaxed form has better convexity properties than the method
discussed above. In certain problems, it may be useful to isolate certain
frequencies which are most important to the clustering solution; this can be
achieved by adding an $\ell_{1}$ or mixed $\ell_{2}/\ell_{1}$ penalty to the
columns of $\bm{U}$, as explored by [9]. If, instead, it is useful to identify
bands of frequencies which behave similarly, a TROUT-based form of convex bi-
clustering [26] could be used to simultaneously cluster both observations and
frequencies into clusters and meaningful bands respectively. Finally, if we
instead seek to cluster $p$-variate time series, our data $\bm{X}$ would now
be an $n\times p\times p$ complex higher-order array (tensor) which we would
seek to cluster along the first mode: in this case, the TROUT problem now
consists of alignment by an arbitrary unitary transform. This problem can be
solved via a unitary Procrustes problem, which has a closed-form solution in
terms of the SVD of $\bm{X}^{H}\bm{U}$ and can be easily embedded in efficient
tensor clustering algorithms [27]. We plan to explore several of these
extensions in future work.

##  6 References

## References

  * [1] Kristiaan Pelckmans, Joseph Brabanter, Bart Moor and Johan Suykens “Convex Clustering Shrinkage” In _PASCAL Workshop on Statistics and Optimization of Clustering_ , 2005
  * [2] Fredrik Lindsten, Henrik Ohlsson and Lennart Ljung “Clustering using sum-of-norms regularization: With application to particle filter output computation” In _SSP 2011: Proceedings of the 2011 IEEE Statistical Signal Processing Workshop_ Nice, France: Curran Associates, Inc., 2011, pp. 201–204 DOI: 10.1109/SSP.2011.5967659
  * [3] Toby Dylan Hocking, Armand Joulin, Francis Bach and Jean-Philippe Vert “Clusterpath: An Algorithm for Clustering using Convex Fusion Penalties” In _ICML 2011: Proceedings of the 28th International Conference on Machine Learning_ Bellevue, Washington, USA: ACM, 2011, pp. 745–752 URL: http://icml-2011.org/papers/419_icmlpaper.pdf
  * [4] Changbo Zhu, Huan Xu, Chenlei Leng and Shuicheng Yan “Convex Optimization Procedure for Clustering: Theoretical Revisit” In _NIPS 2014: Advances in Neural Information Processing Systems 27_ Montréal, Canada: Curran Associates, Inc., 2014, pp. 1619–1627 URL: https://papers.nips.cc/paper/5307-convex-optimization-procedure-for-clustering-theoretical-revisit
  * [5] Kean Ming Tan and Daniela Witten “Statistical Properties of Convex Clustering” In _Electronic Journal of Statistics_ 9.2, 2015, pp. 2324–2347 DOI: 10.1214/15-EJS1074
  * [6] Peter Radchenko and Gourab Mukherjee “Convex Clustering via $\ell_{1}$ Fusion Penalization” In _Journal of the Royal Statistical Society, Series B: Statistical Methodology_ 79.5, 2017, pp. 1527–1546 DOI: 10.1111/rssb.12226
  * [7] Yuliya Marchetti and Qing Zhou “Solution Path Clustering with Adaptive Concave Penalty” In _Electronic Journal of Statistics_ 8.1, 2014, pp. 1569–1603 DOI: 10.1214/14-EJS934
  * [8] Chenyu Liu, Qiang Sun and Kean Ming Tan “Robust convex clustering: How does fusion penalty enhance robustness?” In _ArXiv Pre-Print 1906.09581_ , 2019 URL: https://arxiv.org/abs/1906.09581
  * [9] Minjie Wang and Genevera I. Allen “Integrative Generalized Convex Clustering Optimization and Feature Selection for Mixed Multi-View Data” In _ArXiv Pre-Print 1912.05449_ , 2019 URL: https://arxiv.org/abs/1912.05449
  * [10] Xiaopeng Lucia Sui, Li Xu, Xiaoning Qian and Tie Liu “Convex clustering with metric learning” In _Pattern Recognition_ 81, 2018, pp. 575–584 DOI: 10.1016/j.patcog.2018.04.019
  * [11] Cheolwoo Park et al. “Convex Clustering Analysis for Histogram-Valued Data” In _Biometrics_ 75.2, 2019, pp. 603–612 DOI: 10.1111/biom.13004
  * [12] Saeed Aghabozorgi, Ali Seyed Shirkhorshidi and Ten Ying Wah “Time-Series Clustering - A Decade Review” In _Information Systems_ 53, 2015, pp. 16–38 DOI: 10.1016/j.is.2015.04.007
  * [13] H. Sakoe and S. Chiba “A Dynamic Programming Approach to Continuous Speech Recognition” In _Proceedings of the Seventh International Congress on Acoustics_ 3 - Subjective Acoustics, 1971, pp. 65–69 URL: https://www.icacommission.org/proceedg.html
  * [14] H. Sakoe and S. Chiba “Dynamic programming algorithm optimization for spoken word recognition” In _IEEE Transactions on Acoustics, Speech, and Signal Processing_ 26.1, 1978, pp. 43–49 DOI: 10.1109/TASSP.1978.1163055
  * [15] Chotirat Ann Ratanamahatana and Eamonn Keogh “Three Myths about Dynamic Time Warping Data Mining” In _SDM 2005: Proceedings of the 2005 SIAM International Conference on Data Mining_ , 2005, pp. 506–510 DOI: 10.1137/1.9781611972757.50
  * [16] Joan Serrà and Josp Ll. Arcos “An empirical evaluation of similarity measures for time series classification” In _Knowledge-Based Systems_ 64, 2014 DOI: 10.1016/j.knosys.2014.04.035
  * [17] Anuj Srivastava, Eric Klassen, Shantanu H. Joshi and Ian H. Jermyn “Shape Analysis of Elastic Curves in Euclidean Spaces” In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 33.7, 2011 DOI: 10.1109/TPAMI.2010.184
  * [18] Eric C. Chi and Kenneth Lange “Splitting Methods for Convex Clustering” In _Journal of Computational and Graphical Statistics_ 24.4 Taylor & Francis, 2015, pp. 994–1013 DOI: 10.1080/10618600.2014.948181
  * [19] Michael Weylandt, John Nagorski and Genevera I. Allen “Dynamic Visualization and Fast Computation for Convex Clustering via Algorithmic Regularization” In _Journal of Computational and Graphical Statistics_ 29.1, 2020, pp. 87–96 DOI: 10.1080/10618600.2019.1629943
  * [20] John C. Duchi and Feng Ruan “Stochastic Methods for Composite and Weakly Convex Optimization Problems” In _SIAM Journal on Optimization_ 28.4, pp. 3229–3259 DOI: 10.1137/17M1135086
  * [21] Laurent Condat “Fast Projection onto the Simplex and the $\ell_{1}$ Ball” In _Mathematical Programming_ 158, 2016, pp. 575–585 DOI: 10.1007/s10107-015-0946-6
  * [22] Jason Xu, Eric Chi and Kenneth Lange “Generalized Linear Model Regression under Distance-to-Set Penalties” In _NIPS 2017: Advances in Neural Information Processing Systems 30_ , 2017 URL: http://papers.nips.cc/paper/6737-generalized-linear-model-regression-under-distance-to-set-penalties
  * [23] Heinz H. Bauschke and Patrick L. Combettes “Convex Analysis and Monotone Operator Theory in Hilbert Spaces”, CMS Books in Mathematics Springer, 2017 DOI: 10.1007/978-3-319-48311-5
  * [24] Kenneth Lange “MM Optimization Algorithms”, Other Titles in Applied Mathematics Society for IndustrialApplied Mathematics, 2016 DOI: 10.1137/1.9781611974409
  * [25] Liang Chen, Xudong Li, Defeng Sun and Kim-Chuan Toh “On the equivalence of inexact proximal ALM and ADMM for a class of convex composite programming” In _Mathematical Programming_ , 2019 DOI: 10.1007/s10107-019-01423-x
  * [26] Eric C. Chi, Genevera I. Allen and Richard G. Baraniuk “Convex Biclustering” In _Biometrics_ 73.1 Wiley Online Library, 2017, pp. 10–19 DOI: 10.1111/biom.12540
  * [27] Michael Weylandt “Splitting Methods For Convex Bi-Clustering And Co-Clustering” In _DSW 2019: Proceedings of the 2 nd IEEE Data Science Workshop_ Minneapolis, Minnesota: IEEE, 2019, pp. 237–244 DOI: 10.1109/DSW.2019.8755599

## References

  * [28] Peter H. Schönemann “A generalized solution of the orthogonal Procrustes problem” In _Psychometrika_ 31.1, 1966, pp. 1–10 DOI: 10.1007/BF02289451
  * [29] Lars Eldén and Haesun Park “A Procrustes problem on the Stiefel manifold” In _Numerische Mathematik_ 82.4, 1999, pp. 599–619 DOI: 10.1007/s002110050432

## Supplementary Materials

##  Appendix A Proofs

As discussed in Section 3, the TROUT distance between two vectors can be
computed by solving the complex-valued (Wirtinger) optimization problem

$d_{\text{TROUT}}(\bm{u},\bm{x})=\min_{z\in\mathbb{B}_{\mathbb{C}}(r=1)}\|\bm{u}-z\bm{x}\|_{2}$

If we consider the column vectors $\bm{u}$ and $\bm{x}$ as $p\times 1$
matrices and $z$ as a $1\times 1$ matrix, we note that this is a special case
of the _Unitary_ Procrustes problem:

$\operatorname*{arg\,min}_{\bm{Z}\in\mathfrak{U}_{k}}\|\bm{U}-\bm{X}\bm{Z}\|_{F}^{2}$

where $\mathfrak{U}_{k}$ denotes the unitary group of dimension $k$; that is,
the set of $k\times k$ matrices satisfying
$\bm{Z}^{H}\bm{Z}=\bm{Z}\bm{Z}^{H}=\bm{I}_{k\times k}$. The unitary Procrustes
problem has a closed-form solution given by

$\bm{A}\bm{B}^{H}\text{ where
}\bm{A}\bm{D}\bm{B}^{H}=\textsf{SVD}(\bm{X}^{H}\bm{U})$

(The Procrustes problem is typically presented with $\bm{Z}\bm{X}$ instead of
$\bm{X}\bm{Z}$: in that case, the solution uses the SVD of
$\bm{U}\bm{X}^{H}$.) The unitary Procrustes problem does not often appear in
the literature _per se_ , but analysis of the orthogonal (real-valued)
Procrustes problem translates directly to the complex-domain [28, 29]. With
this in hand, we are able to give a closed-form solution for the TROUT
minimizer.

######  Proposition 1.

Let $\bm{u},\bm{x}\in\mathbb{C}^{p}$ be arbitrary complex vectors. Then

$\angle[\bm{x}^{H}\bm{u}]=\frac{\bm{x}^{H}\bm{u}}{\|\bm{x}\|\|\bm{u}\|}=\operatorname*{arg\,min}_{z\in\mathbb{B}_{\mathbb{C}}(r=1)}\|\bm{u}-z\bm{x}\|_{2}$

and

$d_{\text{TROUT}}(\bm{u},\bm{x})=\min_{\theta\in[0,2\pi)}\|\bm{u}e^{\mathfrak{i}\theta}\bm{x}\|=\left\|\bm{u}-\angle[\bm{x}^{H}\bm{u}]\bm{x}\right\|$

###### Proof.

Using the closed-form solution for the Procrustes problem given above, it
suffices to calculate the SVD of $\bm{x}^{h}\bm{u}$. Note that we have the
analytic SVD

$\bm{x}^{H}\bm{u}=\begin{pmatrix}\frac{\bm{x}}{\|\bm{x}\|}\end{pmatrix}\begin{pmatrix}\|\bm{u}\|\bm{x}\|\end{pmatrix}\begin{pmatrix}\frac{\bm{u}}{\|\bm{u}\|}\end{pmatrix}^{H}$

so the Procrustes solution is given by

$\begin{pmatrix}\frac{\bm{x}}{\|\bm{x}\|}\end{pmatrix}\begin{pmatrix}\frac{\bm{u}}{\|\bm{u}\|}\end{pmatrix}^{H}=\angle[\bm{x}^{H}\bm{u}]$

as claimed. The second claim follows by substitution in the definition of
$d_{\text{TROUT}}$. ∎

##  Appendix B Additional References