Clustering of Series via Dynamic Mode Decomposition and the Matrix Pencil Method
Leonid Pogorelyuk Clarence W. Rowley


Dynamic magnetism in the disordered hexagonal double perovskite BaTi_1/2Mn_1/2O_3
    M.R. Cantarino^1, R.P. Amaral^2, R.S. Freitas^1, J.C.R.
Araújo^2, R. Lora-Serrano^2, H. Luetkens^3, C. Baines^3,
S. Bräuninger^4, V. Grinenko^4, R. Sarkar^4, H.H. Klauss^4,
E. C. Andrade^5, F. A. Garcia^1 
    December 30, 2023
=================================================================================================================================================================================================================




In this paper, a new algorithm for extracting features from sequences of multidimensional observations is presented. The independently developed Dynamic Mode Decomposition and Matrix Pencil methods provide a least-squares model-based approach for estimating complex frequencies present in signals as well as their corresponding amplitudes. Unlike other feature extraction methods such as Fourier Transform or Autoregression which have to be computed for each sequence individually, the least-squares approach considers the whole dataset at once. It invokes order reduction methods to extract a small number of features best describing all given data, and indicate which frequencies correspond to which sequences. As an illustrative example, the new method is applied to regions of different grain orientation in a Transmission Electron Microscopy image.




§ INTRODUCTION


Clustering of series of data points is an unsupervised classification
task which consists of grouping sequences together based on some notion
of similarity between them. In the literature this task is known as “clustering
of time series” or “functional data clustering”, although it
applies to any sequences, not just sequences in time <cit.>.
These methods belong to the field of cluster analysis, which has an
even larger variety of algorithms and techniques <cit.>.

It is common to divide the various time series clustering algorithms
into several approaches <cit.>: The raw-data-based
approaches compute similarity or distance metrics directly on the
time series and then apply more general clustering algorithms (e.g.
k-means <cit.>) to achieve their goal. Feature-based
approaches first compute features using methods such as the Fast Fourier Transform (FFT)
<cit.> or Principal Component Analysis (PCA) <cit.>
on the time series and then apply standard clustering techniques. Finally, model-based
approaches assume some model for the process producing the data series
and compute its parameters.

A common model-based approach assumes some linear relation between
elements of the series, their precursors and process noise <cit.>.
In that case, one fits autoregression coefficients to the series first
<cit.>, and then proceeds by computing a distance metric
between time series and clustering based on those coefficients. Recently,
<cit.> has introduced another metric for time series, which assumes a linear relation between
elements within each series, and relies on a method called Dynamic Mode Decomposition
(DMD) <cit.>.

DMD <cit.> and its variants assume that the data
was generated by a linear dynamical system with measurement noise; the method then estimates the complex frequencies and magnitudes (called “modes”)
corresponding to that system. For sequential data, DMD is often paired with
delay-embedding, a common technique in system identification <cit.>.
It involves defining high dimensional “delayed observables” which
include several adjacent data points from the given sequence, arranging
the data in (generalized) Hankel matrices and then applying standard
DMD <cit.>.

Interestingly, Hankel matrices were employed in the literature of
frequency estimation long before the appearance of DMD (for a review
see <cit.>). Modern algorithms for frequency
estimation include the State Space method <cit.>, ESPRIT
<cit.> and the the Matrix Pencil method <cit.>.
Similarly to DMD, they all compute the Singular Value Decomposition
(SVD) of the Hankel data matrices and, under certain conditions to
be discussed in this paper, give the same results.

Furthermore, the various methods—Matrix Pencil, ESPRIT, and DMD—each provide
a different perspective on the same arrangement of given data series.
Together, they offer an insight on designing a new set of features
for sequential data which will be the main subject of this paper.
Those features, computed on a set of series, are particularly useful for
clustering purposes.

In section <ref>, we give an overview of DMD, Matrix Pencil methods,
and ESPRIT, emphasizing their similarities. Based
on these, an algorithm for extracting features from a set of sequences
in described in section <ref>. Finally, numerical
examples for the use of those features for clustering are presented
in section <ref>.



§ COMPARISON BETWEEN DMD AND MATRIX PENCIL
METHODS


Throughout our discussion, we consider data consisting of measurements
y(t) of a signal x(t) in the presence of noise s(t),
where t=0,1,…,T.  We also assume that the signal is a linear combination of R damped sinusoids, so that it may be written


    x(t)   =j=1R∑v_jλ_j^t
    
      y(t)   =x(t)+s(t),

where x(t),y(t),s(t),v_j∈ℂ^n, λ_j∈ℂ. In this section we discuss two families of methods
which estimate R and λ_j,v_j based on the data { y(t)} _t=0^T.

Earlier frequency estimation methods first approximate the complex frequencies
λ_1,…,λ_R, and leave the task of finding the coefficients
v_1,…,v_R for a standard least squares procedure. Among these
are the ESPRIT <cit.>, the Matrix Pencil <cit.>
and State Space <cit.> methods, which will be discussed
in sections <ref>–<ref>.

The Dynamic Mode Decomposition (DMD) <cit.> and
its variants estimate both the frequencies and the (vector) coefficients v_1,…,v_R
(called “modes”) and are frequently applied to non-sequential
data as well. In the context of sequential data, DMD is often used
with delayed observables <cit.>,
which will be defined in sections <ref>-<ref>.

This section will focus on showing that DMD with delayed observables
is equivalent to the Matrix Pencil methods which historically preceded
it, and will consider a possible insight from the ESPRIT method to
the DMD community.




 §.§ DMD with Delayed Observables and Order Reduction


The Dynamic Mode Decomposition (DMD) was originally
  proposed in the fluid mechanics community, as a method for identifying coherent structures in fluids flows. In DMD setting, each pair of
measurements (y(t),y(t)) from (<ref>)
is called a data snapshot, with the indices t not necessarily
being ordered or sequential.  The general idea behind DMD is to find a linear map such that


    y(t+1)≈ Ly(t),

with L∈ℂ^n× n. However, since L has at most
rank n, this approximation is extremely poor when n<R even if
the measurements are exact, i.e. s(t)=0, ∀ t. To alleviate
this problem in the case when sequential data is available, it is
common to introduce delayed observables <cit.>, as


    z(t)=[   y(t); y(t+1);      ⋮; y(t+d) ]∈ℂ^(d+1)n,

with d>0 delays and (d+1)n≥ R. When no noise is
present (s(t)=0), given enough data, one can fit an autoregressive
model <cit.>, such that


    z(t)=Kz(t),

holds exactly for some K∈ℂ^(d+1)n×(d+1)n.
In this case, R eigenvalues of K are λ_1,…,λ_R, and
the rest are zero.

If R is not precisely known but (d+1)n>R, it is necessary
to get an estimate of the number of frequencies, R̅, and project
the delayed observables z(t) onto an R̅<(d+1)n
dimensional space,


    z̃(t)=Tz(t),

where z̃(t) are the reduced order observables and T∈ℂ^R̅×(d+1)n
has rank R̅. In the DMD literature it is customary to compute the
reduced-order observables using
Proper Orthogonal Decomposition (POD)
<cit.> (also known as Principal Component
Analysis (PCA) <cit.>), in which T is
computed using Singular Value
Decomposition (SVD), as described in the following subsection.

Consequently, considering the data pair (z̃(t),z̃(t+1)),
the fitted linear map


    z̃(t+1)=K̃z̃(t),

with K̃∈ℂ^R̅×R̅ will, assuming
R̅=R and s(t)=0, have exactly the desired eigenvalues
λ_1,…,λ_R.

Besides estimating the frequency, one can decompose the signal into
a sum,


    z(t)=j=1R∑[      v_j;   λ_jv_j;        ⋮; λ_j^dv_j ]λ_j^t+[   s(t); s(t+1);      ⋮; s(t+d) ],

of so called dynamic modes, [      v_j^T        ... λ_j^dv_j^T ]^T, and a remainder, [   s(t)^T        ⋯ s(t+d)^T ]^T. Those DMD modes are eigenvectors of K in (<ref>)
and a discussion on their scaling can be found in <cit.>.

We next present the full implementation details of the DMD method.




 §.§ DMD implementation


Given M pairs of snapshots (z(t),z(t+1)), the data matrices are defined such that their columns consist of the delayed observables
in no particular order,


    X=[   z(t_0)        ⋯ z(t_M-1) ],    
    Y=[   z(t_0+1)          ⋯ z(t_M-1+1) ]
    .

For a sequential signal, one may choose t_i=i and construct the data matrices,


    X=[   y(0)   y(1)      ⋯;   y(1)   y(2);      ⋮      ⋮;   y(d) y(d+1)      ⋯ ]
    , Y=[   y(1)   y(2)      ⋯;   y(2)   y(3);      ⋮      ⋮; y(d+1) y(d+2)      ⋯ ]
    ,

to be used in the following algorithm.

Step 1: Compute the Singular Value Decomposition (SVD) 

    X=UΣ V^*,

and choose the largest R̅ singular values after which there
is a significant drop in their magnitude. The truncated SVD approximation
of the data is then X≈ U_R̅Σ_R̅V_R̅^*
with U_R̅∈ℂ^(d+1)n×R̅,
Σ_R̅∈ℂ^R̅×R̅, V_R̅∈ℂ^M×R̅,
and T=U_R^* is the order reducing transformation.

Step 2: Defining the reduced order data matrices as
X̃=U_R̅^*X, Ỹ=U_R̅^*Y both
in ℂ^R̅× M, we have


    K̃=ỸX̃^+=U_R̅^*YV_R̅Σ_R̅^-1,

where K̃∈ℂ^R̅×R̅.

Step 3: Perform an eigendecomposition of K̃


    K̃=ṼΛ̅Ṽ^-1,

where columns of Ṽ∈ℂ^R̅×R̅ (Ṽ^-1∈ℂ^R̅×R̅)
are the right (left) eigenvectors of K̃ and Λ̅=diag{λ̅_1R̅,...,λ̅_R̅}
consists of the approximated complex frequencies.

Step 4: A rank R̅ approximation of the full dynamics
in (<ref>) is 

    K≈ U_R̅K̃U_R̅^*.

Its right eigenvectors, the columns of U_R̅Ṽ, are
called the DMD modes, while the rows of Ṽ^-1U_R̅^*
are the left eigenvectors and are known as the adjoint DMD modes <cit.>.

Step 5: To find an appropriate scaling for the DMD modes we
first notice that for any t,τ it follows from (<ref>),
(<ref>) and (<ref>)
that


    z(t)≈ K^t-τz(τ)≈j=1R̅∑c_τ,jw_jλ̅_j^t-τ,

