remarkRemark hypothesisHypothesis claimClaim Discriminant Dynamic Mode
DecompositionN. Takeishi, K. Fujii, K. Takeuchi, and Y. Kawahara

# Discriminant Dynamic Mode Decomposition
for Labeled Spatio-Temporal Data Collections††thanks: Preprint. The major part
of this work was done when the first author was working at RIKEN Center for
Advanced Intelligence Project. This work was supported by JSPS KAKENHI Grant
Numbers JP19K21550, JP20H04075, JP19H04941, and JP18H03287; JST PRESTO Grant
Number JPMJPR20C5; JST CREST Grant Number JPMJCR1913; and AMED Grant Number
JP19dm0307009.

Naoya Takeishi RIKEN Center for Advanced Intelligence Project, Tokyo, Japan ()
and University of Applied Sciences and Arts Western Switzerland, Geneva,
Switzerland<EMAIL_ADDRESS>Keisuke Fujii Graduate School of
Informatics, Nagoya University, Aichi, Japan () and RIKEN Center for Advanced
Intelligence Project, Tokyo, Japan<EMAIL_ADDRESS>Koh Takeuchi
Graduate School of Informatics, Kyoto University, Kyoto, Japan () and RIKEN
Center for Advanced Intelligence Project, Tokyo, Japan.
<EMAIL_ADDRESS>Yoshinobu Kawahara Institute of Mathematics for
Industry, Kyushu University, Fukuoka, Japan () and RIKEN Center for Advanced
Intelligence Project, Tokyo, Japan<EMAIL_ADDRESS>
###### Abstract

Extracting coherent patterns is one of the standard approaches towards
understanding spatio-temporal data. Dynamic mode decomposition (DMD) is a
powerful tool for extracting coherent patterns, but the original DMD and most
of its variants do not consider label information, which is often available as
side information of spatio-temporal data. In this work, we propose a new
method for extracting distinctive coherent patterns from labeled spatio-
temporal data collections, such that they contribute to major differences in a
labeled set of dynamics. We achieve such pattern extraction by incorporating
discriminant analysis into DMD. To this end, we define a kernel function on
subspaces spanned by sets of dynamic modes and develop an objective to take
both reconstruction goodness as DMD and class-separation goodness as
discriminant analysis into account. We illustrate our method using a synthetic
dataset and several real-world datasets. The proposed method can be a useful
tool for exploratory data analysis for understanding spatio-temporal data.

###### keywords:

time-series, dynamic mode decomposition, discriminant analysis

68T10, 37M10, 37N99

##  1 Introduction

Spatio-temporal data are ubiquitous in modern science and engineering, but
they are often so complicated and high-dimensional that we can hardly read
useful information from them directly. One of popular approaches to
understanding spatio-temporal data is to extract from them _coherent patterns_
across space and time with which we can summarize the characteristics of
spatio-temporal data in some sense. Such summarization can be useful not only
for understanding underlying dynamics but also for crafting machine learning
features and models, reduced order modeling, and designing controllers.

Extraction of spatio-temporal coherent patterns can be performed in different
ways. Many popular approaches are formulated as dimensionality reduction
techniques including principal component analysis (PCA) and other
matrix/tensor factorization techniques [see, e.g., 23; 26; 11]. The use of PCA
for coherent patterns extraction is indeed prevalent in a wide range of
applications, and there have been countless numbers of use cases. For example,
we may characterize complex neural activities using a small number of spatio-
temporal patterns computed by PCA [10]. Moreover, a method called proper
orthogonal decomposition (POD), which is equivalent to PCA in principle, has
been a standard technique for reduced order modeling in fluid mechanics.
Dynamic mode decomposition (DMD) is also known as a powerful tool to extract
spatio-temporal coherent patterns [see, e.g., 27, and references therein]. DMD
was proposed originally in the field of fluid mechanics [33; 34] and has been
successfully utilized in a broad range of applications such as neuroscience
[6], epidemiology [31], power system analysis [38], and building maintenance
design [18]. DMD is favored mainly because it considers the dynamical
properties of spatio-temporal data, which are usually ignored by PCA and POD.

In this work, we tackle the problem of extracting coherent patterns from
_labeled_ spatio-temporal data collections. In many practices of spatio-
temporal data analysis, label information is available as side information.
For example, biological recordings, such as neural signals, are often
accompanied with additional information of targets (e.g., types of applied
stimuli). Data from social phenomena, such as recordings of traffic and
people, can often be labeled with information of calendar (e.g., days of week
and holidays). Such label information can be a useful clue for extracting
meaningful coherent patterns from spatio-temporal data.

Data summarization considering label information has been addressed also in
several other contexts. For example, Fisher’s linear discriminant [13] is a
classical and well-known method for dimensionality reduction of labeled data.
Supervised PCA [4] is a variant of PCA that takes labels of each data-point
into consideration. In brain-computer interface studies, extraction of common
spatial patterns [32] from a set of labeled neural signals is a popular tool
for classification. However, most existing methods for labeled data
summarization ignore dynamical aspects of data, and thus these methods are not
necessarily suitable for spatio-temporal data.

Figure 1: Construction of the paper. In Section 2, we introduce three
technical building blocks: DMD, kernel discriminant analysis, and Grassmann
kernels. We then define a kernel function on sets of dynamic modes in Section
3.1. Finally, we formulate our proposal of discriminant DMD as a combination
of the DMD-based kernel and kernel discriminant analysis in Section 3.2. We
also show its numerical example in Section 3.3.

In this paper, we present a method to extract coherent patterns from labeled
spatio-temporal data such that they reflect _both_ of the labels and the
dynamics. We build such an extraction method upon DMD. Since the original DMD
algorithm and the most existing variants cannot consider label information, we
propose to incorporate discriminant analysis into DMD’s formulation.
Specifically, we combine an optimization formulation of DMD [9] with a
parameter optimization method for kernel Fisher discriminant analysis [42]. To
this end, we need a kernel function to measure similarity between sets of
dynamic modes extracted from different time-series episodes. We design such a
kernel function for sets of dynamics modes based on the notion of Grassmann
kernels on subspaces [22]. Finally, we formulate an optimization problem to
compute _discriminant DMD_ , with which we can obtain distinctive coherent
patterns that contribute to major differences in spatio-temporal behaviors of
different labels.

The remainder of the paper is structured as follows. In Section 2, we review
three techniques that constitute building blocks. In Section 3, we then
present the discriminant DMD with an illustrative numerical example. Fig. 1
summarizes the construction of Section 2 and Section 3. Afterward, we discuss
the related work in Section 4. We showcase some applied use cases of the
proposed method with real-world datasets in Section 5. Finally, we give
discussion in Section 6 and conclude the paper in Section 7.

##  2 Background

We review three lines of studies that constitute the technical building blocks
of the proposed method, which combines DMD and discriminant analysis. We first
review these two notions in Section 2.1 and Section 2.2, respectively.
Afterward, we review the concept of Grassmann kernels in Section 2.3, which
plays a key role in combining the first two notions.

###  2.1 Modal Decomposition of Dynamical Systems

The operator-theoretic perspective for analyzing dynamical systems is the
first building block toward the purpose of this work. We first introduce some
basic notions of dynamical systems and their operator-theoretic view [29] and
then review the modal decomposition of dynamical systems by DMD [33; 34].

####  2.1.1 Dynamical Systems and Koopman Operators

Variety of physical phenomena can be modeled as _dynamical systems_. In this
paper, we primarily consider discrete-time dynamical systems like

(1) $v_{t+1}=f(v_{t}),$

where $v_{t}\in\mathcal{M}$ is a state vector, $\mathcal{M}$ is a state space,
and $t\in\\{0\\}\cup\mathbb{N}$ is a time index. When map $f$ is nonlinear,
which is often the case, it is highly challenging to analyze the dynamical
system’s behavior directly in state space $\mathcal{M}$.

The operator-theoretic perspectives, such as ones based on the _Koopman
operator_ [see, e.g., 29; 7], are attracting attention as an alternative to
analyze nonlinear dynamical systems. Consider a function on the state space,
namely $g\colon\mathcal{M}\to\mathbb{R}$, which is referred to as an
_observable_. We suppose $g$ is in some Hilbert space $\mathcal{G}$. For the
dynamical system in Eq. 1, the corresponding Koopman operator $\mathcal{K}$ is
defined as

(2) $\mathcal{K}g(v)=g(f(v)).$

Because $\mathcal{G}$ is a vector space, $\mathcal{K}$ is a _linear_ operator
and infinite-dimensional in general.

An advantage of considering the Koopman operator appears as its spectral
decomposition. For simplicity of discussion, let us assume that the spectral
decomposition of $\mathcal{K}$ is well-defined and that $\mathcal{K}$ has only
simple point spectra. Let $\lambda_{j}\in\mathbb{C}$ and
$\varphi_{j}\colon\mathcal{M}\to\mathbb{C}$ be an eigenvalue and an
eigenfunction of $\mathcal{K}$, respectively, i.e.,

(3) $\mathcal{K}\varphi_{j}(v)=\lambda_{j}\varphi_{j}(v),\quad
j\in\mathbb{N}.$

If $g$ is in the span of the (possibly countably many) eigenfunctions of
$\mathcal{K}$, we have a decomposition of the observable:

(4)
$\mathcal{K}g(v)=\sum_{j=1}^{\infty}\lambda_{j}\varphi_{j}(v)\tilde{w}_{j}(g),$

where $\tilde{w}_{j}(g)\in\mathbb{C}$ is a coefficient of the orthogonal
projection of $g$ onto the span of $\varphi_{j}$. By applying $\mathcal{K}$ to
the both sides of Eq. 4 repeatedly (see Eq. 2 and Eq. 3), we obtain

(5)
$g(v_{t})=\sum_{j=1}^{\infty}\lambda_{j}^{t}\big{(}\varphi_{j}(v_{0})\tilde{w}_{j}(g)\big{)},$

where $v_{0}\in\mathcal{M}$ is a initial condition. The expression in Eq. 4 or
Eq. 5 is called the Koopman mode decomposition (KMD) [7] in literature and
have been utilized for discovering meaningful spatio-temporal patterns of
dynamics [see, e.g., 7; 27, and references therein].

####  2.1.2 Dynamic Mode Decomposition

Dynamic mode decomposition (DMD) [33; 34] is a method to compute modal
decomposition of dynamics. DMD and its variants have been utilized in a wide
range of domains [see, e.g., 27, and references therein] partly because of its
mathematical and algorithmic simplicity. It also has a connection to KMD under
some conditions [7; 40; 2; 25]. We briefly review the idea and techniques of
DMD.

##### Eigendecomposition-Based Formulation

We present an abstracted description of a working algorithm of DMD based on
pseudoinverse and eigendecomposition [39]. Suppose we have a time-series
episode $X=(\bm{x}_{1},\dots,\bm{x}_{\tau})$ obtained from some dynamics,
where $\bm{x}_{t}\in\mathbb{R}^{p}$ denotes an observation (also called a
snapshot) at time $t$. We compute the eigendecomposition of a matrix $\bm{A}$
such that

