On Applications of Bootstrap in Continuous Space Reinforcement Learning

TL;DR

This paper explores the use of bootstrap methods in continuous space reinforcement learning, demonstrating their effectiveness in achieving low regret and accurately learning system dynamics.

Contribution

It introduces bootstrap-based policies for continuous RL, showing they attain optimal regret scaling and improve dynamic model learning accuracy.

Findings

Bootstrap policies achieve square root regret scaling.

Bootstrap methods improve learning of system dynamics.

Numerical analysis supports theoretical results.

Abstract

In decision making problems for continuous state and action spaces, linear dynamical models are widely employed. Specifically, policies for stochastic linear systems subject to quadratic cost functions capture a large number of applications in reinforcement learning. Selected randomized policies have been studied in the literature recently that address the trade-off between identification and control. However, little is known about policies based on bootstrapping observed states and actions. In this work, we show that bootstrap-based policies achieve a square root scaling of regret with respect to time. We also obtain results on the accuracy of learning the model's dynamics. Corroborative numerical analysis that illustrates the technical results is also provided.

Figures

Tables

	$A_0$	$=$			(13)
	$Q_x$	$=$			(14)
	$P\left({𝜽}_0\right)$	$=$			(15)

Equations

$$x\left(t+1\right)$$

$$c_t\left(\mathbf{\pi }\right)$$

$${\mathbb{H}}_t=\left(x\left(0\right),\cdots ,x\left(t\right),u\left(0\right),\cdots ,u\left(t-1\right)\right).$$

$$u\left(t\right)=\mathbf{\pi }\left({\mathbb{H}}_t,Q_x,Q_u\right),$$

$$\underset{n\to \infty }{\mathrm{limsup}}\frac{1}{n}\sum _{t=0}^{n-1}{c}_t\left(\mathbf{\pi }\right).$$

$${\Phi }_{\mathbf{\theta }}\left(P\right)=Q_x+A^{\prime }PA-A^{\prime }PB{\left(B^{\prime }PB+Q_u\right)}^{-1}{B}^{\prime }PA.$$

$$G\left(\mathbf{\theta }\right)=-{\left(B^{\prime }P\left(\mathbf{\theta }\right)B+Q_u\right)}^{-1}{B}^{\prime }P\left(\mathbf{\theta }\right)A.$$

$${\mathbf{\pi }}^{\star }:u\left(t\right)=G\left({\mathbf{\theta }}_0\right)x\left(t\right),t=0,1,2,\cdots .$$

$${\mathcal{R}}_n\left(\mathbf{\pi }\right)=\sum _{t=0}^{n-1}\left[c_t\left(\mathbf{\pi }\right)-c_t\left({\mathbf{\pi }}^{\star }\right)\right].$$

$$\left|{\mathbf{\lambda }}_{\max }\left(A_0+B_{0}G\left({\mathbf{\theta }}_0\right)\right)\right|<1.$$

$$F_t={\left[I_p,G_t^{\prime }\right]}^{\prime }\in {\mathbb{R}}^{q\times p},$$

$$x\left(t+1\right)={\mathbf{\theta }}_0{F}_{t}x\left(t\right)+\xi \left(t+1\right),$$

$${\tilde{\mathbf{\theta }}}_n=\arg \underset{\mathbf{\theta }\in {\mathbb{R}}^{p\times q}}{\min }\sum _{t=0}^{n-1}{\left|\left|x\left(t+1\right)-\mathbf{\theta }{F}_{t}x\left(t\right)\right|\right|}^2.$$

$$\zeta \left(t+1\right)=x\left(t+1\right)-{\tilde{\mathbf{\theta }}}_n{F}_{t}x\left(t\right),$$

$${\hat{\mathbb{P}}}_n=\frac{1}{n}\sum _{t=1}^n\mathbf{\delta }\left[\zeta \left(t\right)-{\overline{\zeta }}_n\right],$$

$${\overline{\zeta }}_n=\frac{1}{n}\sum _{t=1}^n\zeta \left(t\right),$$

$$\hat{x}\left(t+1\right)={\tilde{\mathbf{\theta }}}_n{F}_t\hat{x}\left(t\right)+\hat{\xi }\left(t+1\right),$$

$${\hat{\mathbf{\theta }}}_n=\arg \underset{\mathbf{\theta }\in {\mathbb{R}}^{p\times q}}{\min }\sum _{t=0}^{n-1}{\left|\left|\hat{x}\left(t+1\right)-\mathbf{\theta }{F}_t\hat{x}\left(t\right)\right|\right|}^2.$$

$$\underset{n\to \infty }{\mathrm{limsup}}\left(n^{-1/2}{\log }^{-2}n\right){\mathcal{R}}_n\left(\mathbf{\pi }\right)$$

$$\underset{n\to \infty }{\mathrm{limsup}}\left(n^{1/2}{\log }^{-2}n\right){\left|\left|\left|{\tilde{\mathbf{\theta }}}_n-{\mathbf{\theta }}_0\right|\right|\right|}^2$$

$$\underset{n\to \infty }{\mathrm{liminf}}{\sigma }_n^2{\mathcal{R}}_n\left(\mathbf{\pi }\right)>0,$$

$$A_0=\left[\begin{array}{ccc}1.07& 0& -0.37\\ 0.48& -0.89& 0.85\\ 0.44& 0.04& 0\end{array}\right],B_0=\left[\begin{array}{ccc}-0.48& 0.44& -0.30\\ -0.52& 0.59& 0.26\\ 0.30& -0.44& 0\end{array}\right].$$

$$A_0=\left[\begin{array}{ccc}1.07& 0& -1.04\\ 0.48& -0.89& 0.85\\ 0.44& 0.81& 0\end{array}\right],B_0=\left[\begin{array}{ccc}-0.48& 0.44& -0.30\\ -0.52& 0.59& -0.26\\ 0.30& -0.30& 0\end{array}\right].$$

$$\mathbf{\pi }:u\left(t\right)=G_{t}x\left(t\right),t=0,1,2,\cdots ,$$

$$\underset{n\to \infty }{\mathrm{limsup}}\frac{{\mathcal{R}}_n\left(\mathbf{\pi }\right)}{\sum _{t=0}^{n-1}{\left|\left|\left(G\left({\mathbf{\theta }}_0\right)-G_t\right)x\left(t\right)\right|\right|}^2+{\left|\left|\sum _{t=1}^n{\left(A_0+B_{0}G\left({\mathbf{\theta }}_0\right)\right)}^{n-t}\xi \left(t\right)\right|\right|}^2}<\infty .$$

$$\mathcal{S}\left({\mathbf{\theta }}_1\right)=\left\{\mathbf{\theta }\in {\mathbb{R}}^{p\times q}:G\left(\mathbf{\theta }\right)=G\left({\mathbf{\theta }}_1\right)\right\}.$$

$$N^{\prime }P\left({\mathbf{\theta }}_1\right)D_1+B_1^{\prime }Z+B_1^{\prime }\sum _{k=0}^{\infty }{D_1^{\prime }}^k\left(D_1^{\prime }Z+Z^{\prime }{D}_1\right){D_1}^{k+1}=0_{r\times p},$$

$$U_n$$

$$W_n$$

