remarkRemark hypothesisHypothesis claimClaim assumptionAssumption A Model-Free
First-Order Method for LQRC. Ju, G. Kotsalis, and G. Lan
[][nocite]ex_supplement

# A Model-Free First-Order Method for Linear Quadratic Regulator with
$\tilde{O}(1/\varepsilon)$ Sampling Complexity††thanks: Initially released on
arXiv on Dec 1st, 2022. This research was partially supported by NIFA grant
2020-67021-31526 and NSF grant 1909298. CJ is supported by the Department of
Energy Computational Science Graduate Fellowship under Award Number DE-
SC0022158.

Caleb Ju H. Milton Stewart School of Industrial & Systems Engineering, Georgia
Institute of Technology, Atlanta, GA<EMAIL_ADDRESS><EMAIL_ADDRESS>Georgios Kotsalis Amazon Core AI, Seattle, WA ().
<EMAIL_ADDRESS>Guanghui Lan22footnotemark: 2

###### Abstract

We consider the classic stochastic linear quadratic regulator (LQR) problem
under an infinite horizon average stage cost. By leveraging recent policy
gradient methods from reinforcement learning, we obtain a first-order method
that finds a stable feedback law whose objective function gap to the optima is
at most $\varepsilon$ with high probability using $\tilde{O}(1/\varepsilon)$
samples, where $\tilde{O}$ hides polylogarithmic dependence on $\varepsilon$.
Our proposed method seems to have the best dependence on $\varepsilon$ within
the model-free literature without the assumption that all policies generated
by the algorithm are stable almost surely, and it matches the best-known rate
from the model-based literature, up to logarithmic factors. The improved
dependence on $\varepsilon$ is achieved by showing the accuracy scales with
the variance rather than the standard deviation of the gradient estimation
error. Our developments that result in this improved sampling complexity fall
in the category of actor-critic algorithms. The actor part involves a
variational inequality formulation of the stochastic LQR problem, while in the
critic part, we utilize a conditional stochastic primal-dual method and show
that the algorithm has the optimal rate of convergence when paired with a
shrinking multi-epoch scheme.

###### keywords:

linear quadratic regulator, actor-critic, primal-dual methods, variational
inequality

93C05, 65K05

##  1 Introduction

The linear quadratic regulator (LQR) problem is a classic problem in control
theory studied extensively since the work of Kalman [17, 18]. Of interest is
the infinite horizon linear-quadratic stochastic control problem. Techniques
for solving this problem go back to the 1960s. It is well known that when this
system is completely controllable, the optimal feedback law is a linear static
state feedback of $u_{t}^{*}=-K^{*}x_{t}$ for some gain matrix
$K^{*}\in\mathbb{R}^{k\times n}$, where $x_{t}\in\mathbb{R}^{n}$ is the
current state. The matrix $K^{*}$ can be obtained by solving the algebraic
Ricatti equations (ARE), which can be numerically solved by, e.g., iterative
methods or linear matrix inequalities [1, 2, 4, 13]. When the parameters of
the system are unknown, estimates can be obtained using system identification
[27] and robust control [44].

In recent years, methods from machine learning have been used for optimal
control, in particular, for solving LQR. For example, one can estimate the
system parameters using a regularized least-squares estimator. The estimated
system is viewed as the nominal system and then solved by ARE. This is known
as the certainty equivalence principle [7, 30], which is a model-based method.
In contrast, there are model-free methods, which aim to solve the problem in
an “online” manner without explicitly estimating the system parameters. Many
model-free methods are inspired by problems in reinforcement learning (RL)
[36]. RL tackles problems such as LQR using optimization methods that only
require an oracle that can simulate the costs and evolution of the system.
There are a plethora of model-free methods, and most can be categorized as
either approximate dynamic programming [21, 37] or gradient-based (first- or
zeroth-order) methods [24, 25, 35].

We focus on gradient-based methods. It is known LQR is non-convex (the set of
stable feedback laws is non-convex [11]), so generally speaking at best we can
only find locally optimal solutions. The ubiquitous gradient descent has been
used in practice [39] while its convergence to the local optimal solution was
proved in [41]. It was not until the breakthrough result of [11] that showed
gradient descent converges to the global optimal solution by proving the
objective function of LQR satisfies the so-called Polyak-Łojasiewicz-condition
[28, 33], or PŁ-condition.

In this paper, we develop a first-order method to efficiently solve LQR. We
use an actor-critic method [19]. The outer loop, or the actor step, views the
LQR problem as a variational inequality, or VI, that satisfies the generalized
monotonicity condition [9, 20]. Solving the VI requires an estimate of the
control Lyapunov equation, which is solved by the critic. The critic estimates
this quantity to high accuracy via a primal-dual method for solving a
stochastic fixed-point problem. We enhance this result using a shrinking
multi-epoch approach to match the optimal convergence rate of [20] for solving
stochastic VIs under Markovian noise, which includes fixed-point problems.
Overall, our method finds a $K$ such that
$J(K)-\min_{K^{*}}J(K^{*})\leq\varepsilon$ with probability at least
$1-\delta$ and sampling complexity of
$O(\varepsilon^{-1}(\ln\varepsilon^{-1})^{2}(\ln(1/(\varepsilon\delta))^{7/2}))$.
Our main insight to achieving this complexity is to show the accuracy
$\varepsilon$ scales with the variance rather than the standard deviation of
the gradient estimation error.

It seems our result has improved dependence on $\varepsilon$ (up to
logarithmic dependence on $\varepsilon$) when compared to the best rates known
in the literature thus far in the model-free case. The work of [42] used a
similar actor-critic approach to achieve a sampling complexity of
$O(\varepsilon^{-5})$. A matching
$O(\varepsilon^{-1}(\ln\varepsilon^{-1})^{2})$ complexity (up to logarithmic
factors) was shown for convergence in expectation [45], but it relies on the
assumption regarding the almost sure stability of the policy at every
iteration, which may not be realistic since the underlying noise of the linear
dynamics can introduce errors that yield unstable policies with probability
greater than zero. On the other hand, we only assume the initial policy is
stable. We seem to avoid this assumption by obtaining high probability
results. Within the model-based case, a least-squares estimator can estimate
the dynamics using $O(\varepsilon^{-2})$ samples [7]. This was later improved
to $O(\varepsilon^{-1})$ in [30] by noticing the error scales as the square of
the estimation error. Our proposed method shows the model-free case can be
competitive with the model-based case in terms of the dependence on
$\varepsilon$. One may prefer model-free since model-based methods can suffer
from model bias, where insufficient samples and unaccounted-for uncertainties
can lead to poorly fitted models and unstable policies [7, 8].

The paper is organized as follows. In Section 2, we formulate the LQR
optimization problem, explore its properties, and cover related works. In
Section 3, we reformulate LQR as a variational inequality and show the
accompanying algorithm. This algorithm requires an accurate solution to the
control Lyapunov equation, which is solved in Section 4 by using a shrinking
multi-epoch stochastic primal-dual method to achieve the optimal deterministic
and stochastic error rate. We finish with some preliminary numerical
experiments in Section 5.

##  2 Preliminaries, problem formulation, and related works

###  2.1 Notation

We refer to the norm $\|\cdot\|$ as the standard $\ell_{2}$-norm for vectors
and the induced $\ell_{2}$-norm for matrices. The norm $\|\cdot\|_{F}$ is the
Frobenius norm for a given matrix. The notation $\rho(\cdot)$ refers to the
spectral radius of a given matrix. For symmetric matrices $A$ and $B$, we use
the standard Loewner order of $A-B\succeq 0$ to mean that $A-B$ is
semidefinite and $A-B\succ 0$ to mean that it is positive definite. For any
symmetric matrix $X$, let $\mathrm{svec}(X)$ be the vectorization of the upper
triangular part of $X$ with the off-diagonals scaled by $\sqrt{2}$.

###  2.2 Linear systems and quadratic control

We consider a discrete-time, time-invariant linear dynamical system with
disturbances,

(1) $x_{t+1}=Ax_{t}+Bu_{t}+w_{t},\quad t\in\mathbb{Z}_{+},$

where $x_{t}\in\mathbb{R}^{n}$ is the state, $u_{t}\in\mathbb{R}^{k}$ is the
control input, and $w_{t}$ is the stochastic disturbance respectively. It is
assumed that the random variables $\\{x_{0},w_{0},w_{1},\ldots,\\}$ are
jointly independent and Gaussian,
$x_{0}\sim\mathcal{N}(0,\Sigma_{0}),~{}~{}w_{t}\sim\mathcal{N}(0,\Psi),$ where
the covariance matrices $\Sigma_{0}\succ 0,\Psi\succ 0$ are known. The
infinite horizon time-average cost to be minimized over casual policies is

(2)
$J=\lim_{T\rightarrow\infty}\frac{1}{T}\mathbb{E}\bigg{[}\sum_{t=0}^{T-1}x^{\mathrm{T}}_{t}Qx_{t}+u^{\mathrm{T}}_{t}Ru_{t}\bigg{]},$

where $Q\succ 0$ and $R\succ 0$. We assume that the matrices $(A,B)$ form a
controllable pair. This is a classic problem in control theory whose solution
can be reviewed in [6, 22]. The average cost is minimized by a linear state
feedback policy

$u_{t}=-K^{*}x_{t}=-(R+B^{\mathrm{T}}PB)^{-1}B^{\mathrm{T}}PAx_{t},\quad
t\in\mathbb{Z}_{+},$

with $P$ being the unique positive definite solution to the algebraic Riccati
equation

$P=Q+A^{\mathrm{T}}PA-A^{\mathrm{T}}PB(R+B^{\mathrm{T}}PB)^{-1}B^{\mathrm{T}}PA.$

The focus of this paper is on the model-free case, and we will restrict
ourselves to randomized linear state feedback policies of the form

(3) $u_{t}=-Kx_{t}+v_{t}.$

The additive noise to the control input is also Gaussian, with
$v_{t}\sim\mathcal{N}(0,\sigma^{2}I),$ for some sufficiently large known
$\sigma>0$. We also require that the basic random variables
$\\{x_{0},w_{0},v_{0},w_{1},v_{1},\ldots,\\}$ are jointly independent. The
resulting closed loop dynamics under (3) are

(4) $\displaystyle x_{t+1}=(A-BK)x_{t}+w_{t}+Bv_{t}.$

The additional feedback law disturbance, $v_{t}$, is often included to allow
for exploration [7, 21, 42]. Later, we will also see the positive noise level
$\sigma>0$ ensures we can accurately estimate the gradient. We denote the set
of stable policies by

(5) $\mathcal{S}=\\{K\in\mathbb{R}^{k\times n}~{}|~{}\rho(A-BK)<1\\}.$

We will search in an iterative fashion for a feedback gain $K\in\mathcal{S}$
that yields an approximation to the optimum within a prespecified degree of
fidelity. We write for the cost $J(K)$ to emphasize its parametric
dependence111The case $\sigma=0$ is a classical result found in textbooks [6],
while $\sigma>0$ appears to be recently developed from the reinforcement
learning literature. A proof can be found in [42, Proposition 3.1]. and state
the aforementioned policy optimization problem as

(6)
$\min_{K}J(K)=\begin{cases}\mathrm{Tr}[(Q+K^{\mathrm{T}}RK)\Sigma_{K}+\sigma^{2}R]=\mathrm{Tr}[P_{K}(\Psi+\sigma^{2}BB^{\mathrm{T}})],&K\in\mathcal{S}\\\
\infty,&\text{o/w}\end{cases}.$

Here, the matrix $\Sigma_{K}$ is the solution to the Lyapunov equation

(7)
$\Sigma_{K}=(\Psi+\sigma^{2}BB^{\mathrm{T}})+(A-BK)\Sigma_{K}(A-BK)^{\mathrm{T}},$

and it corresponds to the steady state covariance matrix of the zero mean
state vector. The matrix $P_{K}$ is the unique positive definite solution to
the Lyapunov equation

(8) $P_{K}=Q+K^{\mathrm{T}}RK+(A-BK)^{\mathrm{T}}P_{K}(A-BK).$

The following lemma contains an expression for the policy gradient $\nabla
J(K)$.

######  Lemma 2.1.

For a given $K\in\mathcal{S}$ the gradient of the cost function $J(K)$ is

(9) $\nabla
J(K)=2\big{[}(R+B^{\mathrm{T}}P_{K}B)K-B^{\mathrm{T}}P_{K}A\big{]}\Sigma_{K}.$

The proof of this lemma can be traced to [10] with minor adjustments to
account for the fact that we are working in discrete time. For the remainder
of this paper, we write the policy gradient compactly as

$\nabla J(K)=2E_{K}\Sigma_{K},$

where the matrix $E_{K}$, the natural gradient, is given by

(10) $E_{K}=(R+B^{\mathrm{T}}P_{K}B)K-B^{\mathrm{T}}P_{K}A.$

The name natural gradient is from the reinforcement learning literature [11,
16].

If we can accurately estimate the gradient, we can characterize the change in
objective value when we descend along the negative of this direction. The
following result, also known as the performance difference lemma, quantifies
this statement.

######  Lemma 2.2.

For any pair of stable policies $K^{\prime}$ and $K$,

$\displaystyle J(K^{\prime})-J(K)$
$\displaystyle=2\mathrm{Tr}(\Sigma_{K^{\prime}}(K^{\prime}-K)^{T}E_{K})+\mathrm{Tr}(\Sigma_{K^{\prime}}(K^{\prime}-K)^{T}(R+B^{T}P_{K}B)(K^{\prime}-K)).$

A proof can be found in [11, Lemma 6]. The above result is akin to the
smoothness of a function, where the change in $J(K^{\prime})-J(K)$ can be
upper bounded by a linear plus quadratic term. By choosing an appropriate
$K^{\prime}$ w.r.t. to a policy $K$, the above result derives the so-called
PŁ-condition or gradient domination condition [28, 33].

######  Lemma 2.3.

For any stable policy $K_{t}$,

$\displaystyle\frac{\sigma_{\min}(\Psi)}{\|R+B^{T}P_{K_{t}}B\|}\|E_{K_{t}}\|^{2}_{F}\leq
J(K_{t})-J(K^{*})\leq\frac{\|\Sigma_{K^{*}}\|}{\sigma_{\min}(R)}\|E_{K_{t}}\|^{2}_{F}.$

A proof can be readily deduced from [10]. In view of Lemma 2.1, an accurate
estimate of $E_{K}$ requires an estimate of the solution $P_{K}$ from the
Lyapunov equation (8) as well as the system parameters $R$, $B$, and $A$.
Alternatively, one can estimate $E_{K}$ without such knowledge. To do so, we
introduce the differential cost
$V:\mathcal{S}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$,

$V(K,x)=\mathbb{E}[\sum_{t=0}^{\infty}\big{(}x_{t}^{\mathrm{T}}(Q+K^{\mathrm{T}}RK)x_{t}-J(K)\big{)}|x_{0}=x],$

while the corresponding $Q$-function is

$\displaystyle Q(K,x,u)$
$\displaystyle=\mathbb{E}[\sum_{t=1}^{\infty}x_{t}^{\mathrm{T}}\big{(}(Q+K^{\mathrm{T}}RK)x_{t}-J(K)\big{)}|x_{0}=x,u_{0}=u]$

Note that $Q(K,x,-Kx)=V(K,x)$ and the explicit formulas are

$\displaystyle V(K,x)$ $\displaystyle=$ $\displaystyle
x^{\mathrm{T}}P_{K}x-\mathrm{Tr}[P_{K}\Sigma_{K}],$ $\displaystyle Q(K,x,u)$
$\displaystyle=$ $\displaystyle
x^{\mathrm{T}}\Theta(K)_{11}x+x^{\mathrm{T}}\Theta(K)_{12}u+u^{\mathrm{T}}\Theta(K)_{21}x+u^{\mathrm{T}}\Theta(K)_{22}u-\mathrm{Tr}[P_{K}\Sigma_{K}],$

where

(16) $\displaystyle\Theta(K):=\left[\begin{array}[]{c}A^{\mathrm{T}}\\\
B^{\mathrm{T}}\end{array}\right]P_{K}\left[\begin{array}[]{cc}A&B\end{array}\right]+\left[\begin{array}[]{cc}Q&0\\\
0&R\end{array}\right].$

The above decomposition makes evident the fact that $\Theta(K)\succ 0$. The
basis of the online learning algorithm lies in the relationship

(17) $\nabla
J(K)=2E_{K}\Sigma_{K}=2\big{[}\Theta(K)_{22}K-\Theta(K)_{21}\big{]}\Sigma_{K},$

where $\Theta(K)_{ij}$ is the $(i,j)$-th block of $\Theta_{K}$ if we view it
as a $2\times 2$ block matrix. The relative $Q$ function satisfies the Bellman
equation

$\displaystyle J(K)+Q(K,x,u)$ $\displaystyle=$ $\displaystyle
x^{\mathrm{T}}Qx+u^{\mathrm{T}}Ru+\mathbb{E}[Q(K,x_{+},-Kx_{+})|x,u]$
$\displaystyle=$ $\displaystyle
x^{\mathrm{T}}Qx+u^{\mathrm{T}}Ru+\mathbb{E}[V(K,x_{+})|x,u]$

where given $(x,u)$ the random variable $x_{+}$ satisfies the system dynamics

$x_{+}=Ax+Bu+w,\quad w\sim\mathcal{N}(0,\Psi).$

For a given state-action pair $z:=\begin{bmatrix}x&u\end{bmatrix}$, define

(20)
$\displaystyle\phi(z)=\mathrm{svec}(zz^{T}),\quad\theta(K)=\mathrm{svec}(\Theta(K)),\quad\vartheta(K)=\left[\begin{array}[]{c}J(K)\\\
\theta(K)\end{array}\right].$

Denote the stationary measure on the state and control space w.r.t. to the
current gain matrix $K$ as $\Pi_{K}$. We define a variable metric via the
inner product

$\langle\phi_{i}(x,u),\phi_{j}(x,u)\rangle=\mathbb{E}_{\Pi_{K}}[\phi_{i}(x,u)\phi_{j}(x,u)].$

It follows from the Bellman equation for the differential $Q$-function that

(21) $H\vartheta(K)=b,$

where the matrix and vector are defined as

(22) $H=\left[\begin{array}[]{cc}1&0\\\
\mathbb{E}_{\Pi_{K}}[\phi(x,u)]&\Xi_{K}\end{array}\right],\ \
b=\left[\begin{array}[]{c}\mathbb{E}_{\Pi_{K}}[x^{\mathrm{T}}Qx+u^{\mathrm{T}}Ru]\\\
\mathbb{E}_{\Pi_{K}}[(x^{\mathrm{T}}Qx+u^{\mathrm{T}}Ru)\phi(x,u)\end{array}\right],$

and

$\Xi_{K}=\mathbb{E}_{\Pi_{K}}[\phi(x,u)(\phi(x,u)-\phi(x_{+},-Kx_{+}))^{\mathrm{T}}].$

###  2.3 Related Works

Two prominent paradigms of LQR being studied are the random initialization and
noisy dynamics model [29]. The random initialization assumes only the initial
state is random, and the remaining state transitions are deterministic. The
noisy dynamics model, which we study in this paper, considers both a random
initialization and perturbations in each state transition.

Within the noisy dynamics setting, an infinite-horizon average cost was
studied by [42] using an actor-critic method. Around the same time, [21] used
an approximate policy iteration from the RL literature. More recently, [43,
45] utilize a stochastic gradient descent-type analysis to solve LQR. While
they can get similar convergence rates, their results rely on a strong
assumption that $\rho(A-BK_{t})$ is almost surely less than one for every $t$,
where $K_{t}$ is the gain matrix at iteration $t$. The aforementioned methods
are model-free. As for model-based methods, [7] used a least-squares estimator
for the transition dynamics and then use tools from system-level synthesis to
obtain their controller. [30] later improved the sampling complexity by
noticing the function gap scales as the square of the parameter error.

Similar problems include the infinite-horizon discounted objective value. The
PhD thesis of [5] investigated the setting where there was no noise. In the
noisy dynamic setting, [12] developed algorithms that find near-optimal
solutions with high probability guarantees. The authors of [37] studied the
sampling complexity of a least-squares temporal difference in this setting. A
question of minimizing a discount factor of one on the random initialization
setting was first shown to be solvable to global optimality by the
breakthrough paper of [11]. Since their work, improvements in the sampling
complexity have been proposed in the works of [29] and [31]. It should be
mentioned these last two works explored the use of zeroth-order methods. This
setting, while model-free, requires a generative model, where one can start
running the linear dynamics from any initial state.

##  3 The actor: Policy optimization via a variational inequality

Our goal is to solve the optimization problem

$\displaystyle\min_{K\in\mathcal{S}}J(K),$

where recall $\mathcal{S}$ is the set of stable policies (c.f. Eq. 5). Since
it is known the optimal policy $K^{*}$ lies in the interior of the open set
$\mathcal{S}$ [10], a necessary condition of optimality is $\nabla
J(K^{*})=2E_{K^{*}}\Sigma_{K^{*}}=0$. By observing that the covariance matrix
from (7) satisfies $\Sigma_{K}\succ 0$, we can solve the variational
inequality (VI):

$\displaystyle\text{find }K^{*}\text{ s.t. }\langle
E_{K^{*}},K-K^{*}\rangle\geq 0,\hskip 5.0pt\forall K\in\mathbb{R}^{k\times
n}.$

The solution to this VI problem is related to our original LQR problem thanks
to the performance difference lemma. It satisfies a generalized monotonicity
condition [9]. Choosing $K^{\prime}=K^{*}$ in Lemma 2.2 and rearranging terms,

$\displaystyle 2\mathbb{E}_{x\sim\mathcal{N}(0,\Sigma_{K^{*}})}\big{[}\langle
xx^{T}E_{K},K-K^{*}\rangle\big{]}\geq J(K)-J(K^{*})\geq 0.$

This VI shows the connection of our method to the VI formulations in
reinforcement learning [24, 26, 25]. By exploiting the monotonicity condition,
we write an algorithm to solve the VI in Algorithm 1. This method is
equivalent to the natural policy gradient method seen in the reinforcement
learning literature. Each iteration requires an estimate of the natural
gradient $E_{K_{t}}$, which we recall requires the solution $P_{K_{t}}$ to the
Lyapunov equation (8), as defined by the linear system of equations in (21).
In the next section, we provide a primal-dual method to complete this task.

Algorithm 1 Natural Policy Gradient Descent

Input: $K_{0}$, $\eta\in\mathbb{R}_{++}$

for $t=0,\ldots,T-1$ do

Estimate ${E}_{K_{t}}$ with (17) and (20) by solving (21) using Algorithm 3

$K_{t+1}=\underset{K\in\mathbb{R}^{k\times n}}{\mathrm{argmin}}\\{\eta\langle
2E_{K_{t}},K\rangle+\frac{1}{2}\|K-K_{t}\|_{F}^{2}\\}$

end for

return $K_{T}$

Let us now move on to proving the convergence of solving LQR using Algorithm
1. The following result will be helpful to bound the objective function.

######  Lemma 3.1.

When $K$ is a stable policy, then

$\displaystyle\|\Sigma_{K}\|\leq\mathrm{Tr}(\Sigma_{K})\leq\sigma_{\min}(Q)^{-1}J(K)\text{
and }\|P_{K}\|\leq\mathrm{Tr}(P_{K})\leq\sigma_{\min}(\Psi)^{-1}J(K).$

Also, $(1-\rho^{2})^{-1}\leq
J(K)\sigma_{\min}(\Psi)^{-1}\sigma_{\min}(Q)^{-1}$ where
$\rho\equiv\rho(A-BK)<1$.

######  Proof 3.2.

The bounds on $\Sigma_{K}$ and $P_{K}$ follow from Eq. 6 and basic trace
inequalities. The bound on $(1-\rho^{2})^{-1}$ is from the upper bound on
$\mathrm{Tr}(\Sigma_{K})$ and

$\displaystyle\frac{\sigma_{\min}(\Psi)}{1-\rho^{2}}$
$\displaystyle=\sum\limits_{t=0}^{\infty}\sigma_{\min}(\Psi)\rho^{2t}$
$\displaystyle\leq\sum\limits_{t=0}^{\infty}\mathrm{Tr}\big{(}(A-BK)^{t}(\Psi+\sigma^{2}BB^{T})((A-BK)^{t})^{T}\big{)}=\mathrm{Tr}(\Sigma_{K}),$

where
$\mathrm{Tr}(M^{T}M(N+N^{\prime}))\geq\sigma_{\min}(N)\mathrm{Tr}(M^{T}M)\geq\sigma_{\min}(N)\|M\|^{2}\geq\sigma_{\min}(N)\rho(M)^{2}$
for a square matrix $M$ and positive semidefinite matrices $N$ and
$N^{\prime}$.

To ease notation in our analysis, we define the following positive constants:

(23)
$\begin{split}C_{1}&:=2[\|R\|+2\sigma_{\min}^{-1}(\Psi)\|B\|^{2}J(K_{0})]\\\
C_{2}&:=\sigma_{\min}(\Psi)\sigma_{\min}(Q)\sigma_{\min}(R)\\\
C_{3}&:=J(K_{0})/\sigma_{\min}(Q).\end{split}$

We also define the natural gradient error,

(24) $\displaystyle\delta_{t}=E_{K_{t}}^{\star}-{E}_{K_{t}},$

where $E_{K_{t}}^{\star}$ is the true natural gradient and ${E}_{K_{t}}$ is
the stochastic estimate. Note that the solution to $K_{t+1}$ in Algorithm 1
can be explicitly written as $K_{t+1}=K_{t}-2\eta_{t}E_{K_{t}}$. We will now
show a contraction property for the LQR problem.

######  Proposition 3.3.

Let $K_{t}$ be a stable controller and suppose $J(K_{t})\leq 2J(K_{0})$. At
iteration $t$ of Algorithm 1, choose the constant step size

$\displaystyle\eta_{t}\equiv\eta=1/(2C_{1}),$

and let $K_{t+1}$ be the computed (stochastic) policy. If $K_{t+1}$ is stable,
then

(25) $\displaystyle
J(K_{t+1})-J(K^{*})\leq\Big{(}1-\frac{C_{2}}{2J(K^{*})C_{1}}\Big{)}[J(K_{t})-J(K^{*})]+\epsilon_{t},$

where

(26)
$\displaystyle\epsilon_{t}:=\frac{1}{C_{1}}\Big{[}\|\delta_{t}\|_{F}^{2}\|\Sigma_{K_{t+1}}\|+\sigma_{\min}(\Psi)\mathrm{Tr}([E_{K_{t}}^{\star}]^{T}\delta_{t})\Big{]}.$

######  Proof 3.4.

By the performance difference lemma (Lemma 2.2),

$\displaystyle J(K_{t+1})-J(K_{t})$
$\displaystyle=-4\eta_{t}\mathrm{Tr}(\Sigma_{K_{t+1}}{E}_{K_{t}}^{T}\delta_{t})-4\eta_{t}\mathrm{Tr}(\Sigma_{K_{t+1}}{E}_{K_{t}}^{T}{E}_{K_{t}})$
$\displaystyle\hskip
11.0pt+4\eta_{t}^{2}\mathrm{Tr}(\Sigma_{K_{t+1}}E_{K_{t}}^{T}(R+B^{T}P_{K_{t}}B)E_{K_{t}})$
$\displaystyle\leq
2\eta_{t}\|E_{K_{t}}S^{T}\|_{F}^{2}+2\eta_{t}\|\delta_{t}S^{T}\|_{F}^{2}-(4\eta_{t}-4\eta_{t}^{2}\|R+B^{T}P_{K_{t}}B\|)\mathrm{Tr}(\Sigma_{K_{t+1}}E_{K_{t}}^{T}E_{K_{t}})$
$\displaystyle=2\eta_{t}\|\delta_{t}S^{T}\|_{F}^{2}-(2\eta_{t}-4\eta_{t}^{2}\|R+B^{T}P_{K_{t}}B\|)\mathrm{Tr}(\Sigma_{K_{t+1}}E_{K_{t}}^{T}E_{K_{t}})$
$\displaystyle\leq\frac{1}{C_{1}}\|\delta_{t}\|_{F}^{2}\|\Sigma_{K_{t+1}}\|-\frac{1}{2C_{1}}\mathrm{Tr}(\Sigma_{K_{t+1}}E_{K_{t}}^{T}E_{K_{t}})$
$\displaystyle\leq\frac{1}{C_{1}}\|\delta_{t}\|_{F}^{2}\|\Sigma_{K_{t+1}}\|-\frac{\sigma_{\min}(\Psi)}{2C_{1}}\mathrm{Tr}(E_{K_{t}}^{T}E_{K_{t}}),$

where the matrix $S$ satisfies $\Sigma_{K_{t+1}}=S^{T}S$ (since
$\Sigma_{K_{t+1}}$ is symmetric positive semidefinite, or sym-psd). In the
first inequality, we used Cauchy-Schwarz and Young’s inequality to show
$\mathrm{Tr}(M^{T}N)\leq\|M\|_{F}\|N\|_{F}\leq\frac{\|M|\|_{F}^{2}}{2}+\frac{\|N\|_{F}^{2}}{2}$
for arbitrary matrices $M$ and $N$ and
$\mathrm{Tr}(MN)\leq\|M\|\mathrm{Tr}(N)$ when $N$ is sym-psd. In the second
inequality, we used the identity
$\|\delta_{t}S^{T}\|^{2}_{F}=\mathrm{Tr}(\Sigma_{K_{t+1}}\delta_{t}^{T}\delta_{t})$
as well as

(27)
$\displaystyle\|R+B^{T}P_{K_{t}}B\|\leq\|R\|+\|B\|^{2}\|P_{K_{t}}\|\leq\|R\|+\sigma_{\min}^{-1}(\Psi)J(K_{t})\|B\|^{2}\leq
C_{1}/2,$

which is due to $J(K_{t})\leq 2J(K_{0})$ and Lemma 3.1. In the third
inequality, we used $\mathrm{Tr}(MN)\geq\sigma_{\min}(M)\mathrm{Tr}(N)$ when
$N$ is sym-psd and $\Sigma_{K_{t+1}}\succeq\Psi\succ 0$ (c.f. (7)). Noting
$\mathrm{Tr}(\delta_{t}^{T}\delta_{t})\geq 0$, we have
$\mathrm{Tr}(E_{K_{t}}^{T}E_{K_{t}})\geq\mathrm{Tr}([E_{K_{t}}^{\star}]^{T}E_{K_{t}}^{\star})-2\mathrm{Tr}([E_{K_{t}}^{\star}]^{T}\delta_{t})$.
Altogether,

$\begin{split}&J(K_{t+1})-J(K_{t})\\\
&\leq-\frac{\sigma_{\min}(\Psi)}{2C_{1}}\mathrm{Tr}([E_{K_{t}}^{\star}]^{T}E_{K_{t}}^{\star})+\frac{1}{C_{1}}\Big{[}\|\delta_{t}\|_{F}^{2}\|\Sigma_{K_{t+1}}\|+\sigma_{\min}(\Psi)\mathrm{Tr}([E_{K_{t}}^{\star}]^{T}\delta_{t})\Big{]}\\\
&\leq-\frac{C_{2}}{2J(K^{*})C_{1}}[J(K_{t})-J(K^{*})]+\epsilon_{t}.\end{split}$

where in the last line we used the PŁ-condition (Lemma 2.3) followed by the
bound on $\|\Sigma_{K_{*}}\|$ (Lemma 3.1). Adding $J(K_{t})-J(K^{*})$ to both
sides finishes the proof.

The above result has two technical issues. First, we assumed $K_{t+1}$ to be
stable. Second, we need to bound $\|\Sigma_{K_{t+1}}\|\|\delta_{t}\|^{2}_{F}$
by $\epsilon$, but the two multipliers are correlated. We will address both
these issues next using the following lemmas, which establish the stability of
the policy if we knew the true natural gradient and bounds for
$\|\Sigma_{K_{t+1}}\|$. Both results can be found in [11] as Lemma 15 and 16,
respectively.

######  Lemma 3.5.

Let $K_{t}$ be stable and suppose $E_{K_{t}}$ is computed exactly in Algorithm
1. If the step size is not too large, i.e.,

$\displaystyle\eta_{t}\leq\|R+B^{T}P_{K_{t}}B\|^{-1}$

then the subsequent policy $K_{t+1}$ is stable.

Under the assumption $J(K_{t})\leq 2J(K_{0})$ for all $t$, it can be shown our
choice of step size from Proposition 3.3 satisfies this bound thanks to Lemma
3.1.

######  Lemma 3.6.

For any stable controller $K$, let $K^{\prime}$ be an arbitrary controller
where

$\displaystyle\|K^{\prime}-K\|\leq\frac{\sigma_{\min}(Q)\sigma_{\min}(\Sigma_{K})}{4J(K)\|B\|(\|A-BK\|+1)}.$

Then $K^{\prime}$ is stable, and moreover,

$\displaystyle\|\Sigma_{K^{\prime}}-\Sigma_{K}\|\leq
4\Big{(}\frac{J(K)}{\sigma_{\min}(Q)}\Big{)}^{2}\frac{\|B\|(\|A-BK\|+1)}{\sigma_{\min}(\Sigma_{K})}\|K^{\prime}-K\|.$

In the above lemma, we set $K$ as
$K_{t+1}^{\star}=K_{t}+2\eta_{t}E_{K_{t}}^{\star}$ to be the policy with the
correct gradient. It is stable with an appropriate step size. The policy
$K_{t+1}=K_{t}+2\eta_{t}E_{K_{t}}$ is the stochastic estimate produced during
Algorithm 1.

######  Proposition 3.7.

Suppose the assumptions from Proposition 3.3 hold. Also suppose
$\|K_{t+1}-K_{t+1}^{\star}\|$ satisfies the hypothesis of Lemma 3.6. Then the
policy $K_{t+1}$ computed by Algorithm 1 satisfies

$\displaystyle
J(K_{t+1})-J(K^{*})\leq\Big{(}1-\frac{C_{2}}{4J(K^{*})C_{1}}\Big{)}[J(K_{t})-J(K^{*})]+\epsilon_{t}^{\prime},$

where

(28) $\displaystyle\epsilon_{t}^{\prime}$
$\displaystyle=\frac{1}{C_{1}}\bigg{[}\frac{\sigma_{\min}(\Psi)C_{1}J(K^{*})}{2C_{2}}\cdot\|\delta\|_{F}^{2}+2C_{3}\|\delta_{t}\|_{F}^{2}+\frac{32C_{3}^{2}\|B\|}{\sigma_{\min}(\Psi)C_{1}}\|\delta_{t}\|^{3}_{F}\bigg{]}.$

######  Proof 3.8.

We start with the decomposition,
$\|\Sigma_{K_{t+1}}\|\leq\|\Sigma_{K_{t+1}^{\star}}\|+\|\Sigma_{K_{t+1}}-\Sigma_{K_{t+1}^{\star}}\|$.
We can bound the second summand with Lemma 3.6. Thus, using the assumption
$J(K_{t})\leq 2J(K_{0})$ and choice in $C_{3}$, we find

$\displaystyle\|\Sigma_{K_{t+1}}\|$
$\displaystyle\leq\|\Sigma_{K_{t+1}^{\star}}\|+4(2C_{3})^{2}\frac{\|B\|(\|A-BK_{t+1}^{\star}\|+1)}{\sigma_{\min}(\Sigma_{K_{t+1}^{\star}})}\|K_{t+1}-K_{t+1}^{\star}\|$
$\displaystyle\leq
2C_{3}+16C_{3}^{2}\frac{2\|B\|}{\sigma_{\min}(\Psi)C_{1}}\|\delta_{t}\|,$

where in the last line, the bound on the first term comes from combining Lemma
3.1 and $J(K_{t+1}^{\star})\leq J(K_{t})\leq 2J(K_{0})$ (thanks to Proposition
3.3). The bound on the second term comes from recalling the natural gradient
error (24), that $K_{t+1}^{\star}$ is a stable policy, and
$\Sigma_{K}\succeq\Psi\succ 0$ from (7). Thus, we have a bound on
$\|\Sigma_{K_{t+1}}\|$, which is a part of the first term in $\epsilon_{t}$
from (26). It remains to bound the second term in $\epsilon_{t}$:
$\sigma_{\min}(\Psi)\mathrm{Tr}([E_{K_{t}}^{\star}]\delta_{t})/C_{1}$. By
Young’s inequality,

$\displaystyle\mathrm{Tr}([E_{K_{t}}^{\star}]^{T}\delta_{t})$
$\displaystyle\leq\frac{C_{2}}{4J(K^{*})\|R+B^{T}P_{K_{t}}B\|}\cdot\|E_{K_{t}}^{\star}\|_{F}^{2}+\frac{J(K^{*})\|R+B^{T}P_{K_{t}}B\|}{C_{2}}\cdot\|\delta_{t}\|_{F}^{2}$
$\displaystyle\leq\frac{C_{2}}{4\sigma_{\min}(\Psi)J(K^{*})}\cdot[J(K_{t})-J(K^{*})]+\frac{J(K^{*})C_{1}}{2C_{2}}\cdot\|\delta_{t}\|_{F}^{2},$

where in the last line we used the PŁ-condition (Lemma 2.3) and (27) to bound
the two summands. We finish by applying both bounds above back into (25).

We can now prove convergence results for Algorithm 1. To simplify notation,
define the condition number as $\kappa=8J(K^{*})C_{1}/C_{2}$.

######  Theorem 3.9.

Suppose $\kappa\geq 1/2$ and $K_{0}$ is a stable policy. Let
$l:=\lceil\kappa\rceil$ and fix an error tolerance $\varepsilon>0$. If at
every iteration $t=0,\ldots,T-1$ we use the step size from Proposition 3.7 and
have $\|\delta_{t}\|_{F}^{2}\leq\min\\{C_{4},C_{5},C_{6}\\}$ in Algorithm 1,
where

(29)
$\begin{split}C_{4}&:=\Big{(}\frac{\sigma_{\min}(Q)\sigma_{\min}(\Psi)C_{1}}{16J(K_{0})\|B\|}\Big{)}^{2},\hskip
5.0ptC_{5}:=\Big{(}\frac{\sigma_{\min}(\Psi_{\sigma})C_{2}}{3840J(K_{0})\|B\|C_{3}^{2}}\varepsilon\Big{)}^{2/3},\\\
C_{6}&:=\frac{C_{2}}{120J(K_{0})C_{1}}\min\Big{(}\frac{C_{2}}{\sigma_{\min}(\Psi)J(K^{*})},\frac{1}{2C_{3}}\Big{)}\varepsilon,\\\
\end{split}$

then the following holds for every $t=0,\ldots,T$:

  1. 1.

Stability: $K_{t}$ is a stable policy.

  2. 2.

Monotonicity or convergence: $J(K_{t})\leq J(K_{t-1})$ or
$J(K_{t})-J(K^{*})\leq\varepsilon$.

  3. 3.

Linear rate: $J(K_{t})-J(K^{*})\leq 2^{-\lfloor t/l\rfloor}J(K_{0})$ as long
as $t\leq l\cdot\log_{2}(J(K_{0})/\varepsilon)$.

######  Proof 3.10.

We use mathematical induction, where the base case of $t=0$ holds by our
assumptions. Now consider some $t+1\leq T$, and suppose the claim holds for
all iterations up to $t$. We first show $K_{t+1}$ is a stable policy. We have

$\displaystyle\|K_{t+1}-K_{t+1}^{\star}\|$
$\displaystyle\leq\frac{\sigma_{\min}(Q)\sigma_{\min}(\Psi)}{16J(K_{0})\|B\|}\leq\frac{\sigma_{\min}(Q)\sigma_{\min}(\Sigma_{K_{t}})}{4J(K_{t})\|B\|(\|A-BK_{t}\|+1)},$

where the first inequality is from $\|\delta_{t}\|^{2}_{F}\leq C_{4}$ and the
second is by $\sigma_{\min}(\Sigma_{K_{t}})\geq\sigma_{\min}(\Psi)$ from (7)
as well as both $J(K_{t})\leq J(K_{0})+J(K^{*})\leq 2J(K_{0})$ and
$\|A-BK_{t}\|<1$ from the inductive hypothesis. Therefore, the hypothesis for
Lemma 3.6 is satisfied, and $K_{t+1}$ is a stable policy.

Next, we prove monotonicity or convergence. Suppose
$J(K_{t})-J(K^{*})>\varepsilon$. By the inductive hypothesis, $K_{t}$ is
stable and $J(K_{t})\leq 2J(K_{0})$. So invoking Proposition 3.7 and adding
$J(K_{*})-J(K_{t})$,

(30) $\displaystyle J(K_{t+1})-J(K_{t})$
$\displaystyle\leq-\frac{C_{2}}{4J(K^{*})C_{1}}[J(K_{t})-J(K^{*})]+\epsilon_{t}^{\prime}$

where we recall $\epsilon_{t}^{\prime}$ is defined in (28). By our bounds
$\|\delta_{t}\|_{F}^{2}\leq\min\\{C_{5},C_{6}\\}$, we ensure

(31)
$\epsilon_{t}^{\prime}\leq\frac{C_{2}\varepsilon}{40J(K_{0})C_{1}}\leq\frac{C_{2}\varepsilon}{40J(K^{*})C_{1}},$

where upon plugging back into (30) and noticing
$J(K_{t})-J(K^{*})>\varepsilon$, we deduce $J(K_{t+1})\leq J(K_{t})$.

Finally, let us show linear rate of convergence. To simplify our notation, let
us define the constant $Z=\frac{C_{2}}{5C_{1}}$. We already showed (31), which
is equivalent to $\epsilon_{t}^{\prime}\leq\frac{Z\varepsilon}{8J(K_{0})}$.
Recalling that we consider iterations $t\leq
l\cdot\log_{2}(J(K_{0})/\varepsilon)$, this implies $\epsilon_{t}^{\prime}\leq
2^{-(\lfloor t/l\rfloor-2)}(Z/2)$. Thus, using our bound on
$\epsilon_{t}^{\prime}$ and choice of $l$, the recursion (30) can be
simplified into (see, e.g., [23, Lemma 11])

$\displaystyle J(K_{t+1})-J(K^{*})$ $\displaystyle\leq
2^{-\lfloor(t+1)/l\rfloor}\Big{(}J(K_{0})-J(K^{*})+\frac{5l\cdot
Z}{16}\Big{)}\leq 2^{-\lfloor(t+1)/l\rfloor}J(K_{0}),$

where the last inequality is because $\kappa\geq 1/2$, so
$l/2=\lceil\kappa\rceil/2\leq\kappa$. This completes the proof of convergence
and finishes our proof by induction.

The theorem suggests the number of policy optimization iterations is

$\displaystyle\frac{16J(K^{*})[\|R\|+\sigma_{\min}^{-1}(\Psi)\|B\|^{2}J(K_{0})]\ln\big{(}\frac{J(K_{0})}{\varepsilon}\big{)}}{\sigma_{\min}(\Psi)\sigma_{\min}(Q)\sigma_{\min}(R)}$

to ensure $J(K_{T})-J(K^{*})\leq\varepsilon$. An important byproduct of the
theorem is a uniform bound on the spectral radius of $A-BK_{t}$. The following
result is due to Lemma 3.1 combined with the upper bound on $J(K_{t})$ from
the previous theorem. We assume $\varepsilon\leq J(K_{0})-J(K^{*})$ since
otherwise $J(K_{0})-J(K^{*})\leq\varepsilon$.

######  Corollary 3.11.

If the assumptions in Theorem 3.9 hold for iterations $t=0,\ldots,T-1$, then

$\displaystyle\rho_{T}:=\max_{0\leq t\leq
T}\rho(A-BK_{t})\leq\sqrt{1-\frac{\sigma_{\min}(\Psi)\sigma_{\min}(Q)}{J(K_{0})}}.$

Boundedness of $\rho_{T}$ will be used in the next section of policy
evaluation, where we ensure the error
$\|\delta_{t}\|_{F}^{2}\equiv\|E_{K_{t}}-E_{K_{t}}^{\star}\|_{F}^{2}$ is
small, as required by the previous theorem.

##  4 The critic: Policy evaluation via stochastic primal-dual

Given a controller a stable policy $K$, recall from Eq. 21 we can estimate the
gradient by solving the linear system of equations,

(32) $\mathbb{E}_{\Pi_{K}}[\tilde{H}]x=\mathbb{E}_{\Pi_{K}}[\tilde{b}],$

where $\tilde{H}$ and $\tilde{b}$ are the stochastic estimates of $H$ and $b$
as defined in Eq. 22. We re-formulate the above problem into the general min-
max problem of

(33) $\displaystyle\min_{x\in X}\big{\\{}f(x):=\max_{y\in Y}\\{\langle
y,\mathbb{E}_{\Pi_{K}}[\tilde{H}]x-\mathbb{E}_{\Pi_{K}}[\tilde{b}]\rangle-\hat{h}(y)\\}\big{\\}},$

where $X\subseteq\mathbb{R}^{n}$ and $Y\subseteq\mathbb{R}^{n}$ are closed
sets and $\hat{h}:Y\to\mathbb{R}$ is a relatively simple, proper, convex,
lower semicontinuous function that is easy to solve.

To facilitate our analysis, let us define some useful notation. Let
$\omega:X\to\mathbb{R}$ be some strongly convex function of modulus
$\sigma_{\omega}$, which we assume to be 1 by scaling. We define the Bregman
divergence associated with $X$ (or $Y$) and $\omega$ as

$\displaystyle
V(x,y)=\omega(y)-\omega(x)-\langle\omega^{\prime}(x),y-x\rangle,\hskip
5.0pt\forall x,y\in X.$

We define the diameter of $X$ (or $Y$) as $D_{X}^{2}=\max_{x,y\in X}V(x,y)$.
Let the random variables $\xi_{1},\xi_{2},\ldots,\xi_{t}$ (for any $t\geq 1$)
be a stochastic process generated according to some probability measure
$(\Omega,\mathcal{F},\mathbb{P})$ on a state space $\mathcal{S}$, e.g.,
$\mathcal{S}\equiv\mathbb{R}^{n}$. Define
$\mathcal{F}_{t}=\sigma(\xi_{1},\xi_{2},\ldots,\xi_{t})$ as the
$\sigma$-algebra generated by the first $t$ random variables. We may also
write the stochastic estimates as $\tilde{H}_{t}=\tilde{H}(\xi_{t})\text{ and
}\tilde{b}_{t}=\tilde{b}(\xi_{t})$, to indicate their parameterizations by the
random variable at time $t$.

The first assumption generalizes both the bounded and a looser “light-tail”
assumptions [32, 23] on the distribution of the stochastic estimates of the
gradient. There exists
$M_{H},M_{b},\varsigma_{H},\varsigma_{b},\beta\in\mathbb{R}_{+}$ such that for
any $t\geq 1$ and $\delta\in(0,1)$,

(34) $\displaystyle\|\tilde{H}_{t}\|\leq
M_{H}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\text{ and
}\|\tilde{b}_{t}\|\leq M_{b}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\text{
w.p. }1-\delta,$
$\displaystyle\|\tilde{H}_{t}-H\|\leq\varsigma_{H}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\text{
and
}\|\tilde{b}_{t}-b\|\leq\varsigma_{b}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\text{
w.p. }1-\delta.$

This assumption is related to sub-Gaussian random variables with two
differences. First, the numerator inside the natural log above is
$\mathrm{exp}(1)$ (instead of $1$) to ensure a support of positive reals
(e.g., $\|\tilde{H}\|\geq M_{H}$ w.p. 1). Second, the exponent of $\beta$
permits a different moment of the random variable to be sub-Gaussian. Later,
Lemma 4.9 and Lemma 4.12 will specify the constants.

Next, we show we can reduce the bias sufficiently fast. Conditioned on the
event that (34) holds, there exist $C,O_{H},O_{b}\in\mathbb{R}_{+}$ and
$\rho,\delta\in(0,1)$ such that for any $t\in\mathbb{Z}_{++}$ and
$\tau\in\mathbb{Z}_{+}$, we have

$\displaystyle\|H-\mathbb{E}[\tilde{H}_{t+\tau}|\mathcal{F}_{t-1}]\|\leq
CM_{H}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\rho^{\tau}+O_{H}\sqrt{\delta}$
$\displaystyle\|b-\mathbb{E}[b_{t+\tau}|\mathcal{F}_{t-1}]\|\leq
CM_{b}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\rho^{\tau}+O_{b}\sqrt{\delta}.$

We verify this assumption in Lemma 4.13. Although the constant $\delta$ can be
arbitrary, we will set it to the same $\delta$ from Section 4.

To measure progress, we define the $Q$-gap function (where
$\tilde{z}\equiv(\tilde{x},\tilde{y})\in Z:=X\times Y$)

$Q(\tilde{z},z)\equiv Q([\tilde{x};\tilde{y}],[x;y])=[\langle
H\tilde{x}-b,y\rangle-\hat{g}(y)]-[\langle
Hx-b,\tilde{y}\rangle-\hat{g}(\tilde{y})]$

and the corresponding primal-dual gap function

(35) $\displaystyle g(\tilde{z})=\max_{z\in Z}Q(\tilde{z},z).$

In view of this definition, we introduce our last assumption, which is a
regularity condition known as sharpness, as first introduced by Polyak [34].
There exists a $\mu>0$ such that

$\displaystyle\mu\sqrt{V(x,x^{*})}\leq g(z)-g(z^{*}),\hskip 5.0pt\forall
z\equiv(x,y)\in X\times Y.$

###  4.1 A stochastic primal-dual algorithm

We use a primal-dual method, Algorithm 2, to solve the min-max problem of
(33). Besides a dual and primal step, the algorithm separately collects $\tau$
samples to estimate gradients. With some abuse of notation, we refer to the
stochastic estimates for the dual step at iteration $t$ as $\tilde{H}_{t,Y}$
and $\tilde{b}_{t,Y}$, while for the primal step we refer to them as
$\tilde{H}_{t,X}$ and $\tilde{b}_{t,X}$ (see Lines 5 and 7). This method is
similar to the conditional temporal difference methods studied in [20]. The
reason we need separate samples is to reduce the bias in both the primal and
dual steps. We then return a weighted average of all past $k$ points.

Algorithm 2 Conditional stochastic primal-dual (CSPD)

1: Input: $X\subseteq\mathbb{R}^{n}$; $x_{-1}=x_{0}\in X$;
$k,\tau\in\mathbb{Z}_{++}$; $\\{\eta_{t},\lambda_{t},\theta_{t}\\}_{t=1}^{k}$,
Arbitrary $y_{0}\in Y$

2: for $t=1,\ldots,k$ do

3: $z_{t}=x_{t-1}+\theta_{t}(x_{t-1}-x_{t-2})$

4: Collect $\tau$ samples, use $\tau$-th sample to set
$(\tilde{H}_{t,Y},\tilde{b}_{t,Y})=(\tilde{H}_{(2t-1)\tau)},\tilde{b}_{(2t-1)\tau)})$

