# Thompson Sampling Efficiently Learns to Control Diffusion Processes

Mohamad Kazem Shirani Faradonbeh, Mohamad Sadegh Shirani Faradonbeh, Mohsen
Bayati

###### Abstract

Diffusion processes that evolve according to linear stochastic differential
equations are an important family of continuous-time dynamic decision-making
models. Optimal policies are well-studied for them, under full certainty about
the drift matrices. However, little is known about data-driven control of
diffusion processes with uncertain drift matrices as conventional discrete-
time analysis techniques are not applicable. In addition, while the task can
be viewed as a reinforcement learning problem involving exploration and
exploitation trade-off, ensuring system stability is a fundamental component
of designing optimal policies. We establish that the popular Thompson sampling
algorithm learns optimal actions fast, incurring only a square-root of time
regret, and also stabilizes the system in a short time period. To the best of
our knowledge, this is the first such result for Thompson sampling in a
diffusion process control problem. We validate our theoretical results through
empirical simulations with real parameter matrices from two settings of
airplane and blood glucose control. Moreover, we observe that Thompson
sampling significantly improves (worst-case) regret, compared to the state-of-
the-art algorithms, suggesting Thompson sampling explores in a more guarded
fashion. Our theoretical analysis involves characterization of a certain
_optimality manifold_ that ties the local geometry of the drift parameters to
the optimal control of the diffusion process. We expect this technique to be
of broader interest.

##  1 Introduction

One of the most natural reinforcement learning (RL) algorithms for controlling
a diffusion process with unknown parameters is based on Thompson sampling (TS)
[1]: a Bayesian posterior for the model is calculated based on its time
evolution, and a control policy is then designed by treating a sampled model
from the posterior as the truth. Despite its simplicity, guaranteeing
efficiency and whether sampling the actions from the posterior could lead to
unbounded future trajectories is unknown. In fact, the only known such
theoretical result for control of a diffusion process is for an epsilon-greedy
type policy that requires selecting purely random actions at a certain rate
[2].

In this work, we consider a $p$ dimensional state signal
$\left\\{\boldsymbol{x}_{t}\right\\}_{t\geq 0}$ that obeys the (Ito)
stochastic differential equation (SDE)

$\boldsymbol{d}\boldsymbol{x}_{t}=\left(A_{0}\boldsymbol{x}_{t}+B_{0}\boldsymbol{u}_{t}\right)\boldsymbol{d}t+\boldsymbol{d}\mathbb{W}_{t}\,,$
(1)

where the _drift matrices_ $A_{0}$ and $B_{0}$ are unknown,
$\boldsymbol{u}_{t}\in\mathbb{R}^{q}$ is the control action at any time $t\geq
0$, and it is designed based on values of $\boldsymbol{x}_{s}$ for
$s\in[0,t]$. The matrix $B_{0}\in\mathbb{R}^{p\times q}$ models the influence
of the control action on the state evolution over time, while
$A_{0}\in\mathbb{R}^{p\times p}$ is the (open-loop) transition matrix
reflecting interactions between the coordinates of the state vector
$\boldsymbol{x}_{t}$. The diffusion term in (1) consists of a non-standard
Wiener process $\mathbb{W}_{t}$ that will be defined in the next section. The
goal is to study efficient RL policies that can design $\boldsymbol{u}_{t}$ to
minimize a quadratic cost function, defined in the next section, subject to
uncertainties around $A_{0}$ and $B_{0}$.

At a first glance, this problem is similar to most RL problems since the
optimal policy must balance between the two objectives of learning the unknown
matrices $A_{0}$ and $B_{0}$ (exploration) and optimally selecting the control
signals $\boldsymbol{u}_{t}$ to minimize the cost (exploitation). However,
unlike most RL problems that have finite or bounded-support state space,
ensuring _stability_ , that $\boldsymbol{x}_{t}$ stays bounded, is a crucial
part of designing optimal policies. For example, in the discrete-time version
of the problem, robust exploration is used to protect against unpredictably
unstable trajectories [3, 4, 5, 6].

#### Related literature.

The existing literature studies efficiency of TS for learning optimal
decisions in finite action spaces [7, 8, 9, 10, 11, 12]. In this stream of
research, it is shown that, over time, the posterior distribution concentrates
around low-cost actions [13, 14, 15]. TS is also studied in further discrete-
time settings with the environment represented by parameters that belong to a
continuum, and Bayesian and frequentist regret bounds are shown for linear-
quadratic regulators [16, 17, 18, 19, 20, 21, 22, 23, 24]. However,
effectiveness of TS in highly noisy environments that are modeled by diffusion
processes remains unexplored to date, due to technical challenges that will be
described below.

For continuous-time linear time invariant dynamical systems, infinite-time
consistency results are shown under a variety of technical assumptions,
followed by alternating policies that cause (small) linear regrets [25, 26,
27, 28, 29]. From a computational viewpoint, pure exploration algorithms for
computing optimal policies based on multiple trajectories of the state and
action data are studied as well [30, 31, 32], for which a useful survey is
available [33]. However, papers that study exploration versus exploitation,
and provide non-asymptotic estimation rates or regret bounds are limited to a
few recent work about offline RL or stabilized processes [34, 2, 35, 36].

#### Contributions.

This work, first establishes that TS learns to stabilize the diffusion process
(1). Specifically, in Theorem 1 of Section 3, we provide the first theoretical
stabilization guarantee for diffusion processes, showing that the probability
of preventing the state process from growing unbounded grows to 1, at an
exponential rate that depends on square-root of the time length devoted to
stabilization. As mentioned above, for RL problems with finite state spaces,
the process is by definition stabilized, regardless of the policy. However,
for the Euclidean state space of $\boldsymbol{x}_{t}$ in (1), stabilization is
necessary to ensure that the state and the cost do not grow unbounded.

Then, efficiency of TS in balancing exploration versus exploitation for
minimizing a cost function that has a quadratic form of both the state and the
control action is shown. Indeed, we establish in Theorem 2 of Section 4 that
the regret TS incurs, grows as the _square-root of time_ , while the squared
estimation error decays with the same rate. It is also shown that both the
above quantities grow quadratically with the dimension. To the authors’
knowledge, the presented results are the first theoretical analyses of TS for
learning to control diffusion processes.

Additionally, through extensive simulations we illustrate that TS enjoys
smaller average regret and substantially lower worst-case regret than the
existing RL policies, thanks to its informed exploration.

It is important to highlight that theoretical analysis of RL policies for
diffusion processes is highly non-trivial. Specifically, the conventional
discrete-time RL technical tools are not applicable, due to uncountable
cardinality of the random variables involved in a diffusion process, the
unavoidable dependence between them, and the high level of processing and
estimation noise. To address these, we make four main contributions. First,
non-asymptotic and uniform upper bounds for continuous-time martingales and
for Ito integrals are required to quantify the estimation accuracy. For that
purpose, we establish concentration inequalities and show sub-exponential tail
bounds for _double stochastic integrals_. Second, one needs sharp bounds for
the impact of estimation errors on eigenvalues of certain non-linear matrices
of the drift parameters that determine actions taken by TS policy. To tackle
that, we perform a novel and tight _eigenvalue perturbation-analysis_ based on
the approximation error, dimension, and spectrum of the matrices. We also
establish _Lipschitz continuity_ of the control policy with respect to the
drift matrices, by developing new techniques based on matrix-valued curves.
Third, to capture evaluation of both immediate and long-term effects of sub-
optimal actions, we employ _Ito calculus_ to bound the stochastic regret and
specify effects of all problem parameters. Finally, to study learning from
data trajectories that the condition number of their information matrix grows
unbounded, we develop stochastic inequalities for _self-normalized continuous-
time martingales_ , and _spectral analysis_ of non-linear functions of random
matrices.

#### Organization.

The organization of the subsequent sections is as follows. We formulate the
problem in Section 2, while Algorithm 1 that utilizes TS for learning to
stabilize the process and its high-probability performance guarantee are
presented in Section 3. Then, in Section 4, TS is considered for learning to
minimize a quadratic cost function, and the rates of estimation and regret are
established. Next, theoretical analysis are provided in Section 5, followed by
real-world numerical results of Section 6. Detailed proofs and auxiliary
lemmas are delegated to the appendices.

#### Notation.

The smallest (the largest) eigenvalue of matrix $M$, in magnitude, is denoted
by $\underline{\lambda}\left(M\right)$ ($\overline{\lambda}\left(M\right)$).
For a vector $a$,
${\left|\kern-1.50696pt\left|a\right|\kern-1.50696pt\right|}$ is the
$\ell_{2}$ norm, and for a matrix $M$,
${\left|\kern-1.50696pt\left|M\right|\kern-1.50696pt\right|}$ is the operator
norm that is the supremum of
${\left|\kern-1.50696pt\left|Ma\right|\kern-1.50696pt\right|}$ for $a$ on the
unit sphere. $\boldsymbol{N}\left(\mu,\Sigma\right)$ is Gaussian distribution
with mean $\mu$ and covariance $\Sigma$. If $\mu$ is a matrix (instead of
vector), then $\boldsymbol{N}\left(\mu,\Sigma\right)$ denotes a distribution
on matrices of the same dimension as $\mu$, such that all columns are
independent and share the covariance matrix $\Sigma$. In this paper,
transition matrices $A\in\mathbb{R}^{p\times p}$ together with input matrices
$B\in\mathbb{R}^{p\times q}$ are jointly denoted by the $(p+q)\times p$
parameter matrix $\boldsymbol{\theta}=\left[A,B\right]^{\top}$. We employ
$\vee$ ($\wedge$) for maximum (minimum). Finally, $a\lesssim b$ expresses that
$a\leq\alpha_{0}b$, for some fixed constant $\alpha_{0}$.

##  2 Problem Statement

We study the problem of designing provably efficient reinforcement learning
policies for minimizing a quadratic cost function in an uncertain linear
diffusion process. To proceed, fix the complete probability space
$(\Omega,\left\\{\mathcal{F}_{t}\right\\}_{t\geq 0},\mathbb{P})$, where
$\Omega$ is the sample space, $\left\\{\mathcal{F}_{t}\right\\}_{t\geq 0}$ is
a continuous-time filtration (i.e., increasing sigma-fields), and $\mathbb{P}$
is the probability measure defined on $\mathcal{F}_{\infty}$.

The state comprises the diffusion process $\boldsymbol{x}_{t}$ in (1), where
$\boldsymbol{\theta_{0}}=\left[A_{0},B_{0}\right]^{\top}\in\mathbb{R}^{\left(p+q\right)\times
p}$ is the unknown drift parameter. The diffusion term in (1) follows
infinitesimal variations of the $p$ dimensional Wiener process
$\left\\{\mathbb{W}_{t}\right\\}_{t\geq 0}$. That is,
$\left\\{\mathbb{W}_{t}\right\\}_{t\geq 0}$ is a multivariate Gaussian process
with independent increments and with the stationary covariance matrix
$\Sigma_{\mathbb{W}}$, such that for all $0\leq s_{1}\leq s_{2}\leq t_{1}\leq
t_{2}$,

$\begin{bmatrix}\mathbb{W}_{t_{2}}-\mathbb{W}_{t_{1}}\\\
\mathbb{W}_{s_{2}}-\mathbb{W}_{s_{1}}\end{bmatrix}\sim\boldsymbol{N}\left(\begin{bmatrix}0_{p}\\\
0_{p}\end{bmatrix},\begin{bmatrix}(t_{2}-t_{1})\Sigma_{\mathbb{W}}&0_{p\times
p}\\\ 0_{p\times p}&(s_{2}-s_{1})\Sigma_{\mathbb{W}}\end{bmatrix}\right).$ (2)

Existence, construction, continuity, and non-differentiability of Wiener
processes are well-known [37]. It is standard to assume that
$\Sigma_{\mathbb{W}}$ is positive definite, which is a common condition in
learning-based control [33, 34, 2, 35] to ensure accurate estimation over
time.

The RL policy designs the action $\left\\{\boldsymbol{u}_{t}\right\\}_{t\geq
0}$, based on the observed system state by the time, as well as the previously
applied actions, to minimize the long-run average cost

$\limsup\limits_{T\to\infty}\frac{1}{T}\int\limits_{0}^{T}{\left[\boldsymbol{x}_{t}^{\top},\boldsymbol{u}_{t}^{\top}\right]Q\begin{bmatrix}\boldsymbol{x}_{t}\\\
\boldsymbol{u}_{t}\end{bmatrix}}\boldsymbol{d}{t},~{}~{}~{}~{}~{}\text{for}~{}~{}~{}~{}Q=\begin{bmatrix}Q_{x}&Q_{xu}\\\
Q_{xu}^{\top}&Q_{u}\end{bmatrix}.$ (3)

Above, the cost is determined by the positive definite matrix $Q$, where
$Q_{x}\in\mathbb{R}^{p\times p}$, $Q_{u}\in\mathbb{R}^{q\times q}$,
$Q_{xu}\in\mathbb{R}^{p\times q}$. In fact, $Q$ determines the weights of
different coordinates of $\boldsymbol{x}_{t},\boldsymbol{u}_{t}$ in the cost
function, so that the policy aims to make the states small, by deploying small
actions. The cost matrix $Q$ is assumed known to the policy. Formally, the
problem is to minimize (3) by the policy

$\boldsymbol{u}_{t}=\boldsymbol{\widehat{\pi}}\left(Q,\left\\{\boldsymbol{x}_{s}\right\\}_{0\leq
s\leq t},\left\\{\boldsymbol{u}_{s}\right\\}_{0\leq s<t}\right).$ (4)

Without loss of generality, and for the ease of presentation, we follow the
canonical formulation that sets $Q_{xu}=0$; one can simply convert the case
$Q_{xu}\neq 0$ to the canonical form, by employing a rotation to
$\boldsymbol{x}_{t},\boldsymbol{u}_{t}$ [38, 39, 40, 41]. It is well-known
that if, hypothetically, the truth $\boldsymbol{\theta_{0}}$ was known, an
optimal policy $\boldsymbol{\pi_{\mathrm{opt}}}$ could be explicitly found by
solving the continuous-time algebraic Riccati equation. That is, for a generic
drift matrix $\boldsymbol{\theta}=\left[A,B\right]^{\top}$, finding the
symmetric $p\times p$ matrix ${P}\left(\boldsymbol{\theta}\right)$ that
satisfies

$A^{\top}{P}\left(\boldsymbol{\theta}\right)+{P}\left(\boldsymbol{\theta}\right)A-{P}\left(\boldsymbol{\theta}\right)BQ_{u}^{-1}B^{\top}{P}\left(\boldsymbol{\theta}\right)+Q_{x}=0.$
(5)

This means, for the true parameter
$\boldsymbol{\theta_{0}}=\left[A_{0},B_{0}\right]^{\top}$, we can let
${P}\left(\boldsymbol{\theta_{0}}\right)$ solve the above equation, and define
the policy

$\boldsymbol{\pi_{\mathrm{opt}}}:\>\>\>\>\>\boldsymbol{u}_{t}=-Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\boldsymbol{x}_{t},\>\>\>\forall
t\geq 0.$ (6)

It is known that the linear time-invariant policy
$\boldsymbol{\pi_{\mathrm{opt}}}$ minimizes the average cost in (3) [38, 39,
40, 41].

###### Definition 1

The process in (1) is stabilizable, if all eigenvalues of
$\overline{A}={A_{0}+B_{0}K}$ have negative real-parts, for a matrix $K$. Such
$K,\overline{A}$ are called a stabilizer and the stable closed-loop matrix.

We assume that the process (1) with the drift parameter
$\boldsymbol{\theta_{0}}$ is stabilizable. Therefore,
${P}\left(\boldsymbol{\theta_{0}}\right)$ exists, is unique, and can be
computed using continuous-time Riccati differential equations similar to (5),
except that the zero matrix on the right-hand side will be replaced by the
derivative of ${P}\left(\boldsymbol{\theta}\right)$ [38, 39, 40, 41].
Furthermore, it is known that real-parts of all eigenvalues of
$\overline{A}_{0}=A_{0}-B_{0}Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)$
are negative, i.e.,
$\left|\overline{\lambda}\left(\exp\left(\overline{A}_{0}t\right)\right)\right|<1$,
which means the matrix $\exp\left(\overline{A}_{0}t\right)$ decays
exponentially fast as $t$ grows [38, 39, 40, 41]. In the sequel, we use (5)
and refer to the solution ${P}\left(\boldsymbol{\theta}\right)$ for different
stabilizable $\boldsymbol{\theta}$. More details about the above optimal
feedback policy can be found in the aforementioned references.

In absence of exact knowledge of $\boldsymbol{\theta_{0}}$, a policy
$\boldsymbol{\widehat{\pi}}$ collects data and leverages it to approximate
$\boldsymbol{\pi_{\mathrm{opt}}}$ in (6). Therefore, at all (finite) times,
there is a gap between the cost of $\boldsymbol{\widehat{\pi}}$, compared to
that of $\boldsymbol{\pi_{\mathrm{opt}}}$. The cumulative performance
degradation due to this gap is the _regret_ of the policy
$\boldsymbol{\widehat{\pi}}$, that we aim to minimize. Technically, whenever
the control action $\boldsymbol{u}_{t}$ is designed by the policy
$\boldsymbol{\widehat{\pi}}$ according to (4), concatenate the resulting state
and input signals to get the observation
$\boldsymbol{z}_{t}(\boldsymbol{\widehat{\pi}})=\left[\boldsymbol{x}_{t}^{\top},\boldsymbol{u}_{t}^{\top}\right]^{\top}$.
If it is clear from the context, we drop $\boldsymbol{\widehat{\pi}}$.
Similarly, $\boldsymbol{z}_{t}(\boldsymbol{\pi_{\mathrm{opt}}})$ denotes the
observation signal of $\boldsymbol{\pi_{\mathrm{opt}}}$. Now, the regret at
time $T$ is defined by:

$\mathrm{Reg}_{\boldsymbol{\widehat{\pi}}}\left(T\right)=\int\limits_{0}^{T}{\left({\left|\kern-1.50696pt\left|Q^{1/2}\boldsymbol{z}_{t}(\boldsymbol{\widehat{\pi}})\right|\kern-1.50696pt\right|}^{2}-{\left|\kern-1.50696pt\left|Q^{1/2}\boldsymbol{z}_{t}(\boldsymbol{\pi_{\mathrm{opt}}})\right|\kern-1.50696pt\right|}^{2}\right)}\boldsymbol{d}{t}.$

A secondary objective is the learning accuracy of $\boldsymbol{\theta_{0}}$
from the single trajectory of the data generated by
$\boldsymbol{\widehat{\pi}}$. Letting $\boldsymbol{\widehat{\theta}_{t}}$ be
the parameter estimate at time $t$, we are interested in scaling of
${\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}_{t}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}$
with respect to $t$, $p$, and $q$.

##  3 Stabilizing the Diffusion Process

This section focuses on establishing that Thompson sampling (TS) learns to
stabilize the diffusion process (1). First, let us intuitively discuss the
problem of stabilizing unknown diffusion processes. Given that the optimal
policy in (6) stabilizes the process in (1), a natural candidate to obtain a
stable process under uncertainty of the drift matrices $A_{0},B_{0}$, is a
linear feedback of the form $\boldsymbol{u}_{t}=K\boldsymbol{x}_{t}$. So,
letting $\overline{A}=A_{0}+B_{0}K$, the solution of (1) is the
Ornstein–Uhlenbeck process
$\boldsymbol{x}_{t}=e^{\overline{A}t}\boldsymbol{x}_{0}+\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}$
[37]. Thus, if real-part of an eigenvalue of $\overline{A}$ is non-negative,
then the magnitude of $\boldsymbol{x}_{t}$ grows unbounded with $t$ [37].
Therefore, addressing instabilities of this form is important, _prior_ to
minimizing the cost. Otherwise, the regret grows (super) linearly with time.
In particular, if $A_{0}$ has some eigenvalue(s) with non-negative real-
part(s), then it is necessary to employ feedback to preclude instabilities.

In addition to minimizing the cost, the algebraic Riccati equation in (5)
provides a reliable and widely-used framework for stabilization, as discussed
after (6). Accordingly, due to uncertainty about $\boldsymbol{\theta_{0}}$,
one can solve (5) and find ${P}\left(\boldsymbol{\widehat{\theta}}\right)$,
only for an approximation $\boldsymbol{\widehat{\theta}}$ of
$\boldsymbol{\theta_{0}}$. Then, we expect to stabilize the system in (1) by
applying a linear feedback that is designed for the approximate drift matrix
$\boldsymbol{\widehat{\theta}}$. Technically, we need to ensure that all
eigenvalues of
$A_{0}-B_{0}Q_{u}^{-1}\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right)$
lie in the open left half-plane. To ensure that these requirements are met in
a sustainable manner, the main challenges are
(i) fast and accurate learning of $\boldsymbol{\theta_{0}}$ so that after a
short time period, a small error
$\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta_{0}}$ is guaranteed,
(ii) specifying the effect of the error
$\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta_{0}}$, on stability of
$A_{0}-B_{0}Q_{u}^{-1}\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right)$,
and
(iii) devising a remedy for the case that the stabilization procedure fails.
Note that the last challenge is unavoidable, since learning from finite data
can never be perfectly accurate, and so any finite-time stabilization
procedure has a (possibly small) positive failure probability.

Algorithm 1 addresses the above challenges by applying additionally randomized
control actions, and using them to provide a posterior belief $\mathcal{D}$
about $\boldsymbol{\theta_{0}}$. Note that the posterior is _not_ concentrated
at $\boldsymbol{\theta_{0}}$, and a sample $\boldsymbol{\widehat{\theta}}$
from $\mathcal{D}$ approximates $\boldsymbol{\theta_{0}}$, crudely. Still, the
theoretical analysis of Theorem 1 indicates that the failure probability of
Algorithm 1 decays exponentially fast with the length of the time interval it
is executed. Importantly, this small failure probability can shrink further by
repeating the procedure of sampling from $\mathcal{D}$. So, stabilization
under uncertainty is guaranteed, after a limited time of interacting with the
environment.

To proceed, let $\\{w_{n}\\}_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}}$ be a
sequence of independent Gaussian vectors with the distribution
$w_{n}\sim\boldsymbol{N}\left(0,\sigma_{w}^{2}I_{q}\right)$, for some fixed
constant $\sigma_{w}$. Suppose that we aim to devote the time length
$\boldsymbol{\tau}$ to collect observations for learning to stabilize. Note
that since stabilization is performed before moving forward to the main
objective of minimizing the cost functions, the stabilization time length
$\boldsymbol{\tau}$ is desired to be as short as possible. We divide this time
interval of length $\boldsymbol{\tau}$ to $\boldsymbol{\kappa_{\mathrm{}}}$
sub-intervals of equal length, and randomize an initial linear feedback policy
by adding $\left\\{w_{n}\right\\}_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}}$.
That is, for $n=0,1,\cdots,\boldsymbol{\kappa_{\mathrm{}}}-1$, Algorithm 1
employs the control action

$\boldsymbol{u}_{t}=K\boldsymbol{x}_{t}+w_{n},\text{
~{}~{}~{}~{}~{}~{}for~{}~{}~{}~{}~{}~{}~{}
}\frac{n\boldsymbol{\tau}}{\boldsymbol{\kappa_{\mathrm{}}}}\leq
t<\frac{(n+1)\boldsymbol{\tau}}{\boldsymbol{\kappa_{\mathrm{}}}},$ (7)

where $K$ is an initial stabilizing feedback so that all eigenvalues of
$A_{0}+B_{0}K$ lie in the open left half-plane. In practice, such $K$ is
easily found using physical knowledge of the model, e.g., via conservative
control sequence for an airplane [42, 43]. However, note that such actions are
sub-optimal involving large regrets. Therefore, they are only temporarily
applied, for the sake of data collection. Then, the data collected during the
time interval $0\leq t\leq\boldsymbol{\tau}$ will be utilized by the algorithm
to determine the posterior belief $\mathcal{D}_{\boldsymbol{\tau}}$, as
follows. Recalling the notation
$\boldsymbol{z}_{t}^{\top}=\left[\boldsymbol{x}_{t}^{\top},\boldsymbol{u}_{t}^{\top}\right]$,
let $\widehat{\mu}_{0},\widehat{\Sigma}_{0}$ be the mean and the precision
matrix of a prior normal distribution on $\boldsymbol{\theta_{0}}$ (using the
notation defined in Section 1 for random matrices). Nonetheless, if there is
no such prior, we simply let $\widehat{\mu}_{0}=0_{\left(p+q\right)\times p}$
and $\widehat{\Sigma}_{0}=I_{p+q}$. Then, define

$\displaystyle\widehat{\Sigma}_{\boldsymbol{\tau}}=\widehat{\Sigma}_{0}+\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}\boldsymbol{z}_{s}^{\top}}\boldsymbol{d}{s},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\widehat{\mu}_{\boldsymbol{\tau}}=\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\left(\widehat{\Sigma}_{0}\widehat{\mu}_{0}+\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}}\boldsymbol{d}{\boldsymbol{x}_{s}^{\top}}\right).$
(8)

Using
$\widehat{\Sigma}_{\boldsymbol{\tau}}\in\mathbb{R}^{\left(p+q\right)\times\left(p+q\right)}$
together with the mean matrix $\widehat{\mu}_{\boldsymbol{\tau}}$, Algorithm 1
forms the posterior belief

$\mathcal{D}_{\boldsymbol{\tau}}=\boldsymbol{N}\left(\widehat{\mu}_{\boldsymbol{\tau}},\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\right),$
(9)

about the drift parameter $\boldsymbol{\theta_{0}}$. So, as defined in the
notation, the posterior distribution of every column $i=1,\cdots,p$ of
$\boldsymbol{\theta_{0}}$, is an independent multivariate normal with the
covariance matrix $\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}$, while the mean
is the column $i$ of $\widehat{\mu}_{\boldsymbol{\tau}}$. The final step of
Algorithm 1 is to output a sample $\boldsymbol{\widehat{\theta}}$ from
$\mathcal{D}_{\boldsymbol{\tau}}$.

Algorithm 1 : Stabilization under Uncertainty

Inputs: initial feedback $K$, stabilization time length $\boldsymbol{\tau}$

for $n=0,1,\cdots,\boldsymbol{\kappa_{\mathrm{}}}-1$ do

while $n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}\leq
t<(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}$ do

Apply control action $\boldsymbol{u}_{t}$ in (7)

end while

end for

Calculate
$\widehat{\Sigma}_{\boldsymbol{\tau}},\widehat{\mu}_{\boldsymbol{\tau}}$
according to (8)

Return sample $\boldsymbol{\widehat{\theta}}$ from the distribution
$\mathcal{D}_{\boldsymbol{\tau}}$ in (9)

Next, to establish performance guarantees for Algorithm 1, let us quantify the
_ideal_ stability by

$\boldsymbol{\zeta_{0}}=-\log\overline{\lambda}\left(\exp\left[A_{0}-B_{0}Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\right]\right).$
(10)

By definition, $\boldsymbol{\zeta_{0}}$ is positive. In fact, it is the
smallest distance between the imaginary axis in the complex-plane, and the
eigenvalues of the transition matrix
$\overline{A}_{0}=A_{0}-B_{0}Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)$,
under the optimal policy in (6). Since $\boldsymbol{\theta_{0}}$ is
unavailable, it is _not_ realistic to expect that after applying a policy
based on $\boldsymbol{\widehat{\theta}}$ given by Algorithm 1, real-parts of
all eigenvalues of the resulting matrix
$A_{0}-B_{0}Q_{u}^{-1}\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right)$
are at most $-\boldsymbol{\zeta_{0}}$. However, $\boldsymbol{\zeta_{0}}$ is
crucial in studying stabilization, such that stabilizing controllers for
systems with larger _$\boldsymbol{\zeta_{0}}$_ can be learned faster. The
exact effect of this quantity, as well as those of other properties of the
diffusion process, are formally established in the following result.
Informally, the failure probability of Algorithm 1 decays exponentially with
$\boldsymbol{\tau}^{1/2}$.

######  Theorem 1 (Stabilization Guarantee)

For the sample $\boldsymbol{\widehat{\theta}}$ given by Algorithm 1, let
$\mathcal{E}_{\boldsymbol{\tau}}$ be the failure event that
$A_{0}-B_{0}Q_{u}^{-1}\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right)$
has an eigenvalue in the closed right half-plane. Then, if
$\boldsymbol{\kappa_{\mathrm{}}}\gtrsim\boldsymbol{\tau}^{2}$, we have

$\log\mathbb{P}{\left(\mathcal{E}_{\boldsymbol{\tau}}\right)}\lesssim-~{}~{}~{}\frac{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}}{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\vee\sigma_{w}^{2}}~{}~{}\frac{1\wedge\boldsymbol{\zeta_{0}}^{p}}{1\vee{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}^{3}}~{}~{}\sqrt{\frac{\boldsymbol{\tau}}{p^{3}q}}.$
(11)

The above result indicates that more heterogeneity in coordinates of the
Wiener noise renders stabilization harder. Moreover, using (10), the term
$1\wedge\boldsymbol{\zeta_{0}}^{p}$ reflects that less stable diffusion
processes with smaller $\boldsymbol{\zeta_{0}}$, are significantly harder to
stabilize under uncertainty. Also as one can expect, larger dimensions make
learning to stabilize harder. This is contributed by higher number of
parameters to learn, as well as higher sensitivity of eigenvalues for
processes of larger dimensions. Finally, the failure probability decays as
$\boldsymbol{\tau}^{1/2}$, mainly because continuous-time martingales have
sub-exponential distributions, unlike sub-Gaussianity of discrete-time
counterparts [44, 45, 46].

##  4 Thompson Sampling for Efficient Control: Algorithm and Theory

In this section, we proceed towards analysis of Thompson sampling (TS) for
minimizing the quadratic cost in (3), and show that it efficiently learns the
optimal control actions. That is, TS balances the exploration versus
exploitation, such that its regret grows with (nearly) the square-root rate,
as time grows. In the sequel, we introduce Algorithm 2 and discuss the
conceptual and technical frameworks it relies on. Then, we establish
efficiency by showing regret bounds in terms of different problem parameters
and provide the rates of estimating the unknown drift matrices.

In Algorithm 2, first the learning-based stabilization Algorithm 1 is run
during the time period $0\leq t<\boldsymbol{\tau}_{0}$. So, according to
Theorem 1, the optimal feedback of $\boldsymbol{\widehat{\theta}_{0}}$
stabilizes the system with a high probability, as long as
$\boldsymbol{\tau}_{0}$ is sufficiently large. Note that if growth of the
state indicates that Algorithm 1 failed to stabilize, one can repeat sampling
from $\mathcal{D}_{\boldsymbol{\tau}_{0}}$. So, we can assume that the
evolution of the controlled diffusion process remains stable when Algorithm 2
is being executed. On the other hand, the other benefit of running Algorithm 1
at the beginning is that it performs an initial exploration phase that will be
utilized by Algorithm 2 to minimize the regret.

Then, in order to learn the optimal policy $\boldsymbol{\pi_{\mathrm{opt}}}$
with minimal sub-optimality, RL algorithms need to cope with a fundamental
challenge, commonly known as the exploration-exploitation dilemma. To see
that, first note that an acceptable policy that aims to have sub-linear
regret, needs to take near-optimal control actions in a long run;
$\boldsymbol{u}_{t}\approx-
Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\boldsymbol{x}_{t}$.
Although such policies exploit well and their control actions are close to
that of $\boldsymbol{\pi_{\mathrm{opt}}}$, their regret grows large since they
fail to explore. Technically, the trajectory of observations
$\left\\{\boldsymbol{z}_{t}\right\\}_{t\geq 0}$ is not rich enough to provide
accurate estimations, since in
$\boldsymbol{z}_{t}^{\top}=\left[\boldsymbol{x}_{t}^{\top},\boldsymbol{u}_{t}^{\top}\right]$,
the signal $\boldsymbol{u}_{t}$ is (almost) a linear function of the state
signal $\boldsymbol{x}_{t}$, and so does not contribute towards gathering
information about the unknown parameter $\boldsymbol{\theta_{0}}$. Conversely,
for sufficient explorations, RL policies need to take actions that deviate
from those of $\boldsymbol{\pi_{\mathrm{opt}}}$, which imposes large regret
(as quantified in Lemma 7). Accordingly, the above trade-off needs to be
delicately balanced; what we show that TS does.

Algorithm 2 is episodic; the parameter estimates
$\boldsymbol{\widehat{\theta}_{n}}$ are updated only at the end of the
episodes at times $\left\\{\boldsymbol{\tau}_{n}\right\\}_{n=0}^{\infty}$,
while during every episode, actions are taken _as if_
$\boldsymbol{\widehat{\theta}_{n}}=\left[\widehat{A}_{n},\widehat{B}_{n}\right]^{\top}$
is the unknown truth $\boldsymbol{\theta_{0}}$. That is, for
$\boldsymbol{\tau}_{n-1}\leq t<\boldsymbol{\tau}_{n}$, using
${P}\left(\boldsymbol{\widehat{\theta}_{n}}\right)$ in (5), we let
$\boldsymbol{u}_{t}=-Q_{u}^{-1}\widehat{B}_{n}^{\top}{P}\left(\boldsymbol{\widehat{\theta}_{n}}\right)\boldsymbol{x}_{t}$.
Then, for each $n=1,2,\cdots$, at time $\boldsymbol{\tau}_{n}$, we use all the
observations collected so far, to find
$\widehat{\Sigma}_{\boldsymbol{\tau}_{n}},\widehat{\mu}_{\boldsymbol{\tau}_{n}}$
according to (8). Next, we use them to sample
$\boldsymbol{\widehat{\theta}_{n}}$ from the posterior
$\mathcal{D}_{\boldsymbol{\tau}_{n}}$ in (9).

The episodes in Algorithm 2 are chosen such that their end points satisfy

$0<\underline{\alpha}\leq\inf_{n\geq
0}\frac{\boldsymbol{\tau}_{n+1}-\boldsymbol{\tau}_{n}}{\boldsymbol{\tau}_{n}}\leq\sup_{n\geq
0}\frac{\boldsymbol{\tau}_{n+1}-\boldsymbol{\tau}_{n}}{\boldsymbol{\tau}_{n}}\leq\overline{\alpha}<\infty,$
(12)

for some fixed constants $\underline{\alpha},\overline{\alpha}$. Broadly
speaking, (12) lets the episode lengths of Algorithm 2 scale properly to avoid
unnecessary updates of parameter estimates, while at the same time performing
sufficient exploration. To see that, first note that since
$\widehat{\Sigma}_{\boldsymbol{\tau}}$ grows with $\boldsymbol{\tau}$, the
estimation error $\boldsymbol{\widehat{\theta}_{n}}-\boldsymbol{\theta_{0}}$
decays (at best polynomially fast) with $\boldsymbol{\tau}_{n}$. So, until
ensuring that updating the posterior yields to significantly better
approximations, it will not be beneficial to update it, sample from it, and
solve (5). So, the period $\boldsymbol{\tau}_{n+1}-\boldsymbol{\tau}_{n}$ that
the data up to time $\boldsymbol{\tau}_{n}$ is utilized, is set to be as long
as $\underline{\alpha}\boldsymbol{\tau}_{n}$. On the other hand, the above
period cannot be too long, since we aim to improve the parameter estimates
after collecting enough new observations;
$\boldsymbol{\tau}_{n+1}\leq\left(1+\overline{\alpha}\right)\boldsymbol{\tau}_{n}$.
A simple setting is to let $\underline{\alpha}=\overline{\alpha}$, which
yields to exponential episodes
$\boldsymbol{\tau}_{n}=\boldsymbol{\tau}_{0}\left(1+\overline{\alpha}\right)^{n}$.
Note that for TS in continuous time, posterior updates should be limited to
sufficiently-apart time points. Otherwise, repetitive updates are
computationally impractical, and also can degrade the performance by
preventing control actions from having enough time to effectively influence.

