On Mon, Jun 2, 2014 at 9:23 PM, Asaf Weinstein <asafw.at.wharton@gmail.com>
wrote:

> Hello,
>
> On the more theoretical side in lme4:
> I am interested in obtaining the inverse of the sparse matrix M = U*Ut +
> I_q, the PLS matrix for (a permuted version of) which a sparse cholesky
> factor L is computed.
>

You really don't want to do that unless you have a very simple arrangement
such as a single random-effects term or multiple nested terms.  Even when L
is sparse, L^{-1} can be nearly dense.

The first approach should be to look at why you want the inverse and
whether there are alternatives that involve solving systems of equations
instead.  For example, many people assume that it is necessary to evaluate
(X'X)^{-1} to obtain the regression coefficients often written as
(X'X)^{-1}X'y but you don't need to do that.  (I cringe every time I see
the expression

solve(t(X) %*% X)

in R code).  Study the calculation and see if you really need M^{-1}.

If you conclude that there is no way around the M^{-1} calculation you are
probably better off creating M as a dense symmetric matrix (dpoMatrix in
the Matrix package) and inverting that.  On modern computers you can handle
sizes of M up to a few thousand without too much sweat.  Try on smaller
sizes and see how the execution time scales.  Do make sure that you write
the matrix as a symmetric matrix, which, in this case, will be positive
definite, so that you end up using the dense Cholesky factorization.

But first, check whether you really need M^{-1}.  Just because one form of
writing an expression involves that doesn't mean you need to calculate it.


> I wonder if
>
> solve(M, diag(q), system='A')
>
> will do the job as efficiently as
>
> t(P) %*% crossprod(Linv) %*% P
>
> where
>
> Linv=solve(Cholesky(M,LDL=FALSE), diag(q), system='L')
>
>
> Thank you very much,
> Asaf
>
>         [[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]]

