In this paper, we study how the Koopman operator framework can be combined with kernel methods to effectively control nonlinear dynamical systems. While kernel methods have typically large computational requirements, we show how random subspaces (Nyström approximation) can be used to achieve  huge computational savings while preserving accuracy.
Our main technical contribution  is deriving theoretical guarantees on the effect of the Nyström approximation.  More precisely, we study the linear quadratic regulator problem, showing that both the approximated Riccati operator and the regulator objective, for the associated solution of the optimal control problem, converge at the rate $m^{-1/2}$, where $m$ is the random subspace size. Theoretical findings are complemented by numerical experiments corroborating our results.

Koopman operator, kernel methods, Nyström method, linear quadratic regulator, data-driven methods

[auto, node distance=2cm,>=latex']

[chart_node_algorithm] (Simulate) Nonlinear dynamical system $\mathbf x_{k+1} = f(\mathbf x_k, \mathbf u_k)$.;

[chart_node_algorithm, right= of Simulate,xshift=-1.5cm] (Lift) Lift system to RKHS, $\infty$-dim. system $z_{k+1} = \G \begin{bmatrix}z_k\ \mathbf u\end{bmatrix}^T$.;

[chart_node_algorithm, right= of Lift,xshift=-1.5cm] (Sketch) Finite-dim. approximation with Nyström method, $\tilde z_{k+1}
=  \nG \begin{bmatrix}\tilde z_k\  \mathbf u\end{bmatrix}^T$.;

[chart_node_lqr, below=of Sketch, yshift=1.2cm] (Identify) System identification: Learn $\nG$ with sketched kernel ridge regression.;

[chart_node_lqr, left= of Identify, xshift=1.5cm] (LQR) Control: Solve the LQR problem with $\nG$.;

[chart_node_theory, below = of Identify,yshift=1.2cm] (Convergence_dynamics) Fix training set. For any regularization $\gamma$, $\lVert \G - \nG\rVert \lesssim \mathcal O(1/\sqrt m)$.;

[chart_node_theory, below= of LQR,yshift=1.2cm] (Convergence_LQR) The optimal control with $\nG$ is close to the one we would obtain with $\G$.;

(A) [left= of Simulate]  ;
(B) [below=of A.west, anchor=west]  ;
(C) [below = of B.west, anchor=west]  ;

[fit= (A) (Simulate) (Lift) (Sketch), rounded corners, draw, dashed, inner sep=0.05cm,align=left] (modelling)Modelling;
[fit= (B) (Identify) (LQR), rounded corners, draw, dashed, inner sep=0.05cm, align=left] (learning)Learning & control;
[fit= (C) (Convergence_dynamics) (Convergence_LQR), rounded corners, draw, dashed, inner sep=0.05cm, align=left] (theory)Guarantees;

[->] (Simulate) – (Lift);
[->] (Lift) – (Sketch);
[->] (Sketch) – (Identify);
[->] (Identify) – (LQR);
[->] (LQR) – ++ (-4,0) -| node [pos=0.99]  (Simulate) node[midway, xshift=2.9cm, yshift=-0.2cm] (m) Optimal gain for state feedback;

[->] (Identify) – (Convergence_dynamics);
[->] (LQR) – (Convergence_LQR);
[->] (Convergence_dynamics) – (Convergence_LQR);
[dashed, ->, color=Crimson] (Lift) – (LQR) node[midway] (m) Intractable;

A summary of the paper. Given some controls and corresponding state trajectories of a nonlinear dynamical system, we can use kernels to transform the dynamics and build a linear, data-driven model of the system. Kernels in general yield an infinite-dimensional representation of the state space, which we render finite-dimensional by means of the Nyström method. In a fixed design analysis, we prove that, as the number of Nyström landmarks $m$ grows, this approximation converges to the model with exact kernel.  We further show how this fact affects the LQR. Namely, the optimal control sequence computed with the Nyström-based dynamics converges to the optimal control sequence we would obtain when using an exact kernel, i.e., an infinite-dimensional transformation of the nonlinear state. The Riccati gain obtained in this way can therefore be used for a state feedback control scheme considering the true nonlinear system.

Summary of the main notations used in our article.

Variable     Meaning

$\cH_1$     RKHS associated to $k$

$\cH=\cH₁×\bR^{n_u}$     product space for lifted state and control

$\xvec \in\bR^d$,
$\yvec \in\bR^{n_y}$,
$\uvec \in\bR^{n_u}$
state, output, control input

$\uveco_i \in\bR^{n_u},\tuveco_i \in\bR^{n_u},\ i\in\mathbb N$
control sequences for the optimal gains, defined in (<ref>), (<ref>)

$\ldmi_i \in\bR^{d}, \ldmo_i \in\bR^{d},\ i=1, \dots, m$
Nyström input and output landmarks

$z \in\Hspace_1$     lifted state (infinite dimensional), see (<ref>)

$\tilde z \in\Hspace_1$     Nyström lifted state (infinite dimensional), see (<ref>)

$\tildezvec \in\bR^{ m}$     lifted state (finite dimensional, Nyström features), see (<ref>)

$\wvec=[\xvec^T,\uvec^T]^T\in \bR^{d+n_y}$     augmented states (concatenation of state and control input)

$\psi: \bR^d \rightarrow \cH_1$     transformation for the state

$\phi : \bR^{d+n_u} \rightarrow \cH_1 \times \bR^{n_u} $     transformation for state and control input, see (<ref>)

$S:\Hspace\rightarrow\R^n$     sampling on the input with normalization $n^{-1/2}$, see (<ref>)

$\Zx:\cH_1\rightarrow\bR^n$     sampling on the output with normalization $n^{-1/2}$, w/o the control input, see (<ref>)

$\Sx:\cH_1\rightarrow\bR^n$     sampling on the input with normalization $n^{-1/2}$, w/o the control input, see (<ref>)

$\Su:\bR^{n_u}\rightarrow\bR^n$     sampling on the control input with normalization $n^{-1/2}$, see (<ref>)

$\Pin:\Hspace\rightarrow\Hspace$    block diagonal projection for input space (state and control input)

$\Pinx:\Hspace_1\rightarrow\Hspace_1$   orthogonal projection on input state

$\Pout:\Hspace_1\rightarrow\Hspace_1$   orthogonal projection on output state, can be equal to $\Pinx$

$\Stildeout:\Hspace_1\rightarrow \R^m$   sampling operator for Nyström centers (regression output), see (<ref>)

$\G=[\A,\B]:\cH → \cH₁$     estimated (state,input)-state transition operator

$\A :\Hspace_1\rightarrow\Hspace_1$    estimated state-state transition operator, also appearing as $A$

$\B :\R^{n_u}\rightarrow\Hspace_1$    estimated input-state transition operator, also appearing as $B$

$\nG=[\nA, \nB]:\cH → \cH₁$     Nyström (state,input)-state transition operator

$ \nA :\Hspace_1\rightarrow\Hspace_1$    Nyström approximation of state-state transition operator, also appearing as $\Atilde$

$ \nB :\R^{n_u}\rightarrow\Hspace_1$    Nyström approximation input-state transition operator, also appearing as $\Btilde$

$ P :\cH_1\rightarrow\Hspace_1, \Ptilde :\cH_1\rightarrow\Hspace_1$    Riccati operator with exact and Nyströmized kernel, see (<ref>) and (<ref>)

$\K :\cH_1\to\R^{n_u}, \Ktilde :\cH_1\to\R^{n_u}$    LQR optimal gain with exact and Nyströmized kernel, see (<ref>) and (<ref>)

§ INTRODUCTION

Nonlinear dynamical systems are ubiquitous, and pose a great challenge  in terms of system identification and control.

A powerful approach to deal with nonlinear dynamical systems is provided by the Koopman operator framework [27, 10, 34, 11].
In this approach, the nonlinear dynamical system is transformed through a set of nonlinear functions, called observables, so that the dynamics of the transformed states are linear,
and can be used to reconstruct the states of the system.

The Koopman operator, at the basis of this technique, was originally introduced in the context of autonomous dynamical systems [27]. However, as shown in the seminal work of Korda et. al [28], the Koopman operator formulation allows to apply linear control techniques, which are well understood and efficient to compute.
As discussed by [11, 36], the control input can be included in the Koopman framework by either defining a family of Koopman operators (one for each value of the control input) [37], or by directly extending the state space to include the control input as an additional state of the system.  The latter perspective is the one we consider in this paper.

A key challenge  to apply the Koopman operator approach is to choose the suitable space of observable
functions, both with and without  considering the control input  [28].
Common choices include splines [28], polynomial or Fourier bases [1], and neural networks [41, 46, 22], see also  [29, 23, 19].

In this paper, we consider observables in a reproducing kernel Hilbert space (RKHS)  [5].  Band limited functions, splines and Sobolev spaces are then special cases [8].  This choice was already mentioned in [28] and has been recently  analyzed in
[30, 13, 25, 18, 9]
for the  forecasting and analysis of dynamical systems.
Kernel methods are popular in machine learning [40], as an RKHS has the advantage of being a possibly infinite dimensional space and correspond to universal approximators [43], while the associated estimators are nonparametric [44] and their computation reduces to  finite-dimensional numerical problems [39].
When efficiency is needed, further approximations are however required.  In this paper,  we consider the so-called Nyström method, based on considering a random subspace of functions defined by the data [45]. The effect of this approximation has been characterized for supervised machine learning, see, e.g., [38] and references therein, and more recently  for dynamical system identification, see [15, 33].

In this paper, we combine the above ideas and develop an efficient and accurate  control approach for nonlinear dynamical systems, based on using the Koopman operator framework together with kernels and the Nyström method. Our main technical contribution is the analysis of the approximation due to the Nyström method but,  unlike previous works,  we consider a control setting.
We focus on the linear quadratic regulator (LQR) method, which has appealing theoretical properties, such as an analytic form of the optimal solution, and easily allows to deal with multi-input-multi-output systems.   Our main result provides end-to-end theoretical guarantees for the proposed data-driven control framework.
Namely, we quantify how the error due to the Nyström approximation w.r.t. an exact kernel affects the solution of the LQR optimal control problem. Instrumental to our analysis are recent results on  learning and controlling linear systems [42, 32, 14].
We complement our theoretical analysis with some numerical experiments.
The combination of the Koopman approach and the LQR [35]  has been applied to challenging robot learning problems [2, 46], also involving soft robots [21]. Here, we assess the proposed control pipeline on a classic control benchmark, i.e., the Duffing oscillator [29], and on the compelling problem of identifying the dynamics of cloth [12, 31, 4, 47].

To conclude, we observe that concurrent to our work, Driessen [16] also considers kernel-based methods and their Nyström approximation for control purposes, but without any theoretical guarantees, and learning a bilinear [28, 35] data-driven system.

*Contributions As summarized in Fig. <ref>, our main contributions are the following:

* we show how the Koopman operator framework for controlled dynamical systems can be used to design linear control by lifting the state-space representation to an RKHS;

* we show how the Nyström method can be exploited to derive efficient computations;

* we prove finite sample bounds for the convergence of the Nyström data-driven system to the infinite dimensional state-space representation in RKHS;

* when using an LQR, we
show how these error rates due to the Nyström approximation translate to the associated Riccati operator, and to the optimal control sequence computed by solving the LQR;

* we provide a publicly available and open-source implementation of the proposed Nyström-based system identification and optimal control algorithms[Available at <https://github.com/LCSL/nys-koop-lqr>.=-1] that we test on some illustrative examples.

*Outline The structure of this paper is as follows.
In Section <ref>, we report the main technical background of our work, and we detail the notations used throughout the paper. In Section <ref>, we discuss how reproducing kernels and the Nyström method can be combined to obtain linear predictors for nonlinear systems, while in Section <ref> we discuss how such a prediction can be used in an LQR. In Section <ref>, we present the core theoretical contributions of our work, while Section <ref> contains a numerical evaluation of our identification algorithm and the associated LQR problem, on a standard control benchmark and on the challenging task of cloth manipulation.

§ BACKGROUND AND NOTATION

In this paper, we consider
a nonlinear controlled discrete dynamical system,
which is approximated for control purposes by a surrogate dynamical system whose dynamics is linear both in the lifted state and in the input variable;
this setup has been previously considered for instance by Korda et. al [28].
We denote the state vector as $\mathbf x\in \mathbb R^d$, the output as $\mathbf y\in \mathbb R^{n_y}$, and the control input vector as $\mathbf u\in\mathbb R^{n_u}$.
For a sequence of control inputs $(\uvec_t)_{t\in\mathbb N}$, an initial state $\xvec_0$, a function $f:\R^{d + n_u}\to\R^d$, and a matrix $E:\R^{d}\to\R^{n_y}$, we are interested in studying
deterministic nonlinear dynamical systems described by the difference equation

\begin{align}
\mathbf x_{t + 1} &= f(\mathbf x_t, \mathbf u_t)\label{eq:controlled_dynsys},\\
\mathbf y_t &= E\mathbf x_t,\label{eq:dynsys_output}
\end{align}

where $t \in\mathbb N$ denotes the discrete-time instant.
The state and the input can be merged in an augmented finite-dimensional state, denoted by $\wvec = [\mathbf x^T, \mathbf u^T]^T\in \R^{d + n_u}$. We will use both notations interchangeably throughout the article.
For a suitable space $\cA$ of observable functions $\xi$, the Koopman operator $\cK:\cA\to\cA$ can be defined as

\begin{align}
(\mathcal K\xi)(\wvec_t) %\\
\quad \forall \xi∊\cA. \label{eq:koopman_op_with_time}
\end{align}

As discussed by [30], if an invariant measure $\pi$ exists for the Markov transition kernel associated to the augmented dynamical system in $\wvec$,
then $\cA$ can be chosen to be $L_{\pi}^2$, i.e., the space of square integrable functions w.r.t. measure $\pi$. This fact is true also for deterministic systems, as the ones we consider in this paper (see <cit.>).

As we will discuss, a fundamental component of our learning framework is the notion of RKHS [5]. For an input space $\Xcal$, an RKHS $H$ is a Hilbert space of functions for which there exists a $k:\Xcal \times \Xcal \rightarrow \R$, the reproducing kernel, so that, for $\chi \in \Xcal$, and for $f \in H$, $k(\chi, \cdot)\in H$, and $f(\chi) = \innerprod{f, k(\chi, \cdot)}_{H}$. The latter notation denotes the inner product in the RKHS. An RKHS is a potentially infinite-dimensional space with universality properties, as it can be dense in $L_{\pi}^2$. Some examples of reproducing kernels are given by the Gaussian kernel, or by the Matérn family of kernels, which can be linked to Sobolev spaces. As we will show in the sequel, in this work, we will take as space of observables $\cA=L_{\pi}^2$ to ensure that the Koopman operator is well defined, and learn its restriction to an RKHS following [30, 33].

Note that the Koopman operator is linear in the observable $ξ$. It can thus be seen as a convenient way to build a linear approximation of the original dynamical system,
at the cost of manipulating the lifted state $ξ(\mathbf w)$ rather than the state $\mathbf x$.

In the following, we use the notation $\n{·}$ for the operator norm. We denote the Euclidean norm of a vector as $\norm{\cdot}_2$.
For any Hilbert space $\mathbb H$ we denote $\HS(\mathbb H)$ the Hilbert space of Hilbert-Schmidt operators on $\mathbb H$ and $\nHS{·}$ the associated Hilbert-Schmidt norm.
We denote $A^*$ the adjoint of an operator $A$, and $A^\dagger$ the pseudo-inverse of $A$. Besides, $\sigma_{\textrm{min}}(A)$ denotes the smallest singular value of $A$. A comprehensive summary of the notations used in the paper can be found in Table <ref>.

§ KOOPMAN SYSTEM IDENTIFICATION

As discussed previously, we focus on approximations of the original dynamical system which are linear
in the control input.
In the following, we will show how to achieve such a linearity in the control input. We will also detail the regression problems that are solved in order to retrieve the data-driven, kernel-based dynamics, and show how such kernel-based dynamics can be transformed to vector-valued dynamics by leveraging the Nyström method.

§.§ Choosing the Koopman lifting function

In order to define a suitable state and control transformation, we consider the RKHS $\cH₁$ associated to a stationary positive definite kernel $k:\R^d\times\R^d\to\R$. We let $\psi: \R^d\rightarrow \cH_1$, $\psi(x)\de k(x,·)$.	We define $\cH \coloneqq \cH_1\times \R^{n_u}$
and choose as our state and input transformation

\begin{align}
\phi:\R^{d+n_u}\rightarrow\cH,\
\phi(\wvec) &\coloneqq \vvec{\psi(\xvec) \\ \uvec}.
\label{e:def_phi}
\end{align}

The choice of the stationary kernel allows to lift the data through an infinite-dimensional nonlinear transformation. Moreover, $\phi$ is linear in the control input variable, which allows to use linear control techniques on the system of interest. However, $\phi$ is infinite-dimensional for many standard choices of kernel functions.
We will show in Section <ref>
that a related finite-dimensional lifting function can be obtained by
using a Nyström approximation of the kernel $k$.

§.§ Regression problem and corresponding solutions

We now explain how the function $\phi$ introduced in (<ref>)
can be leveraged
to obtain a linear, data-driven surrogate dynamical system to be used in place of the original nonlinear one when designing the control law. The difference equation of this data-driven model can be estimated by least-squares regression, as we detail in the following.

Having a dataset of $n$ training pairs $((\wvec_i, \wvec_{i+1}))_{i=1,\dots,n}$ including control inputs and corresponding states, we define the sampling operators for the system's state ($\Sx$) and control input ($\Su$) as

\begin{align}
S_{\xvec} &:\Hspace_1 \rightarrow\R^n,\
S_{\xvec} l= \tfrac{1}{\sqrt{n}} [l(\xvec_1), \dots, l(\xvec_n)]^T,
\label{e:def_Sx} \\
S_{\mathbf u} &:\R^{n_u} \rightarrow\R^n,\
S_{\uvec} \mathbf u= \tfrac{1}{\sqrt{n}}[\langle \mathbf u_1, \mathbf u\rangle, \dots,\langle \mathbf u_n, \mathbf u\rangle]^T.
\label{e:def_Su}
\end{align}

Moreover, as discussed by [28], when considering Koopman operator regression for controlled systems, the regression output can be limited to be the one-step-ahead state, i.e., we are not interested in forecasting the evolution of the control input variable. Thus, the sampling operator for the output can be written as

\begin{equation}
\Zx  :\Hspace_1 \rightarrow\R^n,\
\Zx g = n^{-1/2}[g(\xvec_{2}), \dots, g(\xvec_{n+1})]^T.
\label{e:def_Zx}
\end{equation}

Note that with these definitions, we can define the following compound sampling operator:

\begin{equation}
S\phi(\wvec) = S_{\xvec} \psi(\xvec)+ S_{\uvec} \uvec.\label{e:def_S}
\end{equation}

*Empirical risk minimization
The dynamics can be estimated by solving the following regularized least-squares problem in $\cH_1$:

\begin{equation}
\label{eq:koopman_regression_full}
\G	\de \arg\min_{W:\mathcal H\rightarrow\mathcal H_1}\mathcal R(W) + \gamma\lVert W\rVert^2_{\HS}
\end{equation}

\begin{equation}
\cR(W) \de \frac1n \sum_{i=1}^{n}\n*{ \psi(\xvec_{i+1})-W ϕ(\wvec_i)}_{\cH_1}²
\label{e:def_risk}.
\end{equation}

Here $γ>0$ is a regularization parameter, and we recall that $ϕ$ is defined in (<ref>). Note that this regression problem is not strictly speaking a Koopman regression problem, as the regression input and output spaces are different.
The objective in (<ref>) is strictly convex and coercive, and thus admits a unique solution.
As shown in Appendix <ref>, the risk function can be rewritten as
$\cR(W)= \nHS{\Zx - SW^*}²$ and thus the solution of (<ref>) can be expressed as

\begin{align}
\G &= \Zx^*(SS^* + \gamma I)^{-1}S.
\label{e:def_G}
\end{align}

For any initial state $\xvec_0∊\bR^d$, this operator defines the following linear dynamics:

\begin{align}[left = \empheqlbrace\,]
z_0 &= \psi(\xvec_0),
\label{e:initial_state}\\
z_{t+1} &=  \G \vvec{z_t \\ \uvec_t}.
\label{eq:dynamics_compound}
\end{align}

*Dynamics are affine in the control input
Note that the operator $\G:\cH→\cH₁$ defined in (<ref>) can be decomposed in two operators $\A:\cH₁ → \cH₁$ and $\B:\R^{n_u}→\cH₁$, controlling respectively the parts of the dynamics due to the state and to the control input.
More precisely, defining

\begin{align}
\A &= \Zx^*(SS^* + \gamma I)^{-1} S_{\xvec},
\label{e:def_A}\\
\B &= \Zx^*(SS^* + \gamma I)^{-1} S_{\uvec},
\label{e:def_B}
\end{align}

the dynamics (<ref>) are equivalent to the autoregressive linear model

\begin{align}
z_{t+1} &= \A z_t+ \B \uvec_t
\label{eq:dynamics}.
\end{align}

§.§ Nyström approximation

Given that the lifted state $z$ is typically infinite-dimensional,
we are now interested in designing a finite-dimensional approximation of the dynamics in (<ref>), which would be more useful for practical control purposes.

In this section, we thus approximate the nonlinear kernel $\kr$ using a Nyström approximation [45].

The approximation is based on the choice of two sets $\ldmi{1}, …,\ldmi{m}$ and $\ldmo{1},…,\ldmo{m}$ of $m$ points, called the input and output landmarks.
Multiple approaches have been studied for landmark selection, but in this paper we consider the simple setting where the landmarks are either all drawn uniformly from the dataset, or the input landmarks are drawn uniformly from the dataset and the output landmarks are then taken one step ahead in time the input ones.

Now, let us define
$\Hin \de \spa\{\psi(\ldmi{1}), …, \psi(\ldmi{m})\}$,
$\Hout\de\spa\{\psi(\ldmo{1}), \dots, \psi(\ldmo{m})\}$,
and let $\Pinx:\Hspace_1→\Hspace_1, \Pout:\Hspace_1\rightarrow\Hspace_1$ be the orthogonal projectors, respectively onto $\Hin$ and $\Hout$.
Moreover, we define $\Pin:\Hspace\rightarrow\Hspace$
as $\Pin \phi(\wvec) = [\Pinx \psi(\xvec) , \uvec]$; as the control variable is already finite-dimensional, we indeed only need to project the lifted state in order to obtain a finite-dimensional approximation.

Following [33], we define the Nyström approximation of $\G$ as follows:

\begin{align}
%[\tilde A,\tilde B]
\nG
%&\de \argmin_{W:\cH\to\cH_1}\frac1n \sum_{i=1}^{n}\n*{ \Pout\psi(\xvec_{i+1})-W \Pin ϕ(\wvec_i)}_{\cH_1}² \\
%&\quad+ \gamma\lVert W\rVert^2_{\HS} \\
&\de \arg\min_{W:\cH→\cH_1} \cR(\Pout W\Pin) + γ\nHS{W} \nonumber\\
&= \Pout\Zx^*S\Pin (\Pin S^*S\Pin + \gamma I)^{-1}\nonumber \\
%&= \Pout\Zx^* (S\Pin\Pin S^* + \gamma I)^{-1} S\Pin\\
&=\Pout\Zx^* (S\Pin S^* + \gamma I)^{-1} S\Pin.
\label{eq:a_b_compound_expr_nystr}
\end{align}

The operator $\nG$ can be used to define linear dynamics approximating the ones from (<ref>),
namely for any initial condition $\xvec_0∊\bR^d$:

\begin{align}[left = \empheqlbrace\,]
\tilde z_0
&= \Pout\psi(\xvec_0),\label{e:n_initial_state}\\
\tilde z_{t+1}
&= \nG \vvec{\tilde z_t\\ \uvec_t}.
\label{eq:nystrom_dynamics_compound}
\end{align}

The projection in the initial condition guarantees that all the states visited during the evolution of the system belong to $\Hout$.

The dynamics in (<ref>) are still linear in the control input, and the operator $\nG$ could be decomposed into operators $(\nA,\nB)$, approximating the operators $(\A,\B)$ defined in (<ref>) and (<ref>).
We choose however to manipulate only $\nG$ in the following to keep expressions more concise.

*Vector-valued state-space representation
The dynamics in (<ref>) are defined for evolving functions: despite being all restricted to a finite-dimensional subspace, the iterates $(\tilde z_t)_{t∊\bN}$ still belong to a functional space.
In order to retrieve a vector-valued state representation that is practically computable,
we will thus look at the dynamics of the coordinates of these evolving functions in a basis of the subspace to which they belong.

More precisely, we can define $\Stildeout$ as

\begin{align}
\Stildeout:\Hspace_1\rightarrow\R^m,\
\Stildeout g &= [g(\ldmo{1}), …,g(\ldmo{m})]ᵀ.
\label{e:def_Stout}
\end{align}

Note that $\Stildeout\Stildeout^*=\Kmout\in\R^{m\times m}$, where $\Kmout$ is the Gram matrix of the stationary kernel $\kr$ computed at the output Nyström centers $\tilde\xvec_i^{\textrm{out}}$.
Defining $\Bout=\Stildeout^* (\Kmout^†)^{1/2}$, it holds $\Bout\Bout^*=\Pout$.

Note that for any initial condition of the form (<ref>), the Nyström states $\tilde z_t$ defined by (<ref>) belong to $\Hout$ for all $t\in\mathbb N$.
defining the following finite-dimensional autoregressive dynamics

\begin{align}[left = \hspace{-0.5cm}\empheqlbrace]
\tildezvec_0
&\de \Bout^*\psi(\xvec_0),
\label{e_n_vec_initial_state}
\\
\tildezvec_{t+1}
&\!=\! \Bout^*\Zx^* (S\Pin S^*\! +\! \gamma I)^{-1} S\Pin\!\vvec{\Bout \tildezvec_t \\ \uvec_t}\!,
\label{eq:nystrom_dynamics_vector_partial}
\end{align}

one can check that

\begin{align}
\tilde z
&\de \Bout \tildezvec
\label{eq:parameterized_nystrom_dynamics}
\end{align}

is equivalent to (<ref>) and (<ref>).

As we detail in Appendix <ref>, the dynamics (<ref>) can be rewritten in terms of matrix products that can be computed efficiently, and only requires to invert an $m×m$ matrix.

Note that the dynamics of $\tilde z$ could similarly be expressed in any orthonormal basis of $\Hout$, however working with $\Bout$ naturally yields expressions where the kernel matrix $\Kmout$ appears, which is convenient for implementation.

*Reconstructing the original state The lifted state $\tildezvec$ can be used to reconstruct the original state $\xvec$, as discussed in <cit.>. This goal can be achieved, for instance, by using a least-squares estimate, or by augmenting the lifted state with the original one. In this work, we consider a regularized least squares estimate for the state reconstruction matrix $C$. For a given regularization parameter $λ > 0$, we define

\begin{equation}
C = \min_{M\in\R^{d\times m}}\sum_{i=1}^{n}\norm{\xvec_{i+1} - M\tildezvec_{i+1}}_2^2 + λ n\norm{M}_{\HS}^2.
\label{eq:reconstruction}
\end{equation}

A closed-form expression for matrix $C$ is also included in Appendix <ref>.

§ KERNELS AND KOOPMAN LQR

Once we have estimated the state-space representation of the dynamical system of interest, we can use linear predictive control techniques, as the Koopman approach transforms a non-linear system in a linear one [28].
In the following, we focus on the LQR,
and more particularly on the infinite horizon LQR [32]
due to its theoretically appealing properties.
The key idea is that the optimization problem at the basis of the LQR is rendered tractable by using a finite dimensional embedding of the state and inputs, which can be achieved by sketching techniques, as discussed in Section <ref>.

*LQR for the the exact dynamics

We first consider using the exact kernel $\kr$ and the transformation $ϕ$.

We consider a control objective that is quadratic in the lifted state and control inputs, via the weighting operators $Q: \mathcal H_1\rightarrow\mathcal H_1$ and $ R:\R^{n_u}\rightarrow\R^{n_u}$.
Then, the LQR strategy requires to solve the following optimal control problem, for which a time-invariant analytical solution is available:

\begin{align}%[left = \empheqlbrace\,]
\label{eq:lqr_koopman_infinite_dimensional}
\min_{\mathbf u_0, \uvec_1,\dots}\lim_{T\rightarrow \infty}&\sum_{i=0}^{T}\langle  z_i, Q z_i\rangle_{\mathcal H_1} +
\langle \mathbf u_i,  R\mathbf u_i\rangle_{\R^{n_u}}\\
\mathrm{s.t.\ }  &\eqref{e:initial_state}, \eqref{eq:dynamics_compound}.%z_{t+1} &= \G \vvec{z_t \\ \uvec_t},
\nonumber
\end{align}

*LQR for the approximated dynamics
When the dynamics are approximated by using the Nyström approach as in (<ref>), the optimal control problem becomes

\begin{align}%[left = \empheqlbrace\,]
\label{eq:lqr_koopman_nystrom_infinite_dimensional}
\min_{\mathbf u_0, \uvec_1,\dots}\lim_{T\rightarrow \infty}\sum_{i=0}^{T}&\langle  \tilde z_i, Q \tilde z_i\rangle_{\mathcal H_1} +
\langle \mathbf u_i,  R\mathbf u_i\rangle_{\R^{n_u}}\\
\nonumber
\mathrm{s.t.\ }  &\eqref{e:n_initial_state}, \eqref{eq:nystrom_dynamics_compound}.%\tilde z_{t+1}
\end{align}

According to (<ref>), the problem in (<ref>) can equivalently be rewritten in the basis $\Bout$ for weighting matrices $\tilde Q\in \R^{m\times m}$,

\begin{equation}
\tilde Q=\Stildeout^*(\Kmout^\dagger)^{1/2} Q(\Kmout^\dagger)^{1/2}\Stildeout,\label{eq:Q_definition}
\end{equation}
and $R$, as follows:

\begin{align}%[left = \empheqlbrace\,]
\label{eq:lqr_koopman_nystrom_finite_dimensional}
\min_{\mathbf u_0, \uvec_1,\dots}\lim_{T\rightarrow \infty}
\sum_{i=0}^{T}&\langle  \tildezvec_i, \tilde Q \tildezvec_i\rangle_{\mathcal \R^m}
+ \ip{ \mathbf u_i,  R\mathbf u_i }_{\R^{n_u}}\\
\nonumber
\mathrm{s.t.\ } &\eqref{e_n_vec_initial_state}, \eqref{eq:nystrom_dynamics_vector_partial}.%
\end{align}

One can typically choose $\tilde Q=C^*Q'C$ where $C$ is the reconstruction matrix defined in (<ref>) in order to define an objective that can be interpreted as a penalization of the states.

*Closed-loop operation
If the problems in (<ref>) and (<ref>) are equivalent, they yield the same optimal control sequence in form of state feedback [20]. Denoting the optimal gain resulting from (<ref>) as $\tilde K: \cH_1\rightarrow \R^{ n_u}$, the optimal control law is

\begin{equation}
\uvec_k = \tilde K\tilde z_k.
\end{equation}

From (<ref>), we have that the state-feedback input can be rewritten as

\begin{equation}
\uvec_k = \tilde K\Bout\tilde \zvec_k,
\end{equation}

with $ \tilde K\Bout:\R^m\rightarrow \R^{n_u}$ being the solution of the LQR problem in (<ref>). When dealing with the true nonlinear dynamics in Section <ref>, we consider a state-feedback control law of the form

\begin{equation}
\uvec_k = \tilde K\Pout\psi(\xvec_k),
\end{equation}

where $\xvec_k$ is the true state of the system, to perform control in closed loop.

*Overall pipeline The content of this section can be summarized in a whole pipeline that can be used for linear identification and control of nonlinear systems. Such a pipeline begins with a system identification phase, in which the state dynamics of the nonlinear system are transformed with the Nyström approach. The linear model identified in this way is then used to solve a linear optimal control problem, and the corresponding state-feedback gain is applied to the real system. A summary of the pipeline can be found in Algorithm <ref>.

Nyström-Koopman system identification & LQR-based closed-loop control strategy
[1]
Training samples, LQR weights $\Qtilde, R$.
Vector-valued, data-driven linear system, optimal state-feedback gain.
Sample Nyström input and output landmarks $\ldmi{1}, …,\ldmi{m}$ and $\ldmo{1},…,\ldmo{m}$ uniformly from the training set.
Compute the operators in (<ref>) to obtain a data-driven linear system.
Use these linear dynamics to solve problem (<ref>) and get optimal state-feedback gain.

§ THEORETICAL ANALYSIS

In this section, we perform a theoretical analysis of the proposed system identification method, and assess its effect on the LQR problem. In particular, we will compare the Nyström-based approach with the data-driven model based on the exact kernel. Overall, we show that the Nyström-based model of (<ref>) is a provably accurate approximation of the dynamics in (<ref>), that can be safely used for control purposes, in place of the intractable infinite-dimensional model of  (<ref>).

*Layout In Section <ref>, we state the assumptions at the basis of our analysis. In Section <ref>, we bound the error, introduced by the Nyström approximation, on operators $\A$ and $\B$ from (<ref>) (compound in the transition operator $\G$). In Section <ref>, we show that the Riccati operator obtained with the Nyström dynamics is close to the one obtained with the exact kernel. Finally, in Section <ref>, we use the aforementioned results to show that, when plugging the optimal control from (<ref>) into the dynamics of (<ref>), the LQR objective function is close to the one obtained when solving (<ref>) directly.

§.§ Hypotheses

In this section, we introduce the hypotheses of our theoretical derivations in the next sections.
[Bounded kernel]

The stationary kernel $k$ is bounded, i.e., $k(\xvec, \xvec')\leq \supk$.

Note that under Assumption <ref>, it holds in particular

\begin{align*}
\norm{\Zx} &\leq κ
\quad\text{and}\quad%\\
\norm{\Sx} \leq κ.
\end{align*}

Indeed, for any $x$, it implies $\n{ψ(x)}=\sqrt{k(x,x)}≤κ$, and thus for any $g∊\cH_1$ with $\n{g}≤1$, it holds

\begin{align*}
\n{\Zx g}_2
&= n^{-1/2}\n{[\ip{g, ψ(\xvec_{2})}, \dots, \ip{g, ψ(\xvec_{n+1})}]ᵀ}_2\\
&= n^{-1/2}\left(\sum_{i=2}^{n+1} |\ip{g, ψ(\xvec_i)}|² \right)^{1/2} \\
&≤ n^{-1/2}\left(\sum_{i=2}^{n+1} \n{g}²κ² \right)^{1/2} %\\
≤ κ,
\end{align*}

and a similar argument holds for $\Sx$.

The following assumptions are standard and guarantee that the LQR problems in (<ref>) and (<ref>) are well-posed and admit an analytical solution in the form of a static state-feedback gain [20]. Minimal assumptions on $\A$, $\B$ from (<ref>) and (<ref>), and $C$ from (<ref>), for stabilizability and detectability to hold, could be derived by leveraging the theoretical deployments reported, e.g., in [7].

The dynamical systems in (<ref>) and (<ref>) are stabilizable, i.e.,

\begin{align}
\exists M&:\Hspace_1\rightarrow \R^{n_u} \textrm{\ such that\ }{\rho(\A + \B M)<1},\nonumber\\
\exists \tilde M&:\Hspace_1\rightarrow \R^{n_u} \textrm{\ such that\ }{\rho(\nA + \nB\tilde M)<1},\nonumber
\end{align}

where $\rho(W)$ is the spectral radius of operator $W$.

Stabilizability assumes that the dynamical systems can be driven to zero with a suitable state-feedback gain. This is a weaker assumption than controllability (i.e., assuming that the system can be driven to any location via the choice of suitable feedback gain), that is often used in the analysis of the LQR algorithm [32]. However, the controllability assumption in $\cH_1$ is not satisfied by the Nyström dynamics, as the $\tilde z$ functions always live in the subspace $\Hout$, as we discussed in Section <ref>.

The dynamical systems in (<ref>) and (<ref>) are detectable, i.e.,

\begin{align}
&\exists t,s,b,d\geq 0 \textrm{\ such that\ }\norm{\A^tz}_{\cH_1}\geq b\norm{z}_{\cH_1}\nonumber\\
&\quad\Rightarrow \innerprod{z, \sum_{i=0}^s\A^{i*}Q\A^iz}_{\cH_1}\geq d\innerprod{z, z}_{\cH_1},\nonumber\\
&\exists t,s,b,d\geq 0 \textrm{\ such that\ }\norm{\nA^t\tilde z}_{\cH_1}\geq b\norm{\tilde z}_{\cH_1}\nonumber\\
&\quad\Rightarrow \innerprod{\tilde z, \sum_{i=0}^s\nA^{i*}Q\nA^i\tilde z}_{\cH_1}\geq d\innerprod{\tilde z, \tilde z}_{\cH_1}.\nonumber
\end{align}

Detectability assumes that, if unstable dynamics happen in the linear system, these must be observed (and taken into account in the design of the LQR objective function).

§.§ Accuracy of the Nyström approximation of the transition operator

We now upper-bound the error (in operator norm) induced by the Nyström approximation on the transition operator $\G$.
Although other works studied sketched estimators of the Koopman operator,
our bound notably differs from [3, 33] who consider different norms (Hilbert-Schmidt and operator norm but on different spaces) and dynamical systems without control.

Under Assumption <ref>, for any $γ>0$ , it holds with probability $1-δ$ that

\begin{align}
\n{\nG-\G}
&≤ \prt*{\frac{κ}{γ}+1} 4\supfmap\sqrt{\frac{3}{m}\log\prt*{\frac{8m}{5δ}}} \nonumber\\
&\quad	+ \frac{48\supfmap^3}{γ^{3/2}} \frac{1}{m}\log\prt*{\frac{8m}{5δ}}.
\label{e:bound_G_nG_opnorm}
\end{align}

The proof for this result is provided in Appendix <ref>.

Note that, if we rewrite the approximate dynamics $\tilde z_{t+1} = \nG \phi(\wvec_t)$ from (<ref>) in the autoregressive form as

\begin{align}
\tilde z_{t+1}
&= \nA ψ(\xvec_t)+ \nB \uvec_t.
\label{e:autoregressive_dynamics}
\end{align}

Theorem <ref> automatically translates in bounds on the approximation of $\A$ and $\B$. Indeed, under the assumptions of Theorem <ref>, the right-hand side of
(<ref>) is also an upper bound on $\n{\A-\nA}$ and $\n{\B-\nB}$.
Formally, the operators $\nA,\nB$ in (<ref>) can be defined from $\nG$ as follows:

\begin{align*}
\nA
&= \nG [I_{\cH_1}, 0_{\cH_1\rightarrow\bR^{n_u}}]^* : \cH_1\rightarrow \cH_1 \\
\nB
&= \nG [0_{\bR^{n_u}\rightarrow \cH_1}, I_{\bR^{n_u}}]^* : \bR^{n_u} \rightarrow \cH_1.
\end{align*}

As a consequence, it holds $\n{\A-\nA}=\n{(\G-\nG) [I_{\cH_1}, 0_{\cH_1\rightarrow\bR^{n_u}}]^*}≤\n{\G-\nG}$, and a similar argument holds for $B$.

In the following, we will use for simplicity the notations $\A\de \A_γ,\B\de \B_γ,\nA\de \nA_γ,\nB \de \nB_γ$. However, all these operators implicitly depend on the choice of the regularization parameter $γ$.

§.§ Convergence analysis for the Riccati operator

In this section, we show that the Riccati operators for problems in (<ref>) and (<ref>) are $\epsilon$-close in operator norm, provided that $\norm{\G - \nG}\leq \epsilon$. Note that, according to Theorem <ref>, we can set $\epsilon$ as function of $m$, namely

\begin{equation*}
\epsilon=\prt*{\frac{κ}{γ}+1} 4\supfmap\sqrt{\frac{3}{m}\log\prt*{\frac{8m}{5δ}}}
+ \frac{48\supfmap^3}{γ^{3/2}} \frac{1}{m}\log\prt*{\frac{8m}{5δ}}.
\end{equation*}

This fact means that we transfer the guarantees for $\G$ from Section <ref> to guarantees on the fundamental building block of the LQR solution, i.e., the Riccati operator.

\begin{align}
F(P, A, B) &= P - A^*[P - PB(R + B^* P B)^{-1}B^* P ]A\nonumber\\
&= P - A^*P(I + B R^{-1}B^*P)^{-1}A - Q.\nonumber
\end{align}

Then, as proven by <cit.>, under the assumptions in Section <ref>, the LQR problems in (<ref>) and (<ref>) admit analytical solutions given by the static state-feedback gains

\begin{align}
\K:\cH_1\to\ \R^{n_u},\ \K &= -(R+\BᵀP\B)^{-1}\B^* P\A
\label{e:def_K},\\
\nK:\cH_1\to\ \R^{n_u},\ \nK &= -(R+\nBᵀ \nP \nB)^{-1} \nB^* \nP \nA
\label{e:def_nK},
\end{align}

where $P,\Ptilde:\cH₁→\cH₁$ are the unique self-adjoint, positive semi-definite operators obtained by solving the following discrete-time algebraic Riccati equations  [20]:

\begin{align}
\label{e:def_P}\\
\label{e:def_nP}.
\end{align}

A fundamental step is now to derive an error rate for the solutions, upper bounding the quantity $\norm{P - \Ptilde}$. To do so, we can state and prove an analogous proposition to <cit.>, based on the fundamental results in [26], and generalized to the case of operator dynamics.

Let the systems in (<ref>) and (<ref>) be stabilizable and detectable. Let $\sigma_{\textit{min}}(P)$ be the smallest singular value of the Riccati operator associated to (<ref>), and let $\Ptilde$ be the Riccati operator for (<ref>). Assume $R$ is positive definite, and $\sigma_{\textrm{min}}(P) \geq 1$. Let $L=A + BK$, where $K$ is the stabilizing Riccati gain defined in (<ref>), let $\rho(L)$ be the spectral radius of $L$, and let

\begin{equation}
\tau(L, \rho)=\sup\{\norm {L^k}\rho^{-k}, k\geq 0\},\nonumber
\end{equation}

with $\rho$ being a constant such that $\rho(L) \leq\rho < 1$.
Then, for $\norm{A - \Atilde}\leq \epsilon$, $\norm{B - \Btilde}\leq \epsilon$, and

\begin{align*}
\epsilon< \min\left\{\norm B,\right.&\left.\frac{1}{12} \frac{1}{\normp{L}^{2} +\normp{P}}\frac{(1 - \rho^2)^2}{\tau(L,\rho)^4}\right.\\
&\cdot\normp A^{-2} \normp P^{-2} \\
&\cdot\left.\normp B^{-3}\cdot\normp {R^{-1}}^{-2}\right\},
\end{align*}

it holds that the error on the Riccati operator due to the Nyström approximation is upper bounded as follows:

\begin{align*}
\norm{P-\tilde P}&\leq 6\epsilon \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp A^2 \normp P^2 \\
&\quad\cdot\normp B \normp {R^{-1}}.
\end{align*}

When talking about stabilizability, we are considering the data-driven dynamics described by (<ref>) and (<ref>). This fact means that we assume these systems are stabilizable in $\cH_1$, not in the original state space $\R^{d}$.

The proof of this result is close to the one of <cit.> which covers the setting where $A$ and $B$ are matrices. Some technical adjustments are required due to the fact that we work with operators. We report a sketch of the proof in Appendix <ref>.

§.§ Convergence analysis of the LQR objective function

In the previous subsections, we showed that $\norm{A - \Atilde}\leq \epsilon$, $\norm{B- \Btilde}\leq \epsilon$ for a sufficiently large number of Nyström landmarks $m$. We further showed that this implies that $\norm{P - \Ptilde}\leq \mathcal O(\epsilon)$. Now, we are interested in assessing how suitable the Nyström dynamics are for solving the LQR problem, instead of using the intractable, infinite-dimensional dynamics associated to an exact kernel. In order to do so, we first solve the LQR problem from (<ref>) based on the Nyström approximation of the dynamics. Then, we plug the optimal control sequence retrieved in this way in the exact kernel dynamics of (<ref>), and compare the LQR objective function with this control sequence, and the control sequence we would obtained by solving the problem in (<ref>). This type of error analysis resembles the one in [32], in the context of certainty equivalence. In order to do so, let us define
two control sequences ($\tuveco_0,\tuveco_1, \dots)$, and ($ \uveco_0, \uveco_1, \dots)$.
Moreover, for a given initial condition $\hat z_0 = z_0$, let

\begin{align}
\cJ &\coloneqq \lim_{T\to\infty}\sum_{i=0}^{T}\langle  z_i,  Q z_i\rangle_{\mathcal \cH_1}
+ \ip{  \uveco_i,  R\uveco_i }_{\R^{n_u}},\nonumber \\
&\quad z_{i+1} = A z_i + B \uveco_i, \uveco_i=\K z_i
\label{eq:obj},\\
\cJhat &\coloneqq \lim_{T\to\infty}\sum_{i=0}^{T}\langle  \hat z_i,  Q \hat z_i\rangle_{\mathcal \cH_1}
+ \ip{  \tuveco_i,  R\tuveco_i }_{\R^{n_u}},\nonumber\\
&\quad \hat z_{i+1} = A \hat z_i + B \tuveco_i, \tuveco_i=\Ktilde \zhat_i
\label{eq:nobj}.
\end{align}

Then, we upper bound the error $\cJhat - \cJ$.

Let $z_0$ be the initial state of the dynamical system of interest. Let $\rho(L)$ be the spectral radius of the closed-loop operator $L=A + BK$, where K is given by (<ref>). Moreover, let $\cJ$ be as in (<ref>), and $\cJhat$ as in (<ref>). Let $\rho$ be a real number such that $\rho(L)\leq \rho < 1$. Lastly, let $\Gamma = 1 + \max\{\norm A, \norm B, \norm P, \norm K\}$.  Under Assumptions <ref>, <ref>, and <ref>, if $\norm{A - \Atilde}\leq \epsilon$, $\norm{B - \Btilde} \leq \epsilon$, $\norm{P - \Ptilde}\leq g(\epsilon)$, with $g(\epsilon)\leq \frac{1 - \rho}{6\norm{B}\tau(L, \rho)\Gamma^2}$, and $\sigma_{\textrm{min}}(R) \geq 1$, we have that

\begin{equation*}
\cJhat - \cJ\leq 36 \sigma_{\textrm{max}}(R)\Gamma^9g(\epsilon)^2\kappa^2\frac{\tau(L, \rho)^2}{1 - \rho^2}.
\end{equation*}

The proof of this result can be obtained by following the same steps of <cit.>. Although deployed for $\A$ and $\B$ being matrices, these steps do not depend on the dimensionality of the state space, i.e., can be used also when considering function-valued dynamics, as we do in our work.

§ SIMULATION RESULTS

§ CONCLUSIONS

In this paper, we have studied the combination of the Koopman operator framework with reproducing kernels and the Nyström method, to design linear control for nonlinear dynamical systems. In a novel theoretical analysis, we have provided end-to-end guarantees, linking novel error rates, due to the Nyström approximation, on the dynamical system's representation with its effect on the linear quadratic regulator problem. Lastly, we have evaluated the effectiveness of the proposed method both on a classic control benchmark, and on the challeging problem of dynamic cloth manipulation.

§ ACKNOWLEDGMENTS

E. Caldarelli, A. Colomé and C. Torras acknowledge support from the project CLOTHILDE (“CLOTH manIpulation Learning from DEmonstrations”), funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Advanced Grant agreement No 741930).
E. Caldarelli acknowledges travel support from ELISE (GA No 951847).
L. Rosasco acknowledges the financial support of the European Research Council (grant SLING 819789), the AFOSR projects FA9550-18-1-7009, FA9550-17-1-0390 and BAA-AFRL-AFOSR-2016-0007 (European Office of Aerospace Research and Development), the EU H2020-MSCA-RISE project NoMADS - DLV-777826, and the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF-1231216.

§ PROOF OF THEOREM <REF>

Via (<ref>) and (<ref>), denoting $\rfK=SS^*+γI$ the (weighted) regularized kernel matrix, it holds

\begin{align*}
\G &= \Zx^*\rfK^{-1}S, \\
%&= \Zx^*S \rfcov^{-1} \\
%\nG &= \Pout\Zx^*S\Pin (\Pin S^*S\Pin + γ I)^{-1} \\
\nG %&= \Pout\Zx^*S g(\fcov) \\
&=\Pout\Zx^* h(S) S\Pin,
\end{align*}

\begin{equation*}
\de (S\Pin S^* + \gamma I)^{-1}.
%\de \Pin (\Pin \fcov \Pin + γ I)^{-1} %\\
%&= \Bfin (\Bfin^* \fcov \Bfin + γI)^{-1} \Bfin^*
\end{equation*}

In particular, we can observe that

\begin{align}
\n{\nG-\G}
&= \norm{\Zx^*\rfK^{-1} S  -  \Pout\Zx^* h(S) S\Pin}\nonumber\\
&≤  \norm{\Zx^*\rfK^{-1} S(I -  \Pin)}\nonumber \\
&\quad+\n{\Zx^* (\rfK^{-1}-h(S)) S\Pin}\nonumber\\%\\
&\quad+\n{(I-\Pout)\Zx^* h(S) S\Pin} .%\\
\label{e:bound_G_nG_dec}
\end{align}

The first addend can be split as follows:

\begin{align}
\norm{\Zx^* \rfK^{-1} S\Pin^\perp}
&≤ \norm{\Zx^* \rfK^{-1}}\norm{S\Pin^\perp}\nonumber\\
&≤ \frac\kappa\gamma\norm{S \Pin^\perp}. %\\
\label{e:bound_G_nG_t1}
\end{align}

By leveraging the identity $A^{-1} - B^{-1} = A^{-1}(B - A) B^{-1}$ for two invertible operators $A$ and $B$, the second addend in (<ref>) can be bounded as

\begin{align*}
&\norm{\Zx^* ( \rfK^{-1} - h(S) ) S\Pin}\nonumber\\
&\quad= \norm{\Zx^* \rfK^{-1} (S\Pin S^* - SS^*) h(S) S\Pin}\\
&\quad≤ \norm{\Zx^* \rfK^{-1}}\n{ (S\Pin S^* - SS^*)}\n{ h(S) S\Pin}\\
&\quad≤ \frac{κ}{γ} \n{\Pin^\perp S^*}² \n{ h(S) S\Pin }. %\\
%&≤ \frac{κ}{γ^{3/2}} \n{\Pin^\perp S^*}² \n{(S\Pin S^* + \gamma I)^{-1/2} S\Pin} %\\
\end{align*}

Note that, using a polar decomposition, $S\Pin = (S\Pin S^*)^{1/2}U$ for some partial isometry $U$, and thus

\begin{align}
\hspace{-0.3cm}\n{h(S) S\Pin }
&≤ \n{(S\Pin S^* + \gamma I)^{-1/2}}\nonumber\\
&\quad\cdot\!	\n{(S\Pin S^* + \gamma I)^{-1/2}(S\Pin S^*)^{1/2} }\nonumber\\
&≤ γ^{-1/2},\!\!\!\!
\label{e:bound_G_nG_hS}
\end{align}

so that

\begin{align}
\norm{\Zx^* ( \rfK^{-1} - h(S) ) S\Pin}
&≤ \frac{κ}{γ^{3/2}} \n{\Pin^\perp S^*}².
\label{e:bound_G_nG_t2}
\end{align}

Eventually, using again (<ref>) the third addend in (<ref>) can be bounded as

\begin{align}
\n{\Pout^\perp \Zx^* h(S) S\Pin} %\\
&≤ \n{\Pout^\perp \Zx^*}\n{ h(S) S\Pin} \nonumber\\
&≤ \n{\Pout^\perp \Zx^*} .%\\
\label{e:bound_G_nG_t3}
\end{align}

Putting together (<ref>), (<ref>) and (<ref>), we get

\begin{align}
\n{\nG-\G}
&≤ \frac{κ}{γ} \norm{\Pin^\perp S^*} + \frac{κ}{γ^{3/2}} \n{\Pin^\perp S^*}²\nonumber  \\
&\quad+ \n{\Pout^\perp \Zx^*}.
\end{align}

We recall that $\Pin = \vvec{\Pinx & 0 \\ 0 & I}$, so that $\n{\Pin^\perp S^*}=\n{\Pinx^\perp \Sx^*}$, and both terms $\n{\Pinx^\perp \Sx^*}$ and $\n{\Pout^\perp \Zx^*}$ can be bounded using Lemma <ref>, which is itself derived from [38]. We thus get

\begin{align}
\n{\nG-\G}
&≤ \prt*{\frac{κ}{γ}+1} 4\supfmap\sqrt{\frac{3}{m}\log\prt*{\frac{8m}{5δ}}} \nonumber\\
&\quad+ \frac{48\supfmap^3}{γ^{3/2}} \frac{1}{m}\log\prt*{\frac{8m}{5δ}}.
\end{align}

§ PROOF SKETCH FOR LEMMA <REF>

Let $L=A + BK$ be the closed-loop operator of the system specified by (<ref>), where $K$ is the stabilizing gain defined in (<ref>). The proof of the result, which closely follows the proof of <cit.>, starts by defining a set $\cS$ of perturbations of the Riccati operator $P$ as

\begin{equation*}
\cS = \{X:\cH_1\to\cH_1, X = X^*, P + X \succcurlyeq 0\}.
\end{equation*}

Then, we can derive a fixed-point operator equation

\begin{equation}
X = \Phi X,\label{e:fixed_point_eq_riccati}
\end{equation}

\begin{align*}
\Phi X &= \cTinv[F(X + P, A, B) - F(X + P, \Atilde, \Btilde) - \mathcal R X],\\
\cT X &= I - L^*XL,\\
\cR X &= L^*X[I + B^*PB(P + X)]^{-1}B^*PBXL.
\end{align*}

Let us define $\CLHo$ as the space of continuous linear operators from $\cH_1$ to $\cH_1$, and $\cL(\CLHo)$ as the space of continuous linear operators from $\CLHo$ to $\CLHo$. The operator $\cT:\CLHo\to\CLHo$ is invertible in $\cL(\CLHo)$. This fact can be proven by showing that $1$ is not in the spectrum of operator $\cD:\CLHo\to\CLHo$, defined as $\cD X = L^*XL$. This is achieved by upper bounding the spectral radius of $\cD$ by means of Gelfand's formula, and by observing that the closed-loop dynamics are stable, i.e., $L$ has spectral radius smaller than 1. Moreover, $\cTinv$ is bounded. This can be proven by showing that the Neumann series $\sum_{k=0}^\infty \cD^k$ is absolutely convergent in operator norm, and therefore equal to $\cTinv$.

Then, it can be shown that the unique solution of (<ref>) belonging to $\cS$ is $P-\Ptilde$, since, by construction of the fixed point equation, it must hold $F(X + P, \Atilde, \Btilde) = 0$.

We can now define the following closed subset of $\cS$

\begin{equation*}
\mathcal S_{\nu} = \{X :\norm{X}\leq \nu, X = X^*, P + X \succcurlyeq 0\}.
\end{equation*}

Let $\norm{\Atilde -A}\leq \epsilon$, $\norm{\Btilde -B}\leq \epsilon$, where $\epsilon$ is the error rate in Theorem <ref>. If we pick

\begin{align*}
\nu= 	\min\left\{\norm{N}^{-1}, \frac12\right., &6\left.\epsilon \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp A^2 \normp P^2\right.\\
&\left. \cdot\normp B \normp {R^{-1}} \right\},\nonumber
\end{align*}

for $X, X_1, X_2\in \cS_\nu$, we can show that $\Phi X\in \cS_\nu$, and $\exists \eta <1$ such that $\norm{\Phi X_1 - \Phi X_2} \leq \eta \norm{X_1- X_2}$, that is, $\Phi$ is a contraction operator. Then, we can apply Banach fixed-point theorem, meaning that the self-adjoint fixed point $\Ptilde - P$ of $\Phi$ belongs to $\cS_{\nu}$, i.e., the error on the Riccati operator due to the Nyström approximation has bounded operator norm as a function of $\epsilon$.

§ CLOSED FORM FOR $\G$

We recall that we defined the risk in (<ref>) as

\begin{align*}
\cR(W) &\de \frac1n \sum_{i=1}^{n}\n*{ \psi(\xvec_{i+1}) - W ϕ(\wvec_i)}_{\cH_1}².
\end{align*}

*Expression of the risk
For any orthonormal basis $(g_j)_{j∊\mathbb N}$ of $\Hspace_1$, it holds

\begin{align*}
\mathcal R&(W)	\\
&= \frac1n \sum_{i=1}^{n} \left[ \sum_{j\in\mathbb N}
\innerprod{g_j, \psi(\xvec_{i+1})-W\phi(\wvec_i)}_{\mathcal H_1}^2\right] \\
%&= \frac1n\sum_{j\in\mathbb N}\left[\sum_{i=1}^{n}\left(\innerprod{g_j, \psi(\xvec_{i+1})}_{\Hspace_1}-\innerprod{g_j, G\phi(\wvec_i)}_{\Hspace_1}\right)^2\right] \\
&= \frac1n\sum_{j\in\mathbb N}\left[\sum_{i=1}^{n}\left(\innerprod{g_j, \psi(\xvec_{i+1})}_{\Hspace_1}-\innerprod{W^*g_j, \phi(\wvec_i)}_{\Hspace}\right)^2\right] \\
&= \frac1n\sum_{j\in\mathbb N}\norm{\sqrt n\Zx g_j-\sqrt nSW^*g_j}^2_{\R^n} \\
&= \nHS{\Zx - SW^*}². %\\
%&= \Tr(\Zx\Zx^* - S^*\Zx G - G^*\Zx^*S  + GS^*S G^*  )
\end{align*}

*Expression of $\G$ Recall $\G=\argmin_{W:\cH→\cH_1}\cR(W)+\gamma\norm{W}_{\HS}$.
Computing the gradient of $\cR(W)$ gives

\begin{equation*}
\pd{\mathcal{R}(W)}{W} = - 2\Zx^*S + 2WS^*S.
\end{equation*}

The first-order optimality condition with additive Tikhonov regularization leads to

\begin{align}
\G (S^*S + \gamma I) &= \Zx^*S, \nonumber\\
\G &= \Zx^*S (S^*S + \gamma I)^{-1} \label{eq:a_b_compound_expr}.%\\
%&= \Zx^*(SS^* + \gamma I)^{-1}S.
\end{align}

§ ERROR INDUCED BY THE NYSTRÖM APPROXIMATION

The following lemma is a restatement of <cit.> with slightly different constants. Albeit being originally written in a random design scenario, the same result holds when conditioning on the data.

Let $x₁,…,x_n∊\bR^d$ and denote $\cC=\tfrac{1}{n} \sum_{i=1}^n ψ(x_i)\kron ψ(x_i)$, where $ψ$ is the canonical feature map associated to a kernel satisfying Assumption <ref>.
Let $\ldm{1},…,\ldm{m}$ be drawn uniformly from all partitions of cardinality $m$ of $\cb{x₁,…,x_n}$.
Denoting $\Pi_m$ the orthogonal projection onto $\spa(ψ(\ldm{1}),…,ψ(\ldm{m}))$,
for any $γ∊]0,\n{\cC}]$, it holds

