Part X
CS229 Lecture notes
Andrew Ng
Factor analysis
When we have data x(i) ∈ Rd that comes from a mixture of several Gaussians, the EM algorithm can be applied to fit a mixture model. In this setting, we usually imagine problems where we have sufficient data to be able to discern the multiple-Gaussian structure in the data. For instance, this would be the case if our training set size n was significantly larger than the dimension d of the data.
Now, consider a setting in which d ≫ n. In such a problem, it might be difficult to model the data even with a single Gaussian, much less a mixture of Gaussian. Specifically, since the n data points span only a low-dimensional subspace of Rd, if we model the data as Gaussian, and estimate the mean and covariance using the usual maximum likelihood estimators,
1 n μ=n x(i)
i=1 1 n
Σ = n (x(i) −μ)(x(i) −μ)T, i=1
we would find that the matrix Σ is singular. This means that Σ−1 does not exist, and 1/|Σ|1/2 = 1/0. But both of these terms are needed in computing the usual density of a multivariate Gaussian distribution. Another way of stating this difficulty is that maximum likelihood estimates of the parameters result in a Gaussian that places all of its probability in the affine space spanned by the data,1 and this corresponds to a singular covariance matrix.
n (i) n
This is the set of points x satisfying x = i=1 αix , for some αi’s so that i=1 α1 =
1.
1
1
2
More generally, unless n exceeds d by some reasonable amount, the max- imum likelihood estimates of the mean and covariance may be quite poor. Nonetheless, we would still like to be able to fit a reasonable Gaussian model to the data, and perhaps capture some interesting covariance structure in the data. How can we do this?
In the next section, we begin by reviewing two possible restrictions on Σ that allow us to fit Σ with small amounts of data but neither will give a satisfactory solution to our problem. We next discuss some properties of Gaussians that will be needed later; specifically, how to find marginal and conditonal distributions of Gaussians. Finally, we present the factor analysis model, and EM for it.
1 Restrictions of Σ
If we do not have sufficient data to fit a full covariance matrix, we may place some restrictions on the space of matrices Σ that we will consider. For instance, we may choose to fit a covariance matrix Σ that is diagonal. In this setting, the reader may easily verify that the maximum likelihood estimate of the covariance matrix is given by the diagonal matrix Σ satisfying
1 n
Σjj = (x(i) −μj)2.
Thus, Σjj is just the empirical estimate of the variance of the j-th coordinate of the data.
Recall that the contours of a Gaussian density are ellipses. A diagonal Σ corresponds to a Gaussian where the major axes of these ellipses are axis- aligned.
Sometimes, we may place a further restriction on the covariance matrix that not only must it be diagonal, but its diagonal entries must all be equal. In this setting, we have Σ = σ2I, where σ2 is the parameter under our control. The maximum likelihood estimate of σ2 can be found to be:
nj i=1
1dn
σ2 = (x(i) −μj)2.
This model corresponds to using Gaussians whose densities have contours that are circles (in 2 dimensions; or spheres/hyperspheres in higher dimen- sions).
nd j j=1 i=1
3
If we are fitting a full, unconstrained, covariance matrix Σ to data, it is necessary that n ≥ d + 1 in order for the maximum likelihood estimate of Σ not to be singular. Under either of the two restrictions above, we may obtain non-singular Σ when n ≥ 2.
However, restricting Σ to be diagonal also means modeling the different coordinates xi, xj of the data as being uncorrelated and independent. Often, it would be nice to be able to capture some interesting correlation structure in the data. If we were to use either of the restrictions on Σ described above, we would therefore fail to do so. In this set of notes, we will describe the factor analysis model, which uses more parameters than the diagonal Σ and captures some correlations in the data, but also without having to fit a full covariance matrix.
2 Marginals and conditionals of Gaussians
Before describing factor analysis, we digress to talk about how to find condi- tional and marginal distributions of random variables with a joint multivari- ate Gaussian distribution.
Suppose we have a vector-valued random variable x=x1 ,
x2
where x1 ∈ Rr, x2 ∈ Rs, and x ∈ Rr+s. Suppose x ∼ N(μ,Σ), where
μ=μ1, Σ=Σ11 Σ12. μ2 Σ21 Σ22
Here, μ1 ∈ Rr, μ2 ∈ Rs, Σ11 ∈ Rr×r, Σ12 ∈ Rr×s, and so on. Note that since covariance matrices are symmetric, Σ12 = ΣT21.
Under our assumptions, x1 and x2 are jointly multivariate Gaussian. What is the marginal distribution of x1? It is not hard to see that E[x1] = μ1, and that Cov(x1) = E[(x1 − μ1)(x1 − μ1)] = Σ11. To see that the latter is true, note that by definition of the joint covariance of x1 and x2, we have
that
Cov(x) = Σ
=Σ11 Σ12 Σ21 Σ22
= E[(x−μ)(x−μ)T]
= Ex1−μ1 x1−μ1 T
x2 − μ2 x2 − μ2
= E (x1 −μ1)(x1 −μ1)T (x1 −μ1)(x2 −μ2)T .
(x2 − μ2)(x1 − μ1)T (x2 − μ2)(x2 − μ2)T
Matching the upper-left subblocks in the matrices in the second and the last lines above gives the result.
Since marginal distributions of Gaussians are themselves Gaussian, we therefore have that the marginal distribution of x1 is given by x1 ∼ N (μ1, Σ11).
Also, we can ask, what is the conditional distribution of x1 given x2? By referring to the definition of the multivariate Gaussian distribution, it can be shown that x1|x2 ∼ N (μ1|2, Σ1|2), where
μ = μ +Σ Σ−1(x −μ), (1) 1|2 1122222
Σ = Σ −Σ Σ−1Σ . (2) 1|2 11122221
When we work with the factor analysis model in the next section, these formulas for finding conditional and marginal distributions of Gaussians will be very useful.
3 The Factor analysis model
In the factor analysis model, we posit a joint distribution on (x, z) as follows, where z ∈ Rk is a latent random variable:
z ∼ N(0,I)
x|z ∼ N(μ+Λz,Ψ).
Here, the parameters of our model are the vector μ ∈ Rd, the matrix Λ ∈ Rd×k, and the diagonal matrix Ψ ∈ Rd×d. The value of k is usually chosen to be smaller than d.
4
5
Thus, we imagine that each datapoint x(i) is generated by sampling a k dimension multivariate Gaussian z(i). Then, it is mapped to a d-dimensional affine space of Rd by computing μ + Λz(i). Lastly, x(i) is generated by adding covariance Ψ noise to μ + Λz(i).
Equivalently (convince yourself that this is the case), we can therefore also define the factor analysis model according to
z ∼ N(0,I) ǫ ∼ N(0,Ψ)
x = μ+Λz+ǫ
where ǫ and z are independent.
Let’s work out exactly what distribution our model defines. Our random
variables z and x have a joint Gaussian distribution xz ∼ N ( μ z x , Σ ) .
We will now find μzx and Σ.
We know that E[z] = 0, from the fact that z ∼ N(0,I). Also, we have
that
E[x] = E[μ+Λz+ǫ]
= μ+ΛE[z]+E[ǫ]
= μ. Putting these together, we obtain
μ z x = ⃗0 μ
Next, to find Σ, we need to calculate Σzz = E[(z − E[z])(z − E[z])T ] (the upper-left block of Σ), Σzx = E[(z − E[z])(x − E[x])T ] (upper-right block), and Σxx = E[(x − E[x])(x − E[x])T ] (lower-right block).
Now, since z ∼ N(0,I), we easily find that Σzz = Cov(z) = I. Also,
E[(z−E[z])(x−E[x])T] = E[z(μ+Λz+ǫ−μ)T] = E[zzT ]ΛT + E[zǫT ]
= ΛT.
In the last step, we used the fact that E[zzT ] = Cov(z) (since z has zero mean), and E[zǫT] = E[z]E[ǫT] = 0 (since z and ǫ are independent, and
6 hence the expectation of their product is the product of their expectations).
Similarly, we can find Σxx as follows:
E[(x−E[x])(x−E[x])T] = E[(μ+Λz+ǫ−μ)(μ+Λz+ǫ−μ)T] = E[ΛzzTΛT +ǫzTΛT +ΛzǫT +ǫǫT]
= ΛE[zzT ]ΛT + E[ǫǫT ] = ΛΛT +Ψ.
Putting everything together, we therefore have that
z ∼ N ⃗0 , I Λ T . ( 3 )
x μΛΛΛT+Ψ
Hence, we also see that the marginal distribution of x is given by x ∼ N(μ,ΛΛT +Ψ). Thus,givenatrainingset{x(i);i=1,…,n},wecanwrite down the log likelihood of the parameters:
l(μ,Λ,Ψ)=logn 1 exp−1(x(i) −μ)T(ΛΛT +Ψ)−1(x(i) −μ). i=1 (2π)d/2|ΛΛT +Ψ|1/2 2
To perform maximum likelihood estimation, we would like to maximize this quantity with respect to the parameters. But maximizing this formula ex- plicitly is hard (try it yourself), and we are aware of no algorithm that does so in closed-form. So, we will instead use to the EM algorithm. In the next section, we derive EM for factor analysis.
4 EM for factor analysis
The derivation for the E-step is easy. We need to compute Qi(z(i)) = p(z(i)|x(i); μ, Λ, Ψ). By substituting the distribution given in Equation (3) into the formulas (1-2) used for finding the conditional distribution of a Gaussian, we find that z(i)|x(i); μ, Λ, Ψ ∼ N (μz(i)|x(i) , Σz(i)|x(i) ), where
μz(i)|x(i) = ΛT (ΛΛT + Ψ)−1(x(i) − μ), Σz(i)|x(i) = I − ΛT (ΛΛT + Ψ)−1Λ.
So, using these definitions for μz(i)|x(i) and Σz(i)|x(i), we have
Q(z(i))= 1 exp−1(z(i) −μ (i) (i))TΣ−1 (z(i) −μ (i)
i (2π)k/2|Σz(i)|x(i) |1/2 2 z |x z(i)|x(i) z |x
(i)).
Let’s now work out the M-step. Here, we need to maximize
n Qi(z(i)) log p(x(i), z(i); μ, Λ, Ψ)dz(i) (4)
i=1 z(i) Qi(z(i))
with respect to the parameters μ, Λ, Ψ. We will work out only the optimiza- tion with respect to Λ, and leave the derivations of the updates for μ and Ψ as an exercise to the reader.
We can simplify Equation (4) as follows:
n i=1 z(i)
n
= Ez(i)∼Qi log p(x(i)|z(i); μ, Λ, Ψ) + log p(z(i)) − log Qi(z(i)) (6)
i=1
Here, the “z(i) ∼ Qi” subscript indicates that the expectation is with respect to z(i) drawn from Qi. In the subsequent development, we will omit this subscript when there is no risk of ambiguity. Dropping terms that do not depend on the parameters, we find that we need to maximize:
n
E log p(x(i)|z(i); μ, Λ, Ψ)
i=1 =
=
i=1
Let’s maximize this with respect to Λ. Only the last term above depends
on Λ. Taking derivatives, and using the facts that tra = a (for a ∈ R), trAB = trBA, and ∇AtrABAT C = CAB + CT ABT , we get:
n 1 (i) (i)T −1 (i) (i) ∇Λ−E 2(x −μ−Λz )Ψ (x −μ−Λz )
Qi(z(i))logp(x(i)|z(i);μ,Λ,Ψ)+logp(z(i))−logQi(z(i))dz(i) (5)
n 1 1(i) (i)T−1(i) (i)
E log (2π)d/2|Ψ|1/2 exp −2(x −μ−Λz ) Ψ (x −μ−Λz ) n 1 n 1 ( i ) ( i ) T − 1 ( i ) ( i )
7
i=1
E −2log|Ψ|−2log(2π)−2(x −μ−Λz ) Ψ (x −μ−Λz )
i=1
= ∇ΛE −tr2z(i) ΛT Ψ−1Λz(i) + trz(i) ΛT Ψ−1(x(i) − μ)
n1TT i=1
n1TT
= ∇ΛE −tr2ΛT Ψ−1Λz(i)z(i) + trΛT Ψ−1(x(i) − μ)z(i)
i=1 n
= E −Ψ−1Λz(i)z(i)T + Ψ−1(x(i) − μ)z(i)T i=1
Setting this to zero and simplifying, we get:
nn
ΛEz(i)∼Qi z(i)z(i)T = (x(i) − μ)Ez(i)∼Qi z(i)T .
i=1 i=1 Hence, solving for Λ, we obtain
n Λ =
Tn T−1
z(i) Ez(i)∼Qi z(i)z(i) . (7)
i=1
(x(i) − μ)Ez(i)∼Qi
It is interesting to note the close relationship between this equation and the
n Λ =
i=1
(x(i) − μ)μTz(i)|x(i)
n − 1 μz(i)|x(i) μTz(i)|x(i) + Σz(i)|x(i) . (8)
i=1
i=1
normal equation that we’d derived for least squares regression, “θT =(yTX)(XTX)−1.”
The analogy is that here, the x’s are a linear function of the z’s (plus noise). Given the “guesses” for z that the E-step has found, we will now try to estimate the unknown linearity Λ relating the x’s and z’s. It is therefore no surprise that we obtain something similar to the normal equation. There is, however, one important difference between this and an algorithm that performs least squares using just the “best guesses” of the z’s; we will see this difference shortly.
To complete our M-step update, let’s work out the values of the expecta- tions in Equation (7). From our definition of Qi being Gaussian with mean μz(i)|x(i) and covariance Σz(i)|x(i) , we easily find
Ez(i)∼Qi z(i)T = μTz(i)|x(i)
Ez(i)∼Qi z(i)z(i)T = μz(i)|x(i) μTz(i)|x(i) + Σz(i)|x(i) .
The latter comes from the fact that, for a random variable Y , Cov(Y ) = E[Y Y T ]−E[Y ]E[Y ]T , and hence E[Y Y T ] = E[Y ]E[Y ]T +Cov(Y ). Substitut- ing this back into Equation (7), we get the M-step update for Λ:
It is important to note the presence of the Σz(i)|x(i) on the right hand side of this equation. This is the covariance in the posterior distribution p(z(i)|x(i)) of z(i) give x(i), and the M-step must take into account this uncertainty
8
9
about z(i) in the posterior. A common mistake in deriving EM is to assume that in the E-step, we need to calculate only expectation E[z] of the latent random variable z, and then plug that into the optimization in the M-step everywhere z occurs. While this worked for simple problems such as the mixture of Gaussians, in our derivation for factor analysis, we needed E[zzT ] aswellE[z];andaswesaw,E[zzT]andE[z]E[z]T differbythequantityΣz|x. Thus, the M-step update must take into account the covariance of z in the posterior distribution p(z(i)|x(i)).
Lastly, we can also find the M-step optimizations for the parameters μ and Ψ. It is not hard to show that the first is given by
1 n μ=n x(i).
i=1
Since this doesn’t change as the parameters are varied (i.e., unlike the update for Λ, the right hand side does not depend on Qi(z(i)) = p(z(i)|x(i); μ, Λ, Ψ), which in turn depends on the parameters), this can be calculated just once and needs not be further updated as the algorithm is run. Similarly, the diagonal Ψ can be found by calculating
1nT T
Φ = n x(i)x(i) −x(i)μTz(i)|x(i) ΛT −Λμz(i)|x(i) x(i) +Λ(μz(i)|x(i) μTz(i)|x(i) +Σz(i)|x(i) )ΛT ,
i=1
and setting Ψii = Φii (i.e., letting Ψ be the diagonal matrix containing only the diagonal entries of Φ).