Perhaps I can clarify some of the issues related to the speed of fitting
mixed-effects models.

At present I am concentrating on Julia development and, in particular, on
linear mixed-effects models with partially-crossed grouping factors fit to
large data sets. This is because developing highly performant code is much
more straightforward in Julia than in R.   In fact, after a couple of
decades of working with R I feel it will be difficult to circumvent some of
the underlying assumptions in R to be able to write highly performant code,
even when using interfaces to compiled languages.

We must use a nonlinear optimizer to determine the ML or REML estimates of
the parameters in the model.  Essentially all the time in fitting complex
models to large data sets is spent in the optimization and all the time in
the optimization is spent in the evaluation of the objective function,

Evaluating the objective function involves solving a large, sparse,
penalized least squares problem, which is done using a sparse Cholesky
decomposition.  In some timings recently I have found that about 5/6 of the
time spend in the objective function evaluation is spent updating the
Cholesky decomposition (1/2 of the time) and solving systems of equations
using this decomposition (1/3 of the time).  This means that the only
possibility for substantial speed gains are in updating and applying this
decomposition.

Both the Matrix package in R and the base linear algebra capabilities in
Julia are based on the SuiteSparse software, which is single-threaded.  It
does allow for a "supernodal" Cholesky decomposition which we have found to
be important in achieving good performance on this type of model.

Some packages allow for a parallel sparse Cholesky.  On modern hardware it
is easiest to use multiple threads with a shared memory model rather than,
say, an MPI distributed memory model.

I am currently experimenting with one such package, PaStiX.  Initial
results are promising but much more development is yet to be done.  I
haven't even shown that a crucial calculation, evaluating the logarithm of
the determinant of the factor, can be done. Of course in theory it can be
done, it is just a matter of how awkward it will be to do it in practice.

An accelerated BLAS will help, especially for the supernodal decomposition,
which is a kind of sparse-dense hybrid.  However, you need to be careful
about using multithreaded BLAS when working with a parallel algorithm that
is itself multithreaded.  If you have 4 cores you don't want to have 4
threads of the algorithm each calling a BLAS routine that itself splits
into 4 threads.  In the end it usually works out best to throttle the
number of threads used by the BLAS and allow the top-level parallel
algorithm to use the threads.




On Sun, Apr 27, 2014 at 1:29 PM, Nathan Doogan <ndoogan@gmail.com> wrote:

> Thanks for the info, Ben.
>
> It sounds like the best way for me to speed things up, then, is to use a
> threaded linear algebra library.
>
> -Nate
>
> On 04/27/2014 02:00 PM, Ben Bolker wrote:
> >    [This is a perfectly reasonable question for r-sig-mixed-models, so
> > I'm forwarding it there]
> >
> > ===================
> > This might be very naive.
> >
> > I assume a very costly part of estimating parameters is evaluating the
> > likelihood -- particularly when the data are large. It seems like it'd
> > be fairly easy to distribute that evaluation across several processes
> > (i.e., multithread it) to speed up the procedure (e.g., mclapply or pvec
> > in R "parallel" package)
> >
> > That said, I'd guess likelihood evaluation in (g)lmer actually happens
> > somewhere else where I am not so comfortable.
> >
> > Any obvious reasons this map-reduce (I think it's called) sort of
> > technique is not in use?
> >
> > Thanks for your time. -Nate
> >
> > --
> > Nathan Doogan, Ph.D.
> > Post Doctoral Researcher
> > The Colleges of Social Work and Public Health
> > The Ohio State University
> >
> >
> >    It is indeed a little naive, but not silly at all.  The problem is
> > that it is *not* "fairly easy" to distribute the evaluation across
> > processors via map-reduce/mclapply etc..  Doug Bates has already done
> > all kind of wizardry to reduce the likelihood evaluation to linear
> > algebra operations that can be done very efficiently. While speeding the
> > process up enormously been doing some careful profiling with his Julia
> > code, and the largest single cost (if I am recalling things correctly)
> > is computing a sparse Cholesky decomposition.  He has been looking
> > *very* recently at the PaStiX library
> > <http://http://pastix.gforge.inria.fr/> as a possible way to parallelize
> > this operation, but it is not completely trivial.
> >    cheers
> >      Ben Bolker
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]

