Year: 2021

Abstract:
This paper is concerned with representing and learning the optimal control law for the linear quadratic Gaussian (LQG) optimal control problem. In recent years, there is a growing interest in re-visiting this classical problem, in part due to the successes of reinforcement learning (RL). The main question of this body of research (and also of our paper) is to approximate the optimal control law without explicitly solving the Riccati equation. For this purpose, a novel simulation-based algorithm, namely an ensemble Kalman filter (EnKF), is introduced in this paper. The algorithm is used to obtain formulae for optimal control, expressed entirely in terms of the EnKF particles. For the general partially observed LQG problem, the proposed EnKF is combined with a standard EnKF (for the estimation problem) to obtain the optimal control input based on the use of the separation principle. The theoretical results and algorithms are illustrated with numerical experiments.

heading: 1 INTRODUCTION
text: This paper is concerned with the problem of reinforcement learning (RL) in continuous-time and continuous (Euclidean) state-space settings. A special case is the linear quadratic Gaussian (LQG) problem where the dynamic model is a linear system, the cost terms are quadratic, and the distributions – of the random initial condition and the noise – are Gaussian.
The LQG problem has a rich and storied history in modern control theory going back to the very origin of the subject (Kalman, 1960). To obtain the optimal control, the bottleneck is to solve the Riccati equation – the differential Riccati equation (DRE) in finite-horizon settings or the algebraic Riccati equation (ARE) in the infinite-horizon settings of the problem. There is a body of literature devoted to the study of these equations (Bittanti et al., 2012; Lancaster & Rodman, 1995) with specialized numerical techniques to compute the solution (Laub, 1991; Benner & Bujanović, 2016).
There are two issues which makes the LQG and related problems a topic of recent research interest: (i) In high-dimensions, the matrix-valued nature of the DRE or ARE means that any algorithm is O(d2) in the dimension d of the state-space; and (ii) the model parameters may not be explicitly available to write down the DRE (or the ARE) let alone solve it. The latter is a concern, e.g., when the model exists only in the form of a black-box numerical simulator.
These two issues have motivated the recent research on the infinite-horizon linear quadratic regulator (LQR) problem (Fazel et al., 2018; Mohammadi et al., 2021b; Tu & Recht, 2019). In LQR, the bottleneck is to solve an ARE. The algorithms studied in the recent papers seek to bypass solving an ARE by directly searching over the space of stabilizing gain matrices. Global convergence rate estimates are established for both discrete-time (Fazel et al., 2018; Dean et al., 2020; Malik et al., 2020; Mohammadi et al., 2021a) and continuous-time (Mohammadi et al., 2020a;b; 2019; 2021b) settings of the LQR problem. Extensions to the H∞ regularized LQR (Zhang et al., 2020) and Markov jump linear systems (Jansch-Porto et al., 2020) have also been carried out.
The motivation and goals of the present work are related to these recent papers, albeit the proposed solution approach is very different. Our inspiration comes from data assimilation (Reich & Cotter, 2015), where in practical applications, e.g., in weather prediction, (i) only simulation-based models are available, and (ii) these are very high-dimensional. The ensemble Kalman filter (EnKF) is an efficient simulation-based algorithm to assimilate sensor data in these applications without the need to explicitly solve a DRE (Evensen, 1994; 2006; Houtekamer & Mitchell, 2001).
Our contribution: A novel extension of the EnKF algorithm is proposed to learn the optimal control law for the general (stochastic and partially observed) setting of the LQG optimal control problem. The proposed algorithm is simulation-based (using particles) and, in particular, avoids the need to solve the DRE. Assuming full-state feedback (e.g, setting of the LQR problem), the optimal control law is directly obtained. For the partially observed case, the proposed EnKF is combined with a standard EnKF (for the estimation problem) to obtain the optimal control input, based on the use of the separation principle.
The algorithm is shown to be exact in its mean-field (when the number of particles N = ∞) limit (Prop. 2). For the finite-N approximation, an error bound is obtained (Eq. (14)). An extensive discussion is included to provide an intuitive explanation of the algorithm, situate the algorithm in the RL landscape (see discussion after Eq. (12) in Sec. 2.2), and to compare and contrast the algorithm with the state-of-the-art (Sec. 2.5). For this purpose, numerical comparisons with competing algorithms are also included (Fig. 2 (c)).
The outline of the remainder of this paper is as follows. The LQG optimal control problem and its simulation-based solution is described in Sec. 2. The algorithms are illustrated with numerical examples in Sec. 3 where comparisons with the state-of-the-art algorithms are also described.

heading: 2 LINEAR QUADRATIC GAUSSIAN (LQG) PROBLEM
text: Problem: The partially observed linear Gaussian model is expressed using the Itô-stochastic differential equations (SDE) as
dXt = AXt dt+BUt dt+ dξt, X0 ∼ N (m0,Σ0) (1a) dZt = HXt dt+ dζt, Z0 = 0 (1b)
where X = {Xt ∈ Rd : t ≥ 0} is the hidden state process, U = {Ut ∈ Rm : t ≥ 0} is the control input, and Z = {Zt ∈ Rp : t ≥ 0} is the observation process. The model parameters A,B,H are matrices of appropriate dimensions, the process noise ξ = {ξt ∈ Rd : t ≥ 0} and the observation noise ζ = {ζt ∈ Rp : t ≥ 0} are Wiener processes (w.p.) with covariance Q and R, respectively. It is assumed that X0, ξ and ζ are mutually independent andR 0. The LQG optimal control objective is to minimize
J(U) = E (∫ T 0 1 2 |CXt|2 + 12U>t RUt dt+ 12X>T PTXT ) (2)
It is assumed that (A,B) is controllable, (A,C) and (A,H) are observable, and the matrices PT 0 and R 0. In the fully observed settings, Ut is allowed to be a function of Xt. In the partially observed settings, Ut is a function of the past observations {Zs : 0 ≤ s ≤ t}. For this purpose, it is convenient to denote Zt := σ({Zs : 0 ≤ s ≤ t}) as the sigma-field generated by the observations, and consider control inputs U which are adapted to the filtration Z := {Zt : t ≥ 0}.
Classical solution: Using the seperation principle1 the LQG controller is obtained in three steps:
Step 1. Filter design: The objective of the filter design step is to compute the causal estimate X̂t := E(Xt|Zt) for t ≥ 0. The evolution equation for {X̂t : t ≥ 0} is the Kalman-Bucy filter. Notably, the optimal gain matrix for the filter is obtained by solving a forward (in time) DRE.
Step 2. Control design: The objective of the control design step is to compute the optimal feedback control law {ut(x) : 0 ≤ t ≤ T, x ∈ Rd} for the fully observed LQG problem. It is well known to be of the linear form ut(x) = Ktx, 0 ≤ t ≤ T where the optimal gain matrix Kt is obtained by solving a backward (in time) DRE.
1The separation principle hinges on the assumption that the control input Ut does not change the observation sigma-field Zt. This is valid under very mild assumptions on the control policy, e.g., Lipschitz with respect to estimate X̂t (Van Handel, 2007, Sec. 7.3), (Georgiou & Lindquist, 2013).
Step 3. Certainty equivalence: The optimal control input is obtained by combining the results from steps 1 and 2:
Ut = ut(X̂t) = KtX̂t, 0 ≤ t ≤ T
The bottleneck is to solve the DREs – the forward DRE for the optimal filter gain and the backward DRE for the optimal control gain2. In the following three sections, we describe a simulation-based algorithm for these three steps which avoids the need to explicitly solve the DREs.

