A Block Coordinate Ascent Algorithm
for Mean-Variance Optimization

Reinforcement Learning (RL) [Sutton and Barto, 1998] is a class of learning problems in which an agent interacts with an unfamiliar, dynamic, and stochastic environment, where the agent’s goal is to optimize some measures of its long-term performance. This interaction is conventionally modeled as a Markov decision process (MDP), defined as the tuple $({\mathcal{S},\mathcal{A},P_{0},P_{ss^{\prime}}^{a},r,\gamma})$, where $\mathcal{S}$ and $\mathcal{A}$ are the sets of states and actions, $P_{0}$ is the initial state distribution, $P_{ss^{\prime}}^{a}$ is the transition kernel that specifies the probability of transition from state $s\in\mathcal{S}$ to state $s^{\prime}\in\mathcal{S}$ by taking action $a\in\mathcal{A}$, $r(s,a):\mathcal{S}\times\mathcal{A}\to\mathbb{R}$ is the reward function bounded by $R_{\max}$, and $0\leq\gamma<1$ is a discount factor. A parameterized stochastic policy $\pi_{\theta}(a|s):\mathcal{S}\times\mathcal{A}\to\left[{0,1}\right]$ is a probabilistic mapping from states to actions, where $\theta$ is the tunable parameter and $\pi_{\theta}(a|s)$ is a differentiable function w.r.t. $\theta$.

One commonly used performance measure for policies in episodic MDPs is the return or cumulative sum of rewards from the starting state, i.e., $R=\sum\nolimits_{k=1}^{{\tau}}{r({s_{k}},{a_{k}})}$, where $s_{1}\sim P_{0}$ and $\tau$ is the first passage time to the recurrent state $s^{*}$ [Puterman, 1994; Tamar et al., 2012], and thus, $\tau\coloneqq\min\{k>0\;|\;{s_{k}}={s^{*}}\}$. In risk-neutral MDPs, the algorithms aim at finding a near-optimal policy that maximizes the expected sum of rewards $J(\theta)\coloneqq\mathbb{E}_{\pi_{\theta}}[R]=\mathbb{E}_{\pi_{\theta}}\big[\sum_{k=1}^{{\tau}}{r({s_{k}},{a_{k}})}\big]$. We also define the square-return $M(\theta)\coloneqq{\mathbb{E}_{{\pi_{\theta}}}}[{R^{2}}]={\mathbb{E}_{{\pi_{\theta}}}}\Big[{\big(\sum_{k=1}^{\tau}{r({s_{k}},{a_{k}})}\big)^{2}}\Big]$. In the following, we sometimes drop the subscript $\pi_{\theta}$ to simplify the notation.

In risk-sensitive mean-variance optimization MDPs, the objective is often to maximize $J(\theta)$ with a variance constraint, i.e.,

where $\text{Var}_{\pi_{\theta}}(R)=M(\theta)-J^{2}(\theta)$ measures the variance of the return random variable $R$, and $\zeta>0$ is a given risk parameter [Tamar et al., 2012; Prashanth and Ghavamzadeh, 2013]. Using the Lagrangian relaxation procedure [Bertsekas, 1999], we can transform the optimization problem (1) to maximizing the following unconstrained objective function:

It is important to note that the mean-variance objective function is NP-hard in general [Mannor and Tsitsiklis, 2011]. The main reason for the hardness of this optimization problem is that although the variance satisfies a Bellman equation [Sobel, 1982], unfortunately, it lacks the monotonicity property of dynamic programming (DP), and thus, it is not clear how the related risk measures can be optimized by standard DP algorithms [Sobel, 1982].

The existing methods to maximize the objective function (3) are mostly based on stochastic approximation that often converge to an equilibrium point of an ordinary differential equation (ODE) [Borkar, 2008]. For example, Tamar et al., 2012 proposed a policy gradient algorithm, a two-time-scale stochastic approximation, to maximize (3) for a fixed value of $\lambda$ (they optimize over $\lambda$ by selecting its best value in a finite set), while the algorithm in Prashanth and Ghavamzadeh, 2013 to maximize (3) is actor-critic and is a three-time-scale stochastic approximation algorithm (the third time-scale optimizes over $\lambda$). These approaches suffer from certain drawbacks: 1) Most of the analyses of ODE-based methods are asymptotic, with no sample complexity analysis. 2) It is well-known that multi-time-scale approaches are sensitive to the choice of the stepsize schedules, which is a non-trivial burden in real-world problems. 3) The ODE approach does not allow extra penalty functions. Adding penalty functions can often strengthen the robustness of the algorithm, encourages sparsity and incorporates prior knowledge into the problem [Hastie et al., 2001].