\[ \n{(I-\Pi_m) (\cC+γI)^{1/2}}^2 ≤ 3γ,\]

with probability at least $1-δ$ provided

\[ m\geq ( 2+ 5 \tfrac{κ²}{γ} )\log\prt*{\tfrac{4κ²}{γδ}}. \]

We apply exactly the same proof as in <cit.> but conditioning on the data (i.e., in particular applying <cit.> with the $v_i$ drawn according to the empirical distribution and $Q$ being the empirical covariance), with $\cN_∞(γ)=κ²/γ$.
Denoting $w=\log\prt*{\tfrac{4\supk}{γδ}}$, we end up with the bound

\begin{equation*}
\n{(I-\Pi_m)(\cC+γI)^{1/2}}² ≤ \frac{γ}{1-β(γ)} ,
\end{equation*}

β(γ) ≤ \frac{2w}{3m} + \sqrt{\frac{2w\supk}{γ m} }.

We now derive slightly better constants than in the original lemma in order to satisfy $β(γ)≤2/3$ (which gives the claimed result). Indeed,

\begin{align}
β(γ) ≤ 2/3
%&⇔ \frac{2w}{3} + \sqrt{\frac{2w\supk}{γ} } \sqrt{m} ≤ 2/3 m \\
%&⇔ w + \frac{3}{2} \sqrt{\frac{2w\supk}{γ} } \sqrt{m} ≤ m \\
&⇔ m - \frac{3}{2} \sqrt{\frac{2w\supk}{γ}} \sqrt{m} - w ≥ 0. %\\
\end{align}

