text:
\section{Introduction}\label{sec:intro}

Consider an infinite-horizon Markov decision process (MDP)
\cite{bellman1957markovian,sutton2018reinforcement} $\mathcal{M} = (S, A, P, r, \gamma)$, where $S$
is a set of states of the system studied, $A$ is a set of actions made by the agent, $P$ is a
transition probability tensor with $P_{st}^{a}$ being the probability of transitioning from state
$s$ to state $t$ when taking action $a$, $r$ is a reward tensor with $r_{s}^a$ being the reward
obtained when taking action $a$ at state $s$, and $0<\gamma<1$ is a discount factor. Throughout the
paper, the state space $S$ and the action space $A$ are assumed to be finite. A policy $\pi$ is a
randomized rule of action-selection where $\pi_s^a$ denotes the probability of choosing action $a$
at state $s$. For a given policy $\pi$, the value function $v_\pi$ is defined as
\begin{equation}
  (v_\pi)_s = \mathbb{E} \sum_{k=0}^\infty \left(\gamma^k r_{s_k}^{a_k}\mid s_0 = s\right),
\end{equation}
which satisfies the Bellman equation:
\begin{equation}\label{eq:vanillabellman}
  (I-\gamma P_\pi)v_\pi = r_\pi,
\end{equation}
where $(P_\pi)_{st} = \sum_a\pi_s^aP_{st}^a$, $(r_\pi)_{s}=\sum_a\pi_s^ar_{s}^a$, and $I$ is the
identity operator.

In order to promote exploration and enhance stability, one often regularizes the problem with a
function $h_\pi$ such as the negative Shannon entropy $(h_\pi)_s = \sum_a \pi_s^a\log\pi_s^a$. With
the regularization $h_\pi$, the original reward $r_\pi$ is replaced with the regularized reward
$r_\pi-\tau h_\pi$ where $\tau>0$ is a regularization coefficient and \eqref{eq:vanillabellman}
becomes
\begin{equation}
  (I-\gamma P_\pi)v_\pi = r_\pi - \tau h_\pi,
\end{equation}
where we overload the notation $v_\pi$ for the regularized value function. Other continuously
differentiable entropy functions can be used as well, as we will show later. Since $\gamma < 1$ and
$P_\pi$ is a transition probability matrix, $(I-\gamma P_\pi)$ is invertible, and
\begin{equation}\label{eq:regularv}
    v_\pi = (I-\gamma P_\pi)^{-1}(r_\pi-\tau h_\pi).
\end{equation}

In a policy optimization problem, we seek a policy $\pi$ that maximizes the value function $v_\pi$. According to the theory of regularized MDPs \cite{geist2019theory}, there is a unique optimal policy $\pi^*$ such that $(v_{\pi^*})_s\geq (v_{\pi})_s$ for any policy $\pi$ and state $s$. It thus suffices to maximize $e^\top v_\pi$ for any positive weight vector $e\in\mathbb{R}_+^{|S|}$. Using \eqref{eq:regularv}, the problem can be stated as
\begin{equation}\label{optim}
    \underset{\pi}{\max} ~ e^\top(I-\gamma P_\pi)^{-1}(r_\pi-\tau h_\pi).
\end{equation}
This problem can be solved by, for example, the policy gradient (PG) method.  However, the vanilla
PG method converges quite slowly. In \cite{agarwal2020optimality}, for instance, the vanilla PG
method is shown to have the $O(T^{-1})$ convergence rate, where $T$ denotes the number of
iterations. For the PG method with entropy regularization and some of its variants, the convergence
rate can be improved to $O(e^{-cT})$, i.e., linear convergence \cite{mei2020global}, which can
still be slow since $c$ is in general close to $0$.  It is also demonstrated in numerical examples
that these algorithms with linear convergence rate can suffer from slow convergence.  For example,
in the example in \cite{zhan2021policy}, thousands of iterations are needed for the algorithm to
converge, even though the model is relatively small and sparse. Therefore, there is a clear need for
designing new methods with faster convergence and one idea is to take the geometry of the problem
into consideration. The Newton method, for example, preconditions the gradient with the Hessian
matrix and obtains second-order local convergence. Since the exact Hessian matrix is usually too
computationally expensive to obtain, the quasi-Newton type methods, which use structurally simpler
approximations of the Hessian instead, are more widely used in generic optimization problems, and
are known to enjoy superlinear convergence \cite{rodomanov2021new, rodomanov2021rates}.

\subsection{Contributions}
In this paper, we investigate the quasi-Newton approach for solving \eqref{optim}. The main
contributions of this paper are the following.
\begin{itemize}
\item First, we present a unified quasi-Newton type method for the policy optimization problem.  The
  main observation is to decompose the Hessian as a sum of a diagonal matrix and a remainder that vanishes at
  the optimal solution. This inspires us to use only the diagonal matrix in the quasi-Newton method. As a result, the proposed method not only leverage the second order information but also enjoys low computational cost due to the diagonal structure of the preconditioner used. When the negative Shannon entropy is used,
  this method reproduces the natural policy gradient (NPG) algorithm. For other forms of entropy regularization, this method results in brand new policy gradient algorithms.

\item Second, we analyze the convergence property of the proposed quasi-Newton algorithms and
  demonstrate local quadratic convergence both theoretically and numerically. By leveraging the
  quasi-Newton framework \cite{dennis1974characterization}, we provide a simple and straightforward
  proof for quadratic convergence near the optimal policy. In the numerical tests, we verify that
  the proposed method leads to fast quadratic convergence even under small regularization and large
  discount rate (close to 1). Even for industrial-size problems with hundreds of thousands states,
  the quasi-Newton method converges in single-digit iterations and within a few minutes on a regular
  laptop. We also provide a proof for the global convergence of the quasi-Newton gradient flow to the optimal solutions.
\end{itemize}

\subsection{Background and related work}
A major workhorse behind the recent success in reinforcement learning (RL) is the large family of
policy gradient (PG) methods \cite{williams1992simple, sutton1999policy}, for
example, the natural policy gradient (NPG) method \cite{kakade2001natural}, the actor-critic method \cite{konda2000actor}, the asynchronous advantage actor-critic (A3C) method
\cite{mnih2016asynchronous}, the deterministic policy gradient (DPG) method
\cite{silver2014deterministic}, the trust region policy optimization (TRPO)
\cite{schulman2015trust}, the generalized advantage estimation (GAE) \cite{schulman2015high}, and
proximal policy optimization (PPO) \cite{schulman2017proximal}, to mention but a few. The NPG method is known to be drastically faster than the original PG method, intuitively because in the NPG method the policy gradient is preconditioned by the Fisher information (an approximation of the Hessian of the KL-divergence) matrix and fits the problem geometry better. This idea is extended in TRPO and PPO where the problem geometry is taken into consideration via trust region constraints (in terms of KL-divergence) and a clipping function of the relative ratio of policies in the objective function, respectively. These implicit ways (in the sense that they do not adjust the gradient by an explicit preconditioner) of adjusting the policy gradient is in essence similar to the mirror descent (MD) method \cite{nemirovskij1983problem} in generic optimization problems.

This similarity in addressing the inherent geometry of the problem is noticed by a line of recent
work including \cite{neu2017unified, geist2019theory, shani2020adaptive, tomar2020mirror, lan2021policy}, and the analysis techniques in MD methods have been adapted to the PG setting. The connection was first built explicitly in \cite{neu2017unified}. The authors consider a linear program formulation where the objective function is the average reward and the domain is the set of stationary state-action distributions, in which case the TRPO method can be viewed as an approximate mirror descent method and the A3C method as an MD method for the dual-averaging \cite{nesterov2009primal} objective. As a complement, \cite{geist2019theory} considers an actor-critic type method where the policy is updated via either a regularized greedy step or an MD step and the value function is updated by a regularized Bellman operator, which also includes TRPO as a special case, and error propagation analysis is provided. In \cite{shani2020adaptive}, an adaptive scaling that naturally arises in the policy gradient is applied to the proximity term of the MD formulation, and sublinear convergence result is proved with a properly decreasing learning rate. In \cite{tomar2020mirror}, the application to the non-tabular setting is enabled by parameterizing the policy and applying MD to the policy parameters, and the corresponding sublinear convergence result is presented.

Regularization, a strategy that considers the modified objective function with an additional penalty term on the policy, is another crucial component in the development of PG type methods. Intuitively, regularization is able to encourage exploration in the policy iteration process and thus avoid local minima. It is also suggested \cite{ahmed2019understanding} that regularization makes the optimization landscape smoother and thus enables possibly faster convergence. Linear convergence result is then established for regularized PG and NPG methods \cite{agarwal2020optimality, mei2020global, cen2020fast}. In these relatively earlier works \cite{agarwal2020optimality, mei2020global, cen2020fast}, the regularization usually takes the form of (negative) entropy or relative entropy. In the more recent work \cite{lan2021policy} and \cite{zhan2021policy} that follow the MD type methods, the regularization is extended to general convex functions with the resulting Bregman divergences different from the KL-divergence and linear convergence is guaranteed as well.

However, most of these algorithms are of either sublinear or linear convergence except the entropy
regularized NPG with full step length (which is a special case of the quasi-Newton method we
propose), and even the linear convergence rate $O(e^{-cT})$ can be slow since $c$ can be close to
zero. This motivates us to invent the quasi-Newton policy gradient method to be introduced in \cref{sec:Newton}.

\section{Quasi-Newton method}\label{sec:Newton}
\subsection{Quasi-Newton method and entropy regularized natural policy gradient}
\begin{sloppypar}
This section derives the quasi-Newton method for the entropy regularized policy optimization
problems. In this paper, we use a more general definition of the term ``quasi-Newton
method''. Unlike classical quasi-Newton methods where the approximate Hessian matrix is constructed
using first-order information, the proposed quasi-Newton method uses second-order information, but
the approximate Hessian matrix remains simple and its inverse is easy to obtain.  We start with the
negative Shannon entropy $(h_\pi)_s = \sum_a \pi_s^a\log\pi_s^a$.
\end{sloppypar}
In what follows, it is assumed that $\pi^*$ is the optimizer of the problem stated in \eqref{optim}.  By
introducing $Z_\pi := I-\gamma P_\pi$, the objective function can be written as
\begin{equation}
  E(\pi) \equiv e^\top(I-\gamma P_\pi)^{-1}(r_\pi-\tau h_\pi) = e^\top Z_\pi^{-1}(r_\pi-\tau h_\pi) = w_\pi^\top(r_\pi-\tau h_\pi),
\end{equation}
where $w_\pi := Z_\pi^{-\top}e$. For any $\epsilon\in \mathbb{R}^{|S|\times|A|}$, introduce
$r_\epsilon\in \mathbb{R}^{|S|}$ and $Z_\epsilon \in \mathbb{R}^{|S|\times|S|}$ such that
\begin{equation}\label{eq:rZexpansion}
(r_\epsilon)_s :=\sum_a\epsilon_s^ar_{s}^a, \quad
(Z_\epsilon)_{st}:=\sum_a\epsilon_s^a(\delta_{st}-\gamma P_{st}^a),
\end{equation}
where $\delta_{st} = 1$ if $s=t$ and $\delta_{st} = 0$ otherwise. Then $r_\epsilon$ and $Z_\epsilon$ are linear with respect to $\epsilon$, which is helpful when expressing the first-order conditions and simplifying the expansion of $E(\pi)$.

Let us first outline the main idea of the quasi-Newton method. The gradient $\nabla_\pi E$ in $\mathbb{R}^{|S||A|}$ of $E(\pi)$ has entries given by
\begin{equation}
  \frac{\partial E}{\partial \pi_s^a} = (r_s^a - \tau (\log\pi_s^a +1) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s)(w_\pi)_s,
\end{equation}
where $c_s$ is a multiplier associated with the constraint $\sum_a\pi_s^a = 1$ that depends on $s$. Our key observation is to decompose the Hessian matrix $D^2E(\pi)$ in $\mathbb{R}^{|S||A|\times |S||A|}$ into two parts
\begin{equation}\label{eq:decomp}
D^2E(\pi) = \Lambda(\pi) + \Delta(\pi),
\end{equation}
where $\Lambda(\pi)$ is a {\em diagonal} matrix given by
$\Lambda_{(sa),(tb)}=-\tau\delta_{\{(sa),(tb)\}}\frac{(w_\pi)_s}{\pi_s^a}$ and $\Delta(\pi)$ is a
remainder that {\em vanishes} at $\pi = \pi^*$, i.e., $\Delta(\pi)=O(\|\pi-\pi^*\|)$ (shown in \Cref{thm: approxnewton}). We emphasize that $\Lambda(\pi)$ is in general not the diagonal part of the Hessian matrix $D^2E(\pi)$, but a diagonal approximation to it. With this decomposition, we can approximate the Hessian matrix $D^2E(\pi)$
by $\Lambda(\pi)$ and obtain the following {\em quasi-Newton flow}:
\[
\begin{aligned}
  \frac{\mathrm{d}\pi_s^a}{\mathrm{d} t} &=-(\Lambda^{-1}\nabla_\pi E)_{sa}= -(\Lambda_{(sa),(sa)})^{-1}\frac{\partial E}{\partial \pi_s^a}\\
  &= \pi_s^a(r_s^a-\tau(\log\pi_s^a+1)-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau,
\end{aligned}
\]
By introducing the parameterization $\theta_s^a = \log \pi_s^a$ and discretizing in time with learning rate $\eta$, we arrive at
\[
\theta_s^a \leftarrow \eta(r_s^a-\tau-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau + (1-\eta)\theta_s^a.
\]
Writing this update back in terms of $\pi_s^a$ leads to the following update rule
\[
\pi_s^a \propto (\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma
  P^av_\pi)_s)/\tau),
\]
which coincides with the NPG scheme with entropy regularization. This result is summarized in the following theorem with the proof given in \cref{sec:approxnewton}.

