# Global Convergence of Receding-Horizon Policy Search
in Learning Estimator Designs

Xiangyuan Zhang Saviz Mowlavi Mouhacine Benosman Tamer Başar X. Zhang and T.
Başar are with the Department of ECE and CSL, University of Illinois
Urbana–Champaign, Urbana, IL 61801 USA (e-mails: {xz7, basar1}@illinois.edu).
S. Mowlavi and M. Benosman are with Mitsubishi Electric Research Laboratories,
Cambridge, MA 02139 USA (e-mails: {mowlavi, benosman}@merl.com).A preliminary
version of this manuscript [1] appeared in the proceedings of the 2023
American Control Conference and was presented on June 2, 2023, in San Diego,
CA.

###### Abstract

We introduce the receding-horizon policy gradient (RHPG) algorithm, the first
PG algorithm with provable global convergence in learning the optimal linear
estimator designs, i.e., the Kalman filter (KF). Notably, the RHPG algorithm
does not require any prior knowledge of the system for initialization and does
not require the target system to be open-loop stable. The key of RHPG is that
we integrate vanilla PG (or any other policy search directions) into a dynamic
programming outer loop, which iteratively decomposes the infinite-horizon KF
problem that is constrained and non-convex in the policy parameter into a
sequence of static estimation problems that are unconstrained and strongly-
convex, thus enabling global convergence. We further provide fine-grained
analyses of the optimization landscape under RHPG and detail the convergence
and sample complexity guarantees of the algorithm. This work serves as an
initial attempt to develop reinforcement learning algorithms specifically for
control applications with performance guarantees by utilizing classic control
theory in both algorithmic design and theoretical analyses. Lastly, we
validate our theories by deploying the RHPG algorithm to learn the Kalman
filter design of a large-scale convection-diffusion model. We open-source the
code repository at https://github.com/xiangyuan-zhang/LearningKF.

##  1 Introduction

In recent years, policy-based reinforcement learning (RL) methods [2, 3, 4, 5]
have gained increasing attention in continuous control applications [6, 7, 8].
While traditional model-based techniques synthesize controller designs in a
case-by-case manner [9, 10], model-free policy gradient (PG) methods promise a
universal framework that learns controller designs in an end-to-end fashion.
The universality of model-free PG methods makes them desired candidates in
complex control applications that involve nonlinear system dynamics and
imperfect state measurements. Despite countless empirical successes, the
theoretical properties of model-free PG methods still need to be thoroughly
investigated in continuous control. Initiated by [11], a recent line of
research has well-analyzed the sample complexity of zeroth-order PG methods in
several linear state-feedback control benchmarks, including linear-quadratic
regulator (LQR) [11, 12, 13, 14, 15, 16, 17], distributed/decentralized LQR
[18, 19], and linear robust control [20, 21, 22]. However, the theoretical
properties of PG methods remain elusive in the output-feedback control
settings, where the state measurement process could be corrupted by
statistical noises and/or other (possibly adversarial) disturbances.

In this work, we study the convergence and sample complexity of PG methods in
the discrete-time infinite-horizon Kalman filtering (KF) problem [23, 9].
Recognized as one of the cornerstones of modern control theory [24], the KF
problem aims to generate optimal estimates of the unknown system states over
time by utilizing a sequence of observed measurements corrupted by statistical
noises. Furthermore, in the linear-quadratic Gaussian (LQG) problem, the
separation principle [25] states that the optimal control law combines KF and
LQR. Thus, KF is a fundamental benchmark for studying the sample complexity of
model-free PG methods beyond state-feedback settings.

Despite being the dual problem to noise-less LQR [26, 25], the KF problem
possesses a substantially more complicated optimization landscape from the
model-free PG perspective since the KF itself is a dynamical system rather
than a static matrix. Specifically, the optimization problem over dynamic
filters might admit multiple suboptimal stationary points, and the optimal KF
possesses a set of equivalent realizations up to similarity transformations
[27, 28]. None of the above challenges appear when using model-free PG to
learn a static LQR policy [11, 12, 13, 14, 15, 16, 17]. As a result of the
challenging landscape the filtering problem presents, only a few papers have
focused on dynamic output-feedback settings. In particular, [27] has analyzed
the optimization landscape of LQG, and [28] has shown that an informativity-
regularized PG method provably converges to an optimal filter in the
continuous-time KF problem, assuming that the model is known. However, [28]
has assumed that the target system is open-loop stable and assumed that the
control engineer has prior knowledge of a filter that satisfies an
informativity condition. It is also unclear if the techniques in [28] can be
directly applied to the model-free setting and result in any sample complexity
guarantees. Thus, obtaining sample complexity of model-free PG methods in the
KF problem has remained a significant challenge.

This work addresses these challenges by introducing a receding-horizon PG
(RHPG) algorithm and establishing its global convergence and sample
complexity. In contrast to direct policy search, the RHPG algorithm integrates
vanilla PG (or any other policy search directions) into a dynamic programming
(DP) outer loop, which iteratively decomposes the infinite-horizon KF problem
that is constrained and non-convex in the policy parameter into a sequence of
static estimation problems that are unconstrained and strongly-convex. Then,
we show that solving the sequence of static estimation problems results in the
global convergence of RHPG toward the KF, which is the optimal linear filter.
We further establish the total sample complexity of the RHPG to be
$\widetilde{\mathcal{O}}(\epsilon^{-2})$ for the learned filter to be
$\epsilon$-close in policy distance to KF, which is the first sample
complexity result of PG methods in the output-feedback control settings.
Notably, the RHPG algorithm does not require any prior knowledge of the system
to generate a valid initialization and does not require the target system to
be open-loop stable. This removes two restrictive assumptions in the previous
work [28]. We validate our theories by learning the KF design of a large-scale
convection-diffusion model.

Compared to the preliminary results included in [1], this work presents a
comprehensive study of RHPG as a model-free RL approach in learning estimator
designs. In particular, our contributions, in addition to those listed in [1],
are three-folded. First, we analyze the optimization landscape in Theorem 4.1,
which clarifies the properties of quadratic programs in RHPG and provides a
theoretical foundation for selecting the algorithmic parameters of PG methods.
Second, we discuss the insights of the RHPG design in Sec. 3.2 and compare
RHPG with standard PG methods regarding parametrization and landscape,
computational efficiencies, and requirements on the simulation oracles. These
discussions provide the intuitions behind the mathematical developments of
RHPG and should benefit future research on improving the algorithm design and
the theoretical analysis. Lastly, we open-source numerical experiments on the
data-driven estimation of a large-scale dynamical system, which corroborates
the theories and demonstrates both the effectiveness and scalability of RHPG.

This work attempts to develop RL algorithms specifically for control and
estimation tasks with performance guarantees, by utilizing classic control
theory in both algorithmic design and theoretical analyses. The dual theories
and implementations of RHPG to the LQR problem have been presented in [17].
Through this line of work, we demonstrate the significant utilization of DP in
overcoming the challenging optimization landscape and streamlining the
analyses when deploying model-free PG methods to linear control and estimation
tasks. Due to the separation principle [25], our results shed light on
applying model-free PG methods in solving the LQG problem through a sequential
design of controller and estimator.

The structure of the paper is as follows. In Sec. 2, we define the infinite-
and finite-horizon settings of the KF problem and formulate them as policy
optimization problems. In Sec. 3, we introduce the RHPG algorithm and provide
the general theory and intuition that backs its design. In Sec. 4, we analyze
the optimization landscape of solving the KF problem using RHPG and establish
the global convergence and sample complexity guarantees of the algorithm.
Lastly, we present the numerical studies on a convection-diffusion model in
Sec. 5. The paper ends with the concluding remarks of Sec. 6, and an appendix
that contains formal proofs of the main results.

###  1.1 Notations

For a square matrix $X$, we denote its trace, spectral norm, condition number,
and spectral radius by $\text{Tr}(X)$, $\|X\|$, $\kappa_{X}$, and $\rho(X)$
resp. We define the $W$-induced norm of $X$ as $\|X\|^{2}_{W}:=\max_{z\neq
0}\frac{z^{\top}X^{\top}WXz}{z^{\top}Wz}$. If $X$ is further symmetric, we use
$X>0$, $X\geq 0$, $X\leq 0$, and $X<0$ to denote that $X$ is positive definite
(pd), positive semi-definite (psd), negative semi-definite (nsd), and negative
definite (nd), resp. We use $x\sim\mathcal{N}(\mu,\Sigma)$ to denote a
Gaussian random vector with mean $\mu$ and covariance $\Sigma$. Lastly, we use
$\bm{I}$ and $\bm{0}$ to denote the identity and zero matrices, resp., with
appropriate dimensions.

##  2 Preliminaries

###  2.1 Infinite-Horizon Kalman Filtering

Consider the discrete-time linear time-invariant system

$\displaystyle x_{t+1}=Ax_{t}+w_{t},\quad y_{t}=Cx_{t}+v_{t},$ (2.1)

where $x_{t}\in\mathbb{R}^{n}$ is the state, $y_{t}\in\mathbb{R}^{m}$ is the
output measurement, and $w_{t}\sim\mathcal{N}(\bm{0},W)$,
$v_{t}\sim\mathcal{N}(\bm{0},V)$ are sequences of i.i.d. zero-mean Gaussian
noises for some $W,V>0$, also independent of each other. The initial state is
also assumed to be a Gaussian random vector such that
$x_{0}\sim\mathcal{N}(\bar{x}_{0},X_{0})$, independent of $\\{w_{t},v_{t}\\}$,
with $\bar{x}_{0}\neq\bm{0}$ and $X_{0}>0$. Additionally, we assume that
$(C,A)$ is observable and note that the condition $W>0$ readily leads to
controllability of $(A,W^{1/2})$, which is a standard condition in KF.

The KF problem aims to generate a sequence of estimated states, denoted by
$\hat{x}_{t}$ for each $t$, that minimizes the infinite-horizon mean-square
error (MSE):

$\displaystyle\mathcal{J}_{\infty}:=\lim_{N\to\infty}\frac{1}{N}\mathbb{E}\bigg{\\{}\sum_{t=0}^{N}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\bigg{\\}}.$
(2.2)

Moreover, each $\hat{x}_{t}$ can only depend on the history and output
measurements up to but not including $t$, i.e., $\\{y_{0},\cdots,y_{t-1}\\}$.
The celebrated result of Kalman [23] showed that the
$\mathcal{J}_{\infty}$-minimizing filter (could also be called $1$-step
predictor), which exists under the controllability and the observability
conditions, has the form of

$\displaystyle\hat{x}_{t+1}^{*}$
$\displaystyle=(A-L^{*}C)\hat{x}_{t}^{*}+L^{*}y_{t},\quad\hat{x}_{0}^{*}=\bar{x}_{0},$
(2.3) $\displaystyle L^{*}$
$\displaystyle=A\Sigma^{*}C^{\top}(V+C\Sigma^{*}C^{\top})^{-1},$ (2.4)

where $L^{*}$ is the Kalman gain and $\Sigma^{*}$ represents the unique pd
solution to the filter algebraic Riccati equation (FARE):

$\displaystyle\Sigma$ $\displaystyle=A\Sigma A^{\top}-A\Sigma
C^{\top}(V+C\Sigma C^{\top})^{-1}C\Sigma A^{\top}+W.$ (2.5)

Hence, without any loss of optimality, we can restrict the search to the class
of filters of the form $\hat{x}_{t+1}=A_{L}\hat{x}_{t}+B_{L}y_{t}$ and then
parametrize the KF problem as a minimization problem over $A_{L}$ and $B_{L}$
subject to a stability constraint111Extending the results in this work to the
setting with instantaneous feedback measurement (i.e., allowing $\hat{x}_{t}$
to depend also on $y_{t}$, and hence replacing $y_{t}$ in (2.6) with
$y_{t+1}$) would be straightforward.

$\displaystyle\min_{A_{L},B_{L}}\quad\mathcal{J}_{\infty}(A_{L},B_{L})\quad\text{s.t.}\quad\hat{x}_{t+1}=A_{L}\hat{x}_{t}+B_{L}y_{t}~{}\
\text{and}~{}\ \rho(A_{L})<1.$ (2.6)

