ZigZagBoomerang

Boomerang(μ, λ) <: ContinuousDynamics

Dynamics preserving the N(μ, Σ) measure (Boomerang) with refreshment time λ

Boomerang1d(Σ, μ, λ) <: ContinuousDynamics

1-d toy boomerang samper. Dynamics preserving the N(μ, Σ) measure with refreshment time λ.

BouncyParticle(λ) <: ContinuousDynamics

Input: argument Γ, a sparse precision matrix approximating target precision. Bouncy particle sampler, λ is the refreshment rate, which has to be strictly positive.

ContinuousDynamics

Abstract type for the deterministic dynamics of PDMPs

ExtendedForm()

Indicates as args[1] that ∇ϕ depends on the extended arguments

∇ϕ(t, x, θ, i, t′, Z, args...)

instead of

∇ϕ(x, i, args...)

Can be used to implement ∇ϕ depending on random coefficients.

FactBoomerang(Γ, μ, λ) <: ContinuousDynamics

Factorized Boomerang dynamics preserving the N(μ, inv(Diagonal(Γ))) measure with refreshment time λ.

Exploits the conditional independence structure of the target measure, in form the argument Γ, a sparse precision matrix approximating target precision. μ is the approximate target mean.

FactTrace

See Trace.

PDMPTrace

See Trace.

SelfMoving()

Indicates as args[1] that ∇ϕ depends only on few coeffients and takes responsibility to call smove_forward!.

Replaced by ExtendedForm.

Trace(t0::T, x0, θ0, F::Union{ZigZag,FactBoomerang})

Trace object for exact trajectory of pdmp samplers. Returns an iterable FactTrace object. Note that iteration iterates pairs t => x where the vector x is modified inplace, so copies have to be made if the x is to be saved. collect applied to a trace object automatically copies x. discretize returns a discretized version.

struct ZigZag(Γ, μ) <: ContinuousDynamics

Local ZigZag sampler which exploits any independence structure of the target measure, in form the argument Γ, a sparse precision matrix approximating target precision. μ is the approximate target mean.

ZigZag1d <: ContinuousDynamics

1-d toy ZigZag sampler, dynamics preserving the Lebesgue measure.

ab(G, i, x, θ, c, Flow)

Returns the constant term a and linear term b when computing the Poisson times from the upper upper bounding rates λᵢ(t) = max(a + b*t)^2. The factors a and b can be function of the current position x, velocity θ, tuning parameter c and the Graph G

conditional_trace(trace::Trace, (p, n))

ConditionalTrace trace on hyperplane H through p with norml n. Returns iterable object iterating pairs t => x such that x ∈ H.

Iteration changes the vector x inplace, collect creates necessary copies.

discretize(trace::Trace, dt)

Discretize trace with step-size dt. Returns iterable object iterating pairs t => x.

Iteration changes the vector x inplace, collect creates necessary copies.

discretize(x::Vector, Flow::Union{ZigZag1d, Boomerang1d}, dt)

Transform the output of the algorithm (a skeleton of points) to a trajectory. Simple 1-d version.

τ = freezing_time(x, θ)

computes the hitting time of a 1d particle with constant velocity θ to hit 0 given the position x

idot(A, j, x) = dot(A[:, j], x)

Compute column-vector dot product exploiting sparsity of A.

idot_moving!(A::SparseMatrixCSC, j, t, x, θ, t′, F)

Compute column-vector dot product exploiting sparsity of A. Move all coordinates needed to their position at time t′

move_forward!(τ, t, x, θ, B::Boomerang)

Updates the position x, velocity θ and time t of the process after a time step equal to τ according to the deterministic dynamics of the Boomerang sampler which are the Hamiltonian dynamics preserving the Gaussian measure: : xt = μ +(x0 − μ)cos(t) + v_0sin(t), vt = −(x0 − μ)sin(t) + v_0cos(t) x: current location, θ: current velocity, t: current time.

move_forward!(τ, t, x, θ, Z::Union{BouncyParticle, ZigZag})

Updates the position x, velocity θ and time t of the process after a time step equal to τ according to the deterministic dynamics of the Bouncy particle sampler (BouncyParticle) and ZigZag: (x(τ), θ(τ)) = (x(0) + θ(0)*t, θ(0)). x: current location, θ: current velocity, t: current time,

move_forward(τ, t, x, θ, B::Boomerang1d)

Updates the position x, velocity θ and time t of the process after a time step equal to τ according to the deterministic dynamics of the Boomerang1d sampler: xt = μ +(x0 − μ)cos(t) + v_0sin(t), vt = −(x0 − μ)sin(t) + v_0cos(t) x: current location, θ: current velocity, t: current time.

