search

marcusps / ExpmV.jl

Julia package to compute the result of expm(t*A)*v when A is a sparse matrix, without computing expm(t*A).
Other
22 stars 12 forks source link

Properly implement `normest1` #1

Open marcusps opened 9 years ago

marcusps commented 9 years ago

The code currently computes the 1 norm using norm(A^m,1) instead of estimating more efficiently. This ends up being the bottleneck in the calculation of not so sparse matrices.

The MATLAB code uses normest1(), a proprietary function. However, the algorithm for normest1 is described in

Nicholas J. Higham and Françoise Tisseur (2000) A block algorithm for matrix 1-norm estimation, with an application to 1-norm pseudospectra. SIAM Journal On Matrix Analysis And Applications, 21 (4). pp. 1185-1201. ISSN 1095-7162 (preprint)

and it should not be too hard to reimplemented it, thereby avoiding this performance bottle neck.

marcusps commented 9 years ago

A small variation of (JuliaLang/julia#12592) largely addresses this issue. Something like it integrated into ExpmV.jl shortly.

carandraug commented 8 years ago

For what it is worth, Octave has a free (GPL licensed) implementation of normest1 also based on Higham's paper which could be used since Matlab's version is proprietary http://hg.savannah.gnu.org/hgweb/octave/file/e35866e6a2e0/scripts/linear-algebra/normest1.m

marcusps commented 8 years ago

Thanks, but there is a branch with a proper implementation of normest1, based on work that was done in Base (Julia's core library). It is not GPL (has the same license as Julia itself).

It is not merged because I haven't gotten around to properly testing it.

thesharpshooter commented 7 years ago

Is this issue still open? @marcusps I would like to work on it if it is still open. I am currently working on Implementation of Julia solver of Exp Runge-Kutta method. The method requires calculation of exponentials of matrix times a vector. expm() performance is poor.

marcusps commented 7 years ago

Yeah, the issue is still open, and I am happy to have some help to finish this up.

Have a look at the normest1 branch, but make sure some write some more extensive tests. I'll push additional tests I wrote when I have a chance, but at one point it seemed like there were some numerical stability issues — I just never got around to tracking those down.

thesharpshooter commented 7 years ago

Other than this I would like to know, What else is left to the current Julia's implementation on Higham's expm and I would like to work on that? Can you refer me some reading materials on it? Thanks :smile:

thesharpshooter commented 7 years ago

@marcusps can you tell me more about this issue? Where does the problem lies, It is still quite unclear to me!

matteoacrossi commented 6 years ago

Any news about this? I'm eager to try ExpmV.jl for my code. I'm currently using Expokit.jl, but I'm really interested in the possibility to use this algorithm to span an interval for t.