$$\underset{n\to \infty }{\mathrm{limsup}}\frac{\left|{\mathbf{\lambda }}_{\max }\left(W_n{W}_n^{\prime }\right)\right|}{\log \left|{\mathbf{\lambda }}_{\max }\left(U_n\right)\right|}<\infty .$$

$$\left({\tilde{\mathbf{\theta }}}_n-{\mathbf{\theta }}_0\right)U_n{\left({\tilde{\mathbf{\theta }}}_n-{\mathbf{\theta }}_0\right)}^{\prime }=W_n{W}_n^{\prime }.$$

Full text

I Introduction

In the theory of reinforcement learning, efficient algorithms with provable theoretical guarantees are established for two canonical settings. The first is finite state Markov decision processes (MDPs) with state spaces of small cardinalities [1]. The second is the continuous space setting of linear quadratic (LQ) systems [2]. In the latter one, the control action and the state both are multidimensional real vectors, and the state evolves according to stochastic linear dynamics determined by the control action. Further, the cost (or negative reward) has a quadratic form in both the state and the control input. Besides being theoretically amenable, LQ models capture a wide range of applications from air conditioning control [3] to portfolio optimization [4]. LQ models also arise when studying the behavior of nonlinear systems around the working equilibrium [5, 6]

In applications where the true system model is not known, data-driven strategies are required for decision making under uncertainty [7]. Then, the learning algorithm has to select actions amongst infinitely many options in order to steer the system toward minimizing the costs incurred. Note that, unlike the finite state MDP case, in LQ systems there is a possible danger of the state vector becoming unbounded [8, 9, 10]. Therefore, the design and analysis of reinforcement learning algorithms for LQ systems involve significantly different conceptual and technical issues to balance exploration (identification) and exploitation (control). For this purpose, one might consider to use upper-confidence bound (UCB) approaches [11, 12, 13, 14] that rely on the optimism in the face of uncertainty (OFU) principle. The UCB approach was historically first developed for finite action bandit problems [15]. While being efficient in the finite action setting, UCB-based approaches have been found to be computationally intractable in more general problems [16].

Recently, various methods for reinforcement learning have been proposed that leverage randomization strategies to guide the learning process. Randomized policy search methods have been studied both empirically, as well as theoretically (see e.g. [17, 18]). For the problem of stabilizing an unknown LQ system, an algorithm leveraging random feedback gains is proposed [19]. There is also work showing the efficiency of achieving the exploration-exploitation trade-off by randomizing the learned model through both posterior sampling [20] and additive randomization [21]. Finally, finite time analysis of Certainty Equivalent policies utilizing input perturbation has led to performance guarantees for both learning [22] and planning [16].

In this paper, we study randomized algorithms that leverage the statistical bootstrap [23] for reinforcement learning in LQ systems. Bootstrap-based exploration has been analyzed in simpler settings, such as bandit problems [24, 25]. There has been a lot of interest recently in using bootstrap-based exploration strategies especially along with deep neural networks [26, 27, 28, 29]. However, results on bootstrap-based reinforcement learning algorithms for LQ models have been limited to primarily numerical analyses for learning the model-misspecification error [30], while rigorous performance guarantees are not currently available.

Further, bootstrap methods are also of practical interest because of their robustness to misspecified models. The amount of exploration in bootstrap-based adaptive control policies is endogenously determined by the history of the system to date. Therefore, the policy adapts its decision-making strategy with possible systematic and/or latent “biases” occurring due to lack of accurate information regarding the system’s dynamics111see the discussion at the end of Section IV for more details. Examples of such “biases” include structural breaks [31], system resets [30], and misspecification of the model dimension [32, 33, 34].

The focus of this work is on the performance of reinforcement learning policies that use the residual bootstrap to balance exploration and exploitation. We show that model-based strategies that use linear regression for learning the model and bootstrapping for policy design, provide a regret that scales as the square root of the total time of interacting with the system. Further, the accuracy of learning the unknown dynamics parameter will be specified. To establish the results, we carefully examine the effect of different converging and diverging quantities involved in the problem, such as the errors in learning the model, the distributions induced by the regression residuals, the correlation within and between the observed input-state signals, and the ongoing learning-planning interactions, that are of independent interest. At the technical level, we leverage results in the literature on the bootstrap [35], martingale central limit [36, 37] and convergence theorems [38].

The remainder of the paper is organized as follows: Section II introduces the mathematical model under consideration and discusses the rigorous formulation of the problem, and also provides some necessary preliminaries. Section III describes the bootstrap procedure and the resulting reinforcement learning algorithm to design the policy. Subsequently, the main result on the performance of the proposed algorithm is presented, together with numerical work showcasing the performance of the algorithm, in Section IV.

Notation

The following notation will be employed throughout the paper. $A^{\prime}$ is the transpose of matrix or vector $A$. The largest and the smallest eigenvalues of square Hermitian matrix $A$ are denoted by $\lambda_{\max}(A)$ and $\lambda_{\min}(A)$, respectively. If $A$ is not Hermitian, the ordering of the eigenvalues is determined by their magnitudes. The norm of the $d$ dimensional vector $v$ is denoted by ${\left|\kern-1.29167pt\left|v\right|\kern-1.29167pt\right|}=\left(\sum\limits_{i=1}^{d}\left|v_{i}\right|^{2}\right)^{1/2}$, and ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\cdot\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$ is used for the operator norm of matrices: ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|A\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}=\sup\limits_{{\left|\kern-0.90417pt\left|v\right|\kern-0.90417pt\right|}=1}{\left|\kern-1.29167pt\left|Av\right|\kern-1.29167pt\right|}$. For atomic measures on Euclidean spaces we use Dirac function $\bm{\delta}{\left[v\right]}$; i.e. it denotes a unit point mass at $v\in\mathbb{R}^{d}$. Finally, the letters $\bm{\pi}$ and $\bm{\theta}$ (or $\widetilde{\bm{\theta}},\widehat{\bm{\theta}}$) are being used for generic reinforcement learning policies and model parameters, respectively, and will be rigorously defined later on.

II Setting and Problem Formulation

The model, denoted by $\mathcal{M}$, consists of multidimensional state and control vectors, parameters that specify its dynamical evolution over time and cost matrices, as defined next. The $p$ dimensional state process $\left\{x\left(t\right)\right\}_{t=0}^{\infty}$ evolves according to an unknown stochastic linear dynamic equation governed by the $r$ dimensional control action $u\left(t\right)$, and the random disturbance (or noise) process $\left\{\xi\left(t\right)\right\}_{t=1}^{\infty}$:

[TABLE]

That is, the current state $x\left(t\right)$ and the input $u\left(t\right)$ determine the next state $x\left(t+1\right)$ through the state transition matrix $A_{0}\in\mathbb{R}^{p\times p}$, and the input influence matrix $B_{0}\in\mathbb{R}^{p\times r}$, respectively.

Definition 1**.**

Henceforth, we will denote the true parameter tuple $\left[A_{0},B_{0}\right]\in\mathbb{R}^{p\times(p+r)}$ by the $p\times q$ dynamics matrix $\bm{\theta}_{0}$, with $q\equiv p+r$. Similarly, we use the parameter $\bm{\theta}\in\mathbb{R}^{p\times q}$ to denote generic dynamics matrices.