move_forward(τ, t, x, θ, ::ZigZag1d)

Updates the position x, velocity θ and time t of the process after a time step equal to τ according to the deterministic dynamics of the ZigZag1d sampler: (x(τ), θ(τ)) = (x(0) + θ(0)*t, θ(0)). x: current location, θ: current velocity, t: current time,

neighbours(G::Vector{<:Pair}, i) = G[i].second

Return extended neighbourhood of i including i. G: graphs of neightbourhoods

normsq(x)

Squared 2-norm.

oscn!(rng, v, ∇ψx, ρ; normalize=false)

Orthogonal subspace Crank-Nicolson step with autocorrelation ρ for standard Gaussian or Uniform on the sphere (normalize = true).

pdmp(∇ϕ, t0, x0, θ0, T, c, F::Union{ZigZag, FactBoomerang}, args..., args) = Ξ, (t, x, θ), (acc, num), c

Outer loop of the factorised samplers, the Factorised Boomerang algorithm and the Zig-Zag sampler. Inputs are a function ∇ϕ giving ith element of gradient of negative log target density ∇ϕ(x, i, args...), starting time and position t0, x0, velocities θ0, and tuning vector c for rejection bounds and final clock T.

The process moves to time T with invariant mesure μ(dx) ∝ exp(-ϕ(x))dx and outputs a collection of reflection points which, together with the initial triple t, x θ are sufficient for reconstructuing continuously the continuous path. It returns a FactTrace (see Trace) object Ξ, which can be collected into pairs t => x of times and locations and discretized with discretize. Also returns the number of total and accepted Poisson events and updated bounds c (in case of adapt==true the bounds are multiplied by factor if they turn out to be too small.)

This version does not assume that ∇ϕ has sparse conditional dependencies, see spdmp.

pdmp(dϕ, ∇ϕ!, t0, x0, θ0, T, c::Bound, flow::BouncyParticle, args...; oscn=false, adapt=false, progress=false, progress_stops = 20, islocal = false, seed=Seed(), factor=2.0)

The first directional derivative dϕ[1] tells me if I move up or down the potential. The second directional derivative dϕ[2] tells me how fast the first changes. The gradient ∇ϕ! tells me the surface I want to reflect on. Refreshes proportional to speed. Keeps only samples at intervals proportional to those times.

 dϕ = function (t, x, v, args...) # two directional derivatives
     u = ForwardDiff.derivative(t -> -ℓ(x + t*v), Dual{:hSrkahPmmC}(0.0, 1.0))
     u.value, u.partials[]
 end
 ∇ϕ! = function (y, t, x, args...)
     ForwardDiff.gradient!(y, ℓ, x)
     y .= -y
     y
 end

The remaining arguments:

d = 25 # number of parameters 
t0 = 0.0
x0 = zeros(d) # starting point sampler
θ0 = randn(d) # starting direction sampler
T = 200. # end time (similar to number of samples in MCMC)
c = 50.0 # initial guess for the bound

# define BouncyParticle sampler (has two relevant parameters) 
Z = BouncyParticle(∅, ∅, # information about graphical structure
    10.0, # momentum refreshment rate 
    0.95, # momentum correlation / only gradually change momentum in refreshment/momentum update
    nothing, # PDMats compatible inverse mass OR
    I # left cholesky factor of momentum precision
) 

trace, final, (acc, num), cs = @time pdmp(
        dneglogp, # return first two directional derivatives of negative target log-likelihood in direction v
        ∇neglogp!, # return gradient of negative target log-likelihood
        t0, x0, θ0, #initial state 
        T, # duration (Real) or number of samples (Int)
        ZZB.LocalBound(c), # use Hessian information 
        Z; # sampler
        oscn=false, # no orthogonal subspace pCR
        adapt=true, # adapt bound c
        adapt_mass=false # adapt PDiag U matrix to fisher information estimate
        progress=true, # show progress bar
)

# to obtain direction change times and points of piecewise linear trace
t, x = ZigZagBoomerang.sep(trace)

pdmp(∇ϕ!, t0, x0, θ0, T, c::Bound, Flow::Union{BouncyParticle, Boomerang}; adapt=false, factor=2.0)

Run a Bouncy particle sampler (BouncyParticle) or Boomerang sampler from time, location and velocity t0, x0, θ0 until time T. ∇ϕ!(y, x) writes the gradient of the potential (neg. log density) into y. c is a tuning parameter for the upper bound of the Poisson rate. If adapt = false, c = c*factor is tried, otherwise an error is thrown.

It returns a PDMPTrace (see Trace) object Ξ, which can be collected into pairs t => x of times and locations and discretized with discretize. Also returns the number of total and accepted Poisson events and updated bounds c (in case of adapt==true the bounds are multiplied by factor if they turn out to be too small.)

pdmp(∇ϕ, x, θ, T, Flow::ContinuousDynamics; adapt=true,  factor=2.0)

Run a piecewise deterministic process from location and velocity x, θ until time T. c is a tuning parameter for the upper bound of the Poisson rate. If adapt = false, c = c*factor is tried, otherwise an error is thrown.

Returns vector of tuples (t, x, θ) (time, location, velocity) of direction change events.

poisson_time(a, b, u)

Obtaining waiting time for inhomogeneous Poisson Process with rate of the form λ(t) = (a + b*t)^+, a,b ∈ R, u uniform random variable

poisson_time(a[, u])

Obtaining waiting time for homogeneous Poisson Process with rate of the form λ(t) = a, a ≥ 0, u uniform random variable

poisson_time((a, b, c), u)

Obtaining waiting time for inhomogeneous Poisson Process with rate of the form λ(t) = c + (a + b*t)^+, where c> 0 ,a, b ∈ R, u uniform random variable

pos(x)

Positive part of x (i.e. max(0,x)).

queue_time!(Q, t, x, θ, i, b, f, Z::ZigZag)

Computes the (proposed) reflection time and the freezing time of the ith coordinate and enqueue the first one. f[i] = true if the next time is a freezing time.

reflect!(∇ϕx, x, θ, F::BouncyParticle, Boomerang)

Reflection rule of sampler F at reflection time. x: position, θ: velocity

    reflect!(i, ∇ϕx, x, θ, F)

Reflection rule of sampler F at reflection time, ∇ϕx: inner product of ∇ϕ and x.i: coordinate which flips sign,x: position,θ: velocity (position and inner product not used for theZigZagandFactBoomerang`.)

discretize(x::Vector, Flow::Union{BouncyParticle, Boomerang}, dt)

Transform the output of the algorithm (a skeleton of points) to a trajectory. multi-dimensional version.

Old version that would not work with the sticky Boomerang sampler not centered in 0

spdmp(∇ϕ, t0, x0, θ0, T, c, [G,] F::Union{ZigZag,FactBoomerang}, args...;
    factor=1.5, adapt=false)
    = Ξ, (t, x, θ), (acc, num), c

Version of spdmp which assumes that i only depends on coordinates x[j] for j in neighbours(G, i).

It returns a FactTrace (see Trace) object Ξ, which can be collected into pairs t => x of times and locations and discretized with discretize. Also returns the number of total and accepted Poisson events and updated bounds c (in case of adapt==true the bounds are multiplied by factor if they turn out to be too small.) The final time, location and momentum at T can be obtained with smove_forward!(t, x, θ, T, F).

spdmp(∇ϕ, t0, x0, θ0, T, c, [G,] F::Union{ZigZag,FactBoomerang}, args...; factor=1.5, adapt=false) = Ξ, (t, x, θ), (acc, num), c

Version of spdmp which assumes that i only depends on coordinates x[j] for j in neighbours(G, i).

It returns a FactTrace (see Trace) object Ξ, which can be collected into pairs t => x of times and locations and discretized with discretize. Also returns the number of total and accepted Poisson events and updated bounds c (in case of adapt==true the bounds are multiplied by factor if they turn out to be too small.) The final time, location and momentum at T can be obtained with smove_forward!(t, x, θ, T, F).

splitpairs(tx) = t, x

Splits a vector of pairs into a pair of vectors.

t, x, θ = ssmove_forward!(G, i, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})

moves forward only the non_frozen particles neighbours of i

t, x, θ = ssmove_forward!(t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})

moves forward only the non_frozen particles

sspdmp_inner!(Ξ, G, G1, G2, ∇ϕ, t, x, θ, Q, c, b, t_old, f, θf, (acc, num),
        F::ZigZag, κ, args...; strong_upperbounds = false, factor=1.5, adapt=false)

Inner loop of the sticky ZigZag sampler. G[i] are indices which have to be moved, G1[i] is the set of indices used to derive the bounding rate λbari and G2 are the indices k in Aj for all j : i in Aj (neighbours of neighbours)

Is assumed that ∇ϕ[x, i] is function of x_i with i in G[i] or that ∇ϕ takes care of moving .

subtrace(tr, J)

Compute the trace of a subvector x[J], returns a trace object.

λ(∇ϕi, i, x, θ, Z::FactBoomerang)

ith Poisson rate of the FactBoomerang sampler

λ(∇ϕ, i, x, θ, Z::ZigZag)

ith Poisson rate of the ZigZag sampler