# Krylov Subspace Method for Nonlinear Dynamical Systems with Random Noise

Yuka Hashimoto<EMAIL_ADDRESS>
NTT Network Technology Laboratories, NTT Corporation
3-9-11, Midori-cho, Musashinoshi, Tokyo, 180-8585, Japan /
Graduate School of Science and Technology, Keio University
3-14-1 Hiyoshi, Kohoku, Yokohama, Kanagawa, 223-8522, Japan Isao Ishikawa
<EMAIL_ADDRESS>
Center for Advanced Intelligence Project, RIKEN
1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan /
Faculty of Science and Technology, Keio University
3-14-1 Hiyoshi, Kohoku, Yokohama, Kanagawa, 223-8522, Japan Masahiro Ikeda
<EMAIL_ADDRESS>
Center for Advanced Intelligence Project, RIKEN
1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan /
Faculty of Science and Technology, Keio University
3-14-1 Hiyoshi, Kohoku, Yokohama, Kanagawa, 223-8522, Japan Yoichi Matsuo
<EMAIL_ADDRESS>
NTT Network Technology Laboratories, NTT Corporation
3-9-11, Midori-cho, Musashinoshi, Tokyo, 180-8585, Japan Yoshinobu Kawahara
<EMAIL_ADDRESS>
Center for Advanced Intelligence Project, RIKEN
1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan /
Institute of Mathematics for Industry, Kyushu University
744, Motooka, Nishi-ku, Fukuoka, 819-0395, Japan

###### Abstract

Operator-theoretic analysis of nonlinear dynamical systems has attracted much
attention in a variety of engineering and scientific fields, endowed with
practical estimation methods using data such as dynamic mode decomposition. In
this paper, we address a lifted representation of nonlinear dynamical systems
with random noise based on transfer operators, and develop a novel Krylov
subspace method for estimating the operators using finite data, with
consideration of the unboundedness of operators. For this purpose, we first
consider Perron-Frobenius operators with kernel-mean embeddings for such
systems. We then extend the Arnoldi method, which is the most classical type
of Kryov subspace method, so that it can be applied to the current case.
Meanwhile, the Arnoldi method requires the assumption that the operator is
bounded, which is not necessarily satisfied for transfer operators on
nonlinear systems. We accordingly develop the shift-invert Arnoldi method for
Perron-Frobenius operators to avoid this problem. Also, we describe an
approach of evaluating predictive accuracy by estimated operators on the basis
of the maximum mean discrepancy, which is applicable, for example, to anomaly
detection in complex systems. The empirical performance of our methods is
investigated using synthetic and real-world healthcare data.

Keywords: Nonlinear dynamical system, Transfer operator, Krylov subspace
methods, Operator theory, Time-series data

##  1 Introduction

Analyzing nonlinear dynamical systems using data is one of the fundamental but
still challenging problems in various engineering and scientific fields.
Recently, operator-theoretic analysis has attracted much attention for this
purpose, with which the behavior of a nonlinear dynamical system is analyzed
through representations with transfer operators such as Koopman operators and
their adjoint ones, Perron-Frobenius operators (Budišić et al., 2012;
Kawahara, 2016). Since transfer operators are linear even if the corresponding
dynamical systems are nonlinear, we can apply sophisticated theoretical
results and useful tools of the operator theory, and access the properties of
dynamics more easily from both theoretical and practical viewpoints. This is
one of the main advantages of using transfer operators compared with other
methods for learning dynamical systems such as using recurrent neural networks
(RNNs) and hidden Markov models. For example, one could consider modal
decomposition of nonlinear dynamics by using the spectral analysis in operator
theory, which provides the global characteristics of the dynamics and is
useful in understanding complex phenomena (Kutz, 2013). This topic has also
been recently discussed in machine learning (Kawahara, 2016; Lusch et al.,
2017; Takeishi et al., 2017b).

However, many of the existing works mentioned above are on deterministic
dynamical systems. Quite recently, the extension of these works to random
systems has been addressed in a few works. The methods for analyzing
deterministic systems with transfer operators are extended to cases in which
dynamical systems are random (Črnjarić-Žic et al., 2017; Takeishi et al.,
2017a). Also, the transfer operator for a stochastic process in reproducing
kernel Hilbert spaces (RKHSs) is defined (Klus et al., 2017), which provides
an approach of analyzing dynamics of random variables in RKHSs.

In this paper, we address a lifted representation of nonlinear dynamical
systems with random noise based on transfer operators, and develop a novel
Krylov subspace method for estimating the operator using finite data, with
consideration of the unboundedness of operators. To this end, we first
consider Perron-Frobenius operators with kernel-mean embeddings for such
systems. We then extend the Arnoldi method, which is the most classical type
of Krylov subspace methods, so that it can be applied to the current case.
However, although transfer operators on nonlinear systems are not necessarily
bounded, the Arnoldi method requires the assumption on the boundedness of an
operator. We accordingly develop the shift-invert Arnoldi method for the
Perron-Frobenius operators to avoid this problem. Moreover, we consider an
approach of evaluating the predictive accuracy with estimated operators on the
basis of the maximum mean discrepancy (MMD), which is applicable, for example,
to anomaly detection in complex systems. Finally, we investigate the empirical
performance of our methods using synthetic data and also apply those to
anomaly detection with real-world healthcare data.

The remainder of this paper is organized as follows. First, in Section 2, we
review transfer operators and Krylov subspace methods. In Section 3, we
consider Perron-Frobenius operators with kernel-mean embeddings for nonlinear
dynamical systems with random noises. In Section 4, we develop Krylov subspace
methods for estimating these operators using data, and in Section 5, we
discuss the connection of our methods to existing methods. In Section 6, we
consider an approach of evaluating the prediction accuracy with estimated
operators. Finally, we empirically investigate the performance of our methods
in Section 7 and conclude the paper in Section 8. Proofs which are not given
after their statements are given in Appendix A.

#### Notations

Standard capital letters and ornamental capital letters denote the infinite
dimensional linear operators. Bold letters denote the matrices (finite
dimensional linear operators) or finite dimensional vectors. Calligraphic
capital letters and itzlicized Greek capital letters denote sets. The inner
product and norm in $\mathcal{H}_{k}$ are denoted as
$\left\langle\cdot,\cdot\right\rangle_{k}$ and $\|\cdot\|_{k}$, respectively.
The operator norm of a bounded linear operator $A$ in $\mathcal{H}_{k}$, which
is defined as $\sup_{v\in\mathcal{H}_{k},\|v\|_{k}=1}\|Av\|_{k}$ is denoted as
$|||A|||_{k}$. While, the Euclid norm in $\mathbb{C}^{S}$ for $S\in\mathbb{N}$
is denoted as $\|\cdot\|$, and the operator norm of a matrix $\mathbf{A}$ is
denoted as $|||\mathbf{A}|||$.

The typical notations in this paper are listed in Table 1.

Table 1: Notation table $(\varOmega,\mathcal{F},P)$ |  A measurable space (sample space) with a probability measure $P$
---|---
$(\mathcal{X},\mathcal{B})$ |  A Borel measurable and locally compact Hausdorff vector space (state space)
$x_{t}$ |  A random variable from $\varOmega$ to $\mathcal{X}$ represents the observation at $t$
$\\{\xi_{t}\\}$ |  An i.i.d. stochastic process corresponds to the random noise, where $\xi_{t}:\varOmega\to\mathcal{X}$
$k$ |  A positive-definite continuous, bounded and $c_{0}$-universal kernel on $\mathcal{X}$
$\phi$ |  The feature map endowed with $k$
$\mathcal{H}_{k}$ |  The RKHS endowed with $k$
$\mathcal{M}(\mathcal{X})$ |  The set of all finite signed measures on $\mathcal{X}$
$\Phi$ |  The kernel mean embedding $\mathcal{M}(\mathcal{X})\to\mathcal{H}_{k}$ defined by $\mu\mapsto\int_{x\in\mathcal{X}}\phi(x)\;d\mu(x)$
$K$ |  A Perron-Frobenius operator
$\operatorname*{\mathcal{D}}(A)$ |  The domain of a linear operator $A$
$\varLambda(A)$ |  The spectrum of an $A$
$\mathcal{W}(A)$ |  The numerical range of an $A$ on $\mathcal{H}$ defined by $\\{\left\langle Av,v\right\rangle\mid\ v\in\mathcal{H},\ \|v\|=1\\}$
$s_{\operatorname{min}}(\mathbf{A})$ |  The minimal singular value of a matrix $\mathbf{A}$ defined by $\min_{\|\mathbf{w}\|_{=}1}\|\mathbf{Aw}\|$
$\gamma$ |  A parameter to transform $K$ to a bounded bijective operator $(\gamma I-K)^{-1}$ which is not in $\varLambda(K)$
$\\{\tilde{x}_{0},\tilde{x}_{1},\ldots\\}$ |  Observed time-series data
$S$ |  A natural number that represents the dimension of the Krylov subspace
$N$ |  A natural number that represents the amount of observed data used for the estimation
$\mu_{t,N}$ |  The empirical measure generated by finite observed data $\\{\tilde{x}_{t},\ldots,\tilde{x}_{t+(S+1)(N-1)}\\}$
$\mu_{t}$ |  The weak limit of $\mu_{t,S}$ in $\mathcal{M}(\mathcal{X})$
$\mathcal{V}(A,v)$ |  The Krylov subspace of a linear operator $A$ and a vector $v$
$Q_{S}$ |  The linear operator from $\mathbb{C}^{S}$ to $\mathcal{H}_{k}$ composed of the orthonormal basis of the Krylov subspace
$\mathbf{R}_{S}$ |  The $S$ times $S$ matrix which transforms the coordinate into the one with the orthonormal basis
$\tilde{\mathbf{K}}_{S}$ |  The estimation of $K$ in an $S$-dimensional Krylov subspace
$\tilde{\mathbf{L}}_{S}$ |  The estimation of $(\gamma I-K)^{-1}$ in a $S$-dimensional Krylov subspace
$a_{t,S}$ |  The abnormality at $t$ computed with $\tilde{\mathbf{K}}_{S}$

##  2 Background

###  2.1 Transfer operators

Consider a deterministic dynamical system

$x_{t+1}=h(x_{t}),$

where $h\colon\mathcal{X}\to\mathcal{X}$ is a map, $\mathcal{X}$ is a state
space and $x_{t}\in\mathcal{X}$. Then, the corresponding Koopman operator
(Koopman, 1931), which is denoted as $\mathscr{K}$, is a linear operator in
some subspace $\mathcal{M}\subseteq\\{g\colon\mathcal{X}\to\mathcal{X}\\}$,
defined by

$\mathscr{K}g:=g\circ h$

for $g\in\mathcal{M}$. From the definition, $\mathscr{K}$ represents the time
evolution of the system as $(\mathscr{K}^{n}g)(x_{0})=g(h(\ldots
h(x_{0})))=g(x_{n})$. Since the Koopman operator is linear even when the
dynamical system $h$ is nonlinear, the operator theory is valid for analyzing
it. And, the adjoint of Koopman operator is called Perron-Frobenius operator.
The concept of the RKHS is combined with transfer operators, and Perron-
Frobenius operators in an RKHS are addressed Kawahara (2016); Ishikawa et al.
(2018). One of the advantages of using transfer operators in RKHSs is that
they can describe dynamical systems defined in non-Euclidean spaces. Let
$\mathcal{H}_{k}$ be the RKHS endowed with a positive definite kernel $k$, and
let $\phi:\mathcal{X}\rightarrow\mathcal{H}_{k}$ be the feature map. Then, the
Perron-Frobenius operator in the RKHS for $h\colon\mathcal{X}\to\mathcal{X}$,
which is denoted by $\mathscr{K}_{\operatorname{RKHS}}$, is a linear operator
in $\mathcal{H}_{k}$ defined by

$\mathscr{K}_{\operatorname{RKHS}}\phi(x):=\phi(h(x))$

for $\phi(x)\in\operatorname{Span}\\{\phi(x)\mid x\in\mathcal{X}\\}$.

Transfer operator has also been discussed for cases in which a dynamical
system is random. Let $(\mathcal{X},\mathcal{B},\mu)$ and
$(\varOmega,\mathcal{F},P)$ be probability spaces. The following random system
is considered (Črnjarić-Žic et al., 2017; Takeishi et al., 2017a):

$x_{t+1}=\pi(t,\omega,x_{t}),$

where $\pi\colon\mathbb{Z}_{\geq
0}\times\varOmega\times\mathcal{X}\to\mathcal{X}$ is a map and
$x_{t}\in\mathcal{X}$. Then, the Koopman operator, which is denoted as
$\bar{\mathscr{K}}_{t}$, is a linear operator in
$\mathcal{L}^{2}(\mathcal{X})$ and defined as

$\bar{\mathscr{K}}_{t}g:=\int_{\omega\in\varOmega}g(\pi(t,\omega,\cdot))\;dP(\omega)$

for $g\in\mathcal{L}^{2}(\mathcal{X})$. Also, Perron-Frobenius operators in
RKHSs for a stochastic process $\\{x_{t}\\}$ on $(\mathcal{X},\mathcal{B})$
whose probability density functions are $\\{p_{t}\\}$ are considered (Klus et
al., 2017). The Perron-Frobenius operator in the RKHS, which is denoted as
$\bar{\mathscr{K}}_{\operatorname{RKHS},t}$, is a linear operator in
$\mathcal{H}_{k}$ and defined by

$\bar{\mathscr{K}}_{\operatorname{RKHS},t}\mathscr{E}p_{t}:=\mathscr{U}p_{t},$

where $\mathscr{E}$ and $\mathscr{U}$ are respectively the embeddings of
probability densities to $\mathcal{H}_{k}$ defined as
$q\mapsto\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}\phi(y)p(y\mid
x)q(x)\;d\mu(y)\;d\mu(x)$ and
$q\mapsto\int_{x\in\mathcal{X}}\phi(x)q(x)\;d\mu(x)$, and $p$ is a function
satisfying $P(x_{t+1}\in A\mid\\{x_{t}=x\\})=\int_{y\in A}p(y\mid
x)\;d\mu(y)$.

