# \[ANN\] ZigZagBoomerang.jl

**URL:** <https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287>\
**Category:** Probabilistic Programming\
**Tags:** package, monte-carlo\
**Created:** [March 16, 2021, 9:53am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287 "2021-03-16T09:53:06Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 16, 2021, 9:53am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/1 "2021-03-16T09:53:06Z")

</div>

Dear all!

We would like to announce our package together with our new paper “Sticky PDMP samplers for sparse and local inference problems” ([Arxiv link](https://arxiv.org/abs/2103.08478))

# [ZigZagBoomerang.jl](https://github.com/mschauer/ZigZagBoomerang.jl)

_Sleek implementations of the ZigZag, Boomerang and other assorted piecewise deterministic Markov processes for Markov Chain Monte Carlo including Sticky PDMPs for variable selection_

## 1. What is this about

ZigZagBoomerang.jl contains a user-friendly and clean implementations of the most prominent piecewise deterministic Monte Carlo methods. These Monte Carlo methods are fast (non-reversible, with momentum) and allow for subsampling of data without creating bias and can be used for sparse and local inference problems.

## 2. What can it do

Within a few lines of code, you can run challenging problems with large data size and in high dimensional space. The only ingredient the user has to add is (the gradient of ) the target log-density and perhaps a rough estimate of the posterior mean.

For example the gradient of a Gaussian density is linear and can be given by a sparse matrix `Γ` and then it is

```julia
# 
Γ = ... # sparse precision matrix

# Define ∇ϕ(x, i, Γ) giving the partial derivative of ϕ(x) with respect to x[i]
∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x) # more efficient that dot(Γ[:, i], x)

# Random initial values
t0 = 0.0
x0 = randn(n*n)
θ0 = rand([-1.0,1.0], n*n)

# Rejection bounds
c = 1.0*ones(length(x0))

# Define ZigZag
Z = ZigZag(Γ, x0*0)

# Run sparse ZigZag for T time units and collect trajectory
T = 20.0
@time trace, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ; adapt=true)
@time traj = collect(discretize(trace, 0.1))

```

![zigzagfield](https://global.discourse-cdn.com/julialang/original/3X/6/a/6a6a92ba3c10835547d1701230855d633fd9a555.gif)

## 4. Using ZigZag and friends as backends in your probabilistic programming environment.

We have written ZigZagBoomerang.jl to be used as sampling engine for probabilistic programming languages. Let’s demo this with [Soss.jl](https://github.com/cscherrer/Soss.jl) where work on Zig-Zag integration is already under way:

```julia-auto
using Soss
m = @model x begin
    α ~ Normal()
    β ~ Normal()
    yhat = α .+ β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], 2.0)
    end
end

# generate data from model
x = randn(20);
obs = -0.1 .+ 2x + 1randn(20); 
posterior = m(x=x) | (y=obs,)

```

```julia-auto
using ZigZagBoomerang

include("zigzag.jl") 

T = 100.0
trace, final, (num, acc) = @time zigzag(posterior, T) # sample Soss posterior

# trace is a continous object, discretize to obtain samples
ts, xs = ZigZagBoomerang.sep(discretize(trace, 0.1))

ts # time points
xs # vector of points

```

(see also [zigzag.jl](https://gist.github.com/mschauer/638f2f50abd6469cec1c38e561374d24))  
Let’s plot

```julia
using Plots
p = plot(ts, getindex.(xs, 1))
plot!(ts, getindex.(xs, 2), color=:red)
p2 = plot(first.(xs), last.(xs))
p

```

 ![tracezigzag](https://global.discourse-cdn.com/julialang/original/3X/5/2/5273618cac41b518d3113e0ec664e579d1ac26a9.png) ![tracezigzag2](https://global.discourse-cdn.com/julialang/original/3X/3/b/3bc248803ea3a1b7443beb9c6e94e719a256206c.png)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 16, 2021, 10:13am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/2 "2021-03-16T10:13:14Z")

</div>

Cool! The arxiv link is not working

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 16, 2021, 10:26am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/3 "2021-03-16T10:26:54Z")

</div>

## So what about sparsity?

The Zig-Zag allows (that is what our paper is about) spike and slab priors for variable selection. Let’s demo this with Soss:

Let’s use MeasureTheory.jl 's `SpikeMixture` (a very spiky spike and slab with a Dirac measure as spike) to indicate that we have ground to believe the coefficients are actually zero:

```julia
m = @model x begin
    α ~ SpikeMixture(Normal(), 0.2) # 0.2*Normal() + 0.8*Dirac(0)
    β ~ SpikeMixture(Normal(), 0.2)
    yhat = α .+ β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], 2.0)
    end
end

x = randn(20);
obs = -0.1 .+ 2x + 1randn(20); 
T = 100.0
posterior = m(x=x) | (y=obs,)

```

sampling with the sticky ZigZag

```julia

include("sparsezigzag.jl")
trace, final, (num, acc) = @time sparse_zigzag(posterior, T, c=50)

ts, xs = ZigZagBoomerang.sep(discretize(trace, 0.1))
xs = xform(posterior).(xs) # make named tuples `NamedTuple{(:β, :α)`

```

now we can answer for example what is the posterior probability that the parameters are zero

```julia-auto
julia> mean(getindex.(xs, :α) .== 0)
0.7211394302848576

julia> mean(getindex.(xs, :β) .== 0)
0.028985507246376812

```

 ![tracesticky](https://global.discourse-cdn.com/julialang/original/3X/9/f/9f9b3ec7c5b649b510259c8334de1b0846c4554f.png)

The nice thing is that this scales very well…

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 16, 2021, 10:36am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/4 "2021-03-16T10:36:16Z")

</div>

For example a cloud of SDE driven _boids_

[![](https://global.discourse-cdn.com/julialang/original/3X/7/1/7116b9e1828a831a67c0ca7bf3e2ce8bd0d9e2d1.jpeg "Animation of boids") ](https://www.youtube.com/watch?v=O1VoURPwVLI)

Each boid is in love with another one and follows them… but whom? It also hates one and avoid it… which one we don’t know. That is a sparse estimation and without sparsity inducing prior there will be not enough signal in the noise to estimate `n*(n-1)/2` love or hate relationships. But with a spike and slab prior

 ![sparseinteractionsticky](https://global.discourse-cdn.com/julialang/original/3X/c/9/c924f7650478af3802814dae357a786092909eb0.png)  
Figure: Thin vertical lines indicate distance to the truth. True zeros are plotted with the symbol ×, other are plotted as points will color gradient corresponding to the ground truth.

we estimate the love and hate matrix with something like 2500 unknowns reasonably well (compared to thresholding) .

 ![Screenshot 2021-03-16 at 11.34.22](https://global.discourse-cdn.com/julialang/original/3X/8/6/86b56df4a0061618ba1a1b117042d4562ab775e0.png)

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [March 16, 2021, 11:53am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/5 "2021-03-16T11:53:07Z")

</div>

Nice package 🙂 at the end of my doctorate I wrote a package for PDMP samplers which is now pretty outdated ([PDSampler.jl](https://github.com/alan-turing-institute/PDSampler.jl)) and which I unfortunately don’t have the energy to maintain anymore. It also implemented the ZZ and BP sampler as well as funkier stuff on factor graphs. It’s nice to see a fresher and better package take over! Good luck!

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 16, 2021, 12:05pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/7 "2021-03-16T12:05:37Z")

</div>

Thank you! I see that I starred your package and then must have forgotten it, sorry for not keeping you in the loop. Do you have something written on the fancy factors? That’s very interesting for us, we thought quite a bit about PDMP on factor graphs. We should also check if there is something to port from your package? @SebaGraz

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [March 16, 2021, 2:18pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/8 "2021-03-16T14:18:34Z")

</div>

No worries I didn’t expect you to (especially given the package is not actively maintained anymore). Re factor graph, you could check [this part of the docs](https://alan-turing-institute.github.io/PDSampler.jl/latest/aboutpdsampler.html#Local-Samplers-1), some references are [[1510.02451] The Bouncy Particle Sampler: A Non-Reversible Rejection-Free Markov Chain Monte Carlo Method](https://arxiv.org/abs/1510.02451) and maybe [[1701.04244] Piecewise Deterministic Markov Processes for Scalable Monte Carlo on Restricted Domains](https://arxiv.org/abs/1701.04244) 🙂

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 16, 2021, 3:10pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/9 "2021-03-16T15:10:17Z")

</div>

We haven’t done yet the BPS on a factor graph, so that goes onto the todo list.

What we can add easily is sampling processes constrained to be positive! Because we already compute the time it takes to hit the coordinate axes and put it in the priority queue, we can also just do a boundary reflection from your [https://arxiv.org/pdf/1701.04244.pdf](https://arxiv.org/pdf/1701.04244.pdf) instead of freezing to the axis as in the sticky version. That’ll be useful for sampling parameters supported on [0, ∞) without reparametrising.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [March 16, 2021, 3:33pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/10 "2021-03-16T15:33:49Z")

</div>

> [@mschauer](#):
>
> What we can add easily is sampling processes constrained to be positive!.. we can also just do a boundary reflection…

This is really interesting. Any thoughts on how far this boundary reflection can be taken? Could this be used to sample arbitrary convex domains?

---

<div class="post-metadata">

**Author:** ![SebaGraz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebagraz/32/12414_2.png) [@SebaGraz](https://discourse.julialang.org/u/SebaGraz)\
**Post date:** [March 16, 2021, 4:58pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/11 "2021-03-16T16:58:03Z")

</div>

The only restriction mentioned on [https://arxiv.org/pdf/1701.04244.pdf](https://arxiv.org/pdf/1701.04244.pdf) is to have an open, pathwise connected domain with Lipschitz boundaries. From a computational point of view, you need to be able to evaluate the boundary and its angle to compute hitting times and reflections at the boundary. One last thing, the process might have some funny recurrent behaviour (I am thinking of the Zig-Zag targeting a uniform density on a square), but they can be easily tackled by adding refreshments.

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [March 16, 2021, 5:51pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/12 "2021-03-16T17:51:34Z")

</div>

Yeah so @SebaGraz nailed it. IIRC from coding it, you can indeed “in theory” have this work for any domain that meets these conditions **but** there’s a bit of a hidden trick in that you have to be able to compute the normal at the incidence point to compute the reflection quickly and efficiently (if you do this with a numerical approximation then you potentially lose unbiasedness properties). @SebaGraz said exactly that, just stressing that it’s important.  
For simple surfaces (basically polygons with few faces) this is easy, for other boundaries it can be pretty hard.

Polygonal boundaries are useful though e.g. for GLMs with things like positivity constraints as mentioned above.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [April 1, 2021, 11:29am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/13 "2021-04-01T11:29:13Z")

</div>

PS:

 ![zigzagtrace2](https://global.discourse-cdn.com/julialang/original/3X/c/a/ca20e821329ced31ac8d123ae492e4e82cc46d8c.png)

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [July 31, 2021, 6:32am UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/14 "2021-07-31T06:32:33Z")

</div>

Fresh from the Youtube press:

[ZigZagBoomerang.jl - parallel inference and variable selection | Moritz Schauer | JuliaCon2021](https://www.youtube.com/embed/wJAjP_I1BnQ)

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [October 13, 2021, 12:11pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/15 "2021-10-13T12:11:02Z")

</div>

Is there code for the boundary reflection somewhere? I’d really like to constrain some parameters to be positive. Thanks

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 13, 2021, 2:43pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/16 "2021-10-13T14:43:16Z")

</div>

I’ve got a sampler I call the White Box sampler. It goes in straight lines until it hits the boundary where the log probability density is below some threshold, then it tries random directions on the hyper sphere until it succeeds in stepping back into the high probability region. It samples essentially uniform from the high probability density region. By using this as an umbrella sampler with a simple metropolis Hastings diffusive sampler you can get quite effective sampling from a true posterior and zero need for gradient calcs.

However I haven’t figured out how to plug this concept into the ecosystem, nor have I written a paper about it. If someone is interested I’d be happy to collaborate. I think the scheme has a lot going for it, especially the gradient free nature.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [October 13, 2021, 2:48pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/17 "2021-10-13T14:48:47Z")

</div>

Sounds a little like Hit-and-Run with a slice sampler. Is that right?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 13, 2021, 2:58pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/18 "2021-10-13T14:58:48Z")

</div>

Haha, I’d have to read about hit and run to be able to answer, but… Sounds plausible. It was motivated by the fact that in higher dimensions the log probability density is near constant, and all the volume is near the surface, so sampling uniformly inside the LP threshold is already nearly exact

Also “white box sampler” is motivated by the idea of a photon bouncing inside a perfectly reflective white painted box. At each reflection it goes towards the inside but in a random direction.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [October 13, 2021, 4:31pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/19 "2021-10-13T16:31:42Z")

</div>

Hi @EvoArt , do you have only the positivity constraint? That’s easy to achieve, we’ll prepare something. @SebaGraz

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [October 13, 2021, 5:00pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/20 "2021-10-13T17:00:46Z")

</div>

Thanks! Yes, just positivity constraint for a subset of parameters.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [October 14, 2021, 12:02pm UTC](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287/22 "2021-10-14T12:02:06Z")

</div>

> [@EvoArt](#):
>
> Thanks! Yes, just positivity constraint for a subset of parameters.

Almost there, do you have a share-able example?

[Next page](https://discourse.julialang.org/t/ann-zigzagboomerang-jl/57287.md?page=2)