4

I want to compute a diffusion kernel, which involves taking exp(b*A) where A is a large matrix. In order to play with values of b, I'd like to diagonalize A (so that exp(A) runs quickly).

My matrix is about 25k x 25k, but is very sparse - only about 60k values are non-zero. Matlab's "eigs" function runs of out memory, as does octave's "eig" and R's "eigen." Is there a tool to find the decomposition of large, sparse matrices?

Dunno if this is relevant, but A is an adjacency matrix, so it's symmetric, and it is full rank.

5 Answers 5

2

Have you tried SVD, svds for sparse matrix in matlab.

EDIT: one more thing, don't do full rank SVD since the dimension is big, use a small rank, say 500, so that your solution fits in the memory. This cuts the small eigenvalues and their vectors out. Thus it does not hurt your accuracy much.

6
  • @ A = USV'. S is a square matrix with eigenvaules on the its diagonal.
    – Yin Zhu
    May 27, 2010 at 0:33
  • 1
    +1 @duffymo: See wikipedia on relationship between svd and eigenvalue decomposition. In any case the OP was asking for diagonalization of A which svd does.
    – vad
    May 27, 2010 at 0:37
  • I tried doing this on mine and it said "maximum variable size exceeded." My computer isn't exactly top-of-the-line though, so it could be a problem there instead of with matlab.
    – Xodarap
    May 27, 2010 at 0:44
  • 1
    @ Xodarap. Matlab's svds is not the best there. You could also try tedlab.mit.edu/~dr/SVDLIBC. When you start, you first make sure SVD does your job by trying small matrices.
    – Yin Zhu
    May 27, 2010 at 0:48
  • @Yin - this is an excellent idea, but do you know how much information is lost? I feel like since this is an adjacency matrix and the rows are really quite independent of each other, finding only 10% of the eigenvalues (or whatever) will lose 90% of the info.
    – Xodarap
    May 27, 2010 at 21:41
1

Have you considered the following property: exp(A*t) = L^(-1) {(sI-A)^(-1)} where L^(-1) the inverse Laplace transform? - provided that you can invert (sI-A)

1
1

If you have access to a 64 bit machine and octave compiled with 64 bit support, you might be able to get around this problem.

Also, I don't know what platform you are running all of this on, but in UNIX based systems you can use ulimit to increase the maximum allowed stack size for user processes.

For example, you can run

ulimit -u unlimited

and this will ensure that there are no memory limits etc on your processes. This is not a good idea in general, since you have have runaway processes that will completely bog down your machine. Try instead

ulimit -s [stacksize]

to increase the stack size limit.

0

Octave has splu which does lu decomposition for sparse matrices. I am not sure whether it can handle 25k x 25k but worth a shot.

Alternatively if your matrix is structured like so: A = [B zeros;zeros C] then you can diagonalize B and C separately and put them together into one matrix. I am guessing you can do something similar for eig.

1
  • Stupid question, but once I do have lu, how do I use this? I realize that det(A) = det(L)*det(U), so there must be some relationship between their diagonalizations, but I'm not smart enough to see it.
    – Xodarap
    May 27, 2010 at 20:22
0

In R you could check igraph package and function arpack which is interface to ARPACK library.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.