Hamiltonian Monte Carlo For Logistic Regression

Why Use Gradients to Explore a Logistic Posterior?

In a Bayesian clinical-trial model, regression coefficients can be strongly related in the posterior. A random-walk proposal that changes one coefficient at a time may then move only a short distance before it leaves a region of substantial posterior density. Hamiltonian Monte Carlo (HMC) uses the posterior’s gradient to guide proposals through that region. This can make a proposal travel much farther than a typical random-walk step while retaining a Metropolis correction that targets the desired posterior distribution.

The key calculation for logistic regression is the gradient of the log-posterior with respect to its coefficient vector. Once that gradient is available, HMC treats the coefficients as positions in a mechanical system, adds an auxiliary momentum, and follows a numerical approximation to the system’s trajectory. HMC is a sampling method, not a change to the likelihood or prior: its accepted draws still target the specified posterior.

Bayesian Logistic Regression and Its Gradient

Suppose there are $n$ conditionally independent binary outcomes $y_i\in\{0,1\}$, with fixed predictor vectors $x_i\in\mathbb{R}^d$. Let $\beta\in\mathbb{R}^d$ be the coefficient vector, and define the linear predictor $\eta_i=x_i^\top\beta$. The logistic model is

P(y_i=1\mid\beta)=p_i=\frac{1}{1+\exp(-\eta_i)}.

Assume a proper multivariate normal prior, $\beta\sim N(m,\Sigma)$, where $\Sigma$ is positive definite. Conditional independence gives the likelihood as the product of Bernoulli probabilities. Ignoring constants that do not depend on $\beta$, the log-posterior is

\ell(\beta) = \sum_{i=1}^n\left[y_i\eta_i-\log\{1+\exp(\eta_i)\}\right] -\frac{1}{2}(\beta-m)^\top\Sigma^{-1}(\beta-m).

Definition: The log-posterior gradient is the vector of partial derivatives of $\ell(\beta)$ with respect to the coefficients. For this model, $$ \nabla\ell(\beta)=X^\top(y-p)-\Sigma^{-1}(\beta-m), $$ where row $i$ of $X$ is $x_i^\top$, $y=(y_1,\ldots,y_n)^\top$, and $p=(p_1,\ldots,p_n)^\top$. The negative log-posterior, or potential energy, is $U(\beta)=-\ell(\beta)$, up to an arbitrary additive constant, and its gradient is $\nabla U(\beta)=X^\top(p-y)+\Sigma^{-1}(\beta-m)$.

The likelihood term has a simple interpretation: observation $i$ contributes $x_i(y_i-p_i)$ to the log-posterior gradient. The contribution points in the direction that would increase the fitted probability when the observed outcome exceeds it, and in the opposite direction when the fitted probability exceeds the observed outcome. The prior adds a term that pulls the coefficients toward $m$, with strength governed by $\Sigma^{-1}$.

To justify the formula, the derivative of $y_i\eta_i-\log\{1+\exp(\eta_i)\}$ with respect to $\eta_i$ is $y_i-p_i$, since the derivative of $\log\{1+\exp(\eta_i)\}$ is $p_i$. Because $\nabla_\beta\eta_i=x_i$, the chain rule gives $x_i(y_i-p_i)$. Summing these contributions gives $X^\top(y-p)$. Differentiating the normal-prior quadratic gives $-\Sigma^{-1}(\beta-m)$. Adding the likelihood and prior derivatives yields the stated gradient.

Worked Example: Calculating a Logistic Posterior Gradient

Consider a hypothetical three-patient dataset with an intercept and a numerical predictor $x$: the rows $(x_i,y_i)$ are $(0,0)$, $(1,1)$, and $(2,1)$. Let the coefficient vector be $\beta=(\beta_0,\beta_1)^\top$, with independent $N(0,1)$ priors, so $\Sigma=I$. Evaluate the log-posterior gradient at $\beta=(0,0)^\top$.

At this coefficient vector, every linear predictor is zero, so each fitted probability is $p_i=1/(1+\exp(0))=0.5$. The three residuals $y_i-p_i$ are $-0.5$, $0.5$, and $0.5$. The intercept component of $X^\top(y-p)$ is their sum, $-0.5+0.5+0.5=0.5$. The predictor component is $0(-0.5)+1(0.5)+2(0.5)=1.5$.

The prior gradient is $-\beta=(0,0)^\top$ at this point. Therefore,

\nabla\ell(0,0)^\top=(0.5,\ 1.5).

The positive components indicate that, at this point, increasing either the intercept or the predictor coefficient locally increases the log-posterior. This is a local derivative statement, not a claim that either coefficient has a positive posterior mean.