5: $y_{t}=\mathrm{argmin}_{y\in Y}\
\\{\langle-\tilde{H}_{t,Y}z_{t}+\tilde{b}_{t,Y},y\rangle+\hat{g}(y)+\lambda_{t}V(y_{t-1},y)\\}$

6: Collect $\tau$ more samples, use $\tau$-th sample to set
$\tilde{H}_{t,X}=\tilde{H}_{2t\tau}$

7: $x_{t}=\mathrm{argmin}_{x\in X}\ \\{\langle
y_{t},\tilde{H}_{t,X}x\rangle+\eta_{t}V(x_{t-1},x)\\}$

8: end for

9:
$(\bar{x}_{k},\bar{y}_{k})=\frac{2}{k(k+1)}\sum\limits_{t=1}^{k}t\cdot(x_{t},y_{t})$

10: return $\bar{x}_{k}$

We provide the general convergence result.

######  Proposition 4.1.

Let $\\{\gamma_{t},\eta_{t},\lambda_{t},\theta_{t}\\}_{t}$ be a set of non-
negative reals satisfying

$\displaystyle\gamma_{t-1}\eta_{t-1}\leq\gamma_{t}\eta_{t},\hskip
3.0pt\gamma_{t-1}\lambda_{t-1}\leq\gamma_{t}\lambda_{t},\gamma_{t}\theta_{t}$
$\displaystyle=\gamma_{t-1},$

and there exist some $p,q\in(0,1)$ for all $t$ satisfying
$\|H\|^{2}/(p\lambda_{t})-q\eta_{t}\leq 0$. Then

$\displaystyle\sum\limits_{t=1}^{k}\gamma_{t}Q(z_{t},z)\leq\gamma_{k}\eta_{k}D_{X}^{2}+\gamma_{k}\lambda_{k}D_{Y}^{2}+\sum\limits_{t=1}^{k}\Lambda_{t}(x,y),$

where

$\displaystyle\Lambda_{t}(x,y)$
$\displaystyle:=-\frac{\gamma_{t}\lambda_{t}}{2}(1-p)\|y_{t}-y_{t-1}\|^{2}-\frac{(1-q)\gamma_{t}\eta_{t}}{2}\|x_{t}-x_{t-1}\|^{2}$
$\displaystyle\hskip
10.0pt+\gamma_{t}\langle(\tilde{H}_{t,Y}-H)z_{t}-(\tilde{b}_{t,Y}-b),y_{t}-y\rangle-\gamma_{t}\langle(\tilde{H}_{t,X}-H)^{T}y_{t},x_{t}-x\rangle.$

A proof can be found in [23]. Next, we bound the stochastic error on a convex
set $U$, which generalizes the primal and dual space, $X$ and $Y$. Due to
space constraints, we only show the main steps of the proof and provide
references to auxiliary parts.

######  Proposition 4.2.

Let $\\{\gamma_{t},\alpha_{t},c_{t}\\}$ be non-negative scalars satisfying
$\gamma_{t-1}\alpha_{t-1}\leq\gamma_{t}\alpha_{t}$. Let $\tilde{G}_{t\tau}$ be
an $\mathcal{F}_{t\tau}$-measurable stochastic vector such that

(36)
$\displaystyle\|\tilde{G}_{t\tau}-\mathbb{E}[\tilde{G}_{t\tau}|\mathcal{F}_{(t-1)\tau}]\|\leq\varsigma\text{
w.p. $1-\delta$},$ (37)
$\displaystyle\|G-\mathbb{E}[\tilde{G}_{t\tau}|\mathcal{F}_{(t-1)\tau}]\|_{*}\leq
CM\rho^{\tau}+O\sqrt{\delta},$

where $C,M,\varsigma,\tau\in\mathbb{R}_{+}$, $\rho,\delta\in[0,1)$, and
$G\in\mathbb{R}^{n}$ are constants. Then for any vectors $\\{u_{t},v_{t}\in
U\\}_{t}$ where $u_{t}$ and $v_{t}$ are $\mathcal{F}_{t\tau}$-measurable, with
probability at least $1-(k+1)\delta$,

