15
$\begingroup$

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.

Improve this question
$\endgroup$
9
  • 6
    $\begingroup$ docs.google.com/… $\endgroup$
    – Steve Huntsman
    Commented May 28, 2010 at 14:02
  • 2
    $\begingroup$ You might find this article of interest cs.cornell.edu/cv/ResearchPDF/19ways+.pdf $\endgroup$
    – Guy Katriel
    Commented May 28, 2010 at 14:26
  • 1
    $\begingroup$ The two links posted so far are identical. $\endgroup$
    – j.c.
    Commented May 28, 2010 at 14:47
  • 1
    $\begingroup$ 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. $\endgroup$
    – Steve Huntsman
    Commented May 28, 2010 at 20:40
  • 4
    $\begingroup$ Xodarap, why are you exponentiating this matrix, and where does the problem come from? In other words, do you want the object $(e^A)$, or are you interested in computing its action on a given vector? These are (at the numerical linear algebra level) somewhat different questions. I'll be happy to point you to some references if you specify what you're trying to do. $\endgroup$
    – Nilima Nigam
    Commented Jun 23, 2011 at 5:26

9 Answers 9

Reset to default
11
$\begingroup$

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

Improve this answer
$\endgroup$
10
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ Al-Mohy and Higham's paper is great is you are dealing with sparse matrices. A preprint of the paper can be found on Higham's website, and he has MATLAB code that implements the algorithm. $\endgroup$
    – Marcus P S
    Commented Jan 15, 2014 at 0:50
8
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – Will Jagy
    Commented May 29, 2010 at 18:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – Xodarap
    Commented May 30, 2010 at 20:54
  • 1
    $\begingroup$ @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$. $\endgroup$
    – Nilima Nigam
    Commented Jun 23, 2011 at 15:02
5
$\begingroup$

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

Improve this answer
$\endgroup$
5
$\begingroup$

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.

Improve this answer
$\endgroup$
3
$\begingroup$

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

Improve this answer
$\endgroup$
2
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Do you know of a good way to diagonalize such a large matrix? Figuring out all 25k eigenvalues seems very time-consuming. $\endgroup$
    – Xodarap
    Commented May 28, 2010 at 15:42
  • $\begingroup$ 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. $\endgroup$
    – Gunnar Þór Magnússon
    Commented May 28, 2010 at 16:33
  • 3
    $\begingroup$ 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. $\endgroup$
    – Mark Meckes
    Commented May 28, 2010 at 16:58
  • $\begingroup$ A full diagonalization will not take advantage of sparsity. $\endgroup$
    – Terry Loring
    Commented Dec 22, 2013 at 18:20
2
$\begingroup$

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.

Improve this answer
$\endgroup$
2
$\begingroup$

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

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .