# A data driven Koopman-Schur decomposition for computational analysis of
nonlinear dynamics††thanks:  This material is based upon work supported by the
DARPA Small Business Innovation Research Program (SBIR) Program Office under
Contract No. W31P4Q-21-C-0007 to AIMdyn, Inc. Any opinions, findings and
conclusions or recommendations expressed in this material are those of the
authors and do not necessarily reflect the views of the DARPA SBIR Program
Office. Distribution Statement A: Approved for Public Release, Distribution
Unlimited.

Zlatko Drmač University of Zagreb, Faculty of Science, Department of
Mathematics, Bijenička cesta 30, 10000 Zagreb, Croatia. The author is also
supported by the Croatian Science Foundation (CSF) grant IP-2019-04-6268
_“Randomized low rank algorithms and applications to parameter dependent
problems”_. Igor Mezić University of California at Santa Barbara and AIMdyn,
Inc., 1919 State St., Ste. 207, Santa Barbara, CA 93101

###### Abstract

This paper introduces a new theoretical and computational framework for a data
driven Koopman mode analysis of nonlinear dynamics. To alleviate the potential
problem of ill-conditioned eigenvectors in the existing implementations of the
Dynamic Mode Decomposition (DMD) and the Extended Dynamic Mode Decomposition
(EDMD), the new method introduces a Koopman-Schur decomposition that is
entirely based on unitary transformations. The analysis in terms of the
eigenvectors as modes of a Koopman operator compression is replaced with a
modal decomposition in terms of a flag of invariant subspaces that correspond
to selected eigenvalues. The main computational tool from the numerical linear
algebra is the partial ordered Schur decomposition that provides convenient
orthonormal bases for these subspaces. In the case of real data, a real Schur
form is used and the computation is based on real orthogonal transformations.
The new computational scheme is presented in the framework of the Extended DMD
and the kernel trick is used.

##  1 Introduction and preliminaries

High resolution measurements (using e.g. particle image velocimetry (PIV),
high speed cameras, magnetometers, microwave radiometers) or high fidelity
numerical simulations of the dynamics under study (given e.g. as a system of
ordinary differential equations obtained in a semi-discretization of partial
differential equations) often provide an abundance of data. After
preprocessing and formatting, the data is represented as a discrete time
sequence of data snapshots $\mathbf{x}_{1},\mathbf{x}_{2},\ldots$ in the state
space that is set to be $\mathbb{C}^{n}$ or $\mathbb{R}^{n}$, and assumed to
be driven by a mapping $\mathbf{x}_{i+1}=\mathbf{F}(\mathbf{x}_{i})$, for some
initial state $\mathbf{x}_{1}$ and with given time step $\delta t$. The
mapping $\mathbf{F}$ may not be known; in a numerical simulation it represents
software implementation of a numerical algorithm. The dimension $n$ ranges,
depending on the application, from moderate to extremely large. For instance,
if the sequence $\mathbf{x}_{1},\mathbf{x}_{2},\ldots$ is obtained by
discretization with time step $\delta t$ of a time dependent partial
differential equation, then $\mathbf{x}_{i}$ is represents the solution on a
spatially discretized domain at the discrete time moment $i$ and the dimension
$n$ (the number of grid points in the discretized domain) can be in hundreds
of thousands or in millions. In such cases the natural structure of the data
is a fourth order tensor (three space dimensions times the time), and the
ambient space $\mathbb{C}^{n}$ or $\mathbb{R}^{n}$ is given implicitly by a
reshaping isomorphism (vectorization).

##### The Dynamic Mode Decomposition (DMD).

Revealing relevant features, latent in the data but valuable for understanding
and control of the underlying processes, requires suitable data mining
techniques. The Dynamic Mode Decomposition (DMD) [45] is such a tool of the
trade. It has been initially motivated by the problems in the computational
fluid dynamics, but its success ranges across a wide spectrum of applications,
such as combustion, [24], robotics [6], power grid [49] and automatic control
[36].

The tacit assumption in the DMD is that the given sequence of snapshots
$\mathbf{x}_{1},\ldots,\mathbf{x}_{m},\mathbf{x}_{m+1}$ is generated by a
linear operator (matrix) $A_{n}$ such that, given $\mathbf{x}_{1}$,
$\mathbf{x}_{i+1}=A_{n}\mathbf{x}_{i}$ for $i=1,\ldots,m$. If we set
$\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{m})$,
$\mathbf{Y}=(\mathbf{x}_{2},\ldots,\mathbf{x}_{m+1})$, then
$\mathbf{Y}=A_{n}\mathbf{X}$. Since the action of $A_{n}$ on $\mathbf{X}$ is
given, we can use the Rayleigh-Ritz procedure to extract from the range of
$\mathbf{X}$ approximate eigenpairs of $A_{n}$. To that end, a POD basis is
computed using the SVD of $\mathbf{X}$, $\mathbf{X}=U_{r}\Sigma_{r}V_{r}^{*}$,
where $r$ is the rank of $\mathbf{X}$. In practice, this will be the numerical
rank, i.e. small singular values of $\mathbf{X}$ will be truncated. Using the
orthonormal basis $U_{r}$, the Rayleigh quotient is computed as
$A=U_{r}^{*}A_{n}U_{r}=U_{r}^{*}\mathbf{Y}V_{r}\Sigma_{r}^{-1}$. If we
diagonalize $A$, $A=G\Lambda G^{-1}$ with
$\Lambda=\mathrm{diag}(\lambda_{i})_{i=1}^{r}$, then the Ritz pairs of $A_{n}$
are $(U_{r}G,\Lambda)$. If we set
$K=U_{r}G=(\mathbf{k}_{1},\ldots,\mathbf{k}_{r})$, then
$A_{n}\mathbf{k}_{i}\approx\lambda_{i}\mathbf{k}_{i}$. Note that there is no
need to form the matrix $A_{n}$ explicitly. A software implementation of the
DMD has been recently included in the LAPACK library [2], [16], [17].

Now, revealing the structure of the snapshots (e.g. discovering coherent
structures in the fluid flow) is based on a spatio-temporal representation of
the snapshots using selected $\ell\leq r$ modes
$\mathbf{k}_{j_{1}},\ldots,\mathbf{k}_{j_{\ell}}$: For given weights
$\omega_{i}\geq 0$, the task is to determine coefficients
$\alpha_{1},\ldots,\alpha_{\ell}$ so that

$\mathbf{x}_{i}\approx\sum_{k=1}^{\ell}\mathbf{k}_{j_{k}}\alpha_{k}\lambda_{j_{k}}^{i-1}\;\;\mbox{and}\;\;\sum_{i=1}^{m+1}\omega_{i}^{2}\|\mathbf{x}_{i}-\sum_{k=1}^{\ell}\mathbf{k}_{j_{k}}\alpha_{k}\lambda_{j_{k}}^{i-1}\|^{2}\longrightarrow\min_{\alpha_{k}}.$
(1)

Note that in the above formula, the representation of the snapshot
$\mathbf{x}_{i+1}=A_{n}\mathbf{x}_{i}$ is obtained simply by raising by one
the powers $\lambda_{j_{k}}^{i-1}$ in the representation of $\mathbf{x}_{i}$.
The key is that the coefficients in the representation (1) of $\mathbf{x}_{i}$
also contain the time stamp. This is a basis for forecasting – extrapolation
into the future by raising the powers of the $\lambda_{j_{k}}$’s in the
representation (1) of the last available (i.e. the present) snapshot. The
representation of the snapshots (1) in terms of selected modes
$\mathbf{k}_{j}$ is most dynamics revealing if it succeeds with small $\ell$,
and such a sparse representation can be obtained using a sparsity promoting
DMD [31].

For more details on the variations of the DMD and its various applications see
e.g. [44], [8], [26], [12], [27], [53], [50], [52], [51], [42], [4], [35].

##### The Koopman connection and the Extended DMD (EDMD).

The theoretical underpinning of the DMD is in close connection to the Koopman
(composition) operator $\mathbb{U}$ associated with the discrete dynamical
system $\mathbf{x}_{i+1}=\mathbf{F}(\mathbf{x}_{i})$. In fact, the leftmost
formula in (1) is a version of the general Koopman Mode Decomposition, first
derived in [37]. $\mathbb{U}$ is defined on a class of functions $\mathcal{F}$
equipped with a Hilbert space structure (typically $L^{2}$ for measure-
preserving systems, see [38] for the Hilbert space construction in the
dissipative case) as $\mathbb{U}f=f\circ\mathbf{F}$, $f\in\mathcal{F}$. Here
$\circ$ denotes the composition of mappings. The functions $f\in\mathcal{F}$
are called the observables of the system. Note that
$(\mathbb{U}f)(\mathbf{x}_{i})=f(\mathbf{F}(\mathbf{x}_{i}))=f(\mathbf{x}_{i+1})$,
so that $\mathbb{U}$ acts similarly as the matrix $A_{n}$ in the DMD, but in a
more general setting. The introduction of $\mathbb{U}$ is actually a
linearization of $\mathbf{F}$, albeit an infinitely dimensional one. The
connection with the DMD is explicitly stated in the Extended DMD [57], which
encompasses the DMD as a special case, and can be understood as a Galerkin
type approximation of the eigenvalues and eigenfunctions of the Koopman
operator, based on its compression to an $N$ dimensional subspace of
$\mathcal{F}$. For an excellent review see e.g. [7].

It is known that, under certain assumptions, when the number of snapshots in
the limit $m\rightarrow\infty$ and the dimension $N\rightarrow\infty$, the
EDMD in several ways (including the convergence of the eigenfunctions and the
eigenvalues along subsequences) converges to the operator $\mathbb{U}$, see
[32]. The use of Hankel-DMD [3] under the assumption of existence of finite-
dimensional invariant subspace leads to stronger results, that were extended
to convergence in pseudospectral sense in [40], without the assumption on
existence of a finite-dimensional invariant subspaces, under the assumption of
discrete spectrum of the Koopman operator.

Here too, the eigenfunctions of $\mathbb{U}$ are numerically approximated
using a discrete compression of $\mathbb{U}$ to an $N$-dimensional subspace
$\mathcal{F}_{N}=\mathrm{span}({\boldsymbol{\psi}}_{1},\ldots,{\boldsymbol{\psi}}_{N})\subset\mathcal{F}$,
where $N\geq n$ and potentially $N\gg n$. Similarly as in the DMD, a Rayleigh
quotient matrix $U_{N}$ is diagonalized and its eigenvectors are used to
construct discrete approximations of certain eigenfunctions of $\mathbb{U}$.
In a special case, the matrix $U_{N}$ and the DMD matrix $A_{n}$ are
transposes of each other,111See Section 3. and the eigenvectors of $A_{n}$ are
the modes in the data snapshot representation (1).

##### A quandary about the eigenvectors.

There is one important detail that is, to the best of our knowledge, not
properly addressed in the literature. If in an application we compute a
spectral decomposition $A=G\Lambda G^{-1}$ of a highly non-normal matrix $A$,
assuming $A$ is diagonalizable, then the matrix $G$ of the eigenvectors can be
severely ill-conditioned and in that case the matrix of the modes $K=U_{r}G$
must also be ill-conditioned. Ill-conditioned means close to singularity and
numerical computation (such as representing the snapshots as linear
combinations of the eigenvectors (1)) becomes prone to large errors.

Even if $A$ is not diagonalizable (i.e. it is defective), a numerical software
will compute spectral information that corresponds to a nearby $A+\delta A$
(this is established by a backward or a mixed error analysis of the algorithm
that is implemented in the software) that is most likely diagonalizable.
(Recall that non-diagonalizable matrices are a nowhere dense subset - the
diagonalizable ones are an open dense set.) As a result, a numerical software
is likely to return a full set of eigenvectors, but alas, not of $A$. A simple
example illustrates this. Let $A=\left(\begin{smallmatrix}1&1\cr
0&1\end{smallmatrix}\right)$, $A+\delta
A=\left(\begin{smallmatrix}1&1\cr\epsilon&1\end{smallmatrix}\right)$, where
$\epsilon>0$ is a small parameter. $A$ is a Jordan block with double
eigenvalue $\lambda=1$ of geometric multiplicity one. The eigenvalues of
$A+\delta A$ are $\lambda_{1,2}=1\pm\sqrt{\epsilon}$ and the spectral
decomposition of $\widetilde{A}=A+\delta A$ reads

$\widetilde{A}=\widetilde{G}\widetilde{\Lambda}\widetilde{G}^{-1},\;\widetilde{G}=\begin{pmatrix}1&1\cr\sqrt{\epsilon}&-\sqrt{\epsilon}\end{pmatrix}\\!,\;\widetilde{\Lambda}=\begin{pmatrix}1+\sqrt{\epsilon}&0\cr
0&1-\sqrt{\epsilon}\end{pmatrix}\\!,\;\widetilde{G}^{-1}=\frac{1}{2}\begin{pmatrix}1&1/\sqrt{\epsilon}\cr
1&-1/\sqrt{\epsilon}\end{pmatrix}.$

Note that $\|\widetilde{G}\|_{2}\|\widetilde{G}^{-1}\|_{2}>1/\sqrt{\epsilon}$
and that no scaling of the eigenvectors can improve this ill-conditioning with
a factor better than $\sqrt{2}$. This is a consequence of a van der Sluis
theorem [56] on optimally scaled matrices and the fact that both columns of
$\widetilde{G}$ are of Euclidean norm $\sqrt{1+\epsilon}\approx 1$.

Hence, theoretical considerations that are based on the nontrivial (i.e. non-
diagonal) Jordan normal form and the generalized eigenvectors are not robust
for practical computation – an arbitrarily small perturbation can change the
Jordan structure entirely and e.g. make the previously identified generalized
eigenvectors to simply disappear and replace them with potentially ill-
conditioned eigenvectors. For a more general result on optimal scaling of the
eigenvectors in the case of multiple eigenvalues and on the Jordan normal form
from a numerical analyst’s perspective we recommend [13].

The analysis of the case when the Koopman operator itself has a degeneracy
that leads to generalized eigenfunctions was presented in [38].222 The notion
of generalized eigenfunctions used in [38] is in the same sense as generalized
eigenvectors of a matrix - i.e.functions that span an invariant subspace of
the Koopman operator of dimension bigger than $1$ that contains only $1$
eigenfunction. This is in contrast with the notion of generalized
eigenfunctions in functional analysis [23] that are objects in Koopman theory
associated with continuous spectrum.  A simple example in which generalized
eigenfunctions arise is the linear system

$\dot{\mathbf{y}}=A\mathbf{y},$ (2)

where $\mathbf{y}=(y_{1},y_{2})^{T}$, $A=\left(\begin{smallmatrix}1&1\cr
0&1\end{smallmatrix}\right)$. Since linear (generalized) eigenfunctions of the
Koopman operator in the linear case are given by
$\left<\mathbf{y},\mathbf{w}\right>,$ where $\mathbf{w}$ is a left
(generalized) eigenvector of $A$ and $\langle\cdot,\cdot\rangle$ is the
complex inner product defined by
$\langle\mathbf{x},\mathbf{y}\rangle=\sum_{i}\mathbf{x}_{i}\cdot\overline{\mathbf{y}_{i}}.$
The eigenvector $\mathbf{w}_{1}$ associated with eigenvalue $1$ is
$\mathbf{w}_{1}=(0,1)^{T}$, and the generalized eigenvector can be chosen to
be $(1,0)^{T}$. therefore $y_{2}$ is an eigenfunction, and $y_{1}$ a
generalized eigenfunction of the system.

The paper [46] contains examples of non-normal operators in fluid flows, that
has substantial impact on treatment of stability of such flows. The most well-
known example of such non-normality is the formulation of the viscous
stability problem for parallel shear flows, in which linearization of the
Navier-Stokes equation leads to the Orr-Sommerfeld linear partial differential
equation whose solutions exhibit non-normal behavior.

##### An alternative proposed in this paper.

We explore a different method to numerically extract spectral information
about $\mathbb{U}$ and to deploy it in applications with the same
functionality as provided by the DMD and the EDMD. The key tool in the new
approach is the Schur triangular form, which is in the framework of the
numerical linear algebra a viable alternative to the diagonalization
$A=G\Lambda G^{-1}$ or computation of the Jordan canonical form.

###### Theorem 1.1

(Schur form) Let $A\in\mathbb{C}^{r\times r}$ and let
$\lambda_{1},\ldots,\lambda_{r}$ be its eigenvalues, listed in an arbitrary
order. Then there exist a unitary $Q$ and an upper triangular $T$ such that
$T_{ii}=\lambda_{i}$, $i=1,\ldots,r$, and $A=QTQ^{*}.$

###### Remark 1.2

If we write $T=\mathrm{diag}(\lambda_{i})+N$, then
$A=Q\mathrm{diag}(\lambda_{i})Q^{*}+QNQ^{*}$ is the sum of a normal and a
nilpotent matrix; $(QNQ^{*})^{r}=\mathbf{0}$. If we let $T_{1:i,1:i}$ denote
the leading $i\times i$ submatrix of $T$, and $Q_{:,1:i}$ the submatrix of the
leading $i$ columns of $Q$, then $AQ_{:,1:i}=Q_{:,1:i}T_{1:i,1:i}$, i.e. the
matrices $Q_{:,1:i}$ span a flag of invariant subspaces of $A$, with the
corresponding compressions (Rayleigh quotients) given by
$T_{1:i,1:i}=Q_{:,1:i}^{*}AQ_{:,i}$. The corresponding orthogonal projectors
$P_{i}=Q_{:,1:i}Q_{:,1:i}^{*}$ are increasing in the Loewner partial order:
$P_{i+1}\succeq P_{i}$, meaning that $P_{i+1}-P_{i}=Q_{:,i}Q_{:,i}^{*}$ is
positive semidefinite.  $\boxtimes$

