Merge pull request #119 from mschauer/mass

test old and new mass matrix

Author: mschauer

src/not_fact_samplers.jl

@@ -158,14 +158,27 @@ function ZigZagBoomerang.reflect!(∇ϕx, t, x, v, F::BouncyParticle{<:Any, <:An
     v .-= (2*dot(∇ϕx, v)/dot(∇ϕx, z)) * z
     v
 end
+function reflect!(∇ϕx, t, x, θ, F::BouncyParticle)
+    θ .-= (2*dot(∇ϕx, θ)/normsq(F.L\∇ϕx))*(F.L'\(F.L\∇ϕx))
+    θ
+end
 function ZigZagBoomerang.refresh!(rng, t, x, v, F::BouncyParticle{<:Any, <:Any, <:Any, <:AbstractPDMat})
-    ρ̄ = sqrt(1-F.ρ^2)
+    ρ̄ = sqrt(1 - F.ρ^2)
     v .*= F.ρ
     s = local_speed(t, x, v, F)
     u = (s*ρ̄)*PDMats.unwhiten(F.U, randn(rng, length(v)))
     v .+= u
     record_rate(v, F)
 end
+function ZigZagBoomerang.refresh!(rng, t, x, v, F::BouncyParticle)
+    ρ̄ = sqrt(1 - F.ρ^2)
+    v .*= F.ρ
+    s = local_speed(t, x, v, F)
+    u = (s*ρ̄)*(F.L'\randn(rng, length(v)))
+    v .+= u
+    record_rate(v, F)
+end
+
 function mass_adapt_init(M::InvChol)
     Cholesky(M.R)
 end

test/maintest.jl

@@ -176,7 +176,7 @@ end
     t0 = 0.0
     θ0 = randn(d)
     x0 = randn(d)
-    M = UpperTriangular(I + 0.4randn(d,d))
+    M = UpperTriangular(I + 0.4randn(d,d)')
     c = 20.0
     B = BouncyParticle(missing, missing, # ignored
         1.0, # momentum refreshment rate 
@@ -206,6 +206,41 @@ end
     @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) < 3/sqrt(length(ts))
 end
 
+@testset "Bouncy Particle Sampler (arbitrary mass matrix 2)" begin
+    Random.seed!(2)
+    t0 = 0.0
+    θ0 = randn(d)
+    x0 = randn(d)
+    M = LowerTriangular(I + 0.4randn(d,d))
+    c = 20.0
+    B = BouncyParticle(missing, missing, # ignored
+        1.0, # momentum refreshment rate 
+        0.9, # momentum correlation / only gradually change momentum in refreshment/momentum update
+        missing, # metric
+        M
+    ) 
+
+    ∇ϕ!(y, t, x, args...) = mul!(y, Γ, x)
+    dϕ(t, x, v, args...) =  dot(v, Γ, x), dot(v, Γ, v)
+    n = 800
+    trace, _, acc, more = @time pdmp(
+        dϕ, # return first two directional derivatives of negative target log-likelihood in direction v
+        ∇ϕ!, # return gradient of negative target log-likelihood
+        t0, x0, θ0, # initial state and duration
+        n, # number of samples
+        ZigZagBoomerang.LocalBound(c), # use Hessian information 
+        B; # sampler
+        adapt=false, # adapt bound c
+        progress=true, # show progress bar
+    )
+    @show more
+    @show acc[1]/acc[2]
+    ts, xs = sep(trace)
+    @show length(ts)
+    @test mean(abs.(mean(xs))) < 3/sqrt(length(ts))
+    @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) < 3/sqrt(length(ts))
+end
+
 @testset "Bouncy Particle Sampler (adapted mass matrix)" begin
     Random.seed!(2)
     t0 = 0.0