\section{Introduction}
The standard formulation of reinforcement learning (RL), which aims to find the optimal policy to optimize the cumulative reward, has been well studied in the recent years. Compared with the model-based algorithms, model-free algorithms do not require the estimation of the transition dynamics and  can be extended to the continuous space. Value function based  algorithms such as Q-learning \cite{watkins1992q,jin2018}, SARSA \cite{rummery1994line}, Temporal Difference (TD) \cite{sutton1988learning} and policy based algorithms such as policy gradient \cite{sutton2000} and natural policy gradient \cite{Kakade2002} have been proposed based on the Bellman Equation, which is a result of the additive structure for the standard RL.

However, many applications require more general non-linear reward functions. As an example, risk-sensitive  objectives have been considered in \cite{mihatsch2002risk}. \cite{Hazan2019} studies the problem of maximizing the entropy of state-action distribution. Further, many realistic applications have multiple objectives, e.g., capacity and power usage in the communication system \cite{aggarwal2017joint}, latency and energy consumption in queueing systems \cite{badita2020optimal}, efficiency and safety in robotic systems \cite{nishimura2020l2b}. 



In this paper, we consider a setting that jointly optimizes a general concave function of the cumulative reward from multiple objectives. With this definition, a non-linear concave function of single objective becomes a special case. Further, fair allocation of resources among multiple users require a non-linear function of the rewards to each user (which correspond to the multiple objectives) \cite{Tian2010}, and is thus a special case of this formulation. Such a setup was first considered in \cite{mridulmb}, where a model-based algorithm was proposed for the problem with provable regret guarantees. However, guarantees for model-free algorithm have not been studied to the best of our knowledge, which we focus on. 




We note that the non-linear objective function looses the additive structure, and thus the Bellman's Equation does not work anymore in this setting \cite{mridulmb,Mengdi2020}. Thus, the value function based algorithm do not directly work in this setup. This paper considers a policy-gradient approach and aim to show the global convergence of such policies. Recently, the authors of \cite{Mengdi2020,Mengdi2021} considered the problem for a single-objective over finite state-action space. However, such a problem is open for continuous state action spaces, and for multiple objectives, which is the focus of this paper. In this paper, we consider a  fundamental policy based algorithm, the vanilla policy gradient, and show the global convergence of this policy based on an efficient estimator of the gradient proposed in this paper. 

We note that in standard reinforcement learning, Policy Gradient Theorem \cite{sutton2000} is used to propose an unbiased gradient estimator such as REINFORCE. However, such an approach can not directly give an unbiased estimator in our setting due to the presence of non-linear function (See Lemma \ref{lem_biased_est}). In this paper, we provide a biased estimator for the policy gradient. This biased estimator is then used to prove the global convergence of the policy gradient algorithm. 


Our contribution can be summarized as follows.
\begin{itemize}[leftmargin=*]
	\item We propose a general biased gradient estimator, which can be applied to both tabular and continuous state-action spaces, and prove that the bias of the estimator decays at order $\mathcal{O}(1/\sqrt{n})$, where $n$ is the number of trajectories sampled (See Remark \ref{rem_bias_decay}).
	\item We prove the policy gradient algorithm with the proposed estimator converges to the global optimal with error $\epsilon$ using  $\mathcal{O}(\frac{M^4\sigma^2}{(1-\gamma)^8\epsilon^4})$ samples, where $M$ is the number of agent, $\sigma^2$ is the variance defined in Assumption \ref{ass_bounded_var} and $\gamma$ is the discount factor. As compared to the number of samples  for standard RL with policy gradient algorithm \cite{Yanli2020}, our result has the same dependence on $\epsilon$. %
\end{itemize} 