###  2.2 Unbounded linear operators

First, we review the definition of a linear operator in a Hilbert space
$\mathcal{H}$.

###### Definition 1

Let $\mathcal{S}$ be a dense subset of $\mathcal{H}$. A linear operator $A$ in
$\mathcal{H}$ is a linear map $A:\mathcal{S}\to\mathcal{H}$. The set
$\mathcal{S}$, which is denoted as $\operatorname*{\mathcal{D}}(A)$, is called
the domain of $A$. If there exists $C>0$ such that the operator norm of $A$,
which is defined as $\|A\|:=\sup_{v\in\mathcal{H},\|v\|=1}\|Av\|$ is bounded
by $C$, then $A$ is called bounded.

For a linear operator $A$, the spectrum and numerical range are defined as
follows:

###### Definition 2

Let $\Gamma(A)$ be the set of $\lambda\in\mathbb{C}$ such that $(\lambda
I-A):\operatorname*{\mathcal{D}}(A)\to\mathcal{H}$ is bijective and $(\gamma
I-A)^{-1}$ is bounded. The spectrum of $A$ is the set
$\mathbb{C}\setminus\Gamma(A)$, which is denoted as $\varLambda(A)$. Moreover,
the numerical range of $A$ is the set $\\{\left\langle
Av,v\right\rangle\in\mathbb{C}\mid v\in\operatorname*{\mathcal{D}}(A),\
\|v\|=1\\}$, which is denoted as $\mathcal{W}(A)$.

If $A$ is bounded, it can be shown that $\varLambda(A)$ is nonempty and
compact (Kubrusly, 2012, Theorem 2.1, Theorem 2.2). Also, by Toeplitz-
Hausdorff theorem, it can be shown that $\mathcal{W}(A)$ is bounded and convex
(McIntosh, 1978). The relation between $\varLambda(A)$ and $\mathcal{W}(A)$ is
characterized by the inclusion
$\varLambda(A)\subseteq\overline{\mathcal{W}(A)}$. However, if $A$ is
unbounded, neither $\varLambda(A)$ nor $\mathcal{W}(A)$ is always bounded.

###  2.3 Krylov subspace methods

Krylov subspace methods are numerical methods for estimating the behavior of a
linear operator by projecting it to a finite dimensional subspace, called
Krylov subspace. Let $A$ be a linear operator in Hilbert space $\mathcal{H}$
and $v\in\mathcal{H}$. Then, the Krylov subspace of $A$ and $v$, which is
denoted by $\mathcal{V}_{S}(A,v)$, is an $S$-dimensional subspace:

$\operatorname{Span}\\{v,Av,\ldots,A^{S-1}v\\}.$

Krylov subspace methods are often applied to compute the spectrum of
$\mathbf{A}$, $\mathbf{A}^{-1}\mathbf{v}$, or $f(\mathbf{A})\mathbf{v}$ for a
given large and sparse $N\times N$ matrix $\mathbf{A}$, vector
$\mathbf{v}\in\mathbb{C}^{N}$ and function $f$ (Krylov, 1931; R Hestenes and
Stiefel, 1952; Saad and Schultz, 1983; Gallopoulos and Saad, 1992; Moret and
Novati, 2004). The theoretical extensions of Krylov subspace methods for
linear operators in infinite dimensional Hilbert spaces are explored in Güttel
(2010); Grimm (2012); Göckler (2014); Hashimoto and Nodera (2018) to deal with
matrices that are finite dimensional approximations of infinite dimensional
linear operators.

The Arnoldi method is a classical and most commonly-used Kryov subspace
method. With the Arnoldi method, the Krylov subspace $\mathcal{V}_{S}(A,v)$ is
first constructed, and $A$ is projected onto it. For a matrix $\mathbf{A}$ and
vector $\mathbf{v}$, since the basis of
$\mathcal{V}_{S}(\mathbf{A},\mathbf{v})$ can be computed only by matrix-vector
products, the projection of $\mathbf{A}$ is also obtained only with matrix-
vector products. Note that the computational cost of the matrix-vector product
is less than or equal to $O(N^{2})$, which is less computationally expensive
than computing the spectrum of $\mathbf{A}$, $\mathbf{A}^{-1}$ or
$f(\mathbf{A})$ directly.

On the other hand, $\mathbf{A}$ is often the matrix approximation of an
unbounded $A$, that is, the spatial discretization of $A$. Theoretically, if
$\mathbf{A}$ is an unbounded operator, $\mathbf{A}^{i}v$ for $i=1,\ldots,S-1$
cannot always be defined, and practically, although $\mathbf{A}$ is a matrix
(bounded), the performance of the Arnoldi method for $\mathbf{A}$ degrades due
to the unboundedness of the original $A$. To overcome this issue, the shift-
invert Arnoldi method, that constructs the Krylov subspace
$\mathcal{V}_{S}((\gamma I-A)^{-1},v)$, where $\gamma$ is not in the spectrum
of $A$, has been investigated. Since $(\gamma I-A)^{-1}$ is bounded, $(\gamma
I-A)^{-i}v$ for $i=1,\ldots,S-1$ is always defined. Thus, the Krylov subspace
$\mathcal{V}_{S}((\gamma I-A)^{-1},v)$ can be constructed. This improve the
performance for matrix $\mathbf{A}$, which is the matrix approximation of
unbounded $A$.

Moreover, the application of the Arnoldi method to estimating transfer
operators has been discussed for the deterministic case
$\mathscr{K}_{\operatorname{RKHS}}$ (Kawahara, 2016) and for the random case
$\bar{\mathscr{K}}_{t}$ (Črnjarić-Žic et al., 2017). An advantage of the
Krylov subspace methods for estimating transfer operators is that they require
one time-series dataset embedded by one observable function or one feature
map, which matches the case of using an RKHS. Meanwhile, the largest
difference between the Krylov subspace methods mentioned in the preceding
paragraphs and those for transfer operators is that the operator to be
estimated is given beforehand or not. That is, calculations appear in Krylov
subspace methods for transfer operators need to be carried out without knowing
the operators.

##  3 Perron-Frobenius Operators with Kernel-Mean Embeddings

Consider the following discrete-time nonlinear dynamical systems with random
noise in $\mathcal{X}$:

$x_{t+1}=h(x_{t})+\xi_{t},$ (1)

where $t\in\mathbb{Z}_{\geq 0}$, $(\varOmega,\mathcal{F})$ is a measurable
space (corresponding to a sample space), $(\mathcal{X},\mathcal{B})$ is a
Borel measurable and locally compact Hausdorff vector space (corresponding to
a state space), $x_{t}$ and $\xi_{t}$ are random variables from sample space
$\varOmega$ to state space $\mathcal{X}$, and
$h\colon\mathcal{X}\to\mathcal{X}$ is a map. Let $P$ be a probability measure
on $\varOmega$. Examples of locally compact Hausdorff space are
$\mathbb{R}^{d}$ and Riemannian manifolds. Assume that $\xi_{t}$ with
$t\in\mathbb{Z}_{\geq 0}$ is an i.i.d. stochastic process and is independent
of $x_{t}$. The $\xi_{t}(\omega)$ corresponds to the random noise in
$\mathcal{X}$. Let $k\colon\mathcal{X}\times\mathcal{X}\to\mathbb{C}$ be a
positive-definite kernel on $\mathcal{X}$. According to the Moore-Aronszajn
theorem, there exists a unique RKHS endowed with a dot product
$\left\langle\cdot,\cdot\right\rangle_{k}$, which is denoted as
$\mathcal{H}_{k}$. The corresponding feature map is denoted by
$\phi:\mathcal{X}\to\mathcal{H}_{k}$. In this paper, we assume that $k$ is
continuous, bounded and $c_{0}$-universal, i.e.,
$\phi(x)\in\mathcal{C}_{0}(\mathcal{X})$ for all $x\in\mathcal{X}$ and
$\mathcal{H}_{k}$ is dense in $\mathcal{C}_{0}(\mathcal{X})$, where
$\mathcal{C}_{0}(\mathcal{X})$ is the space of all continuous functions vanish
at infinity (Sriperumbudur et al., 2011). For example, the Gaussian kernel
$e^{-c\|x-y\|^{2}}$ and the Laplacian kernel $e^{-c\|x-y\|_{1}}$ with $c>0$
for $x,y\in\mathcal{X}$ with $\mathcal{X}=\mathbb{R}^{d}$ are continuous and
bounded $c_{0}$-universal kernels.

Now, we consider the transformation of the random variables in dynamical
system (1) into probability measures to capture the time evolution of the
system starting from several initial states. That is, random variable $x$ is
transformed into probability measure $x_{*}P$, where $x_{*}P$ denotes the push
forward measure of $P$ with respect to $x$, defined by
$x_{*}P(B)=P(x^{-1}(B))$ for $B\in\mathcal{B}$. This transformation replaces
the nonlinear relation $h$ between $x_{t}$ and $x_{t+1}$ with a linear one
between probability measures. Also, let
$\beta_{t}\colon\mathcal{X}\times\varOmega\to\mathcal{X}$ be a map defined by
$(x,\omega)\mapsto h(x)+\xi_{t}(\omega)$. Then, a linear map
$\mu\mapsto{\beta_{t}}_{*}(\mu\otimes P)$ is considered for a probability
measure $\mu$, instead of $h$. Also, we embed the probability measures into
Hilbert space $\mathcal{H}_{k}$, which defines an inner product between
probability measures, to apply the operator theory. Referring to Klus et al.
(2017), this embedding is possible by the kernel mean embedding (Muandet et
al., 2017) as follows. Let $\mathcal{M}(\mathcal{X})$ be the set of all finite
signed measures on $\mathcal{X}$. Then, the kernel mean embedding
$\Phi\colon\mathcal{M}(\mathcal{X})\to\mathcal{H}_{k}$ is defined by
$\mu\mapsto\int_{x\in\mathcal{X}}\phi(x)\;d\mu(x)$.

As a result, the Perron-Frobenius operator for dynamical system (1) is defined
with $\beta_{t}$ and the kernel mean embedding $\Phi$ as follows.

###### Definition 3

The Perron-Frobenius operator for the system (1),
$K:\Phi(\mathcal{M}(\mathcal{X}))\to\mathcal{H}_{k}$, is defined as

$K\Phi(\mu):=\Phi({\beta_{t}}_{*}(\mu\otimes P)).$ (2)

That is, $K$ transfers the measure generated by $x_{t}$ to that by $x_{t+1}$.
In fact, the following lemma holds.

###### Lemma 4

The relation $K\Phi({x_{t}}_{*}P)=\Phi({x_{t+1}}_{*}P)$ holds.

Before discussing the estimation of $K$, we here describe some basic
properties of the kernel mean embedding $\Phi$ and $K$, which are summarized
as follows:

###### Lemma 5

The kernel mean embedding
$\Phi\colon\mathcal{M}(\mathcal{X})\to\mathcal{H}_{k}$ is a linear and
continuous operator.

###### Lemma 6

The Perron-Frobenius operator
$K:\Phi(\mathcal{M}(\mathcal{X}))\to\mathcal{H}_{k}$ does not depend on $t$,
is well-defined and is a linear operator.

Also, the following two propositions show the connections of $K$ to the
existing operators (stated in Section 2.1). We have the following relations of
$K$ with $\bar{\mathscr{K}}_{\operatorname{RKHS},t}$ and with
$\bar{\mathscr{K}}_{t}$:

###### Proposition 7

If the stochastic process $\\{x_{t}\\}$ considered in Klus et al. (2017)
satisfies $x_{t+1}=h(x_{t})+\xi_{t}$, then
$\bar{\mathscr{K}}_{\operatorname{RKHS},t}$ does not depend on $t$ and the
identity
$\bar{\mathscr{K}}_{\operatorname{RKHS}}\mathscr{E}p_{t}=K\Phi({x_{t}}_{*}P)$
holds.

###### Proposition 8

If the random dynamical system $\pi$ satisfies
$\pi(t,\omega,x)=\beta_{t}(\omega,x)=h(x)+\xi_{t}(\omega)$, then the Koopman
operator $\bar{\mathscr{K}}_{t}$ in $\mathcal{H}_{k}$ does not depend on $t$
and is the adjoint operator of $K\mid_{\operatorname{Span}\\{\phi(x)\mid
x\in\mathcal{X}\\}}$.

##  4 Krylov Subspace Methods for Perron-Frobenius Operators in RKHSs

In this section, we describe the estimation problem of the Perron-Frobenius
operator $K$ defined as Eq. (2). For this purpose, we extend Krylov subspace
methods to our case. We first extend the classical Arnoldi method to our case
in Subsection 4.1. Although this method requires $K$ to be bounded for its
convergence, $K$ is not necessarily bounded even for standard situations. For
example, if $k$ is the Gaussian kernel, $h$ is nonlinear and $\xi_{t}\equiv
0$, then $K$ is unbounded (Ikeda et al., 2019). Therefore, we develop a novel
shift-invert Arnoldi method to avoid this issue in Subsection 4.2. Although
these two subsections discuss the ideal situations with infinite length of
time-series data, we consider practical situations with finite ones in
Subsection 4.3.

With both methods, we construct the basis of the Krylov subspace as follows.
Let $S\in\mathbb{N}$ be the dimension of the Krylov subspace constructed using
observed time-series data $\\{\tilde{x}_{0},\tilde{x}_{1},\ldots\\}$, which is
assumed to be generated by dynamical system (1) with sample
$\omega_{0}\in\varOmega$. To generate the elements of the basis of the Krylov
subspace in terms of kernel mean embedding of probability measures, we split
the observed data into $S+1$ datasets as
$\\{\tilde{x}_{0},\tilde{x}_{S^{\prime}},\ldots\\}$,
$\\{\tilde{x}_{1},\tilde{x}_{1+S^{\prime}},\ldots\\}$, $\ldots$,
$\\{\tilde{x}_{S},\tilde{x}_{S+S^{\prime}},\ldots\\}$, where $S^{\prime}=S+1$.
Then we define each element of the basis as the time average of each subset
above in the RKHS.

