Faster CFA without data augmentation

You might consider the SEM approach of modelling the covariance matrix rather than the individual observations. Given sample size n, sample covariance matrix S, and model-implied covariance matrix \Sigma, the log-likelihood is proportional to \log L = -\frac{1}{2} n \left[ \log \lvert\Sigma\rvert + tr(S \Sigma^{-1}) \right]. See Hoyle 2012 “Handbook of Structural Equation Modeling”, page 21.

Alternatively, you could use the Wishart distribution to model the scatter matrix (Q below), from which the equation above is derived.

data{
    int N; // Number of observations
    int P; // Number of indicators
    matrix[N,P] Y;  // Raw data
}
transformed data{
    matrix[N,P] Y_center;
    matrix[P,P] Y_scatter;

    for(p in 1:P){
        Y_center[,p] = Y[,p] - mean(Y[,p]);
    }

    Y_scatter = Y_center' * Y_center;
}
parameters{
    cov_matrix[P] Sigma;  // Model-implied covariance matrix
}
model{
    Y_scatter ~ wishart(N, Sigma);
}

Neither of these model the mean structure, though I suspect it may speed up the estimation of the covariance structure.