How to solve Path Integral Problems with Probabilistic Programming Techniques?

Specific Domains Probabilistic programming
1

Hello everyone, I’m not very familiar with probabilistic programming, so I may not be approaching this problem in the most optimal way. Essentially, I’m trying to solve a problem of the following form:

Z \sim \int\int...\int e^{-U(R_1, R_2 ... R_N)} d R_1 \ldots d R_N

Here, Z is typically considered a density, while the array R_i have the same distribution. The function U(R_1, R_2 ... R_N) can be decomposed into \sum {V(R_i) + T(R_i, R_{i+1}) }.

I’ve tried using DynamicHMC.jl and AdvancedHMC.jl, but the samples produced are almost identical; the differences between the samples (i.e., R_i) are smaller than 1E-3. Furthermore, generating 10,000 samples with approximately 800 dimensions takes about 2 hours.

In many situations, < U(R_1, R_2 ... R_N) > is more important than generating samples. Is it possible to obtain < U(R_1, R_2 ... R_N) > directly with less computational cost?

I would greatly appreciate any advice or suggestions on how to tackle this problem more effectively. Thank you in advance!

2

If the limiting factor is the generation of sample points, and not the evaluation of U, you could try a method that generates samples more quickly but is perhaps less statistically efficient. e.g. if you are integrating over a hypercube, the Sobol package can generate 10000 low-discrepancy samples in 800 dimensions in less than a second on my laptop — just evaluate U at each point and take the average. (And you can also easily parallelize it.) The QuasiMonteCarlo package automates this for you.

3

At first glance, this looks like the kind of problem where you need pruning/enriching to get good statistics. See https://www.jstage.jst.go.jp/article/tjsai/16/2/16_2_279/_pdf for an overview. If you go this route, you actually cannot calculate Z accurately, only the averages of the surviving lineages \langle U \rangle. The U(R_i) will then act as your potential / fitness function at each “time” step from which particles sample according to the conditional distribution governed by T. Not familiar with the Julia versions of this, but not that hard to roll your own.

The structure of your problem is a bit weird though, should it be T(R_{i+1} - R_{i}) like in classical-action discretization? Can we hear a little bit more about these functions?

Edit: the boundary conditions and what you can sample are also quite important and what you have closed distributions for. Sequentially generating independent R_i samples according to the conditional distribution will of course generate completely iid samples from the target distribution.

4

I have heard that before, Quasi Monte Carlo is very interesting. I’m a bit worried that the weights of the generated sample points will be relatively small. After all I’ll still have a try.
Thank you for your advice!

5

This is related to importance-sampled (Quasi-)Monte-Carlo. e.g. the Cuba package has an implementation of this.

6

Thank you! I would I will read it carefully.

Oh, sure! Sorry, I made a mistake, it should be T(R_i, R_i+1). There is a little change for sampling bosons and fermions, which discribed in https://pnas.org/doi/full/10.1073/pnas.1913365116 and http://aip.scitation.org/doi/10.1063/5.0008720. And I am working on using Slater determinant wavefunction to solve it.

Could you please give more information about it? It is similar to Difussion Monte Carlo.

7

I am not sure, but could Cuba give out some samples? Thank you for your reply.

8

Could you please give more information about it? It is similar to Difussion Monte Carlo.

Yes, as you’ve written as a 1-D sum, the target distribution looks like a Markov chain:

p(R_{i+1}|R_{i}) d R_{i+1} \propto e^{U(R_{i+1}) + T(R_{i+1} + R_{i})} d R_{i+1}

so

p(R_{n}|R_{n-1}) ··· p(R_{2}|R_{1}) p(R_{1}|R_{0}) π(R_0) d \vec{R} = p(R_0,R_1,···,R_n)d \vec{R}

So you can sample from the target for calculating averages if you can sample from p(R_{i+1}|R_{i}) .

9

It sounds exciting. The method just like the Difussion Monte Carlo, but different. It make problem easy to solve. I think I can make the problem much simpler.
Thank you! I would always learn a lot.

By the way, it always add p(R_n|R_1) to make it like a circle but it’s not necessary.

Sorry, I just make a mistake, it should be T(R_i, R_{i-1})

1 Like