###  4.1 Arnoldi method for bounded operators

For $t=0,\dots,S$, let
$\mu_{t,N}:=\frac{1}{N}\sum_{i=0}^{N-1}\delta_{\tilde{x}_{t+iS^{\prime}}}$ be
the empirical measure constructed from observed data, where $\delta_{x}$
denotes the Dirc measure centered at $x\in\mathcal{X}$, and
$\Psi_{0,N}:=[\Phi(\mu_{0,N}),\ldots,\Phi(\mu_{S-1,N})]$ with
$N\in\mathbb{N}$. By the definition of $K$, the following relation holds:

$K\Psi_{0,N}=\left[\Phi\left({\beta_{0}}_{*}\left(\mu_{0,N}\otimes
P\right)\right),\ldots,\Phi\left({\beta_{S-1}}_{*}\left(\mu_{S-1,N}\otimes
P\right)\right)\right].$ (3)

The calculation on the right-hand side of the Eq. (3) is possible only if
$\beta_{t}$ is available. However, in practical situations, $\beta_{t}$ is not
available. Therefore, $\Phi({\beta_{t}}_{*}(\mu_{t,N}\otimes P))$ is not
available either. To avoid this problem, we assume the following condition,
which is similar to ergodicity, i.e., for any measurable and integral function
$f$, the following identity holds:

$\begin{split}&\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}\int_{\omega\in\varOmega}f(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t}(\omega))\;dP(\omega)\\\
&\qquad=\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}f(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\eta))\
a.s.\ \eta\in\varOmega.\end{split}$ (4)

Here, while the left-hand side of assumption (4) represents the space average
of $\xi_{t}$, the right-hand side gives its time average. As a result,
$\lim_{N\to\infty}\Phi\left({\beta_{t}}_{*}\left(\mu_{t,N}\otimes
P\right)\right)$ can be calculated without $\beta_{t}$, which is stated as
follows:

###### Proposition 9

Under assumption (4), the following identity holds for $t=0\ldots,S-1$:

$\displaystyle\lim_{N\to\infty}\Phi\left({\beta_{t}}_{*}\left(\mu_{t,N}\otimes
P\right)\right)=\lim_{N\to\infty}\Phi(\mu_{t+1,N}).$

Proof  By the definition of $K$, the identity
$\lim_{N\to\infty}K\Phi(\mu_{t,N})=\lim_{N\to\infty}\Phi\left({\beta_{t}}_{*}(\mu_{t,N}\otimes
P\right))$ holds. Moreover, under assumption (4), the following equalities
hold:

$\displaystyle\lim_{N\to\infty}\Phi\left({\beta_{t}}_{*}\left(\mu_{t,N}\otimes
P\right)\right)=\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}\int_{\omega\in\varOmega}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t}(\omega))\;dP(\omega)$
$\displaystyle\qquad=\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega_{0}))=\lim_{N\to\infty}\Phi\left(\frac{1}{N}\sum_{i=0}^{N-1}\delta_{\tilde{x}_{t+1+iS^{\prime}}}\right)$
$\displaystyle\qquad=\lim_{N\to\infty}\Phi(\mu_{t+1,N}),$

which completes the proof of the proposition.

Assume $\mu_{t,N}$ converges weakly to a finite signed measure $\mu_{t}$.
Then, since $\Phi$ is continuous,
$\lim_{N\to\infty}\Phi(\mu_{t,N})=\Phi(\mu_{t})$ holds. Moreover, if $K$ is
bounded, then $K\Phi(\mu_{t})=\lim_{N\to\infty}K\Phi(\mu_{t,N})$ holds. By
Lemma 9, the limit of the right-hand side of Eq. (3) is represented without
$\beta_{t}$ as $[\Phi(\mu_{1}),\ldots,\Phi(\mu_{S})]$. In addition, that of
the left-hand side becomes
$K\left[\Phi(\mu_{0}),\ldots,\Phi(\mu_{S-1})\right]$. As a result, we have:

$\displaystyle
K\left[\Phi(\mu_{0}),\ldots,\Phi(\mu_{S-1})\right]=\left[\Phi(\mu_{1}),\ldots,\Phi(\mu_{S})\right].$
(5)

Note that the range of the operator $[\Phi(\mu_{0}),\ldots,\Phi(\mu_{S-1})]$
in Eq. (5) is the Krylov subspace $\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ (cf.
Section 2.3) since

$\displaystyle\operatorname{Span}\\{\Phi(\mu_{0}),\ldots,\Phi(\mu_{S-1})\\}=\operatorname{Span}\\{\Phi(\mu_{0}),K\Phi(\mu_{0}),\ldots,K^{S-1}\Phi(\mu_{0})\\}.$

Now, the estimation of $K$ is carried out as follows. First, define $\Psi_{0}$
and $\Psi_{1}$ as

$\Psi_{0}:=[\Phi(\mu_{0})\ldots,\Phi(\mu_{S-1})],~{}\Psi_{1}:=[\Phi(\mu_{1})\ldots,\Phi(\mu_{S})].$

Then, we orthonormally project $K$ to the Krylov subspace
$\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ by QR decomposition. That is, let

$\Psi_{0}=Q_{S}\mathbf{R}_{S},$

be the QR decomposition of $\Psi_{0}$, where $Q_{S}=[q_{0},\ldots,q_{S-1}]$
with an orthonormal basis $q_{0},\ldots,q_{S-1}$ of
$\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ and $\mathbf{R}_{S}$ is an $S\times S$
matrix. Note that since $(Q_{S}Q_{S}^{*})^{2}=Q_{S}Q_{S}^{*}$ and
$(Q_{S}Q_{S})^{*}=Q_{S}Q_{S}^{*}$, $Q_{S}Q_{S}^{*}$ is an orthogonal
projection, where $Q_{S}^{*}$ is the adjacent operator of $Q_{S}$. $Q_{S}$
transforms a vector in $\mathbb{C}^{S}$ into the corresponding vector in
$\mathcal{H}_{k}$, which is the linear combination of the orthonormal basis of
$\mathcal{V}_{S}(K,\Phi(\mu_{0}))$. On the other hand, $Q_{S}^{*}$, the
adjacent operator of $Q_{S}$, projects a vector in $\mathcal{H}_{k}$ to
$\mathbb{C}^{S}$. Moreover, $\mathbf{R}_{S}$ transforms the coordinate with
basis $\Phi(\mu_{0}),\ldots,\Phi(\mu_{S-1})$ into that with
$q_{0},\ldots,q_{S-1}$. By identifying the $S$-dimensional subspace
$\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ with $\mathbb{C}^{S}$, a projection of $K$
to $\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ is represented as the $S\times S$ matrix
$Q_{S}^{*}KQ_{S}$. This matrix gives a numerical approximation of $K$. Let
$\tilde{\mathbf{K}}_{S}:=Q_{S}^{*}KQ_{S}$. As a result, since
$Q_{S}=\Psi_{0}\mathbf{R}_{S}^{-1}$, the following equality is derived using
Eq. (5):

$\displaystyle\tilde{\mathbf{K}}_{S}=Q_{S}^{*}\Psi_{1}\mathbf{R}_{S}^{-1},$

which shows that $\tilde{\mathbf{K}}_{S}$ can be calculated with only observed
time series data $\\{\tilde{x}_{0},\tilde{x}_{1},\ldots\\}$. We give a more
detailed explanation of the QR decomposition for the current case and the
pseudo-code of the above in Appendices B and C, respectively.

Regarding the convergence of $\tilde{\mathbf{K}}_{S}$, we have the following
proposition:

###### Proposition 10

Assume $K$ to be bounded. If the Krylov subspace
$\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ converges to $\mathcal{H}_{k}$, that is, if
the projection $Q_{S}Q_{S}^{*}$ converges strongly to the identity map in
$\mathcal{H}_{k}$ as $S\to\infty$, then for
$v\in\Phi(\mathcal{M}(\mathcal{X}))$, $Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}v$
converges to $Kv$ as $S\to\infty$.

Proof  Since $\tilde{\mathbf{K}}_{S}=Q_{S}^{*}KQ_{S}$ and
$Kv\in\mathcal{H}_{k}$, the following inequality holds:

$\displaystyle\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}v-Kv\|_{k}$
$\displaystyle\leq\|Q_{S}Q_{S}^{*}KQ_{S}Q_{S}^{*}v-Q_{S}Q_{S}^{*}Kv\|_{k}+\|Q_{S}Q_{S}^{*}Kv-
Kv\|_{k}$
$\displaystyle\leq|||Q_{S}Q_{S}^{*}K|||_{k}\|Q_{S}Q_{S}^{*}v-v\|_{k}+\|Q_{S}Q_{S}^{*}Kv-
Kv\|_{k}$
$\displaystyle\leq|||K|||_{k}\|Q_{S}Q_{S}^{*}v-v\|_{k}+\|Q_{S}Q_{S}^{*}Kv-
Kv\|_{k}$ $\displaystyle\to 0,$ (6)

as $S\to\infty$, which completes the proof of the proposition.
Note that $|||K|||_{k}$ in Eq. (6) is not always finite if $K$ is unbounded.
Thus, Proposition 10 is not always true if $K$ is unbounded.

###  4.2 Shift-invert Arnoldi method for unbounded operators

The estimation of $K$ with the Arnoldi method does not always converge to $K$
if $K$ is unbounded. Therefore, in this section, we develop the shift-invert
Arnoldi method for the estimation of $K$ to avoid this issue. With this
method, we fix $\gamma\notin\varLambda(K)$ and consider a bounded bijective
operator $(\gamma I-K)^{-1}:\mathcal{H}_{k}\to\Phi(\mathcal{M}(\mathcal{X}))$,
where $\varLambda(K)$ is the spectrum of $K$ under the assumption
$\varLambda(K)\neq\mathbb{C}$. And, bounded operator $(\gamma I-K)^{-1}$
instead of $K$ is projected to a Krylov subspace.

For the projection of $(\gamma I-K)^{-1}$, we need to calculate the Krylov
subspace of $(\gamma I-K)^{-1}$. However, since $K$ is unknown in the current
case, directly calculating $(\gamma I-K)^{-i}v$ thus, the Krylov subspace is
intractable. Therefore, we construct the Krylov subspace using only data by
setting an appropriate vector $w_{S}\in\mathcal{H}_{k}$, which depends on the
dimension of the Krylov subspace $S$. That is, the following proposition
guarantees a similar identity to Eq. (5):

###### Proposition 11

Define $w_{j}:=\sum_{t=0}^{j}\binom{j}{t}(-1)^{t}\gamma^{j-t}\Phi(\mu_{t})$.
Then, we have:

$\displaystyle(\gamma
I-K)^{-1}\left[w_{1},\ldots,w_{S}\right]=\left[w_{0},\dots,w_{S-1}\right].$

Moreover, space $\operatorname{Span}\\{w_{1},\dots,w_{S}\\}$ is the Krylov
subspace $\mathcal{V}_{S}\left((\gamma I-K)^{-1},w_{S}\right)$.

Proof  Based on Proposition 9, we have:

$\displaystyle\lim_{N\to\infty}(\gamma\Phi(\mu_{t,N})-\Phi(\mu_{t+1,N}))=\lim_{N\to\infty}\left(\gamma\Phi(\mu_{t,N})-\Phi\left(\left({\beta_{t}}_{*}\left(\mu_{t,N}\otimes
P\right)\right)\right)\right)$
$\displaystyle\quad=\lim_{N\to\infty}\left(\gamma\Phi(\mu_{t,N})-K\Phi(\mu_{t,N})\right)=\lim_{N\to\infty}(\gamma
I-K)\Phi(\mu_{t,N}).$ (7)

Since $(\gamma I-K)^{-1}$ is bounded, applying $(\gamma I-K)^{-1}$ to both
sides of Eq. (7) derives the identity $(\gamma
I-K)^{-1}(\gamma\Phi(\mu_{t})-\Phi(\mu_{t+1}))=\Phi(\mu_{t})$. Thus, for
$j=0,\ldots,S-1$, the following identity holds:

$\displaystyle(\gamma
I-K)^{-1}\sum_{t=0}^{j}\binom{j}{t}(-1)^{t}\gamma^{j-t}(\gamma\Phi(\mu_{t})-\Phi(\mu_{t+1}))=\sum_{t=0}^{j}\binom{j}{t}(-1)^{t}\gamma^{j-t}\Phi(\mu_{t}).$
(8)

Since $\binom{j}{t}+\binom{j}{t-1}=\binom{j+1}{t}$, the following identities
also hold:

$\displaystyle\sum_{t=0}^{j}\binom{j}{t}(-1)^{t}\gamma^{j-t}(\gamma\Phi(\mu_{t})-\Phi(\mu_{t+1}))$
$\displaystyle\quad=\sum_{t=0}^{j}\binom{j}{t}(-1)^{t}\gamma^{j+1-t}\Phi(\mu_{t})+\sum_{t=0}^{j}\binom{j}{t}(-1)^{t+1}\gamma^{j-t}\Phi(\mu_{t+1})$
$\displaystyle\quad=\binom{j}{0}\gamma^{j+1}\Phi(\mu_{0})+\sum_{t=1}^{j}\left(\binom{j}{t}+\binom{j}{t-1}\right)(-1)^{t}\gamma^{j+1-t}\Phi(\mu_{t})+\binom{j}{j}(-1)^{j+1}\gamma^{0}\Phi(\mu_{j+1})$
$\displaystyle\quad=\sum_{t=0}^{j+1}\binom{j+1}{t}(-1)^{t}\gamma^{j+1-t}\Phi(\mu_{t}).$
(9)

Since $w_{j}=\sum_{t=0}^{j}\binom{j}{t}(-1)^{t}\gamma^{j-t}\Phi(\mu_{t})$, by
Eqs. (8) and (9), the identity $(\gamma I-K)^{-1}w_{j+1}=w_{j}$ holds. Thus
the following identity holds:

$\displaystyle(\gamma
I-K)^{-1}\left[w_{1},\ldots,w_{S}\right]=\left[w_{0},\dots,w_{S-1}\right],$

and space $\operatorname{Span}\\{w_{1},\dots,w_{S}\\}$ is the Krylov subspace
$\mathcal{V}_{S}((\gamma I-K)^{-1},w_{S})$.
Note that $w_{j}$ can be calculated using only data.

We now describe the estimation procedure. First, define $\Psi_{0}$ and
$\Psi_{1}$ as

$\Psi_{0}:=[w_{1},\ldots,w_{S}],~{}\Psi_{1}:=[w_{0},\ldots,w_{S-1}],$

respectively. And, let $\Psi_{0}=Q_{S}\mathbf{R}_{S}$ be the QR decomposition
of $\Psi_{0}$. Similar to the Arnoldi method, the projection of $(\gamma
I-K)^{-1}$ to $\mathcal{V}_{S}((\gamma I-K)^{-1},w_{S})$ is formulated as

$\tilde{\mathbf{L}}_{S}:=Q_{S}^{*}\Psi_{1}\mathbf{R}_{S}^{-1},$

by using Proposition 11. As a result, $K$ is estimated by transforming the
projected $(\gamma I-K)^{-1}$ back into $K$ as

$\tilde{\mathbf{K}}_{S}:=\gamma\mathbf{I}-\tilde{\mathbf{L}}_{S}^{-1}$

A more detailed explanation of the QR decomposition for the current case and
the pseudo-code are found in Appendices B and C, respectively.

Regarding the convergence of $\tilde{\mathbf{K}}_{S}$, we have the following
proposition:

###### Proposition 12

If the Krylov subspace $\mathcal{V}_{S}((\gamma I-K)^{-1},w_{S})$ converges to
$\mathcal{H}_{k}$, that is, if the projection $Q_{S}Q_{S}^{*}$ converges
strongly to the identity map in $\mathcal{H}_{k}$ as $S\to\infty$, and if the
convergence of $Q_{S}Q_{S}^{*}$ to the identity map is faster than the
increase in $|||Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}(\gamma I-K)^{-1}|||_{k}$
along $S$, i.e., $|||Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}(\gamma
I-K)^{-1}|||_{k}\|Q_{S}Q_{S}^{*}u-u\|_{k}\to 0$ as $S\to\infty$ for arbitrary
$u\in\mathcal{H}_{k}$, then for $v\in\Phi(\mathcal{M}(\mathcal{X}))$,
$Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}v$ converges to $Kv$ as $S\to\infty$.

Proof  Since
$\tilde{\mathbf{K}}_{S}=\gamma\mathbf{I}-\tilde{\mathbf{L}}_{S}^{-1}$ and
$\tilde{\mathbf{L}}_{S}=Q_{S}^{*}(\gamma I-K)^{-1}Q_{S}$, and since
$v\in\Phi(\mathcal{M}(\mathcal{X}))$ can be represented as $v=(\gamma
I-K)^{-1}u$ with some $u\in\mathcal{H}_{k}$ by the bijectivity of $(\gamma
I-K)^{-1}$, the following inequality holds:

$\displaystyle\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}v-Kv\|_{k}$
$\displaystyle\qquad\leq\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}(\gamma
I-K)^{-1}u-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}(\gamma
I-K)^{-1}Q_{S}Q_{S}^{*}u\|_{k}$
$\displaystyle\qquad\qquad\qquad+\|Q_{S}g(\tilde{\mathbf{L}}_{S})Q_{S}^{*}u-g((\gamma
I-K)^{-1})u\|_{k}$
$\displaystyle\qquad\leq|||Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}(\gamma
I-K)^{-1}|||_{k}\|Q_{S}Q_{S}^{*}u-u\|_{k}$
$\displaystyle\qquad\qquad\qquad+|\gamma|\|Q_{S}\tilde{\mathbf{L}}_{S}Q_{S}^{*}u-(\gamma
I-K)^{-1}u\|_{k}+\|Q_{S}Q_{S}^{*}u-u\|_{k},$ (10)

where $g(z):=\gamma z-1$. Since
$u\in\mathcal{H}_{k}=\operatorname*{\mathcal{D}}((\gamma I-K)^{-1})$,
$Q_{S}\tilde{\mathbf{L}}_{S}Q_{S}u$ converges to $(\gamma I-K)^{-1}u$ in the
same manner as Proposition 10. Also, under the assumption of
$|||Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}(\gamma
I-K)^{-1}|||_{k}\|Q_{S}Q_{S}^{*}u-u\|_{k}\to 0$ as $S\to\infty$,
$\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}v-Kv\|_{k}\to 0$ as $S\to\infty$.
Note that since $Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}$ and $(\gamma I-K)^{-1}$
are bounded, $|||Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}(\gamma I-K)^{-1}|||_{k}$
in the first term of the last inequality of Eq. (10) is finite for a fixed
$S$. This situation is completely different from that of the Arnoldi method,
in which case $|||Q_{S}Q_{S}^{*}K|||_{k}$ in Eq. (6) cannot always be defined
when $K$ is unbounded even if $S$ is fixed. Concerning the second term of the
last inequality of Eq. (10), it represents the approximation error of bounded
operator $(\gamma I-K)^{-1}$ by $\tilde{\mathbf{L}}_{S}$, which corresponds to
the approximation error of $K$ by $\tilde{\mathbf{K}}_{S}$ with the Arnoldi
method.

###  4.3 Computation with finite data

In practice, $\mu_{t}$ are not available due to the finiteness of data
$\\{\tilde{x}_{0},\tilde{x}_{1},\ldots\\}$. Therefore, we need $\mu_{t,N}$
instead of $\mu_{t}$. We define $\Psi_{0,N}$ and $\Psi_{1,N}$ as the
quantities that are obtained by replacing $\mu_{t}$ with $\mu_{t,N}$ in the
definitions of $\Psi_{0}$ and $\Psi_{1}$. For example, we define $\Psi_{0,N}$
for the Arnoldi method (described in Subsection 4.1) by
$\Psi_{0,N}:=[\Phi(\mu_{0,N}),\ldots,\Phi(\mu_{S-1,N})]$. Also, we let
$\Psi_{0,N}=Q_{S,N}\mathbf{R}_{S,N}$ be the QR decomposition of $\Psi_{0,N}$,
and $\tilde{\mathbf{K}}_{S,N}$ be the estimator with $\Psi_{0,N}$ and
$\Psi_{1,N}$ that corresponds to $\tilde{\mathbf{K}}_{S}$.

Then, we can show that the above matrices from finite data converge to the
original approximators.

###### Proposition 13

As $N\rightarrow\infty$, the matrix $\tilde{\mathbf{K}}_{S,N}$ converges to
matrix $\tilde{\mathbf{K}}_{S}$, and operator
$Q_{S,N}:\mathbb{C}^{S}\to\mathcal{H}_{k}$ converges to $Q_{S}$ strongly in
$\mathcal{H}_{k}$.

Proof  The elements of $\mathbf{R}_{S,N}\in\mathbb{C}^{S\times S}$ and
$Q_{S,N}^{*}\Psi_{1,N}\in\mathbb{C}^{S\times S}$ are composed of the finite
linear combinations of the inner products between $\Phi(\mu_{t,N})$ in the
RKHS. Since $\lim_{N\to\infty}\mu_{t,N}\to\mu_{t}$ for each
$t\in\\{0,\ldots,S\\}$, and since $\Phi$ is continuous, the identity
$\lim_{N\to\infty}\Phi(\mu_{t,N})\to\Phi(\mu_{t})$ holds. Therefore, by the
continuity of the inner product $\left\langle\cdot,\cdot\right\rangle_{k}$,
$\left\langle\Phi(\mu_{t,N}),\Phi(\mu_{s,N})\right\rangle_{k}$ converges to
$\left\langle\Phi(\mu_{t}),\Phi(\mu_{s})\right\rangle_{k}$ for each
${}^{\forall}t,s\in\\{0,\ldots,S\\}$ as $N\to\infty$. Thus, matrices
$\mathbf{R}_{S,N}$ and $Q_{S,N}^{*}\Psi_{1,N}$ converge to $\mathbf{R}_{S}$
and $Q_{S}^{*}\Psi_{1}$ as $N\to\infty$, respectively. This implies matrix
$\tilde{\mathbf{K}}_{N,S}$ converges to $\tilde{\mathbf{K}}_{S}$ as
$N\to\infty$.

Moreover, by the identity $\lim_{N\to\infty}\Phi(\mu_{t,N})=\Phi(\mu_{t})$,
$\lim_{N\to\infty}\|\Psi_{0,N}v-\Psi_{0}v\|_{k}=0$ holds for all
$v\in\mathbb{C}^{S\times S}$. Since $k$ is bounded, for all
$t\in\\{0,\ldots,S\\}$, there exist $\tilde{C}(t)>0$ such that
$\|\Phi(\mu_{t,N})\|_{k}\leq\tilde{C}(t)$ for all $N\in\mathbb{N}$. Thus,
there exists $C>0$ such that $|||\Psi_{0,N}|||_{k}\leq C$ for all
$N\in\mathbb{N}$. Therefore, for all $v\in\mathbb{C}^{S\times S}$, it is
deduced that

$\displaystyle\|Q_{S,N}v-Q_{S}v\|_{k}=\|\Psi_{0,N}\mathbf{R}_{S,N}^{-1}v-\Psi_{0}\mathbf{R}_{S}^{-1}v\|_{k}$
$\displaystyle\qquad\leq\|\Psi_{0,N}\mathbf{R}_{S,N}^{-1}v-\Psi_{0,N}\mathbf{R}_{S}^{-1}v\|_{k}+\|\Psi_{0,N}\mathbf{R}_{S}^{-1}v-\Psi_{0}\mathbf{R}_{S}^{-1}v\|_{k}$
$\displaystyle\qquad\leq|||\Psi_{0,N}|||_{k}\|\mathbf{R}_{S,N}^{-1}v-\mathbf{R}_{S}^{-1}v\|_{k}+\|\Psi_{0,N}\mathbf{R}_{S}^{-1}v-\Psi_{0}\mathbf{R}_{S}^{-1}v\|_{k}$
$\displaystyle\qquad\leq
C\|\mathbf{R}_{S,N}^{-1}v-\mathbf{R}_{S}^{-1}v\|_{k}+\|\Psi_{0,N}\mathbf{R}_{S}^{-1}v-\Psi_{0}\mathbf{R}_{S}^{-1}v\|_{k}$
$\displaystyle\qquad\to 0,$

as $N\to\infty$. This implies that $Q_{S,N}$ converges to $Q_{S}$ strongly in
$\mathcal{H}_{k}$, which completes the proof of the proposition.
The convergence speeds of $Q_{S,N}\to Q_{S}$ and
$\tilde{\mathbf{K}}_{S,N}\to\tilde{\mathbf{K}}_{S}$ depend on that of
$\Phi(\mu_{t,N})\to\Phi(\mu_{t})$ as described in the proof of this
proposition. And, the following proposition gives the connection of the
convergence of $\Phi(\mu_{t,N})\to\Phi(\mu_{t})$ with the property of noise
$\xi_{t}$:

###### Proposition 14

For all $\epsilon>0$ and for $t=1,\ldots,S$, if $N$ is sufficiently large, the
probability of
$\|\Phi(\mu_{t,N})-\lim_{N\to\infty}\Phi(\mu_{t})\|_{k}\geq\epsilon$ under the
condition of $x_{t-1+iS^{\prime}}=\tilde{x}_{t-1+iS^{\prime}}$ is bounded by
$4\sum_{i=0}^{N-1}\sigma_{t-1,i}^{2}/(N^{2}\epsilon^{2})$, where
$\sigma_{t,i}^{2}:=\int_{\omega\in\varOmega}\left\|\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))-m_{t,i}\right\|_{k}^{2}\;dP(\omega)$
and
$m_{t,i}:=\int_{\omega\in\varOmega}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))\;dP(\omega)$.

Proof  Let
$\bar{m}_{t,N}:=1/N\sum_{i=0}^{N-1}\int_{\omega\in\varOmega}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))\;dP(\omega)$,
$m_{t,i}:=\int_{\omega\in\varOmega}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))\;dP(\omega)$
and
$\bar{\sigma}_{t,N}^{2}:=\int_{\omega\in\varOmega}\|1/N\sum_{i=0}^{N-1}\phi(h(x_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))-\bar{m}_{t,N}\|_{k}^{2}\;dP(\omega)$.
The following identities about $\bar{\sigma}_{t,N}$ hold:

$\displaystyle\bar{\sigma}_{t,N}^{2}$
$\displaystyle=\int_{\omega\in\varOmega}\left\|\frac{1}{N}\sum_{i=0}^{N-1}(\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))-m_{t,i})\right\|_{k}^{2}\;dP(\omega)$
$\displaystyle=\int_{\omega\in\varOmega}\frac{1}{N^{2}}\left(\sum_{i=0}^{N-1}\|\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))-m_{t,i}\|_{k}^{2}\right.$
$\displaystyle\quad\left.+\sum_{\begin{subarray}{c}i,j=0\\\ i\neq
j\end{subarray}}^{N-1}\left\langle\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))-m_{t,i},\phi(h(\tilde{x}_{t+jS^{\prime}})+\xi_{t+jS^{\prime}}(\omega))-m_{t,j}\right\rangle_{k}\right)\;dP(\omega)$
$\displaystyle=\frac{1}{N^{2}}\sum_{i=0}^{N-1}\sigma_{i,t}^{2}.$ (11)

The last equality holds because $\xi_{t+iS^{\prime}}$ and
$\xi_{t+jS^{\prime}}$ are independent if $i\neq j$ and the following equality
holds for $i,j\in\\{0,\ldots,N-1\\},\ i\neq j$ by definition of $m_{t,i}$:

$\displaystyle\int_{\omega\in\varOmega}\left\langle\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))-m_{t,i},\phi(h(\tilde{x}_{t+jS^{\prime}})+\xi_{t+jS^{\prime}}(\omega))-m_{t,j}\right\rangle_{k}\;dP(\omega)$
$\displaystyle=\int_{\omega\in\varOmega}\int_{\eta\in\varOmega}\\!\\!\left\langle\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))-m_{t,i},\phi(h(\tilde{x}_{t+jS^{\prime}})+\xi_{t+jS^{\prime}}(\eta))-m_{t,j}\right\rangle_{k}dP(\omega)dP(\eta)$
$\displaystyle=0.$

Let $\epsilon>0$. By Chebyshev’s inequality and Eq. (11), it is derived that:

$\displaystyle
P\left(\left\|\frac{1}{N}\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t-1+iS^{\prime}})+\xi_{t-1+iS^{\prime}}(\omega))-m_{t-1,N}\right\|_{k}\geq\epsilon\right)\leq\frac{\bar{\sigma}_{t-1,N}^{2}}{\epsilon^{2}}=\frac{1}{N^{2}\epsilon^{2}}\sum_{i=0}^{N-1}\sigma_{t-1,i}^{2}.$

By assumption (4), for sufficiently large $N$,
$P(\|\bar{m}_{t,N}-\lim_{N\to\infty}1/N\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))\|_{k}\leq\epsilon/2)=1$
holds. Thus, we have:

$\displaystyle
P(\|\Phi(\mu_{t,N})-\Phi(\mu_{t})\|_{k}\geq\epsilon\mid\tilde{x}_{t-1+iS^{\prime}}=\tilde{x}_{t-1+iS^{\prime}},\
i=0,1,\ldots)$
$\displaystyle=P\left(\left\|\frac{1}{N}\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega)-\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega)\right\|_{k}\geq\epsilon\right)$
$\displaystyle\leq
P\left(\left\|\frac{1}{N}\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega)-\bar{m}_{t,N}\right\|_{k}\right.$
$\displaystyle\qquad\qquad\qquad\qquad+\left.\left\|\bar{m}_{t,N}-\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega))\right\|_{k}\geq\epsilon\right)$
$\displaystyle\leq
P\left(\left\|\frac{1}{N}\sum_{i=0}^{N-1}\phi(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega)-\bar{m}_{t,N}\right\|_{k}+\frac{\epsilon}{2}\geq\epsilon\right)\leq\frac{4}{N^{2}\epsilon^{2}}\sum_{i=0}^{N-1}\sigma_{t-1,i}^{2},$

which completes the proof of the proposition.
Note that $\sigma_{t,i}$ represents the variance of
$h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega)$ in RKHS and
$\xi_{t}(\omega)$ is the only term that depends on $\omega$ in
$h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega)$. Therefore,
$\sigma_{t,i}$ is small if the variance of $\xi_{t}$ is small. Thus,
Proposition 14 shows if the variance of $\xi_{t}$ is small, the convergence is
fast.

##  5 Connection to Existing Methods

In the previous two sections, we defined a Perron-Frobenius operator for
dynamical systems with random noise based on kernel mean embeddings and
developed the Krylov subspace methods for estimating it. We now summarize the
connection of the methods with the existing Krylov subspace methods for
transfer operators on dynamical systems.

For $\mathcal{L}^{2}(\mathcal{X})$ space and deterministic dynamical systems,
the Arnoldi method for the Krylov subspace with Koopman operator
$\mathscr{K}$,

$\mathcal{V}_{S}(\mathscr{K},g)=\operatorname{Span}\\{g,\mathscr{K}g,\ldots,\mathscr{K}^{S-1}g\\},$

is considered in Kutz (2013), where $g:\mathcal{X}\to\mathbb{C}$ is an
observable function. Let $\\{\tilde{x}_{0},\ldots,\tilde{x}_{S-1}\\}$ be the
sequence generated from deterministic system $x_{t+1}=h(x_{t})$. Then, this
Krylov subspace captures the time evolution starting from many initial values
$\tilde{x}_{0},\ldots,\tilde{x}_{S-1}$ by approximating $g$ with
$[g(\tilde{x}_{0}),\ldots,g(\tilde{x}_{S-1})]$. This idea is extended to the
Krylov subspace with Koopman operator $\bar{\mathscr{K}}$,
$\mathcal{V}_{S}(\bar{\mathscr{K}},g)$, for the case in which the system is
random, by assuming the following ergodicity (Črnjarić-Žic et al., 2017;
Takeishi et al., 2017a): for any measurable function $f$,

$\int_{x\in\mathcal{X}}f(x)\;d\mu(x)=\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}f(\tilde{x}_{i}).$
(12)

The following proposition states the connection with our assumption (4) in
Subsection 4.1 and the assumption (12).

###### Proposition 15

For each $t=0,\ldots,S-1$, if there exists a random variable $y_{t}$ such that
$\mu_{t}={y_{t}}_{*}P$, and $y_{t}$ is independent of $\xi_{t}$, then
assumption (4) is equivalent to assumption (12).

Meanwhile, Kawahara considers a Perron-Frobenius operator for deterministic
systems in an RKHS $\mathcal{H}_{k}$, and projects it to the following Krylov
subspace (Kawahara, 2016):

$\displaystyle\mathcal{V}_{S}(\mathscr{K}_{\operatorname{RKHS}},\phi(\tilde{x}_{0}))=\operatorname{Span}\\{\phi(\tilde{x}_{0}),\mathscr{K}_{\operatorname{RKHS}}\phi(\tilde{x}_{0})\ldots,\mathscr{K}_{\operatorname{RKHS}}^{S-1}\phi(\tilde{x}_{0})\\}.$
(13)

Subspace (13) captures the time evolution starting from a single initial value
$\tilde{x}_{0}$. This prevents the straightforward extension of Krylov
subspace (13) to the subspace that is applicable to the case in which the
dynamics is random. It can be shown that the Krylov subspace with Perron-
Frobenius operator $K$, $\mathcal{V}_{S}(K,\Phi(\mu_{0}))$, which is addressed
in this paper for random systems, is a generalization of the Krylov subspace
for the deterministic systems considered in Kawahara (2016):

###### Proposition 16

The Krylov subspace $\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ generalizes the Krylov
subspace introduced by Kawahara (13) to that for dynamical systems with random
noise.

Note that the framework of the Krylov subspace methods for Perron-Frobenius
operators for random systems has not been addressed in prior works. Also note
that the theoretical analysis for these methods requires the assumption that
the operator is bounded, which is not necessarily satisfied for transfer
operators on discrete-time nonlinear systems (Ikeda et al., 2019).

The shift-invert Arnoldi method is a popular Krylov subspace method discussed
in the field of numerical linear algebra, which is applied to extract
information from given matrices, and some theoretical analyses have been
extended to given unbounded operators (Güttel, 2010; Grimm, 2012; Göckler,
2014; Hashimoto and Nodera, 2018). However, as far as we know, our paper is
the first paper to address the unboundedness of Perron-Frobenius operators for
the estimation problem and apply the shift-invert Arnorldi method to estimate
Perron-Frobenius operators, which are not known beforehand. Since the shift-
invert Arnoldi method was originally investigated for given matrices or
operators, applying it to unknown operators is not straightforward, as
described in Subsection 4.2.

##  6 Evaluation of Prediction Errors with Estimated Operators

In this section, we discuss an approach of evaluating the prediction accuracy
with estimated Perron-Frobenius operators, which is applicable, for example,
to anomaly detection in complex systems.

Consider the prediction of $\phi(\tilde{x}_{t})$ ($\in\mathcal{H}_{k}$) using
estimated operator $\tilde{\mathbf{K}}_{S,N}$ and observation (embedded in
RKHS $\mathcal{H}_{k}$) $\phi(\tilde{x}_{t-1})$. This prediction is calculated
as $Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})$. Thus,
the prediction error can be evaluated as

$\|\phi(\tilde{x}_{t})-Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})\|_{k}.$
(14)

Note that this is the maximum mean discrepancy (MMD) between
$Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})$ and
$\phi(\tilde{x}_{t})$ in the unit disk in $\mathcal{H}_{k}$ Gretton et al.
(2012).

For practical situations such as anomaly detection, we define the degree of
abnormality for prediction at $t$ based on the MMD (14) as follows:

$a_{t,S}:=\frac{\|\phi(\tilde{x}_{t})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}{\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}.$

The $a_{t,S}$ is bounded as follows:

$\displaystyle
a_{t,S}=\frac{\|\phi(\tilde{x}_{t})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}{\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}$
$\displaystyle\qquad\leq\frac{\|\phi(\tilde{x}_{t})-K\phi(\tilde{x}_{t-1})\|_{k}}{\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}+\frac{\|K\phi(\tilde{x}_{t-1})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}{\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}.$
(15)

Concerning the second term of the right-hand side in Eq. (15), the following
proposition is derived directly by Propositions 10 and 12:

###### Proposition 17

Let $\tilde{\mathbf{K}}_{S}$ be the estimation using the shift-invert Arnoldi
method. Under the assumption of Proposition 12,
${\|K\phi(\tilde{x}_{t-1})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}/{\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}$,
the second term of the right-hand side in Eq. (15), converges to $0$ as
$S\to\infty$. For the Arnoldi method, the convergence is attained for the case
in which $K$ is bounded and the Krylov subspace
$\mathcal{V}_{S}(K,\Phi(\mu_{0}))$ converges to $\mathcal{H}_{k}$.

On the other hand, the numerator of the first term of the right-hand side in
Eq. (15) represents the deviation of the observation $\tilde{x}_{t}$ from the
prediction at $t$ under the assumption that $\tilde{x}_{t}$ is generated by
the dynamical system (1) in the RKHS, because the identity
$K\phi(\tilde{x}_{t-1})=\int_{\omega\in\varOmega}\phi(h(\tilde{x}_{t-1})+\xi_{t}(\omega))\;dP(\omega)$
holds by the definition of $K$. And, the following proposition shows that the
denominator indicates how vector $\phi(\tilde{x}_{t-1})$ deviates from the
Krylov subspace.

###### Proposition 18

Let $\tilde{\mathbf{K}}_{S}$ be the estimation with the shift-invert Arnoldi
method and let $g(z):=z/(\gamma z-1)$ for $z\in\mathbb{C}$. If $g$ is
holomorphic in the interior of $\mathcal{W}((\gamma I-K)^{-1})$ and continuous
in $\overline{\mathcal{W}((\gamma I-K)^{-1})}$, then the following inequality
holds:

$\displaystyle\frac{1}{\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}\leq\frac{C}{\|Q_{S}^{*}\phi(\tilde{x}_{t-1})\|},$
(16)

where $C=(1+\sqrt{2})\sup_{z\in\mathcal{W}((\gamma I-K)^{-1})}|g(z)|\geq 0$ is
a constant. For the Arnoldi method, inequality (16) is satisfied with
$g(z)=1/z$ for the case in which $K$ is bounded.

To show Proposition 18, the following lemma by Crouzeix and Palencia (2017) is
used.

###### Lemma 19

Let $\mathbf{A}$ be a matrix. If $f$ is holomorphic in the interior of
$\mathcal{W}(\mathbf{A})$ and continuous in
$\overline{\mathcal{W}(\mathbf{A})}$, then there exists $0<C\leq 1+\sqrt{2}$
such that

$|||f(\mathbf{A})|||\leq C\sup_{z\in\mathcal{W}(\mathbf{A})}|f(z)|.$ (17)

Proof (Proof of Proposition 18) Let
$s_{\operatorname{min}}(\mathbf{A}):=\min_{\|\mathbf{w}\|_{=}1}\|\mathbf{Aw}\|$
be the minimal singular value of a matrix $\mathbf{A}$. Since the relation
$1/s_{\operatorname{min}}(\tilde{\mathbf{L}}_{S})=|||\tilde{\mathbf{L}}_{S}^{-1}|||$
and inclusion $\mathcal{W}(\tilde{\mathbf{L}}_{S})\subseteq\mathcal{W}((\gamma
I-K)^{-1})$ hold, the following inequalities hold:

$\displaystyle\frac{1}{\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}}=\frac{1}{\|\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|}\leq\frac{1}{s_{\operatorname{min}}(\tilde{\mathbf{K}}_{S})\|Q_{S}^{*}\phi(\tilde{x}_{t-1})\|}$
$\displaystyle\qquad\leq\frac{|||\tilde{\mathbf{K}}_{S}^{-1}|||}{\|Q_{S}^{*}\phi(\tilde{x}_{t-1})\|}\leq\frac{(1+\sqrt{2})\sup_{z\in\mathcal{W}(\tilde{\mathbf{L}}_{S})}|g(z)|}{\|Q_{S}^{*}\phi(\tilde{x}_{t-1})\|}$
$\displaystyle\qquad\leq\frac{(1+\sqrt{2})\sup_{z\in\mathcal{W}((\gamma
I-K)^{-1})}|g(z)|}{\|Q_{S}^{*}\phi(\tilde{x}_{t-1})\|}.$

For the Arnoldi method, if $K$ is bounded, the identity
$Q_{S}^{*}KQ_{S}=\tilde{\mathbf{K}}_{S}$ holds. Thus, the inequality
$|||\tilde{\mathbf{K}}_{S}^{-1}|||\leq(1+\sqrt{2})\sup_{z\in\mathcal{W}(K)}|g(z)|$
with $g(z)=1/z$ holds in this case, which deduces the same result as
inequality (16).
If $\phi(\tilde{x}_{t-1})$ deviates from the Krylov subspace, the norm of the
projected vector $Q_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})$, which is equal to
$\|Q_{S}^{*}\phi(\tilde{x}_{t-1})\|$, becomes small. Proposition 18 implies
$a_{t,S}$ becomes large in this case. On the other hand, if
$\phi(\tilde{x}_{t})$ is sufficiently close to the Krylov subspace, that is,
$\min_{u\in\mathcal{V}_{S}}\|\phi(\tilde{x}_{t-1})-u\|_{k}\approx 0$, then we
have:

$\displaystyle 1$
$\displaystyle=\|\phi(\tilde{x}_{t-1})\|_{k}^{2}=\|\phi(\tilde{x}_{t-1})-Q_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}^{2}+\|Q_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}^{2}$
$\displaystyle=\min_{u\in\mathcal{V}_{S}}\|\phi(\tilde{x}_{t-1})-u\|_{k}+\|Q_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}^{2}\approx\|Q_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}^{2}=\|Q_{S}^{*}\phi(\tilde{x}_{t-1})\|^{2},$

As a result, if $\tilde{x}_{t}$ is generated by dynamical system (1), and if
$\phi(\tilde{x}_{t-1})$ is sufficiently close to $\mathcal{V}_{S}$, then
$a_{t,S}$ is bounded by a reasonable value. Conversely, if $\tilde{x}_{t}$ is
unlikely to be generated by dynamical system (1), or $\phi(\tilde{x}_{t-1})$
is not close to the subspace $\mathcal{V}_{S}$, then $a_{t,S}$ becomes large.
In the context of anomaly detection, since both the above cases mean
$\tilde{x}_{t-1}$ or $\tilde{x}_{t}$ deviates from the regular pattern of
times-series $\\{\tilde{x}_{0},\ldots,\tilde{x}_{T-1}\\}$, they should be
regarded as abnormal.

In practice, $Q_{S}$, $\mathbf{R}_{S}$, and $\tilde{\mathbf{K}}_{S}$ are
approximated by $Q_{S,N}$, $\mathbf{R}_{S,N}$ and $\tilde{\mathbf{K}}_{S,N}$,
respectively. Thus, the following empirical value can be used:

$a_{t,S,N}:=\frac{\|\phi(\tilde{x}_{t})-Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})\|_{k}}{\|Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})\|_{k}}.$

By Proposition 13, the following proposition about the convergence of
$a_{t,S,N}$ holds:

###### Proposition 20

The $a_{t,S,N}$ converges to $a_{t,S}$ as $N\to\infty$.

##  7 Numerical Results

We empirically evaluate the behavior of the proposed Krylov subspace methods
in Subsection 7.1 then describe their application to anomaly detection using
real-world time-series data in Subsection 7.2.

###  7.1 Comparative Experiment

The behavior of the Arnoldi and shift-invert Arnoldi methods (SIA in the
figures) were evaluated numerically based on the empirical abnormality. We
used 100 synthetic time-series datasets
$\\{\tilde{x}_{0},\ldots,\tilde{x}_{T}\\}$ randomly generated by the following
dynamical system:

$x_{t+1}=0.99x_{t}\cos(0.1x_{t})+\xi_{t},$ (18)

where $x_{0}=0.5$, and $\\{\xi_{t}\\}$ is i.i.d with
$\xi_{t}\sim\mathcal{N}(0,0.01)$. Using the synthetic data, $K$ was first
estimated, then the empirical abnormalities $a_{t,S,N}$ were computed using
all time-series data with $t=1601$, $N=50,100$ and $S=1,\ldots,12$. We chose
time points $t=1601$ for evaluation because the estimation of $K$ requires
$\\{\tilde{x}_{0},\ldots,\tilde{x}_{N\times(S+1)}\\}$ and $1601>N\times(S+1)$
for all $N=50,100$ and $S=1,\ldots,12$. The Gaussian kernel was used, and
$\gamma=1+1\mathrm{i}$, where $\mathrm{i}$ denotes the imaginary unit, was set
for the shift-invert Arnoldi method.

For evaluating the behavior of each method along with $S$, the values

$|a_{t,S,N}-a_{t,S-1,N}|$

for $S=2,\ldots,12$ were computed with all time-series data, then the averages
of all the time-series data were computed.

The results are shown in Figure 1. If $K$ is bounded, and if the Krylov
subspace $\mathcal{V}(K,\Phi(\mu_{0}))$ converges to $\mathcal{H}_{k}$, that
is, $Q_{S}Q_{S}^{*}$ converges strongly to the identity map in
$\mathcal{H}_{k}$ as $S\to\infty$, then by Proposition 17, $a_{t,S}$ computed
with the Arnoldi method converges to
$\|\phi(x_{t})-K\phi(x_{t-1})\|_{k}/\|K\phi(x_{t-1})\|_{k}$. Therefore, in
this case, $|a_{t,S,N}-a_{t,S-1,N}|$, the difference between the empirical
abnormality with $S$ and that with $S-1$, becomes smaller as $S$ grows.
However, the value $|a_{t,S,N}-a_{t,S-1,N}|$ computed with the Arnoldi method
do not seem to decrease even if $S$ grows. This is due to the unboundedness of
$K$. Meanwhile, we can see $|a_{t,S,N}-a_{t,S-1,N}|$ computed with the shift-
invert Arnoldi method decreases as $S$ grows. This is because $a_{t,S}$
computed with the shift-invert Arnoldi method converges to
$\|\phi(x_{t})-K\phi(x_{t-1})\|_{k}/\|K\phi(x_{t-1})\|_{k}$ even if $K$ is
unbounded. The result indicates that the shift-invert Arnoldi method
counteracts the unboundedness of Perron-Frobenius operators.

Figure 1: Convergence of empirical abnormality along $S$

(a) $p=15$

(b) $p=30$

Figure 2: Accuracy versus false-alarm rate in anomaly-detection experiments
using our shift-invert Arnoldi method and the existing methods: LSTM and AR
model

###  7.2 Anomaly detection with real-world data

We show the empirical results for our shift-invert Arnoldi method in anomaly
detection with real-world healthcare data. We used electrocardiogram (ECG)
data (Keogh et al., 2005).111Available at ‘
http://www.cs.ucr.edu/~eamonn/discords/ ’.
We used ‘chfdb_chf01_275.txt’, ‘chfdb_chf13_45590.txt’ and
‘mitdbx_mitdbx_108.txt’ in the experiment. ECGs are time-series of the
electrical potential between two points on the surface of the body caused by a
beating heart.

The graphs in Figure 2 show the accuracy versus the false-alarm rate for these
datasets. We first computed $\tilde{\mathbf{K}}_{S,N}$ with $S=10$, $N=40$
then computed the empirical abnormality $\hat{a}_{t,N}$ for each $t$ with the
shift-invert Arnoldi method. The Laplacian kernel and $\gamma=1.25$ were used.
To extract the relationship between $\tilde{x}_{t}$ and
$\tilde{x}_{t-p+1},\ldots,\tilde{x}_{t-1}$, we set $x_{t}$ in dynamical system
(1) as $x_{t}:=[y_{t},\ldots,y_{t-p+1}]$ for random variable $y_{t}$ at $t$.
In this example, $p$ was set as $p=15,30$. $\tilde{\mathbf{K}}_{S,N}$ was
computed using the data $\\{\tilde{x}_{0},\ldots,\tilde{x}_{399}\\}$, and then
the empirical abnormalities $\hat{a}_{t,N}$ at $t=430,431,\ldots$ were
computed. Also, the results obtained using long short-term memory (LSTM)
(Malhotra et al., 2015) and autoregressive (AR) model (Takeuchi and Yamanishi,
2006) were evaluated for comparison. LSTM with $15,30$ time-series and $10$
neurons, with the tanh activation function and the AR model
$x_{t+1}=\sum_{i=0}^{p-1}c_{i}x_{t-i}+\xi_{t}$ with $p=15,30$ were used. The
datasets included 12 abnormal parts. As can be seen, our shift-invert Arnoldi
method achieved higher accuracy than LSTM and the AR model while maintaining a
low false-alarm rate for these datasets from complex systems.

##  8 Conclusions

In this paper, we addressed a transfer operator to deal with nonlinear
dynamical systems with random noise, and developed novel Krylov subspace
methods for estimating the transfer operators from finite data. For this
purpose, we first considered the Perron-Frobenius operators with kernel-mean
embeddings for such systems. As for the estimation, we extended the Arnoldi
method so that it can be applied to the current case. Then, we developed the
shift-invert Arnoldi method to avoid the problem of the unboundedness of
estimated operators because transfer operators on nonlinear systems are not
necessarily bounded. We also considered an approach of evaluating the
prediction accuracy by estimated operators on the basis of the maximum mean
discrepancy. Finally, we empirically investigated the performance of our
methods using synthetic and real-world healthcare data.

##  A Proofs

### Proof of Lemma 4

Since $\xi_{t}$ with $t\in\mathbb{Z}_{\geq 0}$ are i.i.d. and independent of
$x_{t}$, the following identities are derived:

$\displaystyle K\Phi({x_{t}}_{*}P)=\Phi({\beta_{t}}_{*}({x_{t}}_{*}P\otimes
P))=\int_{\omega\in\varOmega}\int_{x\in\mathcal{X}}\phi(\beta_{t}(x,\omega))\;d{x_{t}}_{*}P(x)\;dP(\omega)$
$\displaystyle\quad=\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}\phi(h(x)+y)\;d{x_{t}}_{*}(x)\;d{\xi_{t}}_{*}(y)=\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}\phi(x+y)\;dh(x_{t})_{*}(x)\;d{\xi_{t}}_{*}(y)$
$\displaystyle\quad=\int_{x\in\mathcal{X}}\phi(x)\;d(h(X)+\xi_{t})_{*}P(x)=\Phi((h(x_{t})+\xi_{t})_{*}P)=\Phi({x_{t+1}}_{*}P),$

which completes the proof of the lemma.

### Proof of Lemma 5

The linearity of $\Phi$ is verified by the definition of $\Phi$. Next let
$\\{\mu_{N}\\}_{N=1}^{\infty}$ be a sequence in $\mathcal{M}(\mathcal{X})$
such that $\mu=\lim_{N\to\infty}\mu_{N}$ weakly. Then since $k$ is bounded and
continuous, the following relations hold:

$\displaystyle\|\Phi(\mu_{N})-\Phi(\mu)\|_{k}^{2}=\left\langle\Phi\left(\mu_{N}\right),\Phi\left(\mu_{N}\right)\right\rangle_{k}-2\Re\left\langle\Phi\left(\mu_{N}\right),\Phi\left(\mu\right)\right\rangle_{k}+\left\langle\Phi\left(\mu\right),\Phi\left(\mu\right)\right\rangle_{k}$
$\displaystyle\qquad=\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}k(x,y)\;d\mu_{N}(x)\;d\mu_{N}(y)-2\Re\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}k(x,y)\;d\mu(x)\;d\mu_{N}(y)$
$\displaystyle\phantom{\qquad=\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}k(x,y)\;d\mu_{N}(x)\;d\mu_{N}(y)}+\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}k(x,y)\;d\mu(x)\;d\mu(y)$
$\displaystyle\qquad\longrightarrow 0,$

as $N\rightarrow\infty$, which implies
$\lim_{N\to\infty}\Phi(\mu_{N})=\Phi(\mu)$ in $\mathcal{H}_{k}$. This
completes the proof of the lemma.

### Proof of Lemma 6

Since each $\xi_{t}$ for $t\in\mathbb{Z}_{\geq 0}$ is i.i.d. and $K\Phi(\mu)$
is represented as
$K\Phi(\mu)=\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}\phi(h(x)+y)\;d\mu(x)\;d{\xi_{t}}_{*}(y)$,
$K$ does not depend on $t$.

In addition, the identity $\Phi(\delta_{x})=\phi(x)$ holds for any
$x\in\mathcal{X}$, where $\delta_{x}$ is the Dirac measure centered at
$x\in\mathcal{X}$. Thus, the inclusion $\operatorname{Span}\\{\phi(x)\mid
x\in\mathcal{X}\\}\subseteq\Phi(\mathcal{M}(\mathcal{X}))$ holds, which
implies $\Phi(\mathcal{M}(\mathcal{X}))$ is dense in $\mathcal{H}_{k}$.
Moreover, according to Sriperumbudur et al. (2011), $\Phi$ is injective for
$c_{0}$-universal kernel $k$. Therefore, the well-definedness of Perron-
Frobenius operator $K$, defined as Eq. (2) is verified.

Concerning the linearity of $K$, let $c_{1},c_{2}\in\mathbb{C}$ and
$\mu,\nu\in\mathcal{M}(\mathcal{X})$. By the linearity of $\Phi$ and the
definition of $K$, the following identities hold:

$\displaystyle
K\left(c_{1}\Phi(\mu)+c_{2}\Phi(\nu)\right)=K\Phi(c_{1}\mu+c_{2}\nu)=\Phi({\beta_{t}}_{*}((c_{1}\mu+c_{2}\nu)\otimes
P))$ $\displaystyle\qquad=\Phi(c_{1}{\beta_{t}}_{*}(\mu\otimes
P)+c_{2}{\beta_{t}}_{*}(\nu\otimes P))=c_{1}\Phi({\beta_{t}}_{*}(\mu\otimes
P))+c_{2}\Phi({\beta_{t}}_{*}(\nu\otimes P))$
$\displaystyle\qquad=c_{1}K\Phi(\mu)+c_{2}K\Phi(\nu),$

which completes the proof of the lemma.

### Proof of Proposition 7

Let $P_{\mid{x_{t}=x}}$ be a probability measure on $(\varOmega,\mathcal{F})$
satisfying $P_{\mid{x_{t}=x}}(B)=P(B\mid\\{x_{t}=x\\})$ for $B\in\mathcal{F}$.
Since $p_{t}$ is the probability density function of $x_{t}$, the identity
$\int_{x\in B}\;d{x_{t}}_{*}P(x)=\int_{x\in B}p_{t}(x)\;d\mu(x)$ holds for any
$B\in\mathcal{B}$. Moreover, by the definitions of $p$ and
$P_{\mid{x_{t}=x}}$, the equality $\int_{y\in
B}\;d({x_{t+1}}_{*}P_{\mid{x_{t}=x}})(y)=\int_{y\in B}p(y\mid x)\;d\mu(y)$
holds for any $B\in\mathcal{B}$. Thus, the following identities are derived:

$\displaystyle\tilde{\mathscr{K}}_{\operatorname{RKHS}}\mathcal{E}p_{t}$
$\displaystyle=\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}\phi(y)p(y\mid
x)p_{t}(x)\;d\mu(y)\;d\mu(y)$
$\displaystyle=\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}\phi(y)\;d\left({x_{t+1}}_{*}P_{\mid{x_{t}=x}}\right)(y)\;d{x_{t}}_{*}P(y).$

Since $x_{t+1}=h(x_{t})+\xi_{t}$, and $x_{t}$ and $\xi_{t}$ are independent,
the following identities hold for $B\in\mathcal{B}$:

$\displaystyle\int_{y\in
B}\;d\left({x_{t+1}}_{*}P_{\mid{x_{t}=x}}\right)(y)=P\left(\left\\{x_{t+1}\in
B\right\\}\mid\\{x_{t}=x\\}\right)=P\left(\left\\{h(x_{t})+\xi_{t}\in
B\right\\}\mid\\{x_{t}=x\\}\right)$
$\displaystyle\qquad=\frac{P\left(\left\\{h(x)+\xi_{t}\in
B\right\\}\bigcap\\{x_{t}=x\\}\right)}{P(\\{x_{t}=x\\})}=\frac{P\left(\left\\{\xi_{t}\in
B-h(x)\right\\}\bigcap\\{x_{t}=x\\}\right)}{P(\\{x_{t}=x\\})}$
$\displaystyle\qquad=\frac{P\left(\left\\{\xi_{t}\in
B-h(x)\right\\}\right)P\left(\\{x_{t}=x\\}\right)}{P(\\{x_{t}=x\\})}=P\left(\left\\{h(x)+\xi_{t}\in
B\right\\}\right)$ $\displaystyle\qquad=\int_{y\in
B}\;d\left({h(x)+\xi_{t}}_{*}P\right)(y),$

where $B-h(x)$ denotes the set $\\{y=z-h(x)\mid z\in B\\}$. Therefore, by the
definition of $\beta_{t}$, the following identities are derived:

$\displaystyle\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}\phi(y)\;d\left({x_{t+1}}_{*}P_{\mid{x_{t}=x}}\right)(y)\;d{x_{t}}_{*}P(y)$
$\displaystyle\qquad=\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}\phi(y)\;d\left({h(x)+\xi_{t}}_{*}P\right)(y)\;d{x_{t}}_{*}P(y)$
$\displaystyle\qquad=\int_{x\in\mathcal{X}}\phi(x)\;d{\beta_{t}}_{*}({x_{t}}_{*}P\otimes
P)(x)=\Phi\left({\beta_{t}}_{*}({x_{t}}_{*}P\otimes P\right).$

By the definition of $K$, the above identities imply
$\tilde{\mathscr{K}}_{\operatorname{RKHS}}\mathcal{E}p_{t}=K\Phi(x_{t})$,
which completes the proof of the proposition.

### Proof of Proposition 8

Since identity $\phi(x)=\Phi(\delta_{x})$ holds for any $x\in\mathcal{X}$, and
by the definition of $K$, the following identities are derived:

$\displaystyle
K\phi(x)=K\Phi(\delta_{x})=\Phi({\beta_{t}}_{*}(\delta_{x}\otimes
P))=\int_{\omega\in\varOmega}\phi(\pi(t,\omega,x))\;dP(\omega).$

Let $g\in\operatorname{Span}\\{\phi(x)\mid x\in\mathcal{X}\\}$. Then $g$ is
represented as $g=\sum_{i=0}^{N-1}c_{i}\phi(\tilde{x}_{i})$ with some
$N\in\mathbb{N}$, $c_{0},\cdots,c_{N-1}\in\mathbb{C}$ and
$\tilde{x}_{0},\ldots,\tilde{x}_{N-1}\in\mathcal{X}$. Moreover, since
$\phi:\mathcal{X}\rightarrow\mathcal{H}_{k}$ is the feature map, the
reproducing property $\left\langle f,\phi(x)\right\rangle_{k}=f(x)$ holds for
any $f\in\mathcal{H}_{k}$. Therefore, the following identities hold:

$\displaystyle\left\langle\tilde{\mathscr{K}}f,g\right\rangle_{k}=\left\langle\tilde{\mathscr{K}}f,\sum_{i=0}^{N-1}c_{i}\phi(\tilde{x}_{i})\right\rangle_{k}=\sum_{i=0}^{N-1}\overline{c_{i}}\left\langle\int_{\omega\in\varOmega}f(\pi(t,\omega,\cdot))\;dP(\omega),\phi(\tilde{x}_{i})\right\rangle_{k}$
$\displaystyle\qquad=\sum_{i=0}^{N-1}\overline{c_{i}}\int_{\omega\in\varOmega}f(\pi(t,\omega,\tilde{x}_{i}))\;dP(\omega)=\sum_{i=0}^{N-1}\overline{c_{i}}\left\langle
f,\int_{\omega\in\varOmega}\phi(\pi(t,\omega,\tilde{x}_{i}))\;dP(\omega)\right\rangle_{k}$
$\displaystyle\qquad=\sum_{i=0}^{N-1}\overline{c_{i}}\left\langle
f,K\phi(\tilde{x}_{i})\right\rangle_{k}=\left\langle
f,K\sum_{i=0}^{N-1}c_{i}\phi(\tilde{x}_{i})\right\rangle_{k}=\left\langle
f,Kg\right\rangle_{k},$

which implies that $\tilde{\mathscr{K}}$ is the adjoint operator of
$K\mid_{\operatorname{Span}\\{\phi(x)\mid x\in\mathcal{X}\\}}$. This completes
the proof of the proposition.

### Proof of Proposition 15

The right-hand side of assumotion (4) is transformed into:

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}\int_{\omega\in\varOmega}f(h(\tilde{x}_{t+iS^{\prime}})+\xi_{t}(\omega))\;dP(\omega)$
$\displaystyle\quad=\int_{\omega\in\varOmega}\int_{x\in\mathcal{X}}f(h(x)+\xi_{t}(\omega))\;d{Y_{t}}_{*}P\;dP(\omega)$
$\displaystyle\quad=\int_{x\in\mathcal{X}}f(x)\;d{\left(h(Y_{t})+\xi_{t}\right)}_{*}P(x).$

Regarding the left-hand side of assumption (4), since
$h(\tilde{x}_{t+iS^{\prime}})+\xi_{t+iS^{\prime}}(\omega_{0})=\tilde{x}_{t+1+iS^{\prime}}$,
the assumption (4) is equivalent to assumption (12).

### Proof of Proposition 16

If $N=1$, then $\mu_{t,N}$ is represented as
$\mu_{t,N}=\delta_{\tilde{x}_{t}}$. Thus, identity
$\Phi(\mu_{t,N})=\phi(\tilde{x}_{t})$ holds. This implies that in this case,
Krylov subspace
$\mathcal{V}_{S}(K,\Phi(\mu_{0}))=\operatorname{Span}\\{\Phi(\mu_{0}),\ldots,\Phi(\mu_{S-1})\\}$
is equivalent to Krylov subspace (13).

### Proof of Proposition 20

Since $|||Q_{S}|||=|||Q_{S,N}|||=1$, the following inequalities hold:

$\displaystyle\|Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
$\displaystyle\qquad\leq\|Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})-Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
$\displaystyle\qquad\qquad+\|Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S}^{*}\phi(\tilde{x}_{t-1})-Q_{S,N}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
$\displaystyle\qquad\qquad+\|Q_{S,N}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
$\displaystyle\qquad\leq|||\tilde{\mathbf{K}}_{S,N}|||\|Q_{S,N}^{*}\phi(\tilde{x}_{t-1})-Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
$\displaystyle\qquad\qquad+\|\tilde{\mathbf{K}}_{S,N}Q_{S}^{*}\phi(\tilde{x}_{t-1})-\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
$\displaystyle\qquad\qquad+\|Q_{S,N}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
(19)

Since the elements of $Q_{S,N}^{*}\phi(\tilde{x}_{t-1})\in\mathbb{C}^{S}$ are
composed of the finite linear combinations of inner products between
$\Phi(\mu_{t,N})$ and $\phi(\tilde{x}_{t-1})$ in the RKHS, the same discussion
as $\mathbf{R}_{S,N}$ and $Q_{S,N}^{*}\Psi_{1,N}$ in Proposition 13 derives
$Q_{S,N}^{*}\phi(\tilde{x}_{t-1})\to Q_{S}^{*}\phi(\tilde{x}_{t-1})$ as
$N\to\infty$. Thus, the first term of Eq. 19 converges to $0$ as $N\to\infty$.
In addition, by Proposition 13,
$\tilde{\mathbf{K}}_{S,N}\to\tilde{\mathbf{K}}_{S}$ and $Q_{S,N}\to Q_{S}$
strongly in $\mathcal{H}_{k}$ as $N\to\infty$, which implies the second and
third terms of Eq. 19 also converge to $0$ as $N\to\infty$. Therefore,
$Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})$ converges to
$Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})$ as $N\to\infty$.
Since the norm $\|\cdot\|_{k}$ is continuous,
$\|Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})\|_{k}$ and
$\|\phi(\tilde{x}_{t-1})-Q_{S,N}\tilde{\mathbf{K}}_{S,N}Q_{S,N}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
converge to
$\|Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$ and
$\|\phi(\tilde{x}_{t-1})-Q_{S}\tilde{\mathbf{K}}_{S}Q_{S}^{*}\phi(\tilde{x}_{t-1})\|_{k}$
as $N\to\infty$, respectively. This implies $a_{t,S,N}\to a_{t,S}$ as
$N\to\infty$.

##  B Computation of QR decomposition of $\Psi_{1,N}$ and
$\tilde{\mathbf{K}}_{S,N}$

For implementing the Arnoldi method and shift-invert Arnoldi method described
in Section 4, QR decomposition must be computed. In this section, we explain
the method to compute the QR decomposition. The orthonormal basis of
$\operatorname{Span}\\{\Phi(\mu_{0,N}),\ldots,\Phi(\mu_{S-1,N})\\}$ for the
Arnoldi method or
$\operatorname{Span}\\{\sum_{j=0}^{i}\binom{i}{j}(-1)^{j}\gamma^{i-j}\Phi(\mu_{j})\mid
1\leq i\leq S\\}$ for the shift-invert Arnoldi method, which is denoted as
$q_{0,N},\ldots,q_{S-1,N}$, is obtained through QR decomposition. Then,
$Q_{S,N}$ is defined as the operator that maps
$[y_{0},\ldots,y_{S-1}]\in\mathbb{C}^{S}$ to
$y_{0}q_{0,N}+\ldots,y_{S-1}q_{S-1,N}$. The adjacent operator $Q_{S,N}^{*}$
maps $v\in\mathcal{H}_{k}$ to $[\left\langle
v,q_{0,N}\right\rangle_{k},\ldots,\left\langle
v,q_{S-1,N}\right\rangle_{k}]\in\mathbb{C}^{S}$.

First, the QR decomposition for the Arnoldi method is shown. For $t=0$,
$q_{0,N}$ is set as $q_{0,N}:=\Phi(\mu_{0,N})/\|\Phi(\mu_{0,N})\|_{k}$. For
$t>0$, $q_{0,N}$ is computed using $q_{0,N},\ldots,q_{t-1,N}$ as follows:

$\displaystyle\tilde{q}_{t,N}$
$\displaystyle:=\Phi(\mu_{t,N})-\sum_{i=0}^{t-1}\left\langle\Phi(\mu_{t,N}),q_{i,N}\right\rangle_{k}q_{i,N}$
(20) $\displaystyle q_{t,N}$
$\displaystyle:=\tilde{q}_{t,N}/\|\tilde{q}_{t,N}\|_{k}.$

Let the $(i,t)$-element of $\mathbf{R}_{S,N}$ be $r_{i,t}$, where $r_{i,t}$ is
set as $r_{i,t}:=\left\langle\Phi(\mu_{t,N}),q_{i,N}\right\rangle_{k}$ for
$i<t$, $r_{i,t}:=\|\tilde{q}_{t,N}\|_{k}$ for $i=t$, $r_{i,t}=0$ for $i>t$.
Then, by Eq. (20) and the definition of $\mathbf{R}_{S,N}$, $q_{i,N}$ is
represented as
$q_{i,N}=(\Phi(\mu_{i,N})-\sum_{j=0}^{i-1}r_{j,i}q_{j})/r_{i,i}$, and
$\Psi_{0,N}=Q_{S,N}\mathbf{R}_{S,N}$ holds. Therefore, by definition of
$r_{j,t}\ (j=0,\ldots,i-1)$, $r_{t,i}$ is computed as follows for $i<t$:

$\displaystyle r_{i,t}$
$\displaystyle=\left\langle\Phi(\mu_{t,N}),q_{i,N}\right\rangle_{k}=\left\langle\Phi(\mu_{t,N}),\frac{\Phi(\mu_{i,N})-\sum_{j=0}^{i-1}r_{j,i}q_{j}}{r_{i,i}}\right\rangle_{k}$
$\displaystyle=\frac{\left\langle\Phi(\mu_{t,N}),\Phi(\mu_{i,N})\right\rangle_{k}-\sum_{j=0}^{i-1}\overline{r_{j,i}}\left\langle\Phi(\mu_{t,N}),q_{j}\right\rangle_{k}}{r_{i,i}}$
$\displaystyle=\frac{\left\langle\Phi(\mu_{i,N}),\Phi(\mu_{t,N})\right\rangle_{k}-\sum_{j=0}^{i-1}\overline{r_{j,i}}r_{j,t}}{r_{i,i}}.$

Since $\mu_{t,N}=\sum_{j=0}^{N-1}\delta_{\tilde{x}_{t+jS}}$,
$\left\langle\Phi(\mu_{t,N}),\Phi(\mu_{i,N})\right\rangle_{k}$ is computed as
follows:

$\displaystyle\left\langle\Phi(\mu_{t,N}),\Phi(\mu_{i,N})\right\rangle_{k}=\left\langle\sum_{j=0}^{N-1}\delta_{\tilde{x}_{t+jS}},\sum_{j=0}^{N-1}\delta_{\tilde{x}_{i+jS}}\right\rangle_{k}$
$\displaystyle\qquad=\sum_{j,l=0}^{N-1}\left\langle\delta_{\tilde{x}_{t+jS}},\delta_{\tilde{x}_{i+lS}}\right\rangle_{k}=\sum_{j,l=0}^{N-1}\left\langle\phi(\tilde{x}_{t+jS}),\phi(\tilde{x}_{i+lS})\right\rangle_{k}=\sum_{j,l=0}^{N-1}k(\tilde{x}_{t+jS},\tilde{x}_{i+lS}).$

Similarly, by the definition of $r_{j,t}\ (j=0,\ldots,t-1)$, $r_{t,t}$ is
computed, since $\left\langle q_{i,N},q_{j,N}\right\rangle_{k}=1$ for $i=j$
and $\left\langle q_{i,N},q_{j,N}\right\rangle_{k}=0$ for $i\neq j$ as
follows:

$\displaystyle r_{t,t}^{2}$
$\displaystyle=\left\langle\tilde{q}_{t,N}\tilde{q}_{t,N}\right\rangle_{k}$
$\displaystyle=\left\langle\Phi(\mu_{t,N})-\sum_{j=0}^{t-1}\left\langle
q_{j,N},\Phi(\mu_{t,N})\right\rangle_{k}q_{j,N},\Phi(\mu_{t,N})-\sum_{j=0}^{t-1}\left\langle
q_{j,N},\Phi(\mu_{t,N})\right\rangle_{k}q_{j,N}\right\rangle_{k}$
$\displaystyle=\left\langle\Phi(\mu_{t,N})-\sum_{j=0}^{t-1}r_{j,t}q_{j,N},\Phi(\mu_{t,N})-\sum_{j=0}^{t-1}r_{j,t}q_{j,N}\right\rangle_{k}$
$\displaystyle=\left\langle\Phi(\mu_{t,N}),\Phi(\mu_{t,N})\right\rangle_{k}-2\Re{\left\langle\sum_{j=0}^{t-1}r_{j,t}q_{j,N},\Phi(\mu_{t,N})\right\rangle_{k}}+\left\langle\sum_{j=0}^{t-1}r_{j,t}q_{j,N},\sum_{j=0}^{t-1}r_{j,t}q_{j,N}\right\rangle_{k}$
$\displaystyle=\sum_{j,l=0}^{N-1}k(\tilde{x}_{t+jS},\tilde{x}_{t+lS})-2\sum_{j=0}^{t-1}r_{j,t}\overline{r_{j,t}}+\sum_{j=0}^{t-1}r_{j,t}\overline{r_{j,t}}$
$\displaystyle=\sum_{j,l=0}^{N-1}k(\tilde{x}_{t+jS},\tilde{x}_{t+lS})-\sum_{j=0}^{t-1}|r_{j,t}|^{2}.$

The above computations construct $\mathbf{R}_{S,N}$. Then, since the $(i,t)$
element of $Q_{S,N}^{*}\Psi_{1}$ is represented as
$\left\langle\Phi(\mu_{t+1,N}),q_{i,N}\right\rangle_{k}$,
$Q_{S,N}^{*}\Psi_{1,N}$ is computed in the same manner as $\mathbf{R}_{S,N}$.
The $\tilde{\mathbf{K}}_{S,N}$ is obtained by
$Q_{S,N}^{*}\Psi_{1}\mathbf{R}_{S,N}^{-1}$.

For the shift-invert Arnoldi method, by Eq. (8), the projection space is
represented as $\operatorname{Span}\\{w_{1,N},\ldots,w_{S,N}\\}$. Thus,
$\left\langle\Phi(\mu_{t,N}),\Phi(\mu_{i,N})\right\rangle_{k}$ is replaced
with $\left\langle
w_{t+1,N},w_{i+1,N}\right\rangle_{k}=\sum_{l=0}^{i+1}\sum_{j=0}^{t+1}\binom{i+1}{l}(-1)^{j+l}\binom{t+1}{j}\overline{\gamma}^{1+i-l}\gamma^{1+t-j}\left\langle\Phi(\mu_{l,N}),\Phi(\mu_{j,N})\right\rangle_{k}$.

##  C Pseudo-codes of Arnoldi and shift-invert Arnoldi methods

Let $\mathbf{R}_{S:T}$ be the matrix composed of $r_{i,t}\ (S\leq t\leq T,\
0\leq i\leq S-1)$. The pseudo-codes for computing $\tilde{\mathbf{K}}_{S}$
with the Arnoldi method and shift-invert Arnoldi method are shown in Algorithm
1 and 2, respectively.

Algorithm 1 Arnoldi method for Perron-Frobenius operator $K$ in an RKHS

1:$S,N\in\mathbb{N},\ \\{\tilde{x}_{0},\ldots,\tilde{x}_{NS-1}\\}$

2:$\tilde{\mathbf{K}}_{S,N}$

3:for $t=0,\ldots,S$ do

4: for $i=0,\ldots,S-1$ do

5: if $i<t$ then

6:
$r_{i,t}=(\sum_{j,l=0}^{N-1}k(\tilde{x}_{t+lS},\tilde{x}_{i+jS})-\sum_{j=0}^{i-1}\overline{r_{j,i}}r_{j,t})/r_{i,i}$

7: else if $i=t$ then

8:
$r_{t,t}=\sqrt{\sum_{j,l=0}^{N-1}k(\tilde{x}_{t+jS},\tilde{x}_{t+lS})-\sum_{j=0}^{t-1}|r_{j,t}|^{2}}$

9: else

10: $r_{i,t}=0$

11: end if

12: end for

13:end for

14:$\tilde{\mathbf{K}}_{S,N}=\mathbf{R}_{1:S}\mathbf{R}_{0:S-1}^{-1}$

Algorithm 2 Shift-invert Arnoldi method for Perron-Frobenius operator $K$ in
an RKHS

1:$S,N\in\mathbb{N},\ \gamma\notin\varLambda(K),\
\\{\tilde{x}_{0},\ldots,\tilde{x}_{NS-1}\\}$

2:$\tilde{\mathbf{K}}_{S,N}$

3:for $t=0,\ldots,S-1$ do

4: for $i=0,\ldots,t$ do

5: $g_{i,t}=\sum_{j,l=0}^{N-1}k(\tilde{x}_{i+jS},\tilde{x}_{t+lS})$

6: end for

7:end for

8:for $t=0,\ldots,S-1$ do

9: for $i=0,\ldots,S-1$ do

10: if $i<t$ then

11:
$r_{i,t}=(\sum_{l=0}^{i+1}\sum_{j=0}^{t+1}\binom{i+1}{l}\binom{t+1}{j}(-1)^{j+l}\overline{\gamma}^{1+i-l}\gamma^{1+t-j}g_{j,l}-\sum_{j=0}^{i-1}\overline{r_{j,i}}r_{j,t})/r_{i,i}$

12: else if $i=t$ then

13:
$r_{t,t}=\sqrt{\sum_{j,l=0}^{t+1}\binom{t+1}{j}\binom{t+1}{l}(-1)^{j+l}\gamma^{1+t-j}\overline{\gamma}^{1-t-l}g_{j,l}-\sum_{j=0}^{t-1}|r_{j,t}|^{2}}$

14: else

15: $r_{i,t}=0$

16: end if

17: end for

18:end for

19:for $i=0,\ldots,S-1$ do

20:
$\tilde{r}_{i}=(\sum_{l=0}^{i+1}\binom{i+1}{l}(-1)^{j+l}\gamma^{i+1-l}g_{0,l}-\sum_{j=0}^{i-1}r_{j,i}\tilde{r}_{j})/r_{i,i}$

21:end for

22:$\tilde{\mathbf{K}}_{S,N}=[\tilde{r}\
\mathbf{R}_{0:S-2}]\mathbf{R}_{0:S-1}^{-1}$

## References

  * Budišić et al. (2012) M. Budišić, R. Mohr, and I. Mezić.  Applied Koopmanism.  _Chaos (Woodbury, N.Y.)_ , 22:047510, 2012.
  * Črnjarić-Žic et al. (2017) N. Črnjarić-Žic, S. Maćešić, and I. Mezić.  Koopman operator spectrum for random dynamical systems.  _arXiv:1711.03146_ , 2017.
  * Crouzeix and Palencia (2017) M. Crouzeix and C. Palencia.  The numerical range is a $(1+\sqrt{2})$-spectral set.  _SIAM Journal on Matrix Analysis and Applications_ , 38(2):649–655, 2017.
  * Gallopoulos and Saad (1992) E. Gallopoulos and Y. Saad.  Efficient solution of parabolic equations by Krylov approximation methods.  _SIAM Journal on Scientific Statistics_ , 13(5):1236–1264, 1992.
  * Göckler (2014) T. Göckler.  _Rational Krylov Subspace Methods for ’ $\phi$’-functions in Exponential Integrators_.  PhD thesis, Karlsruher Instituts für Technologie, 2014.
  * Gretton et al. (2012) A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola.  A kernel two-sample test.  _Journal of Machine Learning Research_ , 13(1):723–773, 2012.
  * Grimm (2012) V. Grimm.  Resolvent Krylov subspace approximation to operator functions.  _BIT Numerical Mathematics_ , 52(3):639–659, 2012\.
  * Güttel (2010) S. Güttel.  _Rational Krylov Methods for Operator Functions_.  PhD thesis, Technischen Universität Bergakademie Freiberg, 2010\.
  * Hashimoto and Nodera (2018) Y. Hashimoto and T. Nodera.  Shift-invert rational Krylov method for an operator $\phi$-function of an unbounded linear operator.  _Japan Journal of Industrial and Applied Mathematics_ , 36(102):1–13, 2018.
  * Ikeda et al. (2019) M. Ikeda, I. Ishikawa, and Y. Sawano.  Composition operators on reproducing kernel hilbert spaces with analytic positive definite functions.  _arXiv:1911.11992_ , 2019.
  * Ishikawa et al. (2018) I. Ishikawa, K. Fujii, M. Ikeda, Y. Hashimoto, and Y. Kawahara.  Metric on nonlinear dynamical systems with Perron-Frobenius operators.  In _Advances in Neural Information Processing Systems 31_ , pages 2856–2866, 2018.
  * Kawahara (2016) Y. Kawahara.  Dynamic mode decomposition with reproducing kernels for Koopman spectral analysis.  In _Advances in Neural Information Processing Systems 29_ , pages 911–919, 2016.
  * Keogh et al. (2005) E. Keogh, J. Lin, and A. Fu.  Hot sax: efficiently finding the most unusual time series subsequence.  In _Fifth IEEE International Conference on Data Mining_ , 2005.
  * Klus et al. (2017) S. Klus, I. Schuster, and K. Muandet.  Eigendecompositions of transfer operators in reproducing kernel Hilbert spaces.  _arXiv:1712.01572_ , 2017.
  * Koopman (1931) B. O. Koopman.  Hamiltonian systems and transformation in Hilbert space.  _Proceedings of the National Academy of Sciences_ , 17(5):315–318, 1931.
  * Krylov (1931) A. N. Krylov.  On the numerical solution of the equation by which, in technical matters, frequencies of small oscillations of material systems are determined.  _Izvestija Akad. Nauk SSSR_ , 7(4):491–539, 1931\.  (in Russian).
  * Kubrusly (2012) C. S. Kubrusly.  _Spectral Theory of Operators on Hilbert Spaces_.  Birkhäuser$/$Springer, Boston, 2012.
  * Kutz (2013) J. N. Kutz.  _Data-Driven Modeling & Scientific Computation: Methods for Complex Systems & Big Data_.  OUP Oxford, 2013.
  * Lusch et al. (2017) B. Lusch, J. Nathan Kutz, and S. Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  _Nature Communications_ , 9:4950, 2017.
  * Malhotra et al. (2015) P. Malhotra, L. Vig, G. Shroff, and P. Agarwal.  Long short term memory networks for anomaly detection in time series.  In _European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning_ , pages 89–94, 2015.
  * McIntosh (1978) A. McIntosh.  The toeplitz-hausdorff theorem and ellipticity conditions.  _The American Mathematical Monthly_ , 85(6):475–477, 1978.
  * Moret and Novati (2004) I. Moret and P. Novati.  RD-rational approximations of the matrix exponential.  _BIT Numerical Mathematics_ , 44(3):595–615, 2004\.
  * Muandet et al. (2017) K. Muandet, K. Fukumizu, B. Sriperumbudur, and B. Schölkopf.  Kernel mean embedding of distributions: A review and beyond.  _Foundations and Trends in Machine Learning_ , 10(1–2):1–141, 2017.
  * R Hestenes and Stiefel (1952) M. R Hestenes and E. Stiefel.  Methods of conjugate gradients for solving linear systems.  _Journal of Research of the National Bureau of Standards_ , 49:409–436, 1952.
  * Saad and Schultz (1983) Y. Saad and M. H. Schultz.  GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems.  _SIAM Journal on Scientific and Statistical Computing_ , 7(3):856–869, 1983.
  * Sriperumbudur et al. (2011) B. K. Sriperumbudur, K. Fukumizu, and G. R. G. Lanckriet.  Universality, characteristic kernels and rkhs embedding of measures.  _Journal of Machine Learning Research_ , 12:2389–2410, 2011\.
  * Takeishi et al. (2017a) N. Takeishi, Y. Kawahara, and T. Yairi.  Subspace dynamic mode decomposition for stochastic Koopman analysis.  _Physical Review E_ , 96, 2017a.
  * Takeishi et al. (2017b) N. Takeishi, Y. Kawahara, and T. Yairi.  Learning Koopman invariant subspaces for dynamic mode decomposition.  In _Advances in Neural Information Processing Systems 30_ , pages 1131–1141, 2017b.
  * Takeuchi and Yamanishi (2006) J. Takeuchi and K. Yamanishi.  A unifying framework for detecting outliers and change points from time series.  _IEEE Transactions on Knowledge and Data Engineering_ , 18(4):482–492, 2006.  doi: 10.1109/TKDE.2006.1599387.