Video summary
Variational Inference represents a significant evolution in probabilistic machine learning, building upon classic algorithms like K-means and Expectation-Maximization while aiming to induce structural efficiency for faster convergence. Rooted in early 20th-century physics and the work of figures such as Richard Feynman and Edward Teller, this method addresses scenarios where computing exact posteriors or Evidence Lower Bounds is computationally intractable. To overcome these challenges, Variational Inference imposes a factorization assumption on the approximating distribution, effectively treating variables as independent through the Mean Field approximation. This approach transforms complex optimization problems into an iterative process where each variable's approximate distribution is updated based on the expected energy contributions from all other variables, creating a self-consistent loop that simplifies the mathematical landscape.
The historical workflow of developing these models involved manually deriving algorithms by defining joint distributions over data, latent cluster assignments, and parameters using conjugate priors like Dirichlet and Gaussian-Inverse-Gamma distributions. By computing expected values under the mean-field factorization assumption, the mathematics naturally induced a specific structure where latent variables became independent and parameter posteriors emerged as conjugate distributions. Although this manual derivation process was tedious, it yielded highly efficient and bug-free algorithms capable of automatically discovering the optimal number of clusters in a dataset. This technique served as a precursor to modern attention mechanisms, demonstrating how inducing factorization could lead to powerful, automatic model configurations without relying on brute-force computation.
Over time, the machine learning community shifted away from these rigorous derivations in favor of simpler gradient descent methods commonly used in deep learning, inadvertently losing the valuable "induced factorization" structure that made earlier algorithms so effective. The speaker argues that this transition resulted in a loss of functionality, as the concept was only later rediscovered within deep learning architectures under the guise of attention mechanisms. By abandoning the structured approaches of Variational Inference for generic optimization, the field missed opportunities to leverage these inherent mathematical properties for more efficient and interpretable models.
Ultimately, the demonstration of the derived Variational Inference algorithm shows how it iteratively updates responsibilities and parameters to converge on an optimal model configuration, bridging the gap between classical probabilistic methods and modern neural network architectures. This journey highlights the importance of understanding the underlying mathematical structures that drive learning algorithms, rather than treating them as black boxes. The legacy of these early derivations persists today, reminding practitioners that the elegant solutions found in physics and statistics can still inform and improve contemporary machine learning practices, ensuring that future developments continue to benefit from induced structure and efficient convergence.
Read the full video transcript
[Music]
On Tuesday we spoke about K means and
then about EM two classic algorithms. K
means is 100 years old. EM 70 years old.
We're now making our way to 40 years, 50
years old, 20 years old and then uh next
Tuesday we'll be at five years old.
The idea behind these fancy
mathematical ideas that I have had to
compress in these difficult equations is
on the high level that inducing
structure in a probabilistic model can
be used to create algorithmic efficiency
to construct
programming structures in your code that
allow much faster convergence of much
more expressive models.
So the first instance of this from 100
years ago hubush dynos is k means which
I had up there. This is what I spoke
about on Tuesday. It's an algorithm to
fit a mixture model specifically a
gausian mixture model where the gaussian
clusters have no coariance or standard
unit coariance but just a mean that gets
fitted. And we found that this algorithm
which has an iterative structure
alternating between assigning um points
to clusters and then fitting clusters to
the points assigned to them and then
repeating that that amounts to
maximizing
a log likelihood of such a gausian
mixture model. actually maximizing the
log likelihood of
the parameters of the model which are
the assignment of uh points to clusters
and the parameters of the clusters. In
the case of K means by the way this lock
probability for the data given the
clusters is just a gausian with a mean
and a standard unit coariance. But this
could be something else as well for
other data types.
And then we made a step and said, "Oh,
actually there's a refinement of this
algorithm that doesn't maximize the
complete data log likelihood." So the
likelihood involving both theta and zed,
but it maximizes the expected complete
data log likelihood. And that's called
the EM algorithm. So we maximize
this expression but under an expectation
over the latent quantity.
And then it turned out that the right
thing to take the expectation over is
the posterior distribution over zed
because then at the time when we set
this distribution
this posterior to the to the postivia at
theta hat and then optimize with respect
to this theta we are actually closing a
g a b a the gap between a lower bound on
the thing we're trying to maximize the
log evidence for x given the parameters
having marginalized over zed and
the so-called elbow the evidence lower
bound which is actually this expression
here on the left and then because we
close this gap and then raise it we are
uh guaranteed to increase the marginal
lo likelihood of the model in every step
and then therefore get a good model in
the end.
So that is an algorithm. So this you
know um I've introduced it as as having
been invented by Hamato Hartley and um
and Dempster and others in the 50s.
What I want to talk about today is the
the direct continuation of this idea
in in statistics at least throughout the
second half of this of the 20th century.
And um what I want to highlight is also
its sort of story of how it influenced
the development of the Nassen field that
we now call machine learning.
So let's to to do that let me briefly
remind you of how the EM algorithm
worked. So um our goal is we we someone
has given us a model. We want this model
involves a probability for some
observations given some parameters. We'd
like to know what a good what what good
parameters are for this model. We want
to do maximum likelihood inference. So
to do that we invent some new variables
called latent variables zed that we
think somehow um simplifies the
computation um and there is usually some
ingenuity involved in coming up with
those then that means we want to
maximize this expression. This is called
the complete data log likelihood because
there's now this new object that we've
invented. Um and the way we're going to
do this is that we first compute a
postivia over these latent variables.
Hopefully we've chosen them such that we
can actually do that. Um then given a
particular current estimate for the
parameters theta, we set a current
estimate for the distribution over zed
as the posterivia given those
parameters. That means that when we look
at this object called the evidence lower
bound, we're effectively making this
a tight lower bound on this by setting
the K divergence between our
approximation and this particular
posterior to zero. And then we raise
this expression as a function of theta
but only this theta. So the theta in Q
is kept constant. Um, and that's often
possible analytically. And then we're
guaranteed to raise the marginal
likelihood. The first thing I want to
point out is that I said maximum
likelihood inference and this is a
probabilistic learning course. So you
probably should think ah but didn't
shouldn't we do postivia rather than
likelihood because likelihoods can be
dangerous. That's completely true. It
might be better to do this in a postivia
sense rather than in a likelihood sense.
But actually doing that is really
straightforward. So if you want to do
maximum apostroori inference rather than
maximum likelihood inference, the only
thing you need to change on on this
slide even is just this red bit. So you
just add a prior for the parameters
theta. Then you're effectively
maximizing a posterior rather than a
likelihood. And um this is doing this is
very easy because you can see that this
theta does not show up in here. So this
integral is actually just I mean if you
think about what what what we add here
then we just that there's just a term
plus integral over q of z p theta dz p
theta doesn't depend on zed. So the
integral commutes out we just have an p
of theta a log p of theta* integral over
q of zed which is one. So just log p of
theta so we just add a regularizer
basically to the optimization.
Then you keep this running until it
converges. Now you have a maximum
posterior estimate. And I motivated all
of this with this f um like uh well sad
story about um brilliant mathematicians
in well central eastern Europe in the
1930s. And um actually that the path of
these people now leads us as science
tends to be because science tends to
actually happen in very small
communities to across time and space to
a different place where the mathematical
story of the idea actually continues. So
this is a picture from um um a
colloquium. So today we would say like a
weekly talk um at the Los Alamos um
laboratory during the Manhattan project
in 1946. So the second world war had
already ended. The nuclear bomb on
Hiroshima Nagasaki had already been
dropped. And you see here in the um
front row I actually made some notes so
I'm not going to tell you anything
wrong. Um
the uh in the first row Norris Bradbury
he was the director of of uh Los Alamos.
The this person you probably don't know
this is a noble laureate Enrico Farmy
he's uh probably the father of the
nuclear bomb um as such. Does anyone
reci recognize anyone else in the
picture?
Maybe the guy with the black jacket.
Yes.
Yes, this is J. Robert Oenheimer, now
Hollywood famous. In this picture, he's
42 years old. Um, clearly still a young
man. And, um, they are actually watching
a talk on what they called the super
bomb. Now we call them H bombs. Um,
hydrogen fusion bombs. Probably being
that talk is probably given by Edward
Teller who was also an immigr. Uh he
studied in um he was Hungarian, studied
in Germany, immigrated, he was Jewish,
he had to immigrate in 1934 having just
finished his degree with his wife Mitsy
Teller. Um and then he got to work
quickly at Los Alamos and um sort of had
this crazy idea that it should be
possible to build a fusion bomb. He
worked on this together with a colleague
who we remember and met already in Ruf
or Rof or Lmberg last lecture.
Stanislafulam
who was a part of the Scottish cafe.
They uh together did computations using
early Monte Carlo methods for which
Mitsy teller also invented some part of
the algorithm essentially Metropolis
Hastings but they did the computations
with an analog computer called the
ENIAC. was like a little machine that
you had to move around on a drawing to
simulate the path of particles being
scattered inside of a geometry of a
nuclear bomb. They tried to find the
geometry such that fusion would ignite
and they would get very high um
concentration of um of particles. So
there is Monte Carlo methods being used
well for impact let's say because they
created I mean overall this was the
outcome of all of this but in this
picture so this was the very first
hydrogen bomb shot this is IV Mike um in
this picture there's another person a
young guy a group leader today we would
say an independent group leader at one
of our institutions right next to
Oppenheimer does anyone guess who this
is Richard Fineman, later Noble laurate
and he made his name with work
on what we now call variational
inference in physics though he already
gave at this point he's 28 years old he
already has his PhD he's basically a
posttock um he just gave a talk recently
like a few years ago to people like
Albert Einstein and John Fon Noman who
were quite intrigued by his work he was
supposedly an absolutely brilliant
theorist and he's extending work that
was introduced as the under the name of
the calculus of variations my uh various
well old mathematicians from a long time
ago so some names associated with this
are also Isaac Newton very early version
then leonard oiler got interested but
lrange actually provided the first
formal definition of this idea of the
so-called calculus of variations that
oiler then directly took over the idea
of the calculus of variations very
vaguely spoken without showing you any
math but we'll go into the math quickly
enough is that just like we can do
derivatives on functions and then follow
derivatives and gradients in
multivariate calculus to find modes of
functions
right so to find to maximize the value
of a function
it should also be possible to construct
and the equivalent of a derivative in
spaces of functions rather than in real
valued vector spaces in spaces of
functions to maximize what is the
corresponding concept to maximize a
functional. So a function of a function
and for that we need something like a
derivative in function space that we can
then follow to maximize such a
functional an expression and LR was like
totally excited about this idea and
fineman got excited about it as well
because it was a crucial ingredient for
the theory he was working on quantum
field theory where the expressions you
want to maximize for are functionals of
incoming probability distributions over
particles. There is by the way just
before we move on just to not leave it
out also a corresponding equivalent to
um well sort of the idea that you also
know from lagrange multipliers. So
finding directions along which a
functional is invariant right? So just
as just as you can follow a gradient to
find a mode, you can also think about
spaces identified by zero gradient which
are equipotential lines. So
spaces along which something is
invariant and does not change. This is
called um the the theory of of uh
invariances in variations and it's
associated with um the work of emintor
who uh like was probably the like most
famous mathematician of all times which
also made ma major contributions to
physics identifying conservation laws
rather than optimal choices. But today
we're going to focus on the optimization
problem on the the idea of how to how to
maximize a functional. Our elbow is such
a functional. It takes in a distribution
Q,
right? That's a function of Z. And then
we want to find the function Q which
makes this expression maximal.
So this expression is clear. Oh, I
forgot the differential. So it's a
derivative in integral over dz, right?
So it's something that takes in the
input over zed and marginalizes it out,
integrates over it. And now we have a
number or actually we have a function of
x that we can try to maximize
for a particular x.
So what does that mean? So this is the
this is the physics motivation. I'll
talk more about that in a second. Where
does this like what does this come in?
What does this have to do with machine
learning? What does this have to do with
with like with inference and
approximation? Well, machine learning is
the continuation of physics into
everything else. It's the application of
the of mechanisms from science into
every other aspect of human activity.
That's maybe why it's so exciting and
why it's so influential at the moment.
So to try let me try and generalize the
ideas both from the previous lecture and
from this physics background in a few
slides. So let's and for a moment let me
try to be more concrete. So let's say
we're trying to do inference in a model
that assigns probability to some data
points X. So X is the things we get to
see and our model will involve some
parameters theta and it involves some
latent quantities which we need to
describe what's going on. We could think
of them as nuisance variables or maybe
we think of them as the things we're
particularly interested in the variables
of the model zed. So in the algorithms
we've encountered so far in K means and
EM we had to make very specific choices
for zed. So the the the name of the game
was let's find the right latent
variables such that this posterior over
zed given x becomes tractable and we can
just write it out for a particular value
of theta. When we can do that then we
our our elbow h well then we need two
things actually. We need to be able to
get the the posterior attractable. So we
need to be able to evaluate this
function. Not just to evaluate it, it
has to have a closed form. And then in
addition, we need to be able to compute
this quantity where Q is now set to P.
So let me go back up. Why did we do it?
Well, I mean, yeah. So we need to plug
for Q. We need to plug in this postivia
and then be able to do this integral.
Then we can do EN.
But maybe we can't do that integral or
we can't evaluate this posterior. That
seems like a strong constraint and make
and makes em a sort of a very specific
tool that you can only apply when you
can do a when you actually can compute
the posterior and the associated
integral for the elbow.
Maybe sometimes we can get away with
just saying somehow will allow Q to be
something more like more restricted,
something that like that some space of
functions in which we can do tractable
computations
and then we just maximize this elbow
this evidence lower bound within the
space of things in which we can do this
within the space of probability
distributions Q in which we can compute
and when we do that we'll then find some
Q that minimizes KL divergence right
because raising the evidence lower bound
means we get closer to the evidence and
as we do that we narrow the gap the gap
is given by the KL divergence between
our Q and the posterior over Z given X
and then this Q that we find within the
space of tractable distributions will be
the best possible tractable
approximation to the posterior
that we might want to have
and that's actually the entire idea. So
this is an optimization process within
some space of probability distributions
Q and that idea is called variational
inference. Well because because
variational is the word for
optimization or sort of the
generalization of derivatives to
function spaces.
And of course that raises the question
of what do I mean by a space of
tractable distributions? I'll tell you
in a second because it's actually
surprisingly
fluid this definition of of tractable.
But before I do that, let me first make
the connection to physics because I said
this all comes that's all associated
with people like fineman and so on. So
what that you understand why people
might be interested in that and actually
also to make my case that machine
learning is just a continuation of uh
physics
to everything else. So here's how this
works in physics. We um
uh um so these problems show up in
what's called statistical mechanics or
thermodynamics. So this these are this
is the study of systems that involve
many many many particles a large set of
particles or large set of things that uh
all have their degrees of freedom. Let's
say their degrees of freedom are x and
zed whatever this is like the position
the momentum of all the particles in a
gas. spin in some quantum gas and so on.
These are called the micro state of a
system. So it's the set of parameters
that or variables that describe
everything. But we assume that there's
many of them like 10 to the 23, right?
So a very very large number of them and
therefore we can't actually see them.
But we can write down how energy arises
in this system. In principle it's just
too many numbers to track. So we write
down an energy function a total energy E
that's just a function of X and Z
and then we that energy as you now know
already is associated with some
probability distribution. So particles
just take a configuration according to e
to the minus energy. That's actually
sort of a circular statement. The energy
is basically this right? So the energy
is the logarithm of the ne negative
logarithm of the probability.
So the typical challenge now in this
theory like the statistical theory of
mechanics or thermodynamics is that you
can't actually go in and measure the
state of every single particle or every
single element of your system because
it's large. There's like 10 to the 23
molecules and they all behave in a very
complicated way. What you can measure
though is what's called macro states. So
you can measure the temperature, you can
measure the pressure, um maybe you can
measure some chemical potential,
something like this. These are called
macro states. And the way you usually
measure those is by taking the system
and letting it do something. You're
letting it expand to measure its
temperature or its pressure. You measure
the force it exerts onto walls to
measure pressure, something like this.
And um so this means you're using some
energy that's available in the system to
measure the macro state in a very
abstract fashion. And that's why the the
this uh like the the the quantities that
we tend to measure in terms of macro
states are associated with something
called the free energy.
So free energy is really just a word
that physicists and chemists actually
use to describe
properties of a of a microscopic system
of a statistical system that one can use
to do work.
And uh that allows us to now separate
the x's and the zs basically right. So
so sort of as a in a sort of backward
definition we could think of the um the
variables x in our model as the things
that we the aspects of the system that
we get to measure in as a function of
that that we get get to measure a
function of when we measure macro states
like temperature um or uh pressure or
some work that the system can do. And
the zetss are the other variables that
somehow only indirectly affect the
macroates. And so we can um define
something called the free energy as the
log negative logarithm of a marginal
over our probability distribution under
set.
So it's the negative logarithm of the
distribution of x.
And um so x might actually be a very
small set of of variables and zed might
be a huge um remaining set of variables
and there are actually now um well okay
so I'll talk about free energy in a
second again. So what what you now want
might want to do is you might want to
figure out what the remaining degrees of
freedom are what the zs are that you may
still be interested in for some reason
but you can't measure them directly. So
you're going to find some probability
distribution that approximates
the correct postivity on Z, the
conditional you would get under all of
the observations by fitting some
parameters of a model. And that involves
writing down something that um in the in
uh variational
theory of physics is called variational
free energy. And that happens to be the
exa exactly the negative of our elbow
because you're minimizing it. So then
what the physicist does is they try to
find an expression that minimizes
variational free energy which is the
same as maximizing the elbow. So they're
effectively doing approximate inference
to find an approximation to the
posterior on a on uh the distribution
over some interesting latent quantities
of a system that can't be directly
measured through workable energy through
free energy. And then a large part of
the game in uh thermodynamics is that
people have come up with various
different ideas for free energy. And
there's all these different words that
if you've taken a a thermodynamics class
or a statistical physics class, you've
encountered like the helmholds free
energy or the gibs free energy. Um, and
they associated with these people. So
Rudolfph Claus was a German uh physicist
who coined the term entropy first. It
came from him. He was looking for a word
that sounded a bit like energy because
he was thinking about the ability of a
system, the differential of a heat in a
system that can be used to return heat
back and do work with it. He initially
wanted to call it the favandlunks
energy, but that was a weird word that
nobody understood. So he was looking for
a Greek word because he thought ancient
languages survive for everyone and
everyone uses the same ancient languages
at least in the western world in you
know 18 something. Um so he came up with
entropy which is the Greek word for um
transformation
basically. So favandlong's energy um and
as a contemporary said he that he
succeeded in finding a word that meant
the same for everyone namely nothing.
And let's be honest to this day nobody
really knows what entropy actually is
beyond an expectation of a logarithm.
Then he had was the first maybe to come
up with this idea of free energy. He
defined a certain type of free energy
that involves entropy of a system and
it's all about what kind of work the
system can do and how you measure it.
Then Josiah Gibbs had um also a very
like influential English physicist in in
thermodynamics had a a new way of
expanding the idea of um free energy. So
he added an extra term to it. That's
pressure times volume. He calls this um
u by the way is the potential energy
plus pressure times volume. He calls
this the enalpy. Somehow funny that they
all came up with something similar to
energy, right? Entropy, energy,
enthalpy. And then I'll do a jump just
to connect the dots to a contemporary
colleague Dave Dave Bllye from Columbia
University now in in New York. And to my
knowledge, he was the first to came up
with the word elbow for evidence lower
bound. yet another word with an E. Um,
and he's a statistician or machine
learner we would say these days. Um, and
uh um really like phrased this very
clearly as the variational energy being
the negative of a lower bound on the
evidence in a general statistical model
that we can apply way beyond physics and
maximize. He also came up with some
early very structured algorithms for
language modeling long before large
language models.
So how does this optimization in the
space of functions actually work? So
it's going to work into the following. I
now have several slides that try to
explain in like as directly as I managed
to because it's actually quite
challenging how this
idea of optimizing for a a functional
for a function might actually work in
practice. And you can maybe imagine that
this is going to be hard. Why? Because
well for functions over a vector space
we know how to compute gradients and
then the gradients are vectors in the
vector space. So these vectors are
tractable objects right? They're just a
bunch of numbers that you can follow.
You can write out an algorithm that just
follows a gradient. But if you think
about a functional that takes in a
function. So a functional is an integral
and we need to optimize this with
respect to a function. Functions are
infinite dimensional objects. So we
can't expect that we're going to compute
some infinitely long vector that we can
somehow follow on a computer. Right? So
the process is not going to be that we
will compute a derivative in function
space and then do a few steps to
optimize the values of the function.
That doesn't work. And instead we're
going to find implicit statements that
characterize a mode in the space of
functions.
And then those implicit statements will
surprisingly sometimes actually tell us
what the optimal function is. So here
this is how is this how this is going to
work. We start again with our model. We
invent some latent variable Z XR are
actually the given data and we have some
parameters theta. We want to find a good
approximation to the posterior over the
latent variable Z given X. If you have a
model
then the posterior is well defined. It
exists but it might be intractable. we
might not be able to write it down. So
what we're now going to do is we invent
some family of distributions Q of Z.
Um
and then
well we uh
that's not a that's a stupid sentence.
Why did I add this there? No, that's not
actually what I want to write.
This might be again one of the things
that I fixed. Let me just do this.
Well, yeah. Or maybe it's just confusing
and it's what I wanted to write
initially.
But let's actually
keep it simpler and just write it like
this. So we invent some family of
distributions Q of Z. And then we want
to optimize the elbow which is
this thing
as a function of as a functional of Q of
zed. I mean so Q of Z might well depend
on our parameters. Sure. But the goal
will be to find a Q a function of zed
that maximizes this expression. And this
expression is an expectation over this
function under Q of zed. So it's an
integral that gets rid of zed basically.
So here it is again when we do that if
we manage to do that somehow we don't
yet know how we're going to minimize KL
divergence to the true posterior. So
what we find will be a good
approximation in the sense of this
divergence. So if we believe that KL
divergence is a good measure of distance
in the space of probability
distributions then this process will
find the best approximation to the
posterior
within the space of Q's that we
consider. Now if you if you don't expect
any if you don't impose any any
restrictions
on Q if you just say within the space of
all probability distributions find the
optimal Q we know what the extremum is
going to be the optimal value it's going
to be this probability distribution the
posterivia but we know that we can't
compute that or we assume that this is
not tractable because if it were we
would just compute it so we'll probably
have to somehow restrict rict the family
of distributions that we want to
consider. And you might think that this
restriction will be something like, oh,
we're just going to assume that they're
all gausian or something like that. And
that's actually an option. But it turns
out that sometimes you can be much more
uh subtle. You can impose very weak
restrictions on the space of functions
that aren't even of a parametric form
and you still get an actual tractable
distribution out. And I'll show you how
this works after the break. But first
there was a question. Was it done?
>> Ah yes. So we are not so this is this is
the thing like let me let me just make
sure you're not not getting confused.
We're not yet going to assume that we
will that Q of zed will be a parametric
family of distributions. It doesn't have
to be an exponential family even. It
could in principle just be anything for
the moment.
So the restrictions we're going to put
in will be something very subtle. It
will not be oh it's a gausian with a
known variance and we just have to
optimize for the mean. That's actually
also something we can do. This is this
is uh sometimes called fixed form
variational inference and I'm thinking
that I might introduce it at the
beginning of the next lecture also in
another historical anecdote. But um yeah
so well like for now let's actually keep
it at that. So if there are no more
questions then
let's take our customary five minutes.
There was a good question during the
break which was also maybe triggered by
my mistake of like the spelling mistake
here or writing mistake here on the on
the slide. Namely what about this theta?
on Tuesday we spoke about optimizing the
elbow for theta and that's the at the
heart of the EM algorithm in in EM this
maximization of the elbow in terms of
the distribution Q of zed sort of
happens almost implicitly right we just
set Q of Z to the to correct posterior
because we can track the full posterior
and then we optimize in theta so in the
setting that we're now discussing we can
still optimize for theta that's still
fine but this is not actually what I'm
talking about right Now, so the what
what we're now considering is this sort
of maybe like to make this aspect of
finding a good distribution over zed an
explicit part of the thought process in
in EM finding the best approximation Q
is straightforward. It's just P of Z
given theta sorry P of Z given X and
theta the true posterior.
But we're now saying well what if we
can't track the full posterior? What is
another thing we could find? So we're
looking for a function of zed. That
function of zed might still also be a
function of theta. And then we can also
optimize for theta. Fine. But it's just
not important for now, right? So we'll
keep those parameters like out for a
moment and care about this thing as a
function of zed. And we're looking for
the best possible function of zed. And
we know that the one best possible
function of zed is the postivia of zed
given x. That's the thing we'd like to
compute because we want to do Beijian
inference, right? But let's just say
this is not something we can write down.
It's just too complex. I'll come back to
the example of Gaussian mixture models
later on in a few minutes and then we'll
see that there might be generalizations
of that model where this Q of Z given X
is intractable.
And now I said this ominous vague thing
is that we're somehow going to put some
restrictions on the space of functions
probability distributions Q of zed.
And again it might seem like we would
have to put some parametric
descriptions. We'd like to say we are
going to write this function in terms of
a finite number of features. But
actually just like in the gausian
process world we got away sometimes
without actually imposing a functional
form on the space of functions. but
instead using some kernel function.
We're now going to get away without
imposing a functional form but by
imposing a factorization.
So a factorization is a statement like
the one in the first row or above that
says so we're looking for this
probability distribution over these
latent variable set. So of these sets
let's say there is like a 100 right or I
don't know n of them
and what we're going to have a weird
clicker problem. Um, and what we're
going to do is we're going to impose
that Q of Z factorizes into a bunch of
terms that depend on well, let's say
some subsets of Z. So the Z I might
still be multivaried subsets or they
might just be individual scalar
variables. It doesn't really matter
actually. It's just some new variables,
some subp part of this entire set of
variables set I. And you know that
factorization amounts to independence
assumptions, right? So maybe by imposing
independence
under the our approximation to the
posterior so our Q is an approximation
to P of Z given X maybe the computation
becomes tractable and we are not going
to impose anything else we won't say
what the functional form of Q is just
that it has to have a product
over the Z I
so then what happens is the following so
let's look at our variational bound our
elbow that we're trying to maximize
I'm going to put it in um in its
definition. So that's the full
definition. But now I've plugged in the
fact that we assume that Q of Z
factoriizes.
So in our variational bound, remember
that the variational bound is an
expectation over Q over log of joint
over Q, right? Log of P of X and Z
divided by Q of zed. So if you take the
logarithm, we can take this into a
difference. And now we impose the fact
that Q factorizes. So we have a product
here and a sum here. And now we're going
to do a thing namely that we're going to
pick out one particular set
of variables set J. So there's this
product over all the Z I. We're going to
pick one of those. Let's call that set
J. I've actually highlighted this in
this slide in red, but I realize now
that you can barely see the red, right?
So this set J is going to be for a
moment the thing we focus on. Then of
course what's going to happen on the
computer is we're going to iterate over
the J's, right? Going to do this for one
J and then for all the other ones. So if
we do that then we pick out this one
term from the product. Um and then in
here let's of move the rest of the terms
of the product into the into a set of
brackets. And then this joint integral
over all the sets will be an outer
integral over Z J. and then all the
inner integrals over z the all the z i's
minus
um over here we'll get one particular
term because here's a sum over the q i z
i right so in one of those terms set j
is in there so we'll get an integral
over q of zj j log q of zj j that's the
entropy of our distribution q of zj j or
partial distribution or marginal over q
of zj j plus something I'll call a
constant. So why why is that a constant?
Well, because it involves the sum over
the logarithm of q of i.
Um when you integrate out all the set
i's that are in j and then there's an
integral over q of zj j which is just a
one. And so it's just a number that
doesn't depend on zj j. Okay? Because we
want to optimize for q of zj j. That
number doesn't matter.
So now we're going to look at this
expression and we realize that it's sort
of in a functional form. We have
marginalized out the set I we've taken
integrals over Z I all the I's that
aren't J. So what we're left with here
if you think of this as a function then
it's a function of X and of Z J because
all the other Z I not J have been
absorbed. They've all been integrated
out. So we're left with a function that
we can think of as the logarithm of an
unnormalized probability distribution
over x and just z j. So it's a function
of only set j but not the other
variables
minus uh sorry plus the entropy of uh q
of zj j plus constants.
So this is so if uh for a particular
function zj this is just this is just a
constant. So what we have left here is
an x uh sorry uh is of course not right
it's just an entropy but what we have
here is um
this expression which we can uh think of
as an unnormalized probability
distribution over zj
which is equal to the expectation
of the log joint
under all the other approximations. So
now if you think for a moment in terms
of an algorithm what we're going to have
in the end is some some for loop that in
each loop goes through all of our Z eyes
I and when it considers a particular set
J
this object here will be the expectation
under all the other approximations under
all the other Q of Q of Z I not J. So if
we find an algorithm that sort of
inductively always constructs
approximations over all the Z I's then
whenever we consider one of the one of
them Z J all the other ones will be
assigned a distribution and if we're
lucky we might be able to compute this
expected value. Now of course we don't
know yet whether we're going to be able
to do that. We just hope
right and now let's just let's just go
in and see if this might work. How would
we then find a Q a Q of Z j which
maximizes this expression?
So here we now sort of instead of trying
to comput a gradient which we can't in
function space we remember that
maximizing the elbow amounts to
minimizing KL divergence. And here we
now have a reduced elbow. We have an
elbow over Q of Z J.
And that means when we maximize this
elbow, we're going to be minimizing the
K divergence between Q of Z J and this
unnormalized well now hopefully also
normalized probability distribution
implied posterior distribution over Z J
which arises from integrating out all
the other Q of Z I and that object right
we know how to minimize the K divergence
we just set it to zero by setting this
equal to this.
So if we have this
log P of X and Z I this thing sorry Z J
this thing then we can look at this
expression and it might be of a
functional form that we can identify as
a function of Z of J for which we might
know the normalization constant and then
we're done we have a probability
distribution over Z set J that minimizes
the KL divergence maximizes the elbow
and therefore is the best approximation
within this loop that seems crazy at
first let me show you an example and and
you'll maybe get the structure. So we'll
find an algorithm that iterates
repeatedly in a for loop and inside one
iteration of the for loop. It'll iterate
over our variables set i 1 to n in each
case setting the approximate
distribution over z j to this implied
posterior distribution which we get by
computing this expected value of the
logarithm of p of x under and zed under
all the other set js.
and then looking at the expression and
finding hopefully that it happens to be
the logarithm of a distribution that we
can write down. Why might we be able to
write down such a distribution? Well,
because we're going to put in some
assumptions about what P of X and Z
actually looks like, right? We're still
going to define our generative model for
the data for X and Z in some nice
structured form. For example, we might
choose an exponential family. So if p
this joint distribution of a latent
variables and data happens to be an
exponential family parameterized by some
natural parameters then as you know we
can write it like this that's the
standard form of an exponential family
with sufficient statistics natural
parameters log partition function and
then the thing that we need to compute
this expected value under our
approximations
for the logarithm of the joint
distribution log of p of x and zed will
be this expression, right? So now we
look where the zs are, right? So the
zetss don't show up in here. That's just
a log normalization constant. Capital
zed and lowerase zed are different
things. So we're left with just an
expectation over the sufficient
statistics
times natural parameters minus log of
zed. And now we know that we want to set
this expression to zero with respect to
the parameters. Um and that means that
we're going to have to set the gradient
of the log partition function to the
expan um to the ex to this thing right
so in terms of theta we're going to have
to set it to this and we know from our
lecture on exponential families that
this expected value of the sufficient
statistics
sorry that this gradient of the lo
partition function is equal to the
expected value of the sufficient
statistics under our uh p of x and um
zed.
And that might actually be something we
can look up on Wikipedia if it's a known
exponential family. And then we have an
expression here that we can just plug in
there and then we're done.
So
here's a summary and then I'll show you
an example. So this so-called free form
variational inference
is a framework to construct
approximating probability distributions
that approximate a posterior P of Z
given X that might not be analytically
tractable. And we do that by iteratively
finding
um minimizers of the Kel D or overall
finding a minimizer of the Kel
divergence between the approximation and
the true postivia which we do in an
iterative sense by going through subsets
of the variables and imposing some
factorization. This by the way is called
the mean field. I mention I actually had
this on the slide but I sort of rushed
over it. The reason why it's called the
mean field approximation again comes
from physics. That's also a little bit
why I gave you this physics
introduction. Namely, you can think of
this factorization as saying there's
this system that has contains lots and
lots and lots of particles. There's too
many of them for me to track all of
them. And they all interact with each
other all the time across time and
space. And now what I'm going to do is
I'm just going to pretend I'm sitting on
one of these particles that J. And I
want to know what happens to this thing
that sits within the system and I want
to want to see its own distribution.
What's the best local distribution for Q
of Z J that I can think of? Well, for
that I need to consider the expected
value of this overall probability log of
P of X and Z. The overall systems
energy. Remember that log of P of XZ is
a negative energy under the effect of
all the other particles. if everyone
else in the system also pretends that
they are alone.
And then the way that they interact with
each other is solely through this update
step which you can think of as well as
we just derived an expected value
of the overall systems energy under
everyone's individual contribution and
that is a function of myself said J and
that sort of mediation of everyone else
on average to me that's the mean field
right the average contribution ution
from everyone else. That's why this is
called mean field theory. So in practice
this means we'll define a probability
distribution. Then we impose a
factorization. We do a lot of
computations to hopefully find an
expected value of the log of P of X and
Z under all of these approximations.
Look at it. Find that it has an
algebraic form that we might be able to
track. And if we can do this across all
of the variables in our model, then
we're done. And then we can just write
an algorithm that iterates. So I'll show
you now for the rest of the lecture how
this works for our Gaussian mixture
model in an extension of the Gaussian
mixture model. But before I do that, I
want to tell you a little bit also why
I'm doing this in this lecture. Because
every time I bring this up,
I realize that what I'm going to do now
next is going to be going to seem
tedious and complicated and like
mentally painful.
And that's maybe actually true. But
there is a very interesting very high
level lesson hidden behind that. And
what I want you to do now in the next 20
odd minutes is to maybe first keep the
following story in mind that I'm telling
you really like because I think it
contains an interesting insight and it
sounds like old men talking about war.
So when when I did my PhD, I already
mentioned this. What I'm going to show
you now was really the daily work of
machine learning engineers and
statisticians. Namely, it was a very
manual process. We would write down
these graphical models that I'm going to
actually I'm going to show it to you
now. Graphical models that look like
this. This would happen on a blackboard.
You would stand together with your
colleagues, write down models and lot
think a lot about algebraic expressions.
And there was a whole lot of theory
about these graphical models. This
lecture also used to be a lot about
these graphical models. There were all
sorts of concepts like induced
factorizations and um um induced
covariances explaining away structures,
V structures, colliders and so on. We
had words for all of these with the
plates and so on because the work
largely involved writing down math on a
on a piece of paper or a blackboard.
then doing the very tedious thing that
we are now going to do together and then
finding an algorithm arising from all of
this that you then had to implement on a
computer that was very hard to get
right. There was usually a lot of bugs
so they were they were quite hard to
fix. Um and then you had these updates
that we're going to see emerging now
from the computations that would give an
algorithm that actually runs that if you
once you hide it all right worked really
beautifully and then you could write a
cool paper or build some product um and
actually also build products. So this is
how some early tools for like for
machine learning were actually built
like how uh online recommener systems
were built for for large platforms like
Xbox Live and so on. Um
and um that felt like a really tedious
process and it used to be a joke that
that was like how how
PhD students would spend their time and
someone actually said well why don't we
just always use mark of J Monte Carlo
and a running g that by the time you had
implemented this algorithm your mark of
J Monte Carlo algorithm had converged
and you could have just taken the
samples from that so the value of course
was in doing the derivations once and
then having a very efficient algorithm
and because it was so tedious very
precisely I think actually because it
was so tedious and you're going to see
in a moment how tedious it was people
stopped doing it by the time deep
learning arrived and they got excited
about gradient descent and said how
about we just get rid of all of this
nasty derivations and just do gradient
descent from now on I'm going to
highlight this on Tuesday next week and
interestingly In doing that, the whole
community actually moved to a different
type of variational bounds which don't
require these nasty derivations which
are just gradient descent. And in doing
that, we actually lost the plot for a
little bit. It lost some functionality.
And I'm going to show you on Tuesday an
example, a high level example of how I
think this happened. And actually none
of us really realized I think for a
while up until someone discovered that
you need the structure again in some
deep learning architectures and it can
be very powerful and that idea was
called attention. So I'll try and
highlight that story and I'm telling
this to you now giving this spoiler
because I want you to like just keep
your attention for 10 minutes as I'm
doing these nasty derivations. I realize
that they are tedious, but there is a
lesson hidden in them if you squint
through the math for why people might
have stopped doing them. So on Tuesday,
I showed you this um Gausian mixture
model. Here's our data. Remember the
scatter plot with the dots. There's like
we think that there are clusters in
there. We think that they are described
by some Gausian cluster which have a
mean and a variance. And um each cluster
has a weight that is um parameterized by
some probability distribution, a
discrete distribution that uh is
represented by a vector of numbers
called pi k that contains numbers from
between zero and one that sum to one. I
showed you on Tuesday how we can
optimize those parameters
mu sigma pi by inventing these latent
variables set n and k which then induce
a factorization once we know what the
set nks are then it's easy to optimize
for mu and sigma and pi and then when we
have mu and sigma and pi we can again
compute a posterior over z and k that
was called the em algorithm so we
iterate between computing a postivia
over zed optimizing the parameters given
that postivia by optimizing the ex the
the elbow the expected um complete data
log likelihood and then setting again
the posterior over uh z to the the
resulting posterior if you fix mu sigma
and pi we now realize that there was a
special case of variational inference
where we did completely free form
variational inference basically we found
the optimal approximation to the
distribution over Z, namely the true
posterior if we fix the parameters. Now
what we're going to do is to do
generalize even further to a point where
the true posterior is not tractable
anymore. So why might you want to do
that? Well, one problem that this model
has is that we have to say how many
clusters there are. We have to say there
are K clusters
and what we don't know what K is, right?
2 3 4 5 6 7 10 I don't know. You have to
look at a data set and that's difficult
to do. I mean it's easy to do in two
dimensions but difficult in 100
dimensions. So what you might want to do
as a general concept um is a fully
basian gausian mixture model. This was
something that people like Neil Lawrence
and Chris Bishop and David Mai were very
interested in around the sort of early
2000s. Um and they wrote several papers
about them which I forgot to site here.
Maybe I'll do that next Tuesday. So
we're going to do a general Beijian
version of this model. What does it mean
to be basian? Well, we're just going to
replace these parameters pi, mu and
sigma with actual variables. We'll say
we'll be uncertain about what pi mu and
sigma are and compute posteriors over
them.
And how would we do that? Well, first
let's think about the algebraic aspect.
So instead of a p of x and z given mu
sigma and pi we're now going to have a
joint distribution over p of x over x z
pi mu and sigma by um the
uh product rule we can write that as the
likelihood we already know p of x and z
given pi mu and sigma. So that's an
object we already had in our gaussian
mixture model times a prior over those
three parameters P of pi, mu, and sigma.
But we're already going to impose some
generative assumptions. We're going to
say this factorizes into a prior over
pi,
a prior over mu given sigma, and a prior
over sigma. Why? Why? Why might might
you want to do that? Because you've done
a homework exercise about conjugate
prior inference in Gausians and you know
that there is a conjugate prior for
Gausians with unknown mean and variance
called the gamma inverse vishar and you
did some of you did the homework on this
with um Bernardo and we also know what
the conjugate prior for a for a um
categorical probability distribution is
does someone know it's actually on the
slide
so category ical distributions are an
exponential family. They have a
conjugate prior and the conjugate prior
is the
it's the multivariate generalization of
the thing I use to infer who is wearing
glasses.
It's on the slide so you can just point
it out.
It's the durishlay distribution. So the
conjugate prior for a categorical
distribution is called the durishlay and
it's here. So I've just written it down.
So we say and this is sort of what would
use what used to happen in front of a
blackboard, right? Someone said, "Oh,
I've drawn this graph." So this is going
to be a variable. Oh, we've got mu and
sigma. They're going to a gausian. Oh, I
know what the conjugate prior for a
gausian with unknown mean and variance
is. It's a gaus inverse var looks like
this. We're going to put independent
priors over every single component with
a vishar distribution over the inverse
of the coariance matrix the position
matrix and a gausian prior over the mean
and then someone knew from memory that
you have to do this product form where
you have a vishar prior over the
precision and then a gausian prior given
the precision that's how you get a
conjugate prior and then pi is a
probability distribution so everyone
knows what the conjugate prior for a a
categorical probability distribution is
it's a dishly distribution which has
this particular form that's a dishly
distribution
and those things have parameters right
they are themselves exponential families
because the conjugate prior to an
exponential family are themselves
exponential families they have par
parameters they're called alpha for the
dish a vector of numbers and w which is
a symmetric positive definite matrix for
the vishar and mu which is a count
variable you know you might remember
that conjugate prior have these count
variables and um um m and beta for this
Gaus inverse var and you just happen to
know that if you lived in that time and
now we realize okay that posterior over
pi mu and sigma is not going to be
tractable. Actually, we used to be able
to read this off a graph. And if if you
would have taken this class like four
years ago, you would now be able to read
it off the graph as well because you can
see that there's a collider structure
here coming in all these arrows pointing
in. This is a variable. It's sorry, it's
an observable. We've seen this. So
therefore, all of these things become
dependent on each other in the graph and
we cannot expect to find a joint
posterior over them. You don't know this
anymore because I didn't do three
lectures on graphical models and instead
talked about deep learning and atlas
approximations. But you know I can tell
you that that's something you what you
would be able to read off if you are an
expert in graphical models. So now what
people used to do is to say oh we're
going to build a variational bound and
we're going to do a free form
variational bound. So the only thing
we're going to do is we're going to
impose some factorization
over this joint distribution over these
four things.
And the factorization we probably want
to have is that the variables
sorry the parameters of what used to be
our Gaussian mixture model will be one
set of variables and these um cluster
labels these set n case ks those will be
the other variable
why well because that's how em works
right it was somehow very useful to have
this distribution over set. Our
approximation that we're going to impose
is that whatever our whatever our
approximation to the postivio over z pi
mu and sigma will be it will have to
factorized like this. There'll have to
be one distribution for who belongs to
which cluster and then one distribution
for all the variables.
That's it. That's the only only
constraint we're going to put in. And
that seems like a very weak thing to do,
right? Just one little dot between the
distributions. And now we do variational
inference. And this variational
inference with this really like cranking
a big squeaking handle on a big cog. So
you sit down with a with a piece of
paper. Sometimes was like an A3 piece of
paper in like landscape mode and then
you start with the derivation. You say
what was the rule to derive a
variational bound? Well, we will need
under this mean field assumption, this
factorization, we will need to compute
the the well to find this object, we
will need to compute the expected value
of the joint distribution
the loin
under all the other variables in our
distribution. So if we want to find the
the approximation to zed then we need to
take an expected value of the log joint
with respect to pi mu and sigma. And
then afterwards hopefully we might be
able to get an to get a distribution on
pi mu and sigma by taking an expectation
over the same thing under this
distribution over zed that we've now
just found. And then we can iterate.
So actually it might be a good idea to
start with zed. Why? Well, because the
zetss are they are binary variables,
right? So, we're going to clearly
because they're just a set of binary
variables, we will learn probabilities
for them. I mean, the most general
distribution over a set finite set of
binary variables is a probability
distribution, a discrete distribution
that just contains numbers.
Okay, so let's do those first. So our
goal is going to be let's compute an
expected value under whatever the other
distributions are. We don't know yet
what they are over this log joint. And
now we look at the expression for the
log joint. So we go back up again. We
look at it. Okay. This is the thing we
need to plug in, right? That's the
object. And these have actually forms,
right? So P of pi is the dish. P of mu
and sigma is this thing down here which
has this factorization structure inside.
And this thing p of x and z given pi mu
and sigma. Well, it's up there. this
this this thing we're we're just going
to plug that in with a logarithm in
front.
So now we find that this when we plug
this in we'll get one term over p of zed
given pi.
Clearly this only contains pi. There is
no mu and sigma. So our expectation over
pi mu and sigma actually already
simplifies. It's just going to be an
expectation over pi. And then there's a
term over p of x given z mu and uh
sigma. So that contains no pi only mu
and sigma. So that will simplify the
computation as well. We'll only need to
take an expected value over mu and
sigma. And then if you go one slide
back, there's a lot of terms here in the
back that don't actually depend on zed
at all. So we just can forget about
them. Okay? So they won't contribute to
our approximation to q of zed. So those
are the two terms that we need to
compute. And now we plug them in. So
what was p of zed given pi? Well, it was
actually happened to be a product over
the individual Z N and then a product
over the individual component K
times pi K raised to the Z N power. So
if you take the logarithm we get a sums
over N and K Z N K expectation of log of
pi K. So if you take the expectation
under this joint distribution over pi,
mu and sigma, well mu and sigma don't
matter and pi is the only thing left and
it goes through all of these sums and
pass the z and k and we're just left
with an expectation over log of pi k
under whatever our approximating
distribution on pi will be. We don't
know yet what that is but we know now
that we will need to be able to compute
the expectation of log of pi k under
this distribution.
Okay, let's go to the second term. We
just write this in. We leave it like
this. Now we go to the second term. Here
we need to take an expected value over
the log of P of X given zed. So on the
previous slide we had log of um we we
can actually if we know zed right then
we know what um log of P of X given zed
is it's just a product over the
individual XN
over gausian
and specifically the Gausians from
cluster K indicated to us by Z NK. So we
go through a sum over both the n and the
k if you want to do if if we did the
logarithm right all the products turn
into sums over zn and k which is dv
these 01 variables that tell us which
cluster we are in that basically select
an entry in the sum over k and then once
we know which term in the sum we select
we just have to write down the logarithm
of a gausian so that's this right so
gausians are exponentials of minus
quadratic form with mean and coariance
and then there's the log normalization
constant that contains sigma. We need to
take the expected value over this thing
under whatever our approximation to mu
and sigma is which we don't know yet but
that's the thing we're going to need to
compute. So now we know something we'll
need to compute. And now we do this cool
variational inference trick and we say
what do we actually have here as a
function of zed. So this expression here
this will be up to normalization our the
logarithm of our distribution over zed
and that's a function of zed
and we see that what this is is well
it's a sum over n and k ah that means
that there'll be a product in our
probability distribution
that factorizes over all the n over all
the individual samples
and then there'll be a z and k and some
number that we have to compute. We don't
know yet what that number is. It'll have
to be computed from our approximations
over q of distributions over pi and mu
and sigma. But it's just going to be a
number. So let's give a name to that
number. Let's call it the logarithm of
row n k. This is a real number. So it's
it can be the logarithm of something.
And then what we have here if we stare
at this and like squint with our eyes
and think of it as a distribution over
zed we can take the exponential and we
see that what we're going to get is a
product over every individual datim all
the individual components all the
clusters times the number row n k raised
to the z nk power that's not just a
discrete distribution it's a factorizing
discrete distribution. ution. Uhhuh. And
now the only thing we need to make sure
is that it's actually a distribution. So
that means the uh row and k have to sum
to one. And we can do that through the
normalization. So by using this constant
over here which we haven't used for
anything yet to normalize and set um the
actual parameters of our distribution to
RNK the responsibilities from EM um
which are now going to be equal to row
NK divided by well the row sum of the
row NK so that's a um dot sum Xis equal
to one basically operation which will
then give us our approximating
distribution to Q of to the posterior
over zed which is a fully factorizing
probability distribution condition
categorical distribution over the set n
case parameterized by a responsibility
rnk.
Uhhuh. And so what's just happened is we
now have we started by assuming that
there will be some approximating
distribution over the zs and one over
the pi mu and sigma. And what the
variational inference has automatically
told us without us imposing it is that
this bound actually factorizes over all
the n. Every single x will have its own
z n which contains k entries and they'll
be independent of each other. We didn't
impose this. It just happened because
our model has this structure. This is
called induced factorization.
And now we know that this is a discrete
distribution. Discrete distributions are
very easy to work with. the expected
value of a discrete distribution um
sorry the expected value of the random
variable under a discrete distribution
is just the probability
for um
um like that's in the vector right so
the the the the expected frequency with
which you're going to see cluster K in
for sample number n is literally just R
and K. It's just a property of discrete
distributions. Here it is. Yeah. So
discrete distributions and now you could
go to like the Wikipedia page for
discrete distributions and redouble the
properties of it like it's for example
that it's uh it's independent between
the Z KJ and it's like you know like
everything is here okay no so but so far
we only have an approximating
distribution for the zed we need to to
be able to build it we need an
approximating distribution to pi mu and
sigma and then we need to be able to
evaluate this thing this expected value
of this nasty expression here. So now
let's turn our attention to pi, mu and
sigma and see if we can find
approximations for these. Now knowing
that our approximating distribution to
zed will be a fully factorizing
categorical distribution. Okay, here I'm
going to spare you some nasty
computation and I'm just going to say uh
I'm going to motivate some convenient
notation. I'm going to introduce this NK
which are the uh column-wise sums over
the responsibilities. It's basically the
expected number of data points in each
cluster. We're going to compute a
weighted um cluster mean
weighted by the responsibilities and
we're also going to wait uh compute a
weighted uh empirical estimate for the
co-variance in each cluster. We call
that SK. That's just a convenient thing
to have because otherwise the slide
becomes so big that I can't use normal
font size anymore. And now we're going
to try and identify our approximating
distribution over pi mu and sigma.
Remember that we made no further
assumptions. We just said that there
will have to be a distribution over pi
mu and sigma. And the way to find that
is to find its logarithm actually is to
compute the expected value of the log
joint under q of z. That's our operator
operational rule for variational
inference. So let's go back again a few
slides. Write down the log joint that
here it is. I've just copied it over
from the previous slide. We will need to
take the expected value of this thing
under our categorical distribution over
Q over Z. And now we look at this
expression and again we see okay this
stuff does not depend on zed. There's no
zed in here. So that's going to be easy.
We just leave it in there. It's still I
mean it's a function of pi and mu and k.
So we'll need to keep it under
consideration. We can't call it a
constant, but it we don't need to take
an expected value over Q of zed because
zed doesn't show up here. Zed shows up
here and there. Okay, so let's move the
expected value inside. And then
rearrange some terms because we also see
that there is a pi here, but not here
and not there, but there's a pi here.
And there's a mu and sigma here, but not
here and there. And it's a mu and sigma
here and there. So let's rearrange the
sum. Let's move this bit over here and
this bit over there. Now um we are going
to have a term that uh of a bunch of
terms that depend on pi. One of them
only depends on pi or nothing else. One
of them contains zed. So we need to take
an expected value over zed under our
approximation which we now have. And um
then there will be all the terms that
depend on mu. One of them does not
depend on mu. One of them depends h
sorry one of them does not depend on
zed. One of them depends on zed. So
let's take an expected value over zed in
here. But zed actually just shows up
linearly in here in a big sum. So to
take the expected value
over zed n k we literally just replace
this with our r n k our responsibility
that we just found. Okay. So I've just
copied over this expression from above.
It's the exact same line again. And what
we now see is if you stare at this, oh
hang on, there is something over pi
here.
And then there is a sum
over the k
and a sum over the n. So this sum over k
shows up both here and there. So we can
take it outside. And in this q of zed we
um have a factorization over the n's and
the k's as well. So let's plug all those
sums in. We do that here and we will get
I've just rearranged terms. Um ah let's
actually do one one after the other.
Let's first do the one on on pi. So just
this term. Let's just think about it.
This will be a distribution over pi.
Let's write this down. So it's the
logarithm of a distribution over pi. So
it's this term which just contains pi.
That's our derivishly prior log derishly
prior plus an expected value of under
the discrete distribution over log of p
of zed given pi. Remember what is what
is p of z given pi? Well, it's just so
the expected sorry the logarithm of p of
z n k given pi is just
pi k.
So
let's plug that in here. And now we have
from the division prior this term from
this term this we can rearrange the
sums. You can put there's a sum over k
everywhere. So we put it outside.
There's one sum over n left. The log of
pi shows up everywhere here and there.
So we can take it outside of brackets.
And our log probability distribution
over pi will be a sum over k times a
bunch of numbers.
So that means as a distribution over pi
we have to take the exponential. Our
distribution will be a product over k pi
k raised to some power up to
normalization.
And if you've paid attention in the
exponential family class lecture, you
know that that's a dish distribution.
That's the exact form of a dishly
distribution. So without imposing
without requiring that our algorithm
uses a dishly distribution, it turns out
that our approximating distribution is a
dishly
just from the factorization assumption.
The algebraic form emerges naturally.
We'll have a dishly approximation which
contains parameters. So numbers alpha
plus nk where the nk are these things.
It's just a name for them. And then um
so you can look up on a and on on
Wikipedia everything about dishlay
distributions. And on Wikipedia you can
click on this link. You'll find lots and
lots of properties of dish
distributions. In particular you'll find
something about the expected value under
a dish distribution of the log of pi.
Ooh. And this is where the loop closes.
So now we go back and we said well our
problem was initially that to compute
our expected value over uh so to compute
our probability distribution over zed
our approximation we would need an
expected value over the logarithm of pi
of k under whatever our approximation on
pi was. So now we know our approximation
on pi will be a dishly and we know that
this expected value of the logarithm of
pi under a dishly distribution happens
to be something we can implement. It's
that this nasty thing it's called the
digama function and it's available in
scypi.
So now we're going to be able to write a
piece of code that does the first step
to compute the approximation to q of
zed. And now we need to do actually the
same thing again for mu and k. So um I
can I'm not going to do this in detail
here. Right. So so this is the stuff
from before. I've copied it over again.
Let's now focus on mu and sigma. Let's
uh find the uh approximating
distributions for mu and sigma. Well
here we discover that there's a lot of
sums in here again. So everything will
become a product over independent
distributions over mu k and sigma k.
Um so there's again an induced
factorization
and if you do the math it's the
basically the same as for the durish
layer you find that there is a very
tedious long expression that tells us
what the functional form of the
distribution over mu and sigma will be
it'll factoriize over the clusters k and
the individual distributions happen to
be a product of gaus inverse gamas so
they are again conjugate posteriors
And um
then you can go on Wikipedia
and look up see if you can find. So we
the last slide um whether we can find
the terms we need to compute to close
the loop to close the inductive loop and
find the terms that we need to construct
our approximation on um on the postivio
on zed Q of zed. So for that from the
previous slide we know that we need to
be able to compute expected values of
log of pi
we know how to do that it's a division
and then we also need expected values of
the log inverse determinant and a
quadratic form and for that back in the
day people had big books for example uh
there were a set of cool notes by
someone called um Sam Roise that you
could read up on and they would or you
would go to your colleague down the
corridor Tom Minka somewhere like that
and they would tell you a yeah I know
how to compute the expected values of of
quadratic forms other gausian there are
these nasty expressions here they are
you can read them up and you just
implement them and they involve
functions that we know they're nasty
expressions but they're are functions
that we know so that means we can now do
an an update loop that initially
initializes a dish distribution for pi a
gauss inverse gamma distribution for mu
and sigma and the categorical
distribution for zs at random and then
it iteratively updates it first
construct constructs a new approximation
for the the categorical distribution
over Z and K by computing these numbers
which are proportional to this object
which we have read off from the previous
slides and these are all variables that
we have because they are parameters of
the approximating distributions for mu
and sigma and pi for the dishlay and the
gaus inverse gamma and then we do the
update that was on the previous uh slide
to compute conjugate posteriors given
the Zs for pi, mu and sigma. And we keep
doing that. If you want to see how this
works, I actually uploaded the code to
beat up Jupyter Hub, but it's so long
that it doesn't make sense to look at it
now. Instead, I can show you as a
concluding picture how this algorithm
works. So, here is our data set again in
black dots. This this this um
old faithful geyser data set. And now
what I've done is I've initialized this
um uh what's called a dirishly gausian
mixture model or beijian gausian mixture
model with nine clusters. I said like
just try just try nine clusters see what
happens. And this is the initial cut um
um configuration. That's the expected
load of each cluster. It's it's a draw
from a sparse division distribution. And
that those are the initial clusters in
alpha value. to see their weight. So the
the faint ones are very low weight and
the thick ones are high weight and there
are nine of them. They're randomly
initialized. And now I just run the
algorithm for a few steps. Every time it
updates the Rs and then the parameters
of the conjugate posterior to mu, pi,
and sigma. And what this algorithm is
now able to do is it discovers that
they're just two two clusters. And back
then people were super excited about
this. You learn how many things there
are in a data set. Isn't that a cool
thing to be able to do? You show the
algorithm the world and it decides for
itself how many things are in your
world. In this case, two
and also what they look like. Namely,
they have two means and two coariances.
So this framework called variational
inference is maybe it's a more civilized
weapon from a more civilized age. It
used to be the way that we built machine
learning. But this is actually what we
used to call machine learning building
these very structured algorithms that
required a lot of tedious work. In this
case you first write down your
generative model P. Then you decide
you're going to do variational
inference. You impose a factorization on
your approximating distribution over the
unknown latent variables. In the case of
the Gaussian mixture model, there are z,
pi, mu and sigma. Then you go through
the math on a piece of paper like or on
a slide and derive the variational
updates. You find some forms that you
know how to compute. It involves lots of
known integrals, lots of oh, I know how
to compute this under an exponential
family. I know how to compute that
expected value. And then you implement
this as an iterative loop that keeps
doing these updates. In fact, you can
even evaluate the variational bound the
elbow as you do and you can see the
elbow rise. That used to be a buck
fixing test. People did just compute the
elbow and see whether it actually goes
up. If it doesn't go up, you made a
mistake. And um that used to be very
tedious. It was very difficult to
automate. It required a lot of knowledge
and finger
uh you know feeling um and elbow grease
like the ability to sit down for several
hours to do the derivations and it's
what used to be a lot of work for uh
research engineers and on Tuesday I will
briefly tell you the thing that you
probably expect that people stopped
doing this because auto was cool and
then we'll rediscover that in doing that
we actually lost a bit of structure
namely this induced factorization idea
that really helps with algorithmic
efficiency and it was rediscovered
arguably that's going to be my point by
some people in the deep learning
community in the form of attention. So
with that um please give feedback. I'll
see you on Tuesday for the last content
lecture of the of term and then
afterwards there'll just be a summary.
Thank you very much.