$\displaystyle\sum\limits_{t=1}^{k}\gamma_{t}\langle
G-\tilde{G}_{t\tau},u_{t}-u\rangle-\frac{c_{t}\gamma_{t}}{2}\|u_{t}-u_{t-1}\|^{2}\leq\gamma_{k}c_{k}D_{U}^{2}$
$\displaystyle+\sum\limits_{t=1}^{k}\gamma_{t}\bigg{(}(CM\rho^{\tau}+O\sqrt{\delta})D_{U}+\frac{4((CM\rho^{\tau})^{2}+O^{2}\delta+\varsigma^{2})}{c_{t}}\bigg{)}+\varsigma
D_{U}\sqrt{8\sum\limits_{t=1}^{k}\gamma_{t}^{2}\ln\frac{1}{\delta}},$

where $D_{U}^{2}=\max_{x,y\in U}V(x,y)$.

######  Proof 4.3.

With $u_{0}^{v}:=u_{0}$ and $\Delta_{t}:=G-\tilde{G}_{t\tau}$, let
$u^{v}_{t}:=\mathrm{argmin}_{u\in
U}\\{\langle\Delta_{t},u\rangle+\alpha_{t}V(u,u^{v}_{t-1})\\}$. Since $G$ is a
constant and $\tilde{G}_{t\tau}$ is assumed to be
$\mathcal{F}_{t\tau}$-measurable, a simple proof by induction ensures
$u_{t}^{v}$ is $\mathcal{F}_{t\tau}$-measurable. Now, we decompose

$\displaystyle\sum\limits_{t=1}^{k}\gamma_{t}\langle\Delta_{t},u_{t}-u\rangle-\frac{c_{t}\gamma_{t}}{2}\|u_{t}-u_{t-1}\|^{2}$
$\displaystyle=\sum\limits_{t=1}^{k}\underbrace{\gamma_{t}\langle\Delta_{t},u_{t}-u_{t-1}\rangle-\frac{c_{t}\gamma_{t}}{2}\|u_{t}-u_{t-1}\|^{2}}_{A_{t}}+\underbrace{\gamma_{t}\langle\Delta_{t},u_{t-1}-u_{t-1}^{v}\rangle}_{B_{t}}+\underbrace{\gamma_{t}\langle\Delta_{t},u_{t-1}^{v}-u\rangle}_{C_{t}}.$

We have the following bound, which holds with probability $1-k\delta$ by union
bound:

$\displaystyle\sum\limits_{t=1}^{k}A_{t}$
$\displaystyle\leq\sum\limits_{t=1}^{k}\frac{\gamma_{t}}{2c_{t}}\|G-\tilde{G}_{t\tau}\|_{*}^{2}$
$\displaystyle\leq\sum\limits_{t=1}^{k}\frac{2\gamma_{t}}{c_{t}}[\|G-\mathbb{E}[\tilde{G}_{t\tau}|\mathcal{F}_{(t-1)\tau}]\|_{*}^{2}+\|\mathbb{E}[\tilde{G}_{t\tau}|\mathcal{F}_{(t-1)\tau}]-\tilde{G}_{t\tau}\|_{*}^{2}]$
$\displaystyle\leq\sum\limits_{t=1}^{k}\frac{2\gamma_{t}}{c_{t}}((CM\rho^{\tau})^{2}+O^{2}\delta+\varsigma^{2}),$

which is due to Young’s inequality, $(a+b)^{2}\leq 2a^{2}+2b^{2}$, and then
(36) and (37). Next, we can bound

$\displaystyle\sum\limits_{t=1}^{k}B_{t}$
$\displaystyle=\sum\limits_{t=1}^{k}(B_{t}-\mathbb{E}[B_{t}|\mathcal{F}_{(t-1)\tau}])+\mathbb{E}[B_{t}|\mathcal{F}_{(t-1)\tau}]$
$\displaystyle\leq\varsigma
D_{U}\sqrt{8\sum\limits_{t=1}^{k}\gamma_{t}^{2}\ln\frac{1}{\delta}}+[CM\rho^{\tau}+O\sqrt{\delta}]D_{U}\sum\limits_{t=1}^{k}\gamma_{t}\text{
(w.p. $1-\delta$)},$

where we used the Azuma-Hoeffding inequality to bound the martingale-
difference sequence $B_{t}-\mathbb{E}[B_{t}|\mathcal{F}_{(t-1)\tau}]$ (and
used the fact $\Delta_{t}$ is $\mathcal{F}_{t\tau}$-measurable while
$u_{t-1}-u_{t-1}^{v}$ are $\mathcal{F}_{(t-1)\tau}$-measurable as well as (36)
and definition of $D_{U}$; see [38, Lemma 7] for full details) and (37) to
bound $\mathbb{E}[B_{t}|\mathcal{F}_{(t-1)\tau}]$. Finally,

$\sum\limits_{t=1}^{k}C_{t}\leq\gamma_{k}c_{k}D_{U}^{2}+\sum\limits_{t=1}^{k}\frac{\gamma_{t}}{c_{t}}\|\Delta_{t}\|_{*}^{2}\leq\gamma_{k}c_{k}D_{U}^{2}+\sum\limits_{t=1}^{k}\frac{2\gamma_{t}}{c_{t}}((CM\rho^{\tau})^{2}+O^{2}\delta+\varsigma^{2}),$

where the first and second inequality can be shown similarly to [23, Lemma
4.10] and the bound on $A_{t}$, respectively. Combining everything finishes
the proof.

We are now ready to complete the proof of convergence. Recall the gap function
$g(\bar{z})\equiv\max_{z\in Z}Q(\bar{z},z)$ and the weighted average solution
$\bar{z}_{k}\equiv(\bar{x}_{k},\bar{y}_{k})$.

######  Theorem 4.4.

Let Section 4 and Section 4 take place. If the hypothesis of Proposition 4.1
holds, then with probability at least $1-2(k+1)\delta$,

$\begin{split}\big{(}&\sum\limits_{t=1}^{k}\gamma_{t}\big{)}g(\bar{z}_{k})\leq
2\gamma_{k}(\eta_{k}D_{X}^{2}+\lambda_{k}D_{Y}^{2})+(D_{X}\varsigma_{X}+D_{Y}\varsigma_{Y})\sqrt{8\sum\limits_{t=1}^{k}\gamma_{t}^{2}\ln\frac{1}{\delta}}\\\
&+\sum\limits_{t=1}^{k}\gamma_{t}\bigg{[}C(M_{X}D_{X}+M_{Y}D_{Y})\rho^{\tau}+(O_{X}D_{X}+O_{Y}D_{Y})\sqrt{\delta}+\frac{4\varsigma_{X}^{2}}{\eta_{t}(1-q)}+\frac{4\varsigma_{Y}^{2}}{\lambda_{t}(1-p)}\\\
&\hskip
40.0pt+4C^{2}\Big{(}\frac{M_{X}^{2}}{\eta_{t}(1-q)}+\frac{M_{Y}^{2}}{\lambda_{t}(1-p)}\Big{)}\rho^{2\tau}+4\Big{(}\frac{O_{X}^{2}}{\eta_{t}(1-q)}+\frac{O_{Y}^{2}}{\lambda_{t}(1-p)}\Big{)}\delta\bigg{]},\end{split}$

where $\overline{\theta}:=\max_{t}\theta_{t}$,

$\displaystyle\Omega_{Y}$ $\displaystyle:=\max\\{\|y\|:y\in Y\\},$
$\displaystyle\Omega_{X}$
$\displaystyle:=\|x^{*}\|+(1+\overline{\theta})\sqrt{2}D_{X},$ $\displaystyle
M_{X}$
$\displaystyle:=M_{H}\Omega_{Y}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta},$
$\displaystyle M_{Y}$
$\displaystyle:=\big{(}M_{H}\Omega_{X}+M_{b}\big{)}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta},$
$\displaystyle\varsigma_{X}$ $\displaystyle:=2M_{X},$
$\displaystyle\varsigma_{Y}$ $\displaystyle:=2M_{Y},$ $\displaystyle O_{X}$
$\displaystyle:=O_{H}\Omega_{Y},$ $\displaystyle O_{Y}$
$\displaystyle:=O_{H}\Omega_{X}+O_{b}.$

Constants $M_{H},M_{b},O_{H},O_{b}$ are from Section 4 and Section 4.

######  Proof 4.5.

We first bound the error $\sum\limits_{t=1}^{k}\Lambda_{t}(x,y)$ defined in
Proposition 4.1 as

(38)
$\begin{split}&\sum\limits_{t=1}^{k}-\frac{\gamma_{t}\lambda_{t}(1-p)}{2}\|y_{t}-y_{t-1}\|^{2}+\gamma_{t}\langle(\underbrace{\tilde{H}_{t,Y}z_{t}-\tilde{b}_{t,Y}}_{\tilde{G}_{t,Y}^{\tau}})-(\underbrace{Hz_{t}-b}_{G_{Y}}),y_{t}-y\rangle\\\
&\hskip
25.0pt-\frac{(1-q)\gamma_{t}\eta_{t}}{2}\|x_{t}-x_{t-1}\|^{2}-\gamma_{t}\langle\underbrace{\tilde{H}_{t,X}^{T}y_{t}}_{\tilde{G}_{t,X}^{\tau}}-\underbrace{H^{T}y_{t}}_{G_{X}},x_{t}-x\rangle.\end{split}$

We wish to use Proposition 4.2, which requires verification of (36) and (37).
First, because matrix $\tilde{H}_{t,Y}$ and vector $\tilde{b}_{t,Y}$ are
generated using the $((2t-1)\tau)$-th sample in Algorithm 2 and
$z_{t}=x_{t-1}+\theta_{t}(x_{t-1}-x_{t-2})$, then $\tilde{G}^{\tau}_{t,Y}$ is
$\mathcal{F}_{(2t-1)\tau}$-measurable. Second, using Section 4, we have with
probability at least $1-\delta$,

$\displaystyle\|\tilde{G}_{t,Y}^{\tau}-\mathbb{E}[\tilde{G}_{t,Y}|\mathcal{F}_{(2t-2)\tau}]\|_{*}$
$\displaystyle\leq[\|\tilde{H}_{t,Y}\|_{*}+\|\mathbb{E}[\tilde{H}_{t,Y}|\mathcal{F}_{(2t-2)\tau}]\|](\|x^{*}\|+\|x_{t-1}-x^{*}\|+\theta_{t}\|x_{t-1}-x_{t-2}\|)$
$\displaystyle\hskip
10.0pt+\|\tilde{b}_{t,Y}\|_{*}+\|\mathbb{E}[\tilde{b}_{t,Y}|\mathcal{F}_{(2t-2)\tau}]\|_{*}$
$\displaystyle\leq
2M_{H}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\|(\|x^{*}\|+(1+\theta_{t})\sqrt{2}D_{X})+2M_{b}\Big{(}\ln\frac{e}{\delta}\Big{)}^{\beta}\equiv\varsigma_{Y}.$

This verifies (36) up to an appropriate shift in the index $t\tau$, which is
without loss of generality. One can similarly show
$\|{G}_{Y}-\mathbb{E}[\tilde{G}_{t,Y}^{\tau}|\mathcal{F}_{(2t-2)(\tau+1)}]\|_{*}\leq
CM_{Y}\rho^{\tau}+O_{Y}\sqrt{\delta}$ using Section 4, verifying (37). Fixing
$U=Y$, $u_{t}=y_{t}$, $u=y$, and $c_{t}=\lambda_{t}$, Proposition 4.2 bounds
the summation over the first and second summands from (38) by (with
probability $1-(k+1)\delta$)

$\displaystyle\sum\limits_{t=1}^{k}$
$\displaystyle\gamma_{t}\langle\tilde{G}_{t,Y}^{\tau}-G_{Y},y_{t}-y\rangle-\frac{\gamma_{t}\lambda_{t}}{2}(1-p)\|y_{t}-y_{t-1}\|^{2}\leq\varsigma_{Y}D_{Y}\sqrt{8\sum\limits_{t=1}^{k}\gamma_{t}^{2}\ln\frac{1}{\delta}}$
$\displaystyle+\gamma_{k}\lambda_{k}D_{Y}^{2}+\sum\limits_{t=1}^{k}\gamma_{t}\bigg{(}(CM_{Y}\rho^{\tau}+O_{Y}\sqrt{\delta})D_{Y}+\frac{4((CM_{Y}\rho^{\tau})^{2}+O_{Y}^{2}\delta+\varsigma_{Y}^{2})}{\lambda_{t}(1-p)}\bigg{)}.$

A nearly identical argument can be made for the third and fourth terms from
(38). Combining everything back derived so far into (38) yields a bound on
$\sum\limits_{t=1}^{k}\Lambda_{t}(x,y)$.

Now, by definition of the primal-dual gap and convexity of $Q$,
$(\sum\limits_{t=1}^{k}\gamma_{t})\cdot g(\bar{z}_{k})=\max_{z\in
Z}(\sum\limits_{t=1}^{k}\gamma_{t})Q(\bar{z}_{k},z)\leq\max_{z\in
Z}\sum\limits_{t=1}^{k}\gamma_{t}Q(z_{t},z)$. We finish by plugging our bound
on $\sum\limits_{t=1}^{k}\Lambda_{t}(x,y)$ into Proposition 4.1 and dividing
through by the sum over $\gamma_{t}$.

We now fix a step size. The proof is skipped since it is similar to [23,
Corollary 4.3].

######  Corollary 4.6.

Let
$\eta_{t}=\frac{3\sqrt{2}\|H\|D_{Y}t+3\varsigma_{X}t^{3/2}}{2\sqrt{2}D_{X}t}$,
$\lambda_{t}=\frac{3\sqrt{2}\|H\|D_{X}t+3\varsigma_{Y}t^{3/2}}{2\sqrt{2}D_{Y}t}$,
and $\theta_{t}=\frac{t-1}{t}$, where the variance and gradient bound terms
$\varsigma_{I}$, $M_{I}$, and $O_{I}$ for $I\in\\{X,Y\\}$ are defined in
Theorem 4.4. Then with probability at least $1-2(k+1)\delta$,

