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 here, 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... Commented Dec 30, 2017 at 17:59
• Moreover, it would be nice to know what B should represent... Commented Dec 30, 2017 at 18:04
• Maybe the implicit-Euler method (also known as backward-Euler) or the Crank-Nicolson method may help... Commented Dec 30, 2017 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!
– Faz
Commented Dec 30, 2017 at 19:29
• We just know that all the eigenvalues of the matrix $B$ have negative (and zero) real parts.
– Faz
Commented Dec 30, 2017 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[0] == v}, u, {t, 1, 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!
– Faz
Commented Dec 30, 2017 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. Commented Dec 31, 2017 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$.
– Faz
Commented Dec 31, 2017 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?
– Faz
Commented Dec 31, 2017 at 15:09