## 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.

-
docs.google.com/… – Steve Huntsman May 28 2010 at 14:02
Can you fix the title? You did not mean "exponent"... – Mariano Suárez-Alvarez May 28 2010 at 14:04
You might find this article of interest cs.cornell.edu/cv/ResearchPDF/19ways+.pdf – Guy Katriel May 28 2010 at 14:26
The two links posted so far are identical. – jc May 28 2010 at 14:47
In MATLAB you'll want to sparsify explicitly if you haven't already; the "sparse" command does this. Then use "eigs" (not "eig") to return the eigenvectors. Do what everyone else is saying (if your matrix is really that sparse, MATLAB should be up to it on a modern laptop) and then compare the results you obtain with "expm" (if you can). I'd be surprised if the calculation took more than a few minutes. – Steve Huntsman May 28 2010 at 20:40

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.)

-
Glad somebody sees this my way. I googled "diffusion kernel," this problem is very far from being simply about exponentiating matrices. Then I deleted my answer, nobody seemed interested. Paper by Kondor and Lafferty, presentations by Liang Sun and then Bruno Jedynak. – Will Jagy May 29 2010 at 18:26
Yeah, I must admit that when I asked this question I didn't realize it was so unsolved. I thought the answer would be "use the really-big-sparse-matrix add-on to Matlab" or something. That being said, sparse adjacency graphs (e.g. the web, genome mapping, etc.) appear all the time, and so I don't believe that there is no acceptable solution - I will accept that there is no perfect solution, but the problem seems too common for there to be no standard toolkit. – Xodarap May 30 2010 at 20:54
@Xodorap, to be clear: there are scores of excellent algorithms out there for sparse matrix operations. What we need to get from you is a categorical statement like 'I want the matrix exponential itself' or 'I want the solution of the diffusion equation $u_t=Au$ with given data'. There are lots of acceptable approaches in either case, but they are not the same. As Will and Chris point out, it is rare for someone to genuinely need $e^A$ for a large, symmetric and sparse $A$. – Nilima Nigam Jun 23 2011 at 15:02

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).

-

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".

-

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

-

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 Jordon 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.

-
Do you know of a good way to diagonalize such a large matrix? Figuring out all 25k eigenvalues seems very time-consuming. – Xodarap May 28 2010 at 15:42
You don't have to calculate all of the Taylor series. If you let P be the characteristic polynomial of the matrix, then you can write exp(A) = g(A) * P(A) + rest, where g is entire, and Cayley-Hamilton then gives exp(A) = rest (you can divide entire functions of matrices by polynomials). The rest can be calculated by finite differences, if I remember correctly. – Gunnar Magnusson May 28 2010 at 16:33
Xodarap says A is real symmetric, so it is indeed diagonalizable. So as Xodarap points out above, the real question is how to go about diagonalizing. – Mark Meckes May 28 2010 at 16:58

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.

-