The additive noise in the stochastic dynamics (1) satisfies $\mathbb{E}\left[\xi\left(t\right)\right]=0$. For the sake of simplicity, we assume that the sequence of noise vectors are independent, and have a stationary covariance structure: $\mathbb{E}\left[\xi\left(t\right)\xi\left(t\right)^{\prime}\right]=\Sigma$. Further, $\Sigma$ is assumed to be positive definite, and $\sup\limits_{t\geq 1}\mathbb{E}\left[{\left|\kern-1.29167pt\left|\xi\left(t\right)\right|\kern-1.29167pt\right|}^{\alpha}\right]<\infty$, for some fixed $\alpha>4$. As a matter of fact, extensions to more general technical settings such as non-stationary [16] or singular covariance matrices (assuming reachability [9]), as well as conditionally independent processes [13], can be accommodated in a similar manner. Note though that the assumed noise process is not necessarily stationary in the strict sense.

We are interested in finding reinforcement learning policies to minimize the long-term average cost as formally defined next. First, suppose that $Q_{x}$ and $Q_{u}$ are the regulation weight matrices reflecting the effect of the state and the input vectors in the cost function, respectively. Specifically, letting $\bm{\pi}$ be the decision making law (policy) determining the control input $u\left(t\right)$ at every time $t$, define the quadratic instantaneous cost of $\bm{\pi}$ according to

[TABLE]

where $Q_{x}\in\mathbb{R}^{p\times p}$, and $Q_{u}\in\mathbb{R}^{r\times r}$ are symmetric positive definite matrices. Thus, (2) reflects the desire to regulate the state of the system through control actions of small magnitude.

When the dynamics follow (1), and the instantaneous cost is given by (2), we denote the model by $\mathcal{M}\left(\bm{\theta}_{0}\right)=\left(\bm{\theta}_{0},Q_{x},Q_{u}\right)$. Further, the history of the system at time $t$, denoted as $\mathbb{H}_{t}$, consists of the sequence of the control inputs applied so far, and the resulting state vectors:

[TABLE]

A reinforcement learning policy observes the history $\mathbb{H}_{t}$ at time $t$ aiming to control the cost incurred. That is, the policy $\bm{\pi}$ is a (possibly random) mapping which designs the input sequence $\left\{u\left(t\right)\right\}_{t=0}^{\infty}$ according to the history available up to that time;

[TABLE]

so that the average cost is minimized. Thus, the objective is summarized in the following regulation problem:

Problem 1**.**

Find $\bm{\pi}$ to minimize the average cost below, subject to (1), (2), and (3);

[TABLE]

Importantly, according to (3) the true dynamics parameter $\bm{\theta}_{0}$ in (1) is unknown. Therefore, the policy must also employ an exploration procedure to accurately learn the model parameters, thus addressing the following identification problem:

Problem 2**.**

Using (1) and (3), design $\bm{\pi}$ to learn $\bm{\theta}_{0}$, as accurately as possible.

Note that in the above formulation, the true system dynamics are unknown, while the cost matrices are known. It gives rise to a realistic setting, since the decision making algorithm does not know the actual evolution of the underlying system (i.e. $\bm{\theta}_{0}$), but is aware of the criteria according to which a policy to achieve the goal is being assessed (i.e. $Q_{x},Q_{u}$).

Subsequently, we define the regret of a policy, which is the amount of sub-optimality it incurs due to uncertainty about the parameters of the model (1). To do so, we need to introduce an optimal policy $\bm{\pi^{\star}}$ that minimizes the average cost, given full knowledge of the system model $\mathcal{M}\left(\bm{\theta}_{0}\right)$. Then, $\bm{\pi^{\star}}$ will be the baseline for assessing the exploitation performance of the arbitrary reinforcement learning policy $\bm{\pi}$. It is well known that in order to find $\bm{\pi^{\star}}$, an algebraic Riccati equation needs to be solved [39, 40].

To proceed, we introduce some additional notation. First, for an arbitrary $\bm{\theta}=\left[A,B\right]$ define the matrix valued mapping

[TABLE]

Both the domain and the range of $\Phi_{\bm{\theta}}\left(\cdot\right)$ are the set of $p\times p$ matrices. Next, if there is a positive semidefinite matrix $P\left(\bm{\theta}\right)$ satisfying the algebraic Riccati equation $P\left(\bm{\theta}\right)=\Phi_{\bm{\theta}}\left(P\left(\bm{\theta}\right)\right)$, let the feedback gain matrix $G\left(\bm{\theta}\right)$ be

[TABLE]

Furthermore, for $\bm{\theta}_{0}=[A_{0},B_{0}]$ in (1), define the linear time-invariant (LTI) policy

[TABLE]

Finally, using $\bm{\pi^{\star}}$, the regret of $\bm{\pi}$ is naturally defined by:

[TABLE]

It remains to specify settings for which $\bm{\pi^{\star}}$ is well-defined. To that end, the following closed-loop stabilizability condition for the model (1) is necessary and sufficient [39, 40].

Assumption 1**.**

There is a LTI feedback gain $u\left(t\right)=G_{s}x\left(t\right)$, such that $G_{s}\in\mathbb{R}^{r\times p}$ satisfies $\left|\bm{\lambda}_{\max}\left(A_{0}+B_{0}G_{s}\right)\right|<1$.

Note that in general, the stabilizing gain $G_{s}$ mentioned above is only required to exist, and does not need to be known to the decision maker. In other words, to verify that the stabilizability Assumption 1 holds, it suffices to show that a hypothetical omniscient decision maker (who knows the true model $\mathcal{M}\left(\bm{\theta}_{0}\right)$) possessing an omnipotent computational power is able to stabilize the system. However, we briefly outline the available constructive methods to compute $G_{s}$ (as well as $\bm{\pi^{\star}}$). It is shown (see for example [39, 40, 19]) that under Assumption 1 the following statements hold;

The positive definite matrix $P\left(\bm{\theta}_{0}\right)$ uniquely exists. So, both the feedback $G\left(\bm{\theta}_{0}\right)$ and the optimal policy $\bm{\pi^{\star}}$ are well defined. 2. 2.

Letting $P_{0}$ be an arbitrary positive semidefinite $p\times p$ matrix, the recursive formula $P_{k+1}=\Phi_{\bm{\theta}_{0}}\left(P_{k}\right)$ converges exponentially fast to $P\left(\bm{\theta}_{0}\right)$ as $k$ grows. 3. 3.

The feedback matrix $G\left(\bm{\theta}_{0}\right)$ stabilizes the system:

[TABLE] 4. 4.

The minimum of the average cost (4) is achieved by $\bm{\pi^{\star}}$. 5. 5.

In the class of LTI policies (i.e. of the form $u\left(t\right)=Gx\left(t\right)$), the policy $\bm{\pi^{\star}}$ is the only optimal one.