(6) $\bm{A}=\bm{X}^{+}(\bm{X}^{-})^{\dagger},$

where $\bm{X}^{+}$ and $\bm{X}^{-}$ are time-lagged data matrices, i.e.,

(7) $\displaystyle\bm{X}^{-}$
$\displaystyle=\begin{bmatrix}\bm{x}_{1}&\cdots&\bm{x}_{\tau-1}\end{bmatrix}\in\mathbb{R}^{p\times(\tau-1)},$
$\displaystyle\bm{X}^{+}$
$\displaystyle=\begin{bmatrix}\bm{x}_{2}&\cdots&\bm{x}_{\tau}\end{bmatrix}\in\mathbb{R}^{p\times(\tau-1)},$

and $\cdot^{\dagger}$ denotes the Moore–Penrose pseudo inverse. When $p\gg 1$,
we often perform the dimensionality reduction of data to $r$ dimensions
($r<p$) using the singular value decomposition. We then compute the
eigendecomposition of the dimension-reduced version of $\bm{A}$, and finally
the eigenvectors are projected back to the original $p$-dimensional space. See
[39] for the detailed procedures.

DMD’s modal decomposition of a time-series episode is formulated using the
spectra of the $\bm{A}$ matrix. Let $\hat{\bm{w}}_{j}\in\mathbb{C}^{p}$ and
$\bm{z}_{j}\in\mathbb{C}^{p}$ be a pair of the right- and left-eigenvectors of
$\bm{A}$, respectively, corresponding to $\bm{A}$’s non-zero eigenvalue
$\lambda_{j}\in\mathbb{C}$, for $j=1,\dots,p^{\prime}$, where $p^{\prime}$ is
the number of the non-zero eigenvalues. Without loss of generality, we suppose
that they are normalized so that
$\hat{\bm{w}}_{j}^{\mathsf{H}}\bm{z}_{j^{\prime}}=\delta_{j,j^{\prime}}$,
where $\delta_{j,j^{\prime}}=1$ if $j=j^{\prime}$ and
$\delta_{j,j^{\prime}}=0$ otherwise. If $\bm{X}^{-}$ and $\bm{X}^{+}$ are
linearly consistent [39] and $\lambda_{1},\dots,\lambda_{p^{\prime}}$ are
distinct, we have

(8)
$\bm{x}_{t}=\sum_{j=1}^{p^{\prime}}\lambda_{j}^{t-1}(\bm{z}_{j}^{\mathsf{H}}\bm{x}_{1})\hat{\bm{w}}_{j}.$

In Eq. 8, each snapshot is represented as a weighted sum of the vectors,
$\hat{\bm{w}}_{1},\dots,\hat{\bm{w}}_{p^{\prime}}$, which are termed _dynamic
modes_. Note that the weights of each dynamic mode are given by the
corresponding (complex) eigenvalues, $\lambda_{1},\dots\lambda_{p^{\prime}}$.
Hence, the $j$-th dynamic mode, $\bm{w}_{j}$, can be regarded as a spatial
coherent pattern that exhibits oscillation with an angular frequency
$\angle\lambda_{j}$ and a decay/growth rate $|\lambda_{j}|$.

There is a formal resemblance between KMD Eq. 5 and DMD Eq. 8. In fact, we can
find further connections between the two decomposition methods under some
conditions. Let us see one of them. Suppose that a snapshot $\bm{x}$ is
generated by a concatenation of the values of multiple observables
$g_{1},\dots,g_{p}\in\mathcal{G}$, i.e.,

(9)
$\bm{x}_{t}=\bm{g}(v_{t})=\begin{bmatrix}g_{1}(v_{t})&\cdots&g_{p}(v_{t})\end{bmatrix}^{\mathsf{T}}.$

With this interpretation, Tu _et al._ [39] have shown that if $\varphi_{j}$ is
in the span of $\\{g_{1},\dots,g_{p}\\}$, and the data is sufficiently rich,
then pointwise evaluations of the eigenfunction are given by
$\varphi_{j}(v_{t})=\bm{z}_{j}^{\mathsf{H}}\bm{x}_{t}$. Moreover, the
connection between variants of DMD and KMD has been discussed [40; 2; 25]. In
this sense, DMD may be roughly understood as a way to approximate spectral
components of the Koopman operator, though care must be taken when one needs
such an interpretation. We note that DMD is still a useful tool for
dimensionality reduction and pattern extraction even if the connection to the
Koopman operator spectrum is hardly established in practice.

##### Opimitzation-Based Formulation

Apart from the eigendecomposition-based formulation, DMD can also be defined
as an exponential curve fitting to multivariate time-series. This paradigm is
often termed optimized DMD [9] as it is based on nonlinear optimization. Let

$\bm{W}=\begin{bmatrix}\bm{w}_{1}&\cdots&\bm{w}_{r}\end{bmatrix}\in\mathbb{C}^{p\times
r}$

be a matrix comprising $r$ ($\leq p$) dynamic modes (we can think that
$\bm{w}_{j}$ corresponds to
$(\bm{z}_{j}^{\mathsf{H}}\bm{x}_{1})\hat{\bm{w}}_{j}$ in Eq. 8 in its role).
Let $\bm{V}_{\lambda_{1:r}}\in\mathbb{C}^{r\times\tau}$ denote a Vandermonde
matrix parameterized by a set
$\lambda_{1:r}=\\{\lambda_{1},\dots,\lambda_{j}\\}$ as

(10)
$\bm{V}_{\lambda_{1:r}}=\begin{bmatrix}1&\lambda_{1}&\lambda_{1}^{2}&\cdots&\lambda_{1}^{\tau-1}\\\
1&\lambda_{2}&\lambda_{2}^{2}&\cdots&\lambda_{2}^{\tau-1}\\\
\vdots&\vdots&\vdots&&\vdots\\\
1&\lambda_{r}&\lambda_{r}^{2}&\cdots&\lambda_{r}^{\tau-1}\end{bmatrix}\in\mathbb{C}^{r\times\tau}.$

Then, optimized DMD is formulated as an optimization problem

(11)
$\underset{\bm{W},\lambda_{1:r}}{\text{minimize}}~{}~{}\|\bm{X}-\bm{W}\bm{V}_{\lambda_{1:r}}\|_{F}^{2}.$

The objective in Eq. 11 measures the squared error of the reconstruction
$\bm{W}\bm{V}_{\lambda_{1:r}}$ with the dynamic modes. Conceptually it
corresponds to the error between the left- and right-hand sides of Eq. 8 of
the eigendecomposition-based formulation.

The optimization problem in Eq. 11 can be addressed with the strategy called
_variable projection_ [21; 3] because it is a linear least squares problem
with regard to $\bm{W}$. That is, if we fix $\lambda_{1:r}$ at some value
$\lambda_{1:r}^{*}$, the corresponding optimal $\bm{W}^{*}$ is immediately
determined by

(12) $\bm{W}^{*}=\bm{X}\big{(}\bm{V}_{\lambda_{1:r}^{*}}\big{)}^{\dagger}.$

Hence, Equation 11 reduces to another nonlinear problem:

(13)
$\underset{\lambda_{1:r}}{\text{minimize}}~{}~{}\left\|\bm{X}-\bm{X}\bm{V}_{\lambda_{1:r}}^{\dagger}\bm{V}_{\lambda_{1:r}}\right\|_{F}^{2}.$

See [3] for details on this formulation. Later in Section 3, we will show that
formulating DMD as an optimization problem is useful in developing our
proposal to utilize label information.

###  2.2 Discriminant Analysis

Discriminant analysis forms the second building block of our method, and here
we review a basic form and extensions. First, we briefly introduce the
Fisher’s linear discriminant [13] and its nonlinear counterpart using the
kernel trick [30]. Then, we review a method [42] to optimize kernel parameters
of the kernelized discriminant analysis.

####  2.2.1 Fisher’s linear discriminant

Fisher’s linear discriminant [13] is a classical supervised dimensionality
reduction method and finds a linear transform of features so that different
classes are well separated. It is defined as maximization of a ratio of
between- and within-class variances. Suppose we have a set of data points111We
slightly violate the notations in Section 2.1; here, $\bm{x}$ is just a data
point in general and does not necessarily constitute a sequence as a snapshot.
with binary labels, $\\{(\bm{x}_{i},y_{i})\mid i=1,\dots,n\\}$, where
$\bm{x}_{i}\in\mathbb{R}^{p}$ and $y_{i}\in\\{1,2\\}$. Then, Fisher’s linear
discriminant can be formulated as an optimization problem

(14)
$\operatorname*{maximize}_{\bm{b}\in\mathbb{R}^{p}}~{}~{}\frac{\bm{b}^{\mathsf{T}}\bm{S}_{\mathrm{B}}\bm{b}}{\bm{b}^{\mathsf{T}}\bm{S}_{\mathrm{W}}\bm{b}}.$

In Eq. 14, $\bm{S}_{\mathrm{B}}$ and $\bm{S}_{\mathrm{W}}$ denote between- and
within-class scatter matrices defined as

(15) $\displaystyle\bm{S}_{\mathrm{B}}$
$\displaystyle=\frac{1}{n}\sum_{l=1,2}n_{l}(\bm{\mu}_{l}-\bm{\mu}_{\text{all}})(\bm{\mu}_{l}-\bm{\mu}_{\text{all}})^{\mathsf{T}},$
$\displaystyle\bm{S}_{\mathrm{W}}$
$\displaystyle=\frac{1}{n}\sum_{l=1,2}\sum_{i\mid
y_{i}=l}(\bm{x}_{i}-\bm{\mu}_{l})(\bm{x}_{i}-\bm{\mu}_{l})^{\mathsf{T}},$

where $n_{l}$ and $\bm{\mu}_{l}$ denote the number and the average of data
points with label $y=l$ (for $l=1,2$), respectively, i.e.,

(16)
$n_{l}=|\\{i|y_{i}=l,\,i=1,\dots,n\\}|\quad\text{and}\quad\bm{\mu}_{l}=\frac{1}{n_{l}}\sum_{i\mid
y_{i}=l}\bm{x}_{i}.$

$\bm{\mu}_{\text{all}}$ denotes the average of all the data points,
$\bm{\mu}_{\text{all}}=\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}$. Equation 14 can
be solved as a generalized eigenvalue problem.

####  2.2.2 Kernel Fisher Discriminant Analysis

Mika _et al._ [30] formulated a nonlinear extension of Fisher’s linear
discriminant using the kernel trick, namely _kernel Fisher discriminant_
(KFD). Let $k\colon\mathbb{R}^{p}\times\mathbb{R}^{p}\to\mathbb{R}$ be a
kernel function and $\phi\colon\mathbb{R}^{p}\to\mathcal{H}$ be the
corresponding feature map, where $\mathcal{H}$ denotes the reproducing kernel
Hilbert space (RKHS) corresponding to $k$. Then, in a two-class case, KFD is
formulated as the problem:

(17)
$\operatorname*{maximize}_{\bm{b}^{\phi}\in\mathcal{H}}~{}~{}\frac{(\bm{b}^{\phi})^{\mathsf{T}}\bm{S}^{\phi}_{\mathrm{B}}\bm{b}^{\phi}}{(\bm{b}^{\phi})^{\mathsf{T}}\bm{S}^{\phi}_{\mathrm{W}}\bm{b}^{\phi}},$

where $\bm{S}_{\mathrm{B}}^{\phi}$ and $\bm{S}_{\mathrm{W}}^{\phi}$ denote
between- and within-class scatter matrices in feature space $\mathcal{H}$,
respectively. These scatter matrices are defined analogously to Eq. 15 as
follows:

(18) $\displaystyle\bm{S}^{\phi}_{\mathrm{B}}$
$\displaystyle=\frac{1}{n}\sum_{l=1,2}n_{l}(\bm{\mu}^{\phi}_{l}-\bm{\mu}^{\phi}_{\text{all}})(\bm{\mu}^{\phi}_{l}-\bm{\mu}^{\phi}_{\text{all}})^{\mathsf{T}},$
$\displaystyle\bm{S}^{\phi}_{\mathrm{W}}$
$\displaystyle=\frac{1}{n}\sum_{l=1,2}\sum_{i\mid
y_{i}=l}(\phi(\bm{x}_{i})-\bm{\mu}^{\phi}_{l})(\phi(\bm{x}_{i})-\bm{\mu}^{\phi}_{l})^{\mathsf{T}},$

where $\bm{\mu}^{\phi}_{l}$ denote the sample mean of the data points of the
$l$-th class ($l\in[1,c]$) in $\mathcal{H}$, that is,

(19) $\displaystyle\bm{\mu}^{\phi}_{l}$
$\displaystyle=\frac{1}{n_{l}}\sum_{i\mid y_{i}=l}\phi(\bm{x}_{i}).$

$\bm{\mu}^{\phi}_{\text{all}}$ denotes the average of all the data points in
$\mathcal{H}$,
$\bm{\mu}^{\phi}_{\text{all}}=\frac{1}{n}\sum_{i=1}^{n}\phi(\bm{x}_{i})$.

We do not describe the solution of KFD because we do not need to solve it
explicitly. Instead, we use it for providing a criterion of the goodness of
discriminant through a kernel parameter optimization scheme for KFD, which is
reviewed below.

####  2.2.3 Kernel Parameter Optimization for KFD

It is often difficult to predetermine the parameters of a kernel used in KFD.
You _et al._ [42] proposed a method for optimizing kernel parameters of KFD.
Instead of the two-class case, let us consider a $c$-class case in general.
The method of You _et al._ [42] tries to optimize kernel parameters via the
following problem:

(20) $\operatorname*{maximize}~{}~{}Q_{1}Q_{2},$ (21) $\displaystyle Q_{1}$
$\displaystyle=\frac{2}{c(c-1)}\sum_{l=1}^{c-1}\sum_{l^{\prime}=l+1}^{c}\frac{\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{l^{\prime}}^{\phi})}{\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{l}^{\phi})+\operatorname{tr}(\bm{\Sigma}_{l^{\prime}}^{\phi}\bm{\Sigma}_{l^{\prime}}^{\phi})},$
(22) $\displaystyle Q_{2}$
$\displaystyle=\sum_{l=1}^{c-1}\sum_{l^{\prime}=l+1}^{c}\frac{n_{l}n_{l^{\prime}}}{n^{2}}\|\bm{\mu}^{\phi}_{l}-\bm{\mu}^{\phi}_{l^{\prime}}\|_{,}$

where $\bm{\Sigma}^{\phi}_{l}$ denotes the sample covariance matrix of the
data points of the $l$-th class in $\mathcal{H}$, that is,

(23) $\bm{\Sigma}^{\phi}_{l}=\frac{1}{n_{l}}\sum_{i\mid
y_{i}=l}(\phi(\bm{x}_{i})-\bm{\mu}^{\phi}_{l})(\phi(\bm{x}_{i})-\bm{\mu}^{\phi}_{l})^{\mathsf{T}}.$

$Q_{1}$ measures the homoscedasticity of data points of each class in the
feature space, and $Q_{2}$ measures the class-separability in the feature
space. By maximizing the product of these objectives, $Q_{1}Q_{2}$, we can
adjust the kernel parameters so that the data points in the kernel feature
space tend to be linearly separable.

###  2.3 Grassmann Kernels on Subspaces

The third building block of the proposed method is the notion of Grassmann
kernels on subspaces [see, e.g., 22]. Grassmann kernels enable us to compute
similarity between sets of vectors even if we do not have correspondence
between the vectors in different sets. In this work, we specifically use it to
define a similarity between two sets of dynamic modes.

One of the kernels introduced by Hamm and Lee [22] is called the _projection
kernel_. Let $\mathcal{B}$ denote some subspace of Euclidean space, and let
$\mathcal{B}_{1},\mathcal{B}_{2}\subset\mathcal{B}$ be subspaces of
$\mathcal{B}$. The projection kernel between those two subspaces,
$k_{\text{P}}(\mathcal{B}_{1},\mathcal{B}_{2})$, is defined as

$k_{\text{P}}(\mathcal{B}_{1},\mathcal{B}_{2})=\|\bm{B}_{1}^{\mathsf{H}}\bm{B}_{2}\|_{F}^{2},$

where $\bm{B}_{i}$ is the matrix whose columns comprise an orthonormal basis
of $\mathcal{B}_{i}$ for $i=1,2$. In fact, the projection kernel induces a
metric based on principal angles between subspaces, which is called the
projection metric [22]. It is defined as
$d_{\text{P}}(\mathcal{B}_{1},\mathcal{B}_{2})=(\sum_{d}\sin^{2}\alpha_{d})^{1/2}$,
where $\alpha_{d}$ denotes the $d$-th principal angle between
$\mathcal{B}_{1}$ and $\mathcal{B}_{2}$. Hence, we can understand the
projection kernel as a similarity measure based on principal angles between
subspaces.

##  3 Proposed Method

Figure 2: Schematic diagram of the two objectives, $f_{\text{DMD}}$ and
$f_{\text{KFD}}$, which appear in the optimization of the proposed method (see
Eq. 34). (a) We have pairs of time-series $\bm{X}_{i}$ and label $y_{i}$ (for
$i=1,\dots,n$) as a dataset (see Eq. 30–Eq. 33). (b) Given such data, we aim
to compute $\Theta_{i}$ and $\bm{W}_{i}$ for $i=1,\dots,n$. The elements of
$\Theta$ and the columns of $\bm{W}$ correspond to DMD’s eigenvalues and
dynamic modes, respectively. Note that the to-be-optimized variable is only
$\Theta$ because $\bm{W}$ is straightforwardly computed given $\Theta$ (see
Eq. 25). One of the objectives to optimize $\Theta$ is
$f_{\text{DMD}}(\Theta_{i})$, which measures the fitting error between
$\bm{X}_{i}$ and $\bm{W}\bm{V}_{\Theta_{i}}$ (see Eq. 36). (c) As
correspondence between columns of different $\bm{W}$’s cannot be established
in general, we compute the similarity of $\bm{W}$’s via their column space
$\operatorname{range}(\bm{W})$, which can be computed by SVD of $\bm{W}$. We
define a kernel function $k_{\text{DMS}}$ on such subspaces (see Eq. 26). (d)
The corresponding feature map is denoted as $\phi_{\text{DMS}}$. Another
objective to optimize $\Theta$ is $f_{\text{KFD}}(\Theta_{1:n})$, which
measures the homoscedasticity of data points in each class and their class-
separability in an RKHS $\mathcal{H}_{\text{DMS}}$ induced by $k_{\text{DMS}}$
(see Eq. 35). Note that the bold arrow at $f_{\text{KFD}}(\Theta_{1:n})$ in
the figure only depicts the class-separability as the homoscedasticity is more
difficult to depict.

In this section, we describe the details of our proposal of a method for
discovering distinctive coherent patterns from labeled spatio-temporal data
collections. The main idea of the proposed method lies in combining DMD
(reviewed in Section 2.1) and KFD (Section 2.2) using a Grassmann kernel
(Section 2.3). To this end, first we prepare a kernel to measure similarity
between sets of dynamic modes in Section 3.1. Then, we formulate the main
optimization problem Eq. 34 to compute distinctive dynamic modes in Section
3.2. Finally, we show numerical examples on synthetic data in Section 3.3. The
overall schematics of the proposed method is shown in Fig. 2.

###  3.1 Dynamic Mode Subspace Kernels for Time-Series Episodes

In computing dynamic modes, most variants of DMD cannot consider label
information of time-series episodes. We address this issue by developing a new
method combining virtues of DMD and KFD. We would like to compute dynamic
modes that are good both in the sense of DMD (i.e., curve fitting Eq. 13) and
in the sense of KFD (i.e., class homoscedasticity and separability Eq. 20–Eq.
21). To this end, we should prepare a kernel function between _sets of_
dynamic modes. An issue here is that in general, we cannot establish explicit
correspondences between dynamic modes computed from two different datasets.
Hence, simply computing similarity between dynamic mode matrices $\bm{W}$
cannot be a valid approach. Instead, we suggest using the projection kernel
that computes similarity between subspaces spanned by sets of dynamic modes.

First, let us formalize a subspace spanned by dynamic modes, namely a _dynamic
mode subspace_ (DMS):

######  Definition 3.1 (Dynamic mode subspace).

Let $\bm{X}\in\mathbb{C}^{p\times\tau}$ be a data matrix whose columns
comprise $p$-dimensional snapshots of a time-series episode of length $\tau$.
Let $\Theta=\theta_{1:r}=\\{\theta_{j}\in\mathbb{C}\mid j=1,\dots,r\\}$ be a
set of $r$ distinct complex values ($0<r\leq p$). We define the _dynamic mode
subspace_ $\mathcal{W}_{\Theta,\bm{X}}$ of $\bm{X}$ corresponding to
$\Theta=\theta_{1:r}$ as

(24) $\displaystyle\mathcal{W}_{\Theta,\bm{X}}$
$\displaystyle=\operatorname{range}\big{(}\bm{W}\big{)}\subset\mathbb{C}^{p},$
(25) $\displaystyle\bm{W}$
$\displaystyle=\bm{X}\big{(}\bm{V}_{\Theta}\big{)}^{\dagger},$

where $\bm{V}_{\Theta}$ is a Vandermonde matrix defined like in Eq. 10 with
set $\Theta=\theta_{1:r}$.