Algorithm 2 : Thompson Sampling for Efficient Control of Diffusion Processes

Inputs: stabilization time $\boldsymbol{\tau}_{0}$

Calculate sample $\boldsymbol{\widehat{\theta}_{0}}$ by running Algorithm 1
for time $\boldsymbol{\tau}_{0}$

for $n=1,2,\cdots$ do

while $\boldsymbol{\tau}_{n-1}\leq t<\boldsymbol{\tau}_{n}$ do

Apply control action
$\boldsymbol{u}_{t}=-Q_{u}^{-1}\widehat{B}_{n-1}^{\top}{P}\left(\boldsymbol{\widehat{\theta}_{n-1}}\right)\boldsymbol{x}_{t}$

end while

Letting
$\widehat{\Sigma}_{\boldsymbol{\tau}_{n}},\widehat{\mu}_{\boldsymbol{\tau}_{n}}$
be as (8), sample $\boldsymbol{\widehat{\theta}_{n}}$ from
$\mathcal{D}_{\boldsymbol{\tau}_{n}}$ given in (9)

end for

We show next that Algorithm 2 addresses the exploration-exploitation trade-off
efficiently. To see the intuition, consider the sequence of posteriors
$\mathcal{D}_{\boldsymbol{\tau}_{n}}$. The explorations Algorithm 2 performs
by sampling $\boldsymbol{\widehat{\theta}_{n}}$ from
$\mathcal{D}_{\boldsymbol{\tau}_{n}}$, depends on
$\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}$. Now, if hypothetically
$\underline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}\right)$ is
not large enough, then $\mathcal{D}_{\boldsymbol{\tau}_{n}}$ does not
sufficiently concentrate around $\widehat{\mu}_{\boldsymbol{\tau}_{n}}$ and so
$\boldsymbol{\widehat{\theta}_{n}}$ will probably deviate from the previous
samples $\left\\{\boldsymbol{\widehat{\theta}_{i}}\right\\}_{i=1}^{n-1}$. So,
the algorithm explores more and obtains richer data $\boldsymbol{z}_{t}$ by
diversifying the control signal $\boldsymbol{u}_{t}$. This renders the next
mean $\widehat{\mu}_{\boldsymbol{\tau}_{n+1}}$ a more accurate approximation
of $\boldsymbol{\theta_{0}}$, and also makes
$\underline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}_{n+1}}\right)$
grow faster than before. Thus, the next posterior
$\mathcal{D}_{\boldsymbol{\tau}_{n+1}}$ provides a better sample with smaller
estimation error
$\boldsymbol{\widehat{\theta}_{n+1}}-\boldsymbol{\theta_{0}}$. Similarly, if a
posterior is excessively concentrated, in a few episodes the posteriors adjust
accordingly to the proper level of exploration. Hence, TS eventually balances
the exploration versus the exploitation. This is formalized below.

######  Theorem 2 (Regret and Estimation Rates)

Parameter estimates and regret of Algorithm 2, satisfy

$\displaystyle{\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}_{n}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}^{2}$
$\displaystyle\lesssim$
$\displaystyle\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)}\log\left(1+\overline{\alpha}\right)~{}~{}~{}\left(p+q\right)p~{}~{}~{}\boldsymbol{\tau}_{n}^{-1/2}\log\boldsymbol{\tau}_{n}\,,$
$\displaystyle\mathrm{Reg}\left(T\right)$ $\displaystyle\lesssim$
$\displaystyle\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}\right)\boldsymbol{\tau}_{0}+\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)^{2}}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)}\frac{\overline{\alpha}{\left|\kern-1.29167pt\left|{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}^{6}}{\log(\underline{\alpha}+1)\underline{\lambda}\left(Q\right)^{6}}~{}~{}~{}\left(p+q\right)p~{}~{}~{}T^{1/2}\log
T\,.$

In the above regret and estimation rates, and similar to Theorem 1,
$\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)/\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)$
reflects the impact of heterogeneity in coordinates of $\mathbb{W}_{t}$ on the
quality of learning. Also, larger $\log(1+\overline{\alpha})$ corresponds to
longer episodes which compromises the estimation. Further, $p(p+q)$ shows that
larger number of parameters linearly worsens the learning accuracy. In the
regret bound,
${\left|\kern-1.29167pt\left|{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}/\underline{\lambda}\left(Q\right)$
indicates effect of the true problem parameters $\boldsymbol{\theta_{0}},Q$.
Finally,
$\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}\right)\boldsymbol{\tau}_{0}$
captures the initial phase that Algorithm 1 is run for stabilization, which
takes sub-optimal control actions as in (7).

##  5 Intuition and Summary of the Analysis

The goal of this section is to provide a high-level roadmap of the proofs of
Theorems 1 and 2, and convey the main intuition behind the analysis. Complete
proofs and the technical lemmas are provided in Appendices A and B,
respectively.

#### Summary of the Proof of Theorem 1.

The main steps involve analyzing the estimation (Lemma 4), studying its effect
on the solutions of (5) (Lemma 12), and characterizing impact of errors in
entries of parameter matrices on their eigenvalues (Lemma 5). Next, we
elaborate on these steps.

We show that the error satisfies
${\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}\lesssim
p(p+q)^{1/2}\boldsymbol{\tau}^{-1/2}$ (Lemma 4). More precisely, the error
depends mainly on total strength of the observation signals
$\boldsymbol{z}_{t}$, which are captured in the precision matrix
$\widehat{\Sigma}_{\boldsymbol{\tau}}$, as well as total interactions between
the signal $\boldsymbol{z}_{t}$ and the noise $\mathbb{W}_{t}$ in the form of
the stochastic integral matrix
$\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}$.
However, we establish an upper bound
$\overline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\right)\lesssim\boldsymbol{\tau}^{-1}$,
that indicates the concentration rate of the posterior
$\mathcal{D}_{\boldsymbol{\tau}}$ (Lemma 3). Similarly, thanks to the
randomization signal $w_{n}$, the signals $\boldsymbol{z}_{t}$ are diverse
enough to effectively explore the set of matrices
$\boldsymbol{\theta}=\left[A,B\right]^{\top}$, leading to accurate
approximation of $\boldsymbol{\theta_{0}}$ by the posterior mean matrix
$\widehat{\mu}_{\boldsymbol{\tau}}$. Then, to bound the error terms caused by
the Wiener noise $\mathbb{W}_{t}$, we establish the rate
$p(p+q)^{1/2}\boldsymbol{\tau}^{1/2}$ (Lemma 2). Indeed, we show that the
entries of this error matrix are continuous-time martingales, and use
exponential inequalities for quadratic forms and double stochastic integrals
[45, 44] to establish that they have a sub-exponential distribution.

Moreover, the error rate of the feedback satisfies a similar property;
${\left|\kern-1.29167pt\left|\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right)-B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}\lesssim
p(p+q)^{1/2}\boldsymbol{\tau}^{-1/2}$ (Lemma 12). So, letting
$\overline{A}=A_{0}-B_{0}Q_{u}^{-1}\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right)$
and
$\overline{A}_{0}=A_{0}-B_{0}Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)$,
it holds that
${\left|\kern-1.29167pt\left|\overline{A}-\overline{A}_{0}\right|\kern-1.29167pt\right|}\lesssim
p(p+q)^{1/2}\boldsymbol{\tau}^{-1/2}$. Next, to consider the effect of the
errors on the eigenvalues of $\overline{A}$, we compare them to the
eigenvalues of $\overline{A}_{0}$, which are bounded by
$-\boldsymbol{\zeta_{0}}$ in (10). To that end, we establish a novel and tight
perturbation analysis for eigenvalues of matrices, with respect to their
entries and spectral properties (Lemma 5). Using that, we show that the
difference between the eigenvalues of $\overline{A}$ and $\overline{A}_{0}$
scales as
$\left(1\vee\boldsymbol{r}^{1/2}{\left|\kern-1.29167pt\left|\overline{A}-\overline{A}_{0}\right|\kern-1.29167pt\right|}\right)^{1/\boldsymbol{r}}$,
where $\boldsymbol{r}$ is the size of the largest block in the Jordan block-
diagonalization of $\overline{A}_{0}$. Therefore, for stability of
$\overline{A}$, we need
${\left|\kern-1.29167pt\left|\overline{A}-\overline{A}_{0}\right|\kern-1.29167pt\right|}\lesssim
p^{-1/2}\left(1\wedge\boldsymbol{\zeta_{0}}^{p}\right)$, since
$\boldsymbol{r}\leq p$. Note that if $\overline{A}_{0}$ is diagonalizable,
$\boldsymbol{r}=1$ implies that we can replace the above upper bound by
$1\wedge\boldsymbol{\zeta_{0}}$. Putting this stability result together with
the estimation error in the previous paragraph, we obtain (11).

#### Summary of the Proof of Theorem 2.

To establish the estimation rates, we develop multiple intermediate lemmas
quantifying the exact amount of exploration Algorithm 2 performs. First, we
utilize the fact that the bias of the posterior distribution
$\mathcal{D}_{\boldsymbol{\tau}_{n}}$ depends on its covariance matrix
$\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}$, as well as a self-normalized
continuous-time matrix-valued martingale. For the effect of the former, i.e.,
$\overline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}^{-1/2}\right)$,
we show an upper-bound of the order $\boldsymbol{\tau}_{n}^{-1/4}$ (Lemma 9).
To that end, the local geometry of the optimality manifolds that contain drift
parameters $\boldsymbol{\theta}$ that has the same optimal feedback as that of
the unknown truth $\boldsymbol{\theta_{0}}$ in (6) are fully specified (Lemma
6), and spectral properties of non-linear functions of random matrices are
studied. Then, we establish a stochastic inequality for the self-normalized
martingale, indicating that its scaling is of the order
$p(p+q)\log\boldsymbol{\tau}_{n}$ (Lemma 8). Therefore, utilizing the fact
that $\boldsymbol{\widehat{\theta}_{n}}-\widehat{\mu}_{\boldsymbol{\tau}_{n}}$
has the same scaling as the bias matrix
$\widehat{\mu}_{\boldsymbol{\tau}_{n}}-\boldsymbol{\theta_{0}}$, we obtain the
estimation rates of Theorem 2.

Next, to prove the presented regret bound, we establish a delicate and tight
analysis for the dominant effect of the control signal $\boldsymbol{u}_{t}$ on
the regret Algorithm 2 incurs. Technically, by carefully examining the
infinitesimal influences of the control actions at every time on the cost, we
show that it suffices to integrate the squared deviations
${\left|\kern-1.29167pt\left|\boldsymbol{u}_{t}+Q_{u}^{-1}\widehat{B}_{n}^{\top}{P}\left(\boldsymbol{\widehat{\theta}_{n}}\right)\boldsymbol{x}_{t}\right|\kern-1.29167pt\right|}^{2}$
to obtain $\mathrm{Reg}\left(T\right)$ (Lemma 7). We proceed toward specifying
the effect of the exploration Algorithm 2 performs on its exploitation
performance by proving the Lipschitz continuity of the solutions of the
Riccati equation (5) with respect to the drift parameters:
${\left|\kern-1.29167pt\left|{P}\left(\boldsymbol{\widehat{\theta}_{n}}\right)-{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}\lesssim{\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}_{n}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}$
(Lemma 12). This result is a very important property of (5) that lets the
rates of deviations from the optimal action scale the same as the estimation
error, and is proven by careful analysis of integration along matrix-valued
curves in the space of drift matrices, as well as spectral analysis for
approximate solutions of a Lyapunov equation (Lemma 10). Thus, the regret
bound is achieved, using the estimation error result in Theorem 2.

##  6 Numerical Analysis

We empirically evaluate the theoretical results of Theorems 1 and 2 under
three control problems. The first two are for the flight control of X-29A
airplane at 2000 ft [42] and for Boeing 747 [43]. The third simulation is for
blood glucose control [47]. We present the results for X-29A airplane in this
section, and defer the other two examples to the appendix. The true drift
matrices of the X-29A airplane are
$A_{0}=\begin{bmatrix}-0.16\;&0.07\;&-1.00\;&0.04\\\
-15.20\;&-2.60\;&1.11\;&0.00\\\ 6.84\;&-0.10\;&-0.06\;&0.00\\\
0.00\;&1.00\;&0.07\;&0.00\end{bmatrix},~{}~{}~{}~{}B_{0}=\begin{bmatrix}-0.0006\;&0.0007\\\
1.3430\;&0.2345\\\ 0.0897\;&-0.0710\\\ 0.0000\;&0.0000\end{bmatrix}$. Further,
we let $\Sigma_{\mathbb{W}}=0.25~{}I_{p}$, $Q_{x}=I_{p}$, and
$Q_{u}=0.1~{}I_{q}$ where $I_{n}$ is the $n$ by $n$ identity matrix. To update
the diffusion process $\boldsymbol{x}_{t}$ in (1), time-steps of length
$10^{-3}$ are employed. Then, in Algorithm 1, we let
$\sigma_{w}=5,\boldsymbol{\kappa_{\mathrm{}}}=\lfloor\boldsymbol{\tau}^{3/2}\rfloor$,
while $\boldsymbol{\tau}$ varies from $4$ to $20$ seconds. The initial
feedback $K$ is generated randomly. The results for 1000 repetitions are
depicted on the left plot of Figure 1, confirming Theorem 1 that the failure
probability of stabilization, decreases exponentially in $\boldsymbol{\tau}$.

On the right hand side of Figure 1, Algorithm 2 is executed for $600$ second,
for $\boldsymbol{\tau}_{n}=20\times 1.1^{n}$. We compare TS with the
_Randomized Estimate_ algorithm [2] for $100$ different repetitions. Average-
and worst-case values of the estimation error and the regret are reported,
both normalized by their scaling with time and dimension, as in Theorem 2. The
graphs show that (especially the worst-case) regret of TS substantially
outperforms, suggesting that TS explores in a more robust fashion. Simulations
for Boeing 747 and for the blood glucose control in Appendix D corroborate the
above findings.

Figure 1: For the X-29A flight control problem, percentage of stabilization
for 1000 runs of Algorithm 1 is plotted on the left. The graphs on the right
depict the performance of Algorithm 2 (blue) compared to Randomized Estimate
policy (red) [2]. The top graph plots the normalized squared estimation error,
${\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}_{n}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}^{2}$
divided by $p(p+q)\boldsymbol{\tau}_{n}^{-1/2}\log\boldsymbol{\tau}_{n}$,
versus time, while the lower one showcases the regret
$\mathrm{Reg}\left(T\right)$, normalized by $p(p+q)T^{1/2}\log T$. Curves for
the worst-case among $100$ replications are provided for both quantities, as
well as for the averages over all replicates.

##  7 Concluding Remarks and Future Work

We studied Thompson sampling (TS) RL policies to control a diffusion process
with unknown drift matrices. First, we proposed a stabilization algorithm for
linear diffusion processes, and established that its failure probability
decays exponentially with time. Further, efficiency of TS in balancing
exploration versus exploitation for minimizing a quadratic cost function is
shown. More precisely, regret bounds growing as square-root of time and square
of dimensions are established for Algorithm 2. Empirical studies showcasing
superiority of TS over state-of-the-art are provided as well.

As the first theoretical analysis of TS for control of a continuous-time
model, this work implies multiple important future directions. Establishing
minimax regret lower-bounds for diffusion process control problem is yet
unanswered. Moreover, studying the performance of TS for robust control of the
diffusion processes aiming to simultaneously minimize the cost function for a
family of drift matrices, is also an interesting direction for further
investigation. Another problem of interest is efficiency of TS for learning to
control under partial observation where the state is not observed and instead
a noisy linear function of the state is available as the output signal.