\begin{theorem}\label{thm: approxnewton}
  Let $h_\pi \in \mathbb{R}^{|S|}$ be the negative Shannon entropy $(h_\pi)_s=\sum_a\pi_s^a\log\pi_s^a$.

  (a) For any $\epsilon$ with $\sum_a\epsilon_s^a = 0$ and $|\epsilon_s^a|<\pi_s^a$, at $\pi=\pi^*$
  \begin{equation}\label{eq:grad}
    r_\epsilon-\tau Dh_\pi\epsilon-Z_\epsilon Z_\pi^{-1}(r_\pi-\tau h_\pi) = 0,
  \end{equation}
  where $Dh_\pi \in \mathbb{R}^{|S|\times |S||A|}$ is the gradient matrix of $h_\pi$ with respect to
  $\pi$.

  (b) There exists a diagonal approximation $\Lambda(\pi)$ of the Hessian matrix $D^2E(\pi)$ given by $\Lambda_{(sa),(tb)}=-\tau\delta_{\{(sa),(tb)\}}\frac{(w_\pi)_s}{\pi_s^a}$ such that
  \begin{equation}\label{eq:approxH}
      \Lambda(\pi)-D^2E(\pi) = O(\|\pi-\pi^*\|).
  \end{equation}

  (c) The quasi-Newton flow from $\Lambda(\pi)$ is
  \begin{equation}\label{eq:QNflow}
    \frac{\mathrm{d}\pi_s^a}{\mathrm{d} t} = \pi_s^a(r_s^a-\tau(\log\pi_s^a+1)-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau.
  \end{equation}
  With a learning rate $\eta$, the gradient update is
  \begin{equation}\label{eq:NPG}
    \pi_s^a \leftarrow \frac{(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma
      P^av_\pi)_s)/\tau)}{\sum_a(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma P^av_\pi)_s)/\tau)}.
  \end{equation}
\end{theorem}

\begin{remark}
  The policy update scheme \eqref{eq:NPG} is the same as the entropy regularized natural policy
  gradient scheme in, for example, \cite{cen2020fast}.
\end{remark}

\paragraph{Historical note}
The natural gradient methods (including the NPG method) were traditionally developed as a way of
implementing the vanilla gradient descent method with an intrinsic metric that is invariant to the
choice of parameters (Cf. \cite{martens2020new}), and entropy regularization was originally
motivated as a way of encouraging exploration and avoid the suboptimality caused by greedy solvers
(Cf. \cite{neu2017unified}). In this regard, it was more or less a coincidence that the algorithm
combining the two methods -- the regularized NPG obtains a fast quadratic convergence
(Cf. \cite{cen2020fast}). The reason behind this coincidence is that the preconditioner used in the
natural gradient method in fact approximates the second order derivatives introduced by the entropy
regularization in this case, though the fisher information matrix was not designed to approximate
any second-order information in the classical natural gradient literature.

\subsection{Extension to other entropy functions}
\Cref{thm: approxnewton} can be extended to more general entropy functions. It yields brand-new
algorithms with quadratic convergence. Here we consider the entropy functions of the form
\begin{equation}\label{eq:generaletp}
(h_\pi)_s = \sum_a\phi\left(\frac{\pi_s^a}{\mu_s^a}\right)\mu_s^a,
\end{equation}
where $\phi$ is convex on $(0, +\infty)$ and $\phi(1) = 0$, and $\mu_s$ is a prior distribution over $A$ such that $\mu_s^a>0$. The term $(h_\pi)_s$ is also called the ``$f$-divergence'' between $\pi_s$ and $\mu_s$ \cite{renyi1961measures,ali1966general}. If there is no prior knowledge on the policy, one can use the uniform prior, i.e., $\mu_s^a=1/|A|$ for all $a$.
We further assume that $\phi$ is twice continuously differentiable and strongly convex and that $\phi'(x)\rightarrow-\infty$ as $x\rightarrow 0$. Here are some examples:
\begin{itemize}
\item When $\phi(x) = x\log x$, $(h_\pi)_s = \sum_a
  \left(\frac{\pi_s^a}{\mu_s^a}\log\frac{\pi_s^a}{\mu_s^a}\right)\mu_s^a
  =\sum_a\pi_s^a\log\frac{\pi_s^a}{\mu_s^a}$. When the uniform prior is used, we recover the
  negative Shannon entropy regularization $\sum_a \pi_s^a\log\pi_s^a$ used in \Cref{thm:
    approxnewton} after omitting the constant $\log\frac{1}{|A|}$.
\item When $\phi(x) = \frac{4}{1-\alpha^2}(1-x^{\frac{1+\alpha}{2}})~
  (\alpha<1)$, we obtain the $\alpha$-divergence:
  \begin{equation}\label{eq:alpha}
    (h_\pi)_s = \frac{4}{1-\alpha^2}-\frac{4}{1-\alpha^2}\sum_a\mu_s^a\left(\pi_s^a/\mu_s^a\right)^{\frac{1+\alpha}{2}}.
  \end{equation}
  In particular, when $\alpha=0$ we obtain the Hellinger divergence
  $(h_\pi)_s=2-2\sum_a\sqrt{\mu_s^a\pi_s^a}$ after dividing by $2$. When
  $\alpha\rightarrow-1$ we obtain the reverse-KL divergence $(h_\pi)_s =   \sum_a\mu_s^a\log\frac{\mu_s^a}{\pi_s^a}$. Also, when $\alpha\rightarrow1$, we obtain the KL-divergence $(h_\pi)_s=\sum_a\pi_s^a\log\frac{\pi_s^a}{\mu_s^a}$, though the limit of $\phi(x)$ does not exist when $\alpha\rightarrow1$.
\end{itemize}

In the following theorem, we extend the quasi-Newton method in \Cref{thm: approxnewton} to the entropy functions described above. The proof of this theorem can be found in \cref{sec:otheretp}.

\begin{theorem}\label{thm: otheretp}
  Assume that $\pi^*$ is the optimizer of \eqref{optim} where $h_\pi$ is the entropy function defined in \eqref{eq:generaletp}.

  (a) For any $\epsilon$ with $\sum_a\epsilon_s^a = 0$ and $|\epsilon_s^a|<\pi_s^a$, at $\pi=\pi^*$
  \begin{equation}\label{eq:general1stmain}
    r_\epsilon-\tau Dh_\pi\epsilon-Z_\epsilon Z_\pi^{-1}(r_\pi-\tau h_\pi) = 0,
  \end{equation}

  (b) The Hessian matrix $D^2E(\pi)$ can be approximated by a diagonal matrix $\Lambda(\pi)$ given by
  \begin{equation}\label{eq:generalapproxhess}
    \Lambda_{(sa),(tb)} = -\tau\delta_{\{(sa),(tb)\}}\frac{(w_\pi)_s\phi''(\pi_s^a/\mu_s^a)}{\mu_s^a}
  \end{equation}
  near $\pi^*$ such that $\Lambda(\pi)-D^2E(\pi) = O(\|\pi-\pi^*\|)$.

  (c) The quasi-Newton flow from $\Lambda$ is
  \begin{equation}\label{eq:qNflowgen}
  \frac{\mathrm{d}\pi_s^a}{\mathrm{d} t} = \mu_s^a(\phi''(\pi_s^a/\mu_s^a))^{-1}(r_s^a- \tau \phi'(\pi_s^a/\mu_s^a)-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau.
  \end{equation}
  With parameterization $\theta_s^a = \phi'(\pi_s^a/\mu_s^a)$, the quasi-Newton method from $\Lambda(\pi)$ can be expressed as:
  \begin{equation}\label{eq:quasiNTupdate}
    \theta_s^a \leftarrow \eta(r_s^a-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau + (1-\eta)\theta_s^a.
  \end{equation}
  where where $0<\eta\leq1$ is the learning rate and $c_s$ is a multiplier introduced by the constraint $\sum_a\pi_s^a = 1$.

\end{theorem}

For particular choices of $\phi$, the corresponding quasi-Newton update scheme can be obtained directly by plugging $\phi$ into \eqref{eq:quasiNTupdate}.
\begin{itemize}
    \item For the case $\phi(x) = x\log x$ and $(h_\pi)_s=\sum_a\pi_s^a\log\frac{\pi_s^a}{\mu_s^a}$,
one can solve the multipliers $c_s$ explicitly as in \Cref{thm: approxnewton} and obtain the NPG method with prior distribution $\mu$:
\begin{equation}
  \pi_s^a \leftarrow \frac{(\mu_s^a)^\eta(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma P^av_\pi)_s)/\tau)}{\sum_a(\mu_s^a)^\eta(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma P^av_\pi)_s)/\tau)}.
\end{equation}
    \item For the case $\phi(x) = \frac{4}{1-\alpha^2}(1-x^{\frac{1+\alpha}{2}})~(\alpha<1)$ and $h_\pi$ given by the $\alpha$-divergence \eqref{eq:alpha}, we have $\theta_s^a = -\frac{2}{1-\alpha}(\pi_s^a/\mu_s^a)^{\frac{\alpha-1}{2}}$, thus by \eqref{eq:quasiNTupdate} the update scheme is
\begin{equation}\label{eq:alphaupdate}
    \pi_s^a \leftarrow \mu_s^a\left((1-\eta)(\pi_s^a/\mu_s^a)^{\frac{\alpha-1}{2}} + \frac{\alpha-1}{2}\eta(r_s^a-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau\right)^{\frac{2}{\alpha-1}}.
\end{equation}
\end{itemize}

The remaining problem in the update schemes is the determination of the multipliers $c_s$,
since in general they cannot be solved explicitly as in the case of the negative Shannon entropy (Cf. \Cref{thm: approxnewton}). Since $\phi$ is strongly convex, we know that $\phi'$ is strictly increasing, and thus $-\phi'$ is a strictly decreasing function mapping from $(0, +\infty)$ to $(-\sup\phi', +\infty)$ since $\underset{x\rightarrow0+0}{\lim}\phi'(x)=-\infty$. Let $\psi:=(-\phi')^{-1}$, then $\psi: (-\sup\phi', +\infty)\rightarrow (0, +\infty) $ is a strictly decreasing function that satisfies $\underset{x\rightarrow-\sup\phi'+0}{\lim}\psi(x)=+\infty$ and $\underset{x\rightarrow+\infty}{\lim}\psi(x)=0$. From \eqref{eq:quasiNTupdate}, the equation of the multiplier $c_s$ corresponding to $\sum_a\pi_s^a=1$ is:
\begin{equation}\label{eq:eqmultorig}
    \sum_a\mu_s^a\psi\left(-\frac{\eta}{\tau}c_s-(1-\eta)\phi'\left(\pi_s^a/\mu_s^a\right) -\frac{\eta}{\tau}(r_s^a-[(I-\gamma P^a)v_\pi]_s)\right)=1,
\end{equation}
or equivalently,
\begin{equation}\label{eq:eqmult}
    \sum_a\mu_s^a\psi\left(\tilde{c}_s+x_a\right)=1,
\end{equation}
where
\begin{equation}
    \label{eq:deftildecxa}
    \tilde{c}_s=-\frac{\eta}{\tau}c_s, \quad
    x_a = -(1-\eta)\phi'\left(\pi_s^a/\mu_s^a\right) -\frac{\eta}{\tau}(r_s^a-[(I-\gamma P^a)v_\pi]_s).
\end{equation}
We claim that the determination of $\tilde{c}_s$ in equation \eqref{eq:eqmult} (and thus the determination of $c_s$) can be done in a similar way as in \cite{ying2020mirror} based on the following lemma. The proof of this lemma can be found in \cref{sec:mult}.

\begin{lemma}\label{lemma:multipliergen}
  Let $L\in\mathbb{R}\cup\{-\infty\}$. Assume that $\psi: (L, +\infty)\rightarrow (0, +\infty) $ is a strictly decreasing function that satisfies $\underset{x\rightarrow L+0}{\lim}\psi(x)=+\infty$ and $\underset{x\rightarrow+\infty}{\lim}\psi(x)=0$ and assume also that $\mu_i>0$, then for any $x_1, x_2, \ldots, x_k$, there is a unique solution to the equation:
  \eq{\label{eq:lemmagen}
    \mu_1\psi(x+x_1) + \cdots + \mu_k\psi(x+x_k) = 1.
  }
  Moreover, the solution is on the interval
  \begin{equation}\label{bisect_intervalgen}
  \left[\max\left\{L-\underset{1\leq i\leq k}{\min}x_i, \underset{1\leq i\leq k}{\min}\left\{\psi^{-1}\left(\frac{1}{k\mu_i}\right) - x_i\right\}\right\}, \underset{1\leq i\leq k}{\max}\left\{\psi^{-1}\left(\frac{1}{k\mu_i}\right) - x_i\right\}\right].
  \end{equation}
\end{lemma}

Leveraging \Cref{lemma:multipliergen} and the monotonicity of the function $\mu_1\psi(x+x_1) + \cdots + \mu_k\psi(x+x_k) - 1$, many of the established numerical methods (e.g. bisection) for nonlinear equations can be applied to determine the solution for \eqref{eq:lemmagen}. This routine can be used to find $\tilde{c}_s$ in \eqref{eq:eqmult} and thus the multipliers $c_s$ in \eqref{eq:eqmultorig} as stated in \Cref{prop:multiplier} whose proof is given in \cref{sec:multsolve}.

\begin{prop}
\label{prop:multiplier}
  The multipliers $c_s$ in the update scheme \eqref{eq:quasiNTupdate} can be determined uniquely
  such that the updated policy $\pi$ satisfies $\pi_s^a\geq0$ and $\sum_a\pi_s^a =1$.
\end{prop}

When the $\alpha$-divergence is used, we have $\phi=\frac{4}{1-\alpha^2}(1-x^{\frac{1+\alpha}{2}})$ and $\phi'(x) = \frac{2}{\alpha-1}x^{\frac{\alpha-1}{2}}$, then $L=-\sup\phi'=0$ and $\psi(x) = (-\phi')^{-1}(x) = (\phi')^{-1}(-x) = (\frac{1-\alpha}{2}x)^{\frac{2}{\alpha-1}}$. The algorithm proposed in this section is summarized in \Cref{alg:quasi} below.
\begin{algorithm}[htbp]
    \caption{Quasi-Newton method for the regularized MDP}
    \label{alg:quasi}
    \begin{algorithmic}[1]
        \REQUIRE the MDP model $\mathcal{M} = (S, A, P, r, \gamma)$, initial policy $\pi_{\text{init}}$, convergence threshold $\epsilon_{\text{tol}}$,
        regularization coefficient $\tau$, learning rate $\eta$, the regularization type (KL or $\alpha$-divergence).
        \STATE Initialize the policy $\pi = \pi_{\text{init}}$.
        \STATE Set $\xi = 1 + \epsilon_{\text{tol}}$ and $k = |A|$.

        \WHILE {$\xi > \epsilon_{\text{tol}}$}

        \STATE Calculate the regularization term $h_\pi$ by $(h_\pi)_s = \sum_a\mu_s^a\phi(\pi_s^a/\mu_s^a)$.
        \STATE Calculate $P_\pi$ and $r_\pi$ by $(P_\pi)_{st} = \sum_a\pi_s^aP_{st}^a, (r_\pi)_{s} = \sum_a\pi_s^ar_{s}^a$.
        \STATE Calculate $v_\pi$ by \eqref{eq:regularv}, i.e., $v_\pi = (I-\gamma P_\pi)^{-1}(r_\pi-\tau h_\pi)$.

        \IF {the KL divergence is used}
        \STATE $(\pi_{\text{new}})_s^a \leftarrow \frac{(\mu_s^a)^\eta(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma P^av_\pi)_s)/\tau)}{\sum_a(\mu_s^a)^\eta(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma P^av_\pi)_s)/\tau)}$ for $a = 1, 2, \ldots, |A|, ~ s = 1, 2, \ldots, |S|$.
        \ENDIF

        \IF {the $\alpha$-divergence is used}

        \FOR {$s = 1, 2, \ldots, |S|$}

        \STATE Set $L=0$ and $\psi(x) = (\frac{1-\alpha}{2}x)^{\frac{2}{\alpha-1}}$
        \STATE Calculate $x_a = -(1-\eta)\phi'\left(\pi_s^a/\mu_s^a\right) -\frac{\eta}{\tau}(r_s^a-[(I-\gamma P^a)v_\pi]_s), ~ a = 1, \ldots, |A|$.

        \STATE Solve for $\tilde{c}_s=-\frac{\eta}{\tau}c_s$ with the bisection method on the interval described in \eqref{bisect_intervalgen}.
        \STATE Update $(\pi_{\text{new}})_s^a \leftarrow \mu_s^a \psi(\tilde{c}_s+x_a)$ for $a = 1, 2, \ldots, |A|$.
        \ENDFOR
        \ENDIF

        \STATE $\xi = \|\pi_{\text{new}}-\pi\|_F / \|\pi\|_F$.
        \STATE $\pi = \pi_{\text{new}}$

        \ENDWHILE
    \end{algorithmic}
    \vspace{0.5em}
\end{algorithm}
\subsection{Convergence of the quasi-Newton gradient flow}\label{sec:flowconv}
Recall from \eqref{eq:qNflowgen} that the quasi-Newton gradient flow with the general entropy functions is
\[
\frac{\mathrm{d}\pi_s^a}{\mathrm{d} t} = \mu_s^a(\phi''(\pi_s^a/\mu_s^a))^{-1}(r_s^a- \tau \phi'(\pi_s^a/\mu_s^a)-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau,
\]
from which we can obtain the dynamics of the objective function $E$:
\begin{equation}
\begin{aligned}
\frac{\mathrm{d}E(\pi)}{\mathrm{d} t} &= \sum_{sa} \frac{\partial E}{\partial \pi_s^a}\frac{\mathrm{d}\pi_s^a}{\mathrm{d} t}\\
&= \sum_{sa}\bigg[(r_s^a - \tau \phi'(\pi_s^a/\mu_s^a) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s)(w_\pi)_s\\
&\;\;\; \cdot\mu_s^a(\phi''(\pi_s^a/\mu_s^a))^{-1}(r_s^a- \tau \phi'(\pi_s^a/\mu_s^a)-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau \bigg]\\
&= \sum_{sa}\mu_s^a(\tau\phi''(\pi_s^a/\mu_s^a))^{-1}(r_s^a- \tau \phi'(\pi_s^a/\mu_s^a)-[(I-\gamma P^a)v_\pi]_s+c_s)^2(w_\pi)_s\\
&\geq 0,
\end{aligned}
\end{equation}
where we have used the gradient
\begin{equation}
    \frac{\partial E}{\partial \pi_s^a} = (r_s^a - \tau \phi'(\pi_s^a/\mu_s^a) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s)(w_\pi)_s.
\end{equation}
As a result, we have shown that $\frac{\mathrm{d}E(\pi)}{\mathrm{d} t}\geq0$. Since $E(\pi)$ is upper-bounded by $e^\top v_{\pi^*}$, $E(\pi)$ converges. With a closer look, the following theorem states that the limiting policy is exactly the optimal policy $\pi^*$ and the proof is given in \cref{sec:flowconvpf}.
\begin{theorem}\label{thm:flowconv}
The quasi-Newton flow \eqref{eq:qNflowgen} converges globally to the optimal policy $\pi^*$.
\end{theorem}

\section{Quadratic Convergence}\label{sec:analysis}
In this section, we study the quadratic convergence of the quasi-Newton method at learning rate
$\eta=1$, which corresponds to the step size used in the Newton method. Our analysis is inspired by the results in \cite{dennis1974characterization,
  wang2021hessian}. The following theorem states the second-order convergence when $\eta=1$, with
the proof given in \cref{sec:quadconvpf}.

\begin{theorem}\label{thm:quadconv}
Let $f(\pi)_{sa} = -(r_s^a-((I-\gamma P^a)v_\pi)_s)$ and
\[
\Phi(\pi) = \tau\sum_{sa}\mu_s^a\phi\left(\frac{\pi_s^a}{\mu_s^a}\right),
\]
where $\phi$ is twice continuously differentiable and strongly convex and that $\phi'(x)\rightarrow-\infty$ as $x\rightarrow 0$.
Denote the $k$-th policy obtained in \Cref{alg:quasi} by $\pi\sps{k}$. For $\eta=1$, the update scheme in \Cref{alg:quasi} can be summarized as
  \begin{equation}\label{eq:Phiupdate}
  \nabla\Phi(\pi\sps{k+1})-\nabla\Phi(\pi\sps{k}) = -\left(f(\pi\sps{k}) + \nabla\Phi(\pi\sps{k}) -
  B^\top c(\pi\sps{k})\right),
  \end{equation}
where we denote by $B$ the $|S|$-by-$(|S|\times|A|)$ matrix such that $B_{ij} = 1$ for $|A|(i-1)+1\leq j\leq|A|i$ and $B_{ij}=0$ otherwise. Assume that $\theta\sps{k}\rightarrow\theta^*$, where $(\theta\sps{k})_s^a = \phi'((\pi\sps{k})_s^a/\mu_s^a)$, then $\pi\sps{k}$ converges quadratically to $\pi^*$, i.e.,
  \begin{equation}\label{eq:quadconv}
    \norm{\pi\sps{k+1}-\pi^*} \leq C\norm{\pi\sps{k}-\pi^*}^2,
  \end{equation}
for some constant $C$.
\end{theorem}

\paragraph{Connection with mirror descent}
The quasi-Newton algorithm \eqref{eq:Phiupdate} for $\eta=1$ has a deep connection with mirror
descent. The vanilla mirror descent of $-E(\pi)$ with a learning rate $\beta$ and the Bregman
divergence associated with $\Phi$ is given by
\eqs{
  \pi\sps{k+1} &= \underset{\pi}{\arg\min}\{-E(\pi\sps{k})-\nabla E(\pi\sps{k})(\pi-\pi\sps{k}) + \frac{1}{\beta}(\Phi(\pi)-\Phi(\pi\sps{k})-\nabla\Phi(\pi\sps{k})(\pi-\pi\sps{k}))\}\\
    &= \underset{\pi}{\arg\min}\{(\mathrm{diag}(w_{\pi\sps{k}})\otimes I_{|A|})(f(\pi\sps{k})+\nabla\Phi(\pi\sps{k}))(\pi-\pi\sps{k}) + \frac{1}{\beta}(\Phi(\pi)-\nabla\Phi(\pi\sps{k})\pi)\},
}
where $\mathrm{diag}(w_{\pi\sps{k}})$ is the diagonal matrix with the diagonal equal to
$w_{\pi\sps{k}} := (I-\gamma P_{\pi\sps{k}}^\top)^{-1}e$, $\otimes$ denotes the Kronecker product,
and $I_{|A|}$ denotes the $|A|$ by $|A|$ identity matrix. In the last equality, the terms
independent of $\pi$ are dropped and the multiplier term in $\nabla E$ is canceled out using $B\pi =
B\pi\sps{k} = \mathbf{1}_{|S|}$. The first order stationary condition of this minimization problem
reads
\eq{
    \nabla\Phi(\pi\sps{k+1}) - \nabla\Phi(\pi\sps{k}) = -\beta(\mathrm{diag}(w_{\pi\sps{k}})\otimes I_{|A|})(f(\pi\sps{k}) + \nabla\Phi(\pi\sps{k}) - B^\top c(\pi\sps{k})).
}
This suggests that \eqref{eq:Phiupdate} can be reinterpreted as an {\em accelerated} mirror descent method with {\em adaptive} learning rates $\beta_s \equiv 1/(w_{\pi\sps{k}})_s$ that depend on the state $s$ and the current iterate $\pi\sps{k}$.

In \cite{zhan2021policy}, a variant of mirror descent is proposed based on an implicit update scheme
\eq{
    (\nabla\Phi(\pi\sps{k+1}))_{sa} - (\nabla\Phi(\pi\sps{k}))_{sa} = -\beta' \left(f(\pi\sps{k})_{sa} + \nabla\Phi(\pi\sps{k+1})_{sa} - (c(\pi\sps{k}))_s\right),
}
with a manually tuned learning rate $\beta'$. In the next section, we will compare this variant with our quasi-Newton method \eqref{eq:Phiupdate} and show that the quasi-Newton method converges orders of magnitudes faster than the ones in \cite{zhan2021policy}.

\section{Numerical results}\label{sec:numerical}

\subsection{Experiment I}
\label{sec:exptoy}

We first test the quasi-Newton methods derived in \cref{sec:Newton} on the model in
\cite{zhan2021policy}. For the sake of completeness, we include the model description here. The MDP considered has a state space $S$ of size $200$ and an action space $A$ of size $50$. For each state $t$ and action $a$, a subset $S_{t}^a$ of $S$ is uniformly randomly chosen such that $|S_t^a| = 20$, and $P_{tt'}^a = 1/20$ for any $t'\in S_t^a$. The reward is given by $r_s^a = U_s^aU_s$, where $U_{s}^a$ and $U_s$ are independently uniformly chosen on $[0,1]$. The discount rate $\gamma$ is set as $0.99$ and the regularization coefficient $\tau=0.001$.

In the numerical experiment, we implement \Cref{alg:quasi} with the KL divergence, the
reverse KL divergence, the Hellinger divergence and the $\alpha$-divergence with $\alpha=-3$. We adopt the uniform prior $\mu_s^a=1/|A|$ in order to make a fair comparison with the policy mirror descent (PMD) and the general policy mirror descent (GPMD) method in \cite{zhan2021policy}. We set the initial policy as the uniform policy, the convergence threshold as $\epsilon = 10^{-12}$ and the learning rate $\eta$ as $1$. \Cref{fig:toy1} demonstrates that, for these four tests, the quasi-Newton algorithm converges in $7$, $7$, $7$, and $6$ iterations, respectively. In comparison, we apply PMD and GPMD to the same MDP with the same stopping criterion. As also shown in
\Cref{fig:toy1}, many more iterations are needed for GPMD and PMD to reach the same precision:
GPMD converges in $14822$ iterations, and PMD does not reach the desired precision after
$3\times10^5$ iterations. For the implementation of GPMD and PMD, we have already tuned the hyperparameters to optimize their performance.

\begin{figure}[ht]
  \centering
  \subfigure[Relative change of the policy using \Cref{alg:quasi} and methods from \cite{zhan2021policy}.\label{fig:toy1}]{
    \includegraphics[width=0.3\textwidth]{images/model1_conv_in_one} }
  \hspace{0.7em}
  \subfigure[The policy error in the process of training using KL-divergence. \label{fig:model1_KL}]{
    \includegraphics[width=0.3\textwidth]{images/model1_KL_new} }
  \hspace{0.7em}
  \subfigure[The policy error in the process of training using reverse KL-divergence. \label{fig:model1_rKL}]{
    \includegraphics[width=0.3\textwidth]{images/model1_rKL_new} }

  \subfigure[The policy error in the process of training using Hellinger divergence. \label{fig:model1_Hellinger}]{
    \includegraphics[width=0.3\textwidth]{images/model1_Hellinger_new} }
  \hspace{0.7em}
  \subfigure[The policy error in the process of training using $\alpha$-divergence with $\alpha=-3$. \label{fig:model1_alpha}]{
    \includegraphics[width=0.3\textwidth]{images/model1_alpha_new} }
  \caption{Figures for the synthetic medium scale MDP. (a): Relative change of the policy
    $\|\pi_{\text{new}}-\pi\|_F / \|\pi\|_F$ during training of \Cref{alg:quasi} compared
    with PMD and GPMD in \cite{zhan2021policy}, with the logarithmic scale used for both axes.
    Notice that \Cref{alg:quasi} converges in 6-7 iterations to $10^{-12}$ in all cases
    while PMD and GPMD take more than $10^4$ iterations.  (b) - (e): Blue: The convergence of
    $\log\left|\log\|\pi-\pi^*\|_F\right|$ in the training process with the KL divergence, the
    reverse KL divergence, the Hellinger divergence and the $\alpha$-divergence with $\alpha=-3$,
    respectively. Green: A line through the origin with slope $\log2$. Comparison of the convergence plots with the green reference lines shows a clear quadratic convergence for \Cref{alg:quasi}.}
  \label{fig:toy1quad}
\end{figure}

\begin{sloppypar}
In order to verify the quadratic convergence proved in \cref{sec:analysis}, we draw the plots of
$\log\left|\log\|\pi-\pi^*\|\right|$ in \Cref{fig:model1_KL}, \Cref{fig:model1_rKL}, \Cref{fig:model1_Hellinger} and \Cref{fig:model1_alpha}, where $\pi^*$ is the final policy and the norm used is the Frobenius norm. A green reference line with slope $\log2$ through the origin is plotted for comparison. If the error converges exactly at a quadratic rate, the plot of $\log\left|\log\|\pi-\pi^*\|\right|$ shall be parallel to the reference line. The convergence curves approach the reference lines in the end (and are even steeper than the reference lines in the beginning), demonstrating clearly a quadratic convergence for all forms of regularization used here.
\end{sloppypar}

\subsection{Experiment II}

Next, we apply the quasi-Newton methods derived in \cref{sec:Newton} to an MDP model
constructed from the search logs of an online shopping store, with two different ranking strategies.
Each issued query is represented as a state in the MDP. In response to a query, the search can be
done by choosing one of the two ranking strategies (actions) to return a ranked list of products
shown to the customer. Based on the shown products, the customer can refine or update the query,
thus entering a new state. The reward at each state-action pair is a weighted sum of the clicks and
purchases resulting from the action. Based on the data collected from two separate 5-week periods
for both ranking strategies, we construct an MDP with 135k states, and a very sparse transition
tensor $P$ with only $0.01\%$ nonzero entries. The discount rate $\gamma$ is set as $0.99$ and the
regularization coefficient is $\tau=0.001$. In the implementation, we use the uniform prior $\mu_s^a=1/|A|$.

When calculating $v_\pi$ by $v_\pi = (I-\gamma P_\pi)^{-1}(r_\pi-\tau h_\pi)$, we apply the
iterative solver Bi-CGSTAB \cite{van1992bi}, a widely used numerical method with high efficiency
and robustness for solving large sparse nonsymmetric systems of linear equations
\cite{saad2003iterative,de1998comparison}, in order to leverage the sparsity of the transition
tensor.

In the numerical experiment, we implement \Cref{alg:quasi} with the KL divergence, the reverse KL divergence, the Hellinger divergence and the $\alpha$-divergence with $\alpha=-3$. We set the initial policy as the uniform policy, the convergence threshold as $\epsilon = 10^{-12}$ and the learning rate $\eta$ as $1$. All the tests end up with fast convergence as shown in \Cref{fig:real}, where logarithmic scale is used for the vertical axis. More specifically, the quasi-Newton algorithm using the KL divergence, the reverse KL divergence, the Hellinger divergence and the $\alpha$-divergence with $\alpha=-3$ converge in $6, 6, 6, 5$ iterations, respectively. It is worth noticing that even though the size of the state space $S$ here is some magnitudes larger than the examples in \cref{sec:exptoy}, the number of quasi-Newton iterations used is about the same. The comparison with GPMD and PMD is not given for this example since they are intractable to implement due to the high computational cost caused by the large MDP model.

In \Cref{table:BiCGSTAB}, we report the number of BiCGSTAB steps used in the algorithm. In each
quasi-Newton iteration, less than $20$ BiCGSTAB steps are used in order to find $v_\pi$. For all the $4$ regularizers used here, altogether only about $100$ BiCGSTAB steps are needed in the whole training process, thanks to the fast convergence of the quasi-Newton method. As a comparison, the regularized value iteration (a special case for the method in \cite{geist2019theory}) typically needs thousands of matrix-vector multiplication with the MDP transition matrix, since its convergence rate is $O(\gamma^T)$.

\begin{table}[ht]
\centering
\begin{tabular}{ cccc c  }
 \hline\hline
 Regularizer &\quad KL\quad&\quad reverse-KL\quad&\quad Hellinger\quad&$\alpha$-divergence ($\alpha=-3$)\\
 \hline\hline
 Quasi-Newton Iterations   & \quad $6$\quad   & \quad $6$\quad & \quad $6$\quad & \quad $5$\quad \\
 \hline
 Total Bi-CGSTAB steps     & \quad $110$\quad & \quad $109$\quad & \quad $110$\quad & \quad $83$\quad \\
 \hline
 Average Bi-CGSTAB steps   &\quad $18.3$\quad & \quad $18.2$\quad & \quad $18.3$\quad & \quad $16.6$\quad \\
 \hline
 \vspace{1pt}
\end{tabular}
  \caption{Number of quasi-Newton iterations and BiCGSTAB steps used in the training process. } \label{table:BiCGSTAB}
\end{table}

\begin{figure}[ht]
  \centering
  \subfigure[Relative change of the policy $\|\pi_{\text{new}}-\pi\|_F / \|\pi\|_F$ in the training process. \label{fig:real}]{
    \includegraphics[width=0.3\textwidth]{images/modelreal_conv_new.png} }
  \hspace{0.7em}
  \subfigure[The policy error in the process of training using KL-divergence. \label{fig:modelreal_KL}]{
    \includegraphics[width=0.3\textwidth]{images/modelreal_KL_new.png} }

  \subfigure[The policy error in the process of training using reverse KL-divergence. \label{fig:modelreal_rKL}]{
    \includegraphics[width=0.3\textwidth]{images/modelreal_rKL_new.png} }
  \hspace{0.7em}
  \subfigure[The policy error in the process of training using Hellinger divergence. \label{fig:modelreal_Hellinger}]{
    \includegraphics[width=0.3\textwidth]{images/modelreal_Hellinger_new.png} }
  \hspace{0.7em}
  \subfigure[The policy error in the process of training using $\alpha$-divergence with $\alpha=-3$. \label{fig:modelreal_alpha}]{
    \includegraphics[width=0.3\textwidth]{images/modelreal_alpha_new} }
  \caption{Figures for the industrial-size MDP. (a): Relative change of the policy $\|\pi_{\text{new}}-\pi\|_F / \|\pi\|_F$ in the training process of \Cref{alg:quasi}. Logarithmic scale is used for the vertical axis. (b) - (e): Blue: The convergence of $\log\left|\log\|\pi-\pi^*\|_F\right|$ in the training process with the KL divergence, the reverse KL divergence, the Hellinger divergence and the $\alpha$-divergence with $\alpha=-3$, respectively. Green: A line through the origin with slope $\log2$. }
  \label{fig:realquad}
\end{figure}

As in the previous numerical example, in \Cref{fig:modelreal_KL}, \Cref{fig:modelreal_rKL}, \Cref{fig:modelreal_Hellinger} and \Cref{fig:modelreal_alpha} we verify the quadratic convergence by comparing the plot of $\log\left|\log\|\pi-\pi^*\|\right|$ with a green reference line through the origin with slope $\log2$. As the convergence curves are approximately parallel to the reference lines, this verifies that the proposed algorithm converges quadratically with all the regularizations in this example as well.

\section{Proofs}\label{sec:proofs}

\subsection{Proof of \Cref{thm: approxnewton}}
\label{sec:approxnewton}
\begin{proof}
  \textbf{Step 1: expand $E(\pi)$ and prove the first-order condition \eqref{eq:grad}.}
  Since $\pi$ is a policy, $\sum_a\pi_s^a = 1$ for any $s$. Thus
  \begin{equation}\label{eq:linear0}
    (Z_\pi)_{st} = \delta_{st} - \gamma\sum_a\pi_s^aP_{st}^a = \sum_a\pi_s^a(\delta_{st}-\gamma P_{st}^a).
  \end{equation}
  Now consider a policy $\pi+\epsilon$, i.e., $\sum_a\epsilon_s^a = 0$ and
  $\pi_s^a+\epsilon_s^a\geq0$, then thanks to \eqref{eq:linear0} one can obtain:
  \begin{equation}\label{eq:linear}
    Z_{\pi+\epsilon} = Z_\pi+Z_\epsilon, \quad r_{\pi+\epsilon} = r_\pi+r_\epsilon,
  \end{equation}
  where $Z_\epsilon$ and $r_\epsilon$ are defined in \eqref{eq:rZexpansion}, i.e.,
  $(Z_\epsilon)_{st}=\sum_a\epsilon_s^a(\delta_{st}-\gamma P_{st}^a)$. $(r_\epsilon)_s =\sum_a\epsilon_s^ar_{s}^a$.
  Leveraging the linearity \eqref{eq:linear}, we obtain the expansion:
  \begin{equation}\label{eq:expansion}
    \begin{aligned}
      &\quad E(\pi+\epsilon) = e^\top Z_{\pi+\epsilon}^{-1}(r_{\pi+\epsilon}-\tau
      h_{\pi+\epsilon}) =e^\top (Z_{\pi}+Z_{\epsilon})^{-1}(r_{\pi}+r_{\epsilon}-\tau
      h_{\pi+\epsilon})\\ &=e^\top Z_{\pi}^{-1}
      (I-Z_{\epsilon}Z_\pi^{-1}+Z_{\epsilon}Z_\pi^{-1}Z_{\epsilon}Z_\pi^{-1})(r_{\pi}+r_{\epsilon}-\tau
      h_{\pi}-\tau Dh_\pi\epsilon-\frac{1}{2}\epsilon^\top \tau D^2h_\pi\epsilon) + O(\|\epsilon\|^3)\\
      &= E(\pi) + w_\pi^\top\left[-Z_{\epsilon}Z_\pi^{-1}(r_\pi-\tau h_\pi)+(r_\epsilon-\tau Dh_\pi\epsilon)\right]+w_\pi^\top(-\frac{1}{2}\epsilon^\top \tau
        D^2h_\pi\epsilon)\\
      &\quad+w_\pi^\top\left[-Z_{\epsilon}Z_\pi^{-1}(r_\epsilon-\tau
        Dh_\pi\epsilon)+Z_{\epsilon}Z_\pi^{-1}Z_{\epsilon}Z_\pi^{-1}(r_\pi-\tau
        h_\pi)\right] + O(\|\epsilon\|^3),
    \end{aligned}
  \end{equation}
  where $Dh_\pi$ is a second-order tensor that maps from $S\times A$ to $S$, and $D^2h_\pi$ is a
  third-order tensor that maps from $(S\times A)^{\otimes2}$ to $S$. With this expansion, one can see
  that
  \[
    \frac{\partial E}{\partial \pi_s^a} = (r_s^a - \tau (\log\pi_s^a +1) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s)(w_\pi)_s,
  \]
  where $c_s$ is a multiplier that only depends on $s$.
  Then at $\pi = \pi^*$,
  \[
  \frac{\partial E}{\partial \pi_s^a} = (r_s^a - \tau (\log\pi_s^a + 1) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s)(w_\pi)_s = 0.
  \]
  Since $w_\pi = (I-\gamma P_\pi^\top)^{-1}e = e + \sum_{i=1}^\infty \gamma^i(P_\pi^\top)^ie$ and
  all elements of $e$ are positive, we also know that all elements of $w_\pi$ are positive. Thus at
  $\pi = \pi^*$,
  \[
  r_s^a - \tau (\log\pi_s^a + 1) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s = 0.
  \]
  Multiplying the left hand side with $\epsilon_s^a$ and taking the sum over $a$, we obtain:
  \[
  \left(r_\epsilon-\tau Dh_\pi\epsilon-Z_\epsilon v_\pi\right)_s + c_s\sum_a\epsilon_s^a = 0, \quad \forall s, \quad\forall \epsilon.
  \]
  Since $\sum_a\epsilon_s^a = 0$ for any $s$ and $v_\pi = Z_\pi^{-1}(r_\pi-\tau h_\pi)$, we have
  \[
  r_\epsilon-\tau Dh_\pi\epsilon-Z_\epsilon Z_\pi^{-1}(r_\pi-\tau h_\pi) = 0, \quad\forall \epsilon,
  \]
  at $\pi = \pi^*$, which proves \eqref{eq:grad}.

  \textbf{Step 2: Derive the decomposition \eqref{eq:decomp} with the obtained expansion and first-order condition.}
  With \eqref{eq:grad}, one can approximate the second-order term in \eqref{eq:expansion} for $\pi$ near $\pi^*$:
  \[
    \begin{aligned}
      &w_\pi^\top(-\frac{1}{2}\epsilon^\top \tau D^2h_\pi\epsilon)+w_\pi^\top\left[-Z_{\epsilon}Z_\pi^{-1}(r_\epsilon-\tau Dh_\pi\epsilon)+Z_{\epsilon}Z_\pi^{-1}Z_{\epsilon}Z_\pi^{-1}(r_\pi-\tau h_\pi)\right]\\
      &=\frac{1}{2}\epsilon^\top\Lambda(\pi)\epsilon-w_\pi^\top Z_{\epsilon}Z_\pi^{-1}\left(r_\epsilon-\tau Dh_\pi\epsilon-Z_{\epsilon}Z_\pi^{-1}(r_\pi-\tau h_\pi)\right)\\
      &\approx \frac{1}{2}\epsilon^\top\Lambda(\pi)\epsilon.
    \end{aligned}
  \]
  By \eqref{eq:grad} and the twice continuous differentiability of $h$, the approximate Hessian $\Lambda$ converge to the true Hessian as $\pi$ converges to $\pi^*$, and their difference $\Lambda(\pi) - D^2E(\pi) = O(\|\pi-\pi^*\|)$. Hence, the second-order derivatives of $E(\pi)$ can be approximated by
  \begin{equation}
    \frac{\partial^2 E}{\partial \pi_s^a\partial \pi_t^b}\approx \Lambda_{(sa),(tb)} = -\tau \delta_{\{(sa),(tb)\}}\frac{(w_\pi)_s}{\pi_s^a},
  \end{equation}
  from which we have shown that $\Lambda$ is diagonal.

  \textbf{Step 3: Derive the quasi-Newton flow and the policy update scheme with the obtained decomposition.}
  Using this approximate second-order derivative as a preconditioner, $w_\pi$ is canceled out in the
  policy gradient algorithm, which yields the gradient flow:
  \[
  \frac{\mathrm{d}\pi_s^a}{\mathrm{d} t} = \pi_s^a(r_s^a-\tau(\log\pi_s^a+1)-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau.
  \]
  Adopting the parameterization $\pi_s^a = \exp(\theta_s^a)$, we have
  \begin{equation}\label{eq:paramgrad}
    \frac{\mathrm{d}\theta_s^a}{\mathrm{d} t} = (r_s^a-\tau(\theta_s^a+1)-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau.
  \end{equation}
  With a learning rate $\eta$, this becomes
  \begin{equation}
    \theta_s^a \leftarrow \eta(r_s^a-\tau-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau + (1-\eta)\theta_s^a,
  \end{equation}
  which corresponds to
  \begin{equation}
    \pi_s^a \leftarrow (\pi_s^a)^{1-\eta}\exp(\eta(r_s^a-\tau-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau),
  \end{equation}
  and $c_s$ is determined by the condition that $\sum_a\pi_s^a=1$. Equivalently, we have
  \[
  \pi_s^a \leftarrow \frac{(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a-\tau-[(I-\gamma P^a)v_\pi]_s)/\tau)}{\sum_a (\pi_s^a)^{1-\eta}\exp(\eta(r_s^a-\tau-[(I-\gamma P^a)v_\pi]_s)/\tau)} = \frac{(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma P^av_\pi)_s)/\tau)}{\sum_a(\pi_s^a)^{1-\eta}\exp(\eta(r_s^a+(\gamma P^av_\pi)_s)/\tau)},
  \]
  where we cancel out the factors independent of $a$ and obtain \eqref{eq:NPG}. This finishes the proof.
\end{proof}

\subsection{Proof of \Cref{thm: otheretp}}
\label{sec:otheretp}
\begin{proof}

  Similar with the proof of \Cref{thm: approxnewton}, by direct calculations one can get:
  \begin{equation}
    \frac{\partial E}{\partial \pi_s^a} = (r_s^a - \tau \phi'(\pi_s^a/\mu_s^a) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s)(w_\pi)_s,
  \end{equation}
  where $c_s$ is a multiplier that only depends on $s$. Since all elements of $w_\pi$ are positive, at $\pi = \pi^*$,
  \begin{equation}\label{eq:general1st}
    (r_s^a - \tau \phi'(\pi_s^a/\mu_s^a) - \left[(I-\gamma P^a)v_\pi\right]_s + c_s) = 0.
  \end{equation}
  By multiplying \eqref{eq:general1st} with $\epsilon_s^a$ and summing over $a$, one can obtain:
  \[
  r_\epsilon-\tau Dh_\pi\epsilon-Z_\epsilon Z_\pi^{-1}(r_\pi-\tau h_\pi) = 0, \quad\forall \epsilon,
  \]
  at $\pi = \pi^*$, which proves \eqref{eq:general1stmain}.
  Since the only difference between the functional $E(\pi)$ defined here and the $E(\pi)$ in \Cref{thm: approxnewton} lies in the regularizer $h$, one can still obtain the expansion:
  \[
    \begin{aligned}
    &E(\pi+\epsilon) - E(\pi) - w_\pi^\top\left[-Z_{\epsilon}Z_\pi^{-1}(r_\pi-\tau h_\pi)+(r_\epsilon-\tau Dh_\pi\epsilon)\right]\\
    &= w_\pi^\top(-\frac{1}{2}\epsilon^\top \tau D^2h_\pi\epsilon) -w_\pi^\top Z_{\epsilon}Z_\pi^{-1}\left(r_\epsilon-\tau Dh_\pi\epsilon-Z_{\epsilon}Z_\pi^{-1}(r_\pi-\tau h_\pi)\right) + O(\|\epsilon\|^3)\\
    &= \frac{1}{2}\epsilon^\top\Lambda(\pi)\epsilon + O(\|\epsilon\|^2\|\pi-\pi^*\|) + O(\|\epsilon\|^3).
    \end{aligned}
  \]
  Hence we have $D^2E(\pi)-\Lambda(\pi) = O(\|\pi-\pi^*\|)$. Using this expansion, one can derive an approximation for the second-order derivatives:
  \[
    \frac{\partial^2 E}{\partial \pi_s^a\partial \pi_t^b} \approx \Lambda_{(sa),(tb)} = -\tau\delta_{\{(sa),(tb)\}}\frac{(w_\pi)_s\phi''(\pi_s^a/\mu_s^a)}{\mu_s^a},
  \]
  which proves \eqref{eq:generalapproxhess} and shows that $\Lambda$ is diagonal. The quasi-Newton flow thus becomes
  \[
    \frac{\mathrm{d}\pi_s^a}{\mathrm{d} t} = \mu_s^a(\phi''(\pi_s^a/\mu_s^a))^{-1}(r_s^a- \tau \phi'(\pi_s^a/\mu_s^a)-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau,
  \]
  which proves \eqref{eq:qNflowgen}, or equivalently,
  \eq{
    \frac{\mathrm{d}(\phi'(\pi_s^a/\mu_s^a))}{\mathrm{d} t} = (r_s^a- \tau \phi'(\pi_s^a/\mu_s^a)-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau.
  }
  Let $\theta_s^a = \phi'(\pi_s^a/\mu_s^a)$, then
  \[
  \frac{\mathrm{d}\theta_s^a}{\mathrm{d} t} = (r_s^a-\tau\theta_s^a-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau,
  \]
  With a learning rate $\eta$, this becomes
  \[
    \theta_s^a \leftarrow \eta(r_s^a-[(I-\gamma P^a)v_\pi]_s + c_s)/\tau + (1-\eta)\theta_s^a,
  \]
  which proves \eqref{eq:quasiNTupdate}.

\end{proof}

\subsection{Proof of \Cref{lemma:multipliergen}}
\label{sec:mult}

\begin{proof}
  Let
  \[
  g(x) = \mu_1\psi(x+x_1) + \cdots + \mu_k\psi(x+x_k).
  \]
  Since $\psi:(L, +\infty)\rightarrow(0, +\infty)$ is decreasing, $g(x)$ is positive and decreasing on $(L-\underset{1\leq i\leq k}{\min}x_i, +\infty)$. When $x\rightarrow-\underset{1\leq i\leq k}{\min}x_i$
  from the right, $g(x)\rightarrow+\infty$ since at least one of the summands go to $+\infty$. If
  $\underset{1\leq i\leq k}{\min}\left\{\psi^{-1}(\frac{1}{k\mu_i}) - x_i\right\}\geq L-\underset{1\leq i\leq k}{\min}x_i$, when $x
  = \underset{1\leq i\leq k}{\min}\left\{\psi^{-1}(\frac{1}{k\mu_i}) - x_i\right\}$
  \[
  \begin{aligned}
  g(x) &= \sum_{i=1}^k\mu_i\psi\left(\underset{1\leq j\leq k}{\min}\left\{\psi^{-1}\left(\frac{1}{k\mu_j}\right) - x_j\right\}+x_i\right)\\
  &= \sum_{i=1}^k\mu_i\psi\left(\underset{1\leq j\leq k}{\min}\left\{\psi^{-1}\left(\frac{1}{k\mu_j}\right) - x_j\right\}+x_i-\psi^{-1}\left(\frac{1}{k\mu_i}\right)+\psi^{-1}\left(\frac{1}{k\mu_i}\right)\right)\\
  &\geq \sum_{i=1}^k\mu_i\psi\left(\psi^{-1}\left(\frac{1}{k\mu_i}\right)\right) =\sum_{i=1}^k \mu_i\times\frac{1}{k\mu_i} = k \times \frac{1}{k} = 1.
  \end{aligned}
  \]
  Since $\psi^{-1}\left(\frac{1}{k\mu_i}\right)\geq L$, we have $\underset{1\leq i\leq k}{\max}\left\{\psi^{-1}\left(\frac{1}{k\mu_i}\right) - x_i\right\}\geq \underset{1\leq i\leq k}{\max}\left\{L - x_i\right\} = L - \underset{1\leq i\leq k}{\min}x_i$.
  Then when $x = \underset{1\leq i\leq k}{\max}\left\{\psi^{-1}\left(\frac{1}{k\mu_i}\right) - x_i\right\}$,
  \[
  g(x) = \sum_{i=1}^k\mu_i\psi\left(\underset{1\leq j\leq k}{\max}\left\{\psi^{-1}\left(\frac{1}{k\mu_j}\right) - x_j\right\}+x_i\right)\leq \sum_{i=1}^k\mu_i\psi\left(\psi^{-1}\left(\frac{1}{k\mu_i}\right)\right) = 1.
  \]
  By the continuity of $g$, there exists a solution $x$ to \eqref{eq:lemmagen} on
  \[
  \left[\max\left\{L-\underset{1\leq i\leq k}{\min}x_i, \underset{1\leq i\leq k}{\min}\left\{\psi^{-1}\left(\frac{1}{k\mu_i}\right) - x_i\right\}\right\}, \underset{1\leq i\leq k}{\max}\left\{\psi^{-1}\left(\frac{1}{k\mu_i}\right) - x_i\right\}\right],
  \]
  and the solution is unique by the strict monotonicity of $g$ on $(L-\underset{1\leq i\leq k}{\min}x_i,
  \infty)$.
\end{proof}

\subsection{Proof of \Cref{prop:multiplier}}
\label{sec:multsolve}
\begin{proof}
By \Cref{lemma:multipliergen} there is a unique solution $\tilde{c}_s$ to the equation $\sum_a \mu_s^a \psi(\tilde{c}_s + x_a) = 1$, where $x_a$ is defined as in \eqref{eq:deftildecxa}.
Now update the policy by
\[
\pi_s^a\leftarrow \mu_s^a \psi(\tilde{c}_s + x_a) = \sum_a\mu_s^a\psi\left(-\frac{\eta}{\tau}c_s-(1-\eta)\phi'\left(\pi_s^a/\mu_s^a\right) -\frac{\eta}{\tau}(r_s^a-[(I-\gamma P^a)v_\pi]_s)\right)
\]
one ensures that $\pi_s^a\geq0$ and $\sum_a\pi_s^a=1$, and that the multiplier $c_s$ with this property is unique.

\end{proof}

\subsection{Proof of \Cref{thm:flowconv}}
\label{sec:flowconvpf}
\begin{proof}
In \cref{sec:flowconv} we have proved that the quasi-Newton flow:
\[
\frac{\mathrm{d}\pi_s^a}{\mathrm{d} t} = \mu_s^a(\phi''(\pi_s^a/\mu_s^a))^{-1}(r_s^a- \tau \phi'(\pi_s^a/\mu_s^a)-[(I-\gamma P^a)v_\pi]_s+c_s)/\tau
\]
converges globally, so it suffices to show that the limiting policy is the optimal policy. Denote the limiting policy by $\pi^\diamond$. Since $\mu_s^a>0$ and $(\phi''((\pi^\diamond)_s^a/\mu_s^a))^{-1}>0$, we have
\begin{equation}\label{eq:flow1st}
    r_s^a- \tau \phi'((\pi^\diamond)_s^a/\mu_s^a)-[(I-\gamma P^a)v_{\pi^\diamond}]_s+c_s = 0,
\end{equation}
and $c_s$ is a multiplier which ensures $\sum_a(\pi^\diamond)_s^a = 1$. From the theory of regularized MDP (Cf. \cite{geist2019theory}), we know that the optimal policy $\pi^*$ is the unique solution to the Bellman maximal equation:
\begin{equation}\label{eq:regBell}
    v = \underset{\pi}{\max}~ r_\pi + \gamma P_\pi v - \tau h_\pi.
\end{equation}
Since $v_{\pi^\diamond} = (I-\gamma P_{\pi^\diamond})^{-1}(r_{\pi^\diamond}-\tau h_{\pi^\diamond})$, we have $v_{\pi^\diamond}-\gamma P_{\pi^\diamond}v_{\pi^\diamond} = r_{\pi^\diamond}-\tau h_{\pi^\diamond}$, or equivalently
\[v_{\pi^\diamond} = r_{\pi^\diamond} + \gamma P_{\pi^\diamond}v_{\pi^\diamond}-\tau h_{\pi^\diamond}.\]
Thus it now suffices to show that $\pi^\diamond$ is the optimizer of the constrained maximization problem $\underset{\pi}{\max}~ r_\pi + \gamma P_\pi v_{\pi^\diamond} - \tau h_\pi$, or in the component form:
\begin{equation}\label{eq:optprob}
  \underset{\pi}{\max} \sum_a (r_s^a+\gamma (P^av_{\pi^\diamond})_s)\pi_s^a - \tau \sum_a\mu_s^a\phi(\pi_s^a/\mu_s^a).
\end{equation}
Since $\phi$ is convex and $\mu$ is positive, $\tau \sum_a\mu_s^a\phi(\pi_s^a/\mu_s^a)$ is also a convex function in $\pi_s$. By the theory of convex optimization (Cf. \cite{boyd2004convex}, chapter $5$), the Karush-Kuhn-Tucker (KKT) condition is sufficient for the optimality when the objective function is convex, and the KKT condition for the problem \eqref{eq:optprob} is
\[
\begin{aligned}
r_s^a+\gamma (P^av_{\pi^\diamond})_s -\tau \phi'(\pi_s^a/\mu_s^a) + \lambda_s &=0,\\
\sum_a\pi_s^a&=1,\\
\pi_s^a&\geq0,
\end{aligned}
\]
where $\lambda_s$ is the Lagrange multiplier. Now let $\pi=\pi^{\diamond}$ and $\lambda_s = c_s-(v_{\pi^\diamond})_s$. From the first-order condition \eqref{eq:flow1st} one can directly observe that the KKT condition above is satisfied, which make $\pi^\diamond$ the optimizer for \eqref{eq:optprob} and $v_{\pi^\diamond}$ the solution to the Bellman equation \eqref{eq:regBell}. Thus $v_{\pi^\diamond}$ and $\pi^\diamond$ are indeed the optimal value function and the optimal policy, which closes the proof.
\end{proof}

\subsection{Proof of \Cref{thm:quadconv}}
\label{sec:quadconvpf}
\begin{proof}

Since $B_{ij} = 1$ for $|A|(i-1)+1\leq j\leq|A|i$ and $B_{ij}=0$ otherwise, the constraint on $\pi$ can be expressed by $B\pi = \mathbf{1}_{|S|}$. The gradient flow \eqref{eq:qNflowgen} can thus be expressed by
\eq{
  \frac{\mathrm{d}(\nabla\Phi(\pi))}{\mathrm{d} t} = -(f(\pi) + \nabla\Phi(\pi) - B^\top c(\pi)),
}
where $c(\pi)$ is the multiplier. The discrete scheme \eqref{eq:quasiNTupdate} corresponds to:
\eq{\label{eq:quasiupdate}
  \nabla\Phi(\pi\sps{k+1})-\nabla\Phi(\pi\sps{k}) = -\left(f(\pi\sps{k}) + \nabla\Phi(\pi\sps{k}) -
  B^\top c(\pi\sps{k})\right),
}
when setting the learning rate $\eta$ as $1$.

With this update scheme, we have
\begin{equation}\label{eq:replace}
    \begin{aligned}
        &f(\pi\sps{k+1}) - f(\pi\sps{k}) - \nabla f(\pi^*)\left(\pi\sps{k+1}-\pi\sps{k}\right) \\
  =&f(\pi\sps{k+1})+\nabla\Phi(\pi\sps{k+1})-B^\top c(\pi\sps{k})-\nabla f(\pi^*)\left(\pi\sps{k+1}-\pi\sps{k}\right)
    \end{aligned}
\end{equation}
Since $\theta\sps{k} = \nabla\Phi(\pi\sps{k})$ converges to $\theta^*$ and the map from $\theta$ to
$\pi$ is continuous, $\pi$ converges to some $\pi^*$ as well. Since
$\pi\sps{k}\in\{\pi:B^\top\pi=\mathbf{1}_{|S|}, \pi_s^a\geq0\}$ and
$\{\pi:B^\top\pi=\mathbf{1}_{|S|}, \pi_s^a\geq0\}$ is a closed set, we also have
$\pi^*\in\{\pi:B^\top\pi=\mathbf{1}_{|S|}, \pi_s^a\geq0\}$. Moreover, since $\nabla\Phi(\pi)$ goes to
infinity as any entry $\pi_s^a\rightarrow0$, $(\pi^*)_s^a>0$ for any $s$ and $a$. Hence
$\pi\sps{k}$ is contained in a closed set $K$ contained in $\{\pi:B^\top\pi=\mathbf{1}_{|S|},
\pi_s^a>0\}$ for large enough $k$.
Since $\nabla f$ is Lipschitz continuous on $K$,
\eq{\label{eq:Lip}
    \underset{k\rightarrow\infty}{\lim}\frac{f(\pi\sps{k+1}) - f(\pi\sps{k}) - \nabla f(\pi^*)\left(\pi\sps{k+1}-\pi\sps{k}\right)}{\norm{\pi\sps{k+1}-\pi\sps{k}}}=0.
}
Moreover, since $f(\pi)_{sa} = -(r_s^a-((I-\gamma P^a)v_\pi)_s)$ has a similar form with $E(\pi)$, we can directly obtain $\nabla f(\pi)$
\eq{
  (\nabla f(\pi))_{sa, tb} = \lambda_{sa, t}(\pi) \left(-f(\pi)_{tb}+\tilde{c}(\pi)_t-\nabla\Phi(\pi)_{tb}\right),
}
where $\lambda_{sa,t}(\pi) = Z_\pi^{-\top}\tilde{e}_{sa}$ and $\tilde{e}_{sa}$ is the $s$-th row of $I-\gamma P^a$. Notice that at $\pi^*$ we have $\frac{\mathrm{d}(\nabla\Phi(\pi))}{\mathrm{d} t}=0$, so $f(\pi^*)+\nabla\Phi(\pi^*) = B^\top c(\pi^*)$, thus
\begin{equation}
\label{eq:KKT}
\begin{aligned}
    &(\nabla f(\pi^*)(\pi\sps{k+1}-\pi\sps{k}))_{sa}\\
    =& \sum_{tb} \lambda_{sa, t}(\pi^*) \left(-f(\pi^*)_{tb}+\tilde{c}(\pi^*)_t-\nabla\Phi(\pi^*)_{tb}\right)(\pi_{tb}\sps{k+1}-\pi_{tb}\sps{k})\\
    =& \sum_{tb} \lambda_{sa, t}(\pi^*) \left(\tilde{c}(\pi^*)_t-c(\pi^*)_t\right)(\pi_{tb}\sps{k+1}-\pi_{tb}\sps{k})\\
    =& \sum_{t} \left[\Bigg(\lambda_{sa, t}(\pi^*) \left(\tilde{c}(\pi^*)_t-c(\pi^*)_t\right)\Bigg)\left(\sum_b(\pi_{tb}\sps{k+1}-\pi_{tb}\sps{k})\right)\right]\\
    =&~0,
\end{aligned}
\end{equation}
where the last equality results from the fact that $\sum_b\pi_{tb}\sps{k+1}=\sum_b\pi_{tb}\sps{k}=1$ for any $t$.
With \eqref{eq:Lip} and \eqref{eq:KKT}, we arrive at
\eq{\label{eq:supconvergepre}
    \underset{k\rightarrow\infty}{\lim}\frac{f(\pi\sps{k+1})+\nabla\Phi(\pi\sps{k+1})-B^\top c(\pi\sps{k})}{\norm{\pi\sps{k+1}-\pi\sps{k}}}=0.
}
Multiplying the fraction above by the unit vector $\frac{\pi\sps{k+2}-\pi\sps{k+1}}{\norm{\pi\sps{k+2}-\pi\sps{k+1}}}$, we obtain
\begin{equation}\label{eq:transformation}
    \begin{aligned}
    &\frac{\left(f(\pi\sps{k+1})+\nabla\Phi(\pi\sps{k+1})-B^\top c(\pi\sps{k})\right)^\top(\pi\sps{k+2}-\pi\sps{k+1})}{\norm{\pi\sps{k+1}-\pi\sps{k}}\norm{\pi\sps{k+2}-\pi\sps{k+1}}}\\
    &=\frac{\left(f(\pi\sps{k+1})+\nabla\Phi(\pi\sps{k+1})-B^\top c(\pi\sps{k+1})\right)^\top(\pi\sps{k+2}-\pi\sps{k+1})}{\norm{\pi\sps{k+1}-\pi\sps{k}}\norm{\pi\sps{k+2}-\pi\sps{k+1}}}\\
    &=\frac{\left(\nabla\Phi(\pi\sps{k+2})-\nabla\Phi(\pi\sps{k+1})\right)^\top(\pi\sps{k+2}-\pi\sps{k+1})}{\norm{\pi\sps{k+1}-\pi\sps{k}}\norm{\pi\sps{k+2}-\pi\sps{k+1}}},
    \end{aligned}
\end{equation}
where we use the constraint $B\pi\sps{k+1} = B\pi\sps{k+2} = \mathbf{1}_{|S|}$. Then
from \eqref{eq:supconvergepre},
\eq{
    \underset{k\rightarrow\infty}{\lim}\frac{\left(\nabla\Phi(\pi\sps{k+2})-\nabla\Phi(\pi\sps{k+1})\right)^\top(\pi\sps{k+2}-\pi\sps{k+1})}{\norm{\pi\sps{k+1}-\pi\sps{k}}\norm{\pi\sps{k+2}-\pi\sps{k+1}}}=0.
}
By direct calculation of $\nabla^2\Phi$, which is diagonal for all the functions $\Phi$ defined in this theorem, we can see that $\Phi$ is strongly convex since $\phi$ is strongly convex. As a result,
there is some constant $\rho>0$ such that
\[
(\nabla\Phi(\pi)-\nabla\Phi(\tilde{\pi}))^\top(\pi-\tilde{\pi}) \geq \rho\norm{\pi-\tilde{\pi}}^2
\]
for any $\pi$ and $\tilde{\pi}$ in $K$. Then
\eq{
    0=\underset{k\rightarrow\infty}{\lim}\frac{\norm{\pi\sps{k+2}-\pi\sps{k+1}}^2}{\norm{\pi\sps{k+1}-\pi\sps{k}}\norm{\pi\sps{k+2}-\pi\sps{k+1}}}= \underset{k\rightarrow\infty}{\lim}\frac{\norm{\pi\sps{k+2}-\pi\sps{k+1}}}{\norm{\pi\sps{k+1}-\pi\sps{k}}},
}
from which we can conclude that $\pi\sps{k}$ converges to $\pi^*$ superlinearly, i.e., \eqref{eq:suplin} holds.
In fact, for any $\epsilon$ (assume $\epsilon<1/2$ without loss of generality), there is some $k(\epsilon)$ such that for any $k>k(\epsilon)$,
\[
\frac{\norm{\pi\sps{k+2}-\pi\sps{k+1}}}{\norm{\pi\sps{k+1}-\pi\sps{k}}} < \epsilon,
\]
then for any $k>k(\epsilon)$
\eqs{
    \norm{\pi\sps{k+1}-\pi^*}&\leq \sum_{n=k+1}^\infty\norm{\pi\sps{n+1}-\pi\sps{n}}\leq \sum_{n=k+1}^\infty\epsilon^{n-k}\norm{\pi\sps{k+1}-\pi\sps{k}}\\
    &\leq \frac{\epsilon}{1-\epsilon}\norm{\pi\sps{k+1}-\pi\sps{k}}\leq 2\epsilon \norm{\pi\sps{k+1}-\pi\sps{k}}.
}
Then
\eqs{
    \norm{\pi\sps{k}-\pi^*} &\geq\norm{\pi\sps{k+1}-\pi\sps{k}}-\norm{\pi\sps{k+1}-\pi^*}\\
    &\geq (\frac{1}{2\epsilon}-1)\norm{\pi\sps{k+1}-\pi^*}.
}
For any $M>0$, take $\epsilon = 1/(2M+2)$, then for any $k>k(\epsilon)$,
\eq{
    \norm{\pi\sps{k}-\pi^*} \geq (\frac{1}{2\epsilon}-1)\norm{\pi\sps{k+1}-\pi^*} = (M+1)\norm{\pi\sps{k+1}-\pi^*}>M\norm{\pi\sps{k+1}-\pi^*},
}
which shows that $\underset{k\rightarrow\infty}{\lim}\frac{\|\pi\sps{k}-\pi^*\|}{\|\pi\sps{k+1}-\pi^*\|} = +\infty$ and thus
\begin{equation}\label{eq:suplin}
\underset{k\rightarrow\infty}{\lim}\frac{\|\pi\sps{k+1}-\pi^*\|}{\|\pi\sps{k}-\pi^*\|} = 0.
\end{equation}
Now, from \eqref{eq:replace} and \eqref{eq:KKT} we have
\eqs{
&f(\pi\sps{k+1})+\nabla\Phi(\pi\sps{k+1})-B^\top c(\pi\sps{k})\\
=&~f(\pi\sps{k+1}) - f(\pi\sps{k}) - \nabla f(\pi^*)\left(\pi\sps{k+1}-\pi\sps{k}\right)\\
=&\left(\int_0^1\left[\nabla f(\pi\sps{k}+t(\pi\sps{k+1}-\pi\sps{k}))-\nabla f(\pi^*)\right]\mathrm{d} t\right) \left(\pi\sps{k+1}-\pi\sps{k}\right)\\
\leq&~\tilde{C}\norm{\pi\sps{k}-\pi^*}\norm{\pi\sps{k+1}-\pi\sps{k}},
}
for some constant $\tilde{C}$, where we used \eqref{eq:suplin} and the Lipschitz continuity of $\nabla f$ in the last equality. Multiplying both sides by $\frac{\pi\sps{k+2}-\pi\sps{k+1}}{\norm{\pi\sps{k+2}-\pi\sps{k+1}}}$ and by \eqref{eq:transformation} and the strong convexity of $\Phi$, we have
\eqs{
    &\rho\norm{\pi\sps{k+2}-\pi\sps{k+1}}^2\\
    \leq&\left(\nabla\Phi(\pi\sps{k+2})-\nabla\Phi(\pi\sps{k+1})\right)^\top(\pi\sps{k+2}-\pi\sps{k+1})\\
    =&\left(f(\pi\sps{k+1})+\nabla\Phi(\pi\sps{k+1})-B^\top c(\pi\sps{k})\right)^\top(\pi\sps{k+2}-\pi\sps{k+1})\\
    \leq&~\tilde{C}\norm{\pi\sps{k}-\pi^*}\norm{\pi\sps{k+1}-\pi\sps{k}}\norm{\pi\sps{k+2}-\pi\sps{k+1}},
}
which implies that
\eq{\label{eq:quadpre}\norm{\pi\sps{k+2}-\pi\sps{k+1}} \leq \tilde{C}\norm{\pi\sps{k}-\pi^*}\norm{\pi\sps{k+1}-\pi\sps{k}}.}
From \eqref{eq:suplin}, we have
\eq{
    \underset{k\rightarrow\infty}{\lim}\frac{\|\pi\sps{k}-\pi\sps{k+1}\|}{\|\pi\sps{k}-\pi^*\|} = \underset{k\rightarrow\infty}{\lim}\frac{\|\pi\sps{k+1}-\pi\sps{k+2}\|}{\|\pi\sps{k+1}-\pi^*\|} = 1.
}
Combining this with \eqref{eq:quadpre} leads to
\eq{
    \norm{\pi\sps{k+1}-\pi^*} \leq C\norm{\pi\sps{k}-\pi^*}^2,
}
for some constant $C$, which closes the proof.

\end{proof}

\section{Conclusion and Discussion}\label{sec:conclusion}

In this paper, we present a fast quasi-Newton method for the policy gradient algorithm with provable quadratic convergence. The proposed method gives a systematic theory that includes the well-known natural policy gradient algorithm as a special case, and naturally extends to other regularizers such as the reverse KL divergence, the Hellinger divergence and the $\alpha$-divergence.

With a relatively simple proof, we show local quadratic convergence of the proposed quasi-Newton method as well as the global convergence of the quasi-Newton gradient flow to the optimal solution. The quadratic convergence is confirmed numerically on both medium and
large sparse models. In contrast with mirror descent type first-order methods (e.g. \cite{zhan2021policy}) that take up to tens of thousands iterations even with manually tuned learning rate, the proposed quasi-Newton algorithms typically converge in less than $10$ iterations, despite the large discount rate ($\approx 1$) and small regularization coefficient ($\approx 0$).

For future work, we plan to adapt the technique used here to other gradient based algorithms for
solving the MDP problems. Other forms of $f$-divergence can also be included. Another direction is
to consider continuous MDP problems by leveraging function approximation, effective spatial
discretization, or model reduction.

\bibliographystyle{siamplain}

\section{A detailed example}

Here we include some equations and theorem-like environments to show
how these are labeled in a supplement and can be referenced from the
main text.
Consider the following equation:
\begin{equation}
  \label{eq:suppa}
  a^2 + b^2 = c^2.
\end{equation}
You can also reference equations such as \cref{eq:matrices,eq:bb}
from the main article in this supplement.

\lipsum[100-101]

\begin{theorem}
  An example theorem.
\end{theorem}

\lipsum[102]

\begin{lemma}
  An example lemma.
\end{lemma}

\lipsum[103-105]

Here is an example citation: \cite{KoMa14}.

\section[Proof of Thm]{Proof of \cref{thm:bigthm}}
\label{sec:proof}

\lipsum[106-112]

\section{Additional experimental results}
\Cref{tab:foo} shows additional
supporting evidence.

\begin{table}[htbp]
{\footnotesize
  \caption{Example table}  \label{tab:foo}
\begin{center}
  \begin{tabular}{|c|c|c|} \hline
   Species & \bf Mean & \bf Std.~Dev. \\ \hline
    1 & 3.4 & 1.2 \\
    2 & 5.4 & 0.6 \\ \hline
  \end{tabular}
\end{center}
}
\end{table}

\bibliographystyle{siamplain}

redpajama_set_name: RedPajamaArXiv


text:
Working remotely helped FDA's CDER shape a new and improved workforce, director says. What now?
One surprising consequence of the coronavirus pandemic and its work-from-home zeitgeist? Biopharma companies have been compelled to rethink and perhaps even improve how they operate their businesses.
And the same appears to be true at the FDA, according to commentary from the Cowen Washington Research Group.
The United States government's regulatory body that oversees the industry—namely the FDA's Center for Drug Evaluation and Research—had a net gain of more than 100 employees in 2020 and is tracking to reach that number again this year. The two banner years of recruiting and retention follow a 2019 when CDER had a net gain of just 20 employees.
The difference, according to CDER director Patrizia Cavazzoni? The pandemic.
"We think that the greater flexibility with working remotely, not having people all have to move to Silver Spring [Maryland], has made a difference not only in our ability to hire talent, but also to retain talent," Cavazzoni said during an Alliance for a Stronger FDA webinar on Tuesday.
RELATED: Roll call: How Big Pharma is charting a path back to the office in 2021
Cavazzoni said CDER is examining what its "new normal" will look like, with the understanding that employees have "told us loud and clear that they really enjoy the greater flexibility in how they work."
Similar discussions are underway at companies such as Novartis and Amgen, which are looking for ways to leverage remote work for years to come.
CDER employs roughly 5,600 on a full-time basis. It's the division of the FDA that handles applications for new drugs, manages pharmaceutical manufacturing regulations, collects and analyzes safety data, monitors drug advertising and determines which medications require a prescription.
Upgrading CDER's workforce was already underway before the pandemic thanks to the FDA's 21st Century Cures Act, which Cavazzoni called "a real godsend" and "a game-changer." It was signed in 2016 and has given CDER more salary flexibility.
Cavazzoni's move to CDER, in fact, was a product of the legislation. It allowed the FDA to lure her from Pfizer in 2018 and put her second in command at CDER. Cavazzoni said that 20%, or roughly 1,000, of CDER's employees are covered by Cures authority.
"It's not a perfect solution," Cavazzoni said during the webinar. "We are never going to be able to pay talent as much as the private sector does. But we are more competitive and we are also more able to hire faster."
RELATED: A 'springboard' to the future: Genentech, Sanofi and Alnylam execs sound off on pharma's evolving back-to-work plans
Cavazzoni said that she is monitoring the morale of CDER's busy workforce as the organization plots the future "within the limitations of the government rules." As CDER considers remote work in the future, a priority will be to ensure that technology is in place to "accommodate hybrid learning."

redpajama_set_name: RedPajamaCommonCrawl


text:
Sweetgreen To Expand Plant-Powered Salad Menu To Meet Net-Zero Pledge By 2027
Green BusinessClimate ChangeFood & Nutrition
By Sally Ho Last updated Mar 4, 2021
Courtesy sweetgreen
Sweetgreen, the American fast casual salad chain, is the latest food establishment to hop on the sustainability bandwagon, pledging to become carbon neutral within the next six years. Describing the commitment as the "right thing to do", the popular healthy go-to salad spot says that it will be increasing its range of plant-powered options and introducing new sustainably grown ingredients to its menu.
Sweetgreen says that it plans to slash its carbon footprint to zero by 2027, half of which will come from sustainable sourcing decisions and plant-forward menu development. Commissioning climate firm Watershed to assess its suppliers and calculate the footprint of its ingredients, Sweetgreen says that it'll "introduce even more plant-powered salads and soil-friendly ingredients", listing regenerative kelp and sorghum as some of the new items customers can expect to see in future menu rosters.
Sweetgreen announced its new sustainability commitment to go carbon neutral in six years. (Source: Sweetgreen)
With the food system driving 25% of global greenhouse gas emissions, the time for change is now.
Nicolas Jammet, Co-Founder & CCO, Sweetgreen
It stops short, however, to eliminate carbon-intensive animal-based products from the menu, with Sweetgreen saying that it'll move instead to "work closely with poultry and cheese suppliers to implement new strategies for carbon minimisation".
Commenting on its new climate commitment, Sweetgreen co-founder and chief concept officer Nicolas Jammet said: "Simply put, we believe it's the right thing to do for our business and for the planet. With the food system driving 25% of global greenhouse gas emissions, the time for change is now."
Scientific research has shown that the number one driver of food system emissions is animal agriculture, with the industry alone contributing more than 18% of the world's carbon footprint, exceeding that of all forms of global transportation combined.
Sweetgreen says it will introduce more plant-powered salads and soil-friendly ingredients as part of its net-zero plan. (Source: Sweetgreen)
Simply put, we believe it's the right thing to do for our business and for the planet.
A recent G20 food footprint report by Oslo-based nonprofit EAT found that if the world's largest economies switched to a flexitarian diet – lowering the consumption of meat and dairy – as much as 40% of the total global carbon budget for food could be freed up.
Other measures that Sweetgreen plans to take to reach its six-year goal include optimising the building materials used in its restaurants and furniture, investing in renewable energy, and "double clicking into more granular decisions around things like equipment and refrigeration". The chain also says it will invest in meaningful offsetting projects, once it has reduced its footprint "as much as possible".
Watershed, the firm that Sweetgreen has partnered with to assess its carbon footprint, says that at the moment, the chain's menu is already around 30% less carbon-intensive than the average U.S. diet.
"Their commitment to decrease their greenhouse gas intensity by 50% and become carbon neutral is setting a new bar for the industry," commented Watershed co-founder Taylor Francis.
Sweetgreen plans to construct greener restaurants powered by clean energy. (Source: Sweetgreen)
Similar moves have been made by Sweetgreen's rival chain Just Salad, who has pledged to put carbon labels on all of its menu items indicating its footprint, making it the first U.S. restaurant chain to do so. Just Salad has also made moves to remove beef, one of the most carbon-intensive of all animal products, replacing it with plant-based Beyond Meat, and has since committed to add more plant-based cheese options too.
Their commitment to decrease their greenhouse gas intensity by 50% and become carbon neutral is setting a new bar for the industry.
Taylor Francis, Co-Founder, Watershed
While Sweetgreen has only so far revealed that it'll make its menu more plant-centric, adding little details as to whether plant-based alternatives will be added, Jammet reiterated that the chain will be implementing a whole slew of step-by-step measures to make their actions count.
"We know that real change doesn't happen overnight – it's all the steps in between, the little moments that can lead to a big impact. That's why we're making this commitment."
All images courtesy of Sweetgreen.
corporate environmental responsibilitycorporate sustainability initiativesgreen businessplant based foodplant based trendsustainable businessvegan foodvegan food trend
Sally Ho 2721 posts
Sally Ho is Green Queen's resident writer and lead reporter. Passionate about the environment, social issues and health, she is always looking into the latest climate stories in Hong Kong and beyond. A long-time vegan, she also hopes to promote healthy and plant-based lifestyle choices in Asia. Sally has a background in Politics and International Relations from her studies at the London School of Economics and Political Science.
Unsustainable Salmon Farming Cost US$47 Billion In Losses Since 2013, Report Finds
Interview: Vivi The Super Vegan Author 'I Realised That Kids Books Aren't Vegan Friendly'
Food Culture Expert Eve Turow-Paul of Food For Climate League: 'Small Changes To Menu Language Can Dramatically…


The Moon - Earth's only natural satellite.
The Moon makes a complete orbit around the Earth in average every 29.53 days. The Moon is in synchronous rotation, which means it rotates about its axis in about the same time it takes to orbit the Earth. This results in it keeping nearly the same face turned towards the Earth at all times. The diameter of the Moon is 3476 km (2,160 mi) and it's distance from Earth is about 30 Earth diameters or 384.000 km (238,857 mi). With its gravitational force the moon influences the Earth e.g. the tides on the oceans. More about the Moon at Wikipedia.
The following graphics simulates the Sun shining from the right hand side. From 1=New Moon to 5=Full Moon (waxing) and 5 back to 1 again waning there is a monthly cycle of the Moon available to watch from Earth. Here you may find more Infos about the Moon phase.
The following graphic symbol is used to show the state of the current moon cycle.
Moon is about 51 %.
Générer cette page a eu 0.0083 secondes.


Goebel has a long history of producing quality figurines, based on the artwork of many different artists, but are best known for their Hummel figurines. Hummel figurines are based on the artwork of Sister Maria Innocentia Hummel (also known as Berta Hummel). However, Berta Hummel was not the only Hummel that Goebel partnered with.
Eleonore Elisabeth Maria Hummel (Lore Hummel) is the artist who inspired and helped create the Goebel Blumenkinder and Kinderland lines of figurines. Due to their existing relationship with Sister M.I.Hummel and the possibility of confusion, Goebel referred to Lore Hummel simply as "the artist known as Lore".
Prior to her partnership with Goebel, Lore Hummel was already a well-known and respected Bavarian artist and author, creating and illustrating numerous children's books.
Lore was born in 1915 in Karlsruhe. She moved with her parents to Nuremberg in 1923, where she went to school and met her future husband, Franz-Xaver Hummel. She was married in 1937 and in January of 1939 she gave birth to her son Franz, who has become one of the most famous German contemporary composers and pianists.
In 1944, due to bombing by the Allied Forces, Lore and her family moved to Altmannstein in the Altm?hlt valley. Altmannstein remained her home until she passed away in November 1997 at the age of 82.
Lore always had an interest and aptitude for drawing and painting from an early age, but refined her skills through formal training with Professor Joseph (Sep) Ruf at the Academy of Fine Arts in Nuremberg, and later through private tuition with illustrator Professor Wolfgang Felten in Munich.
Lore's trademark style reflected her love of children, music and nature. Her artistry portrayed rosy-cheeked young children engaged in carefree and playful pursuits. As well as writing and illustrating over 90 children's books, Lore has also published anthologies of songs, poems and fairy tales and designed calendars and postcards.
Working together with Goebel Master Sculptor Gerhard Bochmann, Lore designed the wonderful Blumenkinder series of limited edition figuriness and the subsequent Kinderland series.
Blumenkinder (Flower Children) figurines were issued from 1966 to 1982. 58 different figurines exist in this series, produced in seven editions. Each figurine was produced in a limited edition of 2,000 pieces.
The Kinderland (Child Country) series was introduced in 1978 as an open edition series (i.e. not in a limited quantity). This series featured young children in pastoral and nativity settings. 24 figurines were issued in the Kinderland series.
A series of 6 annual Kinderland plates featuring Lore's artistry was introduced by Goebel in 1982 and ran through 1987.
To support the interest and sale of Goebel products, Lore made one trip to the U.S. during the summer of 1982. She spent time at the Goebel Gallery in Tarrytown, NY, where she was happy to meet collectors, demonstrate her talents, answer questions and sign figurines.
A newer series of 14 polyresin figurines, based on Lore Hummel's artwork and introduced in 2006, is currently produced by Fiebiger Floristik GmbH, under license from Kai Hummel. These figurines can be seen, along with Lore's books and postcards at www.lorehummel.de.


The Free the Spirit Festival is Saturday, September 30th but the St.Amant Foundation invites everyone to the Festival Kickoff Lunch Party on Wednesday, September 13th from 11am – 1pm! Habanero Sombrero Food Truck and PimpMyRice Inc will be there! Plus we'll have Mini Donuts and free slushies. It all happens in our parking lot at 440 River Road.
Bring your cash or debit card, have some lunch and then enjoy lawn games, music, snap some photos at our a Polaroid photo booth, purchase your Festival lottery tickets and learn more about the Festival happening at the end of the month.


TORMENTA4UX RUBBYY100 jasmin. EllyJ. loriStyles. TheLadyAnnita.
AmazingShemalexxVanillaRedmarkyhot18CrystalBoxwood .KinkyBreeObedientSluttvalentynnahotcrazylips .SweetAzuremarkyhot18SophieRagazziSiberiancutie .InaHot92Daniellasweet7hotcrazylipsArianaGrant .ClassybischHothornyevelynSweetAnuskaRealHOTCUMMERduo .HotManStrongthebestlatin88ZharyckTSObedientSlutt .SexyPandorradayibigtitsAmiHoneyEllyJ .hotslaySexyPandorraKinkyBreedayibigtits .SEXXRulerALEXAImaniaEllyJhotslay .TeddyNoyzSEXXRulerALEXAHothornyevelynpyrostorm321 .ObedientSluttEmilyhornyxxxVincentFreshKleopatraXHot .
XLilyBad19XAnnaKudikovakatty33FlorenceFay .MissLexxxykatty33dayibigtitspyrostorm321 .ElectraYonkiFlorenceFayFinishCoupleYannickboy18 .KleopatraXHotHOTVENUZLANCEundeniablydiffeHotManStrong .AnnaKudikovahotcrazylipsHonestxLinsyInaHot92 .DaisyMunXpyrostorm321FlorenceFayTallBlondBoy .pyrostorm321Daniellasweet7pyrostorm321DomixJayden .lissabetaFloryDreamsAnnaKudikovaLaciHurst .valentynnaEllyJRealHOTCUMMERduoxPREZYLE .ArianaGrantBeeRoselissabetahotcrazylips .SammyBerryHarmonyBoneShanonMatbigbum .


Click here to be notified by email when Soap & Disposable Towel Set becomes available.
An Oprah Favorite! This is such a cool soap & towel set that we know anyone will love. This comes with 3 soaps and 12 embossed disposable, linen-like towels. These packs are so much fun and perfect for gift giving. What a great set for bridesmaids, hostess gifts, wedding gifts, housewarming gifts, birthday gifts, for Mother's Day & so much more!
Carving Technique: These are hand carved soaps...versus being engraved. You will not believe how much better the quality is when they are done this way - personalization is deep, crisp and defined. Nothing else even comes close.
12 Towels: mbossed in yoru choice of color. Disposable. Size of towel all the way layed out is 13x16.75" - but folded like you will have it displayed is about 4.25"x8.25" - a perfect size.
Personalization: Please type initials in exact order to appear.
Box Description: Exterior is a beautiful ivory & tan floral brocade.
These ship FAST in 2-3 days! Sorry, no gift wrap but we can include a gift note. For custom logos or private label soaps, please contact helpdesk@inthisveryroom.com.
Aqua Mineral: Antique White in color (Safe for Sensitive Skin & Noses – water based fragrance) simple, gentle & clean.


text:
National Weather Service confirms tornado damaged school in Kershaw County
KERSHAW, S.C. — High winds caused damage to buildings and buses at Kershaw County school Saturday evening.
Kershaw County Sheriff Lee Boan said that the damage will take days to repair and could lead to additional off-days for students.
RELATED: Severe storms knock out power to thousands in the Midlands, cause damage
The school district has been notified and will be assessing the damage today.
As of 12:20 p.m. Sunday, a team from the National Weather Service was on site to survey the damage and confirmed the damage is the result of a tornado. The team is still evaluating the track but confirmed a high-end EF-2 Tornado with estimated winds around 130 mph. Radar indicated there were wind gusts as high as 90 miles an hour Saturday night.
The survey team will also be headed to Fairfield County and Lancaster County to examine damage.
The Kershaw County Sheriff's Office is asking for all spectators to stay off North Central High School property until qualified personnel can complete their assessments and make the area safe and secure. KCSO deputies will be patrolling the property.
The South Carolina Emergency Management Division also reported damage in Anderson County with heavy debris and trees blocking the roadway; and in Chester County where trees were reported down, one outbuilding destroyed and a home struck by a tree.
Posted in Charlotte North Carolina News, North Carolina News
New roof, enduring history: Burgwin-Wright House renovated after Florence
Tornado Causes Major Damage At South Carolina High School


text:
Music theory for the guitar made simple. It all starts with the major scale. If you want to get a basic understanding of music theory and guitar theory, then start with the source. The major scale is what all other scales and chords get referred to or compared to. If you go to a piano and start with the c note and play all of the white keys all the way up to the next c note you have played a c major scale. There are seven separate notes in the scale, eight if you count the fact that you start and end on the same note. The thing that can make scales complicated, is if you try and memorize each note in each scale. It would be a big undertaking to remember each note in the (a) major scale and the (a flat) major scale and the (d) major scale and (c sharp) major scale. Instead I will simplify it to get you started, later if you want to be a more well rounded musician you can memorize all the notes. Instead memorize the pattern. When playing notes you have a whole step or a half step. This is why pianos have black keys. The black keys are half steps between each note. If you play a g and go a half step up to the black key you are now playing a g sharp, (or a flat). If you play the g then go a half step down to the next black key you are playing a g flat (or f sharp). Remember in music notes only go from a to g. If you look on the piano you will not see a black key between the (e) and the (f) notes. To go from e to f is only a half step. The same is true with the step from b to c. On the guitar each fret is only a half step. So if you play the top (a) note and slide your finger one fret to the right, (considering your playing a right handed guitar) your now playing an (a) sharp. You would need to slide two frets over to play the (b). So again to make it simple or get down to the basic part of the scale lets look at the pattern. The pattern for any major scale is this, whole step, whole step, half, whole, whole, whole, half. It does not matter which note you start on, If you follow this pattern you are playing a major scale. The first note you start on is the root note or the name of the scale.
Start this pattern on a d note and you are playing a d major scale. Start the pattern on an f# and you are playing an f# major scale. Look at the scales on this web site and look at the pattern for the major scale then apply it to the neck of your guitar. To play a minor scale you will play a pattern close to that of the major scale but with two changes. The third note and the seventh note will be flatten. So your pattern will be this whole, half, whole, whole, whole, half, whole. The first note you play is the root note, so if it is a b, you are playing a b minor and so on.
Now lets turn the scale into a chord. A chord is a group of notes you get from the scale that you can play at the same time or two at the same time or go up and down the chord in an arpeggio. So take a g major scale and take the first, third, and fifth note out of it and you have a g major chord. You will be playing the g, b, d, notes. You may notice that when playing this chord on the guitar you are playing six strings, but the chord only has three notes. The other strings are playing the same notes as well but at different octaves. The first note in the chord is also called the root note.
If we take the g minor scale and turn it into a g minor chord, we still take the first, third, and fifth notes out of it. If you look at the scales or look at the pattern you'll remember the third note gets flatten a half step. This is what makes the minor chord sound different. We are now playing g, b flat, d. Now here is where we tie it all together. You don't have to play a scale from start to end or up and down. Once you have memorized the pattern, play the notes in different order within the pattern and at different lengths. Now you are making up your own melodies. Play along with others and have them strum a g chord while you play a melody you have made up using the g scale pattern. If they play a d minor chord, you can play a melody in the d minor scale pattern. Please post questions or comments. It is my hope to make this more simple for the beginner. I tried learning the hard way and gave up many times as a result.


Costa Rica Vacation – Top 10 Places To Visit!
The Arenal Volcano National Park is one of the 10 most active volcanoes in the world. Visitors can book tours that explore the volcano and surrounding areas.


Read more about Science Museum short film - Why do clouds have names?
Message Impossible - How Would You Redesign The UK's Weather Forecast?
Read more about Message Impossible - How Would You Redesign The UK's Weather Forecast?


text:
Oregon State Hospital Memorial Honored by Public Art Network Year in Review
Americans for the Arts, the nation's leading nonprofit organization for advancing the arts and arts education, honors 31 outstanding public arts projects created in 2014, including the Oregon State Hospital Memorial, through the Public Art Network (PAN) Year in Review program, the only national program that specifically recognizes the most compelling public art. The works were chosen from more than 300 entries across the country and recognized today at Americans for the Arts' 2015 Annual Convention in Chicago.
Read more about Oregon State Hospital Memorial Honored by Public Art Network Year in Review
Artists, supporters and those interested in the creative economy of Clatsop County invited to Mechanisms to Support Visual Artists conversation May 6 in Astoria
The current conditions, needs and aspirations of artists, patrons and the creative economies of Clatsop County will be the focus of an in-depth community conversation hosted by the Oregon Arts Commission from 5:30 to 7:30 p.m. on Wednesday, May 6, at Clatsop Community College.
Read more about Artists, supporters and those interested in the creative economy of Clatsop County invited to Mechanisms to Support Visual Artists conversation May 6 in Astoria
$89,375 awarded for three new 'World of Work' projects
Training as sound technicians at professional concerts; Designing, creating and installing public art; Learning the art of transmedia storytelling – telling a single story via multiple media – at Warm Springs Academy. Dozens of middle- and high-school aged underserved Oregon youth will experience creative industry careers alongside working professionals in three new "Connecting Students to the World of Work" projects announced today by the Oregon Arts Commission.
Read more about $89,375 awarded for three new 'World of Work' projects
Gypsy Prince of Springfield named Oregon's Poetry Out Loud champion
Gypsy Prince, a three-time school champion from Springfield's Academy of Arts and Academics, is the 2015 Oregon Poetry Out Loud state champion. Prince, a 17-year-old senior, was one of nine regional finalists who competed in Saturday's state Poetry Out Loud contest at Willamette University. She advances to the national Poetry Out Loud competition, scheduled for April 27-29 in Washington, D.C.
Read more about Gypsy Prince of Springfield named Oregon's Poetry Out Loud champion
2015 Poetry Out Loud finalists named; state contest set for March 14
Regional Poetry Out Loud competitions held in Beaverton, Medford and Salem March 7 identified the nine Oregon high school students who will compete for the state title at Willamette University on Saturday, March 14. Three finalists were selected from each regional competition; a record 38 schools throughout the state participated.
Read more about 2015 Poetry Out Loud finalists named; state contest set for March 14
By James Bash | Friday, March 6, 2015
Packed like salmon in a tin can, the audience at Astoria's Pier 11 was already overflowing into the lobby on an unusually sunny Saturday afternoon in February. Everyone listened intently as local poet Nancy Cook launched into one of her newest songs about going clamming during low tide. It didn't take much urging from Cook to get to crowd to join her on the refrain, "It's a minus tide this mornin'; I'll meet you there."
Read more about Arts Build Communities: Fisherpoets put the fishing life into verse in Astoria
Eight visual arts institutions receive $95,690 in Art Acquisition Funds
Eight Oregon arts organizations will receive $95,690 in funds from The Ford Family Foundation's Art Acquisition Program, the Oregon Arts Commission announced today. The Art Acquisition Fund, managed by the Arts Commission, helps visual arts institutions with publicly accessible collections acquire seminal works by Oregon visual artists. The effort preserves public access to great works and supports artists and the institutions that sustain their work.
Read more about Eight visual arts institutions receive $95,690 in Art Acquisition Funds
2015-03-13 16:00:00 / Salem
Read more about 2015-03-13 16:00:00 / Salem
Plan now for Oct. 2 Southern Oregon Arts Summit
As previously announced, the 2015 Arts Summit will take place in Southern Oregon. Now scheduled in Ashland on Friday, Oct. 2, the summit will be held at the newly renovated Ashland Hills Hotel & Suites. The agenda will address both statewide and regional topics.
Read more about Plan now for Oct. 2 Southern Oregon Arts Summit
Record 38 Oregon high schools participate in Poetry Out Loud
This year a record 38 Oregon high schools are participating in Poetry Out Loud, a national recitation contest organized by the Oregon Arts Commission in collaboration with the NEA and the Poetry Foundation.
Read more about Record 38 Oregon high schools participate in Poetry Out Loud


text:
"Questioning Justice Barrett's Questions in Fulton v. City of Philadelphia" by Patrick Hornbeck
Dec 9, 2020 Canopy ForumEssaysDecember 2020, Hornbeck, Patrick, Politics
When Amy Coney Barrett took the bench November 4 for her third day of oral arguments as an associate justice of the U.S. Supreme Court, the questions she asked seemed to reveal both her loyalty to her late mentor, Justice Antonin Scalia, and her ongoing interest in abortion rights. The case at hand, Fulton v.…
"When Campaign Attacks Become Theological" by Deirdre Jonese Austin
Nov 23, 2020 Canopy ForumEssaysAustin, Deirdre Jonese, November 2020, Politics
Georgia State Capitol, Atlanta, Georgia. Source: DXR, Wikimedia CC BY-SA 4.0. Currently, all eyes are on Georgia as the state prepares for the runoff elections on January 5th, which will determine which party has control of the Senate. While the battles between Jon Ossoff and David Perdue and Raphael G. Warnock and Kelly Loeffler are…
"Constitutional Grace: Securing the Blessings of Liberty Through Dignity and Forgiveness" by William E. Thro
Nov 18, 2020 Canopy ForumEssaysNovember 2020, Politics, Thro, William E.
We the People … in order … to secure the blessings of liberty to ourselves and our posterity do ordain and establish this Constitution for the United States of America. Because neither the People nor their leaders are angels, the Constitution reflects a Calvinist theological perspective and embodies "obsessive distrust of government — all government…
"Universalist Response to Particularistic Regression: Political Philosophy behind the Pope Francis Encyclical Fratelli tutti" by Mikhail Minakov
Nov 17, 2020 Canopy ForumSeriesMinakov, Mikhail, November 2020, Politics
The modern history of humanity is driven by several contradictions, one of them being the conflict between universalism and particularism. In a nutshell, numerous universalist frameworks are grounded on the idea that norms and values relate equally to all humans, to all living creatures, or — even more widely — to all forms of being.…
"What We Owe the Democracy: Martin Luther King, Jr., the Right to Vote, and the Call to Civic Duty" by Atiba R. Ellis
Nov 3, 2020 Canopy ForumEssaysEllis, Atiba R., November 2020, Politics
The right to vote is a contested concept in American society. The choices made by elected federal and state governments, on the behalf of "We the People" as to who is not included in our democracy both construct American citizenship and reflect American democratic values. The contest around the right to vote, called "the voting…
"How the Divisive Nature of Religion Could Unify our Divided Politics" by Mingyu Jun
Nov 2, 2020 Canopy ForumEssaysJun, Mingyu, November 2020, Politics
In the 2016 presidential election, 81% of white evangelical Christians voted for the GOP candidate Donald Trump, while only 16% gave their votes to Hilary Clinton.This significant margin has the potential to grow further in the 2020 election with the voting polls indicating that 82% of evangelicals would vote for Trump over Biden. To close…
"Trump, Biden, and Religious Claims in a Secular Space" by John E. King
Nov 2, 2020 Canopy ForumEssaysKing, John E., November 2020, Politics
"USA Bible" by Pastor Robert / Flickr / CC BY-ND 2.0 At a campaign event in early August 2020, President Trump made statements about his Democratic rival for the presidency, Joe Biden, that went beyond his typically disparaging remarks about Biden's policies, cognitive abilities, or political record. Instead, Trump's comments were decidedly religious. Trump informed…
"The Cost of Florida's Ballot" by Benjamin Plener Cover
Oct 30, 2020 Canopy ForumEssaysCover, Benjamin Plener, October 2020, Politics
The Leadership Conference on Civil and Human Right / Flickr / CC BY-NC 2.0 Florida is the most populous battleground state, famous for high-stakes contests, razor-thin margins, and trouble with electoral administration. Back in 2000, George W. Bush won Florida by 1,784 votes. In Florida, perhaps more than in any other state, every vote counts.…
"Tax Law, Religion, and Justice: An Exploration of Theological Reflections on Taxation" by Allen Calhoun
Oct 27, 2020 Canopy ForumLiterature HighlightsCalhoun, Allen, October 2020, Politics
Tax Law, Religion, and Justice: An Exploration of Theological Reflections on TaxationAllen Calhoun An Introduction by Allen Calhoun This is an Accepted Manuscript of a book chapter forthcoming from Routledge/CRC Press in Tax Law, Religion, and Justice: An Exploration of Theological Reflections on Taxation on March 8, 2021, available online here. Why does the institution of taxation occupy…
"Religious Questions: Relevant, Legitimate, and Impossible" by Paul Horwitz
Oct 7, 2020 Canopy ForumEssaysHorwitz, Paul, October 2020, Politics
For the writer on law and religion, the "religious test" question is a gift that keeps on giving. Over two centuries after we enshrined the rule against religious tests as a "Qualification to any Office or public Trust under the United States" — and, arguably, decades if not more after the Religious Test Clause had…
Conviction (18)
Literature Highlights (20)
Series (104)