###### Remark 1.3

If $A=QTQ^{*}$, then $A^{*}=QT^{*}Q^{*}=(Q\Pi)(\Pi^{T}T^{*}\Pi)(Q\Pi)^{*}$ is
the Schur form of $A^{*}$, where $\Pi$ is the permutation matrix that reverses
the order of columns, so that $\Pi^{T}T^{*}\Pi$ is upper triangular.
Similarly, the Schur form of the transposed matrix reads
$A^{T}=(Q\Pi)^{*T}(\Pi^{T}T^{T}\Pi)(Q\Pi)^{T}$.  $\boxtimes$

An advantage of this decomposition is that it is numerically computed using
only unitary transformations, and that it does not refer to the Jordan normal
form of $A$. Having $A=QTQ^{*}$, we have in the DMD
$A_{n}(U_{r}Q)\approx(U_{r}Q)T$, and in the EDMD this yields a partial Schur
decomposition of the matrix compression of $\mathbb{U}$.

When using the Schur form, in the general non-normal case we do not have
direct access to the individual eigenvectors333The exception is the first
column of $Q$, $AQ(:,1)=T_{11}Q(:,1)$. (discretized eigenfunctions) anymore –
instead, we have a flag of nearly invariant subspaces and additional work is
needed to get the same functionality as in the EDMD. Our proposed solution is
based on the numerical linear algebra techniques, in particular on the ordered
Schur form. A special case that covers Hermitian and skew-Hermitian DMD (thus
diagonal Schur form) is recently published in [17]. We have already used an
ordered partial Schur decomposition444We called it partial Schur Modal
Decomposition (PSMD) in [15, Section 5.3] to improve the Optimal Mode
Decomposition [58], and we discussed some details of using a partial Schur
decomposition in a DMD computational framework. 555 A similar idea is also
used recently in [54], but only partially. The authors of [54] seemed unaware
of our 2018 paper [15].

For the behaviour of the computed matrix Schur forms in the limit for
$m\rightarrow\infty$ and $N\rightarrow\infty$, we must consider a Schur form
of $\mathbb{U}$. Since the matrix Schur form can be extended to spectral
operators [43], we proceed beyond the numerical procedure and complete the new
proposed framework with the Schur decomposition of the Koopman operator and
the convergence of the numerical method.

##### Outline of the paper.

The rest of the paper is organized as follows. In Section 2 we set the
operator theoretic stage for the Koopman operator and its Schur form – from
the definition of the Koopman semigroup, nonlinear dynamics, to properly
defined function spaces.

In Section 3, we review the EDMD in detail. The Koopman operator based
interpretation is stressed (as opposed to the mere regression) and the details
of a data driven compression of the operator and its Rayleigh quotient are
given. Further, the kernel trick is explained and elements of convergence in
appropriate Hilbert spaces are given. Together with details such as evaluation
of the approximate eigenfunctions on the snapshots, this section serves to
define the framework and set the stage for our new proposed Koopman-Schur
Subspace Modal Decomposition (KS-SSMD), that is described in detail in Section
4. We show that the KS-SSMD has all the functionalities of the EDMD, and, in
particular, that it can be used to reveal coherent structures (Section 4.4)
and for forecasting (Section 4.5). The flexibility of working with only a
subset of selected eigenvalues is achieved using the ordered Schur form, that
is outlined in Section 4.2. The new approach is put to test in Section 5. We
use numerical examples to illustrate subtle numerical issues and the potential
of the KS-SSMD, in particular good forecasting skills. This report concludes
with Section 6, that contains a review of our current work on the further
development of the method.

##  2 Koopman operator and its Schur form

Let $T:M\rightarrow M$ be a dynamical system on a set $M$. Let the Koopman
operator $\mathbb{U}$ associated with $T$, be defined by

$\mathbb{U}f=f\circ T,f\in\cal F$

where $\cal F$ is a space of observables closed under the action of
$\mathbb{U}$ i.e. $\mathbb{U}(\cal F)\subset\cal F.$

We recall the definition of a spectral operator [22]:

###### Definition 2.1

A bounded operator $S$ on a space $\cal F$ is of scalar type if it satisfies

$S=\int_{\mathbb{C}}\lambda dP_{\lambda},$ (3)

where $dP$ is the resolution of the identity for $S$. An operator $U$ is
spectral, if it can be decomposed as

$U=S+Q_{0},$

where $S$ is of scalar type, and $Q_{0}$ is quasi-nilpotent, i.e.

$\lim_{n\rightarrow\infty}(Q_{0}^{n})^{1/n}=0.$

The following theorem is due to Rogers [43].

###### Theorem 2.1

Let U be a spectral operator. Then $U=N+Q_{1}$ where $N$ is normal, $Q_{1}$ is
quasinilpotent, and $\sigma(N)=\sigma(U)$. Also $Q_{1}N$, $NQ_{1}$, and
$NQ_{1}-Q_{1}N$ are quasinilpotent. Furthermore there is a countable nest of
selfadjoint projections ${\cal N}=\\{P_{\alpha}\\}_{\alpha\in\Gamma}$ with the
following properties:

  1. 1.

Each $P_{\alpha}$, reduces $N$ and is invariant under $Q$,

  2. 2.

The operator $N$ is in the $C^{*}$-algebra generated by ${\cal N}$,

  3. 3.

The nest ${\cal N}$ separates $\sigma(U)$ in the sense that given any
$\lambda,\mu\in\sigma(U)$, there exists a $P_{\alpha}$, such that
$\sigma(U_{P_{\alpha}{\cal F}})$ contains exactly one of $\lambda,\mu$.

We work in continuous time $t$, where ${\bf S}^{t}:M\rightarrow M$ is a group
of transformations of $M$. The type system we consider here is the one that
has a global (Milnor) attractor $\mathcal{A}$ of zero Lebesgue measure, where
$\mathcal{A}$ is equipped with a physical measure [9]. The operator $U$ can be
obtained from the family of operators $U^{t}$ defined by $U^{t}f=f\circ{\bf
S}^{t}$ by sampling the evolution of the system at discrete times $n\Delta
t,n\in\mathbb{Z}$, with $T=S^{\Delta t}$. Let ${\cal{D}}$ be a compact forward
invariant set of ${\bf S}^{t}$, and $\mu$ the physical invariant measure on
the (Milnor) attractor $\mathcal{A}\subset{\cal{D}}$, with
$\mu(\mathcal{A})=1$. Let $\mathcal{H}_{\mathcal{A}}$ be the space of square-
integrable observables on $\mathcal{A}$. Note that $U^{t}$ restricted to
$\mathcal{H}_{\mathcal{A}}$ is unitary. Thus the theorem above yields
$U^{t}|_{\mathcal{H}_{\mathcal{A}}}=N^{t}$, i.e. $U^{t}$ is a normal operator
when restricted to $\mathcal{H}_{\mathcal{A}}$. The functions in
$\mathcal{H}_{\mathcal{A}}$ can typically be thought of as being defined on
the whole basin of attraction of $\mathcal{A}$. Let
$\tilde{\mathcal{H}}_{\cal{B}}\subset C({\cal{D}})$ be a Hilbert space of
continuous functions $f:{\cal{D}}\rightarrow\mathbb{C}$ orthogonal to
$\mathcal{H}_{\mathcal{A}}$ with respect to $\mu$, that satisfy

$\int_{\cal{D}}f\phi d\mu=0,\ \forall\phi\in\mathcal{H}_{\mathcal{A}}.$ (4)

We assume that the tensor product space
$\tilde{\mathcal{H}}=\mathcal{H}_{\mathcal{A}}\otimes\tilde{\mathcal{H}}_{\cal{B}}$
is invariant under $U^{t}$. It is clear that the space
$\mathcal{H}_{\mathcal{A}}$ is itself invariant under $U^{t}$. So is
$\tilde{\mathcal{H}}_{\cal{B}}$ [39]:

###### Lemma 2.1

$\tilde{\mathcal{H}}_{\cal{B}}$ is invariant under $U^{t}$.

Define $\mathcal{H}_{\cal{B}}=\tilde{\mathcal{H}}_{\cal{B}}\cup{\mathbf{1}}$,
where ${\mathbf{1}}$ is the constant unit function on ${\cal{D}}$. Define the
tensor product space

$\mathcal{H}=\mathcal{H}_{\mathcal{A}}\otimes\mathcal{H}_{\cal{B}}.$ (5)

Clearly, $U^{t}=U^{t}|_{\mathcal{H}_{\mathcal{A}}}\otimes
U^{t}|_{\mathcal{H}_{\cal{B}}}$ on $\mathcal{H}$. $\boxtimes$

A suitable space $\mathcal{H}_{\cal{B}}$ in this context is the so-called
Modulated Fock Space [39] where principal eigenfunctions of the Koopman
operator are used in the definition of the space. Using the techniques in [39]
the following can be shown:

###### Theorem 2.2

Let $U$ be a discretization of the Koopman family $U^{t}$ for a system with an
$n$-dimensional Milnor attractor $\mathcal{A}$. Let $\mathcal{H}_{\cal{B}}$ be
the Modulated Fock Space. We assume that $U^{t}|_{\mathcal{H}_{\cal{B}}}$ is
spectral. Then $U$ admits the Schur form

$U=N+Q_{1}$ (6)

where $N$ is normal and $Q_{1}$ is quasi-nilpotent.

Proof: We know that $U^{t}|_{\mathcal{H}_{\mathcal{A}}}$ is scalar, as it is
the space of square integrable functions with respect to $\mu$ rendering
$U^{t}|_{\mathcal{H}_{\mathcal{A}}}$ unitary. The fact that
$U^{t}|_{\mathcal{H}_{\cal{B}}}$ is spectral immediately renders
$U^{t}|_{\mathcal{H}_{\mathcal{A}}}\otimes U^{t}|_{\mathcal{H}_{\cal{B}}}$
spectral, and therefore Theorem 2.1 applies. $\boxtimes$

##  3 Extended DMD

In this section, we review the EDMD following the seminal work [57], upon
which we build our discrete Koopman-Schur decomposition. For the sake of
completeness and for the readers’ convenience when comparing the proposed
method with the EDMD, we provide full description of the main ingredients of
the EDMD. For further study of the EDMD the reader is referred to [57], and
for convergence analysis we refer to [32], [11], [10].

In Section 3.1 we first describe the construction of a finite dimensional data
driven compression $U_{N}$ of the Koopman operator to a selected
$N$-dimensional function subspace ${\mathcal{F}}_{N}$, and in Section 3.2 we
analyze the convergence when $N\rightarrow\infty$ and the number of available
snapshot pairs $m\rightarrow\infty$. If $m<N$, then the proper method to
approximate the eigenvalues and eigenvectors is the Rayleigh-Ritz extraction
and in Section 3.3 we review the details. The kernel trick, that allows for an
efficient computation with extremely large dimension $N$ is briefly described
in Section 3.4. The technical details of point-wise evaluation of the
approximate eigenfunctions and representation of the data snapshots as linear
combinations of the computed eigenfunctions are provided in Sections 3.5 and
3.6, respectively. In addition, as a preparation and a motivation for the new
method in Section 4, we discuss some numerical/implementation issues that are
not tackled in [57].

###  3.1 Numerical data driven compression of $\mathbb{U}$

Numerical approximation of some eigenvalues and eigenfunctions of $\mathbb{U}$
begins with a selection of a finite dimensional subspace
${\mathcal{F}}_{N}=\mathrm{Span}({\boldsymbol{\psi}}_{1},\ldots,{\boldsymbol{\psi}}_{N}),$
where ${\boldsymbol{\psi}}_{1},\ldots,{\boldsymbol{\psi}}_{N}$ are suitable
scalar basis functions defined on $\mathbb{C}^{n}$ (or $\mathbb{R}^{n}$ for
real data). In an ideal situation, an abundance of snapshot data
$(\mathbf{x}_{1},\mathbf{y}_{1}),\ldots(\mathbf{x}_{m},\mathbf{y}_{m})$ is
available ($\mathbf{y}_{i}=\mathbf{F}(\mathbf{x}_{i})$ and $m$ is as large as
we can afford) and with a suitable choice of basis functions we can let
$N\rightarrow\infty$ (i.e. we can make $N$ as large as desired – if we take
polynomials in $n$ variables of the total degree of at most $p$, then
$N=\binom{n+p}{n}$). In fact, using a kernel trick with appropriate functions,
we can have even (implicitly) $N=\infty$. Furthermore, the
${\boldsymbol{\psi}}_{i}$’s are preferably chosen from a sufficiently rich
class of functions with good approximation properties in the ambient function
space. We defer the convergence analysis (when $m,N\rightarrow\infty$) to
Section 3.2, and proceed with a review of numerical details of matrix
representation of the data driven compression, using material from [21], [16],
[17], [41].

Data driven setting means that we are supplied with data snapshots pairs
$(\mathbf{x}_{1},\mathbf{y}_{1}=\mathbf{F}(\mathbf{x}_{1}))$, $\ldots,$
$(\mathbf{x}_{m},\mathbf{y}_{m}=\mathbf{F}(\mathbf{x}_{m}))$, and the values
of the basis functions evaluated at these snapshots are conveniently tabulated
in the $m\times N$ matrices ${\boldsymbol{\Psi}}_{x}$,
${\boldsymbol{\Psi}}_{y}$ as

${\boldsymbol{\Psi}}_{x}=\begin{pmatrix}{\boldsymbol{\psi}}_{1}(\mathbf{x}_{1})&\ldots&{\boldsymbol{\psi}}_{N}(\mathbf{x}_{1})\cr\vdots&\ldots&\vdots\cr{\boldsymbol{\psi}}_{1}(\mathbf{x}_{m})&\ldots&{\boldsymbol{\psi}}_{N}(\mathbf{x}_{m})\end{pmatrix},\;\;{\boldsymbol{\Psi}}_{y}=\begin{pmatrix}{\boldsymbol{\psi}}_{1}(\mathbf{y}_{1})&\ldots&{\boldsymbol{\psi}}_{N}(\mathbf{y}_{1})\cr\vdots&\ldots&\vdots\cr{\boldsymbol{\psi}}_{1}(\mathbf{y}_{m})&\ldots&{\boldsymbol{\psi}}_{N}(\mathbf{y}_{m})\end{pmatrix}.$
(7)

If the data are collected from a single trajectory, then
$\mathbf{y}_{i}=\mathbf{x}_{i+1}$. In general, the data may be collected from
several trajectories and the basic assumption is that the data is organized so
that $\mathbf{y}_{i}=\mathbf{F}(\mathbf{x}_{i})$. It is reasonable to assume
that ${\boldsymbol{\psi}}_{1},\ldots,{\boldsymbol{\psi}}_{N}$ are linearly
independent. However, if those functions are visible only on the snapshots
$\mathbf{x}_{i}$, and $m<N$, then at most $m$ of them will be seen as
independent – in other words, the ranks of the matrices
${\boldsymbol{\Psi}}_{x}$, ${\boldsymbol{\Psi}}_{y}$ are at most $m$. In a
data rich scenario with $m>N$, it is reasonable to assume that the rank is
$N$, but a numericall algorithm should check the numerical rank in case of
ill-conditioning.

Once the data matrices (7) are assembled, an operator compression
$\mathbf{P}_{{\mathcal{F}}_{N}}{\left.\kern-1.2pt\mathbb{U}\vphantom{\big{|}}\right|_{{\mathcal{F}}_{N}}}:\mathcal{F}_{N}\longrightarrow\mathcal{F}_{N}$
is computed numerically, where $\mathbf{P}_{{\mathcal{F}}_{N}}$ is an
appropriately defined projection onto ${\mathcal{F}}_{N}$. Since
${\mathcal{F}}_{N}$ is unlikely to be $\mathbb{U}$-invariant, we have

$\mathbb{U}{\boldsymbol{\psi}}_{i}={\boldsymbol{\psi}}_{i}\circ\mathbf{F}=\sum_{j=1}^{N}\mathbf{u}_{ji}{\boldsymbol{\psi}}_{j}+\rho_{i},$

where $\mathbf{u}_{ji}$, $j=1,\ldots,N$, are the coefficients that minimize
the residual function $\rho_{i}$, $\int|\rho_{i}|^{2}d\mu\rightarrow\min$.
Since in a data driven setting the only available information is (7), the only
feasible action is to solve the discrete least squares (LS) problems

$\int\left|\sum_{j=1}^{N}\mathbf{u}_{ji}{\boldsymbol{\psi}}_{j}-{\boldsymbol{\psi}}_{i}\circ\mathbf{F}\right|^{2}d{\boldsymbol{\delta}}_{m}\longrightarrow\min_{\mathbf{u}_{1i},\ldots,\mathbf{u}_{Ni}},\;\;i=1,\ldots,N,$
(8)

i.e. with respect to the empirical measure defined using the sum of the Dirac
measures concentrated at the $\mathbf{x}_{k}$’s,
${\boldsymbol{\delta}}_{m}=(1/m)\sum_{k=1}^{m}{\boldsymbol{\delta}}_{\mathbf{x}_{k}}$.

###### Remark 3.1