Coordinate descent (CD)¹¹ 1 Note that since our problem is maximization, our proposed algorithms are block coordinate ascent. and the more general block coordinate descent (BCD) algorithms solve a minimization problem by iteratively updating variables along coordinate directions or coordinate hyperplanes [Wright, 2015]. At each iteration of BCD, the objective function is (approximately) minimized w.r.t. a coordinate or a block of coordinates by fixing the remaining ones, and thus, an easier lower-dimensional subproblem needs to be solved. A number of comprehensive studies on BCD have already been carried out, such as Luo and Tseng, 1992 and Nesterov, 2012 for convex problems, and Tseng, 2001, Xu and Yin, 2013, and Razaviyayn et al., 2013 for nonconvex cases (also see Wright, 2015 for a review paper). For stochastic problems with a block structure, Dang and Lan, 2015 proposed stochastic block mirror descent (SBMD) by combining BCD with stochastic mirror descent [Beck and Teboulle, 2003; Nemirovski et al., 2009]. Another line of research on this topic is block stochastic gradient coordinate descent (BSG) [Xu and Yin, 2015]. The key difference between SBMD and BSG is that at each iteration, SBMD randomly picks one block of variables to update, while BSG cyclically updates all block variables.

In this paper, we develop mean-variance optimization algorithms based on both nonconvex stochastic BSG and SBMD. Since it has been shown that the BSG-based methods usually have better empirical performance than their SBMD counterparts, the main algorithm we report, analyze, and evaluate in the paper is BSG-based. We report our SBMD-based algorithm in Appendix C and use it as a baseline in the experiments of Section 5. The finite-sample analysis of our BSG-based algorithm reported in Section 4 is novel because although there exists such analysis for convex stochastic BSG methods [Xu and Yin, 2015], we are not aware of similar results for their nonconvex version to the best our knowledge.

In this section, we describe why the vanilla stochastic gradient cannot be used to maximize $J_{\lambda}(\theta)$ defined in Eq. (3). Taking the gradient of $J_{\lambda}(\theta)$ w.r.t. $\theta$, we have

Computing $\nabla_{\theta}{J_{\lambda}}({\theta_{t}})$ in (5) involves computing three quantities: $\nabla_{\theta}J({\theta}),\nabla_{\theta}M({\theta})$, and $J({\theta})\nabla_{\theta}J({\theta})$. We can obtain unbiased estimates of $\nabla_{\theta}J(\theta)$ and $\nabla_{\theta}M(\theta)$ from a single trajectory generated by the policy $\pi_{\theta}$ using the likelihood ratio method [Williams, 1992], as $\nabla_{\theta}J({\theta})=\mathbb{E}[{R_{t}}\omega_{t}(\theta)]$ and $\nabla_{\theta}M({\theta})=\mathbb{E}[R_{t}^{2}\omega_{t}(\theta)]$. Note that $R_{t}$ is the cumulative reward of the $t$-th episode, i.e., $R_{t}=\sum\nolimits_{k=1}^{{\tau_{t}}}r_{k}$, which is possibly a nonconvex function, and $\omega_{t}(\theta)=\sum_{k=1}^{{\tau_{t}}}{\nabla_{\theta}\ln{\pi_{\theta}}({a_{k}}|{s_{k}})}$ is the likelihood ratio derivative. In the setting considered in the paper, an episode is the trajectory between two visits to the recurrent state $s^{*}$. For example, the $t$-th episode refers to the trajectory between the ($t$-1)-th and the $t$-th visits to $s^{*}$. We denote by $\tau_{t}$ the length of this episode.

