text:
\subsection{Proof \cref{thm:structure_W_H}}

\label{sec:p_sructure_W_H}

Let $2\alpha_2-(\alpha_1+\alpha_2)^2 = 0$ and $T(W)=T_d$ is partitioned as in \eqref{eq:structure_T}. Since $T_d$ has to be skew-symmetric (see also \Cref{thm:existence_svd}), the eigenvalues of $T_d$ are purely imaginary. Then condition \eqref{eq:problem} with $T_{12} = -T_{21}^\top = -R\in\mathbb{R}^{n\times n}$, where $R$ is arbitrary, reads

\begin{align}
	 - I &= g'(z)f(z)R^\top-f'(z)g(z)R \nonumber \\
	&+ f'(z)f(z) T_{11} + g'(z)g(z) T_{22}.
	\label{eq_p2:1}
\end{align}

Since, $T_{11} = - T_{11}^\top$, $T_{22} = -T_{22}^\top$, w.l.o.g. we can express $R$ as  $R = I + \tilde{R}$, where $\text{diag}(\tilde{R}) = 0$. Hence, it has to hold
\begin{align}
	g'(z)f(z)-f'(z)g(z) = -1,
	\label{eq_p2:2}
\end{align}
which is satisfied by the generating functions \cite[Theorem 1]{GRUSHKOVSKAYA2018151}
\begin{align}
	g(z) = -f(z)\int \frac{1}{f^2(z)} dz.
	\label{eq_p2:3}
\end{align}
Assume that generating functions $f$ and $g$ as in \eqref{eq_p2:3} satisfy \eqref{eq_p2:2}, then condition \eqref{eq_p2:1} translates into
\begin{align}
    f'(z)f(z) T_{11}\! +\! g'(z)g(z) T_{22} \nonumber \\
	= (f'(z)g(z)+g'(z)f(z))\tilde{R},
	\label{eq_p2:3_1}
\end{align}
implying that $\tilde{R} = -\tilde{R}^\top$ holds.
Next we consider three cases:

\\\\

\textbf{Case 1: $\tilde{R} = 0$}. Hence,

\begin{align}
	f'(z)f(z) T_{11} + g'(z)g(z) T_{22} = 0.
	\label{eq_p2:4}
\end{align}
Clearly, \eqref{eq_p2:4} is satisfied by $T_{11} = T_{22} = 0$ with $f$ arbitrary while satisfying $g$ in \eqref{eq_p2:3}, i.e., \eqref{eq:cases_H1} results.

For the sub-case $a^{-1}T_{11} = b^{-1} T_{22} =: Q$ arbitrary skew-symmetric with $a,b\in \mathbb{R}_{>0}$, $f$ and $g$ have to satisfy $af'(z)f(z)+b g'(z)g(z) = 0$, i.e., w.l.o.g. $af^2(z)+b g^2(z) = 1$. Accordingly, with \eqref{eq_p2:3} and $y'(z) = f^{-2}(z)$ it yields

\begin{align}
	y'(z) = a + b y^2(z).
	\label{eq_p2:5}
\end{align}

The unique solution of \eqref{eq_p2:5} with $\phi \in \mathbb{R}$ is

\begin{align}
	y(z) = \sqrt{\frac{a}{b}}\tan\Big(\sqrt{a b} z + \phi \Big)
	\label{eq_p2:5_2}
\end{align}

and therefore with the definition of $y(z)$ and \eqref{eq_p2:3} one reveals \eqref{eq:cases_H2}.

Repeating the above calculations for $a^{-1} T_{11} = - b^{-1} T_{22} =: Q$ arbitrary skew-symmetric yields \begin{align}
	y(z) = \sqrt{\frac{a}{b}}\tanh\Big(\sqrt{ab} z + \phi \Big)
\end{align}

and therefore with the definition of $y(z)$ and \eqref{eq_p2:3} one reveals \eqref{eq:cases_H3}.

The remaining sub-case $T_{11} \neq \pm a T_{22}$, implies that

\begin{align}
	f'(z)f(z) T_{11} = 0\ \ \text{and}\ \ g'(z)g(z) T_{22} = 0,
	\label{eq_p2:5_1}
\end{align}

since $	f'(z)f(z) T_{11} + g'(z)g(z) T_{22} = 0$ must hold.

If $f'(z)f(z) = 0$ and $T_{11}=Q$ arbitrary skew-symmetric, i.e., $f^2(z) = a$ with $a\in\mathbb{R}_{>0}$, it implies that $f(z) = \sqrt{a}$. Hence, with \eqref{eq_p2:3}, $g'(z)g(z) \neq 0$ for all $z\in\mathbb{R}$ yields $T_{22} = 0$, i.e. \eqref{eq:cases_H4} results.

The same argumentation holds for $g'(z)g(z) = 0$ with arbitrary skew-symmetric $T_{22}$ such that \eqref{eq:cases_H5} is recovered.

The circumstance that  $T_{11} \neq T_{22}$ with $T_{11} \neq 0$ and $T_{22} \neq 0$ is not valid due to \eqref{eq_p2:3} and \eqref{eq_p2:5_1}. Specifically, $f'(z)f(z) = 0$ and $g'(z)g(z) = 0$ has to hold; obviously, based on the above cases, $f(z)=\sqrt{a}$ and $g(z)=\sqrt{a}$ are in conflict with \eqref{eq_p2:2}.
\\\\

\textbf{Case 2:} $f'(z)g(z)+g'(z)f(z) = 0$ for all $z\in \mathbb{R}$. Hence,

\begin{align}
    &f'(z)f(z) T_{11}\! +\! g'(z)g(z) T_{22} = 0 \text{ and }
    \label{eq_p2:6_00}
    \\
    &f'(z)g(z)+g'(z)f(z) = 0
    \label{eq_p2:6_01}
\end{align}

has to be satisfied. Clearly, \eqref{eq_p2:6_00} is satisfied by $T_{11} = T_{22} = 0$, where \eqref{eq_p2:6_01} implies $-af(z)g(z)=1$ with $a \in \mathbb{R}_{>0}$. Accordingly, with \eqref{eq_p2:3} and $y'(z) = f^{-2}(z)$ it yields
\begin{align}
     y'(z) = ay(z) .
     \label{eq_p2:6_1}
\end{align}
 The unique solution of \eqref{eq_p2:6_1} with $c\in\mathbb{R}$ is
\begin{align}
	y(z) = e^{az}+c
	\label{eq_p2:7}
\end{align}
and therefore with the definition of $y(z)$ and \eqref{eq_p2:3} one gets \eqref{eq:cases_H6}.

The sub-cases where $a^{-1}T_{11} = b^{-1} T_{22} =: Q$ arbitrary skew-symmetric or  $T_{11} \neq \pm a T_{22}$ as discussed for Case 1 are not valid. With the same approach as above, i.e., $y'(z) = f^{-2}$ it yields to $y'(z) = 0$ and therefore no solution for $f(z)$ (and $g(z)$) can be found or \eqref{eq_p2:6_00} and \eqref{eq_p2:6_01} are not satisfied as discussed in the last paragraph of Case 1, respectively.

\\\\

\textbf{Case 3: $\tilde{R} = -\tilde{R}^\top \neq 0$}. Hence,
\begin{align}
    f'(z)f(z) T_{11}\! +\! g'(z)g(z) T_{22} \nonumber \\
	= (f'(z)g(z)+g'(z)f(z))\tilde{R}.
	\label{eq_p2:case3}
\end{align}

Clearly, $T_{11} = T_{22} = 0$ is not valid, since $(f'(z)g(z)+g'(z)f(z))\tilde{R} \neq 0$ in this last case.

For the sub-case \corrected{$a^{-1}T_{11} = b^{-1} T_{22} =: Q = -Q^\top \neq 0$ with $a,b \in \mathbb{R}_{>0}$, it has to hold that $a f'(z)f(z) + b g'(z)g(z) = c(f'(z)g(z)+g'(z)f(z))$ and $\tilde{R} = cQ$ with $c \in \mathbb{R}\backslash{\{0\}}$, i.e., w.l.o.g. $\nicefrac{a}{2}f^2(z) + \nicefrac{b}{2} g^2(z) -c f(z)g(z)= 1$. Accordingly, with \eqref{eq_p2:3}} and $y'(z) = f^{-2}(z)$ it yields
\begin{align}
     y'(z) = \frac{a}{2} + cy(z) + \frac{b}{2} y^2(z).
	\label{eq_p2:6}
\end{align}
The unique solution of \eqref{eq_p2:6} with $\phi \in \mathbb{R}$ is
\begin{align}
	y(z) = \frac{\sqrt{ab-c^2}}{b} \tan\Big(\sqrt{ab-c^2} z + \phi \Big)
	\label{eq_p2:7}
\end{align}
and therefore with the definition of $y(z)$ and \eqref{eq_p2:3} one reveals \eqref{eq:cases_H7}.

The remaining sub-case $T_{11} \neq a T_{22}$, implies that

\begin{align}
&f'(z)f(z) T_{11} - (f'(z)g(z)+g'(z)f(z)) \tilde{R} = 0 \quad \text{and}
\nonumber \\
&g'(z)g(z) T_{22} = 0, \quad \text{or}
\label{eq_p2:8}
\\
&g'(z)g(z) T_{22} - (f'(z)g(z)+g'(z)f(z)) \tilde{R} = 0 \quad \text{and}
\nonumber \\
&f'(z)f(z) T_{11} = 0,
\label{eq_p2:9}
\end{align}

since \eqref{eq_p2:case3} must hold. However, we show in the sequel that \eqref{eq_p2:8} and \eqref{eq_p2:9} lead to no new solution or is not valid, respectively.