In the remainder of the paper, we employ reinforcement learning algorithms to address Problem 1, studying the growth rates of $\mathcal{R}_{n}\left(\bm{\pi}\right)$. Similarly, letting $\widetilde{\bm{\theta}}_{n}$ be the learned/estimated parameter at time $n$ (the sample size is $n$ as well), we consider the exploration performance in Problem 2 through the rates of the learning error ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widetilde{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$. Bootstrap is the cornerstone of the proposed algorithms to efficiently randomize the design of the control inputs, and address the trade-off between the learning accuracy and the regret.

III Algorithms

An algorithm needs to address the common dilemma of decision making under uncertainty, as follows. First, if the algorithm makes decisions naively according to the estimated (learned) dynamics parameter, it will presumably fail to provide a small regret. Intuitively, the state $x\left(t\right)$ and the action $u\left(t\right)$ are required to be highly correlated in order to remain close to the optimal strategy $\bm{\pi^{\star}}$ in (6). Because of this correlation, history $\mathbb{H}_{t}$ may fail to accurately learn $\bm{\theta}_{0}$, which can lead to drastically large regret values. Technically, if $u\left(t\right)=Gx\left(t\right)$ for some $r\times p$ feedback gain matrix $G$, then the dimension of the observed history is effectively $p$, while the rows of the matrices $\bm{\theta}$ in the parameter space belong to $\mathbb{R}^{q}$. Therefore, learning can be dramatically misleading. This phenomenon of failing to falsify the imprecise approximations of the true model is extensively discussed in the adaptive control literature [41, 42, 12, 21].

In other words, if the policy fails to sufficiently explore the parameter space, an inaccurate approximation $\widetilde{\bm{\theta}}_{t}$ can falsely be treated as an accurate one. This necessitates an efficient exploration strategy to decrease the aforementioned correlation between the state and the action. Moreover, the above argument reveals the reasoning leading to UCB approaches [13, 14], or statistically independent dither schemes [16, 22], as useful prescriptions to overcome the exploration-exploitation dilemma.

In order to explore, the decision maker needs to deviate from the learned model $\widetilde{\bm{\theta}}_{t}$ prior to using $\mathcal{M}\left(\widetilde{\bm{\theta}}_{t}\right)$ to design the reinforcement learning policy. On the other hand though, the above deviations must be sufficiently small in order to avoid significant deterioration in the exploitation performance (i.e. increase in the regret). The solution we discuss here is to utilize the bootstrap to provide the necessary balance between these two competing objectives.

To this end, the policy $\bm{\pi}$ applies the supposedly optimal control action treating $\mathcal{M}\left(\widehat{\bm{\theta}}_{t}\right)$ as the true model, where $\widehat{\bm{\theta}}_{t}$ is provided by the bootstrap algorithm. It computes the regression residuals for the learned parameter $\widetilde{\bm{\theta}}_{t}$, and bootstraps (i.e. resamples) them to reconstruct a surrogate system. Then, the history of the surrogate system will be the data being used to compute $\widehat{\bm{\theta}}_{t}$. In the first subsection, we explain the least squares estimator for learning the model parameter, as well as the above residual bootstrap procedure.

Subsequently, in the second subsection we present an episodic algorithm which updates the model-based policy at the end of every episode, while the lengths of the episodes grow exponentially fast. Therefore, as the duration of the interaction with the system grows, the number of policy updates scales logarithmically. As a matter of fact, this leads to a significant reduction in the computation of the reinforcement learning policy, by avoiding unnecessary updates before collecting sufficient data, due to the fact that the solution of the algebraic Riccati equation (5) for a hypothetical model is not instantly available. The latter would impose a substantial computational burden, especially for systems whose dimension is fairly large.

III-A Residual Bootstrap

According to the linear dynamical model in (1), a natural procedure to learn $\bm{\theta}_{0}$ through the control input $u\left(t\right)$ and the observed states $x\left(t\right),x\left(t+1\right)$ is based on least squares. In the sequel, we discuss the residual bootstrap method for the least squares learning procedure. Further, we will present the corresponding algorithm which will be used as a subroutine in the reinforcement learning algorithm in the next subsection.

Recall that the LTI policy $\bm{\pi^{\star}}$ in (6) is optimal. Thus, a natural form of the adaptive policies that a reinforcement learning algorithm is expected to provide through planing according to the learned model, is $u\left(t\right)=G_{t}x\left(t\right)$. Assuming so for $t<n$, now the algorithm needs to decide about the action at time $n$. Thus, plugging $u\left(t\right)=G_{t}x\left(t\right)$ in the dynamical model (1), and denoting

[TABLE]

we get the so-called closed-loop evolution of the system by the (possibly time-varying) autoregressive dynamics

[TABLE]

for $0\leq t<n$. Then, Algorithm 1 returns the bootstrapped parameter $\widehat{\bm{\theta}}_{n}$ based on the matrices $\left\{F_{t}\right\}_{t=0}^{n-1}$, as well as the available state observations $\left\{x\left(t\right)\right\}_{t=0}^{n}$. The details are provided below.

First, based on the collected history $\left\{F_{t}\right\}_{t=0}^{n-1},\left\{x\left(t\right)\right\}_{t=0}^{n}$, define the following least square estimate of $\bm{\theta}_{0}$:

[TABLE]

The learning procedure (8) treats the noise vectors $\xi\left(t\right)$ as the errors of a linear regression procedure, based on the dynamical model (1). Therefore, the residuals of the least squares estimate are defined by the difference between the observed response $x\left(t+1\right)$, and the fitted response $\widetilde{\bm{\theta}}_{n}F_{t}x\left(t\right)$. That is,

[TABLE]

for $0\leq t<n$. The residuals $\left\{\zeta\left(t\right)\right\}_{t=1}^{n}$ can conceptually be considered as approximations of the actual regression errors $\left\{\xi\left(t\right)\right\}_{t=1}^{n}$. Using the residuals $\left\{\zeta\left(t\right)\right\}_{t=1}^{n}$, we define the centered empirical distribution

[TABLE]

where $\overline{\zeta}_{n}$, the average of the residuals given by

[TABLE]

is being used for centering the empirical distribution. In fact, $\widehat{\mathbb{P}}_{n}$ is the sample probability measure for the population distribution of the noise process $\left\{\xi\left(t\right)\right\}_{t=1}^{\infty}$. Note that $\widehat{\mathbb{P}}_{n}$ is defined on $\mathbb{R}^{p}$. We then use $\widetilde{\bm{\theta}}_{n}$ and $\widehat{\mathbb{P}}_{n}$ to generate the surrogate state vectors $\left\{\widehat{x}\left(t\right)\right\}_{t=0}^{n}$ by the dynamical model

[TABLE]

where the bootstrap noise vectors $\widehat{\xi}\left(t+1\right)$ are drawn independently from $\widehat{\mathbb{P}}_{n}$. Hence, letting $\widehat{\mathbb{E}}_{n}$ be the expectation with respect to $\widehat{\mathbb{P}}_{n}$, clearly we have $\widehat{\mathbb{E}}_{n}\left[\widehat{\xi}\left(t\right)\right]=0$. Also note that the actual dynamics parameter for the surrogate system $\left\{\widehat{x}\left(t\right)\right\}_{t=0}^{n}$ is the learned parameter $\widetilde{\bm{\theta}}_{n}$ defined in (8). Finally, the algorithm applies the least squares estimator to the generated surrogate states to obtain $\widehat{\bm{\theta}}_{n}$:

[TABLE]

The pseudo-code for the residual Bootstrap explained above is given in Algorithm 1. It will be used later at the heart of Algorithm 2 to design reinforcement learning policies.

Remark 1**.**

If the noise process is parametrized, one can accordingly draw $\widehat{\xi}\left(t\right)$ from the corresponding parametric sample distribution.

To see that, assume we know that the noise vectors belong to a parametric family of stochastic processes. Then, instead of using the nonparametric empirical distribution in (10), one can use the the residuals $\zeta\left(t\right)$ to estimate the parameter of interest. So, letting $\widehat{\mathbb{P}}_{n}$ be the parametric distribution provided by the obtained estimate, the bootstrap noise $\widehat{\xi}\left(t\right)$ can be sampled independently from $\widehat{\mathbb{P}}_{n}$. For example, if we know that $\xi\left(t\right)$ are i.i.d. Gaussian vectors, we can find the sample covariance matrix $\widehat{\Sigma}_{n}=n^{-1}\sum\limits_{t=1}^{n}\zeta\left(t\right)\zeta\left(t\right)^{\prime}-\overline{\zeta}_{n}\>\overline{\zeta}_{n}^{\prime}$, and draw $\widehat{\xi}\left(t\right)$ independently from the centered Gaussian distribution with covariance matrix $\widehat{\Sigma}_{n}$.

Remark 2**.**

In the original version of bootstrap [23], the covariates (i.e. the state vectors) are fixed, and only the residuals are being bootstrapped. In the time series models such as (1), every state vector comprises of the previous noise vectors. Therefore, bootstrapping the residuals automatically leads to new state sequence $\left\{\widehat{x}\left(t\right)\right\}_{t=0}^{\infty}$ for the surrogate system [30].

III-B Policy Design

Next, Algorithm 2 for decision making under uncertainty based on bootstrapping the residuals (Algorithm 1) is discussed. For this purpose, we first define the extended gain matrix $F\left(\bm{\theta}\right)$ based on the optimal feedback $G\left(\bm{\theta}\right)$.

Definition 2**.**

For parameter $\bm{\theta}\in\mathbb{R}^{p\times q}$, using the matrix $G\left(\bm{\theta}\right)$ in (5), define the $q\times p$ matrix $F\left(\bm{\theta}\right)=\left[I_{p},G\left(\bm{\theta}\right)^{\prime}\right]^{\prime}$.

The matrix $F\left(\bm{\theta}\right)$ can be interpreted as an extension of the original feedback gain; applying $u\left(t\right)=G\left(\bm{\theta}\right)x\left(t\right)$, the closed-loop transition matrix takes the form $A_{0}+B_{0}G\left(\bm{\theta}\right)=\bm{\theta}_{0}F\left(\bm{\theta}\right)$.

Recall that the true model is not known, and a reinforcement learning algorithm needs to simultaneously learn the dynamics parameter, and design the control input. To do so, we present an episodic decision making strategy outlined in Algorithm 2. That is, the policy applies control actions during each episode, assuming that the approximation of the model available at the time coincides with the true model. Then, at the end of every episode, the algorithm updates the learned model based on the history collected so far, and continues making decisions as if the new approximate model is the truth. The learning mentioned above is through a linear regression for the dynamics (1), and the approximation consists of bootstrapping (by Algorithm 1) the model estimate obtained by the regression. In the sequel, we explain the details of the above alternating steps of the algorithm.

The reinforcement learning policy is initiated with the history $\mathbb{H}_{0}$ in the first line of Algorithm 2. Then, it chooses an arbitrary stabilizable approximation of $\bm{\theta}_{0}$, denoted by $\widehat{\bm{\theta}}_{0}$, and starts the system by applying the action prescribed by the model $\mathcal{M}\left(\widehat{\bm{\theta}}_{0}\right)$; i.e. $u\left(t\right)=G\left(\widehat{\bm{\theta}}_{0}\right)x\left(t\right)$. Note that selection of $\widehat{\bm{\theta}}_{0}$ is straightforward, since almost all (w.r.t. Lebesgue measure) parameter matrices are stabilizable [20].

The starting time-points of the episodes are determined by the exponents of the reinforcement rate $\beta>1$. That is, at every time $t=\lceil\beta^{m}\rceil$, the approximation $\widehat{\bm{\theta}}_{t}$ will be updated, while for $\beta^{m}\leq t<\beta^{m+1}$ the algorithm freezes $\widehat{\bm{\theta}}_{t}$. In other words, whenever $t=\lceil\beta^{m}\rceil$, Algorithm 2 calls the residual bootstrap Algorithm 1 to get $\widehat{\bm{\theta}}_{t}$ according to the collected history of the control actions and the states. So, for all $\beta^{m}\leq t<\beta^{m+1}$, the matrices $F\left(\widehat{\bm{\theta}}_{t}\right)$ are exactly the same. The efficiency of the policy relies on the idea that the sequence $\left\{\widehat{\bm{\theta}}_{t}\right\}_{t=0}^{\infty}$ will provide finer approximations of the truth $\bm{\theta}_{0}$, as the algorithm proceeds (or more precisely, as $m$ grows).

IV Theoretical Results and Simulations

We start by establishing performance guarantees on the regret and the learning accuracy for bootstrap-based policies, supplemented by numerical examples that illustrate the behavior of Algorithm 2 for both identification and regulation. The following result specifies the growth rate of the regret $\mathcal{R}_{n}\left(\bm{\pi}\right)$ for the policy $\bm{\pi}$ designed by Algorithm 2, as well as the decay rate of the identification error ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widetilde{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$.

Theorem 1**.**

Letting $\bm{\pi}$ be the policy given by Algorithm 2, define the learned parameter $\widetilde{\bm{\theta}}_{n}$ by (8). Then, we have

[TABLE]

The proof is provided in the appendix. Technically, it relies on the careful examination of the effect of Algorithm 1 on the randomization of the feedback gains $G\left(\widehat{\bm{\theta}}_{\lceil\beta^{i}\rceil}\right)$. This randomization in turn diversifies the extended gain matrices $\left\{F\left(\widehat{\bm{\theta}}_{\lceil\beta^{i}\rceil}\right)\right\}_{i=1}^{m}$, so that their superposition efficiently explores the whole parameter space $\mathbb{R}^{p\times q}$, as $m$ grows. To this end, we utilize the state-of-the-art results on the behavior of the algebraic Riccati equation [14], properties of the optimality manifold [43, 44, 21], and results from martingale theory [45, 46, 47], limit distributions of dependent sequences [36, 37], and the bootstrap [35].

The regulation and identification rates of Theorem 1 are modulo logarithmic factors similar to the corresponding rates of the reinforcement learning policies utilizing OFU [13, 14], additive randomization [21], posterior sampling [20], and input perturbation [16]. Moreover, the square root scaling of the regret is efficient for adaptive regulation of LQ systems as discussed next.

Recalling the discussion at the beginning of Section III, an adaptive control policy needs to sufficiently explore the parameter space in order to balance the trade-off of identification and regulation. For falsifying the imprecise approximation $\widetilde{\bm{\theta}}_{n}$ through an exploration procedure, the control signals $\left\{u\left(t\right)\right\}_{t=0}^{n-1}$ need to deviate from the optimal feedback gain $G\left(\bm{\theta}_{0}\right)$. More precisely, for a policy $\bm{\pi}$, let the deviations from the optimal feedback be $\epsilon_{t}={\left|\kern-1.29167pt\left|u\left(t\right)-G\left(\bm{\theta}_{0}\right)x\left(t\right)\right|\kern-1.29167pt\right|}$, for $0\leq t<n$. Then, observing the history $\mathbb{H}_{n}$, the error of estimating the true dynamics parameter $\bm{\theta}_{0}$ is at least $\sigma_{n}$ (modulo a constant factor), where $\sigma_{n}^{-2}=\sum\limits_{t=0}^{n-1}\epsilon_{t}^{2}$ [7, 10]. Hence, if $\bm{\pi}$ aims to falsify $\widetilde{\bm{\theta}}_{n}$, the difference ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widetilde{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$ needs to be in the order of magnitude at least $\sigma_{n}$ [12, 21]. Whenever $\bm{\pi}$ employs $\widehat{\bm{\theta}}_{n}$ for designing control inputs, $\liminf\limits_{n\to\infty}\sigma_{n}^{-1}{\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widehat{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}>0$ holds, since $\widehat{\bm{\theta}}_{n}$ needs to be found through $\widetilde{\bm{\theta}}_{n}$.

On the other hand, for the above deviations we have

[TABLE]

according to the regret specification recently established [21]. Further, applying the adaptive feedback $u(n)=G\left(\widehat{\bm{\theta}}_{n}\right)x(n)$ at time $n$, the increase $\mathcal{R}_{n+1}\left(\bm{\pi}\right)-\mathcal{R}_{n}\left(\bm{\pi}\right)$ in the regret is approximately ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|G\left(\widehat{\bm{\theta}}_{n}\right)-G\left(\bm{\theta}_{0}\right)\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}^{2}$ [16], which is up to a constant factor at least ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|{\widehat{\bm{\theta}}_{n}}-{\bm{\theta}_{0}}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}^{2}$ [14, 21]. Thus, the lower bound $\sigma_{n}$ for ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widehat{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$ implies that $\liminf\limits_{n\to\infty}\sigma_{n}^{-2}\left(\mathcal{R}_{n+1}\left(\bm{\pi}\right)-\mathcal{R}_{n}\left(\bm{\pi}\right)\right)>0$. Putting the latter result and (16) together, we obtain $\liminf\limits_{n\to\infty}\left(\mathcal{R}_{n+1}\left(\bm{\pi}\right)^{2}-\mathcal{R}_{n}\left(\bm{\pi}\right)^{2}\right)>0$, which provides the lower bound $\liminf\limits_{n\to\infty}n^{-1/2}\mathcal{R}_{n}\left(\bm{\pi}\right)>0$. Note that a rigorous proof of the above lower bound argument is beyond the scope of this paper. For more detailed discussions, we refer the reader to the aforementioned references.

IV-A Numerical Illustration

Next, we present numerical analyses employing Algorithm 2 for decision-making under uncertainty. Henceforth, let $\bm{\pi}$ be the reinforcement learning policy provided by Algorithm 2, with reinforcement rate $\beta=1.2$. The true dynamical model and cost matrices are provided in (13) and (14), respectively. Solving the algebraic Riccati equation, we get $P\left(\bm{\theta}_{0}\right),G\left(\bm{\theta}_{0}\right)$ given in (15), which lead to a closed-loop matrix of the spectral radius $\left|\bm{\lambda}_{\max}\left(\bm{\theta}_{0}F\left(\bm{\theta}_{0}\right)\right)\right|=0.26$.

Figure 1 depicts the normalized regret as a function of $n$, for $100$ replicates of the stochastic linear system in (1). The corresponding normalized identification errors are plotted in Figure 2. These figures are in a full agreement with the theoretical result of Theorem 1; both normalized rates $n^{-1/2}\mathcal{R}_{n}\left(\bm{\pi}\right)$ and $n^{1/4}{\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widetilde{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$ are dominated by logarithmic factors of the time index $n$. In Figure 3, we plot the resulting spectral radius of the reinforcement learning policy $\bm{\pi}$ for both the actual system of the dynamics parameter $\bm{\theta}_{0}$, as well as that of the surrogate system of $\widetilde{\bm{\theta}}_{n}$. According to Figure 3, Algorithm 2 fully stabilizes the system, even though in the first few episodes the system is unstable.

The ensuing figures indicate the robustness of Algorithm 2 to structural breaks. Figure 4 shows the performance of $n^{-1/2}\mathcal{R}_{n}\left(\bm{\pi}\right)$ and $n^{1/4}{\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widetilde{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$ for a single break in the model, wherein at time $t=400$ the dynamics matrices suddenly become

[TABLE]

A similar performance analysis while the system incurs two breaks is provided in Figure 5. The first break is similar to the one mentioned above, and occurs at time $t=200$. Then, for the second break at time $t=700$, the true dynamics matrices change to

[TABLE]

According to Figures 4 and 5, Algorithm 2 is robust to remarkably large values of model mis-specifications. Note that the reinforcement learning policy is fully ignorant of the breaks. So, $\bm{\pi}$ needs to adaptively adjust its decision-making law toward the new optimal policies.

The rationale for the exhibited robustness is as follows: when a break occurs, the parameter estimate $\widetilde{\bm{\theta}}$ becomes an inaccurate approximation of the true dynamics matrix $\bm{\theta}_{0}$. This in turn leads the regression residuals $\left\{\zeta\left(t\right)\right\}_{t=1}^{n}$ becoming large. Therefore, the bootstrapped parameter $\widehat{\bm{\theta}}$ computed by Algorithm 1 provides a large randomization, which in turn leads to an increase in the exploration phase. Then, after a few episodes, the resulting enhanced exploration provides more accurate estimates $\widetilde{\bm{\theta}}$, and the above negative feedback procedure proceeds. Thus, as time grows, Algorithm 2 self-tunes to the equilibrium of the suitable amount of exploration.

The above argument intuitively indicates that the aforementioned equilibrium is a stable one. Since the endogenous randomization of the bootstrap procedure consistently assesses the accuracy of the fitted model $\widetilde{\bm{\theta}}$, the resulting adaptive policy automatically adjusts the old decision-making strategy to the new environment. Hence, the algorithm accordingly addresses the unexpected flaw of the sudden and unknown changes in the true model $\mathcal{M}\left(\bm{\theta}_{0}\right)$, as well as the resulting unpredicted deviations in the trajectory of the state sequence $\left\{x\left(t\right)\right\}_{t=0}^{\infty}$.

V Concluding Remarks

We proposed a reinforcement learning algorithm for sequential decision-making for an LQ system with unknown temporal dynamics. The presented model-based policy is based on residual bootstrap, and is shown to be efficient in terms of both identification and regulation. Namely, we establish the rates for the worst-case regret, as well as the learning accuracy. Further, we discussed the robustness of the bootstrap method for handling unexpected changes in the dynamical model.

As the first work on bootstrap-based policies for LQ models, it poses a number of interesting questions. For example, theoretical analysis for addressing the performance of bootstrap method under imperfect observation is a natural direction for future work. Further subjects of interest include design and analysis of fully non-parametric randomization methods such as covariate resampling. Finally, extending the presented framework to model-free algorithms can be considered as another fruitful research direction to examine.

Appendix A Auxiliary Results

Next, we present auxiliary results being used in the proof of Theorem 1 in Appendix B. First, we present the regret bounds of Lemma 1, for which the proof can be found in the work of Faradonbeh et al. [21], using martingale convergence analysis of Lai and Wei [48]. Then, Lemma 2 provides the local specification of the optimality manifold which is established for both full-rank [43, 44] and rank-deficient dynamics matrices [21]. Subsequently, Lemma 3 is presented to study convergence rates of linear regression procedures. The asymptotic [45] and non-asymptotic [46, 47] proofs of Lemma 3 are available in the literature.

Then, we state Lemma 4 which addresses the behavior of the empirical covariance matrix of the state sequence of a stabilized system [8]. Further, we have Lemma 5 that is providing the local Lipschitz continuity property (i.e. in a neighborhood of $\bm{\theta}_{0}$) of the feedback matrix $G\left(\bm{\theta}\right)$ [14]. Finally, we establish Lemma 6 on the population covariance matrix induced by the empirical probability measure $\widehat{\mathbb{P}}_{n}$ defined in (10).

Lemma 1**.**

For the sequence of $r\times p$ matrices $\left\{G_{t}\right\}_{t=0}^{\infty}$, suppose that there is a filtration $\left\{\mathcal{F}_{t}\right\}_{{t}=0}^{\infty}$ such that for all $t\geq 0$, $x\left(t\right),G_{t},\xi\left(t\right)$ are $\mathcal{F}_{t}$-measurable, and \mathbb{E}\left[\xi\left(t+1\right)\big{|}\mathcal{F}_{t}\right]=0. Then, for the regret of the policy

[TABLE]

the following holds:

[TABLE]

Lemma 2**.**

For the stabilizable parameter $\bm{\theta}_{1}=\left[A_{1},B_{1}\right]$, let $\mathcal{S}\left({\bm{\theta}_{1}}\right)$ be the manifold of optimal feedback gains:

[TABLE]

Then, the tangent space of $\mathcal{S}\left({\bm{\theta}_{1}}\right)$ at point $\bm{\theta}_{1}$ consists of matrices $\left[M,N\right]$, such that $M\in\mathbb{R}^{p\times p},N\in\mathbb{R}^{p\times r}$ satisfy

[TABLE]

where $D_{1}=\bm{\theta}_{1}F\left(\bm{\theta}_{1}\right)$, $Z=P\left(\bm{\theta}_{1}\right)\left(M+NG\left(\bm{\theta}_{1}\right)\right)$.

Lemma 3**.**

Consider the dynamical system $x\left(t+1\right)=\bm{\theta}_{0}F_{t}x\left(t\right)+\xi\left(t+1\right)$. Then, define

[TABLE]

where $U_{n}^{-1}$ is the Moore-Penrose inverse. It holds that

[TABLE]

Further, since $\widetilde{\bm{\theta}}_{n}$ provided by (8) satisfies the normal equation $\widetilde{\bm{\theta}}_{n}U_{n}=\sum\limits_{t=0}^{n-1}x\left(t+1\right)x\left(t\right)^{\prime}F_{t}^{\prime}$, we have

[TABLE]

Therefore, we get

[TABLE]

Lemma 4**.**

Suppose that the control feedback matrix $G\left(\widehat{\bm{\theta}}\right)$ is applied to the dynamical model $\mathcal{M}\left(\bm{\theta}_{0}\right)$. Hence, plugging $u\left(t\right)=G\left(\widehat{\bm{\theta}}\right)x\left(t\right)$ in (1), the system evolves according to $x\left(t+1\right)=\bm{\theta}_{0}F\left(\widehat{\bm{\theta}}\right)x\left(t\right)+\xi\left(t+1\right)$. Assuming $\left|\bm{\lambda}_{\max}\left(\bm{\theta}_{0}F\left(\widehat{\bm{\theta}}\right)\right)\right|<1$, for the empirical covariance matrix $V_{n}=\sum\limits_{t=0}^{n}x\left(t\right)x\left(t\right)^{\prime}$ the following holds:

[TABLE]

Lemma 5**.**

Letting $C_{P}<\infty$, there is a constant $C_{G}<\infty$ such that

[TABLE]

Lemma 6**.**

Let $\widehat{\mathbb{E}}_{n}$ be the expectation with respect to $\widehat{\mathbb{P}}_{n}$; the empirical probability measure of the residuals defined in (10). Then, for the Bootstrap covariance matrix $\widehat{\Sigma}_{n}=\widehat{\mathbb{E}}_{n}\left[\widehat{\xi}\left(t\right)\widehat{\xi}\left(t\right)^{\prime}\right]$ we have

[TABLE]

Proof.

First, the definition of $\left\{\zeta\left(t\right)\right\}_{t=1}^{n}$ in (9), in addition to the dynamics (1) yield $\zeta\left(t+1\right)=\left(\bm{\theta}_{0}-\widetilde{\bm{\theta}}_{n}\right)F_{t}x\left(t\right)+\xi\left(t+1\right)$. Defining $U_{n},W_{n}$ similar to Lemma 3, the normal equation $\widetilde{\bm{\theta}}_{n}U_{n}=\sum\limits_{t=0}^{n-1}x\left(t+1\right)x\left(t\right)^{\prime}F_{t}^{\prime}$ implies that

[TABLE]

So, we obtain

[TABLE]

since $\left(\bm{\theta}_{0}-\widetilde{\bm{\theta}}_{n}\right)U_{n}\left(\bm{\theta}_{0}-\widetilde{\bm{\theta}}_{n}\right)^{\prime}=W_{n}W_{n}^{\prime}$. Thus, applying the Law of Large Numbers to the matrices $\xi\left(t\right)\xi\left(t\right)^{\prime}$, we get $\limsup\limits_{n\to\infty}\left|\bm{\lambda}_{\max}\left(\widehat{\Sigma}_{n}\right)\right|\leq\left|\bm{\lambda}_{\max}\left(\Sigma\right)\right|$. Further, by the Law of Large Numbers, $n^{-1}\sum\limits_{t=1}^{n}\xi\left(t\right)$ vanishes as $n$ grows. Therefore, the definitions of $U_{n},W_{n}$ lead to

[TABLE]

Finally, since Lemma 3 implies that $\limsup\limits_{n\to\infty}n^{-1}\left|\bm{\lambda}_{\max}\left(W_{n}W_{n}^{\prime}\right)\right|=0$, we get the desired result on the smallest eigenvalue: $\liminf\limits_{n\to\infty}\left|\bm{\lambda}_{\min}\left(\widehat{\Sigma}_{n}\right)\right|\geq\left|\bm{\lambda}_{\min}\left(\Sigma\right)\right|$. $\blacksquare$

Appendix B Proof of Theorem 1

The following analysis rigorously studies the behavior of both $\mathcal{R}_{n}\left(\bm{\pi}\right)$ and ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widetilde{\bm{\theta}}_{n}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$ as the time of interacting with the system, $n$, grows. In the sequel, we assume that the system is stable. The stabilization problem has been addressed previously in the literature  [19]. In fact, an ad-hoc algorithm for stabilizing the system is presented and analyzed in the work of Faradonbeh et al. [19]. It establishes high probability guarantees for stabilization in finite time, using random feedback gains. The proposed method can be executed a priori, and terminates in a relatively short time period. Moreover, the random feedback framework presented for stabilization algorithms can also be implemented with the bootstrap method of Algorithm 1. Therefore, the subsequent theoretical analysis focuses on establishing regret bounds and learning accuracy of Algorithm 2 after the transient stabilization period.

In the reinforcement learning policy provided by Algorithm 2, the residual bootstrap procedure of Algorithm 1 is being called at the end of every episode. Fixing $i$, for $j<i$ define the following quantities when the algorithm BOOTSTRAP is called at time $t=\lceil\beta^{i}\rceil$:

[TABLE]

Note that the notation is slightly overloaded in the above expressions since in every call, BOOTSTRAP generates a completely new set of surrogate noise and state vectors $\widehat{\xi}\left(\cdot\right),\widehat{x}\left(\cdot\right)$. Similarly, for the original system of dynamics parameter $\bm{\theta}_{0}$ define the matrices

[TABLE]

Using the martingale Central Limit Theorem, $\beta^{-i/2}\bm{\phi}_{i}$ converges in distribution to a Gaussian random matrix. Moreover, for an arbitrary $j<i$, let $\mathcal{F}_{j}$ be the sigma-field generated by $\left\{\xi\left(t\right)\right\}_{t=1}^{\lceil\beta^{j}\rceil}$ and $\left\{\left\{\widehat{\xi}\left(t\right)\right\}_{t=1}^{\lceil\beta^{\ell}\rceil}\right\}_{\ell=1}^{j}$. Then, for the $\mathcal{F}_{j}$-measurable $v\in\mathbb{R}^{p}$, Lemma 4 and Lemma 6 imply that as $k$ grows,

[TABLE]

where “$\Rightarrow$” denotes the convergence in distribution, and $\inf\limits_{{\left|\kern-0.90417pt\left|v\right|\kern-0.90417pt\right|}=1}\left|\bm{\lambda}_{\min}\left(\Sigma_{v}\right)\right|>0$. Next, we study eigenvalues of the matrix $\widehat{U}_{i}$. According to Lemma 4 and Lemma 6, we have:

[TABLE]

In order to establish a lower bound for the smallest eigenvalue of $\widehat{U}_{i}$, compute ${\widehat{U}_{i}}^{-1}$ according to the blocks of

[TABLE]

So,

[TABLE]

where

[TABLE]

Clearly, $\left|\bm{\lambda}_{\min}\left(X_{22}(i)\right)\right|\geq\left|\bm{\lambda}_{\min}\left(\widehat{U}_{i}\right)\right|$. Further, let $v=\left[v_{1}^{\prime},v_{2}^{\prime}\right]^{\prime}\in\mathbb{R}^{q}$ be arbitrary, where $v_{1}\in\mathbb{R}^{p}$, $v_{2}\in\mathbb{R}^{r}$. Then, we have

[TABLE]

So, Lemma 4 and Lemma 6 lead to

[TABLE]

On the other hand, the least squares estimates in (8), (12) lead to

[TABLE]

Therefore, applying Lemma 3, Lemma 4, and Lemma 6, we obtain

[TABLE]

Combine the above two inequalities to get

[TABLE]

Now, we compare $G\left(\widehat{\bm{\theta}}_{\lceil\beta^{i}\rceil}\right)$, $G\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)$ based on the local specification of the optimality manifolds in Lemma 2. To do so, we also need to characterize the right-hand-side $p\times r$ sub-matrix of the randomization matrix $\widehat{\bm{\theta}}_{\lceil\beta^{i}\rceil}-\widetilde{\bm{\theta}}_{\lceil\beta^{i}\rceil}$ provided by the residual bootstrap procedure. Therefore, multiplying both sides of the normal equation $\widehat{\bm{\theta}}_{\lceil\beta^{i}\rceil}-\widetilde{\bm{\theta}}_{\lceil\beta^{i}\rceil}=\widehat{W}_{i}\widehat{U}_{i}^{-1/2}$ by $\left[0_{r\times p},I_{r}\right]^{\prime}$, and using $F\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)^{\prime}\widehat{U}_{i}^{-1}\left[0_{r\times p},I_{r}\right]^{\prime}=\left(-X_{12}(i)+G\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)^{\prime}\right){X_{22}(i)}^{-1}$, we obtain

[TABLE]

where $\bm{\psi}_{i}=\sum\limits_{j=0}^{i-1}\bm{\phi}_{ij}\left(G\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)^{\prime}-X_{12}(i)\right){X_{22}(i)}^{-1/2}$. Letting $v\in\mathbb{R}^{q}$ be an eigenvector of $\widehat{U}_{i}$, decompose it to $v=\left[v_{1}^{\prime},v_{2}^{\prime}\right]^{\prime}$, where $v_{1}\in\mathbb{R}^{p}$, $v_{2}\in\mathbb{R}^{r}$. Then, in order to consider $v^{\prime}\widehat{U}_{i+1}v$, (19) indicates that we need to investigate

[TABLE]

for $j<i$. For this purpose, plug in $\left(\widehat{\bm{\theta}}_{\lceil\beta^{i}\rceil}-\widetilde{\bm{\theta}}_{\lceil\beta^{i}\rceil}\right)\left[0_{r\times p},I_{r}\right]^{\prime}=\bm{\psi}_{i}{X_{22}(i)}^{-1/2}$, as well as (22), in Lemma 2. Using the notation of Lemma 2, ${\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}}F\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)$ is stable, and for the matrix $M+NG\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)$, (22) implies that $\limsup\limits_{j<i;\>\>j\to\infty}i^{-1/2}\beta^{j/2}{\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\left(\widehat{\bm{\theta}}_{\lceil\beta^{i}\rceil}-\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)F\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}<\infty$. Hence