In other words, (8) computes the orthogononal projection of
$\mathbb{U}{\boldsymbol{\psi}}_{i}$ onto $\mathcal{F}_{N}$, computed in the
space $L^{2}({\boldsymbol{\delta}}_{m})$. In the continuous setting in
$L^{2}(\mu)$, it simply reads
$\int\left|\sum_{j=1}^{N}\mathbf{u}_{ji}{\boldsymbol{\psi}}_{j}-{\boldsymbol{\psi}}_{i}\circ\mathbf{F}\right|^{2}d\mu\longrightarrow\min_{\mathbf{u}_{1i},\ldots,\mathbf{u}_{Ni}},\;\;i=1,\ldots,N.$
With properly weighted (quadrature based) algebraic least squares rigorous
error bounds are possible, see e.g. [11], [10].  $\boxtimes$

Since
$(\mathbb{U}{\boldsymbol{\psi}}_{i})(\mathbf{x}_{j})=({\boldsymbol{\psi}}_{i}\circ\mathbf{F})(\mathbf{x}_{j})={\boldsymbol{\psi}}_{i}(\mathbf{y}_{j})$,
the minimization problem (8) reads

$\frac{1}{m}\sum_{j=1}^{m}|\rho_{i}(\mathbf{x}_{j})|^{2}=\frac{1}{m}\left\|\\!\left[{\boldsymbol{\Psi}}_{x}\left(\begin{smallmatrix}\mathbf{u}_{1i}\cr\vdots\cr\mathbf{u}_{Ni}\end{smallmatrix}\right)\\!-\\!\left(\begin{smallmatrix}{\boldsymbol{\psi}}_{i}(\mathbf{y}_{1}))\cr\vdots\cr{\boldsymbol{\psi}}_{i}(\mathbf{y}_{m})\end{smallmatrix}\right)\\!\right]\\!\\!\right\|_{2}^{2}\longrightarrow\min_{\mathbf{u}_{1i},\ldots,\mathbf{u}_{Ni}},\;\;i=1,\ldots,N.$
(9)

To ensure well defined unique matrix representation of
$\mathbf{P}_{{\mathcal{F}}_{N}}{\mathbb{U}}_{|{\mathcal{F}}_{N}}$ (independent
of the rank of ${\boldsymbol{\Psi}}_{x}$), an additional constraint
$\sum_{j=1}^{N}|\mathbf{u}_{ji}|^{2}\rightarrow\min$ is imposed for all $i$,
yielding
$(\mathbf{u}_{ji})_{j=1}^{N}={\boldsymbol{\Psi}}_{x}^{\dagger}({\boldsymbol{\psi}}_{i}(\mathbf{y}_{k}))_{k=1}^{m}={\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}(;,i)$.
Here ${\boldsymbol{\Psi}}_{x}^{\dagger}$ denotes the Moore-Penrose generalized
inverse of ${\boldsymbol{\Psi}}_{x}$, and the additional constraint (built
into the definition of ${\boldsymbol{\Psi}}_{x}^{\dagger}$) means that the
unique shortest vector from the solution manifold is selected.

Hence,
$\mathbf{P}_{{\mathcal{F}}_{N}}{\mathbb{U}}_{|{\mathcal{F}}_{N}}:\mathcal{F}_{N}\longrightarrow\mathcal{F}_{N}$,
defined as

$\mathbf{P}_{{\mathcal{F}}_{N}}{\mathbb{U}}_{|{\mathcal{F}}_{N}}{\boldsymbol{\psi}}_{i}=\sum_{j=1}^{N}\mathbf{u}_{ji}{\boldsymbol{\psi}}_{j},\;\;i=1,\ldots,N,$

is in the basis
$\mathcal{B}=\\{{\boldsymbol{\psi}}_{1},\ldots,{\boldsymbol{\psi}}_{N}\\}$
represented by the uniquelly determined matrix

$\displaystyle\\!\\![\mathbf{P}_{{\mathcal{F}}_{N}}\mathbb{U}_{|\mathcal{F}_{N}}]_{\mathcal{B}}={\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}\equiv
U_{N}=(\mathbf{u}_{ji})_{j,i=1}^{N},$ (10)
$\displaystyle\\!\\!\mathbf{P}_{{\mathcal{F}}_{N}}\mathbb{U}_{|\mathcal{F}_{N}}\\!(\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}\\!\\!\begin{pmatrix}\mathsf{f}_{1}\cr\vdots\cr\mathsf{f}_{N}\end{pmatrix})\\!=\\!\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}(U_{N}\\!\\!\begin{pmatrix}\mathsf{f}_{1}\cr\vdots\cr\mathsf{f}_{N}\end{pmatrix}).$

If we introduce the notation
${\boldsymbol{\psi}}(x)=({\boldsymbol{\psi}}_{1}(x),\ldots,{\boldsymbol{\psi}}_{N}(x))^{T}\in\mathbb{C}^{N}$,
and similarly for the residuals,
${\boldsymbol{\rho}}(x)=(\rho_{1}(x),\ldots,\rho_{N}(x))^{T}\in\mathbb{C}^{N}$,
then the action of $\mathbb{U}$ on ${\mathcal{F}}_{N}$ can be compactly
written as

$\mathbb{U}({\boldsymbol{\psi}}(x)^{T}\mathbf{z})={\boldsymbol{\psi}}(x)^{T}U_{N}\mathbf{z}+{\boldsymbol{\rho}}(x)^{T}\mathbf{z},\;\;x\in\mathbb{C}^{n},\;\mathbf{z}\in\mathbb{C}^{N}.$
(11)

In the sequel, we ease the notation and write
$\mathbb{U}{\boldsymbol{\psi}}(x)^{T}\mathbf{z}$ instead of
$\mathbb{U}({\boldsymbol{\psi}}(x)^{T}\mathbf{z})$, and for a sequence of
functions $\zeta_{1}(x),\zeta_{2}(x),\ldots$ we define
$\mathbb{U}(\zeta_{1}(x),\zeta_{2}(x),\ldots)$ as
$(\mathbb{U}\zeta_{1}(x),\mathbb{U}\zeta_{2}(x),\ldots)$.

###  3.2 Convergence as $m\rightarrow\infty$

In this section we discuss the convergence of
$\mathbf{P}_{{\mathcal{F}}_{N}}{\left.\kern-1.2pt\mathbb{U}\vphantom{\big{|}}\right|_{{\mathcal{F}}_{N}}}:\mathcal{F}_{N}\longrightarrow\mathcal{F}_{N}$
to ${\mathbb{U}}$ as $N\rightarrow\infty$. As shown in [32] it is enough to
show this convergence in order to establish convergence of the EDMD operator
$\mathbf{P}_{{\mathcal{F}}_{N},m}{\left.\kern-1.2pt\mathbb{U}\vphantom{\big{|}}\right|_{{\mathcal{F}}_{N}}}:\mathcal{F}_{N}\longrightarrow\mathcal{F}_{N}$
to $\mathbb{U}$ as $N,m\rightarrow\infty$666The discussion of this issue in
[32] pertains to the functional space being $L^{2}(\mu)$, but can be
established for any separable Hilbert spaces adapted to dynamical systems,
considered in [39]..

For an operator $A$, we denote by $R_{A}(\lambda)=(A-\lambda
I)^{-1},\lambda\in\rho(A),R_{A}:\rho(A)\rightarrow{\cal{B}}(\mathcal{H})$, its
resolvent function, where $\rho(A)=\mathbb{C}/\sigma(A)$ is the resolvent set,
$\sigma(A)$ is the spectrum of $A$, and ${\cal{B}}(\mathcal{H})$ the space of
bounded linear operators on $\mathcal{H}$. We start with the notion of strong
resolvent convergence:

###### Definition 3.1

A sequence of operators $B_{n}\in{\cal{B}}(\mathcal{H})$ converges strongly to
$B\in{\cal{B}}(\mathcal{H})$ provided $B_{n}h\rightarrow Bh,\ \ \forall
h\in\mathcal{H}.$

We assume $\mathbf{P}_{{\mathcal{F}}_{N}}\rightarrow I$ strongly with $N$,
where $I$ is the identity operator on $\mathcal{H}$ (see [32] for some
conditions on the basis functions and $\mathcal{H}$ under which this is true).

###### Definition 3.2

Let $A_{n}\in{\cal{B}}(\mathcal{H}),n\in\mathbb{N}$ and
$A\in{\cal{B}}(\mathcal{H})$. We say that $A_{n}$ converges to $A$ in the
strong resolvent sense if $R_{A_{n}}\rightarrow R_{A}$ strongly for any - and
thus every - $\lambda\in\Gamma$ where

$\Gamma=\mathbb{C}/(\sigma(A)\cup_{n}\sigma(A_{n}))=\mathbb{C}/\Sigma.$ (12)

###### Lemma 3.1

Assume $\mathbb{U}$ admits the Kato spectral expansion

$\mathbb{U}=\sum_{j}\lambda_{j}P^{s}_{j}+D_{j}$ (13)

where $D_{j}$ is quasinilpotent. Then
$\mathbf{P}_{{\mathcal{F}}_{N}}{\left.\kern-1.2pt\mathbb{U}\vphantom{\big{|}}\right|_{{\mathcal{F}}_{N}}}$
converges with $N$ to ${\mathbb{U}}$ in the strong resolvent sense.

Proof: The projections $\mathbb{U}_{N}$ admit the expansion

$\mathbb{U}_{N}=\sum_{j}\lambda_{j}P_{N}P^{s}_{j}P_{N}+P_{N}D_{j}P_{N}$ (14)

where each eigenvalue $\lambda_{j}$ has finite algebraic multiplicity $m_{j}$
and the resolvent operator is

$R_{\mathbb{U}_{N}}(\lambda)=\sum_{j}\frac{P_{N}P^{s}_{j}P_{N}}{\lambda-\lambda_{j}}+\sum_{n=1}^{m_{j}-1}\frac{P_{N}D^{n}_{j}P_{N}}{(\lambda-\lambda_{j})^{-n-1}}.$

Thus, $R_{\mathbb{U}_{N}}$ strongly converges to $R_{\mathbb{U}}(\lambda)$.
$\boxtimes$

###### Proposition 3.1

Assume the spectrum of $\mathbb{U}$ is contained inside the unit circle. Let
$\mathbb{U}_{N}$ converge with $N$ to ${\mathbb{U}}$ in the strong resolvent
sense. then $f(\mathbb{U}_{N})$ converges to $f(\mathbb{U})$ strongly for
every bounded $C^{\infty}$ function $f:\Sigma\rightarrow\mathbb{C}$.

Proof: We need the following Lemma:

###### Lemma 3.2

The spectrum of $\mathbb{U}_{N}$ is inside the unit circle.

Proof: Follows from (14) . $\boxtimes$

Now we finish the proof of the Proposition: There is a $C^{\infty}$ functional
calculus for scalar operators [34], that, together with stability established
in the Lemma 3.2, implies the result. $\boxtimes$

Now we can discuss spectral convergence issues. We have the following:

###### Theorem 3.2

The sequence $\mathbb{U}_{N}$ is spectraly inclusive, i.e.
$\sigma(\mathbb{U})\subseteq\lim_{N\rightarrow\infty}\sigma(\mathbb{U}_{N})$.

Proof: Suppose not. Then there is a $\lambda\in\mathbb{C}$ and $\epsilon>0$
such that $\sigma(\mathbb{U}_{N})\cap B_{\epsilon}(\lambda)=\emptyset$ where
$B_{\epsilon}(\lambda)\subset\mathbb{C}$ is an open ball of radius $\epsilon$
centered at $\lambda$. Consider a $C^{\infty}$ function $f$ that is 1 on
$B_{\epsilon/2}(\lambda)$ and zero outside of $B_{\epsilon}(\lambda)$. Then
$f(\mathbb{U}_{N})=0$ implying $f(\mathbb{U})g=\lim f(\mathbb{U}_{N})g=0$ for
every $g$. Since $\lambda\in\sigma(\mathbb{U})$ we have that the range of the
spectral projection $P(B_{\epsilon/2}(\lambda))$ is not empty, containing a
function $h$ such that $f(\mathbb{U})h=h$, leading to a contradiction.
$\boxtimes$

###  3.3 The Rayleigh quotient of $U_{N}$

The $N\times N$ matrix
$U_{N}={\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}$ is in in
Section 3.1 defined as a matrix representation of the compression
$\mathbf{P}_{{\mathcal{F}}_{N}}{\mathbb{U}}_{|{\mathcal{F}}_{N}}:\mathcal{F}_{N}\longrightarrow\mathcal{F}_{N}$,
where the projection is executed using algebraic least squares – in the
limited data driven setting this is a reasonable approximation. But, if $N>m$,
it is of low rank of at most $m$, and potentially of large dimension.
Moreover, in this case ${\boldsymbol{\Psi}}_{x}$ and ${\boldsymbol{\Psi}}_{y}$
are also of rank of at most $m$, and the data does not contain information to
uniquely define $U_{N}$ in the sense of the formulation in Section 3.1.

Indeed, if $E$ is such that ${\boldsymbol{\Psi}}_{x}E=\mathbf{0}$, then
${\boldsymbol{\Psi}}_{x}(U_{N}+E)={\boldsymbol{\Psi}}_{x}U_{N}={\boldsymbol{\Psi}}_{x}{\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}=\mathbf{P}_{x}{\boldsymbol{\Psi}}_{y}$,
where $\mathbf{P}_{x}$ is the orthogonal projector onto the range of
${\boldsymbol{\Psi}}_{x}$. The formula
$U_{N}={\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}$ is just a
particular way to uniquely select a point in the linear manifold of the
solutions of (8, 9). It is appealing because of the minimal norm property and
the corresponding theory of the Moore-Penrose generalized inverse, but it is
not clear why should this choice be considered intrinsic.

The key observation is that the ambiguity in defining the matrix of the
compressed operator is removed by taking the Rayleigh quotient of $U_{N}$ with
respect to the range of ${\boldsymbol{\Psi}}_{x}$:

$({\boldsymbol{\Psi}}_{x}^{\dagger})^{\dagger}(U_{N}+E){\boldsymbol{\Psi}}_{x}^{\dagger}={\boldsymbol{\Psi}}_{x}U_{N}{\boldsymbol{\Psi}}_{x}^{\dagger}={\boldsymbol{\Psi}}_{x}{\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}{\boldsymbol{\Psi}}_{x}^{\dagger}=\mathbf{P}_{x}{\boldsymbol{\Psi}}_{y}{\boldsymbol{\Psi}}_{x}^{\dagger}.$

For the reasons of numerical stability, the matrix representation of the
Rayleigh quotient should be computed using a unitary basis of the range of
${\boldsymbol{\Psi}}_{x}^{\dagger}$, or of some of its subspaces. Further,
even when $m>N$ and all matrices are of full rank $N$,
${\boldsymbol{\Psi}}_{x}$ can be ill-conditioned so that the numerical
solution of the least squares problem becomes difficult and the numerical rank
must be carefully determined.

To find a unitary basis and determine the correct dimension, a unitary rank
revealing decomposition is used. The SVD is preferred because of its
theoretical and numerical properties, and because of availability of robust
and computationally efficient software implementations. Let
$r=\mathrm{rank}({\boldsymbol{\Psi}}_{x})$ and let the SVD of
${\boldsymbol{\Psi}}_{x}$ be

${\boldsymbol{\Psi}}_{x}=W\Sigma
V^{*},\;\;\Sigma=\mathrm{diag}(\sigma_{i})_{i=1}^{r},\,W\in\mathbb{C}^{m\times
r},\;V\in\mathbb{C}^{N\times r},\;W^{*}W=V^{*}V=\mathbb{I}_{r}.$ (15)

(Here we write the _economy size_ SVD, i.e. the zero singular values and the
corresponding singular vectors are chopped out from the full SVD.)

Using ${\boldsymbol{\Psi}}_{x}^{\dagger}=V\Sigma^{-1}W^{*}$, the Rayleigh
quotient $\widehat{U}=V^{*}U_{N}V$ is the $r\times r$ matrix

$\widehat{U}=V^{*}U_{N}V=V^{*}(V\Sigma^{-1}W^{*}){\boldsymbol{\Psi}}_{y}V=\Sigma^{-1}W^{*}{\boldsymbol{\Psi}}_{y}V.$
(16)

It can be easily checked that $U_{N}V=V\widehat{U}$. We can express $V$ as
$V={\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1}$ and then $\widehat{U}$ as

$\widehat{U}=V^{*}U_{N}V=V^{*}(V\Sigma^{-1}W^{*}){\boldsymbol{\Psi}}_{y}({\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1})=\Sigma^{-1}W^{*}({\boldsymbol{\Psi}}_{y}{\boldsymbol{\Psi}}_{x}^{*})W\Sigma^{-1}.$
(17)

The reason for considering (17) instead of the simpler formula
$\widehat{U}=V^{*}U_{N}V$ will become apparent in Section 3.4.

####  3.3.1 Why is the SVD truncated?

In practical finite precision computation, the rank of a matrix is determined
as the numerical rank, based on the computed singular values and guided by the
Eckart-Young-Mirsky theorem. Since this involves truncation of small (non-
zero) singular values, the numerical rank could differ from the true rank. It
should be stressed that determination of numerical rank is a nontrivial
problem, see e.g. [25].

In the context of the Rayleigh-Ritz procedure of computing the eigenvalues and
eigenvectors based on the available data ${\boldsymbol{\Psi}}_{x}$,
${\boldsymbol{\Psi}}_{y}$, the issue is not about getting a good low rank
approximation of ${\boldsymbol{\Psi}}_{x}$ but about extracting spectral
information based on the pair ${\boldsymbol{\Psi}}_{x}$,
${\boldsymbol{\Psi}}_{y}$. The SVD is used to compute an orthonormal basis so
that the Rayleigh-Ritz algorithm can be numerically more stable. Using the
low-rank approximation means discarding information that is in the range of
${\boldsymbol{\Psi}}_{x}$, so it is better to get as large as possible
orthonormal basis, i.e. more left singular vectors. Unfortunately, the
smallest singular values (and the corresponding singular vectors) are often
computed very inaccurately and the data driven formula (16) for the Rayleigh
quotient is contaminated by large errors if small singular values are used. If
a more accurate SVD algorithm capable of getting even small singular values
accurately (when feasible) is available, then they can be included. Another
important detail is that small singular values are a consequence of
unfavourable scaling of the data matrices. These important issues are
discussed in detail in [15], [14], [19], [18], (For software details see [16],
[17].)

