ZigZagBoomerang
Back to repository: https://github.com/mschauer/ZigZagBoomerang.jl
ZigZagBoomerang.Boomerang
ZigZagBoomerang.Boomerang1d
ZigZagBoomerang.BouncyParticle
ZigZagBoomerang.ContinuousDynamics
ZigZagBoomerang.ExtendedForm
ZigZagBoomerang.FactBoomerang
ZigZagBoomerang.FactTrace
ZigZagBoomerang.PDMPTrace
ZigZagBoomerang.SelfMoving
ZigZagBoomerang.Trace
ZigZagBoomerang.ZigZag
ZigZagBoomerang.ZigZag1d
ZigZagBoomerang.ab
ZigZagBoomerang.conditional_trace
ZigZagBoomerang.discretize
ZigZagBoomerang.discretize
ZigZagBoomerang.freezing_time
ZigZagBoomerang.idot
ZigZagBoomerang.idot_moving!
ZigZagBoomerang.move_forward
ZigZagBoomerang.move_forward
ZigZagBoomerang.move_forward!
ZigZagBoomerang.move_forward!
ZigZagBoomerang.neighbours
ZigZagBoomerang.normsq
ZigZagBoomerang.oscn!
ZigZagBoomerang.pdmp
ZigZagBoomerang.pdmp
ZigZagBoomerang.pdmp
ZigZagBoomerang.pdmp
ZigZagBoomerang.poisson_time
ZigZagBoomerang.poisson_time
ZigZagBoomerang.poisson_time
ZigZagBoomerang.pos
ZigZagBoomerang.queue_time!
ZigZagBoomerang.reflect!
ZigZagBoomerang.reflect!
ZigZagBoomerang.sdiscretize
ZigZagBoomerang.spdmp
ZigZagBoomerang.spdmp
ZigZagBoomerang.splitpairs
ZigZagBoomerang.ssmove_forward!
ZigZagBoomerang.ssmove_forward!
ZigZagBoomerang.sspdmp_inner!
ZigZagBoomerang.subtrace
ZigZagBoomerang.λ
ZigZagBoomerang.λ
ZigZagBoomerang.Boomerang
— TypeBoomerang(μ, λ) <: ContinuousDynamics
Dynamics preserving the N(μ, Σ)
measure (Boomerang) with refreshment time λ
ZigZagBoomerang.Boomerang1d
— TypeBoomerang1d(Σ, μ, λ) <: ContinuousDynamics
1-d toy boomerang samper. Dynamics preserving the N(μ, Σ)
measure with refreshment time λ
.
ZigZagBoomerang.BouncyParticle
— TypeBouncyParticle(λ) <: ContinuousDynamics
Input: argument Γ
, a sparse precision matrix approximating target precision. Bouncy particle sampler, λ
is the refreshment rate, which has to be strictly positive.
ZigZagBoomerang.ContinuousDynamics
— TypeContinuousDynamics
Abstract type for the deterministic dynamics of PDMPs
ZigZagBoomerang.ExtendedForm
— TypeExtendedForm()
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.
ZigZagBoomerang.FactBoomerang
— TypeFactBoomerang(Γ, μ, λ) <: 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.
ZigZagBoomerang.FactTrace
— TypeFactTrace
See Trace
.
ZigZagBoomerang.PDMPTrace
— TypePDMPTrace
See Trace
.
ZigZagBoomerang.SelfMoving
— TypeSelfMoving()
Indicates as args[1]
that ∇ϕ
depends only on few coeffients and takes responsibility to call smove_forward!
.
Replaced by ExtendedForm
.
ZigZagBoomerang.Trace
— MethodTrace(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.
ZigZagBoomerang.ZigZag
— Typestruct 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.
ZigZagBoomerang.ZigZag1d
— TypeZigZag1d <: ContinuousDynamics
1-d toy ZigZag sampler, dynamics preserving the Lebesgue measure.
ZigZagBoomerang.ab
— Methodab(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
ZigZagBoomerang.conditional_trace
— Methodconditional_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.
ZigZagBoomerang.discretize
— Methoddiscretize(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.
ZigZagBoomerang.discretize
— Methoddiscretize(x::Vector, Flow::Union{ZigZag1d, Boomerang1d}, dt)
Transform the output of the algorithm (a skeleton of points) to a trajectory. Simple 1-d version.
ZigZagBoomerang.freezing_time
— Methodτ = freezing_time(x, θ)
computes the hitting time of a 1d particle with constant velocity θ
to hit 0 given the position x
ZigZagBoomerang.idot
— Methodidot(A, j, x) = dot(A[:, j], x)
Compute column-vector dot product exploiting sparsity of A
.
ZigZagBoomerang.idot_moving!
— Methodidot_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′
ZigZagBoomerang.move_forward!
— Methodmove_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.
ZigZagBoomerang.move_forward!
— Methodmove_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,
ZigZagBoomerang.move_forward
— Methodmove_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.
ZigZagBoomerang.move_forward
— Methodmove_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,
ZigZagBoomerang.neighbours
— Methodneighbours(G::Vector{<:Pair}, i) = G[i].second
Return extended neighbourhood of i
including i
. G
: graphs of neightbourhoods
ZigZagBoomerang.normsq
— Methodnormsq(x)
Squared 2-norm.
ZigZagBoomerang.oscn!
— Methodoscn!(rng, v, ∇ψx, ρ; normalize=false)
Orthogonal subspace Crank-Nicolson step with autocorrelation ρ
for standard Gaussian or Uniform on the sphere (normalize = true
).
ZigZagBoomerang.pdmp
— Methodpdmp(∇ϕ, 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 i
th 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 num
ber of total and acc
epted 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
.
ZigZagBoomerang.pdmp
— Methodpdmp(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)
ZigZagBoomerang.pdmp
— Methodpdmp(∇ϕ!, 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 num
ber of total and acc
epted 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.)
ZigZagBoomerang.pdmp
— Methodpdmp(∇ϕ, 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.
ZigZagBoomerang.poisson_time
— Methodpoisson_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
ZigZagBoomerang.poisson_time
— Methodpoisson_time(a[, u])
Obtaining waiting time for homogeneous Poisson Process with rate of the form λ(t) = a, a
≥ 0, u
uniform random variable
ZigZagBoomerang.poisson_time
— Methodpoisson_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
ZigZagBoomerang.pos
— Methodpos(x)
Positive part of x
(i.e. max(0,x)).
ZigZagBoomerang.queue_time!
— Methodqueue_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.
ZigZagBoomerang.reflect!
— Methodreflect!(∇ϕx, x, θ, F::BouncyParticle, Boomerang)
Reflection rule of sampler F
at reflection time. x
: position, θ
: velocity
ZigZagBoomerang.reflect!
— Method 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 the
ZigZagand
FactBoomerang`.)
ZigZagBoomerang.sdiscretize
— Methoddiscretize(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
ZigZagBoomerang.spdmp
— Methodspdmp(∇ϕ, 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 num
ber of total and acc
epted 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)
.
ZigZagBoomerang.spdmp
— Methodspdmp(∇ϕ, 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 num
ber of total and acc
epted 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)
.
ZigZagBoomerang.splitpairs
— Methodsplitpairs(tx) = t, x
Splits a vector of pairs into a pair of vectors.
ZigZagBoomerang.ssmove_forward!
— Methodt, x, θ = ssmove_forward!(G, i, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})
moves forward only the non_frozen particles neighbours of i
ZigZagBoomerang.ssmove_forward!
— Methodt, x, θ = ssmove_forward!(t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})
moves forward only the non_frozen particles
ZigZagBoomerang.sspdmp_inner!
— Methodsspdmp_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 .
ZigZagBoomerang.subtrace
— Methodsubtrace(tr, J)
Compute the trace of a subvector x[J]
, returns a trace object.
ZigZagBoomerang.λ
— Methodλ(∇ϕi, i, x, θ, Z::FactBoomerang)
i
th Poisson rate of the FactBoomerang
sampler
ZigZagBoomerang.λ
— Methodλ(∇ϕ, i, x, θ, Z::ZigZag)
i
th Poisson rate of the ZigZag
sampler
Literature
- Joris Bierkens, Paul Fearnhead, Gareth Roberts: The Zig-Zag Process and Super-Efficient Sampling for Bayesian Analysis of Big Data. The Annals of Statistics, 2019, 47. Vol., Nr. 3, pp. 1288-1320. https://arxiv.org/abs/1607.03188.
- Joris Bierkens, Sebastiano Grazzi, Kengo Kamatani and Gareth Robers: The Boomerang Sampler. ICML 2020. https://arxiv.org/abs/2006.13777.
- Joris Bierkens, Sebastiano Grazzi, Frank van der Meulen, Moritz Schauer: A piecewise deterministic Monte Carlo method for diffusion bridges. 2020. https://arxiv.org/abs/2001.05889.
- https://github.com/jbierkens/ICML-boomerang/ (code accompanying the paper "The Boomerang Sampler")