For \eqref{eq_p2:8}, \corrected{$a^{-1}T_{11} = c^{-1} \tilde{R} =: Q = -Q^\top \neq 0$} with $a \in \mathbb{R}_{>0}$ and $c \in \mathbb{R}\backslash{\{0\}}$, it has to hold that \corrected{$a f'(z)f(z) = c(f'(z)g(z)+g'(z)f(z))$}, i.e., w.l.o.g.

\begin{align}
    \frac{a}{2} f^2(z) -cf(z)g(z)= 1.
    \label{eq_p2:9_1}
\end{align}

Accordingly, with \eqref{eq_p2:3} and $y'(z) = f^{-2}(z)$ it yields

\begin{align}
     y'(z) = \frac{a}{2} + cy(z).
	\label{eq_p2:9_2}
\end{align}
The unique solution of \eqref{eq_p2:9} with $d \in \mathbb{R}$ is
\begin{align}
	y(z) = de^{cz}-\frac{a}{2c}
	\label{eq_p2:10}
\end{align}
and therefore with the definition of $y(z)$ and \eqref{eq_p2:3} one gets $f(z) = (cd)^{-\nicefrac{1}{2}}\text{exp}(-\nicefrac{c}{2}\, z)$ and $g(z) = -(d/c)^{-\nicefrac{1}{2}}\text{exp}(\nicefrac{c}{2}\, z)$ such that $T_{22} = 0$. However, $f,g$ satisfying \eqref{eq_p2:9_1} for all $z\in \mathbb{R}$ only for $a=0$, hence $T_{11}=0$, i.e, the same result as \eqref{eq:cases_H6}.

For \eqref{eq_p2:9}, \corrected{$b^{-1}T_{11} = c^{-1} \tilde{R} =: Q = -Q^\top \neq 0$} with $a \in \mathbb{R}_{>0}$ and $c \in \mathbb{R}\backslash{\{0\}}$, it has to hold that \corrected{$b g'(z)g(z) = c(f'(z)g(z)+g'(z)f(z))$}, i.e., w.l.o.g.

\begin{align}
    \frac{b}{2} g^2(z) -cf(z)g(z)= 1.
    \label{eq_p2:10_1}
\end{align}

Following the procedure above yields $f(z) = (2)^{-\nicefrac{1}{2}} c^{-1} (b \text{exp}(\nicefrac{c}{2}\, z + \nicefrac{d}{2})-\text{exp}(-\nicefrac{c}{2}\, z - \nicefrac{d}{2}))$ and $g(z) = 2^{\nicefrac{1}{2}}b^{-1}\text{exp}(-\nicefrac{c}{2}z - \nicefrac{d}{2})$ with $d\in\mathbb{R}$. However, \eqref{eq_p2:10_1} is only satisfied for $d=-cz-\text{ln}(b)$ which leads to $f(z) = 0$ and therefore no valid solution.

Notice that every feasible structure of the skew-symmetric matrix $T_d$ is discussed above case by case, and the differential equations arising in the analysis are solved uniquely. Hence, we believe that the list of triples $(T_d,f,g)$ in \cref{thm:structure_W_H} for $2\alpha_2-(\alpha_1+\alpha_2)^2=0$ and $T_d$ skew-symmetric is essentially exhaustive, save for some scaled version of the presented cases.

$\hfill \small{\square}$


TITLE: How to factor $y = x^5 + 20x^2 + 5$?
QUESTION [5 upvotes]: How would I factor to solve for x?
$x^5 + 20x^2 + 5=0 $?
Do I use synthetic division?
Is there a faster/easier way?
Do I have to keep plugging in numbers to see if they equal to zero?
Thanks! I'm not asking for full solutions if you don't want to share (but that would be nice) just opinions on what I should do.

REPLY [5 votes]: It’s irreducible over the rational numbers, by the Eisenstein Criterion. So in particular, it doesn’t have a rational root. You can get an approximate root by hand in various ways, and the method I would use is the Newton-Raphson method.


TITLE: Two quadratics with one root opposite to each other
QUESTION [5 upvotes]: Say we have two quadratic equations as follows,
$$x^2-2mx-1=0$$
$$x^2+(m+3)x-4=0$$
where $m$ is real, such that the two equations have one root that is the additive inverse of each other.
(For example, one may have root $x=3$ and the other would have $x=-3$ as one of its roots.)
Find m.

My only thought is using the quadratic formula and compare the two equation's roots, it'll get the job done, but takes a lot of work.
I suppose there's a method to modify one of the equations so that the new equation from modification has one same root as the equation we didn't modify.
That way, we can just subtract them and solve for a much simpler equation.

REPLY [4 votes]: If $\alpha$ is a root of $f(x) = x^2 - 2mx - 1$, then $-\alpha$ is a root of $f(-x)$ or $x^2 + 2mx - 1$.
Now $x^2 + 2mx - 1$ and $x^2 + (m+3)x - 4$ must have a common root, or...


TITLE: Let $F$ be a field. Show that the field of quotients of $F$ is ring isomorphic to $F$. My Attempt Shown
QUESTION [1 upvotes]: Let $F$ be a field. Show that the field of quotients of $F$ is ring isomorphic to $F$.

Attempt: Let $F'$ be the field of Quotients of the field $F$.
Let $\Phi:F \rightarrow F'$ such that $\Phi(x)=x/1$
$x/1$ refers to the equivalence class containing the element $(x,1) ~;~x,1 \in F $ . The equivalence relations satisfies $a/b = c/d$ iff $ad=bc$
Operation Preservation : Then, $\Phi(x+y)= (x/1) + (y/1) = \Phi(x) + \Phi(y)$
$\Phi(xy) = (xy)/1 = (x/1) . (y/1)$
These operations are consistent with the definition of addition and multiplication defined for the elements of the field of quotient and hence, both addition and multiplication are preseerved.
One - One nature : If $x/1$ denotes the equivalence class containing $(x,1)$. Now, how do I prove the one-one nature. I have read about equivalence classes, but, don't possess much expert intuition. So, here's my attempt :
$(x,1) \in x/1 \implies (1,x) \in x/1 \implies (x,x) \in x/1$

Also, the equivalence classes are disjoint, $\implies$ No other equivalence class shall contain elements of the form $(x,1)$ or $(1,x)$ or $(x,x) \implies x/1 =y/1~~ \forall ~~x,y \in F$ . How does this $\implies x=y$ which is necessary for the one-one condition?

Onto Nature : Corresponding to any element $(x/1)$ in the quotient of field, we can find $x \in F$. Hence, the mapping is onto.
Is my attempt correct?
Thank you for your help.

REPLY [1 votes]: There are a few issues here, but good effort.
First of all, it might clarify things to refer to elements of the field of fractions via their representative ordered pairs, rather than as fractions. I know that you have explained your notation, but when you write elements of the fraction field in this way, they look like elements of $F$ which may lead to confusion.
Your proof that $\Phi$ is a ring homomorphism looks good, although you could spell out explicitly how the rules of addition/multiplication in the field of fractions are applied.
Injectivity (one-to-one): you have to show that $\Phi(x)=\Phi(y) \implies x=y$. This is equivalent to proving that if $(x,1)=(y,1)$ then $x=y$.
Surjectivity (onto): you need to show that for any element $(x,y)$ in the field of fractions, there is some element $p \in F$ so that $\Phi(p)=(x,y)$. (This is where it becomes important not to refer to elements of the quotient as fractions -- you're going to need to divide to get the element $p$!).


text:
\begin{document}

\maketitle

\begin{abstract}
The \textit{multi-level}, \textit{multi-disciplinary} and \textit{multi-fidelity} optimization framework developed at Bombardier Aviation has shown great results to explore efficient and competitive aircraft configurations.
This optimization framework has been developed within the Isight software, the latter offers a set of ready-to-use optimizers.
Unfortunately, the computational effort required by the Isight optimizers can be prohibitive with respect to the requirements of an industrial context.
In this paper,  a \textit{constrained Bayesian optimization} optimizer, namely the \textit{super efficient global optimization with mixture of experts}, is used to reduce the optimization computational effort. The obtained results showed significant improvements compared to two of the popular Isight optimizers.
The capabilities of the tested \textit{constrained Bayesian optimization} solver are demonstrated on \textit{Bombardier research aircraft configuration} study cases.
\end{abstract}

\section{Nomenclature}

{\renewcommand\arraystretch{1.0}
\noindent\begin{longtable*}{@{}l @{\quad=\quad} l@{}}
$d$ & the number of design variables \\
$m$ & the number of inequality constraints \\
$\Omega$ & the domain of the design variables, i.e.,  $\Omega \subset \mathbb{R}^d$\\
$f$  & the objective function, i.e., $f: \mathbb{R}^d \mapsto \mathbb{R}$ \\
$\bm{g}$ & the constraints function, i.e.,  $\bm{g}: \mathbb{R}^d \mapsto \mathbb{R}^m$\\
$\mu_s^{(l)}$ & the prediction of the Gaussian process of a given function $s: \mathbb{R}^d \mapsto \mathbb{R}$ build with $l$ samples\\
$\sigma_s^{(l)}$ & the uncertainty of the Gaussian process of a given function  $s: \mathbb{R}^d \mapsto \mathbb{R}$ build with $l$ samples \\
$\alpha_f^{(l)}$ & the acquisition function of the objective function $f: \mathbb{R}^d \mapsto \mathbb{R}$\\
$\alpha_{\bm{g}}^{(l)}$ & the feasibility criterion of the constraint functions $\bm{g}: \mathbb{R}^d \mapsto \mathbb{R}^m$\\
$\Omega^{(l)}_{\bm{g}}$ & the approximated feasible domain defined by the feasibility criterion $\alpha_{\bm{g}}^{(l)}: \mathbb{R}^d \mapsto \mathbb{R}^m$ \\
$C_{L_{max}}$ & the low speed maximum lift coefficient
\end{longtable*}}

\input{introduction.tex}

\input{Bombardier_MDO.tex}

\input{optimizers.tex}

\input{BRAC_optim.tex}

\input{conclusion.tex}

\section*{Acknowledgments}
This work is part of the activities of ONERA - ISAE - ENAC joint research group in a context of  a partnership between ONERA, ISAE-SUPAERO and Bombardier Aviation. In addition, we thank Jasveer Singh for his work on the BRAC CMDO specification.

\clearpage
\bibliography{Biblio}

\end{document}

config: arxiv


text:
\section*{Supplementary Material}

\section{Proof of Lemma \ref{lem:Bounded_Q_var}}

Since the estimate $\hat{Q}$ is unbiased, it suffices to show that
$\mathbb{E}[|\hat{Q}|^{2}]$ is bounded. Indeed, we have
\begin{flalign}
\mathbb{E}[|\hat{Q}|^{2}]
=\mathbb{E}\Bigg[\Bigg|\sum_{t=0}^{T}R(s_{t},a_{t}=\pi_{\theta}(s_{t}))\Bigg|^{2}\Bigg]\nonumber
 & \le\mathbb{E}\Bigg[\Bigg(\sum_{t=0}^{T}|R(s_{t},a_{t}=\pi_{\theta}(s_{t}))|\Bigg)^{2}\Bigg]\nonumber \\
 & \le\mathbb{E}\big[(T+1)^{2}U_{R}^{2}\big]

 =U_{R}^{2}\dfrac{1+\gamma}{(1-\gamma)^{2}}, \nonumber
\end{flalign}
where the last equality follows from the fact that $T$ is a geometric
random variable with probability of success $1-\gamma$. Enough said.	\null\nobreak\hfill\ensuremath{\square}

\section{Proof of Theorem \ref{thm:nesterov_result1}}

	Let Assumptions \ref{asm:bounded_gradient}-\ref{asm:smooth} and also Assumption \ref{asm:nesterov_assumptions} be in effect. Also, recall the definitions of the gradients, namely
	\begin{equation} \nonumber
	\nabla J(\theta) = \E_{s\sim \rho^{\pi_\theta}}\left[\nabla_\theta \pi_\theta(s) \nabla_a Q^{\pi_\theta}(s,a)|_{a = \pi_{\theta}(s)}\right],
	\end{equation}

	and

	\begin{equation} \nonumber
	 \hat \nabla J(\theta) =  \E_{s\sim \rho^{\pi_\theta}} \left[\nabla_\theta \pi_\theta(s) \nabla_a Q_\mu^{\pi_\theta} (s, a) \vert_{a = \pi_\theta(s)  } \right].
	\end{equation}

	By definition, and then by Jensen, we have

	\begin{equation}
	\begin{split}
	\| \nabla J(\theta) &-  \hat \nabla J(\theta)\| \\ \nonumber
	&\hspace{-24pt}= \Big\|
	\E_{s\sim \rho^{\pi_\theta}}\left[\nabla_\theta \pi_\theta(s) \nabla_a Q^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)}\right] - \E_{s\sim \rho^{\pi_\theta}}\left[\nabla_\theta \pi_\theta(s) \nabla_a Q_\mu^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)}\right] \Big\| \\ \nonumber
	&\hspace{-24pt}= \Big\| \E_{s\sim \rho^{\pi_\theta}}\left[\nabla_\theta \pi_\theta(s) \nabla_a Q^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)} -\nabla_\theta \pi_\theta(s) \nabla_a Q_\mu^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)}\right] \Big\|\\ \nonumber
	&\hspace{-24pt} \leq  \E_{s\sim \rho^{\pi_\theta}}\left[ \big\| \nabla_\theta \pi_\theta(s) \nabla_a Q^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)} -\nabla_\theta \pi_\theta(s) \nabla_a Q_\mu^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)}\big\|\right].
	\end{split}
	\end{equation}
	Grouping terms, then applying Cauchy-Schwarz, we obtain
	\begin{equation}
	\begin{split} \nonumber
	\| \nabla J(\theta) - \hat \nabla J(\theta)\| &\leq \E_{s\sim \rho^{\pi_\theta}}\left[ \big\| \nabla_\theta \pi_\theta(s)\big( \nabla_a Q^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)} - \nabla_a Q_\mu^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)}\big)\big\|\right] \\
	&\leq \E_{s\sim \rho^{\pi_\theta}}\left[ \left\| \nabla_\theta \pi_\theta(s) \right\| \cdot \left\| \nabla_a Q^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)} - \nabla_a Q_\mu^{\pi_\theta}(s,a)\big|_{a = \pi_{\theta}(s)}\right\|\right]. \\
	\end{split}
	\end{equation}

	By Assumption \ref{asm:bounded_gradient}, we can bound $\nabla_\theta \pi_\theta(s)$ by $B_\Theta$. As such, the remaining expression can be bounded by invoking Assumption \ref{asm:nesterov_assumptions}, yielding

	\begin{equation}
	\begin{split} \nonumber
	\| \nabla J(\theta) - \hat \nabla J(\theta)\| & \leq B_\Theta\E_{s\sim \rho^{\pi_\theta}}\left[ \left\| \nabla_a \left( Q^{\pi_\theta}(s,a) -  Q_\mu^{\pi_\theta}(s,a)\right)\big|_{a = \pi_{\theta}(s)}\right\|\right] \\
	&\leq B_\Theta\E_{s\sim \rho^{\pi_\theta}}\left[ \left\| \nabla_a \left( Q^{\pi_\theta}(s,a) -   \E_\mathbf{u} \left[Q^{\pi_\theta}(s,a+ \mu \textbf{u})\right]\right)\big|_{a = \pi_{\theta}(s)}\right\|\right]\\
	&\leq  B_\Theta\E_{s\sim \rho^{\pi_\theta}}\left[ \left\| \nabla_a \left(\E_\mathbf{u} \left[ Q^{\pi_\theta}(s,a) -   Q^{\pi_\theta}(s,a+ \mu \textbf{u})\right]\right)\big|_{a = \pi_{\theta}(s)}\right\|\right]\\
	&\leq  B_\Theta\E_{s\sim \rho^{\pi_\theta}}\left[ \left\| \E_\mathbf{u} \left[\nabla_a \left( Q^{\pi_\theta}(s,a) -   Q^{\pi_\theta}(s,a+ \mu \textbf{u})\right]\right)\big|_{a = \pi_{\theta}(s)}\right\|\right]\\
	&\leq  B_\Theta\E_{s\sim \rho^{\pi_\theta}}\left[  \E_\mathbf{u} \left[\left\|\nabla_a \left( Q^{\pi_\theta}(s,a) -   Q^{\pi_\theta}(s,a+ \mu \textbf{u})\right)\big|_{a = \pi_{\theta}(s)}\right\|\right]\right]\\
	&\leq  B_\Theta\E_{s\sim \rho^{\pi_\theta}}\left[  \E_\mathbf{u} \left[G\left\|\mu \mathbf{u}\right\|\right]\right]\\
	&\leq B_\Theta \mu \sqrt{p} G/(1-\gamma),
	\end{split}
	\end{equation}
	where the final inequality comes from \cite[Lemma 1]{nesterov2017random}.

    If, alternatively, Assumption \ref{asm:nesterov_assumptions_H} instead of Assumption \ref{asm:nesterov_assumptions}, holds, then by \cite[Lemma 3]{nesterov2017random}, we obtain
	\begin{equation}
	\begin{split} \nonumber
	\| \nabla J(\theta) - \hat \nabla J(\theta)\| & \leq B_\Theta\E_{s\sim \rho^{\pi_\theta}}\left[ \left\| \nabla_a \left( Q^{\pi_\theta}(s,a) -  Q_\mu^{\pi_\theta}(s,a)\right)\big|_{a = \pi_{\theta}(s)}\right\|\right] \\
	&\leq B_\Theta \mu^2 (p+4)^2 H/(1-\gamma).
	\end{split}
	\end{equation}
    This completes the proof.\null\nobreak\hfill\ensuremath{\square}

\section{Proof of Lemma \ref{lem:bounded_variance}}

First, by definition, we may write
\begin{equation} \nonumber
(1-\gamma)^2\mathbb{E}[\Vert g_{t}\Vert^{2}]=\mathbb{E}\Bigg[\Bigg\Vert\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\dfrac{\hat{Q}_{t}^{+}-\hat{Q}_{t}^{-}}{\mu}{\bf u}\Bigg\Vert^{2}\Bigg]=:\mathbb{E}[\Vert g'_{t}\Vert^{2}],
\end{equation}
where
\begin{align}\nonumber
\hat{Q}_{t}^{\pm} & =Q^{\pi_{\theta_{t}}}(\pm)+\dfrac{1}{N}\sum_{n=1}^{N}w_{t,n}^{\pm}.
\end{align}
By Cauchy-Schwarz, and by exploiting the above, we readily get
\begin{flalign}
\mathbb{E}[\Vert g'_{t}\Vert^{2}] & \le B_{\Theta}^{2}\mathbb{E}\Bigg[\Bigg\Vert\dfrac{Q^{\pi_{\theta_{t}}}(+)-Q^{\pi_{\theta_{t}}}(-)}{\mu}+\dfrac{1}{N}\sum_{n=1}^{N}\dfrac{w_{t,n}^{+}-w_{t,n}^{-}}{\mu}\Bigg\Vert^{2}\Vert{\bf u}\Vert^{2}\Bigg]\nonumber \\
 & =B_{\Theta}^{2}\mathbb{E}\Bigg[\Bigg|\dfrac{Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t})+\mu{\bf u})-Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t}))}{\mu}+\dfrac{1}{N}\sum_{n=1}^{N}\dfrac{w_{t,n}^{+}-w_{t,n}^{-}}{\mu}\Bigg|^{2}\Vert{\bf u}\Vert^{2}\Bigg]\nonumber \\
 & \le2B_{\Theta}^{2}\mathbb{E}\Bigg[\Bigg|\dfrac{Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t})+\mu{\bf u})-Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t}))}{\mu}\Bigg|^{2}\Vert{\bf u}\Vert^{2} \nonumber \\
 &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\Bigg|\dfrac{1}{N}\sum_{n=1}^{N}\dfrac{w_{t,n}^{+}-w_{t,n}^{-}}{\mu}\Bigg|^{2}\Vert{\bf u}\Vert^{2}\Bigg].\label{eq:Bound_1}
\end{flalign}
For the first term on the right-hand side of (\ref{eq:Bound_1}),
we have
\begin{equation}
\mathbb{E}\Bigg[\Bigg|\dfrac{Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t})+\mu{\bf u})-Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t}))}{\mu}\Bigg|^{2}\Vert{\bf u}\Vert^{2}\Bigg]\le L^2\mathbb{E}\big[\Vert{\bf u}\Vert^{4}\big]\le L^2 (p+4)^{2}.\label{eq:Bound_2}
\end{equation}
For the second term on the right-hand side of (\ref{eq:Bound_1}),
we may write
\begin{align} \nonumber
\mathbb{E}\Bigg[\Bigg|\dfrac{1}{N}\sum_{n=1}^{N}\dfrac{w_{t,n}^{+}-w_{t,n}^{-}}{\mu}\Bigg|^{2}\Vert{\bf u}\Vert^{2}\Bigg] & =\dfrac{1}{\mu^{2}N^{2}}\mathbb{E}\Bigg[\mathbb{E}\Bigg[\Bigg|\sum_{n=1}^{N}w_{t,n}^{+}-w_{t,n}^{-}\Bigg|^{2}\Vert{\bf u}\Vert^{2}\Bigg|s_t,{\bf u},\theta_{t},T_{Q}\Bigg]\hspace*{-1bp}\hspace*{-1bp}\Bigg]\nonumber \\
 & =\dfrac{1}{\mu^{2}N^{2}}\mathbb{E}\Bigg[\Vert{\bf u}\Vert^{2}\mathbb{E}\Bigg[\Bigg|\sum_{n=1}^{N}w_{t,n}^{+}-w_{t,n}^{-}\Bigg|^{2}\Bigg|s_t,{\bf u},\theta_{t},T_{Q}\Bigg]\hspace*{-1bp}\hspace*{-1bp}\Bigg]. \nonumber
\end{align}
Naturally, we need to focus on the (inner) conditional expectation
above. By construction, it is true that, for every $t$, the differences
$w_{t,n}^{+}-w_{t,n}^{-},n\in\mathbb{N}_{N}^{+}$ are \textit{conditionally}
independent and identically distributed relative to $s,{\bf u},\theta_{t}$
and $T_{Q}$, and \textit{with zero mean}. This implies that
\begin{flalign}
\mathbb{E}\Bigg[\Bigg|\sum_{n=1}^{N}w_{t,n}^{+}-w_{t,n}^{-}\Bigg|^{2}\Bigg|s_t,{\bf u},\theta_{t},T_{Q}\Bigg] & =\sum_{n=1}^{N}\mathbb{E}\big[\big|w_{t,n}^{+}-w_{t,n}^{-}\big|^{2}\big|s_t,{\bf u},\theta_{t},T_{Q}\big]\nonumber \\
 & =N\mathbb{E}\big[\big|w_{t,1}^{+}-w_{t,1}^{-}\big|^{2}\big|s_t,{\bf u},\theta_{t},T_{Q}\big],\nonumber \\
 & \le N2\mathbb{E}\big[\big(w_{t,1}^{+}\big)^{2}+\big(w_{t,1}^{-}\big)^{2}\big|s_t,{\bf u},\theta_{t},T_{Q}\big]\nonumber \\
 & \le N4U_{R}^{2}\dfrac{1+\gamma}{(1-\gamma)^{2}}=:N\sigma^{2},\nonumber
\end{flalign}
where in the last line we have used Lemma \ref{lem:Bounded_Q_var} twice. As a result, we obtain that
\begin{equation}
\mathbb{E}\Bigg[\Bigg|\dfrac{1}{N}\sum_{n=1}^{N}\dfrac{w_{t,n}^{+}-w_{t,n}^{-}}{\mu}\Bigg|^{2}\Vert{\bf u}\Vert^{2}\Bigg]\le\dfrac{\sigma^{2}p}{\mu^{2}N},\label{eq:Bound_3}
\end{equation}
since $\mathbb{E}\big[\Vert{\bf u}\Vert^{2}\big]\equiv p$. Combining
(\ref{eq:Bound_2}) and (\ref{eq:Bound_3}), we end up with the uniform
bound
\begin{equation} \nonumber
\mathbb{E}[\Vert g_{t}\Vert^{2}]\le
\dfrac{2B_{\Theta}^{2}}{(1-\gamma)^2}
\Big(L^2(p+4)^{2}+\dfrac{\sigma^{2}p}{\mu^{2}N}\Big),
\end{equation}
thus completing the proof.
\null\nobreak\hfill\ensuremath{\square}

\section{Proof of Theorem \ref{thm:main_1}}

First, it is not hard to see that the stochastic ascend direction
$g_{t}$ constructed by Algorithm \ref{alg:ZOPG} is an unbiased estimate of the
quasi-gradient $\hat{\nabla}J$. Indeed, with $s_{t}\sim (1-\gamma)\rho^{\pi_{\theta_{t}}}$,
we may carefully write
\begin{flalign}
(1-\gamma)\mathbb{E}[g_{t}|\theta_{t}] & =\mathbb{E}\Bigg[\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\dfrac{\hat{Q}_{t}^{+}-\hat{Q}_{t}^{-}}{\mu}{\bf u}\Bigg|\theta_{t}\Bigg]\nonumber \\
 & =\mathbb{E}\Bigg[\mathbb{E}\Bigg[\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\dfrac{\hat{Q}_{t}^{+}(s_{t},{\bf u},\theta_{t},T_{Q})-\hat{Q}_{t}^{-}(s_{t},\theta_{t},T_{Q})}{\mu}{\bf u}\Bigg|s_{t},{\bf u},\theta_{t},T_{Q}\Bigg]\hspace*{-1bp}\hspace*{-1bp}\Bigg|\theta_{t}\Bigg]\nonumber \\
 & \equiv\mathbb{E}\Bigg[\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\mathbb{E}\Bigg[\dfrac{\hat{Q}_{t}^{+}(s_{t},{\bf u},\theta_{t},T_{Q})-\hat{Q}_{t}^{-}(s_{t},\theta_{t},T_{Q})}{\mu}\Bigg|s_{t},{\bf u},\theta_{t},T_{Q}\Bigg]\hspace*{-1bp}{\bf u}\Bigg|\theta_{t}\Bigg]\nonumber \\
 & =\mathbb{E}\Bigg[\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\dfrac{Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t})+\mu{\bf u})-Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t}))}{\mu}{\bf u}\Bigg|\theta_{t}\Bigg]\nonumber \\
 & =\mathbb{E}\Bigg[\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\mathbb{E}\Bigg[\dfrac{Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t})+\mu{\bf u})-Q^{\pi_{\theta_{t}}}(s_{t},\pi_{\theta_{t}}(s_{t}))}{\mu}{\bf u}\Bigg|s_{t},\theta_{t}\Bigg]\hspace*{-1bp}\hspace*{-1bp}\Bigg|\theta_{t}\Bigg]\nonumber \\
 & =\mathbb{E}\big[\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\nabla_{a}Q_{\mu}^{\pi_{\theta_{t}}}(s,a)|_{a=\pi_{\theta_{t}}(s_{t})}|\theta_{t}\big]\nonumber \\
 & =
 (1-\gamma)
 \mathbb{E}_{s_{t}\sim\rho^{\pi_{\theta_{t}}}}\big[\nabla_{\theta}\pi_{\theta_{t}}(s_{t})\nabla_{a}Q_{\mu}^{\pi_{\theta_{t}}}(s,a)|_{a=\pi_{\theta_{t}}(s_{t})}|\theta_{t}\big] \nonumber \\
 & =
 (1-\gamma)
 \hat{\nabla}J(\theta_t). \nonumber
\end{flalign}

Next, using the assumption that $J(\theta)$ is $L_J$-smooth, we have

\begin{equation} \nonumber
    J(\theta_{t + 1}) \geq J(\theta_t) + (\theta_{t+1} - \theta_t)^\top \nabla J(\theta_t) - L_J\|\theta_{t+1} - \theta_t\|^2
\end{equation}

From Algorithm \ref{alg:ZOPG}, we have $\theta_{t+1} - \theta_t = \alpha g_t$. Let $\E_t$ denote conditional expectation taken relative to history up to iteration $t-1$. Then, by taking the expectation on both sides, we obtain

\begin{equation} \nonumber
    \E_t\left[J(\theta_{t+1})\right] \geq J(\theta_t) + \alpha \hat \nabla J(\theta_t)^\top \nabla J(\theta_t) - L_J\alpha^2 \E_t\left[ \| g_t  \|^2 \right].

\end{equation}

Now add and subtract $\alpha  \nabla J(\theta_t)^\top \nabla J(\theta_t)$ to the right hand side. By grouping terms, we may write
\begin{equation} \nonumber
    \E_t  \left[J(\theta_{t+1})\right] \ge
    J(\theta_t) + \alpha  \big( \hat \nabla J(\theta_t) - \nabla J(\theta_t)\big)^{\hspace{-1pt}\top} \nabla J(\theta_t) + \alpha  \|\nabla J(\theta_t)\|^2 - L_J\alpha^2\E_t\left[ \| g_t  \|^2 \right].
\end{equation}

Using Cauchy-Schwarz, we bound the term $\big( \hat \nabla J(\theta_t) - \nabla J(\theta_t)\big)^\top \nabla J(\theta_t)$ by $-\| \hat \nabla J(\theta_t) - \nabla J(\theta_t)\|\cdot \|\nabla J(\theta_t)\|$. Further, because $J(\theta)$ is assumed to be $G_J$-Lipschitz, we have

\begin{equation} \nonumber
\E_t \left[J(\theta_{t+1})\right] \geq J(\theta_t) - \alpha  G_J \|\hat \nabla J(\theta_t) - \nabla  J(\theta_t)\| + \alpha  \|\nabla J(\theta_t)\|^2 - L_J\alpha^2 \E_t\left[ \| g_t  \|^2 \right].
\end{equation}

Taking the total expectation, invoking Lemma \ref{lem:bounded_variance}, summing up until $T$, and dividing by $1/(\alpha T)$, we obtain
\vspace{-11pt}
\begin{equation}
\begin{split}\label{equ:before_rearranging}
 \frac{1}{T} \sum_{t = 1}^T \E  \|\nabla J(\theta_t)\|^2 \leq & \hspace{1pt} \frac{\Delta_0}{\alpha T}
 +  L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}(p+4)^2L^2
 + L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}\frac{p \sigma^2}{\mu^2N}  \\
 &+ \frac{ G_J}{T} \sum_{t = 1}^T \E [\|\hat \nabla J(\theta_t) - \nabla J(\theta_t)\|],
 \end{split}
\end{equation}
where $\Delta_0 =  J(\theta^*)-J(\theta_0) $ and $\theta^* \in \arg\max J(\theta)$.

Under Assumption \ref{asm:nesterov_assumptions}, Theorem \ref{thm:nesterov_result1} implies that

\begin{equation}
\frac{1}{T} \sum_{t = 1}^T \E  \|\nabla J(\theta_t)\|^2 \leq \frac{\Delta_0}{\alpha T}

+  L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}(p+4)^2L^2
 + L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}\frac{p \sigma^2}{\mu^2N}
+ G_J\frac{B_\Theta \mu \sqrt{p}G}{1-\gamma}. \nonumber
\end{equation}

If, alternatively, Assumption \ref{asm:nesterov_assumptions_H} holds, then from \eqref{equ:before_rearranging} we use the pertinent result form Theorem \ref{thm:nesterov_result1} to obtain

\begin{equation}
\frac{1}{T} \sum_{t = 1}^T \E  \|\nabla J(\theta_t)\|^2 \leq \frac{\Delta_0}{\alpha T}

+  L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}(p+4)^2L^2
 + L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}\frac{p \sigma^2}{\mu^2N}
+ G_J \frac{B_\Theta \mu^2 (p+4)^2H}{1-\gamma}. \nonumber
\end{equation}

To complete the proof, we set the number of Monte Carlo rollouts ($N$), smoothing parameter ($\mu$), and step size ($\alpha$) in
each of the expressions above, according to Table \ref{sample-table}. Enough said.
\null\nobreak\hfill\ensuremath{\square}

\textbf{Extensions for constrained policy parameterizations:} In such a case, we may consider a projected gradient step in Algorithm \ref{alg:ZOPG}, i.e.,
\begin{equation}
\theta_{t+1}=\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\},\nonumber
\end{equation}
where $\Theta\subset\mathbb{R}^{d}$ is also assumed to be closed
and convex. Below we present two approaches for deriving meaningful convergence rates.

\textit{Noisy projection approach:} In this approach, we may expand the product $(\theta_{t+1}-\theta_{t})^{\top}\nabla J(\theta_{t})$
into two relevant terms as
\begin{align*}
(\theta_{t+1}-\theta_{t})^{\top}\nabla J(\theta_{t}) & =(\theta_{t+1}+\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\theta_{t})^{\top}\nabla J(\theta_{t})\\
 & =(\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\theta_{t})^{\top}\nabla J(\theta_{t})\\
 &\quad\quad+(\theta_{t+1}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\})^{\top}\nabla J(\theta_{t}).
\end{align*}
For the first term, we invoke the key property \cite[Lemma 1]{ghadimi2016mini}
\begin{equation}
(\Pi_{\Theta}\{\theta+\alpha g\}-\theta)^{\top}g\ge\dfrac{1}{\alpha}\Vert\Pi_{\Theta}\{\theta+\alpha g\}-\theta\Vert^{2},\quad\forall(\theta,g)\in\Theta\times\mathbb{R}^{d}, \label{GenProj}
\end{equation}
yielding
\[
(\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\theta_{t})^{\top}\nabla J(\theta_{t})\ge\dfrac{1}{\alpha}\Vert\theta_{t}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}\Vert^{2}.
\]
For the second term, we may further expand as
\begin{flalign*}
 & \hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}(\theta_{t+1}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\})^{\top}\nabla J(\theta_{t})\\
 & =(\theta_{t+1}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\})^{\top}\nabla J(\theta_{t})\\
 & \quad\quad+(\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\})^{\top}\nabla J(\theta_{t}),
\end{flalign*}
which implies that
\begin{flalign*}
 & \hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\mathbb{E}_{t}\{(\theta_{t+1}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\})^{\top}\nabla J(\theta_{t})\}\\
 & =(\mathbb{E}_{t}\{\theta_{t+1}\}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\})^{\top}\nabla J(\theta_{t})\\
 & \quad\quad+(\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\})^{\top}\nabla J(\theta_{t})\\
 & \ge-\Vert\nabla J(\theta_{t})\Vert\Vert\mathbb{E}_{t}\{\theta_{t+1}\}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\Vert\\
 & \quad\quad-\alpha\Vert\nabla J(\theta_{t})\Vert\Vert\hat{\nabla}J(\theta_{t})-\nabla J(\theta_{t})\Vert\\
 & \ge-G_{J}\Vert\mathbb{E}_{t}\{\theta_{t+1}\}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\Vert\\
 & \quad\quad-\alpha G_{J}\Vert\hat{\nabla}J(\theta_{t})-\nabla J(\theta_{t})\Vert.
\end{flalign*}
Consequently, we obtain the lower bound
\begin{flalign*}
\mathbb{E}_{t}\{(\theta_{t+1}-\theta_{t})^{\top}\nabla J(\theta_{t})\} & \ge\dfrac{1}{\alpha}\Vert\theta_{t}-\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}\Vert^{2}\\
 & \quad\;\,\,-\alpha G_{J}\Vert\hat{\nabla}J(\theta_{t})-\nabla J(\theta_{t})\Vert\\
 & \quad\quad-G_{J}\Vert\mathbb{E}_{t}\{\theta_{t+1}\}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\Vert.
\end{flalign*}
Exploiting smoothness as in the unconstrained case discussed above,
we end up with the complexity estimate
\begin{flalign*}
 & \hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\dfrac{1}{T}\sum_{t=1}^{T}\mathbb{E}\dfrac{1}{\alpha^{2}}\Vert \Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\} - \theta_{t}\Vert^{2}\\
 & \le\dfrac{\Delta_{0}}{\alpha T}

 +  L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}(p+4)^2L^2
 + L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}\frac{p \sigma^2}{\mu^2N}
 \\
 & \quad+\dfrac{G_{J}}{T}\sum_{t=1}^{T}\mathbb{E}\big[\Vert\hat{\nabla}J(\theta_{t})-\nabla J(\theta_{t})\Vert\big]
 \\ &  \quad\quad
 +\dfrac{G_{J}}{\alpha T}\sum_{t=1}^{T}\mathbb{E}\big[\Vert\mathbb{E}_{t}[\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}]-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}\Vert\big].
\end{flalign*}
Note that, in the above, whenever it is actually the case that $\theta_{t}+\alpha g_{t}\in\Theta$
\textit{conditioned on} $\theta_{t}$, it follows that (note that $\Theta$
is closed convex)
\begin{flalign*}
\Theta\ni\mathbb{E}_{t}[\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}] & =\mathbb{E}_{t}[\theta_{t}+\alpha g_{t}]\\
 & =\theta_{t}+\alpha\hat{\nabla}J(\theta_{t}),
\end{flalign*}
which implies that
\[
\mathbb{E}_{t}[\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}]-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}=0.
\]
This will happen when $g_{t}$ is (or can be made) uniformly not too
large given $\theta_{t}$, and at the same time the stepsize
$\alpha$ is sufficiently small. Exactly such a behavior has been
observed in our simulations.

\textit{Mini-batch approach:} Alternatively, we may proceed to derive a slightly more pessimistic complexity estimate by writing

\begin{align*}
(\theta_{t+1}-\theta_{t})^{\top}\nabla J(\theta_{t}) & =(\theta_{t+1}-\theta_{t})^{\top}(\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})+\hat{\nabla}J(\theta_{t}))\\
 & =(\theta_{t+1}-\theta_{t})^{\top}(\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t}))+(\theta_{t+1}-\theta_{t})^{\top}\hat{\nabla}J(\theta_{t}).
\end{align*}
For the first term, Cauchy-Schwarz trivially gives
\[
(\theta_{t+1}-\theta_{t})^{\top}(\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t}))\ge-\alpha\Vert g_{t}\Vert\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert.
\]
For the second term, we may further expand as
\begin{flalign*}
(\theta_{t+1}-\theta_{t})^{\top}\hat{\nabla}J(\theta_{t}) & =(\theta_{t+1}-\theta_{t})^{\top}(\hat{\nabla}J(\theta_{t})+g_{t}-g_{t})\\
 & =(\theta_{t+1}-\theta_{t})^{\top}(\hat{\nabla}J(\theta_{t})-g_{t})+(\theta_{t+1}-\theta_{t})^{\top}g_{t},
\end{flalign*}
where again invoking \eqref{GenProj} yields
\[
(\theta_{t+1}-\theta_{t})^{\top}g_{t}\ge\dfrac{1}{\alpha}\Vert\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}-\theta_{t}\Vert^{2},
\]
whereas it also holds that
\begin{flalign*}
(\theta_{t+1}-\theta_{t})^{\top}(\hat{\nabla}J(\theta_{t})-g_{t}) & =(\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}-\theta_{t})^{\top}(\hat{\nabla}J(\theta_{t})-g_{t})\\
 & =(\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\})^{\top}(\hat{\nabla}J(\theta_{t})-g_{t})\\
 & \quad\;\,\,+(\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}-\theta_{t})^{\top}(\hat{\nabla}J(\theta_{t})-g_{t}).
\end{flalign*}
In regard, to the latter expression, it is then true that
\begin{flalign*}
\mathbb{E}_{t}[(\theta_{t+1}-\theta_{t})^{\top}(\hat{\nabla}J(\theta_{t})-g_{t})] & =\mathbb{E}_{t}[(\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\})^{\top}(\hat{\nabla}J(\theta_{t})-g_{t})]\\
 & \quad\quad+0\\
 & \ge-\alpha\mathbb{E}_{t}[\Vert g_{t}-\hat{\nabla}J(\theta_{t})\Vert^{2}].
\end{flalign*}
Combining everything, we get
\begin{align*}
\mathbb{E}_{t}[(\theta_{t+1}-\theta_{t})^{\top}\nabla J(\theta_{t})] & \ge\dfrac{1}{\alpha}\Vert\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}-\theta\Vert^{2}\\
 & \quad\;\,\,-\alpha\mathbb{E}_{t}[\Vert g_{t}\Vert]\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert\\
 & \quad\quad\;\,\,-\alpha\mathbb{E}_{t}[\Vert g_{t}-\hat{\nabla}J(\theta_{t})\Vert^{2}].
\end{align*}
Now, note that that we can write
\begin{align*}
 & \hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\dfrac{1}{\alpha^{2}}\mathbb{E}[\Vert\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\theta_{t}\Vert^{2}]\\
 & =\dfrac{1}{\alpha^{2}}\mathbb{E}[\Vert\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}\\
 & \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}-\theta_{t}\Vert^{2}]\\
 & \le2\mathbb{E}[\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert^{2}]+2\dfrac{1}{\alpha^{2}}\mathbb{E}[\Vert\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}-\theta_{t}\Vert^{2}]\\
 & =2\mathbb{E}[\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert^{2}]\\
 & \quad\;\,\,+2\dfrac{1}{\alpha^{2}}\mathbb{E}[\Vert\Pi_{\Theta}\{\theta_{t}+\alpha\hat{\nabla}J(\theta_{t})\}-\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}\\
 & \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\,\,+\Pi_{\Theta}\{\theta_{t}+\alpha g_{t})\}-\theta_{t}\Vert^{2}]\\
 & \le2\mathbb{E}[\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert^{2}]\\
 & \quad\;\,\,+4\mathbb{E}[\Vert g_{t}-\hat{\nabla}J(\theta_{t})\Vert^{2}]\\
 & \quad\quad\;\,\,+4\dfrac{1}{\alpha^{2}}\mathbb{E}[\Vert\Pi_{\Theta}\{\theta_{t}+\alpha g_{t}\}-\theta_{t}\Vert^{2}].
\end{align*}
Finally, following the usual steps as above, we obtain the complexity
estimate
\begin{flalign*}
 & \hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\dfrac{1}{T}\sum_{t=1}^{T}\mathbb{E}\dfrac{1}{\alpha^{2}}\Vert\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\theta_{t}\Vert^{2}\\
 & \le\dfrac{2}{T}\sum_{t=1}^{T}\mathbb{E}[\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert^{2}]+\dfrac{4}{T}\sum_{t=1}^{T}\mathbb{E}[\Vert g_{t}-\hat{\nabla}J(\theta_{t})\Vert^{2}]\\
 & \quad + \dfrac{\Delta_{0}}{\alpha T}

 +  L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}(p+4)^2L^2
 + L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}\frac{p \sigma^2}{\mu^2N}
 \\
 & \quad\quad+\dfrac{1}{T}\sum_{t=1}^{T}\mathbb{E}[\Vert g_{t}\Vert\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert]+\dfrac{1}{T}\sum_{t=1}^{T}\mathbb{E}[\Vert g_{t}-\hat{\nabla}J(\theta_{t})\Vert^{2}],
\end{flalign*}
or, more compactly,
\begin{flalign*}
 & \hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\hspace*{-1bp}\dfrac{1}{T}\sum_{t=1}^{T}\mathbb{E}\dfrac{1}{\alpha^{2}}\Vert\Pi_{\Theta}\{\theta_{t}+\alpha\nabla J(\theta_{t})\}-\theta_{t}\Vert^{2}\\
 & \le\dfrac{\Delta_{0}}{\alpha T}

 +  L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}(p+4)^2L^2
 + L_J\alpha \frac{2B_\Theta^2}{(1-\gamma)^2}\frac{p \sigma^2}{\mu^2N}
 \\
 & \quad+\dfrac{1}{T}\sum_{t=1}^{T}\mathbb{E}[\Vert g_{t}\Vert\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert]+\dfrac{2}{T}\sum_{t=1}^{T}\mathbb{E}[\Vert\nabla J(\theta_{t})-\hat{\nabla}J(\theta_{t})\Vert^{2}]\\
 & \quad\quad+\dfrac{5}{T}\sum_{t=1}^{T}\mathbb{E}[\Vert g_{t}-\hat{\nabla}J(\theta_{t})\Vert^{2}],
\end{flalign*}

This bound suggests that further variance reduction can be employed
in order to make the term involving $\mathbb{E}[\Vert g_{t}-\hat{\nabla}J(\theta_{t})\Vert^{2}]=\mathbb{E}[\Vert g_{t}-\mathbb{E}[g_{t}|\theta_{t}]\Vert^2]$
small, for instance by exploiting mini-batching, similar in spirit
to \cite{ghadimi2016mini}. However, this comes at the expense of significantly
higher training overhead.

\section{Complementary Numerical Results} \label{supp:results}

