Exponential of large matrices I want to make a diffusion kernel, which involves $e^{\beta A}$, where A is a large matrix (25k by 25k). It is an adjacency matrix, so it's symmetric and very sparse.
Does anyone have a recommendation of a tool to solve this? I use the term "tool" loosely - if you know that transforming it in this way first or whatever is useful then I'd like to know that.

I am going with a hack - since the kernel "diffuses" relatively quickly, I just take only the neighbourhood around the two vertices that I want. This gives me a much reduced adjacency matrix which I can then raise e to without difficulty. 
I'm not familiar enough with the kernel function though to know how severely this is screwing up my results, and it's imperfect at best, so I'm still interested if anyone has a better idea.
 A: This is not an answer, but it's too long for a comment.
First, you need advice from a numerical analyst, not me.  Computing matrix exponentials is a well-studied problem with a large literature.  For one example, the recent book by Higham "Functions of matrices. Theory and computation" devotes a chapter to it. Matlab has a builtin routine for it.
The trick will be to take advantage of the sparseness, which almost certainly rules
out an approach based on diagonalization. Taylor series are not likely to help---try
computing $\exp(100)$ using the series expansion about $0$.
Also, just because you can write down the problem you want to solve using a matrix exponential,
does not guarantee this is the best way to solve it. (To give a crude example, the solution
to the linear system $Ax=b$ is $A^{-1}b$, but no-one in their right mind solves linear systems by computing inverses.)
A: I've asked for some clarification in a comment. In the meanwhile, 
if you're looking for software, I'll assume you've tried PETSc or Trilinos already? Here's a link to the freeware by Jiri Pittner, which links to BLAS routines as well:
http://www.pittnerovi.com/la/
Here's a site from INRIA http://verdandi.gforge.inria.fr/doc/linear_algebra_libraries.pdf
A: You can use the Chebyshev Polynomial expansion
to calculate the effect of  the matrix exponential on a vector. 
Which is a standard technique in quantum chemistry community and the method is extremely stable and fast. This method was developed by Tal-Ezer and Kossloff in an article named An accurate and efficient scheme for propagating the time dependent Schrödinger equation
You can see a Reviews of Modern Physics article by Alexander Wesse which deals with Kernal Polynomial Method (A generalization of the Chebyshev type algorithms).
I assume that to access these references you have the subscription to these scientific journals.
A: If you have a sparse matrix with localized effect (e.g. small valences), fast eigenvalue drop off and are required to compute the full matrix exponential, then you might be interested in 'diffusion wavelets'. While calculating the exponential they are as well calculating a basis where the result is still sparse.
Yet I am not aware of a ready to use implementation.
http://www.math.duke.edu/~mauro/Papers/DiffusionWavelets.pdf
http://en.wikipedia.org/wiki/Diffusion_wavelets
A: If your matrix is diagonalizable, say $A = PDP^-1$, then $\exp(A) = P \exp(D) P^-1$.  If your matrix is not diagonalizable and you need the more general Jordan Canonical Form, this approach may not work.  JCF is not suitable for numerical computation since it forming the JCF is a discontinuous process: arbitrarily close matrices can map to canonical forms that differ by an integer in one entry.
You could calculate $\exp(A)$ directly by its Taylor series.  Then the problem becomes how to efficiently calculate powers of $A$.  Maybe you could take advantage of your particular sparsity structure to calculate these powers.
A: Have a look at a recent paper discussing how matrix sparseness and locality go together:  "Decay Properties of Spectral Projectors with Applications to Electronic Structure" by Benzi et al. in SIAM Review, 55(1), 3--64, (2013).  The paper has applications that go beyond what the title indicates.  Much of the paper covers continuous functions applied to sparse hermitian matrices.
If you have some way of determining a priori which matrix elements will be small, you can compute a polynomial of the matrix quickly.  If your graph is related to a surface, you have an idea of how far apart on the graph two vertices need to be before they can be neglected.
To decide what polynomial to use, I would suggest you get an approximation of the operator norm.  This is fast for a sparse matrix.  In matlab you use normest.  In other languages see:  "Estimating the matrix p-norm" by Nicholas J. Higham, Numerische Mathematik,  62(1), 539--555, (1992).  The code there simplifies in the case $p=2$, which is the case you want.  This norm estimate, rounded up a bit for good measure, tells you where the spectrum of your matrix sits.  Now get (say from a truncated power series) a polynomial that is close enough for your purposes to the actual exponential on the spectrum of your matrix.
Even if you can't figure which matrix elements of the answer you will zero-out, if you can accept a modest error and so deal with a polynomial of relatively small degree, then you are just needing to compute several powers of a sparse matrix.  It is then a question of how-sparse you start with vs. how high a power you need.
I will warn you that I find Matlab does not do so well taking products of sparse matrices.  I think it is optimized for minimizing data storage, not matrix multiplication.
A: Suprised that no one mentioned Expokit, http://www.maths.uq.edu.au/expokit/
It does exactly what was requested, and is available in several different implementations (including Matlab).
A: The book by Higham and the "nineteen dubious ways" paper deal with the dense case only. For the sparse case, the best way to go is using an algorithm that computes the so-called action, i.e., the map $ v \mapsto \exp(A)v$. See e.g. Al-Mohy, http://epubs.siam.org/sisc/resource/1/sjoce3/v33/i2/p488_s1?isAuthorized=no.
The matrix $\exp(A)$ itself is full and unstructured, and generally you do not want to use it. If you really need it, though, check out a series of papers by Benzi and coauthors: they show that the off-diagonal elements of many matrix functions decay exponentially, and thus your matrix might be "nearly banded".
A: I've done 11kx11k matrices in Python with scipy's sparse matrix exponential function in tens of seconds: https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.sparse.linalg.expm.html