heading: 2.1 STEP 1. FILTER DESIGN USING ENKF
text: The filter design objective is to compute the causal estimate X̂t = E(Xt|Zt). In the linear Gaussian settings, the conditional distribution of the Xt is Gaussian whose mean and variance are denoted as mt and Σt, respectively. These evolve according to the Kalman-Bucy filter:
dmt = Amt dt+BUt dt+ Lt( dZt −Hmt dt), m0 = E(X0) (3a) d dt Σt = AΣt + ΣtA > +Q− ΣtH>R−1HΣt, Σ0 = var(X0) (3b)
where Lt = ΣtH>R−1 is the Kalman gain. Note that (3b) is a forward (in time) DRE and its solution Σt is used to compute the optimal Kalman gain Lt.
The EnKF is a simulation-based algorithm to approximate the Kalman filter, that does not require an explicit solution of the DRE (3b). The design of an EnKF proceeds in two steps:
1. Construct a stochastic process, denoted by X̄ := {X̄t ∈ Rd : t ≥ 0}, such that the conditional distribution (given Zt) of X̄t is equal to the conditional distribution of Xt;
2. Simulate N stochastic processes, denoted by {Xit : t ≥ 0, 1 ≤ i ≤ N}, to empirically approximate the distribution of X̄t.
The process X̄ is referred to as the mean-field process. The N processes in the step 2 are referred to as particles. The construction ensures that the EnKF is exact in the mean-field (N =∞) limit. That is, for any bounded and continuous function f ,
E[f(Xt)|Zt] Step 1 = E[f(X̄t)|Zt]︸ ︷︷ ︸
exactness condition
Step 2≈ 1 N N∑ i=1 f(Xit)
The details of the two steps are as follows:
Mean-field process: The mean-field process is constructed as
dX̄t = AX̄t dt+BUt dt+ dξ̄t + L̄t( dZt − HX̄t +Hm̄t
2 dt), X̄0 ∼ N (m0,Σ0) (4)
where ξ̄ := {ξ̄t : t ≥ 0} is an independent copy of the process noise ξ, L̄t := Σ̄tH>R−1 is the Kalman gain and
m̄t := E[X̄t|Zt], Σ̄t := E[(X̄t − m̄t)(X̄t − m̄t)>|Zt] are the conditional mean and the conditional covariance, respectively, of X̄t. The right-hand side of (4) depends upon both the process (X̄t) as well as the statistics of the process (m̄t, Σ̄t). Such an SDE is an example of a McKean-Vlasov SDE. The proof of the following proposition is included in the Appendix A (see also Taghvaei & Mehta (2020, Theorem 1)).
Proposition 1 (Exactness of EnKF). Consider the mean-field EnKF (4) initialized with a Gaussian initial condition X̄0 ∼ N (m0,Σ0). Suppose the control input U is a Z-adapted stochastic process. Then its solution X̄t is a Gaussian random variable whose conditional mean and variance
m̄t = mt, Σ̄t = Σt, a.s., t > 0
evolve the same as the Kalman filter (3). 2In the steady-state or infinite horizon settings (as T → ∞), one may replace the optimal filter gain and the optimal control gain by their steady-state values. These are directly obtained by solving the respective AREs.
Finite-N approximation: The mean-field process is simulated as an interacting particle system:
dXit = AX i t dt+BUt dt+ dξ i t︸ ︷︷ ︸
i-th copy of model (1a)
+L (N) t ( dZt −
HXit +HX̂ (N) t
2 dt)︸ ︷︷ ︸
data assimilation step
, Xi0 i.i.d∼ N (m0,Σ0) (5)
where X̂(N)t := 1 N ∑N i=1X i t is the empirical mean and
L (N) t =
1 N − 1 N∑ i=1 (Xit)(HX i t −HX̂(N)t )>R−1 (6)
is the empirical approximation of the optimal Kalman gain matrix. The system (5) is referred to as the square root form of the EnKF (Bergemann & Reich, 2012, Eq (3.3)). Note that the gain is approximated entirely in terms of particles without solving the DRE (3b).
The EnKF (5) is an example of a simulation-based algorithm in the sense that N copies of the model (1a) are simulated in parallel. The simulations are coupled through a term which is referred to as the data assimilation step. This term has a gain times error feedback control structure.

