Speedup computeSigmaHat and computeSigmaHat.LISREL

[karchjd]

It seems that some algorithmic changes could increase the speed of the two functions.

computeSigmaHat.LISREL

In this function,

LAMBDA..IB.inv <- LAMBDA %*% solve(IB)

A generally faster way to do this would be

LAMBDA..IB.inv <- solve(IB, t(LAMBDA))

For recursive models (which I assume represent the majority of fitted models), IB should be a triangular matrix. In that case, the approach below should be the fastest:

forwardsolve(IB, t(LAMBDA), transpose = TRUE)

I tried the forwardsolve approach on the political democracy model, and the speedup was about 1.3×. For larger B matrices, I would expect the speedup to increase, since the forwardsolve approach has lower complexity than the inverse + matrix multiplication approach.

computeSigmaHat

Here, an eigenvalue decomposition is performed to check whether Sigma-hat is positive definite. I think this should be replaced by a Cholesky decomposition, since it is faster and is used anyway in the success case (positive definite Sigma). The Cholesky decomposition fails if the matrix is not positive definite, so in the success case one eigenvalue decomposition would be saved.

I tested this on a 70×70 sample covariance matrix, and it reduced lavInspect(fit, "timing")$optim from 0.66 seconds to 0.5 seconds.

I see that you are moving to C++, but since these are algorithmic changes, I thought I would report them anyway. I am curious what you think since I expect that you already considered these changes.

[yrosseel]

Indeed. You reached some of the oldest code in lavaan (you can tell that the lav_ naming convention is not even used here), and there are many, many things that need to be improved here from a numerical point of view. (Compare this with, for example, much more recent code in lav_mvnorm.R, where solve(X) is never used standalone).

I welcome pull requests. But I would do this in little steps, because we don't want to break anything that works now ;-)

[yrosseel]

Update:

• the computeSigmaHat() function has been fixed (see Added faster version of computeSigmaHat #442)
• the (old) functions in lav_model_compute() and lav_representation_lisrel() have been renamed (eg computeSigmaHat() is now lav_model_sigma())

Perhaps the next step could be to look at the use of forwardsolve()? The @model@modprop slot contains a check if the structural model (in 'beta') is acyclic (or not). Perhaps this can be useful.

Another long standing idea is to exploit a version of the Woodbury identity for computing the inverse (and determinant) of Sigma. This may speed up things a lot if the number of variables is large (say >50), while the number of factors is modest (say, 5). But it only seems useful if the 'Theta' matrix is diagonal (and contains no zero or negative values), and depending on the version of the Woodbury identity we use, it is not clear how can (efficiently) check if the final Sigma is positive definite.