Correct application of Laplacian operator

Question

Not a physicist, and I'm having trouble understanding how to apply the Laplacian-like operator described in this paper and the original. We let:

$$ \hat{f}(x) = f(x) + \frac{\int H(x,y)\psi(y) dy}{\sqrt{\pi(x)}} $$

Where $H$ is a Hermitian operator (real symmetric, actually), $\pi(x)$ is an un-normalized density function (measurable positive), $\psi(x)$ is arbitrary integrable function, and f(x) is an arbitrary measurable function. $H$ is chosen to ensure that the following condition holds:

$$ \int H(x,y)\sqrt{\pi(x)} dx = 0 \tag{1} $$

This is done so that $\hat{f}(x)$ and $f(x)$ have the same expectation under $\pi(x)$: $\int f(x)\pi(x) dx = \int \hat{f}(x)\pi(x) dx $.

The paper goes on to define a Schrödinger-like $H$ for $\{x \in \mathbb{R}^d\}$:

$$H = -\frac{1}{2}\sum_{i=1}^d \frac{\partial^2}{\partial x_i^2} + V(x)$$

Where $V(x)$ is constructed to meet (1):

$$ V(X) = \frac{1}{2\sqrt{\pi(x)}} \sum_{i=1}^d \frac{\partial^2 \sqrt{\pi(x)}}{\partial x_i^2} $$

Which allows us to verify $(1)$ as follows:

$$ \begin{align} H \sqrt{\pi} =& \int H(x,y)\sqrt{\pi(y)} dy \\ =& \int -\frac{1}{2}\sum_{i=1}^d \frac{\partial^2 \sqrt{\pi(y)}}{\partial y_i^2} + \frac{\sqrt{\pi(y)}}{2\sqrt{\pi(y)}} \sum_{i=1}^d \frac{\partial^2 \sqrt{\pi(y)}}{\partial y_i^2} dy = 0\\ \end{align} $$

I got the right result above, but I'm not sure I applied the H operator correctly...

The paper considers the following form for $\psi(x)$:

$$ \psi(x) = P(x)\sqrt{\pi(x)}$$

The authors then derive the following:

$$ \hat{f}(x) = f(x) -\frac{1}{2} \Delta P(x) + \nabla P(x) \cdot \left(-\frac{1}{2}\nabla \ln \pi(x)\right) $$

Where $\nabla$ denotes the gradient $(\frac{\partial}{\partial x_1},...,\frac{\partial}{\partial x_d})$ and $\Delta$ denotes the Laplacian operator $\sum_{i=1}^d \frac{\partial^2}{\partial x_i^2}$. I'm not able to re-derive this equation. Again, I'm not sure I'm applying the H operator correctly, because it looks like I should end up with something in terms of $x$. My obviously incorrect attempt follows:

$$ \begin{align} \int H(x,y) \psi(y) dy =& \int -\frac{1}{2}\sum_{i=1}^d \frac{\partial^2 \psi(y)}{\partial y_i^2} + \frac{\psi(y)}{2\sqrt{\pi(y)}} \sum_{i=1}^d \frac{\partial^2 \sqrt{\pi(y)}}{\partial y_i^2} dy\\ =& \int -\frac{1}{2}\sum_{i=1}^d \frac{\partial^2 P(y)\sqrt{\pi(y)}}{\partial y_i^2} + \frac{P(y)}{2} \sum_{i=1}^d \frac{\partial^2 \sqrt{\pi(y)}}{\partial y_i^2} dy\\ =& \int -\frac{\sqrt{\pi(y)}}{2} \Delta P(y) - \frac{1}{2\sqrt{\pi(y)}}\nabla P(y) \cdot \nabla \pi(y) dy \end{align} $$

Edit:

I was able to reproduce the result by applying the operator by just swapping $x$ and $y$ and not explicitly solving the integral. here is what I got:

$$ \begin{align} H\psi =& -\frac{1}{2}\sum_{i=1}^d \frac{\partial^2 \psi(x)}{\partial x_i^2} + \frac{\psi(x)}{2\sqrt{\pi(x)}} \sum_{i=1}^d \frac{\partial^2 \sqrt{\pi(x)}}{\partial y_i^2}\\ =& -\frac{\sqrt{\pi(x)}}{2} \Delta P(x) - \frac{1}{2\sqrt{\pi(x)}} \nabla P(x) \cdot \nabla \pi(y) \\ \end{align} $$ Dividing by $\sqrt{\pi(x)}$ yields: $$ = -\frac{1}{2} \Delta P(x) - \frac{1}{2} \nabla P(x) \cdot \nabla \ln \pi(x) $$ Which is the same. Does this mean that the definition of $H$ is describing the solution to the integral? What is $H(x,y)$?

Answer 1

In the original paper it is said that equation (10): $$ H = -\frac{1}{2}\sum_{i=1}^d \frac{\partial^2}{\partial x_i^2} + V(x) \qquad (10) $$ "is written using the standard quantum-mechanical notation for a local Hamiltonian in the x-space realization".

This means that actually $H$ from (10) is the operator and $H(x,y)$ from your eq. (1) is the matrix of this operator in the linear space of integrable functions with the following basis: $$ f_y(\xi) = \delta(\xi-y). $$

So one can write: $$ H\psi(x) = \int H(x,y)\psi(y) \; dy. $$

