In order to learn more about Factor Analysis, I've tried to implement a common model in R by hand, using MLE.

So I simulated data ( data ~ beta_1 + beta_2*x) . I employed PCA for generating starting values for X, and a normal distribution sample for betas. After defining a log likelihood, I started to run the iterations.

However, it takes so long to converge (although it does after one hour). I'm unsure if it is related to some adjustment I overlooked.

Anyway, my question is: is there some way of optimizing these iterations and estimating factor loadings in less time? If so, what are my mistakes here?

# Generating data
nbetas <- 100
nxs <- 36

betas <- matrix(rnorm(2*nbetas), ncol=2)
xs <- rnorm(nxs)
results <- matrix(0, nrow=nxs, ncol=nbetas)
for (i in 1:nxs)
  for (j in 1:nbetas)
    results[i,j] <- betas[j,1] + betas[j,2]*xs[i]

# A simple PCA recovers the xs
plot(cmdscale(dist(results))[,1], xs)

# Log Likelihood function
norm.loglike <- function(betas_h,xs_h,values) {
  soma <- 0
  for (i in 1:nxs)
    for (j in 1:nbetas)
      soma <- soma + dnorm(values[i,j], mean=(betas_h[j] + betas_h[j+nbetas]*xs_h[i])) 
# It seems that optim doesn't take on matrix, so I'm dealing with a vector 

xs_new <- cmdscale(dist(results))[,1] # starting values for x

# First estimation
betas_new <- optim(matrix(rnorm(nbetas*2), ncol=2), norm.loglike, method="L-BFGS-B",xs_h=xs_new, values=results)$par

# Iterations
for (i in 1:30) {
  xs_new <- optim(xs_new, norm.loglike, "L-BFGS-B",betas=betas_new, values=results)$par
  betas_new <- optim(betas_new, norm.loglike, method="L-BFGS-B",xs=xs_new, values=results)$par

# Plotting graphs to check them
plot(betas[,1], betas_new[,1])
plot(betas[,2], betas_new[,2])
  • 1
    $\begingroup$ If I am not mistaken, factorial solutions are undetermined (up to the product by an orthogonal matrix). Therefore, I would think that in order to achieve convergence to a single point you should give some conditions, for otherwise there is an infinity of points with equal likelihood. $\endgroup$ – F. Tusell Sep 12 '16 at 7:38
  • $\begingroup$ Thank you. But what conditions in this case? $\endgroup$ – Guilherme Duarte Sep 12 '16 at 7:42
  • $\begingroup$ Any you see fit. You might choose for instance to have 'beta2' to have a triangle of zeros, that will insure uniqueness of the solution, I think. $\endgroup$ – F. Tusell Sep 12 '16 at 7:55
  • 1
    $\begingroup$ You should be using expectation-maximization (EM) algorithm instead of using general solver to mazimize the likelihood. It should converge in seconds. $\endgroup$ – amoeba Sep 12 '16 at 15:01
  • 1
    $\begingroup$ That is not my point. I'm trying to implement another algorithm in order to show it to my students $\endgroup$ – Guilherme Duarte Sep 12 '16 at 17:50

Here is my programmmatic approach. This function, norm.loglike2.2(), is over 10 times faster than your norm.loglike() on my env (this function is exchangeable for yours because of the same arguments and return values).

norm.loglike2.2 <- function(betas_h, xs_h, values){
  mean_mat <- betas_h[1:nbetas] + matrix(betas_h[(1+nbetas):(2*nbetas)], ncol=1) %*% xs_h
  soma_v <- dnorm(values, t(mean_mat))
  soma <- sum(soma_v)
  • $\begingroup$ As far as I understand from your argument, did you just implement the same function in a matrix form or did you do something else? $\endgroup$ – Guilherme Duarte Sep 12 '16 at 12:46
  • $\begingroup$ @Guilherme Duarte; Yes, I just avoided using for() loop and took a vector and matrix based method to shorten a time for calculation. What it does is the same as yours. $\endgroup$ – cuttlefish44 Sep 12 '16 at 12:56
  • $\begingroup$ Indeed it is way faster. I'm looking for other solutions, though. I implemented with for() for didactic reasons. Anyway, based on my question, if no one answers here, I'll check your answer as correct. $\endgroup$ – Guilherme Duarte Sep 12 '16 at 13:00
  • $\begingroup$ @Guilherme Duarte; You don’t have to do that. I understood my answer isn't essential. $\endgroup$ – cuttlefish44 Sep 12 '16 at 13:06

This is not quite an answer, but too long for a comment. Let the observed multivariate variable be according to the factor model:

$\boldsymbol{X} = A\boldsymbol{F} + \boldsymbol{\epsilon}$

Then the covariance matrix of the $\boldsymbol{X}$ will be $\Sigma = AA' + \Psi$ where $\Psi$ is the diagonal matrix of variances of the unique factors. Clearly, $A$ cannot be unique, for any $A_* = AG$ with $G$ orthogonal will produce exactly the same $\Sigma$:

$\Sigma = AA' + \Psi = A_*A_*' + \Psi$

One way of making $A$ unique is to force it to have (for instance) an upper triangle of zeros. This also reduces the number of parameters over which you have to optimize. A quick surge with Google turned up


In page 18 you can find sketched the idea of how to implement a ML factor estimation.

  • $\begingroup$ Thanks for the tip about the upper triangle. Looks like the link you provided is dead, do you have an alternate resource you can share? $\endgroup$ – andrew_reece Apr 5 '18 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.