Listen on your favorite platform:

Adrian Seyboldt is back for a third appearance, and this time he's brought Eliot Carlson, a research analyst in Bob Carpenter's group at the Flatiron Institute.

Together, they spent the last year working on a very hard problem: how do you precondition Hamiltonian Monte Carlo well, before you've even sampled the thing you're trying to precondition for?

That chicken-and-egg problem, and the paper Adrian and Eliot wrote to attack it, is the spine of the conversation. But the episode also doubles as a state-of-the-report on Nutpie, Adrian's Rust-based NUTS sampler, three years and two prior LBS episodes into its life: what's changed, what's still experimental, and why "just switch your PyMC default sampler" has become genuinely good advice rather than a plug.

Mass Matrix Adaptation Is Just Automatic Reparameterization

But first, what even is max matrix adaptation?? Forget the inner-product interpretation; think of mass matrix adaptation as an automatic, dumb, but very effective reparameterization of your model.

The simplest version, the diagonal mass matrix, just rescales each parameter so its posterior standard deviation becomes one -- exactly the kind of manual rescaling you could do yourself by hand, badly, if you had infinite patience.

Here's where things diverge. Stan and vanilla PyMC only look at one source of information: the posterior standard deviation estimated from warm-up draws. Nutpie uses a second, underused signal: the scores, meaning the gradients of the log density, which HMC is already computing at every leapfrog step to build its trajectory.

Nutpie's diagonal adaptation takes the geometric mean of the two resulting standard deviations rather than relying on draws alone, and while there's no guarantee it always wins, in practice it usually does.

Four Preconditioners, One Practical Default

In practice, you can choose from four versions of mass matrix adaptation in Nutpie:

  • Diagonal: the fast default, essentially free overhead, and the right choice whenever your posterior has no meaningful correlation structure.
  • Dense: mostly a mathematical stepping stone used to derive the third case; a full mass matrix scales too badly with dimension to be practical.
  • Low-rank-plus-diagonal: the one Adrian actually recommends reaching for. It rescales with the diagonal first, then corrects the most significant correlations using a low-rank structure, at a cost of a few extra seconds of tuning time.
  • Normalizing flow (experimental): the option for pathological geometries that no linear mass matrix can fix, like the classic hierarchical-model funnel, where any linear rescaling that tames one side of the funnel makes the other side worse. It's a tool worth having, but rarely the one you'll reach for first.

The paper's headline result is that low-rank-plus-diagonal adaptation beats the standard draws-only approach by a median factor of four across 114 models, drawn primarily from posteriordb.

Adrian and Eliot aren't surprised by the direction of the result: static correlations are everywhere in real models, and low-rank correction targets exactly that.

Learning the Parameterization Itself

A very exciting line of work after this one is to automatically choose between centered and non-centered parameterizations (and every gradation in between) during mass matrix adaptation itself.

Centered and non-centered forms of a normal distribution are mathematically identical but can produce wildly different sampler performance depending on how much data informs a given parameter, and picking correctly by hand is tedious enough that almost nobody does it consistently.

Framed as a single, extremely cheap interpolation parameter per normal distribution, the transformation can be learned automatically during warmup, and on a fine-grained, per-subvariable basis rather than as one binary switch for an entire model.

As someone who builds models on a daily basis, this is an issue I've been facing for years, so I'm genuinely excited about this work!

That's it for this one folks. Remember you can also interact with the episode on NotebookLM -- ask questions, generate flashcards, and more.

Hope you enjoyed it, and see you in two weeks, my dear Bayesians!

Chapters

00:00:00 What is HMC preconditioning?

00:09:03 A more robust low-rank mass matrix

00:11:58 What is mass matrix adaptation?

00:18:06 What does preconditioning HMC mean?

00:20:57 What is normalizing flow adaptation, and when does a linear mass matrix fall short?

00:23:50 When does normalizing flow adaptation actually help, and when is classic mass matrix adaptation enough?

00:27:13 What is Fisher divergence?

00:30:10 Why is HMC's trajectory, not its density, the right target for preconditioning?

00:33:04 What are the diagonal, dense, and low-rank-plus-diagonal versions of mass matrix adaptation?

00:46:25 How much faster is low-rank-plus-diagonal adaptation?

00:51:07 What's the practical recommendation for using Nutpie and its mass matrix adaptation?

00:54:31 Why does low-rank adaptation sometimes fail spectacularly?

01:01:35 Where does this research fit in the bigger picture of HMC?

01:12:12 How could centered vs. non-centered parameterization be chosen automatically?

My guests today are Adrian Seyboldt and Eliot Carlson.

I'm sure you know Adrian by now, but this time he's brought Eliot with him, a research analyst in Bob Carpenter's group at the Flatiron Institute, who spent the last year

co-writing a paper with Adrian.

That paper is about preconditioning, the process of transforming your target distribution into something Hamiltonian Monte Carlo can temple.

Efficiently.

The catch is a chicken and egg problem since preconditioning requires knowing things about your posterior that sampling is supposed to reveal in the first place.

Adrian and Elliott's fix is to use a piece of information samplers already compute, but mostly throw away the gradients of the log density and to formalize the whole problem

using something called Fisher Divergence.

We get into mass matrix adaptation, why Stan and Pime C's default approach leaves information on the table, when to reach for diagonal versus low rank versus normalizing

flow adaptations, and a benchmark results that's hard to ignore.

A median, four times speed-up across 114 models.

This is Learning Bayesian Statistics, episode 163, recorded June 12, 2026.

Show you how to be a good basian change your predictions after taking information.

And if you think it now be less than amazing, let's adjust those expectations.

What's a Bayesian is someone who cares about evidence.

Welcome to Learning Bayesian Statistics, a podcast about Bayesian inference, the methods, the projects, and the people who make it possible.

I'm your host, Alex Andorra.

You can follow me on Twitter at Alex underscore Andorra.

like the country for any info about the show.

LearnBayesStats.com is Lab Plus2B.

Show notes, becoming a corporate sponsor, unlocking Bayesian merch, supporting the show on Patreon.

Everything is in there.

That's LearnBayesStats.com.

If you're interested in one-on-one mentorship, online courses, or statistical consulting, feel free to reach out and book a call at topmate.io slash Alex underscore and Dora.

See you around, folks, and best Bayesian wishes to you all.

Eliot Carlson, welcome to Learning Patient Statistics.

And thanks for having me.

Yeah, yeah, uh, I'm very happy to have you on.

And obviously, Alrien Zabold.

Willkommen now.

Learning patient statistics.

Yeah, always great to have you back, Evan.

Um you're basically a co-host now, so you know you are you're at home.

feel free to interrupt me whenever you need, and and and say whatever you want, you know.

but today we're gonna talk about what you guys have been uh working on because you've been working together.

Uh I think also with uh Bob Carpenter on that project.

and so we're gonna talk about a lot of uh

algorithm nuts, nut by uh rest, th things like that today.

it's gonna be a a fun episode.

But to start with, of course, Eliot, it's your first time on the show, so we need your origin story.

Okay.

What what are you doing nowadays and how did you end up doing that?

Okay, well I'm currently a uh research analyst in Bob Carpenter's group at the Flatiron and also been the last year working

with PyMC Labs, to whom I was introduced via Adrian.

Um so I did a undergrad in in stats and a master's in stats and sort of I think I I arrived at Bayesian stats the reverse way that most people do, where um or excuse me, to

DMC in the reverse way.

So I first learned about MCMC without really knowing much about uh Bayesian stats and was interested in in that.

And um

So I was introduced to Bob and asked if I could uh if he had a project, an MCMC related project for me to work on.

And then he uh introduced me to Adrian.

He said, I have a an algorithm I've been trying to understand for a year now and but I haven't had the time to fully dive into it.

So then that that initiated maybe another a year long process of us writing this paper, which was a lot of fun and and yeah, learned a ton.

This Adrian guy just

never wrote down his algorithm in a proper way, right?

I mean, I understand you, Adrian.

Like writing papers is is not the funniest part, I would say, of of the job.

Uh it's not the one I enjoy most.

So you know, um I think it's much, much more fun when you do it with other uh interesting.

Actually writing it with with Bob and Ali it helped so much.

That was a lot of fun.

So Yeah, yeah, exactly.

Writing alone is just my god.

Yeah, I hate that too.

And so Ali actually so you you answered that a bit, but I'm really curious also to what drew you to patient inference specifically, um, because I think here you you told us how

you end up ended up working on HMC basically and on samplers.

Um is there w were you introduced to patient stats before that or did you discover them also just with this uh with this sample with this sampling work?

Yeah, so I think the chronology was I was first taking a class in stochastic processes and Markov chains and then I learned about Metropolis Hastings and I thought, wow, this is I

don't know why, I just thought, wow, this is really cool.

Uh even though I I didn't really know what it was used for at all.

Um, which is kind of what I meant in like reverse way, 'cause I think most people view sampling as like this ugly thing they have to do when they want their models to fit.

But anyway, so once I

you know, in the past couple of years, now that I've been exposed to more uh you know, of the Bayesian methodology, I I've become more interested in that in and of itself, I think.

Um, especially just seeing how powerful it is for especially like researchers in in, you know, scientific domains who have prior information from other exper experiments and

things like, you know, Adrian and I have worked on a project where this has come up a lot.

Um, where yeah, people have priors or, you know, from other experiments that they need to integrate into their models.

Um, and then but I I would say that Adrian may be the same, though the fun comes in the sampling part.

Yeah, I think I think Adri Adrian will relate to that.

And well actually drainings to third time on the show, so I'll I will refer to

To your previous appearances, uh episode episode seventy-four, I think, was your first one where we talked about your well, that was, you know, I think the first time.

Yeah.

and and that was the first time we talked about Zero Subnormal and also about Nutpie and and the work you did on that and why you started doing that and why Rust is interesting.

So I encourage people to listen to that.

Um episode seventy six also with

Bob Carpenter, uh which is a great background.

In episode one thirty three, where you joined as a guest star with Sean Pickney, uh because we talked about all things um algorithm efficiency, but also Zero Subnormal in

Stan, which Sean implemented, so that was a nice crossover episode.

so I recommend listeners to to to listen to these ones, but um maybe can you tell us

You know, um since then, what what changed or what evolved in the project, especially NutPi and also what you've been working on lately beyond the paper, of course.

Yeah.

So kind of in NutPie itself, kind of from a high level I would say things haven't changed that much.

So the basic mass matrix adaptation algorithm has stayed the same.

There was a little bit of tweaking, like the the window sizes are different and they are growing now, things like that.

Kind of around that though the project has changed quite a bit, I would say.

So got a lot more robust, so has a lot more users, so they definitely notice whenever some whenever whenever I break something, people actually report report it in issues, which is

uh a strange feeling.

So you you notice that people seem to be using it now.

so yeah, just has been a lot more tested, I would say.

I guess one bigger muse

item is that it has an R interface now as well.

So not just a Python interface, but an R interface.

That's actually not written by me.

Um somebody else wrote that and was really happy that that appeared.

So uh people can now use it from R as well if you have stand models or something like that.

And I think that's probably the the biggest change there.

Also on the mass matrix adaptation side, the load rank version, which has been around for a long, long time, it

It works now much better.

So when I tested that years ago before writing the paper, it was like somehow it didn't really work that well.

Then when we wrote the paper, I just ran the benchmarks again and suddenly it was working well.

And never was able to figure out what what really changed.

Maybe my benchmarks were just wrong or I I no idea.

I had the idea it should work well, but it somehow didn't for some reason.

Now it does.

So I'm not sure if that's actually a change in NutPi or just uh

me looking again.

Um but yeah, that's uh definitely worth a look if uh yeah to for anybody who's working on on models that are not trivial to sample.

Yeah, yeah, actually we can we can talk a bit more about that and I think it's a nice segue into your paper.

Um so can you can you tell people what even you know mass matrix adaptation means?

I guess it's a

It's a complicated word, an intimidating word when you haven't heard about that.

So what even is that about?

Why do you why would you even care about it for a for a model?

So they're d dif completely different ways of thinking about the mass matrix.

One way is it defines an inner product and you need an inner product for nuts.

That's kind of the math side.

If you're a mathematician, that maybe tells you something.

Otherwise, maybe that's not too helpful.

Kind of

If you have divergences or something, the sampler always tells you, yeah, you you need to reparameterize.

So you need to find new parameters, new way to write your original model.

And the way I tend to think about the mass matrix adaptation most is as an automatic reparameterization in a sense.

Just a really, really dumb one that doesn't fix much, but it fixes really, really important things.

Namely kind of the diagonal version of mass matrix adaptation.

just rescales your variables.

I would also add that like, yeah, I I think the sometimes the term mass matrix, it sort of drags you into the physical interpretations of HMC, which are really beautiful if you

understand physics, which I do not.

So but but really like you can just think of it as um what we call like preconditioning.

So just changing the geometry of your distribution.

And you can show that

Like using a specific mass matrix is the same thing as just sampling from this uh sampling from a slightly, you know, shrunk or stretched distribution or or a distribution where

that's some linear transformation has been done that transformation B in the mass matrix.

Yeah, thanks, Eliot.

I th I think that's a very intuitive way for for people to understand it.

And actually, so each time each time people use Pine C or Stan or other PPL, there are

doing that kind of things under the hood, even though they may not know.

And and that's the that's exactly the beauty of PPS, right?

It's that you don't have to care about that.

eh And so yeah, HMC has that large, if I understood correctly and you guys reframe me if I'm wrong because you're the expert, but it's it's a ha it has this like large constant

factor that's determined by the curvature of the target density.

And so what

mass matrix adaptation does, as you said earlier, these preconditioning on that to make sure we're able to sample from the posterior afterwards after we've tuned.

but yeah, how does that work currently um in Stan and PyMC, for instance?

And what does Nutpie do differently so that we understand a bit um a bit more the difference and the

The value proposition if you want of NetBuy, and then we'll dive into your newest paper.

So the let's stick with the diagonal mass matrix for a second.

Um so that you can think of it looks at each parameter and tries to just rescale it.

So if it's maybe has a posterior standard deviation of ten to the minus four, then it will just scale it up by ten to the four to make it bigger.

So that kind of

the the this posterior standard deviation of each parameter as the sampler sees it is identical.

And that's exactly what Stan and kind of the default PyMC sampler, what they are doing.

They're just saying I look at the posterior standard deviation and I rescale it such that the posterior standard deviation is one.

And I mean you could do that manually if you really, really wanted, just introduce new parameters that are rescaled versions of the original ones.

You just wouldn't want to.

That's not a good use of your time.

So we do it in the sampler.

m So standard PyMC, they just use um the posterior standard deviation directly.

We look at the second source of information, we look at the at the scores, so the the d uh the derivatives of the lock log density.

Kind of in an ideal world, you want that the standard deviation of your parameter is one, but also that the standard deviation of your

scores so of those derivatives that those are also one.

They are both kind of interacting in in in in some sense.

and NutPi just uses the geometric mean of those two standard deviations.

And then there's a long discussions kind of how you can can kind of say why that is a good idea.

That's what we go to in the paper a bit more.

And kind of a much more general framework of kind of how to

Do good preconditioning and kind of this this particular version just pops out of that as a special case.

But in this special case, it's just the geometric mean of those two two things that we kind of want to want to rescale in a in a sensible way.

And usually, so that's there's no guarantee that you can, if you if you work on it, you can find examples where it's worse.

Um, but usually it's better.

And um

Usually is I think where samplers are concerned, the the best you can pretty much always that's the best what what what you can do.

So yeah.

Aliot first, uh maybe you have you have something to to add or um yeah, anything you're thinking about on that?

Um I would say maybe this is uh jumping ahead, but uh you know, one um basic basically w you know a lot of what the paper is doing is formalizing

some intuitions and general things that we want out of the sampler.

So like one uh one thing I would add is that we we basically know that a standard normal is uh very perform HMC will perform really well on a standard normal.

And you have this nice fact about standard normals relating the covariance to the um

covariance of the draws of the distribution to the covariance of its scores.

So if you think about so if you have like a you know a normal distribution is parameterized by a mean and covariance.

And so if you look at the if you look at the distribution of the scores, the gradients of the log density, their covariance is the inverse uh sigma, the inverse covariance of your

normal, which is kind of this cool

fact.

And so that provides a little bit of intuition as to why.

So so then when we take the uh trick mean, we sort of want to combine those two pieces of information.

Uh so we'd look at the inverse scores and uh the draws.

Okay.

Yeah.

Yeah.

Um actually I think that's a good that's a good moment to um to dive into your your paper, which is in the show notes for people who want to dig deeper because here of course will

We won't make it justice, uh because it's it's a very um thorough and and deep paper.

Um so titled is preconditioning Hamiltonian Monte Carlo, I think all of that we just talked about, by minimizing fissure divergence.

So we'll go into that second part just afterwards.

But first, can you maybe give us the elevator pitch for the paper and and tell us what

what it does intuitively.

So the big problem that Stan faces or any s HMC sampler is that when we're you know preconditioning, basically we want to find a some mass matrix that uh turns our

distribution into something that's friendly for HMC to run on.

And to do that, you need to know something about your distribution.

You want to know maybe it's covariance structure or something.

But the whole point of

Doing sampling is that we want to learn about our distribution.

So, you know, things like its covariance structure.

So that we call that this this codependence, um, the chicken and egg problem.

It's known as.

So, so basically, you know, ordinarily like preconditioning based on draws is how it we sort of initially come up with like a very rough sketch of what the distribution might be

based on just a few warm-up draws and then and then sort of get the motor going.

Um and so that

can often take like a huge amount of the sampling time.

And so basically the so as we said, the, you know, we have this other piece of information in the form of uh scores that samplers like Stan are not using.

And uh that can be like really valuable for preconditioning.

And since we have them, we might as well try to use them.

And so a lot of the work of the paper is formalizing um

you know, what information can these scores give us about uh the the distributions and s and formalizing like how how do we turn this into a how do we best make use of this

information, I guess.

Yeah.

Um yeah, I think that's a great explanation.

Makes it a bit more intuitive for people.

Uh Adrian, do you wanna add anything to that?

Yeah, so I guess what I would add is the idea of the general framework of how to do preconditioning.

So kind of what the paper really does spend most of its time on is a special case of this.

So we say we have this objective, how we want to uh precondition um the model.

And mostly the appendix kind of have such in a very general way written down.

Um and then most of the work is what does that do in a very special case where we just look at the mass matrix that does linear things.

And there we work through kind of what has that what what are the in influences on the condition number, why does this work better?

Kind of but but as I said, just in the linear case.

And what I really like about this framework is that it's not restricted to that.

You can add nonlinear transformations as well.

Um that's kind of where the normalizing flow adaptation that's still experimental in in Nutpie, for instance, is coming in.

Um and I think that's where a lot of the interesting research, more researchy side of things look like, kind of the the

Linear one, I would more say is more in the engineering part where okay, we have the objective, we derive it, we look at look at derive closed form solutions, we look at the

the details of that, the closed form solutions and inspect what it does in that particular simple case.

And on the research end, we kind of have this how can we generalize this?

And what's that what does this tell us in general about how to precondition um HMC?

Actually maybe we can talk a bit about the normalizing flow adaptation because I think it's a pretty intuitive concept and that will that will set us up well to then talk a bit

more about feature divergence and and mass matrix adaptation.

So yeah, can you or Eliot um give us the the idea of why normalizing flow adaptation would be interesting here?

and what it does.

And so you don't need to redefine what normalizing flows are.

I have an episode on that with Marilou Gabrier.

I'm gonna put it in the related episodes.

But yeah, feel free to just like give the elevator pitch for it and then make it in the context like give us in the context of Nutpie and Mass Matrix adaptation what it does and

why it can be interesting.

I think that especially in applied problems, there are a lot of

uh geometries that are induced by certain parameters that create uh that that basically are not solvable with a linear preconditioner at all.

So the most famous example of this is like the funnel, which arises when you have these hierarchical models where there's a variance parameter.

So you get this, you know, big thing that's called a funnel.

And basically uh on either side, any any linear uh preconditioner you try to apply to

sort of um tame one side of it, you know, this very thin over here and very wide over here.

Anything you do to one side will just make the other side worse.

Um so uh normalizing flows since they can sort of uh since it's a a neural network, it the you know obviously very not restricted to be linear.

Um and so you can get kind of sophisticated, more sophisticated preconditioners.

And yeah, as you mentioned, nor normalizing flows, you have the

uh they're designed in such a way that you have the normalizing uh you have the Jacobian, excuse me.

Um and so that that can help you.

Yeah.

So I actually I have not used normalizing flows myself, but I am I would be interested in in trying them out.

Yeah.

Yeah.

Yeah I've tried them I tried the adaptation uh of normalizing flows.

but it was in a

huge GP model and in the end it didn't like it wasn't giving um speed ups for assembly.

Actually the the classic mass matrix adaptation was perfectly fine.

Uh so it really depends on the models and on the geometry uh that you get with the the posterior and if it's especially complicated with well, you know, funnels as you were

saying, then it's gonna help a lot because the what I find super

Intuitive and interesting with the normalizing flow idea is that basically it just um tries with a neural network to learn a transformation of the space that becomes a standard

normal.

and that really ties to what you were talking about at the beginning of the episode where, well, you know, standard normals are super easy to sample in HMC, so why don't we try to

do that?

And

Yeah, the idea is that if you have a very complicated posterior geometry, spending more time in tuning and adapting the mass matrix with the normalizing flow can be very

interesting because then it's gonna pay off a lot afterwards when you start sampling.

Uh obviously the caveat is that, well, you have to learn a neural network, so uh tuning's gonna take much more time.

But basically the idea is a bit like to amortize that.

that tuning time with then really fast and efficient sampling.

Um and that's also why like pairing that adaptation with normalizing flow with GPUs is gonna play extremely well, especially if you are using NutPine to do that.

So yeah, like it's it's interesting to have in your toolbox for sure, even though you're not gonna use that all the time and really not the majority of the time actually.

This kind of I would say like

uh leads into the the Fisher divergence topic.

I don't know if you want to uh talk about that.

Basically Yeah, that's perfect.

Exactly what I was gonna ask you.

Yeah.

Um so just I'll give I'll give like a brief a a simple background, I guess.

So most people know about K L divergence.

It's basically just divergence as being like a way of trying to compare uh I don't I mean I'm gonna use the word metric.

I know they're not

like technically in metrics, I think, but but it's basically a metric of comparing, you know, how close are two probability distributions to one another.

So the KL divergence is the most common one where you you know you're integrating over uh the the support and you know comparing the the log densities.

The Fisher divergence is the slightly different uh metric which where you're comparing the scores and you're basically saying, okay, you know, distribution A and distribution B.

uh when we integrate over those, how close are the scores to one another?

And so the so we can use that as a uh minimization objective for preconditioning.

So how do we get our our uh target distribution, how do we minimize the fissure divergence between our target distribution and a standard normal, which we know works really well.

And so that's

Yeah, that's all that Fisher divergence is sort of a formalization.

And and the nice thing about that is that the the solution to that minimization problem incorporates both draws and scores, which are exactly the two pieces of information that

we have uh during preconditioning.

So it sort of uh works out works out really nicely.

Interesting.

So what is the two two further questions?

what is the difference with KL?

divergence with most of the listeners will know about.

And why is Fisher divergence here a better objective for what we're talking about?

And maybe additional question uh completely random.

Do you know why it's called Fisher Divergence?

Which is uh a like a very difficult, like non intuitive name?

I actually don't.

Adrian, do you know why?

I guess uh because of Fisher?

kind of he was

Doing a lot of work with sc with the score function, so I would assume that's related.

It's also closely related to the Fisher information matrix or Fisher information um in general.

So I would assume it's coming from that, but I never really checked.

It's also really pretty much the same thing as score matching, kind of in different literature.

That's uh it's also called called that.

Um Okay.

Nice.

Well, after that uh deranging of the

of the thought process because of me.

So what is the difference actually with with KL divergence?

so the KL KL divergence really looks at the density and the Fisher divergence looks at the scores, so the derivative of the density.

And it turns out that you can use both kind of if you know the the full if you know know if you have a cop probability distribution and you know what the der derivative looks like

everywhere, then you know the full

the the full the full distribution.

Similarly to how if you know the density, you know the full distribution in a sense.

So that those are really complementary.

You can use either one or the other one.

Um the reason the Fisher divergence in particular or two reasons the Fisher divergence is particularly useful for us.

One is that it just uses the information that we have, which is a nice thing, because we don't really have the density in the first place.

We only have the density up to a constant, but we do have the derivative, the score function directly, because the constant goes away when we take the derivative.

So that's nice.

the other reason is more related to HMC directly.

So HMC really uses the score function to drive you kind of to to find the trajectory through the parameter space.

So that's really the thing that we want to target in a sense.

So we want the trajectory of um our samples to look like the trajectory of a standard normal distribution in a sense.

And that's because that's what makes it efficient, the the trajectory.

It's not so much the density in itself.

Um so for HMC specifically, I think it's the right thing to target.

You could actually add additional terms, but that's a different question.

So you could look at higher order derivatives as well.

I think that's re

interesting research uh direction as well.

Yeah, super thanks, Adrian.

Super clear answer.

So and that ties well with um ties up well with the you know physics interpretation that you were talking about before all Yacht because like yeah okay the the derivat derivatives

like the second order and the like kind of the well if I remember correctly if it were a car that'd be the the speed right or that'd be the acceleration like

I I always mix mix them up.

but yeah, like um okay, so it's like you're interested in the second order the change in the derivative and not necessarily and not the in the density, sorry, and not the density

itself.

And that's what the Fisher divergence tell you.

Um divergence from what though?

Like it's so it's like you look at the density you're trying to sample, uh the the direct derivative of the density you're trying to sample.

And you're looking you're comparing E to something to say that it's diverging?

So you a divergence is just the math term to kind of like a distance like thing for different probability distributions.

It's not really a distance because it's not symmetric, for instance, but it's like a math term for distance for probability distributions.

I'm actually also not entirely sure where that nomenclature, where that name is coming from.

But um Yeah.

So it's the yeah, how different is our standard normal to our transformed posterior?

That's really what the divergence isn't measuring there.

Yeah.

Yeah, and also if I remember correctly, it's like we're talking about divergence and not distance because it can be asymmetric.

and it it doesn't have to be symmetric.

So so that's one of the of the main thing.

I think earlier to you wanted also to to answer my my question.

Yes.

W wait, which which question?

Sorry, I've I forgot.

Divergence to what?

And I think what you were saying is that well we are looking at the difference between like the the distance again in in quotes, so but the divergence between the target

density, well the derivative of the target density and I'm guessing the derivative of a standard domain, I think was your answer.

Yeah.

Yeah.

So I think like underneath the hood uh maybe it's a little in the weeds, but um like at any given point

you sort of want to you know you have the s of your target you have the point and the score at that uh at that point and then you'd want to say, okay, well if this were a

standard normal, what would what what would the score be at this point?

And then we want to change the distributions that those scores match.

Yeah.

Okay.

Yeah.

Yeah.

Yeah.

Um that's that's very clear.

And now actually I'm

So I'm wondering you in the paper you present three version of these of these adaptations and we've also mentioned a few of them here already so it's like diagonal, dense and low

rank plus diagonal uh mass matrix adaptation.

So can you walk us through those and when you would reach for each one?

So the diagonal one is just the default.

That's kind of the simple one that's really, really fast.

So it doesn't really add any overhead.

So if you've got a really nice posterior, no correlations or something, the diagonal is the one because kind of that does the work, doesn't spend much time doing anything else.

Um the dense one is really just a mass math thing that we needed in the paper.

So that's not even implemented.

You could implement it.

There's nothing particularly wrong with it, but I think the low rank plus dense is uh low rank plus diagonal is really kind of the practical thing that that that you want if you

want if you want to have uh a way to fix correlations in your posterior.

So the dense one is just we need that to derive the equation for the third one.

and the third one, the diagonal plus low rank, really is.

Rescale it first with the diagonal mass matrix and then um try to fix some correlations that you can find and in a way that it scales reasonably well with the dimension.

So kind of if you were to look at the full mass matrix, it would get really, really slow, because that's just computationally a lot to do.

Uh in as unless you have a very, very low-dimensional problem.

Um

So we use a low rank version of that.

It still adds some computational overhead, but it's not not crazy.

So it scales reasonably well.

Um it does usually add a couple of seconds or so, maybe to the total sampling time, maybe a bit more.

So if you're really fast that you don't want that.

If if you have a larger problem, it doesn't matter at all.

It's like the that that's a rounding arrow.

Um

And I've actually thought about making the low rank plus uh the low rank version the the default at some point.

That's not the case yet though.

But you can think of it, I think, by now as a reasonable thing to do and and not just this crazy now I need need the the big guns thing.

Like that's that that's the flow the the the flow version.

That's kind of the crazy one.

But the low rank version I think that's mostly em

Mostly quite workable by now.

You might end up in kind of certain situations where something's unstable or so it might not work even though the diagonal one works if something's ill conditioned somewhere.

So that can happen, but I don't think it's particularly common.

Okay.

And and so yeah, actually maybe the normalizing flow adaptation would be the fourth one, right?

The the last one you reach for once you've reached for the the three previous ones and and they are not working for your case.

Yeah.

So the the flow version I'm mostly a bit skeptical of kind of telling people to use it because the implementation is a bit messy right now.

So you might get weird error messages.

It's slower than it could be.

Um so there's a lot kind of to to do on the implementation there.

Math wise, I'm relatively happy with it so far.

Um but just the engineering takes uh I think takes a bit more time.

Also kind of kind of saying a normalizing flow is doing a lot of work there, kind of

what normalizing flow.

There's like a bazillion different ways of doing this precisely.

And which one we kind of would want to be the default is like I don't know how to choose that at this point.

so that is a big question.

Like what particular normalizing flow?

And we have one that's going to change for sure.

So if you have a model that works now, I can no not guarantee that it will also work in a future version if we happen to change the the normalizing flow for instance.

Then that will definitely break some things that worked before and probably hopefully fix m more things that yeah.

It's like that if you change those uh those settings.

Um so that that I would think of an experimental thing that I'm happy if people try it, but do expect a couple of round uh rough uh yeah, roughly.

Rough edges, yeah, for sure, for sure.

Um and actually

So I've I have a question related to that.

But first, Eliot, do you wanna add anything about this uh presentation of the different method of mass matrix adaptations?

Yeah, well I w I was gonna add uh anecdotally, like I I mean, I pretty recently had a GP model that w took like an hour to sample with diagonal and then changing it to low rank,

it uh was like a matter of seconds.

Like the speed up was was crazy.

Um

So I would definitely like you know, if it it's definitely worth trying out.

I especially for for I think models where you have some reason to think that you might have like a very correlated posterior for some reason, either due to like constraints or I

I I didn't go under the hood and really figure out why it was making that GP model so much quicker, but I think um it could have something to do with that it was a

a GP2, maybe it works particularly well for those, I'm not sure.

Um but yeah, so I I would encourage people to try that out if they're uh if uh if it's taking a long time with the diagonal it could it could speed give you an insane speed up.

Yeah.

Yeah, so you mean uh going from the default diagonal adaptation to the uh low rank plus diagonal run, right?

Interesting.

Uh was it was it a vanilla GP or an GP approximation like S HS GP or something else?

Actually don't remember.

It might have been vanilla.

Um Yeah, okay, yeah.

Yeah, that's my intuition too.

Like in my experience, it's gonna it's gonna help with uh vanilla GP, but if you're using an approximation already like HS GP, which was my case actually for the model I mentioned

earlier, it's not.

It's not helping much because the adaptation, like the diagonal adaptation is already working well because the model like the geometry is much easier because you thanks to HSGP

you have it's much more of a almost linear problem.

Um and so that makes for a a much easier geometry.

Which leads to actually my follow-up question on that, um maybe for Adrian first and then Eliot, it's okay, so we have these different ways of

Doing the adaptations.

But the issue is that, and I guess it's gonna be for most users, you know, like not for you, because you have an intuition, and once you start working on these models, you get an

intuition of what's what's gonna work best.

But my question is: well, if I only know what the best adaptation is once I've sampled my posterior model, it's not that useful, right?

Um and that makes the the iteration much longer than it needs to.

So

Is there a way and even is the research and something you guys maybe are working on, which would be like ideally for the users, I'm thinking, well, we do that at adaptation

ourselves, you know, under the hood, and we're able to more or less understand what the posterior geometry is gonna look like and then we just under the hood, not points like,

oh, here I'm gonna use low rank plus diagonal because it's gonna work better without even the user no knowing it, instead of

doing it manually and incrementally.

Um so yeah, like is that something possible and what would you tell people right now to to do in the So I would say i it it's a good question for sure.

Kind of the this idea of users just trying lots and lots of different sampler configurations is not something I'm actually too happy about.

So kind of this stand way of here is the sampler, write the model, if it doesn't work, fix your model.

I think that's in many ways that's the right thing to right way to think about if you're use building statistical models and working with them.

If you start playing with the sampler and spending lots and lots of time of playing with the sampler, at some point it's it's just but much better to fix the model.

Just just fix the model and the sampler will not be the problem anymore.

As a developer, that's of course different because I want to kind of increase the size of things that that work out of the box, right?

So there's a bit of a same.

Yeah.

Um so automatically choosing.

I think that's what we want.

And in a sense, that's what the low rank version already does.

So it will kind of compute the full or not the full, it will compute a mass matrix with a larger number of eigenvalues and then throw away those that look like they are too small

to do anything useful.

So it only keeps those eigenvalues that look like they do useful work.

And it gets rid of the other ones.

Um Okay.

Yeah.

So there's still the additional cost of figuring out what were would would have been the options, right?

There's also the cost of okay, that there's this additional cost, so we don't want to do it every step, for instance, kind of the the diagonal mass matrix we update in every draw

the

Low rank one we don't update in every draw because that's takes time.

So there there is a trade off and I don't don't know a way of getting rid of it entirely.

But as I said, we are trying to kind of make the low rank version kind of the the thing that chooses automatically.

But how well we are there, that's yeah.

Yeah, yeah, yeah.

I I I understand.

But yeah, that'd be

That'd be great.

So and it sounds like at some point you're gonna make the low-rink adaptation the default in NetPy once you're happy with uh what you're working on right now.

That might happen, yeah.

I don't want to kind of set that say that promise that in any way, but I think that might happen.

Mm-hmm.

Okay.

Eliot?

Anything you want to add?

Yeah, I'm I I would add I guess the Yeah, I think for

At least in my experience, I hope I think Adrian would agree, like if for any sampler you can always conceive of a posterior that will, you know, will totally thwart it.

Yeah.

Um, you know, you can always have either like tail behavior that has some totally different covariance structure for some reason or whatever.

There's also just the classic, you know, there could be a mode a million miles away that you'll never find.

So I think it's a matter of

as Adrian said before, of what works on average better.

Yeah.

Um and then yeah, I I would add to what Adrian I I I feel the same way in that there's sort of an interesting uh kind of tension between like you want, you know, the the

software to be as user friendly as possible and and but you know you also don't want i if if a

some posterior is really horrible and um the user just presses go and you know you you'd sort of like there to be some feedback that okay, you know, we we can sample this, but you

should know that, you know, we, you know, we had to whip out this normalizing flow to do it.

So you should also maybe consider, you know, reparameterizing or something.

I know also there's some work uh on automatic like reparameterization before like even pre pre not in the sense of how we've been talking about it, but

pre-sampling, you know, via a PyTensor or something if if you know they see see some pattern in in a model to um just rewrite it, which I I don't know too much about, but I I

I believe that's that's going on, which would help too.

Yeah.

Mm-hmm.

Yeah, yeah, yeah.

Uh that definitely is going on.

I don't know if you're working on any of that, Adrian.

Um if you are we can we can touch on that topic uh a bit later in the in the recording.

But um

I think actually it's uh it's uh good time to also talk a bit about the the benchmarking that you're doing in the paper, uh because you're actually doing exactly that kind of

benchmarking with the the lowering plus diagonal mass matrix adaptation uh in in the paper where you're you're um comparing a lot of models like ha like if I remember closely more

than one hundred models and you're showing that the

low rank plus diagonal approach outperforms the the industry's standard by a big factor.

Don't remember the exact one, but uh it's it's big.

The median difference is like a factor of four.

So it's Right.

Sometimes it's much, much more than four, but uh Yeah.

Yeah.

So were you surprised by by how big the gap was?

And do you have any do you have any idea of why

that adaptation works um so much better and also well addendum uh to that question is you are saying that this adaptation is more or less the um case where Nutpie is trying to

automatically guess the difficulty of the posterior.

So I'm I'm guessing this is this is this is related to that.

But yeah, like can you explain also the intuition behind why you're saying uh that

Yeah.

So what the low rank one does fix that the diagonal one doesn't fix is correlations.

So static correlations.

So you can have correlations that change depending on some other parameter.

It won't will not fix that.

It will also not fix funnels or whatever other kind of weird, crazy geometry you have.

But if you have correlated posteriors, it can fix those.

Up to a point, there's also kind of if

especially if you have a very, very large model with lots and lots of parameters, then kind of estimating what the correlation is gets harder and harder.

So it will fix start fixing things less.

but still correlations just are very common in models.

We know that for various reasons.

For instance, if you use if you have hierarchical linear regressions and you use centered parameterizations, which is kind of usually what I write down first

And you have lots of data, you end up with strong correlations.

sorry, if you use non-centered parameterizations, then you often end up with very strong correlations.

Or if you have over parameterizations, you end up with often very strong correlations.

So there are lots of reasons why many models will have large constant correlations in them.

And fixing those it it makes sense that fixing those helps a lot.

so in in that

that regard it's not not surprising.

Then the question is kind of how common is that actually in models?

And that kind of goes down, get gets back to kind of how difficult it is to actually benchmark samplers.

Because kind of if you look at papers, it's quite common that papers will cite, here are five models.

I tried it on those five models and hey, it's it works.

Yay.

But I mean you don't even need to to be maliciously

cherry picking for that to be a problem because I mean maybe you developed your sampler with those problems in mind and you chose your hyperparameters such as works there.

That happens incredibly easily.

So you need more models to really evaluate it, but where do you get them?

Where do you get kind of what in what what fields different fields use completely different models, have completely different typical posteriors.

What do you use there?

So the best I found the best

set of pot models that we have is posteriordb, I would say, and it's it has big gaps.

It doesn't have that many really big models.

It has kind of lots and lots of relatively small models that are sometimes very easy, sometimes very hard.

And it also has larger models.

It's not not like there are none in there.

But um so

You all still need I mean it's a bigger benchmark, I would say, than most other papers have, but you still need to take it with a grain of salt because like your models might be

different than the models that happen to be in Posterior DP.

But in general, I would say yes, it makes sense that um the low rank version works better in many actual models.

Yeah, I think that's really uh a takeaway for most of the listeners here, right?

It's like first, honestly if you're if you're using

Pine C models just use the Nutpie sampler, honestly, by default.

that's that's what I do in in my own workflow, and and you're gonna get like even out of the box, you're gonna get uh a really sizable speed up on on CPU, and then if you're using

GPU even more, um and you have all these different master matrix adaptations, and then well definitely try the

the diag diagonal plus low rank uh mass matrix adaptation.

I would say even like almost out of the box.

Let's I'll definitely do that now.

each time I have uh I have note pie on my on my environment.

So yeah, that's a great takeaway.

Um Elliott, I know you worked on that too, so I'm guessing you have uh very interesting thoughts also about what we just talked about.

Yeah.

I would I I think I would also add

uh that I I believe posteriordb is open, you know, so if if you, you know, have a model that isn't running well and or or is sampling well, I think you can always add it to

posteriordb.

And I think there's a lot of, you know, m more models.

It it would be great to have more models in there to reflect uh like the diversity of of models that are being used in different in field different fields and whatnot.

Uh so I think that's very much like a

you know, an ongoing thing.

I think there's also one in in Blackjacks, they have inference gym as well.

which I think has some very different models from what I if if I remember correctly.

yeah, and and I know there are some people working on on, you know, building those up so that they reflect, you know, kind of equivalent to what they have in machine learning with

all of these, you know, huge benchmark data sets be great, you know.

To keep adding new posteriors and actually can you we we haven't talked about posterior TB mo mostly on the show, so uh I think a lot of users, um listeners don't know what we're

talking about.

So can you actually define what what this is uh Eliot and it's basically just a repository with a couple of like a hundred and twenty or s something like that.

Mostly stand models.

Yeah.

Uh okay.

I think a couple of them are ported to PyMC.

More are ported in a pull request somewhere, but it's most of the models as I said are in in in Stan.

with the corresponding data sets.

So you can kind of um mix those up and um and and and run those.

So not nothing particularly fancy, but really, really valuable.

Yeah.

Yeah, it's a it's a great effort to try and benchmark the models, basically.

the

Not the models, the samplers.

Uh and as you were saying, it's a very complicated job to make an apples to apples comparison and basically these um this effort is trying to to make that right, to do that

right.

Um and actually I'll uh put a link to Poster Tb in the in the show notes.

I think it's uh interesting for for listeners.

and well Adrian and Eliot, you you both mentioned about that low rank plus diagonal approach.

So

It works incredibly well in a lot of cases, but also on some problems it doesn't work at all.

And it can fail spectacularly.

But you don't have a clear idea and understanding of why right now.

but what are your hypotheses?

You know, what are what are you thinking the problems might be?

Eliot, maybe you can start and then Adrian?

Sure.

I I think yeah

As I said before, like it's always gonna be an issue of averages and and so sometimes uh for specific models, you know, you're always initializing somewhere that can be completely

arbitrary and you know once you're off on the wrong track, you know, things could just get worse and worse.

Uh and so some, you know, like unholy combination of adaptation mode, uh local geometry, initialization points and posterior

Um that would be my hunch.

I th I think probably like it I don't know if it's something that's systematic or just uh it would require like a inspection on a model-to-model basis, but yeah.

I'm not really sure.

What do you think, Adrian?

Yeah, so about the it failed spectacularly.

Uh it's been a while that I've seen one where I would kind of use

way where I would say it like that.

um So we definitely had some in the past where there were just some overflows.

So you could could somehow smuggle a nun through the tests, some kind of through the the cracks somewhere and you ended up with a nun and then everything failed.

um But things where it works less or kind of just maybe the maybe kind of you change from diagonal to low rank

And suddenly your tuning isn't long enough anymore.

And some something didn't actually converge.

So you'll get horrible results or something like that.

That definitely does happen.

And I think that that'll just to be entirely honest, I mean we can fix individual things and hopefully get more of them right, but there'll always be something something that's uh

Right, yeah, yeah.

for sure.

Yeah.

Yeah, as Eliot was saying, it's a it's a matter of like of averages and

I mean not not even only averages, but mostly distribution of cases, right?

If if you're able to cover the bulk of the distribution of cases, that's the goal of NetPy, right?

It's not.

It's not a customized uh package for just one kind of model.

So you you you cannot really zero in on on these kind of issues.

Um something actually Yeah.

What one one thing I've noticed quite a well

Quite a lot as well is kind of if people compare Nutpie and Stan or Nutpie and other other samplers, that they it's relatively common that people kind of compare that, try to

compare cases where both samplers haven't really converged well.

So where the effective sample size is actually still four in some cases for for one parameter or something like that.

And

I think that's an important thing, kind of how how do you compare that?

Because I mean it can still be meaningful, right?

One sampler might get much, much closer to the actual posterior than the other one.

But if you don't know what the actual posterior is, it's really hard to tell which one it is.

So Yeah.

Yeah.

Um So if you compare samplers, be careful, really, really careful if they if none of them have converged yet.

Converged in the sense of

Decent is have effective sample size for all parameters, not just all but one parameter or something, but all parameters.

No divergences and and so on.

Um if you have that, then you can compare the samplers.

If you don't have that, it's still interesting which sample does better, but it gets much, much harder to figure that out.

And to the point of right now I don't know how to do that.

So

In the benchmarks, for instance, we just excluded everything where none of the samples converged.

We just threw those cases out because they are important, but I don't know what to do with them.

Um So I think that's also an an an interesting thing to keep in mind.

Yeah, yeah, definitely.

Okay.

So it's like when you're comparing samplers, you you have to be even more careful of having all the green lines, right?

Because when you're

Working on a model for for a project, you don't have to have all the green lights as long as the parameters that are the most important for your study are are having the green

lights.

And if you have some orange lights for some nuisance parameters, that's that's fine.

But if you're comparing the samplers, you really, really want to be in the best of cases.

And

Elliott, actually I'm I'm curious if you if you know about that and Adrian also for you.

What I'm the intuition I got about the um the diagonal plus low rank adaptation is you do need a bit more time in tuning, right?

Uh for it to be efficient.

Um or is that is that the wrong intuition?

I think that's right.

'Cause you need a an eigen decomposition, uh, which is O of N cubed, if I'm not that's right, I think.

Um yeah, yeah, so there's a little bit of uh additional like work, I think.

But I I mean I think it pays off really well as evidenced in the uh in the paper.

Yeah, for sure.

Yeah.

Yeah, yeah.

The idea here being like kind of amortizing the assembling time.

with uh a bit longer tuning time, uh which in a lot of cases definitely pays off.

It's like yeah it's like yeah eating healthier than you real quick it's not gonna leave longer.

So the if if we're talking about kind of the performance of the low rank one that's not n cubed.

So that's a little bit better because of the low rank structure that we can make use of.

Yes.

Yeah so if you do an eigenvalue decomposition of the full of kind of a full

mass matrix that's n cubed or um the dimensionality of the problem cubed.

Um the low rank one is better than that.

But yeah.

It's still it still costs something.

So yeah, yeah, for sure, for sure.

Okay.

Yeah.

Yeah, but thanks.

Good uh good precision.

Thanks, uh Adrian.

And so I'm like to zoom out a little starting to to play us out here because uh well

Starting to get late for you, Adrian.

Uh so I wanna be mindful of your time, but where where would you see your paper and just uh the the research we just talked about, where does it fit in the bigger picture of

Hamiltonian Monte Carlo research right now?

And other other adaptations or preconditioning ideas that you're watching closely uh and and you think have a lot of potential?

So in general, I think there are lots of ideas floating around that are interesting.

I'm not sure if we have time to list them all, that because there are many, many cu really cool things around, I would say.

But on the preconditioning one, I would say there's not that much from I I I don't want to c maybe I'm missing something important.

but uh so there there was a paper that's definitely important, I should add that, that looked

like s at something that's a bit similar that to the NutPie mass matrix, but just using the scores instead of scores and and draws.

But in general I would say it feels to me like we we kind of have kind of the state of the art there right now.

I would say.

Maybe somebody else disagrees and I don't want to kind of mess that up in any way but um

On other parts of sampling, I think there there are lots and lots of really interesting things happening.

There's the sparse nuts, for instance, that was going on the on the dis in Standos Discord.

There was a bit of a discussion about that a while ago, for instance.

That's more about the idea of using more information about the model.

So right now or most samplers that are used are really black, kind of Stan and Py and C at least, are really black box samplers.

So they

Take the log p function and nothing like log p function and its gradient and nothing else.

But there's lots more information in the model in a sense.

And it would be really, really useful to use that for preconditioning, for sampling, for all manner of ways.

so I think that's a really big topic uh for sampling in general.

Then kind of on the HMC side, there's for instance MCLMC.

kind of the microcanonical sampler, I think that's really cool.

I'm not entirely sure how much better it actually is.

It defin f um I can definitely come up with example where it's better.

um That's really interesting line of work, I would say.

then there's WalNUTS, for instance, kind of the idea of choosing step sizes dynamically, because kind of in HMC as we implemented right away right now, we're always using constant

step sizes.

And in in some ways that's just dumb.

We we we just shouldn't

We d we just should say, okay, here we need small step sizes and here we need need large step sizes, which every ODE solar will do.

But that's hard to do while still keeping the posterior exact.

And WalNUTS is a way of doing that.

I think that's really cool.

also I recently came across a boomerang sampler, which is also kind of com the zigzag sampler, kind of that that family of things, if you've heard that.

Um also has an HMC like variant, for instance.

And I've been thinking that could also be a dynamic step size HMC in a sense.

So that that's cool as well.

Then then there's the whole line of work of using GPUs better.

So kind of having many, many chains instead of only only few.

Yeah.

So um I I'm forgetting things, but the the general gist I would say is of

Things I've been thinking about of uh other people's work I've been thinking about.

But lots of cool things.

Yes.

Yeah, yeah.

That that that's great.

Thanks for this uh for this list, Adrian.

Um I definitely need to put a link to WalNUTS the in the channel, so I'll do that later for sure.

Um and by the way, for those of you looking on YouTube, uh it's amazing now because of Sunset.

Uh Adrian looks like a black shadow.

Um

It's amazing because you look like you're, you know, one of these um anonymous uh testimonies in uh in documentaries.

So I I'm gonna make sure to anonymize your voice for that part of the of the episode where like you're giving us the dark secrets of the of the HMC world.

clouds all day except for right now.

So that's damn.

Uh yeah, Eliot, uh anything you wanna add on that, you know, like ideas that uh that you're paying attention attention to in the adaptation and preconditioning world?

um yeah, I mean I I not so much uh adaptation and and preconditioning.

I know that um diffusion models are like all the all the rave right now.

and I think there's a lot of people working on

how to use those to precondition densities.

I I personally don't know much about about those.

Um nor I I'm also not sure if there are actual like software project, you know, open source projects for wide use, you know, uh with those implemented.

Um yeah, I think I guess broadly what I'm most interested in moving forward is different ways of parallelizing

MCMC, I think, and making it uh you know, scalable to ultra high dimensions.

Uh, you know, because MCMC is inherently like a very sequential algorithm w while meanwhile a lot of the you know computational advances that we're getting in hardware are

for parallel uh computation.

So I think that's kind of an interesting there's there's some interesting ideas on how to bridge that gap.

But I I would agree that for for preconditioning, yeah, it seems seems like not pious state of the art.

Yeah.

Mm-hmm.

Yeah, I was gonna I was gonna ask you actually, uh, you know, what you're gonna what are you gonna focus on in the in the coming month, uh things you wanna learn in particular or

or explore?

Yeah, I I'm particularly interested in it uh a new way of

Parallelizing MCMC over trajectory length.

So there's um so there's a method of uh you basically ch can solve uh an optimization problem or to solve a full MCM2 trajectory instead of computing uh in sequence, you can

make an initial guess for a whole trajectory and then iterate on it, and it provably will converge to

the true trajectory, which uh I think could be really cool for problems where the chain actually is uh somewhat predictable to some extent, but the log p evaluations are uh

really expensive, for instance.

I think it doesn't really work it it's I mean part of what I'm interested in is how that can be done for adaptive uh samplers where it becomes a little more complicated.

So that's one thing I'll be focusing on a little bit.

Yeah.

Yeah.

Yeah, sounds like fun.

And and Adrian, what about you?

Because you always do always do some fun stuff, investigating some new ways of doing things.

Uh so yeah, like I'm I'm curious what's on your mind for uh the coming month.

Yeah, as always I've got more projects running than I can actually finish, or I don't know, so that's that's always a little bit of a problem, kind of focusing on on particular

things that I can actually get done.

But um as I mentioned, MCLMC and kind of the non adjusting methods is something I've been looking into.

For some reason I haven't really gotten it to to work as well as I want it to.

So I'm not sure how much um say maybe it's it's entirely possible I'm just doing something wrong, of course.

So I don't want to bash on the methods at all in any way.

But so far I didn't have too much luck, but that could be me.

so it's a question of how much I I want to try that more.

Um definitely I'm looking into kind of what I mentioned a little bit in the beginning is of kind of using higher order autodif um as well.

So not just the Fisher divergence, but what happens if we also use second order derivatives.

I think there's really cool theory that could could tell us how to do that.

so that could be really, really useful.

interesting.

No no no sh no clue if it's actually will turn up to turn out to be useful, but fun for sure.

and then probably looking more into WalNUTS and this boomerang sampler idea um than that that I mentioned.

Kind of two completely different things, but both I think can be thought of as um the the the WalNUTS foot for sure, but the boomerang as well, kind of an adaptive step size thing.

Um which I think could could really help a lot because I

think cases where that would be really useful are quite, quite common.

um And kind of this idea of using more model information, that's also something I really, really like.

uh I've been looking into that a little bit with with Ricardo, where we now have a little implementation that can choose between centered and non-centered parametrizations

automatically, for instance.

Yeah.

Kind of moves that into the mass matrix adaptation, which I really like.

Um but they have many more ideas and and then Yeah.

Mm-hmm.

Yeah, I think that last one actually is a is a very fun one.

I saw it this morning uh in the Discord uh of PyMC.

So maybe can you can you expand a bit on that?

I mean it's something we've already mentioned on the show with you and and and other uh guests, but that's great to actually see an implementation of it that seems to be working.

Um so yeah, can you can you tell us a bit more about that and and if people can already try it?

Yeah.

Yeah.

So centered non center parameterizations are one of those things that if you've seen them once, kind of it's relatively easy to understand how to how to do that.

So you ni take a normal distribution, you either say it's sigma equals some val it's kind of the some standard deviation that you want, or you set Zigma to one and then you scale

the normal variable up with your standard deviation.

And it's the exact same thing, it doesn't make a difference kind of meaning wise.

It's boring.

But it can really, really influence sampler performance.

So if you do it wrong, and wrong means depends on your data set, right?

It's not like wrong you need to always use non centered or something, but sometimes you need to use centered, sometimes you need to use noncentered.

Basically, if you have lots of information, you want to use centered, so lots of data for that particular parameter, and otherwise you wouldn't want to use non centered.

So that's a mess.

That's no fun at all.

But turns out kind of this general framework that that Nutpie uses to learn transformations, we can use that to choose that automatically.

And not just the extremes of centered or non-centered, but there's a parametric family in between, you can say halfway noncentered and halfway centered, or yeah three quarters,

whatever.

And that's one parameter that you can just think of as a parameter for your parametrization.

um a transformation parameter.

and you can basically treat that as part of a normalizing flow if if you want.

A really, really cheap one.

Really one that's easy to evaluate.

Um so much nicer in many ways than kind of the s the big neural networks or something.

and you can just learn it during sampling.

So the implementation that we have, it's still in a PR, it's a draft, it's not done yet in in

It it works, but it's not uh not deployed in any in any meaningful way.

It just looks at a PyMC model, looks at normal distributions that have standard deviations that are random parameters, and adds kind of parametrization parameters in a

sense that say do we want to have this as a centered or non-centered one?

And then we add that in the mass matrix adaptation and just learn, well, where where where should we place the dial and um automatically choose that.

Yeah.

Yeah, this is this is awesome because basically the idea would be

Well the user writes whatever model they want.

Usually they will write the centered parameterization because it's much easier to uh code and write and understand for humans.

But then once you kick it off once you kick off sampling, you Pine C would give that to to Nutpie and then Nutpie under the hood during the mass matrix adaptation would be doing that

for parameterization between centered and uncentered, which would be like a small

normalizing flow adaptation, so adding a bit of a cost in tuning time of course, but then uh really rewarding you for for sampling time.

And also just making the model much easier to write and and sample from.

Yeah, also kind of if you do it manually often you need to re sample, right?

So I'm not sure kind of ex saying it's taking sampling time is actually true because kind of it's just sampling time you now don't need to spend yes changing your model.

And

Completely.

And you're trying and you're also getting a better scenario where most of the cases, I'm guessing, uh are gonna fall in between the pure centered or pure non-centered.

You know, just pro I don't know yet, but that's an interesting question.

I'm actually curious what what how that turns out.

Yeah.

Yeah, exactly.

I mean, I guess from a prior perspective, I think it's easier to bet on the on the fact that the distribution is probably more like inverse U shape instead of U shaped.

But yeah, that'd be my bet, but I'm uh I'll I'll be happy to hear if I'm wrong.

And most of the cases are like pure-centered or pure non-centered.

Um also something that's useful, I don't know if you have that already in the implementation, but it's the idea that you don't have also to have fully centered for the

whole model or fully non-centered for the whole model, where it would be like some clusters

would be maybe fully centered, some others would be fully noncentered, some others would be in the middle.

Would that be working like that already?

Yeah, it it's doing that.

It's it's doing that already.

So it's kind of on a per per variable basis.

Not even per variable, but finer grained.

So if you have one big normal distribution with some crazy big shape, then it can can learn that okay, some of those individual normal distributions kind of inside the big

normal should be percent and some should be uncentered, non-centered.

Yeah.

Yeah, yeah, okay.

Yeah.

Yeah, which is awesome.

It's basically like the case you want where you're able to infer the whole continuous space of that spectrum instead of artificially reducing it to a binary choice, which is

what we're doing right now, but just because of the limit technology limitation or maybe just like time limitation.

Whereas like if you're modern.

I I've actually used some some versions in between.

It's not too hard to do kind of manually either.

It's just a pain.

So I don't think nobody does it because of that.

Yeah, yeah.

Yeah, because it's just like so many cases to try.

So I think it's it's it's a perfect use case for automaize automatization for sure.

Um awesome.

Well that's that's super cool.

Eliot, anything you wanna add on that before I ask you the last two questions?

Let's let's go ahead.

Okay, awesome.

Well uh

Adrian, you've had already these last two questions.

So if you want, like you can you can answer them again, actually.

uh but uh Eliot, it's your first time, so you know.

Um first one is if you had unlimited time and resources, which problem will you try to solve?

Like any And careful trick question because it doesn't mean just invasion stats you got me with that last time.

Yes, okay.

Well I would say maybe

Maybe something which problem?

I maybe like the clean energy problem.

Like how do we, you know, convert w to one hundred percent clean energy.

That'd be that'd be cool.

I'd I'd probably work on that.

Okay.

Okay.

So I guess the closer closest to that would be fusion, something like nuclear fusion right now or something like that.

Hmm.

Okay.

Yeah.

Well, that sounds that sounds very interesting and important.

So I I love that.

Um Adrian, you wanna you wanna change your answer from last time?

Maybe you've solved already the your first problem, so you you want another one.

So so so the sampler is now actually a solved problem.

Okay, that's that's nice to hear.

Yeah, exactly.

Yeah, I guess somehow fixing the democracy issues that we are facing, but I'm not sure how infinite resources would help with that.

Seems more Yeah.

That some people have infinite resources seems to be a pro more part of the problem, I don't know.

gonna take longer than a few episodes, so I guess when you come back you'll you won't ha you won't have to answer this question again.

and second question, Eliot, if you could have dinner with any great scientific mind, dead, alive or fictional, who'd it be?

Maybe it does can can I say like Kurt Girdle?

I mean I he's not really a scient scientist, I guess mathematician.

Um I I think that that would

Maybe that would be my uh especially since a lot of the work he did was when he was like twenty two years old or something, so I think.

Very um I I would say him.

Yeah.

I mean a great choice and you're the first one to to answer for that, so uh well done.

Um Adrian, you wanna maybe uh have another dinner uh with somebody else.

Yeah, me neither.

But we can go back to your transcript and and and check.

I guess it was probably was Bertrand Russell, so I'm not gonna choose him this time.

That's possible.

Maybe let's let's pick Euclid for instance.

That that that sounds like fun.

If i if we assume that language barriers are not not not not part of the of the thing there, but because otherwise they would be a very very boring dinner if

Okay, so I will check I will check who you chose, but actually you know what?

I think you picked Euclid now that you you Yeah.

I remember you saying something about language barrier or something like that.

Yeah.

We'll see.

But I mean, if you did pick Euclid and then did that again independently, that it means like it's really a re a r um revealed preference.

So you know.

Um we definitely need that dinner to happen.

Hopefully hopefully next time you will have had it.

Yeah.

Time machine is almost done, but still kind of the last the the the last ten percent it's always it's it's always the the last little bit.

Yeah, yeah, you know.

It's the time machine is the engine that's missing some some some tiny piece.

Yeah, yeah, yeah.

Yeah, time er time travel is really an endeavor where uh the prieto effect really doesn't kick in, you know.

It's like eighty percent of the result doesn't it doesn't get you there.

Um

And actually nobody answered time travel for the first question um in the history of the show.

So you know.

Uh I I am surprised actually now that we talk about it.

Maybe it's not that important.

That may be why.

It's all it's only important to me.

Um awesome.

Well, guys, thank you so much for taking the time.

I was pleasure, Adrian.

Always a pleasure to have you on the show.

Um come back any time.

Eliot, really

Great to have you for the first time on the show.

Hope you enjoyed it.

And uh also you're welcome anytime.

Uh thank you to uh Buck Carpenter also for working with you guys on on this paper and all these topics.

uh I know he's he's doing an amazing job uh at the Flatiron Institute and and the whole team.

So, as usual, there will be a lot of links in the show notes.

Uh feel free, guys, to add any any links.

to papers or tutorials or things like that that you think are are worth sharing.

And thank you again for taking the time and being on this show.

Thanks so much.

It's been a blast.

This has been another episode of Learning Bayesian Statistics.

Be sure to rate, review, and follow the show on your favorite podcatcher and visit learnbayesstats.com for more resources about today's topics, as well as access to more

episodes to help you reach true Bayesian state of mind.

That's LearnBayesStats.com.

Our theme music is Good Bayesian by Baba Brinkman.

feat. MC Lars and Mega Ran.

Check out his awesome work at bababrinkman.com.

I'm your host.

Alexandre.

You can follow me on Twitter at Alex underscore Andorra like the country.

You can support the show and unlock exclusive benefits by visiting patreon.com slash learnbayesstats.

Thank you so much for listening and for your support.

You're truly a good baby and change your predictions after taking information.

And if you're thinking I'll be less than amazing, let's adjust those expectations.

Let me show you how to be a good daisy.

Change calculations after taking fresh data.

Those predictions that your brain is making.

Let's get them on a solid foundation.

Key Takeaways

Mass matrix adaptation is best understood as an automatic, fairly dumb, but very effective reparameterization of your model. The simplest version, the diagonal mass matrix, just rescales each parameter so its posterior standard deviation becomes one, which is exactly what you'd do by hand if you had the patience. Every time you sample a PyMC or Stan model, this kind of reparameterization is happening under the hood.

Stan and PyMC's default sampler only use one source of information for diagonal mass matrix adaptation: the posterior standard deviation estimated from warm-up draws. Nutpie also uses the gradients of the log density, which HMC is already computing at every step to build its trajectory. For a standard normal distribution, the covariance of the gradients is exactly the inverse covariance of the draws, so Nutpie takes the geometric mean of the two resulting standard deviations. There's no guarantee it's always better, but in practice it usually is.

Preconditioning HMC means transforming your target distribution into one that's friendly to sample, but doing that well requires knowing things about the distribution, like its covariance, that sampling itself is supposed to discover. This chicken-and-egg problem is usually handled by sketching a rough estimate from a handful of early warm-up draws, which can burn a large share of total sampling time. Adrian and Eliot's paper formalizes how to make better use of a second signal, the score function, that HMC already computes for free but that Stan-style preconditioning ignores.

The diagonal version is the fast default, essentially free overhead, and the right choice when your posterior has no meaningful correlation structure. The dense version is mostly a mathematical stepping stone used to derive the third case and isn't implemented, since a full mass matrix scales too badly with dimension. Low-rank-plus-diagonal is the practical middle ground: it rescales with the diagonal first, then corrects the most significant correlations using a low-rank structure, at the cost of a few extra seconds of tuning. Beyond those, normalizing flow adaptation is the nonlinear, more experimental option for geometries a linear mass matrix simply can't fix, like funnels.

Across 114 models, drawn primarily from posteriordb, low-rank-plus-diagonal adaptation beats the standard draws-only approach by a median factor of four. The result isn't surprising in direction: static correlations are extremely common in real models, often introduced by the modeler's own choices, like centered versus non-centered hierarchical parameterizations or over-parameterization, and low-rank correction targets exactly that. One real example: Eliot had a Gaussian process model that took an hour to sample with diagonal adaptation and dropped to seconds with low-rank-plus-diagonal.

Centered and non-centered forms of a normal distribution are mathematically identical but can produce very different sampler performance depending on how much data informs a given parameter, and picking correctly by hand is tedious enough that few people do it consistently. The new approach frames the choice as a single, extremely cheap interpolation parameter per normal distribution, effectively a minimal one-parameter normalizing flow, that gets learned automatically during mass matrix adaptation. It works on a fine-grained, per-subvariable basis rather than as one global switch for the whole model.

Switch your sampler to Nutpie by default; it typically gives a sizable speedup on CPU and more on GPU, with essentially nothing to lose. On top of that, try the diagonal-plus-low-rank mass matrix adaptation nearly by default too, since it usually helps meaningfully and rarely hurts. Both guests are candid that no adaptation is magic: some posteriors will always defeat a given sampler configuration, and if you find yourself spending more time tweaking sampler settings than improving your model, it's often a signal to fix the model instead.

Related Episodes