JenniNiku / gllvm

Generalized Linear Latent Variable Models
https://jenniniku.github.io/gllvm/
48 stars 19 forks source link

Potential matrix inverse improvements in C++ #83

Closed BertvanderVeen closed 1 year ago

BertvanderVeen commented 1 year ago

In many places in gllvm.cpp there is explicit inversion of matrices followed by a product with a vector, where that might be unnecessary. For example: (A.col(i).matrix()).inverse()*u.row(i).transpose() can instead be calculated by solving the system of equations so that explicit inversion of the matrix can be avoided: (A.col(i).matrix()).lu().solve(u.row(i).transpose()).

This should be less computationally expensive and more numerically stable (see e.g., https://gregorygundersen.com/blog/2020/12/09/matrix-inversion/ for some details). However, it is unclear if it is preferable to solve a system twice in favor of a single explicit inversion, which would often be required due to the nature of the calculations.

In some cases, it is also possible to instead form the inverse of A by the inverse of its cholesky factor. Normally this would come at the expense of having to calculate the factorization, but in gllvm it comes without extra cost due to the parameterization, and the inverse of a cholesky factor is again computationally cheaper than the full inverse.

matrix<Type> Ql(nlvr,nlvr);
matrix>type> Q(nlvr,nlvr);
 matrix<Type> I(nlvr,nlvr);
I.setIdentity();
Ql = A.col(i).matrix().template triangularView<Eigen::Lower>().solve(I);
Q = Ql*Ql.transpose();

edit: r-oldrel does not seem to accept treating triangular matrices this way..

Finally, according to the TMB pages (https://kaskr.github.io/adcomp/group__matrix__functions.html#gacf8e4c3cd2dc2c2859f8c52ad40ccec0), inverse() should be preferred over ``atomic::matinv``` for small matrices. As such, use of the latter might be useful in cases of e.g., n by n covariance matrices as for structured row-effects, or spatial LVs.