We solve this inequality as a second-order equation in $\sqrt{m}$, with discriminant $Δ=w(\frac{9\supk}{2γ}+4)$, which is positive whenever $γ<4\supk δ^{-1}$.
A sufficient condition to satisfy $β(γ)≤2/3$ is thus

\begin{align}
m &≥ \prt*{\frac{3}{4} \sqrt{\frac{2w\supk}{γ}}+ \frac{1}{2} \sqrt{w(\frac{9\supk}{2γ}+4)} }².
\end{align}

Using the identity $2(a²+b²)≥(a+b)²$, we get the sufficient condition

%&≥ 2\frac{9}{16} \frac{2w\supk}{γ} + 2\frac{1}{4} w(\frac{9\supk}{2γ}+4) \\
%&≥ \frac{9}{4} \frac{w\supk}{γ} + \frac{1}{2} w (\frac{9\supk}{2γ}+4) \\
%&≥ w (\frac{9}{2} \frac{\supk}{γ} + 2)
≥ w \prt*{ 5 \frac{\supk}{γ} + 2 }.%\nonumber.

Under Assumption <ref>,
provided that $\Pinx$ and $\Pout$ are built by sampling uniformly $m$ samples respectively from $\cb{\bf x₁,…,\bf x_n}$ and $\cb{\bf x₂,…,\bf x_{n+1}}$, it holds with probability $1-δ$

