Question

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.

Answer 1

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

Answer 2

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

Answer 3

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

Answer 4

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

Answer 5

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.

Answer 6

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.