# find eigenvectors from a PCAmodel

54 views

### Michela Bisighini

Dec 21, 2020, 6:27:41 AM12/21/20
to scalismo

Hi, I'm Michela. I need a help, I'm having problems in calculating the eigenvectors on a model obtained with the PCA. I'm using ScalismoLab and I generated the model in this way:
val PCAmodel = StatisticalMeshModel.createUsingPCA (dc) .get
Later I found the basis Matrix:
val X = PCAmodel.gp.basisMatrix
And finally to get the eigenvectors I tried to calculate the covariance matrix using different methods, but none of them worked:
METHOD 1)
val G = X.t * X
val n = X.cols
val svd.SVD (u, s, _) = breeze.linalg.svd (G * (1.0 / (n-1)))

error: could not find implicit value for parameter op: breeze.linalg.operators.OpMulMatrix.Impl2 [breeze.linalg.DenseMatrix [Float], Double, That]

METHOD 2) In the second method I have used this function:
def decomposeCovarianceMatrix ( X: DenseMatrix [Double], threshold: Double = 1.0e-8 ): (DenseMatrix [Double], DenseVector [Double]) = {
val n = X.cols
val Cov = X * X.t * (1.0 / (n-1))
val svd.SVD (u, s, _) = breeze.linalg.svd (Cov)
val rank = s.toArray.count (_> threshold) (u (::, 0 until rank), s (0 until rank)) }
val Covariance = decomposeCovarianceMatrix (X)

error: type mismatch;
found: breeze.linalg.DenseMatrix [Float]
required: breeze.linalg.DenseMatrix [Double]

Is there a way to solve these problems or another way to find the eigenvectors of the covariance matrix of the model created using PCA?
Best regards
Michela

### Marcel Luethi

Dec 21, 2020, 6:52:31 AM12/21/20
to Michela Bisighini, scalismo
Dear Michaela,

The eigenvalues correspond to the variance of the model. You can access them simply be calling
model.gp.variance

While possible, it is not ideal to access the pcabasis and variance directly. The reason is that there are different ways to represent that basis (normalized or not, including the variance, ...) and how exactly it is represented should remain internal to scalismo. A better way is to access this information using the klBasis:

model.gp.klBasis.map(eigenpair=>  eigenpair .eigenvalue)
model.gp.klBasis.map(eigenpair => eigenpair.eigenfunction)

Best regards,

Marcel

ps. The reason your svd computation did not work is, that you are mixing float and double representations in your computation
pps. ScalismoLab is a toy that we developed to making it easier to start using Scalismo. In case you are working on a real project, you should consider using an IDE, as described here: https://scalismo.org/docs/

--
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 view this discussion on the web visit https://groups.google.com/d/msgid/scalismo/f0c3da65-0134-4805-9549-b4b74fde1269n%40googlegroups.com.

### Michela Bisighini

Dec 21, 2020, 9:53:10 AM12/21/20
to scalismo
Thanks for the advices.
There still remains the problem of finding the eigenfunctions, in fact when I run this code: model.gp.klBasis.map (eigenpair => eigenpair.eigenfunction) gives me as an output a vector composed as follows:
res10: Seq [scalismo.common.DiscreteVectorField [scalismo.geometry._3D, scalismo.geometry._3D]] = Vector (<function1>, <function1>, <function1>, <function1>, <function1>)
I don't understand the reason of this error
Best regards,
Michela

### Marcel Luethi

Dec 21, 2020, 11:14:09 AM12/21/20
to Michela Bisighini, scalismo
Dear Michela,

I just realized that I misspelled you name. I am sorry for that.

I should have been a bit more explicit in my last email. The KLBasis consists of n pairs, where each pair contains the eigenvalue with its associated eigenfunction. The eigenfunction is a function from PointId to Vector, which you can use to obtain the deformation represented by the eigenfunction at the particular point id. The following code illustrates this:

// take first 10 pointIds of the reference mesh
for (pointId <- model.referenceMesh.pointIds.take(10)) {

// loop through all the basis functions
for ((eigenpair, i) <- model.gp.klBasis.zipWithIndex) {
val valueAtId = eigenpair.eigenfunction(pointId)
println(s"value of \$i-th eigenfunction at point with id \$pointId = \$valueAtId")
}
}

In this way you can extract all the information you need.

Best regards,

Marcel

### Michela Bisighini

Dec 21, 2020, 12:24:48 PM12/21/20
to scalismo
So for each point of my model created using PCA there are n eigenfunctions (where n is the number of shapes considered as input to the model)? I expected that for each eigenvalue there was an associated eigenfunction. Probably I'm confusing the meaning of eigenvector and eigenfunction.
Michela

### Marcel Luethi

Dec 21, 2020, 2:48:00 PM12/21/20
to Michela Bisighini, scalismo
Your understanding is completely right. For each eigenvalue there is an associated eigenfunction. The number of eigenfunctions does not have anything to do with the number of points. What I tried to show with the example is that each eigenfunction is a function from PointId to Vector3D and hence can be evaluated at all the points (identified by point ids) of the reference mesh.

I hope this makes things a bit clearer. If not, don't hesitate to ask again. This is one of the points where the interpretation used in Scalismo differs from how PCA is usually interpreted and hence needs some getting used to.

Best regards,

Marcel

### Michela Bisighini

Dec 22, 2020, 3:15:52 AM12/22/20
to scalismo
I apologize for the insistence. I would like in the code of my project to be able to calculate the eigenfunctions associated with each eigenvalue (with this code:
model.gp.klBasis.map(eigenpair => eigenpair.eigenfunction)  it doesn't work).
Find the eigenfunctions for each eigenvalue would be of great help to proceed with statistic analysis.
Michela

### Marcel Luethi

Dec 22, 2020, 6:03:22 AM12/22/20
to Michela Bisighini, scalismo
Hi Michela

No need to apologize - these discussions on the mailing may help others to understand the concept and maybe will prompt us to do better documentation :-)

Unfortunately I am not sore which part you struggle with exactly and therefore can just summarize what the mechanism in Scalismo is:
- the method klBasis of the gp will return n eigenpairs.
- each eigenpair consists of the eigenvalue and the corresponding eigenfunction
- the eigenvector is just a double value, which can be interpreted as the total variance of the data that is represented by the associated eigenfunction
- the eigenfunction is a function mapping PointIds to Vectors. For a point with a given point Id  in your reference mesh, you can obtain the corresponding value of the eigenfunction by calling eigenpair.eigenfunction(pointId).

Can you point out which of the above points (if any) is confusing you?

Best regards,

Marcel

### Michela Bisighini

Dec 22, 2020, 6:40:58 AM12/22/20
to scalismo
The theory behind Scalismo I think is clear. What I would like to print from my model are the n eigenfunctions associated with the n eigenvalues (where n is the number of shapes in input to the model). Being a model created using PCA I expect n eigenfunctions in total and not n eigenfunctions for each point of my shape.
If I use this code:
femursPCAmodel.gp.klBasis
the output is: scalismo.statisticalmodel.DiscreteLowRankGaussianProcess.KLBasis [scalismo.geometry._3D, scalismo.geometry._3D] = Vector (Eigenpair (9000.203, <function1>), Eigenpair (2482.864, <function1>), Eigenpair (78046, function1>), Eigenpair (252.84692, <function1>), Eigenpair (2.7221907E-13, <function1>))
In this case I have considered n=5 , so I have n Eigenpair where each eigenpair contains the i-th eigenvalue and its eigenfunction, unfortunately the eigenfunction is not printed.
I hope my issue is clearer now. Maybe I'm missing some concepts
Thanks
Michela

### Marcel Luethi

Dec 22, 2020, 7:12:03 AM12/22/20
to Michela Bisighini, scalismo
I tried to refactor the previous example a bit to make it more clear what is happening.

def printValuesOfEigenfunction(model : StatisticalMeshModel, eigenFun : DiscreteField[_3D, Vector[_3D]]) : Unit = {
for (pointId <- model.referenceMesh.pointIds.take(10)) { // only take the first 10 points to avoid too much output on the console
println(eigenFun(pointId))
}
}

val model = StatismoIO.readStatismoMeshModel(new File("datasets/bfm.h5")).get

for ((eigenpair, i) <- model.gp.klBasis.zipWithIndex) {
println(s"printing values of the \$i-th eigenfunction")
printValuesOfEigenfunction(model, eigenpair.eigenfunction)
}

This program should loop through all the eigenpairs, and calls the method printValuesOfEigenfunction for each eigenfunction. Within the method, the value of the eigenfunction is printed out for the first 10 points of the reference mesh.

Does this make it a bit clearer?

Best regards,

Marcel

### Michela Bisighini

Dec 23, 2020, 5:27:12 AM12/23/20
to scalismo
Thinking about it for few hours it becomes clearer to me.
Now, having obtained the eigenfunctions and the vector containing the eigenvalues, I would like to generate the possible shapes by considering the average shape of the model and each eigenfunction with its corresponding eigenvalue.
I wrote this code:
val X = socketPCAmodel.gp.basisMatrix
val meanVec = invasaturaPCAmodel.gp.meanVector
val shape = meanVec + X (::, 0) * eigenvalues (0) // I consider the first eigenvalue and the first eigenfunction

Now I have a DenseVector that contains all the points of my shape. What I would like to do is to visualize my shape by convert this denseVector in a Trianglemesh.
Thank you very much for your help
Michela

### Marcel Luethi

Dec 23, 2020, 9:01:10 AM12/23/20
to Michela Bisighini, scalismo
Hi Michela,

Great to hear that things are getting clearer. It always needs some getting used to ...

If I understand you correctly, what you are trying to do is already implemented in Scalismo. The method is called "instance". You can provide a set of shape coefficients, and it gives you the triangle mesh corresponding to these shape coefficients. The following code illustrates this:

val faceModel = StatismoIO.readStatismoMeshModel(new File("datasets/bfm.h5")).get

// create vector of shape coefficients and set first coefficient to 3
val v = DenseVector.zeros[Float](faceModel.rank)
v(0) = 3.0f

// obtain instance corresponding to shape coefficients and visualize it
val mesh = faceModel.instance(v)
show(mesh, "mesh")

When you work with Scalismo, you (almost) never have to manipulate matrices yourself. Scalismo works on a higher abstraction level (Gaussian processes) and hides all the matrix computations from the user. The reason is that the matrices themselves do not carry sufficient information to do the computations. In particular, the matrices do not represent any geometric information. You couldn't tell based on the matrix alone how the points are ordered in the matrix, nor what the neighborhood relationships are. This is the reason why emphasize the interpretation of Gaussian processes so much, as Gaussian processes represent all the geometric information needed.

Best regards,

Marcel