\begin{align}
\max(\n{\Pinx^\perp \Sx^*}, \n{\Pout^\perp \Zx^*})
&≤ 4\supfmap\sqrt{\frac{3}{m}\log\prt*{\frac{8m}{5δ}}}\nonumber .
\end{align}

Using a polar decomposition and applying Lemma <ref> with $m ≥ ( 2+ 5 \tfrac{κ²}{γ} ) \log\prt*{\tfrac{4κ²}{γδ}}$, it holds

\begin{align}
\n{\Pinx^\perp \Sx^*}
&= \n{\Pinx^\perp (\Sx^*\Sx)^{1/2}}  \\
&≤ \n{\Pinx^\perp (\Sx^*\Sx+γI)^{1/2}} %\\
≤ \sqrt{3γ},
\end{align}

and a similar argument holds for $\n{\Pout^\perp \Zx^*}$.

The condition
m %&
≥ ( 2+ 5 \tfrac{κ²}{γ} ) \log\prt*{\tfrac{4κ²}{γδ}}
is in particular satisfied provided that

m/2 ≥ \max(2, 5 \tfrac{κ²}{γ}) \log\prt*{\tfrac{4κ²}{γδ}}

which holds taking $γ=\frac{16κ²}{m}\log(\tfrac{4m}{5δ})$.
This yields the claimed result via a union bound.

§ COMPUTABLE EXPRESSION FOR NYSTRÖM VECTOR DYNAMICS AND $C$

In this appendix, we report the computable expression for the operator

in (<ref>), as well as the expression of the reconstruction matrix $C$ in (<ref>).
They are the formulas used in practice for the implementation of the method. To simplify the notation, for operators $M, N$, we denote a block-diagonal operator as $\bd{M, N} = \vvec{M &0\\ 0&N}$.

We introduce the following additional notation. In the sequel we refer to the nonlinear kernel $k$:

* $K_{nm, \textrm{out}}$ is the kernel between the regression outputs and the output landmarks;

* $K_{nm, \textrm{in}}$ is the kernel between the regression inputs (state only) and the input landmarks;

* $K_{m, \textrm{in}}$ is the kernel at the input landmarks;

* $K_{m, \textrm{in},\textrm{out}}$ is the kernel between the input landmarks and the output landmarks.

Let us define the sampling operator for the input Nyström landmarks as

\begin{align*}
\Stildein:\Hspace_1\rightarrow\R^m,\
%\Stildeout g &= m^{-1/2}[\innerprod{g, \kr(\cdot, \tildexvec_2)}_{\Hspace}, \dots,\innerprod{g, \kr(\cdot, \tildexvec_{m+1})}_{\Hspace}]^T.
\Stildein g &= [g(\ldmi{1}), …,g(\ldmi{m})]ᵀ.
\end{align*}

Note that $\ran(\Stildein^*)=\Hin$.

We denote
$\Stildein = U\Sigma V^*$ its
singular value decomposition, where $U:\R^t\rightarrow\R^m$, $\Sigma:\R^t\rightarrow\R^t$, $V:\R^t\rightarrow \tilde \cH_1$.

Here, $\Sigma$ is diagonal and strictly positive. It holds $\Pin = \bd{VV^*, I}$ and thus

\begin{align*}
%\Bout^*&\nG \bd{\Bout, I}\\
%&= \Bout^* \Pout\Zx^* (S\Pin S^* + \gamma I)^{-1} S\bd{\Pinx & 0 \\ 0 & I } \bd{ \Bout , I } \\
&\Bout^* \Zx^* (S\Pin S^* + \gamma I)^{-1} S\bd{\Pinx\Bout , I } \\
&= \Bout^* \Zx^*\prt*{S \bd{VV^* , I } S^* + \gamma I}^{-1} \\
&\quad\cdot S\bd{VV^* , I } \bd{\Bout , I }\\
&= \Bout^* \Zx^*S\bd{V , I }  \prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}^{-1} \\
&\quad\cdot\bd{V^* , I } \bd{\Bout , I }\\
%	&= \Bout^* \Zx^*S\bd{V\Sigma U^* , I }\bd{U\Sigma^{-1} , I } \\
%	&\quad\cdot \prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}^{-1}\\
%	&\quad\cdot\bd{\Sigma^{-1}U^*, I } \bd{U\Sigma V^* , I } \bd{\Bout , I }\\
&= \Bout^* \Zx^*S\bd{\Stildein^* , I } \bd{U\Sigma^{-1} , I } \\
&\quad\cdot\prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}^{-1}\\
&\quad\cdot \bd{\Sigma^{-1}U^*, I } \bd{\Stildein , I }\bd{\Bout , I }.
\end{align*}

Now, we can observe that $V^*V=I$. Moreover, $\prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}$ is full-rank, $\bd{U\Sigma^{-1} ,I }$ is full column-rank and $\bd{\Sigma^{-1} U^*, I }$ is full row-rank,
then by <cit.> we have that

\begin{align*}
&\bd{U \Sigma^{-1}, I }\prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}^{-1}\bd{\Sigma^{-1}U^*, I }\\
&\quad= \bd{\Sigma U^*, I }^{\dagger} \prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}^{-1}\bd{U\Sigma, I }^\dagger\\
%&\quad= \bd{\Sigma U^*, I }^{\dagger}\left\{ \bd{U\Sigma, I }\prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}\right\}^\dagger\\
&\quad= \left\{ \bd{U\Sigma, I }\prt*{\bd{V^* , I } S^*S \bd{V , I }+ \gamma I}\bd{\Sigma U^*, I }\right\}^\dagger\\
&\quad=  \left\{ \prt*{\bd{U\Sigma V^* , I } S^*S \bd{V\Sigma U^* , I }+ \gamma \bd{U\Sigma\Sigma U^*, I }}\right\}^\dagger\\
&\quad= \prt*{\bd{\Stildein , I } S^*S \bd{\Stildein^* , I }+ \gamma \bd{\Stildein\Stildein^*, I }}^\dagger.
\end{align*}

This yields

\begin{align*}
&\Bout^*\nG \bd{\Bout , I}\\
&\quad= \Bout^* \Zx^*S\bd{\Stildein^* , I }\\
&\qquad \cdot \prt*{\bd{\Stildein , I } S^*S \bd{\Stildein^* , I }+ \gamma \bd{\Stildein\Stildein^*, I }}^\dagger\\
&\qquad \cdot\bd{\Stildein , I }\bd{\Bout , I }\\
&\quad= \Bout^* \Zx^*\vvec{\Sx\Stildein^* &\Su}\\
&\qquad\cdot \prt*{\vvec{\Stildein\Sx^* \\ \Su^T} \vvec{\Sx\Stildein^{*} &\Su}+ \gamma \bd{\Stildein\Stildein^*, I }}^\dagger \\
&\qquad\cdot\bd{\Stildein , I }\bd{\Bout , I }\\
&\quad= (K_{m, \textrm{out}}^{\dagger})^{1/2}K_{mn, \textrm{out}}\vvec{K_{nm, \textrm{in}} &\sqrt n S_{\uvec}}\\
&\qquad\cdot\prt*{\vvec{K_{mn, \textrm{in}}\\ \sqrt n\Suᵀ}
\vvec{ K_{nm, \textrm{in}} &\sqrt n \Su }
\!+\! \gamma n \bd{
K_{m,\textrm{in}} , I } }^{\dagger} \\
&\qquad\cdot\bd{K_{m, \textrm{in},\textrm{out}} (K_{m, \textrm{out}}^†)^{1/2}, I}.
\end{align*}

Moreover, having stacked the regression outputs in a matrix $X_{+1}\in \R^{d\times n}$, the matrix $C$ from (<ref>) takes the form

\begin{align*}
C &= X_{+1}K_{nm, \textrm{out}}\\
&\quad\cdot(K_{mn, \textrm{out}}K_{nm, \textrm{out}} + \lambda n K_{m, \textrm{out}})^{-1}K_{m,\textrm{out}}^{1/2}.
\end{align*}

[1]
Ian Abraham, Gerardo de la Torre, and Todd Murphey.
Model-based control using Koopman operators.
Robotics: Science and Systems XIII, 2017.

[2]
Ian Abraham and Todd D Murphey.
Active learning of dynamics for data-driven control using Koopman
IEEE Transactions on Robotics, 35(5):1071–1083, 2019.

[3]
Tamim El Ahmad, Luc Brogat-Motte, Pierre Laforgue, and Florence d'Alché
Sketch in, sketch out: Accelerating both learning and inference for
structured prediction with kernels.
arXiv preprint arXiv:2302.10128, 2023.

[4]
Fabio Amadio, Juan Antonio Delgado-Guerrero, Adrià Colomé, and Carme
Controlled Gaussian process dynamical models with application to
robotic cloth manipulation.
International Journal of Dynamics and Control, 11:3209–3219,

[5]
Nachman Aronszajn.
Theory of reproducing kernels.
Transactions of the American Mathematical Society,
68(3):337–404, 1950.

[6]
Adi Ben-Israel and Thomas N. E. Greville.
Generalized Inverses: Theory and Applications.
CMS Books in Mathematics. Springer, 2 edition, 2003.

[7]
Alain Bensoussan, Giuseppe Da Prato, Michel C Delfour, and Sanjoy K Mitter.
Representation and control of infinite dimensional systems,
volume 1.
Springer, 2007.

[8]
Alain Berlinet and Christine Thomas-Agnan.
Reproducing kernel Hilbert spaces in probability and
Springer Science & Business Media, 2011.

[9]
Petar Bevanda, Max Beier, Armin Lederer, Stefan Sosnowski, Eyke
Hüllermeier, and Sandra Hirche.
Koopman kernel regression.
Advances in Neural Information Processing Systems, 36, 2024.

[10]
Petar Bevanda, Stefan Sosnowski, and Sandra Hirche.
Koopman operator dynamical models: Learning, analysis and control.
Annual Reviews in Control, 52:197–212, 2021.

[11]
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.

[12]
Franco Coltraro, Jaume Amorós, Maria Alberich-Carramiñana, and Carme
An inextensible model for the robotic manipulation of textiles.
Applied Mathematical Modelling, 101:832–858, 2022.

[13]
Suddhasattwa Das and Dimitrios Giannakis.
Koopman spectra in reproducing kernel hilbert spaces.
Applied and Computational Harmonic Analysis, 49(2):573–607,

[14]
Sarah Dean, Horia Mania, Nikolai Matni, Benjamin Recht, and Stephen Tu.
On the sample complexity of the linear quadratic regulator.
Foundations of Computational Mathematics, 20(4):633–679, 2020.

[15]
Anthony M. DeGennaro and Nathan M. Urban.
Scalable Extended Dynamic Mode Decomposition using Random
Kernel Approximation.
SIAM Journal on Scientific Computing, 41(3):A1482–A1499, 2019.

[16]
Bas JL Driessen.
Koopman operator regression based surrogate modelling for control.
Master's thesis, Eindhoven University of Technology, November 2023.

[17]
Maryam Fazel, Rong Ge, Sham Kakade, and Mehran Mesbahi.
Global convergence of policy gradient methods for the linear
quadratic regulator.
In International Conference on Machine Learning, pages
1467–1476. PMLR, 2018.

[18]
Dimitris Giannakis, Amelia Henriksen, Joel A Tropp, and Rachel Ward.
Learning to forecast dynamical systems from streaming data.
SIAM Journal on Applied Dynamical Systems, 22(2):527–558,

[19]
Andrew J Gibson, Michael L Calvisi, and Xin C Yee.
Koopman linear quadratic regulator using complex eigenfunctions for
nonlinear dynamical systems.
SIAM Journal on Applied Dynamical Systems, 21(4):2463–2486,