Note that by (2.3), there indeed exists a solution to (2.6) where
$(A_{L}^{*},B_{L}^{*})=(A-L^{*}C,L^{*})$. Note also that when the pair $(A,C)$
is known, (2.6) involves an over-parametrization since solving (2.6) is
equivalent to optimizing a single variable $B_{L}$. However, in the model-free
setting where $(A,C)$ is unknown, which is the target setting of our paper, it
is reasonable to parametrize the KF problem as in (2.6). Until now, obtaining
sample complexity of model-free PG methods in solving the KF problem (2.6) has
remained a major challenge.

Figure 1: Algorithm 1, executed forward in time, constructs an $(h+1)$-horizon
KF problem from $t=0$ to $t=h+1$ at each iteration indexed by $h$.

###  2.2 Finite-Horizon Kalman Filtering

We now discuss the finite-$N$-horizon KF problem, also described by the system
dynamics (2.1). Adopting the same parametrization as in (2.6), but this time
allowing time-dependence, and again without any loss of optimality, we
represent the finite-horizon KF problem as a minimization problem over a
sequence of time-varying filter parameters $\\{A_{L_{t}},B_{L_{t}}\\}$, for
all $t\in\\{0,\cdots,N-1\\}$,

$\displaystyle\min_{\\{A_{L_{t}},B_{L_{t}}\\}}$
$\displaystyle\mathcal{J}\big{(}\\{A_{L_{t}},B_{L_{t}}\\}\big{)}:=\mathbb{E}\bigg{\\{}\sum_{t=0}^{N}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\bigg{\\}}$
(2.7)
$\displaystyle\text{s.t.}\quad\hat{x}_{t+1}=A_{L_{t}}\hat{x}_{t}+B_{L_{t}}y_{t},\quad\hat{x}_{0}=\bar{x}_{0}.$

The minimum in (2.7) can be achieved by
$(A^{*}_{L_{t}},B^{*}_{L_{t}})=(A-L^{*}_{t}C,L^{*}_{t})$, where $L^{*}_{t}$ is
the time-varying Kalman gain

$\displaystyle
L_{t}^{*}=A\Sigma_{t}^{*}C^{\top}(V+C\Sigma_{t}^{*}C^{\top})^{-1},\quad\Sigma_{0}^{*}=X_{0}.$
(2.8) $\displaystyle\Sigma^{*}_{t+1}$
$\displaystyle=A\Sigma^{*}_{t}A^{\top}-A\Sigma^{*}_{t}C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}C\Sigma^{*}_{t}A^{\top}+W.$
(2.9)

The solutions $\Sigma^{*}_{t}$, for all $t\in\\{0,\cdots,N-1\\}$, generated by
the filter Riccati difference equation (FRDE) (2.9) always exist and are
unique and pd, due to $V>0$, $W>0$, and the iteration starts with
$\Sigma^{*}_{0}=X_{0}>0$.

##  3 Receding-Horizon Policy Gradient

###  3.1 Kalman Filtering and Dynamic Programming

It is well known that the solution of the FRDE (2.9) converges monotonically
to the stabilizing solution of the FARE (2.5) at an exponential rate [29, 30].
Then, it readily follows that the optimal time-varying filter
$(A_{L_{t}}^{*},B_{L_{t}}^{*})$ to the finite-horizon KF problem (2.7) also
converge monotonically to the time-invariant $(A_{L}^{*},B_{L}^{*})$ as
$N\to\infty$. We present this convergence result in the following theorem,
which plays an important role in our algorithm design.

###### Theorem 3.1

The finite-horizon Kalman gain as in (2.8) converges to the infinite-horizon
Kalman gain defined in (2.4) exponentially fast as $N\to\infty$. Specifically,
using $\|\cdot\|_{*}$ to denote the $\Sigma^{*}$-induced norm and letting

$\displaystyle
N_{0}=\frac{1}{2}\cdot\frac{\log\big{(}\frac{\|X_{0}-\Sigma^{*}\|_{*}\cdot\kappa_{\Sigma^{*}}\cdot\|A_{L}^{*}\|\cdot\|C\|}{\epsilon\cdot\lambda_{\min}(V)}\big{)}}{\log\big{(}\frac{1}{\|A_{L}^{*}\|_{*}}\big{)}}+1.$
(3.1)

where $\|A_{L}^{*}\|_{*}<1$, we have that, for all $N\geq N_{0}$, it holds
that $\|L^{*}_{N-1}-L^{*}\|\leq\epsilon$ for any $\epsilon>0$. If additionally
$X_{0}>\Sigma$ holds, then $L^{*}_{t}$ is stabilizing for all $t\geq 0$ in the
sense that $\rho(A_{L^{*}_{t}})=\rho(A-L^{*}_{t}C)<1$.

The proof of Theorem 3.1 is provided in §A. Theorem 3.1 quantifies how
different system parameters affect the non-asymptotic convergence rate of the
time-varying filters to the time-invariant KF. It further demonstrates that if
$N\sim\mathcal{O}(\log(\epsilon^{-1}))$, then the filter
$(A_{L_{N-1}}^{*},B_{L_{N-1}})$ will be $\epsilon$-close to the infinite-
horizon KF ($A_{L}^{*},B_{L}^{*}$). Furthermore, if $\epsilon$ is sufficiently
small (i.e., smaller than the stability margin of $A_{L}^{*}$), then it holds
that $\rho(A_{L_{N-1}}^{*})<1$. When $X_{0}>\Sigma$, the system is
sufficiently excited, and as a result, the frozen filter at any $t\geq 0$ is
stable in the sense that $\rho(A_{L^{*}_{t}})<1$.

###  3.2 Algorithm Design

Instead of solving the infinite-horizon KF problem (2.7) directly, we
introduce the RHPG algorithm, which first selects a sufficiently large problem
horizon $N$ according to Theorem 3.1, then constructs and solves $N$ static
estimation problems sequentially (see Figure 1 for an illustration) using PG
methods. We describe the procedure of the RHPG algorithm below.

Algorithm 1 Receding-Horizon Policy Gradient (RHPG)

0: Problem horizon $N$

1: Initialize $A_{L_{0}},B_{L_{0}}\leftarrow\bm{0}_{n\times n},\bm{0}_{n\times
m}$

2: for $h\in\\{0,\cdots,N-1\\}$ do

3: Solve (3.2) using PG methods until convergence

4: Use the convergent filter $A_{L_{h}},B_{L_{h}}$ to warm-start PG updates
for the next iteration

5: end for

6: return $A_{L_{N-1}},B_{L_{N-1}}$

We provide detailed implementations of zeroth-order and first-order RHPG in
Algorithms 2 and 3, respectively.

The RHPG algorithm is executed forward in time, and in the first iteration,
the algorithm learns the optimal filter for a one-step static estimation
problem. Then, every subsequent RHPG iteration extends the problem horizon by
adding one additional time step after the initial one. Formally, at each
iteration indexed by $h$, the RHPG algorithm constructs an $(h+1)$-horizon KF
problem from $t=0$ to $t=h+1$, but we fix the filter parameters for all
$t\in\\{0,\cdots,h-1\\}$ as those generated from earlier iterations and only
optimize for latest filter parameters $(A_{L_{h}},B_{L_{h}})$. This renders
each iteration of the RHPG algorithm into solving a static estimation problem
that is quadratic in $(A_{L_{h}},B_{L_{h}})$.

Mathematically, for every $h\in\\{0,\cdots,N-1\\}$, the RHPG algorithm solves
the following minimization problems

$\displaystyle\min_{A_{L_{h}},B_{L_{h}}}\mathcal{J}_{h}:=\mathbb{E}_{x_{0},w_{t},v_{t},\theta_{0}}\Big{\\{}\sum_{t=0}^{h+1}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\Big{\\}}$
(3.2) s.t.
$\displaystyle\hat{x}_{t+1}=A_{L_{t}}^{*}\hat{x}_{t}+B_{L_{t}}^{*}y_{t},\
\forall t\in\\{0,\cdots,h-2\\},\ \hat{x}_{0}=\bar{x}_{0}$ $\displaystyle
x_{t+1}=Ax_{t}+w_{t},\ \forall t\in\\{0,\cdots,h-2\\},\
x_{0}\sim\mathcal{N}(\bar{x}_{0},X_{0})$
$\displaystyle\hat{x}_{h}=A_{L_{h-1}}^{*}\hat{x}_{h-1}+B_{L_{h-1}}^{*}y_{h-1}+\theta_{0}$
(3.3) $\displaystyle x_{h}=Ax_{h-1}+w_{h-1}+\theta_{0}.$ (3.4)

where $\theta_{0}\sim\mathcal{N}(\bm{0},\Theta)\in\mathbb{R}^{n}$ is sampled
independently to $x_{0},w_{t},v_{t}$ and satisfies $\Theta>0$. The purpose of
injecting an additional “small” noise $\theta_{0}$ in (3.3)-(3.4) is to ensure
the strict convexity of the quadratic program (3.2); we will formally justify
it in Sec. 4.

Figure 2: Illustration of the different optimization landscapes and search
directions under vanilla PG and RHPG. *Note that this graph does not imply
that the cost values of $\mathcal{J}_{h}$ are higher than that of
$\mathcal{J}_{\infty}$.

Before presenting the theoretical analyses, we provide a few comments and
discussions regarding the algorithmic designs of RHPG. From the optimization
landscape perspective and compared to the vanilla PG update that is trapped by
suboptimal first-order stationary points, RHPG provides a new policy search
direction that points toward the global optimum by learning the KF design
step-by-step, where each subproblem admits a unique global optimum due to the
quadratic landscape. See Figure 2 for an illustration. Our algorithmic design
shares a similar flavor as the curriculum learning literature [31, 32, 33],
where an agent (controller/estimator) evolves by mastering the simplest tasks
first and then gradually conquering tasks that are more and more challenging.
In our case, the RHPG algorithm first learns the filter capable of predicting
only the immediate next state from scratch but then keeps adapting/evolving to
handle new filtering tasks with longer and longer problem horizons (cf., the
warm-start step in Line $4$ of Algorithm 1). When the problem horizon becomes
sufficiently large (as characterized by Theorem 3.1), the filter converges
globally, and the behaviors of the infinite-horizon KF, such as closed-loop
stability, begin to emerge. Lastly, since RHPG starts by learning the simplest
static estimation task and every subproblem is unconstrained, RHPG does not
require any specific filter for initializations. In other words, it suffices
to initialize arbitrarily for searching a static estimator.

On the computational side, it may seem at first glance that RHPG is less
efficient compared to vanilla PG since it solves $N$ optimization problems
instead of $1$. This, however, turns out not to be true. When applying
(sampled-based) vanilla PG to the infinite-horizon objective
$\mathcal{J}_{\infty}$ directly, the rollout length is typically a very large
finite number so that an accurate PG estimate can be obtained. The rollout
length for RHPG is $1$ in its first iteration, which coincides with the
iteration that needs the largest number of PG steps since it learns an optimal
one-step static estimator _from scratch_. The rollout length remains very
small in the first few iterations, where the time-varying filters (2.8) are
distinct across time. When the rollout length becomes moderate, only a few PG
updates are needed to fine-tune the filter; see Figure 2. Due to the much
shorter rollout trajectories, the computational efficiency of RHPG is
comparable, if not better, to vanilla PG. We provide the complexity analysis
in Sec. 4.

Lastly, we discuss the requirements for the simulation oracle. To sample the
gradients, we require the standard assumption that the user has access to a
simulator such that for any input filter $(A_{L_{h}},B_{L_{h}})$, the
simulator can return an empirical value of the objective function (3.2). This
requires the simulator to generate exact state trajectories of the simulated
model, but it only reveals a noisy scalar objective value to the learning
algorithm. The requirement is reasonable for the offline learning setting
since the algorithm does not use any system information directly. However,
building a simulator naturally requires knowledge of the system model, which
could be exact, approximate, or simplified. Transferring the simulated
policies to a real system (a.k.a., Sim2Real) that might exhibit different
dynamics requires a _provable_ robustness guarantee of the learned
controller/estimator, which is beyond the scope of the present paper, but it
is an important current and future research topic [34, 21, 35, 22, 36, 37].