$\displaystyle g(\bar{z}_{k})$
$\displaystyle\leq\frac{12\|H\|D_{X}D_{Y}}{k+1}+\frac{(48+3\sqrt{2}+\frac{4\sqrt{2}}{\sqrt{3}}\sqrt{\ln
1/\delta})(D_{X}\varsigma_{X}+D_{Y}\varsigma_{Y})}{\sqrt{k}}$
$\displaystyle\hskip
10.0pt+\frac{48}{\sqrt{k}}\Big{[}C^{2}\big{(}\frac{D_{X}M_{X}^{2}}{\varsigma_{X}}+\frac{D_{Y}M_{Y}^{2}}{\varsigma_{Y}}\big{)}\rho^{2\tau}+\big{(}\frac{D_{X}O_{X}^{2}}{\varsigma_{X}}+\frac{D_{Y}O_{Y}^{2}}{\varsigma_{Y}}\big{)}\delta\Big{]}$
$\displaystyle\hskip
10.0pt+C\big{(}D_{X}M_{X}+D_{Y}M_{Y}\big{)}\rho^{\tau}+(D_{X}O_{X}+D_{Y}O_{Y})\sqrt{\delta}.$

Recall that $\varsigma_{X}$, $\varsigma_{Y}$, $M_{X}$, and $M_{Y}$, defined in
Theorem 4.4, are related to $\ln(e/\delta)$. Thus, $\delta$ must be
judiciously chosen to minimize the overall error. In view of Section 4 and
assuming $g(z^{*})=0$, the above result implies we can find a $\bar{z}_{k}$
where $V(\bar{z}_{k},z^{*})\leq\epsilon$ with high probability in
$O(\frac{\|H\|^{2}}{\mu^{2}k^{2}}+\frac{\mathrm{var}}{\mu^{2}k})$ iterations,
where “$\mathrm{var}$” bounds statistical errors. However, this is suboptimal
when $\mathrm{var}=0$, i.e., the deterministic case.

To achieve the optimal deterministic error rate, we introduce a shrinking
multi-epoch algorithm in Algorithm 3. It repeatedly calls Algorithm 2 and uses
the solution to the previous call to warm-start the next call. We name each
call an epoch. Each epoch $s$ will use an updated feasible set for the primal
variable $X_{s}$ that is shrunk from the previous epoch. The feasible set for
the dual variable $Y$ remains unchanged. A key feature is that $x^{*}\in
X_{s}$ for every epoch. Since $X_{s}$ gets smaller, we also update the
diameter $D_{X_{s}}$. The choice of step sizes still adheres to those
prescribed in Corollary 4.6 up to the difference in the diameter $D_{X_{s}}$,
while the number of iterations $k_{s}$ increases each epoch. The mixing time
and failure rate, respectively $\tau_{s}$ and $\delta_{s}$, remain constant.

Algorithm 3 Shrinking multi-epoch CSPD

Input: $X\subseteq\mathbb{R}^{n}$; $p_{0}\in X$; $D_{0}\text{ s.t.
}V(p_{0},x^{*})\leq D_{0}^{2}$; $k,\tau,S\in\mathbb{Z}_{++}$

for $s=1,\ldots,S$ do

$D_{s}^{2}:=2^{-(s-1)}D_{0}^{2}$

$X_{s}:=\\{x\in X:V(x,p_{s-1})\leq D_{s}^{2}\\}$

$p_{s}\leftarrow$ Algorithm 2 with $X_{s}$, $p_{s-1}$, $k_{s}$, $\tau_{s}$,
$\delta_{s}$, $\\{\eta_{t},\lambda_{t},\theta_{t}\\}_{t}$ specified by
Corollaries 4.6 and 4.7

end for

return $p_{S}$

We now establish the main convergence result. Recall the terms $M_{I}$,
$\varsigma_{I}$, and $O_{I}$ for $I\in\\{X,Y\\}$ are defined as in Theorem
4.4. In the following result, we assume
$D_{0},M_{X},M_{Y},\varsigma_{X},\varsigma_{Y}\geq 1$ and
$\varsigma_{X}\geq\|H\|D_{Y}$. We skip the proof, which follows similarly to
[23, Proposition 4.5] (as well as applying the previous result and Section 4).

######  Proposition 4.7.

Suppose Section 4 holds. Let
$\\{\eta_{t},\lambda_{t},\theta_{t}\\}_{t=1}^{k_{s}}$ be defined as in
Corollary 4.6 and the number of iterations, mixing time, and failure rate
satisfy

$\displaystyle k_{s}$ $\displaystyle\geq
5000\max\Big{\\{}\frac{\|H\|D_{Y}}{\mu},\frac{225+\ln(1/\delta)}{\mu^{2}}\big{(}\varsigma_{X}^{2}+\frac{D_{Y}^{2}\varsigma_{Y}^{2}}{4D_{0}^{2}}\cdot
2^{s-1}\big{)}\Big{\\}}$ $\displaystyle\tau_{s}$
$\displaystyle\equiv\tau\geq\frac{\ln\big{(}10C(D_{0}M_{X}+D_{Y}M_{Y})\max\\{\frac{1}{\mu\epsilon},\frac{1}{\sqrt{\epsilon}}\\}\big{)}}{\ln(1/\rho)}$
$\displaystyle\delta_{s}$
$\displaystyle\equiv\delta\leq\frac{\epsilon}{200}\min\big{\\{}\frac{\mu^{2}}{D_{0}^{2}O_{X}^{2}+D_{Y}^{2}O_{Y}^{2}},\frac{1}{D_{0}O_{X}^{2}+D_{Y}O_{Y}^{2}}\big{\\}}.$

Choosing the number of epochs as
$S=\big{\lceil}\log_{2}\frac{2D_{0}^{2}}{\epsilon}\big{\rceil}$, then
$\|p_{S}-x^{*}\|^{2}\leq\epsilon$ with probability at least
$1-2(N_{S}+S)\delta$, and the total iterations is at most $N_{S}$, where

$\displaystyle
N_{S}:=O(1)\max\Big{\\{}\frac{\|H\|D_{Y}}{\mu}\ln\Big{(}\frac{D_{0}^{2}}{\epsilon}\Big{)},\frac{\ln(e/\delta)}{\mu^{2}}\max\\{\varsigma_{X}^{2}\ln\Big{(}\frac{D_{0}}{\epsilon}\Big{)},\frac{D_{Y}^{2}\varsigma_{Y}^{2}}{\epsilon}\\}\Big{\\}},$

while the sample complexity is $2N_{S}\tau$. Here, $O(1)$ is some absolute
constant.

We have achieved the optimal deterministic and stochastic error rate,
improving upon Corollary 4.6. Next, we apply this generic primal-dual method
to the setting of estimating the natural gradient matrix $E_{K_{t}}$.

###  4.2 Sampling complexity of policy optimization and evaluation

We use the setup of $X=\mathbb{R}^{n\cdot k}$, $h=0$, and
$Y=\\{y:\|y\|_{2}\leq 1\\}$ in (33). Under the assumption
$H=\mathbb{E}_{\Pi_{K}}[\tilde{H}]\succ 0$, this ensures the primal-dual gap
function from (33) is $g([x,y])=\|Hx-b\|_{2}$ for any feasible primal-dual
solution $[x,y]$. We also choose the strongly convex function
$\omega=\frac{1}{2}\|\cdot\|_{2}^{2}$ so that the Bregman divergence is
$V(x,x^{\prime})=\frac{1}{2}\|x^{\prime}-x\|_{2}^{2}$.

For the remainder of this section, we will denote $x_{t}$ as the state at time
$t$ of the linear dynamical system (1) while $u_{t}=-Kx_{t}+v_{t}$ is the
feedback, where $v_{t}$ is the policy disturbance introduced in (3). We let
the stochastic process be
$\\{\xi_{t}:\xi_{t}=[x_{t},u_{t},x_{t+1},u_{t+1}]\\}_{t}$, i.e., the current
and next state-action pair. It remains to show Section 4, Section 4, and
Section 4.

First, we show the gradient is bounded with high probability. The proof, which
is skipped due to space constraints, follows from the Hanson-Wright inequality
(e.g., [3, Example 2.12]) and a change of variables.

######  Lemma 4.8.

Let $\ell\sim\mathcal{N}(\mu,\Sigma)$. For any $\delta\in(0,1)$, we have

$\displaystyle\mathbb{P}\bigg{\\{}\|\ell\|^{2}>\frac{(4\|\Sigma\|\sqrt{\ln(1/\delta)}+\|\Sigma\|_{F})^{2}-\|\Sigma\|_{F}^{2}}{8\|\Sigma\|}+\mathrm{Tr}(\Sigma)+2\|\mu\|^{2}\bigg{\\}}\leq\delta.$

Define the following covariance matrices

(39)
$\displaystyle\widetilde{\Sigma}_{K}^{(t)}=\begin{bmatrix}\Sigma_{K}^{(t)}&-\Sigma_{K}^{(t)}K^{T}\\\
-K\Sigma_{K}^{(t)}&K\Sigma_{K}^{(t)}K^{T}+\sigma^{2}I_{k}\end{bmatrix},\Sigma_{K}^{(t)}=\sum\limits_{p=0}^{t-1}(A-BK)^{p}\Psi[(A-BK)^{T}]^{p},$

where $I_{k}$ is the identity matrix of size $k$. The matrix
$\widetilde{\Sigma}_{K}^{(t)}$ is the covariance matrix of the state-action
pair $\begin{bmatrix}x_{t}&u_{t}\end{bmatrix}$ with $u_{t}=-Kx_{t}$
conditioned on $x_{0}$, and similarly, $\Sigma_{K}^{(t)}$ is the covariance
matrix of the state $x_{t}$. Here, $\Psi$ is the covariance matrix for the
uncertainty in the dynamics of (1).

In view of the lemma, we construct the event for the $t$-th state-action pair,

(40)
$\displaystyle\mathcal{E}_{t}=\Big{\\{}\|x_{t}\|^{2}+\|u_{t}\|^{2}\leq\big{(}4\|\widetilde{\Sigma}^{(t)}\|+\mathrm{Tr}(\widetilde{\Sigma}^{(t)})+2(\|x_{0}\|^{2}+\|u_{0}\|^{2})\big{)}\sqrt{\ln\frac{e}{\delta}}\Big{\\}},$

This event occurs with probability (w.p.) at least $1-\delta$. Since we
collect every $(\tau\cdot t)$ and $(\tau\cdot t+1)$-th state-action pair
(i.e., $\tau$-th sample) in Algorithm 2 (see Lines 5 and 7), we define the
intersection of the following events:

$\displaystyle\mathcal{E}=\bigcap_{1\leq t\leq
T}(\mathcal{E}_{t\cdot(\tau+1)-1}\cap\mathcal{E}_{t\cdot(\tau+1)}),$

where $T\geq 1$. By union bound, $\mathbb{P}(\mathcal{E})\geq 1-2T\delta$.

Now we are ready to show Section 4. Recall the feedback lies in
$u_{t}\in\mathbb{R}^{k}$ and that $\sigma>0$ is the added policy disturbance
level (c.f. (3)) that is a free parameter.

######  Lemma 4.9.

Let $K$ be the current controller, assumed to be stable, and the events
$\mathcal{E}_{t}$ and $\mathcal{E}_{t+1}$ take place. Then Section 4 is
satisfied with $\beta=1$ and

$\displaystyle M_{H}$
$\displaystyle:=17(1+\|K\|)^{4}\max\Big{\\{}1,\Big{(}\frac{5J(K)}{\sigma_{\min}(Q)}+2\|x_{0}\|^{2}+(4+k)\sigma^{2}\Big{)}^{2}\Big{\\}},$
$\displaystyle M_{b}$ $\displaystyle:=M_{H}\max\\{\|Q\|,\|R\|\\}.$

A proof can be found in Section A.1. To bound the variance-like term in
Section 4, we now bound the norm of the expected matrix and vector,
$H=\mathbb{E}[\tilde{H}]$ and $b=\mathbb{E}[\tilde{b}]$. First, we will need a
couple of technical lemmas.

######  Lemma 4.10 (Section 7.1 P5 [14]).

For symmetric positive semidefinite matrix $S$, $\|S\|_{F}\leq\mathrm{Tr}(S)$.

The next lemma is due to [42, Eqn 5.4] and Lemma 4.10.

######  Lemma 4.11.

Let $K$ be the current stable controller and define the radius
$R_{*}=(\|Q\|_{F}+\|R\|_{F})+(\|A\|_{F}^{2}+\|B\|_{F}^{2})\cdot\frac{J(K)}{\sigma_{\min}(\Psi)}$.
Then the optimal solution to (32) satisfies $\|x^{*}\|\leq R_{*}$.

We can now bound $\|H\|$, which follows similarly to [42, Lemma B.2] as well
as $\|b\|$, which follows from the equality $b=Hx^{*}$ and the previous lemma.

######  Lemma 4.12.

Matrix $H=\mathbb{E}_{\Pi_{K}}[\tilde{H}]$ and vector
$b=\mathbb{E}_{\Pi_{K}}[\tilde{b}]$ have norms

$\displaystyle\|H\|$ $\displaystyle\leq
L_{H}:=1+4(1+\|K\|^{2}_{F})\frac{J(K)}{\sigma_{\min}(Q)}+\sigma^{2}(k+2)\text{
and }\|b\|\leq L_{H}R_{*}.$

By triangle inequality, $\|\tilde{H}_{t}-H\|\leq\|\tilde{H}_{t}\|+\|H\|$ and
$\|\tilde{b}_{t}-b\|\leq\|\tilde{b}_{t}\|+\|b\|$. This combined with the last
three lemmas completes verification of Section 4.

Next, we supply constants to bound the bias of the gradient estimation error.

######  Lemma 4.13.

