Log-posterior calculation involve PDE numerical analysis


#1

Hello, I have a quick question. I need to fit 25 parameters within a set of PDEs. (Partial Differential Equation).
The way to fit is abstracted as follow:

Random sampling the 25 parameters
Plug into the PDEs then do numerical analysis(iteration)
Calculate a set of result which can be considered as probabilities of binomial distribution
Calculate the likelihood along with the collected data

I’m asking if it’s possible to use STAN and the HMC algo to sample the 25 parameters till the likelihood reaches the maximum (MLE) ?
Since the PDE solving(also log-posterior) includes iterative numerical methods, I don’t know if the STAN auto-differentiation algorithm will work or not.

Because the log-posterior includes a bunch of codes, I come to ask if it’s doable before I start to translate the original MATLAB code into the STAN’s functions.


#2

I’m asking if it’s possible to use STAN and the HMC algo to sample the 25 parameters till the likelihood reaches the maximum (MLE)

So HMC is gonna generate samples from your posterior. If it’s an MLE you want, then you can just run an optimizer on that.

Since the PDE solving(also log-posterior) includes iterative numerical methods, I don’t know if the STAN auto-differentiation algorithm will work or not.

It’s possible it can, but the autodiff for something like this is probably a little slow. If you want to generate a couple thousand samples from your posterior and NUTS needs to take 32 HMC steps for each one (which is pretty efficient sampling), then you’re looking at 64000 model evaluations/gradient calculations. With a PDE or something you probably want to think about custom autodiff.

That said there are numerous pretty complicated functions in the Stan math library that are autodiffed. The eigenvalue solver is one.

Random sampling the 25 parameters

Hmm, Stan doesn’t really directly draw samples from it’s posterior distribution. It’s less of a guess-and-check process and more of an exploring process. Check out the Michael Betancourt or Radford Neal Intro HMC papers. There are some cool animations here: https://chi-feng.github.io/mcmc-demo/app.html#HamiltonianMC,standard

Plug into the PDEs then do numerical analysis(iteration)

Is this like an iterative sparse solver?

Calculate a set of result which can be considered as probabilities of binomial distribution

The results of the PDE are a distribution? Or get fed into a distribution as data? Can you write out what a hypothetical Stan model would look like? Do you have one that has maybe just a few parameters – 3-4 maybe instead of 25?


#3

Also amusingly enough I learned FEM from a guy named Yang Xu (UCSB).


#4

For inverse problem like this, for each HMC step the computing effort is in the inversion of the Hessian of the Tikhonov functional. If your parameter-to-observable mapping is linear and using a Gussian prior as well as likelihood, this is equivalent to the posterior cov matrix. I think stan can handle this. On the other hand, I’d be concerned about the amount of computing. Since you opt an iterative solver, I assume the linear system from PDE discretization is large. This could make it impractical to use MCMC. Another possible issue is the convergence tolerance of your Krylov solver, as if the solution is not accurate enough the number of effective samples can decrease.

So maybe you can profile it by hard coding a low dimensional matrix inversion in stan and examining its performance.

This is a problem where parallel could really help.