# Speeding up the matrix exponential computation for sparse matrices

I want to compute the action of matrix exponential on a vector. My matrix $B$ is of a very large size and it is sparse, e.g., it is written as follows:

size = 350;
vector0 = Range[0., 300., 300./(size - 1)];
B0 = NDSolveFiniteDifferenceDerivative[2, vector0,
"DifferenceOrder" -> 20]["DifferentiationMatrix"];

B = SparseArray@KroneckerProduct[B0, B0];
vector = Flatten@KroneckerProduct[vector0, vector0];
Dimensions[B]


Based on the documentation and the discussion in https://mathematica.stackexchange.com/questions/83730/possible-method-for-matrixexp/83732?noredirect=1#comment431557_83732, I keep going simply as comes next:

sol2 = MatrixExp[B, vector, Method -> "Krylov"]; // AbsoluteTiming


This took 31 seconds in my office laptop. It computes the action successfully, BUT it is getting slower and slower when the dimension grows. In the above example, the dimension of $B$ is $122500\times 122500$, while in my real practice I wish to compute this action for matrices of dimension up to million.

So, I wish to know that is it possible to speed up such a computation by imposing a Compile or by a tolerance or something like that? Since, I am not looking for a result which is correct up to 15 digits! A lower accuracy, maybe up to 6 digits is enough.

Any hints would be welcomed and thanks to Mr. Wizard for encouraging asking this question.

• Is this matrix built correctly? It is not symmetric. While that is okay mathematically, it might be easier if it was symmetric... – Henrik Schumacher Dec 30 '17 at 17:59
• Moreover, it would be nice to know what B should represent... – Henrik Schumacher Dec 30 '17 at 18:04
• Maybe the implicit-Euler method (also known as backward-Euler) or the Crank-Nicolson method may help... – Henrik Schumacher Dec 30 '17 at 18:50
• In general such a matrix is constructed in the process of solving high dimensional PDEs and after imposing the boundaries, which are of mixed type, the matrix does not admit a precise structure! $B$ is used as the coefficient matrix for a system of large ODEs after discretization. The above $B$ is just a simple example. I do not know what you mean with the backward-Euler here! – Fazlollah Dec 30 '17 at 19:29
• We just know that all the eigenvalues of the matrix $B$ have negative (and zero) real parts. – Fazlollah Dec 30 '17 at 19:35

Let me propose a simple solution using the integration of PDE. Maybe, you try to avoid it, but I did not find such request in the question.

size = 100;
vector0 = Range[0., 300., 300./(size - 1)];
B0 = NDSolveFiniteDifferenceDerivative[2, vector0,
"DifferenceOrder" -> 20]["DifferentiationMatrix"];

B = SparseArray@KroneckerProduct[B0, B0];
vector = Flatten@KroneckerProduct[vector0, vector0];
Dimensions[B]
(* {10000, 10000} *)

sol2 = MatrixExp[B, vector, Method -> "Krylov"]; // AbsoluteTiming
(* {0.205314, Null} *)

exp[m_, v_] := Module[{f, u},
f[x_List] := m.v;
NDSolveValue[{u'[t] == -f[u[t]], u == v}, u, {t, 1, 1}]
]

sol3 = exp[B, vector]; // AbsoluteTiming
Norm[sol2 - sol3]/Norm[sol2]
(* {0.060318, Null} *)
(* 4.45283*10^-8 *)


I have reduced size because the result of MatrixExp seems to be unstable for large values of size. The proposed solution uses all cores of CPU since the multiplication of SparseArray and a vector is parallelized.

• This is a nice response. However, the absolute error is not that accurate, but it is OK. I am just worrying that when the size grows, the MatrixExp gives a better response for me. I mean in practice, it is more reliable rather than this NDSolveValue approach. Please, write more about what is happening in the background of this approach. I still think that the solution based on your approach loses much accuracy and stability for larger values of size! – Fazlollah Dec 30 '17 at 22:14
• @Fazlollah I'm sorry, I did a mistake because sol2 is almost the same as vector in this case. There should be f[x_List] := m.x; and u'[t] == f[u[t]]. As a result, the code becomes much slower. I have an another idea, but it takes a more time to think about it. Do you have a more reliable example of B? Moreover, it seems to me that some eigenvalues of B have positive real part. – ybeltukov Dec 31 '17 at 9:48
• I think with correcting this typo, the final results would be better, but the point is in this way, MatrixExp[matrix,vector] is faster than your approach! I try to edit the question by giving a practical $B$. – Fazlollah Dec 31 '17 at 11:03
• Another question is that what sort of "Krylov" method does the MMA use internally for the action in the question? GMRES, BiCGStab or something else? – Fazlollah Dec 31 '17 at 15:09