[20]
William W Hager and Larry L Horowitz.
Convergence and stability properties of the discrete Riccati
operator equation and the associated optimal control and filtering problems.
SIAM Journal on Control and Optimization, 14(2):295–312, 1976.

[21]
David A Haggerty, Michael J Banks, Ervin Kamenar, Alan B Cao, Patrick C Curtis,
Igor Mezić, and Elliot W Hawkes.
Control of soft robots with inertial dynamics.
Science Robotics, 8(81):eadd6864, 2023.

[22]
Wenjian Hao, Bowen Huang, Wei Pan, Di Wu, and Shaoshuai Mou.
Deep Koopman learning of nonlinear time-varying systems.
Automatica, 159:111372, 2024.

[23]
Eurika Kaiser, J Nathan Kutz, and Steven L Brunton.
Data-driven discovery of Koopman eigenfunctions for control.
Machine Learning: Science and Technology, 2(3):035023, 2021.

[24]
Sham Kakade and John Langford.
Approximately optimal approximate reinforcement learning.
In International Conference on Machine Learning, pages
267–274. PMLR, 2002.

[25]
Mohammad Khosravi.
Representer theorem for learning Koopman operators.
IEEE Transactions on Automatic Control, 68(5):2995–3010, 2023.

[26]
Michail M Konstantinov, P Hr Petkov, and Nikolai D Christov.
Perturbation analysis of the discrete Riccati equation.
Kybernetika, 29(1):18–29, 1993.

[27]
Bernard O Koopman.
Hamiltonian systems and transformation in Hilbert space.
Proceedings of the National Academy of Sciences,
17(5):315–318, 1931.

[28]
Milan Korda and Igor Mezić.
Linear predictors for nonlinear dynamical systems: Koopman operator
meets model predictive control.
Automatica, 93:149–160, 2018.

[29]
Milan Korda and Igor Mezić.
Optimal construction of Koopman eigenfunctions for prediction and
IEEE Transactions on Automatic Control, 65(12):5114–5129,

[30]
Vladimir Kostic, Pietro Novelli, Andreas Maurer, Carlo Ciliberto, Lorenzo
Rosasco, and Massimiliano Pontil.
Learning dynamical systems via Koopman operator regression in
reproducing kernel Hilbert spaces.
In Advances in Neural Information Processing Systems,
volume 35, pages 4017–4031, 2022.

[31]
Adrià Luque, David Parent, Adrià Colomé, Carlos Ocampo-Martinez, and Carme
Model predictive control for dynamic cloth manipulation: Parameter
learning and experimental validation.
IEEE Transactions on Control Systems Technology, pages 1–17,

[32]
Horia Mania, Stephen Tu, and Benjamin Recht.
Certainty equivalence is efficient for linear quadratic control.
In Advances in Neural Information Processing Systems,
volume 32, pages 10154–10164, 2019.

[33]
Giacomo Meanti, Antoine Chatalic, V. Kostic, P. Novelli, M. Pontil, and
L. Rosasco.
Estimating Koopman operators with sketching to provably learn
large scale dynamical systems.
In Advances in Neural Information Processing Systems, 2023.

[34]
Igor Mezić.
Koopman operator, geometry, and learning of dynamical systems.
Not. Am. Math. Soc., 68(7):1087–1105, 2021.

[35]
Joseph Moyalan, Hyungjin Choi, Yongxin Chen, and Umesh Vaidya.
Data-driven optimal control via linear transfer operators: A convex
Automatica, 150:110841, 2023.

[36]
Samuel E. Otto and Clarence W. Rowley.
Koopman Operators for Estimation and Control of
Dynamical Systems.

[37]
Sebastian Peitz and Stefan Klus.
Koopman operator-based model reduction for switched-system control of
Automatica, 106:184–191, 2019.

[38]
Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco.
Less is more: Nyström computational regularization.
In Advances in Neural Information Processing Systems,
volume 28, pages 1657–1665, 2015.

[39]
Bernhard Schölkopf, Ralf Herbrich, and Alex J Smola.
A generalized representer theorem.
In International conference on computational learning theory,
pages 416–426. Springer, 2001.

[40]
Bernhard Schölkopf and Alexander J Smola.
Learning with kernels: support vector machines, regularization,
optimization, and beyond.
MIT press, 2002.

[41]
Haojie Shi and Max Q-H Meng.
Deep Koopman operator with control for nonlinear systems.
IEEE Robotics and Automation Letters, 7(3):7700–7707, 2022.

[42]
Max Simchowitz, Horia Mania, Stephen Tu, Michael I Jordan, and Benjamin Recht.
Learning without mixing: Towards a sharp analysis of linear system
In Conference on Learning Theory, pages 439–473. PMLR, 2018.

[43]
Ingo Steinwart and Andreas Christmann.
Support vector machines.
Springer Science & Business Media, 2008.

[44]
Larry Wasserman.
All of nonparametric statistics.
Springer Science & Business Media, 2006.

[45]
Christopher Williams and Matthias Seeger.
Using the Nyström method to speed up kernel machines.
In Advances in Neural Information Processing Systems,
volume 13, pages 682–688, 2000.

[46]
Hang Yin, Michael C Welle, and Danica Kragic.
Embedding Koopman optimal control in robot policy learning.
In IEEE/RSJ International Conference on Intelligent Robots and
Systems (IROS), pages 13392–13399. IEEE, 2022.

[47]
Ce Xu Zheng, Adrià Colomé, Luis Sentis, and Carme Torras.
Mixtures of controlled Gaussian processes for dynamical modeling of
deformable objects.
In Learning for Dynamics and Control Conference, pages
415–426. PMLR, 2022.

§ PROOF OF LEMMA <REF>

Note that in our case the operators $Q$ and $R$ are the same both for the exact and Nyströmized dynamics, and are assumed to be positive definite. The structure of the proof is as follows. First, we write down a fixed-point operator equation, $X=\Phi X$, for a suitable operator $\Phi$, that has a unique solution, given by the error on the Riccati operator due to the Nyström approximation, $\Ptilde - P$. Then, we define a set of perturbations $X$ with bounded norm, and we show that $\Phi$ maps this set to itself, and is a contraction. This implies that the fixed point of $\Phi$, which is shown to be $\Ptilde-P$, lies in the set we defined, and has bounded norm. We begin by defining

\begin{align}
N:\cH_1\to\cH_1,\ N &= B R^{-1}B^*,
\label{e:def_N} \\
\Ntilde:\cH_1\to\cH_1,\ \Ntilde &= \Btilde R^{-1}\Btilde^*.
\label{e:def_nN}
\end{align}

Given $P$ as the self-adjoint positive semi-definite solution of (<ref>), let us consider the following set of self-adjoint perturbations

\begin{equation}
\cS = \{X:\cH_1\to\cH_1, X = X^*, P + X \succcurlyeq 0\}.
\end{equation}

Let $X \in \cS$. Let  $L = A + BK$, where $K$ is defined in (<ref>). Then, it holds

\begin{equation}
F(P + X, A, B) = \underbrace{X - L^*XL }_{\eqqcolon \mathcal T X}+\underbrace{ L^*X[I + N(P + X)]^{-1}NXL}_{\eqqcolon \mathcal R X}.\nonumber\label{eq:perturbed_riccati}
\end{equation}

This statement can be found, in the matrix setting, in [32, 26]. We report here the derivation for completeness. By definition of $P$ it holds $F(P,A,B)=0$ thus:

\begin{align}
&X - L^*XL + L^*X[I + N(P + X)]^{-1}NXL\nonumber \\
&= X - L^*XL + L^*X[I + N(P + X)]^{-1}NXL + F(P,A,B)\\
&= X - L^*XL + L^*X[I + N(P + X)]^{-1}NXL
+ P-A^*P(I+NP)^{-1}A -Q \\
&= (P+X)-Q + \brk*{ L^*X[I + N(P + X)]^{-1}NXL
- L^*XL
-A^*P(I+NP)^{-1}A } \\
&= (P+X)-Q + \brk*{
L^*X\prt*{[I + N(P+X)]^{-1}NX-I}L
} \\
&= (P+X)-Q + \Big[
L^*X\prt*{[I + N(P+X)]^{-1}NX-[I + N(P+X)]^{-1}[I + N(P+X)]}L \\
&\quad -A^*P(I+NP)^{-1}A
\Big] \nonumber\\
&= (P+X)-Q - \brk*{
L^*X\prt*{[I + N(P+X)]^{-1}[I + NP]}L
} \\
&= (P+X)-Q - A^*\brk*{
\\
&= (P+X)-Q - A^*\brk*{
\\
&= (P+X)-Q - A^*\brk*{
\\
&= (P+X)-Q - A^*(P+X)(I+N(P+X))^{-1}A\\
&= F(P+X, A, B),
\end{align}

which concludes the proof.

We can now study the operator $\mathcal T$ appearing in Proposition <ref>.

Let $\CLHo$ denote the set of continuous linear operators from $\cH₁$ to itself.
The operator $\mathcal T: \CLHo→\CLHo$ appearing in Proposition <ref> is invertible.
Under assumptions <ref> and <ref>, the fixed point equation

\begin{equation}
X = \mathcal T^{-1} [F(X + P, A, B) - F(X + P, \Atilde, \Btilde) - \mathcal R X] \nonumber
\end{equation}

has a unique solution in $\cS$, namely $X = P-\Ptilde$.

The operator $\mathcal T$ can be rewritten as

\begin{equation}
\mathcal T X =( I - \mathcal D)X
\end{equation}

where $\mathcal D X = L^*XL:\CLHo\to\CLHo$ is linear ($\mathcal D (X + Y) = L^*XL + L^*YL,  \mathcal D (\alpha X) =\alpha L^*XL$).
Let $\cL(\CLHo)$ denote the set of continuous linear operators from $\CLHo$ to itself.
Then $\mathcal T∊\cL(\CLHo)$ (because $A, B, Q, R$, and $K$, defined in (<ref>), appearing in the expression of $L$ are bounded).
We denote the spectrum of $\mathcal D$ as

\begin{equation}
\sigma(\mathcal D) = \{\lambda : \mathcal D - \lambda I \text{ is not  invertible in $\cL(\CLHo)$}\},
\end{equation}

and its spectral radius $\rho(\mathcal D) = \sup_{\lambda \in \sigma(\mathcal D)} |\lambda|$.
The operator $\mathcal T$ is invertible in $\cL(\CLHo)$ iff 1 is not in the spectrum $\sigma(\mathcal D)$, and we now show that this holds by proving $ρ(\cD)<1$.

We can observe that since $K$ in (<ref>) is a stabilizing gain by assumption, it holds <cit.> $\rho(L) <1$. The spectral radius of $\cD $ can be computed with the Gelfand's formula,
using the fact that $L^k, (L^*)^k$ are bounded for any $k$:

\begin{align}
\rho(\mathcal D) &= \lim_{k\to\infty}\norm{\mathcal D^k}^{1/k} \\
&=  \lim_{k\to\infty}\sup_{\norm X = 1}\norm{\mathcal D^k X}^{1/k} \\
%&=  \lim_{k\to\infty}\sup_{\norm X = 1}\norm{ \underbrace{L^*\dots L^*}_{k \text{ times}}
%X\underbrace{L\dots L}_{k \text{ times}}}^{1/k} \\
&= \lim_{k\to\infty}\sup_{\norm X = 1}\norm{ {L^*}^k	XL^k}^{1/k} \\
&\leq \lim_{k\to\infty}\norm{ {L^*}^k}^{1/k}\norm{L^k}^{1/k}\\
&= \rho(L)^2 \\
&< 1.
\end{align}

Now, we can consider the following equation:

\begin{equation}
\label{eq:perturbed_riccati_with_nystrom}
F(X + P, A, B) - F(X + P, \Atilde, \Btilde) = \mathcal T X + \mathcal R X.
\end{equation}

By Proposition <ref> we have

\begin{equation}
F(X + P, A, B )= \cT X + \cR X,
\end{equation}

and thus it must hold

\begin{align}
F(X + P, \Atilde, \Btilde) &= 0.
\label{e:prop_nP}
\end{align}

By definition of $\nP$ in (<ref>),
$X=\nP-P$ satisfies this equation,
but it is also known <cit.> that under Assumptions <ref> and <ref>, (<ref>) admits a unique self-adjoint positive semi-definite solution. Moreover as both $P$ and $\nP$ are self-ajdoint, $X=\nP-P$ is self-adjoint as well and thus the unique solution of (<ref>) in $\cS$.

Also, note that (<ref>), for a suitable operator $\Phi:\cS\to\CLHo$, can be rewritten as

\begin{equation}
X = \underbrace{\mathcal T^{-1} [F(X + P, A, B) - F(X + P, \Atilde, \Btilde) - \mathcal R X]}_{\eqqcolon \Phi X},
\label{e:def_Phi}
\end{equation}

so that $\nP-P$ is the unique fixed point of $Φ$.

We have shown that the operator $\Phi$ is well defined and has a unique fixed point in $\cS$ equal to $\Ptilde - P$, and we will now show that the latter is upper bounded by $\epsilon$ in operator norm, where $\epsilon$ is the error on the matrices of the dynamical system, similarly to [32].
As mentioned at the beginning of this section, in the following, we will define a set $\cS_{\nu}$ of perturbations $X\in\cS$ with bounded norm. Furthermore, we will show that the function $\Phi$ is a mapping from $\cS_ν$ to itself, and that it is a contraction (i.e. $\alpha$-Lipschitz, $\alpha < 1$). Therefore, since $\cS_\nu$ is a closed subset of $\cS$, by the Banach fixed point theorem, $\Phi$ has a fixed point in $\cS_{\nu}$. Since the fixed point of $\Phi$ in $\cS$ has been proven to be $ \Ptilde - P$, this means that the error on the $P$ operator due to the Nyström approximation is bounded, and the bound depends on the error rate of the system's matrices.
More formally, let us define:

\begin{equation}
\mathcal S_{\nu} = \{X :\norm{X}\leq \nu, X = X^*, P + X \succcurlyeq 0\}.
\end{equation}

We can now prove some technical bounds that are used in [32]. We can define

\begin{align}
\Delta A &\coloneqq \nA - A,\\
\Delta B &\coloneqq \nB - B,\\
\Delta N &\coloneqq \nN - N,\\
\Delta P &\coloneqq \nP - P.
\end{align}

where we recall that $N,\nN$ are defined in <ref>, <ref>.

Assume $\norm{\Delta A}\leq \epsilon$,  and $\norm {\Delta B} \leq \epsilon$. Then, for $\epsilon\leq \norm{B}$, the following bound holds:

\begin{equation}
\norm{\Delta N}\leq 3\epsilon \norm{R^{-1}}\norm{B}.\nonumber
\end{equation}

We can observe that:

\begin{align}
\norm{\Delta N}& = \norm{BR^{-1} B^* - \Btilde R^{-1} \Btilde^*}\\
&=  \norm{BR^{-1} B^* - \Btilde R^{-1} B^* + \Btilde R^{-1} B^*-\Btilde R^{-1} \Btilde^*}\\
&= \norm{(B - \Btilde)R^{-1} B^*  + \Btilde R^{-1} (B^*-\Btilde^*)}\\
&\leq \norm{(B - \Btilde)R^{-1} B^*}  + \norm{\Btilde R^{-1} (B^*-\Btilde^*)}\\
&\leq \norm{B- \Btilde}\norm{R^{-1}}(\norm{B} + \norm{\Btilde})\\
&\leq \epsilon \norm{R^{-1}}(\norm{B} + \norm{\Btilde}).
\end{align}

We can observe that, by choosing $\epsilon\leq \norm B$, by the reversed triangle inequality:

\begin{align}
\norm{\Btilde} - \norm{B}
& \leq \norm{\Btilde - B}
\leq \epsilon
\leq \norm{B},
\end{align}

meaning that

\begin{equation}
\norm{\Delta N}\leq 3\epsilon \norm{R^{-1}}\norm{B},
\end{equation}

as also shown in [32].

Let $L$, $\cT$ be as defined in Proposition <ref>. Moreover, let

\begin{equation}
\tau(L, \rho)=\sup\{\norm {L^k}\rho^{-k} : k≥ 0\}.\nonumber
\end{equation}

By assumption $K$ stabilizes the system, $ρ(L)<1$.
For any $\rho(L)\leq \rho<1$ we have that

\begin{equation}
\norm{\cTinv} \leq \frac{ \tau(L, \rho)^2}{1-\rho^2}.\nonumber
\end{equation}

By definition we have that operator $\cT$ can be expressed as $\cT= I - \cD$,
where $\cD X = L^*XL$.
Note that $L$ is bounded because $A,B$ and $K$ are bounded,
and thus $\cD:\CLHo→\CLHo$ is also bounded as

\begin{align}
\norm{\cD} &= \sup_{\norm X = 1}\norm{\cD X}%\\
= \sup_{\norm X = 1}\norm{L^* XL}%\\
\leq \norm L^2.
\end{align}

Now, observe that, since the spectral radius of $L$ is strictly smaller than 1 (stability), we can choose $\rho(L)< \rho<1$ so that the quantity $\tau(L, \rho)$ is finite [32].
Indeed according to Gelfand's formula, we are guaranteed that there exists $k_0$ such that for any $k\geq k_0$, $\norm{L^k}\leq \rho^k$.
Hence, we can choose

\begin{equation}
\tau'(L, \rho) =\max \cb{1} \cup \cb{ \norm{L^k}/\rho^k : k∊\irange{0,…,k_0}}
\end{equation}

which involves a finite sequence of values, so that for any $k∊\bN$ it holds $\n{L^k}≤\tau'(L, \rho) ρ^k$, and $τ(L,ρ)≤τ'(L,ρ)<∞$.

The Neumann series $\sum_{k=0}^{\infty} \cD^k$ is (absolutely) convergent in operator norm, as

\begin{align}
\sum_{k=0}^{\infty}\norm{\cD^k}
&=\sum_{k=0}^{\infty} \sup_{\norm X=1}\norm{{L^*}^k X L^k}\\
&\leq \sum_{k=0}^{\infty} \norm{L^k}^2\\
&\leq \sum_{k=0}^{\infty}\left[\tau(L, \rho)\rho^k\right]^2\\
&=\frac{ \tau(L, \rho)^2}{1-\rho^2},
\end{align}

since we have a geometric series with $\rho^2 < 1$.
Therefore, we can express $\cTinv$ using the Neumann series as

\begin{align}
\cTinv &= (I - \cD)^{-1}
= \sum_{k=0}^{\infty}\cD^k.
\end{align}

This yields the claimed result by the triangle inequality.

Let $\cR$, $L$, $N$, $X$ be the operators defined in Proposition <ref>. Then, we have

\begin{equation}
\norm{\cR X}\leq \norm{L}^2\norm{X}^2\norm{N}.\nonumber
\end{equation}

The first part of the proof is the same as the one of Lemma 7 in [32]. Note that $N$ and $P+X$ are self-adjoint, positive semi-definite operators.
Let us now consider $\alpha >0$. We have that

\begin{align}
0 \preccurlyeq
(N+\alpha I)[I + (P + X+\alpha I)(N+\alpha I)]^{-1}&=[(N+\alpha I)^{-1} + P + X+\alpha I]^{-1}\\
&\preccurlyeq(N+\alpha I).
\end{align}

This implies the following norm inequality:

\begin{equation}
\norm{(N+\alpha I)[I + (P + X+\alpha I)(N+\alpha I)]^{-1}}\leq\norm{N+\alpha I}.
\end{equation}

By definition, since limits preserve inequalities and by continuity of the norm, we have that

\begin{align}
\norm{N[I + (P + X)N]^{-1}}
&=\norm{\lim_{\alpha\to0} (N+\alpha I)[I + (P + X+\alpha I)(N+\alpha I)]^{-1}} \\
&=\lim_{\alpha\to0}\norm{(N+\alpha I)[I + (P + X+\alpha I)(N+\alpha I)]^{-1}}\\
&\leq \lim_{\alpha\to0}\norm{N+\alpha I}\\
\end{align}

\begin{align}
\norm{\cR X}&=\norm{L^*X[I + N(P + X)]^{-1}NXL} \\
&= \norm{L^*XN[I + (P + X)N]^{-1}XL}\\
&\leq \norm{L}^2\norm{X}^2\norm{N}.
\end{align}

With the same notations and assumptions of Proposition <ref>, for $X\in \cS_ν$, we have that

\begin{equation}
\norm{F(P+X, \Atilde, \Btilde) - F(P+X,A,B)}\leq \norm{A}^2\norm{P+X}^2\norm{\Delta N}+2\norm{A}\norm{P+X}\epsilon + \norm{P+X}\epsilon^2.\nonumber
\end{equation}

Let $P_X$ be a shorthand for $P + X$.
We can consider the following decomposition [32]:

\begin{align}
F(P_X, \Atilde, \Btilde) - F(P_X,A,B) =& A^*\PX(I + N\PX)^{-1}A- \Atilde^*\PX(I + \tilde N\PX)^{-1}\Atilde\\
=& A^*\PX(I + N\PX)^{-1}A -  A^*\PX(I + \Ntilde\PX)^{-1}A\nonumber\\
&+  A^*\PX(I + \Ntilde\PX)^{-1}A -  A^*\PX(I + \Ntilde\PX)^{-1}\Atilde\nonumber\\
&+  A^*\PX(I + \Ntilde\PX)^{-1}A -  \Atilde^*\PX(I + \Ntilde\PX)^{-1}A\nonumber\\
&+  A^*\PX(I + \Ntilde\PX)^{-1}\Atilde  +\Atilde^*\PX(I + \Ntilde\PX)^{-1}A\nonumber \\
&-A^*\PX(I + \Ntilde\PX)^{-1}A  -\Atilde^*\PX(I + \Ntilde\PX)^{-1}\Atilde\\
=& A^*\PX\left[(I + N\PX)^{-1} -  (I + \Ntilde\PX)^{-1}\right]A\nonumber\\
&-  A^*\PX(I + \Ntilde\PX)^{-1}\Delta A\nonumber\\
&-  \Delta A^*\PX(I + \Ntilde\PX)^{-1}A\nonumber\\
&-  \Delta A^*\PX(I + \Ntilde\PX)^{-1}\Delta A\\
=&A^*\PX(I + N\PX)^{-1} \Delta N \PX  (I + \Ntilde\PX)^{-1}A\nonumber\\
&-  A^*\PX(I + \Ntilde\PX)^{-1}\Delta A\nonumber\\
&-  \Delta A^*\PX(I + \Ntilde\PX)^{-1}A\nonumber\\
&-  \Delta A^*\PX(I + \Ntilde\PX)^{-1}\Delta A.
\end{align}

By <cit.>, we get

\begin{align}
\norm{F(P_X, \Atilde, \Btilde) - F(P_X,A,B)}&\leq \norm{A}^2\norm{\PX}^2\norm{\Delta N}+2\norm{A}\norm{\PX}\epsilon + \norm{\PX}\epsilon^2.
\end{align}

With the technical results introduced above, we are now ready to show that $\Phi$ maps $\mathcal S_{\nu}$ to itself, and is a contraction. We can firstly compute suitable upper bounds, merging the results above.

With the same notations and assumptions of Proposition <ref>, for $X\in\cS_ν$, $\nu\leq 1/2$, $\epsilon\leq \norm{B}$, and $\rho(L) \leq\rho < 1$, it holds that

\begin{equation}
\norm{\Phi X} \leq  \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}\nu^2 +3\epsilon(\norm{P} + 1)^2(\norm{A} + 1)^2(\norm{R^{-1}} + 1)(\norm{B} + 1) \right].\nonumber
\end{equation}

We already know from Proposition <ref> that $\norm{\Delta N}\leq 3\epsilon \norm{R^{-1}}\norm{B}$. Moreover, $\norm{\PX}\leq\norm{P} + \nu ≤ \n{P}+1$ by the triangle inequality and given that $X∊\cS_ν$ and $ν≤1/2$. According to Proposition <ref>, <ref>, <ref> and using $\epsilon\leq \norm{B}$:

\begin{align}
\norm{\Phi X} &= \norm{\mathcal T^{-1} [F(X + P, A, B) - F(X + P, \Atilde, \Btilde) - \mathcal R X]}\\
&\leq \norm{\mathcal T^{-1}}\left[\norm{F(X + P, A, B) - F(X + P, \Atilde, \Btilde)} + \norm{ \mathcal R X}\right]\\
&\leq \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}\nu^2 +\norm{A}^2\norm{\PX}^23\epsilon \norm{R^{-1}}\norm{B}+2\norm{A}\norm{\PX}\epsilon + \norm{\PX}\epsilon^2 \right]\\
&\leq  \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}\nu^2 +\epsilon\norm{\PX}\left(3\norm{A}^2\norm{\PX} \norm{R^{-1}}\norm{B}+2\norm{A} +\epsilon\right) \right]\\
&\leq  \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}\nu^2 +\epsilon\norm{\PX}\left(3\norm{A}^2\norm{\PX} \norm{R^{-1}}\norm{B}+2\norm{A} +\norm{B}\right) \right].%\\
%&\leq  \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}\nu^2 +3\epsilon(\norm{P} + 1)^2(\norm{A} + 1)^2(\norm{R^{-1}} + 1)(\norm{B} + 1) \right].
\end{align}