Further, even for the case when there is a non-linear function of a single objective, the approach and results are novel, and have not been considered in the prior works for continuous state-action spaces. 
\section{Policy Gradient Method for Joint Optimization}
	Policy gradient algorithm aims to update the parameter with the iteration
	\begin{equation}
		\theta^{k+1}=\theta^k+\eta\nabla_\theta f(\bm{J}^{\pi_{\theta^k}})
	\end{equation}
	where $\eta$ is the step size. However, it is impossible to compute the true gradient because the transition dynamics is unknown in practice. Thus, an estimator for the true gradient is necessary. In this section, we first give the form of the true gradient. Then, we propose a biased estimator and bound the bias. The policy-gradient algorithm is also formally described based on the estimator. Finally, we analyze some properties of the objective function, which will be used in the proof of the main result.
	\subsection{Computation of the Gradient}
	Starting from the Chain Rule
	\begin{equation}
		\nabla_{\theta}f(\bm{J}^{\pi_\theta})=\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_\theta})}{\partial J_m^{\pi_\theta}}\nabla_\theta J_m^{\pi_\theta}
	\end{equation}
	%
	Define $\tau=(s_0,a_1,s_1,a_1,s_2,a_2\cdots)$ as a trajectory, whose distribution induced by policy $\pi_\theta$ is $p(\tau\vert \theta)$ that can be expressed as
	\begin{equation}\label{eq:traj_dis}
		p(\tau\vert \theta)=\rho(s_0)\prod_{t=0}^{\infty}\pi_\theta(a_t\vert s_t)P(s_{t+1}\vert s_t,a_t)
	\end{equation}
	Define $R_m(\tau)=\sum_{t=0}^{\infty}\gamma^tr_m(s_t,a_t)$ as the cumulative reward for $m^{th}$ agent following the trajectory $\tau$. Then, the expected return $J_m^\pi(\theta)$ can also be expressed as
	\begin{equation*}
		J_m^{\pi_\theta}=\mathbf{E}_{\tau\sim p(\tau\vert\theta)}[R_m(\tau)]
	\end{equation*}
	and the gradient can be calculated as
	\begin{equation*}
		\begin{aligned}
			\nabla_\theta J_m^\pi(\theta)&=\int_\tau R_m(\tau)p(\tau\vert \theta)d\tau =\int_{\tau}R_m(\tau)\frac{\nabla_\theta p(\tau\vert \theta)}{p(\tau\vert \theta)}p(\tau\vert \theta)d\tau=\mathbf{E}_{\tau\sim p(\tau\vert \theta)}\big[\nabla_\theta \log p(\tau\vert\theta)R_m(\tau)\big]
		\end{aligned}
	\end{equation*}
	Notice that $\nabla_\theta \log p(\tau\vert \theta)$ is independent of the transition dynamics because 
	\begin{equation*}
		\nabla_\theta \log p(\tau\vert \theta)=\nabla_\theta\bigg[\log \rho(s_0)+\sum_{t=0}^{\infty}\big[\log\pi_\theta(a_t\vert s_t)+\log P(s_{t+1}\vert s_t,a_t)\big]\bigg]=\sum_{t=0}^\infty\nabla_\theta \log\pi_\theta(a_t\vert s_t)
	\end{equation*}
	and thus the gradient for the objective function is
	\begin{equation}\label{eq:gradient_origin}
		\nabla_\theta f(\bm{J}^{\pi_\theta})=\mathbf{E}_{\tau\sim p(\tau\vert \theta)}\bigg[\bigg(\sum_{t=0}^\infty \nabla_\theta \log \pi_\theta(a_t\vert s_t)\bigg)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{t=0}^\infty \gamma^tr_m(s_t,a_t)\big)\bigg)\bigg]
	\end{equation}
	Notice that removing the past reward from the return doesn't change the expectation value \cite{Peters2008}. Thus, we can rewrite Eq. \eqref{eq:gradient_origin} as
	\begin{equation}\label{eq:gradient_simplified}
		\nabla_\theta f(\bm{J}^{\pi_\theta})=\mathbf{E}_{\tau\sim p(\tau\vert \theta)}\bigg[\sum_{t=0}^\infty \nabla_\theta \log(\pi_\theta(a_t\vert s_t))\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^\infty \gamma^hr_m(s_h,a_h)\big)\bigg)\bigg]
	\end{equation}
	
	
	\if 0
	\begin{lemma}\label{lem_grad_simplified}
		Eq. \eqref{eq:gradient_origin} and Eq. \eqref{eq:gradient_simplified} are the same
	\end{lemma}
	\begin{proof}
		Notice that the difference between Eq. \eqref{eq:gradient_origin} and \eqref{eq:gradient_simplified} can be expressed as below
		\begin{equation}
			\mathbf{E}_{\tau\sim p(\tau\vert \theta)}\bigg[\sum_{t=0}^\infty \nabla_\theta \log(\pi_\theta(a_t\vert s_t))\bigg(\sum_{m=1}^{M}\big(\frac{1}{N_2}\sum_{j=1}^{N_2}\frac{\partial f}{\partial J_m^\pi}\big)\big(\sum_{h=0}^{t-1} \gamma^hr_m(s_h,a_h)\big)\bigg)\bigg]
		\end{equation}
		Thus, it is sufficient to prove the above equation is $\bm{0}$. Divide the trajectory $\tau$ into two parts $\tau_{1}=(s_0,a_0,\cdots,s_{t-1},a_{t-1})$ and $\tau_{2}=(s_t,a_t,\cdots)$. Then, by law of total expectation,
		\begin{equation}
			\begin{aligned}
			&\mathbf{E}_{\tau\sim p(\tau\vert \theta)}\bigg[\sum_{t=0}^\infty \nabla_\theta \log(\pi_\theta(a_t\vert s_t))\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\big)\big(\sum_{h=0}^{t-1} \gamma^hr_m(s_h,a_h)\big)\bigg)\bigg]\\
			&=\sum_{t=0}^\infty\mathbf{E}_{\tau_1}\bigg\{\mathbf{E}_{\tau_2}\bigg[\nabla_\theta \log(\pi_\theta(a_t\vert s_t))\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\big)\big(\sum_{h=0}^{t-1} \gamma^hr_m(s_h,a_h)\big)\bigg)\bigg]\bigg|\tau_1\bigg\}\\
			&=\sum_{t=0}^\infty\mathbf{E}_{\tau_1}\bigg\{\mathbf{E}_{\tau_2}\bigg[\nabla_\theta \log(\pi_\theta(a_t\vert s_t))\bigg]\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\big)\big(\sum_{h=0}^{t-1} \gamma^hr_m(s_h,a_h)\big)\bigg)\bigg|\tau_1\bigg\}\\
			&=\sum_{t=0}^\infty\mathbf{E}_{\tau_1}\bigg\{\bm{0}\cdot\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\big)\big(\sum_{h=0}^{t-1} \gamma^hr_m(s_h,a_h)\big)\bigg)\bigg|\tau_1\bigg\}=\bm{0}
			\end{aligned}
		\end{equation}
	\end{proof}
	\fi
	
	
	\subsection{Proposed Estimator}
	From Eq. \eqref{eq:gradient_simplified}, an estimator for $\nabla_\theta f(\bm{J}^{\pi_\theta})$ can be directly derived as
	\begin{equation}\label{eq:est_origin}
		g(\tau_i,\tau_{j=1:N_2}\vert \theta)=\sum_{t=0}^\infty \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_m^{\pi}}\big)\big(\sum_{h=t}^\infty \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)
	\end{equation}
	where
	\begin{equation}
		\hat{J}_{m}^{\pi}=\frac{1}{N_2}\sum_{j=1}^{N_2}\sum_{t=0}^{\infty}\gamma^tr_m(s_t^j,a_t^j)
	\end{equation}
	and $N_2$ is the number of trajectories of $\tau_j$ that we need to sample to estimate $\frac{\partial f}{\partial J_m^\pi}$. Notice that the trajectories $\tau_i=(s_0^i,a_0^i,s_1^i,a_1^i,\cdots)$ and $\tau_j=(s_0^j,a_0^j,s_1^j,a_1^j,\cdots)$ are sampled independently from the distribution $p(\tau\vert \theta)$. However, notice that the proposed estimator is not unbiased due to the concavity of the function $f$ (See Lemma \ref{lem_biased_est} in Appendix \ref{sec_app_bias}). Moreover, the estimator in Eq. \eqref{eq:est_origin} is unachievable because it requires a sum over infinite range of $t$. Thus, we define a truncated version of Eq. \eqref{eq:est_origin} as
	
	\begin{equation}\label{eq:est_truncated}
		g(\tau_i^H,\tau_{j=1:N_2}^H\vert \theta)=\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_{m,H}^{\pi}}\big)\big(\sum_{h=t}^{H-1} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)
	\end{equation}
	where
	\begin{equation}\label{eq:def_Jm}
		\hat{J}_{m,H}^{\pi}=\frac{1}{N_2}\sum_{j=1}^{N_2}\sum_{t=0}^{H-1}\gamma^tr_m(s_t^j,a_t^j)
	\end{equation}
	In the remaining part of this paper, we denote $g(\tau_i^H,\tau_{j=1:N_2}^H\vert \theta)$ as $g(\tau_i^H,\tau_j^H\vert \theta)$ for simplicity. With this truncated estimator, the proposed algorithm is in Algorithm \ref{alg:PG}. In each iteration of policy gradient ascent, $N_2$ trajectories are sampled in line 3 and used to estimate the value function for each agent. Line 4 samples another $N_1$ trajectories independent of $N_2$ and uses Eq. \eqref{eq:est_truncated} to calculate the gradient estimator. Line 5 and 6 perform one-step gradient descent using the gradient estimator.
	 	
	\begin{algorithm*}[tb]
		\caption{Policy Gradient for Joint Optimization of Multi-Objective RL}
		\label{alg:PG}
		\begin{algorithmic}[1]
			\State Initialize $\theta^0$ and step size $\eta=\frac{1}{4L_J}$
			\For{episode $k=0,...,K-1$} 
			\State Sample $N_2$ trajectories $\tau_j$ under policy $\theta^k$ of length $H$ and compute $\hat{J}_{m,H}^{\pi}$ by Eq. \eqref{eq:def_Jm}.
			\State Sample $N_1$ trajectories $\tau_i$ under policy $\theta^k$ of length $H$ and for each trajectory compute the gradient estimator $g(\tau_i^H,\tau_j^H\vert \theta^k)$ by Eq. \eqref{eq:est_truncated}
			\State Compute the gradient update direction $\omega^k=\frac{1}{N_1}\sum_{i=1}^{N_1}g(\tau_i^H\vert \theta^k)$
			\State Update the parameter $\theta^{k+1}=\theta^k+\eta\omega^k$
			\EndFor 
		\end{algorithmic}
	\end{algorithm*}
	
	\subsection{Bounding the Bias of the Truncated Estimator}
	To bound the bias of the proposed truncated estimator, we define three auxiliary functions.
	\begin{equation}
		\tilde{g}(\tau_i,\tau_j\vert \theta)=\sum_{t=0}^{\infty} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\big)\big(\sum_{h=t}^{\infty} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)
	\end{equation}
	\begin{equation}\label{eq:finite_auxi}
		\tilde{g}(\tau_i^H,\tau_j\vert \theta)=\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\big)\big(\sum_{h=t}^{H-1} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)
	\end{equation}
	\begin{equation}\label{eq:truncated_auxi}
		\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)=\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=J_{m,H}^{\pi}}\big)\big(\sum_{h=t}^{H-1} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)
	\end{equation}
	where $J_{m,H}^{\pi}=\mathbf{E}\bigg[\sum_{t=0}^{H-1}\gamma^tr_m(s_t,a_t)\bigg]$.
	
	It should be noticed that Eq. \eqref{eq:finite_auxi} and \eqref{eq:truncated_auxi} are different because the value function used in the partial derivatives are truncated in \eqref{eq:truncated_auxi} but not in \eqref{eq:finite_auxi}.	Moreover, Eq. \eqref{eq:truncated_auxi} and the proposed estimator in Eq. \eqref{eq:est_truncated} are also different because \eqref{eq:est_truncated} uses the empirical value for trajectories $\tau_j^H$ while Eq. \eqref{eq:truncated_auxi} uses the expected value. We note that $\tilde{g}(\tau_i,\tau_j\vert \theta)$ is an unbiased estimator for $\nabla_\theta f(\bm{J}^{\pi_\theta})$. Thus, the bias of the truncated estimator Eq. \eqref{eq:est_truncated} can be decomposed as
	\begin{equation}
		\begin{split}
		\mathbf{E}[g(\tau_i^H,\tau_j^H\vert \theta)]-\nabla_\theta f(\bm{J}^{\pi_\theta})&=\mathbf{E}\underbrace{[g(\tau_i^H,\tau_j^H\vert \theta)-\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)]}_{(I)}\\
		&+\mathbf{E}\underbrace{[\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)-\tilde{g}(\tau_i^H,\tau_j\vert \theta)]}_{(II)}+\mathbf{E}\underbrace{[\tilde{g}(\tau_i^H,\tau_j\vert \theta)-\tilde{g}(\tau_i,\tau_j\vert \theta)]}_{(III)}
		\end{split}
	\end{equation}
	which means the bias includes three parts: (I) denotes the bias coming from the finite samples of trajectories $\tau_j$. (II) and (III) denote the bias due to the truncation of trajectories $\tau_j$ and $\tau_i$, respectively. In the following, we give three lemmas to bound each of them. The detailed proofs are provided in Appendix \ref{sec_app_bound_bias}.
	\begin{lemma}\label{lem_bound_bias1}
		For any $\epsilon'>0$ and $p\in(0,1)$, with probability at least $1-p$, if the number of samples for $\tau_j$ satisfies,
		\begin{equation}
			N_2\geq \frac{M(1-\gamma^H)^2}{2(1-\gamma)^2\epsilon'^2}\log(\frac{2MH}{p})
		\end{equation} 
		then for each trajectory $\tau_i$, the first part of bias for the proposed truncated estimator, Eq. \eqref{eq:est_truncated}, is bounded by
		\begin{equation}
			\Vert g(\tau_i^H,\tau_j^H\vert \theta)-\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)\Vert\leq MGL_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^2}\epsilon'
		\end{equation} 
	\end{lemma}
	\begin{lemma}\label{lem_bound_bias2} 
	For each trajectory $\tau_i$, the second part of bias for the proposed truncated estimator, Eq. \eqref{eq:est_truncated}, is bounded by
		\begin{equation}
			\Vert\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)-\tilde{g}(\tau_i^H,\tau_j\vert\theta)\Vert\leq M^{3/2}GL_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^3}\gamma^H
		\end{equation}
	\end{lemma}
	\begin{lemma}\label{lem_bound_bias3} 
		For each trajectory $\tau_i$, the third part of bias for the proposed truncated estimator, Eq. \eqref{eq:est_truncated}, is bounded by
		\begin{equation}
			\Vert\tilde{g}(\tau_i^H,\tau_j\vert \theta)-\tilde{g}(\tau_i,\tau_j\vert\theta)\Vert\leq MGC\frac{\gamma^H(1+H(1-\gamma))}{(1-\gamma)^2}
		\end{equation}
	\end{lemma}
	\begin{remark}\label{rem_bias_decay}
		Combining the Lemmas \ref{lem_bound_bias1}, \ref{lem_bound_bias2}, and \ref{lem_bound_bias3}, it is found that if the length of sampled trajectories is long enough, the bias of the proposed estimator decays as $\mathcal{O}(\frac{1}{\sqrt{N_2}})$.
	\end{remark}
	\subsection{Properties of the Objective Function}
	Similar to the truncated estimator, we define a truncated version for the objective function as follows
	\begin{equation*}
		f(\bm{J}_H^{\pi_\theta})=f(\mathbf{E}[\sum_{t=0}^{H-1}\gamma^tr_1(s_t,a_t)],\cdots,\mathbf{E}[\sum_{t=0}^{H-1}\gamma^tr_M(s_t,a_t)])
	\end{equation*}
	In this subsection, we will give some properties of $f(\bm{J}^{\pi_\theta})$ and $f(\bm{J}_H^{\pi_\theta})$. The detailed proofs are provided in  Appendix \ref{sec_app_property}.
	The following lemma shows the smoothness property for $f(\bm{J}^{\pi_\theta})$ and $f(\bm{J}_H^{\pi_\theta})$.
	\begin{lemma}\label{lem_smooth}
		Both the objective function $f(\bm{J}^{\pi_\theta})$ and the truncated version $f(\bm{J}_H^{\pi_\theta})$ are $L_J$-smooth w.r.t. $\theta$, where
		\begin{equation*}
			L_J=\frac{MCB}{(1-\gamma)^2}
		\end{equation*}
	\end{lemma}
	
	It is reasonable to expect that the truncated objective function and the original one can be arbitrary close when the length of horizon is long enough, and the next lemma bounds the gap between original and truncated objective function.
	\begin{lemma}\label{lem_bound_truncated}
		The difference between the gradient of objective function and that of truncated version is bounded by
		\begin{equation}
			\Vert \nabla_\theta f(\bm{J}^{\pi_\theta})-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert\leq \frac{MG\gamma^H}{(1-\gamma)^2}\bigg[ \sqrt{M}L_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{1-\gamma}+C[1+H(1-\gamma)]\bigg]
		\end{equation}
	\end{lemma}
	
	In order to introduce the following result, it is helpful to define the state visitation measure
	\begin{equation}
		d_{\rho}^{\pi}:=(1-\gamma)\mathbf{E}_{s_0\sim\rho}\bigg[\sum_{t=0}^{\infty}\gamma^t\mathbf{Pr}^{\pi}(s_t=s\vert s_0)\bigg]
	\end{equation}
	where $\mathbf{Pr}^{\pi}(s_t=s\vert s_0)$ denotes the probability that $s_t=s$ with policy $\pi$ starting from $s_0$. In the theoretical analysis of policy gradient for standard reinforcement learning, one key result is the performance difference lemma. In the multi-objective setting, a similar performance lemma is  derived as follows.
	\begin{lemma}\label{lem_performance_diff}
		The difference in the performance for  any policies $\pi_\theta$ and $\pi_{\theta'}$is bounded as follows
		\begin{equation}
			(1-\gamma)[f(\bm{J}^{\pi_\theta})-f(\bm{J}^{\pi_{\theta'}})]\leq \sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_{\theta'}})}{\partial J_m^{\pi_{\theta'}}}\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta(\cdot\vert s)}\big[A_m^{\pi_{\theta'}}(s,a)\big]
		\end{equation}
	\end{lemma}
	
\section{Related Work}
{\bf Policy Gradient with Cumulative Return: } As the core result for policy based algorithms, Policy Gradient Theorem \cite{sutton2000} provides a method to obtain the gradient ascent direction for the standard reinforcement learning with the policy parameterization. However, in general, the objective in the reinforcement learning is non-convex with respective to the parameters \cite{Alekh2020}. Thus, the research on policy gradient algorithm focuses on the first order stationary point guarantees for a long time \cite{Tianbing2017,Xu2019,xu2020}. Recently, there is a line of interest on the global convergence result for reinforcement learning. \cite{Kaiqing2019} utilizes the idea of escaping saddle points in the policy gradient and shows the convergence to the second order stationary points, which gives the locally optimal. \cite{Alekh2020} provides provable global convergence result for direct parameterization and softmax parameterization in the tabular case. For the restrictive parameterization, they propose a variant of NPG, Q-NPG and analyze the global convergence result with the function approximation error for both NPG and Q-NPG. \cite{Mei2020} improves the convergence rate for policy gradient with softmax parameterization from $\mathcal{O}(1/\sqrt{t})$ to $\mathcal{O}(1/t)$ and shows a significantly faster linear convergence rate $\mathcal{O}(\exp(-t))$ for the entropy regularized policy gradient. With actor-critic method \cite{konda2000actor}, \cite{Lingxiao2019} establishes the global optimal result for neural policy gradient method. \cite{bhandari2020global} identifies the structure properties which shows that there are no sub-optimal stationary points for reinforcement learning. \cite{Yanli2020} proposes a general framework of the analysis for policy gradient type of algorithms and gives the sample complexity for PG, NPG and the variance reduced version of them. However, all of the above research have been done on the standard reinforcement learning, where the objective function is the direct summation of the reward. This paper focuses on a joint optimization of multi-objective problem, where multiple objectives are combined with a concave function. 

{\bf Policy Gradient with General Objective Function: } Even though standard reinforcement learning has been widely studied, there are few results on the policy gradient algorithm with a general objective function. Some special examples are variance-penalty \cite{YingHuang1994} and maximizing entropy \cite{Hazan2019}. Very recently, \cite{Mengdi2020,Mengdi2021} study the global convergence result of the policy gradient with general utilities. They consider the setting that the objective is a concave function of the state-action occupancy measure, which is similar to our setting. By the method of convex conjugate, \cite{Mengdi2020} proposed a variational policy gradient theorem to obtain the gradient ascent direction and gives the global convergences result of PG with general utilities. Despite enjoying a rate of $\mathcal{O}(1/t)$ in terms of iterations, their algorithm requires an additional saddle point problem to fulfill the gradient update and thus introduce extra computation complexity. \cite{Mengdi2021} further proposes the SIVR-PG algorithm and improves the convergence rate in the same setting. However, the SIVR-PG algorithm requires the estimation of state-action occupancy measure, which means that the algorithm can only be applied to the tabular setting.  We note that  our method does not have such limitation and thus can be applied even if the state and action space is large or continuous. Finally, note that \cite{Mengdi2020,Mengdi2021} improve the previous convergence rate for policy gradient by exploring the hidden convexity of the proposed problem. However, in order to utilize such convexity, they require the assumption that the inverse mapping of visitation measure $\lambda:\Theta\rightarrow \lambda(\Theta)$ exists and the Lipschitz property of such inverse mapping is assumed. It has been shown that such assumption holds for direct parameterization. However, such assumptions for continuous state-action space or other types of parameterization may not be valid. 


{\bf Multi-Objective Reinforcement Learning: } Similar to our setting, multi-objective reinforcement learning also considers the problem including several different objective functions. \cite{Multiobjective_overview1,Multiobjective_overview2} give a comprehensive overview of the research in multi-objective reinforcement learning. Two lines of methods have been studied, single-policy approach and multi-policy approach. In the multi-policy method, the goal is to achieve the Pareto Optimal solution, where the vector-valued utilities are used. Single policy method proposes some scalarization function to transform to problem back into single-objective MDP, which is similar to our setting. Several scalarization function such as weighted sum \cite{multiobjective_linear1,multiobjective_linear2}, W-learning \cite{multiobjective_W}, AHP \cite{multiobjective_AHP}, ranking \cite{multiobjective_ranking1,multiobjective_ranking2} have been proposed. However, none of these approaches  work for a combination of multiple objectives through a general non-linear concave function. Recently, such a problem has been investigated in \cite{mridulmb}, where regret guarantees for a model-based algorithm have been derived. In our work, we aim to provide guarantees for a model-free policy-gradient based algorithm. 

\section{Proof for the bias of Estimator in Eq. \eqref{eq:est_origin}}\label{sec_app_bias}
\begin{lemma}\label{lem_biased_est}
	The proposed estimator, Eq. \eqref{eq:est_origin}, is biased w.r.t $\nabla_\theta f(\bm{J}^{\pi_\theta})$
\end{lemma}
\begin{proof}
	By the law of total expectation
	\begin{equation}
		\begin{aligned}
			\mathbf{E}_{\tau_i,\tau_{j=1:N_2}}[g(\tau_i\vert \theta)]	&=\mathbf{E}_{\tau_i,\tau_{j=1:N_2}}\bigg[\sum_{t=0}^\infty \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_{m}^{\pi}}\big)\big(\sum_{h=t}^\infty \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\bigg]\\
			&=\mathbf{E}_{\tau_i}\bigg\{\mathbf{E}_{\tau_{j=1:N_2}}\bigg[\sum_{t=0}^\infty \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_{m}^{\pi}}\big)\big(\sum_{h=t}^\infty \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\bigg]\bigg|\tau_i\bigg\}\\
			&=\mathbf{E}_{\tau_i}\bigg\{\sum_{t=0}^\infty \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\mathbf{E}_{\tau_{j=1:N_2}}\bigg[\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_{m}^{\pi}}\bigg]\big(\sum_{h=t}^\infty \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\bigg|\tau_i\bigg\}\\
			&\overset{(*)}\neq\mathbf{E}_{\tau_i}\bigg\{\sum_{t=0}^\infty \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^\infty \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\bigg\}\\
			&=\nabla_\theta f(J_1^\pi(s),J_2^\pi(s),\cdots,J_M^\pi(s))
		\end{aligned}
	\end{equation}
	Notice that the key step (*) holds because
	\begin{equation}
		\begin{aligned}
			\mathbf{E}_{\tau_{j=1:N_2}}\bigg[\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_{m}^{\pi}}\bigg]&=\mathbf{E}_{\tau_{j=1:N_2}}\bigg[\frac{\partial f}{\partial J_m^\pi}(\frac{1}{N_2}\sum_{j=1}^{N_2}\sum_{t=0}^{\infty}\gamma^tr_1(s_t^j,a_t^j),\cdots, \frac{1}{N_2}\sum_{j=1}^{N_2}\sum_{t=0}^{\infty}\gamma^tr_M(s_t^j,a_t^j))\bigg]\\
			&\leq \frac{\partial f}{\partial J_m^\pi}(\mathbf{E}_{\tau_{j=1:N_2}}\bigg[\frac{1}{N_2}\sum_{j=1}^{N_2}\sum_{t=0}^{\infty}\gamma^tr_1(s_t^j,a_t^j)\bigg],\cdots,\mathbf{E}_{\tau_{j=1:N_2}}\bigg[\frac{1}{N_2}\sum_{j=1}^{N_2}\sum_{t=0}^{\infty}\gamma^tr_M(s_t^j,a_t^j)\bigg])\\
			&=\frac{\partial f}{\partial J_m^\pi}(J_1^\pi,\cdots,J_M^\pi)
		\end{aligned}
	\end{equation}
	where the inequality holds by the Assumption \ref{ass_concave_obj}