## References

  * [1] W. R. Thompson, “On the likelihood that one unknown probability exceeds another in view of the evidence of two samples,” _Biometrika_ , vol. 25, no. 3/4, pp. 285–294, 1933.
  * [2] M. K. S. Faradonbeh and M. S. S. Faradonbeh, “Efficient estimation and control of unknown stochastic differential equations,” _arXiv preprint arXiv:2109.07630_ , 2021.
  * [3] P. A. Ioannou and J. Sun, _Robust adaptive control_. PTR Prentice-Hall Upper Saddle River, NJ, 1996, vol. 1.
  * [4] F. L. Lewis, L. Xie, and D. Popa, _Optimal and robust estimation: with an introduction to stochastic control theory_. CRC press, 2017.
  * [5] A. Subrahmanyam and G. P. Rao, _Identification of Continuous-time Systems: Linear and Robust Parameter Estimation_. CRC Press, 2019.
  * [6] J. Umenberger, M. Ferizbegovic, T. B. Schön, and H. Hjalmarsson, “Robust exploration in linear quadratic reinforcement learning,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [7] S. Agrawal and N. Goyal, “Analysis of thompson sampling for the multi-armed bandit problem,” in _Conference on learning theory_. JMLR Workshop and Conference Proceedings, 2012, pp. 39–1.
  * [8] ——, “Further optimal regret bounds for thompson sampling,” in _Artificial intelligence and statistics_. PMLR, 2013, pp. 99–107.
  * [9] A. Gopalan and S. Mannor, “Thompson sampling for learning parameterized markov decision processes,” in _Conference on Learning Theory_. PMLR, 2015, pp. 861–898.
  * [10] M. J. Kim, “Thompson sampling for stochastic control: The finite parameter case,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 12, pp. 6415–6422, 2017.
  * [11] M. Abeille and A. Lazaric, “Linear thompson sampling revisited,” in _Artificial Intelligence and Statistics_. PMLR, 2017, pp. 176–184.
  * [12] N. Hamidi and M. Bayati, “On worst-case regret of linear thompson sampling,” _arXiv preprint arXiv:2006.06790_ , 2020.
  * [13] D. Russo and B. Van Roy, “Learning to optimize via posterior sampling,” _Mathematics of Operations Research_ , vol. 39, no. 4, pp. 1221–1243, 2014\.
  * [14] ——, “An information-theoretic analysis of thompson sampling,” _The Journal of Machine Learning Research_ , vol. 17, no. 1, pp. 2442–2471, 2016.
  * [15] D. Russo, B. Van Roy, A. Kazerouni, I. Osband, and Z. Wen, “A tutorial on thompson sampling,” _arXiv preprint arXiv:1707.02038_ , 2017.
  * [16] M. Abeille and A. Lazaric, “Improved regret bounds for thompson sampling in linear quadratic control problems,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1–9.
  * [17] Y. Ouyang, M. Gagrani, and R. Jain, “Posterior sampling-based reinforcement learning for control of unknown linear systems,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 8, pp. 3600–3607, 2019.
  * [18] M. K. S. Faradonbeh, A. Tewari, and G. Michailidis, “Finite-time adaptive stabilization of linear systems,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 8, pp. 3498–3505, 2018.
  * [19] ——, “Randomized algorithms for data-driven stabilization of stochastic linear systems,” in _2019 IEEE Data Science Workshop (DSW)_. IEEE, 2019, pp. 170–174.
  * [20] ——, “On applications of bootstrap in continuous space reinforcement learning,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 1977–1984.
  * [21] ——, “On adaptive linear–quadratic regulators,” _Automatica_ , vol. 117, p. 108982, 2020.
  * [22] ——, “Input perturbations for adaptive control and learning,” _Automatica_ , vol. 117, p. 108950, 2020.
  * [23] ——, “Optimism-based adaptive regulation of linear-quadratic systems,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 4, pp. 1802–1808, 2020.
  * [24] S. Sudhakara, A. Mahajan, A. Nayyar, and Y. Ouyang, “Scalable regret for learning to control network-coupled subsystems with unknown dynamics,” _arXiv preprint arXiv:2108.07970_ , 2021.
  * [25] P. Mandl, “Consistency of estimators in controlled systems,” in _Stochastic Differential Systems_. Springer, 1989, pp. 227–234.
  * [26] T. E. Duncan and B. Pasik-Duncan, “Adaptive control of continuous-time linear stochastic systems,” _Mathematics of Control, signals and systems_ , vol. 3, no. 1, pp. 45–60, 1990.
  * [27] P. Caines, “Continuous time stochastic adaptive control: non-explosion, $\varepsilon$-consistency and stability,” _Systems & control letters_, vol. 19, no. 3, pp. 169–176, 1992.
  * [28] T. E. Duncan, L. Guo, and B. Pasik-Duncan, “Adaptive continuous-time linear quadratic gaussian control,” _IEEE Transactions on automatic control_ , vol. 44, no. 9, pp. 1653–1662, 1999.
  * [29] P. E. Caines and D. Levanony, “Stochastic $\varepsilon$-optimal linear quadratic adaptation: An alternating controls policy,” _SIAM Journal on Control and Optimization_ , vol. 57, no. 2, pp. 1094–1126, 2019.
  * [30] S. A. A. Rizvi and Z. Lin, “Output feedback reinforcement learning control for the continuous-time linear quadratic regulator problem,” in _2018 Annual American Control Conference (ACC)_. IEEE, 2018, pp. 3417–3422.
  * [31] K. Doya, “Reinforcement learning in continuous time and space,” _Neural computation_ , vol. 12, no. 1, pp. 219–245, 2000.
  * [32] H. Wang, T. Zariphopoulou, and X. Y. Zhou, “Reinforcement learning in continuous time and space: A stochastic control approach.” _J. Mach. Learn. Res._ , vol. 21, pp. 198–1, 2020.
  * [33] Z.-P. Jiang, T. Bian, and W. Gao, “Learning-based control: A tutorial and some recent results,” _Foundations and Trends® in Systems and Control_ , vol. 8, no. 3, 2020.
  * [34] M. Basei, X. Guo, A. Hu, and Y. Zhang, “Logarithmic regret for episodic continuous-time linear-quadratic reinforcement learning over a finite-time horizon,” _Available at SSRN 3848428_ , 2021.
  * [35] L. Szpruch, T. Treetanthiploet, and Y. Zhang, “Exploration-exploitation trade-off for continuous-time episodic reinforcement learning with linear-convex models,” _arXiv preprint arXiv:2112.10264_ , 2021.
  * [36] M. K. S. Faradonbeh and M. S. S. Faradonbeh, “Bayesian algorithms learn to stabilize unknown continuous-time systems,” _arXiv preprint arXiv:2112.15094_ , 2021.
  * [37] I. Karatzas and S. Shreve, _Brownian motion and stochastic calculus_. Springer Science & Business Media, 2012, vol. 113.
  * [38] G. Chen, G. Chen, and S.-H. Hsu, _Linear stochastic control systems_. CRC press, 1995, vol. 3.
  * [39] J. Yong and X. Y. Zhou, _Stochastic controls: Hamiltonian systems and HJB equations_. Springer Science & Business Media, 1999, vol. 43.
  * [40] H. Pham, _Continuous-time stochastic control and optimization with financial applications_. Springer Science & Business Media, 2009, vol. 61.
  * [41] S. P. Bhattacharyya and L. H. Keel, _Linear Multivariable Control Systems_. Cambridge University Press, 2022\.
  * [42] J. T. Bosworth, _Linearized aerodynamic and control law models of the X-29A airplane and comparison with flight data_. National Aeronautics and Space Administration, Office of Management …, 1992, vol. 4356.
  * [43] T. Ishihara, H.-J. Guo, and H. Takeda, “A design of discrete-time integral controllers with computation delays via loop transfer recovery,” _Automatica_ , vol. 28, no. 3, pp. 599–603, 1992.
  * [44] P. Cheridito, H. M. Soner, and N. Touzi, “Small time path behavior of double stochastic integrals and applications to stochastic control,” _The Annals of Applied Probability_ , vol. 15, no. 4, pp. 2472–2495, 2005.
  * [45] B. Laurent and P. Massart, “Adaptive estimation of a quadratic functional by model selection,” _Annals of Statistics_ , pp. 1302–1338, 2000.
  * [46] J. A. Tropp, “User-friendly tail bounds for sums of random matrices,” _Foundations of computational mathematics_ , vol. 12, no. 4, pp. 389–434, 2012.
  * [47] T. Zhou, J. L. Dickson, and J. Geoffrey Chase, “Autoregressive modeling of drift and random error to characterize a continuous intravascular glucose monitoring sensor,” _Journal of Diabetes Science and Technology_ , vol. 12, no. 1, pp. 90–104, 2018.
  * [48] P. Hartman and A. Wintner, “The spectra of toeplitz’s matrices,” _American Journal of Mathematics_ , vol. 76, no. 4, pp. 867–882, 1954.
  * [49] L. Reichel and L. N. Trefethen, “Eigenvalues and pseudo-eigenvalues of toeplitz matrices,” _Linear algebra and its applications_ , vol. 162, pp. 153–185, 1992.
  * [50] R. Gondhalekar, E. Dassau, and F. J. Doyle III, “Periodic zone-mpc with asymmetric costs for outpatient-ready safety of an artificial pancreas to treat type 1 diabetes,” _Automatica_ , vol. 71, pp. 237–246, 2016.

## Organization of Appendices

This paper has four appendices. First, we prove Theorem 1 in Appendix A,
together with multiple lemmas that the proof of the theorem relies on, and
their statements and proofs are provided in Appendix A as well. Similarly, the
proof for Theorem 2 together with intermediate steps, all are presented in
Appendix B. Then, Appendix C consists of statements and proofs of other
results that are used for establishing both theorems. Finally, empirical
simulations beyond those presented in Section 6 are presented in Appendix D.

##  Appendix A Proof of Theorem 1

For analyzing the estimation error, we establish Lemma 4, which under the
condition $\boldsymbol{\kappa_{\mathrm{}}}\gtrsim\boldsymbol{\tau}^{2}$
provides that with probability at least $1-\delta$,

${\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}\lesssim\frac{p\left(p+q\right)^{1/2}}{\boldsymbol{\tau}^{1/2}}\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\vee\sigma_{w}^{2}}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}}\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right)^{3}\log\left(\frac{pq\boldsymbol{\kappa_{\mathrm{}}}}{\delta}\right).$

Note that in the proof of Lemma 4, results of Lemmas 1, 2, and 3 are used.

Therefore, Lemma 12 implies that for solutions of (5), with probability at
least $1-\delta$, it holds that

${\left|\kern-1.29167pt\left|{P}\left(\boldsymbol{\widehat{\theta}}\right)-{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}\lesssim\frac{p\left(p+q\right)^{1/2}}{\boldsymbol{\tau}^{1/2}}\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\vee\sigma_{w}^{2}}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}}\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right)^{3}\log\left(\frac{pq\boldsymbol{\kappa_{\mathrm{}}}}{\delta}\right).$

Note that we get the same expression as the right-hand-side above, as an upper
bound for
${\left|\kern-1.29167pt\left|\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right)-B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}$.
So, letting

$\overline{A}=A_{0}-B_{0}Q_{u}^{-1}\widehat{B}^{\top}{P}\left(\boldsymbol{\widehat{\theta}}\right),~{}~{}~{}~{}~{}~{}~{}~{}\overline{A}_{0}=A_{0}-B_{0}Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right),$

we obtain

${\left|\kern-1.29167pt\left|\overline{A}-\overline{A}_{0}\right|\kern-1.29167pt\right|}\lesssim\sqrt{\frac{p^{2}q}{\boldsymbol{\tau}}}\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\vee\sigma_{w}^{2}}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}}\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right)^{3}\log\left(\frac{pq\boldsymbol{\kappa_{\mathrm{}}}}{\delta}\right),$
(13)

with probability at least $1-\delta$.

Next, to consider the effect of the above errors on the eigenvalues of
$\overline{A}$, we compare them to that of $\overline{A}_{0}$. Note that real-
parts of all eigenvalues of $\overline{A}_{0}$ are at most
$-\boldsymbol{\zeta_{0}}$, as defined in (10). So, using the result and the
notation of Lemma 5, for all eigenvalues of $\overline{A}$ being on the open
left half-plane it suffices to have
$\Delta_{\overline{A}_{0}}\left(\overline{A}-\overline{A}_{0}\right)\leq\boldsymbol{\zeta_{0}}$.
Also, in lights of Lemma 5, suppose that $\boldsymbol{r}$ is the size of the
largest block in the Jordan block-diagonalization of $\overline{A}_{0}$. So,
(36) implies that if

${\left|\kern-1.29167pt\left|\overline{A}-\overline{A}_{0}\right|\kern-1.29167pt\right|}\lesssim\boldsymbol{r}^{-1/2}\left(1\wedge\boldsymbol{\zeta_{0}}^{\boldsymbol{r}}\right),$

then Algorithm 1 successfully stabilizes the diffusion process in (1). Thus,
(13) shows that the failure probability of Algorithm 1;
$\mathbb{P}{\left(\mathcal{E}_{\boldsymbol{\tau}}\right)}$, satisfies

$\sqrt{\frac{p^{2}q}{\boldsymbol{\tau}}}\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\vee\sigma_{w}^{2}}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}}\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right)^{3}\log\left(\frac{pq\boldsymbol{\kappa_{\mathrm{}}}}{\mathbb{P}{\left(\mathcal{E}_{\boldsymbol{\tau}}\right)}}\right)\gtrsim\boldsymbol{r}^{-1/2}\left(1\wedge\boldsymbol{\zeta_{0}}^{\boldsymbol{r}}\right).$

Finally, $\boldsymbol{r}\leq p$ together with
$\log\left(pq\boldsymbol{\kappa_{\mathrm{}}}\right)\lesssim\boldsymbol{\tau}^{1/2}$,
lead to the desired result.

In the remainder of this section, technical lemmas above are stated and their
proofs will be provided.

###  A.1 Bounding cross products of state and randomization

###### Definition 2

For a set $\mathcal{S}$, let $\mathds{1}{\left\\{\mathcal{S}\right\\}}$ be the
indicator function that is $1$ on $\mathcal{S}$, and vanishes outside of
$\mathcal{S}$.

###### Lemma 1

In Algorithm 1, for $t\geq 0$, define the piecewise-constant signal $v(t)$
below according to the randomization sequence $w_{n}$:

$v(t)=\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\mathds{1}{\left\\{\frac{n\boldsymbol{\tau}}{\boldsymbol{\kappa_{\mathrm{}}}}\leq
t<\frac{(n+1)\boldsymbol{\tau}}{\boldsymbol{\kappa_{\mathrm{}}}}\right\\}}w_{n}.$
(14)

Then, with probability at least $1-\delta$, we have

$\displaystyle{{\left|\kern-1.29167pt\left|\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{x}_{s}v(s)^{\top}}\boldsymbol{d}{s}-\frac{\boldsymbol{\tau}^{2}}{2\boldsymbol{\kappa_{\mathrm{}}}^{2}}B_{0}\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}w_{n}w_{n}^{\top}\right|\kern-1.29167pt\right|}}$
$\displaystyle\lesssim$
$\displaystyle{\left(\sigma_{w}^{2}+\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\right)\left(1+\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right)}\left(pq^{1/2}\boldsymbol{\tau}^{1/2}\log\frac{pq}{\delta}+q\left[1+\frac{\boldsymbol{\tau}^{2}}{\boldsymbol{\kappa_{\mathrm{}}}^{2}}\right]\frac{\boldsymbol{\tau}}{\boldsymbol{\kappa_{\mathrm{}}}^{1/2}}\log^{3/2}\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta}\right).$

Proof. First, after plugging the control signal $\boldsymbol{u}_{t}$ in (1)
and solving the resulting stochastic differential equation, we obtain

$\boldsymbol{x}_{t}=e^{\overline{A}t}\boldsymbol{x}_{0}+\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}+\int\limits_{0}^{t}{e^{\overline{A}(t-s)}B_{0}v(s)}\boldsymbol{d}{s}.$
(15)

This implies that

$\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{x}_{t}v(t)^{\top}}\boldsymbol{d}{t}=\Phi_{1}+\Phi_{2}+\Phi_{3},$

where

$\displaystyle\Phi_{1}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}t}\boldsymbol{x}_{0}v(t)^{\top}}\boldsymbol{d}{t}=\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}t}}\boldsymbol{d}{t}\right)\boldsymbol{x}_{0}w_{n}^{\top},$
$\displaystyle\Phi_{2}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}v(t)^{\top}}\boldsymbol{d}{t}=\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}}\boldsymbol{d}{t}\right)w_{n}^{\top},$
$\displaystyle\Phi_{3}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{\int\limits_{0}^{t}{e^{\overline{A}(t-s)}B_{0}v(s)}\boldsymbol{d}{s}v(t)^{\top}}\boldsymbol{d}{t}.$

To analyze $\Phi_{1}$, we use the fact that every entry of $\Phi_{1}$ is a
normal random variable with mean zero and variance at most

$\sigma_{w}^{2}\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}{\left|\kern-1.29167pt\left|\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}t}}\boldsymbol{d}{t}\right|\kern-1.29167pt\right|}^{2}{\left|\kern-1.50696pt\left|\boldsymbol{x}_{0}\right|\kern-1.50696pt\right|}^{2}\leq\sigma_{w}^{2}\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right)^{2}{\left|\kern-1.50696pt\left|\boldsymbol{x}_{0}\right|\kern-1.50696pt\right|}^{2}.$

Therefore, with probability at least $1-\delta$, it holds that

${\left|\kern-1.29167pt\left|\Phi_{1}\right|\kern-1.29167pt\right|}\lesssim\sigma_{w}\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right){\left|\kern-1.50696pt\left|\boldsymbol{x}_{0}\right|\kern-1.50696pt\right|}\sqrt{pq\log\left(\frac{pq}{\delta}\right)}.$
(16)

Furthermore, to study $\Phi_{2}$, Fubini Theorem [37] gives

$\displaystyle\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}}\boldsymbol{d}{t}$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\left(\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(t-n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1})}}\boldsymbol{d}{t}\right)e^{\overline{A}(n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}-s)}}\boldsymbol{d}{\mathbb{W}_{s}}$
$\displaystyle+$
$\displaystyle\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\left(\int\limits_{s}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(t-s)}}\boldsymbol{d}{t}\right)}\boldsymbol{d}{\mathbb{W}_{s}}$
$\displaystyle=$ $\displaystyle
F\sum\limits_{m=1}^{n}e^{\overline{A}(n-m)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}\int\limits_{(m-1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{m\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(m\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}-s)}}\boldsymbol{d}{\mathbb{W}_{s}}$
$\displaystyle+$
$\displaystyle\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{G_{s}}\boldsymbol{d}{\mathbb{W}_{s}}$