heading: 2.2 STEP 2. CONTROL DESIGN USING ENKF
text: This section contains the main contribution of this paper. Assuming full-state feedback, the optimal control law for the LQG problem (1a)-(2) is
ut(x) = Ktx where Kt = −R−1B>Pt is the optimal gain matrix and Pt is a solution of the backward (in time) DRE
− d dt Pt = A >Pt + PtA+ C >C − PtBR−1B>Pt, PT (given) (7)
It is known that Pt 0 for 0 ≤ t ≤ T whenever PT 0 (Brockett, 2015, Sec. 24). Therefore, St = P −1 t is well-defined. It is straightforward to verify that St also solves a backward DRE
d dt St = ASt + StA > −BR−1B> + StC>CSt, ST = P−1T (8)
The bottleneck is to solve the DRE (7). In the following, an EnKF is proposed to obtain a simulationbased approximation of the optimal control law ut. As before, the construction proceeds in two steps: (i) definition of an exact mean-field process and (ii) its finite-N approximation.
Mean-field process: The objective is to construct a stochastic process, denoted Ȳt ∈ Rd at time t, whose variance equals Pt, the solution of the DRE (7). This is done by constructing Ȳ = {Ȳt ∈ Rd : 0 ≤ t ≤ T} as a solution of the following backward (in time) McKean-Vlasov SDE:
dȲt = AȲt dt+B d η̄t +
1 2 S̄tC >(CȲt + Cn̄t) dt, ȲT ∼ N (0, ST ) (9)
and defining Ȳt := S̄ −1 t (Ȳt − n̄t)
where η̄ = {η̄t ∈ Rm : t ≥ 0} is a w.p. with covariance matrix R−1, d η̄ in (9) denotes the backward Itô-integral (Nualart & Pardoux, 1988, Sec. 4.2), and
n̄t := E[Ȳt], S̄t := E[(Ȳt − n̄t)(Ȳt − n̄t)>] (10)
The proof of the following proposition is included in Appendix B.
Proposition 2. Consider the mean-field EnKF (9) initialized with a Gaussian initial condition ȲT ∼ N (0, ST ). Then its solution Ȳt is a Gaussian random variable whose mean and variance
n̄t = 0, S̄t = St, 0 ≤ t ≤ T
Consequently, Ȳt is also a Gaussian random variable with mean and variance
E(Ȳt) = 0, E(ȲtȲ > t ) = Pt, 0 ≤ t ≤ T
and therefore (i) If the matrix B is explicitly known then the optimal gain matrix
Kt = −R−1B>E(ȲtȲ >t ) or else (when it is not) then (ii) define the Hamiltonian3 (or the Q-function) (Mehta & Meyn, 2009)
H(x, a, t) := 12 |Cx|2 + 12a>Ra+ x>E(ȲtȲ >t )(Ax+Ba) from which the optimal control law is obtained as ut(x) = arg mina∈Rm H(x, a, t).
Finite-N approximation: The mean-field process is empirically approximated by simulating a system of interacting particles {Yit ∈ Rd : 0 ≤ t ≤ T, i = 1, . . . , N} according to
dYit = AYit dt+B d η i t︸ ︷︷ ︸ i-th copy of model (1a) +S (N) t C
>( CYit + Cn(N)t
2 ) dt︸ ︷︷ ︸
RL step
, YiT i.i.d∼ N (0, P−1T ) (11)
Y it = (S (N) t ) −1(Yit − n(N)t ) where ηi is a copy of η̄, and
n (N) t = N −1 N∑ i=1 Yit , S(N)t = 1 N − 1 N∑ i=1 (Yit − n(N)t )(Yit − n(N)t )>. (12)
Relationship to RL: The following intuitive explanation is included to situate the proposed simulation-based algorithm in the RL landscape:
Representation: In designing any RL algorithm, the first issue is representation – how does one represent the unknown value function (Pt in the linear case)? Our novel idea – the first key innovation of this paper – is to represent Pt is in terms of statistics (variance) of the particles. Such a representation is fundamentally distinct from representing the value function, or its proxies, such as the Q function, in terms of a set of basis functions (Bradtke et al., 1994; Devraj et al., 2020; Maei et al., 2010; Fujimoto et al., 2018; Melo et al., 2008).
The algorithm is entirely simulation based: N copies of the model (1a) are simulated in parallel where the terms on the righthand-side of (11) have the following intuitive explanation:
Dynamics: The first term “AYit dt” on the right-hand side of (11) is simply a copy of uncontrolled dynamics in the model (1a).
Control: The second term is the control input “BUt dt” implemented as “B dηit”. That is, the control input U for the i-th particle is a white noise process with covariance R−1. One may interpret this as an exploration step whereby the cheaper control directions are explored more.
In summary, for the i-th particle, the dynamics and control are the same as any RL algorithm (white noise is used for exploration). The difference arises due to the third term (which is the second key innovation of this paper).
RL step: The third term indicated as the RL step engenders a particle flow that effectively implements the value iteration step of the RL. There are several points to be made:
1. The RL step is a function of the state cost term in (2). This is most easily seen by writing the equation for the empirical mean
dn (N) t = (A+ S (N) t C >C)n (N) t dt+B d η
(N) t , n (N) T = 0
where η(N)t = N −1∑ i η i t. Noting that the state cost term is x
>C>Cx, the RL step implements a gradient with S(N)t C
> as the gain matrix. 3The Hamiltonian H(x, a, t) is in the form of an oracle because (Ax + Ba) is the right-hand side of the
simulation model (1a).
2. The RL step has a linear feedback control structure and serves to couple the particles. Without the RL step, the particles are independent of each other.
3. This form of the RL step is possible only if one has access to a simulator and an ability to add additional terms outside the control channel (same as the data assimilation step in estimation). For example, this is not possible when the system exists only as an experiment.
Arrow of simulation time: The particles are simulated backward – from terminal time t = T to initial time t = 0. This is consistent with the dynamic programming (DP) equation which also proceeds backward in time. However, because the focus of the RL is on the infinite-horizon problem, this important property of DP is not considered in algorithm design. It is noted that the infinite-horizon problem is easily handled in our formulation by taking T as suitably large (which is how the infinitehorizon limit is defined in DP).
It remains to use the finite-N system to approximate the optimal control law.
Optimal control: There are two cases to consider: (i) If the matrix B is explicitly known then4
K (N) t = −
1 N − 1 N∑ i=1 R−1(B>Y it )(Y i t ) > (13)
or else (when it is not) then (ii) approximate the Hamiltonian as
H(N)(x, a, t) := 12 |Cx|2 + 12a>Ra+ 1 N − 1 N∑ i=1 (x>Y it )(Y i t ) > (Ax+Ba)︸ ︷︷ ︸
model (1a)
from which the optimal control law is obtained as
u (N) t (x) = arg min a∈Rm H(N)(x, a, t)
There are several zeroth-order approaches to solve the minimization problem, e.g., by constructing 2-point estimators for the gradient. Since the objective function is quadratic and the matrix R is known, m queries of H(N)(x, ·, t) are sufficient to compute u(N)t (x).

heading: 2.3 ENKF-BASED ALGORITHM FOR THE LQG CONTROLLER
text: By combining the result of steps 1 and 2, the optimal control input
Ut = K (N) t X̂ (N) t if B is known
u (N) t (X̂ (N) t ) o.w.
The overall algorithm including steps 1, 2 and 3 is tabulated in the Appendix C where additional remarks on numerical approximation of backward SDEs also appear.

heading: 2.4 CONVERGENCE AND ERROR ANALYSIS
text: The mean-field process (9) represents the mean-field limit of the finite-N system (11), as the number of particles N → ∞. The convergence analysis is a delicate matter based on the propagation of chaos (Bishop & Del Moral, 2018). In Appendix D, under certain additional assumptions on system matrices, the following error bound is derived:
E[‖S(N)t − S̄t‖2F ] ≤ 3‖S̄T ‖2F N e−(4µ− 1 N )(T−t) + cvar N , 0 ≤ t ≤ T (14)
where µ and cvar are positive constants and ‖ · ‖F is the Frobenius (matrix) norm. The computational complexity of an ENKF is O(Nd). EnKF is a workhorse in applications such as weather prediction where models are simulation-based (Evensen, 2006; Reich & Cotter, 2015). In these applications, the number of simulations N << d.
4Eq. (13) for the optimal control gain matrix K(N)t is dual to the eq. (6) for optimal filter gain matrix L (N) t .

heading: 2.5 COMPARISON WITH RELATED WORK
text: Classical RL algorithms for the LQR problem are based on a linear function approximation, using quadratic basis functions, of the value function or its surrogate the Q-function (Bradtke et al., 1994). Convergence guarantees typically require (i) a persistence of excitation condition, and (ii) use of the on policy methods whereby the parameters are learned for a given fixed policy (which is subsequently improved). For the deterministic LQR problem, the persistence of excitation condition is difficult to justify theoretically, and in practice can lead to poor performance related to slow convergence rates. These limitations have spurred recent research on the LQR problem.
Our algorithm is compared to the model-free policy optimization based methods of Mohammadi et al. (2021a) and Fazel et al. (2018). In policy optimization, the value J(U) is minimized over the search space of stabilizing gain matrices, using a gradient-descent approach, starting from a stabilizing controller (which is not straightforward to obtain in a model-free setting). For each perturbation of the gain matrix, the gradient is estimated by simulating N particles (trajectories) over a time-horizon [0, T ]. Multiple iterations are needed to converge to the minimizer.
In contrast, using the EnKF, the matrix Pt is approximated by simulating N particles once over the time-horizon [0, T ]. The optimal control is obtained from an application of the minimum principle that only requiresm evaluations at each t (the minimization of Hamiltonian is an online calculation).
The trade-off between the two approaches is as follows: While policy optimization methods require multiple iterations with a small number of particles, the EnKF requires only a single iteration with relatively larger number of particles. As illustrated with the aid of numerical examples in Sec. 3, this can lead to an order of magnitude gain in the computation time.
As a final point, our paper provides a simulation-based algorithm for the most general class of linear quadratic problem: fully or partially observed, finite or infinite-horizon, stochastic or deterministic, and furthermore does not require an initial stabilizing controller.
Duality: In Appendix E, it is shown that the Gaussian density of the random variable Ȳt equals a log transformation of the value function (at time t) of the LQG optimal control problem. Additional remarks are also included to expound the duality roots of the proposed algorithm, including a survey of the relevant literature on this beautiful subject.

