7

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 = NDSolve`FiniteDifferenceDerivative[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.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Faz
  • 1,817
  • 1
  • 14
  • 19
  • 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! – Faz Dec 30 '17 at 19:29
  • We just know that all the eigenvalues of the matrix $B$ have negative (and zero) real parts. – Faz Dec 30 '17 at 19:35
  • MatrixExp[1. B].vector is the solution of the parabolic flow associated to B(operator semi group generated by B) (similar to heat flow). Implicit Euler and Crank-Nicolson are two methods for computing the solution MatrixExp[t. B].vector iteratively (for all t between, say 0 and 1. It depends however a bit on the operator B: If it "diffusion dominated" (i.e., the symmetric part is dominant and semidefinite), then it might work very well; for "convection dominated" problems (e.g. significant nonsymmetric part), this may work not so well. – Henrik Schumacher Dec 30 '17 at 20:15
  • Moreover: What's "high-dimensional" and what is the PDE you try to solve. That's very important to know... And by the way: Where is the mass matrix? It should also show up if I am not mistaken... – Henrik Schumacher Dec 30 '17 at 20:18
  • Here "high-dimensional" comes from a 3D time-dependent parabolic PDE with variable coefficients. I think for computing MatrixExp[t.B].vecotr, a fast code was given in the answer of the folllowing: https://mathematica.stackexchange.com/questions/78278/ndsolve-in-mathematica-wont-use-all-the-cores-avaiable/78578#78578. – Faz Dec 30 '17 at 20:28
  • I know this is an old post but here is a code for Matlab but that could be applied to Mathematica. https://www.mathworks.com/matlabcentral/fileexchange/84058-fastexpm – Rob Dec 09 '20 at 14:07

1 Answers1

7

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 = NDSolve`FiniteDifferenceDerivative[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.

ybeltukov
  • 43,673
  • 5
  • 108
  • 212
  • 1
    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 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 09:48
  • 1
    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 Dec 31 '17 at 11:03
  • 2
    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 Dec 31 '17 at 15:09