\end{proof}

\section{Bound the Bias for the Proposed Estimator}\label{sec_app_bound_bias}
\subsection{Proof for Lemma \ref{lem_bound_bias1}}
\begin{proof}
	By the triangle inequality, Assumptions \ref{ass_bound_reward} and \ref{ass_score}, we have
	\begin{equation}\label{eq:bound_bias1}
	\begin{split}
	\Vert g(\tau_i^H,\tau_j^H\vert \theta)-\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)\Vert&=\bigg\Vert\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_{m,H}^\pi}-\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=J_{m,H}^\pi}\big)\big(\sum_{h=t}^{H-1} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\bigg\Vert\\
	&\leq\frac{G}{1-\gamma}\bigg|\sum_{t=0}^{H-1}(\gamma^t-\gamma^H) \bigg(\sum_{m=1}^{M}\big(\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=\hat{J}_{m,H}^\pi}-\frac{\partial f}{\partial J_m^\pi}\bigg|_{J_m^\pi=J_{m,H}^\pi}\big)\bigg)\bigg|\\
	&\leq G\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^2} \sum_{m=1}^{M}\bigg|\frac{\partial f(\hat{J}_{m,H}^\pi)}{\partial J_m^\pi}-\frac{\partial f(J_{m,H}^\pi)}{\partial J_m^\pi}\bigg|\\
	&\leq GML_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^2} \Vert \hat{\bm{J}}_H^\pi-\bm{J}_H^\pi\Vert\\
	\end{split}
	\end{equation}
	where the last step follows from Assumption \ref{ass_partial_grad}. Moreover, an entry  in the difference $\hat{\bm{J}}_H^\pi-\bm{J}_H^\pi$ can be bounded as 
	\begin{equation}\label{eq:bound_JmH}
	\begin{split}
	\vert \hat{J}_{m,H}^\pi-J_{m,H}^\pi\vert&=\bigg|\frac{1}{N_2}\sum_{j=1}^{N_2}\sum_{t=0}^{H-1}\gamma^tr_m(s_t,a_t)-\mathbf{E}\big[\sum_{t=0}^{H-1}\gamma^tr_m(s_t,a_t)\big]\bigg|\\
	&\leq \sum_{t=0}^{H-1}\gamma^t\bigg|\frac{1}{N_2}\sum_{j=1}^{N_2}r_m(s_t,a_t)-\mathbf{E}[r_m(s_t,a_t)]\bigg|
	\end{split}
	\end{equation}
	By Hoeffding Lemma, if we have $N_2\geq \frac{M(1-\gamma^H)^2}{2(1-\gamma)^2\epsilon'^2}\log(\frac{2MH}{p})$, then
	\begin{equation}
	\begin{split}
	P\bigg(\bigg|\frac{1}{N_2}\sum_{j=1}^{N_2}r_m(s_t,a_t)-\mathbf{E}[r_m(s_t,a_t)]\bigg|\geq \frac{(1-\gamma)\epsilon'}{(1-\gamma^H)\sqrt{M}}\bigg)&\leq 2\exp(-\frac{2N_2^2\frac{(1-\gamma)^2\epsilon'^2}{(1-\gamma^H)^2 M}}{\sum_{j=1}^{N_2}(1-0)^2})\leq\frac{p}{MH}
	\end{split}
	\end{equation}
	Finally, by using an union bound, with probability at least $1-p$, we have
	\begin{equation}\label{eq:hoeffding_rm}
	\bigg|\frac{1}{N_2}\sum_{j=1}^{N_2}r_m(s_t,a_t)-\mathbf{E}[r_m(s_t,a_t)]\bigg|\leq \frac{(1-\gamma)\epsilon'}{(1-\gamma^H)\sqrt{M}}\quad \forall m\in[M],\forall t\in[0,H-1]
	\end{equation}
	Substituting Eq. \eqref{eq:hoeffding_rm} back into \eqref{eq:bound_JmH}, we have $\vert \hat{J}_{m,H}^\pi-J_{m,H}^\pi\vert\leq\frac{\epsilon'}{\sqrt{M}}$ and thus $\Vert \bm{J}_H^\pi-\hat{\bm{J}}_H^\pi\Vert_2\leq \epsilon'$, which gives the result in the statement of the Lemma.	
\end{proof}

\subsection{Proof for Lemma \ref{lem_bound_bias2}}
\begin{proof}
	Similar to Eq. \eqref{eq:bound_bias1}, we have
	\begin{equation}\label{eq:bound_bias2}
	\Vert \tilde{g}(\tau_i^H,\tau_j^H\vert \theta)-\tilde{g}(\tau_i^H,\tau_j\vert \theta)\Vert\leq GML_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^2} \Vert \bm{J}_H^\pi-\bm{J}^\pi\Vert\\
	\end{equation}
	By triangle inequality, the element of $\bm{J}_H^\pi-\bm{J}^\pi$ can be bounded by
	\begin{equation}\label{eq:bound_Jm}
	\begin{split}
	\vert J_{m,H}^\pi-J_m^\pi\vert&\leq \bigg|\mathbf{E}\big[\sum_{t=0}^{\infty}\gamma^tr_m(s_t,a_t)\big]-\mathbf{E}\big[\sum_{t=0}^{H-1}\gamma^tr_m(s_t,a_t)\big]\bigg|\\
	&\leq\sum_{t=H}^{\infty}\gamma^t\bigg|\mathbf{E}[r_m(s_t,a_t)]\bigg|\leq \frac{\gamma^H}{1-\gamma}
	\end{split}
	\end{equation}
	where the last step holds by Assumption \ref{ass_bound_reward}. Substituting Eq \eqref{eq:bound_Jm} back into \eqref{eq:bound_bias2} gives the result in the statement of the Lemma.
\end{proof}

\subsection{Proof for Lemma \ref{lem_bound_bias3}}
\begin{proof}
	By the triangle inequality,
	\begin{equation}
	\begin{aligned}
	&\Vert \tilde{g}(\tau_i^H,\tau_j\vert \theta)-g(\tau_i,\tau_j\vert \theta)\Vert\\
	&=\Vert\sum_{t=0}^\infty \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^{\infty} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)-\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^{\infty} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\\
	&+\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^{\infty} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)-\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^{H-1} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\Vert\\
	&\leq \Vert\sum_{t=H}^\infty \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^{\infty} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\Vert\\
	&+\Vert\sum_{t=0}^{H-1} \nabla_\theta \log\pi_\theta(a_t^i\vert s_t^i)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=H}^{\infty} \gamma^hr_m(s_h^i,a_h^i)\big)\bigg)\Vert\\
	&\leq \frac{MGC\gamma^H}{(1-\gamma)^2}+\frac{MGCH\gamma^H}{(1-\gamma)}=MGC\frac{\gamma^H(1+H(1-\gamma))}{(1-\gamma)^2}
	\end{aligned}
	\end{equation}
	where the last inequality holds by  Lemma \ref{lem_bound_parital} and Assumption \ref{ass_score}.