However, it is not possible to compute an unbiased estimate of $J(\theta)\nabla_{\theta}J(\theta)$ without having access to a generative model of the environment that allows us to sample at least two next states $s^{\prime}$ for each state-action pair $(s,a)$. As also noted by Tamar et al., 2012 and Prashanth and Ghavamzadeh, 2013, computing an unbiased estimate of $J(\theta)\nabla_{\theta}J(\theta)$ requires double sampling (sampling from two different trajectories), and thus, cannot be done using a single trajectory. To circumvent the double-sampling problem, these papers proposed multi-time-scale stochastic approximation algorithms, the former a policy gradient algorithm and the latter an actor-critic algorithm that uses simultaneous perturbation methods [Bhatnagar et al., 2013]. However, as discussed in Section 2.1, multi-time-scale stochastic approximation approach suffers from several weaknesses such as no available finite-sample analysis and difficult-to-tune stepsize schedules. To overcome these weaknesses, we reformulate the mean-variance objective function and use it to present novel algorithms with in-depth analysis in the rest of the paper.

In this section, we present a new formulation for $J_{\lambda}(\theta)$ that is later used to derive our algorithms and do not suffer from the double-sampling problem in estimating $J(\theta)\nabla_{\theta}J(\theta)$. We begin with the following lemma.

This is a special case of the Lengendre-Fenchel duality [Boyd and Vandenberghe, 2004] that has been used in several recent RL papers (e.g., Liu et al., 2015; Du et al., 2017; Liu et al., 2018). Let ${F_{\lambda}}(\theta)\coloneqq{\left({J(\theta)+\frac{1}{{2\lambda}}}\right)^{2}}-M(\theta)$, which follows ${F_{\lambda}}(\theta)=\frac{{{J_{\lambda}}(\theta)}}{\lambda}+\frac{1}{{4{\lambda^{2}}}}-\zeta$. Since $\lambda>0$ is a constant, maximizing $J_{\lambda}(\theta)$ is equivalent to maximizing ${F_{\lambda}}(\theta)$. Using Lemma 1 with $z=J(\theta)+\frac{1}{2\lambda}$, we may reformulate ${F_{\lambda}}(\theta)$ as

Using (6), the maximization problem $\mathop{\max}_{\theta}{F_{\lambda}}(\theta)$ is equivalent to

Our optimization problem is now formulated as the standard nonconvex coordinate ascent problem (7). We use three stochastic solvers to solve (7): SBMD method [Dang and Lan, 2015], BSG method [Xu and Yin, 2015], and the vanilla stochastic gradient ascent (SGA) method [Nemirovski et al., 2009]. We report our BSG-based algorithm in Section 3.3 and leave the details of the SBMD and SGA based algorithms to Appendix C. In the following sections, we denote by $\beta_{t}^{\theta}$ and $\beta_{t}^{y}$ the stepsizes of $\theta$ and $y$, respectively, and by the subscripts $t$ and $k$ the episode and time-step numbers.

We now present our main algorithm that is based on a block coordinate update to maximize (7). Let ${g}_{t}^{\theta}$ and ${g}_{t}^{y}$ be block gradients and $\tilde{g}_{t}^{\theta}$ and $\tilde{g}_{t}^{y}$ be their sample-based estimations defined as

The block coordinate updates are

To obtain unbiased estimates of $g^{y}_{t}$ and $g^{\theta}_{t}$, we shall update $y$ (to obtain $y_{t+1}$) prior to computing $g^{\theta}_{t}$ at each iteration. Now it is ready to introduce the Mean-Variance Policy Gradient (MVP) Algorithm 1.

Before presenting our theoretical analysis, we first introduce the assumptions needed for these results.

Assumption 1 is standard in nonconvex coordinate descent algorithms [Xu and Yin, 2015; Dang and Lan, 2015]. We also need the following assumption that is standard in the policy gradient literature.

In practice, we can choose either Option I with the result of the final iteration as output or Option II with the result of a randomly selected iteration as output. In what follows in this section, we report an asymptotic convergence analysis of MVP with Option I, and in Section 4, we derive a finite-sample analysis of MVP with Option II.

The proof of Theorem 1 follows from the analysis in Xu and Yin, 2013. Due to space constraint, we report it in Appendix A.

Algorithm 1 is a special case of nonconvex block stochastic gradient (BSG) methods. To the best of our knowledge, no finite-sample analysis has been reported for this class of algorithms. Motivated by the recent papers by Nemirovski et al., 2009, Ghadimi and Lan, 2013, Xu and Yin, 2015, and Dang and Lan, 2015, in Section 4, we provide a finite-sample analysis for general nonconvex block stochastic gradient methods and apply it to Algorithm 1 with Option II.

In this section, we provide a finite-sample analysis of the general nonconvex block stochastic gradient (BSG) method, where the problem formulation is given by

$\xi$ is a random vector, and $F(\cdot,\xi):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is continuously differentiable and possibly nonconvex for every $\xi$. The variable $x\in\mathbb{R}^{n}$ can be partitioned into $b$ disjoint blocks as $x=\{x^{1},x^{2},\dotsc,x^{b}\}$, where $x^{i}\in\mathbb{R}^{n_{i}}$ denotes the $i$-th block of variables, and $\sum_{i=1}^{b}n_{i}=n$. For simplicity, we use $x^{<i}$ for $(x_{i},\dots,x_{i-1})$, and $x^{\leq i}$,$x^{>i}$, and $x^{\geq i}$ are defined correspondingly. We also use $\nabla_{x^{i}}$ to denote $\frac{\partial}{\partial x^{i}}$ for the partial gradient with respect to $x^{i}$. $\Xi_{t}$ is the sample set generated at $t$-th iteration, and $\mathbf{\Xi}_{[t]}=(\Xi_{1},\dots,\Xi_{t})$ denotes the history of sample sets from the first through $t$-th iteration. $\{\beta_{t}^{i}:i=1,\cdots,b\}_{t=1}^{\infty}$ are denoted as the stepsizes. Also, let $\beta_{t}^{\max}=\max_{i}\beta_{t}^{i}$, and $\beta_{t}^{\min}=\min_{i}\beta_{t}^{i}$. Similar to Algorithm 1, the BSG algorithm cyclically updates all blocks of variables in each iteration, and the detailed algorithm for BSG method is presented in Appendix B.

Without loss of generality, we assume a fixed update order in the BSG algorithm. Let $\Xi_{t}=\{\xi_{t,1},\dotsc,\xi_{t,m_{t}}\}$ be the samples in the $t$-th iteration with size $m_{t}\geq 1$. Therefore, the stochastic partial gradient is computed as $\tilde{g}_{t}^{i}=\frac{1}{m_{t}}\sum_{l=1}^{m_{t}}\nabla_{x^{i}}F(x_{t+1}^{<i},x_{t}^{\geq i};\xi_{t,l}).$ Similar to Section 3, we define $g_{t}^{i}=\nabla_{x^{i}}f(x_{t+1}^{<i},x_{t}^{\geq i})$, and the approximation error as $\Delta_{t}^{i}=\tilde{g}_{t}^{i}-g_{t}^{i}$. We assume that the objective function $f$ is bounded and Lipschitz smooth, i.e., there exists a positive Lipschitz constant $L>0$ such that $\|\nabla_{x^{i}}f(x)-\nabla_{x^{i}}f(y)\|_{2}\leq L\|x-y\|_{2}$, $\forall i\in\{1,\dots,b\}$ and $\forall x,y\in\mathbb{R}^{n}$. Each block gradient of $f$ is also bounded, i.e., there exist a positive constant $G$ such that $\|\nabla_{x^{i}}f(x)\|_{2}\leq G$, for any $i\in\{1,\dots,b\}$ and any $x\in\mathbb{R}^{n}$. We also need Assumption 1 for all block variables, i.e., $\mathbb{E}[\|\Delta_{t}^{i}\|_{2}]\leq\sigma$, for any $i$ and $t$. Then we have the following lemma.

The proof of Lemma 2 is in Appendix B. It should be noted that in practice, it is natural to take the final iteration’s result as the output as in Algorithm 1. However, a standard strategy for analyzing nonconvex optimization methods is to pick up one previous iteration’s result randomly according to a discrete probability distribution over $\{1,2,\dotsc,N\}$ [Nemirovski et al., 2009; Ghadimi and Lan, 2013; Dang and Lan, 2015]. Similarly, our finite-sample analysis is based on the strategy that randomly pick up $\bar{x}_{N}=x_{z}$ according to

Now we provide the finite-sample analysis result for the general nonconvex BSG algorithm as in [Xu and Yin, 2015].

As a special case, we discuss the convergence rate with constant stepsizes $\mathcal{O}(1/\sqrt{N})$ in Corollary 1, which implies that the sample complexity $N=\mathcal{O}(1/\varepsilon^{2})$ in order to find $\varepsilon$-stationary solution of problem (18).

We present the major theoretical results of this paper, i.e., the finite-sample analysis of Algorithm 1 with Option II. The proof of Theorem 3 is in Appendix A.

The portfolio domain Tamar et al., 2012 is composed of the liquid and non-liquid assets. A liquid asset has a fixed interest rate $r_{l}$ and can be sold at any time-step $k\leq{\tau}$. A non-liquid asset can be sold only after a fixed period of $W$ time-steps with a time-dependent interest rate $r_{\text{nl}}(k)$, which can take either $r_{\text{nl}}^{\text{low}}$ or $r_{\text{nl}}^{\text{high}}$, and the transition follows a switching probability $p_{\text{switch}}$. The non-liquid asset also suffers a default risk (i.e., not being paid) with a probability $p_{\text{risk}}$. All investments are in liquid assets at the initial time-step $k=0$. At the $k$-th step, the state is denoted by $x(k)\in\mathbb{R}^{W+2}$, where $x_{1}\in[0,1]$ is the portion of the investment in liquid assets, $x_{2},\cdots,x_{W+1}\in[0,1]$ is the portion in non-liquid assets with time to maturity of $1,\cdots,W$ time-steps, respectively, and $x_{W+2}(k)=r_{\text{nl}}(k)-\mathbb{E}[r_{\text{nl}}(k)]$. The investor can choose to invest a fixed portion $\eta$ ($0<\eta<1$) of his total available cash in the non-liquid asset or do nothing. More details about this domain can be found in Tamar et al., 2012. Figure 1(a) shows the results of the algorithms. PG has a large variance and the Tamar’s method has the lowest mean return. The results indicate that MVP yields a higher mean return with less variance compared to the competing algorithms.

An American-style option Tamar et al., 2014 is a contract that gives the buyer the right to buy or sell the asset at a strike price $W$ at or before the maturity time ${\tau}$. The initial price of the option is $x_{0}$, and the buyer has bought a put option with the strike price $W_{\text{put}}<x_{0}$ and a call option with the strike price $W_{\text{call}}>x_{0}$. At the $k$-th step ($k\leq\tau$), the state is $\{x_{k},k\}$, where $x_{k}$ is the current price of the option. The action $a_{k}$ is either executing the option or holding it. $x_{k+1}$ is $f_{u}x_{k}$ w.p. $p$ and $f_{d}x_{k}$ w.p. $1-p$, where $f_{u}$ and $f_{d}$ are constants. The reward is $0$ unless an option is executed and the reward for executing an option is ${r_{k}}=\max(0,{W_{{\rm{put}}}}-{x_{k}})+\max(0,{x_{k}}-{W_{{\rm{call}}}})$. More details about this domain can be found in Tamar et al., 2014. Figure 1(b) shows the performance of the algorithms. The results suggest that MVP can yield a higher mean return with less variance compared to the other algorithms.

The optimal stopping problem [Chow and Ghavamzadeh, 2014; Chow et al., 2018] is a continuous state domain. At the $k$-th time-step ($k\leq\tau$, $\tau$ is the stopping time), the state is $\{x_{k},k\}$, where $x_{k}$ is the cost. The buyer decide either to accept the present cost or wait. If the buyer accepts or when $k=T$, the system reaches a terminal state and the cost $x_{k}$ is received, otherwise, the buyer receives the cost $p_{h}$ and the new state is $\{x_{k+1},k+1\}$, where $x_{k+1}$ is $f_{u}x_{k}$ w.p. $p$ and $f_{d}x_{k}$ w.p. $1-p$ ($f_{u}>1$ and $f_{d}<1$ are constants). More details about this domain can be found in Chow and Ghavamzadeh, 2014. Figure 1(c) shows the performance of the algorithms. The results indicate that MVP is able to yield much less variance without affecting its mean return. We also summarize the performance of these algorithms in all three risk-sensitive domains as Table 1, where Std is short for Standard Deviation.