Bug in expm

[John W. Eaton]

On 22-Apr-2008, Marco Caliari wrote:

| Dear maintainers,
| 
| the bug "solved" here 
| http://www.cae.wisc.edu/pipermail/octave-maintainers/2008-January/005841.html
| has not been solved. It seems that the patch was not applied. In 
| practice, the trailing permutations should be performed before the leading 
| ones. I'm sorry I didn't notice that in Octave 3.0.0+.
| I enclosed the files used for checking: expmdemo1.m is the m-file 
| implementation of the built-in expm.

Can someone create a patch for the expm function in liboctave that is
based on this code?

Also, you might want to look at the current expm in R.  It was
originally derived from the code in Octave, but then there was a brief
discussion (unfortunately not on this or apparently any other mailing
list) earlier in the year about a bug and with a fix.  I don't know
whether the bug reported against R is the same or similar to the ones
reported for Octave.  If you are interested in fixing this, then
please see the code below and ask Martin about the details.

I'd consider a patch but do not have time to generate one from this
code or the function you sent.

  Subject: Re: Found and fixed bug !! {"Matrix exponential in R"}
  From: Martin Maechler <address@hidden>
  Reply-To: Martin Maechler <address@hidden>
  To: "John W. Eaton" <address@hidden>
  Cc: Martin Maechler <address@hidden>, John Eaton <address@hidden>,
          address@hidden
  Date: Thu, 28 Feb 2008 17:42:42 +0100
  Message-ID: <address@hidden>

  Hi John,

  >>>>> "JWE" == John W Eaton <address@hidden>
  >>>>>     on Thu, 28 Feb 2008 03:54:18 -0500 writes:

      JWE> On 27-Feb-2008, Martin Maechler wrote:
      JWE> | Gentlemen,
      JWE> | 
      JWE> | I am  *very*  happy to report
      JWE> | that thanks  to my progress  with the R-forge package 'expm',
      JWE> | with an added R-interface to the LAPACK  dgebal  routine
      JWE> | reading  Ward(1977),
      JWE> | and using the R version  of  
      JWE> |     "scaling + Pade + squaring"  ( =: s.P.s )
      JWE> |     the original of which Stig wrote,
      JWE> | I have been able to confirm that 
      JWE> | 
      JWE> | - it's not Ward(1977) which would be flawed
      JWE> |   {Ward does a bit more than s.P.s : he proposed two additional
      JWE> |    pre-conditioning steps to the "s + s" one; 
      JWE> |    and that's what the octave code was also aiming at},
      JWE> | 
      JWE> |   but rather
      JWE> | 
      JWE> | - it's the code we originally took from octave which has a bug :
      JWE> |    
      JWE> |   When reverting the permutation (from pre-conditioning step 2),
      JWE> |   the code does things in a wrong order.
      JWE> | 
      JWE> | I've replaced the (somewhat complicated) wrong code by much
      JWE> | easier to understand correct code, and have found that several
      JWE> | examples which gave wrong results previously no longer do so.
      JWE> | 
      JWE> | This is very good news, and I will also update the
      JWE> | expm() C version of Matrix with a non-buggy one.

      JWE> Hi,

      JWE> I didn't write the expm function in Octave and I'm not an expert in
      JWE> this area.  There was recently a discussion on the Octave lists about
      JWE> how Octave's expm fails for various matrices, but no fix has been
      JWE> proposed.  What are the changes that you've made?  Would you be
      JWE> willing to also prepare a patch for the current Octave code?  

  Can you tell me which source file I have to change?
  Then I'll try {if I can build octave from source Quickly,
  without contortions}

      JWE> If not, where can I find the changes you made to the matrix 
exponential
      JWE> function in R?

  The new code is on R-forge, and e.g. visible in
  
http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/src/expm.c?rev=22&root=expm&view=markup

  However it's easier to show you what change fixed the bug:

  
--------old-buggy-code---------------------------------------------------------
          /* Preconditioning 2: Inversion of 'dgebal()' :
           * ------------------ Note that dgebak() seems *not* applicable */

          /* Step 2 a)  apply inverse scaling -- TODO speedup: only inside 
{ilo:ihi} */
          for (j = 0; j < n; j++)
              for (i = 0; i < n; i++)
                  z[i + j * n] *= scale[i]/scale[j];

          /* 2 b) Inverse permutation  (if not the identity permutation) */

          if (ilo != 1 || ihi != n) {
              /* inverse permutation vector: */
              int *invperm = (int *) R_alloc(n, sizeof(int));

              /* balancing permutation vector */
              for (i = 0; i < n; i++)
                  invperm[i] = i;       /* identity permutation */

              /* leading permutations applied in forward order */
              for (i = 0; i < (ilo - 1); i++)
              {
                  int p_i = (int) (perm[i]) - 1;
                  int tmp = invperm[i]; invperm[i] = invperm[p_i]; invperm[p_i] 
= tmp;
              }

              /* trailing permutations applied in reverse order */
              for (i = n - 1; i >= ihi; i--)
              {
                  int p_i = (int) (perm[i]) - 1;
                  int tmp = invperm[i]; invperm[i] = invperm[p_i]; invperm[p_i] 
= tmp;
              }

              /* construct inverse balancing permutation vector */
              Memcpy(pivot, invperm, n);
              for (i = 0; i < n; i++)
                  invperm[pivot[i]] = i;

              /* apply inverse permutation */
              Memcpy(work, z, nsqr);
              for (j = 0; j < n; j++)
                  for (i = 0; i < n; i++)
                      z[i + j * n] = work[invperm[i] + invperm[j] * n];
          }

  
----END-old-buggy-code---------------------------------------------------------

  and here is the fixed ... much shorter! ... version of that

  
--------new-correct-code---------------------------------------------------------
          /* Preconditioning 2: Inversion of 'dgebal()' :
           * ------------------ Note that dgebak() seems *not* applicable */

          /* Step 2 a)  apply inverse scaling */
          for (j = 0; j < n; j++)
              for (i = 0; i < n; i++)
                  z[i + j * n] *= scale[i]/scale[j];

          /* 2 b) Inverse permutation  (if not the identity permutation) */

          if (ilo != 1 || ihi != n) {
              /* ---- new code by Martin Maechler ----- */

  #define SWAP_ROW(I,J) F77_CALL(dswap)(&n, &z[(I)], &n, &z[(J)], &n)

  #define SWAP_COL(I,J) F77_CALL(dswap)(&n, &z[(I)*n], &i1, &z[(J)*n], &i1)

  #define RE_PERMUTE(I)                         \
              int p_I = (int) (perm[I]) - 1;    \
              SWAP_COL(I, p_I);                 \
              SWAP_ROW(I, p_I)

              /* reversion of "leading permutations" : in reverse order */
              for (i = (ilo - 1) - 1; i >= 0; i--) {
                  RE_PERMUTE(i);
              }

              /* reversion of "trailing permutations" : applied in forward 
order */
              for (i = (ihi + 1) - 1; i < n; i++) {
                  RE_PERMUTE(i);
              }
          }

  
----END-new-correct-code---------------------------------------------------------

      JWE> Also, how long ago was it that the Octave code was adapted for R?

  I'd guess about two years.  Doug Bates did that.

      JWE> Thanks,
      JWE> jwe

  You're welcome,
  Martin




jwe