where w_j=(U_R̅Ṽ)_j∈ℂ^(d+1)n
is the jth DMD mode and c_τ,j=(Ṽ^-1U_R̅^*z(τ))_j∈ℂ
are scaling coefficients based on the τth snapshot.

In general, since M>R̅, one has c_τ_1,jλ̅_j^-τ_1 c_τ_2,jλ̅_j^-τ_2
(unlike in the underconstrained case discussed in <cit.>),
hence we suggest a scaling for the jth DMD mode based on an average
of its scaling coefficients at all time snapshots,


    c_j=1/Mτ∈{ t_1,...,t_M}∑(Ṽ^-1U_R̅^*z(τ))_jλ̅_j^-τ.

This gives an approximation of the dynamics in delayed observables
(see (<ref>)),


    z(t)≈j=1R̅∑c_jw_jλ̅_j^t,

where c_jw_j∈ℂ^(d+1)n are the scaled
dynamic modes of the delayed observable. The coefficients v_j
in (<ref>) are then estimated via


    v̅_j=1/d+1k=0d∑(c_jw_j)_kn:(k+1)nλ̅_j^-k,

where (c_jw_j)_kn:(k+1)n∈ℂ^n
are the n consequent elements of c_jw_j beginning at the
knth elements.




 §.§ Matrix Pencil Methods

The authors of the Matrix Pencil <cit.> and the State Space <cit.> methods originally considered the problem of retrieving parameters of sinusoidal processes (with frequencies close to one another) from noisy measurements. There, the data matrices are defined as 


    X=[      ⋯   y(1)   y(0);          y(2)   y(1);             ⋮      ⋮;      ⋯ y(d+1)   y(d) ], Y=[      ⋯   y(2)   y(1);          y(3)   y(2);             ⋮      ⋮;      ⋯ y(d+2) y(d+1) ],

which correspond to a reversed ordering with respect to (<ref>). The Matrix Pencil method proceeds by finding the generalized eigenvectors
of the matrix pencil (X,Y), i.e. λ̅∈ℂ
and p∈ℂ^M, q∈ℂ^(d+1)n in
the row and column spaces of X respectively, such that


    (Y-λ̅X)p=0,



    q^*(Y-λ̅X)=0.

Similarly to DMD, this is accomplished by computing the truncated
SVD of X, X≈ U_R̅Σ_R̅V_R̅^*.
The generalized eigenvalues λ̅ are estimated by constructing
the matrix <cit.>


    L̃=Σ_R̅^-1U_R̅^*YV_R̅

and computing its eigenvalues.

Indeed, if L̃=W̃Λ̅^#W̃^-1
where W̃,Λ̅^#∈ℂ^R̅×R̅
and Λ̅^# is diagonal, the matrix pencil (X,Y)
can be approximated by


    Y-λ X≈ U_R̅Σ_R̅(L̃-λ I)V_R̅^*.

Therefore, the columns of V_R̅W̃ are its right generalized
eigenvectors and the rows of W̃^-1Σ_R̅^-1U_R̅^*
are its left generalized eigenvectors.

When comparing to the DMD matrix in (<ref>) and (<ref>),
we conclude that L̃=Σ_R̅K̃Σ_R̅^-1,
Λ̅^#=Λ̅ and W̃=Σ_R̅^-1Ṽ.
In other words, DMD and Matrix Pencil produce the same estimates for
eigenvalues and the reduced order matrices are similar (with Σ_R̅
being the similarity transformation).

Finally we find that the estimates of both the adjoint DMD modes and
the left generalized eigenvectors of (X,Y) are the same,


    Ṽ^-1U_R̅^*=W̃^-1Σ_R̅^-1U_R̅^*,

up to scaling and ordering.

Among closely related approaches are the State Space Method <cit.>
and ESPRIT <cit.> which construct the same Hankel
matrices and are equivalent to first order in the noise terms to the
Matrix Pencil method <cit.>. Other related approaches
are the SVD based Prony <cit.> and the Prony Koopman
Mode Decomposition <cit.> methods, which also yield
the same frequencies as DMD or the Matrix Pencil method for purely
sinusoidal data.

So far, our treatment of the X and Y data matrices was asymmetrical.
In our least squares approach we implicitly minimized the noise contribution
of the X matrix alone when computing L̃=ỸX̃^+.
However, Total Least Squares (TLS) variants of the above methods (TLS
Matrix Pencil <cit.>, TLS ESPRIT <cit.>,
TLS DMD <cit.>), have a symmetrical and
mutually similar treatment of the errors in both X and Y data
matrices.

Another asymmetry arises due to our treatment of the data series as
going forward in time, even though this direction should not be preferred
over backward time propagation at least when the eigenvalues are all
on the unit circle. An interested reader may refer to the forward-and-backward
(FB) Matrix Pencil Method <cit.> or the FB DMD <cit.>
for an unbiased treatment of the data, although it is worth noticing
that the two methods are not the same. FB Matrix Pencil extends the
observable to include data which propagates both forward and backwards
in time, which is appropriate for sequential time series. On the other
hand, FB DMD first computes two propagation matrices, in forward and
backward time, and then computes the geometric average of the former
with the inverse of the later. At any rate, the TLS and FB versions
of the above algorithms all require at least twice the computational
effort of standard variants.




 §.§ ESPRIT