\begin{figure*}[!t]
	\centering
	\begin{tabular}{cc}
		\hspace{5.2mm}\includegraphics[trim = {5cm, 10cm, 5cm, 10cm}, width=.415\linewidth]{small_noisy_10_corrected}   &\hspace{8mm}
		\includegraphics[trim = {5cm, 10cm, 5cm, 10cm},width=.415\linewidth]{noisy_10_corrected} \\ \vspace{.5cm}
		~&~\\
		\hspace{5.2mm}\includegraphics[trim = {5cm, 10cm, 5cm, 10cm}, width=.415\linewidth]{small_noisy_1_corrected}   &\hspace{8mm}
		\includegraphics[trim = {5cm, 10cm, 5cm, 10cm},width=.415\linewidth]{noisy_1_corrected} \\ \vspace{.4cm}
		~&~\\
	\end{tabular}
	\caption{Average reward per episode of length $T = 20$ with confidence bounds over 50 trials.  }\vspace{4mm}
	\label{fig:LC_20}
	\vspace{10pt}
\end{figure*}

\begin{figure*}[!t]
	\centering
	\begin{tabular}{cc}
		\hspace{5.2mm}\includegraphics[trim = {5cm, 10cm, 5cm, 10cm}, width=.415\linewidth]{small_noisy_eval_100}   &\hspace{8mm}
		\includegraphics[trim = {5cm, 10cm, 5cm, 10cm},width=.415\linewidth]{noisy_eval_100} \\ \vspace{.5cm}
		~&~\\
	\end{tabular}
	\caption{Average reward per episode of length $T = 100$ with confidence bounds over 50 trials.}\vspace{2mm}
	\label{fig:LC_100}
	\vspace{-10pt}
\end{figure*}

\begin{figure*}[!t]
	\centering
	\begin{tabular}{ccc}
		\hspace{3mm}\includegraphics[trim = {5cm, 10cm, 5cm, 10cm}, width=.274\linewidth]{True_baseline20}   &\hspace{3.5mm}
		\includegraphics[trim = {5cm, 10cm, 5cm, 10cm},width=.274\linewidth]{True_baseline200} &\hspace{3.5mm}
		\includegraphics[trim = {5cm, 10cm, 5cm, 10cm},width=.274\linewidth]{True_baseline600} \\ \vspace{.4cm}
		~&~&~\\
		\small (a) & \small (b) & \small (c)
	\end{tabular}
	\caption{Average reward per episode of length $T = 20$ with confidence bounds over 50 trials. System noise $\omega = 0.01$ with no MC variance reduction ($N = 1$). Variance reduction is applied \textit{only} on value-function estimates to approximate the true (deterministic) baseline exclusively for PG-B, with (a) $20$, (b) $200$, and (c) $600$ rollouts.}
	\label{fig:true_baseline}
\end{figure*}
In section \ref{sec:Numericals}, we presented the learning curves for our algorithm implemented on a navigation problem with and without system noise. We showed that either with or without system noise, ZDPG and ZDPG-S enjoy a substantial advantage over PG with baseline (PG-B), and of course standard PG, especially without MC variance reduction, which is of importance in various practical settings.

The main purpose of this section is to confirm that the advantage of ZDPG does not disappear as the noise in the system further increases; in fact, it becomes even more drastic. We show this by comparing the learning curves on our navigation example for system noise variance $\omega = 0.01$ and $\omega = 0.1$, for different values of Monte-Carlo variance reduction $(N)$, and policy evaluation rollout length $(T)$. In this set of simulations, we only compare ZDPG-S with PG-B, since these methods seem to be empirically superior to ZDPG and PG, respectively.

In Figure \ref{fig:LC_20}, the evaluation rollout length is set to $T = 20$. The smoothing parameter for ZDPG-S is $\mu  = 0.5$ for the $\omega = 0.01$ system noise case and $\mu = 0.6$ for the $\omega = 0.1$ system noise case. Analogously, for PG-B, we select action variance as $0.025$ for both levels of system noise as this value showed the most consistent performance throughout this line of experiments (by trial-and-error). For all simulations and for both methods, the learning rate was set to $\alpha = 10^{-7}$. Figure \ref{fig:LC_20} (top) shows that, indeed, ZDPG-S remains advantageous over and consistently outperforms PG-B when the system noise is increased. More interestingly, ZDPG-S shows significant advantage over PG-B in terms of noise resilience especially when there is no MC variance reduction at all (i.e. $N = 1$), which, as also stated above, is a case of interest in most practical scenarios.

Figure \ref{fig:LC_100} shows the learning curves for both methods when the evaluation length is set to $T = 100$. While the performance of both ZDPG-S and PG-B is more sensitive over such a longer evaluation horizon, the plots show that ZDPG-S experiences much more stability  in solving the navigation problem for both levels of system noise. In contrast, PG-B becomes significantly unstable, resulting in trajectories that either do not necessarily reach the target (see Figure \ref{fig:traj}), or reach the target after a phase of prolonged, suboptimal wandering.

Because PG-B relies on the estimate ($\hat Q(s,a) - \tilde{V}(s)$),

MC variance reduction (whenever $N>1$) is essentially applied to both parts of the involved two-point difference for the results presented in Figure \ref{fig:LC_20} (and analogously for the two-point gradient estimates in ZDPG). In that way, the comparison between ZDPG and PG-B is fair for every $N\ge1$. Still, it would be interesting to compare ZDPG-S with PG-B, when in the latter

a more accurate version of the corresponding true value function is employed, resembling a standard (non-stochastic) PG with baseline. To do this, we conducted the same experiment as in Figure \ref{fig:LC_20} for the low system noise case ($\omega = 0.01$) and with no MC reduction ($N=1$), but where only the baseline $\tilde{V}(s)$ is constructed as a \textit{multiple-rollout} estimate of the associated value function at state $s$ (i.e., an average).

We note, though, that this procedure results in a PG-B algorithm which is much more wasteful in terms of resources, whereas at the same time its comparison with ZDPG is extremely unfair, in favor of ZDPG.
This resulting plots are shown in Figure \ref{fig:true_baseline}. It very interesting to see that, even though the performance of PG-B slightly improves with the consideration of multiple-rollouts for the corresponding value function estimate, it still performs rather poorly compared with ZDPG-S, which of course demands just two system rollouts per episode.

Again, we would like to reiterate that the success of ZDPG comes from the non-heuristic use of deterministic policies, in contrast to random policies essential in log-trick-based standard (baseline) PG methods. At least with regards to our simulation setting, removing this level of randomness from the evolution of the system during training allows the agent to reach a solution quickly and consistently.


TITLE: Showing a certain series equals a certain trig function in complex analysis
QUESTION [0 upvotes]: I was reading a a complex analysis proof that was showing that for a fixed $\alpha$ s.t. $Im(\alpha)>0$ we have
$\sum \limits_{m=-\infty}^{\infty} \frac{1}{(\alpha +n)^2}=(-4\pi^2)\sum \limits_{m=1}^{\infty} m e^{2\pi i m\alpha}$.
I understood this proof.  However, I need help with two other related statements.  First, the proof notes at the end this means that the series above on the left equals $\frac{\pi^2}{\sin^2(\pi\alpha)}$.  In other words I am suppose to see that the one on the right above is equal to $\frac{\pi^2}{\sin^2(\pi\alpha)}$.  I don't know the correct identities to show this so I would appreciate if someone would explain that.  Also, after the proof there is sort of a food for thought question that asks the reader if this identity was true for all complex numbers other than integers.  This seems valid but I would appreciate if someone would let me know that as well.
Thanks!

REPLY [1 votes]: With $\operatorname{Im} \alpha > 0$, we have $\lvert e^{2\pi i\alpha}\rvert = e^{-2\pi \operatorname{Im}\alpha} < 1$, so the series is basically a differentiated geometric series,
$$\sum_{m=1}^\infty m\cdot x^m = x\sum_{m=1}^\infty mx^{m-1} = x\frac{d}{dx}\left(\frac{1}{1-x}\right) = \frac{x}{(1-x)^2}.$$
Inserting $x = e^{2\pi i\alpha}$, we obtain
$$-4\pi^2 \sum_{m=1}^\infty me^{2\pi im\alpha} = \frac{-4\pi^2e^{2\pi i\alpha}}{(1- e^{2\pi i\alpha})^2},$$
moving the $e^{2\pi i\alpha}$ to the denominator, we then obtain
$$\sum_{n=-\infty}^\infty \frac{1}{(\alpha+n)^2} =\frac{-4\pi^2}{(e^{\pi i \alpha} - e^{-\pi i\alpha})^2} = \frac{\pi^2}{\sin^2 \pi\alpha}.$$
Since the sum on the left hand side converges locally uniformly for $\alpha \notin \mathbb{Z}$, it defines a holomorphic function there. By the identity theorem, that holomorphic function is $\frac{\pi^2}{\sin^2 \pi\alpha}$ on all of $\mathbb{C}\setminus\mathbb{Z}$.


TITLE: Relations between $M_k(\Gamma_0(N))$ and $M_k(\Gamma_0(L,M))$, and Eisenstein Series.
QUESTION [2 upvotes]: Consider the following congruence subgroups
$$\Gamma_0(N):=\Big\{ \begin{pmatrix} a & b \\ c & d \end{pmatrix}\in\text{SL}_2(\mathbb{Z})\; \Big| \; c\equiv0\text{ (mod }N)\Big\}$$
and
$$\Gamma_0(L,M):=\Big\{ \begin{pmatrix} a & b \\ c & d \end{pmatrix}\in\text{SL}_2(\mathbb{Z})\; \Big| \; c\equiv0\text{ (mod }L)\;\text{ and }
  \;b\equiv0\text{ (mod }M) \Big\}.$$
Is there some sort of relationship between modular forms over these groups? Say, for example, if $N=LM$ or $N=\text{lcm}(N,M)$?
I am asking this because I was reading Miyake's book (titled "Modular Forms"), at chapter 7, for two characters $\chi$ and $\psi$ of modulus $L$ and $M$, he defines the Eisenstein Series $E_k(z,s,\chi,\psi)$. He then proves (for $s=0$ and under some conditions), that they are modular over $\Gamma_0(L,M)$ with nebentypus $\chi\bar{\psi}$.
However, everywhere else that I see these Eisensteins series they are understood as being modular over $\Gamma_0(LM)$. For example https://wstein.org/books/modform/modform/eisenstein.html (see the comments after equation (4)).
I am confused. Are these the same series, or do they require some slight modification to be modular over $\Gamma_0(LM)$?

REPLY [3 votes]: It's not true that there is a convenient relationship between $\Gamma_0(L, M)$ and $\Gamma_0(N)$. But there is a relationship between the Eisenstein series defined on them. The Eisenstein series in Stein and the Eisenstein series in Miyake aren't quite the same Eisenstein series.
Stein's Eisenstein series appears earlier in Miyake's book. It's the one described in Theorem 4.7.1 (pg 177) of Miyake, which is shown to be modular on $\Gamma_0(LM)$, where $L$ and $M$ are the conductors of the Dirichlet characters (where I use the notation from chapter 7 and your question, which differs from the notation in the book in chapter 4).
Below Theorem 7.1.3 of Miyake, in equation 7.1.13 (pg 271), Miyake states that one gets the Stein-type Eisenstein series via
$$ E_k(Mz; \chi, \psi) = c E_{\mathrm{Stein}}(z; \chi, \overline{\psi}) $$
for a constant $c = c(k, \psi, M)$. Here, $E_{\mathrm{Stein}}$ is my notation for the Stein-type Eisenstein series.
Thus the conversion is to consider $E_k(Mz; \chi, \psi)$, which is modular on $\Gamma_0(LM)$. Stein translates results from Miyake through this identity for his book.
Additional note: a typo in Stein
There is a typo that make this comparison more confusing. Equation 4 in https://wstein.org/books/modform/modform/eisenstein.html (or equivalently equation 5.3.1 in the book version), $\psi(n)$ should be replaced by $\overline{\psi}(n)$. And in Theorem 5.9, the condition should either read
$\chi = \varepsilon \psi$ or $\chi \overline{\psi} = \varepsilon$.
In fact, this was the source of a problem in sagemath for a bit. Here is the source of that trac ticket


TITLE: Calculus implicit differential
QUESTION [0 upvotes]: So i am here with my question , Question asks $(x-y)^2 = x + y - 1$.
$2 (x-y) (1- y') = 1 + y'$
I am a little confused when it comes to right here, can anyone clarify a little,
I know the step that should come next which is $2 (x-y) - 2 (x-y) y' = 1 + y'$
but how did they get $-2(x-y)y'$?

REPLY [1 votes]: What the authors probably do is suppose that $y$ is a function of $x$. That way, you have that $$(x+y-1)' = (x)' + (y(x))' - (1)' = 1 + y' - 0$$
on the right side. The left side is best seen as a function of $x$ that maps $x$ to $(x-y(x))^2$. Differenting this gives $$2(x-y(x)) \cdot (x-y(x))'= 2(x-y(x))(1-y')$$


text:
\begin{document}

\title[\it Stationary states]{Stationary states of two-level open quantum systems}

\author{Bart{\l}omiej Gardas}
\address{Institute of Physics, University of Silesia, PL-40-007 Katowice, Poland}
\ead{bartek.gardas@gmail.com}

\author{Zbigniew Pucha{\l}a}
\address{Institute of Theoretical and Applied Informatics, Polish Academy
of Sciences, Ba{\l}tycka 5, 44-100 Gliwice, Poland}
\ead{z.puchala@iitis.pl}

\begin{abstract}
A problem of finding stationary states of  open  quantum systems is addressed.
We focus our attention on a generic type of open system: a qubit coupled to its
environment. We apply the theory of block operator matrices and find  stationary
states of two--level open  quantum systems under certain conditions applied both
on the qubit and the surrounding.
\end{abstract}

\pacs{03.65.Yz, 03.67.-a, 02.30.Tb, 03.65.Db}

\maketitle

\section{Introduction}

In quantum mechanics, the density operator $\rho$ of a quantum system is called a stationary state if
$[H,\rho]=0$, where $H$ is a given time-independent Hamiltonian of the system. Since $\rho$ satisfies
the Liouville-von Neumann equation $i\partial_t\rho=[H,\rho]$, it is clear that stationary states are invariant
with respect to the transformation $\rho \mapsto U_t\rho U_t^{\dagger}$, where $U_t=\exp(-\rmi Ht)$ is the
time evolution operator. In other words, stationary states do not change during the time evolution.