where the matrix
$F=\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(t-n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1})}}\boldsymbol{d}{t}$,
does not depend on $s$ or $n$, and the matrix
$G_{s}=\int\limits_{s}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(t-s)}}\boldsymbol{d}{t}$,
does not depend on $n$, since
$n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}\leq
s\leq(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}$.

So, letting $\mathrm{e}_{i},i=1,\cdots,p$, be the standard basis of the
Euclidean space, conditioned on the Wiener process
$\left\\{\mathbb{W}_{s}\right\\}_{s\geq 0}$, for every $j=1,\cdots,q$, the
coordinate $j$ of $\mathrm{e}_{i}^{\top}\Phi_{2}$ is a mean zero normal random
variable. Thus, given $\left\\{\mathbb{W}_{s}\right\\}_{s\geq 0}$, with
probability at least $1-\delta$, it holds that

$\left(\mathrm{e}_{i}^{\top}\Phi_{2}\mathrm{e}_{j}\right)^{2}\lesssim\mathrm{var\left(\mathrm{e}_{i}^{\top}\Phi_{2}\mathrm{e}_{j}\Big{|}\mathcal{F}\left(\mathbb{W}_{0:\boldsymbol{\tau}}\right)\right)}\log\frac{1}{\delta}.$

Now, to calculate the conditional variance, we can write

$\frac{\mathrm{var\left(\mathrm{e}_{i}^{\top}\Phi_{2}\mathrm{e}_{j}\Big{|}\mathcal{F}\left(\mathbb{W}_{0:\boldsymbol{\tau}}\right)\right)}}{\sigma_{w}^{2}}=\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left[\mathrm{e}_{i}^{\top}\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}}\boldsymbol{d}{t}\right]^{2}\lesssim\sum\limits_{n=1}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left[\left(\sum\limits_{m=1}^{n}\beta_{m,n}\right)^{2}+\alpha_{n}^{2}\right],$

where

$\displaystyle\beta_{m,n}$ $\displaystyle=$
$\displaystyle\mathrm{e}_{i}^{\top}Fe^{\overline{A}(n-m)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}\int\limits_{(m-1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{m\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(m\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}-s)}}\boldsymbol{d}{\mathbb{W}_{s}},$
$\displaystyle\alpha_{n}$ $\displaystyle=$
$\displaystyle\mathrm{e}_{i}^{\top}\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{G_{s}}\boldsymbol{d}{\mathbb{W}_{s}}.$

To proceed, define the matrix $H=\left[H_{n,m}\right]$, where for $1\leq
m,n\leq\boldsymbol{\kappa_{\mathrm{}}}-1$, every block
$H_{n,m}\in\mathbb{R}^{1\times p}$ is

$H_{n,m}=\mathrm{e}_{i}^{\top}Fe^{\overline{A}(n-m)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}},$

for $m\leq n$, and is $0$ for $m>n$. Then, denote

$\Gamma=\begin{bmatrix}\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}-s)}}\boldsymbol{d}{\mathbb{W}_{s}}\\\
\int\limits_{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{2\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(2\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}-s)}}\boldsymbol{d}{\mathbb{W}_{s}}\\\
\vdots\\\
\int\limits_{(\boldsymbol{\kappa_{\mathrm{}}}-1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{\boldsymbol{\tau}}{e^{\overline{A}(\boldsymbol{\tau}-s)}}\boldsymbol{d}{\mathbb{W}_{s}}\end{bmatrix}\in\mathbb{R}^{p(\boldsymbol{\kappa_{\mathrm{}}}-1)\times
1},$

to get

$\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(\sum\limits_{m=1}^{n}\beta_{m,n}^{2}\right)={\left|\kern-1.50696pt\left|H\Gamma\right|\kern-1.50696pt\right|}^{2}\leq\overline{\lambda}\left(H^{\top}H\right){\left|\kern-1.50696pt\left|\Gamma\right|\kern-1.50696pt\right|}^{2}.$

Now, for the matrix $H$, we have [48, 49]:

$\overline{\lambda}\left(H^{\top}H\right)\lesssim\left(\sum\limits_{n=1}^{\boldsymbol{\kappa_{\mathrm{}}}-1}{\left|\kern-1.50696pt\left|H_{n,1}\right|\kern-1.50696pt\right|}\right)^{2}\lesssim\left(\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}\sum\limits_{n=1}^{\boldsymbol{\kappa_{\mathrm{}}}}e^{\overline{A}n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}\right)^{2}\lesssim\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right)^{2}.$

Note that thanks to the independent increments of the Wiener process, the
blocks of $\Gamma$ are statistically independent. Further, by Ito Isometry
[37], every block of $\Gamma$ is a mean-zero normally distributed vector with
the covariance matrix

$\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}(\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}-s)}\Sigma_{\mathbb{W}}e^{\overline{A}^{\top}(\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}-s)}}\boldsymbol{d}{s}.$

So, according to the exponential inequalities for quadratic forms of normally
distributed random variables [45], it holds with probability at least
$1-\delta$, that

${\left|\kern-1.50696pt\left|\Gamma\right|\kern-1.50696pt\right|}^{2}\lesssim
p\boldsymbol{\kappa_{\mathrm{}}}\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\left(\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}\right)\log\frac{1}{\delta}.$

Thus, with probability at least $1-\delta$, we have

$\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(\sum\limits_{m=1}^{n}\beta_{m,n}^{2}\right)\lesssim\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right)^{2}p\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\boldsymbol{\tau}\log\frac{1}{\delta}.$

Similarly, the bound above can be shown for
$\sum\limits_{n=1}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\alpha_{n}^{2}$. Hence,
we obtain the corresponding high probability bound for a single entry
$\mathrm{e}_{i}^{\top}\Phi_{2}\mathrm{e}_{j}$ of $\Phi_{2}$, which together
with a union bound, implies that

${\left|\kern-1.29167pt\left|\Phi_{2}\right|\kern-1.29167pt\right|}\lesssim\sigma_{w}pq^{1/2}\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right)\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)^{1/2}\boldsymbol{\tau}^{1/2}\log\left(\frac{pq}{\delta}\right),$
(17)

with probability at least $1-\delta$.

Next, according to Fubini Theorem, $\Phi_{3}$ can also be written as

$\displaystyle\Phi_{3}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{\int\limits_{0}^{s}{e^{\overline{A}(s-t)}B_{0}v(t)v(s)^{\top}}\boldsymbol{d}{t}}\boldsymbol{d}{s}=\int\limits_{0}^{\boldsymbol{\tau}}{\int\limits_{t}^{\boldsymbol{\tau}}{e^{\overline{A}(s-t)}B_{0}v(t)v(s)^{\top}}\boldsymbol{d}{s}}\boldsymbol{d}{t}.$

Thus, we have

$2\Phi_{3}=\int\limits_{0}^{\boldsymbol{\tau}}{\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}\left|t-s\right|}B_{0}v(t\wedge
s)v(s\vee t)^{\top}}\boldsymbol{d}{t}}\boldsymbol{d}{s}.$

Recall that the signal $v(t)$ in (14) is piecewise-constant, with values
determined by the randomization sequence $w_{n}$. So, the above double
integral can be written as a double sum

$\displaystyle 2\Phi_{3}$ $\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\sum\limits_{m=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{m\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(m+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}\left|t-s\right|}}\boldsymbol{d}{s}}\boldsymbol{d}{t}\right)B_{0}w_{m\wedge
n}w_{m\vee n}^{\top}$ $\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\sum\limits_{m=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(e^{\overline{A}\left|m-n\right|\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}\left|t-s\right|}}\boldsymbol{d}{s}}\boldsymbol{d}{t}\right)B_{0}w_{m\wedge
n}w_{m\vee n}^{\top}.$

Thus, we have

$2\Phi_{3}-\frac{\boldsymbol{\tau}^{2}}{\boldsymbol{\kappa_{\mathrm{}}}^{2}}B_{0}\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}w_{n}w_{n}^{\top}=\Phi_{4}+\Phi_{5},$
(18)

for

$\displaystyle\Phi_{4}$ $\displaystyle=$
$\displaystyle\left(\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}\left|t-s\right|}}\boldsymbol{d}{s}}\boldsymbol{d}{t}-\boldsymbol{\tau}^{2}\boldsymbol{\kappa_{\mathrm{}}}^{-2}I_{q}\right)B_{0}\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}w_{n}w_{n}^{\top},$
$\displaystyle\Phi_{5}$ $\displaystyle=$ $\displaystyle
2\left(\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}\left|t-s\right|}}\boldsymbol{d}{s}}\boldsymbol{d}{t}\right)\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\sum\limits_{m=n+1}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(e^{\overline{A}\left(m-n\right)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}B_{0}w_{n}w_{m}^{\top}\right).$

To proceed, we use the following concentration inequality for random matrices
with martingale difference structures, titled as Matrix Azuma inequality [46].

###### Theorem 3

Let $\left\\{\Psi_{n}\right\\}_{n=1}^{k}$ be a $d_{1}\times d_{2}$ martingale
difference sequence. That is, for some filtration
$\left\\{\mathcal{F}_{n}\right\\}_{n=0}^{k}$, the matrix $\Psi_{n}$ is
$\mathcal{F}_{n}$-measurable, and
$\mathbb{E}\left[\Psi_{n}\Big{|}\mathcal{F}_{n-1}\right]=0$. Suppose that
${\left|\kern-1.29167pt\left|\Psi_{n}\right|\kern-1.29167pt\right|}\leq\sigma_{n}$,
for some fixed sequence $\left\\{\sigma_{n}\right\\}_{n=1}^{k}$. Then, with
probability at least $1-\delta$, we have

${\left|\kern-1.29167pt\left|\sum\limits_{n=1}^{k}\Psi_{n}\right|\kern-1.29167pt\right|}^{2}\lesssim\left(\sum\limits_{n=1}^{k}\sigma_{n}^{2}\right)\log\frac{d_{1}+d_{2}}{\delta}.$

So, to study $\Phi_{4}$, we apply Theorem 3 to the random matrices
$\Psi_{n}=w_{n}w_{n}^{\top}-\sigma_{w}^{2}I_{q}$, using the trivial filtration
and the high probability upper-bounds for
${\left|\kern-1.29167pt\left|\Psi_{n}\right|\kern-1.29167pt\right|}\leq{\left|\kern-1.50696pt\left|w_{n}\right|\kern-1.50696pt\right|}^{2}+\sigma_{w}^{2}$;

${\left|\kern-1.29167pt\left|\Psi_{n}\right|\kern-1.29167pt\right|}\leq\sigma_{n}=\sigma_{w}^{2}\left(1+q\log\frac{q\boldsymbol{\kappa_{\mathrm{}}}}{\delta}\right),$

as well as the fact

${\left|\kern-1.29167pt\left|\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\left(e^{\overline{A}\left|t-s\right|}-I_{q}\right)}\boldsymbol{d}{s}}\boldsymbol{d}{t}\right|\kern-1.29167pt\right|}\lesssim\boldsymbol{\tau}^{3}\boldsymbol{\kappa_{\mathrm{}}}^{-3},$

to obtain the following bound, which holds with probability at least
$1-\delta$:

${\left|\kern-1.29167pt\left|\Phi_{4}\right|\kern-1.29167pt\right|}\lesssim{\left|\kern-1.29167pt\left|B_{0}\right|\kern-1.29167pt\right|}\sigma_{w}^{2}\boldsymbol{\tau}^{3}\boldsymbol{\kappa_{\mathrm{}}}^{-2}\left(1+\frac{q}{\boldsymbol{\kappa_{\mathrm{}}}^{1/2}}\log^{3/2}\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta}\right).$
(19)

On the other hand, to establish an upper-bound for $\Phi_{5}$, consider the
random matrices

$\Psi_{n}=\sum\limits_{m=n+1}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\left(e^{\overline{A}\left(m-n\right)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}B_{0}w_{n}w_{m}^{\top}\right),$

subject to the natural filtration they generate, and apply Theorem 3, using
the bounds

${\left|\kern-1.29167pt\left|\Psi_{n}\right|\kern-1.29167pt\right|}\leq\sigma_{n}\lesssim\boldsymbol{\tau}^{-1}\boldsymbol{\kappa_{\mathrm{}}}\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right){\left|\kern-1.29167pt\left|B_{0}\right|\kern-1.29167pt\right|}\sigma_{w}^{2}q\log\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta},$

together with

${\left|\kern-1.29167pt\left|\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{\int\limits_{0}^{\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}\left|t-s\right|}}\boldsymbol{d}{s}}\boldsymbol{d}{t}\right|\kern-1.29167pt\right|}\lesssim\boldsymbol{\tau}^{2}\boldsymbol{\kappa_{\mathrm{}}}^{-2}.$

Therefore, Theorem 3 indicates that with probability at least $1-\delta$, it
holds that

$\Phi_{5}\lesssim\frac{\boldsymbol{\tau}}{\boldsymbol{\kappa_{\mathrm{}}}^{1/2}}\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right){\left|\kern-1.29167pt\left|B_{0}\right|\kern-1.29167pt\right|}\sigma_{w}^{2}q\log^{3/2}\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta}.$
(20)

Finally, put (16), (17), (18), (19), and (20) together, to get the desired
result.

$\blacksquare$

###  A.2 Bounding cross products of state and Wiener process

###### Lemma 2

In Algorithm 1, with probability at least $1-\delta$, we have

${\left|\kern-1.29167pt\left|\int\limits_{0}^{t}{\boldsymbol{x}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}\right|\kern-1.29167pt\right|}\lesssim\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right)\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}\right)p\left(p+q\right)^{1/2}\boldsymbol{\tau}^{1/2}\log\left(\frac{pq}{\delta}\right).$

Proof. First, according to (15), we can write

$\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{x}_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}=\Phi_{1}+\Phi_{2}+\Phi_{3},$

where

$\displaystyle\Phi_{1}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}t}\boldsymbol{x}_{0}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}},$
(21) $\displaystyle\Phi_{2}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{\int\limits_{0}^{t}{e^{\overline{A}(t-s)}B_{0}v(s)}\boldsymbol{d}{s}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}},$
(22) $\displaystyle\Phi_{3}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}.$
(23)

Now, according to Ito Isometry [37], similar to (16), we have

${\left|\kern-1.29167pt\left|\Phi_{1}\right|\kern-1.29167pt\right|}\lesssim\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)^{1/2}\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right){\left|\kern-1.50696pt\left|\boldsymbol{x}_{0}\right|\kern-1.50696pt\right|}\sqrt{pq\log\left(\frac{pq}{\delta}\right)},$
(24)

with probability at least $1-\delta$. Moreover, in a procedure similar to the
one that lead to (17), one can show that with probability at least $1-\delta$,
it holds that

${\left|\kern-1.29167pt\left|\Phi_{2}\right|\kern-1.29167pt\right|}\lesssim\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}t}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{t}\right)\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)^{1/2}\sigma_{w}pq^{1/2}\boldsymbol{\tau}^{1/2}\log\left(\frac{pq}{\delta}\right).$
(25)

Therefore, we need to find a similar upper-bound for $\Phi_{3}$. To that end,
Ito formula provides

$\boldsymbol{d}\left(e^{-\overline{A}s}\mathbb{W}_{s}\right)=-\overline{A}e^{-\overline{A}s}\mathbb{W}_{s}\boldsymbol{d}s+e^{-\overline{A}s}\boldsymbol{d}\mathbb{W}_{s}.$

Therefore, integration gives

$\int\limits_{0}^{t}{e^{-\overline{A}s}}\boldsymbol{d}{\mathbb{W}_{s}}=e^{-\overline{A}t}\mathbb{W}_{t}+\overline{A}\int\limits_{0}^{t}{e^{-\overline{A}s}\mathbb{W}_{s}}\boldsymbol{d}{s},$

which after rearranging and letting
$\Psi_{t}=\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}$,
leads to

$\Psi_{t}\mathbb{W}_{t}^{\top}=\left(\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}\right){\mathbb{W}_{t}^{\top}}={\mathbb{W}_{t}}{\mathbb{W}_{t}^{\top}}+\overline{A}{\left(\int\limits_{0}^{t}{e^{\overline{A}(t-s)}\mathbb{W}_{s}}\boldsymbol{d}{s}\right)}{\mathbb{W}_{t}^{\top}}.$

Now, since $\boldsymbol{d}\Psi_{t}=\boldsymbol{d}\mathbb{W}_{t}$, Ito Isometry
[37] implies that
$\boldsymbol{d}\Psi_{t}\boldsymbol{d}\mathbb{W}_{t}^{\top}=\Sigma_{\mathbb{W}}\boldsymbol{d}t$.
So, apply integration by part and use the above equation to get

$\displaystyle\Phi_{3}=\int\limits_{0}^{\boldsymbol{\tau}}{\Psi_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{}\boldsymbol{d}{\left(\Psi_{t}\mathbb{W}_{t}^{\top}\right)}-\left(\int\limits_{0}^{\boldsymbol{\tau}}{\mathbb{W}_{t}}\boldsymbol{d}{\Psi_{t}^{\top}}\right)^{\top}-\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{d}\Psi_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}$
$\displaystyle=$
$\displaystyle{\Psi_{\boldsymbol{\tau}}}{\mathbb{W}_{\boldsymbol{\tau}}^{\top}}-\left(\int\limits_{0}^{\boldsymbol{\tau}}{\mathbb{W}_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}\right)^{\top}-\Sigma_{\mathbb{W}}\boldsymbol{\tau}$
$\displaystyle=$
$\displaystyle{\mathbb{W}_{\boldsymbol{\tau}}}{\mathbb{W}_{\boldsymbol{\tau}}^{\top}}+\overline{A}{\left(\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}(\boldsymbol{\tau}-s)}\mathbb{W}_{s}}\boldsymbol{d}{s}\right)}{\mathbb{W}_{\boldsymbol{\tau}}^{\top}}-\left(\int\limits_{0}^{\boldsymbol{\tau}}{\mathbb{W}_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}\right)^{\top}-\Sigma_{\mathbb{W}}\boldsymbol{\tau}.$

