Structured Policy Iteration for Linear Quadratic Regulator

The general method for solving these problems is dynamic programming (DP) (Bellman et al., 1954; Bertsekas et al., 1995). To overcome the issues of DP such as intractability and computational cost, the common alternatives are approximated dynamic programming (ADP) (Bertsekas & Tsitsiklis, 1996; Bertsekas & Shreve, 2004; Powell, 2007; De Farias & Van Roy, 2003; O’Donoghue et al., 2011) or Reinforcement Learning (RL) (Sutton et al., 1998) including policy gradient (Kakade, 2002; Schulman et al., 2017) or Q-learning based methods (Watkins & Dayan, 1992; Mnih et al., 2015).

LQR.   LQR is a classical problem in control that is able to capture problems in continuous state and action space, pioneered by Kalman in the late 1950’s (Kalman, 1964). Since then, many variations have been suggested and solved such as jump LQR, random LQR, (averaged) infinite horizon objective, etc (Florentin, 1961; Costa et al., 2006; Wonham, 1970). When the model is (assumed to be) known, Arithmetic Riccati Equation (ARE) (Hewer, 1971; Laub, 1979) can be used to efficiently compute the optimal value function and policy for generic LQR. Alternatively, we can use eigendecomposition (Lancaster & Rodman, 1995) or transform it into a semidefinite program (SDP) (Balakrishnan & Vandenberghe, 2003).

LQR under model uncertainty.   When the model is unknown, one class is model-based approaches where the transition dynamics is attempted to be estimated. For LQR, in particular, (Abbasi-Yadkori & Szepesvári, 2011; Ibrahimi et al., 2012) developed online algorithms with regret bounds where linear dynamic (and cost function) are learned and sampled from some confidence set, but without any computational demonstrations. Another line of work is to utilize robust control with system identification (Dean et al., 2017; Recht, 2019) with sample complexity analysis (Tu & Recht, 2017). Certainty Equivalent Control for LQR is also analyzed with suboptimality gap (Mania et al., 2019). The other class is model-free approaches where policy is directly learned without estimating dynamics. Regarding discrete time LQR, in particular, Q-learning (Bradtke et al., 1994) and policy gradient (Fazel et al., 2018) were developed together with lots of mathematically machinery therein. but with little empirical demonstrations.

LQR with structured policy.   For continuous time LQR, many work in control literature (Wang & Davison, 1973; Sandell et al., 1978; Lau et al., 1972) has been studied for sparse and decentralized feedback gain, supported by theoretical demonstrations. Recently, (Wytock & Kolter, 2013; Lin et al., 2012) developed fast algorithms for a large-scale system that induce the sparse feedback gain utilizing Alternating Direction Method of Multiplier (ADMM), Iterative Shrinkage Thresholding Algorithm (ISTA), Newton method, etc. And, most of these only consider sparse feedback gain under continuous time LQR.

However, to the best of our knowledge, there are few algorithms for discrete time LQR, which take various structured (linear) policies such as sparse, low-rank, or proximity to some reference policy into account. Furthermore, regarding learning such a structured policy, there is no theoretical work that demonstrates a computational complexity, sample complexity, or the dependency of stable stepsize on algorithm parameters either, even in known-model setting (and model-free setting).

From standard LQR in Eq. (2), we restrict policy to linear class, i.e., $u_{t}=Kx_{t}$, and add a regularizer on the policy to induce the policy structure. We formally state a regularized LQR problem as

for a nonnegative parameter $\lambda\geq 0$. Here $f(K)$ is the (averaged) cost-to-go under policy $K$, and $r:\mathbf{R}^{n\times m}\rightarrow\mathbf{R}$ is a nonnegative convex regularizer inducing the structure of policy $K$.

Sensitivity to initial and fixed stepsize choice. Note that we use the initial stepsize $\eta=\mathcal{O}(1/\lambda)$ as a rule of thumb whereas the typical initial stepsize (Barzilai & Borwein, 1988; Wright et al., 2009; Park et al., 2020) does not depend on the regularization parameter $\lambda$. This stepsize choice is motivated by theoretical analysis (see Lemma 7 in Section 2.3) as well as empirical demonstration (see Fig. 1 in Section 4). This order of stepsize automatically scales well when experimenting over various $\lambda$s, alleviating iteration counts and leading to a faster algorithm in practice. It turns out that proximal gradient step is very sensitive to stepsizes, often leading to an unstable policy $K$ with $\rho(A+BK)\geq 1$ or requiring a large number of iteration counts to converge. Therefore, we utilize linesearch over fixed stepsize choice.