Another important issue is the possibility of losing (through truncation) low
energy unstable modes, see e.g. a discussion in [1].

###  3.4 The EDMD with the kernel trick

The key issue is practical computation when the data snapshots are mapped in
ever higher dimensional space, i.e. when $N\rightarrow\infty$. The matrices
$W$ and $\Sigma$ in (15), (17) can be computed from the spectral decomposition
of the $m\times m$ cross-product
$C_{xx}={\boldsymbol{\Psi}}_{x}{\boldsymbol{\Psi}}_{x}^{*}$, since
$C_{xx}=W\Sigma^{2}W^{*}$. The cross-product matrix
$C_{yx}={\boldsymbol{\Psi}}_{y}{\boldsymbol{\Psi}}_{x}^{*}$, used in (17), is
$m\times m$ and its elements, as those of $C_{xx}$, are inner products of
vectors from $N$-dimensional space. Hence, the only place where the
potentially extremely large dimension $N$ incurs computational cost is in
computing $C_{xx}$ and $C_{yx}$.

Note that $C_{xx}=C_{xx}^{*}$, and that for all $i,j$ the elements of these
matrices can be expressed using the Euclidean inner product
$\langle\cdot,\cdot\rangle_{\mathbb{C}^{N}}$ as

$\displaystyle(C_{xx})_{ij}$ $\displaystyle=$
$\displaystyle({\boldsymbol{\Psi}}_{x}{\boldsymbol{\Psi}}_{x}^{*})_{ij}=\sum_{k=1}^{N}{\boldsymbol{\psi}}_{k}(\mathbf{x}_{i})\overline{{\boldsymbol{\psi}}_{k}(\mathbf{x}_{j})}=\langle{\boldsymbol{\psi}}(\mathbf{x}_{i}),{\boldsymbol{\psi}}(\mathbf{x}_{j})\rangle_{\mathbb{C}^{N}},$
(18) $\displaystyle(C_{yx})_{ij}$ $\displaystyle=$
$\displaystyle({\boldsymbol{\Psi}}_{y}{\boldsymbol{\Psi}}_{x}^{*})_{ij}=\sum_{k=1}^{N}{\boldsymbol{\psi}}_{k}(\mathbf{y}_{i})\overline{{\boldsymbol{\psi}}_{k}(\mathbf{x}_{j})}=\langle{\boldsymbol{\psi}}(\mathbf{y}_{i}),{\boldsymbol{\psi}}(\mathbf{x}_{j})\rangle_{\mathbb{C}^{N}}.$
(19)

This is where the kernel trick comes to the rescue. The trick is to work with
the basis functions ${\boldsymbol{\psi}}_{k}$ implicitly by using a suitable
function ${\boldsymbol{f}}$ such that

${\boldsymbol{f}}(x,y)=\langle{\boldsymbol{\psi}}(x),{\boldsymbol{\psi}}(y)\rangle_{\mathbb{C}^{N}},\;\;x,y\in\mathbb{C}^{n},$
(20)

so that $(C_{xx})_{ij}={\boldsymbol{f}}(\mathbf{x}_{i},\mathbf{x}_{j})$,
$(C_{yx})_{ij}={\boldsymbol{f}}(\mathbf{y}_{i},\mathbf{x}_{j})$. Since the
construction starts with the function ${\boldsymbol{f}}$, its defining
property must be the definiteness in the sense of the Mercer’s theorem. With a
broad palette of kernel functions, so that computation of the Rayleigh
quotient (17) is feasible for large dimensions, even when
$N\rightarrow\infty$. For more details we refer to [57].

###### Remark 3.3

The kernel trick is an excellent example of trade-off between computational
efficiency and numerical accuracy. The key of the accuracy of the computed SVD
of ${\boldsymbol{\Psi}}_{x}$ is the condition number
$\kappa_{2}({\boldsymbol{\Psi}}_{x})=\|{\boldsymbol{\Psi}}_{x}\|_{2}\|{\boldsymbol{\Psi}}_{x}^{\dagger}\|_{2}$.
On the other hand, the condition number of
$C_{xx}={\boldsymbol{\Psi}}_{x}{\boldsymbol{\Psi}}_{x}^{*}$ is squared,
$\kappa_{2}(C_{xx})=\kappa_{2}({\boldsymbol{\Psi}}_{x})^{2}$.  $\boxtimes$

####  3.4.1 Connection with the DMD

The DMD is a special case of EDMD with $N=n$ and with the coordinate functions
${\boldsymbol{\psi}}_{i}(x)=e_{i}^{T}x$ as basis functions. The kernel
function is ${\boldsymbol{f}}(x,y)=\langle x,y\rangle_{\mathbb{C}^{n}}$.

In the DMD, the snapshots are collected column-wise into the matrices
$\mathbf{X}=(\mathbf{x}_{1}\;\ldots\;\mathbf{x}_{m})$,
$\mathbf{Y}=(\mathbf{y}_{1}\;\ldots\;\mathbf{y}_{m})$. Note that
$\mathbf{X}={\boldsymbol{\Psi}}_{x}^{T}$,
$\mathbf{Y}={\boldsymbol{\Psi}}_{y}^{T}$. The DMD matrix is
$A_{n}=\mathbf{Y}\mathbf{X}^{\dagger}={\boldsymbol{\Psi}}_{y}^{T}{\boldsymbol{\Psi}}_{x}^{T\dagger}$,
hence $A_{n}=A_{N}=U_{N}^{T}$. The SVD of $\mathbf{X}$ (using (15)) reads
$\mathbf{X}=V^{*T}\Sigma W^{T}$, so that the $r\times r$ Rayleigh quotient
$\widehat{A}$ used in the DMD is

$\widehat{A}=V^{T}A_{N}V^{*T}=(V^{*}U_{N}V)^{T}=\widehat{U}^{T}.$

Here too, in numerical computations, $r$ is the numerical rank, and
$\widehat{A}$ is computed using the formula

$\widehat{A}=V^{T}A_{N}V^{*T}=V^{T}\mathbf{Y}W^{T*}\Sigma^{-1}V^{T}V^{*T}=V^{T}\mathbf{Y}W^{T*}\Sigma^{-1}.$
(21)

In DMD, a spectral decomposition of $\widehat{A}$, $\widehat{A}=G\Lambda
G^{-1}$, is computed and the approximate eigenvectors of $A_{N}$ are the
columns of $V^{*T}G$.

For a numerically robust implementation of the DMD see [15], [16] and for a
special case of Hermitian DMD see [17].

###### Remark 3.4

The DMD is often in the literature interpreted as a mere linear regression
problem, i.e. $A_{n}=\mathbf{Y}\mathbf{X}^{\dagger}$ is taken as a mere
solution of the linear least squares problem
$\|\mathbf{Y}-A_{n}\mathbf{X}\|_{F}\rightarrow\min$, with the solution
$A_{n}=\mathbf{Y}\mathbf{X}^{\dagger}$, and ignoring the connection with the
Koopman linearization of the underlying dynamics and the fact that $A_{n}$ is,
in the case of tall data matrices, just one out of many solutions of the LS
problem. See [19], [18] for more detailed discussion.777See also [16], [17].
$\boxtimes$

###  3.5 Evaluation of the approximate eigenfunctions of $\mathbb{U}$

If we compute the $r$ eigenpairs $(\omega_{i},s_{i})$ of $\widehat{U}$,
$\widehat{U}s_{i}=\omega_{i}s_{i}$, $i=1,\ldots,r$, then

$U_{N}Vs_{i}={\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}Vs_{i}=V\Sigma^{-1}W^{*}{\boldsymbol{\Psi}}_{y}{\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1}s_{i}=\omega_{i}Vs_{i},$

so that $(\omega_{i},Vs_{i})$ are the eigenpairs of $U_{N}$. The approximate
eigenfunctions of $\mathbb{U}$ are

${\boldsymbol{\phi}}_{i}(x)=\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&{\boldsymbol{\psi}}_{2}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}Vs_{i}=\sum_{j=1}^{N}{\boldsymbol{\psi}}_{j}(x)(Vs_{i})_{j},\;\;i=1,\ldots,r.$
(22)

This is understood in the sense that

$\displaystyle\mathbb{U}{\boldsymbol{\phi}}_{i}(x)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&{\boldsymbol{\psi}}_{2}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}U_{N}Vs_{i}+{\boldsymbol{\rho}}(x)^{T}Vs_{i}$
(23) $\displaystyle=$
$\displaystyle\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&{\boldsymbol{\psi}}_{2}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}\omega_{i}Vs_{i}+{\boldsymbol{\rho}}(x)^{T}Vs_{i}$
(24) $\displaystyle=$
$\displaystyle\omega_{i}{\boldsymbol{\phi}}_{i}(x)+{\boldsymbol{\rho}}(x)^{T}Vs_{i}\backsimeq\omega_{i}{\boldsymbol{\phi}}_{i}(x).$
(25)

(Recall that $\backsimeq$ is a shorthand that means equal up to a residual
that is minimized over the available data in the sense of least squares.)

In a practical computation, the ${\boldsymbol{\phi}}_{i}$’s are not
numerically evaluated using the definition (22), as it involves the dimension
$N$.

####  3.5.1 Evaluation at the snapshots $\mathbf{x}_{j}$

Evaluating the eigenfunctions ${\boldsymbol{\phi}}_{i}$ at the snapshots
$\mathbf{x}_{j}$ means tabulating the values

${\boldsymbol{\Phi}}_{x}=\begin{pmatrix}{\boldsymbol{\phi}}_{1}(\mathbf{x}_{1})&\ldots&{\boldsymbol{\phi}}_{r}(\mathbf{x}_{1})\cr\vdots&\ddots&\vdots\cr{\boldsymbol{\phi}}_{1}(\mathbf{x}_{m})&\ldots&{\boldsymbol{\phi}}_{r}(\mathbf{x}_{m})\end{pmatrix}={\boldsymbol{\Psi}}_{x}VS\in\mathbb{C}^{m\times
r},$

where the relation ${\boldsymbol{\Phi}}_{x}={\boldsymbol{\Psi}}_{x}VS$ follows
directly from the definition of the matrix ${\boldsymbol{\Psi}}_{x}$ and the
eigenfunctions (22). Next, observe that the SVD
${\boldsymbol{\Psi}}_{x}=W\Sigma V^{*}$ eliminates $V$ so that we have
${\boldsymbol{\Phi}}_{x}=W\Sigma S$. In practice, the SVD of
${\boldsymbol{\Psi}}_{x}$ may not be feasible and the kernel trick is used, so
that the same conclusion follows from the spectral decomposition of
${\boldsymbol{\Psi}}_{x}{\boldsymbol{\Psi}}_{x}^{*}$ and the relation
$V={\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1}$ :

${\boldsymbol{\Phi}}_{x}={\boldsymbol{\Psi}}_{x}{\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1}S=W\Sigma^{2}W^{*}W\Sigma^{-1}S=W\Sigma
S.$

####  3.5.2 Evaluation at an arbitrary $\mathbf{x}\in\mathbb{C}^{n}$

To evaluate ${\boldsymbol{\phi}}_{i}(\mathbf{x})$ at an arbitrary
$\mathbf{x}\in\mathbb{C}^{n}$ we need

${\boldsymbol{\phi}}_{i}(\mathbf{x})={\boldsymbol{\psi}}(\mathbf{x})^{T}Vs_{i}={\boldsymbol{\psi}}(\mathbf{x})^{T}{\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1}s_{i},$

where, using the kernel function ${\boldsymbol{f}}(\cdot,\cdot)$,

$({\boldsymbol{\psi}}(\mathbf{x})^{T}{\boldsymbol{\Psi}}_{x}^{*})_{j}=\sum_{k=1}^{N}{\boldsymbol{\psi}}_{k}(\mathbf{x})\overline{{\boldsymbol{\psi}}_{k}(\mathbf{x}_{j})}=\langle{\boldsymbol{\psi}}(\mathbf{x}),{\boldsymbol{\psi}}(\mathbf{x}_{j})\rangle_{\mathbb{C}^{N}}={\boldsymbol{f}}(\mathbf{x},\mathbf{x}_{j}).$

Hence

${\boldsymbol{\phi}}_{i}(\mathbf{x})=\begin{pmatrix}{\boldsymbol{f}}(\mathbf{x},\mathbf{x}_{1})&{\boldsymbol{f}}(\mathbf{x},\mathbf{x}_{2})&\ldots&{\boldsymbol{f}}(\mathbf{x},\mathbf{x}_{m})\end{pmatrix}W\Sigma^{-1}s_{i}.$

###  3.6 Snapshot representation

Next computational task is to represent the snapshots
$\mathbf{x}_{1},\ldots,\mathbf{x}_{m}$ using the computed eigenfunctions
${\boldsymbol{\phi}}_{1},\ldots,{\boldsymbol{\phi}}_{r}$. More generally,
given $\ell$ scalar functions $g_{1},\ldots,g_{\ell}$, we seek an
approximation

$\begin{pmatrix}g_{1}(x)&\ldots&g_{\ell}(x)\end{pmatrix}\backsimeq\begin{pmatrix}{\boldsymbol{\phi}}_{1}(x)&\ldots&{\boldsymbol{\phi}}_{r}(x)\end{pmatrix}\Gamma,\;\;\Gamma\in\mathbb{C}^{r\times\ell}.$

The approximation is computed in the least squares sense, column-wise, e.g.

$\Gamma(:,i)=\mathrm{argmin}_{h\in\mathbb{C}^{r}}\|\left(\begin{smallmatrix}g_{i}(\mathbf{x}_{1})\cr\vdots\cr
g_{i}(\mathbf{x}_{m})\end{smallmatrix}\right)-{\boldsymbol{\Phi}}_{x}h\|_{2}={\boldsymbol{\Phi}}_{x}^{\dagger}\left(\begin{smallmatrix}g_{i}(\mathbf{x}_{1})\cr\vdots\cr
g_{i}(\mathbf{x}_{m})\end{smallmatrix}\right),\;\;i=1,\ldots,\ell.$ (26)

For snapshot reconstruction, $\ell=n$ and $g_{i}(x)=e_{i}^{T}x$ is the $i$-th
coordinate function, and in the sequel we consider this special case. If we
define $\mathbf{X}=(\mathbf{x}_{1}\;\ldots\;\mathbf{x}_{m})$ then the column-
wise definition (26) of $\Gamma$ gives (in the least squares sense)
$\mathbf{X}^{T}\approx{\boldsymbol{\Phi}}_{x}\Gamma$ and the solution matrix
$\Gamma$ can be compactly written as

$\Gamma={\boldsymbol{\Phi}}_{x}^{\dagger}\mathbf{X}^{T}=(W\Sigma
S)^{\dagger}\mathbf{X}^{T}=S^{-1}\Sigma^{-1}W^{*}\mathbf{X}^{T}.$

Hence, we have the representation of the snapshots (in a least squares sense)

$\mathbf{X}=(\mathbf{x}_{1}\;\ldots\;\mathbf{x}_{m})\approx\Gamma^{T}{\boldsymbol{\Phi}}_{x}^{T},\;\;\Gamma^{T}=\mathbf{X}W^{*T}\Sigma^{-1}S^{-T}=(\gamma_{1}\;\ldots\;\gamma_{r}).$

Hence, each snapshot has a modal decomposition

$\mathbf{x}_{i}\approx\sum_{j=1}^{r}\gamma_{j}{\boldsymbol{\phi}}_{j}(\mathbf{x}_{i}),\;\;i=1,\ldots,m.$

Note that $\widehat{U}=S\Omega S^{-1}$ yields $\widehat{U}^{T}=S^{-T}\Omega
S^{T}$, i.e. the columns of $S^{-T}$ are the eigenvectors of
$\widehat{U}^{T}$. (In the case of DMD, $\widehat{U}^{T}=\widehat{A}$; see
Section 3.4.1.)

###### Remark 3.5

Note that both $S$ and its inverse appear in the formuals – for the values
${\boldsymbol{\phi}}_{k}(\mathbf{x}_{i})$ (the matrix
${\boldsymbol{\Phi}}_{x}$) we need $S$, and for the modes $\gamma_{j}$ (the
matrix $\Gamma$) we need $S^{-T}$. There is an important caveat here: $S$ may
be severely ill-conditioned (even numerically singular) and the method
requires using $S^{-1}$. This important issue is addressed in Section 4.
$\boxtimes$

####  3.6.1 Connection to the Schmid’s DMD

In the case of DMD (recall Section 3.4.1),
${\boldsymbol{\Psi}}_{x}=\mathbf{X}^{T}=W\Sigma V^{*}$,
$\mathbf{X}=V^{*T}\Sigma W^{T}$ and

$\Gamma^{T}=V^{*T}\Sigma W^{T}W^{*T}\Sigma^{-1}S^{-T}=V^{*T}S^{-T}=V^{*T}G.$

##  4 Koopman-Schur decomposition

In the EDMD and the DMD (see Sections 3.4.1, 3.5) computational frameworks it
is tacitly assumed that the Rayleigh quotient matrix $\widehat{U}$
($\widehat{A}$ in the case of DMD) is diagonalizable, and the whole
construction is build on the decomposition