Therefore, every entry of $\Phi_{3}$ is a quadratic function of the normally
distributed random vectors
$\mathbb{W}_{\boldsymbol{\tau}},\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}(\boldsymbol{\tau}-s)}\mathbb{W}_{s}}\boldsymbol{d}{s}$.
Note that we used the fact that

$\mathbb{W}_{\boldsymbol{\tau}}\mathbb{W}_{\boldsymbol{\tau}}^{\top}=\int\limits_{0}^{\boldsymbol{\tau}}{}\boldsymbol{d}{\left(\mathbb{W}_{t}\mathbb{W}_{t}^{\top}\right)}=\left(\int\limits_{0}^{\boldsymbol{\tau}}{\mathbb{W}_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}\right)^{\top}+\left(\int\limits_{0}^{\boldsymbol{\tau}}{\mathbb{W}_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}\right)+\Sigma_{\mathbb{W}}\boldsymbol{\tau}.$

Thus, exponential inequalities for quadratic forms of normal random vectors
[45] imply that for all $i,j=1,\cdots,p$, it holds that

$\left(\mathrm{e}_{i}^{\top}\Phi_{3}\mathrm{e}_{j}\right)^{2}\lesssim
p\mathbb{E}\left[\left(\mathrm{e}_{i}^{\top}\Phi_{3}\mathrm{e}_{j}\right)^{2}\right]\log^{2}\frac{1}{\delta},$
(26)

since $\mathbb{E}\left[\mathrm{e}_{i}^{\top}\Phi_{3}\mathrm{e}_{j}\right]=0$.
So, it suffices to find the expectation in (26). For that purpose, we use Ito
Isometry [37] to obtain:

$\displaystyle\mathbb{E}\left[\left(\mathrm{e}_{i}^{\top}\Phi_{3}\mathrm{e}_{j}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\int\limits_{0}^{\boldsymbol{\tau}}{\mathrm{e}_{i}^{\top}\Psi_{t}\mathrm{e}_{j}^{\top}\Sigma_{\mathbb{W}}^{1/2}}\boldsymbol{d}{\left(\Sigma_{\mathbb{W}}^{-1/2}\mathbb{W}_{t}\right)}\right)^{2}\right]=\mathbb{E}\left[\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.50696pt\left|\mathrm{e}_{i}^{\top}\Psi_{t}\Sigma_{\mathbb{W}}^{1/2}\mathrm{e}_{j}\right|\kern-1.50696pt\right|}^{2}}\boldsymbol{d}{t}\right]$
$\displaystyle\leq$
$\displaystyle\mathrm{e}_{j}^{\top}\Sigma_{\mathbb{W}}\mathrm{e}_{j}\mathbb{E}\left[\int\limits_{0}^{\boldsymbol{\tau}}{\left(\mathrm{e}_{i}^{\top}\Psi_{t}\right)^{2}}\boldsymbol{d}{t}\right]=\mathrm{e}_{j}^{\top}\Sigma_{\mathbb{W}}\mathrm{e}_{j}\mathbb{E}\left[\int\limits_{0}^{\boldsymbol{\tau}}{\left(\mathrm{e}_{i}^{\top}\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}\right)^{2}}\boldsymbol{d}{t}\right].$

To proceed with the above expression, apply Fubini Theorem [37] to interchange
the expected value with the integral, and then use Ito Isometry again:

$\displaystyle\mathbb{E}\left[\int\limits_{0}^{\boldsymbol{\tau}}{\left(\mathrm{e}_{i}^{\top}\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}\right)^{2}}\boldsymbol{d}{t}\right]$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{\mathbb{E}\left[\left(\mathrm{e}_{i}^{\top}\int\limits_{0}^{t}{e^{\overline{A}(t-s)}\Sigma_{\mathbb{W}}^{1/2}}\boldsymbol{d}{\left(\Sigma_{\mathbb{W}}^{-1/2}\mathbb{W}_{s}\right)}\right)^{2}\right]}\boldsymbol{d}{t}$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{\boldsymbol{\tau}}{\mathrm{e}_{i}^{\top}\left(\int\limits_{0}^{t}{e^{\overline{A}(t-s)}\Sigma_{\mathbb{W}}e^{\overline{A}^{\top}(t-s)}}\boldsymbol{d}{s}\right)\mathrm{e}_{i}}\boldsymbol{d}{t}$
$\displaystyle\leq$
$\displaystyle\mathrm{e}_{i}^{\top}\left(\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}s}\Sigma_{\mathbb{W}}e^{\overline{A}^{\top}s}}\boldsymbol{d}{s}\right)\mathrm{e}_{i}\boldsymbol{\tau}.$

Therefore, (26) yields to

$\displaystyle{\left|\kern-1.29167pt\left|\Phi_{3}\right|\kern-1.29167pt\right|}^{2}\leq\sum\limits_{i,j=1}^{p}\left(\mathrm{e}_{i}^{\top}\Phi_{3}\mathrm{e}_{j}\right)^{2}$
$\displaystyle\lesssim$
$\displaystyle\sum\limits_{i,j=1}^{p}\left[\mathrm{e}_{j}^{\top}\Sigma_{\mathbb{W}}\mathrm{e}_{j}\mathrm{e}_{i}^{\top}\left(\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}s}\Sigma_{\mathbb{W}}e^{\overline{A}^{\top}s}}\boldsymbol{d}{s}\right)\mathrm{e}_{i}\right]\boldsymbol{\tau}p\log^{2}\frac{p}{\delta}$
(27) $\displaystyle=$
$\displaystyle\boldsymbol{\mathrm{tr}}\left(\Sigma_{\mathbb{W}}\right)\boldsymbol{\mathrm{tr}}\left(\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}s}\Sigma_{\mathbb{W}}e^{\overline{A}^{\top}s}}\boldsymbol{d}{s}\right)p\boldsymbol{\tau}\log^{2}\frac{p}{\delta}$
$\displaystyle\lesssim$
$\displaystyle\boldsymbol{\mathrm{tr}}\left(\Sigma_{\mathbb{W}}\right)^{2}\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|e^{\overline{A}s}\right|\kern-1.29167pt\right|}}\boldsymbol{d}{s}\right)^{2}p\boldsymbol{\tau}\log^{2}\frac{p}{\delta}.$

Finally, putting (24), (25), and (27) together, we obtain the desired result.

$\blacksquare$

###  A.3 Concentration of normal posterior distribution in Algorithm 1

###### Lemma 3

In Algorithm 1, letting $\overline{A}=A_{0}+B_{0}K$, suppose that

$\displaystyle\boldsymbol{\tau}$ $\displaystyle\gtrsim$
$\displaystyle\left(\int\limits_{0}^{\boldsymbol{\tau}}{{\left|\kern-1.29167pt\left|\exp(\overline{A}s)\right|\kern-1.29167pt\right|}^{2}}\boldsymbol{d}{s}\right)\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}{\left|\kern-1.29167pt\left|B_{0}\right|\kern-1.29167pt\right|}^{2}\right)(p+q)\log\frac{1}{\delta},$
(28) $\displaystyle\frac{\boldsymbol{\kappa_{\mathrm{}}}}{\boldsymbol{\tau}}$
$\displaystyle\gtrsim$
$\displaystyle\frac{\sigma_{w}^{2}}{\sigma_{w}^{2}\wedge\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)}~{}~{}{\left|\kern-1.29167pt\left|B_{0}\right|\kern-1.29167pt\right|}\left(1\vee{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right)q\log\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta}.$
(29)

Then, for the matrix $\widehat{\Sigma}_{\boldsymbol{\tau}}$ in (8), with
probability at least $1-\delta$ we have

$\displaystyle\underline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}}\right)\gtrsim\boldsymbol{\tau}\left(\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}\right)\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}^{2}\right)^{-1}.$

Proof. First, we can write the control action in (7) as
$\boldsymbol{u}_{t}=K\boldsymbol{x}_{t}+v(t)$, for the piecewise-constant
signal $v(t)$ in (14). Then, the dynamics in (1) provides

$\boldsymbol{d}\boldsymbol{x}_{t}=\left(\overline{A}\boldsymbol{x}_{t}+B_{0}v(t)\right)\boldsymbol{d}t+\boldsymbol{d}\mathbb{W}_{t}.$

Therefore, similar to (15), one can solve the above stochastic differential
equation to get

$\boldsymbol{x}_{t}=e^{\overline{A}t}\boldsymbol{x}_{0}+\int\limits_{0}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}+\int\limits_{0}^{t}{e^{\overline{A}(t-s)}B_{0}v(s)}\boldsymbol{d}{s}.$

So, using the exponential inequalities for quadratic forms [45], with
probability at least $1-\delta$, it holds that

${\left|\kern-1.50696pt\left|\boldsymbol{x}_{\boldsymbol{\tau}}-e^{\overline{A}\boldsymbol{\tau}}\boldsymbol{x}_{0}\right|\kern-1.50696pt\right|}^{2}\lesssim\overline{\lambda}\left(\int\limits_{0}^{\boldsymbol{\tau}}{e^{\overline{A}s}\Sigma_{\mathbb{W}}e^{\overline{A}^{\top}s}}\boldsymbol{d}{s}+\sigma_{w}^{2}\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}J_{n}B_{0}B_{0}^{\top}J_{n}^{\top}\right)\left(p+p^{1/2}\log\frac{1}{\delta}\right),$
(30)

where

$J_{n}=\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{e^{\overline{A}s}}\boldsymbol{d}{s}.$

Furthermore, an application of Ito calculus [37] leads to
$\boldsymbol{d}\boldsymbol{x}_{t}\boldsymbol{d}\boldsymbol{x}_{t}^{\top}=\boldsymbol{d}\mathbb{W}_{t}\boldsymbol{d}\mathbb{W}_{t}^{\top}=\Sigma_{\mathbb{W}}\boldsymbol{d}t$.
Now, by defining the matrix valued processes

$\Phi_{t}=\int\limits_{0}^{t}{\boldsymbol{x}_{s}\boldsymbol{x}_{s}^{\top}}\boldsymbol{d}{s},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}M_{t}=\int\limits_{0}^{t}{\boldsymbol{x}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}+\int\limits_{0}^{t}{\boldsymbol{x}_{s}v(s)^{\top}B_{0}^{\top}}\boldsymbol{d}{s},$

we obtain

$\displaystyle\boldsymbol{d}\left(\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{\top}\right)$
$\displaystyle=$
$\displaystyle\boldsymbol{x}_{t}\boldsymbol{d}\boldsymbol{x}_{t}^{\top}+\boldsymbol{d}\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{\top}+\boldsymbol{d}\boldsymbol{x}_{t}\boldsymbol{d}\boldsymbol{x}_{t}^{\top}$
$\displaystyle=$
$\displaystyle\boldsymbol{x}_{t}\left(\left(\overline{A}\boldsymbol{x}_{t}+B_{0}v(t)\right)\boldsymbol{d}t+\boldsymbol{d}\mathbb{W}_{t}\right)^{\top}$
$\displaystyle+$
$\displaystyle\left(\left(\overline{A}\boldsymbol{x}_{t}+B_{0}v(t)\right)\boldsymbol{d}t+\boldsymbol{d}\mathbb{W}_{t}\right)\boldsymbol{x}_{t}^{\top}+\Sigma_{\mathbb{W}}\boldsymbol{d}t$
$\displaystyle=$
$\displaystyle\boldsymbol{d}\Phi_{t}\overline{A}^{\top}+\overline{A}\boldsymbol{d}\Phi_{t}+\boldsymbol{d}M_{t}+\boldsymbol{d}M_{t}^{\top}+\Sigma_{\mathbb{W}}\boldsymbol{d}t.$

Thus, after integrating both sides of the above equality, we obtain

$\Phi_{t}\overline{A}^{\top}+\overline{A}\Phi_{t}+M_{t}+M_{t}^{\top}+t\Sigma_{\mathbb{W}}+\boldsymbol{x}_{0}\boldsymbol{x}_{0}^{\top}-\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{\top}=0.$

Because all eigenvalues of $\overline{A}$ are in the open left half-plane, we
can solve the above equation for $\Phi_{t}$, to get

$\Phi_{t}=\int\limits_{0}^{\infty}{\exp\left(\overline{A}s\right)\left[M_{t}+M_{t}^{\top}+t\Sigma_{\mathbb{W}}+\boldsymbol{x}_{0}\boldsymbol{x}_{0}^{\top}-\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{\top}\right]\exp\left(\overline{A}^{\top}s\right)}\boldsymbol{d}{s}.$
(31)

Next, putting Lemma 1, Lemma 2, and (30) together, as long as (28) holds, with
probability at least $1-\delta$ we have

$\underline{\lambda}\left(M_{\boldsymbol{\tau}}+M_{\boldsymbol{\tau}}^{\top}+\boldsymbol{\tau}\Sigma_{\mathbb{W}}+\boldsymbol{x}_{0}\boldsymbol{x}_{0}^{\top}-\boldsymbol{x}_{\boldsymbol{\tau}}\boldsymbol{x}_{\boldsymbol{\tau}}^{\top}\right)\gtrsim\boldsymbol{\tau}\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right).$

Thus, (31) implies that
$\underline{\lambda}\left(\Phi_{\boldsymbol{\tau}}\right)\gtrsim\boldsymbol{\tau}\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)$.
To proceed, consider the matrix $\widehat{\Sigma}_{\boldsymbol{\tau}}$ in (8),
which comprises two signals $\boldsymbol{x}_{t},v(t)$. The empirical
covariance matrix of the state signal is studied above, while for the
piecewise-constant randomization signal $v(t)$ in (14), we have

$\int\limits_{0}^{t}{v(s)v(s)^{\top}}\boldsymbol{d}{s}=\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}\int\limits_{n\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}^{(n+1)\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}}{w_{n}w_{n}^{\top}}\boldsymbol{d}{s}=\boldsymbol{\tau}\boldsymbol{\kappa_{\mathrm{}}}^{-1}\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}w_{n}w_{n}^{\top}.$

Thus, according to Theorem 3, similar to (19) we have

${\left|\kern-1.29167pt\left|\sum\limits_{n=0}^{\boldsymbol{\kappa_{\mathrm{}}}-1}w_{n}w_{n}^{\top}-\boldsymbol{\kappa_{\mathrm{}}}\sigma_{w}^{2}I_{q}\right|\kern-1.29167pt\right|}\lesssim\boldsymbol{\kappa_{\mathrm{}}}^{1/2}\sigma_{w}^{2}q\log^{3/2}\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta},$

with probability at least $1-\delta$, which for

$H_{\boldsymbol{\tau}}=\int\limits_{0}^{\boldsymbol{\tau}}{\begin{bmatrix}0_{p}\\\
v(s)\end{bmatrix}\begin{bmatrix}0_{p}\\\
v(s)\end{bmatrix}^{\top}}\boldsymbol{d}{s}-\boldsymbol{\tau}\sigma_{w}^{2}\begin{bmatrix}0_{p\times
p}&0_{p\times q}\\\ 0_{q\times p}&I_{q}\end{bmatrix},$

leads to

${\left|\kern-1.29167pt\left|H_{\boldsymbol{\tau}}\right|\kern-1.29167pt\right|}\lesssim\sigma_{w}^{2}q\log^{3/2}\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta},$
(32)

because $\boldsymbol{\kappa_{\mathrm{}}}\gtrsim\boldsymbol{\tau}^{2}$.

Next, using
$\boldsymbol{z}_{s}=\left[\boldsymbol{x}_{s}^{\top},\boldsymbol{x}_{s}^{\top}K^{\top}+v(s)^{\top}\right]^{\top}$,
the matrix $\widehat{\Sigma}_{\boldsymbol{\tau}}$ can be written as

$\widehat{\Sigma}_{\boldsymbol{\tau}}=\begin{bmatrix}I_{p}\\\
K\end{bmatrix}\Phi_{\boldsymbol{\tau}}\begin{bmatrix}I_{p}\\\
K\end{bmatrix}^{\top}+\boldsymbol{\tau}\sigma_{w}^{2}\begin{bmatrix}0_{p\times
p}&0_{p\times q}\\\ 0_{q\times
p}&I_{q}\end{bmatrix}+F_{\boldsymbol{\tau}}+H_{\boldsymbol{\tau}},$ (33)

where

$F_{\boldsymbol{\tau}}=\int\limits_{0}^{\boldsymbol{\tau}}{\left(\begin{bmatrix}I_{p}\\\
K\end{bmatrix}\boldsymbol{x}_{s}\begin{bmatrix}0_{p}\\\
v(s)\end{bmatrix}^{\top}+\begin{bmatrix}0_{p}\\\
v(s)\end{bmatrix}\boldsymbol{x}_{s}^{\top}\begin{bmatrix}I_{p}\\\
K\end{bmatrix}^{\top}\right)}\boldsymbol{d}{s}.$

However, Lemma 1 and
$\boldsymbol{\kappa_{\mathrm{}}}\gtrsim\boldsymbol{\tau}^{2}$ give a high
probability upper-bound for the above matrix:

${\left|\kern-1.29167pt\left|F_{\boldsymbol{\tau}}\right|\kern-1.29167pt\right|}\lesssim\left(1\vee{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right)\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}\right)\left(pq^{1/2}\boldsymbol{\tau}^{1/2}\log\frac{pq}{\delta}+q\log^{3/2}\frac{\boldsymbol{\kappa_{\mathrm{}}}q}{\delta}\right).$
(34)

In the sequel, we show that with probability at least $1-\delta$, it holds
that