ESPRIT was originally developed for estimating the direction of arrival of signals generated by a set of narrowband emitters and collected by an array of sensors <cit.>. As mentioned earlier, it is equivalent to the Matrix Pencil Method and  DMD with delayed observables when applied to the same sequential
data. However, an observation was made by the authors of ESPRIT, that when the emitted signals are uncorrelated their power can be estimated. This observation provided an inspiration for us to extend DMD to a clustering procedure for multiple times series
which will be described in the next section.

In the ESPRIT approach the data matrices are interpreted as


    X   =j=1R∑a_jb_j^T=AB^T
    
    Y   =j=1R∑λ_ja_jb_j^T=AΛ B^T,

where the columns of X are outputs of a first array of sensors, and
the columns of Y are outputs of a second array—identical, but
shifted relative to the first. Also, the jth column of A=[ a_1   ⋯ a_R ] corresponds to direction-of-arrival vectors of (or sensors response
to) the jth signal, the jth columns of B=[ b_1   ⋯ b_R ] consists of samples of the complex envelopes of the jth signal,
and Λ=diag{λ_j}∈ℂ^R× R
corresponds to a phase shift between the first array of sensors and
the second one <cit.>.

Although presented in a different setup, those definition of the data
matrices are equivalent to the DMD matrices when constructed from
delayed observables. Indeed, from (<ref>), (<ref>)
and (<ref>) we see that


    X=j=1R∑[   v_jλ_j^t_1            ⋯   v_jλ_j^t_M; v_jλ_j^t_1+1            ⋯ v_jλ_j^t_M+1;            ⋮                         ⋮; v_jλ_j^t_1+d            ⋯ v_jλ_j^t_M+d ], Y=j=1R∑[   v_jλ_j^t_1+1              ⋯   v_jλ_j^t_M+1;   v_jλ_j^t_1+2              ⋯   v_jλ_j^t_M+2;              ⋮                             ⋮; v_jλ_j^t_1+d+1              ⋯ v_jλ_j^t_M+d+1 ],

and hence may choose


    A=[ v_1/‖ v_1‖λ_1^0               ⋯ v_R/‖ v_R‖λ_R^0; v_1/‖ v_1‖λ_1^1               ⋯ v_R/‖ v_R‖λ_R^1;               ⋮                               ⋮; v_1/‖ v_1‖λ_1^d               ⋯ v_R/‖ v_R‖λ_R^d ], B=[ ‖ v_1‖λ_1^t_1             ⋯ ‖ v_1‖λ_1^t_M;             ⋮                           ⋮; ‖ v_R‖λ_R^t_1             ⋯ ‖ v_R‖λ_R^t_M ]

to put the matrices in the form of (<ref>).

In <cit.> the authors show that when the sources
are uncorrelated their “power” can be estimated. In other words
one can estimate B^*B if it is diagonal.  Unfortunately, while in the ESPRIT framework it is possible for the sources to be uncorrelated, in the formulation consistent with DMD (Eq. (<ref>)), B^*B can never be diagonal because each data snapshot is always affected by all of the frequencies.
However,
in the next section we will consider data consisting of multiple signals,
each exhibiting a subset of frequencies with different magnitudes,
i.e.,


    B=[ b_1,t_1λ_1^t_1              ⋯ b_1,t_Mλ_1^t_M;              ⋮                             ⋮; b_R,t_1λ_R^t_1              ⋯ b_R,t_Mλ_R^t_M ]],

where b_j,t_k∈ℝ^+. In that case DMD (or Matrix
Pencil method) might find the coefficients in B and thus characterize
each time snapshot according to the prevalent frequencies it exhibits.
Furthermore, it might be possible to group “similar” snapshots
together as will be discussed next.




§ CLUSTERING OF TIME SERIES


In this section we consider the prospects of DMD and Matrix Pencil
analysis to unsupervised classification of sequential data. Specifically,
we present a novel method for computing features of data series to
allow clustering them into groups of similar underlying dynamics,
that is, similar values of λ_j in (<ref>).

In section <ref> the underlying assumption
was that all of the data comes from a single system, or a group of
independent systems to which DMD or Matrix Pencil methods are to be
applied separately. If we suspect that the set of time series constituting
the data comes from a smaller set of possible dynamical systems, we
may wish to group the series based on similar underlying dynamics.

Among the many approaches to clustering of time series, we focus on
ones that assume a model behind the data <cit.>.
As an example, fitting autoregression parameters, defining some associated
metric and using them as features is one common approach <cit.>.
However, to the best of our knowledge, all contemporary clustering
techniques compute the model parameters (or features) for each time
series separately. Inspired by the ESPRIT algorithm (see discussion
in <ref>), we introduce a model-based approach which considers
all the times series at once, and uses elements of the DMD modes (or
generalized eigenvectors) as features for clustering purposes.




 §.§ DMD Clustering of Time Series


Consider N series each consisting of d+2 data points, { y_i(t)}⊂ℂ^n ,i=1,…,N, t=0,…,d+1.
We assume that each time series is a sum of damped sinusoids as in
section <ref> (Eq. (<ref>)).  Each time
series was generated by one of several different dynamical systems, and the
number of different systems (which we call P) is significantly smaller than the
number of time series N. Formally, we assume that the time series are split
into P partitions


    I_1∪⋯∪ I_P   ={ 1,⋯,N}
    
    I_1∩⋯∩ I_P   =∅,

