[SciPy-User] Diagonalization of a tridiagonal, symmetric, sparse matrix with scipy.sparse.linalg.eigsh()

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] Diagonalization of a tridiagonal, symmetric, sparse matrix with scipy.sparse.linalg.eigsh()

Simone Bolognini

Hi everyone! Here's my problem.

I have an NxN symmetric and tridiagonal matrix computed by a Python code and I want to diagonalize it.

In the specific case I'm dealing with N = 6000, but the matrix can become larger. Since it is sparse, I assumed the best way to diagonalize it was to use the algorithm scipy.sparse.linalg.eigsh(), which performed extremely good with other sparse and symmetric matrices (not tridiagonal ones, though) I worked with. In particular, since I need only the low lying part of the spectrum, I'm specifying k=2 and which='SM' in the function.

However, in this case this algorithm seems not to work, since after approximately 20 minutes of computation I get the following error:

ArpackNoConvergence: ARPACK error -1: No convergence (60001 iterations, 0/2 eigenvectors converged)

Why is this happening? Is it a problem related to some properties of tridiagonal matrices? What other scipy routine can I use in order to diagonalize my matrix in an efficient way?

Here's a minimal code to reproduce my error:

import scipy.sparse.linalg as sl
import numpy as np

dim = 6000
a = np.empty( dim - 1 )
a.fill( 1. )
diag_up = np.diag( a, 1 )
diag_bot = np.diag( a, -1 )

b = np.empty( dim )
b.fill( 1. )

mat = np.diag( b ) + diag_up + diag_bot
v, w = sl.eigsh(mat, 2, which = 'SM')

On my pc the construction of the matrix mat takes 364ms, while only the last line in which the diagonalization is performed gives the reported error.


Thanks a lot for the support.

Simone


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Diagonalization of a tridiagonal, symmetric, sparse matrix with scipy.sparse.linalg.eigsh()

Luis JM Amoreira
Hi
Try scipy.linalg.solve_banded.

You have to store the matrix in diagonal ordered form. Take a random
matrix example:

In [1]: import scipy.linalg

In [2]: m=scipy.rand(3,6000)

In [3]: m[0,0]=0.; m[-1,-1]=0.

In [4]: b=scipy.rand(6000)

In [5]: %timeit scipy.linalg.solve_banded((1,1),m,b)
The slowest run took 49.27 times longer than the fastest. This could
mean that an intermediate result is being cached.
1000 loops, best of 3: 262 µs per loop

In steps 2 and 3 I create the diagonal ordered form. Each diagonal of
the coefficient matrix (I considered a random tridiagonal 6000x6000
matrix) is stored in each line of matrix m; zero is taken for the first
element of the upper diagonal and the last element of the lower diagonal.

Hope this helps, but your problem may be harder than my ranom example...

Regards
Ze Amoreira

On 09/22/2017 08:52 AM, Simone Bolognini wrote:

> Hi everyone! Here's my problem.
>
> I have an NxN symmetric and tridiagonal matrix computed by a Python code
> and I want to diagonalize it.
>
> In the specific case I'm dealing with |N = 6000|, but the matrix can
> become larger. Since it is sparse, I assumed the best way to diagonalize
> it was to use the algorithm |scipy.sparse.linalg.eigsh()|, which
> performed extremely good with other sparse and symmetric matrices (not
> tridiagonal ones, though) I worked with. In particular, since I need
> only the low lying part of the spectrum, I'm specifying |k=2| and
> |which='SM'| in the function.
>
> However, in this case this algorithm seems not to work, since after
> approximately 20 minutes of computation I get the following error:
>
>     ArpackNoConvergence: ARPACK error -1: No convergence (60001
>     iterations, 0/2 eigenvectors converged)
>
> Why is this happening? Is it a problem related to some properties of
> tridiagonal matrices? What other scipy routine can I use in order to
> diagonalize my matrix in an efficient way?
>
> Here's a minimal code to reproduce my error:
>
> |importscipy.sparse.linalg assl importnumpy asnp dim =6000a =np.empty(dim
> -1)a.fill(1.)diag_up =np.diag(a,1)diag_bot =np.diag(a,-1)b =np.empty(dim
> )b.fill(1.)mat =np.diag(b )+diag_up +diag_bot v,w =sl.eigsh(mat,2,which
> ='SM')|
>
> On my pc the construction of the matrix mat takes 364ms, while only the
> last line in which the diagonalization is performed gives the reported
> error.
>
>
> Thanks a lot for the support.
>
> Simone
>
>
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> https://mail.python.org/mailman/listinfo/scipy-user
>
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Diagonalization of a tridiagonal, symmetric, sparse matrix with scipy.sparse.linalg.eigsh()

Pauli Virtanen-3
In reply to this post by Simone Bolognini
pe, 2017-09-22 kello 09:52 +0200, Simone Bolognini kirjoitti:
[clip]
> However, in this case this algorithm seems not to work, since after
> approximately 20 minutes of computation I get the following error:
>
> ArpackNoConvergence: ARPACK error -1: No convergence (60001
> iterations, 0/2
> eigenvectors converged)
>
> Why is this happening? Is it a problem related to some properties of
> tridiagonal matrices?

It is not a question of matrix structure (tridiagonal, sparse) but the
eigenvalue spectrum. For your matrix, there are large number of
eigenvalues clustered to -1 which is the spectrum bottom, and this is
bad for Krylov methods.

> What other scipy routine can I use in order to
> diagonalize my matrix in an efficient way?

scipy.linalg.eig_banded
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user