###  3.3 Bias of Model-Free Receding-Horizon Filtering

Figure 3: As illustrated, Theorem 3.1 bounds the distance between
$L^{*}_{N-1}$ and $L^{*}$ and Theorem 3.2 analyzes the forward propagation of
the optimization errors from each iteration of RHPG. Combining two theorems
upper-bounds the total policy gap to the infinite-horizon KF.

The RHPG algorithm is backed by Bellman’s principle of optimality, which
requires solving each iteration exactly. However, iterative algorithms such as
PG can only return an $\epsilon$-accurate solution in a finite time. We
generalize the dynamic programming principle in the following theorem to
analyze how computational errors accumulate in the forward DP process and
provide an optimality guarantee for the filter that the RHPG algorithm
returns.

###### Theorem 3.2

Choose the problem horizon $N$ following Theorem 3.1 and assume that one can
compute, for all $h\in\\{0,\cdots,N-1\\}$ and some $\epsilon>0$, filter
$(\widetilde{A}_{L_{h}},\widetilde{B}_{L_{h}})$ that satisfies

$\displaystyle\big{\|}\widetilde{A}_{L_{h}}-\widetilde{A}_{L_{h}}^{*}\big{\|},\big{\|}\widetilde{B}_{L_{h}}-\widetilde{B}_{L_{h}}^{*}\big{\|}\sim\mathcal{O}(\epsilon\cdot\texttt{poly}(\text{system
parameters})),$

where $(\widetilde{A}_{L_{h}}^{*},\widetilde{B}_{L_{h}}^{*})$ is the unique
minimizing solution of $\mathcal{J}_{h}$ in (3.2), after setting the filters
for all $t\in\\{0,\cdots,h-2\\}$ to be those computed in the previous
iterations and are $\epsilon$-close to the minimum of
$\mathcal{J}_{0},\cdots,\mathcal{J}_{h-1}$, respectively. Then, the RHPG
algorithm outputs $(\widetilde{A}_{L_{N-1}},\widetilde{B}_{L_{N-1}})$ that
satisfies $\big{\|}[\widetilde{A}_{L_{N-1}}\
\widetilde{B}_{L_{N-1}}]-[A^{*}_{L}\ B^{*}_{L}]\big{\|}\leq\epsilon$, where
$(A^{*}_{L},B^{*}_{L})$ represents the infinite-horizon KF. If further
$\epsilon$ is sufficiently small such that $\epsilon<1-\|A_{L}^{*}\|_{*}$,
then $\widetilde{A}_{L_{N-1}}$ satisfies $\rho(\widetilde{A}_{L_{N-1}})<1$.

We illustrate Theorem 3.2 in Figure 3 and defer its proof to §B. Theorem 3.2
guarantees that if every iteration of the DP is solved to an
$\mathcal{O}(\epsilon)$-accuracy, then the convergent filter after completing
the $N$-step DP procedure is at most $\epsilon$-away from the exact infinite-
horizon KF. We note that the RHPG algorithm utilizes two layers of
approximations. First, the solution to the infinite-horizon KF problem (2.6)
is approximated by the solution to a finite-horizon KF problem (2.7), where we
choose $N\sim\mathcal{O}(\log(\epsilon^{-1}))$ due to the exponential
attraction of the Riccati equation. Then, we solve the finite-horizon KF
problem (2.7) by integrating forward DP with model-free policy search.
Combining the two steps addresses the infinite-horizon KF task with a provable
global convergence guarantee using only samples of system trajectories.

##  4 Optimization Landscape, Convergence, and Sample Complexity

We first present the optimization landscape of the static estimation problem
(3.2) in the following theorem.

###### Theorem 4.1

For every $h\in\\{0,\cdots,N-1\\}$, the quadratic objective $\mathcal{J}_{h}$
defined in (3.2) is twice continuously differentiable, and its Hessian matrix
can be represented as

$\displaystyle H_{h}$
$\displaystyle=\mathbb{E}_{x_{h},\hat{x}_{h},v_{h}}\begin{bmatrix}\hat{x}_{h}\hat{x}_{h}^{\top}&\hat{x}_{h}y_{h}^{\top}\\\
y_{h}\hat{x}_{h}^{\top}&y_{h}y_{h}^{\top}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\mu_{\hat{x}_{h}}\mu_{\hat{x}_{h}}^{\top}+\Theta&(\mu_{\hat{x}_{h}}\mu_{x_{h}}^{\top}+\Theta)C^{\top}\vspace{0.2em}\\\
C(\mu_{x_{h}}\mu_{\hat{x}_{h}}^{\top}+\Theta)&C(\mu_{x_{h}}\mu_{x_{h}}^{\top}+\Theta)C^{\top}+V\end{bmatrix}>0,$

where we have used $\mu_{\hat{x}_{h}}$ and $\mu_{x_{h}}$ to denote
$\mathbb{E}[\hat{x}_{h}]$ and $\mathbb{E}[x_{h}]$, respectively, and
$\Theta>0$ is the covariance matrix of the zero-mean Gaussian random vector
$\theta_{0}$ in (3.3)-(3.4). Moreover, the objective function
$\mathcal{J}_{h}$ is strongly convex with constant $\lambda_{\min}(H_{h})$ and
smooth with constant $\lambda_{\max}(H_{h})$. Lastly, introducing the
additional random vector $\theta_{0}$ is without any loss of optimality in the
sense that the time-varying KF characterized by (2.8)-(2.9) represents the
unique minimum of (3.2).

The proof of Theorem 4.1 is deferred to Sec. C, where an extended discussion
on the effect of $\theta_{0}$ is also provided. In short, introducing an
additional “small” $\theta_{0}$ in (3.3)-(3.4) ensures the strict convexity of
the quadratic objective $\mathcal{J}_{h}$ with respect to $A_{L_{h}}$, while
$\mathcal{J}_{h}$ is strictly convex in $B_{L_{h}}$ with or without
$\theta_{0}$ due to the condition $V>0$.

Denote $\pi_{t}=\big{[}A_{L_{t}}\mid B_{L_{t}}\big{]}$. We define the analytic
vanilla PG of $\mathcal{J}_{h}$ for every $h\in\\{0,\cdots,N-1\\}$ to
be222Note that one can use any other policy search directions such as natural
PG [3] or least-squares policy iteration [38] to replace vanilla PG in RHPG.

$\displaystyle\nabla_{\pi_{h}}\mathcal{J}_{h}(\pi_{h})=2\Big{[}\pi_{h}(\Psi_{h}+\Delta)-(G_{h}+\Xi)\Big{]},$
(4.1)

where

$\displaystyle\Delta=\bigg{[}\begin{array}[]{c|c}\Theta&\Theta C^{\top}\\\
\hline\cr C\Theta&C\Theta C^{\top}\end{array}\bigg{]},\
\Xi=\big{[}\begin{array}[]{c|c}A\Theta&A\Theta C^{\top}\end{array}\big{]}$
(4.5) $\displaystyle\Psi_{t}=\bigg{[}\begin{array}[]{c|c}\operatorname{{\rm
Var}}(\hat{x}_{t})&\operatorname*{\rm
Cov}(x_{t},\hat{x}_{t})^{\top}C^{\top}\\\ \hline\cr C\operatorname*{\rm
Cov}(x_{t},\hat{x}_{t})&C\operatorname{{\rm
Var}}(x_{t})C^{\top}+V\end{array}\bigg{]}$ (4.8) $\displaystyle
G_{t}=\big{[}\begin{array}[]{c|c}A\operatorname*{\rm
Cov}(x_{t},\hat{x}_{t})&A\operatorname{{\rm
Var}}(x_{t})C^{\top}\end{array}\big{]}$ (4.10)
$\displaystyle\operatorname{{\rm
Var}}(\hat{x}_{t+1})=\pi_{t}\Psi_{t}\pi_{t}^{\top},\ \operatorname{{\rm
Var}}(x_{t+1})=A\operatorname{{\rm Var}}(x_{t})A^{\top}+W$
$\displaystyle\operatorname*{\rm
Cov}(x_{t+1},\hat{x}_{t+1})=A\big{[}\operatorname*{\rm
Cov}(x_{t},\hat{x}_{t})\mid\operatorname{{\rm
Var}}(x_{t})C^{\top}\big{]}\pi_{t}^{\top}$ $\displaystyle\operatorname{{\rm
Var}}(\hat{x}_{0})=\operatorname*{\rm
Cov}(x_{0},\hat{x}_{0})=\bar{x}_{0}\bar{x}_{0}^{\top},\ \operatorname{{\rm
Var}}(x_{0})=\bar{x}_{0}\bar{x}_{0}^{\top}+X_{0}.$

We next define the vanilla PG update as

$\displaystyle\pi_{h}^{\prime}=\pi_{h}-\eta_{h}\cdot\nabla_{\pi_{h}}\mathcal{J}_{h}(\pi_{h}).$
(4.11)

where $\eta_{h}>0$ is a constant stepsize. When the exact PG in (4.1) is not
available, it can be estimated from samples of system trajectories using (two-
point) zeroth-order optimization techniques as described in Algorithm 2. Due
to the landscape properties listed in Theorem 4.1, global convergence and
sample complexity of the PG update (4.11) and its zeroth-order implementation
naturally follow. We present them in the following propositions.

Algorithm 2 Two-Point Zeroth-Order Oracle

0: Iteration index $h$, filter $\pi_{h}$, smoothing radius $r_{h}$

1: Uniformly sample $U_{h}$ from the surface of a unit sphere; Compute
$\pi_{h}^{+}\leftarrow\pi_{h}+r_{h}U_{h}$ and
$\pi_{h}^{-}\leftarrow\pi_{h}-r_{h}U_{h}$

2: Sample $x_{0}$ and compute $x_{h}$ and $\hat{x}_{h}$ by (3.3)-(3.4) using
convergent filters $\\{\pi_{t}\\}_{t\in\\{0,\cdots,h-1\\}}$

3: Generate $x_{h+1}$ and $y_{h}$ using a simulator.

4: Compute estimates of $x_{h+1}$ using $\pi_{h}^{+}$ and $\pi_{h}^{-}$,
denoted as $\hat{x}_{h+1}^{+}$ and $\hat{x}_{h+1}^{-}$, respectively

5: Calculate objective values $J_{h}(\pi_{h}^{+})$ and $J_{h}(\pi_{h}^{-})$ of
(3.2)

6: return
$\widetilde{\nabla}_{\pi_{h}}\mathcal{J}_{h}(\pi_{h})\leftarrow\frac{n(m+n)}{2r_{h}}\cdot\big{[}J_{h}(\pi_{h}^{+})-J_{h}(\pi_{h}^{-})\big{]}\cdot
U_{h}$

###### Proposition 4.2

For all $h\in\\{0,\cdots,N-1\\}$ and a fixed $\epsilon>0$, choose a constant
stepsize $\eta_{h}\leq 1/\lambda_{\max}(H_{h})$. Then, the PG update (4.11)
converges linearly to the unique minimum of (3.2). That is, such that
$\|\pi_{h}^{T_{h}}-\widetilde{\pi}_{h}^{*}\big{\|}\leq\epsilon$ after a total
number of $T_{h}\sim\mathcal{O}(\log(\epsilon^{-1}))$ iterations, where
$\widetilde{\pi}_{h}^{*}$ is the unique minimizing solution of (3.2) after
setting filters for all $t\in\\{0,\cdots,h-2\\}$ to be those computed in the
previous iterations and are $\epsilon$-close to the unique minimum,
respectively.

###### Proposition 4.3