$\widehat{U}=S\Omega
S^{-1},\;\;\Omega=\mathrm{diag}(\omega_{i})_{i=1}^{r},\;\;S=(s_{1},\ldots,s_{r}),\;\widehat{U}s_{i}=\omega_{i}s_{i}.$
(27)

Here $\widehat{U}=V^{*}U_{N}V$, where $U_{N}$ represents a compression of
$\mathbb{U}$ onto the range of the orthonormal matrix $V$; see (17) in Section
3.3.

In Section 1, we used a $2\times 2$ example to illustrate the problem of non-
existence of the full system of eigenvectors. Furthermore, in practical
computation, although diagonalizability of a matrix is not a priori
guaranteed, a numerical algorithm will most probably return the maximal number
of $r$ eigenvectors, even if there exists only one. The following example will
illustrate.

###### Example 4.1

We use the Matlab’s function eig to compute the eigenvalues and eigenvectors
of a Jordan block of dimension $10$. The computed eigenvector matrix $S$ will
be singular, without any warning/error message, as shown by the following
simple code

    >> Uhat=diag(ones(10,1),1);
    >> [S,~]=eig(Uhat); [G,~]=eig(Uhat’);
    >> [ cond(S) cond(G) ]
    ans =
       Inf   Inf

More generally, we can think of $\widehat{U}$ being equal or close to a matrix
whose Jordan normal form contains several (Jordan) blocks of different sizes,
corresponding to same and/or different eigenvalues. As shown in Sections 3.4.1
and 3.6, the eigenvector matrices $S$ and $G$ are the ingredients of the
matrix of the Koopman modes, which is in this case of numerical rank one (in
the general case possibly highly ill-conditioned and numerically rank
deficient) and as such it has a rather limited potential for discovering
latent properties.

Recall that matrices with non-diagonal Jordan form are a nowhere dense set, so
it is rather unlikely that our $\widehat{U}$ is such. However, even if
diagonalizable, $\widehat{U}$ can be highly non-normal and the matrix $S$ can
be numerically ill-conditioned, and the eigenvectors can be so volatile that
using them numerically may not have sense at all. On the other hand, those
eigenvectors are the key ingredients in the modal decomposition in the EDMD –
recall Remark 3.5.

It is also possible that some individual eigenvectors themselves are sensitive
to perturbations, but a spectral subspace that encompasses them is much more
robust. Multiple or closely clustered eigenvalues that make a cluster well
separated from the rest of the spectrum have such eigenvectors.

These nontrivial numerical difficulties have been the main motivation and a
driving force for the development of an alternative approach that is based on
the ordered partial Schur decomposition. The key idea is to avoid the spectral
decomposition (diagonalization, Jordan normal form) of the compressed operator
$\mathbb{U}$ (see Section 3.1) and to use numerically more robust Schur
decomposition. From a numerical linear algebra perspective, this seems
plausible and it only remains to adapt this approach to the DMD/EDMD and to
ensure the same functionality in the applications. In this section we provide
the details of the new approach.

###  4.1 Partial Schur decomposition of $\mathbb{U}$ and $U_{N}$

Following the above discussion, let us, instead of the diagonalization (27),
compute the Schur form (Theorem 1.1) of $\widehat{U}$:

$\widehat{U}=QTQ^{*},\;\;T=\left(\begin{smallmatrix}t_{11}&\ldots&t_{1r}\cr&\ddots&\vdots\cr&&t_{rr}\end{smallmatrix}\right),\;\;t_{ii}=\omega_{i},$
(28)

where $Q$ is unitary and $T$ is upper triangular with the eigenvalues of
$\widehat{U}$, in some order, on the diagonal. Then

$U_{N}(VQ)={\boldsymbol{\Psi}}_{x}^{\dagger}{\boldsymbol{\Psi}}_{y}VQ=V\Sigma^{-1}W^{*}{\boldsymbol{\Psi}}_{y}{\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1}Q=V\widehat{U}Q=(VQ)T$
(29)

represents a partial Schur form of $U_{N}$ – the matrix $VQ$ is $N\times r$
with $N>r$ (including the case of $N=\infty\gg r$). Its range is an
$r$-dimensional $U_{N}$-invariant subspace and the corresponding compression
of $U_{N}$ is $T$. If we set $Z=VQ$, then

$U_{N}Z=ZT,\;\;\mbox{and, more
generally,}\;\;U_{N}^{k}Z=ZT^{k},\;\;k=1,2,\ldots$ (30)

In fact, $f(U_{N})Z=Zf(T)$, where $f$ is any analytic function defined in a
neighborhood of the spectrum of $U_{N}$. Note that
$Z^{*}Z=Q^{*}V^{*}VQ=\mathbb{I}_{r}$.

The relation $U_{N}Z=ZT$ can be illustrated as

$\left(\begin{smallmatrix}\square&\square&\square&\square\cr\square&\square&\square&\square\cr\square&\square&\square&\square\cr\square&\square&\square&\square\cr\end{smallmatrix}\right)\left(\begin{smallmatrix}\blacksquare&\blacksquare\cr\blacksquare&\blacksquare\cr\blacksquare&\blacksquare\cr\blacksquare&\blacksquare\cr\end{smallmatrix}\right)=\left(\begin{smallmatrix}\blacksquare&\blacksquare\cr\blacksquare&\blacksquare\cr\blacksquare&\blacksquare\cr\blacksquare&\blacksquare\cr\end{smallmatrix}\right)\left(\begin{smallmatrix}\boxtimes&\boxtimes\cr\mathbf{0}&\boxtimes\end{smallmatrix}\right),$

which explains the term _partial Schur decomposition_ of $U_{N}$ – the
triangular factor $T$ is of smaller dimension than $U_{N}$ and only a portion
of the eigenvalues is computed.

####  4.1.1 Schur functions of $\mathbb{U}$

On the operator level, the relation (30) becomes

$\displaystyle\\!\\!\\!\\!\\!\mathbb{U}(\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}Z)$
$\displaystyle=$
$\displaystyle\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}U_{N}Z+{\boldsymbol{\rho}}(x)^{T}Z$
(31) $\displaystyle=$
$\displaystyle\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}ZT+{\boldsymbol{\rho}}(x)^{T}Z$
$\displaystyle\backsimeq$
$\displaystyle\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}ZT.$
(32)

If we define a new sequence of functions as

$\begin{pmatrix}\zeta_{1}(x)&\zeta_{2}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}=\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}Z,$
(33)

then we have, using (31),

$\displaystyle\\!\\!\\!\\!\mathbb{U}\begin{pmatrix}\zeta_{1}(x)&\zeta_{2}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}$
$\displaystyle=$
$\displaystyle\begin{pmatrix}\zeta_{1}(x)&\zeta_{2}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}T+{\boldsymbol{\rho}}(x)^{T}Z$
(34) $\displaystyle\backsimeq$
$\displaystyle\begin{pmatrix}\zeta_{1}(x)&\zeta_{2}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}T.$

Since $T$ is upper triangular, relation (34) contains a nested sequence of
analogous partial triangulations

$\mathbb{U}\begin{pmatrix}\zeta_{1}&\zeta_{2}&\ldots&\zeta_{i}\end{pmatrix}\backsimeq\begin{pmatrix}\zeta_{1}&\zeta_{2}&\ldots&\zeta_{i}\end{pmatrix}T(1:i,1:i),\;\;i=1,\ldots,r.$
(35)

In particular,
$\mathbb{U}\;\mathrm{Span}\\!\begin{pmatrix}\zeta_{1}&\zeta_{2}&\ldots&\zeta_{i}\end{pmatrix}\precsim\mathrm{Span}\\!\begin{pmatrix}\zeta_{1}&\zeta_{2}&\ldots&\zeta_{i}\end{pmatrix},\;\;i=1,\ldots,r,$
where the relation $\precsim$ is used as a shorthand to be interpreted as
follows: any function in the subspace on the left hand side can be written as
a linear combination of the $\zeta_{i}$’s plus a residual that has certain
minimality property over the available data (See Section 3.1.).

###### Remark 4.2

If the basis functions ${\boldsymbol{\psi}}_{i}$ are an orthonormal set,888If
the basis function are not normalized (but are mutually orthogonal), we simply
redefine them to
${\boldsymbol{\psi}}_{i}\leftarrow{\boldsymbol{\psi}}_{i}/\|{\boldsymbol{\psi}}_{i}\|_{L^{2}}$,
to enforce
$\langle{\boldsymbol{\psi}}_{i},{\boldsymbol{\psi}}_{j}\rangle_{L^{2}}={\boldsymbol{\delta}}_{ij}$.
$\langle{\boldsymbol{\psi}}_{i},{\boldsymbol{\psi}}_{j}\rangle_{L^{2}}={\boldsymbol{\delta}}_{ij}$,
then

$\langle\zeta_{i},\zeta_{j}\rangle_{L^{2}}=\langle\sum_{k}(z_{i})_{k}{\boldsymbol{\psi}}_{k},\sum_{\ell}(z_{j})_{\ell}{\boldsymbol{\psi}}_{\ell}\rangle_{L^{2}}=\sum_{k}(z_{i})_{k}\overline{(z_{j})_{k}}=\langle
z_{i},z_{j}\rangle={\boldsymbol{\delta}}_{ij}.$ (36)

Otherwise, we can apply the Gram-Schmidt orthogonalization (in $L^{2}$)

$\begin{pmatrix}\widehat{\zeta}_{1}&\widehat{\zeta}_{2}&\ldots&\widehat{\zeta}_{r}\end{pmatrix}=\begin{pmatrix}\zeta_{1}&\zeta_{2}&\ldots&\zeta_{r}\end{pmatrix}R,$

where
$\langle\widehat{\zeta}_{i},\widehat{\zeta}_{j}\rangle_{L^{2}}={\boldsymbol{\delta}}_{ij}$,
and $R$ is upper triangular. Then, $\widehat{T}=R^{-1}TR$ is upper triangular
and we have

$\displaystyle\\!\\!\\!\\!\mathbb{U}\begin{pmatrix}\widehat{\zeta}_{1}(x)&\widehat{\zeta}_{2}(x)&\ldots&\widehat{\zeta}_{r}(x)\end{pmatrix}$
$\displaystyle=$
$\displaystyle\begin{pmatrix}\widehat{\zeta}_{1}(x)&\widehat{\zeta}_{2}(x)&\ldots&\widehat{\zeta}_{r}(x)\end{pmatrix}\widehat{T}+{\boldsymbol{\rho}}(x)^{T}VQR$
(37) $\displaystyle\backsimeq$
$\displaystyle\begin{pmatrix}\widehat{\zeta}_{1}(x)&\widehat{\zeta}_{2}(x)&\ldots&\widehat{\zeta}_{r}(x)\end{pmatrix}\widehat{T}.$

$\boxtimes$

###  4.2 Ordered partial Schur form

The partial Schur form can be truncated at any index $\ell<r$, resulting again
in a partial Schur form of lower order, see (35). It may be of interest to
have such a low order partial Schur form that contains only some selected
eigenvalues that are not necessarily the leading $i$ values on the diagonal of
$T$. On the other hand, when we compute the Schur form numerically, the order
of the eigenvalues on the diagonal of $T$ is in general not known in advance –
it depends on the convergence of the numerical algorithm/software for a given
input $\widehat{U}$.

A given Schur form can be reordered in the sense that an efficient algorithm
determines a unitary matrix $\Theta$ such that
$\widetilde{T}=\Theta^{*}T\Theta$ is again upper triangular with diagonal
entries corresponding to the eigenvalues in any given order; see e.g. [5]. The
new partial Schur form of $U_{N}$ becomes

$U_{N}(VQ\Theta)=(VQ\Theta)\widetilde{T},\;\;\widetilde{T}=\Theta^{*}T\Theta,$

and we replace (34) with

$\mathbb{U}\begin{pmatrix}\zeta_{1}(x)&\zeta_{2}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}\Theta\backsimeq\begin{pmatrix}\zeta_{1}(x)&\zeta_{2}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}\Theta(\Theta^{*}T\Theta),$
(38)

i.e. the new functions are generated using $\widetilde{Z}=Z\Theta=VQ\Theta$,

$\begin{pmatrix}\widetilde{\zeta}_{1}(x)&\ldots&\widetilde{\zeta}_{r}(x)\end{pmatrix}=\begin{pmatrix}\zeta_{1}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}\Theta=\begin{pmatrix}{\boldsymbol{\psi}}_{1}(x)&\ldots&{\boldsymbol{\psi}}_{N}(x)\end{pmatrix}\widetilde{Z}.$
(39)

Since $\widetilde{Z}^{*}\widetilde{Z}=\mathbb{I}_{r}$, the orthogonality
similar to (36) remains true. The new partial Schur decomposition
$U_{N}\widetilde{Z}=\widetilde{Z}\widetilde{T}$ can be truncated as in (35) at
some index $i$, with the desired eigenvalues on the diagonal of
$\widetilde{T}(1:i,1:i)$, and the corresponding Schur functions
$\widetilde{\zeta}_{1},\ldots,\widetilde{\zeta}_{i}$.

###### Remark 4.3

It should be noted that the reordering of the Schur form works on the
triangular factor that is of modest dimension, and computing the additional
reordering transformation $\Theta$ is not prohibitively expensive. LAPACK
library provides functions STRSEN, DTRSEN, CTRSEN, ZTRSEN that implement Schur
form reordering for real and complex matrices, single and double precision.
Its application to $Q$ or $VQ$ is a BLAS 3 operation that can be executed very
efficiently, or (depending on the application of the decomposition) the Schur
vectors can be kept in factored form. Matlab package contains the function
ordschur that implements a method by Kressner [33].  $\boxtimes$

###  4.3 Evaluation of the Schur functions

The kernel trick evaluation of the functions $\zeta_{i}(\cdot)$ follows the
lines of Section 3.5. At an arbitrary $\mathbf{x}\in\mathbb{C}^{n}$ we have
(see Section 3.5.2)
$\zeta_{i}(x)={\boldsymbol{\psi}}(x)^{T}VQe_{i}={\boldsymbol{\psi}}(x)^{T}{\boldsymbol{\Psi}}_{x}^{*}W\Sigma^{-1}Qe_{i}$
and thus

$\begin{pmatrix}\zeta_{1}(\mathbf{x})&\ldots&\zeta_{r}(\mathbf{x})\end{pmatrix}=\begin{pmatrix}{\boldsymbol{f}}(\mathbf{x},\mathbf{x}_{1})&{\boldsymbol{f}}(\mathbf{x},\mathbf{x}_{2})&\ldots&{\boldsymbol{f}}(\mathbf{x},\mathbf{x}_{m})\end{pmatrix}W\Sigma^{-1}Q.$

Note that a kernel trick can be used to compute $\zeta_{i}(\mathbf{x})$. At
the snapshots, we have the tabulated values computed as

${\boldsymbol{\zeta}}_{x}=\begin{pmatrix}\zeta_{1}(\mathbf{x}_{1})&\ldots&\zeta_{r}(\mathbf{x}_{1})\cr\vdots&\ddots&\vdots\cr\zeta_{1}(\mathbf{x}_{m})&\ldots&\zeta_{r}(\mathbf{x}_{m})\end{pmatrix}={\boldsymbol{\Psi}}_{x}VQ=W\Sigma^{2}W^{*}W\Sigma^{-1}Q=W\Sigma
Q\in\mathbb{C}^{m\times r}.$ (40)

Note that ${\boldsymbol{\zeta}}_{x}$ is conveniently given factored in a SVD
decomposition.

###  4.4 Representation of observables using Schur functions

As in Section 3.6, we seek representation of given scalar functions
$g_{1},\ldots,g_{\ell}$, but this time in terms of the Schur functions
$\zeta_{1},\ldots,\zeta_{r}$,

$\begin{pmatrix}g_{1}(x)&\ldots&g_{\ell}(x)\end{pmatrix}\backsimeq\begin{pmatrix}\zeta_{1}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}\Xi,\;\;\Xi\in\mathbb{C}^{r\times\ell}.$
(41)

In the restrictive data driven framework, the representation (41) is feasible
only at the data snapshots $\mathbf{x}_{j}$, i.e. at the data matrix
$G_{ji}=g_{i}(\mathbf{x}_{j})$. The approximation is computed in the least
squares sense, column-wise, e.g.

$\Xi(:,i)=\mathrm{argmin}_{h\in\mathbb{C}^{r}}\|\left(\begin{smallmatrix}g_{i}(\mathbf{x}_{1})\cr\vdots\cr
g_{i}(\mathbf{x}_{m})\end{smallmatrix}\right)-{\boldsymbol{\zeta}}_{x}h\|_{2}={\boldsymbol{\zeta}}_{x}^{\dagger}\left(\begin{smallmatrix}g_{i}(\mathbf{x}_{1})\cr\vdots\cr
g_{i}(\mathbf{x}_{m})\end{smallmatrix}\right),\;\;i=1,\ldots,\ell.$

Hence, $G\approx{\boldsymbol{\zeta}}_{x}\Xi$ with

$\Xi={\boldsymbol{\zeta}}_{x}^{\dagger}G=\begin{pmatrix}\zeta_{1}(\mathbf{x}_{1})&\ldots&\zeta_{r}(\mathbf{x}_{1})\cr\vdots&\ddots&\vdots\cr\zeta_{1}(\mathbf{x}_{m})&\ldots&\zeta_{r}(\mathbf{x}_{m})\end{pmatrix}^{\dagger}\begin{pmatrix}g_{1}(\mathbf{x}_{1})&\ldots&g_{\ell}(\mathbf{x}_{1})\cr\vdots&\ddots&\vdots\cr
g_{1}(\mathbf{x}_{m})&\ldots&g_{\ell}(\mathbf{x}_{m})\end{pmatrix}=Q^{*}\Sigma^{-1}W^{*}G,$

