Machine Learning

Computing with Latent Variables

Rajesh Ranganath

1
Formal definition

Probabilistic latent variable models:

■ Data: x

■ Hidden Structure (latent variables): z

■ Model: p(x, z) = p(z)p(x | z) - prior p(z) and likelihood p(x | z)

■ Posterior: p(z | x) - probability of the hidden structure

2
 DSTN @ CMU . EDU
al Laboratory DJSCHLEGEL @ LBL . GOV
ratoryAstrophysics PRABHAT @ LBL . GOV

del of op-
h a varia-
xel inten-
able, with
properties
erties are
r distribu-
data sets.
mages. We
y survey,
he current
stial bod-

Figure 1. An image from the Sloan Digital Sky Survey (SDSS,

nerative model 2015) of a galaxy from the constellation Serpens, 100 million
[Regier+ 2015]
light years from Earth, along with several other galaxies and many
odel to be em-
The work we stars from our own galaxy.
3
dams, Harvard University RPA @ SEAS . HARVARD . EDU
offman, Adobe Research MDHOFFMA @ CS . PRINCETON . EDU
Lang, Carnegie Mellon University DSTN @ CMU . EDU

Astrophysics
chlegel, Lawrence Berkeley National Laboratory
, Lawrence Berkeley National Laboratory
DJSCHLEGEL @ LBL . GOV
PRABHAT @ LBL . GOV

Abstract
present a new, fully generative model of op-
l telescope image sets, along with a varia-
nal procedure for inference. Each pixel inten-
is treated as a Poisson random variable, with
ate parameter dependent on latent properties
stars and galaxies. Key latent properties are
mselves random, with scientific prior distribu-
ns constructed from large ancillary data sets.
check our approach on synthetic images. We
o run it on images from a major sky survey,
ere it exceeds the performance of the current
e-of-the-art method for locating celestial bod-
and measuring their colors.

oduction
Figure 1. An image from the Sloan Digital Sky Survey (SDSS,
Why a latent variable
2015)model?
of a galaxy from the constellation Serpens, 100 million
er presents Celeste, a new, fully generative model
omical image sets—the first such model to be em- light years from Earth, along with several other galaxies and many
investigated, to our knowledge. The work we stars from our own galaxy.
an encouraging example of principled statistical
e applied successfully to a science domain under-
y the machine learning community. It is unfortu- from a particular celestial body or from background at-
astronomy and cosmology receive comparatively mospheric noise—that pass through a telescope’s lens dur-
our attention: the scientific questions are funda- ing an exposure. Multiple celestial bodies may contribute
there are petabytes of data available, and we as a photons to a single image (e.g. Figure 1), and even to a
lysis community have a lot to offer the domain sci- single pixel of an image. Locating and characterizing the
One goal in reporting this work is to raise the profile imaged celestial bodies is an inference problem central to
problems for the machine-learning audience and astronomy. To date, the algorithms proposed for this in-
at much interesting research remains to be done. ference problem have been primarily heuristic, based on
finding bright regions in the images (Lupton et al., 2001;
w to the science. Stars and galaxies radiate photons. Stoughton et al., 2002).
nomical image records photons—each originating
Generative models are well-suited to this problem—for
ngs of the 32nd International Conference on Machine three reasons. First, to a good approximation, photon
Lille, France, 2015. JMLR: W&CP volume 37. Copy- counts from celestial objects are independent Poisson pro-
5 by the author(s).
cesses: each star or galaxy has an intrinsic brightness that
4
dams, Harvard University RPA @ SEAS . HARVARD . EDU
offman, Adobe Research MDHOFFMA @ CS . PRINCETON . EDU
Lang, Carnegie Mellon University DSTN @ CMU . EDU

Astrophysics
chlegel, Lawrence Berkeley National Laboratory
, Lawrence Berkeley National Laboratory
DJSCHLEGEL @ LBL . GOV
PRABHAT @ LBL . GOV

Abstract
present a new, fully generative model of op-
l telescope image sets, along with a varia-
nal procedure for inference. Each pixel inten-
is treated as a Poisson random variable, with
ate parameter dependent on latent properties
stars and galaxies. Key latent properties are
mselves random, with scientific prior distribu-
ns constructed from large ancillary data sets.
check our approach on synthetic images. We
o run it on images from a major sky survey,
ere it exceeds the performance of the current
e-of-the-art method for locating celestial bod-
and measuring their colors.

oduction
Figure 1. An image from the Sloan Digital Sky Survey (SDSS,
Prior knowledge about
2015) of a galaxy from the constellation Serpens, 100 million
er presents Celeste, a new, fully generative model
omical image sets—the first such model to be em- light years from Earth, along with several other galaxies and many
investigated, to our knowledge. The work we stars from our own galaxy.

p(x : light-measurement | zi = galaxy) = Known from physics

an encouraging example of principled statistical
i
e applied successfully to a science domain under-
y the machine learning community. It is unfortu- from a particular celestial body or from background at-
astronomy and cosmology receive comparatively mospheric noise—that pass through a telescope’s lens dur-
our attention: the scientific questions are funda- ing an exposure. Multiple celestial bodies may contribute
there are petabytes of data available, and we as a photons to a single image (e.g. Figure 1), and even to a
lysis community have a lot to offer the domain sci- single pixel of an image. Locating and characterizing the
One goal in reporting this work is to raise the profile imaged celestial bodies is an inference problem central to
problems for the machine-learning audience and astronomy. To date, the algorithms proposed for this in-
at much interesting research remains to be done. ference problem have been primarily heuristic, based on
finding bright regions in the images (Lupton et al., 2001;
w to the science. Stars and galaxies radiate photons. Stoughton et al., 2002).
nomical image records photons—each originating
Generative models are well-suited to this problem—for
ngs of the 32nd International Conference on Machine three reasons. First, to a good approximation, photon
Lille, France, 2015. JMLR: W&CP volume 37. Copy- counts from celestial objects are independent Poisson pro-
5 by the author(s).
cesses: each star or galaxy has an intrinsic brightness that
4
dams, Harvard University RPA @ SEAS . HARVARD . EDU
offman, Adobe Research MDHOFFMA @ CS . PRINCETON . EDU
Lang, Carnegie Mellon University DSTN @ CMU . EDU

Astrophysics
chlegel, Lawrence Berkeley National Laboratory
, Lawrence Berkeley National Laboratory
DJSCHLEGEL @ LBL . GOV
PRABHAT @ LBL . GOV

Abstract
present a new, fully generative model of op-
l telescope image sets, along with a varia-
nal procedure for inference. Each pixel inten-
is treated as a Poisson random variable, with
ate parameter dependent on latent properties
stars and galaxies. Key latent properties are
mselves random, with scientific prior distribu-
ns constructed from large ancillary data sets.
check our approach on synthetic images. We
o run it on images from a major sky survey,
ere it exceeds the performance of the current
e-of-the-art method for locating celestial bod-
and measuring their colors.

oduction
Figure 1. An image from the Sloan Digital Sky Survey (SDSS,
Prior knowledge about
2015) of a galaxy from the constellation Serpens, 100 million
er presents Celeste, a new, fully generative model
omical image sets—the first such model to be em- light years from Earth, along with several other galaxies and many
investigated, to our knowledge. The work we stars from our own galaxy.

p(x : light-measurement | zi = galaxy) = Known from physics

an encouraging example of principled statistical
i
e applied successfully to a science domain under-
y the machine learning community. It is unfortu- from a particular celestial body or from background at-
astronomy and cosmology receive comparatively mospheric noise—that pass through a telescope’s lens dur-
Posterior over mixtureing component “classifies”
an exposure. Multiple celestial
our attention: the scientific questions are funda- bodies may contribute
photons to a single image (e.g. Figure 1), and even to a
there are petabytes of data available, and we as a
lysis community have a lot to offer the domain sci- single pixel of an image. Locating and characterizing the
imaged celestial bodies is an inference problem central to
p(z = galaxy | x )
One goal in reporting this work is to raise the profile
i
problems for the machine-learning audience and i astronomy. To date, the algorithms proposed for this in-
at much interesting research remains to be done. ference problem have been primarily heuristic, based on
p(z = galaxy)p(x : light-measurement | zi = galaxy)
w to the science. Stars and galaxies radiatei photons.
finding bright regions in the images (Lupton et al., 2001;
Stoughton et al.,i 2002).
=P
nomical image records photons—each originating
ngs of the 32nd type∈{galaxy,
International Conference on Machine planet,...}
three
p(x : “light′′ |problem—for
Generative models are well-suited to this
zi =
reasons. First, toi a good approximation,
type)p(zi = type)
photon
Lille, France, 2015. JMLR: W&CP volume 37. Copy- counts from celestial objects are independent Poisson pro-
5 by the author(s).
cesses: each star or galaxy has an intrinsic brightness that
4
Uses of Latent Variable Models

■ Encoding prior knowledge

■ Combining simple distributions to create a more complex one
■ Uncovering hidden structure

5
Combining Simple Distributions

Take a categorical distribution

Categorical(1...K)

Take a normal distribution

Normal(µ, σ)

Combine

zi ∼ Categorical(1...K)
xi ∼ Normal(µzi , σzi )

Get a mixture of Gaussians, which is far more flexible

6
Q + A from Last Year
■ Encoding prior knowledge
■ Combining simple distributions to create a more complex one
■ Uncovering hidden structure

7
Q + A from Last Year
■ Encoding prior knowledge
■ Combining simple distributions to create a more complex one
■ Uncovering hidden structure
Do Latent Variables Help Predict x

What does it mean to predict x?

■ In some sense yes

p(x) = p(x(1) )p(x(2) | x(1) )p(x(3) |x(1) , x(2) )...

Predictions

7
Q + A from Last Year
■ Encoding prior knowledge
■ Combining simple distributions to create a more complex one
■ Uncovering hidden structure
How do I know my latent variable is correct? Can I used cross-
validation like when we predict Y?

This cannot be checked without assumptions because

■ The data generating distribution F(x) is all you get from the data

■ F(x) can be modeled without latent variables (though maybe slow)

■ For a fixed class, predictions and predictive checks help

7
Q + A from Last Year
■ Encoding prior knowledge
■ Combining simple distributions to create a more complex one
■ Uncovering hidden structure
How can we create graphs?

■ Based on prior knowledge

■ Based on computational considerations

■ Based on the hidden structure useful for a problem

7
Q + A from Last Year
■ Encoding prior knowledge
■ Combining simple distributions to create a more complex one
■ Uncovering hidden structure
Where does one use latent variables? Can we have some more
real-world examples?

■ We’ll have one more in a second

7
Q + A from Last Year
■ Encoding prior knowledge
■ Combining simple distributions to create a more complex one
■ Uncovering hidden structure
Are latent variables about the noise and getting an accurate data
generating distribution?

Yes/No
■ Noise variables could be latent variables

■ Noise is part of the unpredictability of x. How can the first

dimension of x be predicted?

p(x) = p(x(1) )p(x(2) | x(1) )p(x(3) |x(1) , x(2) )...

■ Latent variables can exist in concert with noise variables (Gaussian

noise)

■ Latent variables could also be structure that’s unobserved

7
Uncovering hidden structure: Finding Topics in
Documents

Data
■ M number of documents
■ V number of words
■ W: M × V matrix of words

Hidden Structure
■ A group of topics that describe the documents
■ Each document contains a distribution over topics

8
Data
■ M number of documents
■ V number of words
■ W: A M × V matrix of words
How do we describe the topics?

9
Data
■ M number of documents
■ V number of words
■ W: A M × V matrix of words
How do we describe the topics?
A distribution over the words called β k

9
■ M number of documents
■ V number of words
■ W: A M × V matrix of words
s How do we describe the document’s topic composition?
A distribution over the topics called θ i

10
Have two distributions
■ β k distributions over words for each topic
■ θ i distribution over topics for each document
Priors?

11
Have two distributions
■ β k distributions over words for each topic
■ θ i distribution over topics for each document
Priors? Dirichlet distribution

11
Still need a likelihood for data
■ M number of documents
■ V number of words
■ W: A M × V matrix of words

with hidden structure

■ β k distributions over words for each topic
■ θ i distribution over topics for each document
Assume all documents have same length?

12
For word m in document l
1. Draw word’s topic from zm,l ∼ Categorical(θ i )
2. Draw topic from for wm,l ∼ Categorical(β zm,l )

13
Topic Model
For each topic:
1. Draw distribution over words from Dirichlet(α)

For each document:

1. Draw distribution over topics from Dirichlet(κ)

For word m in document l:

1. Draw word’s topic from zm,l ∼ Categorical(θ i )
2. Draw word from topic for wm,l ∼ Categorical(β zm,l )

14
Topic Model
For each topic:
1. Draw distribution over words from Dirichlet(α)

For each document:

1. Draw distribution over topics from Dirichlet(κ)

For word m in document l:

1. Draw word’s topic from zm,l ∼ Categorical(θ i )
2. Draw word from topic for wm,l ∼ Categorical(β zm,l )
Compute posterior

p(β, θ , z, w)
p(β, θ , z | w) =
p(w)

14
 L ATENT D IRICHLET A LLOCATION

α θ z w N
M

aphical model representation of LDA. The boxes are “plates” representing

e outer plate represents documents, while the inner plate represents the repe
opics and words within a document.

15
 Figure 1: Posterior topics from the hierarchical Dirichlet process topic model on two large data sets.
[Hoffman+ 2013]These posteriors were approximated using stochastic variational inference with 1.8M ar- 16
Computing the posterior. Is it easy?

17
Computing the posterior is hard

The posterior distribution

p(x, z)
p(z | x) =
p(x)

The model p(x, z) is given; the challenge lies in computing p(x)

Z
p(x) = p(x, z)dz

This is a high dimensional integral (sum) which in general is

(analytically) intractable

18
Computing the posterior is hard
Bayesian Mixture of Gaussians

µk ∼ Normal(0, 1)
zi ∼ Categorical(1...K)
xi ∼ Normal(µzi , 1)

Marginal likelihood is
K
Z Y N X
Y K
p(x) = p(µk ) p(zi = j)p(xi | µj )dµ1 ...dµk
k=1 i=1 j=1

Swapping
K Z
XY N
Y
p(x) = p(µk ) p(zi = k)p(xi | µk )dµk
z k=1 i:z[i]=k

19
Show the equality

K
Z Y N X
Y K
p(x) = p(µk ) p(zi = j)p(xi | µj )dµ1 ...dµk
k=1 i=1 j=1
K
Z Y N
XY
= p(µk ) p(zi = z[i])p(xi | µz[i] )dµ1 ...dµk
k=1 z i=1
K
Z XY YN
= p(µk ) p(zi = z[i])p(xi | µz[i] )dµ1 ...dµk
z k=1 i=1
XZ YK YN
= p(µk ) p(zi = z[i])p(xi | µz[i] )dµ1 ...dµk
z k=1 i=1
K Z
XY YN
= p(µk ) p(zi = k)p(xi | µk )dµk
z k=1 i:z[i]=k
R R R
Uses a,b
f (a)g(b) = a
f (a) b
g(b)

20
Variational Inference

p.z j x/
⇤
[Link] / jj p.z j x//
[Link] / ⇤

init

Posit a family of distributions q(z; λ) indexed with parameter λ

21
Variational Inference

p.z j x/
⇤
[Link] / jj p.z j x//
[Link] / ⇤

init

Find λ such that q is close to p(z | x)

21
Variational Inference

p.z j x/
⇤
[Link] / jj p.z j x//
[Link] / ⇤

init

Closeness measured by the KL divergence

21
 KL(q(z; λ)||p(z | x)) = Eq [log q(z; λ) − log p(z | x)]
= Eq [log q(z; λ) − log p(z, x) + log p(x)]
= Eq [log q(z; λ) − log p(z, x)] + log p(x)

Equivalently,

log p(x) = Eq [log p(z, x) − log q(z; λ)] + KL(q(z; λ)||p(z | x))
≥ Eq [log p(z, x) − log q(z; λ)]

A lower bound on the evidence log p(x)

22
The Evidence Lower Bound

L (λ) = Eq [log p(x, z)] − Eq [log q(z; λ)]

■ KL is intractable; VI optimizes the evidence lower bound (ELBO)

□ It is a lower bound on log p(x)
□ Maximizing the ELBO is equivalent to minimizing the KL

■ The ELBO trades off two terms

□ The first term prefers q(·) to place its mass on the MAP estimate
□ The second term encourages q(·) to be diffuse

■ Approximation q is chosen to match types

23
The Evidence Lower Bound

L (λ) = Eq [log p(x | z)] − KL(q(z; λ)||p(z))

■ KL is intractable; VI optimizes the evidence lower bound (ELBO)

□ It is a lower bound on log p(x)
□ Maximizing the ELBO is equivalent to minimizing the KL

■ The ELBO trades off two terms

□ The first term prefers q(·) to maximize the likelihood
□ The second term regularizes q(·) to the prior

■ Approximation q is chosen to match types

24
The Recipe

p(x, z) Z [Link] /

(· · · )q(z; )dz r
q(z; )

25
 L ATENT D IRICHLET A LLOCATION

α θ z w N
M

aphical model representation of LDA. The boxes are “plates” representing

e outer plate represents documents, while the inner plate represents the repe
opics and words within a document.

26
VI for LDA B LEI , N G , AND J ORDAN B LEI , N G , AND J ORDAN

A.1 Computing E[log(θi | α)] Finally, we expand Eq. (14) in terms of the model parameters (α, β) and the variational parameters
(γ, φ). Each of the five lines below expands one of the five terms in the bound:
The need to compute the expected value of the log of a single probability component under the
Dirichlet arises repeatedly in deriving the inference and parameter estimation procedures for LDA.
This value can be easily computed from the natural parameterization of the exponential family k k

representation of the Dirichlet distribution.

L (γ, φ; α, β) = log Γ ∑kj=1 α j ∑ log Γ(αi ) + ∑ (αi 1) Ψ(γi ) Ψ ∑kj=1 γ j
i=1 i=1
Recall that a distribution is in the exponential family if it can be written in the form: N k
T +∑ ∑ φni Ψ(γi ) Ψ ∑kj=1 γ j
p(x | η) = h(x) exp η T (x) A(η) , n=1 i=1
N k V
where η is the natural parameter, T (x) is the sufficient statistic, and A(η) is the log of the normal- +∑ ∑ ∑ φni wnj log βi j (15)
ization factor. n=1 i=1 j=1
We can write the Dirichlet in this form by exponentiating the log of Eq. (1): k k
log Γ ∑kj=1 γ j + ∑ log Γ(γi ) ∑ (γi 1) Ψ(γi ) Ψ ∑kj=1 γ j
p(θ | α) = exp ∑ki=1 (αi 1) log θi + log Γ ∑ki=1 αi ∑ki=1 log Γ(αi ) . i=1 i=1
N k
From this form, we immediately see that the natural parameter of the Dirichlet is ηi = αi 1 and ∑ ∑ φni log φni ,
the sufficient statistic is T (θi ) = log θi . Furthermore, using the general fact that the derivative of n=1 i=1
the log normalization factor with respect to the natural parameter is equal to the expectation of the
sufficient statistic, we obtain: where we have made use of Eq. (8).
L ATENT D IRICHLET A LLOCATION
In the following two sections, we show how to maximize this lower bound with respect to the
E[log θi | α] = Ψ(αi ) Ψ ∑kj=1 α j
variational parameters φ and γ.
where Ψ is the digamma function, the first derivative of the log Gamma function.
A.3.2 VARIATIONAL D IRICHLET A.3.1 VARIATIONAL MULTINOMIAL
A.2 Newton-Raphson methods for a Hessian with special structure
Next, we maximize Eq. (15) with respect to γi , the ith component of the posterior Dirichlet param- We first maximize Eq. (15) with respect to φni , the probability that the nth word is generated by
In this
eter. section
The terms we describeγiaare:
containing linear algorithm for the usually cubic Newton-Raphson optimization latent topic i. Observe that this is a constrained maximization since ∑ki=1 φni = 1.
method. This method is used for maximum likelihood estimation of the Dirichlet distribution (Ron-
We form the Lagrangian by isolating the terms which contain φni and adding the appropriate
ning, 1989, Minka, 2000).
k N
Lagrange multipliers. Let βiv be p(wvn = 1 | zi = 1) for the appropriate v. (Recall that each wn is
L[γ] = ∑ (αi optimization
The Newton-Raphson Ψ ∑kj=1 γ finds
1) Ψ(γi ) technique j +a∑stationary ) Ψof∑akj=1
φni Ψ(γipoint γj
function by iterating:
i=1 n=1 a vector of size V with exactly one component equal to one; we can select the unique v such that
αnew = αold H(αkold ) 1 g(αold ) wvn = 1):
log Γ ∑kj=1 γ j + log Γ(γi ) ∑ (γi 1) Ψ(γi ) Ψ ∑kj=1 γ j .
i=1
where H(α) and g(α) are the Hessian matrix and gradient respectively at the point α. In general,
3 L[φni ] = φni Ψ(γi ) Ψ ∑kj=1 γ j + φni log βiv φni log φni + λn ∑kj=1 φni 1 ,
this algorithm scales as O(N ) due to the matrix inversion.
This simplifies to:
If the Hessian matrix is of the form:
k where we have dropped the arguments of L for simplicity, and where the subscript φni denotes that
T
L[γ] = ∑ Ψ(γi ) Ψ ∑kj=1 γ j H =
αi diag(h) φni1z1γi,
+ ∑Nn=1 + log Γ ∑kj=1 γ j + log Γ(γi ). (10) we have retained only those terms in L that are a function of φni . Taking derivatives with respect to
i=1 φni , we obtain:
where diag(h) is defined to be a diagonal matrix with the elements of the vector h along the diagonal,
thentake
We we the
canderivative
apply the with
matrix inversion
respect to γi :lemma and obtain:
∂L
diag(h) 1 11T diag(h) 1 = Ψ(γi ) Ψ ∑kj=1 γ j + log βiv log φni 1 + λ.
∂L H 1 = diag(h) 1 k ∂φni
= Ψ0 (γi ) αi + ∑Nn=1 φni γi Ψ0 ∑ k1
z j=1 ∑h jα1 j + ∑Nn=1 φn j
+γ∑j kj=1 γj .
∂γi j=1
Multiplying by the gradient, we obtain the ith component: Setting this derivative to zero yields the maximizing value of the variational parameter φni (cf. Eq. 6):
Setting this equation to zero yields a maximum at:
gi c
(H 1 g)i = φni ∝ βiv exp Ψ(γi ) Ψ ∑kj=1 γ j . (16)
γi = αi + ∑Nn=1hφi ni . (17)

Since Eq. (17) depends on the variational 1018

multinomial φ, full variational inference requires 1020
alternating between Eqs. (16) and (17) until the bound converges.

A.4 Parameter estimation

In this final section, we consider the problem of obtaining empirical Bayes estimates of the model 27
The Variational Inference Recipe
Start with a model:

p(z, x)

28
The Variational Inference Recipe
Choose a variational approximation:

q(z; λ)

28
The Variational Inference Recipe
Write down the ELBO:

L (λ) = Eq(z;λ) [log p(x, z) − log q(z; λ)]

28
The Variational Inference Recipe
Compute the expectation(integral):

Example: L (λ) = xλ2 + log λ

28
The Variational Inference Recipe
Take derivatives:

1
Example: ∇λ L (λ) = 2xλ +
λ

28
The Variational Inference Recipe
Optimize:

λt+1 = λt + ρt ∇λ L

28
The Variational Inference Recipe

p(x, z) Z [Link] /

(· · · )q(z; )dz r
q(z; )

28
Example: Bayesian Logistic Regression

■ Data pairs yi , xi
■ xi are covariates
■ yi are label
■ z is the regression coefficient
■ Generative process

p(z) ∼ N(0, 1)
p(yi | xi , z) ∼ Bernoulli(σ(zxi ))

29
VI for Bayesian Logistic Regression

Assume:
■ We have one data point (y, x)
■ The approximating family q is the normal; λ = (µ, σ2 )
The ELBO is

L (µ, σ2 ) = Eq [log p(z) + log p(y | x, z) − log q(z)]

30
VI for Bayesian Logistic Regression

L (µ, σ2 )
= Eq [log p(z) − log q(z) + log p(y | x, z)]

31
VI for Bayesian Logistic Regression

L (µ, σ2 )
= Eq [log p(z) − log q(z) + log p(y | x, z)]
1 1
= − (µ2 + σ2 ) + log σ2 + Eq [log p(y | x, z)] + C
2 2

31
VI for Bayesian Logistic Regression

L (µ, σ2 )
= Eq [log p(z) − log q(z) + log p(y | x, z)]
1 1
= − (µ2 + σ2 ) + log σ2 + Eq [log p(y | x, z)] + C
2 2
1 2 1
= − (µ + σ2 ) + log σ2 + Eq [yxz − log(1 + exp(xz))]
2 2

31
VI for Bayesian Logistic Regression

L (µ, σ2 )
= Eq [log p(z) − log q(z) + log p(y | x, z)]
1 1
= − (µ2 + σ2 ) + log σ2 + Eq [log p(y | x, z)] + C
2 2
1 2 1
= − (µ + σ2 ) + log σ2 + Eq [yxz − log(1 + exp(xz))]
2 2
1 2 1
= − (µ + σ ) + log σ2 + yxµ − Eq [log(1 + exp(xz))]
2
2 2

31
VI for Bayesian Logistic Regression

L (µ, σ2 )
= Eq [log p(z) − log q(z) + log p(y | x, z)]
1 1
= − (µ2 + σ2 ) + log σ2 + Eq [log p(y | x, z)] + C
2 2
1 2 1
= − (µ + σ2 ) + log σ2 + Eq [yxz − log(1 + exp(xz))]
2 2
1 2 1
= − (µ + σ ) + log σ2 + yxµ − Eq [log(1 + exp(xz))]
2
2 2
We are stuck.
1. We cannot analytically take that expectation.
2. The expectation hides the objectives dependence on the variational
parameters. This makes it hard to directly optimize.

31
Options?

■ Derive a model specific bound:

[Jordan and Jaakola; 1996], [Braun and McAuliffe; 2008], others

■ More general approximations that require model-specific analysis:

[Wang and Blei; 2013], [Knowles and Minka; 2011]

32
Nonconjugate Models

■ Discrete Choice Models

■ Nonlinear Time series Models
■ Bayesian Neural Networks
■ Deep Latent Gaussian Models
■ Deep Exponential Families
■ Models with Attention
(e.g. Sparse Gamma or Poisson)
(such as DRAW)
■ Correlated Topic Model
■ Generalized Linear Models
(including nonparametric
(Poisson Regression)
variants)
■ Stochastic Volatility Models
■ Sigmoid Belief Network

We need a solution that does not entail model specific work

33
 Black Box Variational Inference (BBVI)
Black box variational inference

REUSABLE
REUSABLE
REUSABLE MASSIVE
VARIATIONAL
VARIATIONAL
VARIATIONAL DATA
FAMILIES
FAMILIES
FAMILIES

ANY MODEL
BLACK BOX
p.ˇ; z j x/
VARIATIONAL
INFERENCE

I Sample from q. /
I Form noisy gradients without model-specific computation

I Use stochastic optimization 34

The Problem in the Classical VI Recipe

p(x, z) Z [Link] /

(· · · )q(z; )dz r
q(z; )

35
The New VI Recipe

p(x, z) Z [Link] /

r (· · · )q(z; )dz
q(z; )

Use stochastic optimization!

36
Computing Gradients of Expectations
■ Define

g(z, λ) = log p(x, z) − log q(z; λ)

■ What is ∇λ L
Z
∇λ L = ∇λ q(z; λ)g(z, λ)dz
Z
= ∇λ q(z; λ)g(z, λ) + q(z; λ)∇λ g(z, λ)dz
Z
= q(z; λ)∇λ log q(z; λ)g(z, λ) + q(z; λ)∇λ g(z, λ)dz

= Eq(z;λ) [∇λ log q(z; λ)g(z, λ) + ∇λ g(z, λ)]

∇λ q
Using ∇λ log q = q

37
Roadmap

■ Score Function Gradients

■ Reparameterization Gradients

38
Score Function Gradients of the ELBO

39
Score Function Estimator

L (λ) = Eq(z;λ) [log p(x, z) − log q(z; λ)]

Recall

∇λ L = Eq(z;λ) [∇λ log q(z; λ)g(z, λ) + ∇λ g(z, λ)]

Simplify:

Eq [∇λ g(z, λ)] = Eq [∇λ log q(z; λ)] = 0

Gives the gradient:

∇λ L = Eq(z;λ) [∇λ log q(z; λ)(log p(x, z) − log q(z; λ))]

Sometimes called likelihood ratio or REINFORCE gradients

[Glynn 1990; Williams, 1992; Wingate+ 2013; R+ 2014; Mnih+ 2014]
40
Noisy Unbiased Gradients

Gradient: Eq(z;λ) [∇λ log q(z; λ)(log p(x, z) − log q(z; λ))]

Noisy unbiased gradients with Monte Carlo!

S
1X
∇λ log q(zs ; λ)(log p(x, zs ) − log q(zs ; λ)),
S s=1
where zs ∼ q(z; λ)

41
Basic BBVI

Algorithm 1: Basic Black Box Variational Inference

Input : Model log p(x, z),

Variational approximation q(z; λ)
Output : Variational Parameters: λ

while not converged do

z[s] ∼ q // Draw S samples from q
ρ = t-th value
PSof a Robbins Monro sequence
λ = λ + ρ 1S s=1 ∇λ log q(z[s]; λ)(log p(x, z[s]) − log q(z[s]; λ))
t=t+1
end

42
The requirements for inference

The noisy gradient:

S
1X
∇λ log q(zs ; λ)(log p(x, zs ) − log q(zs ; λ)),
S s=1
where zs ∼ q(z; λ)

To compute the noisy gradient of the ELBO we need

■ Sampling from q(z)
■ Evaluating ∇λ log q(z; λ)
■ Evaluating log p(x, z) and log q(z)

There is no model specific work: black box criteria are

satisfied

43
 Black Box Variational Inference
Black box variational inference

REUSABLE
REUSABLE
REUSABLE MASSIVE
VARIATIONAL
VARIATIONAL
VARIATIONAL DATA
FAMILIES
FAMILIES
FAMILIES

ANY MODEL
BLACK BOX
p.ˇ; z j x/
VARIATIONAL
INFERENCE

I Sample from q. /
I Form noisy gradients without model-specific computation

I Use stochastic optimization 44

Discussion

Variational Inference for Bayesian Mixtures of Gaussians

45
Problem: Basic BBVI doesn’t work

Variance of the gradient can be a problem

Varq(z;ν) = Eq(z;ν) [(∇ν log q(z; ν)(log p(x, z) − log q(z; ν)) − ∇ν L )2 ].

2.0
PDF
1.5 Abs Mu Score
1.0

0.5

0.0
2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0
Intuition:
Sampling rare values can lead to large scores and thus high
variance

46
Solution: Control Variates
Replace with f with f̂ where E[f̂ (z)] = E[f (z)]. General such
class:
f̂ (z) ≜ f (z) − a(h(z) − E[h(z)])

6
PDF
5 f = x + x2
fˆ; h = x2
4
fˆ; h = f
3

−1
−2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0

■ h is a function of our choice

■ a is chosen to minimize the variance
■ Good h have high correlation with the original function f

47
Solution: Control Variates
Replace with f with f̂ where E[f̂ (z)] = E[f (z)]. General such
class:
f̂ (z) ≜ f (z) − a(h(z) − E[h(z)])

6
PDF
5 f = x + x2
fˆ; h = x2
4
fˆ; h = f
3

−1
−2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0

■ For variational inference we need functions with known q

expectation
■ Set h as ∇λ log q(z; λ)
■ Simple as Eq [∇λ log q(z; λ)] = 0 for any q
47
Solution: Control Variates
Replace with f with f̂ where E[f̂ (z)] = E[f (z)]. General such
class:
f̂ (z) ≜ f (z) − a(h(z) − E[h(z)])

6
PDF
5 f = x + x2
fˆ; h = x2
4
fˆ; h = f
3

−1
−2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0

Many of the other techniques from Monte Carlo can help:

■ Importance Sampling, Quasi Monte Carlo, Rao-Blackwellization

[Ruiz+ 2016; Ranganath+2014; Titsias+2015; Mnih+2016]

47
Nonconjugate Models

■ Discrete Choice Models

■ Nonlinear Time series Models
■ Bayesian Neural Networks
■ Deep Latent Gaussian Models
■ Deep Exponential Families
■ Models with Attention
(e.g. Sparse Gamma or Poisson)
(such as DRAW)
■ Correlated Topic Model
■ Generalized Linear Models
(including nonparametric
(Poisson Regression)
variants)
■ Stochastic Volatility Models
■ Sigmoid Belief Network

We can design models based on data rather than inference.

48
More Assumptions?

The current black box criteria

■ Sampling from q(z)
■ Evaluating ∇λ log q(z; λ)
■ Evaluating log p(x, z) and log q(z)
Can we make additional assumptions that are not too restrictive?

49
Pathwise Gradients of the ELBO

50
Reparameterization Estimator

Assume
1. z = t(ε, λ) for ε ∼ s(ε) implies z ∼ q(z; λ)
Example:

ε ∼ Normal(0, 1)
z = εσ + µ
→ z ∼ Normal(µ, σ2 )

2. log p(x, z) and log q(z) are differentiable with respect to z

51
Reparameterization Estimator
Recall

∇λ L = Eq(z;λ) [∇λ log q(z; λ)g(z, λ) + ∇λ g(z, λ)]

Rewrite using using z = t(ε, λ)

∇λ L = Es(ε) [∇λ log s(ε)g(t(ε, λ), λ) + ∇λ g(t(ε, λ), λ)]

To differentiate:

∇L (λ) = Es(ε) [∇λ g(t(ε, λ), λ)]

= Es(ε) [∇z [log p(x, z) − log q(z; λ)]∇λ t(ε, λ) − ∇λ log q(z; λ)]
= Es(ε) [∇z [log p(x, z) − log q(z; λ)]∇λ t(ε, λ)]

This is also known as the pathwise gradient.

[Glasserman 1991; Fu 2006; Kingma+ 2014; Rezende+ 2014; Titsias+ 2014]

52
Variance Comparison

[R+ 2018]

53
What’s an example problem that might have gradients where one
is better than the other?

54
Score Function Estimator vs. Reparameterization
Estimator

Pathwise
Score Function
■ Differentiates the function
■ Differentiates the density
∇z [log p(x, z) − log q(z; λ)]
∇λ q(z; λ)
■ Requires differentiable
■ Works for discrete and
models
continuous models
■ Requires variational
■ Works for large class of
approximation to have form
variational approximations
z = t(ε, λ)
■ Variance can be a big
■ Generally better behaved
problem
variance

55
How do we use both estimators at the same time?

56
Parameter Learning with Variational Inference

Train model pθ (x, z) by maximum likelihood

Z
log pθ (x) = log pθ (x, z)dz

Hard to compute the integral. If posterior was known,

pθ (x, z)
log pθ (x) = log
pθ (z | x)

Posterior is hard because of integration (unknown p(x)).

57
Parameter Learning with Variational Inference

Maximize lower bound on the likelihood

log pθ (x) = Eq [log p(z, x) − log q(z; λ)] + KL(q(z; λ)||p(z | x))
≥ Eq [log p(z, x) − log q(z; λ)] := L (λ, θ )

A lower bound on the evidence log pθ (x)

58
Parameter Learning with Variational Inference

Maximize lower bound on the likelihood

log pθ (x) = Eq [log p(z, x) − log q(z; λ)] + KL(q(z; λ)||p(z | x))
≥ Eq [log p(z, x) − log q(z; λ)] := L (λ, θ )

A lower bound on the evidence log pθ (x)

Optimize L (λ, θ ) using gradients
■ Use standard gradients for θ

■ Use score/reparameterization gradients for λ

58
Parameter Learning with Variational Inference

Maximize lower bound on the likelihood

log pθ (x) = Eq [log p(z, x) − log q(z; λ)] + KL(q(z; λ)||p(z | x))
≥ Eq [log p(z, x) − log q(z; λ)] := L (λ, θ )

A lower bound on the evidence log pθ (x)

Could instead maximize θ with λ fixed and vice-versa
■ λ = arg max L (λ, θ
t λ t−1 )

■ θ t = arg maxθ L (λt , θ )

Called coordinate ascent

58
Parameter Learning with Variational Inference

Maximize lower bound on the likelihood

log pθ (x) = Eq [log p(z, x) − log q(z; λ)] + KL(q(z; λ)||p(z | x))
≥ Eq [log p(z, x) − log q(z; λ)] := L (λ, θ )

A lower bound on the evidence log pθ (x)

Could instead maximize θ and choose optimal q = pθ (z | x)
■ Compute optimal q = p
t θ t−1 (z | x)

■ θ t = arg maxθ L (qt , θ )

Called the Expectation Maximization Algorithm

58