For low-dimensional closed systems, the stationary states can be obtained relatively easily~\cite{sakurai,galindo1}.
It is a common situation, that a small quantum system is immersed in other, mostly large, system called the
environment~\cite{Leggett}. Such an open system does not evolve unitarily in time. An analysis of open quantum
systems~\cite{davies,alicki,breuer} is much more complicated as they are a stage of a variety of physical
phenomena~\cite{dajka2,dajka4,dajka5}. The famous decoherence process~\cite{zurek} may serve as an example. In open
quantum systems character of potentially existing stationary states is not obvious.

There are various physical problems related to the properties of open quantum systems, which has already been addressed
and intensively discussed (see \eg, ~\cite{dajka1,dajka3}). Nevertheless, the procedure of deriving stationary states is
not one of them. The existence and properties of stationary states have significant importance in quantum information
processing and quantum theory itself. One can pose natural questions:

\begin{enumerate}
 \item do the stationary states exist for a given open system?
 \item what features of a given model are responsible for existence of such states?
 \item how such states can be constructed?
\end{enumerate}

The answers to the above questions are still incomplete. For example, it is  known that the stationary states
exist for completely positive (CP)~\cite{cp} evolution of the open system, this fact follows directly from
Schauder fixed point theorem~\cite{Chuang}. However, this is only an existential result and so far there are
no available methods to determine explicit form of stationary states. Furthermore, the very existence of the
stationary states in general case is an open problem \eg, in the presence of initial system-environment correlations~\cite{NCP}.

The purpose of this paper is to propose the method of calculating the stationary states in the case of two-dimensional
open quantum system. The theory of block operator matrices~\cite{Vadim,bom,spectral} is adapted to achieve this goal.
In particular, we use the \emph{Riccati operator equation}~\cite{RiccEq} to solve the eigenproblem for the total Hamiltonian.
It is shown how to derive the stationary states by using the solution of the equation.

\section{Block operator matrix approach}

We begin with a brief review of the block operator matrices approach to the problem of decoherence in the case of a single
qubit~\cite{mgr,gardas,gardas2}. Let $\mf{H}$ be the Hamiltonian of the total system. We will assume that it has the following form

\begin{equation}
 \label{total}
\mf{H} = H_{\text{Q}}\ot\idE+\idQ\ot H_{\text{E}} + \mf{H}_{\text{int}},
\end{equation}

where $H_{\text{Q}}$ and $H_{\text{E}}$ represent the Hamiltonian of the qubit and the environment, respectively,
while $\mf{H}_{\text{int}}$ specifies the interaction between the systems. The Hamiltonian $\mf{H}$ acts on
the Hilbert space $\mc{H}_{\text{tot}}=\mb{C}^{2}\ot\mc{H}_{\text{E}}$, where $\mc{H}_{\text{E}}$ is the Hilbert
space (possibly infinite-dimensional) related to the environment. $\idQ$ and $\idE$ are the identity operators on
$\mb{C}^2$ and $\mc{H}_{\text{E}}$, respectively.

Since the isomorphism $\mb{C}^2\ot\mc{H}_{\text{E}}\simeq\mc{H}_{\text{E}}\oplus\mc{H}_{\text{E}}$ holds true,
the Hamiltonian~(\ref{total}) admits the block operator matrix representation~\cite{spectral}:

\begin{equation}
 \label{bom}
\mf{H} =
  \begin{bmatrix}
  H_{\p} & V \\
  V^{\dagger} & H_{\m}
  \end{bmatrix}
  \quad\text{on}\quad
  \mathcal{D}(\mf{H})=\left(\mathcal{D}(H_{\p})\cap\mathcal{D}(V^{\dagger})\right)\oplus
\left(\mathcal{D}(V)\cap\mathcal{D}(H_{\m})\right).
\end{equation}

All the entries of~(\ref{bom}) are operators acting on $\mc{H}_{\text{E}}$. Moreover, the diagonal entries, \ie, $H_{\spm}$
are self-adjoint. In this paper, we will focus on the case in which $V$ is bounded, thus $V^{\dagger}$ is bounded as well;
however, no assumption on boundedness of $H_{\spm}$ is made. Under these circumstances we have
$\mathcal{D}(\mf{H})=\mathcal{D}(H_{\p})\oplus\mathcal{D}(H_{\m})$, where domains $\mc{D}(H_{\spm})$ are assumed to be dense in $\mc{H}_{\text{E}}$.

The generally accepted procedure to obtain the reduced time evolution of the open system, the so-called reduced dynamics, reads

\begin{equation}
 \label{reduced}
 \rho_t  = \mbox{Tr}_{\text{E}}[\mf{U}_t\Phi(\rho_0)\mf{U}_t^{\dagger}] \equiv T_t(\rho_0).
\end{equation}

Above, $\rho_0$ specifies the state of the open system at $t=0$. The map $\Phi$ assigns to each initial state $\rho_0$ a single
state $\Phi(\rho_0)$ of the total system. The assignment map must be chosen properly so that $T_t$ can be well-defined~\cite{reduced1,
reduced2,reduced3}. For instance, if no correlations between the systems are initially present, then $\Phi(\rho_0)=\rho_0\otimes\omega$,
for some initial state of the environment $\omega$. It is worth mentioning that, if the initial state cannot be factorized, the definition
of $\Phi$ is not accessible~\cite{korelacje,erratum}. The unitary operator $\mf{U}_t=\exp(-\rmi\mf{H}t)$ describes the time evolution of
the total system.

The map $\mbox{Tr}_{\text{E}}(\cdot):\mathcal{T}(\mathcal{H}_{\text{E}}\oplus\mathcal{H}_{\text{E}})\rightarrow M_2(\mathbb{C})$
denotes the so-called partial trace:

\begin{equation}
 \label{partial}
  \mbox{Tr}_{\text{E}}
  \begin{bmatrix}
   M_{11} & M_{12} \\
   M_{21} & M_{22}
   \end{bmatrix}
   =
   \begin{pmatrix}
   \mbox{Tr}M_{11} & \mbox{Tr}M_{12} \\
   \mbox{Tr}M_{21} & \mbox{Tr}M_{22}
   \end{pmatrix}.
 \end{equation}

$\mbox{Tr}(\cdot)$ refers to the usual trace operation on $\mc{H}_{\text{E}}$,
$\mathcal{T}(\mathcal{H}_{\text{E}}\oplus\mathcal{H}_{\text{E}})$ denotes the
Banach space of trace classes operator with the trace norm: $\|A\|_{1}=\mbox{Tr}(\sqrt{AA^{\dagger}})$, whereas
$M_2(\mathbb{C})$ is the Banach space of $2\times 2$ complex matrices. Note, the partial trace is a linear operation
transforming the block operator matrices (square brackets) to the ordinary matrices (round brackets).

\section{Main results}
\label{main}

Fixed point theorems like Banach or Schauder indicate the existence of stationary states for a given evolution $T_t$.
However, there is no general analytical procedure to obtain explicit form of such states. In this section we propose
a method of deriving stationary states for two-level open quantum systems. The generalization to the higher-dimensions
seems to be possible. However, we will not deal with this issue in this paper. We begin with some definitions.

\begin{deff}
The density matrix $\rho$ is said to be a stationary state if it is invariant
with respect to reduced evolution, $T_t(\rho)=\rho$.
\end{deff}

\begin{deff}
Let $X$ be an operator acting on the Hilbert space $\mc{H}_{\rm{E}}$. The subset $\Gamma_X$ of $\mc{H}_{\rm{E}}\oplus\mc{H}_{\rm{E}}$
defined as

   \begin{equation}
      \label{space}
      \Gamma_X:= \left\{
         \begin{bmatrix}
              \ket{\psi}\\
              X\ket{\psi}
         \end{bmatrix}
         :\ket{\psi}\in\mathcal{D}(X)\subset\mc{H}_{\text{E}}
              \right\}
    \end{equation}
is said to be the graph of $X$.
\end{deff}

The graph of a linear and closed operator is a subset of the Hilbert space, which is a Hilbert space itself equipped with the inner product

\begin{equation}
\label{inner}
\braket{\Psi_1}{\Psi_2} = \braket{\psi_1}{\psi_2}+\braket{\phi_1}{\phi_2},\quad
\ket{\Psi_i}=
\begin{bmatrix}
 \ket{\psi_i}\\
 \ket{\phi_i}
\end{bmatrix}\in\Gamma_X \quad\text{(}i=1,2\text{)}.
\end{equation}

$\braket{\psi}{\phi}$ is an inner product on $\mc{H}_{\text{E}}$. It is a known fact (see Lemma $5.3$ in~\cite{graph}) that the graph $\Gamma_X$
is $\mf{H}-$invariant, that is $\mf{H}\left(\Gamma_X\cap\mathcal{D}(\mf{H})\right)\subset\Gamma_X$ if and only if $X$ is a bounded solution (with
$\text{Ran}(X|_{\mathcal{D}(H_{\p})})\subset\mathcal{D}(H_{\m})$) of the \emph{Riccati} equation:

\begin{equation}
 \label{ricc}
       XVX+XH_{\p}-H_{\m}X-V^{\dagger}=0\quad \text{on}\quad\mathcal{D}(H_{\p}).
\end{equation}

Along with the equation above we introduce the dual Riccati equation, namely

\begin{equation}
 \label{ricc2}
       YV^{\dagger}Y+YH_{\m}-H_{\p}Y-V=0\quad \text{on}\quad\mathcal{D}(H_{\m}).
\end{equation}

It is proved in~\cite{graph} that $Y=-X^{\dagger}$ is a solution (with $\text{Ran}(X^{\dagger}|_{\mathcal{D}(H_{\m})}\subset\mathcal{D}(H_{\p})$)
of~(\ref{ricc2}) if and only if the orthogonal complement of $\Gamma_X$, \ie, the subspace

\begin{equation}
      \label{space2}
      \Gamma_X^{\bot}= \left\{
         \begin{bmatrix}
              -X^{\dagger}\ket{\psi}\\
              \ket{\psi}
         \end{bmatrix}
         :\ket{\psi}\in\mathcal{D}(X^{\dagger})\subset\mc{H}_{\text{E}}
              \right\}
    \end{equation}

is $\mf{H}$-invariant. It is straightforward to see that a bounded operator $X$ solves~(\ref{ricc}) if and only if $-X^{\dagger}$
is a solution of~(\ref{ricc2}). Therefore, $\Gamma_X$ as well as $\Gamma_X^{\bot}$ are $\mf{H}-$invariant if and only if $X$ is a
bounded solution of~(\ref{ricc}). In other words, $\Gamma_X$ is reducing subspace of $\mf{H}$ if and only if $X$ is a bounded solution
of~(\ref{ricc}). From considerations above follow also that $\Gamma_X$ and $\Gamma_X^{\bot}$ are $\mf{U}_t-$invariant.

\begin{deff}
\label{def}
Elements from the graph and its orthogonal complement are denoted by $\ket{X_{\psi}}$ and $\ket{X^{\psi}}$, respectively. The Riccati
states are defined as $\rho_{\psi}=\rm{Tr}_{\rm{E}}(\varrho_{\psi})$ and $\rho^{\psi}=\rm{Tr}_{\rm{E}}(\varrho^{\psi})$, where
$\varrho_{\psi}:=\ket{X_{\psi}}\bra{X_{\psi}}$ and $\varrho^{\psi}:=\ket{X^{\psi}}\bra{X^{\psi}}$.
\end{deff}

The vectors $\ket{X_{\psi}}$ and $\ket{X^{\psi}}$ are \emph{not} normalized with respect to the norm induced by the inner
product~(\ref{inner}). Moreover, the states $\varrho_{\psi}$ and $\varrho^{\psi}$ are not factorisable (\ie, correlations occur),
unless $X\ket{\psi}\sim\ket{\psi}$ and $X^{\dagger}\ket{\psi}\sim\ket{\psi}$, respectively. However, they are
$\mf{U}_t(\cdot)\mf{U}_t^{\dagger}-$invariant, which is obvious because the vectors $\ket{X_{\psi}}$ and $\ket{X^{\psi}}$ are
$\mf{U}_t-$invariant. As a consequence, the Riccati states $\rho_{\psi}$ and $\rho^{\psi}$ are $T_t-$invariant, where the map
$T_t$ has been defined in~(\ref{reduced}). Therefore, the set of all the Riccati states is invariant under the time evolution.
Nevertheless, the Riccati states are not the stationary states, in general. However, we show that the latter can be found among
the Riccati states. To be specific, we will prove the following

\begin{thh}
  \label{th1}
Let $X$ be a bounded solution of the Riccati equation~\eref{ricc}. Then

\begin{enumerate}

\item[i)] the Riccati state $\rho_{\psi}$ is a stationary state if the vector
          $\ket{\psi}$ is a eigenvector of the operator $Z_{\p}\equiv H_{\p}+VX :\mathcal{D}(H_{\p})\rightarrow\mc{H_{\text{E}}}$,

\item[ii)] the Riccati state $\rho^{\phi}$ is a stationary state if the vector
          $\ket{\phi}$ is a eigenvector of the operator $Z_{\m}\equiv H_{\m} -
          V^{\dagger}X^{\dagger}:\mathcal{D}(H_{\m})\rightarrow\mc{H_{\text{E}}}$.

\end{enumerate}
\end{thh}