[TABLE]

By leveraging (18), (19), and (20), we get

[TABLE]

Note that since the system is stabilized, $\left|\bm{\lambda}_{\min}\left(Q_{u}\right)\right|>0$ implies that $\sup\limits_{1\leq k}{\left|\kern-1.29167pt\left|\kern-1.29167pt\left|G\left(\widehat{\bm{\theta}}_{\lceil\beta^{k}\rceil}\right)\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}<\infty$. Further, according to Lemma 4 we have $\limsup\limits_{i\to\infty}{\left|\kern-1.29167pt\left|\kern-1.29167pt\left|X_{12}(i)\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}<\infty$. Thus, using (21), as long as $\limsup\limits_{i\to\infty}i\beta^{-i}\left|\bm{\lambda}_{\min}\left(\widehat{U}_{i}\right)\right|=0$, we have

[TABLE]

Therefore, (23) yields

[TABLE]

In addition, since the extended feedback matrices $F\left(\widehat{\bm{\theta}}_{\lceil\beta^{j}\rceil}\right)$ in the definitions of $U_{i}$ and $\widehat{U}_{i}$ are the same, we get a similar result for the actual system:

[TABLE]

Next, applying Lemma 3 to the least squares estimates in (8), (12), and using (24), (25), we obtain

[TABLE]

Combining the above two inequalities, we get the desired result regarding the learning accuracy. Further, Lemma 5 implies that

[TABLE]

So, according to Lemma 4 we get

[TABLE]

On the other hand, the moment condition $\sup\limits_{t\geq 1}\mathbb{E}\left[{\left|\kern-1.29167pt\left|\xi\left(t\right)\right|\kern-1.29167pt\right|}^{\alpha}\right]<\infty$ implies that for all $\gamma>1/\alpha$, $\epsilon>0$,

[TABLE]

which by the Borel-Cantelli Lemma leads to $\limsup\limits_{t\to\infty}t^{-1/4}{\left|\kern-1.29167pt\left|\xi\left(t\right)\right|\kern-1.29167pt\right|}=0$, since $\alpha>4$. By stability of the closed-loop matrix $A_{0}+B_{0}G\left(\bm{\theta}_{0}\right)$, it holds that

[TABLE]

Putting (26), (27), and Lemma 1 together, the desired result on the growth rate of the regret is established, which finishes the proof.