each corresponding to the dynamics


    x_i(t)=k=1l_j∑v_ijkλ_jk^t∈ℂ^n, i∈ I_j

where l_j is the number of oscillators in the j-th partition,
{λ_jk} _k=1^l_j⊂ℂ are
the corresponding complex frequencies for that partition and x_iv_ijk,∈ℂ^n.

To state the general problem: we wish to estimate the complex frequencies
λ_jk as well as the partitions I_j from data given
by 

    y_i(t)=x_i(t)+s_i(t),

where x_i(t) is given by (<ref>) and s_i(t)
is a noise term.

First, consider the left generalized eigenvectors of the matrix pencil
(X,Y), where X and Y are data matrices


    X=[ y_1(0)      ⋯ y_N(0);      ⋮             ⋮; y_1(d)      ⋯ y_N(d) ], Y=[   y_1(1)        ⋯   y_N(1);        ⋮                 ⋮; y_1(d+1)        ⋯ y_N(d+1) ],

of the dimensions nd× N, which is the standard arrangement for DMD with delayed observables (similarly to (<ref>)).
The following property of the matrix
pencil in the noise-free case will be useful when determining the
partitions I_j. It requires that the number of series in each partition, |I_j| (here, |I_j| denotes the cardinality of the set I_j), is greater than nl_j, the dimensionality of the dynamic modes of the jth system.

Proposition 1: For j∈{1,…,P}, let Ξ_j be a matrix with
|I_j| columns given by

    Ξ_j=[   v_ij1;       ⋮; v_ijl_j ]_i∈ I_j∈ℂ^nl_j×|I_j|,
         i∈ I_j.

The matrix Ξ_j
thus has nl_j rows and consists of all the coefficients (modes) of all systems
in (<ref>) belonging to the jth partition and stacked
on top of each other. Assume that the matrices Ξ_j have full column
rank for all j∈{1,…,P} (and in particular |I_j|≥ nl_j=rank(Ξ_j), ∀ j),
that d≥∑ l_j and that no noise is present, s_i(t)=0, ∀ i,t.
Then the following hold:


  * The matrix pencil (X,Y) has nl linearly independent
left generalized eigenvectors where

    l=|{λ_jk : 1≤ j≤ P, 1≤ k≤ l_j}|

is the number of distinct eigenvalues among all systems in all partitions
(l≤∑ l_j).

  * Let p_λ be a generalized right eigenvector of the matrix
pencil (X,Y) corresponding to the eigenvalue λ.
If the ith data snapshot corresponding to the jth partition
(i∈ I_j) doesn't exhibit λ in its dynamics (λ∉{λ_jk} _k=1^l_j),
then the ith element of p_λ is zero. 

Proof: See Appendix A.

In (<ref>) one views the data as snapshots taken
at different times, similar to the approach in section <ref>
(Eq. (<ref>)). Alternatively, one can
arrange the data as


    X̂=[ y_1(0)      ⋯ y_1(d);      ⋮             ⋮; y_N(0)      ⋯ y_N(d) ], Ŷ=[   y_1(1)        ⋯ y_1(d+1);        ⋮                 ⋮;   y_N(1)        ⋯ y_N(d+1) ]

with X̂,Ŷ∈ℂ^nN× d, viewing all the
time series as consequent time snapshots of one large system. The
matrix pencil (X̂,Ŷ) has similar properties
to (X,Y) but requires “less” data and is therefore more useful for a numerical algorithm as stated in the next proposition and its discussion.

Proposition 2: For j∈{1,…,P}, let Ξ_j be a matrix with
n|I_j| rows, given by

    Ξ_j = [   v_ij1       ⋯ v_ijl_j ],
         i∈ I_j.∈ℂ^n|I_j|× l_j

Thus, the matrix Ξ_j has l_j columns and consists of all the coefficients (modes) of all systems
in (<ref>) belonging to the jth partition and placed
next to each other. Assume that Ξ_j has full row rank
for each j (i.e. n|I_j|≥ l_j=rank(Ξ_j)),
that d≥∑ l_j and that no noise is present, s_i(t)=0, ∀ i,t.
Then the following hold:


  * The matrix pencil (X̂,Ŷ) has l linearly
independent left generalized eigenvectors where l≤∑ l_j
is the number of unique eigenvalues (see (<ref>)). 

  * Let q_λ be a generalized left eigenvector of the matrix
pencil (X̂,Ŷ) corresponding to the eigenvalue
λ. If the ith data snapshot corresponding to the jth partition
(i∈ I_j) doesn't exhibit λ in its dynamics (λ∉{λ_jk} _k=1^l_j),
the n elements of q_λ beginning at i· n are all
zero. 

Proof: This follows directly from Proposition 1 for the matrix
pencil (X̂^T,Ŷ^T) when considering each
n dimensional system as n one dimensional systems.

We note that the assumptions in Proposition 1 require more snapshots
and yield more generalized eigenvectors compared to Proposition 2.
For this reason we chose the arrangement in (<ref>)
over (<ref>) in the example presented in section
<ref>.

The presence of zero elements in the generalized eigenvectors
corresponding to different data series can be exploited in order to find
the partitions in (<ref>). In other words, we propose
using the elements of the generalized eigenvectors (or adjoint DMD
modes) as features for clustering the data series.

