Orthonormality of PCA basis matrix after LowRankGaussianProcess.approximateGP

Christophe Van Dijck

unread,

Dec 21, 2015, 8:15:14 AM12/21/15

to scalismo

Hi

I am using scalismo to expand my SSMs with GP modes. I would like to be able to plugin the resulting extended SSM back into my framework. However, I saw that the PCA basis of the extended SSM is not orthonormal..

My code :

val model: StatisticalMeshModel = StatismoIO.readStatismoMeshModel(ssmFile).get
    
val gausKernel = GaussianKernel[_3D](stdDev) * magnitude
val matrixValuedKernel : MatrixValuedPDKernel[_3D,_3D] = UncorrelatedKernel[_3D](gausKernel)
val combKernel = model.gp.interpolateNearestNeighbor.cov + matrixValuedKernel
    
val zeroMean = VectorField(RealSpace[_3D], (pt:Point[_3D]) => Vector(0,0,0)) 
val gp = GaussianProcess (zeroMean,combKernel)
    
val sampler = RandomMeshSampler3D(model.mean, 500, 42)
val lowrankGP : LowRankGaussianProcess[_3D, _3D] = LowRankGaussianProcess.approximateGP(gp, sampler, nModes) 
    
val extModel : StatisticalMeshModel = StatisticalMeshModel(model.mean, lowrankGP)

It's possible I'm not completely grasping the low rank GPs, but I would think they are orthonormalised as well when converting to a statistical mesh model...

Best wishes,

Christophe

Marcel Luethi

unread,

Dec 21, 2015, 3:12:18 PM12/21/15

to Christophe Van Dijck, scalismo

Hi Christophe,

Thanks for pointing this out. This is a problem in scalismo, which we were not yet aware of.

The problem is a bit subtle and occurs because the method we use for computing the eigenvectors / eigenfunctions in the case of PCA is not the same as the one we use for the more general low-rank approximation, which is needed in the case where we combine kernels.

It will take a bit of time to fix this problem in a conceptually sound way. However, in the meantime you can orthonormalize the statistical model yourself. You will find a function that does that here
https://gist.github.com/anonymous/051dc37d491812678dae

I also noticed that there is are two potential problems in your code: You should always sample points from the reference (model.referenceMesh) and not from the mean to compute the low rank approximation. Furthermore, the new Gaussian process is not zero mean in general, but the mean should be model.gp.mean (i.e. the mean that was estimated from the data). In case you used the mean as a reference for building your model, the code would be correct, but in the general case it might not work as expected.

You find these corrections also in the gist.

Best regards

Marcel

--
You received this message because you are subscribed to the Google Groups "scalismo" group.
To unsubscribe from this group and stop receiving emails from it, send an email to scalismo+u...@googlegroups.com.
To post to this group, send email to scal...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/scalismo/214e7805-954d-45db-b3c3-df6f912baf65%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Marcel Lüthi, PhD
Department of Mathematics and Computer Science
University of Basel
Tel: +41 61 267 04 42

Christophe Van Dijck

unread,

Dec 22, 2015, 1:38:11 AM12/22/15

to scalismo, christop...@gmail.com, marcel...@unibas.ch

Hi Marcel,

Thanks for the quick reply! Indeed, the mean was taken as reference. However, I'll change my code to be fully conform :) .

I'll try the code you send me, thanks for that, and will be awaiting the update in scalismo.