\begin{proof}
Let $Z_{\p}\ket{\psi}=\lambda\ket{\psi}$ for $\lambda\in\mb{C}$ and $\ket{\psi}\in\mathcal{D}(H_{\p})$.
From~(\ref{ricc}) we obtain that $V^{\dagger}+H_{\m}X=XZ_{\p}$, hence in view of~(\ref{bom}) the last equality
leads to $\mf{H}\ket{X_{\psi}} = \lambda\ket{X_{\psi}}$. Thus, the vector state $\ket{X_{\psi}}$ is the
eigenvector of the total Hamiltonian with the corresponding eigenvalue $\lambda$. Since $\mf{H}$ is self-adjoint
we have $\lambda\in\mb{R}$ and in consequence $\mf{U}_t\varrho_{\psi}\mf{U}_t^{\dagger}=\varrho_{\psi}$,
which ultimately leads to $T_t(\rho_{\psi})=\rho_{\psi}$.

In a comparable manner we have $\mf{H}\ket{X^{\phi}}=\xi\ket{X^{\phi}}$ for $\xi\in\mb{R}$ and
$\ket{\phi}\in\mc{D}(H_{\m})$ so that $Z_{\m}\ket{\phi}=\xi\ket{\phi}$. Just as before
$\mf{T}_t(\varrho^{\phi})=\varrho^{\phi}$, therefore $T_t(\rho^{\phi})=\rho^{\phi}$.
\end{proof}

At this point, some remarks, regarding theorem given above, should be made.

\begin{Rem}
The question whether all stationary states are Riccati states or if it is possible that stationary states exist that
are not Riccati states is still open.
\end{Rem}

\begin{Rem}
 \label{re1}
Since the space $\Gamma_{X}$ is closed, we have the following decomposition $\mc{H}_{\text{tot}}=
\Gamma_{X}\oplus\Gamma_{X}^{\bot}$. Thus the total Hamiltonian is similar to certain block diagonal
operator  matrix, $\mf{S}^{-1}\mf{H}\mf{S}=\mf{H}_{\text{d}}$, where

\begin{equation}

 \mf{H}_{\text{d}}=
      \begin{bmatrix}
        Z_{\p} & 0 \\
        0 & Z_{\m}
      \end{bmatrix}
\quad\text{with}\quad\mathcal{D}(Z_{\spm})=\mathcal{D}(H_{\spm})\quad\text{and}\quad
 \mf{S}=
      \begin{bmatrix}
     \idE & -X^{\dagger} \\
        X & \idE
    \end{bmatrix}.
\end{equation}

This, implies that $\sigma(\mf{H})=\sigma(Z_{\p})\cup\sigma(Z_{\m})$. Therefore, the eigenvalues of $Z_{\spm}$ are exactly
the eigenvalues of the Hamiltonian $\mf{H}$.
\end{Rem}

\begin{proof}
Let $X$ be a bounded solution of~(\ref{ricc}). $V$ is assumed to be bounded as well. From the definition of $Z_{\spm}$ we have
$\mathcal{D}(Z_{\spm})=\mathcal{D}(H_{\spm})$, and thus $\mc{D}(\mf{H})=\mc{D}(\mf{H}_{\text{d}})$. Since $X$ solves the Riccati
equation~(\ref{ricc}), it is clear that $\mf{H}\mf{S}=\mf{S}\mf{H}_{\text{d}}$. To prove $\mf{H}\sim\mf{H}_{\text{d}}$ we will
show that $\mf{S}$ is invertible and $\mf{S}^{-1}$ is bounded. Indeed, $\mf{S}=\mf{I}+\mf{X}$, where

   \begin{equation}
       \mf{X} =
       \begin{bmatrix}
        0 & -X^{\dagger} \\
        X & 0
        \end{bmatrix}.
       \end{equation}

Since $\mf{X}^{\dagger}=-\mf{X}$, the spectrum of $\mf{X}$ is a subset of the imaginary axis. In particular, $-1\not\in\sigma(\mf{X})$
thus $0\not\in\sigma(\mf{S})$ and, hence, $\mf{S}$ has a bounded inverse.

\end{proof}
\begin{Rem}
  \label{Rem4}
The stationary states $\rho_{\psi}$, $\rho^{\phi}$ indicated in theorem~\ref{th1} are given by

 \begin{equation}
  \label{explicit}
  \rho_{\psi}=A
      \begin{pmatrix}
          1 & \langle X\rangle^{*}_{\psi} \\
           \langle X\rangle_{\psi} & \|X\psi\|^2
      \end{pmatrix}\quad\text{and}\quad
        \rho^{\phi}=B
      \begin{pmatrix}
          \|X^{\dagger}\phi\|^2 & -\langle X\rangle^{*}_{\phi} \\
          -\langle X\rangle_{\phi} & 1
      \end{pmatrix},
 \end{equation}

where $\ket{\psi}\in\mathcal{D}(H_{\p})$ and $\ket{\phi}\in\mathcal{D}(H_{\m})$ are the eigenvectors of $Z_{\m}$ and $Z_{\p}$,
respectively. $A=\mbox{\emph{Tr}}(\rho_{\psi})$, $B=\mbox{\emph{Tr}}(\rho^{\phi})$ are the normalization constants and $\langle
X\rangle_{\varphi}=\bra{\varphi}X\ket{\varphi}$.
\end{Rem}

\begin{proof}
Since $\mbox{Tr}\ket{\psi}\bra{\phi}=\braket{\phi}{\psi}$, the equations~(\ref{explicit}) can be obtained directly from the
definition~(\ref{def}) and the formula~(\ref{partial}).
\end{proof}

\section{Examples}

\subsection{Spin-boson model}

In this subsection we will demonstrate an application of the presented method to a non-trivial example, namely, the paradigmatic
spin-boson model~\cite{SB_paradigm, SB_Parity}. Assume that the Hamiltonian of the qubit (spin-half) and its environment (boson)
are in the following forms

\begin{equation}
 \label{eq:QandE}
    H_Q = \beta\sigma_z+\alpha\sigma_x \quad\text{and}\quad
    H_E = \omega a^{\dagger}a,
\end{equation}

respectively. For the sake of simplicity, we consider only the case of a single boson. The interaction between the systems reads

\begin{equation}
 \mf{H}_{\text{int}} = \sigma_z\ot (g^{\ast}a+ga^{\dagger})
                       \equiv\sigma_z\ot V.
\end{equation}

In the above description, $\sigma_x$, $\sigma_z$ are the standard Pauli matrices and $\alpha, \beta \in \mb{R}$. The creation
$a^{\dagger}$ and annihilation $a$ operators obey the canonical commutation relation (CCR) $[a,a^{\dagger}]=\idE$~\cite{sakurai}.
Parameters $\omega>0$ and $g\in\mathbb{C}$ represent the energy of the boson and the coupling constant between the qubit and the
environment, respectively.

If $\alpha=0$ (no energy exchange between the systems), the model can be solved, \ie, the reduced dynamics can be obtained
exactly~\cite{SB_Alicki,Luczka}. The solution describes the physical phenomena known as the pure decoherence or dephasing
~\cite{dajka_cat}. On the other hand, when $\alpha\not=0$ the exact solution in not known. The objective is to estimate the
stationary states for the latter case.

To proceed, we must clarify some technical aspects (\eg, domains of $H_{\spm}$). Clearly, the operators $a$ and $a^{\dagger}$ cannot
both be bounded since the trace of their commutator does not vanish~\cite{reed}. Therefore, the CCR holds only on some dense subspace
$\mc{D}_2$ of $\mc{H}_{\text{E}}$. Let $\mc{D}_1$ be the dense domain of both $a$ and $a^{\dagger}$, on which they are mutually adjoint,
that is, $(a^{\dagger})^{*}=a$, $a^{*}=a^{\dagger}$. At this point, it is not obvious that the sets $\mc{D}_1$, $\mc{D}_2$, having desire
properties, exist. The detailed construction can be found in~\cite{sz0} and the right choice is given by

\begin{equation}
 \mc{D}_k = \left\{\ket{\psi}\in\mc{H}_{\text{E}}: \sum_{n=0}^{\infty}n^k|\braket{\psi}{\phi_n}|^2<\infty\right\},\quad k =1,2.
\end{equation}

On $\mc{D}_1$ the creation and annihilation operators can be defined explicitly as (see also~\cite{sz1,sz2,sz3})

\begin{equation}
 \label{aa}
   a\ket{\phi} = \sum_{n=1}^{\infty}\sqrt{n}\braket{\phi_n}{\phi}\ket{\phi_{n-1}}, \quad
   a^{\dagger}\ket{\phi} = \sum_{n=0}^{\infty}\sqrt{n+1}\braket{\phi_n}{\phi}\ket{\phi_{n+1}},
   \quad\ket{\phi}\in\mc{D}_1.
  \end{equation}

$\{\ket{\phi_n}\}_{n=0}^{\infty}$ is an orthonormal basis in $\mc{H}_{E}$. From~(\ref{aa}) follows that
$a^{\dagger}\ket{\phi_n}=\sqrt{n+1}\ket{\phi_{n+1}}$ and $a\ket{\phi_n}=\sqrt{n}\ket{\phi_{n-1}}$, which in most books
on quantum mechanics is a definition of the creation and annihilation operators. However, the operators defined
in such a way are not closed; nevertheless, they are closable and their closures are given by~(\ref{aa}).

Since $a$ is closed, $H_E$ is a positive self-adjoint operator and $\mc{D}_2\subset\mc{D}_1$ is a core of $a$ (see \eg, Theorem
4.2.1 in~\cite{blank}). Henceforward, we assume that the basis $\{\ket{\phi_n}\}_{n=0}^{\infty}$ is composed with the eigenvectors
of $H_E$. In this case we have

\begin{equation}
H_E\ket{\phi} = \sum_{n=0}^{\infty}\omega n\braket{\phi_n}{\phi}\ket{\phi_n}, \quad\ket{\phi}\in\mc{D}_2.
\end{equation}

By choosing a suitable coupling constant $g$, it is possible to make $H_{\spm} = H_E \pm V$ self-adjoint on $D_{\spm}:=\mc{D}_1\cap\mc{D}_2=\mc{D}_2$.
To see this, let us first note that $(g^*a+ga^{\dagger})^{*}\supset g^*a+ga^{\dagger}$, thus $V$ is Hermitian (symmetric). Since $H_E$ is
self-adjoint and $V$ Hermitian, it is sufficient to show that $V$ is relatively bounded with respect to $H_E$ ($H_E$ bounded) and with
$H_E$-bound less than one. The fundamental result of perturbation theory, known as the Kato-Rellich theorem~\cite{reed} assure Hermiticity
of $H_{\spm}$ in this case.

Recall that $B$ is $A$ bounded if i) $\mc{D}(A)\subset\mc{D}(B)$ and ii) $\|B\ket{\phi}\|^2\le a\|A\ket{\phi}\|^2 +
b\|\ket{\phi}\|^2$, for all $\ket{\phi}\in\mc{D}(A)$ and some nonnegative constants $a$, $b$. The infimum of all $a$ for which a corresponding
$b$ exists such that the last inequality holds is called $A-$bound of $B$. Note, sometimes it is convenient to replace the condition ii) by
the equivalent one: $\|B\ket{\phi}\|\le a\|A\ket{\phi}\| + b\|\ket{\phi}\|$.

It is not difficult to see that if two operators $B_i$, $i=1,2$ are bounded with respect to the same operator $A$ and their relative
bound are less than $b_1$ and $b_2$, respectively; then $a_1B_1+a_2B_2$ is also $A$ bounded and its relative bounded is less than
$|a_1|b_1+|a_2|b_2$ (see Lemma 6.1 in~\cite{gt}). In other words, the set of all $A$ bounded operator form a linear space. Therefore,
to see that $V$ is $H_E$ bounded it is sufficient to prove that both $a$ and $a^{\dagger}$ are $H_E$ bounded. To finish this note

\begin{eqnarray}
 \|a\ket{\phi}\|^2 &=& \sum_{n=1}^{\infty}n|\braket{\phi_{n-1}}{\phi}|^2 = \sum_{n=0}^{\infty}(n+1)|\braket{\phi_n}{\phi}|^2 \\\nonumber
 &\le& \sum_{n=0}^{\infty}n^2|\braket{\phi_n}{\phi}|^2 +\sum_{n=0}^{\infty}|\braket{\phi_n}{\phi}|^2
  = \omega^{-1}\|H_E\ket{\phi}\|^2 + \|\ket{\phi}\|^2.
\end{eqnarray}

In comparable manner one can also verify that $\|a^{\dagger}\ket{\phi}\|^2\le \omega^{-1}\|H_E\ket{\phi}\|^2$. Since $H_E-$bound of both
$a$ and $a^{\dagger}$ is less than one, the $H_E-$bound of $V$ is also less than one for $|z|<1/2$.

The block operator matrix representation of the spin-boson Hamiltonian reads

 \begin{equation}
  \label{eq:SBbom}
 \mf{H} =
   \begin{bmatrix}
    H_{\p} & \alpha \\
    \alpha & H_{\m}
   \end{bmatrix},
    \quad\text{where}\quad H_{\spm} = H_E\pm V
    \quad\text{and}\quad\mc{D}(\mf{H})=\mc{D}_2\oplus\mc{D}_2;
\end{equation}

the quantity $\alpha$ is understood as $\alpha\idE$. For the sake of simplicity we have set $\beta=0$, the example remains
non-trivial because $[H_{Q}\ot\idE,\mf{H}_{\text{int}}]\not=0$. The corresponding Riccati equation takes the form

\begin{equation}
 \label{eq:SBriccati}
\alpha X^2 + XH_{\p}-H_{\m}X-\alpha =0\quad\text{on}\quad\mc{D}_2.
\end{equation}

In order to solve this equation we define an operator $P$ as

\begin{equation}
 \label{parity}
 P\ket{\psi} = \sum_{n=0}^{\infty}e^{i\pi n}\braket{\phi_n}{\psi}\ket{\phi_n},\quad\ket{\psi}\in\mc{H}_{\text{E}}.
\end{equation}

Directly from~(\ref{parity}) we have $P=P^{\dagger}$ and $P^2=\idE$, thus $P$ is both self-adjoint and unitary. Formally, $P$
can be written as $P=\exp(i\pi H_{E})$, however, unlike $H_{E}$, $P$ is everywhere defined. The Hellinger-Toeplitz theorem
guaranties that $P$ is bounded, which can also be seen directly. Indeed, form unitarity we obtain $\|P\psi\|=\|\psi\|$, for
$\ket{\psi}\in\mc{H}_{\text{E}}$, hence $\|P\|=1$. $P$ is, in fact, the bosonic parity operator~\cite{Bender}. We will show
that $X=P$ solves~(\ref{eq:SBriccati}). Since $P^2=\idE$ it is sufficient to show that

\begin{equation}
 \label{syl}
 PH_{\m}-H_{\p}P = 0
  \quad\text{or equivalently}\quad
  PH_{\spm}P=H_{\smp}\quad\text{on}\quad\mc{D}_2.
\end{equation}

In order to prove~(\ref{syl}) let us first note that $P\ket{\psi}\in\mc{D}_2$, for $\ket{\psi}\in\mc{D}_2$, which means
$\mbox{Ran}(P|_{\mc{D}_2})\subset\mc{D}_2$. This follows from $|e^{i\pi n}|=1$. Furthermore, $PH_EP=H_E$ and $PVP=-V$.
The first equality is obvious, while the second one follows from $PaP=-a$ and $Pa^{\dagger}P=-a^{\dagger}$. And as a result we
obtain~(\ref{syl}).

The stationary states in this example can be written as

 \begin{equation}
  \label{eq:SBstate}
\rho_{\spm} = \frac{1}{2}
   \begin{pmatrix}
    1 & r_{\spm} \\
     r^{*}_{\spm} & 1
   \end{pmatrix},
   \quad\text{where}\quad
   r_{\spm}=\bra{\psi_{\spm}}P\ket{\psi_{\spm}}
\end{equation}

and $\ket{\psi_{\spm}}$ are the eigenvectors of $Z_{\spm}=H_{\spm}\pm\alpha P$. Unfortunately, the solution of this eigenproblem is not known
for $\alpha\not=0$. However, one can determine certain bounds on $r_{\spm}$ using properties of $P$ and $\rho_{\spm}$. First, $r_{\pm}$ are
real numbers because $P$ is self-adjoint. From the non-negativity of $\rho_{\spm}$ we obtain that $r_{\spm}\in[-1,1]$.

Equation~(\ref{eq:SBstate}) provides estimation of the stationary states of the spin-half immersed within the bosonic bath. This result has been
obtained without any approximations. Of course, $r_{\spm}$ can be computed approximately with the use of known methods. It is important to stress
that to obtain the exact reduced dynamics for the model in question one needs to resolve an eigenvalue problem for $Z_{\spm}$.

\subsection{Commuting environment}

In the second example we consider the Hamiltonian in the following form

\begin{equation}
 \label{spin}
   \mf{H} = \alpha\sigma_x\ot\idE+\idQ\ot H_0+\sigma_z\ot H_1, \quad\alpha\not=0.
\end{equation}

Here we assume that linear, bounded operators $H_0$, $H_1$ commute. In this example we impose restriction to the spectra of $H_0$ and $H_1$, \ie{}
$\sigma(H_0)$, $\sigma(H_1)$ are discrete and non-degenerated. The Hamiltonian~(\ref{spin}) describes a qubit in contact with an environment and
in the presence of the magnetic field $\vec{B}=B\hat{e}_{x}$, where $B\sim\alpha$. Examples of such systems occur in the literature,
\eg~\cite{gardas3,spinStar1,spinStar2}. The block operator matrix representation of~(\ref{spin}) is given by~(\ref{eq:SBbom}) with
$H_{\spm}=H_0\pm H_1$ and $\mc{D}(\mf{H})=\mc{H}_{\text{E}}$. The corresponding Riccati equation reads~(\ref{eq:SBriccati}).

Using the fact that $H_0$ and $H_1$ commute, so they have a common set of eigenvectors, we write

\begin{equation}
H_0\ket{\phi_n}=\lambda_n\ket{\phi_n}\quad\text{and}\quad H_1\ket{\phi_n}=\xi_n\ket{\phi_n},
\end{equation}

where $\lambda_n\in\sigma(H_0)$, $\xi_n\in\sigma(H_1)$ and $\braket{\phi_n}{\phi_m}=\delta_{nm}$, for $n,m\in\mb{N}$.
The Riccati equation has a positive and self-adjoint solution $X=f(H_1)$, where the function $f$ is given by

\begin{equation}
 f(x)=\frac{\sqrt{x^2+\alpha^2}-x}{\alpha},\quad\text{for}\quad x\in\sigma(H_1).
\end{equation}

Unlike to the spin-boson model, in this case the eigenproblem for $Z_{\spm}$ can be readily solved. Indeed, we have

\begin{equation}
Z_{\spm}\ket{\phi_n} = \left[\lambda_n\pm\xi_n f(\xi_n)\right]\ket{\phi_n}.
\end{equation}

According to the Remark~\ref{Rem4} we obtain

\begin{equation}
  \label{final}
 \rho_n^{\p} = C_n
   \begin{pmatrix}
   1 & f(\xi_n) \\
   f(\xi_n)^{\ast} & |f(\xi_n)|^2
   \end{pmatrix}
   \quad\text{and}\quad
   \rho_n^{\m} = C_n
   \begin{pmatrix}
   |f(\xi_n)|^2 & -f(\xi_n)^{\ast} \\
   -f(\xi_n) & 1
   \end{pmatrix},
\end{equation}
where $ C_n=(1+|f(\xi_n)|^2)^{-1}$.

In this case there are no initial correlations between the systems since $X\ket{\phi_n}\sim\ket{\phi_n}$. We wish to emphasize that there
may exist other solutions of the Riccati equation. For an explicit example see Ref.~\cite{gardas2}.

\subsection{Sylvester equation}

In the case $V=0$, the Hamiltonian~(\ref{bom}) is already in a block diagonal form. One can notice that the Riccati equation simplifies to
the Sylvester equation:

\begin{equation}
 XH_{\p}-H_{\m}X = 0\quad\text{on}\quad\mc{D}(H_{\p}).
\end{equation}

There exist at least one solution, namely $X=0$. The corresponding stationary states are given by the projections
$P_0=\mbox{diag}(0,1)=\ket{0}\bra{0}$ and $P_1=\mbox{diag}(1,0)=\ket{1}\bra{1}$.

\subsection{All unbounded entries}

In the last example we consider an interesting example in which all the entries of $\mf{H}$ are unbounded, but still the solution of the
Riccati equation exists as a bounded operator. To see this, let us choose $H_{\spm}=H_0$ with domain $\mathcal{D}(H_0)$ and let us assume
that $V$ is self-adjoint with domain $\mathcal{D}(V)$. Then, the Riccati equation

\begin{equation}
  \label{cap}
  XVX+XH_0-H_0X-V = 0
  \quad\text{on}\quad
  \mathcal{D}(H_0)\cap\mathcal{D}(V)
\end{equation}

has at least two bounded solutions, that is $X_{\spm}=\pm\idE$. The stationary states read

\begin{equation}
 \rho_{\spm} =
 \frac{1}{2}
      \begin{pmatrix}
       1 & \varepsilon_{\spm} \\
       \varepsilon_{\spm} & 1
      \end{pmatrix},\quad\text{where}\quad\varepsilon_{\spm}=\pm 1.
\end{equation}

An example in which all the entries of a block operator matrix are bounded and the solution of the Riccati equation is unbounded has been
provided in~\cite{Vadim}.

\section{Summary.}

In this paper we have proposed the method of calculating the stationary states for two-level open quantum system. The theory of block operator
matrices has been adapted to achieve this purpose. In particular, we have related the solution of the algebraic Riccati equation to stationary
states. In the presented method, the stationary states are generated from the stationary states of the total system by tracing out the environment.
Our investigation includes the case when the initial system-environment correlations occur. In fact, this case is embedded in the method since the
eigenstates of the total Hamiltonian are entanglement.

At the end, we want to stress that the method cannot be used when the total Hamiltonian is not known. Such a situation arises \eg, when the details
of the interaction between the systems are not accessible. We hope that despite aforementioned weaknesses the results of the paper may serve as a
starting point for further investigations.

\ack
The authors are deeply grateful to the referee for very careful reading
of the original manuscript and valuable suggestions.
B. Gardas would like to thank Jerzy Dajka for his suggestions and comments.
This work was financed from the Polish science budget resources in the years
2010-2013 as a research projects: grant number N N519 442339 and project number
IP 2010 0334 70.

\section*{References}


TITLE: Prove that for any given positive integer $n$,there exists a number having digits $0,1$ which is divisible by $n$.
QUESTION [0 upvotes]: Prove that for any given positive integer $n$,there exists a number having digits $0,1$ which is divisible by $n$.

Let that number be of the general form: $x=\overline {b_kb_{k-1}...b_1b_0}$ where $b_i\in \{0,1\}$.How can we construct $x$ to be divisible by the given $n$?

REPLY [3 votes]: Let $$a_k = \underbrace{111...11}_{k} = {10^k-1\over 9}$$
So, among $a_1,a_2,...a_{n+1}$ two must have the same remainder modulo $n$, say $a_i$ and $a_j$ and $i>j$. Thus $$n\mid a_i-a_j ={10^i-10^j\over 9} =10^j{10^{i-j}-1\over 9}= \underbrace{111...11}_{i-j}\underbrace{000...00}_{j}$$

REPLY [2 votes]: Note that $10^{k \, \varphi(n)} \equiv 1,$ where $\varphi$ is the Euler phi function and $k\in \mathbb{N}$, and then see that
$$
\sum_{k=1}^n 10^{k \, \varphi(n)} \equiv 0 \pmod n
$$
Thus the integer whose decimal expansion has a 1 at positions $k \, \varphi(n)$ for $k=1,\dots,n$ and zeros everywhere else is divisible by $n.$
Alas, this only works if $n$ is relative prime to 10, and needs to be modified otherwise.


TITLE: Usefulness of Variance
QUESTION [5 upvotes]: I've had a look for intuitive explanations of the variance of an RV (e.g. Intuitive explanation of variance and moment in Probability) but unfortunately for me, I still don't feel comfortable with the concept. Why would you opt to use variance over the standard deviation (which usefully is in the same units as the expectation)?
Also, if the expectation,
$E(X) = \Sigma_{i=0}^n i P(X = i)$,
what is $E(X^2)$? Is it simply
$E(X^2) = \Sigma_{i=0}^n i P(X^2 = i)$?

REPLY [0 votes]: Variance is a natural parameter for the normal distribution.


TITLE: Show that for every $\epsilon > 0$ there is $\delta > 0$ such that$\int \limits _{E_0}|f| < \epsilon$...
QUESTION [1 upvotes]: Assume $E \subset R_n$ and $f \in L(E)$. Show that for every $\epsilon > 0$ there is $\delta > 0$ such that $\int \limits _{E_0}|f| < \epsilon$ for every measurable $E_0 \subset E$ such that $|E_0| < \delta$. We are given that $|f| = \lim \limits_{k = \infty} f_k$ where each $f_k : E \rightarrow [0,\infty]$ is simple, Lebesgue measurable, and $f_k \leq f_{k+1}$.
This is confusing off the bat because the only $\epsilon - \delta$ definition that I know is: for $\epsilon > 0$ there is $\delta > 0$ such that $|f(x) - L| < \epsilon$ and $|x - x_0|< \delta$.
With that being said, my approach to this problem was to manipulate $\int \limits _{E_0}|f|$ and get some form of $a*|E_0|$, where $a$ is a constant. I know that $\int\limits _{E_0} a = a |E_0|$. The problem is I can't seem to get $\int \limits _{E_0}|f|$ to be the integral of a constant. Any ideas/suggestions?
EDIT: I think that I should start with $|\int \limits_{E_0}f(x) - L| <\epsilon$ and manipulate that first to get $\int \limits _{E_0}|f|<\epsilon$. Is this right?
EDIT (Attempted simple case soln, see user 15464 answer):
For $\epsilon >0$ there is $\delta > 0$ such that $|E| < \delta$ if $\int \limits_E \phi < \epsilon$
$$
\int \limits_E \phi < \epsilon \Rightarrow \int \limits_{E_i} a_i < \epsilon
$$
where $a_i$ is a value of $\phi$ on $E_i$ which is a disjoint subset of $E$.
$$
\int \limits_{E_i} a_i < \epsilon \Rightarrow a_i |E_i| < \epsilon \Rightarrow |E_i| < \frac {\epsilon}{a_i }
$$
then I can take the sum of all the subsets ($m \in N$)
$$
\sum \limits_{i=1}^{m} |E_i| = |E| < \sum \limits_{i=1}^{m} \frac {\epsilon}{a_i } = \frac {\epsilon}{\sum \limits_{i=1}^{m} a_i } = \epsilon \sum \limits_{i=1}^{m} \frac {1}{a_i }
$$
finally $$\delta = \epsilon \sum \limits_{i=1}^{m} \frac {1}{a_i } > 0$$

REPLY [1 votes]: Notice that if $\phi$ is a simple function, then $\phi$ takes on only finitely many values, so in particular $|\phi|$ is bounded. Can you use this to prove the result in the case of a simple function?
Once you have the result for a simple function, and you want to extend it to the general $f$, choose a simple function $\phi$ with $\int |f - \phi| < \epsilon/2$. Then choose $\delta$ such that if $|E| < \delta$, then $\int_E |\phi|  < \epsilon/2$. Write
$$\int_E |f| \leq \int_E |f - \phi| + \int_E |\phi|$$
and bound each term separately.