views:

267

answers:

5

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.

A: 

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.

Anon
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
+2  A: 

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.

Yin Zhu
@ A = USV'. S is a square matrix with eigenvaules on the its diagonal.
Yin Zhu
+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.
Anon
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
@ Xodarap. Matlab's svds is not the best there. You could also try http://tedlab.mit.edu/~dr/SVDLIBC/. When you start, you first make sure SVD does your job by trying small matrices.
Yin Zhu
@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
@Xodarap. Very little. http://en.wikipedia.org/wiki/Singular_value_decomposition#Matrix_approximationYou can also do some experiments on small matrices to see the difference between approximation and the exact value.
Yin Zhu
+1  A: 

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)

gd047
I think inverse laplace runs in O(n^4), which is very slow. It's also prone to errors. See http://www.cs.cornell.edu/cv/ResearchPDF/19ways+.pdf
Xodarap
A: 

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

Marek
+1  A: 

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.

Marcus P S

related questions