Hyperplane Disk Density¶
In this experiment, we consider the target distribution $\nu$ given by (unnormalized) density
$$ \varrho_{\nu}(x) = \exp\left( - \left( \textstyle \sum_{i=1}^d x_i \right)^2 - \Vert x \Vert^2 \right) , \qquad x = (x_1,...,x_d)^T \in \mathbb{R}^d , $$
in dimension $d=200$. It is easily seen that the points $x$ at which the sum term vanishes constitute a hyperplane of $\mathbb{R}^d$, i.e.~a $(d-1)$-dimensional subspace. Moreover, on rays emanating orthgonally from the hyperplane, the sum term increases like the squared Euclidean norm. It therefore leads $\nu$ to be concentrated around the hyperplane with Gaussian tails. By itself, this would not yet constitute an unnormalized density (i.e. an integrable function), so the second term is introduced to further concentrate $\nu$ on a region near the coordinate origin, with Gaussian tails on all rays emanating from the origin and running within the hyperplane. Together, the two terms produce a distribution that is concentrated on a zero-centered spherical subset of the hyperplane, which can be viewed as a (hyper-)disk.
This target serves as an example of a distribution with highly non-trivial covariance structure. Accordingly, GPSS outperforming its competitors here is meant to demonstrate that it is relatively robust against deviations of the target covariance from the identity. We compute the usual performance metrics and generate radius trace plots and step size histograms below.
import sys
try:
if "pyodide" in sys.modules:
import piplite
await piplite.install('fastprogress')
except:
print("piplite not found, ignore this if you are not using Jupyterlite")
import scipy
import slice_sampling as ss
import mcmc_utils as mcu
import plotting_functions as pfs
import numpy as np
import numpy.linalg as alg
import pandas as pd
# Set parameters
d = 200 # dimension, should be > 1
itnum = int(1e4) # number of chain iterations
path_prefix = "./plots/hyperplane_disk_"
Construct the Target¶
tde_cnt = 0 # target density evaluation count
def log_density(x):
global tde_cnt
tde_cnt += 1
return -alg.norm(x)**2 - np.sum(x)**2
Run the Samplers¶
# initialize chains in an area of high probability mass to avoid burn-in
x_0 = np.ones(d)
x_0[-1] = 1-d
x_0 *= 10.0 / alg.norm(x_0) # WARNING: the 10.0 is hand-adjusted to dim d=200, for d=100 it should be 7.0
w = 20.0 # WARNING: also hand-adjusted to dim d=200
samples_gpss = ss.gibbsian_polar_ss(log_density, x_0, w, itnum)
tdes_gpss = tde_cnt
tde_cnt = 0
samples_hruss = ss.hit_and_run_uniform_ss(log_density, x_0, w, itnum)
tdes_hruss = tde_cnt
tde_cnt = 0
samples_ess_u = ss.naive_elliptical_ss(log_density, x_0, itnum)
tdes_ess_u = tde_cnt
tde_cnt = 0
# compute empirical covariance to feed to ESS for tuned performance
from sklearn.covariance import EmpiricalCovariance
tuned_cov = EmpiricalCovariance().fit(samples_gpss).covariance_
print(tuned_cov[0,0], np.sum(tuned_cov[0,1:]))
# near-perfectly tuned ESS
log_likelihood = lambda x: -ss.log_prior(tuned_cov, x) + log_density(x)
samples_ess_t = ss.elliptical_ss(tuned_cov, log_likelihood, x_0, itnum)
tdes_ess_t = tde_cnt
tde_cnt = 0
Collect Outputs, Compute Performance Metrics¶
snames = ["Gibbsian Polar Slice Sampling", "Hit-and-Run Uniform Slice Sampling", "Untuned Elliptical Slice Sampling", "Tuned Elliptical Slice Sampling"]
samples = [samples_gpss, samples_hruss, samples_ess_u, samples_ess_t]
tdes_per_it = np.array([tdes_gpss, tdes_hruss, tdes_ess_u, tdes_ess_t]) / itnum # target density evals per iteration
msamples = [sams[:,0] for sams in samples]
radii = mcu.get_radii_list(samples)
steps = np.array(mcu.get_steps_list(samples))
iats = mcu.iat_list(radii)
pd.DataFrame(
np.array([snames, np.round(tdes_per_it, 2), np.round(iats, 2), np.round(np.mean(steps, axis=1), 2)]).T,
columns = ["Sampler", "TDEs/it", "IAT", "MSS"]
)
Plotting¶
pfs.plot_trace_and_step_hists(radii, steps, snames, dpi=250, filepath = path_prefix + "radii_and_steps.png")