Suppose the event $\mathcal{E}$ occurs with some $\delta\in(0,1]$. Then
Section 4 is satisfied with $M_{H}$ and $M_{b}$ from Lemma 4.9, $\beta=5/4$,
$\rho=\rho(A-BK)$,

$\displaystyle C$
$\displaystyle:=\frac{1}{2}\Big{(}M_{H}+\sqrt{\frac{n}{1-\rho}}\Big{)}$
$\displaystyle O_{H}$
$\displaystyle:=65(n+k)\max\\{\|\widetilde{\Sigma}_{K}\|^{2},\|\widetilde{\Sigma}_{K}\|^{4}\\}$
$\displaystyle O_{b}$
$\displaystyle:=41\max\\{\|Q\|,\|R\|\\}(n+k)^{2}\max\\{\|\widetilde{\Sigma}_{K}\|^{2},\|\widetilde{\Sigma}_{K}\|^{4}\\}.$

The proof is in Section A.2. This verifies Section 4.

Our final assumption is Section 4, which addresses the sharpness condition.
When solving $Hx=b$ where $H$ is positive definite, the sharpness constant is
the smallest singular value since $\|Hx-b\|\geq\sigma_{\min}(H)\|x-x^{*}\|$.
Here, we see the reason for having $\sigma$, the added disturbance level from
(3), be positive.

######  Lemma 4.14.

Let $K$ be the stable policy and denote $H=\mathbb{E}_{\Pi_{K}}[\tilde{H}]$
with the induced spectrum of $\rho\equiv\rho(A-BK)$. Then

$\displaystyle\sigma_{\min}(H)$
$\displaystyle\geq\frac{1}{2}\min\Big{\\{}1,\big{(}1+(1+\|K\|)^{2}\frac{J(K)}{\sigma_{\min}(Q)}+\sigma^{2}\cdot
k\big{)}^{-1}(1-\rho^{2})$ $\displaystyle\hskip
10.0pt\cdot\min\\{\max(0,\sigma^{2}-\sigma_{\min}(\Psi)\|K\|^{2}),\sigma_{\min}(\Psi)\\}\Big{\\}}.$

The proof is given in Section A.3. Thus, we have verified all assumptions.

The following corollary is a result of combining Theorem 3.9 and Proposition
4.7 as well as the constants provided in the results above. We assume the
initial diameter is $D_{0}=O(R_{*})$, where $R_{*}$ is from Lemma 4.11, and
let $\rho_{T}:=\max_{0\leq t\leq T}\rho(A-BK_{t})$.

######  Corollary 4.15.

Consider using the parameters from Theorem 3.9 to run Algorithm 1 and
Proposition 4.7 to run Algorithm 3, where the initial solution $x_{0}$
satisfies $\|x_{0}\|=O(R_{*})$, and the error tolerance and failure rate are

$\displaystyle\epsilon$
$\displaystyle=O(1)\min\\{1,\varepsilon/J(K_{0})^{2},\varepsilon/J(K_{0})^{3},\varepsilon^{2/3}/J(K_{0})^{8/3}\\}$
$\displaystyle\delta$
$\displaystyle=O(1)\min\big{\\{}\frac{\varepsilon}{n^{4}J(K_{0})^{13}},\frac{\delta^{*}}{J(K_{0})^{6}\ln(1/\varepsilon)(\ln
1/\delta^{*})^{7/2}}\big{\\}}$

for some $\delta^{*}\in(0,1)$. If $\delta\leq O(1)J(K_{*})^{9}$, then it
suffices to run Algorithm 3 for

$\displaystyle N_{\text{in}}$
$\displaystyle:=O(1)\Big{(}\frac{J(K_{0})^{9}}{\varepsilon(1-\rho_{T}^{2})^{2}}\big{(}\ln\frac{1}{\varepsilon\delta^{*}}\big{)}^{7/2}\ln\frac{1}{\varepsilon}\Big{)}$

iterations to estimate a $E_{K_{t}}$ for the hypothesis of Theorem 3.9 to hold
w.p. at least $1-\delta^{*}O(1)(\ln 1/\varepsilon)^{-1}$. Putting everything
together, under the assumption $K_{0}$ is stable and $\|K_{t}\|$ is bounded by
a constant for all $t$, the sampling complexity to ensure
$J(K_{T})-J(K^{*})\leq\varepsilon$ w.p. at least $1-\delta^{*}$ is
$O(1)N_{\text{in}}\ln(1/\varepsilon)/\ln(1/\rho_{T})$, or equivalently,

$\displaystyle
O\Big{\\{}\frac{J(K_{0})^{9}}{\varepsilon(1-\rho^{2}_{T})^{3}}\Big{(}\ln\frac{1}{\varepsilon\delta^{*}}\Big{)}^{7/2}\Big{(}\ln\frac{1}{\varepsilon}\Big{)}^{2}\Big{\\}}.$

The upper bound on $\delta$ is mild since we can choose a smaller
$\varepsilon$ and/or $\delta^{*}$, while the assumption $\|K_{t}\|$ is bounded
seems to be common in the literature [42, 45], and it can be relaxed at the
expense of having the sampling complexity depend on $\max_{t}\|K_{t}\|$. The
initial objective value $J(K_{0})$ is explicitly shown since it is related to
the mixing rate $\rho_{T}=\max_{0\leq t\leq T}\rho(A-BK_{t})$. In particular,
recall Corollary 3.11 upper bounds $\rho_{T}$ by a constant strictly less than
one that depends only on problem-dependent constants and $J(K_{0})$.
Therefore, we do not need the assumption that $\rho_{T}$ is strictly less than
one almost surely like in [43, 45]. While our work studies the average reward
setting, our work should be amenable to the infinite-horizon discounted
setting as well, since they both share some properties such as the PŁ-
condition and the natural gradient being partially defined by a solution to a
linear system of equations.

##  5 Numerical Experiments

We consider a simple scalar linear dynamical system ($n=k=1$) where the
transition matrices $A$ and $B$ are identity matrices and the cost matrices
$Q$ and $R$ are the identity matrix scaled by 100. We set the system noise
level to $\Psi=0.1^{2}$ and controller noise to $\sigma=0.1$. Figure 1 shows
the absolute error when using natural policy gradient (Algorithm 1) with a
step size of $\eta=0.01$ for 50 iterations. For the policy evaluation, we ran
the shrinking multi-epoch primal-dual algorithm (Algorithm 3) using 3 epochs.
We tuned the number of iterations, samples, and step sizes. In Figure 1, we
see the algorithm on average improves the objective function every iteration.
There are some iterations where the objective function increases, and we
observed this is due to instances where the injected noise was large. The
error decreases steadily at a near linear rate, as suggested by Theorem 4.4.

Figure 1: Error vs. sampling complexity of NPG. The shaded regions is the 10th
to 90th percentile over $20$ random seeds, and the solid blue line is the
median performance.

##  6 Conclusion

Using a variational inequality formulation for the actor step and conditional
stochastic primal-dual method in the critic step to estimate the natural
gradient $E_{K_{t}}$, we obtain a sampling complexity of
$\tilde{O}(1/\varepsilon)$ to ensure $J(K_{t})-J(K^{*})\leq\varepsilon$ with
high probability. Our analysis exploits the fact that the accuracy of the
overall algorithm scales with the variance rather than the standard deviation
for estimating $E_{K_{t}}$. This motivates us to formulate the critic step as
a non-smooth problem to obtain the optimal rate of convergence.

## Acknowledgments

We thank Sihan Zeng for pointing us to [43, 45] and discussing the results
therein. Disclaimer: This report was prepared as an account of work sponsored
by an agency of the United States Government. Neither the United States
Government nor any agency thereof, nor any of their employees, makes any
warranty, express or implied, or assumes any legal liability or responsibility
for the accuracy, completeness, or usefulness of any information, apparatus,
product, or process disclosed, or represents that its use would not infringe
privately owned rights. Reference herein to any specific commercial product,
process, or service by trade name, trademark, manufacturer, or otherwise does
not necessarily constitute or imply its endorsement, recommendation, or
favoring by the United States Government or any agency thereof. The views and
opinions of authors expressed herein do not necessarily state or reflect those
of the United States Government or any agency thereof.

## References

  * [1] B. D. Anderson and J. B. Moore, Optimal control: linear quadratic methods, Courier Corporation, 2007.
  * [2] M. Athans and P. L. Falb, Optimal control: an introduction to the theory and its applications, Courier Corporation, 2013.
  * [3] S. Boucheron, G. Lugosi, and P. Massart, Concentration inequalities: A nonasymptotic theory of independence, Oxford university press, 2013.
  * [4] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear matrix inequalities in system and control theory, SIAM, 1994.
  * [5] S. J. Bradtke, Incremental dynamic programming for on-line adaptive optimal control, PhD thesis, Citeseer, 1994.
  * [6] P. E. Caines, Linear Stochastic Systems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2018.
  * [7] S. Dean, H. Mania, N. Matni, B. Recht, and S. Tu, On the sample complexity of the linear quadratic regulator, Foundations of Computational Mathematics, 20 (2020), pp. 633–679.
  * [8] M. Deisenroth and C. E. Rasmussen, Pilco: A model-based and data-efficient approach to policy search, in Proceedings of the 28th International Conference on machine learning (ICML-11), 2011, pp. 465–472.
  * [9] F. Facchinei and J.-S. Pang, Finite-dimensional variational inequalities and complementarity problems, Springer, 2003.
  * [10] I. Fatkhullin and B. Polyak, Optimizing static linear feedback: Gradient method, SIAM Journal on Control and Optimization, 59 (2021), pp. 3887–3911.
  * [11] M. Fazel, R. Ge, S. Kakade, and M. Mesbahi, Global convergence of policy gradient methods for the linear quadratic regulator, in International Conference on Machine Learning, PMLR, 2018, pp. 1467–1476.
  * [12] C.-N. Fiechter, Pac adaptive control of linear systems, in Proceedings of the tenth annual conference on Computational learning theory, 1997, pp. 72–80.
  * [13] G. Hewer, An iterative technique for the computation of the steady state gains for the discrete optimal regulator, IEEE Transactions on Automatic Control, 16 (1971), pp. 382–384.
  * [14] R. A. Horn and C. R. Johnson, Matrix analysis, Cambridge university press, 2012.
  * [15] L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika, 12 (1918), pp. 134–139.
  * [16] S. M. Kakade, A natural policy gradient, in Advances in Neural Information Processing Systems, T. Dietterich, S. Becker, and Z. Ghahramani, eds., vol. 14, MIT Press, 2001.
  * [17] R. E. Kalman, On the general theory of control systems, in Proceedings First International Conference on Automatic Control, Moscow, USSR, 1960, pp. 481–492.
  * [18] R. E. Kalman et al., Contributions to the theory of optimal control, Bol. soc. mat. mexicana, 5 (1960), pp. 102–119.
  * [19] V. Konda and J. Tsitsiklis, Actor-critic algorithms, Advances in neural information processing systems, 12 (1999).
  * [20] G. Kotsalis, G. Lan, and T. Li, Simple and optimal methods for stochastic variational inequalities, ii: Markovian noise and policy evaluation in reinforcement learning, SIAM Journal on Optimization, 32 (2022), pp. 1120–1155.
  * [21] K. Krauth, S. Tu, and B. Recht, Finite-time analysis of approximate policy iteration for the linear quadratic regulator, Advances in Neural Information Processing Systems, 32 (2019).
  * [22] P. R. Kumar and P. Varaiya, Stochastic Systems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015.
  * [23] G. Lan, First-order and stochastic optimization methods for machine learning, Springer, 2020.
  * [24] G. Lan, Policy mirror descent for reinforcement learning: Linear convergence, new sampling complexity, and generalized problem classes, Mathematical programming, (2022), pp. 1–48.
  * [25] T. Li, F. Wu, and G. Lan, Stochastic first-order methods for average-reward markov decision processes, arXiv preprint arXiv:2205.05800, (2022).
  * [26] Y. Li, T. Zhao, and G. Lan, Homotopic policy mirror descent: Policy convergence, implicit regularization, and improved sample complexity, arXiv preprint arXiv:2201.09457, (2022).
  * [27] L. Ljung, System identification, in Signal analysis and prediction, Springer, 1998, pp. 163–173.
  * [28] S. Lojasiewicz, Une propriété topologique des sous-ensembles analytiques réels, Les équations aux dérivées partielles, 117 (1963), pp. 87–89.
  * [29] D. Malik, A. Pananjady, K. Bhatia, K. Khamaru, P. Bartlett, and M. Wainwright, Derivative-free methods for policy optimization: Guarantees for linear quadratic systems, in The 22nd international conference on artificial intelligence and statistics, PMLR, 2019, pp. 2916–2925.
  * [30] H. Mania, S. Tu, and B. Recht, Certainty equivalence is efficient for linear quadratic control, Advances in Neural Information Processing Systems, 32 (2019).
  * [31] H. Mohammadi, A. Zare, M. Soltanolkotabi, and M. R. Jovanović, Convergence and sample complexity of gradient methods for the model-free linear–quadratic regulator problem, IEEE Transactions on Automatic Control, 67 (2021), pp. 2435–2450.
  * [32] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro, Robust stochastic approximation approach to stochastic programming, SIAM Journal on optimization, 19 (2009), pp. 1574–1609.
  * [33] B. T. Polyak, Gradient methods for the minimisation of functionals, USSR Computational Mathematics and Mathematical Physics, 3 (1963), pp. 864–878.
  * [34] B. T. Polyak, Sharp minima. institute of control sciences lecture notes, moscow, ussr, 1979, in IIASA workshop on generalized Lagrangians and their applications, IIASA, Laxenburg, Austria, 1979.
  * [35] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, Proximal policy optimization algorithms, arXiv preprint arXiv:1707.06347, (2017).
  * [36] R. S. Sutton and A. G. Barto, Reinforcement learning: An introduction, MIT press, 2018.
  * [37] S. Tu and B. Recht, Least-squares temporal difference learning for the linear quadratic regulator, in International Conference on Machine Learning, PMLR, 2018, pp. 5005–5014.
  * [38] Y. Wang, W. Chen, Y. Liu, Z.-M. Ma, and T.-Y. Liu, Finite sample analysis of the gtd policy evaluation algorithms in markov setting, Advances in Neural Information Processing Systems, 30 (2017).
  * [39] D. Wilson, R. Sharp, and S. Hassan, The application of linear optimal control theory to the design of active automotive suspensions, Vehicle System Dynamics, 15 (1986), pp. 105–118.
  * [40] C. S. Withers, The moments of the multivariate normal, Bulletin of the Australian Mathematical Society, 32 (1985), pp. 103–107.
  * [41] W.-Y. Yan, K. L. Teo, and J. B. Moore, A gradient flow approach to computing lq optimal output feedback gains, Optimal Control Applications and Methods, 15 (1994), pp. 67–75.
  * [42] Z. Yang, Y. Chen, M. Hong, and Z. Wang, Provably global convergence of actor-critic: A case for linear quadratic regulator with ergodic cost, Advances in neural information processing systems, 32 (2019).
  * [43] S. Zeng, T. T. Doan, and J. Romberg, A two-time-scale stochastic optimization framework with applications in control and reinforcement learning, arXiv preprint arXiv:2109.14756, (2021).
  * [44] K. Zhou, J. Doyle, and K. Glover, Robust and optimal control, Prentice hall, 1996.
  * [45] M. Zhou and J. Lu, Single time-scale actor-critic method to solve the linear quadratic regulator with convergence guarantees, arXiv preprint arXiv:2202.00048, (2022).

##  Appendix A Proofs from Section 4

###  A.1 Proof of Lemma 4.9

######  Proof A.1 (Proof of Lemma 4.9).

Recall from (20) that
$\phi_{t}\equiv\phi(x_{t},-Kx_{t})=\mathrm{svec}(z_{t}z_{t}^{T})$ and
$z_{t}=\left[\begin{array}[]{cc}x_{t}&-Kx_{t}\end{array}\right]$. We have

$\displaystyle\|\tilde{H}_{t}\|=\bigg{\|}\begin{bmatrix}1&0^{T}\\\
\phi_{t}&\phi_{t}[\phi_{t}-\phi_{t+1}]^{T}\end{bmatrix}\bigg{\|}$
$\displaystyle\leq 2+2\|z_{t}\|^{4}+\|z_{t}\|^{2}\|z_{t+1}\|^{2}$
$\displaystyle\leq 2+4\max\\{\|z_{t}\|^{4},\|z_{t+1}\|^{4}\\}.$

By definition of $\mathcal{E}_{t}$, the covariance matrix
$\widetilde{\Sigma}_{K}^{(t)}$ from (39), then for any $t$,

(41)
$\begin{split}\|z_{t}\|_{2}^{2}&\leq\big{[}4\|\widetilde{\Sigma}_{K}^{(t)}\|+\mathrm{Tr}(\widetilde{\Sigma}_{K}^{(t)})+2(\|x_{0}\|^{2}+\|u_{0}\|^{2})\big{]}\sqrt{\ln(e/\delta)}\\\
&\leq
4\big{[}((1+\|K\|)^{2}\|\Sigma_{K}^{(t)}\|+\sigma^{2})+((1+\|K\|^{2})\mathrm{Tr}(\Sigma_{K}^{(t)})+\sigma^{2}k)\\\
&\hskip 10.0pt+2(\|x_{0}\|^{2}+\|u_{0}\|^{2})\big{]}\sqrt{\ln(e/\delta)}\\\
&\leq\Big{[}5(1+\|K\|)^{2}\big{(}\frac{5J(K)}{\sigma_{\min}(Q)}+2\|x_{0}\|^{2}+(4+k)\sigma^{2}\big{)}\Big{]}\sqrt{\ln(e/\delta)},\end{split}$

where the last line used
$\|\Sigma_{K}^{(t)}\|\leq\mathrm{Tr}(\Sigma_{K}^{(t)})\leq\mathrm{Tr}(\Sigma_{K})\leq
J(K)/\sigma_{\min}(Q)$ thanks in part to Lemma 3.1. Combining the last two
inequalities gets us the upper bound on $\|\tilde{H}_{t}\|$. One can similarly
bound $\|\tilde{b}_{t}\|$ by recalling
$c(x_{t},u_{t})=x_{t}^{T}Qx_{t}+u_{t}^{T}Ru_{t}$.

###  A.2 Proof of Lemma 4.13

We will first need the geometrically fast mixing property of a linear
dynamical system. In the context of our system (4), we have $F=A-BK$,
$\Lambda=\Psi+\sigma^{2}I$ and $\Gamma=1$ in the following lemma.

######  Proposition A.2 (Proposition 3.1 [37]).

Take the linear dynamical system $X_{t+1}=FX_{t}+w_{t}$, where
$\omega_{t}\sim\mathcal{N}({0},\Lambda)$ and $\Lambda\succ 0$. Suppose
$\left\|F^{k}\right\|\leq\Gamma\rho^{k}$ for all $k\geq 0$, where $\Gamma>0$
and $\rho\in(0,1)$. Let $\mathbb{P}_{X_{k}}\left(\cdot\mid X_{0}=x\right)$
denote the conditional distribution of $X_{k}$ given $X_{0}=x$. Then for all
$k\geq 0$ and any distribution $\nu_{0}$ over $x$,

$\displaystyle\mathbb{E}_{x\sim\nu_{0}}\left[\left\|\mathbb{P}_{X_{k}}\left(\cdot\mid
X_{0}=x\right)-\nu_{\infty}\right\|_{\mathrm{tv}}\right]\leq\frac{\Gamma}{2}\sqrt{\mathbb{E}_{\nu_{0}}\left[\|x\|^{2}\right]+\frac{n}{1-\rho^{2}}}\rho^{k}.$

We need to bound higher moments of a mean-zero Gaussian random vector. The
result is a corollary of Isserlis’ Theorem [15].

######  Lemma A.3.

Let $X\sim\mathcal{N}(0,\Sigma)$, where $\Sigma\in\mathbb{R}^{n\times n}$. Let
$I=\\{i_{j}\in[n]\\}_{j=1}^{2p}$ be a set of (possibly non-distinct) indices
for some positive integer $p$. Then

$\displaystyle\mathbb{E}[\prod_{j=1}^{2p}X_{i_{j}}]\leq\big{(}3\cdot
5\cdot\ldots\cdot(2p-1)\big{)}\|\Sigma\|^{p}.$

######  Proof A.4.

Since $|I|=2p$ is even, Isserlis’ Theorem [15, 40] derives us the closed-form
expression
$\mathbb{E}[\prod_{j=1}^{2p}X_{i_{j}}]=\sum\limits_{\sigma\in\Pi(I)}\prod_{(i,j)\in\sigma}\mathrm{Cov}(X_{i}X_{j})\leq|\Pi(I)|\cdot\|\Sigma\|^{p}$,
where $\Pi(I)$ denotes the set of all partitions of $I$ into (disjoint) pairs.
Here, it is known the cardinality of $\Pi(I)$ is
$|\Pi(I)|=(2p)!/(2^{p}p!)=1\cdot 3\cdot\ldots\cdot(2p-1)$.

######  Proof A.5 (Proof of Lemma 4.13).

Let $z_{t}=\begin{bmatrix}x_{t}&u_{t}\end{bmatrix}$ be the state-action pair
and $\mathrm{d}P_{t^{\prime}|t}$ be the probability density function for
$z_{t^{\prime}}$ conditioned on $z_{t}$. Recall
$H=\mathbb{E}_{\xi\sim\Pi_{K}}[\tilde{H}(\xi)]$, and denote $\tilde{H}_{1}$ as
the stochastic estimate of $\tilde{H}$ without the constant value “1” in the
leftmost column of the first row. Denoting $\mathcal{E}^{c}_{t+\tau}$ as the
complement of $\mathcal{E}_{t+\tau}$, then

$\displaystyle\|H-\mathbb{E}[\tilde{H}_{t+\tau}|\mathcal{F}_{t-1},\mathcal{E}_{t+\tau}]\|_{*}$
$\displaystyle=\Big{\|}\int_{\xi\in\mathcal{E}_{t+\tau}}\tilde{H}_{1}(\xi)[\mathrm{d}\Pi_{K}(\xi)-\mathrm{d}P_{t+\tau|t}(\xi)]+\int_{\xi\not\in\mathcal{E}_{t+\tau}}\tilde{H}_{1}(\xi)\mathrm{d}\Pi_{K}(\xi)\Big{\|}_{*}$
$\displaystyle\leq\int_{\xi\in\mathcal{E}_{t+\tau}}\|\tilde{H}(\xi)\|_{*}|\mathrm{d}\Pi_{K}(\xi)-\mathrm{d}P_{t+\tau|t}(\xi)|+\int_{\xi\not\in\mathcal{E}_{t+\tau}}\|\tilde{H}_{1}(\xi)\|_{*}\mathrm{d}\Pi_{K}(\xi)$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}\Big{(}\frac{M_{H}}{2}\sqrt{\|z_{t}\|^{2}+\frac{n}{1-\rho^{2}}}\Big{)}\ln\frac{e}{\delta}\Big{)}\rho^{\tau}+\sqrt{\mathbb{E}_{\Pi_{K}}\|\tilde{H}_{1}\|^{2}}\sqrt{\mathrm{Pr}\\{\mathcal{E}_{t+\tau}^{c}\\}}$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}\frac{M_{H}}{2}\Big{(}M_{H}+\sqrt{\frac{n}{1-\rho}}\Big{)}\Big{(}\ln\frac{e}{\delta}\Big{)}^{5/4}\rho^{\tau}+\sqrt{\mathbb{E}_{\Pi_{K}}\|\tilde{H}_{1}\|^{2}}\sqrt{\mathrm{Pr}\\{\mathcal{E}_{t+\tau}^{c}\\}}$

where the bound on the first term in (i) is by Proposition A.2 and Lemma 4.9
and the second term by Cauchy-Schwarz. Inequality (ii) is by
$\|z_{t}\|^{2}\leq M_{H}\sqrt{\ln(e/\delta)}$, where we used (41) and
$M_{H}\geq 1$. By definition of $\tilde{H}(\xi)$ in (22) and denoting
$\mathbb{E}=\mathbb{E}_{\Pi_{K}}$,

$\displaystyle\mathbb{E}\|\tilde{H}_{1}\|^{2}\leq\mathbb{E}\|\tilde{H}_{1}\|^{2}_{F}$
$\displaystyle\leq
10(n+k)^{2}\max_{i,j,u,v}\Big{\\{}\mathbb{E}([z_{t}]_{i}[z_{t}]_{j})^{2},\mathbb{E}([z_{t}]_{i}[z_{t}]_{j}[z_{t}]_{u}[z_{t}]_{v})^{2},\mathbb{E}([z_{t}]_{i}[z_{t}]_{j}[z_{t+1}]_{i}[z_{t+1}]_{v})^{2}\Big{\\}}$
$\displaystyle\leq
40(n+k)^{2}\max_{i}\\{\mathbb{E}([z_{t}]_{i})^{4},\mathbb{E}([z_{t}]_{i})^{8},\mathbb{E}([z_{t+1}]_{i})^{8}\\}$
$\displaystyle\leq\big{(}65(n+k)\max\\{\|\widetilde{\Sigma}_{K}\|^{2},\|\widetilde{\Sigma}_{K}\|^{4}\\}\big{)}^{2}\equiv
O_{H}^{2},$

where the last line is by Lemma A.3. Combining the last two results yields the
bound on the bias for $\tilde{H}$. By noting the similarity in definition
between $\tilde{H}$ and $\tilde{b}$ in (22), one can similarly bound
$\|b-\mathbb{E}[b_{t+\tau}|\mathcal{F}_{t-1}]\|_{*}$.

###  A.3 Proof of Lemma 4.14

######  Proof A.6 (Proof of Lemma 4.14).

In [42, Lemma B.2], it was shown that

$\displaystyle\|H^{-1}\|$ $\displaystyle\leq
1+(1+\|\widetilde{\Sigma}_{K}\|_{F})\|\Xi_{K}^{-1}\|$ $\displaystyle\leq
1+(1+\|\widetilde{\Sigma}_{K}\|_{F})(1-\rho^{2})^{-1}\max\Big{\\{}\frac{1}{\max\\{0,\sigma^{2}-\sigma_{\min}(\Psi)\|K\|^{2}\\}},\frac{1}{\sigma_{\min}(\Psi)}\Big{\\}}.$

Here, $\widetilde{\Sigma}_{K}=\lim_{t\to\infty}\widetilde{\Sigma}_{K}^{(t)}$
is the steady state covariance matrix for the state-action pair
$\begin{bmatrix}x&u\end{bmatrix}$, defined according to (39). We can further
refine this bound by observing $\widetilde{\Sigma}_{K}$ is symmetric and
positive definite, hence

$\displaystyle\|\widetilde{\Sigma}_{K}\|_{F}$
$\displaystyle\leq(1+\|K\|)^{2}\|\Sigma_{K}\|_{F}+\sigma^{2}\cdot
k\leq(1+\|K\|)^{2}\frac{J(K)}{\sigma_{\min}(Q)}+\sigma^{2}\cdot k,$

where the second inequality is by combining Lemma 3.1 and Lemma 4.10.
Combining the two bounds above, we conclude

$\displaystyle\sigma_{\min}(H)$
$\displaystyle\geq\frac{1}{2}\min\Big{\\{}1,\big{(}1+(1+\|K\|)^{2}\frac{J(K)}{\sigma_{\min}(Q)}+\sigma^{2}\cdot
k\big{)}^{-1}(1-\rho^{2})$ $\displaystyle\hskip
10.0pt\cdot\min\\{\max(0,\sigma^{2}-\sigma_{\min}(\Psi)\|K\|^{2}),\sigma_{\min}(\Psi)\\}\Big{\\}}.$