unibas-gravis / scalismo

Scalable Image Analysis and Shape Modelling
https://scalismo.org
Apache License 2.0
247 stars 65 forks source link

Feat/gp map #422

Closed madsendennis closed 9 months ago

madsendennis commented 9 months ago

MR related to https://github.com/unibas-gravis/scalismo/issues/421

Provides handle to directly calculate the .posterior.mean with .MAP without calculating the posterior kernel/covariance matrix.

madsendennis commented 9 months ago

I noted that in the regression function in GaussianProcess.scala, the inverse calculation is what takes up the majority of the time: val invKernel: DenseMatrix[Double] = breeze.linalg.inv(K). Does anyone have an idea how we can calculate this faster?

madsendennis commented 9 months ago

I added an attempt to remove the inv from the mean and kernel calculation with the Cholesky decomposition. The tests pass, and the mean values seem to be the same in the instances I tested, but there seems to be some difference in the kernel calculation. Maybe @JonathanAellen can help have a look, especially on the last commit?

madsendennis commented 9 months ago

I did some more test to check the difference for the inverse vs cholesky. They seem to provide the same output while cholesky is a lot faster to compute.

I see that that the MultivariateNormalDistribution case class has a private regularizedCholesky function. Should we move this out to a separate file and use instead in the GP file?

Andreas-Forster commented 9 months ago

I noted that in the regression function in GaussianProcess.scala, the inverse calculation is what takes up the majority of the time: val invKernel: DenseMatrix[Double] = breeze.linalg.inv(K). Does anyone have an idea how we can calculate this faster?

Maybe you can have a look at using native libraries like blas under the hood (see https://github.com/luhenry/netlib)

Andreas-Forster commented 9 months ago

For speed, there is a paper about cholesky and matrix inversion here: https://arxiv.org/abs/1111.4144

marcelluethi commented 9 months ago

@Andreas-Forster How about calling it posteriorMean instead of Map? This makes it clear that it is just a point estimate and also addresses the possible confusion with Scala's map function.