SMAC Talk

# .orange[HaRAM]: .orange[**Ha**]miltonian .orange[**R**]epelling .orange[**A**]ttracting .orange[**M**]etropolis

### .bolder[Siddharth Vishwanath & Hyungsuk Tak]

---
class: center, middle
count:false

.footnote.tiny.pull-left[

 

 
Courtesy the [template by Chi Feng](https://github.com/chi-feng/mcmc-demo)
]

---
class: inverse, center, middle, animated slideInUp slideOutRight
count: false
class: left
# .left[Sampling]

<div style="display:none">
$$
\newcommand{\defeq}{\stackrel{\small\bullet}{=}}
\newcommand{\ra}{\rangle}
\newcommand{\la}{\langle}
\newcommand{\norm}[1]{\left\|#1\right\|}
\newcommand{\abs}[1]{\left\lvert#1\right\rvert}
\newcommand{\Abs}[1]{\Bigl\lvert#1\Bigr\rvert}
\newcommand{\pr}{{\mathbb P}}
\newcommand{\qr}{{\mathbb Q}}
\newcommand{\muv}{{\boldsymbol{mu}}}
\newcommand{\xv}{{\boldsymbol{x}}}
\newcommand{\av}{{\boldsymbol{a}}}
\newcommand{\bv}{{\boldsymbol{b}}}
\newcommand{\cv}{{\boldsymbol{c}}}
\newcommand{\dv}{{\boldsymbol{d}}}
\newcommand{\ev}{{\boldsymbol{e}}}
\newcommand{\fv}{{\boldsymbol{f}}}
\newcommand{\gv}{{\boldsymbol{g}}}
\newcommand{\hv}{{\boldsymbol{h}}}
\newcommand{\nv}{{\boldsymbol{n}}}
\newcommand{\sv}{{\boldsymbol{s}}}
\newcommand{\tv}{{\boldsymbol{t}}}
\newcommand{\uv}{{\boldsymbol{u}}}
\newcommand{\vv}{{\boldsymbol{v}}}
\newcommand{\wv}{{\boldsymbol{w}}}
\newcommand{\zerov}{{\mathbf{0}}}
\newcommand{\onev}{{\boldsymbol{1}}}
\newcommand{\phiv}{{\boldsymbol{\phi}}}
\newcommand{\cc}{{\check{C}}}
\newcommand{\xv}{{\boldsymbol{x}}}
\newcommand{\Xv}{{\boldsymbol{X}\!}}
\newcommand{\yv}{{\boldsymbol{y}}}
\newcommand{\Yv}{{\boldsymbol{Y}}}
\newcommand{\zv}{{\boldsymbol{z}}}
\newcommand{\Zv}{{\boldsymbol{Z}}}
\newcommand{\Iv}{{\boldsymbol{I}}}
\newcommand{\Jv}{{\boldsymbol{J}}}
\newcommand{\Cv}{{\boldsymbol{C}}}
\newcommand{\Ev}{{\boldsymbol{E}}}
\newcommand{\Fv}{{\boldsymbol{F}}}
\newcommand{\Gv}{{\boldsymbol{G}}}
\newcommand{\Hv}{{\boldsymbol{H}}}
\newcommand{\alphav}{{\boldsymbol{\alpha}}}
\newcommand{\epsilonv}{{\boldsymbol{\epsilon}}}
\newcommand{\betav}{{\boldsymbol{\beta}}}
\newcommand{\deltav}{{\boldsymbol{\delta}}}
\newcommand{\gammav}{{\boldsymbol{\gamma}}}
\newcommand{\etav}{{\boldsymbol{\eta}}}
\newcommand{\piv}{{\boldsymbol{\pi}}}
\newcommand{\thetav}{{\boldsymbol{\theta}}}
\newcommand{\tauv}{{\boldsymbol{\tau}}}
\newcommand{\muv}{{\boldsymbol{\mu}}}
\newcommand{\phiinv}{\Phi^{-1}}
\newcommand{\Fiinv}{F^{-1}}
\newcommand{\giinv}{g^{-1}}
\newcommand{\fhat}{\hat{f}}
\newcommand{\ghat}{\hat{g}}
\newcommand{\ftheta}{f_\theta}
\newcommand{\fthetav}{f_{\thetav}}
\newcommand{\gtheta}{g_\theta}
\newcommand{\gthetav}{g_{\thetav}}
\newcommand{\ztheta}{Z_\theta}
\newcommand{\xtheta}{\Xv_\theta}
\newcommand{\ytheta}{\Yv_\theta}
\newcommand{\p}{\partial}
\newcommand{\f}{\frac}
\newcommand{\cf}{\cfrac}
\newcommand{\e}{\epsilon}
\newcommand{\indep}{\perp\kern-5pt \perp}
\newcommand{\inner}[1]{\langle#1\rangle}
\newcommand{\pa}[1]{\left(#1\right)}
\newcommand{\pb}[1]{\left\{#1\right\}}
\newcommand{\pc}[1]{\left[#1\right]}
\newcommand{\pA}[1]{\Big(#1\Big)}
\newcommand{\pB}[1]{\Big\{#1\Big\}}
\newcommand{\pC}[1]{\Big[#1\Big]}
\newcommand{\ty}[1]{\texttt{#1}}
\newcommand{\borel}[1]{\mathscr{B}\pa{#1}}
\newcommand{\scr}{\mathcal}
\newcommand{\scrb}{\mathscr}
\newcommand{\argmin}{\mathop{\text{arg}\ \!\text{min}}}
\newcommand{\arginf}{\mathop{\text{arg}\ \!\text{inf}}}
\newcommand{\argmax}{\mathop{\text{arg}\ \!\text{max}}}
\newcommand{\argsup}{\mathop{\text{arg}\ \!\text{sup}}}
\newcommand{\bigo}[1]{\mathcal{O}_{p}\!\left(#1\right)}
\newcommand{\f}{\frac}
\newcommand{\e}{\epsilon}
\newcommand{\inv}{^{-1}}
\newcommand{\phiinv}{\Phi^{-1}}
\newcommand{\Fiinv}{F^{-1}}
\newcommand{\giinv}{g^{-1}}
\newcommand{\fhat}{\hat{f}}
\newcommand{\ghat}{\hat{g}}
\newcommand{\ftheta}{f_\theta}
\newcommand{\fthetav}{f_{\thetav}}
\newcommand{\gtheta}{g_\theta}
\newcommand{\gthetav}{g_{\thetav}}
\newcommand{\ztheta}{Z_\theta}
\newcommand{\xtheta}{\Xv_\theta}
\newcommand{\ytheta}{\Yv_\theta}
\newcommand{\absdet}[1]{\abs{\det\pa{#1}}}
\newcommand{\jac}[1]{\Jv_{#1}}
\newcommand{\absdetjx}[1]{\abs{\det\pa{\Jv_{#1}}}}
\newcommand{\absdetj}[1]{\norm{\Jv_{#1}}}
\newcommand{\sint}{sin(\theta)}
\newcommand{\cost}{cos(\theta)}
\newcommand{\sor}[1]{S\mathcal{O}(#1)}
\newcommand{\ort}[1]{\mathcal{O}(#1)}
\newcommand{\A}{{\mathcal A}}
\newcommand{\C}{{\mathbb C}}
\newcommand{\E}{{\mathbb E}}
\newcommand{\F}{{\mathcal{F}}}
\newcommand{\N}{{\mathbb N}}
\newcommand{\R}{{\mathbb R}}
\newcommand{\Q}{{\mathbb Q}}
\newcommand{\Z}{{\mathbb Z}}
\newcommand{\X}{{\scr{X}}}
\newcommand{\Y}{{\mathbb{Y}}}
% \newcommand{\G}{{\mathcal{G}}}
\newcommand{\M}{{\mathcal{M}}}
\newcommand{\betaequivalent}{\beta\text{-equivalent}}
\newcommand{\betaequivalence}{\beta\text{-equivalence}}
\newcommand{\Mb}{{\boldsymbol{\mathsf{M}}}}
\newcommand{\Br}{{\mathbf{\mathsf{Bar}}}}
\newcommand{\dgm}{{\mathbf{\mathsf{Dgm}}}}
\newcommand{\Db}{{\mathbf{\mathsf{D}}}}
\newcommand{\Img}{{\mathbf{\mathsf{Img}}}}
\newcommand{\mmd}{{\mathbf{\mathsf{MMD}}}}
\newcommand{\Xn}{{\mathbb{X}_n}}
\newcommand{\Yn}{{\mathbb{Y}_n}}
\newcommand{\Xb}{{\mathbb{X}}}
\newcommand{\Yb}{{\mathbb{Y}}}
\newcommand{\s}{{{\sigma}}}
\newcommand{\fnsbar}{{\bar{f}^n_\s}}
\newcommand{\fns}{{f^n_\s}}
\newcommand{\fs}{{f_\s}}
\newcommand{\fsbar}{{\bar{f}_\s}}
\newcommand{\barfn}{{{f}^n_\sigma}}
\newcommand{\barfnm}{{{f}^{n+m}_\sigma}}
\newcommand{\barfo}{{{f}_\sigma}}
\newcommand{\fn}{{f^n_{\rho,\sigma}}}
\newcommand{\fnm}{{f^{n+m}_{\rho,\sigma}}}
\newcommand{\fo}{{f_{\rho,\sigma}}}
\newcommand{\K}{{{K_{\sigma}}}}
\newcommand{\barpn}{{\bar{p}^n_\sigma}}
\newcommand{\barpo}{{\bar{p}_\sigma}}
\newcommand{\pn}{{p^n_\sigma}}
\newcommand{\po}{{p_\sigma}}
\newcommand{\J}{{\mathcal{J}}}
\newcommand{\B}{{\mathcal{B}}}
\newcommand{\pt}{{\tilde{\mathbb{P}}}}
\newcommand{\Winf}{{W_{\infty}}}
\newcommand{\HH}{{{\scr{H}_{\sigma}}}}
\newcommand{\D}{{{\scr{D}_{\sigma}}}}
\newcommand{\Ts}{{T_{\sigma}}}
\newcommand{\Phis}{{\Phi_{\sigma}}}
\newcommand{\nus}{{\nu_{\sigma}}}
\newcommand{\Qs}{{\mathcal{Q}_{\sigma}}}
\newcommand{\ws}{{w_{\sigma}}}
\newcommand{\vs}{{v_{\sigma}}}
\newcommand{\ds}{{\delta_{\sigma}}}
\newcommand{\fp}{{f_{\pr}}}
\newcommand{\prs}{{\widetilde{\pr}_{\sigma}}}
\newcommand{\qrs}{{\widetilde{\qr}_{\sigma}}}
\newcommand{\Inner}[1]{\Bigl\langle#1\Bigr\rangle}
\newcommand{\innerh}[1]{\langle#1\rangle_{\HH}}
\newcommand{\Innerh}[1]{\Bigl\langle#1\Bigr\rangle_{\HH}}
\newcommand{\normh}[1]{\norm{#1}_{\HH}}
\newcommand{\norminf}[1]{\norm{#1}_{\infty}}
\newcommand{\gdelta}{{\G_{\delta}}}
\newcommand{\supgdelta}{{\sup\limits_{g\in\gdelta}\abs{\Delta_n(g)}}}
\newcommand{\id}{\text{id}}
\newcommand{\supp}{\text{supp}}
\newcommand{\cech}{\v{C}ech}
\newcommand{\Zz}{{\scr{Z}}}
\newcommand{\psis}{\psi_\s}
\newcommand{\phigox}{\Phis(\xv)-g}
\newcommand{\phigoy}{\Phis(\yv)-g}
\newcommand{\fox}{{f^{\epsilon,{\xv}}_{\rho,\sigma}}}
\newcommand{\prx}{{\pr^{\epsilon}_{\xv}}}
\newcommand{\pro}{{\pr_0}}
\newcommand{\dotfo}{\dot{f}_{\!\!\rho,\s}}
\newcommand{\phifo}{{\Phis(\yv)-\fo}}
\newcommand{\phifox}{{\Phis(\xv)-\fo}}
\newcommand{\kinf}{{\norm{\K}_{\infty}}}
\newcommand{\half}{{{\f{1}{2}}}}
\newcommand{\Jx}{\J_{\epsilon,{\xv}}}
\newcommand{\dpy}{\text{differential privacy}}
\newcommand{\edpy}{$\epsilon$--\text{differential privacy}}
\newcommand{\eedpy}{$\epsilon$--edge \text{differential privacy}}
\newcommand{\dpe}{\text{differentially private}}
\newcommand{\edpe}{$\epsilon$--\text{differentially private}}
\newcommand{\eedpe}{$\epsilon$--edge \text{differentially private}}
\newcommand{\er}{Erdős-Rényi}
\newcommand{\krein}{Kreĭn}
% \newcommand{\grdpg}{\mathsf{gRDPG}}
% \newcommand{\rdpg}{\mathsf{RDPG}}
% \newcommand{\eflip}{{\textsf{edgeFlip}}}
\newcommand{\grdpg}{\text{gRDPG}}
\newcommand{\rdpg}{\text{RDPG}}
\newcommand{\eflip}{{\text{edgeFlip}}}
\newcommand{\I}{{\mathbb I}}
\renewcommand{\pa}[1]{\left(#1\right)}
\renewcommand{\pb}[1]{\left\{#1\right\}}
\renewcommand{\pc}[1]{\left[#1\right]}
%
%
%
\newcommand{\N}{\mathcal{N}}
\newcommand{\D}{\boldsymbol{\nabla}}
\newcommand{\d}[1]{{\partial}_{#1}}
\newcommand{\dd}[1]{\f{d}{d #1}}
% \newcommand{\ddd}[2]{\f{\partial #2}{\partial #1}}
\newcommand{\ddd}[2]{\f{d #2}{d #1}}
\newcommand{\q}{\boldsymbol{q}}
\newcommand{\q}{\boldsymbol{q}}
\renewcommand{\tr}{^{\top}}
\newcommand{\e}{\epsilon}
\newcommand{\qo}{{\boldsymbol{q}_0}}
\newcommand{\p}{{\boldsymbol{p}}}
\newcommand{\tp}{{\widetilde{\boldsymbol{p}}}}
\newcommand{\z}{{\boldsymbol{z}}}
\newcommand{\tz}{\widetilde{\boldsymbol{z}}}
\newcommand{\init}{\textup{init}}
\newcommand{\half}{{\f12}}
\newcommand{\vhalf}{{1/2}}
\newcommand{\F}{{\boldsymbol{\Psi}}}
\newcommand{\P}{{\boldsymbol{\Phi}}}
\newcommand{\flip}{\mathbf{F}}
\newcommand{\vec}[1]{\begin{bmatrix}#1\end{bmatrix}}
\newcommand{\S}{\Sigma}
\newcommand{\Ov}{\mathbf{O}}
\newcommand{\Om}{\boldsymbol{\Omega}}
\newcommand{\G}{\boldsymbol{\Gamma}}
% \newcommand{\inv}[1]{{}^{#1}}
$$
</div>

Generate samples `$\pb{ \q_1, \q_2, \dots, \q_n } \sim \pi$` from a target distribution 
.content-box-blue[
$$
  \pi(\q) \propto e^{- U(\q) }
$$
]

### .orange[Methods]
1. .large[Rejection sampling ❌ ] 
2. .large[Random-walk Metropolis ]✅ 
--

3. .large[Metropolis adjusted Langevin algorithm ]😮 
4. .large[Hamiltonian Monte Carlo ]😲

--
5. .large[Stein Variational gradient descent ]🤯 
6. .large[Wasserstein gradient flows ]🤯 
7. .large[Normalizing flows ]🤯 
8. ...

---
class: left
count: false

# .left[Objective]

.Large.content-box-orange[<body> Can we design a .bolder.dblue.text-shadow[sampler] which: 
1. Can efficiently sample from .bolder.dblue.text-shadow[multi-modal distributions]? 
2. .bolder.dblue.text-shadow[Preserves all the nice properties] of existing mainstays?
</body>]

---
class: inverse, center, middle
count: false

# Background
## (⚡️ Speed)

<!-- ---
class: left

### .dblue[Robust TDA]

.large[&nbsp;&nbsp;&nbsp;&nbsp; 2. Efficient and Robust Topological Inference]
   
### .dblue[Statistical Inference using TDA]

### .dblue[Differential Privacy under TDA Lens]
.large[&nbsp;&nbsp;&nbsp;&nbsp; 4. The Shape of Edge Differential Privacy]

### .dblue[Geometry in Multimodal Sampling and Non-convex Optimization]

---

.pull-left[
.bold.text-shadow.red[Random Walk Metropolis]
<body> 
1) Given the current state $\q_n$, propose a new state
$$ 
\q_{n+1} \sim \kappa(\q|\q_n)
$$

2) Accept / Reject with probability
$$
\alpha(\q_{n+1} | \q_n) = 1 \wedge \f{\kappa(\q_{n}|\q_{n+1})\pi(\q_{n+1})}{\kappa(\q_{n+1}|\q_{n})\pi(\q_{n})}
$$
</body>]
--

.content-box-orange.center[ Wastes too much time in high dimensions. ]

---
class: left
# More intelligent diffusions

.pull-left[
.bold.text-shadow.red[Langevin Monte Carlo]
<body> 
1) Given the current state $\q_n$, propose a new state
$$ 
\q_{n+1} \sim \q_n + {h}\D\log\pi(\q_n) + 2h \N(\zerov, \Sigma)
$$

2) Accept / Reject with probability
$$
\alpha(\q_{n+1} | \q_0) = 1 \wedge \f{\kappa(\q_{n}|\q_{n+1})\pi(\q_{n+1})}{\kappa(\q_{n+1}|\q_{n})\pi(\q_{n})}
$$
</body>]
--

.content-box-orange.center[ Only works for small step sizes! ]

---

# Hamiltonian Monte Carlo

---
class: center, middle

---
class: left
# Why walk when you can flow?

1. Augment state space with auxiliary variables `$\p \sim \N(\boldsymbol{0}, \S)$`
   
2. Joint distribution `$(\q,\p) \sim \exp(-H(\q,\p))$` where `$H(\q,\p)$` is given by
`\begin{align}
H(\q,\p) = U(\q) + \f{1}{2} \p^{\top}\S^{-1}\p
\end{align}`

3. Treat `$H(\q,\p)$` as the Hamiltonian of a system and generate trajectories using Hamiltonian dynamics

<body>$$\dd{t}\begin{bmatrix}\q_t\\ \p_t\end{bmatrix} = \begin{bmatrix}\Ov & \I\\ -\I & \Ov\end{bmatrix}\begin{bmatrix}\D_{\!\q}H(\q_t, \p_t)\\ \D_{\!\p}H(\q_t, \p_t)\end{bmatrix} = \begin{bmatrix}\S\inv\p_t\\ - \D U(\q_t)\end{bmatrix}.$$
</body>

1. Accept/reject state `$(\q_t, \p_t)$` with probability 
.center.content-box-blue[<body>$$\alpha(\q_t, \p_t | \q_0, \p_0) = \min\pb{1, \f{\kappa(\q_0, \p_0 | \q_t, \p_t) e^{-H(\q_t,\p_t)}}{\kappa(\q_t, \p_t | \q_0, \p_0) e^{-H(\q_0,\p_0)}} \Bigg| \ddd{(\q_0, \p_0)}{(\q_t, \p_t)} \Bigg| }.$$
</body>]

---
class: center
count: false
# .left[Why walk when you can flow?]

---
class: center
count: false
# .left[Why walk when you can flow?]

---
class: center
count: false
# .left[Why walk when you can flow?]

---
class: center
count: false
# .left[Why walk when you can flow?]

---

# The Four Pillars of Hamiltonian Monte Carlo

.theorem[<body> 
For every trajectory $\pb{(\q_t, \p_t) : t > 0}$ satisfying the Hamiltonian dynamics
$$
\dd{t}H(\q_t, \p_t) = 0.
$$
</body>]
]
.pull-right[
<img src=images/gifs/wo-friction2.gif height="270" />
]

.content-box-blue.center[<body>
$H(\q_t, \p_t) = H(\q_0, \p_0)$
</body>]

]

.theorem[<body>Let: 
1. $\flip: (\q, \p) \mapsto (\q, -\p)$ denote a .bolder.text-shadow[momentum flip] 
2. $\F_t: (\q_0, \p_0) \mapsto (\q_t, \p_t)$ be the .bolder.text-shadow[flow operator]

Then
$$
\flip \circ \F_t \circ \flip \circ \F_t = \id.
$$
</body>]
]
.pull-right[
<img src=images/gifs/reversible.gif height="270"/>
]

.content-box-blue.center[<body>
Detailed balance condition is satisfied, i.e., $\kappa(\q_t, \p_t | \q_0, \p_0) = \mathbf{1}.$
</body>]

]

.theorem[<body>
For every trajectory $\pb{(\q_t, \p_t) : t > 0}$$$\ddd{(\q_0, \p_0)}{(\q_t, \p_t)} = \Bigg| \begin{bmatrix} \ddd{\q_0}{\q_t} & \ddd{\q_t}{\p_0} \\ \ddd{\p_t}{\q_0} & \ddd{\p_t}{\p_0} \end{bmatrix} \Bigg| = 1.$$
</body>]
]
.pull-right[
![](https://statmodeling.stat.columbia.edu/wp-content/uploads/2017/01/catpendulum.png)
]

.content-box-blue.center[<body>$$\alpha(\q_t, \p_t | \q_0, \p_0) = \min\pb{1, \f{\kappa(\q_0, \p_0 | \q_t, \p_t) e^{-H(\q_t,\p_t)}}{\kappa(\q_t, \p_t | \q_0, \p_0) e^{-H(\q_0,\p_0)}} \Bigg| \ddd{(\q_0, \p_0)}{(\q_t, \p_t)} \Bigg| } = 1.$$
</body>]

]
.panel[.panel-name[Symplecticity]

.theorem[<body>Let $\omega(\cdot, \cdot)$ be the .bolder.text-shadow[symplectic 2-form]
$$
 \omega(\q,\p) = \sum_{i=1}^d dq_i \wedge dp_i.
$$
Then $\omega(\q_t, \p_t) = \omega(\q_0, \p_0)$.
</body>]
]
.pull-right[
![](https://www.ias.edu/sites/default/files/styles/wysiwyg_half/public/media-assets/Nelson_symplectic_form.jpg)
]

.content-box-blue.center[<body>
Enables efficient numeric approximations. 
</body>]

]

---
class: left

# Symplectic Integration

* The solution `$\F_t: (\q_0, \p_0) \mapsto (\q_t, \p_t)$` to the system of .bolder.text-shadow.dblue[differential equations]
$$
\ddd{t}{\q_t} = \S\inv\p_t, \quad\quad \ddd{t}{\p_t} = - \D U(\q_t) 
$$
is rarely available in practice. It can be numerically solved using the .bolder.red.text-shadow[leapfrog scheme]. 
--

* Let `$\e \approx dt$` be a small .bolder.red[step-size], then `$\F_{dt} \approx \F_\e: (\q_t, \p_t) \mapsto (\q_{t+\e}, \p_{t+\e})$` is given by

.content-box-red.center[<body>$$\begin{cases}
 \p_{t+\f \e 2} = \p_t - \f \e 2 \D U(\q_t) \\[5pt]
 \q_{t+\e} = \q_t + \e \S\inv \p_{t+\f \e 2}\\[5pt]
 \p_{t+\e} = \p_{t+\f \e 2} - \f \e 2 \D U(\q_{t+\e}) \end{cases}
 $$
</body>]
--

.content-box-orange.center[<body>For time $T$ and $L = \lfloor{T/\e}\rfloor$, we have $\F_T \approx \F_{\e, L} = (\F_\e)^{\otimes L}$.
</body>]

---
class: center

# .left[Why symplectic integration?]

.content-box-orange[<body> `$|H\left( \F_{T}(\q,\p) \right) - H\left( \F_{\e, L}(\q,\p) \right)| \approx O(\e^3)$`. &nbsp;&nbsp;&nbsp;&nbsp; Therefore, `$\alpha(\q_T, \p_T | \q_0, \p_0) \approx e^{O(\e^3)}.$`</body>]

---
class: left

# .left[Is HMC the panacea?]

.theorem-blue[<body>For $b \in \R$ and $\muv = b \mathbf{1}_d$ consider the multimodal target
$$
\pi(\q) \sim \half \mathcal{N}(\q | -\muv, \I_d) + \half \mathcal{N}(\q | +\muv, \I_d),
$$
and let $\mathcal{E}(b, d)$ denote the event that for initial state $\q_0=-\muv$
$$
\mathcal{E}(b, d) = \pb{ \norm{\P_T(\q_0) - \muv} \le \norm{\P_T(\q_0) + \muv}}
$$
Then for all $b > \sqrt{2/d}$
$$
\pr\pa{\mathcal{E}(b, d)} \le \exp\bigg(-\f{1}{8} \pa{\f{b^2d - 2}{b}}^2\bigg)
$$
</body>]

.center[<body>Markov chain mixes well only when `$U(\q)$` is strongly convex. Terrible for multimodal distributions.</body>]

---
class: center
count: false

# .left[Is HMC the panacea?]

Markov chain mixes well only when `$U(\q)$` is strongly convex. Terrible for multimodal distributions.

---
class: center
count: false

# .left[Is HMC the panacea?]

Markov chain mixes well only when `$U(\q)$` is strongly convex. Terrible for multimodal distributions.

---

# HaRAM
## Hamiltonian Repelling-Attracting Metropolis

---
class: left

# Observation #1: .orange.text-shadow[Friction dissipates energy]

Consider the trajectory of a particle `$(\q_t, \p_t)$` on a rough surface with friciton `$\gamma > 0$`

$$
\ddd{t}{\q_t} = \S\inv\p_t, \quad\quad \ddd{t}{\p_t} = - \D U(\q_t) - \gamma\p_t
$$

This can be rewritten as

<body>$$\dd{t}\begin{bmatrix}\q_t\\ \p_t\end{bmatrix} = \underbrace{\begin{bmatrix}\Ov & \I\\ -\I & \Ov\end{bmatrix}}_{=\Om}\underbrace{\begin{bmatrix}\D_{\!\q}H(\q_t, \p_t)\\ \D_{\!\p}H(\q_t, \p_t)\end{bmatrix}}_{=\D H(\q_t, \p_t)} - \underbrace{\begin{bmatrix}\Ov & \Ov\\ \Ov & \gamma\I\end{bmatrix}}_{=\G}\begin{bmatrix}\q_t\\ \p_t\end{bmatrix}.$$
</body>

Therefore,

.content-box-blue.center[<body>Eq. (A) is called a .text-shadow.red[conformal Hamiltonian system]. </body>]

---
class: left
count: false

# Observation #1: .orange.text-shadow[Friction dissipates energy]

---
class: left
count: false

# Observation #1: .orange.text-shadow[Friction dissipates energy]

---
class: left

# Observation #2: .orange.text-shadow[Negative friction accumulates energy]

If we just flip the sign of the friction parameter `$\gamma$` we get

<body>$$\dd{t}\begin{bmatrix}\q_t\\ \p_t\end{bmatrix} = \underbrace{\begin{bmatrix}\Ov & \I\\ -\I & \Ov\end{bmatrix}}_{=\Om}\underbrace{\begin{bmatrix}\D_{\!\q}H(\q_t, \p_t)\\ \D_{\!\p}H(\q_t, \p_t)\end{bmatrix}}_{=\D H(\q_t, \p_t)} + \underbrace{\begin{bmatrix}\Ov & \Ov\\ \Ov & \gamma\I\end{bmatrix}}_{=\G}\begin{bmatrix}\q_t\\ \p_t\end{bmatrix}.$$
</body>

i.e.,

.content-box-blue.center[<body> Eq. (B) is also a .text-shadow.red[conformal Hamiltonian system.] </body>]

---
class: left
count: false

# Observation #1: .orange.text-shadow[Negative friction accumulates energy]

---
class: left
count: false

# Observation #1: .orange.text-shadow[Negative friction accumulates energy]

.content-box-blue.center[<body>But as `$t \uparrow$`, the particle magically gains energy (हराम) </body>]

---
class: center, middle
count: false

# HaRAM
## / हराम /

.Large[Proscribed or forbidden by Law.]

---

# Hamiltonian Repelling-Attracting Metropolis

1. Choose a hypothetical friction parameter `$\gamma \in [0, \infty)$`, and integration time `$T$`
   
2. Let `$\Om = \begin{bmatrix}\boldsymbol{0} & I\\ -I & \boldsymbol{0}\end{bmatrix}$` and `$\G = \begin{bmatrix} \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \gamma \end{bmatrix}$`
   
3. For time `$t \in [0, T/2]$` generate .dblue[conformal Hamiltonian dynamics] using .red[negative friction]
`\begin{align}
\f{d}{d t} (\q_t, \p_t) = \Om \D H(\q_t, \p_t) + \G (\q_t, \p_t)
\end{align}`
--

1. For time `$t \in [T/2, T]$` generate .dblue[conformal Hamiltonian dynamics] using .green[positive friction]
`\begin{align}
\f{d}{d t} (\q_t, \p_t) = \Om \D H(\q_t, \p_t) - \G (\q_t, \p_t)
\end{align}`
--

1. Accept/reject state `$(\q_T, \p_T)$` with MH probability
<body>$$\alpha(\q_T, \p_T | \q_0, \p_0) = \min\pb{1, \f{\kappa(\q_0, \p_0 | \q_T, \p_T) e^{-H(\q_T,\p_T)}}{\kappa(\q_T, \p_T | \q_0, \p_0) e^{-H(\q_0,\p_0)}} \Bigg| \ddd{(\q_0, \p_0)}{(\q_t, \p_t)} \Bigg| }$$
</body>

---
count: false
# Hamiltonian Repelling Attracting Metropolis

---
count: false
# Hamiltonian Repelling Attracting Metropolis

---
count: false
# Hamiltonian Repelling Attracting Metropolis

---
class: left

# Properties of HaRAM

.left.content-box-purple[<body> 
For the Haram trajectory $\pb{(\q_t, \p_t) : t > 0}$$$\dd{t}H(\q_t, \p_t) \neq 0,$$but$$\int\limits_0^T\dd{t}H(\q_t, \p_t) dt = 0.$$
</body>]

.content-box-green.center[<body>
This implies that $H(\q_t, \p_t) = H(\q_0, \p_0)$ only when $t=T$.</body>]
]

.left.content-box-purple[<body> 
Let .red[<body>$\F_t$</body>] and .green[<body>$\P_t$</body>] be the conformal .bolder.text-shadow[Hamiltonian flow operators] w.r.t. the .red[repelling] and .green[attracting] dynamics, and let $\flip$ be the .bold[momentum-flip] operator. Then for all $t \le T/2$$$\flip \circ \P_t \circ \flip = \F_{-t}$$

In particular,$$\flip \circ \pa{\P_T \circ \F_T} \circ \flip \circ \pa{\P_T \circ \F_T} = \id.$$
</body>]

.content-box-green.center[<body>
This implies that $\kappa(\q_T, \p_T | \q_0, \p_0) = \onev$.</body>]
]
.pull-right[
 <img src="images/gifs/reversible.gif" height="270"/>
]

]

.left.content-box-purple[<body> 
.red[<body>$\F_t$ expands volume</body>] and .green[<body>$\P_t$ shrinks volume</body>]. But, over the entire time interval $[0,T]$:$$\ddd{(\q_0, \p_0)}{(\q_T, \p_T)} = 1.$$
</body>]

.content-box-green.center[<body>
This implies that the Jacobian correction term $\Bigg|\ddd{(\q_0, \p_0)}{(\q_T, \p_T)}\Bigg|$ is not needed.</body>]

]

.left.content-box-purple[<body>The deformation of the .bolder.text-shadow[symplectic 2-form] is given by
$$
 \omega(\q_t, \p_t) = \begin{cases}
 e^{\gamma t}\omega(\q_0, \p_0), & 0 \le t \le T/2\\ \\
 e^{\gamma (T-t)}\omega(\q_0, \p_0), & T/2 \le t \le T
 \end{cases}.
$$

In particular, HaRAM preserves the .bolder.text-shadow[symplectic 2-form], i.e.,$$\omega(\q_T, \p_T) = \omega(\q_0, \p_0).$$</body>]

]

For `$\e \approx dt$` and the .red[repelling] conformal Hamiltonian flow `$\F_\e \approx \F_{dt}$` the .red[conformal symplectic leapfrog algorithm] is

$$
`\begin{cases}
    \tp_{t+\f \e 2} = e^{\gamma\e/2}\p_t \\[5pt]
    \p_{t+\f \e 2} = \tp_{t+\f \e 2} - \f \e 2 \D U(\q_t) \\[5pt]
    \q_{t+\e} = \q_t + \e\S\inv \tp_{t+\f \e 2}\\[5pt]
    \tp_{t+\e} = \tp_{t+\f \e 2} - \f \e 2 \D U(\q_{t+\e})\\[5pt]
    \p_{t+\e} = e^{\gamma\e/2}\p_{t+\e}
\end{cases}`
$$
and for `$L = \lfloor T/\e \rfloor$`, `$\F_T \approx \F_{\e, L} = (\F_\e)^{\otimes L}$`. Similarly, `$\P_T \approx \P_{\e, L} = (\P_\e)^{\otimes L}$` by changing the sign of `$\gamma$`.

.content-box-purple.center[<body> `$|H\left( \P_{T} \circ \F_{T}(\q,\p) \right) - H\left( \P_{\e, L} \circ \F_{\e, L}(\q,\p) \right)| \approx O(\e^3)$`</body>]

]

---

# Experiments

---
class: left

# Bimodal Gaussian

.pull-left[
For `$\muv = 5 \onev_d \in \R^d$` the target is
$$
\pi(\q) \sim \half \mathcal{N}(\muv, \Sigma_1) + \half \mathcal{N}(-\muv, \Sigma_2)
$$
where
* `$\Sigma_1(i, j) = 0.5^{|i-j|}$` 
* `$\Sigma_2 = R \Sigma_1 R^\top$` for `$R \in SO(d)$`

.pull-right.center[
 <img src="images/plots/gaussian3/d/plt.svg">
]

| Method | Time/L(d=2) | Time/L(d=10) | Time/L(d=50) | W2 metric (d=2) | W2 metric (d=10) | W2 metric (d=50) |
|-------:|:-----------:|:------------:|:------------:|:---------------:|:----------------:|:----------------:|
|  HaRAM |   0.2697 s  |   0.7125 s   |   2.9816 s   |       0.51      |       7.93       |       7.78       |
|    HMC |   0.2555 s  |   0.6602 s   |   2.9375 s   |      17.16      |       45.93      |       48.04      |
|    RAM |   0.5364 s  |   0.5847 s   |   33.7778 s  |       0.9       |       33.88      |       51.52      |
|   MHMC |   0.5211 s  |   1.3530 s   |      ---     |      18.41      |       89.50      |       47.84      |

]

d=2
]]

d=10
]]

]

.pull-left.center[
 <img src="images/plots/gaussian3/d/plt2_tr.svg">
]

.pull-right.center[
 <img src="images/plots/gaussian3/d/plt10_tr.svg">
]
]

d=2
]

.pull-right.center[
 <img src="images/plots/gaussian3/d/scatter10.svg">

d=10
]
]

]

]

]
]

---
class: left

# Benchmark Dataset

]

.pull-left.center[

| Method |    Time/L   |    W2 metric    |
|-------:|:-----------:|:---------------:|
|  HaRAM |   0.2697 s  |      133.733    |
|    HMC |   0.2555 s  |      2989.539   |
|    RAM |   0.5364 s  |      309.421    |
|   MHMC |   0.5211 s  |      1916.740   |

]

.pull-right.center[
 <img src="figures/sims/gaussian3/acf-d2.svg">
]
]

]

<img src="figures/sims/gaussian3/scatter-d2.svg">
]
]

]

---
class: left

# Vanilla Gaussian Distribution

.pull-left.center[
 <img src="figures/sims/gaussian1/d2.svg">
]

.pull-right.center[
 <img src="figures/sims/gaussian1/acf-d2.svg">
]

---
class: left

# Autotuning

For `$\muv = \frac{5}{\sqrt{d}} \onev_d \in \R^d$` the target is
$$
\pi(\q) \sim \half \mathcal{N}(\muv, I) + \half \mathcal{N}(-\muv, I)
$$

]
]

---
class: left

# Bayesian Neural Network

.pull-right.center[
<img src="images/plots/nn/nn.png">
]

]
]

---
class: left

# Time Delay Estimation

$$
`\begin{align}
dX(t)&=-\frac{1}{\tau}(X(t)-\mu)dt+\sigma dB(t)\\
Y(t)&=X(t-\Delta)+\beta
\end{align}`
$$
]

$$
`\begin{align}
x_i\mid X(t_i) &\sim N(X(t_i),~ \sigma^2_{x, i})\\
y_i\mid Y(t_i) &\sim N(Y(t_i),~ \sigma^2_{y, i})~~\textrm{and}\\
y_i\mid X(t_i-\Delta), \Delta, \beta &\sim N(X(t_i-\Delta) + \beta,~ \sigma^2_{y, i}),
\end{align}`
$$

$$
`\begin{align}
\Delta &\sim \textrm{Unif}(-1178.939,~ 1179.939)\\
\beta &\sim \textrm{Unif}(-60,~ 60)\\
\mu &\sim \textrm{Unif}(-30,~ 30)\\
\sigma^2 &\sim \textrm{inverse-Gamma}(1, 2\times 10^{-7})\\
\tau &\sim \textrm{inverse-Gamma}(1, 1).
\end{align}`
$$

]

]

]
]

---
class: left

# Code

x = randn(1000)
y = x .+ randn(1000) .* e

@model function lr(x, y)
    σ2 ~ Truncated(Normal(), 1e-6, Inf)
    b0 ~ Cauchy(2.0)
    b1 ~ Normal(0.0, 1.0)
    y ~ MvNormal(b0 .+ b1 .* x, σ2 * I)
end

samples, accepts = mcmc(
  DualAverage(λ=10, δ=0.8)
  HaRAM();
  lr(x, y),
  n=1e4, n_burn=1e4
);
```
]

.content-box-blue.center[<body> github.com/sidv23/haram </body>]

]

```

---

# References

1. Arnold, V. I. (2013). *Mathematical methods of classical mechanics*
2. Betancourt, M. (2017). *A Conceptual Introduction to Hamiltonian Monte Carlo*
3. Betancourt, M., Byrne, S., Livingstone, S. & Girolami, M. (2017). *The geometric foundations of Hamiltonian monte carlo*
4. Duane, S., Kennedy, A. D., Pendleton, B. J. & Roweth, D. (1987). *Hybrid Monte Carlo*
5. Hairer, E., Lubich, C. & Wanner, G. (2006). *Geometric Numerical Integration*
6. McLachlan, R. & Perlmutter, M. (2001). *Conformal Hamiltonian systems*
7. Neal, R. (2011). *MCMC Using Hamiltonian Dynamics*
8. Tak, H., Meng, X.-L. & van Dyk, D. A. (2018). *A repelling-attracting Metropolis algorithm for multimodality.*
9. Tierney, L. (1994). *Markov chains for exploring posterior distributions.*
10. Tripuraneni, N., Rowland, M., Ghahramani, Z. & Turner, R. (2017). *Magnetic Hamiltonian Monte Carlo*
---

# Thanks for listening!
## Questions?