$\underline{\lambda}\left(\begin{bmatrix}I_{p}\\\
K\end{bmatrix}\Phi_{\boldsymbol{\tau}}\begin{bmatrix}I_{p}\\\
K\end{bmatrix}^{\top}+\boldsymbol{\tau}\sigma_{w}^{2}\begin{bmatrix}0_{p\times
p}&0_{p\times q}\\\ 0_{q\times
p}&I_{q}\end{bmatrix}\right)\gtrsim\boldsymbol{\tau}\left(\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}\right)\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}^{2}\right)^{-1},$

which, according to (32), (33), and (34), implies the desired result. To show
the above least eigenvalue inequality, we use
$\underline{\lambda}\left(\Phi_{\boldsymbol{\tau}}\right)\gtrsim\boldsymbol{\tau}\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)$
to obtain

$\underline{\lambda}\left(\begin{bmatrix}I_{p}\\\
K\end{bmatrix}\Phi_{\boldsymbol{\tau}}\begin{bmatrix}I_{p}\\\
K\end{bmatrix}^{\top}+\boldsymbol{\tau}\sigma_{w}^{2}\begin{bmatrix}0_{p\times
p}&0_{p\times q}\\\ 0_{q\times
p}&I_{q}\end{bmatrix}\right)\gtrsim\boldsymbol{\tau}\left(\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}\right)\underline{\lambda}\left(\begin{bmatrix}I_{p}&K^{\top}\\\
K&KK^{\top}+I_{q}\end{bmatrix}\right).$

However, block matrix inversion gives

$\underline{\lambda}\left(\begin{bmatrix}I_{p}&K^{\top}\\\
K&KK^{\top}+I_{q}\end{bmatrix}\right)=\overline{\lambda}\left(\begin{bmatrix}I_{p}&K^{\top}\\\
K&KK^{\top}+I_{q}\end{bmatrix}^{-1}\right)^{-1}=\overline{\lambda}\left(\begin{bmatrix}K^{\top}K+I_{p}&-K^{\top}\\\
-K&I_{q}\end{bmatrix}\right)^{-1},$

that is clearly at least
$\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}^{2}\right)^{-1}$,
apart from a constant factor. Therefore, we get the desired result.
$\blacksquare$

###  A.4 Approximation of true drift parameter by Algorithm 1

###### Lemma 4

Suppose that $\boldsymbol{\widehat{\theta}}$ is given by Algorithm 1. Then,
with probability at least $1-\delta$, we have

${\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}\lesssim\frac{p\left(p+q\right)^{1/2}}{\boldsymbol{\tau}^{1/2}}\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\vee\sigma_{w}^{2}}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\wedge\sigma_{w}^{2}}\left(1+{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right)^{3}\log\left(\frac{pq\boldsymbol{\kappa_{\mathrm{}}}}{\delta}\right).$
(35)

Proof. First, consider the mean matrix of the Gaussian posterior distribution.
Using the data generation mechanism
$\boldsymbol{d}\boldsymbol{x}_{t}=\boldsymbol{\theta_{0}}^{\top}\boldsymbol{z}_{t}\boldsymbol{d}t+\boldsymbol{d}\mathbb{W}_{t}$,
we have

$\widehat{\mu}_{\boldsymbol{\tau}}=\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}}\boldsymbol{d}{\boldsymbol{x}_{s}^{\top}}=\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\left(\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}\boldsymbol{z}_{s}^{\top}}\boldsymbol{d}{s}\boldsymbol{\theta_{0}}+\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}\right)=\boldsymbol{\theta_{0}}-\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\left(\boldsymbol{\theta_{0}}-\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}\right),$

where we used the definition of $\widehat{\Sigma}_{\tau}$ in (8). Now, the
sample $\boldsymbol{\widehat{\theta}}$ from $\mathcal{D}_{\boldsymbol{\tau}}$
can be written as
$\boldsymbol{\widehat{\theta}_{\boldsymbol{\tau}}}=\widehat{\mu}_{\boldsymbol{\tau}}+\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1/2}\Phi$,
where $\Phi\sim\boldsymbol{N}\left(0_{(p+q)\times p},I_{p+q}\right)$ is a
standard normal random matrix, as defined in the notation. So, for the error
matrix, it holds that

${\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}\leq{\left|\kern-1.29167pt\left|\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\right|\kern-1.29167pt\right|}\left({\left|\kern-1.29167pt\left|\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}+{\left|\kern-1.29167pt\left|\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}\right|\kern-1.29167pt\right|}\right)+{\left|\kern-1.29167pt\left|\widehat{\Sigma}_{\boldsymbol{\tau}}^{-1}\right|\kern-1.29167pt\right|}^{1/2}{\left|\kern-1.29167pt\left|\Phi\right|\kern-1.29167pt\right|}.$

To proceed towards bounding the above error matrix, use

$\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}=\int\limits_{0}^{\boldsymbol{\tau}}{\left(\begin{bmatrix}I_{p}\\\
K\end{bmatrix}\boldsymbol{x}_{s}+\begin{bmatrix}0\\\
v(s)\end{bmatrix}\right)}\boldsymbol{d}{\mathbb{W}_{s}^{\top}},$

to obtain

${\left|\kern-1.29167pt\left|\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{z}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}\right|\kern-1.29167pt\right|}\leq\left(1\vee{\left|\kern-1.29167pt\left|K\right|\kern-1.29167pt\right|}\right){\left|\kern-1.29167pt\left|\int\limits_{0}^{\boldsymbol{\tau}}{\boldsymbol{x}_{s}}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}\right|\kern-1.29167pt\right|}+{\left|\kern-1.29167pt\left|\int\limits_{0}^{\boldsymbol{\tau}}{v(s)}\boldsymbol{d}{\mathbb{W}_{s}^{\top}}\right|\kern-1.29167pt\right|},$

To proceed, note that with probability at least $1-\delta$, we have

${\left|\kern-1.29167pt\left|\Phi\right|\kern-1.29167pt\right|}^{2}\lesssim
p(p+q)\log\frac{p(p+q)}{\delta}.$

Now, by putting this together with the results of Lemma 2, Lemma 3, and (25),
we get the desired result.

$\blacksquare$

###  A.5 Eigenvalue ratio bound for sum of two matrices

###### Lemma 5

Suppose that $M,{E}$ are $p\times p$ matrices, and let
$M=\Gamma^{-1}\Lambda\Gamma$ be the Jordan diagonalization of $M$. So, for
some positive integer $k$, we have $\Lambda\in\mathbb{C}^{p\times
p}=\mathrm{diag}\left(\Lambda_{1},\cdots,\Lambda_{k}\right)$, where the blocks
${\Lambda_{1},\cdots,\Lambda_{k}}$ are Jordan matrices of the form

$\displaystyle\Lambda_{i}=\begin{bmatrix}\lambda_{i}&1&0&\cdots&0&0\\\
0&\lambda_{i}&1&0&\cdots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&\cdots&0&\lambda_{i}&1\\\
0&0&0&\cdots&0&\lambda_{i}\end{bmatrix}\in\mathbb{C}^{\boldsymbol{r}_{i}\times\boldsymbol{r}_{i}}.$

Then, let $\boldsymbol{r}=\max\limits_{1\leq i\leq k}\boldsymbol{r}_{i}\leq
p$, and define $\Delta_{M}({E})$ as the difference between the largest real-
part of the eigenvalues of $M+{E}$ and that of $M$. Then, it holds that

$\Delta_{M}({E})\leq\left(1\vee\boldsymbol{r}^{1/2}{\left|\kern-1.29167pt\left|{E}\right|\kern-1.29167pt\right|}\mathrm{cond}\left(\Gamma\right)\right)^{1/\boldsymbol{r}},$
(36)

where $\mathrm{cond}\left(\Gamma\right)$ is the condition number of $\Gamma$:
$\mathrm{cond}\left(\Gamma\right)=\overline{\lambda}\left(\Gamma^{\top}\Gamma\right)^{1/2}\underline{\lambda}\left(\Gamma^{\top}\Gamma\right)^{-1/2}$.

Proof. Since the expression on the right-hand-side of (36) is positive, it is
enough to consider an eigenvalue $\lambda$ of $M+{E}$ which is not an
eigenvalue of $M$, and show that
$\Re\left(\lambda\right)-\log\overline{\lambda}\left(\exp\left(M\right)\right)$
is less than the expression on the RHS of $\eqref{EigPerturb}$. So, for such
$\lambda$, the matrix $M-\lambda I_{p}$ is non-singular, while $M+{E}-\lambda
I_{p}$ is singular. Let the vector $v\neq 0$ be such that $\left(M+{E}-\lambda
I_{p}\right)v=0$, which by Jordan diagonalization above implies that

$v=-\Gamma^{-1}\left(\Lambda-\lambda I\right)^{-1}\Gamma{E}v.$ (37)

Then, $\Lambda=\mathrm{diag}\left(\Lambda_{1},\cdots,\Lambda_{k}\right)$
indicates that $\Lambda-\lambda I$ and $\left(\Lambda-\lambda I\right)$ are
block diagonal, the latter consisting of the blocks
$\mathrm{diag}\left(\left(\Lambda_{1}-\lambda
I_{\boldsymbol{r}_{1}}\right)^{-1},\cdots,\left(\Lambda_{k}-\lambda
I_{\boldsymbol{r}_{k}}\right)^{-1}\right)$.

Now, multiplications show that

$\left(\Lambda_{i}-\lambda
I_{\boldsymbol{r}_{i}}\right)^{-1}=-\begin{bmatrix}\left(\lambda-\lambda_{i}\right)^{-1}&\left(\lambda-\lambda_{i}\right)^{-2}&\cdots&\left(\lambda-\lambda_{i}\right)^{-\boldsymbol{r}_{i}}\\\
0&\left(\lambda-\lambda_{i}\right)^{-1}&\cdots&\left(\lambda-\lambda_{i}\right)^{-\boldsymbol{r}_{i}+1}\\\
\vdots&\vdots&\vdots&\vdots\\\
0&\cdots&0&\left(\lambda-\lambda_{i}\right)^{-1}\end{bmatrix}.$

Therefore, according to the definition of matrix operator norms in Section 1,
we obtain

${\left|\kern-1.29167pt\left|\left(\Lambda_{i}-\lambda
I_{\boldsymbol{r}_{i}}\right)^{-1}\right|\kern-1.29167pt\right|}^{2}\leq\boldsymbol{r}\left(1\vee\left|\lambda-\lambda_{i}\right|^{-\boldsymbol{r}}\right)^{2}.$

Putting these bounds for the blocks of $\left(\Lambda-\lambda I\right)^{-1}$
together, (37) leads to

$\displaystyle 1$ $\displaystyle\leq$
$\displaystyle{\left|\kern-1.29167pt\left|\left(\Lambda-\lambda
I\right)^{-1}\right|\kern-1.29167pt\right|}{\left|\kern-1.29167pt\left|\Gamma\right|\kern-1.29167pt\right|}{\left|\kern-1.29167pt\left|\Gamma^{-1}\right|\kern-1.29167pt\right|}{\left|\kern-1.29167pt\left|{E}\right|\kern-1.29167pt\right|}$
$\displaystyle\leq$
$\displaystyle\boldsymbol{r}^{1/2}\mathrm{cond}\left(\Gamma\right){\left|\kern-1.29167pt\left|{E}\right|\kern-1.29167pt\right|}\max\limits_{1\leq
i\leq
k}\left(1\wedge\left|\lambda-\lambda_{i}\right|^{\boldsymbol{r}}\right)^{-1}$
$\displaystyle\leq$
$\displaystyle\boldsymbol{r}^{1/2}\mathrm{cond}\left(\Gamma\right){\left|\kern-1.29167pt\left|{E}\right|\kern-1.29167pt\right|}\left(1\wedge\left(\Re(\lambda)-\log\overline{\lambda}\left(\exp\left(M\right)\right)\right)^{\boldsymbol{r}}\right)^{-1}.$

To see the last inequality above, note that if
$\Re(\lambda)-\log\overline{\lambda}\left(\exp\left(M\right)\right)$ is
positive, then it is larger than all the terms
$\left|\lambda-\lambda_{i}\right|$, for $i=1,\cdots,k$. Thus, for

$\Re\left(\lambda\right)=\log\overline{\lambda}\left(\exp\left(M+{E}\right)\right),$

we obtain (36). $\blacksquare$

##  Appendix B Proof of Theorem 2

To establish the rates of exploration Algorithm 2 performs, we utilize Lemma
8, which indicates that

${\left|\kern-1.29167pt\left|\widehat{\mu}_{\boldsymbol{\tau}_{n}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}\lesssim{\left|\kern-1.29167pt\left|\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}^{-1/2}\right|\kern-1.29167pt\right|}{\left|\kern-1.29167pt\left|\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}^{-1/2}\int\limits_{0}^{\boldsymbol{\tau}_{n}}{\boldsymbol{x}_{t}}\boldsymbol{d}{\mathbb{W}_{t}^{\top}}\right|\kern-1.29167pt\right|}\lesssim\underline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}\right)^{-1/2}\left(p(p+q)\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\log\overline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}\right)\right)^{1/2}.$

Now, (51) gives
$\log\overline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}\right)\lesssim\log\boldsymbol{\tau}_{n}$,
while Lemma 9 provides
$\underline{\lambda}\left(\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}\right)\gtrsim\boldsymbol{\tau}_{n}^{1/2}\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)$.
Moreover, since
$\widehat{\Sigma}_{\boldsymbol{\tau}_{n}}^{1/2}\left(\boldsymbol{\widehat{\theta}_{n}}-\widehat{\mu}_{\boldsymbol{\tau}_{n}}\right)$
is a standard normal $(p+q)\times p$ matrix, we have

${\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}_{n}}-\widehat{\mu}_{\boldsymbol{\tau}_{n}}\right|\kern-1.29167pt\right|}\lesssim\boldsymbol{\tau}_{n}^{-1/4}\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)^{-1/2}\left(p(p+q)\log(pq)\right)^{1/2}.$

Thus, we obtain the desired result for the estimation error.

To proceed toward establishing the regret bound, Lemma 7 shows that we need to
integrate
${\left|\kern-1.29167pt\left|\boldsymbol{u}_{t}+Q_{u}^{-1}\widehat{B}_{n}^{\top}{P}\left(\boldsymbol{\widehat{\theta}_{n}}\right)\boldsymbol{x}_{t}\right|\kern-1.29167pt\right|}^{2}$
over the stabilized period of Algorithm 2: $\boldsymbol{\tau}_{0}\leq t\leq
T$:

$\mathrm{Reg}\left(T\right)\lesssim\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}\right)\boldsymbol{\tau}_{0}+\int\limits_{\boldsymbol{\tau}_{0}}^{T}{{\left|\kern-1.50696pt\left|\boldsymbol{u}_{t}+Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\boldsymbol{x}_{t}\right|\kern-1.50696pt\right|}^{2}}\boldsymbol{d}{t}.$

Further, according to (51), for
$\boldsymbol{\tau}_{n-1}<T\leq\boldsymbol{\tau}_{n}$, we have

$\int\limits_{\boldsymbol{\tau}_{0}}^{T}{{\left|\kern-1.50696pt\left|\boldsymbol{u}_{t}+Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\boldsymbol{x}_{t}\right|\kern-1.50696pt\right|}^{2}}\boldsymbol{d}{t}\lesssim\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)\sum\limits_{i=0}^{n-1}\left(\boldsymbol{\tau}_{i+1}-\boldsymbol{\tau}_{i}\right){\left|\kern-1.29167pt\left|{K}\left(\boldsymbol{\widehat{\theta}_{i}}\right)-{K}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}^{2}.$

On the other hand, Lemma 12 implies that

${\left|\kern-1.29167pt\left|{K}\left(\boldsymbol{\widehat{\theta}_{i}}\right)-{K}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}\lesssim\frac{{\left|\kern-1.29167pt\left|{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}^{3}}{\underline{\lambda}\left(Q_{x}\right)\underline{\lambda}\left(Q_{u}\right)^{2}}{\left|\kern-1.29167pt\left|\boldsymbol{\widehat{\theta}_{i}}-\boldsymbol{\theta_{0}}\right|\kern-1.29167pt\right|}.$

Thus, we have

$\mathrm{Reg}\left(T\right)\lesssim\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}\right)\boldsymbol{\tau}_{0}+\frac{\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)^{2}}{\underline{\lambda}\left(\Sigma_{\mathbb{W}}\right)}\frac{{\left|\kern-1.29167pt\left|{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}^{6}}{\underline{\lambda}\left(Q_{x}\right)^{2}\underline{\lambda}\left(Q_{u}\right)^{4}}p(p+q)\sum\limits_{i=0}^{n-1}\left(\boldsymbol{\tau}_{i+1}-\boldsymbol{\tau}_{i}\right)\frac{\log\boldsymbol{\tau}_{i}}{\boldsymbol{\tau}_{i}^{1/2}}.$

Thus, according to (12), we obtain the desired regret bound result in Theorem
2.

###  B.1 Geometry of drift parameters and optimal policies

###### Lemma 6

For the drift parameter $\boldsymbol{\theta_{1}}$, and for
$X\in\mathbb{R}^{p\times p},Y\in\mathbb{R}^{p\times q}$, define

$\displaystyle\Delta_{\boldsymbol{\theta_{1}}}(X,Y)={P}\left(\boldsymbol{\theta_{1}}\right)Y+\int\limits_{0}^{\infty}{e^{\overline{A}_{1}^{\top}t}\left[M\left(X,Y\right)^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)+{P}\left(\boldsymbol{\theta_{1}}\right)M\left(X,Y\right)\right]e^{\overline{A}_{1}t}B_{1}}\boldsymbol{d}{t},$