Choose the smoothing radius of the zeroth-order PG to satisfy
$r_{h}\sim\mathcal{O}(\epsilon)$ and the stepsize
$\eta_{h}\sim\mathcal{O}(\epsilon^{2})$. Then, the zeroth-order PG update in
Algorithm 2 converges after
$T_{h}\sim\widetilde{\mathcal{O}}(\epsilon^{-2}\log(\frac{1}{\delta\epsilon^{2}}))$
iterations in the sense that
$\|\pi_{h}^{T_{h}}-\widetilde{\pi}_{h}^{*}\big{\|}\leq\epsilon$ with a
probability of at least $1-\delta$.

Proposition 4.2 is standard, and Proposition 4.3 follows from the proof of
Proposition 3.3 in [17]. Combining Theorem 3.2 with Proposition 4.3, we
conclude that if we spend $\widetilde{O}(\epsilon^{-2})$ samples in solving
every one-step KF problem to $\mathcal{O}(\epsilon)$-accuracy with a
probability of $1-\delta$, for all $h\in\\{0,\cdots,N-1\\}$, then Algorithm 1
is guaranteed to output $\pi_{N-1}$ that is $\epsilon$-close to the infinite-
horizon KF with a probability of at least $1-N\delta$. The total sample
complexity of the RHPG algorithm is thus $\widetilde{O}(\epsilon^{-2})\cdot
O(\log(\epsilon^{-1}))\sim\widetilde{O}(\epsilon^{-2})$. This complexity
result matches the complexity of applying RHPG to the LQR task [17].

##  5 Numerical Experiment: Estimation of the Convection-Diffusion Model

We conducted numerical experiments to design state estimators for the one-
dimensional convection-diffusion linear PDE333We open-source the code
repository at https://github.com/xiangyuan-zhang/LearningKF. The convection-
diffusion equation models physical phenomena involving the transfer of
particles, energy, or other quantities within a system due to convection and
diffusion. These quantities are described by a continuous concentration
function $c(\mathsf{x},t):\Omega\times\mathbb{R}^{+}\rightarrow\mathbb{R}$,
where $\mathsf{x}$ and $t$ represent spatial and temporal coordinates,
respectively, and $\Omega\subset\mathbb{R}$ is the spatial domain of interest.
The one-dimensional convection-diffusion equation can then be expressed as

$\displaystyle\frac{\partial{c}}{\partial
t}=\nu\frac{\partial^{2}c}{\partial\mathsf{x}^{2}}-v\frac{\partial
c}{\partial\mathsf{x}},$ (5.1)

where $\nu$ is the diffusion coefficient and $v$ is the convection velocity;
these scalar physical parameters characterize the strength of convection and
diffusion, respectively. When $v=0$, the convection-diffusion equation (5.1)
reduces to the heat equation. As with any PDE, the convection-diffusion
equation must be accompanied by initial and boundary conditions. Here, we
consider the domain $\Omega=[0,1]$ with periodic boundary conditions, while
the initial condition will be defined shortly.

The convection-diffusion equation can be solved numerically by discretizing
space and time, resulting in a linear state-space model of the form (2.1). To
do so, we define a state vector $x_{t}\in\mathbb{R}^{n}$ that contains the
values of $c$ at $n$ equally-spaced points in $\Omega$ and at time $t=k\Delta
t$, where $n$ is even, $k=0,1,\dots$, and $\Delta t$ is the discrete time
step. The dynamics governed by the convection-diffusion equation can then be
approximated by a state-space model $x_{t+1}=Ax_{t}$ with

$\displaystyle A=\frac{1}{n}\cdot D^{\dagger}{\rm
diag}(e^{-ivk_{\mathsf{x}}-\nu k_{\mathsf{x}}^{2}}\Delta t)D,$ (5.2)

where $k_{\mathsf{x}}=2\pi[0,\dots,n/2-1,0,-n/2+1,\dots,-1]$ is the vector of
spatial wavenumbers, $i$ is the imaginary unit, $D$ is the discrete Fourier
transform (DFT) matrix defined by $D_{pq}=e^{-2\pi i(p-1)(q-1)/n}$, and its
scaled conjugate transpose $D^{\dagger}/n$ is the inverse discrete Fourier
transform (IDFT) matrix [39]. The matrix $A$ combines a spectral
discretization of the spatial derivatives in the convection-diffusion
equation, which takes advantage of the periodicity of the spatial domain with
an exact temporal integration of its continuous-time dynamics. Such spectral
evaluation of the derivatives enjoys exponential convergence properties [40].
As a result, even a small state dimension $n$ yields a state-space model that
faithfully reproduces the dynamics of the convection-diffusion equation.

Figure 4: Eigen-spectrum of $A$ constructed in (5.2), where $\rho(A)=1$ since
the integrated concentration over the domain is conserved by the convection-
diffusion equation.

To set up the numerical experiments, we chose $\Delta t=0.05$, and set the
dimensions of the state vector $x_{t}$ and measurement/observation vector to
$n=200$ and $m=5$, respectively. The five sensors were evenly distributed
across the physical domain $\Omega$, where each sensor measures the (unscaled)
values of the state at the corresponding location, subject to additive zero-
mean Gaussian white noises. Moreover, we chose the diffusion coefficient to
$\nu=2\times 10^{-3}$, the convection velocity to $v=5\times 10^{-2}$, the
distribution of $x_{0}$ to be

$\displaystyle x_{0}\sim\mathcal{N}\Big{(}$ $\displaystyle
c(\mathsf{x},t=0),\frac{1}{16}\cdot\sin(2\pi\mathsf{x})(\sin(2\pi\mathsf{x}))^{\top}\Big{)}$
(5.3) $\displaystyle c(\mathsf{x},t=0)=\mathrm{sech}(10(\mathsf{x}-1/2)).$
(5.4)

Furthermore, we set the covariance matrix of the measurement noise to
$V=10^{-1}\cdot\bm{I}$, and the covariance matrix of the process to
$W=10^{-9}\cdot\bm{I}$. Lastly, we set the covariance matrix of the additional
noise $\theta_{0}$ in (3.3)-(3.4) to $\Theta=10^{-2}\cdot\bm{I}$.

Algorithm 3 First-Order RHPG with Adam

0: Problem horizon $N$, $\Theta$

1: Compute $\Delta$, $\Xi$ by (4.5)

2: for $h\in\\{0,\cdots,N-1\\}$ do

3: Compute $\Psi_{h}$, $G_{h}$ by (4.8)-(4.10)

4: $i_{h}\leftarrow 0$, $m_{h}\leftarrow\bm{0}_{n\times(n+m)}$,
$v_{h}\leftarrow\bm{0}_{(n+m)\times(n+m)}$

5: while $\|\nabla_{\pi_{h}}\|\geq 10^{-4}$ do

6: $\nabla_{\pi_{h}}\leftarrow
2\big{[}\pi_{h}(\Psi_{h}+\Delta)-(G_{h}+\Xi)\big{]}$

7:
$\pi_{h},m_{h},v_{h}\leftarrow\textsc{Adam}(\pi_{h},\nabla_{\pi_{h}},m_{h},v_{h},10^{-3},i_{h})$

8: $i_{h}\leftarrow i_{h}+1$

9: end while

10: if $h\neq N-1$ then

11: $\pi_{h+1}\leftarrow\pi_{h}$ // Warm start the next iteration

12: end if

13: end for

14: return $\pi_{N-1}$

15:

16:
Adam$(P,\nabla,m,v,\eta,i,\beta_{1}=0.9,\beta_{2}=0.999,\varsigma=10^{-8})$

17: $m\leftarrow\beta_{1}m+(1-\beta_{1})\nabla$,
$v\leftarrow\beta_{2}v+(1-\beta_{2})(\nabla^{\top}\nabla)$

18: $\hat{m}\leftarrow m/(1-\beta_{1}^{i})$, $\hat{v}\leftarrow
v/(1-\beta_{2}^{i})$

19: $P\leftarrow P-\eta\cdot\hat{m}(\hat{v}^{1/2}+\varsigma\bm{I})^{-1}$

20: return $P,m,v$

After setting up the PDE environment, we applied the RHPG algorithm to learn
the KF design with several different problem horizons ranging from $N=1$ to
$N=101$. We implemented the inner loop of the RHPG algorithm using the first-
order vanilla PG update (4.11) with the stepsize selected based on the Adam
rule [41]. We detail the procedure of the first-order RHPG algorithm in
Algorithm 3.

Figure 5: Visualizing the estimated trajectories by RL-based filters through
RHPG and the model-based infinite-horizon KF (2.3)-(2.4), along with the
ground truth system trajectory. The $x$, $y$, and $z$ axes of the three-
dimensional plots represent time ($700$ total steps), space ($200$
dimensions), and the normalized concentration of the convection-diffusion
equation. With large $N$, the learned filter converges globally to the
infinite-horizon KF.

In the first experiment, we generated a ground truth system trajectory with
the length of $700$ discrete time steps and with a deterministic initial
condition being $x_{0}=c(\mathsf{x},t=0)$ in (5.4). Then, we visualized the
estimated state trajectories computed using the convergent filters from RHPG
with several different input horizons $N\in\\{1,2,3,6,11,21,31,51,101\\}$),
where for all filters the initial state estimate was set to
$\hat{x}_{0}=c(\mathsf{x},t=0)$. In Figure 5, we provide a comparison of the
estimated trajectories by the RL-based filters with the ground truth
trajectory and the estimated trajectory generated by the infinite-horizon KF
(2.3)-(2.4). One can observe from the first row of Figure 5 that when $N$ is
small, the convergent filters of RHPG are myopic and can only predict system
states over a short period since the convergences of RHPG in the initial
iterations are toward the time-varying KFs in (2.8). As we increase the number
of iterations of RHPG, the filter adapts to new tasks with longer and longer
problem horizons (cf., the illustration in Figure 2). Lastly, one can observe
from the third row of Figure 5 that when the problem horizon becomes
sufficiently large, the RHPG algorithm converges globally to the infinite-
horizon KF, which corroborates the theories developed in the paper.

Figure 6: Average state estimation error $\|x_{t}-\hat{x}_{t}\|_{2}$ of RHPG
filters and the model-based KF over $100$ different system trajectories. Each
trajectory starts from a random initial condition and is affected by randomly
sampled process and measurement noises.

In the second experiment, we generated $100$ random system trajectories of the
convection-diffusion equation for $700$ time steps, each starting from a
randomly sampled initial condition according to (5.3). We then applied the
convergent filters of RHPG as well as the model-based KF to estimate the $100$
system trajectories and plotted the average estimation error (i.e.,
$\|x_{t}-\hat{x}_{t}\|_{2}$) over time. As shown in Figure 6, the RHPG filters
with a small problem horizon $N$ are myopic, resulting in lower estimation
costs over the short time period. However, these myopic filters are not able
to regulate the state estimation error in the asymptotic regime. When we run
RHPG for a sufficiently large number of steps (e.g., $51$ and $101$), the
learned filter performs well in the asymptotic regime due to the global
convergence of RHPG to the infinite-horizon KF.

##  6 Conclusion

We have introduced the RHPG algorithm and provided rigorous analyses for its
convergence and sample complexity in learning the infinite-horizon KF. RHPG is
the first model-free PG algorithm with provable global convergence in learning
the KF design, and it does not require any prior knowledge of the system for
initialization and does not require the target system to be open-loop stable.
We have validated our theories in extensive numerical experiments on a large-
scale convection-diffusion model. On a higher level, our work has proposed
developing RL algorithms specifically for control applications by utilizing
classic control theory in the algorithmic design. The proposed paradigm
enables certifying provable performance guarantees in a purely model-free
setting, overcoming the nonconvex optimization landscape. Following this work
and the dual theory to LQR [17], several ongoing and future research
directions include designing and analyzing RHPG-type algorithms for LQG and
$\mathcal{H}_{\infty}$-robust filtering.

## Acknowledgment

X. Zhang and T. Başar were supported in part by the US Army Research
Laboratory (ARL) Cooperative Agreement W911NF-17-2-0181, and in part by the
Army Research Office (ARO) MURI Grant AG285. S. Mowlavi and M. Benosman were
supported solely by Mitsubishi Electric Research Laboratories. X. Zhang
acknowledges helpful discussions with Bin Hu of UIUC in the early stage of the
project and with Arvind Raghunathan of MERL. X. Zhang and T. Başar acknowledge
anonymous reviewers of ACC ’23 for their helpful comments.

## References

  * [1] X. Zhang, B. Hu, and T. Başar, “Learning the Kalman filter with fine-grained sample complexity,” in _American Control Conference_ , 2023, pp. 4549–4554.
  * [2] R. S. Sutton, D. A. McAllester, S. P. Singh, and Y. Mansour, “Policy gradient methods for reinforcement learning with function approximation,” in _Advances in Neural Information Processing Systems_ , 2000, pp. 1057–1063.
  * [3] S. M. Kakade, “A natural policy gradient,” in _Advances in Neural Information Processing Systems_ , 2002, pp. 1531–1538.
  * [4] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust region policy optimization,” in _International Conference on Machine Learning_ , 2015, pp. 1889–1897.
  * [5] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [6] J. Schulman, P. Moritz, S. Levine, M. Jordan, and P. Abbeel, “High-dimensional continuous control using generalized advantage estimation,” _arXiv preprint arXiv:1506.02438_ , 2015.
  * [7] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.
  * [8] B. Recht, “A tour of reinforcement learning: The view from continuous control,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 2, pp. 253–279, 2019.
  * [9] B. D. Anderson and J. B. Moore, _Optimal Filtering_. Prentice-Hall, 1979.
  * [10] ——, _Optimal Control: Linear Quadratic Methods_. Prentice-Hall, Inc., 1990.
  * [11] M. Fazel, R. Ge, S. M. Kakade, and M. Mesbahi, “Global convergence of policy gradient methods for the linear quadratic regulator,” in _International Conference on Machine Learning_ , 2018, pp. 1467–1476.
  * [12] 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_ , vol. 67, no. 5, pp. 2435–2450, 2021.
  * [13] D. Malik, A. Pananjady, K. Bhatia, K. Khamaru, P. L. Bartlett, and M. J. Wainwright, “Derivative-free methods for policy optimization: Guarantees for linear quadratic systems,” _Journal of Machine Learning Research_ , vol. 21, no. 21, pp. 1–51, 2020.
  * [14] B. Hambly, R. Xu, and H. Yang, “Policy gradient methods for the noisy linear quadratic regulator over a finite horizon,” _SIAM Journal on Control and Optimization_ , vol. 59, no. 5, pp. 3359–3391, 2021.
  * [15] J. C. Perdomo, J. Umenberger, and M. Simchowitz, “Stabilizing dynamical systems via policy gradient methods,” in _Advances in Neural Information Processing Systems_ , 2021, pp. 29 274–29 286.
  * [16] C. Ju, G. Kotsalis, and G. Lan, “A model-free first-order method for linear quadratic regulator with $\widetilde{O}(1/\varepsilon)$ sampling complexity,” _arXiv preprint arXiv:2212.00084_ , 2022.
  * [17] X. Zhang and T. Başar, “Revisiting LQR control from the perspective of receding-horizon policy gradient,” _IEEE Control Systems Letters_ , vol. 7, pp. 1664–1669, 2023.
  * [18] Y. Li, Y. Tang, R. Zhang, and N. Li, “Distributed reinforcement learning for decentralized linear quadratic control: A derivative-free policy optimization approach,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 12, pp. 6429–6444, 2021.
  * [19] L. Furieri, Y. Zheng, and M. Kamgarpour, “Learning the globally optimal distributed LQ regulator,” in _Learning for Dynamics and Control_ , 2020, pp. 287–297.
  * [20] B. Gravell, P. M. Esfahani, and T. Summers, “Learning robust controllers for linear quadratic systems with multiplicative noise via policy gradient,” _arXiv preprint arXiv:1907.03680_ , 2019.
  * [21] K. Zhang, B. Hu, and T. Başar, “Policy optimization for $\mathcal{H}_{2}$ linear control with $\mathcal{H}_{\infty}$ robustness guarantee: Implicit regularization and global convergence,” _SIAM Journal on Control and Optimization_ , vol. 59, no. 6, pp. 4081–4109, 2021.
  * [22] K. Zhang, X. Zhang, B. Hu, and T. Başar, “Derivative-free policy optimization for linear risk-sensitive and robust control design: Implicit regularization and sample complexity,” in _Advances in Neural Information Processing Systems_ , 2021, pp. 2949–2964.
  * [23] R. E. Kalman, “A new approach to linear filtering and prediction problems,” _Journal of Basic Engineering_ , vol. 82, no. 1, pp. 35–45, 1960.
  * [24] T. Başar, _Control Theory: Twenty-five Seminal Papers (1931-1981)_. IEEE Press New York, 2001\.
  * [25] K. J. Astrom, _Introduction to Stochastic Control Theory_. Elsevier, 1971.
  * [26] R. E. Kalman, “On the general theory of control systems,” in _Proceedings First International Conference on Automatic Control, Moscow, USSR_ , 1960, pp. 481–492.
  * [27] Y. Tang, Y. Zheng, and N. Li, “Analysis of the optimization landscape of linear quadratic Gaussian (LQG) control,” _Mathematical Programming_ , 2023. [Online]. Available: https://doi.org/10.1007/s10107-023-01938-4
  * [28] J. Umenberger, M. Simchowitz, J. C. Perdomo, K. Zhang, and R. Tedrake, “Globally convergent policy search over dynamic filters for output estimation,” _arXiv preprint arXiv:2202.11659_ , 2022.
  * [29] S. Chan, G. Goodwin, and K. Sin, “Convergence properties of the Riccati difference equation in optimal filtering of nonstabilizable systems,” _IEEE Transactions on Automatic Control_ , vol. 29, no. 2, pp. 110–118, 1984\.
  * [30] B. Hassibi, A. H. Sayed, and T. Kailath, _Indefinite-Quadratic Estimation and Control: A Unified Approach to $H_{2}$ and $H_{\infty}$ Theories_. SIAM, 1999.
  * [31] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, “Curriculum learning,” in _International Conference on Machine Learning_ , 2009, pp. 41–48.
  * [32] X. Wang, Y. Chen, and W. Zhu, “A survey on curriculum learning,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 44, no. 9, pp. 4555–4576, 2021.
  * [33] P. Soviany, R. T. Ionescu, P. Rota, and N. Sebe, “Curriculum learning: A survey,” _International Journal of Computer Vision_ , vol. 130, no. 6, pp. 1526–1565, 2022.
  * [34] L. Pinto, J. Davidson, R. Sukthankar, and A. Gupta, “Robust adversarial reinforcement learning,” in _International Conference on Machine Learning_ , 2017, pp. 2817–2826.
  * [35] K. Zhang, B. Hu, and T. Başar, “On the stability and convergence of robust adversarial reinforcement learning: A case study on linear quadratic systems,” in _Advances in Neural Information Processing Systems_ , 2020, pp. 22 056–22 068.
  * [36] X. Zhang, R. K. Velicheti, and T. Başar, “Learning minimax-optimal terminal state estimators and smoothers,” in _22nd IFAC World Congress_ , 2023, pp. 12 391–12 396.
  * [37] L. Cui, T. Başar, and Z.-P. Jiang, “A reinforcement learning look at risk-sensitive linear quadratic Gaussian control,” in _Learning for Dynamics and Control Conference_ , 2023, pp. 534–546.
  * [38] M. G. Lagoudakis and R. Parr, “Least-squares policy iteration,” _The Journal of Machine Learning Research_ , vol. 4, pp. 1107–1149, 2003.
  * [39] K. R. Rao and P. C. Yip, _The Transform and Data Compression Handbook_. CRC Press, 2018.
  * [40] L. N. Trefethen, “Finite difference and spectral methods for ordinary and partial differential equations,” 1996. [Online]. Available: http://people.maths.ox.ac.uk/trefethen/pdetext.html
  * [41] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [42] C. De Souza, “Monotonicity and stabilizability results for the solutions of the Riccati difference equation,” in _Workshop on the Riccati Equation in Control, Systems and Signals_ , 1989, pp. 38–41.
  * [43] R. R. Bitmead, M. R. Gevers, I. R. Petersen, and R. J. Kaye, “Monotonicity and stabilizability-properties of solutions of the Riccati difference equation: Propositions, lemmas, theorems, fallacious conjectures and counterexamples,” _Systems & Control Letters_, vol. 5, no. 5, pp. 309–315, 1985.
  * [44] R. A. Horn and C. R. Johnson, _Matrix Analysis_. Cambridge University Press, 2012.

##  Appendix A Proof of Theorem 3.1

We first present a technical lemma due to [42].

###### Lemma A.1

Consider two RDEs

$\displaystyle\Sigma_{t+1}^{1}=W+A\Sigma_{t}^{1}A^{\top}-A\Sigma_{t}^{1}C^{\top}(V+C\Sigma_{t}^{1}C^{\top})^{-1}C\Sigma_{t}^{1}A^{\top},$
$\displaystyle\Sigma_{t+1}^{2}=W+A\Sigma_{t}^{2}A^{\top}-A\Sigma_{t}^{2}C^{\top}(V+C\Sigma_{t}^{2}C^{\top})^{-1}C\Sigma_{t}^{2}A^{\top}.$

Then, the difference between the two solutions,
$\widetilde{\Sigma}_{t}:=\Sigma^{2}_{t}-\Sigma^{1}_{t}$, for all $t$,
satisfies

$\displaystyle\widetilde{\Sigma}_{t+1}=\overline{A}_{t}\widetilde{\Sigma}_{t}\overline{A}^{\top}_{t}-\overline{A}_{t}\widetilde{\Sigma}_{t}C^{\top}(\widetilde{V}_{t}+C\widetilde{\Sigma}_{t}C^{\top})^{-1}C\widetilde{\Sigma}_{t}\overline{A}^{\top}_{t},$
(A.1)

where $\widetilde{V}_{t}=V+C\Sigma^{1}_{t}C^{\top}$ and
$\overline{A}_{t}:=A-A\Sigma^{1}_{t}C^{\top}\widetilde{V}_{t}^{-1}C$.

Next, identify $\Sigma^{1}_{t}$ with $\Sigma^{*}$ and $\Sigma^{2}_{t}$ with
$\Sigma^{*}_{t}$ in Lemma A.1. Then,
$\widetilde{V}_{t}=V+C\Sigma^{*}C^{\top}=:\widetilde{V}$ and
$\overline{A}_{t}=A-A\Sigma^{*}C^{\top}\widetilde{V}^{-1}C=:\overline{A}$ for
all $t$. Invoking Lemma A.1 leads to

$\displaystyle\widetilde{\Sigma}_{t+1}$
$\displaystyle=\overline{A}\widetilde{\Sigma}_{t}\overline{A}^{\top}-\overline{A}\widetilde{\Sigma}_{t}C^{\top}(\widetilde{V}+C\widetilde{\Sigma}_{t}C^{\top})^{-1}C\widetilde{\Sigma}_{t}\overline{A}^{\top}$
(A.2)
$\displaystyle=\overline{A}\widetilde{\Sigma}_{t}^{1/2}\big{[}\bm{I}+\widetilde{\Sigma}^{1/2}_{t}C^{\top}\widetilde{V}^{-1}C\widetilde{\Sigma}_{t}^{1/2}\big{]}^{-1}\widetilde{\Sigma}_{t}^{1/2}\overline{A}^{\top}$
$\displaystyle\leq\big{[}1+\lambda_{\min}(\widetilde{\Sigma}^{1/2}_{t}C^{\top}\widetilde{V}^{-1}C\widetilde{\Sigma}_{t}^{1/2})\big{]}^{-1}\overline{A}\widetilde{\Sigma}_{t}\overline{A}^{\top}=:\mu_{t}\overline{A}\widetilde{\Sigma}_{t}\overline{A}^{\top},$
(A.3)

