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

Quantum mechanics takes place in infinite dimensional Hilbert space. Naturally, any numerical approach to solving the equations of quantum mechanics -- or any other physical system with a continuum description -- involves a finite truncation of this space. When we turn to \emph{many-body} quantum mechanics, the dimension of Hilbert space necessarily grows exponentially with the number of particles relative to the finite truncation used for a single particle. To be specific, if a single particle is described by a wavefunction $\psi(\vec{r})$ defined on a real-space grid of linear size $L$ (with $L^3$ points in three dimensions), the wavefunction of $N$ particles $\Psi(\vec{r}_1,\ldots,\vec{r}_N)$ is defined on a grid in $3N$ dimensions of $L^{3N}$ points.

Traditionally, this challenge has been dealt with by considering many body wavefunctions of restricted form. For example, the Hartree--Fock method employs factorized wavefunctions\footnote{We leave aside the issue of the statistics of indistinguishable particles for the moment.}
$$
\Psi(\vec{r}_1,\ldots,\vec{r}_N)=\psi_1(\vec{r}_1)\ldots \psi_N(\vec{r}_N).
$$
This reduces the memory cost to linear in the number of particles,\footnote{The computational complexity of the Hartree--Fock method scales as the cube of the number of basis functions used to represent the $\psi_i(\vec{r})$.} but represents a drastic simplification that performs especially poorly when the interaction between particles is strong. Over the years many \emph{post--Hartree--Fock} hand-crafted improvements of the many-body wavefunction have been introduced, including Jastrow factors, and the coupled cluster and configuration interactions methods (\cite{Foulkes:2001aa}).

The exponential growth in the complexity of many-body quantum mechanics closely parallels the curse of dimensionality encountered in computer vision and other traditional applications of machine learning. Therefore it is natural that neural methods that have recently proven successful in the latter domain be applied to quantum mechanical calculations.

\subsection{Deep Learning Approaches to Quantum Ground States} \label{sec:dlq}

Beginning with the work of \cite{Carleo:2017aa}, the past few years have seen numerous attempts to leverage the expressive power of deep networks to solve the quantum many-body problem.\footnote{There is earlier work on few-body quantum mechanics from before the deep learning era (\cite{Lagaris:1997aa}). A separate strand of work is concerned with efficient parameterization of interatomic potentials that have been calculated by other means, see e.g. (\cite{Behler:2016aa}).}

This initial work dealt with lattice models of spins that are a mainstay of quantum condensed matter physics, representing the wavefunction using a Restricted Boltzmann Machine. Since then neural representations of many body wavefunctions have multiplied to include lattice systems of fermions (\cite{Nomura:2017aa}) and bosons (\cite{Saito:2017aa}), particles in continuous space (\cite{Ruggeri:2018aa,Kessler:2019aa}), systems with symmetry (\cite{Choo:2018aa}), and the many-electron problem (\cite{Han:2019aa,Pfau:2019aa,Hermann:2019aa}).

All of the above-mentioned applications of neural networks to quantum mechanics have used the Schr\"odinger picture, in which a representation is sought of the wavefunction. There are several mathematically equivalent formulations of quantum mechanics that provide alternatives to the Schr\"odinger equation, and hence potential alternative routes for the application of machine learning. In this work we show how \emph{reinforcement learning} may be used to solve problems in quantum mechanics via the \emph{path integral} representation.

\subsection{Feynman--Kac Representation of the Ground State}

Feynman's path integral formulation of quantum mechanics starts from the Schr\"odinger equation (\cite{Feynman:1965aa})
\begin{equation}
  i\frac{\partial\psi(\vec{r},t)}{\partial t} = \left[H\psi\right](\vec{r},t),
\end{equation}
where for a single particle the Hamiltonian operator describing motion in a potential $V(\vec{r})$ has the form
\begin{equation}
  \left[H\phi\right](\vec{r}) = -\frac{1}{2}\nabla^2\phi(\vec{r}) + V(\vec{r})\phi(\vec{r}).
\end{equation}
The wavefunction at time $t_2>t_1$ can be expressed in terms of the wavefunction at time $t_1$ by
\begin{equation}
  \psi(\vec{r}_2,t_2) = \int d\vec{r}_1 \mathcal{K}(\vec{r}_2,t_2;\vec{r}_1,t_1)\psi(\vec{r}_1,t_1),
\end{equation}
where the \emph{propagator} $\mathcal{K}(\vec{r}_2,t_2;\vec{r}_1,t_1)$ is the Green's function of the Schr\"odinger equation
\begin{equation}
  \left(i\frac{\partial}{\partial t_2} - H_{\vec{r}_2}\right)\mathcal{K}(\vec{r}_2,t_2;\vec{r}_1,t_1) =i\delta(\vec{r}_1-\vec{r}_2)\delta(t_1-t_2).
\end{equation}
$\mathcal{K}(\vec{r}_2,t_2;\vec{r}_1,t_1)$ has the physical interpretation of the \emph{probability amplitude} to move from $\vec{r}_1$ at time $t_1$ to $\vec{r}_2$ at time $t_2$. The \emph{Born rule} then states that the square modulus of the probability amplitude gives the probability for this event.

The path integral is a representation of the propagator as a `sum over paths', written formally as
\begin{equation}\label{eq:real-time}
  \mathcal{K}(\vec{r}_2,t_2;\vec{r}_1,t_1) = \int_{\vec{r}(t_1)=\vec{r}_1 \atop \vec{r}(t_2)=\vec{r}_2} \mathcal{D}\vec{r}(t)\exp\left(i\int_{t_1}^{t_2} L(\vec{r},\dot{\vec{r}})dt\right)
\end{equation}
where $L(\vec{r},\vec{v}) = \frac{1}{2}\vec{v}^2 - V(\vec{r})$ is the classical Lagrangian function of the system, and the domain of the integral is all paths satisfying the stated endpoint conditions.

Despite being successfully wielded by physicists for decades Feynman's original idea has never found rigorous mathematical formulation due to the difficulties of defining a suitable measure on the path space. \cite{Kac:1949aa} discovered, however, that a fully rigorous path integral formula exists for the heat-type equations
\begin{equation}\label{eq:im-time}
  \frac{\partial\psi(\vec{r},t)}{\partial t} = -\left[H\psi\right](\vec{r},t),
\end{equation}
also known as the imaginary time Schr\"odinger equation. Moving to imaginary time makes the exponent in \eqref{eq:real-time} real. Kac observed that the part of the exponent arising from the kinetic energy could be interpreted as a measure on Brownian paths, leading to the \emph{Feynman--Kac (FK) formula} for the solution of \eqref{eq:im-time}
\begin{equation}\label{eq:FK}
  \psi(\vec{r}_2,t_2) =  \E_{\vec{r}(t_2)=\vec{r}_2}\left[\exp\left(-\int_{t_1}^{t_2}V(\vec{r}(t))dt\right)\psi(\vec{r}(t_1),t_1)\right],
\end{equation}
where the expectation is over Brownian paths finishing at $\vec{r}_2$ at time $t_2$. In this way quantum mechanics is brought into the realm of stochastic processes, albeit in `imaginary time'. Whilst this formulation is therefore not of direct utility in studying quantum dynamics, it provides a very useful tool for studying \emph{ground states}. This is because the propagator $K(\vec{r}_2,t_2;\vec{r}_1,t_1)$ for \eqref{eq:im-time} has a spectral representation in terms of the eigenfunctions $\varphi_n$ and eigenenergies $E_n$ of the time independent Schr\"odinger equation $H\varphi_n = E_n\varphi_n$ as
\begin{align}\label{eq:long-time}
  K(\vec{r}_2,t_2;\vec{r}_1,t_1) &= \sum_n \varphi_n(\vec{r}_2)\varphi^*_n(\vec{r}_1)e^{-E_n(t_2-t_1)}\\
  &\longrightarrow \varphi_0(\vec{r}_2)\varphi^*_0(\vec{r}_1)e^{-E_0(t_2-t_1)} \qquad \text{ as } t_2-t_1\to\infty.
\end{align}
Thus, as $t_2-t_1\to\infty$, only the ground state contributes.

The FK formula defines a new path measure $\mathbb{P}_\text{FK}$ that differs from the Brownian measure $\mathbb{P}_0$ by the Radon--Nikodym derivative
\begin{equation}\label{eq:RN}
  \frac{d\mathbb{P}_\text{FK}}{d\mathbb{P}_0} = \mathcal{N}\exp\left(-\int_{t_1}^{t_2}V(\vec{r}(t))dt\right)
\end{equation}
where $\mathcal{N}$ is a normalization factor. Intuitively, \eqref{eq:RN} describes paths that spend more time in the attractive regions of the potential ($V(\vec{r})<0$) and less time in the repulsive regions ($V(\vec{r})>0$). It is natural to conjecture that if $t_2=-t_1=T/2$ with $T\to\infty$, the distribution of $\vec{r}(0)$ under this measure  coincides with the ground state probability distribution  $|\varphi_0(\vec{r})|^2$ from the Born rule. To see that this is the case, consider a path that passes through $(\vec{r}_-,-T/2)$, $(\vec{r},0)$ and $(\vec{r}_+,T/2)$ for some arbitrary initial and final point $\vec{r}_\pm$. The overall propagator is then
\begin{equation}
  K(\vec{r}_+,T/2;\vec{r},0)K(\vec{r},0;\vec{r}_-,-T/2;)\sim  |\varphi_0(\vec{r})|^2\varphi_0(\vec{r}_+)\varphi^*_0(\vec{r}_-)e^{-E_0T}.
\end{equation}
Apart from a normalization factor that depends on $\vec{r}_\pm$ and $T$,\footnote{In particular, we see that $\mathcal{N}\sim e^{E_0 T}$ in the long time limit.} this is just the expected ground state distribution $|\varphi_0(\vec{r})|^2$. Thus, the ability to sample from the FK measure for long trajectories would also allow us to sample from the ground state distribution.

\subsection{Quantum Mechanics and Optimal Control}

Long before Feynman, the path measure was studied for finite $T$ and the case of a free particle ($V(\vec{r})=0$) by \cite{Schrodinger:1931aa,Schrodinger:1932aa} in an early exploration of the connection between his equation and the heat equation. Schr\"odinger sought the path measure that interpolates between given marginal distributions at the initial and final times. The generalizations of this question are now known as the \emph{Schr\"odinger problem}: see \cite{Leonard:2014aa} for an insightful recent review with a historical survey. One of the major conclusions of these works is that the path measure is \emph{Markov} (\cite{Jamison:1974aa}) and that paths $\vec{r}(t)$ satisfy a stochastic differential equation (SDE)
\begin{equation}\label{eq:sde}
  d\vec{r}_t = d\vec{B}_t + \vec{v}(\vec{r}_t,t)dt,
\end{equation}
where $\vec{B}_t\in \mathbb{R}^d$ is a standard Brownian motion and $\vec{v}(\vec{r}_t,t)$ is a drift that is determined by the potential $V(\vec{r})$, as well as the initial and final conditions. The problem of finding the drift can be given an optimal control formulation, with $\vec{v}(\vec{r}_t,t)$ achieving the minimum of the cost function (\cite{Holland:1977aa,Fleming:1977aa})
\begin{equation}\label{eq:hollandcost}
  C_T[\vec{v}] = \frac{1}{T}\E\left[\int_0^T\left[\frac{1}{2}(\vec{v}(\vec{r}_t,t))^2 + V(\vec{r}_t)\right]dt\right],
\end{equation}
where the expectation is over the process \eqref{eq:sde}. When this cost is minimized the path measure of this process coincides with the FK measure. For the infinite horizon $T\to\infty$ case the optimal drift has no explicit time dependence and we obtain the ground state energy as $E_0 = \lim_{T\to\infty} \min_{\vec{v}} C_T[\vec{v}]$.

The existence of a variational principle allows the tools of deep learning to be brought to bear. In this work the drift function will be parameterized by a neural network with parameters $\theta$: $\vec{v}(\vec{r})=\vec{v}_\theta(\vec{r})$. Finding the optimal drift can be regarded as a \emph{reinforcement learning} problem with the cost \eqref{eq:hollandcost} as our (negative) reward function. Although the expectation in \eqref{eq:hollandcost} is intractable, a Monte Carlo estimate may be made by generating a batch of solutions of the SDE using standard discretizations (\cite{Kloeden:2013aa}). Optimization of this estimate is then possible through automatic differentiation of the estimated cost with respect to $\theta$. Derivatives with respect to the expectation are straightforward as the SDE increments are written in terms of standard Brownian increments, analogous to the reparameterization trick used in stochastic backpropagation through deep latent Gaussian models (\cite{Rezende:2014aa}) or variational autoencoders (\cite{Kingma:2013aa}).

As well as providing a route to the evaluation of the ground state energy, \eqref{eq:sde} provides a sampler for the ground state probability distribution, as the stationary distribution of the SDE with the optimal drift coincides with $|\varphi_0|^2$.

The imaginary time Feynman path integral (or FK formula) is the basis of other numerical techniques in physics, notably the path integral Monte Carlo method (\cite{Ceperley:1995aa}). In this approach, the state of the Monte Carlo simulation corresponds to an entire Feynman trajectory, which is updated according to a Markov chain that ensures the paths sample the FK distribution in the stationary state. Our approach by contrast learns to sample the trajectories in the optimal way. We intend our method to be a drop-in replacement for path integral Monte Carlo as far as ground state properties are concerned.

\subsection{The Many-Body Problem}

All of the above considerations extend straightforwardly to systems of many identical particles, with one important caveat. The wavefunction of a system of identical particles must either be completely \emph{symmetric} or \emph{antisymmetric} under exchange of particles, corresponding to identical bosons or fermions (the latter described by the Pauli exclusion principle). The overall ground state of an identical particle Hamiltonian is symmetric under very general conditions and therefore corresponds to the bosonic case (\cite{Feynman:1998aa}). Apart from some two-electron problems where the spin state of the electrons is antisymmetric, allowing the spatial wavefunction to be symmetric (see \Cref{sec:exp}), the many-fermion problem is out of reach of our method for now.\footnote{The same restriction applies to the path integral Monte Carlo method, where it is known as the \emph{sign problem}.}

\subsection{Outline}

The outline of the remainder of this paper is as follows. In \Cref{sec:oc} we discuss the connection between quantum mechanics and optimal control in more detail, focusing on the variational principle that can be used to learn the optimal drift. In \Cref{sec:rl} we introduce our neural representation of the drift, with particular attention given to permutation equivariance which is a feature of systems of indistinguishable particles. We also describe the  algorithm used to learn the drift. \Cref{sec:exp} describes experiments on some simple physical systems. Finally, in \Cref{sec:con} we summarize our findings and provide an outlook to future work.

\section{Quantum Mechanics and Optimal Control}\label{sec:oc}

\subsection{Fokker--Planck to Schr\"odinger}\label{sec:fp}

The most direct way to establish a link between stochastic processes and quantum mechanics is via a standard mapping between the Fokker--Planck and Schr\"odinger equations. This is discussed for example in \cite{Risken:1996aa,Pavliotis:2014aa}, but we repeat the main points here. Starting from a solution $ p({\vec{r}},t)$ of the Fokker--Planck (FP) equation
\begin{equation} \label{eq:fp}
  \frac{\partial p({\vec{r}},t)}{\partial t} = \frac{1}{2} \nabla^2 p({\vec{r}},t) + \nabla \cdot \left ( p({\vec{r}},t) \nabla U({\vec{r}}) \right ),
\end{equation}
with a drift $\vec{v}(\vec{r}) = - \nabla U({\vec{r}})$ given in terms of some potential function $U({\vec{r}})$, the stationary state $\pi({\vec{r}})$ of this FP equation has the form of a Boltzmann distributio
\begin{equation}
  \pi({\vec{r}}) \propto \exp \left(-2 U({\vec{r}}) \right).
\end{equation}
If we define the function
\begin{equation}\label{eq:psi_def}
  \psi({\vec{r}},t) = \frac{p({\vec{r}},t)}{\sqrt{\pi({\vec{r}})}},
\end{equation}
then $\psi({\vec{r}},t)$ satisfies the (imaginary time) Schr\"odinger equation
\begin{equation}
  \frac{\partial \psi({\vec{r}},t)}{\partial t} = -H \psi({\vec{r}},t)
\end{equation}
with Hamiltonian $H=\frac{1}{2} \nabla^2+V_U(\vec{r})$,
where the potential $V_U(\vec{r})$ has the form
\begin{equation}\label{eq:FP-pot}
V_U(\vec{r})\equiv \frac{1}{2} \left[-\nabla^2 U + (\nabla U)^2\right]=\frac{1}{2}\left[\nabla\cdot \vec{v}+\vec{v}^2\right].
\end{equation}
The zero energy ground state wavefunction of this Hamiltonian is
\begin{equation}
  \varphi_0({\vec{r}}) = \sqrt{\pi({\vec{r}})}.
\end{equation}
The Fokker--Planck equation \eqref{eq:fp} describes the evolution of the probability density of a stochastic process described by the stochastic differential equation (SDE)
\begin{equation}\label{eq:SDE2}
  d{\vec{r}}_t = d\vec{B}_t +\vec{v}({{\vec{r}}_t})dt,
\end{equation}
where $\vec{B}_t$ is a standard Brownian motion.
We emphasize that the quantum probability in the ground state $|\varphi_{0}({\vec{r}})|^2$ coincides with the classical stationary distribution $\pi_0$ of this process.

\subsubsection{Example: Calogero--Sutherland and Dyson Brownian Motion}

An instructive example of a many body problem and the associated stochastic process is provided by the Calogero--Sutherland model describing particles in one dimension in a harmonic potential interacting by an inverse square potential (\cite{Sutherland:1972aa})
$$
H = \sum_i \frac{1}{2}\left[-\frac{\partial^2}{\partial x_i^2}+x_i^2\right] + \lambda(\lambda-1)\sum_{i<j} \frac{1}{(x_i-x_j)^2}.
$$
In this case, the ground state is known exactly and has the form of the Gaussian ground state wavefunction of the harmonic oscillator multiplied by a pairwise Jastrow factor
$$
\Phi_0(x_1,\ldots x_N) = \prod_{i<j}|x_i-x_j|^{\lambda}\exp\left(-\frac{1}{2}\sum_i x_i^2\right).
$$
The drift of the associated SDE is
\begin{equation}\label{eq:CS_drift}
v_i =\partial_i \log\Phi_0= - x_i + \lambda \sum_{j\neq i} \frac{1}{x_i-x_j}
\end{equation}
corresponding to Dyson's Brownian motion (\cite{Dyson:1962aa}, see \Cref{fig:cs-sim}).

\begin{figure}
\floatconts
  {fig:cs-sim
  {\caption{Trajectories of 15 particles evolving under the Dyson process \eqref{eq:CS_drift}.}
  {\includegraphics[width=0.5\textwidth]{15_particles.png}}
\end{figure}

\subsection{The Variational Principle}\label{sec:var}

\begin{comment}
We start from the time independent Schr\"odinger equation for the ground state wavefunction $\varphi_0(\vec{r})$
\begin{equation}
  [H\varphi_0](\vec{r}) =\left[-\frac{1}{2}\nabla^2+V(\vec{r})\right]\varphi_0(\vec{r}) = E_0 \varphi_0(\vec{r})
\end{equation}
Since the ground state wavefunction is nodeless it may be chosen real and positive, and hence expressed as
\begin{equation}
  \varphi_0(\vec{r}) = \exp\left[-U(\vec{r})\right],
\end{equation}
where $U(\vec{r})$ satisfies
\begin{equation}
  \frac{1}{2}\nabla^2 U -\frac{1}{2}(\nabla U)^2 + V(\vec{r}) = E_0.
\end{equation}
We can write this as a Hamilton--Jacobi--Bellman (HJB) equation
\begin{equation}\label{eq:aux}
  \frac{1}{2}\nabla^2 U+\min_{\vec{v}}\left[{\vec{v}}\cdot\nabla U + \frac{1}{2}|{\vec{v}}|^2\right] + V(\vec{r}) = E_0,
\end{equation}
where the minimum is taken over all drifts $\vec{v}(\vec{r})$. The HJB equation may be interpreted in terms of the SDE
\begin{equation}\label{eq:holland_sde}
  d\vec{r}_t = d\vec{B}_t + \vec{v}(\vec{r}_t)dt,
\end{equation}
with $\vec{r}_0 = \vec{r}$. Define the cost function
\begin{equation}\label{eq:hollandcost}
  C[\vec{v}] = \lim_{T\to\infty}\E\left[\frac{1}{T} \int_0^T dt~\left(\frac{1}{2}(\vec{v}(\vec{r}_t))^2 + V(\vec{r}_t)\right)\right].
\end{equation}
\cite{Holland:1977aa} shows that the minimum of the functional $C[{\vec{v}}]$ is  $E_0$, which is achieved for
\begin{equation}
  {\vec{v}}_0 = \frac{\nabla \varphi_0}{\varphi_0} = -\nabla U.
\end{equation}
Although Holland's analysis provides a cost function that will be basis of our approach, the above formal manipulations don't provide much intuition.
\end{comment}

In order to turn the above connection into a calculational tool, we exploit a connection between the path measure $\mathbb{P}_{\vec{v}}$ of the SDE \eqref{eq:SDE2} and the corresponding quantity $\mathbb{P}_\text{FK}$ from the FK formula. For a FP measure that relates the ground state distribution $\pi_0(\vec{r}_{0,T})$ at $t=0$ and $t=T$, the log likelihood ratio (log Radon--Nikodym derivative) of the two measures obeys
\begin{align}\label{eq:KLlim}
  \log\left(\frac{d\mathbb{P}_{\vec{v}}}{d\mathbb{P}_\text{FK}}\right) &=\ell_T - E_0 T+\log\left(\frac{\varphi_0(\vec{r}_0)}{\varphi_0(\vec{r}_T)}\right)\nonumber\\
   \ell_T&\equiv \int \vec{v}(\vec{r}_t)
  \cdot d\vec{B}_t+\int dt\left(\frac{1}{2}|\vec{v}(\vec{r}_t)|^2+V(\vec{r}_t)\right),
\end{align}
See Appendix \ref{sec:prob-proof} for a heuristic proof. The  Kullback–Leibler (KL) divergence of the two measures is then
\begin{equation}\label{eq:full-KL}
    D_{\text{KL}}\left(\mathbb{P}_{\vec{v}} \middle\| \mathbb{P}_\text{FK}\right)=\E_{\mathbb{P}_{\vec{v}}}\left[\ell_T-E_0 T+\log\left(\frac{\varphi_0(\vec{r}_0)}{\varphi_0(\vec{r}_T)}\right)\right].
\end{equation}
Note that \eqref{eq:full-KL} is true for any $\vec{r}_0$. If we additionally average over $\vec{r}_0$ drawn from the stationary distribution of the SDE, then the distributions of $\vec{r}_T$ and $\vec{r}_0$ coincide, and the final term vanishes. Because the KL divergence is positive $D_{\text{KL}}\left(\mathbb{P}_{\vec{v}} \middle\| \mathbb{P}_\text{FK}\right)\geq 0$ we then have $\E_{\mathbb{P}_{\vec{v}}}\left[\ell_T\right]\geq E_0 T$, with equality when the two measures match. This is the variational principle that forms the basis of our approach. For $T\to\infty$ we can ignore the final term so that
\begin{equation}
    \lim_{T\to\infty} \frac{\E_{\mathbb{P}_{\vec{v}}}\left[\ell_T\right]}{T}\geq E_0,
\end{equation}
irrespective of the initial state distribution (as $T\to\infty$ the final state distribution will be the stationary state of the SDE, assuming ergodicity). Using the fact that Brownian increments have zero expectation recovers Holland's cost \eqref{eq:hollandcost}.

\section{Reinforcement Learning}\label{sec:rl}

We have identified a variational principle for the drift with a time integrated cost function. Learning the optimal drift may be regarded as an entropy regularized reinforcement learning (RL) problem. The usual Markov decision process formulation of RL involves Markov dynamics in a state space together with a set of actions that can influence the dynamics within this space. We seek a \emph{policy} -- a distribution of actions to be taken at each step given the current state -- optimal with respect to a time integrated reward function. In the present case, optimizing over all possible drift functions corresponds to complete controllability of the state, instead of a restriction to a (possibly finite) set of actions. As we have seen in \Cref{sec:var}, however, the optimization problem is regularized by a KL divergence that originates from the kinetic energy, and penalizes deviations of the stochastic process from the `passive' Brownian dynamics. In RL the resulting cost gives rise to \emph{linearly solvable Markov decision problems} (\cite{Todorov:2007aa,Todorov:2008aa}). \footnote{In RL this entropy regularization has a practical utility in encouraging diverse policies (\cite{Levine:2018aa}).}

In the remainder of this section we introduce a concrete algorithm to learn the optimal drift for a given quantum many-body problem.

\subsection{Neural representation of the drift}\label{sec:drift}

In practice, one cannot optimize over arbitrary drifts. We introduce a neural parameterization $\vec{v}_\theta(\vec{r})$ of the drift function, so the quality of our approximation is limited by the expressiveness of the network. Thus we are performing a variational approximation to the optimal control problem, with the cost bounded from below by the ground state energy.

For $N>1$ particles, the drift of each is a function of all the particle coordinates: $\vec{v}_{i,\theta}(\vec{r}_1,\ldots,\vec{r}_N)$ $i=1,\ldots, N$. For identical particles, this drift must be \emph{permutation equivariant}, meaning that for any permutation $P$ of the particles
\begin{equation}
    \vec{v}_{i,\theta}(\vec{r}_1,\ldots,\vec{r}_N) = \vec{v}_{P(i),\theta}(\vec{r}_{P(1)},\ldots,\vec{r}_{P(N)}).
\end{equation}
We now introduce two ways in which equivariance may be achieved.

\subsubsection{DeepSets}\label{sec:deepsets}

\cite{Guttenberg:2016aa,Zaheer:2017aa} identified deep neural architectures that respect permutation equivariance, building on the much earlier work of \cite{Shawe-Taylor:1989aa}. We use a variant of the \textsc{DeepSets} architecture of \cite{Zaheer:2017aa}, in which an equivariant hidden state is based on the sum-pooling operation
$$
\mathbf{h}_i = \phi\left(\boldsymbol{\Gamma}\vec{r}_i + \boldsymbol{\Lambda}\sum_j\vec{r}_j\right)
$$
(biases omitted for clarity) for weight matrices $\boldsymbol{\Gamma}, \boldsymbol{\Lambda}\in \mathbb{R}^{H\times d}$, where $d$ is the spatial dimension and $H$ the number of hidden units. $\phi$ is a nonlinear activation function. These layers can then be stacked to produce a more expressive representation.

\subsubsection{PairDrift}\label{sec:pairdrift}

Alternatively, we may use a representation of single particle and pair features (c.f. \cite{Pfau:2019aa})
\begin{align}
\mathbf{h}_i &= \boldsymbol{\sigma}_1(\mathbf{r}_i) + \sum_j \boldsymbol{\pi}_1(\mathbf{r}_i-\mathbf{r}_j)\\
\mathbf{h}_{ij} &= \boldsymbol{\Pi}_1(\mathbf{r}_i-\mathbf{r}_j).
\end{align}
where $\boldsymbol{\sigma}, \boldsymbol{\pi}:\mathbb{R}^d\to \mathbb{R}^H$ and $\boldsymbol{\Pi}:\mathbb{R}^d\to \mathbb{R}^{H\times H}$ are neural networks. The single particle hidden state $\mathbf{h}_i$ contains information about single particle and pair locations. This process may be continued
\begin{align}
\tilde{\mathbf{h}}_i &= \boldsymbol{\sigma}_2(\mathbf{h}_i) + \sum_j \boldsymbol{\pi}_2(\mathbf{h}_{ij})\\
\tilde{\mathbf{h}}_{ij} &= \boldsymbol{\Pi}_2(\mathbf{h}_{ij}).
\end{align}
Finally, the permutation equivariant drift function is written in terms of the last hidden state as
\begin{align}
\vec{v}_i &= \boldsymbol{\sigma}_3(\tilde{\mathbf{h}}_i) + \sum_j \boldsymbol{\pi}_3(\tilde{\mathbf{h}}_{ij}).
\end{align}
Note that if we had represented the drift simply as
$$
\vec{v}_i &= \boldsymbol{\sigma}(\mathbf{r}_i) + \sum_j \boldsymbol{\pi}(\mathbf{r}_i-\mathbf{r}_j),
$$
this would correspond to a pairwise Jastrow factor in the many-body wavefunction via the formula $\vec{v}_i=\nabla_i \log\Psi$. Repeatedly merging the pair features into the single particle features enriches the representation. In our experiments, we find it sufficient to include a single hidden layer.

\subsection{Integration of the SDE} \label{sec:int}

From \eqref{eq:KLlim}, we see that to evaluate the cost function we need to solve the system of SDEs
\begin{align}\label{eq:sde_system}
  d\vec{r}_t &= d\vec{B}_t + \vec{v}(\vec{r}_t)dt,\nonumber\\
  d\ell_t &= \vec{v}(\vec{r}_t)\cdot d\vec{B}_t + \left(\frac{1}{2}|\vec{v}(\vec{r}_t)|^2 + V(\vec{r}_t)\right)dt.
\end{align}
As $T\to\infty$, $\ell_T$ satisfies $\E_{\mathbb{P}_{\vec{v}}}\left[\ell_T\right]\geq E_0 T$. Numerous techniques exist for the numerical solution of such problems, beginning with the Euler--Maruyama scheme
\begin{align}\label{eq:em}
    \vec{r}_{t+\Delta t} &= \vec{r}_{t} + \Delta\vec{B}_t + \vec{v}(\vec{r}_t)\Delta t,\nonumber\\
    \ell_{t+\Delta t} &= \ell_{t} + \vec{v}(\vec{r}_t)\cdot \Delta\vec{B}_t +\left(\frac{1}{2}|\vec{v}(\vec{r}_t)|^2 + V(\vec{r}_t)\right)\Delta t.
\end{align}
where $\Delta\vec{B}_t\sim \mathcal{N}(0,t)$ are iid Gaussian increments.

Higher order methods offering improved accuracy with longer time steps $\Delta t$ are described in \cite{Kloeden:2013aa}. For most of our experiments we used the SOSRA method from \cite{Rackauckas:2018aa}. Implementing these methods in a framework that supports automatic differentiation allows the cost to be optimized by gradient descent with respect to the parameters $\theta$ describing the drift.

As emphasized in recent work on neural ordinary differential equations \cite{Chen:2018aa}, one can regard integration schemes like \eqref{eq:em} as a (recurrent) residual network. This affords a clean separation of the model (the drift $\vec{v}_\theta(\vec{r})$) from the implementation and allows
us to adapt the integration scheme on the fly for different epochs of training, or at evaluation time.

We evolve a batch of trajectories for $T/\Delta t$ steps of size $\Delta t$, starting from an initial state equal to the final state of the previous batch. This is analogous to training a recurrent network by backpropagation through time on infinite sequences by unrolling for a finite number of steps, stopping the gradients from propagating into the infinite past. Starting each step from the final state of the previous step means that the batch tracks closely the stationary distribution of the stochastic process, so no burn-in is required before the integration of the cost function can begin.

This procedure leads to the following issue, arising from the presence of the boundary term in \eqref{eq:KLlim}. Its expectation is zero for a stationary process where the marginal distributions of the initial and terminal points coincide. The gradients of this expectation with respect to the drift parameters are however nonzero, as changing these parameters changes the endpoint $\vec{r}_T$. Evaluating this gradient to get the gradient of the KL in \eqref{eq:full-KL} is intractable because the ground state wavefunction is not known, but it may be estimated from the current drift as
\begin{equation}\label{drift-correct}
\nabla_\theta \log\left(\frac{\varphi_0(\vec{r}_0)}{\varphi_0(\vec{r}_T)}\right) = -\nabla_\theta \vec{r}_T \cdot \vec{v}(\vec{r}_T)\approx -\nabla_\theta \vec{r}_T \cdot \vec{v}_\theta(\vec{r}_T),
\end{equation}
an approximation that improves as the true drift $\vec{v}(\vec{r})=\nabla\log\varphi_0(\vec{r})$ is approached. Although this correction may be neglected in the $T\to\infty$ limit, we find it significantly improves convergence in our experiments.

\subsection{Stochastic Backpropagation}

Recall that the cost function \eqref{eq:full-KL} is the expectation of a time integral over the stochastic process. Although the expectation is intractable, a batch of $B$ simulated trajectories $\vec{r}^{(b)}_{t}$ can be used to make a Monte Carlo estimate as follows
\begin{equation}\label{eq:cost-est}
    \ell_T[\vec{v}_\theta] \approx \frac{1}{B T} \sum_{b,t}\left[\vec{v}_\theta\left(\vec{r}^{(b)}_t\right)\cdot \Delta\vec{B}^{(b)}_t + \left(\frac{1}{2}\left(\vec{v}_\theta\left(\vec{r}^{(b)}_t\right)\right)^2 + V\left(\vec{r}^{(b)}_t\right) \right)\Delta t \right ].
\end{equation}
Because $\vec{r}^{b}_{t}$ have been generated by an SDE discretization such as \eqref{eq:em}, backpropagating through this estimated expectation is analogous to the reparameterization trick in deep latent Gaussian models (\cite{Rezende:2014aa}) or variational autoencoders (\cite{Kingma:2013aa}). This similarity has been emphasized in recent work on neural stochastic differential equations (\cite{Tzen:2019aa}). The memory cost of backpropagation is proportional to $BT/\Delta t$, and is the main bottleneck in our approach.

Though the exact $T\to\infty$ value of the cost is an upper bound on the ground state energy, we have to bear in mind that the approximations of finite time, finite batch size, and any discretization error in the integration of the SDE mean that in practice the upper bound may be violated during training, particularly as optimization converges.

\subsection{Related Work}

We have already mentioned wavefunction based neural approaches to quantum mechanics in \Cref{sec:dlq}. Here we summarize work that is technically related to ours.

Learning SDEs using stochastic optimization was introduced in \cite{Ryder:2018aa} to infer SDE parameters from trajectories and in \cite{Tzen:2019aa} as a generative model to map a simple initial distribution of $\vec{r}_{-T/2}$ to a given empirical distribution of $\vec{r}_{T/2}$.

The path integral point of view has been applied to RL in \cite{Kappen:2007aa,Theodorou:2010aa}. Also from the RL side, \cite{Todorov:2007aa} introduced \emph{linearly solvable Markov decision problems}. Their reward function -- as in the quantum case -- is a sum of a cost on the state space and the KL divergence between the passive and controlled dynamics. Todorov showed that for these problems the nonlinear Bellman equation could be transformed to a linear equation by an exponential transformation. As we have seen in \Cref{sec:oc}, for the quantum cost this linear equation is just the Schr\"odinger equation.

Finally, \cite{Han:2018aa} solved semilinear parabolic PDEs in high dimension using a representation in terms of a backward SDE with a neural representation, but with a cost function based on matching to a given terminal condition.

\section{Experiments} \label{sec:exp}

Our experiments use \textsc{PairDrift} architectures described in \Cref{sec:drift}, which we found to perform best, with a single hidden layer and HardTanh activation. To ensure the trajectories remain close to the origin from the outset, these networks are initialized to give zero drift, but we add a skip connection to give a restoring drift of simple form. For example, a linear skip connection $\mathbf{W}_\text{skip}\vec{r}_i$ to the drift of the $i^\text{th}$ particle ensures that the SDE is initially a 3D Ornstein--Uhlenbeck process.

Optimization of the estimated cost \eqref{eq:cost-est} uses Adam (\cite{Kingma:2015ab}). The gradient is corrected to account for the boundary terms as described in \Cref{sec:int}
\begin{equation}
    \nabla_\theta \ell_T \to \nabla_\theta \ell_T-\nabla_\theta \vec{r}_T \cdot \vec{v}_\theta(\vec{r}_T)
\end{equation}
All models were implemented in PyTorch. For hyperparameters see Appendix \ref{sec:hyperparameters}. The code is available at \url{https://github.com/AustenLamacraft/QuaRL}

\subsection{Hydrogen Atom}

Hydrogen is an electrically neutral atom containing a nucleus consisting of a single proton, and a single electron which is bound to the nucleus via the Coulomb force. Consider a single non-relativistic electron subject to an isotropic Coulomb interaction force provided by the nucleus. In atomic units\footnote{\emph{Atomic units} denotes energies measured in Hartrees (27.211 eV) and distances measured in Bohr radii (0.529 $\AA$).} the explicit form of the Hamiltonian reads
\begin{equation}\label{eq:hydrogen_Hamiltonian}
H = -\frac{1}{2}\nabla^{2} - \frac{1}{|\vec{r}|}
\end{equation}
where $\vec{r}$ is the position of the electron with respect to the center of the nucleus. The Schr\"odinger equation for hydrogen can be solved explicitly via separation of variables, yielding the ground state wave-function
\begin{equation}\label{eq:hydrogen_wave}
\varphi_{0}(\vec{r}) = \frac{1}{\sqrt{\pi}}e^{-|\vec{r}|}.
\end{equation}
As the ground state of hydrogen has an exact analytic form, it is also possible to derive an exact expression for the electron's drift. The ground state drift velocity of the electron is given by
\begin{equation}\label{eq:hydrogen_drift}
\vec{v} = \nabla\log \varphi_0(\vec{r}) = -\hat{\vec{r}}.
\end{equation}
The exact ground state energy of Hydrogen is $E_{0}= -0.5$. We used a skip connection with $\mathbf{W}_\text{skip}=-\mathsf{1}_3$, and found an energy of $-0.497$, which is within $0.6 \%$  of the exact value. The learned drift is illustrated in \Cref{fig:h-and-he}. From \eqref{eq:hydrogen_drift} and Figure \ref{fig:h-and-he}, we see that the Hydrogen ground state drift has a cusp at the position of the proton.

\begin{figure}[htbp]
\floatconts
  {fig:h-and-he
  {\caption{(Left) The learned drift field for the Hydrogen atom in the $x-y$ plane, colorscaled by its magnitude. (Right) The electron density distribution in the stretched $\text{H}_2$ molecule at proton separation $R=2.8$.}
  {\includegraphics[height=0.43\textwidth]{hydrogen_vector_field.jpg}\hfill\includegraphics[height=0.43\textwidth]{H2_density.jpg}}
\end{figure}

\subsection{Helium atom and Hydrogen molecule}
The Helium atom consists of two electrons bound to a nucleus of charge $+2e$. The non-relativistic two-body Hamiltonian reads
\begin{equation}\label{eq:helium_atom_Hamiltonian} H = - \frac{1}{2}    \left(\nabla_{\vec{r}_1}^{2} + \nabla_{\vec{r}_2}^{2} \right ) -  \frac{2}{|\vec{r}_1|} - \frac{2}{|\vec{r}_2|} + \frac{1}{|\vec{r}_{1} - \vec{r}_{2}|}
\end{equation}
where $\vec{r}_{1}$ and $\vec{r}_{2}$ denote the positions of the electrons with respect to the  nucleus. An analytic solution to the Schr\"odinger equation does not exist for Helium. Accurate numerical calculations may be found in the recent paper \cite{Baskerville:2019aa}, together with a comparison with the Hartree--Fock method. In the Hartree--Fock (HF) approximation the wavefunction is taken as the product of two single particle wave-functions \footnote{The antisymmetry of the wavefunction demanded by the Pauli principle is achieved for the ground state by putting the spin states of the two electrons in a spin singlet. Spin triplet states -- requiring an antisymmetric spatial wavefunction -- are not accessible by our method.}
\begin{equation}\label{eq:helium_wave}
\Psi_{\text{HF}}\left(\vec{r}_1 , \vec{r}_2\right) = \psi_{\text{HF}}\left(r_1\right)\psi_{\text{HF}}\left(r_2\right),
\end{equation}
where the single particle wavefunction $\psi_\text{HF}(r)$ is optimized to minimize the energy. \cite{Baskerville:2019aa} find a correlation energy -- the difference between the best HF value and the numerically exact values -- of $0.042$, or $1.4 \%$ of the exact value.

We used the \textsc{PairDrift} network to model the Helium atom. The ground state energy was found to be within $0.2\%$ of the exact value.

For both the Helium atom and the Hydrogen molecule, we used a skip connection taking advantage of Kato's cusp conditions (\cite{Kato:1957aa}). These are conditions on the cusps that appear in the wavefunction whenever two particles with a Coulomb interaction are at the same position. Kato's cusp conditions for wave-functions translate directly into cusp conditions for drifts: they should be proportional to $\hat{\vec{r}}$ in the limit of zero separation. For electron-electron interactions, the proportionality factor is $\frac{1}{2}$ and the drift is repulsive, while for electron-nucleus interactions, the proportionality factor is the charge of the nucleus and the drift is attractive. For our experiments on the Helium atom and Hydrogen molecule, we therefore use a skip connection that is a sum of terms $\hat{\vec{r}}$, i.e. one term for each Coulomb interaction.

Although we could have incorporated other conditions arising from rotational symmetries, we chose not to do so because these symmetries do not generalize to other continuum systems.

The Hydrogen molecule consists of two electrons orbiting a pair of bound protons at separation $R$, which repel each other. With the protons at positions $\pm R\hat{\vec{z}}/2$, the Hamiltonian becomes
\begin{equation}\label{eq:hydrogenMol_Hamiltonian}
H = - \frac{1}{2}   \left(\nabla_{\vec{r}_1}^{2} + \nabla_{\vec{r}_2}^{2} \right)  + \frac{1}{R} + \frac{1}{|\vec{r}_{1} - \vec{r}_{2}|} -  \frac{1}{|\vec{r}_{1} - R\hat{\vec{z}}/2|} -  \frac{1}{|\vec{r}_{1} + R\hat{\vec{z}}/2|}
-  \frac{1}{|\vec{r}_{2} - R\hat{\vec{z}}/2|} -  \frac{1}{|\vec{r}_{2} + R\hat{\vec{z}}/2|}.
\end{equation}
The correlation energy of the Hydrogen molecule is -0.044 Hartrees, corresponding to 3.8\% of the numerically exact value.

To model the Hydrogen molecule, we also used the  \textsc{PairDrift} network with the skip connection describing the electron-electron and electron-nucleus cusps. We find a ground state energy within $0.3 \%$ of the exact value at the equilibrium separation of $R=1.401$. For a stretched Hydrogen molecule at $R=2.800$ (i.e. twice the equilibrium separation), our result is similarly $0.3 \%$ from the exact value.

All our results for the atomic and molecular systems are given in \Cref{tab:results}.

\subsection{Bosons in a Harmonic Trap}

As an example of a many-particle system we study a system of $N$ identical bosons trapped in a two-dimensional isotropic harmonic potential, with the atom-atom interactions approximated by a finite-sized Gaussian contact potential (\cite{Mujal:2017aa}).  The relevant Hamiltonian reads
\begin{equation}\label{eq:harmonic_Hamiltonian}
H =  \sum_{i = 1}^{N} \left[ - \frac{1}{2}  \nabla_{i}^{2} + \frac{1}{2}\vec{\vec{r}_{i}}^{2} \right] +  \frac{g}{\pi s^{2}}\sum_{i < j }^{N}\exp\left[-\left(\vec{r}_{i} - \vec{r}_{j}\right)^2/s^2\right]
\end{equation}
where $g$ and $s$ are respectively the strength and range of the pairwise interaction. \break

We simulated the ground state for $N = 2,3,$ and $4$ bosons for $s=0.5$ and $g=0-15$ using the \textsc{PairDrift} neural architecture. The energies derived via the neural network are compared with those of \cite{Mujal:2017aa} (obtained by exact diagonalization) in \Cref{fig:boson-drift-and-results}. An energy correspondence within 1\% of the exact diagonalization results was observed for all parameter combinations.
\break

\begin{figure}[htbp]
\floatconts
  {fig:boson-drift-and-results

 \includegraphics[height=0.32\textwidth]{bosons_vector_field.jpg}\hfill\includegraphics[height=0.32\textwidth]{boson_energies.jpeg}

 {\caption{(Left) The learned drift field for interaction strength $g=15$ of one boson (orange) in the presence of the others (green), together with the restoring drift arising from the harmonic potential pushing particles to the center (cross). Note the resultant drift field tends to keep the bosons apart. (Right) Comparison of ground state energies of $N=2$, $N=3$, and $N=4$ bosons (blue, orange and green, respectively) with results of \cite{Mujal:2017aa}.}

\end{figure}

\begin{table}[]
    \centering
    \begin{tabular}{ p{3cm} | p{2.3cm} p{2.3cm} p{2.3cm} p{2.3cm}}
              & H atom    & He atom   & $\text{H}_2$ molecule & $\text{H}_2$ molecule ($R=2.8$)\\
 \hline
  Numerically exact  & -0.5      & -2.903    &  -1.173   & -1.071\\
 Hartree-Fock   & N/A        & -2.862 (1.4\%)   &  -1.129 (3.8\%)   & \\
 \textbf{Ours}  &-0.497 (0.6\%)  & -2.898 (0.2\%) &  -1.169 (0.3\%)    & -1.068 (0.3\%)
\end{tabular}
    \caption{Comparison of our approach with Hartree-Fock and numerically exact values for the ground state energies. The percentages within parentheses are the relative errors of an energy estimate compared to the numerically exact value.}
    \label{tab:results}
\end{table}

\hfill \break

\section{Conclusion}\label{sec:con}

In this work we have presented a novel neural approach to many-body quantum mechanics that leverages the optimal control formulation of \cite{Holland:1977aa} to approximate the Feynman--Kac path measure. This contrasts with earlier deep learning approaches to the many-body problem that start from the Schr\"odinger picture. We have demonstrated the utility of our approach for some simple many-body problems.

Future work will extend this approach in the following directions:

\begin{enumerate}
    \item Many models studied in condensed matter physics are defined on a \emph{lattice} rather than in continuous space. In this case identical particles are handled by working with \emph{occupation numbers} -- the count of particles on each lattice site -- rather than by introducing a fictitious index for each particle and then demanding invariance under labelling. The control theoretic perspective is well suited to such discrete state spaces as described in \cite{Todorov:2007aa}.
    \item Rotating systems or systems with finite angular momentum may be handled with a modified cost function that
    leads to a rotational component to the drift. This is how time reversal symmetry is broken in our formalism, corresponding to complex wavefunctions in the Schr\"odinger picture.

    \item The memory cost of backpropagating through SDE trajectories is a major limitation of our approach that could be addressed using checkpointing  \cite{Martens:2012aa,Chen:2016aa,Gruslys:2016aa}. \cite{Chen:2018aa} suggested that neural ODE models be trained at constant memory cost by solving the linear adjoint equation backwards in time. The difficulties associated with implementing this method in the SDE case are discussed in \cite{Tzen:2019aa}, though very recent work has addressed this issue \cite{li2020scalable}.

\end{enumerate}

\acks{We thank Alex Gaunt for useful discussions. We acknowledge support from a National Science Foundation Graduate Research Fellowship under Grant No. DGE-174530 (ARB) and EPSRC Grant No. EP/P034616/1 (AL).}

\section{Introduction}

Finding the ground state of a quantum mechanical system involves finding the lowest eigenvalue of the Hamiltonian operator that describes the system. For many body systems, this is a problem that grows exponentially harder as the number of particles increases, and finding tractable approaches has been the principal computational challenge of quantum mechanics since its inception.

When studying a non-relativistic system of electrons and nuclei, a natural starting point is a realistic Hamiltonian in continuous space with Coulomb interactions. Often, however, \emph{model Hamiltonians} are used instead: for example, the singularity of the Coulomb potential describing the electron-ion attraction is often replaced with a less singular \emph{pseudopotential} in quantum chemistry calculations. Another type of simplification, more common in strongly correlated physics, is to study \emph{lattice models} where particle locations are discretized. Although less realistic, this simplification often allows the treatment of larger systems where qualitative different thermodynamic phases become apparent.

In recent years, progress in deep learning has opened up new fronts to attack the many body problem by parameterizing the wavefunction describing the ground state using neural networks (e.g.\ \cite{carleo_solving_2017,choo_two-dimensional_2019, pfau_ab-initio_2019}). In this work we introduce a new neural approach to lattice models using an alternative formulation based on Reinforcement Learning (RL), in which the optimal policy matches a stochastic process describing the ground state probability distribution. \cite{barr_quantum_2020} applied the same perspective to continuum models using the Feynman--Kac formula and tools of stochastic calculus. We will present the analogous theory for the lattice case, and introduce algorithms to learn the optimal policy in this setting.

Our formulation draws together concepts from many-body physics, control theory, and reinforcement learning. The next Section introduces the necessary background, before \Cref{sec:rl4qm} presents the RL formulation of quantum states. \Cref{sec:neural} describes the application of this formulation to the neural representation of quantum states, before we present our conclusions in \Cref{sec:conclusion}.

\section{Theoretical Background}

In this Section we will introduce the models to be studied and two stochastic representations: in continuous time using a generalized Feynman--Kac formula, and in discrete time using the Schr\"odinger equation directly. We then describe the class of linearly solvable Markov decision problems introduced in \cite{todorov2006linearly}, which form the basis of our mapping to RL. The connections between RL and quantum physics are well established in the case of continuous state spaces: see \cite{barr_quantum_2020}
 and references therein.

\subsection{Models}

For illustration purposes, we restrict ourselves to the simplest family of lattice models describing a system of $N$ spins taking values $Z_i=\pm 1$. \footnote{The term \emph{lattice gas} is sometimes used, where the interpretation is that sites in the lattice are either occupied or empty.} The Hamiltonians act on the $2^N$-dimensional space $\mathbb{C}^2\otimes \mathbb{C}^2 \otimes \cdots \mathbb{C}^2$ and are expressed in terms of Pauli spin matrices $\{X_i, Y_i, Z_i\}$
\begin{equation}
    X_i = \begin{pmatrix}
    0 & 1 \\
    1 & 0
    \end{pmatrix}_i\qquad
    Y_i = \begin{pmatrix}
    0 & -i \\
    i & 0
    \end{pmatrix}_i\qquad
    Z_i = \begin{pmatrix}
    1 & 0 \\
    0 & -1
    \end{pmatrix}_i
\end{equation}
acting on spin $i$. We consider two examples, the \textbf{transverse field Ising model}
\begin{equation} \label{eq:ising}
      H_\text{Ising} = -J \sum_{\braket{i,j}} Z_{i} Z_{j}
        - h \sum_{i} X_i,
\end{equation}
and the \textbf{XXZ model}
\begin{equation} \label{eq:xxz}
      H_{XXZ} = -\sum_{\braket{i,j}} \left[J Z_{i} Z_{j} + J_\perp\left(X_iX_j + Y_iY_j\right)\right],
\end{equation}
where $\braket{i,j}$ denotes a sum over nearest neighbors on the lattice. The case $J=J_\perp$ is known as the XXX or Heisenberg model \citep{auerbach2012interacting}.

\subsection{Connection to stochastic dynamics}

Particular limits of the above models can be described in terms of stochastic dynamics. We illustrate this connection in the case of the XXX model by observing that in the space of two neighboring spins
\begin{equation}
    \left[Z_{i} Z_{j} + X_iX_j + Y_iY_j\right]-1= \begin{pmatrix}
    0 & 0 & 0 & 0 \\
    0 & -1 & 1 & 0 \\
    0 & 1 & -1 & 0 \\
    0 & 0 & 0 & 0
    \end{pmatrix}_{{ij} }.
\end{equation}
This may be interpreted as a transition rate matrix for a Markov process with transitions $\circ \bullet\longleftrightarrow \bullet \circ$, where filled and empty circles represent spin up and spin down states. Apart from an additive constant, the operator $-H_\text{XXX}$ for $J>0$ is then a transition rate matrix of the \textbf{symmetric simple exclusion process}, describing a lattice gas of diffusing particles (\Cref{fig:isingandxy}).
\begin{equation}
-H_\text{XXX} = {\overbrace{J\sum_{\braket{i,j}}\left[Z_{i} Z_{j} + X_iX_j + Y_iY_j-1\right]}^{\equiv\Gamma_{\text{SEP}}}} + JN.
\end{equation}
A valid transition matrix has columns summing to zero; by symmetry of the Hamiltonian this is also true for the rows. The ground state of $H_\text{XXX}$ corresponds to the dominant (maximal) eigenvector of the transition matrix with eigenvalue 0. This eigenvector describes a uniform distribution of particle configurations.

When $J\neq J_\perp$ the Hamiltonian no longer has a simple stochastic interpretation, having the form
\begin{equation}\label{eq:stoq}
 H_\text{XXZ} = -\Gamma_{\text{SEP}} + V,
\end{equation}
where $V$ is diagonal but not constant. Hamiltonians of this form are called \emph{stoquastic}, and we will refer to the Markov process defined by $\Gamma$ as the \emph{passive dynamics}, for reasons that will become clear. In terms of the matrix elements of the Hamiltonian, the potential is $V_{ss'}=V(s)\delta_{ss'}$, where
\begin{equation} \label{eq.potential.from.H}
  V(\vec{s}) = H_{\vec{s}\bs} + \sum_{\vec{s}' \neq \vec{s}} H_{\vec{s} \vec{s}'}.
\end{equation}
Moreover, note that \eqref{eq:stoq} implies that $H_{\vec{s}\bs'}\leq0$ for $\vec{s}\neq\vec{s}'$.

As a second example, consider $H_\text{Ising}$, which is stoquastic for $h>0$ with very simple transition rates
\begin{equation}
 \Gamma_{{\vec{s}}\to{\vec{s}'}} = h,
\end{equation}
if the two spin configurations $s$ and $s'$ differ at only one site (\Cref{fig:isingandxy}). From \Cref{eq.potential.from.H}, the potential is
\begin{align*}
V(\vec{s})
    &= -\sum_{{\vec{s}'}\neq{\vec{s}}} \Gamma_{{\vec{s}}\to{\vec{s}'}} + H_{ss}\\
    &= -hN -  J \sum_{\braket{i,j}} {s_i} {s_j}.
\end{align*}
\begin{figure}
    \centering
    \includegraphics[height=0.43\textwidth]{ising_and_xy.png}
    \caption{Transitions in the simple exclusion process (above) and the Ising model (below). In the simple exclusion process, a transition is a move of a particle to a non-occupied (white) lattice site. For the Ising model, a transition is a flip of spin up (black) to spin down (white) or the other way around.}
    \label{fig:isingandxy}
\end{figure}

A key property of stoquastic Hamiltonians is the following: if the passive dynamics is ergodic, then the ground state wavefunction is positive \citep{bravyi_monte_2015}. This is a consequence of the Perron-Frobenius theorem, and explains why stoquastic Hamiltonians are `free of the sign problem': the ground state wavefunction does not change sign.

\subsection{Feynman--Kac formula}\label{sec:fk}

For a stoquastic Hamiltonian, a representation of the ground state in terms of the passive dynamics exists, in the form of the Feynman--Kac formula. Although this is more familiar in the context of continuous state spaces where the underlying stochastic dynamics is Brownian motion, we will require a discrete analog.

First, note that the ground state $\varphi_0$ of the Hamiltonian is the asymptotic solution of the imaginary time Schr\"odinger equation
\begin{equation} \label{eq.imaginary.schrodinger}
    \partial_t \varphi(t) = - H \varphi(t)
\end{equation}
$\varphi(t)\to e^{-E_{0}  t}\varphi_0$ as $t\to\infty$, where $E_0$ is the ground state energy \citep{goldberg_integration_1967}.

The Feynman--Kac formula gives the solution of the imaginary time Schr\"odinger equation with a \emph{terminal} condition at some time $t+T$.
For a stoquastic Hamiltonian with passive rates $\Gamma$ and potential $V$, the Feynman--Kac formula states that the solution at the earlier time $t$ can be expressed as an expectation \citep{rogers_diffusions_2000}
\begin{equation}
 \varphi({\vec{s}_t}, {t}) = \E_{\mathbb{P}}
 \left[\exp\left(-\int_{t}^{t+T} V({\vec{s}_{t'}}) dt'\right)\varphi({\vec{s}_{t+T}}, {t+T})\right].
\end{equation}
Here $\varphi({\vec{s}_{t+T}}, {t+T})$ is the terminal condition, and the expectation is over trajectories following the passive dynamics $\Gamma$ starting at $s_t$ at time $t$.
Since $\varphi({\vec{s}_t}, {t}) = e^{-E_0 t}\varphi_0({\vec{s}_t})$ is an exact solution of the imaginary time Schr\"odinger equation, the Feynman--Kac formula implies
\begin{equation} \label{eq.fk.ground}
\varphi_0(\vec{s}_t) = \E_{\mathbb{P}}
 \left[\exp\left(-\int_{t}^{t+T} \left (V({\vec{s}_{t'}}) {+} E_0 \right ) dt' \right)\varphi_0({\vec{s}_{t+T}})\right].
\end{equation}
In \eqref{eq.fk.ground} the role of the potential is to modify the measure over trajectories relative to that of the passive dynamics. It is perhaps not obvious, but nonetheless true, that the modified measure describes the trajectories of a different Markov process. In \Cref{sec:transform.fk}, we formulate the problem of learning the Markov process that correctly describes the trajectories in the Feynman--Kac formula as a Reinforcement Learning problem.

\begin{comment}
\section{Stoquastic Ground States}
The models described in the introduction are all examples of stoquastic Hamiltonians. Stoquastic Hamiltonians, which we will define shortly, are sometimes also described as being `free of the sign problem'. For this reason, they are more amenable to Monte Carlo approaches and their complexity is somewhere between quantum systems with the sign problem and classical systems.

For our purposes, it is important that the ground state of a stoquastic Hamiltonian can be characterized by stochastic equations. One of these stochastic equations is the well-known Feynman--Kac formula, but we will also derive stochastic representations of the Schr\"odinger equation.

\subsection{Stoquastic Hamiltonians}
\label{sec:stoquastic}

Strictly speaking, stoquasticity is a statement about a Hamiltonian matrix that depends on the basis used for the matrix representation. Given a matrix representation of a Hamiltonian in a certain basis, the Hamiltonian matrix is stoquastic if the Hamiltonian matrix has real and non-positive off-diagonal elements. Nevertheless, we will also call the Hamiltonian operator itself stoquastic if there is a known basis for which the Hamiltonian matrix is stoquastic.

Stoquasticity for a Hamiltonian matrix $H$ implies that we can decompose the Hamiltonian matrix as
\begin{equation}
 H = -\Gamma + V
\end{equation}
where $\Gamma$ is a transition rate matrix and $V$ is diagonal. The interpretation is as follows. $\Gamma$ is the kinetic term and defines a continuous time Markov chain. We will call this the `passive dynamics'. The second term is the potential $V$. In terms of the Hamiltonian, the passive transition rates are given by
\begin{equation} \label{eq.passive.dynamics}
 \Gamma_{\vec{s}\to \vec{s}'} =
\begin{cases}
    - H_{\vec{s} \vec{s}'} &\vec{s}\neq \vec{s} \\
    \sum_{\vec{s}' \neq \vec{s}} H_{\vec{s} \vec{s}'} & \vec{s}=\vec{s}' \\
\end{cases}
\end{equation}
and the potential is given by
\begin{equation} \label{eq.potential.from.H}
  V(\vec{s}) = H_{\vec{s}\bs} + \sum_{\vec{s}' \neq \vec{s}} H_{\vec{s} \vec{s}'}.
\end{equation}
So the passive transition rates \eqref{eq.passive.dynamics} describe a continuous time Markov chain. The transitions of this Markov chain are between states $\vec{s}\neq\vec{s}'$ that satisfy $H_{\vec{s} \vec{s}'} \neq 0 $. We will call such states `adjacent'.

Another important property of stoquastic Hamiltonians is the following: if the passive dynamics \eqref{eq.passive.dynamics} is ergodic, then the ground state wavefunction can be chosen positive in the stoquastic basis \citep{crosson_discussion_2017}, i.e.\
\begin{equation} \label{eq:no.node.discrete}
 \varphi_0(\vec{s}) > 0.
\end{equation}
This is a consequence of the Perron-Frobenius theorem, and also explains why stoquastic Hamiltonians are `free of the sign problem': the ground state wavefunction does not change sign.

\subsubsection{Example: the transverse Ising model}
To illustrate, we go back to the transverse field Ising model \eqref{eq:ising}. In the z-spin basis (i.e. a basis state $\vec{s} = (s_1, s_2, \dots, s_N)$ measures the spin $s_i$ of each site $i$ in the z-direction) the Hamiltonian becomes
\begin{equation}
    \braket{\vec{s}'|H|\vec{s}} = -J \sum_{\braket{i,j}} {s_i} {s_j} \braket{\vec{s}'|\vec{s}}
        -h \sum_i \braket{\vec{s}'|\sigma^+_i + \sigma^-_i|\vec{s}}.
\end{equation}
Therefore, the transverse field Ising model is stoquastic: in the z-spin basis, we get $H_{\vec{s}\bs'}\le 0$ if $\vec{s}\neq\vec{s}'$. Furthermore, states $\vec{s}$ and $\vec{s}'$ are adjacent if $s_i=s'_i$ for each site $i$ except one, which is illustrated in Figure \ref{fig:isingandxy}. The passive transition rates are
\begin{equation}
 \Gamma_{{\vec{s}}\to{\vec{s}'}} = h
\end{equation}
if $\vec{s}$ and $\vec{s}'$ are adjacent and zero otherwise. Following \eqref{eq.potential.from.H}, the potential becomes
\begin{align*}
V(\vec{s})
    &= -\sum_{{\vec{s}'}\neq{\vec{s}}} \Gamma_{{\vec{s}}\to{\vec{s}'}} + H_{ss}\\
    &= -hN -  J \sum_{\braket{i,j}} {s_i} {s_j}.
\end{align*}

\begin{figure}[htbp]
\floatconts
  {fig:ising.action
  {\caption{Two adjacent Ising spin states: one spin can be flipped (red icon) to obtain the new state (on the right).}
  {\includegraphics[height=0.43\textwidth]{ising_action.png}\hfill\includegraphics[height=0.43\textwidth]{ising_new_state.png}}
\end{figure}

\subsection{Feynman--Kac representation}
We can use the passive dynamics \eqref{eq.passive.dynamics} of a stoquastic Hamiltonian to give a stochastic representation of its ground state. This representation is known as the Feynman--Kac formula. Although the Feynman--Kac formula is more well-known in continuous state spaces, it also applies to discrete state spaces, as we will describe in this section.

To start with, the Feynman--Kac formula is a stochastic solution of the imaginary time Schr\"odinger equation
\begin{equation} \label{eq.imaginary.schrodinger}
    \partial_t \varphi(t) = - H \varphi(t).
\end{equation}
The imaginary time Schr\"odinger equation is literally the time-dependent Schr\"odinger equation with the substitution $t\to i t$. This is useful because the solution $\varphi(t)$ of the imaginary time Schr\"odinger equation approaches the ground state $\varphi_0$ as $t\to\infty$ \citep{goldberg_integration_1967}.

The Feynman--Kac formula is a stoquastic solution of the imaginary time Schr\"odinger equation with a \emph{terminal} condition at some time $t+T$.
For a stoquastic Hamiltonian with passive rates $\Gamma$ and potential $V$, the Feynman--Kac formula states that the solution at the earlier time $t$ can be expressed as an expectation
\begin{equation}
 \varphi({\vec{s}_t}, {t}) = \E_{\mathbb{P}}
 \left[\exp\left(-\int_{t}^{t+T} V({\vec{s}_{t'}}) dt'\right)\varphi({\vec{s}_{t+T}}, {t+T})\right].
\end{equation}
Here $\varphi({\vec{s}_{t+T}}, {t+T})$ is the \emph{terminal} condition, and the expectation is over trajectories following the passive dynamics $\Gamma$ starting at $s_t$ at time $t$.

Although this may not seem immediately useful because it goes backwards in time, it also implies another characterization of the ground state wavefunction. Since $\varphi({\vec{s}_t}, {t}) = e^{-E_0 t}\varphi_0({\vec{s}_t})$ is an exact solution of the imaginary time Schr\"odinger equation, the Feynman--Kac formula implies
\begin{equation} \label{eq.fk.ground}
\varphi_0(\vec{s}_t) = \E_{\mathbb{P}}
 \left[\exp\left(-\int_{t}^{t+T} \left (V({\vec{s}_{t'}}) - E_0 \right ) dt' \right)\varphi_0({\vec{s}_{t+T}})\right].
\end{equation}

\end{comment}

\subsection{Stochastic representations of the Schr\"odinger equation} \label{sec:stoch}

The Feynman--Kac formula gives a natural stochastic representation of the ground state in terms of a \emph{continuous} time Markov chain. As discussed in \cite{neirotti1996monte}, the ground state of a stoquastic Hamiltonian can also be characterized by the stationary state of a \emph{discrete} time Markov chain. Inspired by this result, we give two stochastic representations of the time-independent Schr\"odinger equation
\begin{equation} \label{eq:schrodinger}
  \sum_{\vec{s}'} H_{\vec{s} \vec{s}'} \varphi_0(\vec{s}') = E_0 \varphi_0(\vec{s}).
\end{equation}
For stoquastic Hamiltonians, the ground state wavefunction $\varphi_0$ is the unique positive eigenfunction.

For the first representation, fix $C > \max_{\vec{s}} H_{\vec{s}\bs}$. For any such constant, we can rewrite the Schr\"odinger equation \eqref{eq:schrodinger} as
\begin{equation}
  \varphi_0(\vec{s}) = \frac{1}{C-E_0} \left [ (C - H_{\vec{s}\bs}) \varphi_0(\vec{s}) - \sum_{\vec{s}'\neq\vec{s}} H_{\vec{s} \vec{s}'} \varphi_0(\vec{s}') \right].
\end{equation}
Since $C > \max H_{\vec{s}\bs}$ and $H_{\vec{s}\bs'}<0$ for $\vec{s}'\neq\vec{s}$, we can turn this into a stochastic representation
\begin{equation} \label{eq:stochastic.schrodinger.1}
 \varphi_0(\vec{s}) = \frac{Z_1(\vec{s})}{C-E_0} \E_{\vec{s}'\sim p_1(\cdot|\vec{s})} \left [ \varphi_0(\vec{s}') \right].
\end{equation}
Here the transition probabilities $p_1(\vec{s}'|\vec{s})$ are
\begin{equation}
    p_1(\vec{s}'|\vec{s}) =
    \begin{cases}
        \left ( C - H_{\vec{s}\bs} \right) / Z_1(\vec{s}) & \vec{s}' = \vec{s} \\
        - H_{\vec{s}\bs'} / Z_1(\vec{s}) & \vec{s}' \neq \vec{s}
    \end{cases}
\end{equation}
with $Z_1(\vec{s}) = C - H_{\vec{s}\bs} - \sum_{\vec{s}'\neq\vec{s}} H_{\vec{s} \vec{s}'}$.

For the second representation, we rewrite the Schr\"odinger equation slightly differently as
\begin{equation}
   \varphi_0(\vec{s}) = -\frac{1}{H_{\vec{s}\bs} - E_0} \sum_{\vec{s}'\neq\vec{s}} H_{\vec{s}\bs'} \varphi_0(\vec{s}')
\end{equation}
This leads to the stochastic representation
\begin{align}\label{eq:stochastic.schrodinger.2}
 \varphi_0(\vec{s}) &= \frac{Z_2(\vec{s})}{H_{\vec{s}\bs} - E_0} \E_{\vec{s}'\sim p_2(\cdot|\vec{s})} \left [ \varphi_0(\vec{s}') \right]\\
 p_2(\vec{s}'|\vec{s}) &=
            \begin{cases}
        0 &  \vec{s}' = \vec{s} \\
        - H_{\vec{s}\bs'} / Z_2(\vec{s}) & \vec{s}' \neq \vec{s}
    \end{cases}
\end{align}
with $Z_2(\vec{s}) = -\sum_{\vec{s}'\neq\vec{s}} H_{\vec{s}\bs'}$. Note that in this second formulation the system always changes state at each step of the process. In both cases, the formulation requires knowledge of the ground state energy $E_0$.

\subsection{Linearly solvable reinforcement learning problems}

The stochastic formulation of the Schr\"odinger equation given in the previous section may be reinterpreted as a linearly solvable Markov decision problem of the type introduced in \cite{todorov2006linearly}.

Todorov considers the maximum entropy reinforcement learning paradigm, where stochastic policies are encouraged by complementing the reward function $r(\vec{s}, \vec{a})$ with the entropy of the policy relative to some reference policy $p$.\footnote{For details see Appendix \ref{sec:maxent.rl}.} As is usual in RL, a problem setting can be summarized and studied by its Bellman equation. Two examples of (soft) Bellman equations are
\begin{subequations}
\begin{equation} \label{eq.u.soft.bellman}
 U^*(\vec{s}) = R^* + \log \E_{\vec{a} \sim p(\cdot|\vec{s})} \left [\exp \left ( r(\vec{s},\vec{a}) + U^*(\vec{a}(\vec{s})) \right )\right ]
\end{equation}
\begin{equation} \label{eq.u.soft.bellman.terminal}
 U^*(\vec{s}) = \log \E_{\vec{a} \sim p(\cdot|\vec{s})} \left [\exp \left ( r(\vec{s},\vec{a}) + U^*(\vec{a}(\vec{s})) \right )\right ].
\end{equation}
\end{subequations}
Here $r(\vec{s}, \vec{a})$ is the reward following the action $\vec{a}$ in the state $\vec{s}$ and $R^*$ is a constant describing the average reward of the optimal policy. $U^*$ is the optimal state-value function and $p$ is the reference policy. Information implicitly contained in these soft Bellman equations is the objective function and time horizon. While \eqref{eq.u.soft.bellman} describes an infinite time horizon, \eqref{eq.u.soft.bellman.terminal} describes a problem with terminal states. Both consider an undiscounted setting where the objective function is simply the sum (or average) of received rewards.

\cite{todorov_efficient_2009} showed that for action-independent rewards, both soft Bellman optimality equations \eqref{eq.u.soft.bellman} and \eqref{eq.u.soft.bellman.terminal} can be transformed into linear equations. Consider action-independent rewards $r(\vec{s},\vec{a}) = r(\vec{s})$. The key is an exponential transformation to the `desirability' $z$ of a state
\begin{equation}
 z(\vec{s}) = \exp(U^*(\vec{s})).
\end{equation}
By simply exponentiating the soft Bellman optimality equations \eqref{eq.u.soft.bellman} and \eqref{eq.u.soft.bellman.terminal}, we obtain
\begin{subequations}
\begin{equation} \label{eq.todorov.desirability}
 z(\vec{s}) = e^{r(\vec{s}) + R^*} \E_{\vec{a}\sim p(\cdot|\vec{s})} \left [ z(\vec{a}(\vec{s})) \right],
\end{equation}
\begin{equation} \label{eq.todorov.desirability.2}
 z(\vec{s}) = e^{r(\vec{s})} \E_{\vec{a}\sim p(\cdot|\vec{s})} \left [ z(\vec{a}(\vec{s})) \right],
\end{equation}
\end{subequations}
which are indeed linear equations for the desirability $z$. { Note that the introduction of discounting results in nonlinear equations, breaking the connection with the Schr\"odinger equation.}

\section{Reinforcement learning for quantum ground states} \label{sec:rl4qm}

Our first contribution is to show how to transform the stochastic representations of \Cref{sec:fk,sec:stoch} to reinforcement learning problems. The key step is a logarithmic transformation that transforms the stochastic representation into a Bellman equation, and represents the inverse of Todorov's approach.

\subsection{Transforming the Feynman-Kac representation}
\label{sec:transform.fk}

The Feynman--Kac representation of the ground state \Cref{eq.fk.ground} may be formulated as a continuous time RL problem. In this section, we give an idea of the proof, see \Cref{sec:fkbellman} for details.

Consider a small time step $\Delta t$ for the Feynman--Kac characterization of the ground state \eqref{eq.fk.ground}
\begin{equation}
\varphi_0(\vec{s}_{t})
    = e^{-V(\vec{s})\Delta t+E_0 \Delta t} \E_{{\mathbb{P}}} \left [  \varphi_0 (\vec{s}_{t+\Delta t})) \right] +{\mathcal{O}}(\Delta t^2).
\end{equation}
Comparison with \eqref{eq.todorov.desirability} gives a reward $r(\vec{s}) = - V(\vec{s}) \Delta t$. Furthermore, the ground state wavefunction takes the role of the desirability $z$, so that we can make the identification
\begin{equation} \label{eq:wave.from.u}
    \log \varphi_0(\vec{s}) = U^*(\vec{s}).
\end{equation}
This identification is visualized in \Cref{fig:bellman} as a Bellman backup diagram. Taking the continuous time limit, this leads to the following maximum entropy RL problem. The control is provided by transition rates $\Gamma^{\theta}$ that drive a continuous time Markov chain for the state $\vec{s}_t$ of the system. We seek the optimal rates $\Gamma^*$
\begin{align}
\Gamma^* &= \mathop{\text{argmax}}_{\Gamma^\theta} R\\
\label{eq:holland.reward.discrete}
  R &= -\lim_{T\to\infty} \frac{1}{T} \E_{s_{0:T}} \left[\int_0^T \left(V(\vec{s}_t) +    \mathcal{H} (\vec{s}_t) \right) dt \right],
\end{align}
where the entropy term $\mathcal{H}$ is given by
\begin{equation} \label{eq:entropy}
    \mathcal{H} (\vec{s}_t) = \sum_{\vec{s}'\neq \vec{s}_t} \left(
    {\Gamma_{\vec{s}_t\to \vec{s}'}}
    - \Gamma^{\theta}_{\vec{s}_t\to \vec{s}'}
    + \Gamma^{\theta}_{\vec{s}_t\to \vec{s}'}
        \log\left(\frac{\Gamma^{\theta}_{\vec{s}_t\to \vec{s}'}}{\Gamma_{\vec{s}_t\to \vec{s}'}}\right) \right),
\end{equation}
and $\Gamma_{s\to s'}$ denotes the passive dynamics of the system
\begin{equation} \label{eq.passive.dynamics}
 \Gamma_{\vec{s}\to \vec{s}'} =
\begin{cases}
    - H_{\vec{s} \vec{s}'} & \vec{s}\neq \vec{s} \\
    \sum_{\vec{s}' \neq \vec{s}} H_{\vec{s} \vec{s}'} & \vec{s}=\vec{s}' \\
\end{cases}
\end{equation}
The expectation in the objective \eqref{eq:holland.reward.discrete} is over trajectories $s_{0:T}$ following the parameterized dynamics $\Gamma^\theta$.

The optimal policy strikes a balance between steering towards states $\vec{s}$ with low potential energy $V(\vec{s})$ and staying close to the passive dynamics $\Gamma$ of \Cref{eq.passive.dynamics}. In physical terms, the objective \eqref{eq:holland.reward.discrete} minimizes the sum of potential energy $V(\vec{s})$ and kinetic energy ${\mathcal{H}}$. { The optimal rates $\Gamma^*$ reproduce the Feynman--Kac measure on the space of trajectories . Note that \eqref{eq:holland.reward.discrete} involves an expectation over the parameterized Markov process. In the case of continuous state spaces, optimizing this expectation is straightforward using the reparameterization trick. The discrete case requires a different approach as described below.}

\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{bellman_backup.png}
    \caption{The Feynman--Kac formula seen as a Bellman backup diagram.}
    \label{fig:bellman}
\end{figure}

\subsection{Transforming the Schrodinger representations}

In a similar way, the stochastic representations of the Schr\"odinger equation from \Cref{sec:stoch} can be mapped to RL in discrete time. The linear equations for the desirability \Cref{eq.todorov.desirability,eq.todorov.desirability.2} will be identified with the stochastic representations \Cref{eq:stochastic.schrodinger.1,eq:stochastic.schrodinger.2} of the Schr\"odinger equation.

The first representation \Cref{eq:stochastic.schrodinger.1} maps to the soft Bellman equation \Cref{eq.u.soft.bellman} i.e.\
\begin{equation}
     U^*(\vec{s}) = R^* + \log \E_{\vec{a} \sim p_1(\cdot|\vec{s})} \left [\exp \left ( r(\vec{s}) + U^*(\vec{a}(\vec{s})) \right )\right ].
\end{equation}
The identification is $\log \varphi_0(\vec{s}) = U^*(\vec{s})$ as before, but now with reward function
\begin{equation}
   r(\vec{s}) = \log Z_1(\vec{s})
\end{equation}
and $R^* = -\log(C - E_0)$.
Therefore, the ground state problem is also equivalent to a discrete time RL problem with an infinite time horizon.

The second representation \eqref{eq:stochastic.schrodinger.2} maps to the soft Bellman equation \eqref{eq.u.soft.bellman.terminal}
\begin{equation}
 U^*(\vec{s}) =  \log \E_{\vec{a}\sim p_2(\cdot|\vec{s})} \left [\exp \left ( r(\vec{s}) + U^*(\vec{a}(\vec{s})) \right )\right ]
\end{equation}
with reward
\begin{equation}
    r(\vec{s}) = \log \left ( \frac{Z_2(\vec{s})}{H_{\vec{s}\bs}-E_0} \right).
\end{equation}
This Bellman equation describes a RL problem with terminal states. In this type of RL problem, there is no fixed time horizon, but exploration stops when one of terminal states is reached. The reward is the total received until the terminal state is reached, which means terminal states necessarily have state-value $U^*=0$.

We are free to choose the terminal states. Choosing a terminal state, i.e.\ setting $U^*(\vec{s})$ to zero, is equivalent to forcing $\varphi_0(\vec{s})=1$ because of the logarithmic transformation $U^*= \log\varphi_0$. Therefore, a choice of terminal state entails a choice of normalization for the wavefunction. A natural candidate for the terminal state is the classical ground state, i.e.\ the state $\vec{s}$ for which $H_{\vec{s}\bs}$ is minimal. Another choice that could be advantageous computationally is a state $\vec{s}$ with minimum symmetry. This is because states $\vec{s}$ and $\vec{s}'$ that are related by a symmetry (such as translation) must obey $\varphi_0(\vec{s}) = \varphi_0(\vec{s}')$, so that choosing a minimum symmetry state gives you many other terminal states for free.

\section{Application: neural quantum states} \label{sec:neural}

In this section, we describe how our mapping of the ground state problem to RL leads to new optimization techniques for neural quantum states. Usually, the quantum wavefunction is parameterized by a neural network $\varphi(\vec{s}) =\mathsf{NN}(\vec{s})$. In contrast, almost all deep reinforcement learning algorithms obtain a neural approximation of either the state-value function $U^*(\vec{s})$ or the action-value function $Q^*(\vec{s}, \vec{a})$. Given a neural approximation of the state-value function $U^*(\vec{s})$, we immediately obtain a neural approximation of the ground state wavefunction
\begin{equation}
 \varphi_0(\vec{s}) = \exp(U^*(\vec{s}))
\end{equation}
because this is the transformation \eqref{eq:wave.from.u} we used to derive the reinforcement learning formulations.

Given a neural approximation of the action-value function $Q^*(\vec{s}, \vec{a})$, we can also obtain a neural quantum state. For this we use the following relation between the state-value function $U^*$ and the action-value function $Q^*$
\begin{equation}
 U^*(\vec{s}) = \log \E_{\vec{a} \sim p(\cdot|\vec{s})} \left [\exp \left (  Q^*(\vec{s}, \vec{a}) \right) \right ].
\end{equation}
Therefore, a neural approximation of $Q^*$ gives a neural approximation of the ground state wavefunction
\begin{equation} \label{eq:wave.from.q}
 \varphi_0(\vec{s}) = \E_{\vec{a} \sim p(\cdot|\vec{s})} \left [\exp \left (  Q^*(\vec{s}, \vec{a}) \right) \right ].
\end{equation}

{

}

\subsection{Advantages of reinforcement learning neural quantum states}
Most papers on Neural Quantum States use either Variational Monte Caro (e.g.\ \cite{han_solving_2019, hermann_deep_2020, kessler_artificial_2019, saito_solving_2017}) or Stochastic Reconfiguration (e.g.\ \cite{choo_two-dimensional_2019, nomura_restricted_2017, pfau_ab-initio_2019}) to optimize neural quantum states. Our optimization method using reinforcement learning could be preferable for the following reasons.

Firstly, our approach needs less data for one update step. For example, Variational Monte Carlo needs to calculate the variational energy for each update step. This means that, given a parametrized wavefunction $\varphi$, samples $\vec{s}\propto \varphi^2$ need to be produced and the so-called local energy
\begin{equation}
    \frac{H\varphi}{\varphi}(\vec{s}) = \frac{1}{\varphi(\vec{s})} \sum_{\vec{s}'\neq\vec{s}} H_{\vec{s}\bs'} \varphi(\vec{s}')
\end{equation}
needs to be averaged over these samples. Therefore, for each update step, new samples need to be generated, and for each sample, the neural wavefunction $\varphi$ needs to be evaluated ${\mathcal{O}}(N)$ times, where $N$ is the number of spins in the system. In contrast, in a Q-learning approach, the action-value function $Q$ can be optimized using off-policy pairs of states $\vec{s}$ and actions $\vec{a}$. The fact that only pairs of states and actions are needed means that the neural action-value function $Q$ only needs to evaluated ${\mathcal{O}}(1)$ times per sample $\vec{s}$. The fact that off-policy data can be used means that generated samples can be reused for different update steps. So one update step can be completed more efficiently using our reinforcement learning approach.

Secondly, our approach can more efficiently produce samples from the neural ground state. To see this, note that such samples are usually generated using the Metropolis--Hastings algorithm. In the Metropolis--Hastings algorithm, the proposal distribution of moves greatly influences its speed of convergence. All reinforcement learning algorithms produce a policy that prefers states of high state-value. In our formulation, that translates to a policy that points towards states of high ground state probability. { In the next section, we show how this leads to a more efficient proposal distribution.
}

{
\subsection{Experiments}
To demonstrate the application to neural quantum states, we consider the two-dimensional transverse Ising model (\Cref{eq:ising}). Training details can be found in \Cref{sec:hparams} and the code at \ifanonsubmission \href{https://github.com/[anonymized]}{https://github.com/[anonymized]}\else \href{https://github.com/WillemGispen/Lattice-QuaRL}{https://github.com/WillemGispen/Lattice-QuaRL} \fi.

We use Soft Q-learning \citep{haarnoja_soft_2018} with discrete space adaptations \citep{christodoulou_soft_2019} to approximate the optimal action-value function $Q^*$ with a convolutional neural network. As described before, this leads via \eqref{eq:wave.from.q} to a neural representation of the ground state wavefunction. One advantage of Christodoulou's discrete space adaptations is that all action-values $Q^*(\vec{s}, \vec{a})$ for a specific state $\vec{s}$ are obtained with one neural network call, so that the neural quantum state \eqref{eq:wave.from.q} may be calculated simply by weighting the sum of this output.

For square lattices up to $6\times 6$, we can compare the variational energy of our optimized wavefunctions with exact diagonalization \citep{hamer_finite-size_2000}. The numerically exact ground state energy of the $6\times 6$ square Ising model, at $h=1$ and $J=0.32758$, is $E_0 = -1.06375$ \citep{hamer_finite-size_2000}. All our methods come within $0.1\%$ of this value, obtaining $-1.0628 (0.09\%)$, $-1.0632 (0.05\%)$, and $-1.0633 (0.04\%)$ for the continuous, non-terminal, and terminal states formulations, respectively. These errors should not be viewed as the limits of each method, but purely as a demonstration that each method can obtain the ground state with reasonable accuracy.

To produce samples from the ground state, we use Metropolis--Hastings with a variety of proposal distributions. The standard proposal distribution is uniform: a single random spin flip is proposed per Metropolis--Hastings step. Using our learned action-values, we can instead propose single spin flips according to $\exp(Q(\vec{s}, \vec{a}))$. Finally, we can sample multiple spin flips from $\exp(Q(\vec{s}, \vec{a}))$ and propose to flip those all at once.

To quantify how efficiently the different proposal distributions sample our optimized $6\times6$ Ising wavefunctions, we measure the autocorrelation of the potential energy. This autocorrelation decays approximately exponentially $\propto \exp(-t/\tau)$. Sampling single spin flips from $\exp(Q(\vec{s}, \vec{a}))$ instead of uniformly approximately halves $\tau$. Sampling $6$ spin flips from $\exp(Q(\vec{s}, \vec{a}))$ results in a $6\times$ speedup of $\tau$. Sampling $10$ spin flips from $\exp(Q(\vec{s}, \vec{a}))$ for a $10\times10$ Ising model (trained with the terminal states method) results in a $10\times$ speedup of $\tau$. These `speedups' refer to comparisons with a single uniformly generated spin flip and are comparable for all our methods. Although our measurements of $\tau$ are only rough indications for the true autocorrelation or mixing times, they suggest the following: sampling $\sqrt{N}$ spin flips from $\exp(Q(\vec{s}, \vec{a}))$ for an Ising system with $N$ spins may reduce autocorrelation/mixing times from ${\mathcal{O}}(N)$ to ${\mathcal{O}}(\sqrt{N})$. Anyhow, it seems clear that using $\exp(Q(\vec{s}, \vec{a}))$ to propose spin flips can signicantly reduce autocorrelation/mixing times. This could be particularly beneficial when autocorrelation/mixing times are long, such as for large systems or when random actions have very low acceptance rates.

}

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

In this paper we have introduced three different reinforcement learning formulations of the problem of finding the ground state of a quantum lattice model, as summarized in \Cref{tab:overview}.

Even before doing {extensive} experiments, we can note some relative advantages and disadvantages of the different RL formulations. Firstly, discrete time RL settings are more commonly studied than continuous time ones, so it is easier to find efficient RL algorithms to solve the discrete time formulations. A disadvantage of the third formulation is that the reward contains the ground state energy $E_0$. Still, an estimate for the ground state energy may be obtained from the variational energy, which must be calculated anyway for validation purposes. Moreover, sparse terminal states are challenging, although they can be dealt with \citep{agostinelli2019solving}.

\begin{table}[]
\begin{tabular}{l|lll}
Basis                  & Time       & Time horizon    & Reward                      \\
\hline
Feynman-Kac formula    & Continuous & Infinite        & $V(s)$                      \\
Schr\"odinger equation & Discrete   & Infinite        & $\log Z_1 (\vec{s}) $           \\
Schr\"odinger equation & Discrete   & Terminal states & $\log \left ( \frac{Z_2(\vec{s})}{H_{\vec{s}\bs} - E_0} \right)$
\end{tabular}
\caption{Overview of the three different reinforcement learning formulations of the ground state problem.}
\label{tab:overview}
\end{table}

Future experiments will explore the performance of our three approaches relative to each other and existing neural representations of quantum states. As well as being the basis of methods in their own right, the learned process described by the rates $\Gamma_\theta$ could be used for importance sampling of the Feynman--Kac formula, which will be optimal (i.e. with zero variance) for $\Gamma_\theta=\Gamma^*$.

\acks{We thank Ariel Barr for early discussions of this work and collaboration on a related project. AL acknowledges support from EPSRC Grant No. EP/P034616/1.}

redpajama_set_name: RedPajamaArXiv


text:
Renée Graham
For Trump, impeachment means civil war
By Renée Graham ,Updated October 2, 2019, 5:00 a.m.
Illustration by Lesley Becker/Globe Staff | Adobe | Globe file photo
If the Civil War had truly ended more than 150 years ago, there would be no President Trump.
The same seething white resentment that doomed the post-war Reconstruction era ushered Trump into the White House. Now he's promoting a hate pastor's rubbish that America will face a "Civil War like fracture in the Nation from which our Country will never heal," if he is removed from office.
Keep this in mind as this increasingly unhinged and desperate president faces a deepening impeachment inquiry. We already know — from Charlottesville to Pittsburgh, from Christchurch to El Paso — Trump's incendiary rhetoric has garnered a body count.
Because he's such a grotesque caricature of a man and a president, it's easy to laugh at Trump's sheer ridiculousness. After Trump tweet-quoted Christian extremist Robert Jeffress's dire, hyperbolic prediction made (where else?) on "Fox & Friends Weekend," #CivilWarSignup began burning up Twitter with endless jokes about imaginary battles, fighting prowess (or lack thereof), and even appropriate attire for a civil war.
None of this is funny.
Not when white hate groups are talking about a "full blown 'hot' civil war" and calling the left "domestic enemies of the Constitution bent on the destruction of the Republic." Never mind that the Constitution has no greater enemy, in this or any nation, than the man who took an oath to defend it when he was sworn in as the 45th president.
And not when Trump, a merchant of lies and discord, hints at potential violence if he is impeached. He is rallying his supporters with the belief that an attack on him is also an attack on them.
Back in 2017, just four months into his presidency when the "i-word" was already bubbling up after he fired FBI director James Comey, Trump sent out a fund-raising e-mail claiming "the media is out to get us." It continued: President Trump is waging "an uphill battle — and we need to be prepared to go into the trenches to FIGHT BACK."
Of course, this has always been Trump's message to supporters — if he goes down, they go down. Politicians love the all-for-one, one-for-all inclusiveness of plural pronouns, but Trump has weaponized them. It's not the first time he has implied that his supporters would "demand" to keep him in office, even beyond a second term — if he lasts that long.
Like the coward he is, Trump hides behind tough-guy poses and talk. He often rhapsodizes about a time when "troublemakers" were dealt with "a little differently than we do now." That's what he said about the whistle-blower — a true American patriot — who exposed information about Trump's improper phone conversation with the president of Ukraine. His implied violence was very clear.
Naturally, it's never Trump who gets his hands dirty. When violence erupts — and hate crimes have soared since he became president — he feigns innocence even as his own invective indicts him. His hate speech is prominent in manifestos by several mass shooters, including the alleged murderer of 22 people at an El Paso Walmart in August.
Trump's words are like IEDs detonating in the violent minds of fragile white men. And now, finally facing consequences for his presidency's pervasive corruption, he's fomenting notions of war in America. It's so outrageous even Republican Representative Adam Kinzinger of Illinois called Trump's comment "beyond repugnant."
If only Trump's comments were merely repugnant, and not potentially lethal.
With 47 percent of Americans now favoring his impeachment and removal from office, Trump is again mining white resentment, the oxygen of his campaign and presidency, as his lifeline. Convinced it will keep him in office, he will drag this nation through mayhem, terror, and grief beyond imagination.
Renée Graham can be reached at renee.graham@globe.com. Follow her on Twitter @reneeygraham.

redpajama_set_name: RedPajamaCommonCrawl


text:
Study: Closing tax breaks for retirement would do little to help budget
By Lori Montgomery
Lori Montgomery
Business Editor
Email Bio Follow
As Washington considers closing a profusion of tax breaks in its quest to reduce record budget deficits, a study suggests that some targets might produce far less cash than is commonly thought.
Tax preferences for retirement plans, such as 401(k) accounts, rank among the biggest money losers for the federal government, with official projections showing a loss of about $600 billion to federal coffers over the next five years. But a study set for release this week suggests that the accounts would have to be closed, or contributions sharply limited, to realize much in the way of budget savings — a move that could undermine some of the nation's most popular retirement tools.
"Given the questions about the future of Social Security, this is no time to be making bad policy decisions about private savings," said Brian Graff, executive director of the American Society of Pension Professionals and Actuaries, which conducted the study.
The analysis highlights the murky math surrounding hundreds of tax breaks that collectively deprive the U.S. Treasury of more than $1 trillion each year. Known as "tax expenditures," the list includes dozens of popular credits, deductions and other policies that benefit many special interests, but also millions of ordinary taxpayers.
The most expensive tax breaks from the government's standpoint are the tax-free treatment of employer-provided health benefits and the mortgage-interest deduction for homeowners. Those provisions allow taxpayers to keep nearly $300 billion a year that would otherwise be collected in taxes, according to administration estimates.
Although that figure gauges the value of the policies to taxpayers, it says nothing about how much of that cash could be recovered through repeal or reform, tax experts say. For example, if the mortgage interest deduction were eliminated, people would probably shift their investments to other tax-preferred vehicles, thereby denying the tax man a portion of his expected rewards, said Ed Kleinbard, a former director of the congressional Joint Committee on Taxation, which, along with the Treasury, estimates the value of federal tax expenditures.
The estimates are "an accurate measure of the extent to which the preference is being used by taxpayers," said Kleinbard, who teaches law at the University of Southern California. "But that's a completely different question than how much revenue could be raised if the preference were eliminated."
Tax preferences for retirement savings are particularly tricky, Kleinbard said, because they are not typically straightforward tax breaks, whose entire value is realized each year. In the case of 401(k)s, for example, taxpayers defer taxation on their contributions until they retire. In the meantime, they accumulate interest on their contributions tax free.
To determine the loss to the government, the Joint Committee on Taxation and the Treasury look at taxes lost each year because of new 401(k) contributions, as well as the interest earned on existing contributions, and then subtract taxes paid on distributions from 401(k)s. The difference ranks 401(k) retirement plans as the third most expensive provision on the administration's list of tax expenditures, projected to deprive the Treasury of an estimated $67 billion in 2012.
Because 401(k)s have only been around for a few decades, however, the value of the contributions and interest earnings still far exceeds the value of distributions to retirees, Graff said. That means traditional methods dramatically overstate the amount of revenue that ultimately will be lost.
The study says that a more accurate measure of retirement savings incentives would estimate the value of the tax breaks over the lifetime of the accounts. By that measure, the study shows, the loss to the government would be as much as 77 percent lower than official estimates.
"If the system were mature, then it might not make much difference. But the system is still growing," said Judy Miller, the society's director of regulatory affairs. Because "the taxes being deferred are far more significant than the amounts being paid out," she said, the official estimates are "very distorted."
Graff noted that President Obama's fiscal commission recommended slashing the total amount people could contribute to retirement plans as part of its proposal to rein in the national debt. Lawmakers, he said, have also told the group that everything is on the table.
Kleinbard said he has seen no move in Congress to limit retirement incentives such as the 401(k), the nation's biggest retirement savings program with 70 million participants. "I do think they are being a bit paranoid," he said of the society.
Still, Roberton Williams, a senior fellow at the nonpartisan Tax Policy Center, said there are good arguments for reexamining tax breaks for retirement savings, which disproportionately benefit the wealthy.
"You could design these things in ways that would encourage people to save for retirement, but in a much more equitable way that would actually impact behavior," Williams said. "The rich are saving anyway. It's the poor guy that needs the incentive."
Lori Montgomery Lori Montgomery oversees coverage of business, the economy and technology. Before becoming an editor, she was a reporter covering U.S. economic policy, with a focus on efforts to tame the national debt. Follow


text:
Q: mvc partial views loses track of the images folder when using jquery This is what I have and it works:
$(function(){
$('.slide-out-div').tabSlideOut({
        tabHandle: '.handle',                     //class of the element that will become your tab
        pathToTabImage: 'http://mhmiisdev2/images/contact_tab.gif', //path to the image for the tab //Optionally can be set using css
        imageHeight: '122px',                     //height of tab image           //Optionally can be set using css
        imageWidth: '40px',                       //width of tab image            //Optionally can be set using css
        tabLocation: 'right',                      //side of screen where tab lives, top, right, bottom, or left
        speed: 300,                               //speed of animation
        action: 'click',                          //options: 'click' or 'hover', action to trigger animation
        topPos: '200px',                          //position from the top/ use if tabLocation is left or right
        leftPos: '20px',                          //position from left/ use if tabLocation is bottom or top
        fixedPosition: true                      //options: true makes it stick(fixed position) on scroll
    });
});

This is what I want and it doesnt work when I change from one controller to another controller. NOTICE THE IMAGE PATH IS NOT ABSOLUTE
$(function(){
    $('.slide-out-div').tabSlideOut({
        tabHandle: '.handle',                     //class of the element that will become your tab
        pathToTabImage: '/images/contact_tab.gif', //path to the image for the tab //Optionally can be set using css
        imageHeight: '122px',                     //height of tab image           //Optionally can be set using css
        imageWidth: '40px',                       //width of tab image            //Optionally can be set using css
        tabLocation: 'right',                      //side of screen where tab lives, top, right, bottom, or left
        speed: 300,                               //speed of animation
        action: 'click',                          //options: 'click' or 'hover', action to trigger animation
        topPos: '200px',                          //position from the top/ use if tabLocation is left or right
        leftPos: '20px',                          //position from left/ use if tabLocation is bottom or top
        fixedPosition: true                      //options: true makes it stick(fixed position) on scroll
    });
});

the html for completeness....

    Content

Medical Variance Reports

Individual Medicines

I suspect somehow pathToTabImage: /images/contact_tab.gif' loses its context when browsing throught controllers. Help me understand...

A:
I suspect somehow pathToTabImage: /images/contact_tab.gif' loses it
  context when browsing throught controllers.

Oh yes, you are correct. Always use url helpers when dealing with urls in ASP.NET MVC:
pathToTabImage: '@Url.Content("~/images/contact_tab.gif")'

Never hardcode urls as you did in your code. This will ensure that your application will work when you deploy it in IIS under a virtual directory.

A: Just use the the helper Url.Content() function inside your JS like this
pathToTabImage: @Url.Content("~/images/contact_tab.gif")
If your JS is inside an external Javascript file you get just refer to an variable that is global.
At the top of the view:

var myPath = @Url.Content("~/images/");

Inside external JS file:
pathToTabImage: myPath + "contact_tab.gif"

A: Here is what I did
 alert(UrlContent("/hdd/images.gif")); //for debugging purposes
 alert(UrlContent("~/hdd/images.gif"));
 alert(UrlContent("hdd/images.jpg");
 //they all return the correct url

function UrlContent(url) {

// first lets take care of users that are smart

url = $.trim(url.replace("~/", "")); //this should take care of smart users!
if (url.charAt(0) == "/") {
    url=url.substring(1,url.length)
}
//now that we have taken care of smart users lets form the url

var UrlPath = location.protocol + '//' + location.hostname +location.port+'/'+url

return UrlPath

}

redpajama_set_name: RedPajamaStackExchange


text:
Het seizoen 1962/63 van de Belgische Derde Klasse ging van start op 1 september 1962 en eindigde op 16 juni 1963. KSV Waregem won in Derde Klasse A, K. Boom FC won in Derde Klasse B.

Gedegradeerde teams
Deze teams waren gedegradeerd uit de Tweede Klasse voor de start van het seizoen:
 K. Sint-Niklase SK (rechtstreeks)
 RRC Tournai (rechtstreeks)

Gepromoveerde teams
Deze teams waren gepromoveerd uit de Vierde Klasse voor de start van het seizoen:
 VC Zwevegem Sport (kampioen 4A)
 KFC Vigor Hamme (kampioen 4B)
 R. Stade Waremmien FC (kampioen 4C)
 FC Eendracht Houthalen (kampioen 4D)

Promoverende teams
Deze teams promoveerden naar Tweede Klasse op het eind van het seizoen:
 KSV Waregem (kampioen 3A)
 K. Boom FC (kampioen 3B)

Degraderende teams
Deze teams degradeerden naar Vierde Klasse op het eind van het seizoen:
 RAA Louviéroise (rechtstreeks uit 3A)
 RCS Brainois (rechtstreeks uit 3A)
 R. Jeunesse Arlonaise (rechtstreeks uit 3B)
 K. Tubantia Borgerhout FC (rechtstreeks uit 3B)

Eindstand

Derde Klasse A

Derde Klasse B

Titelwedstrijd Derde Klasse
De twee teams die kampioen werden in hun reeks, KSV Waregem en K. Boom FC, speelden een heen -en terugwedstrijd om te bepalen wie de winnaar van Derde Klasse zou worden.

Referenties
 Belgian Soccer Database
 RSSSF

Derde klasse voetbal België
Belgie 3

redpajama_set_name: RedPajamaWikipedia


12 So his father agreed to divide his wealth between his sons.
37 As Jesus was speaking, one of the Pharisees invited him home for a meal. So he went in and took his place at the table.
1 One Sabbath day Jesus went to eat dinner in the home of a leader of the Pharisees, and the people were watching him closely.
20 And while he was still a long way off, his father saw him coming. Filled with love and compassion, he ran to his son, embraced him, and kissed him.
Through this parable Jesus challenges what nearly everyone has ever thought about God, sin, and salvation. Jesus's purpose is not warm our hearts but to shatter our categories. Jesus is redefining everything we thought we knew about connecting to God. He redefining sin, what it means to be lost, and what it means to be saved. pp.10, 28.
9 The Lord isn't really being slow about his promise, as some people think. No, he is being patient for your sake. He does not want anyone to be destroyed, but wants everyone to repent.
29 And people will come from all over the world—from east and west, north and south—to take their places in the Kingdom of God.
42 Many people did believe in him [Jesus], however, including some of the Jewish leaders. But they wouldn't admit it for fear that the Pharisees would expel them from the synagogue.
38 Afterward Joseph of Arimathea, who had been a secret disciple of Jesus (because he feared the Jewish leaders), asked Pilate for permission to take down Jesus' body. When Pilate gave permission, Joseph came and took the body away.
7 So God's message continued to spread. The number of believers greatly increased in Jerusalem, and many of the Jewish priests were converted, too.
8 Crispus, the leader of the synagogue, and everyone in his household believed in the Lord.
I have been a member of the Pharisees, the strictest sect of our religion. 6 Now I am on trial because of my hope in the fulfillment of God's promise made to our ancestors.
28 The older brother was angry and wouldn't go in. His father came out and begged him.