We use a rough upper-bounded on the second term of this last expression to get claimed result, which is more convenient to work with:

\begin{align}
3\epsilon(\norm{P} + 1)&^2(\norm{A} + 1)^2(\norm{R^{-1}} + 1)(\norm{B} + 1)\\
& =3\epsilon(\norm{P} + 1)^2(\norm{A} + 1)^2(\norm{R^{-1}}\norm{B} +\norm{R^{-1}} + \norm{B} + 1)\\
&=3\epsilon(\norm{P} + 1)^2(\norm{A}^2 + 2\norm{A} + 1)(\norm{R^{-1}}\norm{B} +\norm{R^{-1}} + \norm{B} + 1)\\
&=3\epsilon(\norm{P}^2 + 2\norm{P} + 1)(\norm{A}^2 + 2\norm{A} + 1)(\norm{R^{-1}}\norm{B} +\norm{R^{-1}} + \norm{B} + 1)\\
&\geq 3\epsilon\norm{A}^2(\norm{P}+1)^2\norm{R^{-1}}\norm{B} + 12\norm{A}(\norm{P} + 1)\epsilon + 6(\norm{P} + 1)\norm{B}\epsilon)\\
&\geq \norm{A}^2\norm{\PX}^23\epsilon \norm{R^{-1}}\norm{B}+2\norm{A}\norm{\PX}\epsilon + \norm{\PX}\epsilon^2.
\end{align}

With the same assumptions and notations of Proposition <ref>, if $X_1,X_2\in\mathcal S_{\nu}$, $\nu\leq \min\{1/2, \norm{N}^{-1}\}$, $‖ΔA‖≤ε$, $‖ΔB‖≤ε$, $\epsilon \leq \norm{B}$, and $\rho(L) \leq\rho < 1$, it holds that

\begin{align}
\norm{\Phi X_1 - \Phi X_2} \leq 32 \frac{\tau(L, \rho)^2}{1 - \rho^2}[&(\norm{A} + 1)^2(\norm{P} + 1)^3( \norm{B} + 1)^3(\norm R^{-1} + 1)^2 \norm{X_1 - X_2}\epsilon\nonumber \\
&+ \norm{L}^2\nu\norm{N}\norm{X_1 - X_2}].\nonumber
\end{align}

Let us choose $\nu\leq\min\{1/2, \norm{N}^{-1}\}$, and observe that:

\begin{align}
\norm{\cR X_1 - \cR X_2} \leq& \norm{L^*X_1 N[I + (P + X_1)N]^{-1}X_1L - L^*X_2 N[I + (P + X_2)N]^{-1}X_2L}\\
\leq& \norm{L}^2\norm{X_1 N[I + (P + X_1)N]^{-1}X_1 - X_1 N[I + (P + X_1)N]^{-1}X_2\nonumber\\
&+X_1 N[I + (P + X_1)N]^{-1}X_2-X_2 N[I + (P + X_1)N]^{-1}X_2\nonumber\\
&+ X_2 N[I + (P + X_1)N]^{-1}X_2
- X_2 N[I + (P + X_2)N]^{-1}X_2}\\
\leq& \norm{L}^22\nu\norm{N}\norm{X_1 - X_2} \nonumber\\
&+ \norm{X_2 N[I + (P + X_1)N]^{-1}\left(X_2N - X_1N\right)[I + (P + X_2)N]^{-1}X_2}\\
=&  \norm{L}^22\nu\norm{N}\norm{X_1 - X_2} \nonumber\\
&+ \norm{X_2 N[I + (P + X_1)N]^{-1}\left(X_2 - X_1\right)N[I + (P + X_2)N]^{-1}X_2}\\
\leq &\norm{L}^2\left\{ 2\nu\norm{N} +\nu^2\norm{N}^2\right\}\norm{X_1 - X_2}\\
\leq & 3\norm{L}^2\nu\norm{N}\norm{X_1 - X_2}.
\end{align}

Again, let $\PX$ be a shorthand for $ P + X$. Having defined for convenience $\cG X = F(\PX, \Atilde, \Btilde) - F(\PX, A, B)$,
and using the definition of $\Phi$ in (<ref>), it holds

\begin{align}
\norm{\Phi X_1 - \Phi X_2} =& \norm{\cTinv}\norm{\cG X_1 - \cG X_2 + \cR X_2 - \cR X_1}\\
\leq&  \frac{\tau(L, \rho)^2}{1 - \rho^2}[\norm{\cG X_1 - \cG X_2} + 3\norm{L}^2\nu\norm{N}\norm{X_1 - X_2}].
\end{align}

Now, we can study $\cG$. Note that, as in [32] $\norm{(I + N\PX)^{-1}}\leq 2\norm{\PX}$, since $\nu\leq \frac12$. Then,

\begin{align}
\norm{\cG X_1 - \cG X_2}=&\norm{A^*{\PX}_1(I + N{\PX}_1)^{-1}A- \Atilde^*{\PX}_1(I + \tilde N{\PX}_1)^{-1}\Atilde \nonumber\\
&- A^*{\PX}_2(I + N{\PX}_2)^{-1}A+ \Atilde^*{\PX}_2(I + \tilde N{\PX}_2)^{-1}\Atilde}\\
=&\lVert A^*{\PX}_1(I + N{\PX}_1)^{-1} \Delta N {\PX}_1  (I + \Ntilde{\PX}_1)^{-1}A- A^*{\PX}_2(I + N{\PX}_2)^{-1} \Delta N {\PX}_2 (I + \Ntilde{\PX}_2)^{-1}A\label{e:cerulean}\\
&-  A^*{\PX}_1(I + \Ntilde{\PX}_1)^{-1}\Delta A+  A^*{\PX}_2(I + \Ntilde{\PX}_2)^{-1}\Delta A\label{e:orange}\\
&-  \Delta A^*{\PX}_1(I + \Ntilde{\PX}_1)^{-1}A+  \Delta A^*{\PX}_2(I + \Ntilde{\PX}_2)^{-1}A\label{e:purple}\\
&-  \Delta A^*{\PX}_1(I + \Ntilde{\PX}_1)^{-1}\Delta A+  \Delta A^*{\PX}_2(I + \Ntilde{\PX}_2)^{-1}\Delta A\rVert\label{e:sepia}.
\end{align}

To control the difference in (<ref>), we can use the following inequality. Note that, by definition, $\norm N\leq \norm B^2 \norm {R^{-1}}$. Then, leveraging <cit.>:

\begin{align}
\norm{A^*{\PX}_1(I + &N{\PX}_1)^{-1} \Delta N {\PX}_1  (I + \Ntilde{\PX}_1)^{-1}A - A^*{\PX}_2(I + N{\PX}_2)^{-1} \Delta N {\PX}_2 (I + \Ntilde{\PX}_2)^{-1}A}\nonumber\\
=& \norm{A^*{\PX}_1(I + N{\PX}_1)^{-1} \Delta N {\PX}_1  (I + \Ntilde{\PX}_1)^{-1}A - A^*{\PX}_1(I + N{\PX}_2)^{-1} \Delta N {\PX}_1(I + \Ntilde{\PX}_1)^{-1}A\nonumber\\
&+ A^*{\PX}_1(I + N{\PX}_2)^{-1} \Delta N {\PX}_1  (I + \Ntilde{\PX}_1)^{-1}A - A^*{\PX}_2(I + N{\PX}_2)^{-1} \Delta N {\PX}_1 (I + \Ntilde{\PX}_1)^{-1}A\nonumber\\
&+A^*{\PX}_2(I + N{\PX}_2)^{-1} \Delta N {\PX}_1  (I + \Ntilde{\PX}_1)^{-1}A - A^*{\PX}_2(I + N{\PX}_2)^{-1} \Delta N {\PX}_2 (I + \Ntilde{\PX}_1)^{-1}A\nonumber\\
&+A^*{\PX}_2(I + N{\PX}_2)^{-1} \Delta N {\PX}_2  (I + \Ntilde{\PX}_1)^{-1}A - A^*{\PX}_2(I + N{\PX}_2)^{-1} \Delta N {\PX}_2 (I + \Ntilde{\PX}_2)^{-1}A}\\
\leq& \norm{A}^2(\norm{P} + 1)\norm{\Delta N}\norm{P_{X_1}(I +N{\PX}_1 )^{-1} - P_{X_1}(I +N{\PX}_2 )^{-1}}\nonumber\\
&+ \norm{A}^2\norm{\Delta N}(\norm P + 1)^3 \norm{X_1 - X_2}\nonumber\\
&+ \norm{A}^2\norm{\Delta N}(\norm P + 1)^3 \norm{X_1 - X_2}\nonumber\\
&+ \norm{A}^2(\norm{P} + 1)\norm{\Delta N}\norm{P_{X_2}(I +N{\PX}_1 )^{-1} - P_{X_2}(I +N{\PX}_2 )^{-1}}\\
\leq& \norm{A}^2(\norm{P} + 1)\norm{\Delta N}\norm{P_{X_1}(I +N{\PX}_1 )^{-1} (I + NP_{X_2} - I - NP_{X_1})(I +N{\PX}_2 )^{-1}}\nonumber\\
&+ \norm{A}^2\norm{\Delta N}(\norm P + 1)^3 \norm{X_1 - X_2}\nonumber\\
&+ \norm{A}^2\norm{\Delta N}(\norm P + 1)^3 \norm{X_1 - X_2}\nonumber\\
&+ \norm{A}^2(\norm{P} + 1)\norm{\Delta N}\norm{P_{X_2}(I + NP_{X_2})^{-1}(I + NP_{X_2} - I - NP_{X_1})(I +N{\PX}_1 )^{-1} }\\
\leq& \norm{A}^22(\norm{P} + 1)^3\norm{\Delta N}\norm{N}\norm{X_2-X_1}\nonumber\\
&+ \norm{A}^2\norm{\Delta N}(\norm P + 1)^3 \norm{X_1 - X_2}\nonumber\\
&+ \norm{A}^2\norm{\Delta N}(\norm P + 1)^3 \norm{X_1 - X_2}\nonumber\\
&+ \norm{A}^22(\norm{P} + 1)^3\norm{\Delta N}\norm{N}\norm{X_2 -X_1}\\
\leq& 12\norm{A}^2(\norm{P} + 1)^3\norm{B}^3\norm{R^{-1}}^2\norm{X_2-X_1}\epsilon\nonumber\\
&+ 6\norm{A}^2\norm{B}\norm{R^{-1}}(\norm P + 1)^3 \norm{X_1 - X_2}\epsilon.
\end{align}

With a similar reasoning, we can compute an upper bound for the difference in (<ref>), as follows:

\begin{align}
\norm{A^*{\PX}_2(I +& \Ntilde{\PX}_2)^{-1}\Delta A - A^*{\PX}_1(I + \Ntilde{\PX}_1)^{-1}\Delta A }\nonumber\\
=&\norm{A^*{\PX}_2(I + \Ntilde{\PX}_2)^{-1}\Delta A - A^*{\PX}_2(I + \Ntilde{\PX}_1)^{-1}\Delta A\nonumber \\
&+A^*{\PX}_2(I + \Ntilde{\PX}_1)^{-1}\Delta A - A^*{\PX}_1(I + \Ntilde{\PX}_1)^{-1}\Delta A }\\
\leq&4\norm{A}(\norm{P} + 1)^3\norm{X_2 - X_1}\norm{B}^2\norm{R^{-1}}\epsilon\nonumber\\
&+ 2\norm{A}\norm{X_2 - X_1}(\norm P  + 1)^2\epsilon.
\end{align}

The difference in (<ref>) is the adjoint of the difference in (<ref>), and can be upper bounded by the same factor. An upper bound on the difference in (<ref>) can be computed as follows:

\begin{align}
\norm{\Delta A^*{\PX}_2(I +& \Ntilde{\PX}_2)^{-1}\Delta A - \Delta A^*{\PX}_1(I + \Ntilde{\PX}_1)^{-1}\Delta A }\nonumber\\
=&\norm{\Delta A^*{\PX}_2(I + \Ntilde{\PX}_2)^{-1}\Delta A - \Delta A^*{\PX}_2(I + \Ntilde{\PX}_1)^{-1}\Delta A \nonumber\\
&+\Delta A^*{\PX}_2(I + \Ntilde{\PX}_1)^{-1}\Delta A - \Delta A^*{\PX}_1(I + \Ntilde{\PX}_1)^{-1}\Delta A }\\
\leq&4(\norm{P} + 1)^3\norm{X_2 - X_1}\norm{B}^2\norm{R^{-1}}\epsilon^2\nonumber\\
&+ 2\norm{X_2 - X_1}(\norm P  + 1)^2\epsilon^2\\
\leq & 4(\norm{P} + 1)^3\norm{X_2 - X_1}\norm{B}^3\norm{R^{-1}}\epsilon\nonumber\\
&+ 2\norm{X_2 - X_1}(\norm P  + 1)^2\norm{B}\epsilon,
\end{align}

where we used on the last line the assumption $ε≤‖B‖$.
The overall upper bound has 6 terms in which $\norm A$ appears at most with degree 2, $\norm P + 1$ with degree 3, $\norm B$ with degree 3, $\norm{R^{-1}}$ with degree 2. With the same reasoning as for the other upper bound, we can write

\begin{equation}
\norm{ \cG X_1 - \cG X_2} \leq 12 (\norm{A} + 1)^2(\norm{P} + 1)^3( \norm{B} + 1)^3(\norm R^{-1} + 1)^2 \norm{X_1 - X_2}\epsilon.
\end{equation}

To conclude,

\begin{align}
\norm{\Phi X_1 - \Phi X_2} \leq 12 \frac{\tau(L, \rho)^2}{1 - \rho^2}[&(\norm{A} + 1)^2(\norm{P} + 1)^3( \norm{B} + 1)^3(\norm R^{-1} + 1)^2 \norm{X_1 - X_2}\epsilon\nonumber \\
&+ \norm{L}^2\nu\norm{N}\norm{X_1 - X_2}].
\end{align}

We are now ready to show that $\Phi$ is a contraction from $\mathcal S_{\nu}$ to itself. First, we show that it maps $\cS_{\nu}$ to $\cS_{\nu}$ (Lemma <ref>). Then, we show that it is a contraction on $\cS_{\nu}$ (Lemma <ref>).

Let us choose $\nu$ as a function of $\epsilon$, namely,

\begin{equation}
\nu= 	\min\left\{6\epsilon \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp A^2 \normp P^2 \normp B \normp {R^{-1}}, \norm{N}^{-1}, \frac12\right\}.\nonumber
\end{equation}

Moreover, let $\epsilon$ be small enough, namely,

\begin{equation}
\epsilon \leq\min\left\{ \frac{1}{12}\normp{L}^{-2} \frac{(1 - \rho^2)^2}{\tau(L, \rho)^4}\normp A^{-2} \normp P^{-2} \normp B^{-3} \normp {R^{-1}}^{-2}, \norm B\right\}.\nonumber
\end{equation}

Lastly, let $\sigma_{\textrm{min}}(P)\geq 1$, and $\rho(L) \leq\rho < 1$. Note that the condition on the singular values of $P$ can be achieved by rescaling $R$ and $Q$ accordingly, as discussed by [32].
Then, for $X\in\mathcal S_{\nu}$, we have that

\begin{equation}
\Phi X\in \mathcal {S}_{\nu}.\nonumber
\end{equation}

In order to prove the lemma, we need to first show that $\Phi X$ has norm upper bounded by $\nu$, for $X\in\cS_ν$ and a suitable choice of $\nu$. Then, we need to show that $P + \Phi X$ is self-adjoint and positive semi-definite. Starting from Proposition <ref>, we see that:

\begin{align}
\norm{\Phi X} &\leq \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}\nu^2 +3\epsilon\normp{P}^2\normp{A}^2\normp{R^{-1}}\normp{B}  \right]\\
&\leq	\frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}	\left(6\epsilon \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp A^2 \normp P^2 \normp B \normp {R^{-1}}\right)^2\right.\nonumber\\
&\quad\left.+3\epsilon\normp{P}^2\normp{A}^2\normp{R^{-1}}\normp{B}  \right]\\
&= \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^2\norm{N}12\epsilon \frac{\tau(L, \rho)^4}{(1 - \rho^2)^2}\normp A^2 \normp P^2 \normp B \normp {R^{-1}} + 1 \right]\nonumber\\
&\quad\cdot3\epsilon\normp{P}^2\normp{A}^2\normp{R^{-1}}\normp{B} \\
&\leq  \frac{ \tau(L, \rho)^2}{1-\rho^2}\left[\norm{L}^212\epsilon \frac{\tau(L, \rho)^4}{(1 - \rho^2)^2}\normp A^2 \normp P^2 \normp B^3 \normp {R^{-1}}^2 + 1 \right]\nonumber\\
&\quad\cdot3\epsilon\normp{P}^2\normp{A}^2\normp{R^{-1}}\normp{B} .
\end{align}

For the choice of $\epsilon$ in the statement of this lemma, the bracket term is bounded by 2, meaning that:

\begin{align}
\norm{\Phi X}
\leq  \frac{ \tau(L, \rho)^2}{1-\rho^2}6\epsilon\normp{P}^2\normp{A}^2\normp{R^{-1}}\normp{B} .
\end{align}

Note that, for our choice of $\epsilon$,

\begin{align}
\norm{\Phi X}
%\frac{ \tau(L, \rho)^2}{1-\rho^2}6\epsilon&\normp{P}^2\normp{A}^2\normp{R^{-1}}\normp{B} \\
&\leq \frac{1}{2}\frac{1-\rho^2}{ \tau(L, \rho)^2}\normp{R^{-1}}^{-1}\normp{B}^{-2}\normp{L}^{-2}\\
&\leq\min\left\{\frac 12, \norm{N}^{-1}\right\},
\end{align}

since $‖N‖=‖BR^{-1}B^*‖≤‖B‖²‖R^{-1}‖$ and $\frac{1-\rho^2}{ \tau(L, \rho)^2}\leq 1$. This can be noted by observing that $\tau(L, \rho)$ is a decreasing function of $\rho$ and is equal to 1 for $\rho\geq\norm{L}$ (see <cit.>), and $1 - \rho^2 < 1$ by our choice of $\rho$.
All in  all, this means that $\norm{\Phi X}\leq \nu$.
To conclude the proof, we can observe that, since $P$ is positive semi-definite and $\sigma_{\textit{min}}(P)\geq 1$ by assumption, for $a\in \cH_1$,

\begin{align}
\innerprod{(P + \Phi X)a, a}_{\cH_1}&= \innerprod{Pa, a}_{\cH_1} + \innerprod{\Phi Xa, a}_{\cH_1}\\
&\geq \innerprod{Pa, a}_{\cH_1} - \norm{\Phi X}\norm{a}^2\\
&\geq \innerprod{Pa, a}_{\cH_1} -\nu\norm{a}^2\\
&= \norm{P^{1/2}a}^2 -\nu\norm{a}^2\\
&\geq [\sigma_{\textit{min}}(P) - \nu]\norm a^2\\
&\geq 0.
\end{align}

Hence, $P + \Phi X\succcurlyeq 0$.
Moreover, $ΦX$ is self-adoint because $\cT,\cR$ and $X\mapsto F(P+X,A,B)$ preserve self-adjointness.

Let $\nu $ be defined as in Lemma <ref>, and $\rho(L) \leq\rho < 1$. Let

\begin{equation}
12\frac{\tau(L, \rho)^4}{(1 - \rho^2)^2}\left[\normp{P}+\normp{L}^2\right]\normp A^2 \normp P^2\normp B^3\normp{R^{-1}}^2\epsilon< 1.\nonumber
\end{equation}

Then, $\exists\eta < 1$ such that, $\forall X_1, X_2 \in \cS_{\nu}$,

\begin{equation}
\norm{\Phi X_1 - \Phi X_2}\leq\eta\norm{X_1 - X_2}.\nonumber
\end{equation}

Let us plug our choice of $\nu$ from Lemma <ref> in the bound of Proposition <ref>, to get

\begin{align}
12\frac{\tau(L, \rho)^2}{1 - \rho^2}&\left[ \normp{A}^2\normp{P}^3\normp{B}^3\normp{ R^{-1}}^2 \norm{X_1 - X_2}\epsilon+ \norm{L}^2\nu\norm{N}\norm{X_1 - X_2}\right]\nonumber\\
&=	 12\frac{\tau(L, \rho)^2}{1 - \rho^2}\left[\normp{P}\normp{B}^2\normp{ R^{-1}} + \frac{\tau(L, \rho)^2}{1 - \rho^2}\norm{L}^2\norm{N}\right]\nonumber\\
&\quad\cdot\normp A^2 \normp P^2\normp B\normp{R^{-1}}\epsilon\norm{X_1 - X_2}\\
&\leq 12\frac{\tau(L, \rho)^2}{1 - \rho^2}\left[\normp{P}\normp{B}^2\normp{ R^{-1}} + \frac{\tau(L, \rho)^2}{1 - \rho^2}\norm{L}^2\normp B^2 \normp{R^{-1}}\right]\nonumber\\
&\quad\cdot\normp A^2 \normp P^2\normp B\normp{R^{-1}}\epsilon\norm{X_1 - X_2}\\
&=12 \frac{\tau(L, \rho)^2}{1 - \rho^2}\left[\normp{P}+ \frac{\tau(L, \rho)^2}{1 - \rho^2}\norm{L}^2\right]\nonumber\\
&\quad\cdot\normp A^2 \normp P^2\normp B^3\normp{R^{-1}}^2\epsilon\norm{X_1 - X_2}\\
&\leq 12 \frac{\tau(L, \rho)^2}{1 - \rho^2}\left[ \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp{P}+ \frac{\tau(L, \rho)^2}{1 - \rho^2}\norm{L}^2\right]\nonumber\\
&\quad\cdot\normp A^2 \normp P^2\normp B^3\normp{R^{-1}}^2\epsilon\norm{X_1 - X_2}\\
&\leq 12 \frac{\tau(L, \rho)^2}{1 - \rho^2}\left[ \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp{P}+ \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp{L}^2\right]\nonumber\\
&\quad\cdot\normp A^2 \normp P^2\normp B^3\normp{R^{-1}}^2\epsilon\norm{X_1 - X_2}.
\end{align}

Again, in the last step we used the fact that $\tau(L, \rho)$ is a decreasing function of $\rho$ and is equal to 1 for $\rho\geq\norm{L}$ (see <cit.>), and $1 - \rho^2 < 1$ by our choice of $\rho$.
The operator $\Phi$ is a contraction if  $\epsilon$ is small enough, namely,

\begin{equation}
12\frac{\tau(L, \rho)^4}{(1 - \rho^2)^2}\left[\normp{P}+\normp{L}^2\right]\normp A^2 \normp P^2\normp B^3\normp{R^{-1}}^2\epsilon< 1.
\end{equation}

Lemmas <ref>, <ref> imply that, for $\epsilon$ small enough, the unique fixed point of $\Phi$ is in $\mathcal S_{\nu}$. According to Proposition <ref>, that fixed point is exactly the error on the Riccati operator due to the Nyström approximation, whose norm scales linearly in $\epsilon$. To conclude the proof of Lemma <ref>, we can observe that the conditions stated in Lemmas <ref>, <ref> yield an upper bound on $\epsilon$:

\begin{equation}
\epsilon< \frac{1}{12} \frac{1}{\normp{L}^{2} +\normp{P}}\frac{(1 - \rho^2)^2}{\tau(L, \rho)^4}\normp A^{-2} \normp P^{-2} \normp B^{-3} \normp {R^{-1}}^{-2} ,
\end{equation}

\begin{equation}
\norm{P-\tilde P}\leq 6\epsilon \frac{\tau(L, \rho)^2}{1 - \rho^2}\normp A^2 \normp P^2 \normp B \normp {R^{-1}}.
\end{equation}

§ PROOF OF THEOREM <REF>

The proof follows the steps in <cit.>, adapted to the operator setting considered in this work. We begin by computing an upper bound on the error for the Riccati gain $\norm{K-\Ktilde}$. We then show that when this error is small enough, $\Ktilde$ stabilizes the system with exact kernel defined in (<ref>). Lastly, we use <cit.> to upper bound the error on the objective functions.
We begin by re-stating two technical lemmas that appear in <cit.> generalized to the case in which the LQR objective contains operators.

Let $f_1$, $f_2$ be two $\mu$-strongly convex twice differentiable functions on a Hilbert space. Let $\xvec_1=\argmin_{\xvec} f_1(\xvec)$ and $\xvec_2 = \argmin_{\xvec}=f_2(\xvec)$. If $\norm{\frac{\partial f_1(\xvec)}{\partial \xvec}\rvert_{\xvec=\xvec_2}}\leq \epsilon$, then $\norm{\xvec_1 - \xvec_2}\leq \frac\epsilon\mu$.

Let $A_i:\cH_1\to\cH_1$, $B_i:\R^{n_u}\to\cH_1$, $i=1,2$.
Let us define $f_i:\R^{n_u}\times\cH_1\rightarrow \R$, $f_i(\uvec, z) = \frac12\innerprod{\uvec, R\uvec}_{\R^{n_u}} + \frac12\innerprod{A_iz + B_i\uvec, P_i(A_iz + B_i\uvec)}_{\cH_1}$ where $R$, $P_i$ are positive definite operators, $i=1, 2$, $z\in \cH_1$, $\uvec\in\R^{n_u}$. Let $K_i:\cH_1\to\R^{n_u}$ be the unique operator s.t. $\uvec_i = \argmin_{\uvec}f_i(\uvec, z)=K_i z$ for any $z$. Let us define the quantity $\Gamma = 1 +\max\{\norm {A_1}, \norm {B_1}, \norm {P_1}, \norm {K_1}\}$, and $\norm{A_1 - A_2}\leq \epsilon$, $\norm{B_1 - B_2}\leq \epsilon$, and $\norm{P_1 - P_2}\leq \epsilon$, for $\epsilon\in[0, 1)$. Then,

\begin{equation}
\norm{K_1 - K_2} \leq \frac{3\epsilon\Gamma^3}{\sigma_{\textrm{min}}(R)}.
\end{equation}

Each $f_i$ in the statement of this lemma is the so-called state-action value function in the language of [17]. It corresponds to an LQR objective when picking an arbitrary control input $\uvec$ at time 0, and then proceeding with the optimal control policy (static state feedback gain) for the next time steps. It is known that, if at timestep 1 we start using the LQR optimal control, the LQR objective from that timestep onwards is equal to $\frac{1}{2}\innerprod{z_1, Pz_1}_{\cH_1}$.

The gains $K_1$ (resp. $K_2$) in the statement of this lemma are the Riccati gains, obtained with the dynamics given by $A_1$, $B_1$ (resp. $A_2$, $B_2$).

The structure of the proof is as follows. We begin by showing that our choice of $f_1, f_2$ fulfills the hypotheses of <ref>. Then we apply such a lemma to upper bound the error of interest.
Computing the gradient of $f_i$ w.r.t. $\uvec$ yields:

\begin{align}
\frac{\partial f_i (\uvec, z)}{\partial \uvec} &= R\uvec + \frac{\partial}{\partial \uvec}\left\{\frac12\innerprod{A_iz + B_i\uvec, P_i(A_iz + B_i\uvec)}_{\cH_1}\right\} \\
&= R\uvec +  \frac{\partial}{\partial \uvec}\left\{\frac12\innerprod{A_iz, P_iA_iz }_{\cH_1} + \innerprod{A_iz, P_iB_i\uvec}_{\cH_1} +  \frac12\innerprod{B_i\uvec, P_iB_i\uvec}_{\cH_1}\right\}\\
&= R\uvec + B_i^*P_iA_iz + B_i^*P_i^{1/2}P_i^{1/2}B_i\uvec\\
&= B_i^*P_iA_iz + (B_i^*P_iB_i + R)\uvec.\label{eq:derivative_objective}
\end{align}

After having derived the gradient expression above, we can bound two differences that will appear in the remainder of the proof. In the same way as in [32], we can observe that, having defined $\Lambda=\max\{\norm{P_1}, \norm{B_1}, \norm{A_1}\}$, and considering the assumption that $\epsilon\in[0, 1)$:

\begin{align}
\norm{B_1^*P_1B_1 - B_2^*P_2B_2} =& \norm{(B_1^* - B_2^*)P_1B_1 + B_2^*(P_1 - P_2)B_1 + B_2^*P_2(B_1 - B_2)}\\
=&\norm{(B_1^* - B_2^*)P_1B_1 + B_2^*(P_1 - P_2)B_1 + B_2^*P_2(B_1 - B_2)\nonumber\\
&-B_1^*(P_1 - P_2)B_1 + B_1^*(P_1 - P_2)B_1\nonumber\\
&-B_1^*P_2(B_1 - B_2) + B_1^*P_2(B_1 - B_2)}\\
=&\norm{(B_1^* - B_2^*)P_1B_1 + (B_2^* - B_1^*)(P_1 - P_2)B_1 + (B_2^* - B_1^*)P_2(B_1 - B_2)\nonumber\\
&+ B_1^*(P_1 - P_2)B_1+ B_1^*P_2(B_1 - B_2)\nonumber\\
&- (B_2^* - B_1^*)P_1(B_1 - B_2) + (B_2 ^*- B_1^*)P_1(B_1 - B_2)\nonumber\\
&-B_1^*P_1(B_1 - B_2) + B_1^*P_1(B_1 - B_2)}\\
=& \norm{(B_1^* - B_2^*)P_1B_1 + (B_2^* - B_1^*)(P_1 - P_2)B_1 \nonumber\\
&+ (B_2^* - B_1^*)(P_2 - P_1)(B_1 - B_2)\nonumber\\
&+ B_1^*(P_1 - P_2)B_1+ B_1^*(P_2 - P_1)(B_1 - B_2)\nonumber\\
&+ (B_2^* - B_1^*)P_1(B_1 - B_2)\nonumber\\
&+ B_1^*P_1(B_1 - B_2)}\\
\leq& \epsilon \norm{P_1}\norm{B_1} + \epsilon^2\norm{B_1} + \epsilon^3+ \norm{B_1}^2\epsilon+ \norm{B_1}\epsilon^2+ \epsilon^2\norm{P_1}+ \norm{B_1}\norm{P_1}\epsilon\\
\leq& \epsilon\Lambda^2 + \epsilon^2\Lambda + \epsilon^3 + \epsilon\Lambda^2 + \epsilon\Lambda + \epsilon^2 \Lambda + \Lambda^2\epsilon\\
\leq& \epsilon (3\Lambda^2 + 3\Lambda + 1)\\
\leq& 3\epsilon\Gamma^2.
\end{align}