where the pseudoinverse ${\boldsymbol{\zeta}}_{x}^{\dagger}$ is computed
directly from the readily available SVD (40) of ${\boldsymbol{\zeta}}_{x}$.

In a DMD framework, the data snapshots are usually stored column-wise in a
matrix, so that the column index represents the discrete time. For such a
representation, write $G^{T}\approx\Xi^{T}{\boldsymbol{\zeta}}_{x}^{T}$,
introduce column partition
$\Xi^{T}=\begin{pmatrix}\eta_{1}&\ldots&\eta_{r}\end{pmatrix}$ and write the
representation (41) at the snapshots $\mathbf{x}_{i}$ as

$\begin{pmatrix}g_{1}(\mathbf{x}_{1})&g_{1}(\mathbf{x}_{2})&\ldots&g_{1}(\mathbf{x}_{m})\cr
g_{2}(\mathbf{x}_{1})&g_{2}(\mathbf{x}_{2})&\ldots&g_{2}(\mathbf{x}_{m})\cr\vdots&\vdots&\cdots&\vdots\cr
g_{\ell}(\mathbf{x}_{1})&g_{\ell}(\mathbf{x}_{2})&\ldots&g_{\ell}(\mathbf{x}_{m})\cr\end{pmatrix}\approx\begin{pmatrix}\eta_{1}&\ldots&\eta_{r}\end{pmatrix}\\!\\!\begin{pmatrix}\zeta_{1}(\mathbf{x}_{1})&\zeta_{1}(\mathbf{x}_{2})&\ldots&\zeta_{1}(\mathbf{x}_{m})\cr\zeta_{2}(\mathbf{x}_{1})&\zeta_{2}(\mathbf{x}_{2})&\ldots&\zeta_{2}(\mathbf{x}_{m})\cr\vdots&\vdots&\cdots&\vdots\cr\zeta_{r}(\mathbf{x}_{1})&\zeta_{r}(\mathbf{x}_{2})&\ldots&\zeta_{r}(\mathbf{x}_{m})\cr\end{pmatrix}.$
(42)

The special case $G^{T}=\mathbf{X}$ (the case of full state observables) is
described in the next two propositions.

###### Proposition 4.1

For the data snapshots representation, we specify the observables $g_{i}$ as:
$\ell=n$, $g_{i}(x)=\mathfrak{e}_{i}^{T}x$. Then $G=\mathbf{X}^{T}$, and we
have

$\Xi={\boldsymbol{\zeta}}_{x}^{\dagger}\mathbf{X}^{T}=(W\Sigma
Q)^{\dagger}\mathbf{X}^{T}=Q^{*}\Sigma^{-1}W^{*}\mathbf{X}^{T},$

i.e.

$\mathbf{X}\approx\Xi^{T}{\boldsymbol{\zeta}}_{x}^{T},\;\;\mbox{where}\;\;\Xi^{T}=\mathbf{X}W^{*T}\Sigma^{-1}Q^{*T},\;\;{\boldsymbol{\zeta}}_{x}^{T}=Q^{T}\Sigma
W^{T},$ (43)

and, with the column partition
$\Xi^{T}=\begin{pmatrix}\eta_{1}&\ldots&\eta_{r}\end{pmatrix}$,

$\mathbf{x}_{i}\approx\sum_{j=1}^{r}\eta_{j}\zeta_{j}(\mathbf{x}_{i}).$

This representation of the data snapshot matrix $\mathbf{X}$ simplifies if the
basis functions are ${\boldsymbol{\psi}}_{i}(x)=\mathfrak{e}_{i}^{T}x$, i.e.
when ${\boldsymbol{\Psi}}_{x}=\mathbf{X}^{T}$.

###### Proposition 4.2

In the case ${\boldsymbol{\Psi}}_{x}=\mathbf{X}^{T}=W\Sigma V^{*}$ and
${\boldsymbol{f}}(x,y)=\langle x,y\rangle_{\mathbb{C}^{n}}$, the
representation (43) of the data snapshot matrix $\mathbf{X}$ in Proposition
4.1 holds with

$\Xi^{T}=\mathbf{X}W^{*T}\Sigma^{-1}Q^{*T}=V^{*T}Q^{*T}=Z^{*T}=\begin{pmatrix}\eta_{1}&\ldots&\eta_{r}\end{pmatrix}.$

Note that in this case $\Xi^{T}=Z^{*T}$ is orthonormal, where $Z=VQ$ is the
orthonormal matrix from the partial Schur decomposition (29). Further,
${\boldsymbol{\zeta}}_{x}^{T}=Z^{T}\mathbf{X}$.
(If the data matrices are scaled as in [15], then the scaling must be
accordingly applied to ${\boldsymbol{\zeta}}_{x}^{T}$.)

For the readers convenience, the case of canonical basis functions (the DMD
framework) is summarized in Algorithm 1. (Note that the algorithm is written
in an EDMD framework, to express the DMD as a special case of the EDMD. A
usual DMD-style formulation works on $\mathbf{X}$ and $\mathbf{Y}$, and the
Rayleigh quotient is $\widehat{U}^{T}$.)

Algorithm 1
$(T,Z,{\boldsymbol{\zeta}}_{x}^{T})=\textsf{Koopman\\_Schur\\_SSMD}(\mathbf{X},\mathbf{Y})$

1: Compute the SVD $\mathbf{X}^{T}=W\Sigma V^{*}$ {This step includes
numerical rank determination.}

2: $\widehat{U}=\Sigma^{-1}W^{*}(\mathbf{Y}^{T}V)$ {The Rayleigh quotient. See
Section 3.4.1}

3: Compute the Schur decomposition $\widehat{U}=QTQ^{*}$ ;

4: $Z=VQ$ ; {The Schur vectors.}

5: ${\boldsymbol{\zeta}}_{x}^{T}=Q^{T}\Sigma W^{T}$ {The Schur functions
values tabulated at the snapshots (40).}

####  4.4.1 Snapshot representation using a subset of modes

High fidelity representation of the data snapshots using only a selected
subset of the computed modes allows for a model order reduction and revealing
of underlying structures in the dynamics (e.g. identification and tracking of
coherent structures in a fluid flow). The sparsity promoting DMD [30], [31] is
such a technique.

Our goal is to develop a similar snapshot reconstruction scheme, but for the
Koopman-Schur framework. We start with a general result on a structured least
squares problem that generalizes a Khatri-Rao product based formulation of the
snapshot reconstruction problem [20].

###### Theorem 4.4

Let $\mathbf{X}\in\mathbb{C}^{n\times m}$, $B\in\mathbb{C}^{n\times\ell}$,
$C\in\mathbb{C}^{m\times\ell}$, ${\boldsymbol{\Omega}}\in\mathbb{C}^{m\times
m}$. If the Khatri-Rao product $C\odot B$ is of full column rank and
${\boldsymbol{\Omega}}$ is nonsingular, then the optimal solution
${\boldsymbol{\alpha}}=(\alpha_{1},\ldots,\alpha_{\ell})$ of the least squares
problem

$\|(\mathbf{X}-B\left(\begin{smallmatrix}\alpha_{1}&&\cr&\cdot&\cr&&\alpha_{\ell}\end{smallmatrix}\right)C^{T}){\boldsymbol{\Omega}}\|_{F}\longrightarrow\min_{\alpha_{1},\ldots,\alpha_{\ell}}$
(44)

is unique and it can be compactly written as

${\boldsymbol{\alpha}}=[(B^{*}B)\circ{(C^{*}{\boldsymbol{\Omega}}^{T*}{\boldsymbol{\Omega}}^{T}C)}]^{-1}[{C^{*}}\circ(B^{*}\mathbf{X}{\boldsymbol{\Omega}}{\boldsymbol{\Omega}}^{*})]\vec{\mathbf{1}},\;\;\vec{\mathbf{1}}=(1,\ldots,1)^{T}\in\mathbb{R}^{m},$
(45)

where $\circ$ denotes the Hadamard matrix product.

Proof: The linear least squares problem (44) can be rewritten into

$vec((\mathbf{X}-B\left(\begin{smallmatrix}\alpha_{1}&&\cr&\cdot&\cr&&\alpha_{\ell}\end{smallmatrix}\right)C^{T}){\boldsymbol{\Omega}})=({\boldsymbol{\Omega}}^{T}\otimes\mathbb{I}_{n})(vec(\mathbf{X})-(C\odot
B)\left(\begin{smallmatrix}\alpha_{1}\cr\vdots\cr\alpha_{\ell}\end{smallmatrix}\right)),$

where $\otimes$ and $\odot$ denote the Kronecker and the Khatri-Rao matrix
products, respectively, and $vec(\cdot)$ is the column-wise reshaping
isomorphism between the vector spaces $\mathbb{C}^{n\times m}$ and
$\mathbb{C}^{n\cdot m}$. Further, if we define

$S_{\boldsymbol{\omega}}=({\boldsymbol{\Omega}}^{T}\otimes\mathbb{I}_{n})(C\odot
B)=({\boldsymbol{\Omega}}^{T}C)\odot B,$

then the minimization problem (44) reads

$\|S_{w}\left(\begin{smallmatrix}\alpha_{1}\cr\vdots\cr\alpha_{\ell}\end{smallmatrix}\right)-({\boldsymbol{\Omega}}^{T}\otimes\mathbb{I}_{n})vec(\mathbf{X})\|_{F}\longrightarrow\min_{\alpha_{1},\ldots,\alpha_{\ell}}.$

The assumptions of the theorem guarantee that $S_{\boldsymbol{\omega}}$ has
full column rank and the unique solution can be expressed using the Moore-
Penrose pseudo-inverse as
${\boldsymbol{\alpha}}=S_{\boldsymbol{\omega}}^{\dagger}({\boldsymbol{\Omega}}^{T}\otimes\mathbb{I}_{n})vec(\mathbf{X})=(S_{\boldsymbol{\omega}}^{*}S_{\boldsymbol{\omega}})^{-1}S_{\boldsymbol{\omega}}^{*}({\boldsymbol{\Omega}}^{T}\otimes\mathbb{I}_{n})vec(\mathbf{X})$.
The factors in this formula can be expressed, using the well known relations
between the Hadamard, the Kronecker and the Khatri-Rao matrix products, as

$\displaystyle S_{\boldsymbol{\omega}}^{*}S_{\boldsymbol{\omega}}$
$\displaystyle=$ $\displaystyle(({\boldsymbol{\Omega}}^{T}C)\odot
B)^{*}(({\boldsymbol{\Omega}}^{T}C)\odot
B)=[({\boldsymbol{\Omega}}^{T}C)^{*}({\boldsymbol{\Omega}}^{T}C)]\circ(B^{*}B)$
(46) $\displaystyle=$
$\displaystyle{(C^{*}{\boldsymbol{\Omega}}^{T*}{\boldsymbol{\Omega}}^{T}C)}\circ(B^{*}B)=(B^{*}B)\circ{(C^{*}{\boldsymbol{\Omega}}^{T*}{\boldsymbol{\Omega}}^{T}C)}$
$\displaystyle
S_{\boldsymbol{\omega}}^{*}({\boldsymbol{\Omega}}^{T}\otimes\mathbb{I}_{n})vec(\mathbf{X})$
$\displaystyle=$ $\displaystyle(C\odot
B)^{*}({\boldsymbol{\Omega}}^{T*}{\boldsymbol{\Omega}}^{T}\otimes\mathbb{I}_{n})vec(\mathbf{X})=[{C^{*}}\circ(B^{*}X{\boldsymbol{\Omega}}{\boldsymbol{\Omega}}^{*})]\vec{\mathbf{1}}$

which yields (45). $\boxtimes$

###### Remark 4.5

The condition that $C\odot B$ is of full column rank is satisfied if e.g. one
of the two matrices is of full column rank and the other one has no zero
column. This is a consequence of a Schur product theorem, see e.g. [29,
Section 5.2].  $\boxtimes$

###### Remark 4.6

Clearly, with an additional nonsigular matrix ${\boldsymbol{\Omega}}_{1}$, the
problem can be formulated as

$\|{\boldsymbol{\Omega}}_{1}(\mathbf{X}-B\left(\begin{smallmatrix}\alpha_{1}&&\cr&\cdot&\cr&&\alpha_{\ell}\end{smallmatrix}\right)C^{T}){\boldsymbol{\Omega}}\|_{F}\longrightarrow\min_{\alpha_{1},\ldots,\alpha_{\ell}}$
(47)

with the solution
${\boldsymbol{\alpha}}=[(B^{*}{\boldsymbol{\Omega}}_{1}^{*}{\boldsymbol{\Omega}}_{1}B)\circ{(C^{*}{\boldsymbol{\Omega}}^{T*}{\boldsymbol{\Omega}}^{T}C)}]^{-1}[{C^{*}}\circ(B^{*}{\boldsymbol{\Omega}}_{1}^{*}{\boldsymbol{\Omega}}_{1}\mathbf{X}{\boldsymbol{\Omega}}{\boldsymbol{\Omega}}^{*})]\vec{\mathbf{1}}.$
In this case ${\boldsymbol{\Omega}}_{1}$ can be e.g. the Cholesky factor of
the inverse noise covariance matrix. Throughout this paper, for the sake of
brevity, we have ${\boldsymbol{\Omega}}_{1}=\mathbb{I}_{n}$. See [20] for more
details.  $\boxtimes$

The formula (45) is in particular interesting in our case where
$C={\boldsymbol{\zeta}}_{x}$ and $B=\Xi^{T}=Z^{*T}$ is orthonormal and thus
$B^{*}B=\mathbb{I}_{\ell}$.

###### Corollary 4.7

If in Theorem 4.4 the matrix $B$ is orthonormal ($B^{*}B=\mathbb{I}_{\ell}$),
then the solution (45) of the least squares problem (44) reads

${\boldsymbol{\alpha}}=\mathrm{diag}(1/\|{\boldsymbol{\Omega}}^{T}C(:,j)\|_{2}^{2})_{j=1}^{\ell}[{C^{*}}\circ(B^{*}\mathbf{X}{\boldsymbol{\Omega}}{\boldsymbol{\Omega}}^{*})]\vec{\mathbf{1}}.$

###### Remark 4.8

The coefficients in Corollary 4.7 are actually the Fourier coefficients and
the formula for ${\boldsymbol{\alpha}}$ follows directly from the projection
theorem. Namely, if $b_{i}$ and $c_{i}$ denote the $i$th column of $B$ and
$C$, respectively, then the approximation problem reads

$\mathbf{X}{\boldsymbol{\Omega}}\approx\sum_{i=1}^{\ell}\alpha_{i}b_{i}c_{i}^{T}{\boldsymbol{\Omega}},\;\;\mbox{where}\;\;\langle
b_{i}c_{i}^{T}{\boldsymbol{\Omega}},b_{j}c_{j}^{T}{\boldsymbol{\Omega}}\rangle_{F}={\boldsymbol{\delta}}_{ij}\|b_{i}\|_{2}^{2}\|{\boldsymbol{\Omega}}^{T}c_{i}\|_{2}^{2}.$

Here ${\boldsymbol{\delta}}_{ij}$ is the Kroneecker symbol and
$\langle\cdot,\cdot\rangle_{F}$ is the Frobenius inner product. The optimal
coefficients are thus

$\alpha_{i}\\!=\\!\frac{\langle\mathbf{X}{\boldsymbol{\Omega}},b_{i}c_{i}^{T}{\boldsymbol{\Omega}}\rangle_{F}}{\langle
b_{i}c_{i}^{T}{\boldsymbol{\Omega}},b_{i}c_{i}^{T}{\boldsymbol{\Omega}}\rangle_{F}}\\!=\\!\frac{\mathrm{Trace}({\boldsymbol{\Omega}}^{*}c_{i}^{T*}b_{i}^{*}\mathbf{X}{\boldsymbol{\Omega}})}{\|b_{i}\|_{2}^{2}\|{\boldsymbol{\Omega}}^{T}c_{i}\|_{2}^{2}}\\!=\\!\frac{\mathrm{Trace}(b_{i}^{*}\mathbf{X}{\boldsymbol{\Omega}}{\boldsymbol{\Omega}}^{*}c_{i}^{T*})}{\|{\boldsymbol{\Omega}}^{T}c_{i}\|_{2}^{2}}\\!=\\!\frac{[c_{i}^{*}\\!\circ\\!(b_{i}^{*}\mathbf{X}{\boldsymbol{\Omega}}{\boldsymbol{\Omega}}^{*})]\vec{\mathbf{1}}}{\|{\boldsymbol{\Omega}}^{T}c_{i}\|_{2}^{2}},$

precisely as stated in Corollary 4.7.  $\boxtimes$

The role of the matrix ${\boldsymbol{\Omega}}$ can be to weigh the importance
of individual snapshots in the overall reconstruction error. This is achieved
by choosing ${\boldsymbol{\Omega}}$ as diagonal positive definite.999Actually,
${\boldsymbol{\Omega}}$ can be semidefinite, where
${\boldsymbol{\omega}}_{i}=0$ means that the $i$th snapshot is simply
excluded. After removing all zero-weighted snapshots, we proceed with the
problem with positive weights.

###### Corollary 4.9

