CODEWITHSUNDEEP
pythonnumpyscipy
Σ baryon
Σ baryon

Reputation: 246

How to optimse eigenvalue/vector computation of tridiagonal matrix

I am trying to optimise my computation for the eigenvalue/vectors of a tridiagonal matrix. My code currently takes about 30s to run and I need it to run much faster, at around 1s. I am currently using the scipy.sparse.linalg.eigsh method for finding the eigenvalue/vectors but I believe that there should be other ways which are much faster. I will post my code below so it is a little easier to see what I am trying to do. Any suggestions are greatly welcome and if there is anything unclear please let me know.

    H = (dx**-2)*diags([-1, 2, -1], [-1, 0, 1], shape=(N, N))
    H = sp.lil_matrix(H)
    H[0,0]=0.5*H[0,0]
    H[N-1,N-1]=0.5*H[N-1,N-1]
    lam, phi = eigsh(H, 400, which="SM")

Upvotes: 2

Views: 377

Answers (1)

senderle
senderle

Reputation: 150977

You're unlikely to find a simple way to speed up this calculation by 30x. As you can see from the timings in this post, even Google's impressively optimized jax library only manages a bit less than a 3x speedup over the base SVD implementation in scipy. Since SVD is equivalent to eigsh for symmetric matrices, we can expect that the speedup for eigsh will be about the same at best.

So even Google can't speed this calculation up that much, at least in the traditional way.

However, for very large, sparse matrices, there are specialized stochastic algorithms that can speed things up by much larger factors under the right circumstances. One of them has an sklearn implementation: sklearn.utils.extmath.randomized_svd. (I believe it's based on the algorithm described here, and it may also be the same as sklearn.decomposition.TruncatedSVD when algorithm='randomized'.)

Here's a slightly modified version of your example, showing how to use it:

from scipy.sparse import diags, lil_matrix
from scipy.sparse.linalg import eigsh

from sklearn.utils.extmath import randomized_svd

N = 3200
k = 400
dx = 1.0

H = (dx**-2)*diags([-1, 2, -1], [-1, 0, 1], shape=(N, N))
H = lil_matrix(H)
H[0, 0] = 0.5*H[0, 0]
H[N-1, N-1] = 0.5*H[N-1, N-1]
# lam, phi = eigsh(H, k, which="SM")

U, s, V = randomized_svd(H, k)

Here, s contains the singular values (equivalent here to eigenvalues), and U and V contain the singular vectors (equivalent here to eigenvectors). In theory, if H is a symmetric matrix, U should be the same as V.T. I didn't find that to be the case for this example, which puzzles me a bit... but I'm going ahead and posting this anyway, because this does in fact run in about a second, and so might be a solution for you.

However, there is still one more big caveat. You're passing which="SM" to eigsh, which as I understand it means you are asking for the 400 smallest eigenvalues. Is that really what you want? In almost all applications I am aware of, you'd want the 400 largest eigenvalues. If in fact you really do want the 400 smallest eigenvalues, then this won't work, because it relies on the fact that the largest eigenvalues are easier to coax out using random matrix projections.

The upshot is, you will want to test this to see if it actually gives results you can use. But it's an interesting possible solution, given what you've said so far about your problem.

Upvotes: 2

Related Questions