Re: the competition's expm vs ours

[Jaroslav Hajek]

On Tue, Dec 9, 2008 at 6:47 PM, John W. Eaton <address@hidden> wrote:
> On  9-Dec-2008, Jaroslav Hajek wrote:
>
> | On Sat, Nov 22, 2008 at 9:51 AM, Marco Caliari <address@hidden> wrote:
> | > Dear all,
> | >
> | > sorry for the delay. I sent an email but probabily it went lost. In Octave
> | > expm there are first some "reductions" (trace reduction and balancing) and
> | > then Pade' approximation (8,8) with scaling and squaring.
> | > On the other hand, Matlab does not make any reductions and uses Pade'
> | > approximation (6,6) with scaling and squaring.
> | >
> | > I have an m function implemeting exactly what Octave does with expm (like
> | > expmdemo1.m for Matlab). I can send it out of the list.
> | >
> | > Best regards,
> | >
> | > Marco
> |
> | Hi all,
> |
> | I've just written an m-file version of expm based on Marco Caliari's
> | demo above, and compared
> | with the built-in expm using recent tip (note that balance was
> | modified to be Matlab compatible).
> | Attached.
> |
> | A quick test on a 1000x1000 random matrix shows 5.5 seconds for the
> | built-in version, 3.8 seconds for the m-file version.
> | Anyone has objections against replacing?
>
> I don't object to replacing it, but I find the speed difference a
> little surprising.  What accounts for the speedup?
>

See my last reply to Marco.

> | I guess an m-file version is going to be much easier to play with.
>
> Yes.
>
> Why not write
>
>  r ^= (2^s);
>
> instead of
>
>  # Undo scaling by repeated squaring
>  for k = 1:s
>    r ^= 2;
>  endfor
>
> Hmm.  If I make this change the timings are once again approximately
> equal (so that answers my question above).  And I thought that
> internally, we were doing essentially the same thing when performing
> an operation like x^N but apparently not. Well, there's another thing
> to optimize if anyone is interested (look at src/xpow.cc).
>

xpow implements an O(log N) algorithm, but it does not minimize the
number of multiplications. That's an NP class problem, I think, so it
would be very hard, for small matrices and large exponents even harder
than the work itself :)
In r^(2^s) we know the optimal path, so we do it explicitly.

> BTW, will you please use ## for comments that are indented with the
> code?  That way, Emacs will automatically indent them properly.
>

I don't basically object, but is it a good idea to tailor coding style
specifically for Emacs?
Not everyone uses Emacs (I use ViM). Using ## for comments is OK, I'm
more annoyed by those weird trailing lisp-like comment blocks in C++
sources. I always thought Emacs was even better customizable than ViM,
so can't such settings be configured in Emacs configuration files?


-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz