{'fraction_non_alphanumeric': 0.038900544038687195, 'fraction_numerical': 0.011911958183693063, 'mean_word_length': 4.485906563932508, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning robust control for LQR systems with multiplicative noise via policy gradient

The linear quadratic regulator (LQR) problem has reemerged as an important theoretical benchmark for reinforcement learning-based control of complex dynamical systems with continuous state and action spaces. In contrast with nearly all recent work in this area, we consider multiplicative noise models, which are increasingly relevant because they explicitly incorporate inherent uncertainty and variation in the system dynamics and thereby improve robustness properties of the controller. Robustness is a critical and poorly understood issue in reinforcement learning; existing methods which do not account for uncertainty can converge to fragile policies or fail to converge at all. Additionally, intentional injection of multiplicative noise into learning algorithms can enhance robustness of policies, as observed in ad hoc work on domain randomization. Although policy gradient algorithms require optimization of a non-convex cost function, we show that the multiplicative noise LQR cost has a special property called gradient domination, which is exploited to prove global convergence of policy gradient algorithms to the globally optimum control policy with polynomial dependence on problem parameters. Results are provided both in the model-known and model-unknown settings where samples of system trajectories are used to estimate policy gradients.


Introduction
Reinforcement learning-based control has recently achieved impressive successes in games [47,48] and simulators [41]. But these successes are significantly more challenging to translate to complex physical systems with continuous state and action spaces, safety constraints, and non-negligible operation and failure costs that demand data efficiency. An intense and growing research effort is creating a large array of models, algorithms, and heuristics for approaching the myriad of challenges arising from these systems. To complement a dominant trend of more computationally focused work, the canonical linear quadratic regulator (LQR) problem in control theory has reemerged as an important theoretical benchmark for learning-based control [45,17]. Despite its long history, there remain fundamental open questions for LQR with unknown models, and a foundational understanding of learning in LQR problems can give insight into more challenging problems.
Almost all recent work on learning in LQR problems has utilized either deterministic or additive noise models [45,17,19,12,20,1,35,51,2,53,38,54], but here we consider multiplicative noise models. In control theory, multiplicative noise models have been studied almost as long as their deterministic and additive noise counterparts [56,16], although this area is somewhat less developed and far less widely known. We believe the study of learning in LQR problems with multiplicative noise is important for three reasons. First, this class of models is much richer than deterministic or additive noise while still allowing exact solutions when models are known, which makes it a compelling additional benchmark [55,4,8]. Second, they explicitly incorporate model uncertainty and inherent stochasticity, thereby improving robustness properties of the controller. Robustness is a critical and poorly understood issue in reinforcement learning; existing methods which do not account for uncertainty can converge to fragile policies or fail to converge at all [4,34,9]. Additionally, intentional injection of multiplicative noise into learning algorithms is known to enhance robustness of policies from ad hoc work on domain randomization [49]. Third, in emerging difficult-to-model complex systems where learning-based control approaches are perhaps most promising, multiplicative noise models are increasingly relevant; examples include networked systems with noisy communication channels [3,25], modern power networks with large penetration of intermittent renewables [15,40], turbulent fluid flow [37], and neuronal brain networks [13].

Related literature
Multiplicative noise LQR problems have been studied in control theory since the 1960s [56]. Since then a line of research parallel to deterministic and additive noise has developed, including basic stability and stabilizability results [55], semidefinite programming formulations [18,10,36], robustness properties [16,8,27,6,24], and numerical algorithms [7]. This line of research is less widely known perhaps because much of it studies continuous time systems, where the heavy machinery required to formalize stochastic differential equations is a barrier to entry for a broad audience. Multiplicative noise models are well-poised to offer data-driven model uncertainty representations and enhanced robustness in learning-based control algorithms and complex dynamical systems and processes. A related line of research which has seen recent activity is on learning optimal control of Markovian jump linear systems with unknown dynamics and noise distributions [46,29], which under certain assumptions form a special case of the multiplicative noise system we analyze in this work.
In contrast to classical work on system identification and adaptive control, which has a strong focus on asymptotic results, more recent work has focused on non-asymptotic analysis using newly developed mathematical tools from statistics and machine learning. There remain fundamental open problems for learning in LQR problems, with several addressed only recently, including non-asymptotic sample complexity [17,51], regret bounds [1,2,38], and algorithmic convergence [19]. Alternatives to reinforcement learning include other data-driven model-free optimal control schemes [23,5] and those leveraging the behavioral framework [39,42]. Subspace identification methods offer a model-based generalization to the output feedback setting [30].

Our contributions
In §2 we establish the multiplicative noise LQR problem and motivate its study via a connection to robust stability. We then give several fundamental results for policy gradient algorithms on linear quadratic problems with multiplicative noise. Our main contributions are as follows, which can be viewed as a generalization of the recent results of Fazel et al. [19] for deterministic LQR to multiplicative noise LQR: • In §3 we show that although the multiplicative noise LQR cost is generally non-convex, it has a special property called gradient domination, which facilitates its optimization (Lemmas 3.1 and 3.3). • In particular, in §4 the gradient domination property is exploited to prove global convergence of three policy gradient algorithm variants (namely, exact gradient descent, "natural gradient descent, and Gauss-Newton/policy iteration) to the globally optimum control policy with a rate that depends polynomially on problem parameters (Theorems 4.1, 4.2, and 4.3). • Furthermore, in §5 we show that a model-free policy gradient algorithm, where the gradient is estimated from trajectory data ("rollouts") rather than computed from model parameters, also converges globally (with high probability) with an appropriate exploration scheme and sufficiently many samples (polynomial in problem data) (Theorem 5.1).
In comparison with the deterministic dynamics studied by [19], we make the following novel technical contributions: • We quantify the increase in computational burden of policy gradient methods due to the presence of multiplicative noise, which is evident from the bounds developed in Appendices B and C. The noise acts to reduce the step size and thus convergence rate, and increases the required number of samples and rollout length in the model-free setting.
• A covariance dynamics operator F K is established for multiplicative noise systems with a more complicated form than the deterministic case. This necessitated a more careful treatment and novel proof by induction and term matching argument in the proof of Lemma B.4. • Several restrictions on the algorithmic parameters (step size, number of rollouts, rollout length, exploration radius) which are necessary for convergence are established and treated. • An important restriction on the support of the multiplicative noise distribution, which is naturally absent in [19], is established in the model-free setting. • A matrix Bernstein concentration inequality is stated explicitly and used to give explicit bounds on the algorithmic parameters in the model-free setting in terms of problem data. • Discussion and numerical results on the use of backtracking line search is included.
• When the multiplicative variances α i , β j are all zero, the assertions of Theorems 4.1, 4.2, 4.3, 5.1 recover the same step sizes and convergence rates of the deterministic setting reported by [19].
Thus, policy gradient algorithms for the multiplicative noise LQR problem enjoy the same global convergence properties as deterministic LQR, while significantly enhancing the resulting controllers robustness to variations and inherent stochasticity in the system dynamics, as demonstrated by our numerical experiments in §6.
To our best knowledge, the present paper is the first work to consider and obtain global convergence results using reinforcement learning algorithms for the multiplicative noise LQR problem. Our approach allows the explicit incorporation of a model uncertainty representation that significantly improves the robustness of the controller compared to deterministic and additive noise approaches.

Optimal Control of Linear Systems with Multiplicative Noise and Quadratic Costs
We consider the infinite-horizon linear quadratic regulator problem with multiplicative noise (LQRm) where x t ∈ R n is the system state, u t ∈ R m is the control input, the initial state x 0 is distributed according to P 0 with covariance Σ 0 := E x0 [x 0 x 0 ], Σ 0 0, and Q 0 and R 0. The dynamics are described by a dynamics matrix A ∈ R n×n and input matrix B ∈ R n×m and incorporate multiplicative noise terms modeled by the i.i.d. (across time), zero-mean, mutually independent scalar random variables δ ti and γ tj , which have variances α i and β j , respectively. The matrices A i ∈ R n×n and B i ∈ R n×m specify how each scalar noise term affects the system dynamics and input matrices. Alternatively, supposeĀ andB are zero-mean random matrices with a joint covariance structure 1 over their entries governed by the covariance matrices Σ A := E[vec(Ā)vec(Ā) ] ∈ R n 2 ×n 2 and Σ B := E[vec(B)vec(B) ] ∈ R nm×nm . Then it suffices to take the variances α i and β j and matrices A i and B j as the eigenvalues and (reshaped) eigenvectors of Σ A and Σ B , respectively, after a projection onto a set of orthogonal real-valued vectors [57]. The goal is to determine a closed-loop state feedback policy π * with u t = π * (x t ) from a set Π of admissible policies which solves the optimization in (1).
We assume that the problem data A, B, α i , A i , β j , and B j permit existence and finiteness of the optimal value of the problem, in which case the system is called mean-square stabilizable and requires mean-square stability of the closed-loop system [33,55]. The system in (1) is called mean-square stable if lim t→∞ E x0,δ,γ [x t x t ] = 0 for any given initial covariance Σ 0 , where for brevity we notate expectation with respect to the noises E {δti},{γtj } as E δ,γ . Mean-square stability is a form of robust stability, implying stability of the mean (i.e. lim t→∞ Ex t = 0 ∀ x 0 ) as well as almost-sure stability (i.e. lim t→∞ x t = 0 almost surely) [55].
Mean-square stability requires stricter and more complicated conditions than stabilizability of the nominal system (A, B) [55], which are discussed in the sequel. This essentially can limit the size of the multiplicative noise covariance [4], which can be viewed as a representation of uncertainty in the nominal system model or as inherent variation in the system dynamics.

Control Design with Known Models: Value Iteration
Dynamic programming can be used to show that the optimal policy π * is linear state feedback u t = π * (x t ) = K * x t , where K * ∈ R m×n denotes the optimal gain matrix. When the control policy is linear state feedback u t = π(x t ) = Kx t , with a very slight abuse of notation the cost becomes Dynamic programming further shows that the resulting optimal cost is quadratic in the initial state, i.e. C(K * ) = E x0 x 0 P x 0 = Tr(P Σ 0 ), where P ∈ R n×n is a symmetric positive definite matrix [9]. Note that the optimal controller does not need to directly observe the noise variables δ ti , γ tj . When the model parameters are known, there are several ways to compute the optimal feedback gains and corresponding optimal cost. The optimal cost is given by the solution of the generalized algebraic Riccati equation (GARE) This is a special case of the GARE for optimal static output feedback given in [8] and can be solved via the value iteration with P 0 = Q, or via semidefinite programming formulations [10,18,36], or via more exotic iterations based on the Smith method and Krylov subspaces [22,28]. The associated optimal gain matrix is It was verified in [55] that existence of a positive definite solution to the GARE (2) is equivalent to mean-square stabilizability of the system, which depends on the problem data A, B, α i , A i , β j , and B j ; in particular, mean-square stability generally imposes upper bounds on the variances α i and β j [4], but may be infinite depending on the structure of A, B, A i , and B j [55]. At a minimum, uniqueness and existence of a solution to the GARE (2) requires the standard conditions for uniqueness and existence of a solution to the standard ARE, namely of (A, B) stabilizable and (A, Q 1/2 ) detectable.
Although (approximate) value iteration can be implemented using sample trajectory data, policy gradient methods have been shown to be more effective for approximately optimal control of high-dimensional stochastic nonlinear systems e.g. those arising in robotics [43]. This motivates our following analysis of the simpler case of stochastic linear systems wherein we show that policy gradient indeed facilitates a data-driven approach for learning optimal and robust policies.

Control Design with Known Models: Policy Gradient
Consider a fixed linear state feedback policy u t = Kx t . Defining the stochastic system matrices the deterministic nominal and stochastic closed-loop system matrices and the closed-loop state-cost matrix the closed-loop dynamics become A gain K is mean-square stabilizing if the closed-loop system is mean-square stable. Denote the set of mean-square stabilizing K as K. If K ∈ K, then the cost can be written as where P K is the unique positive semidefinite solution to the generalized Lyapunov equation We define the state covariance matrices and the infinite-horizon aggregate state covariance matrix as If K ∈ K then Σ K also satisfies a dual generalized Lyapunov equation Vectorization and Kronecker products can be used to convert (3) and (4) into systems of linear equations. Alternatively, iterative methods have been suggested for their solution [22,28]. The state covariance dynamics are captured by two closed-loop finite-dimensional linear operators which operate on a symmetric matrix X: Thus F K (without an argument) is a linear operator whose matrix representation is The Σ t evolve according to the dynamics We define the t-stage of F K (X) as F t K (X) := F K (F t−1 K (X)) with F 0 K (X) = X, which gives the natural characterization and in particular We then have the following lemma: Lemma 2.1 (Mean-square stability). A gain K is mean-square stabilizing if and only if the spectral radius ρ(F K ) < 1.
Proof. Mean-square stability implies lim t→∞ E[x t x t ] = 0, which for linear systems occurs only when Σ K is finite, which by (5) is equivalent to ρ(F K ) < 1.
Recalling the definition of C(K) and (4), along with the basic observation that K / ∈ K induces infinite cost, gives the following characterization of the cost: The evident fact that C(K) is expressed as a closed-form function, up to a Lyapunov equation, of K leads to the idea of performing gradient descent on C(K) (i.e., policy gradient) via the update K ← K − η∇C(K) to find the optimal gain matrix. However, two properties of the LQR cost function C(K) complicate a convergence analysis of gradient descent. First, C(K) is extended valued since not all gain matrices provide closed-loop mean-square stability, so it does not have (global) Lipschitz gradients. Second, and even more concerning, C(K) is generally non-convex in K (even for deterministic LQR problems, as observed by Fazel et al. [19]), so it is unclear if and when gradient descent converges to the global optimum, or if it even converges at all. Fortunately, as in the deterministic case, we show that the multiplicative LQR cost possesses further key properties that enable proof of global convergence despite the lack of Lipschitz gradients and non-convexity.

From Stochastic to Robust Stability
Additional motivation for designing controllers which stabilize a stochastic system in mean-square is to ensure robustness of stability of a nominal deterministic system to model parameter perturbations. Here we state a condition which guarantees robust deterministic stability for a perturbed deterministic system given mean-square stability of a stochastic single-state system with multiplicative noise where the noise variance and parameter perturbation size are related.
Example 2.2 (Robust stability). Suppose the stochastic closed-loop system where a, x t , δ t are scalars with E[δ 2 t ] = α is mean-square stable. Then, the perturbed deterministic system is stable for any constant perturbation |φ| ≤ √ a 2 + α − |a|.
Although this is a simple example, it demonstrates that the robustness margin increases monotonically with the multiplicative noise variance. We also see that when α = 0 the bound collapses so that no robustness is guaranteed, i.e., when |a| → 1. This result can be extended to multiple states, inputs, and noise directions, but the resulting conditions become considerably more complex [8,24]. We now proceed with developing methods for optimal control.

Gradient Domination and Other Properties of the Multiplicative Noise LQR Cost
In this section, we demonstrate that the multiplicative noise LQR cost function is gradient dominated, which facilitates optimization by gradient descent. Gradient dominated functions have been studied for many years in the optimization literature [44] and have recently been discovered in deterministic LQR problems by [19].

Multiplicative Noise LQR Cost is Gradient Dominated
First, we give the expression for the policy gradient of the multiplicative noise LQR cost. 2 Define Lemma 3.1 (Policy Gradient Expression). The policy gradient is given by Proof. Substituting the RHS of the generalized Lyapunov equation (3) into the cost C(K) = Tr(P K Σ 0 ) yields Taking the gradient with respect to K and using the product rule and rules for matrix derivatives we obtain We include a factor of 2 on the gradient expression that was erroneously dropped in [19]. This affects the step size restrictions by a corresponding factor of 2.
where the tilde on K and overbar on K are used to denote the terms being differentiated. Applying this gradient formula recursively to the last term in the last line (namely ∇K Tr(PKΣ 1 )), and recalling the definition of Σ K we obtain For brevity the gradient is implied to be with respect to the gains K in the rest of this work, i.e., ∇ K denoted by ∇. Now we must develop some auxiliary results before demonstrating gradient domination. Throughout Z and Z F are the spectral and Frobenius norms respectively of a matrix Z, and σ(Z) and σ(Z) are the minimum and maximum singular values of a matrix Z. The value function of V K (x) for x = x 0 is defined as which relates to the cost as C(K) = E x0 V K (x 0 ). The advantage function is defined as where the expectation is with respect to A and B inside the parentheses of V K ( Ax + Bu). The advantage function can be thought of as the difference in cost ("advantage") when starting in state x of taking an action u for one step instead of the action generated by policy K. We also define the state, input, and cost sequences Throughout the proofs we will consider pairs of gains K and K and their difference ∆ := K − K.
Lemma 3.2 (Value difference). Suppose K and K generate the (stochastic) state, action, and cost sequences Then the value difference and advantage satisfy Proof. The proof follows the "cost-difference" lemma in [19] exactly substituting versions of value and cost functions, etc. which take expectation over the multiplicative noise. By definition we have so we can write the value difference as We expand the following value function difference as where the last equality is valid by noting that the first term in sequence {x t } K ,x is x. Continuing the value difference expression we have where the fifth equality holds since For the second part of the proof regarding the advantage expression, we expand and substitute in definitions: where the third equality follows from all of the δ ti and γ tj being zero-mean and mutually independent. Substituting and continuing, We also have the following expression from the recursive relationship for P K Substituting, we cancel the V K (x) term which leads to the result after rearrangement: which completes the proof.
Next, we see that the multiplicative noise LQR cost is gradient dominated.

Lemma 3.3 (Gradient domination).
The LQR-with-multiplicative-noise cost C(K) satisfies the gradient domination condition Proof. We start with the advantage expression Next we rearrange and complete the square: with equality only when ∆ = −R −1 K E K . Let the state and control sequences associated with the optimal gain K * be {x t } K * ,x and {u t } K * ,x respectively. We now obtain an upper bound for the cost difference by writing the cost difference in terms of the value function as Using the first part of the value-difference Lemma 3.2 and negating we obtain where the second step used the advantage inequality in (8). Now using | Tr(Y Z)| ≤ Y | Tr(Z)| we obtain where the first and second inequalities will be used later in the Gauss-Newton and gradient descent convergence proofs respectively. Combining which will be used later in the natural policy gradient descent convergence proof. Now we rearrange and substitute in the policy gradient expression 1 where the last step used the definition and submultiplicativity of spectral norm. Using completes the proof.
The gradient domination property gives the following stationary point characterization.
In other words, so long as Σ K is full rank, stationarity is both necessary and sufficient for global optimality, as for convex functions. Note that it is not sufficient to just have multiplicative noise in the dynamics with a deterministic initial state x 0 to ensure that Σ K is full rank. To see this, observe that if x 0 = 0 and Σ 0 = 0 then Σ K = 0, which is clearly rank deficient. By contrast, additive noise is sufficient to ensure that Σ K is full rank with a deterministic initial state x 0 , although we will not consider this setting. Using a random initial state with Σ 0 0 ensures rank(Σ K ) = n and thus ∇C(K) = 0 implies K = K * .
Although the gradient of the multiplicative noise LQR cost is not globally Lipschitz continuous, it is locally Lipschitz continuous over any subset of its domain (i.e., over any set of mean-square stabilizing gain matrices). The gradient domination is then sufficient to show that policy gradient descent will converge to the optimal gains at a linear rate (a short proof of this fact for globally Lipschitz functions is given in [32]). We prove this convergence of policy gradient to the optimum feedback gain by bounding the local Lipschitz constant in terms of the problem data, which bounds the maximum step size and the convergence rate.

Additional Setup Lemmas
Following [19] we refer to Lipschitz continuity of the gradient as (C 1 -)smoothness, and so this section deals with showing that the LQR cost satisfies an expression that is almost of the exact form of a Lipschitz continuous gradient.
Lemma 3.5 (Almost-smoothness). The LQR-with-multiplicative-noise cost C(K) satisfies the almost-smoothness expression Proof. As in the gradient domination proof, we express the cost difference in terms of the advantage by taking expectation over the initial states to obtain From the value difference lemma for the advantage we have Noting that {u t } K ,x = K x we substitute to obtain Using the definition of Σ K completes the proof.
Remark 3.6. For small deviations K − K the equation in the almost-smoothness lemma exactly describes a Lipschitz continuous gradient. The naming should not be taken to imply that the LQRm cost is not smooth, but rather that the equation as stated does not immediately yield a Lipschitz constant; indeed the Lipschitz constant is what much of the later proofs go towards bounding (implicitly) i.e. by bounding higher-order terms which must be accounted for when K = K.
To be specific, a Lipschitz continuous gradient to C(K) implies there exists a Lipschitz constant L such that for all K , K. This is the quadratic upper bound which is used e.g. in Thm. 1 of [32] to prove convergence of gradient descent on a gradient dominated objective function. The "almost"-smoothness condition is that which is exactly of the form of the Lipschitz gradient condition with Lipschitz constant 2 Σ K R K . Note this is not a global Lipschitz condition since 2 Σ K R K becomes unbounded as K becomes mean-square destabilizing, but rather a local Lipschitz condition since 2 Σ K R K is bounded on any sublevel set of C(K).
Lemma 3.7 (Cost bounds). We always have Proof. The proof follows that in [19] exactly. The cost is lower bounded as which gives the first inequality. The cost is also lower bounded as which gives the second inequality.

Global Convergence of Policy Gradient in the Model-Based Setting
In this section we show that the policy gradient algorithm and two important variants for multiplicative noise LQR converge globally to the optimal policy. In contrast with [19], the policies we obtain are robust to uncertainties and inherent stochastic variations in the system dynamics. We analyze three policy gradient algorithm variants:

Ks
The more elaborate natural gradient and Gauss-Newton variants provide superior convergence rates and simpler proofs. A development of the natural policy gradient is given in [19] building on ideas from [31]. The Gauss-Newton step with step size 1 2 is in fact identical to the policy improvement step in policy iteration (a short derivation is given shortly) and was first studied for deterministic LQR in [26]. This was extended to a model-free setting using policy iteration and Q-learning in [12], proving asymptotic convergence of the gain matrix to the optimal gain matrix. For multiplicative noise LQR, we have the following results.

Derivation of the Gauss-Newton step from policy iteration
We start with the policy improvement expression for the LQR problem: Stationary points occur when the gradient is zero, so differentiating with respect to u we obtain Setting (11) to zero and solving for u gives Differentiating (11) with respect to u we obtain confirming that the stationary point is indeed a global minimum.
Thus the policy iteration gain matrix update is This can be re-written in terms of the gradient as so: Parameterizing with a step size gives the Gauss-Newton step

Gauss-Newton Descent
Theorem 4.1 (Gauss-Newton convergence). Using the Gauss-Newton step Ks with step size 0 < η ≤ 1 2 gives global convergence to the optimal gain matrix K * at a linear rate described by Proof. The next-step gain matrix difference is Using the almost-smoothness Lemma 3.5 and substituting in the next-step gain matrix difference we obtain Ks E Ks . Recalling and substituting in (9) we obtain Adding C(K s ) − C(K * ) to both sides and rearranging completes the proof.

Natural Policy Gradient Descent
Theorem 4.2 (Natural policy gradient convergence). Using the natural policy gradient step with step size 0 < η ≤ c npg where gives global convergence to the optimal gain matrix K * at a linear rate described by Proof. First we bound the one-step progress, where the step size depends explicitly on the current gain K s . Using the update (12), the next-step gain matrix difference is Ks . Using Lemma 3.5 and substituting we obtain Recalling and substituting (10) we obtain Adding C(K s ) − C(K * ) to both sides and rearranging gives the one step progress bound Next, using the cost bound in Lemma 3.7, the triangle inequality, and submultiplicativity of spectral norm we have 1 .
Accordingly, choosing the step size as 0 < η ≤ c npg ensures (13) holds at the first step. This ensures that which allows (13) to be applied at the next step as well. Proceeding inductively by applying (13) at each successive step completes the proof.

Policy Gradient Descent
Theorem 4.3 (Policy gradient convergence). Using the policy gradient step with step size 0 < η ≤ c pg gives global convergence to the optimal gain matrix K * at a linear rate described by Proof. The proof is developed in Appendix B.
The proofs for these results explicitly incorporate the effects of the multiplicative noise terms δ ti and γ tj in the dynamics. For the policy gradient and natural policy gradient algorithms, we show explicitly how the maximum allowable step size depends on problem data and in particular on the multiplicative noise terms. Compared to deterministic LQR, the multiplicative noise terms decrease the allowable step size and thereby decrease the convergence rate; specifically, the state-multiplicative noise increases the initial cost C(K 0 ) and the norms of the covariance Σ K * and cost P K , and the input-multiplicative noise also increases the denominator term B 2 + q j=1 β j B j 2 . This means that the algorithm parameters for deterministic LQR in [19] may cause failure to converge on problems with multiplicative noise. Moreover, even the optimal policies for deterministic LQR may actually destabilize systems in the presence of small amounts of multiplicative noise uncertainty, indicating the possibility for a catastrophic lack of robustness; observe the results of the example in Section 6.1. The results and proofs also differ from that of [19] because the more complicated mean-square stability must be accounted for, and because generalized Lyapunov equations must be solved to compute the gradient steps, which requires specialized solvers.

Global Convergence of Policy Gradient in the Model-Free Setting
The results in the previous section are model-based; the policy gradient steps are computed exactly based on knowledge of the model parameters. In the model-free setting, the policy gradient is estimated to arbitrary accuracy from sample trajectories with a sufficient number of sample trajectories n sample of sufficiently long horizon length using gain matrices randomly selected from a Frobenius-norm ball around the current gain of sufficiently small exploration radius r. We show for multiplicative noise LQR that with a finite number of samples polynomial in the problem data, the model-free policy gradient algorithm still converges to the globally optimal policy, despite small perturbations on the gradient.
In the model-free setting, the policy gradient method proceeds as before except that at each iteration Algorithm 1 is called to generate an estimate of the gradient via the zeroth-order optimization procedure described by Fazel et al. [19].

Algorithm 1 Model-Free policy gradient estimation
Input: Gain matrix K, number of samples n sample , rollout length , exploration radius r 1: for i = 1, . . . , n sample do 2: Generate a sample gain matrix K i = K + U i , where U i is drawn uniformly at random over matrices with Frobenius norm r 3: Generate a sample initial state x Simulate the closed-loop system for steps according to the stochastic dynamics in (1) starting from x

5:
Collect the empirical finite-horizon cost estimate Theorem 5.1 (Model-free policy gradient convergence). Let and µ be a given small tolerance and probability respectively and N be the number of gradient descent steps taken. Suppose that the distribution of the initial states is bounded such that x 0 ∼ P 0 implies x i 0 ≤ L 0 almost surely for any given realization x i 0 of x 0 . Suppose additionally that the distribution of the multiplicative noises is bounded such that the following inequality is satisfied almost surely for any given realized sequence x i t of x t with a positive scalar z ≥ 1: under the closed-loop dynamics with any gain such that C(K) ≤ 2C(K 0 ). Suppose the step size η is chosen according to the restriction in Theorem 4.3 and at every iteration the gradient is estimated according to the finite-horizon procedure in Algorithm 1 where the number of samples n sample , rollout length , and exploration radius r are chosen according to the fixed polynomials of the problem data A, B, α i , β j , A i , B j , Q, R, Σ 0 , K 0 , L 0 and z which are all defined in the proofs in Appendix C. Then, with high probability of at least 1 − µ, performing gradient descent results in convergence to the global optimum over all N steps: at each step, either progress is made at the linear rate or convergence has been attained with C(K s ) − C(K * ) ≤ .
Proof. The proof is developed in Appendix C.
From a sample complexity standpoint, it is notable that the number of samples n sample , rollout length , and exploration radius r in Theorem 5.1 are polynomial in the problem data A, B, α i , β j , A i , B j , Q, R, Σ 0 , C(K 0 ). The constant z imposes a bound on the multiplicative noise, which is naturally absent in [19]. Note that z ≥ 1 since any upper bound of a scalar distribution with finite support must be equal to or greater than the mean. In general, this implicitly requires the noises to have bounded support. Such an assumption is qualitatively the same as the condition imposed on the initial states. These assumptions are reasonable; in a practical setting with a physical system the initial state and noise distributions will have finite support. There is no restriction on how large the support is, only that it not be unbounded. Also note that the rate is halved compared with the model-based case of Theorem 4.3; this is because the "other half" is consumed by the error between the estimated and true gradient.

Numerical Experiments
In this section we present results for three systems: A. Shows that "optimal" control that ignores actual multiplicative noise can lead to loss of mean-square stability, B. Shows the efficacy of the policy gradient algorithms on a networked system, C. Shows the increased difficulty of estimating the gradient from sample data in the presence of multiplicative noise.
All systems considered permitted a solution to the GARE (2). The bounds on the step size, number of rollouts, and rollout length given by the theoretical analysis can be rather conservative. For practicality, we selected the constant step size, number of rollouts, rollout length, and exploration radius according to a grid search over reasonable values. Additionally, we investigated the use of backtracking line search to adaptively select the step size; see e.g. [11]. Throughout the simulations, we computed the baseline optimal cost C(K * ) by solving the GARE (2) to high precision via value iteration. Python code which implements the algorithms and generates the figures reported in this work can be found in the GitHub repository at https://github.com/TSummersLab/polgrad-multinoise/. The code was run on a desktop PC with a quad-core Intel i7 6700K 4.0GHz CPU, 16GB RAM; no GPU computing was utilized.

Importance of Accounting for Multiplicative Noise
We first considered an open-loop mean-square unstable system with four states and one input representing an active two-mass suspension converted from continuous to discrete time using a standard bilinear transformation, with parameters: We performed model-based policy gradient descent; at each iteration gradients were calculated by solving generalized Lyapunov equations (3) and (4) using the problem data. The gains K m and K represent iterates during optimization of ("training" on) the LQRm and LQR cost (with the multiplicative noise variances set to zero), respectively. We performed the optimization starting from the same feasible initial gain, which was generated by perturbing the exact solution of the generalized algebraic Riccati equation such that the LQRm cost under the initial control was approximately 10 times that of the optimal control. The step size was chosen via backtracking line search. The optimization stopped once the Frobenius norm of the gradient fell below a small threshold. The plot in Fig. 1 shows the "testing" cost of the gains at each iteration evaluated on the LQRm cost (with multiplicative noise). From this figure, it is clear that K m minimized the LQRm as desired.
When there was high multiplicative noise, the noise-ignorant controller K actually destabilized the system in the mean-square sense; this can be seen as the LQRm cost exploded upwards to infinity after iteration 10. In this sense, the multiplicative noise-aware optimization is generally safer and more robust than noise-ignorant optimization, and in examples like this is actually necessary for mean-square stabilization.

Policy Gradient Methods Applied to a Network
Many practical networked systems can be approximated by diffusion dynamics with losses and stochastic diffusion constants (edge weights) between nodes; examples include heat flow through uninsulated pipes, hydraulic flow through leaky pipes, information flow between processors with packet loss, electrical power flow between generators with resistant electrical power lines, etc. A derivation of the discrete-time dynamics of this  vs. iteration during policy gradient descent on the 4-state, 1-input suspension example system. system is given in [57]. We considered a particular 4-state, 4-input system and open-loop mean-square stable with the following parameters:  This system is open-loop mean-square stable, so we initialized the gains to all zeros for each trial. We performed policy optimization using the model-free gradient, and the model-based gradient, model-based natural gradient, and model-based Gauss-Newton step directions on 20 unique problem instances using two step size schemes: Backtracking line search: Step sizes η were chosen adaptively at each iteration by backtracking line search with parameters α = 0.01, β = 0.5 (see [11] for a description), except for Gauss-Newton which used the optimal constant step-size of 1/2. Model-free gradients and costs were estimated with 100,000 rollouts per iteration. We ran a fixed number, 20, of iterations chosen such that the final cost using model-free gradient descent was no more than 5% worse than optimal. Constant step size: Step sizes were set to constants chosen as large as possible without observing infeasibility or divergence, which on this problem instance was η = 5 × 10 −5 for gradient, η = 2 × 10 −4 for natural gradient, and η = 1/2 for Gauss-Newton step directions. Model-free gradients were estimated with 1,000 rollouts per iteration. We ran a fixed number, 20,000, of iterations chosen such that convergence was achieved with all step directions.
In both cases sample gains were chosen for model-free gradient estimation with exploration radius r = 0.1 and the rollout length was set to = 20. The plots in Fig. 2 show the relative cost over the iterations; for the model-free gradient descent, the bold centerline is the mean of all trials and the shaded region is between the 10 th and 90 th percentile of all trials. Using backtracking line search, it is evident that in terms of convergence the Gauss-Newton step was extremely fast, and both the natural gradient and model-based gradient were slightly slower, but still quite fast. The model-free policy gradient converged to a reasonable neighborhood of the minimum cost quickly, but stagnated with further iterations; this is a consequence of the inherent gradient and cost estimation errors that arise due to random sampling and the multiplicative noise. Using constant stepsizes, we were forced to take small steps due to the steepness of the cost function near the initial gains, slowing overall convergence using the gradient and natural gradient methods. Here we observed that Gauss-Newton again converged most quickly, followed by natural gradient and lastly the gradient methods. The smaller step size also allowed us to use far fewer samples in the model-free setting, where we observed somewhat faster initial cost decrease with eventual stagnation around 10 −2 , or 1%, relative error, which represents excellent control performance. All algorithms exhibited convergence to the optimum, confirming the asserted theoretical claims.  vs. iteration during policy gradient methods on a 4-state, 4-input lossy diffusion network with multiplicative noise using a) backtracking line search and b) constant step sizes.

Gradient Estimation
Multiplicative noise can significantly increase the variance and sample complexity of cost gradient estimates relative to the noiseless case, which is novelly reflected in the theoretical analysis for the number of rollouts and rollout length. To demonstrate this empirically, we evaluated the relative gradient estimation error vs. number of rollouts for the system with K = 0, Q = Σ 0 = I 2 , R = 1, δ t ∼ N (0, 0.1), rollout length l = 40, exploration radius r = 0.2, averaged over 10 gradient estimates. The results are plotted in Figure 3. To achieve the same gradient estimate error of 10%, the system with multiplicative noise required 200× the number of rollout samples (10 8 ) as when there was no noise (5 × 10 5 ).

Conclusions
We have shown that policy gradient methods in both model-based and model-free settings give global convergence to the globally optimal policy for LQR systems with multiplicative noise. These techniques are directly applicable for the design of robust controllers of uncertain systems and serve as a benchmark for data-driven control design. Our ongoing work is exploring ways of mitigating the relative sample inefficiency of model-free policy gradient methods by leveraging the special structure of LQR models and Nesterov-type acceleration, and exploring alternative system identification and adaptive control approaches. We are also investigating other methods of building robustness through H ∞ and dynamic game approaches. Another extension relevant to networked control systems is enforcing sparse structure constraints on the gain matrix via projected policy gradient as suggested in [14].