Linesearch.   We adopt a backtracking linesearch (See (Parikh et al., 2014)). Given $\eta_{i}$, $K^{i}$, $\nabla f(K^{i})$, and the potential next iterate $K^{i+1}$, it check if the following criterion (the stability and the decrease of the objective) is satisfied,

where $G_{\eta_{i}}(K)=\frac{1}{\eta_{i}}(K-\mathbf{prox}_{r,\lambda\eta_{i}}(K-\eta_{i}\nabla f(K)))$ and $\rho(\cdot)$ is the spectral radius. Otherwise, it shrinks the stepsize $\eta_{i}$ by a factor of $\beta<1$ and check it iteratively until Eq. (7) is satisfied.

Stabilizing sequence of policy.   We can start with a stable policy $K^{0}$, meaning $\rho(A+BK^{0})<1$. For example, under the standard assumptions on $A,B,Q,R$, Riccati recursion provides a stable policy, the solution of standard LQR in Eq. (2). Then, satisfying the linesearch criterion in Eq. (7) subsequently, the rest of the policies $\{K^{i}\}$ are a stabilizing sequence.

Computational cost.   The major cost incurs when solving the Lyapunov equations in the policy (and covariance) evaluation step. Note that if $A+BK$ is stable, so does $(A+BK)^{T}$ since they share the same eigenvalues. Under the stability, each Lyapunov equation in Eq. (4) has a unique solution with the computational cost $O(n^{3})$ (Jaimoukha & Kasenally, 1994; Li & White, 2002). Additionally, we can solve a sequence of Lyapunov equations with less cost, via using iterative methods with adopting the previous one (warm-start) or approximated one.

Recall that $\left\|\cdot\right\|$, $\left\|\cdot\right\|_{F}$, $\sigma_{\mathrm{min}}(\cdot)$ is $\ell_{2}$ matrix norm, Frobenius norm, and smallest singular value, as mentioned in Section 1.1. First, we start with (local) smoothness and strong convexity around a stable policy. Here we regard a policy $K$ is stable if $\rho(A+BK)<1$.

Assumptions.  $\rho(A+BK^{0})<1$, $\Sigma_{0}=\mathbf{E}x_{0}x_{0}^{T}\succ 0$, $\left\|K^{0}-K^{\star}\right\|\leq\Delta$, and $\left\|K^{\mathrm{ref}}-K^{\star}\right\|\leq\Delta$.

Then, we provide a proper stepsize that guarantees one iteration of proximal gradient step is still inside of the stable and (local) smooth ball $\mathcal{B}(K;\rho_{K})$.

Next proposition describes that next iterate policy has the decrease in function value and is stable under sufficiently small stepsize.

We also derive the bound on stepsize in Lemma 7, not dependent on iteration numbers.

Now we prove that the S-PI in Algorithm 1 converges to the stationary point linearly.

Note that (the global bound on) stepsize $\eta_{\mathrm{min}}$ is inversely proportional to $\lambda$, which motivates the initial stepsize in linesearch.

Note that, in model-free setting, model parameters $A,B,Q$ and $R$ cannot be directly accessed, which hinders the direct computation of $P$, $\Sigma$, and $\nabla f(K)$ accordingly. Instead, we adopt a smoothing procedure to estimate the gradient based on samples.

Model-free S-PI in Algorithm 3 consists of two steps: (1) policy evaluation step and (2) policy improvement step. In (perturbed) policy evaluation step, perturbation $U^{j}$ is uniformly drawn from the surface of the ball with radius $r$, $\mathbb{S}_{r}\subset\mathbf{R}^{n\times m}$. These data are used to estimate the gradient via a smoothing procedure for the policy improvement step. With this approximate gradient, proximal gradient subroutine tries to decrease the objective while inducing the structure of policy. Comparing to (known-model) S-PI in Algorithm 1, one important difference is its usage of a fixed stepsize $\eta$, rather than an adaptive stepsize from a backtracking linesearch that requires to access function value $f(K)=\mathbf{Tr}(\Sigma_{0}P)$ explicitly.