\end{proof}
\section{First Order Stationary Result for Policy Gradient}\label{sec_app_first_order}
\begin{lemma}\label{lem_fisrt_order}
	The policy gradient algorithm can achieve first-order stationary. More formally, if we choose the step size $\eta=\frac{1}{4L_J}$ and
	\begin{equation}
		N_1=\mathcal{O}(\frac{\sigma^2}{\epsilon})\quad
		N_2=\mathcal{O}(\frac{M^3}{(1-\gamma)^6\epsilon})\quad
		K=\mathcal{O}(\frac{M}{(1-\gamma)^2\epsilon})
		%
	\end{equation}
	then,
	\begin{equation}
		\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]\leq\epsilon
	\end{equation}
\end{lemma}
\begin{proof}
	Recall the definition of $\omega^k$ and $\tilde{\omega}^k$ in Eq. \eqref{eq:def_omega^k} and \eqref{eq:def_tilde_omega^k}, respectively. By Lemma \ref{lem_smooth}, we have
	\begin{equation}
		\begin{split}
		f(\bm{J}_H(\theta^{k+1}))&\geq f(\bm{J}_H(\theta^{k}))+\left<\nabla_\theta f(\bm{J}_H(\theta^k)),\theta^{k+1}-\theta^k\right>-\frac{L_J}{2}\Vert \theta^{k+1}-\theta^k\Vert_2^2\\
		&=f(\bm{J}_H(\theta^{k}))+\eta\left<\nabla_\theta f(\bm{J}_H(\theta^k)),\omega^k\right>-\frac{L_J\eta^2}{2}\Vert \omega^k\Vert_2^2\\
		&\overset{(a)}=f(\bm{J}_H(\theta^{k}))+\eta\left<\nabla_\theta f(\bm{J}_H(\theta^k)),\omega^k-\nabla_\theta f(\bm{J}_H(\theta^k))+\nabla_\theta f(\bm{J}_H(\theta^k))\right>\\
		&\quad-\frac{L_J\eta^2}{2}\Vert \omega^k-\nabla_\theta f(\bm{J}_H(\theta^k))+\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2\\
		&\overset{(b)}\geq f(\bm{J}_H(\theta^{k}))+\eta\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2-\eta\vert\left<\nabla_\theta f(\bm{J}_H(\theta^k)),\omega^k-\nabla_\theta f(\bm{J}_H(\theta^k))\right>\vert\\
		&\quad-L_J\eta^2\bigg(\Vert \omega^k-\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2+\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2\bigg)\\
		&\geq f(\bm{J}_H(\theta^{k}))+\eta\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2-\frac{\eta}{2}\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2-\frac{\eta}{2}\Vert \omega^k-\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2\\
		&\quad-L_J\eta^2\bigg(\Vert \omega^k-\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2+\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2\bigg)\\
		&= f(\bm{J}_H(\theta^{k}))+(\frac{\eta}{2}-L_J\eta^2)\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2-(\frac{\eta}{2}+L_J\eta^2)\Vert \omega^k-\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2\\
		&\overset{(c)}\geq f(\bm{J}_H(\theta^{k}))+(\frac{\eta}{2}-L_J\eta^2)\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2-(\eta+2L_J\eta^2)\Vert \omega^k-\tilde{\omega}^k\Vert_2^2\\
		&\quad -(\eta+2L_J\eta^2)\Vert \tilde{\omega}^k-\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2\\
		&\overset{(d)}\geq f(\bm{J}_H(\theta^{k}))+(\frac{\eta}{2}-L_J\eta^2)\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2-(\eta+2L_J\eta^2)\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2\\
		&-(\eta+2L_J\eta^2)\Vert \tilde{\omega}^k-\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2
		\end{split}
	\end{equation}
	where the step (a) holds by $\theta^{k+1}=\theta^k+\eta \omega^k$. Step (b) and (c) holds by Cauchy-Schwarz Inequality. Step (d) holds by Lemma \ref{lem_bound_bias1}. Then, take expectation with respect to the trajectories $\tau_i,\tau_j$ (Recall that $\theta^k,\theta^{k+1}$ is a function of $\tau_i,\tau_j$), we have
	\begin{equation}\label{eq:first_order_telescope}
		\begin{split}
		\mathbf{E}[f(\bm{J}_H(\theta^{k+1}))]&\geq \mathbf{E}[f(\bm{J}_H(\theta^{k}))]+(\frac{\eta}{2}-L_J\eta^2)\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]-(\eta+2L_J\eta^2)\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2\\
		&\quad-(\eta+2L_J\eta^2)\mathbf{E}[\Vert \tilde{g}^k-\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]\\
		&\geq \mathbf{E}[f(\bm{J}_H(\theta^{k}))]+(\frac{\eta}{2}-L_J\eta^2)\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]-(\eta+2L_J\eta^2)\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2\\
		&\quad-(\eta+2L_J\eta^2)\frac{\sigma^2}{N_1}
		\end{split}
	\end{equation}
	where the last step holds by Assumption \ref{ass_bounded_var}. Notice that in Eq. \eqref{eq:first_order_telescope}, $\mathbf{E}[f(\bm{J}_H(\theta^{k+1}))]$ and $\mathbf{E}[f(\bm{J}_H(\theta^{k}))]$ give a recursive form. Thus, telescoping from $k=0$ to $k=K-1$, we have
	\begin{equation}\label{eq:first_order}
		\frac{\mathbf{E}[f(\bm{J}_H(\theta^{K}))-f(\bm{J}_H(\theta^0))]}{K}\geq (\frac{\eta}{2}-L_J\eta^2)\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]-(\eta+2L_J\eta^2)[\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+\frac{\sigma^2}{N_1}]
	\end{equation}
	and thus
	\begin{equation}
		\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]\leq \frac{\frac{\mathbf{E}[f(\bm{J}_H(\theta^K))-f(\bm{J}_H(\theta^0))]}{K}+(\eta+2L_J\eta^2)[\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+\frac{\sigma^2}{N_1}]}{\frac{\eta}{2}-L_J\eta^2}
	\end{equation}
	Taking $\eta=\frac{1}{4L_J}$ and letting $N_1=\frac{18\sigma^2}{\epsilon}$, $K=\frac{48L_J\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^K))-\nabla_\theta f(\bm{J}_H(\theta^0))\Vert_2^2]}{\epsilon}$ and $\epsilon'=\frac{(1-\gamma)^2}{MGL_f}\sqrt{\frac{\epsilon}{6}}$, we have
	\begin{equation}
		\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]\leq\epsilon
	\end{equation}
	Recalling the definition of $N_2$ in the statement of Lemma \ref{lem_bound_bias1}, we have
	\begin{equation}
		N_2=\frac{6M^3G^2L_f^2(1-\gamma^H)^2}{(1-\gamma)^6\epsilon}\log(\frac{2MH}{p})
	\end{equation}
	\if 0
	\begin{equation}
		H=-\log_\gamma(\frac{2MGL_f}{(1-\gamma)^3}\sqrt{\frac{3M}{\epsilon}})
	\end{equation}
	\fi 
	Also, by the definition of $L_J$ in the lemma \ref{lem_smooth}
	\begin{equation}
		K=\frac{48MCB}{(1-\gamma)^2\epsilon}\mathbf{E}[\Vert\nabla_\theta f(\bm{J}_H(\theta^k))\Vert_2^2]
	\end{equation}
\end{proof}
\section{Proof of Lemma \ref{lem_framework}}\label{sec_app_framwork}
\begin{proof}
	Starting with the definition of KL divergence,
	\begin{equation}
	\begin{aligned}
	&\mathbf{E}_{s\sim d_\rho^{\pi^*}}[KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^k}(\cdot\vert s))-KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^{k+1}}(\cdot\vert s))]\\
	=&\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}\bigg[\log\frac{\pi_{\theta^{k+1}(a\vert s)}}{\pi_{\theta^k}(a\vert s)}\bigg]\\
	\overset{(a)}\geq&\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot(\theta^{k+1}-\theta^k)]-\frac{M}{2}\Vert\theta^{k+1}-\theta^k\Vert^2\\
	=&\eta\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot\omega^k]-\frac{M\eta^2}{2}\Vert\omega^k\Vert^2\\
	=&\eta\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot\omega_*^k]+\eta\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot(\omega^k-\omega_*^k)]-\frac{M\eta^2}{2}\Vert\omega^k\Vert^2\\
	=&\eta[f(\bm{J}^{\pi^*})-f(\bm{J}^{\pi_{\theta^k}})]+\eta\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot\omega_*^k]-\eta[f(\bm{J}^{\pi^*})-f(\bm{J}^{\pi_{\theta^k}})]\\
	&+\eta\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot(\omega^k-\omega_*^k)]-\frac{M\eta^2}{2}\Vert\omega^k\Vert^2\\		\overset{(b)}=&\eta[f(\bm{J}^{\pi^*})-f(\bm{J}^{\pi_{\theta^k}})]+\frac{\eta}{1-\gamma}\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}\bigg[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot(1-\gamma)\omega_*^k-\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_{\theta^k}})}{\partial J_m^{\pi_{\theta^k}}}A_m^{\pi_{\theta^k}}(s,a)\bigg]\\
	&+\eta\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}[\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot(\omega^k-\omega_*^k)]-\frac{M\eta^2}{2}\Vert\omega^k\Vert^2\\
	\overset{(c)}\geq&\eta[f(\bm{J}^{\pi^*})-f(\bm{J}^{\pi_{\theta^k}})]\\
	&-\frac{\eta}{1-\gamma}\sqrt{\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}\bigg[\bigg(\nabla_\theta\log\pi_{\theta^k}(a\vert s)\cdot(1-\gamma)\omega_*^k-\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_{\theta^k}})}{\partial J_m^{\pi_{\theta^k}}}A_m^{\pi_{\theta^k}}(s,a)\bigg)^2\bigg]}\\
	&-\eta\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}\Vert\nabla_\theta\log\pi_{\theta^k}(a\vert s)\Vert_2\Vert(\omega^k-\omega_*^k)\Vert_2-\frac{M\eta^2}{2}\Vert\omega^k\Vert^2\\
	\overset{(d)}\geq&\eta[f(\bm{J}^{\pi^*})-f(\bm{J}^{\pi_{\theta^k}})]-\frac{\eta\sqrt{\epsilon_{bias}}}{1-\gamma}-\eta G\Vert(\omega^k-\omega_*^k)\Vert_2-\frac{M\eta^2}{2}\Vert\omega^k\Vert^2\\
	\end{aligned}	
	\end{equation}
	where the step (a) holds by Assumption \ref{ass_score} and step (b) holds by Lemma \ref{lem_performance_diff}. Step (c) uses the convexity of the function $f(x)=x^2$. Finally, step (d) comes from the Assumption \ref{ass_transfer_error}. Rearranging items, we have
	\begin{equation}
	\begin{split}
	f(\bm{J}^{\pi^*})-f(\bm{J}^{\pi_{\theta^k}})&\leq \frac{\sqrt{\epsilon_{bias}}}{1-\gamma}+ G\Vert(\omega^k-\omega_*^k)\Vert_2+\frac{M\eta}{2}\Vert\omega^k\Vert^2\\
	&+\frac{1}{\eta}\mathbf{E}_{s\sim d_\rho^{\pi^*}}[KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^k}(\cdot\vert s))-KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^{k+1}}(\cdot\vert s))]
	\end{split}
	\end{equation}
	Summing from $k=0$ to $K-1$ and dividing by $K$, we get the desired result.