We note that, in Definition 3.1, the columns of
$\bm{W}=\bm{X}(\bm{V}_{\Theta})^{\dagger}\in\mathbb{C}^{p\times r}$ are the
dynamic modes computed from data $\bm{X}$ with $\Theta$ being the set of DMD
eigenvalues (in the sense of optimized DMD [9]), using the variable projection
technique in Eq. 12. In other words, dynamic modes appear in the definition of
DMS only implicitly, and we can parametrize a DMS only with a set of $r$
complex scalars, $\Theta$, instead of a dynamic mode matrix $\bm{W}$ of size
$p\times r$. This property of DMS plays a key role in the optimization problem
in Section 3.2. We then define the projection kernel on DMS as follows.

######  Definition 3.2 (Projection kernel on DMS).

Let $\bm{X}_{1}$ and $\bm{X}_{2}$ be two data matrices, and let $\Theta_{1}$
and $\Theta_{2}$ be two sets of $r$ distinct complex values. Let
$\mathcal{W}_{\Theta_{1},\bm{X}_{1}}$ and
$\mathcal{W}_{\Theta_{2},\bm{X}_{2}}$ be their DMS’s corresponding to
$\Theta_{1}$ and $\Theta_{2}$, respectively. The projection kernel [22]
between the two subspaces, $\mathcal{W}_{\Theta_{1},\bm{X}_{1}}$ and
$\mathcal{W}_{\Theta_{2},\bm{X}_{2}}$, is

(26) $\displaystyle
k_{\text{DMS}}(\mathcal{W}_{\Theta_{1},\bm{X}_{1}},\mathcal{W}_{\Theta_{2},\bm{X}_{2}})$
$\displaystyle=\|\bm{B}_{\Theta_{1},\bm{X}_{1}}^{\mathsf{H}}\bm{B}_{\Theta_{2},\bm{X}_{2}}\|_{F}^{2},$

where the columns of $\bm{B}_{\Theta_{1},\bm{X}_{1}}$ and
$\bm{B}_{\Theta_{2},\bm{X}_{2}}$ comprise an orthonormal basis of
$\mathcal{W}_{\Theta_{1},\bm{X}_{1}}$ and
$\mathcal{W}_{\Theta_{1},\bm{X}_{1}}$, respectively.

We denote the feature map corresponding to kernel function $k_{\text{DMS}}$ by

$\phi_{\text{DMS}}:\mathcal{W}_{\Theta,\bm{X}}\mapsto\phi_{\text{DMS}}(\mathcal{W}_{\Theta,\bm{X}}).$

Moreover, we denote the RKHS induced by $k_{\text{DMS}}$ by
$\mathcal{H}_{\text{DMS}}$. We can compute $\bm{B}_{\Theta,\bm{X}}$ by the
singular value decomposition (SVD) of
$\bm{W}=\bm{X}(\bm{V}_{\Theta})^{\dagger}$; that is, $\bm{B}_{\Theta,\bm{X}}$
can be a matrix comprising left singular vectors of
$\bm{X}(\bm{V}_{\Theta})^{\dagger}$.

The derivative of $k_{\text{DMS}}$ with regard to an element of $\Theta_{1}$
or $\Theta_{2}$ can be computed as follows. Let $\theta_{1,j}$ denote the
$j$-th element of $\Theta_{1}$. Then222When $\bm{Y}$ is an $n\times m$ matrix,
we define both $\partial x/\partial\bm{Y}$ and $\partial\bm{Y}/\partial x$ as
an $n\times m$ matrix.,

(27) $\frac{\partial
k_{\text{DMS}}(\mathcal{W}_{\Theta_{1},\bm{X}_{1}},\mathcal{W}_{\Theta_{2},\bm{X}_{2}})}{\partial\theta_{1,j}}=\operatorname{sum}\left(\frac{\partial
k_{\text{DMS}}}{\partial\bm{V}_{\Theta_{1}}}\circ\frac{\partial\bm{V}_{\Theta_{1}}}{\partial\theta_{1,j}}\right),$

where $\circ$ is the Hadamard product, and $\operatorname{sum}(\bm{Y})$ means
the summation over all the elements of matrix $\bm{Y}$. Suppose that
$\bm{B}_{\Theta_{1},\bm{X}_{1}}$ is obtained by the SVD of
$\bm{X}_{1}(\bm{V}_{\Theta_{1}})^{\dagger}$. Then, we have

(28) $\left(\frac{\partial
k_{\text{DMS}}}{\partial\bm{V}_{\Theta_{1}}}\right)^{\mathsf{T}}=(\bm{V}_{\Theta_{1}})^{\dagger}(\bm{V}_{\Theta_{1}}^{\mathsf{H}})^{\dagger}\bm{C}_{1}^{\mathsf{H}}(\bm{I}-\bm{V}_{\Theta_{1}}(\bm{V}_{\Theta_{1}})^{\dagger})^{\mathsf{H}}\\\
+(\bm{I}-(\bm{V}_{\Theta_{1}})^{\dagger}\bm{V}_{\Theta_{1}})^{\mathsf{H}}\bm{C}_{1}^{\mathsf{H}}(\bm{V}_{\Theta_{1}}^{\mathsf{H}})^{\dagger}(\bm{V}_{\Theta_{1}})^{\dagger}-(\bm{V}_{\Theta_{1}})^{\dagger}\bm{C}_{1}(\bm{V}_{\Theta_{1}})^{\dagger},$

where

(29)
$\bm{C}_{1}=(\bm{X}_{1}(\bm{V}_{\Theta_{1}})^{\dagger})^{\dagger}\bm{B}_{\Theta_{2},\bm{X}_{2}}\bm{B}_{\Theta_{2},\bm{X}_{2}}^{\mathsf{H}}(\bm{I}-\bm{B}_{\Theta_{1},\bm{X}_{1}}\bm{B}_{\Theta_{1},\bm{X}_{1}}^{\mathsf{H}})\bm{X}_{1}.$

The above discussion analogously applies to the computation of $\partial
k_{\text{DMS}}/\partial\theta_{2,j}$.

###  3.2 Discriminant Dynamic Mode Decomposition

This section describes our problem setting more formally and the main part of
the proposed method.

####  3.2.1 Problem Setting

We suppose that we have a _dataset_ $\mathcal{D}$:

(30) $\mathcal{D}=\\{(X_{1},y_{1}),\dots,(X_{n},y_{n})\\},$

whose element is a pair of an _episode_ $X_{i}$ and a _label_ $y_{i}$. An
episode $X_{i}$ is a sequence of observation vectors, i.e.,

(31) $X_{i}=(\bm{x}_{i,1},\dots,\bm{x}_{i,\tau_{i}}),$

where $\bm{x}_{i,t}\in\mathbb{R}^{p}$ ($t=1,\dots,\tau_{i}$) is an observation
at time-step $t$ and also called a _snapshot_. $\tau_{i}$ is the length of the
episode $X_{i}$, and $p$ is the dimensionality of a snapshot. We may denote
the episode $X_{i}$ also by a matrix

(32)
$\bm{X}_{i}=\begin{bmatrix}\bm{x}_{i,1}&\cdots&\bm{x}_{i,\tau_{i}}\end{bmatrix}\in\mathbb{R}^{p\times\tau_{i}}.$

In a dataset, each episode is associated with a label

(33) $y_{i}\in\\{1,\ldots,c\\}.$

The label of an episode represents some distinctive property (within $c$
classes) of the episode.

Given such a collection of labeled time-series episodes as a dataset, we would
like to obtain sets of complex values
$\Theta_{i}=\\{\theta_{i,1},\dots,\theta_{i,r}\\}$ (for $i=1,\dots,n$) such
that:

  1. 1.

The labeled episodes are well separated in KFD’s sense; and

  2. 2.

The snapshots of each episode are well fitted in the sense of the optimized
DMD’s sense.

Here we note that we only deal with the sets of time-evolution parameters
$\Theta_{i}$ (for $i=1,\dots,n$) as optimization variables because once they
are determined, we can immediately compute the corresponding dynamic modes
$\bm{W}$ by substituting $\theta_{i,1:r}$ to $\lambda^{*}_{1:r}$ in Eq. 12.

####  3.2.2 Formulation

We formulate the proposed method, _discriminant DMD_ , as the following
optimization problem:

(34)
$\begin{gathered}\underset{\Theta_{1:n}}{\text{minimize}}~{}~{}\frac{n^{-1}\sum_{i=1}^{n}f_{\text{DMD}}(\Theta_{i})}{f_{\text{KFD}}(\Theta_{1:n})^{\alpha}+\epsilon},\end{gathered}$

where $\alpha\geq 0$ is a hyperparameter that balances the importance of
$f_{\text{KFD}}$ and $f_{\text{DMD}}$, and $\epsilon\geq 0$ is a small number
for numerical stability. The two terms in Eq. 34 are to achieve the purposes
of the discriminant DMD listed above and defined as follows.

##### KFD term

For letting the DMS projection kernel well separate data points (i.e., labeled
episodes), we try to maximize

(35) $f_{\text{KFD}}(\Theta_{1:n})=Q_{1}Q_{2},$

with respect to $\Theta_{1:n}$. We define $Q_{1}$ and $Q_{2}$ like Eq. 21 and
Eq. 22, respectively, with $\phi$ being $\phi_{\text{DMS}}$, the feature map
corresponding to $k_{\text{DMS}}$. Hence, $Q_{1}$ and $Q_{2}$ depend on
$\Theta_{1:n}$ via the arguments of $k_{\text{DMS}}$. Again let us emphasize
that a data point of KFD here corresponds to an episode (_not_ a snapshot), so
kernel function $k_{\text{DMS}}$ will be called $O(n^{2})$ times.

##### DMD term

For making the elements of $\Theta_{i}$ good for fitting $\bm{X}_{i}$ in the
optimized DMD’s sense (i.e., Eq. 13), for $i=1,\dots,n$, we try to minimize

(36)
$f_{\text{DMD}}(\Theta_{i})=\frac{1}{\tau_{i}}\big{\|}\bm{X}_{i}-\bm{X}_{i}(\bm{V}_{\Theta_{i}})^{\dagger}\bm{V}_{\Theta_{i}}\big{\|}_{F}^{2}.$

It is exactly the same with the one used in the optimized DMD with the
variable projection [3].

Algorithm 1 Discriminant DMD

0: Dataset $\mathcal{D}=\\{(X_{i},y_{i})\mid i=1,\dots,n\\}$ (i.e., a set of
pairs of an episode $X_{i}=(\bm{x}_{i,1},\dots,\bm{x}_{i,\tau_{i}})$ and a
label $y_{i}$) and hyperparameter $\alpha$

0: Sets of time-evolution parameters $\\{\Theta_{i}\mid i=1,\dots,n\\}$ (i.e.,
what corresponds to DMD’s eigenvalues) and sets of distinctive dynamic modes
$\\{\bm{W}_{i}\mid i=1,\dots,n\\}$ (each column of $\bm{W}_{i}$ is a dynamic
mode)