The elements of the matrix can be calculated as follows: $$ H(x,y) = H_{x,y} = (f_x, Hf_y) = \int f_x^*(\xi) H(\xi) f_y(\xi) \; d\xi = $$ $$ \int f_x^*(\xi) \left(-\frac{1}{2}\sum_{i=1}^d \frac{\partial^2}{\partial \xi_i^2} + V(\xi)\right) f_y(\xi) \; d\xi = $$ $$ \int \delta(\xi-x) \left(-\frac{1}{2}\sum_{i=1}^d \delta''_i(\xi-y) + V(\xi) \delta(\xi-y) \right) \; d\xi. $$ So $$ H(x,y) = -\frac{1}{2}\sum_{i=1}^d \delta''_i(x-y) + V(x) \delta(x-y), $$ where $$ \delta''_i(r) = \frac{\partial^2}{\partial r_i^2} \delta(r). $$

Answer 2

While the previous answer is complete for what you were asking (you were only missing the delta-function factor in the matrix element of a local operator), I would like to motivate the arbitrary looking factors of $\sqrt{\pi}$ used in the method, and explain it in a different, perhaps more illuminating way, because I think it is very interesting.

When you have a quantum Hamiltonian of the Schrodinger form $H=\nabla^2 + V$, the ground state is a positive definite real function $\psi_0(x)= e^{-{W(x)\over 2}}$, where the previous equation defines W, and by adding a constant to the energy, you can make the ground state energy exactly zero. Under these conditions, time evolution keeps the ground state the ground state, with no phases in time.

The ground state itself can be interpreted as a probability measure for the imaginary time continuation of the quantum mechanical Hamiltonian. The imaginary time Schrodinger equation is a diffusion process, and the Markov Chain version is given by the following transformation on the Hamiltonian (I believe originally described this way, possibly with factor of 2 or squaring errors, I didn't check this recently, in the Wikipedia page for "quantum superposition"):

$$ T(x,y) = {1\over \psi_0(x)} H(x,y) \psi_0(y) $$

While H is symmetric, T is not. Where H has the property that it kills the ground state, T has the property that it is stochastic for continuous time evolution:

$$ \sum_y T(x,y) = 0 $$

The transformation between the two is a pointwise rescaling of the linear space of amplitudes by the diagonal matrix whose x,x entry is $\psi_0(x)$.

The stochastic nature of the T matrix means that the equation

$$ \dot{\rho} = \sum_x \rho(x) T(x,y) $$

has a probability interpretation--- it describes a Markov Chain in imaginary time. Notice now that the stationary distribution of this Markov Chain is:

$$ \pi(x) = \psi_0^2(x) = e^{-W(x)}$$

So that the factors of $\sqrt{\pi}$ are easy to understand--- they are transforming back and forth between the quantum mechanics and the imaginary time stochastic mechanics using the ground state wavefunction. The quantity W is called the SUSY superpotential (this is not true--- everyone other than me calls its derivative the superpotential, but I think this is a stupid convention and it might not be too late to change it), and the expression for the potential V from W is well known from SUSY QM (the Schrodinger map between stochastic Markov-chain evolution and quantum mechanics is known to have a supersymmetry relating the forward and backward version of the Markov chain, as described well in a book by Junkers)

$$ V(x) = {1\over 2} |\nabla W|^2 + {1\over 2} \nabla^2 W $$

These relations are used in this work without explaining where they appear most naturally. You can become completely familiar with the transformation between unitary QM and stochastic imaginary time Markov processes (or Langevin equations when space is continuous), so that the factors become obvious, and the formulas are no longer opaque.

The basic idea

The method in the paper is as follows: given an operator f(x), you want to find its expected value in the stationary distribution of a Markov chain,

$$\int f(x)\mu(x) dx $$

You consider the function $f(x)$ as a function of the operator X, and map the Markov chain to a quantum mechanics problem, so that you are after

$$ \langle 0 | f(x) |0\rangle $$

The issue with evaluating this by Monte-Carlo is that f(x) produces states other than the ground state when acting on the ground state, in an energy basis, $|E_i\rangle$, there are nonzero matrix elements

$$ f_{i0} = \langle E_i | f(x) |0\rangle $$

You want to zero out all the $f_{i0}$, or make them as small as possible. You then note that the Hamiltonian H kills $|0\rangle$ (on either side), so it acts as a projection operator to the complement of $|0\rangle$, with additional rescaling. So if you find an operator A with the property that

$$ \hat{f}_{ij} = \langle E_i | f(x) - H A |0 \rangle = 0 $$

for all i not equal to 0, you have made the matrix element of $\hat{f}$ in the vacuum equal to the mean of $f$, so that the monte carlo sampling will always sample the mean. This is an exact reformulation of the problem of computing the mean, so there is nothing gained. But if you make any ansatz at all for A, like taking A to be $\psi(x)$ with some specific variational formula for $\psi$ (this $\psi$ has no relation to the wavefunction, it's an unfortunate notation in the paper), you can reduce the matrix elements of f between the vacuum and other states (the variance of f in the vacuum, the variance of f in the measure) by varying the constants in $\psi$. This is the main idea of the method.

I think that formulating it in this way makes the conceptual relation to other work more apparent, although the formulations in the original paper you cite is self contained and clear, of course.