\end{proof}
\section{Proof for Properties of the Objective Function}\label{sec_app_property}
\subsection{Proof for Lemma \ref{lem_smooth}}
\begin{proof}
	In order to show the smoothness, it is sufficient to bound  $\Vert\nabla^2_\theta f(\bm{J}^{\pi_\theta})\Vert$ and $\Vert\nabla^2_\theta f(\bm{J}_H^{\pi_\theta})\Vert$. By  Eq. \eqref{eq:gradient_simplified}, we have
	\begin{equation}
	\begin{split}
	\Vert\nabla_\theta^2 f(\bm{J}^{\pi_\theta})\Vert&=\Vert\mathbf{E}_{\tau\sim p(\tau\vert \theta)}\bigg[\sum_{t=0}^{\infty}\nabla_\theta^2\log\pi_\theta(a_t\vert s_t)\bigg(\sum_{m=1}^{M}\frac{\partial f}{\partial J_m^\pi}\big(\sum_{h=t}^{\infty}\gamma^hr_m(s_h,a_h)\big)\bigg)\bigg]\Vert\\
	&\leq \frac{MC}{(1-\gamma)}\sum_{t=0}^{\infty}\gamma^t\Vert\nabla_\theta^2\log\pi_\theta(a_t\vert s_t)\Vert\leq \frac{MCB}{(1-\gamma)^2}
	\end{split}
	\end{equation}
	where the last inequality holds by the Assumption \ref{ass_score}. The smoothness property for the truncated version $f(\bm{J}_H^{\pi_\theta})$ can be proved similarly.
\end{proof}
\subsection{Proof for Lemma \ref{lem_bound_truncated}}
\begin{proof}
	Notice that $\tilde{g}(\tau_i,\tau_j\vert \theta)$ is an unbiased estimator for $\nabla_\theta f(\bm{J}^{\pi_\theta})$. Moreover, $\tilde{g}(\tau_i^H,\tau_j^H)$ is an unbiased estimator for $\nabla_\theta f(\bm{J}_H^{\pi_\theta})$. Thus,
	\begin{equation}
		\begin{split}
		\Vert\nabla_\theta f(\bm{J}^{\pi_\theta})-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert&\overset{(a)}=\Vert\mathbf{E}[\tilde{g}(\tau_i,\tau_j\vert \theta)-\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)]\Vert\leq \mathbf{E}\Vert\tilde{g}(\tau_i,\tau_j\vert \theta)-\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)\Vert\\
		&\overset{(b)}\leq \mathbf{E}\Vert\tilde{g}(\tau_i,\tau_j\vert \theta)-\tilde{g}(\tau_i^H,\tau_j\vert \theta)\Vert+\mathbf{E}\Vert\tilde{g}(\tau_i^H,\tau_j\vert \theta)-\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)\Vert\\
		&\overset{(c)}\leq M^{3/2}GL_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^3}\gamma^H+ MGC\frac{\gamma^H[1+H(1-\gamma)]}{(1-\gamma)^2}
		\end{split}
	\end{equation}
	where the step (a) and (b) hold by the triangle inequality. Step (c) holds by the Lemma \ref{lem_bound_bias2} and \ref{lem_bound_bias3}
\end{proof}	

\subsection{Proof for Lemma \ref{lem_performance_diff}}
\begin{proof}
	By the concavity of the function $f$, we have
	\begin{equation}
	\begin{aligned}
	f(\bm{J}^{\pi_\theta})&\leq  f(\bm{J}^{\pi_{\theta'}})+\nabla_{\bm{J}^{\pi_{\theta'}}} f(\bm{J}^{\pi_{\theta'}})^T(\bm{J}^{\pi_\theta}-\bm{J}^{\pi_{\theta'}})\\
	&=f(\bm{J}^{\pi_{\theta'}})+\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_{\theta'}})}{\partial J_m^{\pi_{\theta'}}}(J_m^{\pi_\theta}-J_m^{\pi_{\theta'}})\\
	&=f(\bm{J}^{\pi_{\theta'}})+\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_{\theta'}})}{\partial J_m^{\pi_{\theta'}}}\frac{1}{1-\gamma}\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta(\cdot\vert s)}\big[A_m^{\pi_{\theta'}}(s,a)\big]
	\end{aligned}
	\end{equation}
	where the last step comes from the policy gradient theorem \cite{sutton2000} for the standard reinforcement learning. Finally,  we get the desired result by rearranging terms.
\end{proof}
\section{Proof for Theorem \ref{thm1}}\label{sec_app_thm}
	In this part, we prove the Theorem \ref{thm1} by bounding the three terms on the right hand side of Eq. \eqref{eq:general_bound}. These terms are:  the difference between the update direction $\frac{G}{K}\sum_{k=0}^{K-1}\Vert(\omega^k-\omega_*^k)\Vert_2$, norm of estimated gradient $\frac{M\eta}{2K}\sum_{k=0}^{K-1}\Vert\omega^k\Vert^2$, and the term about KL divergence $\frac{1}{\eta K}\mathbf{E}_{s\sim d_\rho^{\pi^*}}[KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^0}(\cdot\vert s))]$
	
	\subsection{Bounding the difference between the update direction}
	Recall the estimated policy gradient update direction is
	\begin{equation}\label{eq:def_omega^k}
	\omega^k=\frac{1}{N_1}\sum_{i=1}^{N_1}g(\tau_i^H,\tau_j^H\vert \theta)
	\end{equation}
	and the true natural policy gradient update direction is
	\begin{equation}
	\omega_*^k=F_\rho(\theta_k)^{\dagger}\nabla_\theta f(\bm{J}^{\pi_\theta})
	\end{equation}
	We define an auxiliary update direction as
	\begin{equation}\label{eq:def_tilde_omega^k}
		\tilde{\omega}^k=\frac{1}{N_1}\sum_{i=1}^{N_1}\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)
	\end{equation}
	Thus, we can decompose the difference as 
	\begin{equation}
	\begin{split}
	&\bigg(\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}\Vert\omega^k-\omega_*^k\Vert_2\bigg)^2\leq\frac{1}{K}\sum_{k=0}^{K-1}\bigg(\mathbf{E}\Vert\omega^k-\omega_*^k\Vert_2\bigg)^2\leq \frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}\bigg[\Vert\omega_k-\omega_*^k\Vert_2^2\bigg]\\
	&=\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}\bigg[\Vert(\omega^k-\tilde{\omega}^k)+(\tilde{\omega}^k-\nabla_\theta f(\bm{J}_H^{\pi_\theta}))+(\nabla_\theta f(\bm{J}_H^{\pi_\theta})-\nabla_\theta f(\bm{J}^{\pi_\theta}))+(\nabla_\theta f(\bm{J}^{\pi_\theta})-F_\rho(\theta^k)^\dagger\nabla_\theta f(\bm{J}^{\pi_\theta}))\Vert_2^2\bigg]\\
	&\leq \frac{4}{K}\sum_{k=0}^{K-1}\mathbf{E}\bigg[\Vert\omega^k-\tilde{\omega}^k\Vert_2^2\bigg]+\frac{4}{K}\sum_{k=0}^{K-1}\mathbf{E}\bigg[\Vert\tilde{\omega}^k-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert_2^2\bigg]+\frac{4}{K}\sum_{k=0}^{K-1}\mathbf{E}\bigg[\Vert\nabla_\theta f(\bm{J}^{\pi_\theta})-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert_2^2\bigg]\\
	&+\frac{4}{K}\sum_{k=0}^{K-1}\mathbf{E}\bigg[\Vert\nabla_\theta f(\bm{J}^{\pi_\theta})-F_\rho(\theta^k)^\dagger\nabla_\theta f(\bm{J}^{\pi_\theta})\Vert_2^2\bigg]
	\end{split}
	\end{equation}
The different terms in the above are bounded as follows:

	\begin{itemize}
		\item Bounding $\mathbf{E}\bigg[\Vert\omega^k-\tilde{\omega}^k\Vert_2^2\bigg]$:  By  Lemma \ref{lem_bound_bias1}, with $N_2$ large enough, for any $\tau_i$ and $\theta$, we have
		\begin{equation}
			\Vert 		g(\tau_i^H\vert\theta)-\tilde{g}(\tau_i^H\vert\theta)\Vert_2\leq MGL_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^2}\epsilon'
		\end{equation}
		Thus,
		\begin{equation}
		\begin{split}
		\Vert\omega^k-\tilde{\omega}^k\Vert_2&=\Vert\frac{1}{N_1}\sum_{i=1}^{N_1}(g(\tau_i^H\vert \theta^k)-\tilde{g}(\tau_i^H\vert \theta^k))\Vert_2\leq \frac{1}{N_1}\sum_{i=1}^{N_1}\Vert(g(\tau_i^H\vert \theta^k)-\tilde{g}(\tau_i^H\vert \theta^k))\Vert_2\\
		&\leq MGL_f\frac{1-\gamma^H-H\gamma^H(1-\gamma)}{(1-\gamma)^2}\epsilon'\leq \frac{MGL_f}{(1-\gamma)^2}\epsilon'
		\end{split}
		\end{equation}
		Thus,
		\begin{equation}
			\mathbf{E}\bigg[\Vert\omega^k-\tilde{\omega}^k\Vert_2^2\bigg]\leq \frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2
		\end{equation}
		\item Bounding $\mathbf{E}\bigg[\Vert\tilde{\omega}^k-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert_2^2\bigg]$: 	Notice that $\tilde{g}(\tau^H\vert\theta)$ is an unbiased estimator for $\nabla_\theta f(\bm{J}_H^{\pi_\theta})$ and thus by  Assumption \ref{ass_bounded_var}, we have $
		\mathbf{E}\bigg[\Vert\omega^k-\tilde{\omega}^k\Vert_2^2\bigg]\leq \frac{\sigma^2}{N_1}$
		\item Bounding $\mathbf{E}\bigg[\Vert\nabla_\theta f(\bm{J}^{\pi_\theta})-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert_2^2\bigg]$: 
		By Lemma \ref{lem_bound_truncated}, we have
		\begin{equation}
			\mathbf{E}\bigg[\Vert\nabla_\theta f(\bm{J}^{\pi_\theta})-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert_2^2\bigg] \leq \frac{M^2G^2\gamma^{2H}}{(1-\gamma)^4}\bigg[ \sqrt{M}L_f+C[1+H(1-\gamma)]\bigg]^2
		\end{equation}
		\item Bounding $\mathbf{E}\bigg[\Vert\nabla_\theta f(\bm{J}^{\pi_\theta})-F_\rho(\theta^k)^\dagger\nabla_\theta f(\bm{J}^{\pi_\theta})\Vert_2^2\bigg]$:	By  Assumption \ref{ass_pd}, we have
		\begin{equation}
			\begin{split}
			&\mathbf{E}\bigg[\Vert\nabla_\theta f(\bm{J}^{\pi_\theta})-F_\rho(\theta^k)^\dagger\nabla_\theta f(\bm{J}^{\pi_\theta})\Vert_2^2\bigg]\leq (1+\frac{1}{\mu_F})^2\mathbf{E}[\Vert \nabla_\theta f(\bm{J}^{\pi_k})\Vert_2^2]\\
			&\leq (1+\frac{1}{\mu_F})^2\bigg(2\mathbf{E}[\Vert \nabla_\theta f(\bm{J}_H^{\pi_k})\Vert_2^2]+2\mathbf{E}[\Vert \nabla_\theta f(\bm{J}^{\pi_k})-\nabla_\theta f(\bm{J}_H^{\pi_k})\Vert_2^2]\bigg)\\
			&\leq (1+\frac{1}{\mu_F})^2\bigg(2\mathbf{E}[\Vert \nabla_\theta f(\bm{J}^{\pi_k})]+\frac{2M^2G^2\gamma^{2H}}{(1-\gamma)^4}\bigg[ \sqrt{M}L_f+C[1+H(1-\gamma)]\bigg]^2\bigg)
		\end{split}
		\end{equation}
	\end{itemize}
	Finally, we obtain the bound
	\begin{equation}
	\begin{split}
	&\bigg(\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}\Vert\omega^k-\omega_*^k\Vert_2\bigg)^2\leq 4\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+4\frac{\sigma^2}{N_1}+4\frac{M^2G^2\gamma^{2H}}{(1-\gamma)^4}\bigg[ \sqrt{M}L_f+C[1+H(1-\gamma)]\bigg]^2\\
	&+4(1+\frac{1}{\mu_F})^2\bigg(\frac{2}{K}\sum_{k=0}^{K-1}\mathbf{E}[\Vert \nabla_\theta f(\bm{J}^{\pi_k})]+\frac{2M^2G^2\gamma^{2H}}{(1-\gamma)^4}\bigg[ \sqrt{M}L_f+C[1+H(1-\gamma)]\bigg]^2\bigg)\\
	&\overset{(a)}=(1+2(1+\frac{1}{\mu_F})^2)4\frac{M^2G^2\gamma^{2H}}{(1-\gamma)^4}\bigg[ \sqrt{M}L_f+C[1+H(1-\gamma)]\bigg]^2+4\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+4\frac{\sigma^2}{N_1}\\
	&+8(1+\frac{1}{\mu_F})^2\frac{\frac{\mathbf{E}[f(\bm{J}_H(\theta^K))-f(\bm{J}_H(\theta^0))]}{K}+(\eta+2L_J\eta^2)[\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+\frac{\sigma^2}{N_1}]}{\frac{\eta}{2}-L_J\eta^2}\\
	&=(1+2(1+\frac{1}{\mu_F})^2)4\frac{M^2G^2\gamma^{2H}}{(1-\gamma)^4}\bigg[ \sqrt{M}L_f+C[1+H(1-\gamma)]\bigg]^2+(1+6(1+\frac{1}{\mu_F})^2)4\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2\\
	&+(1+6(1+\frac{1}{\mu_F})^2)4\frac{\sigma^2}{N_1}+128(1+\frac{1}{\mu_F})^2L_J\frac{\mathbf{E}[f(\bm{J}_H(\theta^K))-f(\bm{J}_H(\theta^0))]}{K}
	\end{split}
	\end{equation}
	where the step (a) requires the first-order stationary property Eq. \eqref{eq:first_order} and it is proved in the Lemma \ref{lem_fisrt_order} in the Appendix \ref{sec_app_first_order}. Given the fixed $\epsilon$, choose the value for $H,\epsilon',N_1,K$ as follows,
	\begin{equation}
	\frac{1}{4}\bigg(\frac{\epsilon^2}{3G^2}\bigg)\geq (1+2(1+\frac{1}{\mu_F})^2)\frac{4M^2G^2\gamma^{2H}}{(1-\gamma)^4}\bigg[ \sqrt{M}L_f+C[1+H(1-\gamma)]\bigg]^2
	\end{equation}
	\begin{equation}
		\epsilon'^2\leq\frac{1}{4(1+6(1+\frac{1}{\mu_F})^2)}\frac{(1-\gamma)^4}{M^2G^2L_f^2}\cdot\frac{1}{4}\bigg(\frac{\epsilon^2}{3G^2}\bigg)
	\end{equation}
	\begin{equation}
		N_1\geq\frac{(1+6(1+\frac{1}{\mu_F})^2)4\sigma^2}{\frac{1}{4}\bigg(\frac{\epsilon^2}{3G^2}\bigg)}
	\end{equation}
	\begin{equation}
		K\geq\frac{128(1+\frac{1}{\mu_F})^2L_J\mathbf{E}[f(\bm{J}_H(\theta^K))-f(\bm{J}_H(\theta^0))]}{\frac{1}{4}\bigg(\frac{\epsilon^2}{3G^2}\bigg)}
	\end{equation}
	then we have
	\begin{equation}
	\frac{G}{K}\sum_{k=0}^{K-1}\mathbf{E}[\Vert\omega^k-\omega_*^k]\Vert_2\leq \frac{\epsilon}{3}
	\end{equation}
	Given the choice of $H,N_1,\epsilon',K$, the dependence of $N_1,N_2,K$ and $H$ on $\sigma, \epsilon, 1-\gamma$ are as follows.
	\begin{equation}\label{eq_require1}
		N_1=\mathcal{O}(\frac{\sigma^2}{\epsilon^2})\quad N_2=\mathcal{O}(\frac{M^3}{(1-\gamma)^6\epsilon^2})\quad K=\mathcal{O}(\frac{M}{(1-\gamma)^2\epsilon^2})\quad
		H=\mathcal{O}(\log\frac{M}{(1-\gamma)\epsilon})
	\end{equation}
	\subsection{Bounding the norm of estimated gradient}
	\begin{equation}
	\begin{split}
	\frac{M\eta}{2K}\sum_{k=0}^{K-1}\Vert \omega^k\Vert_2^2&\leq\frac{M\eta}{2}\bigg[\frac{3}{K}\sum_{k=0}^{K-1}\Vert \omega^k-\tilde{\omega}^k\Vert_2^2+\frac{3}{K}\sum_{k=0}^{K-1}\Vert \tilde{\omega}^k-\nabla_\theta f(\bm{J}_H^{\pi_\theta})\Vert_2^2+\frac{3}{K}\sum_{k=0}^{K-1}\Vert \nabla_\theta f(\bm{J}_H^{\pi_\theta}))\Vert_2^2\bigg]\\
	&\leq \frac{M\eta}{2}\bigg[3\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+3\frac{\sigma^2}{N_1}+3\frac{\frac{\mathbf{E}[f(\bm{J}_H(\theta^K))-f(\bm{J}_H(\theta^0))]}{K}+(\eta+2L_J\eta^2)[\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+\frac{\sigma^2}{N_1}]}{\frac{\eta}{2}-L_J\eta^2}\bigg]\\
	&=M\eta\bigg[6\frac{M^2G^2L_f^2}{(1-\gamma)^4}\epsilon'^2+6\frac{\sigma^2}{N_1}+24L_J\frac{\mathbf{E}[f(\bm{J}_H(\theta^K))-f(\bm{J}_H(\theta^0))]}{K}\bigg]
	\end{split}
	\end{equation}
	Given the fixed $\epsilon$, choose the value for $\epsilon',N_1,K$ as follows,
	\begin{equation}
	\epsilon'^2\leq\frac{(1-\gamma)^4}{M^2G^2L_f^2}\cdot\frac{1}{6M\eta}\bigg(\frac{\epsilon}{9}\bigg)
	\end{equation}
	\begin{equation}
	N_1\geq\frac{54\sigma^2}{\epsilon}
	\end{equation}
	\begin{equation}
	K\geq\frac{216L_J\mathbf{E}[f(\bm{J}_H(\theta^K))-f(\bm{J}_H(\theta^0))]}{\epsilon}
	\end{equation}
	then we have
	\begin{equation}
	\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{E}[\Vert\omega^k]\Vert_2^2\leq \frac{\epsilon}{3}
	\end{equation}
	Given the choice of $\epsilon',N_1,K$, the dependence of $N_1,N_2,K$ and $H$ on $\sigma, \epsilon, 1-\gamma$ are as follows.
	\begin{equation}\label{eq_require2}
	N_1=\mathcal{O}(\frac{\sigma^2}{\epsilon})\quad N_2=\mathcal{O}(\frac{M^3}{(1-\gamma)^6\epsilon})\quad K=\mathcal{O}(\frac{M}{(1-\gamma)^2\epsilon})\quad
	H=\mathcal{O}(\log\frac{M}{(1-\gamma)\epsilon})
	\end{equation}
	\subsection{Bounding the KL divergence}
	It is obvious if we choose
	\begin{equation}
	K\geq\frac{3\mathbf{E}_{s\sim d_\rho^{\pi^*}}[KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^0})]}{\eta\epsilon(\cdot\vert s)}
	\end{equation}
	then
	\begin{equation}
	\frac{1}{\eta K}\mathbf{E}_{s\sim d_\rho^{\pi^*}}[KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^0})]\leq \frac{\epsilon}{3}
	\end{equation}
	In other word, the dependence of $K$ on $\epsilon$ is
	\begin{equation}\label{eq_require3}
		K=\mathcal{O}(\frac{M}{\epsilon})
	\end{equation}

\section{Evaluation Details}\label{app:simulations}
We now describe the details of the environment and the simulation setup used. We begin by describing the simulation setup. We use softmax parameterization for implementing our policies. Further, we use PyTorch version 1.0.1 to implement the policies and perform gradient ascent. The experiments are run on a machine with Intel i9 processor with $36$ logical cores running at $3.00$ GHz each. The machines are equipped with Nvidia GeForce RTX 2080 GPU. Each of the $10$ independent runs for both environment took about $500$ seconds to finish.

We now explain the details of the environments:

\textbf{1. Wireless Scheduler Environment} \\
In our wireless scheduler environment, we have $4$ users connected to a base station. At time step $t$, each of the user can be in two states, \textit{good} or \textit{bad}. The scheduler has access to the states of all the users. After observing the states of the users, the scheduler selects a user and the user can download some data with data rate which depends on the state of the user. If the user is not selected by the user, it observe $0$ data rate. The system as $2^4 = 16$ states and $4$ actions. The data rate observed by a user acts as reward for the system. The reward matrix for the users is given in Table \ref{tab:rate_table}. Further, the channel conditions of the users are time varying. After every time step the state of a user can toggle with probability $0.1$ and remain the same with probability $0.9$.

\begin{table}[htbp]
\centering
\begin{tabular}{|c |c |c| c| c|} 
 \hline
 Agent state & $r_{1,t}$ & $r_{2,t}$ & $r_{3,t}$ & $r_{4,t}$\\
 \hline
 $good$ &1.5 & 2.25& 1.25 & 1.5\\ 
 \hline
 $bad$ &0.768 & 1.0& 0.384 & 1.12\\
 \hline