Unfortunately, in presence of noise we do not expect any of the elements
of the generalized eigenvectors to be exactly zero as stated in Propositions
1 and 2. However, when using truncated SVD with the Matrix Pencil
method, the resulting generalized eigenvectors are continuous functions
of the data <cit.>. Consequently, we expect the elements
which should theoretically be zero, to remain close to zero for
low levels of noise; this will be verified in an example in
section <ref>.  We first discuss some details of
the implementation.




 §.§ Implementation


As mentioned earlier, the choice of arranging the data as in (<ref>)
has the benefits of giving a lower dimensional system and requiring
less data. We therefore will use the matrix pencil (X̂,Ŷ)
and its generalized right eigenvectors as features for clustering
the given time series.

Step 1: Compute the Singular Value Decomposition (SVD): 

    X̂=UΣ V^*

and choose a truncation value l̅ for the singular values.
This gives X̂≈ U_l̅Σ_l̅V_l̅^*
where U_l̅∈ℂ^nN×l̅, Σ_l̅∈ℂ^l̅×l̅
and V∈ℂ^d×l̅.

Step 2: Construct K̂=U_l̅^*ŶV_l̅Σ_l̅^-1∈ℂ^l̅×l̅
and find its left eigenvectors ŴK̂=Λ̂Ŵ,
where Λ̂,Ŵ∈ℂ^l̅×l̅
and Λ̂ is diagonal.

Step 3: Construct 

    Q=ŴU_l̅^*∈ℂ^l̅× nN,

whose rows are generalized right eigenvectors of the matrix pencil
(X̂,Ŷ).

Step 4: We define DMD features as sub-matrices q_i∈ℂ^l̅× n, 1≤ i≤ N
of Q, whose elements are defined as 

    (q_i)_j,k=|(Q)_k,in+j|, 1≤ j≤l̅, 1≤ k≤ n,

and |·| stands for the norm which is applied element-wise.
Furthermore there are many ways to define a metric on the above features,
for example


    d(q_i,q_j)=‖ q_i-q_j‖ _F,

where ‖·‖ _F is the Frobenius norm.

Step 5: Invoke any of the standard clustering methods on the
{ q_i} _i=1^N features which have a one-to-one
correspondence with the N time series.

We note that the matrix of right eigenvectors may also be expressed
as 

    Q=ŴΣ_l̅^-1V_l̅^*X̂^*.

Consequently, each of the DMD features can be seen as a transformation
of the appropriate data snapshot,


    q_i=|ŴΣ_l̅^-1V_l̅^*[ y_i^*(0);        ⋮; y_i^*(d) ]|,

where |·| stands for applying absolute value to
each elements of the matrix.

Moreover, as suggested previously in <cit.>, DMD
features can be used during an unsupervised learning stage and later
applied to never-seen-before data. In our case, the “training”
phase would consist of finding Ŵ, Σ_l̅, V_l̅
and the centers of clusters in the feature space, ℝ_+^l̅× n.
Then, new data { y(t)} _t=0^d can be mapped into
the feature space via 

    q=|ŴΣ_l̅^-1V_l̅^*[ y^*(0);      ⋮; y^*(d) ]|,

and assigned to a previously “learned” cluster.




§ NUMERICAL EXAMPLES





 §.§  A Toy Example - Multiple 1D Sinusoidal Signals


For illustration purposes we consider a set of one dimensional signals consisting of one or two sinusoids:


    x_1,...,6(t)   =α_1,...,6exp(iω_At)
    
      x_7,...,12(t)   =α_7,...,12exp(iω_Bt)
    
      x_13,...,23(t)   =α_13,...,23exp(iω_Ct)+β_13,...,23exp(iω_Dt),

where ω_A=1.0, ω_B=1.7, ω_C=0.8, ω_D=1.5 and the coefficients α_i,β_i∈ℂ are uniformly distributed on {.1≤|z|≤2|z∈ℂ}.
The measurements y(t)=Re{ x(t)} +s(t) were taken at t∈{ 0,1,...,19}
and the noise term s was normally distributed with σ=0.1. 

We now apply the technique presented in Section <ref> to recover the frequencies from the data generated by (<ref>) and cluster the signals based on those frequencies.
We arrange the data as in (<ref>) into X̂,Ŷ∈ℝ^23×19,
where N=23 is the number of signals.



The singular values of the X̂ are shown in figure <ref>(a) and exhibit a sharp drop in magnitude after the first eight. This implies that the data contains four superimposed sinusoids. To estimate the frequencies we compute the DMD as suggested in section <ref> with l̅=8 modes. 

Figure <ref>(b) shows the estimated eigenvalues of the matrix pencil (X̂,Ŷ) which match the frequencies corresponding to the discrete measurements – exp(± iω_A,B,C,D). However, the singular values and the eigenvalues alone are not sufficient to determine which signals correspond to which frequencies. To do that, we compute the generalized right eigenvectors per (<ref>)

Figure <ref>(c) and (d) show magnitudes  of four eigenvectors. Two of them (fig. <ref>(c)) correspond to frequencies close to ω_C and ω_D and have elements with magnitudes close to zero at indices 1,...,12. With proposition 2 in mind, this suggests that the other sequences (13,...,23) must exhibit those frequencies, which is indeed the case for the data generated by (<ref>). Similarly, the other two eigenvectors (fig. <ref>(d)), corresponding to frequencies close to ω_A and ω_B, have elements close to zero at all indices except 1,...,6 and 7,...,12 respectively. Again, we have correctly identified the signals those frequencies belong to. 