Similarly, we have

\begin{align}
\norm{B_1^*P_1A_1 - B_2^*P_2A_2} =& \norm{(B_1^* - B_2^*)P_1A_1 + B_2^*(P_1 - P_2)A_1 + B_2^*P_2(A_1 - A_2)}\\
=&\norm{(B_1^* - B_2^*)P_1A_1 + B_2^*(P_1 - P_2)A_1 + B_2^*P_2(A_1 - A_2)\nonumber\\
&-B_1^*(P_1 - P_2)A_1 + B_1^*(P_1 - P_2)A_1\nonumber\\
&-B_1^*P_2(A_1 - A_2) + B_1^*P_2(A_1 - A_2)}\\
=&\norm{(B_1^* - B_2^*)P_1A_1 + (B_2^* - B_1^*)(P_1 - P_2)A_1 + (B_2^* - B_1^*)P_2(A_1 - A_2)\nonumber\\
&+ B_1^*(P_1 - P_2)A_1+ B_1^*P_2(A_1 - A_2)\nonumber\\
&- (B_2^* - B_1^*)P_1(A_1 - A_2) + (B_2 ^*- B_1^*)P_1(A_1 - A_2)\nonumber\\
&-B_1^*P_1(A_1 - A_2) + B_1^*P_1(A_1 - A_2)}\\
=& \norm{(B_1^* - B_2^*)P_1A_1 + (B_2^* - B_1^*)(P_1 - P_2)A_1 \nonumber\\
&+ (B_2^* - B_1^*)(P_2 - P_1)(A_1 - A_2)\nonumber\\
&+ B_1^*(P_1 - P_2)A_1+ B_1^*(P_2 - P_1)(A_1 - A_2)\nonumber\\
&+ (B_2^* - B_1^*)P_1(A_1 - A_2)\nonumber\\
&+ B_1^*P_1(A_1 - A_2)}\\
\leq& \epsilon \norm{P_1}\norm{A_1} + \epsilon^2\norm{A_1} + \epsilon^3\nonumber\\
&+ \norm{B_1}\norm{A_1}\epsilon+ \norm{B_1}\epsilon^2+ \epsilon^2\norm{P_1}+ \norm{B_1}\norm{P_1}\epsilon\\
\leq& \epsilon\Lambda^2 + \epsilon^2\Lambda + \epsilon^3 + \epsilon\Lambda^2 + \epsilon\Lambda + \epsilon^2 \Lambda + \Lambda^2\epsilon\\
\leq& \epsilon (3\Lambda^2 + 3\Lambda + 1)\\
\leq& 3\epsilon\Gamma^2.
\end{align}

We can then consider the following difference. For arbitrary $\uvec$, $z$, by using the gradient expression in (<ref>) and the triangle inequality,

\begin{equation}
\left\lVert	\frac{\partial f_1 (\uvec, z)}{\partial \uvec} - 	\frac{\partial f_2 (\uvec, z)}{\partial \uvec}\right\rVert\leq 3\epsilon\Gamma^2(\norm \uvec + \norm z).\label{eq:grad_diff}
\end{equation}

In particular, for any $z$ $‖z‖≤1$, this inequality applied to $\uvec_1=K_1 z$ yields by the first-order optimality condition of $\uvec_1$:

\begin{align}
\left\lVert	\frac{\partial f_1 (\uvec, z)}{\partial \uvec} \rvert_{\uvec=\uvec_1}- 	\frac{\partial f_2 (\uvec, z)}{\partial \uvec}\rvert_{\uvec=\uvec_1}\right\rVert
&=	\left\lVert	\frac{\partial f_2 (\uvec, z)}{\partial \uvec}\rvert_{\uvec=\uvec_1}\right\rVert %\\
\leq 3\epsilon\Gamma^2(\norm {\uvec_1} + 1).
\end{align}

Hence, still for any $z$ with $‖z‖≤1$, applying Lemma <ref> to the functions $f₁,f₂$ defined in the statement, which are $\mu$-strongly convex with $\mu\geq \sigma_{\textrm{min}}(R)$, we get

\begin{align}
\n{K₁z-K₂z}
&≤ \sigma_{\textrm{min}}(R)^{-1} 3εΓ^2(\norm{\uvec_1} +1) \\
&\leq \sigma_{\textrm{min}}(R)^{-1}3εΓ^2(\norm{K_1}\norm{z} +1)\\
&\leq \sigma_{\textrm{min}}(R)^{-1}3εΓ^2(\norm{K_1} +1)\\
&\leq \sigma_{\textrm{min}}(R)^{-1}3εΓ^3.
\end{align}

This yields the claimed result given that

\begin{align}
\norm{K_1 - K_2}
&= \sup_{\norm{z}\leq 1} \norm{(K_1 - K_2)z}. %\\
%&= \norm{(K_1 - K_2)z_\textup{sup}}.
\end{align}

Having proved the lemma above, we now specialize it to the case of exact vs. Nyström kernel-based gains.

Let $\epsilon>0$ be s.t. $\norm{ A - \Atilde}\leq \epsilon$, $\norm{ B - \Btilde}\leq \epsilon$, and $\norm{P - \Ptilde}\leq g(\epsilon)$ with $g(\epsilon)\geq \epsilon$. Then, assuming $R$ and $Q$ are positive definite and $\sigma_{\textrm{min}}(R)\geq1$, we have

\begin{equation}
\label{eq:bound_gain_error}
\norm{\Ktilde - K} \leq 3\Gamma^3g(\epsilon).
\end{equation}

Moreover, choose $\rho$ such that $\rho(L) \leq \rho<1$. If $\epsilon$ is small enough s.t. the r.h.s. of (<ref>) is smaller than $\frac{1 - \rho}{2\tau(L, \rho)}$, then

\begin{equation}
\tau\left(A + B\Ktilde, \frac{1 + \rho}{2}\right)\leq \tau(L, \rho).
\end{equation}

The first part of the proposition is a corollary of Lemma <ref> based on the condition $g(\epsilon)> \epsilon$ and $\sigma_{\textrm{min}}(R) \geq 1$. The second part of the proposition follows from applying the first inequality of <cit.> with $M=A+BK$ and $Δ=B(\tilde K- K)$, which guarantees that, if

\begin{equation}
\norm{\Ktilde - K}\leq \frac{1 - \rho}{2\tau(L, \rho)\norm B},
\end{equation}

\begin{align}
\norm{(A + B\Ktilde)^k}&= \norm{(A + BK + B\Ktilde - BK)^k}\\
&= \norm{(A + BK + B(\Ktilde - K))^k}\\
&\leq\tau(L, \rho)\left(\frac{1 - \rho}{2} + \rho\right)^k\\
&= \tau(L, \rho)\left(\frac{1 + \rho}{2}\right)^k.
\end{align}

Now we are ready to conclude the proof. Let us consider the following decomposition of the error on the LQR objective, which corresponds to the well known performance difference lemma in the reinforcement learning literature [24]. The following lemma is a restatement of <cit.>, with our notation and setup (i.e., with operators).

Let $\cJ$ and $\cJhat$ be as of (<ref>) and (<ref>). Then, the following decomposition holds:

\begin{equation}
\cJhat - \cJ = \lim_{T\to\infty}\sum_{i=0}^{T}\ \innerprod{(\Ktilde - K)z_i, (R + B^*PB)(\Ktilde - K)z_i}_{\cH_1}
\end{equation}

We begin by using the same telescoping argument applied by [17].
Let $A$ and $B$ be defined as in (<ref>), (<ref>). Let us define, for an arbitrary initial condition $\hat z_0\in\cH_1$:

\begin{equation}
\hat z_{i+1} = A\hat z_i + B\tuveco_i,\tuveco_i=\Ktilde\zhat_i, i\geq 0.
\end{equation}

Similarly to (<ref>), let us define, for an initial condition, at timestep $i$, $\hat z_i\in\cH_1$:

\begin{align}
\hat z_{i,j+1}' &= A\hat z_{i,j}' + B\uveco_j, \uveco_j = K \zhat_j',  j\geq i,\\
\hat z_{i,i}' &= \hat z_i.
\end{align}

We can decompose the error of interest as follows, noting that, by the definition in (<ref>),

\begin{equation}
\cJ=\sum_{j=0}^T \left(\innerprod{\hat z_{0, j}', Q\hat z_{0, j}'}_{\cH_1}+ \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right):
\end{equation}

\begin{align}
\cJhat - \cJ =& \lim_{T\to\infty}\left\{\sum_{i=0}^T \left[\innerprod{\hat z_i, Q\hat z_i}_{\cH_1} + \innerprod{\tuveco_i, R\tuveco_i}_{\R^{n_u}}\right] - \cJ\right\} \\
=& \lim_{T\to\infty}\left\{\sum_{i=0}^T \left[\innerprod{\hat z_i, Q\hat z_i}_{\cH_1} + \innerprod{\tuveco_i, R\tuveco_i}_{\R^{n_u}} \right.\right.\nonumber\\
& +\sum_{j=i}^T \left(\innerprod{\hat z_j', Q\hat z_j'}_{\cH_1} + \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right)\nonumber\\
&- \left.\left.\sum_{j=i}^T \left(\innerprod{\hat z_j', Q\hat z_j'}_{\cH_1}+ \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right)\right] - \cJ\right\}\\
=&\lim_{T\to\infty}\left\{\sum_{i=0}^T \left[\innerprod{\hat z_i, Q\hat z_i}_{\cH_1} + \innerprod{\tuveco_i, R\tuveco_i}_{\R^{n_u}} +\sum_{j=i+1}^T \left(\innerprod{\hat z_{i,j}', Q\hat z_{i,j}'}_{\cH_1} + \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right)\right.\right.\nonumber\\
&- \left.\left.\sum_{j=i}^T \left(\innerprod{\hat z_{i,j}', Q\hat z_{i,j}'}_{\cH_1}+ \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right)\right]\right\}.
\end{align}

Now we can consider a single addend in the outer-most sum.
Consider that by definition

\begin{equation}
\tuveco_i = \Ktilde \zhat_i.
\end{equation}

Moreover, note that

\begin{align}
\zhat_{i, i+1}' &= A \zhat_{i, i}' + B \Ktilde\zhat_{i, i}'\\
&=A \zhat_i + B \Ktilde\zhat_i.
\end{align}

Lastly, note that, subject to the Riccati-optimal state feedback law $\uveco$, the following simplification is allowed, for $T\to \infty$:

\begin{align}
\sum_{j=i+1}^T \left(\innerprod{\hat z_{i,j}', Q\hat z_{i,j}'}_{\cH_1}
+ \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right)&= \innerprod{\zhat_{i, i+1}' , P \zhat_{i, i+1}' }_{\cH_1} \\
&=  \innerprod{(A + B\Ktilde)\hat z_i, P (A + B\Ktilde) \hat z_i}_{\cH_1}.
\end{align}

Similarly,  observe that, for $T\to \infty,$

\begin{align}
\sum_{j=i}^T \left(\innerprod{\hat z_{i,j}', Q\hat z_{i,j}'}_{\cH_1}
+ \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right)&= \innerprod{\zhat_{i, i}' , P \zhat_{i, i}' }_{\cH_1} \\
&=  \innerprod{\hat z_i, P \hat z_i}_{\cH_1}.
\end{align}

Hence, we can rewrite the addend as the following, for $T\to\infty$:

\begin{align}
\innerprod{\hat z_i, Q\hat z_i}_{\cH_1}& + \innerprod{\tuveco_i, R\tuveco_i}_{\R^{n_u}}
+\sum_{j=i+1}^T \left(\innerprod{\hat z_{i,j}', Q\hat z_{i,j}'}_{\cH_1}
+ \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right) \nonumber\\
&-\sum_{j=i}^T \left(\innerprod{\hat z_{i,j}', Q\hat z_{i,j}'}_{\cH_1}+ \innerprod{\uveco_j, R\uveco_j}_{\R^{n_u}}\right)\nonumber\\
=& \innerprod{\hat z_i, (Q + \Ktilde^*R\Ktilde) \hat z_i}_{\cH_1} + \innerprod{(A + B\Ktilde)\hat z_i, P (A + B\Ktilde) \hat z_i}_{\cH_1} -\innerprod{\hat z_i, P\hat z_i}_{\cH_1}\\
%=& \innerprod{\hat z_i, (Q + \Ktilde^*R\Ktilde) z_i}_{\cH_1} + \innerprod{(A + B\Ktilde)\hat z_i, P (A + B\Ktilde)z_i}_{\cH_1} -\innerprod{\hat z_i, P\hat z_i}_{\cH_1}\\
=& \innerprod{\hat z_i, (Q + (\Ktilde - K + K)^*R(\Ktilde- K + K) \hat z_i}_{\cH_1} \nonumber\\
&+ \innerprod{(A + B\Ktilde -BK + BK)\hat z_i, P (A + B\Ktilde - BK + BK)\hat z_i}_{\cH_1} -\innerprod{\hat z_i, P\hat z_i}_{\cH_1}\\
=& \innerprod{\zhat_i, Q\zhat_i}_{\cH_1} + \innerprod{\zhat_i, (\Ktilde - K)^*R(\Ktilde - K)\zhat_i}_{\cH_1} + \innerprod{\zhat_i, K^*R(\Ktilde - K)\zhat_i}_{\cH_1}\nonumber \\
&+ \innerprod{\zhat_i, (\Ktilde - K)^*RK\zhat_i}_{\cH_1} +
\innerprod{\zhat_i, K^*RK\zhat_i}_{\cH_1}\nonumber\\
&+\innerprod{B(\Ktilde - K)\hat z_i, PB (\Ktilde - K)\hat z_i}_{\cH_1} + \innerprod{(A + BK)\hat z_i, P (A + BK)\hat z_i}_{\cH_1}\nonumber\\
&+\innerprod{B(\Ktilde - K)\hat z_i, P (A + BK)\hat z_i}_{\cH_1}+ \innerprod{(A + BK)\hat z_i, PB (\Ktilde - K)\hat z_i}_{\cH_1}\nonumber\\
&-\innerprod{\zhat_i, P\zhat_i}_{\cH_1}\\
=& \innerprod{\zhat_i, (\Ktilde - K)^*(R + B^*PB)(\Ktilde -K)\zhat_i}_{\cH_1}\nonumber\\
&+ 2\innerprod{\zhat_i, (\Ktilde - K)^*[RK  + B^*P(A+ BK)]\zhat_i}_{\cH_1}\\
=&\innerprod{\zhat_i, (\Ktilde - K)^*(R + B^*PB)(\Ktilde -K)\zhat_i}_{\cH_1}\nonumber\\
&+ 2\innerprod{\zhat_i, (\Ktilde - K)^*[(R + B^* P B)K  + B^*PA]\zhat_i}_{\cH_1}\\
=&\innerprod{\zhat_i, (\Ktilde - K)^*(R + B^*PB)(\Ktilde -K)\zhat_i}_{\cH_1}\nonumber\\
&+ 2\innerprod{\zhat_i, (\Ktilde - K)^*[(R + B^* P B)(-1)(R + B^*PB)^{-1}B^*PA  + B^*PA]\zhat_i}_{\cH_1}\\
=& \innerprod{\zhat_i, (\Ktilde - K)^*(R + B^*PB)(\Ktilde -K)\zhat_i}_{\cH_1}.
\end{align}

To conclude the proof of the theorem, we can apply Cauchy-Schwarz inequality, Lemma <ref>, and the fact that, since $\sigma_{\textrm{min}}(R)\geq 1$, it holds

\begin{align}
\norm{R + B^*PB}&\leq \norm {R} + \norm{B^*PB}\\
&\leq \norm {R} + \norm R\norm{B^*PB}\\
&\leq \sigma_{\textrm{max}}(R)\Gamma^3,
\end{align}

where we recall that $Γ\de 1 + \max(‖A‖,‖P‖,‖K‖,‖B‖)$ in the statement of Theorem <ref>.
We have that

\begin{align}
\cJhat - \cJ &= \lim_{T\to\infty}\sum_{i=0}^T\innerprod{\zhat_i, (\Ktilde - K)^*(R + B^*PB)(\Ktilde -K)\zhat_i}_{\cH_1}\\
&\leq \lim_{T\to\infty}\sum_{i=0}^T\norm{\zhat_i}_{\cH_1}^2\norm{\Ktilde - K}^2\norm{R + B^*PB}\\
&\leq 9 \sigma_{\textrm{max}}(R)\Gamma^9g(\epsilon)^2\lim_{T\to\infty}\sum_{i=0}^T\norm{\zhat_i}_{\cH_1}^2\\
&=9 \sigma_{\textrm{max}}(R)\Gamma^9g(\epsilon)^2\lim_{T\to\infty}\sum_{i=0}^T\norm{(A + B\Ktilde)^iz_0}_{\cH_1}^2\\
&\leq 9 \sigma_{\textrm{max}}(R)\Gamma^9g(\epsilon)^2\norm{z_0}_{\cH_1}^2\lim_{T\to\infty}\sum_{i=0}^T\norm{(A + B\Ktilde)^i}^2\\
&\leq 9 \sigma_{\textrm{max}}(R)\Gamma^9g(\epsilon)^2\norm{z_0}_{\cH_1}^2\lim_{T\to\infty}\sum_{i=0}^T\tau(A + BK, \rho)^2\left(\frac{1 + \rho}{2}\right)^{2i}\\
&\leq 9 \sigma_{\textrm{max}}(R)\Gamma^9g(\epsilon)^2\norm{z_0}_{\cH_1}^2\tau(A + BK, \rho)^2\frac{1}{1 - \left(\frac{1 + \rho}{2}\right)^2}\\
&\leq 36 \sigma_{\textrm{max}}(R)\Gamma^9g(\epsilon)^2\norm{z_0}_{\cH_1}^2\frac{\tau(A + BK, \rho)^2}{1 - \rho^2}.
\end{align}

To conclude, note that, according to Assumption <ref>,

\begin{align}
\norm{z_0}^2_{\cH_1}&=‖k(x_0,·)‖²_{\cH_1}\\
&=\ip{k(x_0,·), k(x_0,·)}_{\cH_1}\\
\end{align}