Let
$B=\Xi^{T}(:,[i_{1},\ldots,i_{\ell}])=(\eta_{i_{1}},\ldots,\eta_{i_{\ell}})$,
$C={\boldsymbol{\zeta}}_{x}(:,[i_{1},\ldots,i_{\ell}])$, and
${\boldsymbol{\Omega}}=\mathrm{diag}({\boldsymbol{\omega}}_{i})_{i=1}^{m}$,
${\boldsymbol{\omega}}_{i}>0$ for all $i=1,\ldots,m$. Then the optimal
reconstruction $\mathbf{X}\approx B{\boldsymbol{\alpha}}C^{T}$ in the sense of
Theorem 4.4 is obtained with

${\boldsymbol{\alpha}}=\mathrm{diag}(\frac{1}{\sum_{k=1}^{m}{\boldsymbol{\omega}}_{k}^{2}|\zeta_{i_{j}}(\mathbf{x}_{k})|^{2}})_{j=1}^{\ell}\;[({\underbrace{(\overline{\zeta_{i_{j}}(\mathbf{x}_{k})})_{j,k=1,1}^{\ell,m}}_{C^{*}}})\circ(Q(:,[i_{1},\ldots,i_{\ell}])^{T}V^{T}\mathbf{X}{\boldsymbol{\Omega}}^{2})]\vec{\mathbf{1}}.$

The computational complexity of the above formula (given its ingredients) is
independent of the state space dimension $n$. In a Koopman-Schur-ESSMD
setting, the kernel trick is used, as outlined in Section 4.3. If
${\boldsymbol{\Psi}}_{x}=\mathbf{X}^{T}$, as in Proposition 4.2, then
${\boldsymbol{\alpha}}=(1,\ldots,1)^{T}$.

Proof: Note that
$B^{*}\mathbf{X}=(Z^{T*}(:,[i_{1},\ldots,i_{\ell}]))^{*}X=Q(:,[i_{1},\ldots,i_{\ell}])^{T}V^{T}\mathbf{X}.$
If ${\boldsymbol{\Psi}}_{x}=\mathbf{X}^{T}$, then we have (from the SVD of
${\boldsymbol{\Psi}}_{x}$) the SVD of $\mathbf{X}$ as well, and

$B^{*}\mathbf{X}=Q(:,[i_{1},\ldots,i_{\ell}])^{T}V^{T}V^{*T}\Sigma
W^{T}=Q(:,[i_{1},\ldots,i_{\ell}])^{T}\Sigma
W^{T}={\boldsymbol{\zeta}}_{x}(:,[i_{1},\ldots,i_{\ell}])^{T}.$

Since $B^{*}B=\mathbb{I}_{\ell}$, Corollary 4.7 yields

${\boldsymbol{\alpha}}=\mathrm{diag}(\frac{1}{\sum_{k=1}^{m}{\boldsymbol{\omega}}_{k}^{2}|\zeta_{i_{j}}(\mathbf{x}_{k})|^{2}})_{j=1}^{\ell}\;[({\underbrace{(\overline{\zeta_{i_{j}}(\mathbf{x}_{k})})_{j,k=1,1}^{\ell,m}}_{C^{*}}})\circ(\underbrace{Q(:,[i_{1},\ldots,i_{\ell}])^{T}\Sigma
W^{T}}_{C^{T}}{\boldsymbol{\Omega}}^{2})]\vec{\mathbf{1}}.$

It is easily checked that the above expression actually reads
${\boldsymbol{\alpha}}=(1,\ldots,1)^{T}$. $\boxtimes$

####  4.4.2 Reordering trick

A reduced order representation of the snapshots may be based on a selection of
eigenvalues
$\lambda_{i_{1}}=T_{i_{1}i_{1}},\ldots,\lambda_{i_{\ell}}=T_{i_{\ell}i_{\ell}}$.
In the standard DMD/EDMD approaches, one simply takes the corresponding modes
and solves a least squares problem for the coefficients. In the Koopman-Schur
setting, the basis for the representation is an orthonormal basis that spans
the spectral subspace encompassing the selected eigenvalues. In a Schur form,
such a basis is in the leading columns of the unitary factor, provided the
eigenvalues are at the corresponding leading positions (upper left corner) of
the triangular factor; the decomposition is then truncated and only selected
values are used in the next steps. Another situation where reordering is
needed is when spurious eigenvalues appear (e.g. that are far outside unit
circle) – reordering can move them to the bottom right part of the triangular
form and then truncated from the decomposition.

The reordering trick, that enables orthonormal representation with any subset
of the eigenvalues (Section 4.4.1) is to reorder the Schur form (Section 4.2)
so that in the new triangular Schur factor $\widetilde{T}$ we have
$\widetilde{T}_{11}=\lambda_{i_{1}},\ldots,\widetilde{T}_{\ell\ell}=\lambda_{i_{\ell}}$,
and then to truncate it at index $\ell$.

###  4.5 Using the Schur functions for forecasting

The Koopman mode decomposition is a powerful tool for revealing
coherent/latent structure of the dynamics under study, model order reduction
and for forecasting the future values of the observables of interest. For the
latter, the key is the accuracy of the numerically computed eigenfunctions and
the eigenvalues.

The proposed Koopman-Schur decomposition can also be adapted for forecasting.
The forecasting mechanism works as follows. Once we have the representation
(41), the relation (34) implies

$\begin{pmatrix}(\mathbb{U}g_{1})(x)&\ldots&(\mathbb{U}g_{\ell})(x)\end{pmatrix}\backsimeq\begin{pmatrix}\zeta_{1}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}T\Xi,$

and, inductively,

$\begin{pmatrix}(\mathbb{U}^{k}g_{1})(x)&\ldots&(\mathbb{U}^{k}g_{\ell})(x)\end{pmatrix}\backsimeq\begin{pmatrix}\zeta_{1}(x)&\ldots&\zeta_{r}(x)\end{pmatrix}T^{k}\Xi,\;\;k=1,2,\ldots.$
(48)

Hence, pushing the observables $g_{i}$ forward
($\mathbb{U}^{k}g_{i}=g_{i}\circ\mathbf{F}\circ\cdots\circ\mathbf{F}$) is, in
the representation (41), realized by raising the power of the triangular
factor $T$. This includes the forecasting, i.e. extrapolating the
representation of the present snapshot into the future. In terms of the values
at the snapshots (see (42)) the relations (48) read

$((\mathbb{U}^{k}g_{i})(\mathbf{x}_{j}))_{i,j=1}^{\ell,m}\approx\Xi^{T}(T^{T})^{k}{\boldsymbol{\zeta}}_{x}^{T},\;\;k=1,2,\ldots.$
(49)

Of course, with repeated applications of $\mathbb{U}$, the residual will build
up (see (31)) and the approximations (48) will deteriorate in the sense of
$\backsimeq$ . Nevertheless, we can use (49) to approximate, for
$j=2,3,\ldots$

$\Xi^{T}\\!\begin{pmatrix}\zeta_{1}(\mathbf{x}_{j})\cr\zeta_{2}(\mathbf{x}_{j})\cr\vdots\cr\zeta_{r}(\mathbf{x}_{j})\end{pmatrix}\\!\\!\approx\\!\\!\begin{pmatrix}g_{1}(\mathbf{x}_{j})\cr
g_{2}(\mathbf{x}_{j})\cr\vdots\cr
g_{\ell}(\mathbf{x}_{j})\end{pmatrix}\\!\\!=\\!\\!\begin{pmatrix}(g_{1}\circ\mathbf{F})(\mathbf{x}_{j-1})\cr(g_{2}\circ\mathbf{F})(\mathbf{x}_{j-1})\cr\vdots\cr(g_{\ell}\circ\mathbf{F})(\mathbf{x}_{j-1})\end{pmatrix}\\!\\!=\\!\\!\begin{pmatrix}(\mathbb{U}g_{1})(\mathbf{x}_{j-1})\cr(\mathbb{U}g_{2})(\mathbf{x}_{j-1})\cr\vdots\cr(\mathbb{U}g_{\ell})(\mathbf{x}_{j-1})\end{pmatrix}\\!\approx\\!\Xi^{T}T^{T}\\!\begin{pmatrix}\zeta_{1}(\mathbf{x}_{j-1})\cr\zeta_{2}(\mathbf{x}_{j-1})\cr\vdots\cr\zeta_{r}(\mathbf{x}_{j-1})\end{pmatrix},$

i.e.

$\begin{pmatrix}\zeta_{1}(\mathbf{x}_{j})\cr\zeta_{2}(\mathbf{x}_{j})\cr\vdots\cr\zeta_{r}(\mathbf{x}_{j})\end{pmatrix}\approx
T^{T}\\!\begin{pmatrix}\zeta_{1}(\mathbf{x}_{j-1})\cr\zeta_{2}(\mathbf{x}_{j-1})\cr\vdots\cr\zeta_{r}(\mathbf{x}_{j-1})\end{pmatrix},\;\;j=2,3,\ldots,m,m+1.$
(50)

The quality of the approximation $\approx$ in (50) is an indirect indicator of
the success/failure of the entire discretization process from Section 3.1.

Recall that for an arbitrary $x\in\mathbb{C}^{n}$

$\begin{pmatrix}\zeta_{1}(x)\cr\zeta_{2}(x)\cr\vdots\cr\zeta_{r}(x)\end{pmatrix}=Q^{T}\Sigma^{-1}W^{T}\begin{pmatrix}{\boldsymbol{f}}(x,\mathbf{x}_{1})\cr{\boldsymbol{f}}(x,\mathbf{x}_{2})\cr\vdots\cr{\boldsymbol{f}}(x,\mathbf{x}_{m})\end{pmatrix}=Q^{T}\Sigma^{-1}W^{T}\mathbf{X}^{*}x=\underbrace{Q^{T}\Sigma^{-1}W^{T}W^{*T}\Sigma
V^{T}}_{Z^{T}}x,$

where the last two equalities hold in the case
${\boldsymbol{\Psi}}_{x}=\mathbf{X}^{T}$, ${\boldsymbol{f}}(x,y)=\langle
x,y\rangle_{\mathbb{C}^{n}}=y^{*}x$. In particular, for $x=\mathbf{x}_{m+1}$,

$\begin{pmatrix}\zeta_{1}(\mathbf{x}_{m+1})\cr\zeta_{2}(\mathbf{x}_{m+1})\cr\vdots\cr\zeta_{r}(\mathbf{x}_{m+1})\end{pmatrix}\\!=Q^{T}\Sigma^{-1}W^{T}\\!\begin{pmatrix}\langle\mathbf{x}_{m+1},\mathbf{x}_{1}\rangle_{\mathbb{C}^{n}}\cr\langle\mathbf{x}_{m+1},\mathbf{x}_{2}\rangle_{\mathbb{C}^{n}}\cr\vdots\cr\langle\mathbf{x}_{m+1},\mathbf{x}_{m}\rangle_{\mathbb{C}^{n}}\end{pmatrix}=\overline{Z}^{*}\mathbf{x}_{m+1}=Z^{T}\mathbf{x}_{m+1},$

as it should be from the least squares approximation theory (Proposition 4.2).
Note that (50) represents an explicit time stepping that can be deployed for
the purposes of forecasting beyond the index $j=m+1$.

###  4.6 The real Koopman-Schur SSMD

If the data is real, it is advantageous to perform as much as possible of the
computation using only the real computer arithmetic. Since the eigenvalues are
in general complex, we can reformulate the algorithm using the real Schur form
[28, Sec. 2.3]. This means that in the Schur form (28) of a real $\widehat{U}$
the matrix $Q$ is real orthogonal and $T$ is block upper triangular with
diagonal blocks $1\times 1$ for each real eigenvalue and $2\times 2$ for each
complex conjugate pairs of eigenvalues. For instance, this real form is the
default in the Matlab function schur and the complex form is obtained with the
input option ’complex’. The LAPACK subroutines SGEES, DGEES are designed to
compute the real Schur form of real matrices.

Using the real Schur form for real data is recommended and it is implemented
in our software. The required changes (as compared with the complex form) are
merely technical and require most additional effort during software
development. Some care is needed when defining the nested invariant subspaces
(see (35)) in the sense that complex conjugate pairs should always be taken
together, with the corresponding Schur vectors. Further, if a reordering of
the Schur form is needed, the reordering algorithm should be adapted to the
real form. Our LAPACK-based implementation of the Koopman-Schur SSMD method (a
separate work) is designed to use the real Schur form for real data, similarly
as we have recently done for DMD implementations [16], [17].

###  4.7 QR compressed Koopman-Schur SSMD

The QR compressed DMD [15], with implementation details [16, Section 3.7] can
considerably increase the run time efficiency of the DMD in case of high
dimensional data snapshots, and it can be used for fast updates (adding new
data snapshots, single or in batches) and down-dates (discarding oldest data
snapshots, while receiving new ones) for the streaming DMD. The technique is
easily transferred to the Hermitian DMD in [17].

It is possible to apply the QR compression to the KS-SSMD Algorithm 1 and to
formulate a streaming version. We will provide details in a forthcoming work.

##  5 Numerical evaluation

We now illustrate the Koopman-Schur framework using numerical examples.
Software prototypes of the presented algorithms are developed using Matlab and
the numerical examples are selected to check the validity and the
functionality of the new proposed approach. The experiments provide insights
and show differences between the methods, both row-wise (DMD vs. EDMD) and
column-wise (diagonalization vs. Schur form) in Table 1.

| Dynamic Mode Decomp. | Extended Dynamic Mode Decomp.
---|---|---
diagonalization | DMD | EDMD
Schur form | Koopman-Schur-SSMD | Koopman-Schur-ESSDMD
Table 1: The four methods used in the numerical experiments in this section.
In the EDMD-based methods (second column) the kernel-trick formulation is
used.

We use two kernel functions

  1. 1.

$k_{1}(x,y)=y^{*}x=\langle x,y\rangle_{\mathbb{C}^{n}}$. With this choice EDMD
and DMD are theoretically equivalent but EDMD is numerically more difficult
because of higher condition numbers of the matrices used in the algorithm.

  2. 2.

$k_{2}(x,y)=e^{-\frac{\|x-y\|_{2}^{2}}{2\sigma^{2}}}$. The Gaussian kernel is
popular and powerful kernel function, and for good performance proper choice
of the hyper-parameter $\sigma$ is critical. We will not tackle the question
of choosing $\sigma$; instead we have empirically selected a value that works
well on our examples.

###  5.1 What is being tested

The main goal of the experiments presented here is to illustrate the proposed
framework. The themes of the test are as follows:

##### Eigenvalues.

The eigenvalues computed by all four methods are compared. If in some test
case we have the eigenvalues of the matrix that generates the data, we use
them as ground truth reference values. Such synthetic examples are useful for
the first phase of implementation and debugging of the code. For the data from
simulations of physical processes we study the numerical behaviour of the
computed eigenvalues along a sequence of sliding data windows.

##### Snapshot reconstruction.

The supplied data snapshots are represented using the Schur functions (see
Section 4.4) and the residuals of the representation are computed. At each
time step, the maximal reconstruction error over all data snapshots in the
active windows is computed and recorded. At the end of each test, the maximal
errors over all windows are displayed for all four methods.

##### Consistency.

We use relation (50) as an implicit test of the correctness of the
implementation by computing the residuals

$\left\|\left(\begin{smallmatrix}\zeta_{1}(\mathbf{x}_{j})\cr\zeta_{2}(\mathbf{x}_{j})\cr\vdots\cr\zeta_{r}(\mathbf{x}_{j})\end{smallmatrix}\right)-T^{T}\\!\left(\begin{smallmatrix}\zeta_{1}(\mathbf{x}_{j-1})\cr\zeta_{2}(\mathbf{x}_{j-1})\cr\vdots\cr\zeta_{r}(\mathbf{x}_{j-1})\end{smallmatrix}\right)\right\|_{2},\;\;j=2,3,\ldots,m,m+1.$

These residuals check whether, with given data, the data driven compression of
the Koopman operator, outlined in Section 3.1, was numerically computed to
sufficient accuracy. Further, if these residuals are large (meaning that the
known past and the present cannot be reproduced) then the forecasting is
doomed to fail. Maximal residuals are computed for each active window and
displayed for the methods in the second row of Table 1.

##### Forecasting skill.

Selected data sets are used to show that the new method has the prediction
skill as the standard DMD and EDMD implementations.

We should also emphasize what is not included in this presentation of the new
method. High performance software implementations, numerical robustness in
terms of the core numerical linear algebra, and adaptations to particular
applications (such as dynamic adaptation for forecasting dynamics with rapid
changes, or revealing coherent structures) are challenging themes that are the
subjects of separate work. A LAPACK style implementation that addresses more
numerical details, following our recent work [16], [17], is in progress and
will be soon available.

###  5.2 Examples

We have tested the new proposed scheme on many examples. Here we present the
result of two test cases and discuss the key elements of the experiment
outlined in Section 5.1.

###### Example 5.1

In the first example, we use a simulated data of a 2D laminar flow around a
cylinder. The $1201$ data snapshots of dimension $201\times 101$ are reshaped
into a $20301\times 1201$ matrix and used as test data, divided into
$\mathbf{X}$ and $\mathbf{Y}$ as explained in Section 1. The sliding data
window is set of fixed size $w=100$, and simulation is run for $100$ time
steps. At each step a new snapshot is received, the oldest snapshot is
discarded and the DMD, EDMD, Koopman-Schur-SSMD and Koopman-Schur-ESSMD
decompositions are computed. The Koopman-Schur-SSMD and the Koopman-Schur-
ESSMD are used for forecasting, to check the framework outlined in Section
4.5. The forecasting scheme is kept simple, following the formulas derived in
Section 4.5, and no attempt is made to dynamically resize the active window.
The forecasting horizon is set to $\tau=40$. Both kernel functions are tested,
and for the Gaussian kernel the parameter $\sigma$ is set to $300$ without any
fine tuning.