In this example one could easily conclude that a single frequency close to ω_A is present in sequences 1,...6, ω_B in 7,...12 and exactly two frequencies, ω_C,ω_D, are present in sequences 13,...23. For a large multidimensional dataset, this identification task requires a clustering algorithm as will be illustrated in the next example.




 §.§ Clustering of Lattice Regions in TEM
Image


For a more practical example, we consider a Transmission Electron Microscopy
(TEM) image of a lattice of gold atoms from <cit.>, shown
in figure <ref>(a). We consider the variation in brightness
of the image as one transverses its pixels from left to right (x axis)
or top to bottom (y axis). Figure <ref>(b) shows those
variations over 51 pixels centered about a single pixel in a particular
region of the lattice. In that region, the hexagonal close-packed
(HCP) lattice is oriented (almost) parallel to the x axis of the image.
The brightness variation in the x direction has a period of about
19 pixels - the closest neighbor distance in the lattice. The variation
in the y axis, in this region, has a period of exactly 1.5 times
the period in the x axis (about 28 pixels). One period in the y axis
however includes two rows of the lattice, hence the variations
in the y direction in that region exhibit two frequencies,
ω_y=1.5ω_x≈ 1/28
and 2ω_y.




In various regions in figure <ref>(a), the grain has
different orientations and therefore different frequencies of variations
in brightness in x and y directions. In this example will use DMD features and constrained hierarchical clustering <cit.>
to split the pixels of the image based on these variations.

We construct the data matrices by combining the brightness data
(figure <ref>(b)) for each pixel in each direction x
and y (n=2).  Letting P_ij denote the brightness of pixel (i,j), we
arrange the data as


    { y_ij(t)} _t=0^d={[ P_i-d/2,j; P_i,j-d/2 ],...,
        [ P_i-1,j; P_i,j-1 ],
        [ P_i,j; P_i,j ],
        [ P_i+1,j; P_i,j+1 ],...,
        [ P_i+d/2,j; P_i,j+d/2 ]}

The number of data series (N) is the number of pixels for which
brightness variations over 51 pixels (d=50) were collected. The
partitions in (<ref>) therefore correspond to regions
where grain orientation remains the same or regions where no lattice is present.

The generalized left eigenvectors of (X̂,Ŷ)
in (<ref>) have 2N elements (per <ref>)
which we rearrange back into the shape of the original image, once
for the x and once for the y direction. Figure <ref>
shows the element-wise absolute value of a generalized left eigenvector
corresponding to the frequency closest to 2ω_y≈ 1/14.
As discussed earlier, this frequency is present only in brightness
variations in the y direction in the triangular shaped region where
the lattice is parallel to the x axis (marked in figure <ref>(a)).
Indeed figure <ref>(a) shows the part of
the eigenvector corresponding to the y direction, which exhibits large
magnitudes in the triangular area where this frequency occurs. On
the other hand, the elements corresponding to the 2ω_y frequency
but the x direction (figure <ref>(b)) are close
to zero in the same area as expected per Proposition 2 in
<ref>.




Having computed the feature vectors (Eq. (<ref>)),
we may proceed with a clustering scheme of our choice. Since in this
particular case the feature space alone doesn't contain any spatial
information (i.e., pixels relative position), we speculate
that a contiguity-constrained clustering <cit.>
approach would be appropriate. Specifically, a constrained version
<cit.> of the Ward hierarchical agglomerative
clustering algorithm <cit.> seems to give satisfactory
results. In this approach, each pixel begins as a singleton cluster.
Then clusters are iteratively merged in a greedy manner, such that
each two newly merged clusters minimize the Ward criterion <cit.>
with respect to the DMD features. However, clusters can be merged only
if they have adjacent pixels, i.e., if they satisfy the spatial connectivity
constraints. This reduces computational time while ensuring that regions
of the image corresponding to similar dynamics remain connected.
An implementation of the Ward agglomerative clustering scheme is available
with Scikit-learn <cit.>, a Python machine
learning library.

The results are shown on figure <ref>, where the pixels
are split into 6 regions based on their DMD features computed as specified
in section <ref>. Since hierarchical clustering produces
a dendogram which has separate pixels as leaves and a single cluster
containing the whole image as a root, any number of clusters between
those two extremes is a priori valid. However, once the dendogram is
complete, it requires little additional computational effort to retrieve
the partitions for any possible number of clusters. In the example
presented here, we find that a choice of 6 partitions seem to capture
the different regions in figure <ref> (a) correctly.
Indeed in figure <ref> we see four partitions where the
lattice orientation causes different frequencies in variation in brightness
in x and y direction. Two more partitions have no lattice present
and form two spatially separated regions.







§ CONCLUSIONS


We have described connections between Dynamic Mode Decomposition (DMD) and a family of frequency detection
methods (namely, the Matrix Pencil, State Space, and ESPRIT methods) that were
formerly known to be equivalent to one another <cit.>. Furthermore, it
was shown that, for sequential data, DMD with delayed observables
yields exactly the same eigenvalues as the Matrix Pencil method, while
the generalized eigenvectors were found to be the adjoint DMD modes.
Based on these similarities, and inspired by the formulation of the
ESPRIT method, a novel approach for clustering of time series was
proposed.

