Re: the competition's expm vs ours

[John W. Eaton]

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?

| 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).

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

Thanks,

jwe