Appendix A. Standard matrix expressions
Before proceeding with the proof of the main results of this study, we first review several basic expressions that will be used later throughout the section. In this section we let A, B, C, M i be generic matrices ∈ R n×m , a, b be generic vectors, and s be a generic scalar.

Spectral norm:
We denote the matrix spectral norm as A = σ max (A) which clearly satisfies

Frobenius norm:
We denote the matrix Frobenius norm as A F whose square satisfies Frobenius norm ≥ spectral norm: For any matrix A the Frobenius norm is greater than or equal to the spectral norm: Inverse of spectral norm inequality: Invariance of trace under cyclic permutation: Trace-spectral norm inequalities: Sub-multiplicativity of spectral norm: Positive semidefinite matrix inequality: Suppose A 0 and B 0. Then Vector self outer product positive semidefiniteness: This is the same inequality given in [19]. See [50] for the exact scale constants and a proof.
Appendix B. Model-based policy gradient descent The proof of convergence using gradient descent proceeds by establishing several technical lemmas, bounding the infinite-horizon covariance Σ K , then using that bound to limit the step size, and finally obtaining a one-step bound on gradient descent progress and applying it inductively at each successive step. We begin with a bound on the induced operator norm of T K : Lemma B.1. (T K norm bound) The following bound holds for any mean-square stabilizing K: Proof. The proof follows that given in [19] using our definition of T K .
Lemma B.2. (F K perturbation) Consider a pair of mean-square stabilizing gain matrices K and K . The following F K perturbation bound holds: Proof. Let ∆ = −∆. For any matrix X we have The operator norm F K − F K is Applying submultiplicativity of spectral norm to (16) and noting that ∆ = ∆ completes the proof.
Lemma B.3 (T K perturbation). If K and K are mean-square stabilizing and Proof. The proof follows [19] using our modified definitions of T K and F K .
We have a generic inequality for a sum of n matrices M i : where the last step is due to the triangle inequality. Recalling the definitions of F t K (I) and F t K we see they are of the form of the LHS and RHS in (18) with all terms matched between F t K (I) and F t K so that the inequality in (18) holds; this can be seen by starting with t = 1 and incrementing t up by 1 which will give (1 + p + q) t terms which are all matched. Thus, Continuing from (17) we have By hypothesis ρ(F K ) < 1, and taking the sum of the geometric series completes the proof.
then the associated state covariance matrices satisfy Proof. First we show that K is mean-square stabilizing and ∆ ≤ h ∆ (K) then K is also mean-square stabilizing. This follows from an analogous argument in [19] by characterizing mean-square stability in terms of ρ(F K ) rather than ρ(A K ) and using Lemma B.4. Let K be distinct from K with ρ(F K ) < 1 and Since K and K are mean-square stabilizing Lemma B.5 holds so we have where we define Using Lemma B.4 we have Rearranging and substituting for Γ, Now we construct the proof by contradiction. Suppose there is a K with ρ(F K ) > 1 satisfying the perturbation restriction Since spectral radius is a continuous function (see [52]) there must be a point K on the path between K and K such that ρ(F K ) = 1 − < 1. Since K and K are mean-square stabilizing Lemma B.5 holds so we have and rearranging However since K is mean-square stabilizing Lemma B.4 holds so we have which is a contradiction. Therefore no such mean-square unstable K satisfying the hypothesized perturbation restriction can exist, completing the first part of the proof.
The rest of the proof follows [19] by using the condition on ∆ , Σ K ≥ σ(Σ 0 ), and Lemmas 3.7, B.1, and B.3. The condition on ∆ directly implies where the last step is due to the combination of Lemma 3.7 with Σ K ≥ σ(Σ 0 ). By Lemma B.2 we have Combining this with Lemma B.1 we have By the condition on ∆ we have 2 which allows us to use Lemma B.3 by which we have where the last step used (19).
Using Lemma 3.7 completes the proof.
Now we bound the one step progress of policy gradient where we allow the step size to depend explicitly on the current gain matrix iterate K s . Lemma B.6 (Gradient descent, one-step). Using the policy gradient step update K s+1 = K s − η∇C(K s ) with step size gives the one step progress bound Proof. The gradient update yields ∆ = −2ηE Ks Σ Ks . Putting this into Lemma 3.5 gives where the last step is due to σ(Σ 0 ) ≤ σ(Σ Ks ) and Lemma 3.3. Note that the assumed condition on the step size ensures the gain matrix difference satisfies the condition for Lemma B.5 as follows: where the last inequality is due to Lemma 3.7. Thus we can indeed apply Lemma B.5, by which we have where the last inequality is due to using the substitution ∆ = η ∇C(K s ) and the hypothesized condition on η. Using this and Lemma 3.7 we have where the second-to-last inequality used the hypothesized condition on η. Therefore Adding 1 to both sides completes the proof.
Lemma B.7 (Cost difference lower bound). The following cost difference inequality holds: Proof. The proof follows that for an analogous condition located in the gradient domination lemma in [19]. Let K and K generate the (stochastic) state and action sequences {x t } K,x , {u t } K,x and {x t } K ,x , {u t } K ,x respectively. By definition of the optimal gains we have C(K * ) ≤ C(K ). Then by Lemma 3.2 we have