where
$\overline{A}_{1}=A_{1}-B_{1}Q_{u}^{-1}B_{1}^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)$
and
$M\left(X,Y\right)=X-YQ_{u}^{-1}B_{1}^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)$.
Then, $\Delta_{\boldsymbol{\theta_{1}}}(X,Y)$ is the directional derivative of
$B^{\top}{P}\left(\boldsymbol{\theta}\right)$ at $\boldsymbol{\theta_{1}}$ in
the direction $\left[X,Y\right]$. Importantly, the tangent space of the
manifold of matrices $\boldsymbol{\theta}\in\mathbb{R}^{p\times(p+q)}$ that
satisfy
$B^{\top}{P}\left(\boldsymbol{\theta}\right)=B_{1}^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)$
at $\boldsymbol{\theta_{1}}$ contains all matrices $X,Y$ that
$\Delta_{\boldsymbol{\theta_{1}}}(X,Y)=0$.

Proof. First, note that according to the Lipschitz continuity of
${P}\left(\boldsymbol{\theta}\right)$ in Lemma 12, the directional derivative
exists and is well-defined, as long as
${\left|\kern-1.29167pt\left|{P}\left(\boldsymbol{\theta_{1}}\right)\right|\kern-1.29167pt\right|}<\infty$.
However, Lemma 11 provides that ${P}\left(\boldsymbol{\theta_{1}}\right)$ is
finite in a neighborhood of $\boldsymbol{\theta_{0}}$, and so the required
condition holds. Below, we start by establishing the second result to identify
the tangent space, and then prove the general result on the directional
derivative.

To proceed, let
$\boldsymbol{\theta}=\boldsymbol{\theta_{1}}+\epsilon\left[X,Y\right]^{\top}$
be such that
$B^{\top}{P}\left(\boldsymbol{\theta}\right)=B_{1}^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)$,
and denote
${K}\left(\boldsymbol{\theta_{1}}\right)=-Q_{u}^{-1}B_{1}^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)$.
So, the directional derivative of ${P}\left(\boldsymbol{\theta_{1}}\right)$
along the matrix $\left[X,Y\right]^{\top}$ can be found as follows. First,
denoting the closed-loop transition matrix by
$\overline{A}=A-BQ_{u}^{-1}B^{\top}{P}\left(\boldsymbol{\theta}\right)$, since

$\overline{A}^{\top}{P}\left(\boldsymbol{\theta}\right)+{P}\left(\boldsymbol{\theta}\right)\overline{A}+Q_{x}+{K}\left(\boldsymbol{\theta}\right)^{\top}Q_{u}{K}\left(\boldsymbol{\theta}\right)=0,$

we have

$\displaystyle\left(\overline{A}_{1}+\epsilon X+\epsilon
Y{K}\left(\boldsymbol{\theta_{1}}\right)\right)^{\top}{P}\left(\boldsymbol{\theta}\right)+{P}\left(\boldsymbol{\theta}\right)\left(\overline{A}_{1}+\epsilon
X+\epsilon Y{K}\left(\boldsymbol{\theta_{1}}\right)\right)$ $\displaystyle=$
$\displaystyle-
Q_{x}-{K}\left(\boldsymbol{\theta_{1}}\right)^{\top}Q_{u}{K}\left(\boldsymbol{\theta_{1}}\right)=\overline{A}_{1}^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)+{P}\left(\boldsymbol{\theta_{1}}\right)\overline{A}_{1}.$

For the matrix ${E}=\lim\limits_{\epsilon\to
0}\epsilon^{-1}\left({P}\left(\boldsymbol{\theta}\right)-{P}\left(\boldsymbol{\theta_{1}}\right)\right)$,
the latter result implies that

$\overline{A}_{1}^{\top}{E}+{E}\overline{A}_{1}+\left(X+Y{K}\left(\boldsymbol{\theta_{1}}\right)\right)^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)+{P}\left(\boldsymbol{\theta_{1}}\right)\left(X+Y{K}\left(\boldsymbol{\theta_{1}}\right)\right)=0.$

Then, since all eigenvalues of ${\overline{A}_{1}}$ are in the open left half-
plane, the above Lyapunov equation for ${E}$ leads to the integral form

${E}=\int\limits_{0}^{\infty}{e^{\overline{A}_{1}^{\top}t}\left(\left(X+Y{K}\left(\boldsymbol{\theta_{1}}\right)\right)^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)+{P}\left(\boldsymbol{\theta_{1}}\right)\left(X+Y{K}\left(\boldsymbol{\theta_{1}}\right)\right)\right)e^{\overline{A}_{1}t}}\boldsymbol{d}{t}.$

On the other hand,
${K}\left(\boldsymbol{\theta}\right)=-Q_{u}^{-1}B^{\top}{P}\left(\boldsymbol{\theta}\right)$
gives

$\displaystyle 0=\lim\limits_{\epsilon\to
0}\frac{1}{\epsilon}\left(B^{\top}{P}\left(\boldsymbol{\theta}\right)-B_{1}^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)\right)=\lim\limits_{\epsilon\to
0}\frac{1}{\epsilon}\left[\left(B^{\top}-B_{1}^{\top}\right){P}\left(\boldsymbol{\theta}\right)-B_{1}^{\top}\left({P}\left(\boldsymbol{\theta_{1}}\right)-{P}\left(\boldsymbol{\theta}\right)\right)\right],$

which, according to the definitions of ${E},M(X,Y)$, implies the desired
result about the tangent space of the manifold under consideration.

Next, to establish the more general result on the directional derivative, we
use the directional derivative of ${P}\left(\boldsymbol{\theta}\right)$ in
(65):

$\int\limits_{0}^{\infty}{e^{\overline{A}_{1}^{\top}t}\left({P}\left(\boldsymbol{\theta_{1}}\right)\left[X+Y{K}\left(\boldsymbol{\theta_{1}}\right)\right]+\left[X+Y{K}\left(\boldsymbol{\theta_{1}}\right)\right]^{\top}{P}\left(\boldsymbol{\theta_{1}}\right)\right)e^{\overline{A}_{1}t}}\boldsymbol{d}{t}.$

Finally, since the directional derivative for $B^{\top}$ is $Y$, for
$B^{\top}{P}\left(\boldsymbol{\theta}\right)$, by the product rule it is
$\Delta_{\boldsymbol{\theta_{1}}}(X,Y)$, which finishes the proof.
$\blacksquare$

###  B.2 Regret bounds in terms of deviations in control actions

###### Lemma 7

Let $\boldsymbol{u}_{t}$ be the action that Algorithm 2 takes at time $t$.
Then, for the regret of Algorithm 2, it holds that

$\displaystyle\mathrm{Reg}\left(T\right)$ $\displaystyle\lesssim$
$\displaystyle\left(\overline{\lambda}\left(\Sigma_{\mathbb{W}}\right)+\sigma_{w}^{2}\right)\boldsymbol{\tau}_{0}{\left|\kern-1.29167pt\left|K+Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\right|\kern-1.29167pt\right|}^{2}$
$\displaystyle+$
$\displaystyle\int\limits_{\boldsymbol{\tau}_{0}}^{T}{{\left|\kern-1.50696pt\left|\boldsymbol{u}_{t}+Q_{u}^{-1}B_{0}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)\boldsymbol{x}_{t}\right|\kern-1.50696pt\right|}^{2}}\boldsymbol{d}{t}+{x^{*}_{T}}^{\top}{P}\left(\boldsymbol{\theta_{0}}\right)x^{*}_{T},$

where $x^{*}_{T}$ is the terminal state under the optimal trajectory
$\boldsymbol{\pi_{\mathrm{opt}}}$ in (6).

Proof. First, denote the optimal linear feedback of
$\boldsymbol{\pi_{\mathrm{opt}}}$ in (6) by
$\boldsymbol{u}_{t}={K}\left(\boldsymbol{\theta_{0}}\right)\boldsymbol{x}_{t}$,
where
${K}\left(\boldsymbol{\theta_{0}}\right)=-Q_{u}^{-1}B_{0}{P}\left(\boldsymbol{\theta_{0}}\right)$.
According to the episodic structure of Algorithm 2, for
$\boldsymbol{\tau}_{n}\leq t<\boldsymbol{\tau}_{n+1}$, denote

$K_{t}=-Q_{u}^{-1}\widehat{B}_{n}^{\top}{P}\left(\boldsymbol{\widehat{\theta}_{n}}\right).$

We first consider the regret of Algorithm 2 after finishing stabilization by
running Algorithm 1; i.e., for $\boldsymbol{\tau}_{0}\leq t\leq T$. Fix some
small $\epsilon>0$, that we will let decay later. We proceed by finding an
approximation of the regret through sampling at times
$\boldsymbol{\tau}_{0}+k\epsilon$, for non-negative integers $k$. To do that,
denote $N=\lceil(T-\boldsymbol{\tau}_{0})/\epsilon\rceil$, and define the
sequence of policies
$\left\\{\boldsymbol{\widehat{\pi}}_{i}\right\\}_{i=0}^{N}$ according to

$\boldsymbol{\widehat{\pi}}_{i}=\begin{cases}\boldsymbol{u}_{t}=K_{t}\boldsymbol{x}_{t}&t<\boldsymbol{\tau}_{0}+i\epsilon\\\
\boldsymbol{u}_{t}={K}\left(\boldsymbol{\theta_{0}}\right)\boldsymbol{x}_{t}&t\geq\boldsymbol{\tau}_{0}+i\epsilon\end{cases}.$

That is, the policy $\boldsymbol{\widehat{\pi}}_{i}$ switches to the optimal
feedback at time $\boldsymbol{\tau}_{0}+i\epsilon$. So, the zeroth policy
$\boldsymbol{\widehat{\pi}}_{0}$ corresponds to applying the optimal policy
$\boldsymbol{\pi_{\mathrm{opt}}}$ after stabilization at time
$\boldsymbol{\tau}_{0}$, while the last one $\boldsymbol{\widehat{\pi}}_{N}$
is nothing but the one in Algorithm 2, that we denote by
$\boldsymbol{\widehat{\pi}}$, for the sake of brevity. As such, we have
$\mathrm{Reg}_{\boldsymbol{\widehat{\pi}}_{0}}\left(T\right)=0$, and the
telescopic summation below holds true:

$\mathrm{Reg}_{\boldsymbol{\widehat{\pi}}}\left(T\right)=\sum\limits_{i=0}^{N-1}\left(\mathrm{Reg}_{\boldsymbol{\widehat{\pi}}_{i+1}}\left(T\right)-\mathrm{Reg}_{\boldsymbol{\widehat{\pi}}_{i}}\left(T\right)\right).$
(38)

Now, to consider the difference
$\mathrm{Reg}_{\boldsymbol{\widehat{\pi}}_{i+1}}\left(T\right)-\mathrm{Reg}_{\boldsymbol{\widehat{\pi}}_{i}}\left(T\right)$,
for a fixed $i$ in the range $0\leq i<N$, denote
$t_{1}=\boldsymbol{\tau}_{0}+i\epsilon$ and let
$\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}},\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i+1}}$
be the state trajectories under
$\boldsymbol{\widehat{\pi}}_{i},\boldsymbol{\widehat{\pi}}_{i+1}$,
respectively. By definition, we have
$\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}=\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i+1}}$,
for all $t\leq t_{1}$. So, we drop the policy superscript and use
$\boldsymbol{x}_{t_{1}}$ to refer to the states of both of them at time
$t_{1}$. Therefore, as long as $t_{1}\leq t<t_{1}+\epsilon$, similar to (15),
the solutions of the stochastic differential equation are

$\displaystyle\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}$
$\displaystyle=$ $\displaystyle
e^{\overline{A}_{0}(t-t_{1})}\boldsymbol{x}_{t_{1}}+\int\limits_{t_{1}}^{t}{e^{\overline{A}_{0}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}},$
$\displaystyle\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i+1}}$
$\displaystyle=$ $\displaystyle
e^{\overline{A}(t-t_{1})}\boldsymbol{x}_{t_{1}}+\int\limits_{t_{1}}^{t}{e^{\overline{A}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}},$

where $\overline{A}_{0}=A_{0}+B_{0}{K}\left(\boldsymbol{\theta_{0}}\right)$
and $\overline{A}=A_{0}+B_{0}K_{t_{1}}$ are the closed-loop transition
matrices under $\boldsymbol{\widehat{\pi}}_{i}$ and
$\boldsymbol{\widehat{\pi}}_{i+1}$, respectively. To work with the above two
state trajectories, we define some notations for convenience:

$\displaystyle M_{0}$ $\displaystyle=$ $\displaystyle
Q_{x}+{K}\left(\boldsymbol{\theta_{0}}\right)^{\top}Q_{u}{K}\left(\boldsymbol{\theta_{0}}\right),$
$\displaystyle M_{1}$ $\displaystyle=$ $\displaystyle
Q_{x}+K_{t_{1}}Q_{u}K_{t_{1}},$ $\displaystyle y_{t}$ $\displaystyle=$
$\displaystyle\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i+1}}-\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}},$
$\displaystyle{E}_{t}$ $\displaystyle=$ $\displaystyle
e^{\overline{A}(t-t_{1})}-e^{\overline{A}_{0}(t-t_{1})}.$

Thus, letting

$Z_{t}=\int\limits_{t_{1}}^{t}{\left[e^{\overline{A}(t-s)}-e^{\overline{A}_{0}(t-s)}\right]}\boldsymbol{d}{\mathbb{W}_{s}},$

it holds that
$y_{t}={E}_{t}\boldsymbol{x}_{t_{1}}+Z_{t}+{O}\left(\epsilon^{2}\right)$.
Further, for the observation signal $\boldsymbol{z}_{t}$ and the cost matrix
$Q$ defined in Section 2, we have

$\displaystyle\int\limits_{t_{1}}^{t_{1}+\epsilon}{\left({\left|\kern-1.50696pt\left|Q^{1/2}\boldsymbol{z}_{t}\left(\boldsymbol{\widehat{\pi}}_{i+1}\right)\right|\kern-1.50696pt\right|}^{2}-{\left|\kern-1.50696pt\left|Q^{1/2}\boldsymbol{z}_{t}\left(\boldsymbol{\widehat{\pi}}_{i}\right)\right|\kern-1.50696pt\right|}^{2}\right)}\boldsymbol{d}{t}$
(39) $\displaystyle=$
$\displaystyle\int\limits_{t_{1}}^{t_{1}+\epsilon}{\left[\left(\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}+y_{t}\right)^{\top}M_{1}\left(\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}+y_{t}\right)-{\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}}^{\top}M_{0}\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}\right]}\boldsymbol{d}{t}$
$\displaystyle=$
$\displaystyle\int\limits_{t_{1}}^{t_{1}+\epsilon}{\left[{\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}}^{\top}S\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}+2y_{t}^{\top}M_{1}\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}+y_{t}^{\top}M_{1}y_{t}\right]}\boldsymbol{d}{t},$

where
$S=M_{1}-M_{0}=K_{t_{1}}^{\top}Q_{u}K_{t_{1}}-{K}\left(\boldsymbol{\theta_{0}}\right)^{\top}Q_{u}{K}\left(\boldsymbol{\theta_{0}}\right)$.

On the other hand, for $t\geq t_{1}+\epsilon$, the evolutions of the state
vectors are the same for the two policies and we have

$\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}=e^{\overline{A}_{0}(t-t_{1}-\epsilon)}\boldsymbol{x}_{t_{1}+\epsilon}^{\boldsymbol{\widehat{\pi}}_{i}}+\int\limits_{t_{1}+\epsilon}^{t}{e^{\overline{A}_{0}(t-s)}}\boldsymbol{d}{\mathbb{W}_{s}}.$

Therefore, the difference signal becomes

$y_{t}=e^{\overline{A}_{0}(t-t_{1}+\epsilon)}\left[\boldsymbol{x}_{t_{1}+\epsilon}^{\boldsymbol{\widehat{\pi}}_{i+1}}-\boldsymbol{x}_{t_{1}+\epsilon}^{\boldsymbol{\widehat{\pi}}_{i}}\right]=e^{\overline{A}_{0}(t-t_{1}+\epsilon)}y_{t_{1}+\epsilon}=e^{\overline{A}_{0}(t-t_{1}+\epsilon)}\left[{E}_{t_{1}+\epsilon}\boldsymbol{x}_{t_{1}}+Z_{t_{1}+\epsilon}\right],$

and we obtain

$\displaystyle\int\limits_{t_{1}+\epsilon}^{T}{\left({\left|\kern-1.50696pt\left|Q^{1/2}\boldsymbol{z}_{t}\left(\boldsymbol{\widehat{\pi}}_{i+1}\right)\right|\kern-1.50696pt\right|}^{2}-{\left|\kern-1.50696pt\left|Q^{1/2}\boldsymbol{z}_{t}\left(\boldsymbol{\widehat{\pi}}_{i}\right)\right|\kern-1.50696pt\right|}^{2}\right)}\boldsymbol{d}{t}$
(40) $\displaystyle=$
$\displaystyle\int\limits_{t_{1}+\epsilon}^{T}{\left[\left(\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}+y_{t}\right)^{\top}M_{0}\left(\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}+y_{t}\right)-{\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}}^{\top}M_{0}\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}\right]}\boldsymbol{d}{t}$
$\displaystyle=$
$\displaystyle\int\limits_{t_{1}+\epsilon}^{T}{\left[2y_{t}^{\top}M_{0}\boldsymbol{x}_{t}^{\boldsymbol{\widehat{\pi}}_{i}}+y_{t}^{\top}M_{0}y_{t}\right]}\boldsymbol{d}{t}.$

Now, after doing some algebra, the expressions in (39) and (40) lead ro the
following for small $\epsilon$: