brain-slam / slam

Surface anaLysis And Modeling
MIT License
13 stars 24 forks source link

Scipy for fast eigen components computation #97

Closed alexpron closed 1 year ago

alexpron commented 1 year ago

Scipy appears to be the most optimized solutions and easy to access (already part of the dependencies ) solution to reproduce eigen components computation of the discretized laplace beltami operator over a mesh.

We would like to find scipy function and parameters values that:

  1. Replicate former Matlab code
  2. In the fastest manner

A former intern suggested to use the eigsh function from the sparse module whose signature is described here (https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html)

As a first try I suggest to use these parameters that seems close to the matlab signature:

eigVal, eigVects = eigsh(lap.tocsr(), nb_eig, M=lap_b.tocsr(), sigma=None, which='SM')

where SM indicate the nb_eig smallest values are considered only.

Note that previously the parameter were :

eigVal, eigVects = eigsh(lap.tocsr(), nb_eig, M=lap_b.tocsr(), sigma=1e-6)

I performed some tests using the example data mesh and the first 100 eigen components extracted in matlab using suggested syntax Visually results are identical but this has to be further tested. Since the precision seems to be 10^⁻6 some other options may be considered like tol=10^-6 (truncation at some point) or sigma=10^-6 as exact computations appears to be super slow (around 13s for 2000 points mesh and the first 100 eigen components)

alexpron commented 1 year ago

Image

GerardMJuan commented 1 year ago

Hi,

I tried the code on my machine and the eigenvalues and eigenvectors obtained from MATLAB/Python seem to be identical within a 1e-6 margin:

Image

Both assertions are true. For the eigenvectors, some signs are flipped, but values are the same.

Code for computing them takes around 1.1s in MATLAB, compared to 20s in Python. Python code probably can be sped up with the right libraries.

alexpron commented 1 year ago

https://docs.scipy.org/doc/scipy/tutorial/arpack.html for speeding up the code