Efficient Cosmological Parameter Estimation with Hamiltonian Monte Carlo

1
Efficient Cosmological Parameter Estimation with
Hamiltonian Monte Carlo

• Amir Hajian
• Cosmo06 September 25, 2006

Astro-ph/0608679
2
Parameter estimation
NASA/WMAP science team
Fig. M. White 1997
3
The Problem

• Power Spectrum calculation takes a long time for
  large l
• Likelihood takes time too
• Lengthy chains are needed specially for
• Curved distributions
• Non-Gaussian distributions
• High dimensional parameter spaces

4
Possible Solutions

• Speed up the calculations
• Parallel computation
• Power Spectrum
• CMBWarp, Jimenez et al (2004)
• Pico, Fendt Wandelt (2006)
• CosmoNet , Auld et al (2006)
• Likelihood
• Improve MCMC method
• Reparametrization, e.g. Verde et al (2003)
• Optimized step-size, e.g. Dunkley et al (2004)
• Parallel chains
• Use more efficient MCMC algorithms,
• e.g. CosmoMC, Cornish et al (2005), HMC.

5
Traditional (Random Walk) Metropolis Algorithm
Current position
p(x)
6
Traditional (Random Walk) Metropolis Algorithm
Proposed position
p(x)
p(x)
p(x) gt p(x) accept the step
7
Traditional (Random Walk) Metropolis Algorithm
Proposed position
p(x)
p(x)
p(x) lt p(x) accept the step with probability
p(x)/p(x) Otherwise take another sample at x
8
Traditional (Random Walk) Metropolis Algorithm
9
Issues with MCMC

• Long burn-in time
• Correlated samples
• Low efficiency in high dimensions
• Low acceptance rate

10
Hamiltonian Monte Carlo

• Proposed by
• Duan et al, Phys. Lett. B, 1987
• Used by
• condensed matter physicists,
• particle physicists and
• statisticians.
• Uses Hamiltonian dynamics to perform big
  uncorrelated jumps in the parameter space.

11
Hamiltonian Monte Carlo
p(x)
x
Define the potential energy U(x) -Log(p(x))
12
Hamiltonian Monte Carlo
U(x)
x
13
Hamiltonian Monte Carlo
U(x)
Give it an initial momentum
u(x)
x
Total energy H(x)U(x)1/2u2
14
Hamiltonian Monte Carlo
U(x)
Evolve the system for a given time Hamiltonian
dynamics
u(x)
u(x)
x
H(x) U(x) K(x)
15
Hamiltonian dynamics
H conserved, only if done accurately
16
Hamiltonian dynamics (in practice)

• Discretized time-steps
• Leapfrog method

Total energy may not remain conserved
Accept the proposed position according to the
Metropolis rule
?/2
?/2
u(t?) x(t?)
u(t) x(t)
u(t?/2)
?
17
Extended Target Density
Sample from H(x,u) Marginal distribution of x is
p(x)
18
How does it work?

• Assume Gaussian distribution
• Trajectories in the phase space
• Randomizing the momentum in the beginning of each
  leapfrog guarantees the coverage of the whole
  space

Fig. K. Hanson, 2001
19
Hamiltonian Monte Carlo
20
Important questions

• Are we sampling from the distribution of
  interest?
• Are we seeing the whole parameter space?
• How many samples do we need to estimate the
  parameters of interest to a desired precision?
• How efficient is our algorithm?

21
Convergence Diagnostics
Autocorrelation
22
Convergence Diagnostics

• Xi
• P(k) ?k2

FFT
?k
23
Convergence Diagnostics

• Power spectrum P(k)
• Averaged

Ideal sampler
Flat
24
Efficiency of MCMC sequence

• ratio of the number of independent draws from the
  target pdf to the number of MCMC iterations
  required to achieve the same variance in an
  estimated quantity.
• For a Gaussian distribution
• Where P0P(k0)

See Dunkley et al (2004) for more details
25
Example Gaussian PDFSampled with different
chains
Better efficiency
Low efficiency
26
Example

• Simplest example Gaussian distribution
• Energy

27
Comparison Acceptance Rate
HMC 100
MCMC 25
28
Comparison Correlations
29
Comparison distributions
30
Comparison Efficiency
Compare to 1/D behavior of the efficiency of
traditional MCMC methods
31
(No Transcript)
32
Cosmological Applications
33
Flat 6-parameter LCDM model

• 0th approximation
• Approximate
• the Lnlikelihood by
• Estimate the fit parameters from an exploratory
  MCMC run.
• Evaluate Gradients,
• Run HMC.

34
Result

• Acceptance rate boosted up to 81 while reducing
  the correlation in the chain.
• Good improvement, but can do better!

35
Better approximation for gradients

• Modified likelihood routine of Pico (Fendt and
  Wandelt, 2006) to evaluate the gradient.

36
Lico (Likelihood routine of Pico)
F(x)
x
Cut the parameter space into pieces and fit a
different function to each piece.
37
The Gradient
38
Flat 6-parameter LCDM model
Acceptance rate 98
39
Correlation lengths
40
Summary

• HMC is a simple algorithm that can improve the
  efficiency of the MCMC chains dramatically
• HMC can be easily added to popular parameter
  estimation softwares such as CosmoMC and
  AnalyzeThis!
• HMC can be used along with methods of speeding up
  power spectrum and likelihood calculations,
• HMC is ideal for curved, non-Gaussian and
  hard-to-converge distributions,
• Approximations made in evaluating the gradient
  just reduce the acceptance rate, but dont get
  propagated into the results of parameter
  estimation.
• It is easy to get a non-optimized HMC, but hard
  to get a wrong answer!