From a Posterior Gradient to Hamiltonian Dynamics

HMC introduces an auxiliary momentum vector $r\in\mathbb{R}^d$ with distribution $N(0,M)$, where $M$ is a positive-definite mass matrix. At each iteration, momentum is freshly drawn independently of the current coefficient vector. The joint target density for $(\beta,r)$ is proportional to the posterior density of $\beta$ times the momentum density. Its negative log-density, up to a constant, is the Hamiltonian

H(\beta,r)=U(\beta)+K(r), \qquad K(r)=\frac{1}{2}r^\top M^{-1}r.

Here $U$ is potential energy and $K$ is kinetic energy. Under exact Hamiltonian dynamics, the position and momentum evolve according to

\frac{d\beta}{dt}=M^{-1}r, \qquad \frac{dr}{dt}=-\nabla U(\beta).

The first equation moves the coefficients according to momentum; the second changes momentum according to the negative log-posterior gradient. Exact dynamics conserve $H$. Thus a trajectory can move a long way in coefficient space without systematically moving to lower posterior density: changes in potential energy are balanced by changes in kinetic energy.

For logistic regression, the potential is differentiable, and the formula above supplies the force needed in the second equation. HMC usually uses a numerical approximation rather than solving these differential equations exactly. A standard choice is leapfrog integration, which alternates half momentum steps and full position steps. For step size $\epsilon$, one leapfrog step is

\begin{aligned} r_{1/2}&=r_0-\frac{\epsilon}{2}\nabla U(\beta_0),\\ \beta_1&=\beta_0+\epsilon M^{-1}r_{1/2},\\ r_1&=r_{1/2}-\frac{\epsilon}{2}\nabla U(\beta_1). \end{aligned}

Repeated steps produce a proposed endpoint. The numerical trajectory does not conserve the Hamiltonian exactly, so HMC applies a Metropolis correction. If the initial and proposed states are $(\beta_0,r_0)$ and $(\beta_L,r_L)$, the acceptance probability is

\min\left\{1,\exp\big[-H(\beta_L,r_L)+H(\beta_0,r_0)\big]\right\}.

If the proposal is rejected, the coefficient vector remains at its current value. The leapfrog map is reversible when combined with momentum reversal and preserves volume; these properties, together with the acceptance correction, allow the method to preserve the target distribution. The momentum is auxiliary: posterior summaries concern the accepted coefficient draws, not the momentum draws.

Worked Example: One Leapfrog Step for a Hypothetical Trial

Suppose two hypothetical trial participants have predictor values $x_1=1$, $x_2=2$, and outcomes $y_1=1$, $y_2=0$. Use the model $\operatorname{logit}(p_i)=x_i\beta$ and prior $\beta\sim N(0,1)$. Start at $\beta_0=0$ with momentum $r_0=0.6$, mass $M=1$, and step size $\epsilon=0.1$. The initial fitted probabilities are both $0.5$. Thus the negative log-posterior gradient is

\nabla U(0)=\sum_i x_i(p_i-y_i)+0 =1(0.5-1)+2(0.5-0)=0.5.

The half momentum step is $r_{1/2}=0.6-(0.1/2)(0.5)=0.575$. The position step gives $\beta_1=0+0.1(0.575)=0.0575$. At this endpoint, the linear predictors are $0.0575$ and $0.115$, giving fitted probabilities $p_1\approx0.514371$ and $p_2\approx0.528718$.

The endpoint negative log-likelihood is

\log\{1+\exp(-0.0575)\}+\log\{1+\exp(0.115)\} \approx 1.417110.

The prior contribution to $U(\beta_1)$ is $0.0575^2/2=0.001653$. Therefore $U(\beta_1)\approx1.418763$. The endpoint gradient is $(p_1-1)+2p_2+\beta_1\approx0.629308$, so the final half momentum step gives $r_1=0.575-(0.05)(0.629308)\approx0.543535$.

Initially $U(\beta_0)=2\log 2\approx1.386294$, and the kinetic energy is $0.6^2/2=0.18$. Hence $H(\beta_0,r_0)\approx1.566294$. At the endpoint, $K(r_1)=0.543535^2/2\approx0.147715$, giving $H(\beta_1,r_1)\approx1.566478$. The acceptance probability is approximately $\exp[-(1.566478-1.566294)]\approx0.9998$. This small numerical energy increase makes rejection unlikely in this particular step; it does not establish that every proposal will be accepted.

Worked Example: Energy Error in a Simple Gaussian Target