Figure 1 and Figure 2 show the reconstruction errors for all four methods and
the consistency errors for the Koopman-Schur based methods, respectively. The
eigenvalues computed by all four methods and the condition number of the key
matrices, using the kernel $k_{1}(\cdot,\cdot)$, are shown in Figure 3. The
prediction skills of the KS-SSMD and KS-EDMD (with the kernels
$k_{1}(\cdot,\cdot)$ and $k_{2}(\cdot,\cdot)$ ) are shown in Figure 4. The
forecast values are computed by direct implementation of the formulas from
Section 4.5, without any additional modifications (such as dynamic window
resizing, removing spurious modes, curbing ill-conditioning etc.), and the
results are encouraging and promising.

Figure 5 shows the eigenvalues of all four methods as computed from the last
active window (with the kernel $k_{2}(\cdot,\cdot)$) and the condition numbers
of the relevant matrices, As expected, there are two groups of the
eigenvalues: the one computed by DMD and KS-SSMD and the other by EDM and KS-
ESSMD. Since the two groups of methods operate using different function
dictionaries, they not necessarily return the same eigenvalues. But if the
values from those two groups get close together, it is an indication of
numerical convergence. We do not have analytical result but empirical evidence
is interesting, We have collected all computed Ritz values (from all $100$
windows) ans plotted them all as in the left panel in Figure 6. Then we
identified in the left panel in Figure 5 the eigenvalues that had been
computed nearly the same by all methods, and zoomed around them in the left
panel of Figure 6. Two particular selection are shown in the middle and the
right panel. Tightly clustered values computed by different methods are easily
identified by visual inspection. A more systematic empirical approach is to
use, say, machine learning techniques to cluster and follow the eigenvalues
computed from a sequence of sliding data windows,101010We tackle this in a
separate work.

Figure 1: (Example 5.1) The reconstruction error and the consistency of the
computed Schur functions. The kernel function is $k_{1}(\cdot,\cdot)$. The
yellow horizontal line indicates the level of machine precision eps times the
state space dimension $n$.

Figure 2: (Example 5.1) The reconstruction error and the consistency of the
computed Schur functions. The kernel function is $k_{2}(\cdot,\cdot)$. The
yellow horizontal line indicates the level of machine precision eps times the
state space dimension $n$.

Figure 3: (Example 5.1) _Left_ : The eigenvalues (Ritz values) computed in the
last active window. _Right_ : The relevant condition numbers. The kernel
function is $k_{1}(\cdot,\cdot)$.

Figure 4: (Example 5.1) The forecast error with the horizon $\tau=40$. _Left_
: Koopman-Schur-SSMD. _Middle_ : Koopman-Schur-ESSMD with the kernel function
$k_{1}(\cdot,\cdot)$. The KS-SSMD and the KS-ESSMD method are equivalent and
the difference is only due to the finite precision arithmetic. _Right_ :
Koopman-Schur-ESSMD with the kernel function $k_{2}(\cdot,\cdot)$.

Figure 5: (Example 5.1) _Left_ : The eigenvalues (Ritz values) computed in the
last active window. _Right_ : The relevant condition numbers. The kernel
function is $k_{2}(\cdot,\cdot)$.

Figure 6: (Example 5.1) _Left_ : All Ritz values computed throughout $100$
time steps. _Middle and right_ : Zoomed in around values identified in Figure
5 (all four methods computing nearly identical values) to show how the
sequence of Ritz values build tight clusters. This motivates using clustering
method throughout iterations and techniques based e.g. on diffusion maps for
tracking the Ritz values to reveal possible numerical convergence. (See Figure
8 below.)

$\boxtimes$

Example 5.1 is well-conditioned and it is a good test case to check the
concept and the software implementation. It shows good prediction skill of the
method based on the Schur decomposition.

###### Example 5.2

The next example uses vorticity data of a quasi-two-dimensional Kolmogorov
flow [55]. This is a difficult example, for discussions and details we refer
to [47], [55], [48]. (The methods tested here (Table 1) are, as in Example
5.1, entirely data driven and completely oblivious to the underlying nature of
the data.) The experiment is performed with a similar setting as in Example
5.1. Figure 7 shows the reconstruction errors and the condition numbers of two
runs, with the window sizes $w=70$ and $w=100$. The impact of the condition
number of the matrix of the modes to the reconstruction of the snapshots is
apparent.

Figure 7: (Example 5.2) _First row_ : Active window width $w=70$. _Second row_
: Active window width $w=100$. The kernel function is $k_{2}(\cdot,\cdot)$.

Figure 8: (Example 5.2) The eigenvalues from $100$ sliding windows with the
kernel $k_{2}$. _Left_ : $w=70$. _Middle_ : The eigenvalues with $w=100$ (red)
added to the plot with the eigenvalues with $w=100$ (blue). _Right_ : Zoomed
in the middle figure. The most outer blue and red eigenvalues agree to two or
three digits.

Figure 9: (Example 5.2) _Left_ : The relative errors (in norm) with the
prediction horizon $\tau=20$. _Middle_ : A detail of one forecast for $5$
steps ahead, zoomed to a part of the snapshot. _Right_ : A detail of one
forecast for $10$ steps ahead, zoomed to a part of the snapshot.

##  6 Conclusions and ongoing and future work

We have presented a new theoretical and computational framework for Koopman
operator based data driven analysis of nonlinear dynamics. Preliminary testing
has shown good performance and, most importantly, the KS-SSMD is based on the
unitary Schur decomposition and it is thus independent of the
diagonalizability of the EDMD matrix, or the (potentially high) condition
number of the eigenvectors.

The ongoing and planned work, based on this report, includes

  * •

Further operator theoretic analysis of the approximation error.

  * •

A LAPACK style implementation of the presented algorithms, following our
recent implementations of the DMD and the Hermitian DMD [16], [19]. (The
Hermitian DMD [17], [18] is a special case of the KS-SSMD, where the Schur
form is actually diagonal.)

  * •

A QR compressed implementation, including a streaming KS-SSMD based on fast
updating/down-dating of the QR compressed representation of the snapshots.

  * •

Algorithmic details and software implementation for the applications:

    * –

Identification of the coherent structures and their representations using
flags of invariant subspaces.

    * –

Forecasting with adaptation to rapid changes in the dynamics.

## References

  * [1] S. Ahuja and C. W. Rowley.  Feedback control of unstable steady states of flow past a flat plate using reduced-order estimators.  Journal of Fluid Mechanics, 645:447–478, 2010.
  * [2] E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. J. Dongarra, J. Du Croz, S. Hammarling, A. Greenbaum, A. McKenney, and D. Sorensen.  LAPACK Users’ Guide (Third Ed.).  Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1999.
  * [3] Hassan Arbabi and Igor Mezic.  Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the koopman operator.  SIAM Journal on Applied Dynamical Systems, 16(4):2096–2126, 2017.
  * [4] T. Askham and J. Kutz.  Variable projection methods for an optimized dynamic mode decomposition.  SIAM Journal on Applied Dynamical Systems, 17(1):380–416, 2018.
  * [5] Zhaojun Bai and James W. Demmel.  On swapping diagonal blocks in real Schur form.  Linear Algebra and its Applications, 186:75–95, 1993.
  * [6] E. Berger, M.Sastuba, D. Vogt, B. Jung, and H. B Amor.  Estimation of perturbations in robotic behavior using dynamic mode decomposition.  Advanced Robotics, 29(5):331–343, 2015.
  * [7] Steven L. Brunton, Marko Budišić, Eurika Kaiser, and J. Nathan Kutz.  Modern koopman theory for dynamical systems.  SIAM Review, 64(2):229–340, 2022.
  * [8] 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(6):887–915, April 2012.
  * [9] Vaughn Climenhaga, Stefano Luzzatto, and Yakov Pesin.  The geometric approach for constructing sinai–ruelle–bowen measures.  Journal of Statistical Physics, 166(3-4):467–493, 2017.
  * [10] Matthew J. Colbrook, Lorna J. Ayton, and Máté Szőke.  Residual dynamic mode decomposition: robust and verified koopmanism.  Journal of Fluid Mechanics, 955:A21, 2023.
  * [11] Matthew J. Colbrook and Alex Townsend.  Rigorous data-driven computation of spectral properties of Koopman operators for dynamical systems.  arXiv e-prints, page arXiv:2111.14889, November 2021.
  * [12] S. T. M. Dawson, M. S. Hemati, M. O. Williams, and C. W. Rowley.  Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition.  Experiments in Fluids, 57(3):42, 2016.
  * [13] J. W. Demmel.  A Numerical Analyst’s Jordan Canonical Form.  PhD thesis, Center for Pure and Applied Mathematics, University of California, Berkeley, 1983.
  * [14] Z. Drmač.  Dynamic Mode Decomposition-A Numerical Linear Algebra Perspective, pages 161–194.  Springer International Publishing, Cham, 2020.
  * [15] Z. Drmač, I. Mezić, and R. Mohr.  Data driven modal decompositions: analysis and enhancements.  SIAM Journal on Scientific Computing, 40(4):A2253–A2285, 2018.
  * [16] Z. Drmač.  A LAPACK implementation of the Dynamic Mode Decomposition I.  Technical report, Department of Mathematics, University of Zagreb, Croatia, and AIMdyn Inc. Santa Barbara, CA, October 2022.  LAPACK Working Note 298.
  * [17] Z. Drmač.  A LAPACK implementation of the Dynamic Mode Decomposition II.  Technical report, Department of Mathematics, University of Zagreb, Croatia, and AIMdyn Inc. Santa Barbara, CA, November 2022.  LAPACK Working Note 300.
  * [18] Zlatko Drmač.  Hermitian dynamic mode decomposition – numerical analysis and software solution.  ACM Trans. Math. Soft. (in revision), 00:72–83, 2023.
  * [19] Zlatko Drmač.  A LAPACK implementation of the Dynamic Mode Decomposition.  ACM Trans. Math. Soft. (in revision), 00:??–??, 2023.
  * [20] Zlatko Drmač, Igor Mezić, and Ryan Mohr.  On least squares problems with certain Vandermonde–Khatri–Rao structure with applications to DMD.  SIAM Journal on Scientific Computing, 42(5):A3250–A3284, 2020.
  * [21] Zlatko Drmač, Igor Mezić, and Ryan Mohr.  Identification of nonlinear systems using the infinitesimal generator of the koopman semigroup – a numerical implementation of the mauroy-goncalves method.  Mathematics, 9(17), 2021.
  * [22] Nelson Dunford et al.  Spectral operators.  Pacific Journal of Mathematics, 4(3):321–354, 1954.
  * [23] Izrail Moiseevich Gelfand and Georgiĭ Evgenevich Shilov.  Generalized Functions, Volume 2: Spaces of Fundamental and Generalized Functions, volume 261.  American Mathematical Soc., 2016.
  * [24] S. Ghosal, V. Ramanan, S. Sarkar, S. Chakravarthy, and S. Sarkar.  Detection and analysis of combustion instability from hi-speed flame images using dynamic mode decomposition.  In ASME. Dynamic Systems and Control Conference, Volume 1, 2016.
  * [25] G. H. Golub, V. C. Klema, and G. W. Stewart.  Rank degeneracy and least squares problems.  Technical Report CS-TR-76-559, Stanford, CA, USA, 1976.
  * [26] M. S. Hemati, C. W. Rowley, E. A. Deem, and L. N. Cattafesta.  De-biasing the dynamic mode decomposition for applied Koopman spectral analysis.  ArXiv e-prints, February 2015.
  * [27] M. S. Hemati, M. O. Williams, and C. W. Rowley.  Dynamic mode decomposition for large and streaming datasets.  Physics of Fluids, 26(11):111701, November 2014.
  * [28] R. A. Horn and C. R. Johnson.  Matrix Analysis.  Cambridge University Press, 1990.
  * [29] R. A. Horn and C. R. Johnson.  Topics in Matrix Analysis.  Cambridge University Press, 1991.
  * [30] M. R. Jovanović, P. J. Schmid, and J. W. Nichols.  Low–rank and sparse dynamic mode decomposition.  In Annual Research Briefs 2012, pages 139–152. Center for Turbulence Research, Stanford University, 2012.
  * [31] M. R. Jovanović, P. J. Schmid, and J. W. Nichols.  Sparsity-promoting dynamic mode decomposition.  Phys. Fluids, 26(2):024103 (22 pages), 2014.
  * [32] Milan Korda and Igor Mezić.  On convergence of Extended Dynamic Mode Decomposition to the Koopman operator.  Journal of Nonlinear Science, 28:687–710, April 2018.
  * [33] Daniel Kressner.  Block algorithms for reordering standard and generalized schur forms.  ACM Trans. Math. Softw., 32(4):521–532, dec 2006.
  * [34] Fumi-Yuki Maeda.  Function of generalized scalar operators.  Journal of Science of the Hiroshima University, Series AI (Mathematics), 26(2):71–76, 1962.
  * [35] Alexandre Mauroy, Igor Mezić, and Yoshihiko Susuki.  The Koopman Operator in Systems and Control: Concepts, Methodologies, and Applications.  01 2020.
  * [36] Alexandre Mauroy, Y Susuki, and I Mezić.  Koopman operator in systems and control.  Springer, 2020.
  * [37] Igor Mezić.  Spectral properties of dynamical systems, model reduction and decompositions.  Nonlinear Dynamics, 41:309–325, 2005.
  * [38] Igor Mezić.  Spectrum of the koopman operator, spectral expansions in functional spaces, and state-space geometry.  Journal of Nonlinear Science, 30(5):2091–2145, 2020.
  * [39] Igor Mezić.  Spectrum of the koopman operator, spectral expansions in functional spaces, and state-space geometry.  Journal of Nonlinear Science, 30(5):2091–2145, 2020.
  * [40] Igor Mezić.  On numerical approximations of the koopman operator.  Mathematics, 10(7):1180, 2022.
  * [41] Igor Mezic, Zlatko Drmac, Nelida Crnjaric-Zic, Senka Macesic, Maria Fonoberova, Ryan Mohr, Allan Avila, Iva Manojlovic, and Aleksandr Andrejcuk.  A Koopman Operator-Based Prediction Algorithm and its Application to COVID-19 Pandemic.  arXiv e-prints, page arXiv:2304.13601, April 2023.
  * [42] J. Proctor, S. Brunton, and J. Kutz.  Dynamic mode decomposition with control.  SIAM Journal on Applied Dynamical Systems, 15(1):142–161, 2016.
  * [43] Robert R Rogers.  Triangular form for bounded linear operators.  Journal of functional analysis, 88(1):135–152, 1990.
  * [44] P. J. Schmid, L. Li, M. P. Juniper, and O. Pust.  Applications of the dynamic mode decomposition.  Theoretical and Computational Fluid Dynamics, 25(1):249–259, 2011.
  * [45] Peter Schmid and Jörn Sesterhenn.  Dynamic Mode Decomposition of numerical and experimental data.  Journal of Fluid Mechanics, 656, 11 2008.
  * [46] Peter J Schmid.  Nonmodal stability theory.  Annual review of fluid mechanics, 39(1):129–162, 2007.
  * [47] Balachandra Suri.  Quasi-two-dimensional Kolmogorov flow: Bifurcations and exact coherent structures.  PhD thesis, Georgia Institute of Technology, 2017.
  * [48] Balachandra Suri, Jeffrey Tithof, Roman O. Grigoriev, and Michael F. Schatz.  Forecasting fluid flows using the geometry of turbulence.  Phys. Rev. Lett., 118:114501, Mar 2017.
  * [49] Yoshihiko Susuki, Igor Mezic, Fredrik Raak, and Takashi Hikihara.  Applied koopman operator theory for power systems technology.  Nonlinear Theory and Its Applications, IEICE, 7(4):430–459, 2016.
  * [50] N. Takeishi, Y. Kawahara, Y. Tabei, and T. Yairi.  Bayesian dynamic mode decomposition.  In Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17, pages 2814–2821, 2017.
  * [51] N. Takeishi, Y. Kawahara, and T. Yairi.  Learning Koopman invariant subspaces for dynamic mode decomposition.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 1130–1140. Curran Associates, Inc., 2017.
  * [52] N. Takeishi, Y. Kawahara, and T. Yairi.  Sparse nonnegative dynamic mode decomposition.  In 2017 IEEE International Conference on Image Processing (ICIP), pages 2682–2686, Sep. 2017.
  * [53] N. Takeishi, Y. Kawahara, and T. Yairi.  Subspace dynamic mode decomposition for stochastic Koopman analysis.  Phys. Rev. E, 96:033310, Sep 2017.
  * [54] Makhin Thitsa, Maison Clouatre, Erik Verriest, Samuel Coogan, and Clyde Martin.  A numerically stable dynamic mode decomposition algorithm for nearly defective systems.  IEEE Control Systems Letters, 5(1):67–72, 2021.
  * [55] Jeffrey Tithof, Balachandra Suri, Ravi Kumar Pallantla, Roman Grigoriev, and Michael Schatz.  Bifurcations in a quasi-two-dimensional kolmogorov-like flow.  Journal of Fluid Mechanics, 828:837–866, 10 2017.
  * [56] A. van der Sluis.  Condition numbers and equilibration of matrices.  Numerische Mathematik, 14:14–23, 1969.
  * [57] 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(6):1307–1346, 2015.
  * [58] A. Wynn, D. S. Pearson, B. Ganapathisubramani, and P. J. Goulart.  Optimal mode decomposition for unsteady flows.  Journal of Fluid Mechanics, 733:473–503, 2013.