where $\widetilde{\Sigma}_{t}^{1/2}$ denotes the unique psd square root of the
psd matrix $\widetilde{\Sigma}_{t}$, $0<\mu_{t}\leq 1$ for all $t$, and
$\overline{A}$ is the closed-loop matrix of the optimal infinite-horizon KF
that has all its eigenvalue inside the unit circle (i.e.,
$\rho(\overline{A})<1$). Next, we use $\|\cdot\|_{*}$ to represent the
$\Sigma^{*}$-induced matrix norm defined as $\|X\|^{2}_{*}:=\max_{z\neq
0}\frac{z^{\top}X^{\top}\Sigma^{*}Xz}{z^{\top}\Sigma^{*}z}$. Then, we invoke
Theorem 14.4.1 of [30], where our $\widetilde{\Sigma}_{t}$, $\overline{A}$ and
$\Sigma^{*}$ correspond to $P_{i}-P^{*}$, $F_{p}$ and $W$ in [30],
respectively. By Theorem 14.4.1 of [30] and (A.3), we obtain
$\|\overline{A}\|_{*}<1$ and given that $\mu_{t}\leq 1$,
$\|\widetilde{\Sigma}_{t+1}\|_{*}\leq\|\overline{A}\|^{2}_{*}\cdot\|\widetilde{\Sigma}_{t}\|_{*}$.
Therefore, the convergence rate is exponential in the sense that
$\|\widetilde{\Sigma}_{t}\|_{*}\leq\|\overline{A}\|_{*}^{2t}\cdot\|\widetilde{\Sigma}_{0}\|_{*}$.
Next, recall the condition number of a matrix $X$ is defined as
$\kappa_{X}:=\sigma_{\max}(X)/\sigma_{\min}(X)$. Then, the convergence of
$\widetilde{\Sigma}_{t}$ to the zero matrix in spectral norm can be
characterized as

$\displaystyle\|\widetilde{\Sigma}_{t}\|\leq\kappa_{\Sigma^{*}}\cdot\|\widetilde{\Sigma}_{t}\|_{*}\leq\kappa_{\Sigma^{*}}\cdot\|\overline{A}\|_{*}^{2t}\cdot\|\widetilde{\Sigma}_{0}\|_{*}.$

In other words, to ensure $\|\widetilde{\Sigma}_{N}\|\leq\epsilon$, it
suffices to require

$\displaystyle
N\geq\frac{1}{2}\cdot\frac{\log\big{(}\frac{\|\widetilde{\Sigma}_{0}\|_{*}\cdot\kappa_{\Sigma^{*}}}{\epsilon}\big{)}}{\log\big{(}\frac{1}{\|\overline{A}\|_{*}}\big{)}}.$
(A.4)

Furthermore, since $(A,W^{1/2})$ is controllable and $(C,A)$ is observable,
and if further $X_{0}>\Sigma^{*}$, the closed-loop system at any time $t\geq
0$ is exponentially asymptotically stable such that the time-invariant
(frozen) filter satisfies
$\rho(A-A\Sigma^{*}_{t}C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}C)<1$ [43, 42].
Lastly, we show that the (monotonic) convergence of the filter gain to the
time-invariant Kalman gain follows from the convergence of $\Sigma^{*}_{t}$ to
$\Sigma^{*}$, which can be verified through:

$\displaystyle L^{*}_{t}-L^{*}$
$\displaystyle=A\Sigma^{*}_{t}C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}-A\Sigma^{*}C^{\top}(V+C\Sigma^{*}C^{\top})^{-1}$
$\displaystyle=A\Sigma^{*}C^{\top}\big{[}(V+C\Sigma^{*}_{t}C^{\top})^{-1}-(V+C\Sigma^{*}C^{\top})^{-1}\big{]}+A(\Sigma^{*}_{t}-\Sigma^{*})C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}$
$\displaystyle=A\Sigma^{*}C^{\top}(V+C\Sigma^{*}C^{\top})^{-1}C(\Sigma^{*}-\Sigma^{*}_{t})C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}-A(\Sigma^{*}-\Sigma^{*}_{t})C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}$
$\displaystyle=(L^{*}C-A)(\Sigma^{*}-\Sigma^{*}_{t})C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}.$
(A.5)

Hence, we have
$\|L^{*}_{t}-L^{*}\|\leq\frac{\|\overline{A}\|\cdot\|C\|}{\lambda_{\min}(V)}\cdot\|\Sigma^{*}_{t}-\Sigma^{*}\|$.
Substituting $\epsilon$ in (A.4) with
$\frac{\epsilon\cdot\lambda_{\min}(V)}{\|\overline{A}\|\cdot\|C\|}$ and
identifying that $\overline{A}$ is exactly $A_{L}^{*}$ completes the proof.

##  Appendix B Proof of Theorem 3.2

To prove $\big{\|}[\widetilde{A}_{L_{N-1}}\
\widetilde{B}_{L_{N-1}}]-[A^{*}_{L}\ B^{*}_{L}]\big{\|}\leq\epsilon$, it
suffices to bound the error between the approximated filter gain
$\widetilde{L}_{N-1}$ and the exact Kalman gain $L^{*}$ as in (2.4). First,
according to Theorem 3.1, we select

$\displaystyle
N=\frac{1}{2}\cdot\frac{\log\big{(}\frac{2\|X_{0}-\Sigma^{*}\|_{*}\cdot\kappa_{\Sigma^{*}}\cdot\|A_{L}^{*}\|\cdot\|C\|}{\epsilon\cdot\lambda_{\min}(V)}\big{)}}{\log\big{(}\frac{1}{\|A_{L}^{*}\|_{*}}\big{)}}+1.$
(B.1)

which ensures that $L^{*}_{N-1}$ is stabilizing and
$\|L^{*}_{N-1}-L^{*}\|\leq\epsilon/2$. Then, it remains to show that Algorithm
1 returns a filter $\widetilde{L}_{N-1}$ such that
$\|\widetilde{L}_{N-1}-L^{*}_{N-1}\|\leq\epsilon/2$.

Recall that the FRDE is the following forward iteration starting with
$\Sigma^{*}_{0}=X_{0}>0$:

$\displaystyle\Sigma^{*}_{t+1}$
$\displaystyle=A\Sigma^{*}_{t}A^{\top}-A\Sigma^{*}_{t}C^{\top}(V+C\Sigma^{*}_{t}C^{\top})^{-1}C\Sigma^{*}_{t}A^{\top}+W$
(B.2) $\displaystyle=(A-L^{*}_{t}C)\Sigma^{*}_{t}A^{\top}+W$ (B.3)
$\displaystyle=(A-L^{*}_{t}C)\Sigma^{*}_{t}(A-L^{*}_{t}C)^{\top}+L^{*}_{t}V(L^{*}_{t})^{\top}+W.$
(B.4)

Moreover, for an arbitrary $L_{t}$, it holds that:

$\displaystyle\Sigma_{t+1}=(A-L_{t}C)\Sigma_{t}(A-L_{t}C)^{\top}+L_{t}VL_{t}^{\top}+W.$
(B.5)

Furthermore, for clarity of the proof, we define/recall:

$\displaystyle L^{*}_{t}\text{: Exact Kalman gain at time $t$ defined in
\eqref{eqn:kalman_gain_finite}}$ $\displaystyle\widetilde{L}_{t}^{*}\text{:
Optimal gain of the current cost-to-come function,}$ absorbing errors in prior
steps $\displaystyle\widetilde{L}_{t}\text{: An approximation of
$\widetilde{L}_{t}^{*}$ obtained by applying \eqref{eqn:KF_update}}$
$\displaystyle\delta_{t}:=\widetilde{L}_{t}-\widetilde{L}_{t}^{*}\text{:
Policy optimization error at time $t$}$
$\displaystyle\widetilde{\Sigma}^{*}_{t+1}\text{: Solution generated by
\eqref{eqn:filter_RDE_Lya} with $L^{*}_{t}=\widetilde{L}^{*}_{t}$ and
$\Sigma^{*}_{t}=\widetilde{\Sigma}_{t}$.}$

We argue that $\|\widetilde{L}_{N-1}-L^{*}_{N-1}\|\leq\epsilon/2$ can be
achieved by carefully controlling $\delta_{t}$ for all $t$. At $t=N-1$, it
holds that

$\displaystyle\|\widetilde{L}_{N-1}-L^{*}_{N-1}\|$
$\displaystyle\leq\|\widetilde{L}^{*}_{N-1}-L^{*}_{N-1}\|+\|\delta_{N-1}\|,$

where from (A.5) we have

$\displaystyle\widetilde{L}^{*}_{N-1}-L^{*}_{N-1}=(L^{*}_{N-1}C-A)(\Sigma^{*}_{N-1}-\widetilde{\Sigma}_{N-1})C^{\top}\cdot(V+C\widetilde{\Sigma}_{N-1}C^{\top})^{-1}.$

As a result, it holds that

$\displaystyle\|\widetilde{L}^{*}_{N-1}-L^{*}_{N-1}\|\leq\frac{\|A_{L_{N-1}}^{*}\|\|C\|}{\lambda_{\min}(V)}\cdot\|\Sigma^{*}_{N-1}-\widetilde{\Sigma}_{N-1}\|.$
(B.6)

Define the helper constant

$\displaystyle
C_{1}:=\frac{\varphi\|C\|}{\lambda_{\min}(V)}>0,\quad\varphi:=\max_{t\in\\{0,\cdots,N-1\\}}\|A^{*}_{L_{t}}\|.$

Next, require $\|\delta_{N-1}\|\leq\frac{\epsilon}{4}$ and
$\|\widetilde{L}^{*}_{N-1}-L^{*}_{N-1}\|\leq\frac{\epsilon}{4}$ to fulfill
$\|\widetilde{L}_{N-1}-L^{*}_{N-1}\|\leq\frac{\epsilon}{2}$. By (B.6), this is
equivalent to requiring

$\displaystyle\|\Sigma^{*}_{N-1}-\widetilde{\Sigma}_{N-1}\|\leq\frac{\epsilon}{4C_{1}}.$
(B.7)

Subsequently, by (B.5), we have

$\displaystyle\Sigma^{*}_{N-1}-\widetilde{\Sigma}_{N-1}=(\Sigma^{*}_{N-1}-\widetilde{\Sigma}^{*}_{N-1})+(\widetilde{\Sigma}^{*}_{N-1}-\widetilde{\Sigma}_{N-1}),$
(B.8)

where the first difference term on the RHS of (B.8) is

$\displaystyle\Sigma^{*}_{N-1}-\widetilde{\Sigma}^{*}_{N-1}=(A-L^{*}_{N-2}C)\Sigma^{*}_{N-2}A^{\top}-(A-\widetilde{L}^{*}_{N-2}C)\widetilde{\Sigma}_{N-2}A^{\top}$
$\displaystyle=(A-L^{*}_{N-2}C)(\Sigma^{*}_{N-2}-\widetilde{\Sigma}_{N-2})A^{\top}+(\widetilde{L}^{*}_{N-2}-L^{*}_{N-2})C\widetilde{\Sigma}_{N-2}A^{\top}.$
(B.9)

Moreover, the second term on the RHS of (B.8) is

$\displaystyle\hskip
9.00002pt\widetilde{\Sigma}^{*}_{N-1}-\widetilde{\Sigma}_{N-1}$
$\displaystyle=(A-\widetilde{L}^{*}_{N-2}C)\widetilde{\Sigma}_{N-2}(A-\widetilde{L}^{*}_{N-2}C)^{\top}+\widetilde{L}^{*}_{N-2}V(\widetilde{L}^{*}_{N-2})^{\top}-(A-\widetilde{L}_{N-2}C)\widetilde{\Sigma}_{N-2}(A-\widetilde{L}_{N-2}C)^{\top}-\widetilde{L}_{N-2}V(\widetilde{L}_{N-2})^{\top}$
$\displaystyle=A\widetilde{\Sigma}_{N-2}C^{\top}\widetilde{L}_{N-2}^{\top}-A\widetilde{\Sigma}_{N-2}C^{\top}(\widetilde{L}^{*}_{N-2})^{\top}+\widetilde{L}_{N-2}C\widetilde{\Sigma}_{N-2}A^{\top}$
$\displaystyle\hskip
9.00002pt+\widetilde{L}^{*}_{N-2}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)(\widetilde{L}^{*}_{N-2})^{\top}-\widetilde{L}_{N-2}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)\widetilde{L}_{N-2}^{\top}-\widetilde{L}^{*}_{N-2}C\widetilde{\Sigma}_{N-2}A^{\top}$
$\displaystyle=\widetilde{L}^{*}_{N-2}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)(\widetilde{L}^{*}_{N-2})^{\top}-A\widetilde{\Sigma}_{N-2}C^{\top}(\widetilde{L}^{*}_{N-2})^{\top}-\widetilde{L}^{*}_{N-2}C\widetilde{\Sigma}_{N-2}A^{\top}+A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1}C\widetilde{\Sigma}_{N-2}A^{\top}$
$\displaystyle-(\widetilde{L}_{N-2}-A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1})(C\widetilde{\Sigma}_{N-2}C^{\top}+V)\cdot(\widetilde{L}_{N-2}-A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1})^{\top}$
(B.10)
$\displaystyle=(A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1}-\widetilde{L}^{*}_{N-2})(C\widetilde{\Sigma}_{N-2}C^{\top}+V)\cdot(A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1}-\widetilde{L}^{*}_{N-2})^{\top}$
$\displaystyle-(\widetilde{L}_{N-2}-A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1})(C\widetilde{\Sigma}_{N-2}C^{\top}+V)\cdot(\widetilde{L}_{N-2}-A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1})^{\top},$
(B.11)

where (B.10) is due to completion of squares. Substituting
$\widetilde{L}^{*}_{N-2}=A\widetilde{\Sigma}_{N-2}C^{\top}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)^{-1}$
into (B.11) leads to

$\displaystyle\widetilde{\Sigma}^{*}_{N-1}-\widetilde{\Sigma}_{N-1}=-\delta_{N-2}(C\widetilde{\Sigma}_{N-2}C^{\top}+V)\delta_{N-2}^{\top}.$
(B.12)

Thus, combining (B.8), (B.9), and (B.12) yields

$\displaystyle\hskip 10.00002pt\|\Sigma^{*}_{N-1}-\widetilde{\Sigma}_{N-1}\|$
$\displaystyle\leq\|\Sigma^{*}_{N-2}-\widetilde{\Sigma}_{N-2}\|\varphi\|A\|+\|\widetilde{L}^{*}_{N-2}-L^{*}_{N-2}\|\|C\|\|\widetilde{\Sigma}_{N-2}\|\|A\|+\|\delta_{N-2}\|^{2}\|C\widetilde{\Sigma}_{N-2}C^{\top}+V\|$
$\displaystyle\leq\|A\|\cdot[\varphi+C_{1}\cdot\|C\|\cdot\|\widetilde{\Sigma}_{N-2}\|]\cdot\|\Sigma^{*}_{N-2}-\widetilde{\Sigma}_{N-2}\|+\|\delta_{N-2}\|^{2}\|C\widetilde{\Sigma}_{N-2}C^{\top}+V\|,$
(B.13)

where the last inequality follows from (B.6). Now, require

$\displaystyle\|\Sigma^{*}_{N-2}-\widetilde{\Sigma}_{N-2}\|$
$\displaystyle\leq\frac{\epsilon}{4C_{1}C_{2}},\quad\|\delta_{N-2}\|\leq\frac{1}{2}\sqrt{\frac{\epsilon}{C_{1}C_{3}}},$
(B.14)

where $C_{2}$ and $C_{3}$ are positive constants defined as

$\displaystyle
C_{2}:=2\|A\|\cdot\big{[}\varphi+C_{1}\cdot\|C\|\cdot\big{(}\|X_{0}\|+\|\Sigma^{*}\|\big{)}\big{]}>0,\quad
C_{3}:=2\big{[}\|V\|+\|C\|^{2}\big{(}\|X_{0}\|+\|\Sigma^{*}\|\big{)}\big{]}>0.$

Then, condition (B.14) is sufficient for (B.7) (and thus for
$\|\widetilde{L}_{N-1}-L^{*}_{N-1}\|\leq\epsilon/2$) to hold. Subsequently, we
can propagate the requirements in (B.14) backward in time. Specifically, we
iteratively apply the arguments in (B.13) (i.e., by plugging quantities with
subscript $t$ into the LHS of (B.13) and plugging quantities with subscript
$t-1$ into the RHS of (B.13)) to obtain the result that if at all
$t\in\\{1,\cdots,N-2\\}$, we require

$\displaystyle\|\Sigma^{*}_{t}-\widetilde{\Sigma}_{t}\|\leq\frac{\epsilon}{4C_{1}C_{2}^{N-t-1}},\
\|\delta_{t}\|\leq\frac{1}{2}\sqrt{\frac{\epsilon}{C_{1}C_{2}^{N-t-2}C_{3}}},$
(B.15)

then (B.14) holds true and therefore (B.7) is satisfied.

We now compute the required accuracy for $\delta_{0}$. As illustrated in
Figure 3, we have $\Sigma^{*}_{1}=\widetilde{\Sigma}^{*}_{1}$ because

$\displaystyle\widetilde{\Sigma}^{*}_{1}$
$\displaystyle=(A-\widetilde{L}^{*}_{0}C)\Sigma^{*}_{0}(A-\widetilde{L}^{*}_{0}C)^{\top}+\widetilde{L}^{*}_{0}V(\widetilde{L}^{*}_{0})^{\top}+W=(A-L^{*}_{0}C)\Sigma^{*}_{0}(A-L^{*}_{0}C)^{\top}+L^{*}_{0}V(L^{*}_{0})^{\top}+W=\Sigma^{*}_{1},$

where the second equality is due to $\widetilde{L}^{*}_{0}=L^{*}_{0}$ since
there are no prior computational errors yet at $t=0$. By (B.13), the distance
between $\Sigma^{*}_{1}$ and $\widetilde{\Sigma}_{1}$ can be bounded as

$\displaystyle\|\Sigma^{*}_{1}-\widetilde{\Sigma}_{1}\|=\|\widetilde{\Sigma}^{*}_{1}-\widetilde{\Sigma}_{1}\|\leq\|\delta_{0}\|^{2}\cdot
C_{3}.$

To fulfill the requirement (B.15) for $t=1$, which is
$\|\Sigma^{*}_{1}-\widetilde{\Sigma}_{1}\|\leq\frac{\epsilon}{4C_{1}C_{2}^{N-2}}$,
it suffices to let

$\displaystyle\|\delta_{0}\|\leq\frac{1}{2}\sqrt{\frac{\epsilon}{C_{1}C_{2}^{N-2}C_{3}}}.$
(B.16)

Lastly, we analyze the worst-case complexity of the proposed algorithm by
computing, at the most stringent case, the required size of $\|\delta_{t}\|$.
When $C_{2}\leq 1$, the most stringent dependence of $\|\delta_{t}\|$ on
$\epsilon$ happens at $t=N-1$, which is of the order $\mathcal{O}(\epsilon)$,
and the dependences on system parameters (through the dependence on constants
$C_{1},C_{2}$ and $C_{3}$) are polynomial. We argue that if $C_{2}>1$, then
the requirement on $\|\delta_{N-1}\|$ is still the most stringent one. This is
because $\|\delta_{0}\|\leq\|\delta_{t}\|$ for all $t\in\\{1,\cdots,N-2\\}$
and by (B.16), we have

$\displaystyle\|\delta_{0}\|\sim\mathcal{O}\Big{(}\sqrt{\frac{\epsilon}{C_{2}^{N-2}}}\Big{)}.$
(B.17)

Since we require $N$ to satisfy (B.1), the dependence of $\|\delta_{0}\|$ on
$\epsilon$ in (B.17) becomes
$\|\delta_{0}\|\sim\mathcal{O}(\epsilon^{\frac{3}{4}})$, which is milder than
that of $\|\delta_{N-1}\|$. Therefore, it suffices to require the most
stringent error bound for all $t$, which is
$\|\delta_{t}\|\sim\mathcal{O}(\epsilon)$, to reach the
$\epsilon$-neighborhood of the infinite-horizon KF. Lastly, for
$\widetilde{A}_{L_{N-1}}$ to be stable, it suffices to let $\epsilon$ to be
sufficiently small such that
$\epsilon<1-\|A_{L}^{*}\|_{*}\Longrightarrow\|\widetilde{A}_{L}\|_{*}<1$. This
completes the proof.

##  Appendix C Proof of Theorem 4.1

We first introduce a standard result in linear algebra, see Eq. (7.7.5) and
Theorem 7.7.7 on pages 495-496 of [44].

###### Lemma C.1

For any symmetric matrix

$\displaystyle X=\begin{bmatrix}X_{11}&X_{12}\\\
X_{12}^{\top}&X_{22}\end{bmatrix}\in\mathbb{R}^{(n+m)\times(n+m)},$

$X$ is pd if and only if $X_{11}>0$ and
$X_{22}-X_{12}^{\top}(X_{11})^{-1}X_{12}>0$, where
$X_{11}\in\mathbb{R}^{n\times n}$, $X_{22}\in\mathbb{R}^{m\times m}$, and
$X_{12}\in\mathbb{R}^{n\times m}$.

We start by analyzing the first iteration of the RHPG algorithm. The objective
function for the one-step KF problem is defined as

$\displaystyle\mathcal{J}_{0}:=\mathbb{E}_{x_{0},w_{0},v_{0}}\Big{\\{}\sum_{t=0}^{1}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\Big{\\}}.$

The one-step static estimation problem, when formulated as a policy
optimization problem, can be represented as

$\displaystyle(\textbf{P1})\ \min_{A_{L_{0}},B_{L_{0}}}\
\mathcal{J}_{0},\quad\text{ s.t. }\
x_{1}=Ax_{0}+w_{0},\quad\hat{x}_{1}=A_{L_{0}}\hat{x}_{0}+B_{L_{0}}y_{0},\quad
x_{0}=\mathcal{N}(\bar{x}_{0},X_{0}),\quad\hat{x}_{0}=\bar{x}_{0}$

The Hessian matrix for the quadratic program $(\textbf{P1})$ is

$\displaystyle
H_{(\textbf{P1})}=\begin{bmatrix}\frac{\partial^{2}\mathcal{J}_{0}}{A_{L_{0}}^{2}}&\frac{\partial^{2}\mathcal{J}_{0}}{A_{L_{0}}\cdot
B_{L_{0}}}\vspace{0.2em}\\\ \frac{\partial^{2}\mathcal{J}_{0}}{B_{L_{0}}\cdot
A_{L_{0}}}&\frac{\partial^{2}\mathcal{J}_{0}}{B_{L_{0}}^{2}}\end{bmatrix}=\begin{bmatrix}\bar{x}_{0}\bar{x}_{0}^{\top}&\bar{x}_{0}\bar{x}_{0}^{\top}C^{\top}\\\
C\bar{x}_{0}\bar{x}_{0}^{\top}&C(\bar{x}_{0}\bar{x}_{0}^{\top}+X_{0})C^{\top}+V\end{bmatrix}.$

Since $\bar{x}_{0}\bar{x}_{0}^{\top}$ is not pd (only one eigenvalue is
positive) and according to Lemma C.1, the quadratic program $(\textbf{P1})$ is
not strictly convex, there exist multiple (in fact, infinite number of)
stationary points that are equally good, and all of these stationary point
policies attain the minimum value of (P1). However, we know that the one-step
filter, denoted as $(A-L_{0}^{*}C,L_{0}^{*})$, consists of one of such
stationary points.

To convexify (P1), we sample $\theta_{0}\sim\mathcal{N}(\bm{0},\Theta)$ and
set $\hat{x}_{0}=\bar{x}_{0}+\theta_{0}$, where $\theta_{0}$ is independent to
$x_{0}$, $w_{0}$, and $v_{0}$. Concretely, we define the second quadratic
program (P2) as follows:

$\displaystyle(\textbf{P2})\ \min_{A_{L_{0}},B_{L_{0}}}\
\mathbb{E}_{x_{0},w_{0},v_{0},\theta_{0}}\Big{\\{}\sum_{t=0}^{1}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\Big{\\}}$
s.t. $\displaystyle
x_{1}=Ax_{0}+w_{0},\quad\hat{x}_{1}=A_{L_{0}}\hat{x}_{0}+B_{L_{0}}y_{0},\quad
x_{0}=\mathcal{N}(\bar{x}_{0},X_{0}),\quad\hat{x}_{0}=\bar{x}_{0}+\theta_{0}$

In other words, the initial state estimate now satisfies
$\hat{x}_{0}\sim\mathcal{N}(\bar{x}_{0},\Theta)$. Then, (P2) can be shown to
be equivalent to

$\displaystyle\min_{A_{L_{0}},B_{L_{0}}}\
\mathbb{E}_{x_{0},\theta_{0}}\Big{\\{}(x_{0}-\bar{x}_{0}-\theta_{0})^{\top}(x_{0}-\bar{x}_{0}-\theta_{0})\Big{\\}}$
$\displaystyle\hskip
10.00002pt+\mathbb{E}_{x_{0},\theta_{0},w_{0},v_{0}}\Big{\\{}(x_{1}-A_{L_{0}}(\bar{x}_{0}+\theta_{0})-B_{L_{0}}y_{0})^{\top}(x_{1}-A_{L_{0}}(\bar{x}_{0}+\theta_{0})-B_{L_{0}}y_{0})\Big{\\}}$
$\displaystyle\equiv\min_{A_{L_{0}},B_{L_{0}}}\
\mathcal{J}_{0}+\mathbb{E}_{\theta_{0}}\Big{\\{}\theta_{0}^{\top}(\bm{I}+A_{L_{0}}^{\top}A_{L_{0}})\theta_{0}\Big{\\}}\equiv\min_{A_{L_{0}},B_{L_{0}}}\
\mathcal{J}_{0}+\text{Tr}\Big{(}(\bm{I}+A_{L_{0}}^{\top}A_{L_{0}})\Theta\Big{)}.$
(C.1)

Notice that $\text{Tr}((\bm{I}+A_{L_{0}}^{\top}A_{L_{0}})\Theta)$ in (C.1) can
be identified as an additional regularization term. The Hessian matrix becomes

$\displaystyle
H_{(\textbf{P2})}=\begin{bmatrix}\bar{x}_{0}\bar{x}_{0}^{\top}+\Theta&\bar{x}_{0}\bar{x}_{0}^{\top}C^{\top}\\\
C\bar{x}_{0}\bar{x}_{0}^{\top}&C(\bar{x}_{0}\bar{x}_{0}^{\top}+X_{0})C^{\top}+V\end{bmatrix}.$

Since $\bar{x}_{0}\bar{x}_{0}^{\top}+\Theta>0$ and

$\displaystyle
C(\bar{x}_{0}\bar{x}_{0}^{\top}+X_{0})C^{\top}+V-C\bar{x}_{0}\bar{x}_{0}^{\top}(\bar{x}_{0}\bar{x}_{0}^{\top}+\Theta)^{-1}\bar{x}_{0}\bar{x}_{0}^{\top}C^{\top}$
$\displaystyle\geq
C(\bar{x}_{0}\bar{x}_{0}^{\top}+X_{0})C^{\top}+V-C\bar{x}_{0}\bar{x}_{0}^{\top}(\bar{x}_{0}\bar{x}_{0}^{\top})^{\dagger}\bar{x}_{0}\bar{x}_{0}^{\top}C^{\top}=CX_{0}C^{\top}+V>0,$

where we have used $\dagger$ to denote the Moore–Penrose matrix inverse, by
Lemma C.1, adding $\theta_{0}$ to the initial estimate $\hat{x}_{0}$ helps
ensuring that (P2) is strictly convex. However, the solution to (P2) will not
be identical to the solution to (P1) due to the extra regularization term.

To achieve unbiased convexification, we inject $\theta_{0}$ to both
$\hat{x}_{0}$ and $x_{0}$, which results in the following optimization problem

$\displaystyle(\textbf{P3})\ \min_{A_{L_{0}},B_{L_{0}}}\
\mathbb{E}_{x_{0},w_{0},v_{0},\theta_{0}}\Big{\\{}\sum_{t=0}^{1}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\Big{\\}}$
s.t. $\displaystyle
x_{1}=Ax_{0}+w_{0},\quad\hat{x}_{1}=A_{L_{0}}\hat{x}_{0}+B_{L_{0}}y_{0},\quad
x_{0}=\mathcal{N}(\bar{x}_{0},X_{0})+\theta_{0},\quad\hat{x}_{0}=\bar{x}_{0}+\theta_{0}$

Let $\Lambda_{0}:=A-B_{L_{0}}C-A_{L_{0}}$. Then, we can show that (P3) is
equivalent to

$\displaystyle\min_{A_{L_{0}},B_{L_{0}}}\mathbb{E}_{x_{0},\theta_{0}}\Big{\\{}(x_{0}-\bar{x}_{0}+\theta_{0}-\theta_{0})^{\top}(x_{0}-\bar{x}_{0}+\theta_{0}-\theta_{0})\Big{\\}}$
$\displaystyle\hskip
15.00002pt+\mathbb{E}\Big{\\{}(Ax_{0}+w_{0}-A_{L_{0}}\bar{x}_{0}-B_{L_{0}}(Cx_{0}+v_{0})+\Lambda_{0}\theta_{0})^{\top}(Ax_{0}+w_{0}-A_{L_{0}}\bar{x}_{0}-B_{L_{0}}(Cx_{0}+v_{0})+\Lambda_{0}\theta_{0}\Big{\\}}$
$\displaystyle\equiv\min_{A_{L_{0}},B_{L_{0}}}\mathcal{J}_{0}+\mathbb{E}_{\theta_{0}}\Big{\\{}\theta_{0}^{\top}\Lambda_{0}^{\top}\Lambda_{0}\theta_{0}\Big{\\}}\equiv\min_{A_{L_{0}},B_{L_{0}}}\mathcal{J}_{0}+\text{Tr}(\Lambda_{0}^{\top}\Lambda_{0}\Theta).$
(C.2)

The Hessian matrix of the quadratic program (P3) is

$\displaystyle
H_{(\textbf{P3})}=\begin{bmatrix}\bar{x}_{0}\bar{x}_{0}^{\top}+\Theta&(\bar{x}_{0}\bar{x}_{0}^{\top}+\Theta)C^{\top}\\\
C(\bar{x}_{0}\bar{x}_{0}^{\top}+\Theta)&C(\bar{x}_{0}\bar{x}_{0}^{\top}+\Theta+X_{0})C^{\top}+V\end{bmatrix}.$

Injecting $\theta_{0}$ into $x_{0}$ leads to a different “regularization term”
as can be seen in (C.2). By Lemma C.1, it is immediate to check that (P3) is
strictly convex, and hence has a unique minimum, and is also smooth.
Subsequently, we know that the one-step filter $(A-L_{0}^{*}C,L_{0}^{*})$
attains the minimum of $\mathcal{J}_{0}$; it also achieves a zero value for
the regularization term $\text{Tr}(\Lambda_{0}^{\top}\Lambda_{0}\Theta)$ in
(C.2), which is non-negative for any $\Theta>0$. Hence,
$(A-L_{0}^{*}C,L_{0}^{*})$ corresponds to the unique minimum of
$(\textbf{P3})$. This proves that injecting $\theta_{0}$ into both $x_{0}$ and
$\hat{x}_{0}$ results in the strictly convex quadratic problem $(\textbf{P3})$
whose unique minimum is exactly the one-step filter
$(A-L_{0}^{*}C,L_{0}^{*})$.

We proceed with the induction step. In the $\tau$-th iteration of the RHPG
algorithm, for an arbitrary $\tau>0$, assume that we have computed time-
varying filters
$\\{A_{L_{0}}^{*},\cdots,A_{L_{\tau-1}}^{*},B_{L_{0}}^{*},\cdots,B_{L_{\tau-1}}^{*}\\}$.
The quadratic program for the $\tau$-th iteration is constructed such that
$\theta_{0}$ is injected into both $\hat{x}_{\tau}$ and $x_{\tau}$, but not
$\hat{x}_{0}$ and $x_{0}$. We define

$\displaystyle\mathcal{J}_{\tau}:=\mathbb{E}_{x_{0},w_{t},v_{t}}\Big{\\{}\sum_{t=0}^{\tau+1}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\Big{\\}}.$

An associated quadratic program (Pτ) can be represented as

$\displaystyle(\textbf{P}_{\tau})$
$\displaystyle\min_{A_{L_{\tau}},B_{L_{\tau}}}\mathbb{E}_{x_{0},w_{t},v_{t},\theta_{0}}\Big{\\{}\sum_{t=0}^{\tau+1}(x_{t}-\hat{x}_{t})^{\top}(x_{t}-\hat{x}_{t})\Big{\\}},$
s.t.
$\displaystyle\hat{x}_{t+1}=A_{L_{t}}^{*}\hat{x}_{t}+B_{L_{t}}^{*}y_{t},\quad\forall
t\in\\{0,\cdots,\tau-2\\},\quad\hat{x}_{0}=\bar{x}_{0}$ $\displaystyle
x_{t+1}=Ax_{t}+w_{t},\quad\forall t\in\\{0,\cdots,\tau-2\\},\
x_{0}\sim\mathcal{N}(\bar{x}_{0},X_{0})$
$\displaystyle\hat{x}_{\tau}=A_{L_{\tau-1}}^{*}\hat{x}_{\tau-1}+B_{L_{\tau-1}}^{*}y_{\tau-1}+\theta_{0},\quad
x_{\tau}=Ax_{\tau-1}+w_{\tau-1}+\theta_{0}.$

where the expectation in $\mathcal{J}_{\tau}$ is taken over $x_{\tau}$,
$\hat{x}_{\tau}$, $w_{\tau}$, and $v_{\tau}$. The Hessian matrix of the
quadratic program (Pτ) is

$\displaystyle H_{(\textbf{P}_{\tau})}$
$\displaystyle=\mathbb{E}_{x_{\tau},\hat{x}_{\tau},v_{\tau}}\begin{bmatrix}\hat{x}_{\tau}\hat{x}_{\tau}^{\top}&\hat{x}_{\tau}y_{\tau}^{\top}\\\
y_{\tau}\hat{x}_{\tau}^{\top}&y_{\tau}y_{\tau}^{\top}\end{bmatrix}=\begin{bmatrix}\mu_{\hat{x}_{\tau}}\mu_{\hat{x}_{\tau}}^{\top}+\Theta&(\mu_{\hat{x}_{\tau}}\mu_{x_{\tau}}^{\top}+\Theta)C^{\top}\vspace{0.2em}\\\
C(\mu_{x_{\tau}}\mu_{\hat{x}_{\tau}}^{\top}+\Theta)&C(\mu_{x_{\tau}}\mu_{x_{\tau}}^{\top}+\Theta)C^{\top}+V\end{bmatrix},$

where we have used $\mu_{\hat{x}_{\tau}}$ and $\mu_{x_{\tau}}$ to denote
$\mathbb{E}[\hat{x}_{\tau}]$ and $\mathbb{E}[x_{\tau}]$, respectively. By
Lemma C.1, the Hessian matrix $H_{(\textbf{P}_{\tau})}$ is pd for any
$\tau>0$, which proves that $(\textbf{P}_{\tau})$ is strictly convex and
smooth. Furthermore, we can apply a derivation similar to (C.2) to show that
$(\textbf{P}_{\tau})$ is equivalent to

$\displaystyle\min_{A_{L_{\tau}},B_{L_{\tau}}}\mathcal{J}_{\tau}+\text{Tr}(\Lambda_{\tau}^{\top}\Lambda_{\tau}\Theta),\quad\Lambda_{\tau}:=A-B_{L_{\tau}}C-A_{L_{\tau}}.$

The unique minimum to $(\textbf{P}_{\tau})$ is, therefore,
$(A-L^{*}_{\tau}C,L^{*}_{\tau})$, which is the time-varying filter at
$t=\tau$. By induction, we have proved Theorem 4.1.