[Notes]0908_MCMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feng 2018-09-08\n",
    "\n",
    "This revisit of the topic is initially motivated by the poundering of inverse Ising problems, with potential applications to a generic clustering algorithm that is purely based on pairwise distance matrix. The initial optimisation attempts indicate that gradient-based methods are capable of inferring a small coupling matrix from the un-normalised likelihood function, provided that temperature is low. The partition function is harder to evaulate, but would be dominate by the lowest energy state, so should be well approximated by its energy and redundancy. \n",
    "\n",
    "But given the lattice/discrete nature of the Ising model, markov-chain monte-carlo may actually represents a better approach. This paper [Lee and Daniels 2018](https://www.researchgate.net/publication/322713743_Convenient_Interface_to_Inverse_Ising_ConIII_A_Python_package_for_solving_maximum_entropy_models) summarised existing algorithms for inverse Ising, and alluded to me the minimal probability flow MPF which is introduced in this [PRL letter] (https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.107.220601). MPF seeks to optimise a model-dependent MCMC dynamics wrt model parameters. Or more specifically, minimises the KL divergence between the data distribution and the model's asymptotic/stationary distribution (not sure which direction it is but model is definitely more general than data so model is zero-forcing). The key here is to replace the asymptotic distribution with the instanteous probability flow at the data distribution. The main issue would be choosing the connectivity matrix/proposal distribution. \n",
    "\n",
    "Initially I thought the MPF might help with the gradient-based dynamics in variational methods, and remembered there was a thing called Langevin dynamics. Upon reading this is simply injecting random Gaussian noises into the stochastic gradients so that it does not collapse into a point estimate [Wellling and Teh 2011](https://www.ics.uci.edu/~welling/publications/papers/stoclangevin_v6.pdf). Stochastic gradient Langevin Dynamics (SGLD) is attractive in it allows seamless transition from optimisation to sampling posterior. Nevertheless, Langevin dynamics is in fact just naive random walk [video](https://www.youtube.com/watch?v=qBf5EBdEw7Q) which makes it rather inefficient in sampling.\n",
    "\n",
    "Recall that pyMC3 concerns MCMC, I looked up its tutorial and apparently it used no-U-turn-sampler (NUTS) to speed up MC in continuous space [pyMC3](https://docs.pymc.io/notebooks/getting_started.html). In fact, [NUTS](https://chi-feng.github.io/mcmc-demo/app.html#NaiveNUTS,banana) is a n intuitive extension of Hamiltonian Monte Carlo [HMC](http://arogozhnikov.github.io/2016/12/19/markov_chain_monte_carlo.html). Other popular MC algorithms are visualised [here](https://chi-feng.github.io/mcmc-demo/) (beware of memory crash!)\n",
    "\n",
    "Upon rethinking of how HMC is implemented using conservation of (kinetic energy + potential energy), it is actually precisely the missing piece I was looking for realising the \"surfer\", an algorithm that traverse the equipotential surface rather than slipping down on it. In fact, HMC is a much more principled way than naively fiddling with the discretised updating equations, and all I needs to do is to replace the log-prob potential with a function that measures distance to the equipotential surface $|logP - E|$, and then simply leapfrog the resultant Hamiltonian dynamics with suitable initialisation. This is exactly the opposite of NUTS in a way: we need to finish a whole cycle to traverse and hence measure the equipotential surface. Potentially we will need to deal with the local curvature with a rescaling matrix, or by reparametrisation, but using insight from AdaDelta should be enough."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The study of equipotential surface is obviously interesting: \n",
    "\n",
    "1. The MAP estimator is a zero-dimensional equipotential surface.\n",
    "\n",
    "1. The shape of equipotential surface near MAP could inform the volume of the posterior distribution near it."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}