We have considered data in the form of a large number of time series
with a smaller number of possible underlying dynamics (frequencies)
for each series. It was shown that, when arranging all the data together,
the generalized eigenvectors (adjoint DMD modes) of the corresponding
Hankel matrices have an interesting property: their elements corresponding
to certain series and frequencies tend to have magnitudes close to
zero when those series do not exhibit those particular frequencies.
We therefore conclude that projecting a sequence onto DMD modes computed
from a larger set of given sequences, is a viable feature extraction
method given a large amount of unlabeled series.

While existing approaches to clustering of time series require extracting
features from each sequence first, DMD and the Matrix Pencil methods
estimate the number of required features and extract them across all
series at once. Various metrics can then be defined on the elements
of the adjoint DMD modes and a range of clustering algorithms can
be applied on the time series accordingly. This method has been illustrated
by clustering regions of distinct patterns in an image based on the
variance in brightness of spatially close pixels.




§ APPENDIX A


Proof of Proposition 1:

Without loss of generality, we assume throughout the proof that the
columns of X and Y defined in (<ref>) are
grouped together based on the partitions, i.e. the first |I_1|
columns are data from systems in the first partition etc. This makes
the proof more tractable but doesn't affect our conclusion regarding
the elements of the generalized eigenvectors.

First, consider the case when n=1, |I_j|=l_j
and Ξ_j=I_l_j× l_j is the identity matrix. The data
matrices take the following form:


    X^#=[        1        1        ⋯        1        1        ⋯        1;     λ_11     λ_12        ⋯   λ_1l_1     λ_21        ⋯   λ_Pl_P;        ⋮        ⋮                 ⋮        ⋮                 ⋮;   λ_11^d   λ_12^d        ⋯ λ_1l_1^d   λ_21^d        ⋯ λ_Pl_P^d ]∈ℂ^d×∑ l_j,

where the columns are increasing powers of all the complex frequencies
of all systems arranged in order of their partitions.

From the statement of the theorem d≥ l, where l≤∑ l_j
is the number of unique exponents among all partitions, thus the rank
of X^# is l. Let e_jk be vectors of the standard basis
of ℂ^∑ l_j such that the non-zero element of e_jk
corresponds to the column of X^# where the powers of λ_jk
appear, and consider the matrix pencil


    λ X^#-Y^#=[             λ-λ_11                  ⋯           λ-λ_1l_1             λ-λ_21                  ⋯;                  ⋮                                     ⋮                  ⋮;     λ_11^d(λ-λ_11)                  ⋯ λ_1l_1^d(λ-λ_1l_1)     λ_21^d(λ-λ_21)                  ⋯ ],

where Y^# is defined similarly to X^# with all powers increased
by one.

Clearly, if λ is one of the exponents λ_jk, it
is a generalized eigenvalue of the matrix pencil (X^#,Y^#),
since in that case rank(Y^#-λ X^#)=rank(X^#)-1.
The corresponding generalized right eigenvector is


    p_λ^#=λ_jk=λ∑e_jk

since it can be easily shown that p_λ^# in the row space
of X^#. We have therefore found all the l generalized eigenvectors
of (X^#,Y^#). Moreover, generalized eigenvectors
that belong to distinct eigenvalues are orthogonal in this case (whether
they belong to the same partition or not).

Remaining in the one dimensional case (n=1), we now relax the assumptions
on the number of series and modes per partition, allowing |I_j|≥ l_j
and Ξ_j of any form as long as satisfies rank(Ξ_j)=l_j.
The data matrices may now be expressed as


    X=[ [        1        ⋯        1;        ⋮                 ⋮;   λ_11^d        ⋯ λ_1l_1^d ]                                                                                      ⋯ [        1        ⋯        1;        ⋮                 ⋮;   λ_P1^d        ⋯ λ_Pl_P^d ] ][ Ξ_1       0;       ⋱;   0     Ξ_P ],

and the matrix pencil as


    Y-λ X=(Y^#-λ X^#)[ Ξ_1       0;       ⋱;   0     Ξ_P ]



Since Ξ_j∈ℂ^l_j×|I_j|, ∀ j
have full column rank,


    p_λ=[ Ξ_1       0;       ⋱;   0     Ξ_P ]^+p_λ^#=λ_jk=λ∑[ Ξ_1^+           0;           ⋱;     0       Ξ_P^+ ]e_jk

is clearly a generalized eigenvector of (X,Y) when λ
is one of the l unique exponents. Since X has rank l, we
have found all the generalized right eigenvectors. Moreover the structure
of e_jk and the matrix on the right hand side of (<ref>),
reveals that p_λ may only have non-zero elements at indices
belonging to partitions which exhibit λ in their dynamics.
Thus, the proposition is proved for this case.

Finally, we consider the case n>1 in which Ξ_j∈ℂ^nl_j×|I_j|
has rank nl_j. In this case, the data matrices can be expressed
as


    X=(X^#⊗ I_n× n)[ Ξ_1       0;       ⋱;   0     Ξ_P ], Y=(Y^#⊗ I_n× n)[ Ξ_1       0;       ⋱;   0     Ξ_P ]

where ⊗ denotes the Kronecker product and I_n× n
and an n by n identity matrix.

The argument in the proof of the n=1 case applies to each dimension
separately, resulting in n times as many generalized eigenvalues
and eigenvectors, and the conclusion still holds.

 abbrv


