Optimizing Factor Analysis implementation in R - Maximum Likelihood

Question

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) {
  print(class(betas_h))
  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 
  print(-soma)
  return(-soma)
}


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) {
  print(i)
  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(xs,xs_new)
plot(betas[,1], betas_new[,1])
plot(betas[,2], betas_new[,2])

Answer 1

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){
  print(class(betas_h))
  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)
  print(-soma)
  return(-soma)
}

Answer 2

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

https://stat.ethz.ch/~maathuis/teaching/fall08/Notes5.pdf

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