The outline of proof is as following: We first claim that for proper parameters (perturbation, horizon number, numbers of trajectory), the gradient estimate from the smoothing procedure is close enough to actual gradient with high probability. Next we demonstrate that approximate proximal gradient still converges linearly with high probability.

In these experiments, we use the unstable Laplacian system  (Recht, 2019).

Large Laplacian dynamics.   $A\in\mathbf{R}^{n\times n}$ where

$B=Q=I_{n}\in\mathbf{R}^{n\times n}$ and $R=1000\times I_{n}\in\mathbf{R}^{n\times n}$.

Synthetic system parameters.   For the Laplacian system, we regard $(n,m)=(3,3)$, $(n,m)=(20,20)$, and $(n,m)=(10^{3},10^{3})$ dimension as small, medium, and large size of system. In addition, we experiment with Lasso regularizer over various $\lambda=10^{-2}\sim 10^{6}$.

Based on our theoretical results (Lemma 7) and sensitivity experiments that empirically show the dependency of stepsize $\eta$ on $\lambda$ (in Fig. 1), we set the initial stepsize $\eta=\frac{1}{\lambda}$. For the backtracking linesearch, we set $\beta=\frac{1}{2}$ and the convergence tolerance $\epsilon_{\mathrm{tol}}=10^{-6}$. We start with an initial policy $K^{0}$ from Riccati recursion.

We denote the stationary point from S-PI (at each $\lambda>0$) as $K^{\star}$ and the solution of the (unregularized) LQR as $K^{lqr}$ .

Dependency of stepsize $\eta$ on $\lambda$.   Under the same problem but with different choices of weight $\lambda$, the largest fixed stepsize $\eta$ that demonstrates the sequence of stable systems, i.e., $A+BK^{i}<1$ actually varies, as Lemma 7 implies. Fig. 1 shows that the largest stepsize diminishes as $\lambda$ increases. This motivates the choice of the initial stepsize (in linesearch) to be inversely proportional to $\lambda$.

Convergence behavior and stepsize sensitivity.   S-PI with a linesearh strategy converges very fast, within 2-3 iterations for most of the Laplacian system with dimension $n=3\sim 10^{3}$ and weight $\lambda=10^{-2}\sim 10^{6}$. However, S-PI with fixed stepsize may perform with subtlety even though the fixed stepsize choice is common in typical optimization problems. In Fig. 2, S-PI with fixed stepsize for the small Laplacian system $(n,m)=(3,3)$ gets very close to optimal but begins to deviate after a certain amount of iterations. Note that it was performed over long iterations with small enough stepsize $\eta_{\mathrm{fixed}}=3\times 10^{-5}$. Please refer to supplementary materials for more detailed results over various stepsizes. This illustrates that a linesearch strategy can be essential in S-PI (or even for the existing policy gradient method for LQR (Fazel et al., 2018)), even though the convergence analysis was shown with some fixed stepsize (but difficult to compute in practice).

Trade off between LQR performance and structure $K$.   In Fig. 3 for a medium size Laplacian system, we show that as $\lambda$ increases, the LQR cost $f(K^{\star})$ increases whereas cardinality decreases (sparsitiy is improved). Note that the LQR performance barely changes (or is slightly worse) for $\lambda\leq 615$ but the sparsity is significantly improved by more than $50\%$. In Fig. 4, we show the sparsity pattern (location of non-zero elements) of the policy matrix with $\lambda=600$ and $\lambda=620$.

Scalability & runtime performance.   In Fig. 5, we report the runtime for Laplacian system of $n=10,\ldots,500$. Notably, it shows the scalability of S-PI, as it takes less than 2 minutes to solve a large system with $n=500$ dimensions, where we used a MacBook Air (with a 1.3 GHz Intel Core i5 CPU) for experiments. These results demonstrate its applicability to large-scale problems in practice.