BAMresearch / bayem

Implementation and derivation of "Variational Bayesian inference for a nonlinear forward model." [Chappell et al. 2008] for arbitrary, user-defined model errors.
MIT License
2 stars 1 forks source link

Use only exp kernel. Rewrite vb with C_inv #82

Closed TTitscher closed 2 years ago

TTitscher commented 2 years ago

The only practical kernel (for now?) is the exponential kernel as we have the analytic inverse. So I removed the rest. Also, I adapted the simple_vb_with_c to use the inverse correlation matrix. The only annoying place where the sparse C_inv is turned to dense is the

sign, logdet = np.linalg.slogdet(C_inv)

as this method is not implemented for sparse matrices. A possible remedy is to use the identity

log(det(M)) = trace(log(M))

which would translate to

C_inv_log = scipy.linalg.logm(C_inv)
logdet = C_inv_log.trace()

But that fails with a error that C_inv is not square which IMO makes no sense...

ajafarihub commented 2 years ago

Concerning the det. issue (of the inverse matrix), maybe a way is to facilitate the method that returns the correlation matrix with an extra output; namely the det. of the C_itself or C_inv . This way, we would probably more easily compute the det. value. Also, at the end if we need to convert it to dense matrix and compute the inverse, I think it is not a big deal, since we do need to do it only once. So, maybe we can just change the method for free energy with an extra input:

free_energy(m, m0, L, L0, L_inv, s, s0, c, c0, k, J, C_inv, C_det)

TTitscher commented 2 years ago

There is a formula to get the determinant of a tridiagonal matrix in O(n) on wikipedia. This should be it!

TTitscher commented 2 years ago

So the formula on wikipedia works, but the determinant gets huge for large matrices. Thus the logdet unavoidable. As I found no trivial way to use the wiki formula to compute the logdet, I used a calculation based on a LU decomposition. At dim 10^6 that takes under a second and I didn't go further.

ajafarihub commented 2 years ago

Nice! I would say it is now a good point to merge the branch.