\end{tabular}
\caption{Agent rate (in Mbps) based on agent state. Rate values are practically observable data rates over a wireless network such as 4G-LTE.}
\label{tab:rate_table}
\end{table}

For the wireless scheduler environment setup, we consider the the joint objective function $f$ defined as
\begin{align*}
f(J_1^\pi, \cdots, J_K^\pi) = \sum_{k=1}^K \frac{-1}{J_k^\pi}.
\end{align*}
For this joint objective function, the gradient becomes,
\begin{align*}
\nabla_\theta f(J_1^\pi, \cdots, J_K^\pi) = \sum_{k=1}^K \frac{1}{\left(J_k^\pi\right)^2}\nabla_\theta J_k^\pi
\end{align*}

For the gradient ascent of objective, we used PyTorch's Adam \cite{kingma2015adam} optimization with learning rate of $0.01$.


\textbf{2. Queuing Environment}\\
For our queuing environment, we consider a server serving customers coming from $4$ queues. Each queue follows Poisson arrivals with different arrival rates given in Table \ref{tab:arrival_rates}. The server has access to the length of the queues. On observing the length of the queue, the server selects a queue to process. If the a customer from a queue is served, the queue gets a reward of $1$ unit.

\begin{table}[ht]
    \centering
    \begin{tabular}{|c|c|c|c|}
    \hline
        $\lambda_1$ &$\lambda_2$ & $\lambda_3$& $\lambda_4$  \\
    \hline
        $0.08$ & $0.16$ & $0.24$ & $0.32$\\
    \hline
    \end{tabular}
    \caption{Arrival rates of the multiple queues for Queuing system environment}
    \label{tab:arrival_rates}
\end{table}

For the queuing environment setup, we consider the the joint objective function $f$ defined as
\begin{align*}
f(J_1^\pi, \cdots, J_K^\pi) = \sum_{k=1}^K \log{J_k^\pi}.
\end{align*}
For this joint objective function, the gradient becomes,
\begin{align*}
\nabla_\theta f(J_1^\pi, \cdots, J_K^\pi) = \sum_{k=1}^K \frac{1}{J_k^\pi}\nabla_\theta J_k^\pi
\end{align*}

For the gradient ascent of objective, we used PyTorch's Adam \cite{kingma2015adam} optimization with learning rate of $0.005$.
\section{Conclusion}\label{sec_conclusion}
In this paper, we formulate a problem which  optimizes a general concave function of multiple objectives. We propose a policy-gradient based approach for the problem, where an estimator for the gradient is used.  We analyze the bias of the policy gradient estimator and show the global convergence result with a vanilla policy gradient algorithm. Extension of the proposed approach to evaluate the convergence rate guarantees of the Natural Policy Gradient and the variance reduced algorithms is an important future direction. 


\section{Evaluations}\label{evaluations}
To validate the understanding of our analysis, we perform evaluations using two different environments with multiple objectives and concave utility functions. We study the impact of the number of trajectory used for gradient estimation. We keep the number of trajectories $N_1 = N_2 = N$ and vary $N$ from $1, 4, 16$, and $64$ and observe the convergence rates for a softmax policy parameterization. Implementation details are provided in Appendix \ref{app:simulations}.

The first environment is wireless scheduler to which $4$ users are connected. Each user can exist in two states, \textit{good} or \textit{bad}. The action is the user to which the scheduler allocates the resource. This system has $16$ states with $4$ actions. At time $t$, each user $k$ achieves different rates $r_{k,t}$ based on their states and resource allocation. We let the length of episode be $H=500$ steps. The joint objective function is $\alpha$-concave utility defined as:
\begin{align}
f(\sum_t r_{1,t}, \cdots, \sum_t r_{K,t}) = -\sum\nolimits_{k=1}^K H/ \left(\sum_t r_{k,t}\right)
\end{align}
The second environment is a server serving $4$ queues with Poisson arrivals with different arrival rates. The system state is $4$ dimensional vector of the length of the $4$ queues. The action at each time is the queue which the server serves. At time $t$, each queue $k$ achieves a reward of $1$ unit if a customer from this queue is served.  The joint objective function is sum-logarithmic utility defined as:
\begin{align}
f(\sum_t r_{1,t}, \cdots, \sum_t r_{K,t}) = \sum\nolimits_{k=1}^K \log \left(\sum_t r_{k,t}/{H}\right),
\end{align}
\begin{figure*}[th]
    \centering
    \subfigure[Wireless Scheduler Environment]{
        \includegraphics[trim=0.15in 0.05in 0.5in 0.45in, clip, width=0.47\textwidth]{wireless_experiment.png}
        \label{fig:wireless_convergence}
    }\hspace{-0.1in}
    \subfigure[Queuing Environment]{
        \includegraphics[trim=0in 0.05in 0.5in 0.45in, clip, width=0.49\textwidth]{queueing_experiment.png}
        \label{fig:queuing_convergence}
    }
    \caption{ Convergence plot for the joint objective policy gradient algorithms on (a) wireless scheduler environment and (b) queuing environment. As the number of trajectories used for sampling gradient of the function increase, the convergence becomes steeper.
    }
    \label{fig:convergence_plots}
    \vspace{-.1in}
\end{figure*}

We plot the behavior of the policy gradient for joint optimization for different values of $N$ in Figure \ref{fig:convergence_plots}. We  run $10$ independent iterations and plot mean in solid lines and the shaded region is $\pm$ standard deviation. For both the setups, we find that increasing $N$, the number of trajectories used for sampling gradient of the function, leads to faster convergence of the joint reward objective. We observe that with single trajectory the convergence towards a good policy is not good for both the enviroments. For the wireless environment with $16$ state, even with $4$ trajectories, we start observing convergence \ref{fig:wireless_convergence}. We infer that for joint optimization of multiple objectives it is absolutely necessary to increase the number trajectories as the environments become complex.
\section{Formulation}\label{sec_formulation}
	We consider an infinite horizon discounted Markov Decision Process (MDP) $\mathcal{M}$ defined by the tuple $(\mathcal{S},\mathcal{A},\mathbb{P},r_1,r_2,\cdots,r_M,\gamma,\rho)$, where $\mathcal{S}$ and $\mathcal{A}$ denote the  state and  action space, respectively. $\mathbb{P}: \mathcal{S}\times\mathcal{A}\rightarrow\Delta^{\mathcal{S}}$ (where $\Delta^{\mathcal{S}}$ is a probability simplex over $\mathcal{S}$) denotes the transition probability distribution from a state-action pair to another state. $M$ denotes the number of agents and $r_m: \mathcal{S}\times\mathcal{A}\rightarrow \mathbb{R}$ denotes the reward for the $m^{th}$ agents. $\gamma\in(0,1)$ is the discounted factor and $\rho: \mathcal{S}\rightarrow \Delta^{\mathcal{S}}$ is the distribution for initial state. In this paper, we make following assumption.
	\begin{assumption}\label{ass_bound_reward}
		The absolute value of the reward functions $r_m,m\in[M]$ is bounded by some constant. Without loss of generality, we assume $r_m\in[0,1],\forall m\in[M]$.
	\end{assumption} 
	Define a joint stationary policy $\pi:\mathcal{S}\rightarrow\Delta^{\mathcal{A}}$ that maps a state $s\in\mathcal{S}$ to a probability distribution of actions with a probability assigned to each action $a\in\mathcal{A}$. At the beginning of the MDP, an initial state $s_0\sim\rho$ is given and all agents together make a decision $a_0\sim\pi(\cdot\vert s_0)$. Each agent receives its reward $r_m(s_0,a_0)$ and together they transits to a new state $s_1\sim\mathbb{P}(\cdot\vert s_0,a_0)$. We define the value function $J_m^{\pi}$ for the $m^{th}$ agent following policy $\pi$ as a discounted sum of reward over infinite horizon. By  Assumption \ref{ass_bound_reward},
	\begin{equation}
		J_m^{\pi}=\mathbf{E}_{\rho,\pi,\mathbb{P}}\bigg[\sum_{t=0}^{\infty}\gamma^tr_m(s_t,a_t)\bigg]
	\end{equation} 
	where $s_0\sim\rho$, $a_t\sim\pi(\cdot\vert s_t)$ and $s_{t+1}\sim\mathbb{P}(\cdot\vert s_t,a_t)$.
	The agents collaboratively aim to maximize the joint objective function $f:\mathbb{R}^M\rightarrow\mathbf{R}$, which is a function of the long-term discounted reward of individual agent. Formally, the problem is written as
	\begin{equation}\label{eq:origin_problem}
		\max_{\pi}f(J_1^{\pi},J_2^{\pi},\cdots,J_M^{\pi})
	\end{equation}
	We consider a policy-gradient based algorithm on this problem and parameterize the policy $\pi$ as $\pi_\theta$ for some parameter $\theta\in\Theta$ such as softmax parameterization or a deep neural network. Commonly, the log-policy function $\log\pi_\theta(a\vert s)$ is called log-likelihood function and we make the following assumption.
	\begin{assumption}\label{ass_score}
		The log-likelihood function is $G$-Lipschitz and $B$-smooth. Formally,
		\begin{equation}
			\Vert \nabla_\theta\log\pi_\theta(a\vert s)\Vert\leq G\quad\forall \theta\in\Theta,\forall (s,a)\in\mathcal{S}\times\mathcal{A}
		\end{equation}
		\begin{equation}
		\Vert \nabla_\theta\log\pi_{\theta_1}(a\vert s)-\nabla_\theta\log\pi_{\theta_2}(a\vert s)\Vert\leq B\Vert \theta_1-\theta_2\Vert\quad\forall \theta_1,\theta_2 \in\Theta,\forall (s,a)\in\mathcal{S}\times\mathcal{A}
		\end{equation}
	\end{assumption}
	\begin{remark}
		The Lipschitz and smoothness properties for the log-likelihood are quite common in the field of policy gradient algorithm \cite{Alekh2020,Mengdi2021,Yanli2020}. Such properties can also be verified for simple parameterization such as Gaussian policy.   
	\end{remark}
	Define the value function vector $\bm{J}^{\pi_\theta}=(J_1^{\pi_\theta},J_2^{\pi_\theta},\cdots,J_M^{\pi_\theta})$. The original problem, Eq. \eqref{eq:origin_problem}, can be rewritten as
	\begin{equation}\label{eq:parameterized_problem}
		\max_{\theta\in\Theta}f(\bm{J}^{\pi_\theta})
	\end{equation}
	We make the following assumptions on the objective function $f$:
	\begin{assumption}\label{ass_concave_obj}
		The objective function $f$ is jointly concave. Hence for any arbitrary distribution $\mathcal{D}$, the following holds.
		\begin{equation}
			f(\mathbf{E}_{\bm{x} \sim\mathcal{D}}[\bm{x}])\geq \mathbf{E}_{\bm{x} \sim\mathcal{D}}[f(\bm{x})])\quad\forall \bm{x}\in\mathbb{R}^M
		\end{equation}
	\end{assumption}
	\begin{remark}(Non-Concave Optimization)
		It is worth noticing that the above problem is a non-concave optimization problem despite the above joint-concave assumption on the objective function. This  is because the parameterized value function $\bm{J}_m^{\pi_\theta}$ is non-concave with respect to $\theta$ (See Lemma 3.1 in \cite{Alekh2020}). Thus, the standard theory from convex optimization can't be directly applied to this problem. Moreover, such assumption is common in the literature, and has also been adopted in \cite{Mengdi2020,Mengdi2021}.
	\end{remark}
	\begin{assumption}\label{ass_partial_grad}
		All partial derivatives of function $f$ are assumed to be locally $L_f$-Lipschitz functions. Formally,
		\begin{equation}
		\vert \frac{\partial f}{\partial x_i}(\bm{y}_1)-\frac{\partial f}{\partial x_i}(\bm{y}_2)\vert\leq L_f\Vert \bm{y}_1-\bm{y}_2\Vert_2\quad \forall\bm{y}_1,\bm{y}_2\in[0,\frac{1}{1-\gamma}]^M,\forall i\in[M]
		\end{equation}
	\end{assumption}
	\begin{remark}
		By Assumption \ref{ass_bound_reward},  $J_m^{\pi_\theta}$ is bounded in $[0,\frac{1}{1-\gamma}]$. Thus, it is enough to assume the locally Lipschitiz property for the partial derivatives of the objective. Such an assumption has also been adopted widely for the general objective function \cite{Mengdi2020,Mengdi2021}.
	\end{remark}
	\noindent Finally, based on the Assumption. \ref{ass_partial_grad}, we derive the following result for the objective function.
	\begin{lemma}\label{lem_bound_parital}
		All partial derivative functions of $f$ are locally bounded by a constant. Formally,
		\begin{equation}
			\bigg| \frac{\partial f}{\partial x_i}(\bm{y})\bigg|\leq C\quad \forall \bm{y}\in[0,\frac{1}{1-\gamma}]^M,\forall i\in[M]
		\end{equation}
	\end{lemma}
	\begin{proof}
		By Assumption \ref{ass_partial_grad}, the partial derivative function is locally Lipschitz and thus is continuous on the set $[0,\frac{1}{1-\gamma}]^M$, which is compact. Since a continuous function with a compact set is bounded, the result follows.
	\end{proof}

