How to compute efficiently A^(-1)*(1-exp(-A*h))?
2 views (last 30 days)
Show older comments
I'm experimenting with exponential integrators for ODEs and SDEs, and try to avoid the built-in matlab solvers. As defined here,
I need to compute the object
efficiently, where
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/253639/image.png)
- A is a large (up to 1e4 x1e4) sparse, complex and symmetric matrix.
- h is the small numerical timestep with typical value 1e-4.
What would be the best way do do so?
- Explicit ' A\(1-expm(A*h))' takes about 2000 sec on my desktop. (and I believe it is unable to handle singlular A)
- Diagonalizing ([V,D]=eig(full(A)))) and working with 'V*D^(-1)(1-exp(D*dt))*V^(-1)' is an other approach. I tried it before and don't remamber the computing time, but it won't be too much faster than the first alternative. In this method, singular eigenvalues can also be regularized manually with l'Hopital's rule.
Which approach is recommended? Is there perhaps a better alternative? Ideally, the approach should be able to handle singular A, but even without this feature, other approaches would be interesting.
I somehow have the feeling there must be a faster way, because the series expansion of the expression of interest only contains positive matrix powers, so there should be no necessity to invert matrices? It seems to me as such that the complexity of this object should inherently not be worse than that of
() alone.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/253640/image.png)
2 Comments
Steven Lord
on 11 Dec 2019
I'm experimenting with exponential integrators for ODEs and SDEs, and try to avoid the built-in matlab solvers.
I'm curious why you're trying to avoid the built-in ODE solvers. What do they do that you don't need or want? Or do they not do what you need or want?
Accepted Answer
Christine Tobler
on 9 Dec 2019
For most sparse matrices expm(A) will be dense, so that should be expected to be expensive with a 1e4-by-1e4 matrix. If you are computing A^(-1)*(1-exp(-A*h)) with the goal of multiplying it with a vector, as I remember there are Krylov subspace methods for computing expm(A)*x directly without building the matrix.
To avoid the singularities, you may be able to use the Taylor expansion explicitly, in a function handle passed to funm. There is a helper function expm1(x) which computes (exp(x) - 1) with smaller round-off error. However, there is no such helper for (exp(x)-1)/x. An alternative to a Taylor expansion would be to check for cases of x being close to zero, and just modifying the value for that case.
3 Comments
Christine Tobler
on 11 Dec 2019
Dear Wouter,
Yes, the combination of function names expm and expm1 can be really confusing, unfortunately.
I'm wondering about the goal of keeping the resulting matrix sparse: For some special matrices this could be close to sparse, but in general I would expect there to be very few elements that are close to zero, just like usually hapens for the inverse or the eigenvectors of a sparse matrix. So before spending too much time on dropping small elements to zero, I would suggest computing the value for a moderately-sized dense matrix similar to your example and checking how many elements are close to zero.
Using a Taylor expansion of the function (1-exp(h*x))/x seems like a good approach, it will certainly remove the problems with singularities.
An excellent book on the topic of computing matrix functions is "Functions of Matrices: Theory and Computation" by N. J. Higham. You can also type "edit expm" or "edit funm" to see how these are implemented in MATLAB code.
More Answers (0)
See Also
Categories
Find more on Operating on Diagonal Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!