heading: 3 NUMERICAL SIMULATIONS
text: In all the following three numerical examples, we consider the infinite-horizon fully observed deterministic LQR problem. Although our algorithms are more generally applied, such a choice allows us to compare and contrast with the recent RL work on the LQR problem. For the LQR problem, the solution of the ARE is denoted P∞ and the associated optimal feedback gain is denoted K∞.
In each of the three numerical examples, a time-horizon T is fixed. Over the time-horizon, the finiteN EnKF algorithm (11) is run to obtain an empirical approximation {P (N)t ∈ Rd×d : 0 ≤ t ≤ T}. For the sake of comparison, the exact {Pt ∈ Rd×d : 0 ≤ t ≤ T} is obtained by numerically solving the backward DRE (7). Typically, we chose PT = I , the identity matrix.

heading: 3.1 EXPONENTIAL CONVERGENCE
text: An attractive feature of the proposed EnKF algorithm is that the mean-field limit is exact. This means the EnKF recovers the properties of the DRE in the limit. One such property is the exponential convergence: For any fixed t, Pt → P∞ as T → ∞, starting from any initialization PT (Ocone & Pardoux, 1996, Remark 2.1). Moreover, because of the error formula (14), P (N)t → P∞ with a small additional bias that decays as O( 1N ).
For a d = 2 dimensional system, Fig. 1 depicts this exponential convergence of the four entries of the symmetric P (N)t matrix (using N = 100 particles). As part of Appendix F.1, additional figures are included to depict exponential convergence of the 100 entries of the 10×10 matrix P (N)t for a d = 10 dimensional system (using N = 1000 particles). In these evaluations, the system (A,B) is chosen to be in its controllable canonical form where the appropriate entries (last row) of the A matrix are randomly generated. Because these are randomly generated, the matrix A is
t shows exponential convergence, with
a small (random) error, to the ARE limit P∞ which is also depicted together with the exact DRE solution. (Note the x-axis of the plots (b)-(d) is T − t so convergence is easy to see).
often unstable (see also Fig. 1 (a)) which does not pose any problem with our algorithm. Additional details on the numerics appear as part of Appendix F.1.

heading: 3.2 EVALUATION AND COMPARISON WITH LITERATURE: MASS SPRING DAMPER SYSTEM
text: In order to evaluate and compare the performance of the proposed algorithm, we describe next the results for the coupled mass-spring system studied in Mohammadi et al. (2019). For each mass, there are two states, position and velocity. A system with d/2 masses is d-dimensional. To evaluate the performance of the finite-N algorithm, we consider the following error metric:
MSE := 1
T E (∫ T 0 ‖Pt − P (N)t ‖2F ‖Pt‖2F dt ) The expectation is approximated empirically by averaging over 100 simulation runs. See Appendix F.2 for details on modeling, parameter values, and the numerical discretization. Fig. 2 depicts the results of the numerical experiments showing the O( 1N ) decay of MSE as N increases (for d fixed).
Part (c) of the figure depicts a comparison with the two state-of-the-art algorithms, namely Mohammadi et al. (2021b) and Fazel et al. (2018). In the figure, computation time is plotted against the relative error in approximating the LQR gain K∞ (plots are qualitatively similar using other error metrics). The computation time is obtained using the Python process time() function from the time library for measuring execution time. The result depicted in the figure is for d = 10 but qualitatively similar results are also obtained for other values of d. See Appendix F.3 for details on numerical implementation of the algorithms. We also implemented the classical RL algorithm in Bradtke et al. (1994) but its convergence was not reliable because of the persistence of excitation issues.
Based on the preliminary analysis, the main reason for the order of magnitude improvement in computational time is as follows: An EnKF requires only a single iteration over a fixed time-horizon [0, T ]. Although the number of particles (N ) is an order of magnitude larger for the EnKF algorithm, certain vectorisation features of numpy yield significant gains in computational time. Because the other algorithms require multiple iterations which must necessarily be carried out serially, these computations are slower. In our comparisons, we used the same time-horizon T and discretization time-step δt for all the algorithms. It is certainly possible that some of these parameters can be optimized to improve the performance of the other algorithms. In particular, one may consider shorter or longer time-horizon T or use parallelization to speed up the gradient calculation.

heading: 3.3 EVALUATION OF ENKF FOR A NONLINEAR CART-POLE SYSTEM
text: It is a time-honored practice to design the optimal control law for a (linearized) LQR model and then implement this control law on a nonlinear system. By the principle of linearization, the control law works well (in theory) for small perturbations and often (in practice) for even large ones. We consider the nonlinear conservative cart pole model (Tedrake; Rawlik et al., 2013). The control acts as external force applied to the cart. The four-dimensional state for the system is (θ, x, ω, v), where θ ∈ S1 (the circle) is the angle of the pole (pendulum) as measured from the stable equilibrium, x ∈ R is the displacement of cart along the horizontal, and (ω, v) := (θ̇, ẋ) ∈ R2 is the velocity vector. The control objective is to balance the pole – stabilize the system at the inverted equilibrium (π, 0, 0, 0), assuming full state feedback. For the purposes of control design, the nonlinear system is first linearized at the desired equilibrium and an LQR problem is formulated (see Appendix F.4 for details). For the purposes of evaluation, the optimal control is applied to the nonlinear model. Figure 3 depict the results of numerical experiments for three different choices of N . It is seen that as few as N = 10 particles are sufficient to stabilize the equilibrium. Using N = 1000 particles, the closed-loop trajectories are virtually indistinguishable from the DRE-based solution.

heading: 4 CONCLUSIONS
text: In recent years, there has been a concerted effort to revisit the LQR problem in the context of RL (Fazel et al., 2018; Mohammadi et al., 2021b; Tu & Recht, 2019; Dean et al., 2020; Malik et al., 2020; Zhang et al., 2020; Jansch-Porto et al., 2020). The effort is largely in response to the well known issues that arise when the state and action space are continuous (Lu et al., 2021).
In this paper, we present a new paradigm for RL. There are two key innovations: (i) the representation of the unknown value function in terms of the statistics (variance) of the particles; and (ii) design of interactions between simulations to solve the optimal control problem. For the LQR problem, this is shown to yield a learning rate that closely approximates the exponential rate of convergence of the solution of the DRE. In numerical examples, this property is shown to lead to an order of magnitude better performance than the state-of-the-art algorithms. Given the enormous success of EnKF in data assimilation (Evensen, 2006; Reich & Cotter, 2015), the contributions of this paper potentially open up new opportunities for RL. It is our hope that the paper will engender new synergies between the data assimilation and the RL communities.

heading: 5 REPRODUCIBILITY STATEMENT
text: The implementable pseudo-code of the algorithm is tabulated in Appendix C. For the three examples in Sec. 3, the numerical parameters are tabulated in Appendix F. The codes for all numerical simulations are provided as part of the supplementary material. All the codes have a README file associated with it. Table 1 provides a list of the Appendix containing the simulation parameters and location of codes for each numerical example.

heading: A PROOF OF PROP. 1
text: It is first shown that the conditional mean m̄t and the conditional covariance Σ̄t evolve according to the Kalman filter equations (3). The SDE for the mean is obtained by taking the conditional expectation of (4), dm̄t = Am̄t dt+BUt dt+ L̄t( dZt −Hm̄t dt) where we used the fact that Ut is assumed to be Zt measurable. To obtain the equation for the covariance, define the error process et = X̄t − m̄t which evolves according to
det = (A− 1
2 L̄tH)et dt+ dξ̄t
Then, upon an application of the Itô rule
d(ete > t ) = (A−
1 2 L̄tH)(ete > t ) dt+ (ete > t )(A− 1 2 L̄tH) > dt+Qdt+ dξ̄te > t + et dξ̄ > t
and taking the conditional expectation, the equation for the conditional covariance Σ̄t = E[ete>t |Zt] is obtained as
d dt Σ̄t = AΣ̄t + Σ̄tA > +Q− Σ̄tH>HΣ̄t where we used the definition L̄t = Σ̄tH>.
The equation for the conditional covariance Σ̄t is identical to the DRE (3b). Therefore, if the initial condition Σ̄0 = Σ0 then Σ̄t = Σt for all t > 0. This also implies L̄t = Lt, which in turn implies that the SDE for the conditional mean m̄t is identical to Kalman filter equation for mt. If m̄0 = m0 then m̄t = mt for all t > 0.
By replacing the mean-field terms m̄t and Σ̄t by exogenous processes mt and Σt, the McKeanVlasov SDE (4) simplifies to an Ornstein-Uhlenbeck SDE. Because the distribution of the initial condition X̄0 is Gaussian, the distribution of X̄t is also Gaussian and given by N (mt,Σt).

heading: B PROOF OF PROP. 2
text: The proof is similar to the proof of Prop. 1. The equation for the mean n̄t is obtained by taking the expectation of SDE (9),
dn̄t = (A+ S̄tC >C)n̄t dt
Because n̄T = 0, we have n̄t = 0 for all t ∈ [0, T ]. The equation for the covariance S̄t is obtained by writing the SDE for the error et := Ȳt − n̄t:
det = (A+ 1
2 S̄tC
>C)et dt+B d η̄t,
Using the Itô rule for ete>t ,
d(ete > t ) = (A+
1 2 S̄tC >C)(ete > t ) dt+(ete > t )(A+ 1 2 S̄tC >C)>−BR−1B>+B d η̄te>t +et(B d η̄t) >
The Itô correction term appears with a negative sign because the SDE involves a backward Wiener process
η̄t (Nualart & Pardoux, 1988, Sec. 4.2). Taking an expectation yields the following equation
for S̄t: d
dt S̄t = (A+
1 2 S̄tC >C)S̄t + S̄t(A+ 1 2 S̄tC >C)> −BR−1B>
The SDE is identical to the SDE for St. Because S̄T = ST , we have S̄t = St for all t ∈ [0, T ]. The conclusion that Ȳt is Gaussian follows from the fact that with n̄t = nt and S̄t = St, the SDE for Ȳt is an Ornstein-Uhlenbeck SDE with a Gaussian terminal condition.
The proof for the rest of proposition is straightforward. By definition,
E[Ȳt] = E[S̄ −1 t (Ȳt − n̄t)] = S̄−1t (n̄t − n̄t) = 0
E[ȲtȲ > t ] = E[S̄ −1 t (Ȳt − n̄t)(Ȳt − n̄t)>S̄−1t ] = S̄−1t = S−1t = Pt
Since E[ȲtȲ >t ] = Pt, the optimal gain matrix Kt = −R−1B>E[ȲtȲ >t ]. For a given x ∈ Rd, the optimal control ut(x) = Ktx = −R−1B>E[ȲtȲ >t ]x is the unique minimizer of the Hamiltonian. This is referred to as the minimum principle of optimal control.

heading: C DETAILS OF THE ALGORITHMS TO SOLVE THE PARTIALLY OBSERVED LQG PROBLEM IN SEC. 2
text: In this section, the implementation details of the EnKF-based algorithms are presented to numerically solve the partially observed LQG problem (1)-(2) introduced in the main body of the paper. For better readability, the overall algorithm is broken down into three separate algorithms:
1. Algorithm 1 is an offline algorithm. It is based on the finite-N approximation of the (control) EnKF as described in Sec. 2.2 of the main body of the paper. The algorithm is run offline to obtain an approximation of the {Pt : 0 ≤ t ≤ T}.
2. Algorithm 2 is an online algorithm. It is based on the finite-N approximation of the (filter) EnKF as described in Sec. 2.1 of the main body of the paper. The algorithm is run online in a real-time manner. At each time step, it processes the sensor measurements from the true system (plant) and computes the optimal control input by calling Algorithm 3.
3. Algorithm 3 computes the optimal control input. It is based upon minimization of the Hamiltonian function. In an online implementation, the optimal control input is applied to the plant at each time step.
The input structure of each of the three algorithms is clearly delineated. In particular, the algorithms require only the simulator in the form of the function evaluator f(x, u) = Ax+Bu for the dynamics, c(x) = Cx for the cost function, and h(x) = Hx for the observation function. The algorithms do not require solution of the DRE. Remark 1. In a numerical implementation of step 1 and step 2, there are two sources of error: (i) error on account of finite-N approximation; and (ii) error on account of time discretization which depends upon the step size ∆t. The finite-N approximation error has been studied in the EnKF literature where it is shown that the error in approximating ‖Pt − P (N)t ‖ (in step 1) and the error in approximating the gain the mean ‖X̂t − X̂(N)t ‖ (in step 2) converges to zero with the rate O( 1√
N ) Del Moral & Tugaut (2018); Bishop et al. (2019); Bishop & Del Moral (2018); Taghvaei & Mehta (2018); Taghvaei (2019). The time discretization error is also expected to be bounded and converges to zero with the rate O(∆t) when the system is stable Kloeden & Platen (1999).
Algorithm 1 [offline] EnKF algorithm to approximate {Pt : 0 ≤ t ≤ T} Input: Simulation time T , simulation step-size ∆t, number of particles N , simulator f(x, u) =
Ax+Bu, terminal cost PT , cost function c(x) = Cx, and control cost matrix R. Output: {P (N)k : k = 0, 1, 2, . . . , T∆t}.
1: TF = T ∆t 2: P (N)TF = PT 3: Initialize {YiTF }Ni=1 i.i.d∼ N (0, P−1T )
4: calculate n(N)TF = N −1∑N
i=1 YiTF 5: for k = TF to 1 do 6: Calculate ĉ(N)k = N −1∑N i=1 c(Yik)
7: Calculate M (N)k = (N − 1)−1 ∑N i=1(Yik − n (N) k )(c(Yik)− ĉ (N) k ) > 8: for i = 1 : N do 9: ∆ηik
i.i.d∼ N (0, 1∆tR−1) 10: ∆Yik = f(Yik,∆ηik)∆t+ 12M (N) k (c(Yik) + ĉ (N) k )∆t 11: Yik−1 = Yik −∆Yik 12: end for 13: Calculate n(N)k−1 = N −1∑N i=1 Yik−1
14: Calculate S(N)k−1 = (N − 1)−1 ∑N i=1(Yik−1 − n (N) k−1)(Yik−1 − n (N) k−1) > 15: P (N)k−1 = (S (N) k−1) −1 16: end for
Algorithm 2 [online] EnKF algorithm to approximate state estimate X̂t and optimal control ut Input: Simulation time T , simulation step-size ∆t, number of particles N , simulator f(x, u) =
Ax + Bu, initial distribution N (m0,Σ0), process noise covariance R, observation function h(x) = Hx, {P (N)k : k = 0, 1, 2, . . . , T∆t} from the offline algorithm F.1. Output: estimate {X̂(N)k : k = 0, 1, 2, . . . , T∆t} and optimal control input {u (N) k : k =
0, 1, 2, . . . , T∆t − 1}. 1: Define TF := T∆t 2: Initialize particles {Xi0}Ni=1
i.i.d∼ N (m0,Σ0) 3: for k = 0 to TF − 1 do 4: Calculate X̂(N)k = 1 N ∑N i=1X i k 5: Calculate u(N)k = arg minaH(X̂ (N) k , P (N) k X̂ (N) k , a) from algorithm 3 6: Apply control u(N)k to the true system and obtain the observation ∆Zk = Z(k+1)∆t − Zk∆t 7: Calculate ĥ(N)k = 1 N ∑N i=1 h(X i k) 8: Calculate L(N)k = 1 N−1 ∑N i=1(X i k − X̂ (N) k )(h(X i k)− ĥ (N) k )
>R−1 9: for i = 1 to N do

heading: 10: ∆ξik
text: i.i.d∼ N (0, Q∆t)
11: ∆Xik = f(X i k, u (N) k )∆t+ ∆ξ i k + L (N) k (∆Zk −
h(Xik)+ĥ (N) k
2 ∆t)
12: Xik+1 = X i k + ∆X i k 13: end for 14: end for
Algorithm 3 Computation of optimal control Input: state x, momentum y, control cost matrix R,
HamiltonianH(x, y, α) = yT (a(x) + b(x)α) + 12 |c(x)|2 + 12α>Rα. Output: optimal control u = arg minαH(x, y, α).
1: if b(x) is known then 2: u = −R−1b(x)>y 3: else 4: for k = 1 to m do 5: ek = [0, . . . , 0, 1︸︷︷︸
k-th component
, 0, . . . , 0] ∈ Rm
6: uk = H(x, y,R−1ek)−H(x, y, 0)− 12 (R−1)kk 7: end for 8: end if

heading: D ERROR ANALYSIS
text: The objective of the error analysis is to study the convergence of the empirical variance of the particles S(N)t to its mean-field limit S̄t as N →∞. According to the Proposition 2, the mean-field limit S̄t = St, hence it evolves backward in time according to
d dt S̄t = AS̄t + S̄tA > + S̄tC >CS̄t −BR−1B>, S̄T = PT (15)
In order to study the error S(N)t − S̄t, we obtain a stochastic differential equation for the evolution of S(N)t .
Lemma 1. The empirical covariance matrix S(N)t , defined in (12), evolves according to the stochastic differential equation
dS (N) t =
( AS
(N) t + S (N) t A > + S (N) t C >CS (N) t −BR−1B> ) dt+ dMt, (16)
where Mt is a Martingale given by
dMt = 1 N − 1 N∑ i=1 eit(B d η i t) > +B d η i t(e i t) >, eit := Yit − n(N)t
with quadratic variation
d〈M〉t = 1 N − 1 [ Tr(BR−1B>)S(N)t +BR −1B>Tr(S(N)t ) +BR −1B>S (N) t + S (N) t BR −1B> ]
Proof. The evolution for the error process eit := Yit − n(N)t is obtained by subtracting the evolution for the empirical mean n(N)t from (11):
deit = Ae i t +B( d
η i t − 1
N N∑ j=1 B d η j t ) + S (N) t Ce i t dt
The evolution for the empirical covariance S(N)t = 1 N−1 ∑N i=1 e i t(e i t) T is obtained by application of the Itô rule.
Note that the first term in the evolution of S(N)t in (16) is exactly identical to that of S̄t in (15), while the second term is an additional stochastic term that converges to zero as N → ∞. Even though the fluctuations scale as O(N− 12 ), the analysis is challenging as has been noted in literature (see the remark after Theorem 3.1 in (Del Moral & Tugaut, 2018)). Error analysis of the (stochastic) ensemble Kalman filter has been carried out in the literature (Del Moral & Tugaut, 2018; Bishop & Del Moral, 2018) under additional assumptions which we also make below: Assumption 1. The matrix C>C is identity. And the control matrix B is full rank.
The assumption regarding the cost matrix C can be relaxed to CTC being non-singular by a change of coordinates. The assumption regarding the control matrix is strong. It is an open problem in the literature regarding stability of EnKF to carry out the error analysis under the more natural assumption that the system is controllable and observable.
Under assumption 1, we can prove the main result stated in the following proposition. Proposition 3. Let S̄t be the mean-field covariance defined in (10) and SNt be the empirical covariance of the particles defined in(16). Then, under Assumption 1, the error between S(N)t and S̄t satisfies the upper-bound
E[‖S(N)t − S̄t‖2F ] ≤ 3cµ‖PT ‖2F
N e−(4µ∞− 1 N )(T−t) + cvar N , 0 ≤ t ≤ T (17)
where µ∞, cµ, cvar are positive constants that are independent of time t.
In the remainder of this section, we prove the proposition 3. To simplify the presentation, we use the time-reversed quantitative Ω(N)t := S (N) T−t and Ωt := S̄T−t which evolve forward in time according to
dΩt = (−AΩt − ΩtA> − ΩtC>CΩt +BR−1B>)dt (18) dΩ
(N) t = (−AΩ(N)t − Ω(N)t A> − Ω(N)t C>CΩ(N)t +BR−1B>) dt+ dMt (19)
where Mt is defined in Lemma 1.
First, we show some basic results regarding the asymptotic properties of the Ricatti equation that governs Ωt. These results are straightforward application of the classical stability theory of the Ricatti equation. Lemma 2. Consider the Ricatti equation (18) that governs Ωt. Then,
1. There exists a positive definite matrix solution Ω∞ to the ARE
0 = −AΩ∞ − Ω∞A> − Ω∞C>CΩ∞ +BR−1B>
such that −A− Ω∞C>C is Hurwitz. 2. Ωt converges to Ω∞ exponentially fast
3. Let µt denote the minimum eigenvalue for the matrix Ft := A + 12C >CΩt. Then, there
exists a positive constant µ∞ > 0 and cµ such that
e− ∫ t 0 µτ dτ < cµe −µ∞t, ∀t ≥ 0 (20)
Proof. Under Assumption 1, the pair (A,BR− 12 ) is controllable and the pair (A,C) is observable. Part (1) follows from (Kwakernaak & Sivan, 1972, Theorem 3.7) and part (2) follows from (Ocone & Pardoux, 1996, Lemma 2.2 and Theorem 2.3). The proof of part (3) follows from the following two facts (i) (−A− 12Ω∞C>C) is Hurwitz because
−BR−1B> = (−A− 1 2 Ω∞C >C)Ω∞ + Ω∞(−A− 1 2 Ω∞C >C)>
is a Lyapunov equation, with BR−1B> positive definite since B is assumed to be full rank; and (ii) Ωt converges to Ω∞ exponentially fast. Therefore, there exists a time τ1 > 0 such that (−A − 1 2ΩtC
>C is Hurwitz for t ≥ τ1. As a result, the minimum eigenvalue µt of Ft = A+ 12ΩtC>C is lower-bounded µt ≥ µ∞ > 0 for t ≥ τ1. This concludes the bound (20) where cµ = e− ∫ τ1 0 µτ dτ
Next, we obtain a bound for E[Tr(Ω(N)t )]. Upon taking the expectation and trace of (19)
d dt E[tr(Ω(N)t )] = −tr(AE[Ω(N)t ])− tr(E[Ω(N)t ]A>)− tr(E[Ω(N)t C>CΩ(N)t ]) + tr(BR−1B>)
≤ −tr(AE[Ω(N)t ])− tr(E[Ω(N)t ]A>)− tr(E[Ω(N)t ]C>CE[Ω(N)t ]) + tr(BR−1B>)
where we used E[Tr(X>X)] ≥ Tr(E[X]E[X]>) forX = Ω(N)t C> to conclude the inequality. Then, by inspecting the equation for Tr(Ωt), it follows that
E[Tr(Ω(N)t )] ≤ Tr(Ωt) ≤ E0 (21) where E0 is a uniformly in time bounded constant, because Ωt converges to Ω∞.
Next, we obtain the bound for the error Ω(N)t − Ωt. Subtracting (18) from (19) yields d(Ω
(N) t − Ωt) = (−A(Ω(N)t − Ωt)− (Ω(N)t − Ωt)A> − (Ω(N)t − Ωt)C>C(Ω(N)t − Ωt)
− (Ω(N)t − Ωt)C>CΩt − ΩtC>C(Ω(N)t − Ωt))dt+ dMt Then, by application of the Itô rule to (Ω(N)t − Ωt)>(Ω(N)t − Ωt) = (Ω(N)t − Ωt)2 and taking the trace,
dTr((Ω(N)t − Ωt)2) = 2tr((Ω(N)t − Ωt)2)(−A−A> − C>CΩt − ΩtC>C))dt − 2tr(C>C(Ω(N)t − Ωt)3))dt+ tr((dMt + dM>t )(Ω(N)t − Ωt))
+ 2
N
( tr(Ω(N)t BR −1B>) + tr(Ω(N)t )tr(BR −1B>) ) dt
Upon taking the expectation and defining Θt = E[tr((Ω (N) t − Ωt)2)],
d dt Θt ≤ 2E
[ tr(−A−A> − C>CΩt)(Ω(N)t − Ωt)2) ] − 2E [ tr(C>CΩ(N)t (Ω (N) t − Ωt)2) ] + 1
N
( Θt + 2tr((BR−1B>)2) + tr(Ω2t ) + 2tr(BR −1B>)tr(Ωt))
where we used the bound (21), and Tr(XY ) ≤ 12 Tr(X>X) + 12 Tr(Y >Y ) for X = Ω (N) t − Ωt and Y = BR−1B>.
In order to continue with proving the bound, we use the assumption C>C = I to conclude d
dt Θt ≤ 2E
[ tr(−A−A> − Ωt)(Ω(N)t − Ωt)2) ] + 1
N
( Θt + 2tr((BR−1B>)2) + 2tr(Ωt)tr(BR−1B>) + tr(Ω2t ))
because tr(Ω(N)t (Ω (N) t − Ωt)2) ≥ 0. Morevoer, due to exponential convergence of Ωt to Ω∞, tr(Ω2t ) ≤ E1 uniformly. Thus, d
dt Θt ≤ 4µtΘt +
1
N (Θt + 2tr((BR−1B>)2) + 2E0tr(BR−1B>) + E1)
where µt is the minimum eigenvalue for A+ 12Ωt. Therefore,
Θt ≤ Θ0e−4 ∫ t 0 µt′dt
′+ 1N t + cvar N
≤ Θ0cµe−4µ∞t+ 1 N t + cvar N
where cvar = 2tr((BR−1B>)2) + 2E0tr(BR−1B>) + E1 and we used (20). The proof of the proposition follows by substitution S(N)T−t − St = Ω (N) t − Ωt and the bound Θ0 ≤ 3‖PT ‖FN .

heading: E DUALITY ROOTS OF THE PROPOSED ALGORITHM
text: For t ∈ (0, T ), the LQG value function is defined as follows
vt(x) := min {Us:t≤s≤T} E (∫ T t ( 1 2 |CXs|2 + 12U>s RUs ) ds+X>T PTXT ∣∣∣∣Xt = x )
(22)
From the DP optimality principle, it is easily shown Liberzon (2012) that
• The value function vt(x) = 12x>Ptx+ (constant) is quadratic where • The matrix-valued process {Pt ∈ Rd×d : 0 ≤ t ≤ T} solves the DRE with terminal
condition PT at time t = T .
Log transformation: The log transformation is a manifestation of the duality between optimal control and optimal estimation Mitter & Newton (2003); Todorov (2008). Define the probability density pt(x) ∝ e−vt(x) Fleming & Mitter (1982); Fleming (1977). For the LQG problem, vt is quadratic and therefore, pt is a Gaussian density with variance P−1t = St. Indeed, the Gaussian density of the random variable Ȳt equals pt. The following proposition gives the precise connection: Proposition 4. Suppose pt ∝ e−vt where vt is the value function defined in (22). Let p̄t be the density of the random variable Ȳt defined in (9). Then, provided pT = p̄T ,
pt(x) = p̄t(x), ∀x ∈ Rd, 0 ≤ t ≤ T
Proof. By Proposition 2, Ȳt is a Gaussian random variableN (0, St) where St = P−1t by definition. Therefore,
− log(p̄t(x)) = 1
2 x>Ptx+ (constant) = vt(x) + (constant) = − log(pt(x)) + (constant)
Therefore, since pt(x) and p̄t(x) are both normalized to one, the constant is zero and p̄t(x) = pt(x).
In the past, duality has been used to:
(i) Obtain linearly solvable sampling algorithms to solve optimal control problems Kappen & Ruiz (2016); Rawlik et al. (2013); Theodorou et al. (2010); Schütte et al. (2012); Todorov (2009); and
(ii) Set up optimal control problems for the purposes of estimation and simulation Sutter et al. (2016); Hartmann & Schütte (2012); Van Handel (2006); Kim & Mehta (2020).
(iii) A parallel approach to duality involves posing optimal control problems as problems of probabilistic inference, and using methods like expectation maximisation (Toussaint & Storkey, 2006), message passing (Hoffmann & Rostalski, 2017; Watson & Peters, 2021; Watson et al., 2021) or minimisation of KL divergence (Kappen et al., 2012; Rawlik et al., 2013; Watson et al., 2020).
Although novel and distinct from the algorithms described in these earlier works – which instead rely on the use of a Feynman-Kac representation – the proposed algorithm is an example of a sampling algorithm to solve the optimal control problem. In a related work, we have extended the results presented here to a class of nonlinear optimal control problems. These nonlinear generalization more clearly reveals the duality roots of the proposed algorithm. The EnKF algorithm arises as a special case in the LQG settings. The nonlinear generalization will be a subject of future publication.

heading: F DETAILS OF THE NUMERICAL EXAMPLES IN SEC. 3
text: Notation: In ∈ Rn×n denotes the identity matrix and 1d ∈ Rd denotes a vector with all entries equal to 1.
F.1 EXPONENTIAL CONVERGENCE
The EnKF algorithm is implemented for the linear system (1) to approximate the optimal gain matrix. The model matrices A and B are in the control canonical form
A =  0 1 0 0 . . . 0 0 0 1 0 . . . 0 ...
... a1 a2 a3 a4 . . . ad
 , B = 
0 0 ... 1  where the entries (a1, . . . , ad) ∈ Rd is selected randomly from N (1.651d, Id). Additional model and simulation parameters are summarized in Table 2 and Table 3.
The numerical result for a single realization of the algorithm, for d = 2, is depicted in Figure 1. Part-(a) depicts the spectral properties of the open-loop (ẋ = Ax) and the closed-loop (ẋ = (A + BK
(N) 0 )x) dynamics (recall that due to the backward in time nature of the EnKF, K0 is the terminal
gain which the algorithm yields). It is observed that the algorithm learns the gain matrix that make the closed-loop dynamics stable, while the open loop dynamics A is unstable.
The entry-wise convergence of the matrix P (N)t as t goes from T to 0 is depicted in part (b)-(c)-(d). The result is compared with the exact solution to the DRE Pt and the asymptomatic limit, as the time horizon T →∞, which is the solution to the ARE.
Additional numerical results, for the case where d = 10, are depicted in Fig. 4 for the open and closed-loop eigenvalues and in Fig. 5 for the convergence of the 100 entries of the 10 × 10 matrix P
(N) t .
F.2 COUPLED MASS SPRING DAMPER SYSTEM
This system is taken from Mohammadi et al. (2019). The matrices A and B are as follows:
A = [ 0ds×ds Ids −T −T ] , B = [ 0ds×ds Ids ] where ds = d2 is the number of masses and T ∈ Rds×ds is a Toeplitz matrix with 2 on the main diagonal and−1 on the first sub-diagonal and first super-diagonal. Additional model and simulation parameters are listed in Table 4 and 5, respectively.
F.3 COMPARISON BETWEEN ENKF AND POLICY-GRADIENT METHODS
The performance of the EnKF algorithm is compared with policy gradient algorithms presented in Mohammadi et al. (2021c) (denoted as [M21]) and Fazel et al. (2018) (denoted as [F18]). These two algorithms are based on gradient free optimization for the gain matrix, with the difference that [M21] is for continuous-time systems while [F18] is for discrete-time systems. We used the discretized linear system xt+∆t = (I + A∆t)xt + But∆t to implement [F18]. The two algorithms involve similar set of hyper-parameters that are selected optimally for each experiment for fair comparison.
The three algorithms, EnKF, [M21], and [F18] are implemented for the mass-spring system of Appendix F.2 for d = 2, 4, 10. The comparison is made by studying the relationship between the computational time and the error in approximating the optimal gain matrix. The error is measured with the formula
error =  ‖K(N)0 −K∞‖F ‖K∞‖F ; for EnKF
‖K(N)T −K∞‖F ‖K∞‖F ; for [M21] and [F18]
where K(N)0 is the terminal gain output by EnKF, K (N) T is the terminal gain output by the [M21] and [F18] algorithms, and K∞ is the analytical gain obtained through the ARE.
In order to study the error vs computational time relationship for EnKF algorithm, the number of particles N is varied, while for [M21] and [F18] the number of iterations is varied. The simulation time horizon T = 10, and the step-size ∆t = 0.01 is the same for all algorithms. The hyperparameters of the [M21], and [F18] algorithms are presented in Table 6. The initial guess K0 = 0, initial distribution D0 = N (0, Id), and gradient descent step α = 0.0001 for both [M21] and [F18]. The values of the other hyper parameters are in Table 7. The numerical results for d = 10 are depicted in Figure 2 and for d = 2, 4 in Figure 6. Additionally, Figure 7 shows comparison when the error is measured as
error =  c (N) 0 −c∞ c (N) T −c∞ ; for EnKF
c (N) T −c∞ c (N) 0 −c∞ ; for [M21] and [F18]
where c(N)0 is the optimal control cost corresponding to the terminal gainas per EnKF, and by the initial gain as per [M21] and [F18] algorithms; and c(N)T is the optimal control cost corresponding to the initial gain as per EnKF, and terminal gain as per the [M21] and [F18] algorithms; and c∞ is the optimal control cost achieved by analytical gain obtained through the ARE.
The simulations are implemented in Python 3 on a Intel Xeon E3-1240 V2 3.40 Ghz CPU, and the process time() function from the time module is used to evaluate the execution time.
F.4 CART-POLE SYSTEM
The nonlinear model is taken from (Tedrake, Chapter 3.2.1):
θ̇ = ω
ω̇ = 1
l(M +m sin2(θ))
( −F cos(θ)−mlω2 cos(θ) sin(θ)− (m+M)g sin(θ) ) ẋ = v
v̇ = 1
M +m sin2(θ)
( F +m sin(θ)(lω2 + g cos(θ)) )
For the LQG control design, we linearize the equations about the desired inverted equilibrium (π, 0, 0, 0). The associated A and B matrices are as follows:
A =  0 0 1 0
0 0 0 1 (M+m)g Ml 0 0 0 mg M 0 0 0
 , B =  0 1 Ml
0 1 M  The model parameters are listed in Table 8 and the simulation parameters are in Table 9.

Review summaries:
The paper is about designing a simulation-based ensemble Kalman filter algorithm for learning the optimal control policy for the Linear Quadratic Gaussian (LQG) control problem. The main contribution is that the paper extends the existing ensemble Kalman filter algorithm for filter (state estimation) design to the optimal control policy learning for the LQG control setting. And it is shown to be exact in its mean-field limit (N=infinity).