holds with equality as
Lemma B.8. The following inequalities hold: where h 0 (K), h 1 (K), h 2 (K) are the polynomials Proof. The proof follows [19] with R K defined here. From the policy gradient expression we have Using Lemma 3.7 we have Tr(E K E K ), and using Lemma B.7 we have Taking square roots completes the proof of the first part of the lemma. For the second part of the lemma, we have where the second line is due to Weyl's inequality for singular values [19]. Using Lemma B.7 again on the Tr(E K E K ) term completes the proof.
We now give the parameter and proof of global convergence of policy gradient descent in Theorem 4.3.
Theorem B.9 (Policy gradient convergence). Consider the assumptions and notations of Theorem 4.3 and define Then the claim of Theorem 4.3 holds.
Proof. We have by Weyl's inequality for singular values (15), submultiplicativity of spectral norm, and Lemma B.8 that Thus by choosing 0 < η ≤ c pg we satisfy the requirements for Lemma B.6 at s = 1, which implies that progress is made at s = 1, i.e., that C(K 1 ) ≤ C(K 0 ) according to the rate in Lemma B.6. Proceeding inductively and applying Lemma B.6 at each step completes the proof.
Remark B.10. The quantities h 1 , h 2 , and R K may be upper bounded by quantities that depend only on problem data and C(K 0 ) e.g. using the cost bounds in Lemma 3.7, which we omit for brevity, so a conservative minimum step size η may be computed exactly.
Appendix C. Model-free policy gradient descent The overall proof technique of showing high probability convergence of model-free policy gradient proceeds by showing that the difference between estimated gradient and true gradient is bounded by a sufficiently small value that the iterative descent progress remains sufficiently large to maintain a linear rate. We use a matrix Bernstein inequality to ensure that enough sample trajectories are collected to estimate the gradient to high enough accuracy.
We begin with a lemma shows that C(K) and Σ K can be estimated with arbitrarily high accuracy as the rollout length increases.
Lemma C.1 (Approximating C(K) and Σ K with infinitely many finite horizon rollouts). Let be an arbitrary small constant. Suppose K gives finite C(K). Define the finite-horizon estimates Σ ( ) where expectation is with respect to x 0 , {δ ti }, {γ tj }. Then the following hold: Proof. The proof follows [19] exactly using suitably modified definitions of C(K), T K , F K .
Next we bound cost and gradient perturbations in terms of gain matrix perturbations and problem data. Using the same restriction as in Lemma B.5 we have Lemmas C.2 and C.3.
Lemma C.2 (C(K) perturbation). If ∆ ≤ h ∆ (K), then the cost difference is bounded as Proof. The proof follows [19] using suitably modified definitions of C(K), T K , F K , however compared with [19] we terminate the proof bound earlier so as to avoid a degenerate bound in the case of K = 0, and we also correct typographical errors. Note that ∆ has a more restrictive upper bound due to the multiplicative noise.
Lemma C.3 (∇C(K) perturbation). If ∆ ≤ h ∆ (K), then the policy gradient difference is bounded as Proof. The proof generally follows [19] using Lemmas B.5, C.2, and B.7 with R K and E K modified appropriately.
Recalling ∇C(K) = 2E K Σ K , and using the triangle inequality, First we bound the second term of (20). By Lemma B.7 Since ∆ ≤ h ∆ (K) Lemma B.5 holds and we have Therefore the second term is bounded as Now we bound the first term of (20). Since ∆ ≤ h ∆ (K) Lemma B.5 holds and by the reverse triangle inequality Rearranging and using Lemma 3.7 gives By an intermediate step in the proof of Lemma C.2 we have Expanding the difference E K − E K gives Using the reverse triangle inequality with the ∆ bound yields Substituting in (21) and P K ≤ C(K) σ(Σ0) from Lemma 3.7: Using (E K − E K ) Σ K ≤ E K − E K Σ K and adding the two terms of (20) completes the proof.
We now discuss smoothing of the cost in the context of model-free gradient descent. As in [19], in the model-free setting we apply Frobenius-norm ball smoothing to the cost. Let S r be the uniform distribution over all matrices with Frobenius norm r (the boundary of the ball), and B r be the uniform distribution over all matrices with Frobenius norm at most r (the entire ball). The smoothed cost is where U is a random matrix with the same dimensions as K and Frobenius norm r. The following lemma shows that the gradient of the smoothed function can be estimated just with an oracle of the function value.
Lemma C.4 (Zeroth-order gradient estimation). The gradient of the smoothed cost is related to the unsmoothed cost by Proof. The result is proved in Lemma 2.1 in [21]. For completeness, we prove it here. By Stokes' Theorem we have By definition we have C r (K) = Br C(K + U )dU vol d (B r ) .
We also have E U ∼Sr [C(K + U )U ] = rE U ∼Sr C(K + U ) U r = r δSr C(K + U ) U r dU vol d−1 (S r ) and the ratio of ball surface area to volume vol d−1 (S r ) vol d (B r ) = d r .
Combining, we have Lemma C.4 shows that the gradient of the smoothed cost can be found exactly with infinitely many infinite-horizon rollouts. Much of the remaining proofs goes towards showing that the error between the gradient of the smoothed cost and the unsmoothed cost, the error due to using finite-horizon rollouts, and the error due to using finitely many rollouts can all be bounded by polynomials of the problem data. As noted by [19] the reason for smoothing in a Frobenius norm ball rather than over a Gaussian distribution is to ensure stability and finiteness of the cost of every gain within the smoothing domain, although now in the multiplicative noise case we must be even more restrictive about our choice of perturbation on K because we require not only mean stability, but also mean-square stability.
We now give a Bernstein inequality for random matrices and some simple variants; this allows us to bound the difference between the sample average of a random matrix and its expectation.
Lemma C.5 (Matrix Bernstein inequality [50]). Let {Z i } N i=1 be a set of N independent random matrices of dimension d 1 × d 2 with d + = d 1 + d 2 , E[Z i ] = 0, Z i ≤ R almost surely, and maximum variance Proof. The lemma and proof follow [50] exactly. Proof. First note that K − K F = ∆ F = U F = r. We break the difference between estimated and true gradient∇C(K) − ∇C(K) into two terms as ∇C r (K) − ∇C(K) + ∇ C(K) − ∇C r (K) .
Since r ≤ h ∆ (K) we see that Lemmas C.2 and C.3 hold. By enforcing the bound r ≤ 1 hcost(K) , by Lemma C.2 and noting that ∆ ≤ ∆ F we have This ensures stability of the system under the perturbed gains so that C(K + U ) is well-defined. For the first term ∇C r (K) − ∇C(K), by enforcing r ≤ 2h grad (K) , by Lemma C.3 we have Since ∇C r (K) is the expectation of ∇C(K + U ), by the triangle inequality we have For the second term∇C(K) − ∇C r (K), we work towards using the matrix Bernstein inequality and adopt the notation of the associated lemma. First note that by Lemma C.4 we have Z := ∇C r (K) = E[∇C(K)]. Each individual sample Z i := mn r 2 C(K + U i )U i has the bounded Frobenius norm Proof. Similar to before, we break the difference between estimated and true gradient into three terms as Thus the choice of r and n sample satisfy the conditions of Lemma C.8, so with high probability of at least 1 − µ For the second term, by using the choices ≥ h r 4mn and C(K + U i ) ≤ 2C(K), Lemma C.1 holds and implies that By the triangle inequality, submultiplicativity, and U i F ≤ r 1 n sample For the first term, Applying this to the cost, summing over time and using the assumed restriction on the multiplicative noise we have Using this and an argument identical to Lemma C.8, each sample Z i := mn r 2 Ĉ (K + U i )U i has bounded Frobenius norm .