1: for $i=1\,\ldots\,n$ do

2: $\Theta_{i},\bm{W}_{i}\leftarrow\text{StandardDMD}(X_{i})$ {[39] or [9; 1];
see Section 2.1}

3: end for

4:
$\Theta_{1},\dots,\Theta_{n}\leftarrow\arg\min_{\Theta_{1}:n}\frac{n^{-1}\sum_{i=1}^{n}f_{\text{DMD}}(\Theta_{i})}{f_{\text{KFD}}(\Theta_{1:n})^{\alpha}+\epsilon}$
{See Eq. 34, Eq. 35, and Eq. 36}

5: for $i=1\,\ldots\,n$ do

6: $\bm{W}_{i}\leftarrow$ compute Eq. 12 with $\lambda^{*}$ being $\Theta_{i}$

7: end for

####  3.2.3 Optimization

Local solutions of Eq. 34 can be found using gradient methods. In the
numerical experiments below, we used a quasi-Newton optimizer with initial
values computed by the standard DMD algorithm based on eigendecomposition
[39]. We present the details on the gradient of the objective in Appendix A.
The overall algorithm of the discriminant DMD is summarized in Algorithm 1.

###  3.3 Numerical Example

|  |
---|---|---
$\bm{w}_{\mathrm{d},y=1}$ | $\bm{w}_{\mathrm{d},y=2}$ | $\bm{w}_{\mathrm{c}}$
(a)

(b)

standard DMD
[e.g., 33]

supervised DMD
[14]

discriminant DMD
(ours)

(c)

Figure 3: (a) True dynamic modes used to generate data as (37). (b) Values of
the two objective functions achieved by the optimization with different
$\alpha$’s. (c, upper row) Averages of the computed dynamic modes for each
label. (c, lower row) Visualization of episodes using classical
multidimensional scaling with the distance function induced from
$k_{\text{DMS}}$. Each point corresponds to each episode. Blue circles and red
squares denote $y=1$ and $y=2$, respectively.

We show a numerical example with synthetic data.

##### Data Generation

We use a synthetic dataset comprising $n=20$ episodes of length $\tau=100$,
each of which is generated as

(37)
$\bm{x}_{i,t}=\lambda_{\mathrm{d},i}^{t}\bm{w}_{\mathrm{d},y_{i}}+\lambda_{\mathrm{c},i}^{t}\bm{w}_{\mathrm{c}}+\bm{e}_{i,t}\in\mathbb{C}^{p},\quad\text{for}\quad
i=1,\dots,n\quad\text{and}\quad t=1,\dots,\tau,$

where the $i$-th episode is defined as $(\bm{x}_{i,1},\dots,\bm{x}_{i,\tau})$,
and $y_{i}\in\\{1,2\\}$ is the label of the $i$-th episode. We set $y_{i}=1$
for $i=1,\dots,10$ and $y_{i}=2$ for $i=11,\dots,20$. Equation 37 says that
each episode is generated by two dynamic modes and noise. One of the true
dynamic modes, $\bm{w}_{\mathrm{d},y}\in\mathbb{R}^{p}$, depends on label $y$,
whereas the other does not. Our interest is to extract this
$\bm{w}_{\mathrm{d},y}$ as a _distinctive_ dynamic mode. In this example, we
defined $\bm{w}_{\mathrm{d},y=1}$ and $\bm{w}_{\mathrm{d},y=2}$ as the left
two images shown in Fig. 3a; each image is of size $10\times 10$, and thus
$p=100$. The corresponding time-evolution parameter,
$\lambda_{\mathrm{d},i}=\exp(-\gamma_{\mathrm{d}}+\sqrt{-1}\omega_{\mathrm{d},i})\in\mathbb{C}$
was randomly generated for each episode by drawing $\omega_{\mathrm{d},i}$
from the uniform distribution on $[0,1]$ with $\gamma_{\mathrm{d}}=0.1$.
Vector $\bm{w}_{\mathrm{c}}$ is a non-distinctive dynamic mode used commonly
for both $y=1$ and $y=2$ (defined as a Gaussian-like image shown in Fig. 3a),
and the corresponding time-evolution parameter, $\lambda_{\mathrm{c},i}$, was
randomly generated for each episode similarly to $\lambda_{\mathrm{d},i}$.
Finally, $\bm{e}\in\mathbb{C}^{p}$ denotes a noise vector whose elements were
independently drawn from the zero-mean circularly-symmetric complex normal
distribution with standard deviation $0.05$.

##### Results

We applied the eigendecomposition-based standard DMD algorithm [39],
supervised DMD ([14]; we review it in Section 4), and the proposed
discriminant DMD to the synthetic dataset. Each method was set so that it
would compute a dynamic mode (i.e., $r=1$)333The true value of $r$ is
obviously $2$, but we set $r=1$. We intentionally misspecified $r$ because if
$r$ had been set to a true value, the problem would have been too easy.. Let
$\hat{\bm{w}}_{i}$ be the dynamic mode computed for the $i$-th episode. In
Fig. 3c, we show the averages of the dynamic modes computed by the three
methods for each label (i.e., $\frac{1}{10}\sum_{i=1}^{10}\hat{\bm{w}}_{i}$
and $\frac{1}{10}\sum_{i=11}^{20}\hat{\bm{w}}_{i}$). We can observe that the
modes computed by the discriminant DMD are less occluded by the non-
distinctive mode. Moreover, we visualize the relation between the 20 episodes
using the classical multidimensional scaling (MDS) with the distance function
induced from $k_{\text{DMS}}$. We can see that the episodes with different
labels are well separated in the feature space induced from $k_{\text{DMS}}$.
Furthermore, we show the values of the two objectives of the discriminant DMD,
$f_{\text{DMD}}$ and $f_{\text{KFD}}$, achieved by solving Eq. 34 with
different $\alpha$’s ($\alpha=0,0.2,0.4,\dots,1.2$), in Fig. 3b. We can
observe that the value of $f_{\text{KFD}}$ improves (i.e., increases), trading
off the reconstruction error, $f_{\text{DMD}}$.

##  4 Related Work

###  4.1 Extraction of Distinctive Patterns

Informative data patterns are often obtained as a by-product of dimensionality
reduction, and supervised dimensionality reduction has been studied in various
contexts. Representative examples include Fisher discriminant [13; 30], common
spatial patterns [32], partial least squares [41], neighborhood component
analysis [20], sufficient dimensionality reduction [16], and supervised PCA
[4]. These methods are promising tools for extracting informative patterns
from labeled non-sequential data, but are not appropriate for our problem,
where each data point is a sequence with a label.

Dimensionality reduction on collections of labeled sequences was addressed in
the line of studies by Su _et al._ [35; 36; 37], in which they compute
discriminative ordered templates of multivariate time-series. However, such
templates do not necessarily match what we would like to compute in this work,
that is, spatio-temporal coherent patterns.

Supervised DMD [14] is seemingly very similar but is fundamentally different
from our method in its purpose and technique. First, whereas supervised DMD
aims to extract _common_ dynamic modes shared by episodes of same labels, our
method extracts _distinctive_ dynamic modes that well distinguish episodes of
different labels. More importantly, there is a critical technical difference;
whereas supervised DMD modifies coefficients multiplied to dynamic modes (and
thus dynamic modes remain unchanged), our method directly modifies dynamic
modes (via eigenvalues) in accordance with label information. In other words,
supervised DMD adapts the _usage_ of dynamic modes to label information,
whilst our method adapts the _shape_ of dynamic modes itself.

###  4.2 Spectral Analysis of Koopman Operator

Recall that the Koopman operator is infinite-dimensional in general, and
consequently, the corresponding modal decomposition is taken into infinitely-
many terms (see Eq. 4 and Eq. 5). DMD is sometimes regarded as a “finite-
dimensional approximation” of such a decomposition in some sense, but it is
not straightforward to characterize which $r$ modes are actually targeted by
DMD out of the infinite number of modes of the Koopman operator. Consequently,
there is no definitive way of determining the number of modes computed by DMD,
$r$, nor choosing good (in some sense) dynamic modes out of computed ones.

A working practice is to determine $r$ by the effective rank of a dataset via
SVD and then choose representative modes according to their energy [33; 34;
39]. A more systematic approach is to consider additional coefficients
multiplied to dynamics modes and to optimize them with sparsity regularization
[24; 14] so that we can choose a relatively small number of modes effectively
used in the decomposition.

Apart from DMD-like methods, other types of approaches to approximating the
Koopman operator also equip strategy to pick up specific spectral components.
A method called generalized Laplace analysis [7] is a rigorous way to
approximate Koopman modes. It is inherently spared for the need of mode
selection because we give, in advance, eigenvalues for which modes are to be
computed. Giannakis [19] and Das and Giannakis [12] suggest computing the
spectra of Koopman operators by a Galerkin method on eigenspaces of the
Laplace–Beltrami operator of data space. In this method, they refer to the
Dirichlet energy of eigenfunctions of Koopman operator, which characterize
their “roughness,” to select valid eigenvalues.

Given the various approaches to selecting good spectral components, we believe
that our method in this work provides a new insight to this end. That is, the
proposed method is informed from label information for modifying (not
selecting, though) eigenvalues and dynamic modes. We have not rigorously
characterized how the proposed method has connection to the Koopman operator
theory yet, but this is an interesting research direction.

##  5 Experiments

We show the application of the discriminant DMD to four real-world datasets.
The first two experiments in this section are mainly for understanding the
behavior of the proposed method; we observe that the method can extract
patterns that well reconstruct data _and_ distinguish different labels. The
last two experiments are rather for demonstrating other possible application
scenes of the proposed method.

###  5.1 House Temperature Data

####  5.1.1 Dataset

We applied the proposed discriminant DMD to measurements of temperature in a
house, which exhibit particular yet varying patterns each day. From the
original dataset444github.com/LuisM78/Appliances-energy-prediction-data [8],
we extracted the temperature measured in the eight rooms (kitchen, living
room, laundry room, office room, bathroom, ironing room, teenager room, and
parent room) of a house and at the nearby weather station. We subtracted the
measurements at the weather station from ones in the eight rooms, so a
snapshot is $p=8$ dimensional. We partitioned the original sequence by day
into $n=136$ episodes (i.e., $136$ days). We performed 6-point moving average
and subsampling by 1/3 on the original data, so each episode is of length
$\tau=48$ with the measurement interval being $30$ minutes. We labeled each
episode as $y=\text{``weekday''}$ (i.e., neither weekend nor holiday) or
$y=\text{``holiday''}$ (i.e., Saturdays, Sundays, and national holidays).
There were $n_{\text{weekday}}=94$ weekdays and $n_{\text{holiday}}=42$
holidays.

####  5.1.2 Configuration

