[SciPy-User] 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
> SciPy-User at python.org
> https://mail.python.org/mailman/listinfo/scipy-user
>