[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Bug in expm
From: |
John W. Eaton |
Subject: |
Bug in expm |
Date: |
Tue, 22 Apr 2008 14:14:29 -0400 |
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
- Bug in expm, Marco Caliari, 2008/04/22
- Bug in expm,
John W. Eaton <=
- Re: Bug in expm, Marco Caliari, 2008/04/23
- Re: Bug in expm, Jaroslav Hajek, 2008/04/24
- Re: Bug in expm, Marco Caliari, 2008/04/24
- Re: Bug in expm, Jaroslav Hajek, 2008/04/24
- Re: Bug in expm, Jaroslav Hajek, 2008/04/26
- Re: Bug in expm, Marco Caliari, 2008/04/26
- Re: Bug in expm, Jaroslav Hajek, 2008/04/26
- Re: Bug in expm, Marco Caliari, 2008/04/27
- Re: Bug in expm, John W. Eaton, 2008/04/29