We applied the discriminant DMD with $r$ (the number of dynamic modes) being
$r=6$ and $\alpha$ (the hyperparameter that balances reconstruction and class
separation) varying from $\alpha=0$ to $\alpha=1$. Setting $\alpha=0$
corresponds to the standard optimized DMD algorithm [9], whereas $\alpha>0$
realizes our proposal for incorporating label information. As a further
baseline, we also show the results obtained by PCA on each sequence.

MDS All DMD eigenvalues Median of dominant modes

Figure 4: Results on the house temperature dataset. (_left_) Visualization of
episodes via MDS with the projection kernel. Each point corresponds to each
episode of the dataset. (_center_) DMD eigenvalues. All the eigenvalues
computed from all episodes are plotted altogether. (_right_) Visualization of
average dominant modes overlaid on the floor plan. The letters denote: (A)
kitchen, (B) living room, (C) laundry room, (D) office room, (E) bathroom, (F)
ironing room, (G) teenager room, and (H) parent room. Best viewed in color.

####  5.1.3 Results

Figure 5: Median of NRMSE of the reconstruction by the discriminant DMD with
different values of $\alpha$, for the house temperature dataset. $\alpha=0$
corresponds to the standard DMD, and $\alpha>0$ corresponds to the
discriminant DMD. The reconstruction errors by PCA with one, two, or three PCs
are also shown for comparison.

In the left column of Fig. 4, we show the two-dimensional embedding of the
episodes by MDS with the distance between episodes computed via the projection
kernel. The projection kernel was computed between sets of dynamic modes for
DMDs and sets of principal components (PCs) for PCA555Projection kernel on PCs
is similar to its original usage [22].. We can observe that the episodes with
different labels are well separated with the discriminant DMD ($\alpha>0$) as
expected from the formulation.

In the center column of Fig. 4, we plot the values of $\lambda$ computed on
all the episodes. This quantity, $\lambda$, appears both in standard DMD and
our discriminant DMD (see Eq. 8 and Eq. 11), and we refer to it in both cases
simply as DMD eigenvalues. The DMD eigenvalues of the two labels distribute
quite similarly when $\alpha=0$ (i.e., no label information is reflected),
whereas the distributions are distinctively different between the two labels
when $\alpha>0$.

In the right column of Fig. 4, we visualize the dynamic modes666As we compute
multiple dynamic modes for all the episodes, we examined summary statistics of
the dynamic modes for each label. First, we chose the dominant mode with the
largest norm $\|\bm{w}\|$ for each episode. Then, we calculated the element-
wise median of the dominant modes over all the episodes for each label. by
painting each room of a floor plan of the house according to the dynamic
modes. With the discriminant DMD (i.e., $\alpha>0$), the difference between
$y=\text{``weekday''}$ and $y=\text{``holiday''}$ are consistently emphasized
in some rooms such as B (living room) and F (ironing room). We can utilize
such information, for example, for designing a controller of air conditioning
systems.

In Fig. 5, we show the median of the normalized root mean square error (NRMSE)
between the original time-series sequences and the reconstructed ones by the
discriminant DMD with different $\alpha$ as well as by PCA. We can observe
that the reconstruction by the discriminant DMD achieves the similar
reconstruction errors as PCA with two or three PCs, regardless of the value of
$\alpha$. Recall that we set $r=6$ for the discriminant DMD. Because a DMD-
like decomposition usually results in (nearly) complex conjugate pairs of
modes for real-valued data, it is reasonable to expect that the reconstruction
capability of DMD-like decomposition with $r=6$ modes (i.e., three pairs)
roughly corresponds to that of PCA with three PCs. The result of $\alpha=0$
(i.e., standard DMD) in Fig. 5 fits this expectation almost completely.
Moreover, we can see the reconstruction capability with $\alpha>0$ (i.e.,
discriminant DMD) is not severely sacrificed, which supports the validity of
the proposed method in the sense that it can compute distinctive _and_
representative dynamic modes.

###  5.2 Motion Capture Data

####  5.2.1 Dataset

Figure 6: Seven angles used as features of the motion capture data.

We also applied the discriminant DMD to measurements of human locomotion,
which has already been a target of analysis of standard DMD techniques [15].
We prepared a dataset from the CMU Graphics Lab Motion Capture
Database777mocap.cs.cmu.edu. We chose $n_{\text{walk}}=40$ episodes of “walk”
motion and $n_{\text{run/jog}}=40$ episodes of “run/jog” motion and
constructed a dataset with these $n=80$ episodes. We only used the first
$\tau=128$ timesteps (about one second) of each sequence. We transformed the
original data into the angles of trunk, thigh, shank, and foot elevation in
the sagittal plane (schematically shown in Fig. 6) for extracting
intersegmental motion coherence [17]. Consequently, each episode is a $p=7$
dimensional time-series sequence. As preprocessing, we normalized each feature
and applied a low-pass filter at $6$ [Hz].

####  5.2.2 Configuration

The configuration is similar to that in the previous experiment; we applied
the discriminant DMD with $r=4$, and the hyperparameter $\alpha$ was varied
from $\alpha=0$ to $\alpha=1$. We also show the results obtained by PCA on
each sequence.

####  5.2.3 Results

In the left column of Fig. 7, we show the two-dimensional embedding of the
episodes computed similarly to the previous experiment, that is, via MDS with
projection kernels. The overlap of the episodes of different labels is
alleviated by the proposed method with $\alpha>0$. In the center and right
columns of Fig. 7, we show all the DMD eigenvalues and the medians of dominant
modes, respectively, similarly to the previous experiment. The distinction
between $y=\text{``walk''}$ and $y=\text{``run/jog''}$ is less obvious
visually in this case, but such information, possibly combined with domain
knowledge, may play an important role in gait analysis.

In Fig. 8, we show the median of the NRMSE between the original sequences and
the reconstructed ones, with which we can obtain similar observations as in
the previous experiment.

MDS All DMD eigenvalues Med. of dominant modes

feature number

Figure 7: Results on the motion capture dataset. (_left_) Visualization of
episodes via MDS with the projection kernel. Each point corresponds to each
episode of the dataset. (_center_) DMD eigenvalues. All the eigenvalues
computed from all episodes are plotted altogether. (_right_) Visualization of
average dominant modes. The feature numbers correspond to ones shown in Fig.
6.

Figure 8: Median of NRMSE of the reconstruction by the discriminant DMD with
different values of $\alpha$, for the motion capture dataset. $\alpha=0$
corresponds to standard DMD, and $\alpha>0$ corresponds to the discriminant
DMD. The reconstruction errors by PCA with one, two, or three PCs are also
shown for comparison.

###  5.3 Taxi Trip Data

---
(a) weekday

---
(b) weekend

Figure 9: Examples of the distinctive coherent patterns extracted from the
taxi trip data. Areas in green, light blue, yellow, and blue indicate
Manhattan, Brooklyn, Queens, and Bronx, respectively. (_upper_) Magnitude of
the coherent patterns. (_lower_) Temporal profiles of corresponding dynamics.

####  5.3.1 Dataset

Certain kinds of spatio-temporal data have characteristic structures we can
anticipate based on knowledge of social and biological rhythms (e.g., days,
months, and years). The discriminant DMD will be remarkably useful for such
data because we can examine specific coherent structures focusing on their
dynamical property, that is, frequencies of rhythms. As an example of such a
usage, we applied the discriminant DMD to records of taxi trips. The dataset
is a record of the numbers of taxies that departed and arrived between four
boroughs of New York City (Manhattan, Brooklyn, Queens, and Bronx)888We used
trip records of boro taxies that are allowed to pick up passengers in outer
boroughs. www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page and thus
comprises $p=4^{2}$ dimensional time-series sequences. We used the original
sequence for ten weeks and partitioned it into $y=\text{``weekday''}$ and
$y=\text{``weekend''}$ episodes for each week, which resulted in a dataset
comprising $n=20$ episodes (i.e., $n_{\text{weekday}}=n_{\text{weekend}}=10$).
As the data record the numbers of departures and arrivals within each 30
minutes, the length of each episode is $\tau_{\text{weekday}}=240$ or
$\tau_{\text{weekend}}=96$. As preprocessing, we performed 24-point (i.e.,
half day) moving average and subtracted the mean value from each feature of
each episode.

####  5.3.2 Results

We set $r=4$ and $\alpha=1$. We examined the dynamic modes whose eigenvalues
$\lambda$ roughly coincide with daily periodicity, which is one of the
expected dynamical characteristics of the traffic. In Fig. 9, we exemplify
such dynamic modes extracted by the discriminant DMD from a pair of weekday
and weekend episodes. We also illustrate the corresponding temporal profiles,
$|\varphi_{t}|=|\lambda^{t-1}(\bm{z}^{\mathsf{H}}\bm{x}_{1})|$ (see Eq. 8) in
Fig. 9. In the weekend (Fig. 9b), almost every traffic except between Brooklyn
and Bronx exhibits large amplitudes with the daily periodicity. In contrast,
in the weekdays (Fig. 9a), the movements between Brooklyn and Queens and ones
from them to Bronx have larger amplitudes. Though directly interpreting such
observations requires intensive domain knowledge, they can be useful, for
example, for understanding the mechanism of the traffic and for making road-
building plans.

###  5.4 EEG Data

---
(a) left-hand

---
(b) right-hand Figure 10: Examples of the distinctive coherent patterns and
their dynamics extracted from the EEG data. Representative ones are shown for
each of the two labels: (a) left- and (b) right-hand motor imagery. The
patterns in the upper row similarly follow the slightly decaying oscillations,
and ones in the lower row follow the rapidly converging dynamics.

####  5.4.1 Dataset

We applied the discriminant DMD to electroencephalography (EEG) data to show a
possibility of application in such a biological domain. We used the dataset
publicly provided in BCI Competition IV 2a999www.bbci.de/competition/iv/. The
original data was recorded using twenty-two Ag/AgCl electrodes with sampling
rate of 250 [Hz]. We eliminated noises by applying the band-pass filter
between 0.5 [Hz] and 100 [Hz] and the spatial filter with the surface
Laplacian [28]. We used ten episodes of the same subject behaving two motor
imagery tasks: imagination of movement of the left hand or the right hand.
Hence, the dataset comprises $n=20$ episodes (i.e., $n_{\text{left-
hand}}=n_{\text{right-hand}}=10$). For each sequence, we used only a part of
length $\tau=250$ (corresponding to $1$ [sec]) from one second after the task
started.

####  5.4.2 Results

We set $r=4$ and $\alpha=1$. As with the previous experiment, we show
representative examples of distinctive dynamic modes and their temporal
profiles in Fig. 10b. While such information alone is not necessarily
sufficient for analyzing EEG signals, it may be informative when combined with
existing domain knowledge. For example, we may find some similarity between
the patterns in Fig. 10b and the ones obtained in the previous work [5] using
the method of common spatial pattern.