\section{Main Result}\label{sec_result}
	Before stating the convergence result for the policy gradient algorithm, we describe the following assumptions which will be needed for the main result. 
	
	\begin{assumption}\label{ass_bounded_var}
		The auxiliary estimator $\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)$ defined in Eq. \eqref{eq:truncated_auxi} has bounded variance. Formally,
		\begin{equation}
			Var(\tilde{g}(\tau_i^H,\tau_j^H\vert \theta):=\mathbf{E}[\Vert \tilde{g}(\tau_i^H,\tau_j^H\vert \theta)-\mathbf{E}[\tilde{g}(\tau_i^H,\tau_j^H\vert \theta)]\Vert^2]\leq \sigma^2
		\end{equation}
		for any $\theta$ and $\tau_i^H,\tau_j^H\sim p^H(\cdot\vert\theta)$, where $p^H(\cdot\vert\theta)$ is a truncated version of $p(\cdot\vert\theta)$ defined in Eq. \eqref{eq:traj_dis}.
	\end{assumption}
	\begin{remark}
		In the standard reinforcement learning problem, it is common to assume that variance of the estimator is bounded \cite{Yanli2020}, \cite{Xu2019} and \cite{xu2020}. Such assumption has been verified for Gaussian policy \cite{tingting2011} and \cite{Pirotta2013}. By  Lemma \ref{lem_bound_parital}, it can be verified similarly in the  multi-objective setting.
	\end{remark}
	\begin{assumption}\label{ass_pd}
		For all $\theta\in\mathbb{R}^d$, the Fisher information matrix induced by policy $\pi_{\theta}$ and initial state distribution $\rho$ satisfies
		\begin{equation}
		F_\rho(\theta)=\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta}[\nabla_{\theta}\log\pi_\theta(a|s)\nabla_\theta\log\pi_\theta(a|s)^T]\succeq\mu_F\cdot \mathbf{I}_d
		\end{equation}
		for some constant $\mu_F>0$
	\end{assumption}
	\begin{remark}
		The positive definiteness assumption is standard in the field of policy gradient based algorithms \cite{Kakade2002,Peters2008,Yanli2020,Kaiqing2019}. A  common example which satisfies such assumption is Gaussian policy with mean parameterized linearly (See Appendix B.2 in \cite{Yanli2020}).
	\end{remark}
	\begin{assumption}\label{ass_transfer_error}
		Define the transferred function approximation error as below
		\begin{equation}\label{eq:transfer_error}
			L_{d_\rho^{\pi^*},\pi^*}(\omega_*^\theta,\theta)=\mathbf{E}_{s\sim d_\rho^{\pi^*}}\mathbf{E}_{a\sim\pi^*(\cdot\vert s)}\bigg[\bigg(\nabla_\theta\log\pi_{\theta}(a\vert s)\cdot(1-\gamma)\omega_*^\theta-\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_\theta})}{\partial J_m^{\pi_\theta}}A_m^{\pi_\theta}(s,a)\bigg)^2\bigg]
		\end{equation}
		We assume that this error satisfies $L_{d_\rho^{\pi^*},\pi^*}(\omega_*^\theta,\theta)\leq \epsilon_{bias}$ for any $\theta\in\Theta$, where $\omega_*^{\theta}$ is given as
		\begin{equation}\label{eq:NPG_direction}
			\omega_*^{\theta}=\arg\min_{\omega}\mathbf{E}_{s\sim d_\rho^{\pi_{\theta}}}\mathbf{E}_{a\sim\pi_{\theta}(\cdot\vert s)}\bigg[[\nabla_\theta\log\pi_{\theta}(a\vert s)\cdot (1-\gamma)\omega-\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_{\theta}})}{\partial J_m^{\pi_{\theta}}}A_m^{\pi_{\theta}}(s,a)]^2\bigg]
		\end{equation}
		It can be shown that $\omega_*^\theta$ is the exact Natural Policy Gradient (NPG) update direction.
	\end{assumption}
	\begin{remark}
		By  Eq. \eqref{eq:transfer_error} and \eqref{eq:NPG_direction}, the transferred function approximation error expresses an approximation error with distribution shifted to $(d_\rho^{\pi^*},\pi^*)$. With the softmax parameterization or linear MDP structure \cite{Chi2019}, it has been shown that $\epsilon_{bias}=0$ \cite{Alekh2020}. When parameterized by the restricted policy class, $\epsilon>0$ due to $\pi_\theta$ not containing all policies. However, for a rich neural network parameterization, the $\epsilon_{bias}$ is small \cite{Lingxiao2019}. Similar assumption has been adopted in \cite{Yanli2020} and \cite{Alekh2020}. 
	\end{remark}
	\if 0
	\begin{equation}
		\begin{aligned}
		&\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta(\cdot\vert s)}\bigg[[\nabla_\theta\log\pi_{\theta}(a\vert s)\cdot(1-\gamma) \omega_*^\theta-\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_\theta})}{\partial J_m^{\pi_\theta}}A_m^{\pi_\theta}(s,a)][\nabla_\theta\log\pi_\theta(a\vert s)]\bigg]=0\\
		\Rightarrow&\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta}[\nabla_\theta\log\pi_\theta(a\vert s)\nabla_\theta\log\pi_\theta(a\vert s)^T] \omega_*^\theta=\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta}[\frac{1}{1-\gamma}\nabla_\theta\log\pi_\theta(a\vert s)\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_\theta})}{\partial J_m^{\pi_\theta}}A_m^{\pi_\theta}(s,a)]\\
		\Rightarrow&\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta}[\nabla_\theta\log\pi_\theta(a\vert s)\nabla_\theta\log\pi_\theta(a\vert s)^T] \omega_*^\theta=\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta}[\frac{1}{1-\gamma}\nabla_\theta\log\pi_\theta(a\vert s)\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_\theta})}{\partial J_m^{\pi_\theta}}A_m^{\pi_\theta}(s,a)]\\
		\Rightarrow&F_\rho(\theta)^{\dagger} \omega_*^\theta=\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_\theta})}{\partial J_m^{\pi_\theta}}\mathbf{E}_{s\sim d_\rho^{\pi_\theta}}\mathbf{E}_{a\sim\pi_\theta}[\frac{1}{1-\gamma}\nabla_\theta\log\pi_\theta(a\vert s)A_m^{\pi_\theta}(s,a)]\\
		\Rightarrow&F_\rho(\theta) \omega_*^\theta=\sum_{m=1}^{M}\frac{\partial f(\bm{J}^{\pi_\theta})}{\partial J_m^{\pi_\theta}}\nabla_\theta J_m^{\pi_\theta}\Rightarrow 			F_\rho(\theta) \omega_*^\theta=\nabla_\theta f(J_1^{\pi_\theta},\cdots,J_M^{\pi_\theta})\Rightarrow \omega_*^\theta=F_\rho(\theta)^{\dagger}\nabla_\theta f
		\end{aligned}
	\end{equation}
	\fi
	
	\subsection{Global Convergence in Multi-Objective Setting}
	
	
	Inspired by the global convergence analysis framework for policy gradient in \cite{Yanli2020}, we present a general framework for convergence analysis of non-linear multi-objective policy gradient in the following. 
	
	
	\begin{lemma}\label{lem_framework} (Generalization of Proposition 4.5 in \cite{Yanli2020})
		Suppose a general gradient ascent algorithm updates the parameter in the way
		\begin{equation}
			\theta^{k+1}=\theta^k+\eta\omega^k
		\end{equation}
		When Assumptions \ref{ass_score} and \ref{ass_transfer_error} hold, we have
		\begin{equation}\label{eq:general_bound}
			\begin{split}
			f(\bm{J}^{\pi^*})-\frac{1}{K}\sum_{k=0}^{K-1}f(\bm{J}^{\pi_{\theta^k}})&\leq \frac{\sqrt{\epsilon_{bias}}}{1-\gamma}+ \frac{G}{K}\sum_{k=0}^{K-1}\Vert(\omega^k-\omega_*^k)\Vert_2+\frac{M\eta}{2K}\sum_{k=0}^{K-1}\Vert\omega^k\Vert^2\\
			&+\frac{1}{\eta K}\mathbf{E}_{s\sim d_\rho^{\pi^*}}[KL(\pi^*(\cdot\vert s)\Vert\pi_{\theta^0}(\cdot\vert s))]		\end{split}
		\end{equation}
		where $\omega_*^k:=\omega_*^{\theta^k}$ and is defined in Eq. \eqref{eq:NPG_direction}
	\end{lemma}
	
	\begin{proof}
		We generalize the Proposition 4.5 in \cite{Yanli2020} by using the Lemma \ref{lem_performance_diff} and propose the framework of global convergence analysis in the joint optimization for multi-objective setting. Thus, the framework proposed in the Proposition 4.5 in \cite{Yanli2020} can be considered as a special case. The detailed proof is provided in Appendix \ref{sec_app_framwork}.
	\end{proof}
	\noindent Now, we provide the main result of global convergence for the policy gradient algorithm with multi-objective setting (with detailed proof  in Appendix \ref{sec_app_thm}).
	\begin{theorem}\label{thm1}
		For any $\epsilon>0$, in the Policy Gradient Algorithm \ref{alg:PG} with the proposed estimator in Eq. \eqref{eq:est_truncated}, if step-size $\eta=\frac{1}{4L_J}$, the number of iteration $K=\mathcal{O}(\frac{M}{(1-\gamma)^2\epsilon})$, the number of samples $N_1=\mathcal{O}(\frac{\sigma^2}{\epsilon})$ and  $N_2=\mathcal{O}(\frac{M^3}{(1-\gamma)^6\epsilon})$ achieves the following bound
		\begin{equation}
			f(\bm{J}^{\pi^*})-\frac{1}{K}\sum_{k=0}^{K-1}f(\bm{J}^{\pi_{\theta^k}})\leq \frac{\sqrt{\epsilon_{bias}}}{1-\gamma}+\epsilon
		\end{equation}
		In other words, policy gradient algorithm needs $\mathcal{O}\big(\frac{M^4\sigma^2}{(1-\gamma)^8\epsilon^4}\big)$ trajectories.
	\if 0	\begin{proof}
			The detailed proof is provided in Appendix \ref{sec_app_thm}.
		\end{proof}\fi 
	\end{theorem}