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.