Another interesting observation here is that those pairs of distinctive
patterns have different dynamical properties. That is, in Fig. 10b, the upper-
row patterns similarly follow slightly decaying oscillations, and the lower-
row patterns have the rapidly converging temporal profiles. One of the
advantages of the proposed method is that it can extract dynamics of spatial
coherent patterns as shown in this example.

##  6 Discussion

The numerical results in the previous section show the utility of the proposed
method, discriminant DMD, for extracting sets of distinctive coherent patterns
(i.e., oscillating spatial patterns) that well reconstruct time-series signals
_and_ distinguish labels. Such information extraction is useful in
applications where we want to understand data characteristics, craft features
for classification and regression, and design controllers. For example, the
distinctive patterns of house temperature (Section 5.1) can be utilized in
designing the location and controllers of air conditioners. Analysis of human
motion capture data with label information (Section 5.2) can help more precise
data-driven understanding of human locomotion (see, e.g., [15]). Nonetheless,
the discriminant DMD itself can hardly be sufficient for _solving_ some task
completely, as so are any other data analysis methods such as PCA and DMD. We
need to combine the discriminant DMD with other subsequent methodologies for
more practical utility, and investigating such possibility in various domains
is a promising direction of future study.

One of notable limitations of the current discriminant DMD is that it does not
provide a principled way to determine which dynamic modes are _the most_
important both for representing sequences and distinguishing labels. As we
define the similarity between episodes via a kernel function on _sets_ of
dynamic modes, it is usually challenging to tell exactly which pairs of
dynamic modes are solely to be focused out of $r$ modes computed by the
method. Hence, one of the primal usages of the current method will be in
exploratory data analysis, where we can interpret the results based on our
domain knowledge. For example, in Section 5.3, we picked up the pair of
dynamic modes based on the knowledge that the data should indicate daily
periodicity as its characteristics. However, it will be more useful if there
is a principled way to rank importance of dynamic modes for data analysis with
less prior knowledge.

##  7 Conclusion

In this work, we developed a method for discovering spatio-temporal coherent
patterns from labeled data collections, namely _discriminant DMD_. The
discriminant DMD computes distinctive coherent patterns that contribute to
major differences of dynamics with different labels by optimizing the
objective that takes the reconstruction goodness of DMD and the class-
separation goodness of discriminant analysis into account. We have
demonstrated applications of the discriminant DMD using four different types
of real-world datasets, with which we empirically validated that the proposed
method extracts spatial patterns that well reconstruct data _and_ distinguish
different labels. Such pattern extraction is useful for exploratory data
analysis towards understanding spatio-temporal data. Important directions of
future research include the investigation of the utility of the discriminant
DMD in further subsequent tasks such as feature crafting and controller
design.

##  Appendix A Gradient of Objective Function

The gradient of the objective function in (34), $f$, with regard to
$\theta_{i,j}$ (for $i=1,\dots,n$ and $j=1,\dots,r$), is obtained as follows.
The gradient is computed from derivatives of $f$ as

(38) $\nabla_{\theta_{i,j}}f=2\frac{\partial
f}{\partial\bar{\theta}_{i,j}}=2\overline{\frac{\partial
f}{\partial\theta_{i,j}}},$

where $\bar{\cdot}$ means the complex conjugate. The derivative is

(39) $\frac{\partial
f}{\partial\theta_{i,j}}=\frac{\partial}{\partial\theta_{i,j}}\left(\frac{f_{\text{DMD}}(\Theta_{1:n})}{f_{\text{KFD}}(\Theta_{1:n})^{\alpha}+\epsilon}\right)=\frac{(f_{\text{KFD}}^{\alpha}+\epsilon)\frac{\partial
f_{\text{DMD}}}{\partial\theta_{i,j}}-\alpha
f_{\text{DMD}}f_{\text{KFD}}^{\alpha-1}\frac{\partial
f_{\text{KFD}}}{\partial\theta_{i,j}}}{\big{(}f_{\text{KFD}}^{\alpha}+\epsilon\big{)}^{2}}.$

The first of the two derivatives in (39), $\frac{\partial
f_{\text{DMD}}}{\partial\theta_{i,j}}$ is

(40) $\frac{\partial
f_{\text{DMD}}}{\partial\theta_{i,j}}=\operatorname{sum}\left(\frac{\partial
f_{\text{DMD}}}{\partial\bm{V}_{\Theta_{i}}}\circ\frac{\partial\bm{V}_{\Theta_{i}}}{\partial\theta_{i,j}}\right),$

where $\operatorname{sum}(\bm{A})$ denotes the summation of all the elements
of matrix $\bm{A}$, and

(41) $\frac{\partial
f_{\text{DMD}}}{\partial\bm{V}_{\Theta_{i}}}=\frac{1}{n\tau_{i}}\Big{(}(\bm{V}_{\Theta_{i}}^{\dagger}\bm{V}_{\Theta_{i}}-\bm{I})\bm{X}_{i}^{\mathsf{H}}\bm{X}_{i}\bm{V}_{\Theta_{i}}^{\dagger}\Big{)}^{\mathsf{T}},$

and

(42)
$\frac{\partial\bm{V}_{\Theta_{i}}}{\partial\theta_{i,j}}=\begin{bmatrix}0&0&0&\cdots&0\\\
\vdots&\vdots&\vdots&&\vdots\\\
0&1&2\theta_{j}&\cdots&(\tau_{i}-1)\theta_{j}^{\tau_{i}-2}\\\
\vdots&\vdots&\vdots&&\vdots\\\ 0&0&0&\cdots&0\end{bmatrix}.$

The other derivative in (39), $\frac{\partial
f_{\text{KFD}}}{\partial\theta_{i,j}}$, is obviously

(43) $\frac{\partial f_{\text{KFD}}}{\partial\theta_{i,j}}=\frac{\partial
Q_{1}}{\partial\theta_{i,j}}Q_{2}+\frac{\partial
Q_{2}}{\partial\theta_{i,j}}Q_{1}.$

Below we basically follow the deviation by You _et al._ [42] to complete this
calculation.

As for $Q_{1}$ in (43), we have

(44) $\frac{\partial
Q_{1}}{\partial\theta_{i,j}}=\frac{2}{c(c-1)}\sum_{l=1}^{c-1}\sum_{m=l+1}^{c}\left\\{\frac{\frac{\partial\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{m}^{\phi})}{\partial\theta_{i,j}}\Big{(}\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{l}^{\phi})+\operatorname{tr}(\bm{\Sigma}_{m}^{\phi}\bm{\Sigma}_{m}^{\phi})\Big{)}}{\Big{(}\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{l}^{\phi})+\operatorname{tr}(\bm{\Sigma}_{m}^{\phi}\bm{\Sigma}_{m}^{\phi})\Big{)}^{2}}\right.\\\
\left.-\frac{\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{m}^{\phi})\Big{(}\frac{\partial\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{l}^{\phi})}{\partial\theta_{i,j}}+\frac{\partial\operatorname{tr}(\bm{\Sigma}_{m}^{\phi}\bm{\Sigma}_{m}^{\phi})}{\partial\theta_{i,j}}\Big{)}}{\Big{(}\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{l}^{\phi})+\operatorname{tr}(\bm{\Sigma}_{m}^{\phi}\bm{\Sigma}_{m}^{\phi})\Big{)}^{2}}\right\\},$

where $\bm{\Sigma}_{l}^{\phi}$ denotes the sample covariance matrix of the
$l$-th class ($l\in[1,c]$) in the feature space as in section 2.2.

Recall that we gave an index to the sequences in a dataset $\mathcal{D}$ from
$1$ to $n$ (see (30)), which was indexed by subscript $i$. Suppose a subset of
$\mathcal{D}$ comprising sequences with label $l$, that is,

(45) $\mathcal{D}_{l}=\\{(X_{i},y_{i})\mid i=1,\dots,n,\,y_{i}=l\\},$

and let $n_{l}=|\mathcal{D}_{l}|$ denote the number of sequences with label
$l$. Now suppose to give an arbitrary order to the elements of
$\mathcal{D}_{l}$. We denote the _original index_ (in $\mathcal{D}$) of the
$p$-th element of $\mathcal{D}_{l}$ by $i_{l,p}$. In other words, $i_{l,p}$
denotes the index of the $p$-th (in an arbitrary order) sequence with label
$l$. Given such definitions, we denote the value of the kernel function
$k_{\text{DMS}}$ in (26) between $(\Theta_{i_{l,p}},\bm{X}_{i_{l,p}})$ and
$(\Theta_{i_{m,q}},\bm{X}_{i_{m,q}})$ by $K_{l,m}^{p,q}$, that is,

(46)
$K_{l,m}^{p,q}=k_{\text{DMS}}(\mathcal{W}_{\Theta_{i_{l,p}},\bm{X}_{i_{l,p}}},\mathcal{W}_{\Theta_{i_{m,q}},\bm{X}_{i_{m,q}}}).$

Moreover, we let $\bm{K}_{l,m}$ be the kernel matrix whose $(p,q)$-element is
$K_{l,m}^{p,q}$.

Given the above quantities, we can compute the quantities related to
$\bm{\Sigma}^{\phi}$ as follows:

(47)
$\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{m}^{\phi})=\operatorname{tr}(\underbrace{\bm{K}_{l,m}(\bm{I}-\bm{O}_{n_{m}})}_{\tilde{\bm{K}}_{l,m}}\underbrace{\bm{K}_{m,l}(\bm{I}-\bm{O}_{n_{l}})}_{\tilde{\bm{K}}_{m,l}})=\sum_{p=1}^{n_{l}}\sum_{q=1}^{n_{m}}\tilde{K}_{l,m}^{p,q}\tilde{K}_{m,l}^{q,p},$

and

(48)
$\frac{\partial\operatorname{tr}(\bm{\Sigma}_{l}^{\phi}\bm{\Sigma}_{m}^{\phi})}{\partial\theta_{i,j}}=\sum_{p}\sum_{q}\Bigg{\\{}\Bigg{(}\frac{\partial
K_{l,m}^{p,q}}{\partial\theta_{i,j}}-\frac{1}{n_{m}}\sum_{q}\frac{\partial
K_{l,m}^{p,q}}{\partial\theta_{i,j}}\Bigg{)}\tilde{K}_{m,l}^{q,p}\\\
\qquad\qquad+\tilde{K}_{l,m}^{p,q}\Bigg{(}\frac{\partial
K_{m,l}^{q,p}}{\partial\theta_{i,j}}-\frac{1}{n_{l}}\sum_{p}\frac{\partial
K_{m,l}^{q,p}}{\partial\theta_{i,j}}\Bigg{)}\Bigg{\\}},$

where $\bm{O}_{n_{l}}$ denotes an $n_{l}\times n_{l}$ matrix with all elements
being $1/n_{l}$.

The remaining part of (43) is computed by

(49) $\frac{\partial
Q_{2}}{\partial\theta_{i,j}}=\sum_{l=1}^{c-1}\sum_{m=l+1}^{c}\frac{n_{l}n_{m}}{n^{2}}\Bigg{(}\frac{1}{n_{l}^{2}}\sum_{q}\sum_{q^{\prime}}\frac{\partial
K_{l,l}^{q,q^{\prime}}}{\partial\theta_{i,j}}\\\
-\frac{2}{n_{l}n_{m}}\sum_{q}\sum_{p}\frac{\partial
K_{m,l}^{q,p}}{\partial\theta_{i,j}}+\frac{1}{n_{m}^{2}}\sum_{p}\sum_{p^{\prime}}\frac{\partial
K_{m,m}^{p,p^{\prime}}}{\partial\theta_{i,j}}\Bigg{)}.$

The final missing piece is the derivative of kernel, $\frac{\partial
K_{l,m}^{p,q}}{\theta_{i,j}}$, which is computed as

(50) $\frac{\partial
K_{i,k}^{p,q}}{\partial\theta_{i,j}}=\operatorname{sum}\left(\frac{\partial
K_{i,k}^{p,q}}{\partial\bm{V}_{s}}\circ\frac{\partial\bm{V}_{s}}{\partial\theta_{i,j}}\right).$

As the second term is already given in (42), we focus on the first term. Let
$\bm{A}_{i}\in\mathbb{C}^{p\times r}$ and $\bm{B}_{i}\in\mathbb{C}^{p\times
r}$ denote the matrices whose columns comprise left- and right-singular
vectors of $\bm{X}_{i}\bm{V}_{\Theta_{i}}^{\dagger}$. Let $\bm{S}_{i}$ be a
diagonal matrix comprising the singular values of
$\bm{X}_{i}\bm{V}_{\Theta_{i}}^{\dagger}$ (in the same order with $\bm{A}_{i}$
and $\bm{B}_{i}$). Then,

(51) $\frac{\partial
K_{l,m}^{p,q}}{\partial\bm{V}_{\Theta_{i}}}=\Big{(}\bm{V}_{\Theta_{i}}^{\dagger}(\bm{V}_{\Theta_{i}}^{\dagger})^{\mathsf{H}}(\bm{C}_{l,m}^{p,q})^{\mathsf{H}}(\bm{I}-\bm{V}_{\Theta_{i}}\bm{V}_{\Theta_{i}}^{\dagger})^{\mathsf{H}}\\\
+(\bm{I}-\bm{V}_{\Theta_{i}}^{\dagger}\bm{V}_{\Theta_{i}})^{\mathsf{H}}(\bm{C}_{l,m}^{p,q})^{\mathsf{H}}(\bm{V}_{\Theta_{i}}^{\dagger})^{\mathsf{H}}\bm{V}_{\Theta_{i}}^{\dagger}-\bm{V}_{\Theta_{i}}^{\dagger}\bm{C}_{l,m}^{p,q}\bm{V}_{\Theta_{i}}^{\dagger}\Big{)}^{\mathsf{T}},$

where

(52)
$\bm{C}_{l,m}^{p,q}=(\bm{B}_{i_{l,p}}\bm{S}_{i_{l,p}}^{-1}\bm{A}_{i_{l,p}}^{\mathsf{H}})\bm{A}_{i_{m,q}}\bm{A}_{i_{m,q}}^{\mathsf{H}}(\bm{I}-\bm{A}_{i_{l,p}}\bm{A}_{i_{l,p}}^{\mathsf{H}})\bm{Y}_{i_{l,p}}.$

## References

  * [1] I. Abraham, G. D. L. Torre, and T. D. Murphey, Model-based control using Koopman operators, in Robotics: Science and Systems 2017 Proceedings, 2017.
  * [2] H. Arbabi and I. Mezić, Ergodic theory, dynamic mode decomposition and computation of spectral properties of the Koopman operator, SIAM Journal on Applied Dynamical Systems, 16 (2017), pp. 2096–2126.
  * [3] T. Askham and J. N. Kutz, Variable projection methods for an optimized dynamic mode decomposition, SIAM Journal on Applied Dynamical Systems, 17 (2018), pp. 380–416.
  * [4] E. Barshan, A. Ghodsi, Z. Azimifar, and M. Z. Jahromi, Supervised principal component analysis: Visualization, classification and regression on subspaces and submanifolds, Pattern Recognition, 44 (2011), pp. 1357–1371.
  * [5] B. Blankertz, R. Tomioka, S. Lemm, M. Kawanabe, and K.-R. Muller, Optimizing spatial filters for robust EEG single-trial analysis, IEEE Signal Processing Magazine, 25 (2008), pp. 41–56.
  * [6] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, Extracting spatial-temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition, Journal of Neuroscience Methods, 258 (2016), pp. 1–15.
  * [7] M. Budišić, R. Mohr, and I. Mezić, Applied Koopmanism, Chaos, 22 (2012), 047510, p. 047510.
  * [8] L. M. Candanedo, V. Feldheim, and D. Deramaix, Data driven prediction models of energy use of appliances in a low-energy house, Energy and Buildings, 140 (2017), pp. 81–97.
  * [9] K. K. Chen, J. H. Tu, and C. W. Rowley, Variants of dynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses, Journal of Nonlinear Science, 22 (2012), pp. 887–915.
  * [10] M. Churchland, J. Cunningham, M. Kaufman, J. Foster, P. Nuyujukian, S. Ryu, and K. Shenoy, Neural population dynamics during reaching, Nature, 487 (2012), pp. 51–56.
  * [11] A. Cichocki, R. Zdunek, A. H. Phan, and S. Amari, Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation, John Wiley & Sons, 2009.
  * [12] S. Das and D. Giannakis, Delay-coordinate maps and the spectra of Koopman operators, Journal of Statistical Physics, 175 (2019), pp. 1107–1145.
  * [13] R. A. Fisher, The use of multiple measurements in taxonomic problems, Annals of Eugenics, 7 (1936), pp. 179–188.
  * [14] K. Fujii and Y. Kawahara, Supervised dynamic mode decomposition via multitask learning, Pattern Recognition Letters, 122 (2019), pp. 7–13.
  * [15] K. Fujii, N. Takeishi, B. Kibushi, M. Kouzaki, and Y. Kawahara, Data-driven spectral analysis for coordinative structures in periodic human locomotion, Scientific Reports, 9 (2019), p. 16755.
  * [16] K. Fukumizu, F. R. Bach, and M. I. Jordan, Kernel dimension reduction in regression, The Annals of Statistics, 37 (2009), pp. 1871–1905.
  * [17] T. Funato, S. Aoi, H. Oshima, and K. Tsuchiya, Variant and invariant patterns embedded in human locomotion through whole body kinematic coordination, Experimental Brain Research, 205 (2010), pp. 497–511.
  * [18] M. Georgescu, B. Eisenhower, and I. Mezić, Creating zoning approximations to building energy models using the Koopman operator, in IBPSA-USA SimBuild 2012, 2012, pp. 40–47.
  * [19] D. Giannakis, Data-driven spectral decomposition and forecasting of ergodic dynamical systems, Applied and Computational Harmonic Analysis, 42 (2019), pp. 338–396.
  * [20] J. Goldberger, G. E. Hinton, S. T. Roweis, and R. R. Salakhutdinov, Neighbourhood components analysis, in Advances in Neural Information Processing Systems, vol. 17, 2005, pp. 513–520.
  * [21] G. H. Golub and V. Pereyra, The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate, SIAM Journal on Numerical Analysis, 10 (1973), pp. 413–432.
  * [22] J. Hamm and D. D. Lee, Grassmann discriminant analysis: a unifying view on subspace-based learning, in Proceedings of the 25th International Conference on Machine Learning, 2008, pp. 376–383.
  * [23] I. T. Jolliffe, Principal Component Analysis, Springer, 2nd ed., 2002\.
  * [24] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, Sparsity-promoting dynamic mode decomposition, Physics of Fluids, 26 (2014), p. 024103.
  * [25] M. Korda and I. Mezić, On convergence of extended dynamic mode decomposition to the Koopman operator, Journal of Nonlinear Science, 28 (2018), pp. 687–710.
  * [26] Y. Koren, R. Bell, and C. Volinsky, Matrix factorization techniques for recommender systems, Computer, 42 (2009), pp. 30–37.
  * [27] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems, SIAM, 2016.
  * [28] D. J. McFarland, L. M. McCane, S. V. David, and J. R. Wolpaw, Spatial filter selection for EEG-based communication, Electroencephalography and Clinical Neurophysiology, 103 (1997), pp. 386–394.
  * [29] I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics, 41 (2005), pp. 309–325.
  * [30] S. Mika, G. Rätsch, J. Weston, B. Schölkopf, and K.-R. Müller, Fisher discriminant analysis with kernels, in Proceedings of the 9th IEEE Workshop on Neural Networks for Signal Processing, 1999, pp. 41–48.
  * [31] J. L. Proctor and P. A. Eckhoff, Discovering dynamic patterns from infectious disease data using dynamic mode decomposition, International Health, 7 (2015), pp. 139–145.
  * [32] H. Ramoser, J. Müller-Gerking, and G. Pfurtscheller, Optimal spatial filtering of single trial EEG during imagined hand movement, IEEE Transactions on Rehabilitation Engineering, 8 (2000), pp. 441–446.
  * [33] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, Spectral analysis of nonlinear flows, Journal of Fluid Mechanics, 641 (2009), pp. 115–127.
  * [34] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics, 656 (2010), pp. 5–28.
  * [35] B. Su, X. Ding, C. Liu, H. Wang, and Y. Wu, Discriminative transformation for multi-dimensional temporal sequences, IEEE Transactions on Image Processing, 26 (2017), pp. 3579–3593.
  * [36] B. Su, X. Ding, H. Wang, and Y. Wu, Discriminative dimensionality reduction for multi-dimensional sequences, IEEE Transactions on Pattern Analysis and Machine Intelligence, 40 (2018), pp. 77–91.
  * [37] B. Su and Y. Wu, Learning low-dimensional temporal representations, in Proceedings of the 35th International Conference on Machine Learning, 2018, pp. 4761–4770.
  * [38] Y. Susuki and I. Mezić, Nonlinear Koopman modes and power system stability assessment without models, IEEE Transactions on Power Systems, 29 (2014), pp. 899–907.
  * [39] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, On dynamic mode decomposition: Theory and applications, Journal of Computational Dynamics, 1 (2014), pp. 391–421.
  * [40] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science, 25 (2015), pp. 1307–1346.
  * [41] S. Wold, M. Sjöström, and L. Eriksson, PLS-regression: a basic tool of chemometrics, Chemometrics and Intelligent Laboratory Systems, 58 (2001), pp. 109–130.
  * [42] D. You, O. C. Hamsici, and A. M. Martinez, Kernel optimization in discriminant analysis, IEEE Transactions on Pattern Analysis and Machine Intelligence, 33 (2011), pp. 631–638.