dsuess / mpnum

Matrix Product Representation library for Python
BSD 3-Clause "New" or "Revised" License
54 stars 17 forks source link

Custom svd function #55

Open NicolaLorenzoni opened 3 years ago

NicolaLorenzoni commented 3 years ago

Hi everyone, I've just noticed that, in mparray.py, and specifically in compress_svd_l and compress_svd_r, with relerr we use svd instead of svdfunc, which does not enable to use a different svd function respect numpy.linalg.svd.

Screenshot from 2021-02-03 12-40-42

dsuess commented 3 years ago

The idea behind svdfunc was to provide different implementations for the truncated SVD. The most naive implementation just computes the full SVD and truncates at given rank. However, you can do much better using e.g. this implementation or randomized algorithms.

However, for relerr=True, we need to compute all singular values. Do you need to change the implementation of the SVD in this case as well?

NicolaLorenzoni commented 3 years ago

I see, many thanks for the answer. I stumbled in this because in my case the SVD was not reaching convergence. Searching more in the library, this is due the fact that in the relerr = True case we use numpy.linalg.svd, while in the relerr = None one scipy.linalg.svd. When the numpy doesn't convergence, using the lapack driver 'gsvd' within the scipy one gives no problem, since it is more stable and precise. So, the error I was stumbling into is not directly related to the truncated or the randomized method, but with the different svd we use in the relerr = None or relerr = True case.

dsuess commented 3 years ago

Ahh, interesting. I thought the implementations in scipy and numpy would be the same. I'll need to benchmark whether gsvd is slower than the current implementation and, if not, just change it out.