Symbolic Cholesky factorization (& other things)

Just a quick update. I'm working on developing a version of evol.vcv() that can be constrained in various ways. evol.vcv() implements the method of Revell & Collar (2009), in which we fit two or more different VCV matrices for the Brownian process to different parts of a phylogenetic tree. The type of constraints we might like to impose include, say, forcing rates (diagonal elements) to be equal among different regimes; or rates to differ and covariances to be the same; etc.

Unfortunately, implementing this has been somewhat more complicated than I originally appreciated for a few reasons, not the least among them being that in the function evol.vcv() (code here), I perform ML optimization on the Cholesky factorized rate matrices, rather than on the original matrices. This is hard to avoid because the likelihood function is undefined where the VCV matrices are not positive definite, and this causes problems for the R optimizers optim() if ignored.

I haven't had as much time to address this as I'd hoped, but as a first pass I am working on solving the problem for 3 × 3 matrices and two rate regimes. Hopefully a general solution appropriate for any size matrix or number of regimes will follow.

I'm trying a couple of different tacks. First, I thought I would symbolically factorize a symmetric matrix and then perhaps I could optimize with respect to the parameters of the symbolic functions in the Cholesky factorization. Unfortunately, on closer examination, this is no more straightforward than the original problem! This is because the parameters of the symbolic functions have the same constraints as the original matrices! (Duh.)

My newest attempt is cheaper & simpler. I am now trying to optimize the original matrices, but just returning logL=-Inf for any set of proposed values that fail a test for positive-definiteness. I am also trying some of the different optimization methods implemented in optim() to see if this has an effect. So far this simpler approach is showing some promise. I will report back when I have some further results.