Properly implement `normest1`
Opened this issue · 8 comments
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.
A small variation of (JuliaLang/julia#12592) largely addresses this issue. Something like it integrated into ExpmV.jl shortly.
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
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.
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.
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.
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 😄
@marcusps can you tell me more about this issue? Where does the problem lies, It is still quite unclear to me!
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.