To isolate the role of the Metropolis correction, take a hypothetical one-dimensional target with $U(q)=q^2/2$, mass $M=1$, initial position $q_0=1$, momentum $r_0=0$, and step size $\epsilon=0.2$. This is a standard normal target. The initial Hamiltonian is $H_0=1/2$. Since $\nabla U(q)=q$, the first half momentum update is $r_{1/2}=0-(0.1)(1)=-0.1$, and the position update is $q_1=1+(0.2)(-0.1)=0.98$.

The final half step gives $r_1=-0.1-(0.1)(0.98)=-0.198$. Thus $U(q_1)=0.98^2/2=0.4802$, $K(r_1)=0.198^2/2=0.019602$, and $H_1=0.499802$. The energy decreased by $0.000198$, so the acceptance probability is $\min\{1,\exp(0.000198)\}=1$. In other trajectories the numerical energy can increase, and the correction then gives a probability below one. The example shows why leapfrog integration is paired with an acceptance step even when its energy error is small.

Assumptions, Diagnostics, and Alternatives

The derivation assumes a correctly specified Bernoulli logistic likelihood conditional on the predictors, a differentiable log-posterior, and a proper target distribution. The stated normal prior is proper when $\Sigma$ is positive definite. For ordinary logistic regression with this prior, the posterior is proper because the likelihood is bounded and the prior is integrable. HMC also requires unconstrained coordinates for the numerical trajectory; a constrained parameter must first be represented on an unconstrained scale, with the appropriate change-of-variables adjustment to the log density when needed.

Before trusting a calculation, check the dimensions and signs: $X^\top(y-p)$ has one component per coefficient, and the prior term in the log-posterior gradient is negative when $\beta$ is above its prior mean. A useful numerical check is to compare a calculated directional derivative with a finite difference of the log-posterior for a small perturbation. If the likelihood code, coefficient ordering, or prior precision is wrong, the resulting force can send the trajectory in an inappropriate direction.

HMC can be inefficient or unreliable in practice if the trajectory is poorly scaled, if the posterior has strong curvature that the numerical steps do not resolve, or if the gradient is incorrect or unstable. Large energy errors and frequent rejections are warnings that the numerical proposals are not accurately following the target geometry. Very short trajectories can behave much like local proposals and fail to make useful progress. These signs call for checking the model and gradient and reconsidering the numerical setup; an acceptance probability alone is not proof that the chain explores all relevant posterior regions.

Approach	How proposals or updates are formed	Trade-off for logistic regression
Random-walk Metropolis	Adds a random perturbation to the current coefficient vector.	Simple and needs no gradient, but can make slow progress when coefficients are correlated.
Gibbs sampling	Draws parameters from their full conditional distributions.	Useful when those distributions can be sampled directly; a logistic model may not provide convenient direct conditional draws.
Hamiltonian Monte Carlo	Uses gradients and auxiliary momentum to generate a trajectory, then applies a Metropolis correction.	Can make distant proposals, but requires gradients and careful numerical integration.

The main benefit of HMC is not that it removes posterior correlation or guarantees rapid convergence. Rather, its gradient-guided trajectories can traverse a posterior region more effectively than an unstructured random walk. Its limitations include computational cost per proposal, dependence on smooth gradients, and sensitivity to the geometry and numerical resolution of the target. As with any MCMC method, accepted draws require appropriate assessment before they are used for posterior summaries.

Common Mistake

Common mistake: Treating the leapfrog endpoint as an exact draw from the posterior because it follows a gradient-informed path. Leapfrog integration only approximates Hamiltonian dynamics and generally has some energy error. The Metropolis correction is essential: it adjusts for that numerical error so the intended target is preserved under the method’s reversibility and volume-preservation conditions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For a Bernoulli logistic model with a normal prior, write the log-posterior gradient and identify the likelihood and prior contributions.
In the two-participant example above, why is the negative log-posterior gradient at $\beta=0$ equal to $0.5$, rather than $-0.5$?
What is the role of momentum in HMC, and why is the Metropolis correction still needed after leapfrog integration?
Name two conditions or checks that matter when applying HMC to logistic regression.

Solutions. (1) The gradient is $X^\top(y-p)-\Sigma^{-1}(\beta-m)$; the first term is from the Bernoulli likelihood and the second from the normal prior. (2) The log-posterior gradient is $-0.5$ at zero, but the negative log-posterior gradient has the opposite sign, $+0.5$. (3) Momentum supplies a direction and speed for movement through coefficient space, allowing a trajectory to travel beyond a local random-walk step. Leapfrog integration is approximate, so the Metropolis correction accounts for the endpoint’s Hamiltonian error. (4) Examples include using a proper, normalizable posterior and a differentiable log density; checking gradient signs and dimensions; and checking for large energy errors or frequent rejections.