[Rd] eigen in beta
Martin Maechler
maechler at stat.math.ethz.ch
Wed Apr 11 17:53:17 CEST 2007
>>>>> "PaulG" == Paul Gilbert <pgilbert at bank-banque-canada.ca>
>>>>> on Wed, 11 Apr 2007 10:51:22 -0400 writes:
PaulG> Hmmm. It is a bit disconcerting that make check passes and I can get a
PaulG> fairly seriously wrong answer. Perhaps this test could be added to make
PaulG> check.
Yes. If I set
EV <- function(x) eigen(x, symmetric = FALSE, only.values = TRUE)$values
dDet <- function(x) prod(EV(x))
something like
(pp <- dDet(z) * dDet(solve(z)))
all.equal(pp, 1 + 0i)
Can you see the problem also with e.g.,
z. <- round(z, 4)
or even
z. <- round(z * 100)
that would make the example code slightly nicer "to the eye"
...
PaulG> check. Whether or not the one answer is correct may be questionable, but
PaulG> there is no question that
PaulG> prod(eigen(z, symmetric = FALSE, only.values = TRUE)$values ) *
PaulG> prod(eigen(solve(z), symmetric = FALSE, only.values = TRUE)$values )
PaulG> [1] 1.01677-0i
PaulG> is wrong. (The product of the determinants should equal the determinant
PaulG> of the product, and the determinant of I is 1.)
yes, indeed (included in my code above)
PaulG> On this machine I am using GNU Fortran (GCC 3.2.3 20030502 (Red Hat
PaulG> Linux 3.2.3-54)) 3.2.3 20030502 (Red Hat Linux 3.2.3-54), which is a
PaulG> bit old but not really old. (I am using gcc 4.1.1 on my home machine,
PaulG> but the failing machine is suppose to be fairly new and "supported."
PaulG> (There is an OS upgrade planned.)
I'd say your OS version (using a 2.4 kernel) is quite old indeed.
PaulG> Should I really be thinking of this as an old compiler?
maybe not old, but buggy?
PaulG> Is the compiler the most likely problem or is it
PaulG> possible I have a bad BLAS configuration, or
PaulG> something else? Previous versions of R have compiled
PaulG> without problems on this machine. (I am never very
PaulG> sure where to find all the information to report for
PaulG> a problem like this. Is there a simple way to get all
PaulG> the relevant information?)
"simple" yes: the config.log file in your build directory,
but that's typically much more than you'd want in a
specific case, i.e. contains much more than just the
relevant information.
Regards,
Martin
PaulG> Paul Gilbert
PaulG> Prof Brian Ripley wrote:
>> All the systems I tried this on give the 'correct' answer, including
>>
>> x86_64 Linux, FC5 (gcc 4.1.1)
>> i686 Linux, FC5
>> ix86 Windows (both gcc 3.4.5 and gcc pre-4.3.0)
>> Sparc Solaris, with gcc3, gcc4 and SunPro compilers.
>>
>> Mainly with R 2.5.0 beta, some with R-devel (where the code is
>> unchanged).
>>
>> We have seen problems specific to RHEL's Fortran compilers on x86_64
>> several times before. I would strongly recommend compiler updates.
>>
>>
>> On Tue, 10 Apr 2007, Peter Dalgaard wrote:
>>
>>> Paul Gilbert wrote:
>>>> Here is the example. Pehaps others could check on other platforms.
>>>> It is
>>>> only the first eigenvalue that is different. I am relatively sure the
>>>> old values are correct, since I compare with an alternate calculation
>>>> using the expansion of a polynomial determinant.
>>>>
>>>>
>>>> z <- t(matrix(c(
>>>> 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0064083373167516857,
>>>> -0.14786612501440565826, 0.368411802235074137,
>>>> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0568624483195125444,
>>>> 0.08575928008564302762, -0.101993668348446601,
>>>> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0039684327579889069,
>>>> -0.00002857482925046247, 0.202241897806646448,
>>>> 1, 0, 0, 0, 0, 0, 0, 0, 0, -0.0222834092601282285,
>>>> -0.09126708346036176145, 0.644249961695308682,
>>>> 0, 1, 0, 0, 0, 0, 0, 0, 0, -0.0032676036920228878,
>>>> 0.16985862929849462888, 0.057282326361118636,
>>>> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0.0148488735227452068,
>>>> -0.06175528918915401677, 0.109566197834008949,
>>>> 0, 0, 0, 1, 0, 0, 0, 0, 0, -0.0392756265125193960,
>>>> 0.04921079262665441212, 0.078176878215115805,
>>>> 0, 0, 0, 0, 1, 0, 0, 0, 0, -0.0013937451966661973,
>>>> 0.02009823693764142133, -0.207228935136287512,
>>>> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.0273358858605219357,
>>>> 0.03830466468488327725, 0.224426004034737836,
>>>> 0, 0, 0, 0, 0, 0, 1, 0, 0, -0.1456426235151105919,
>>>> 0.28688029213315069388, 0.326933845656016908,
>>>> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0.0164670122082246559,
>>>> -0.21966261349875662590, 0.036404179329694988,
>>>> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.0146156940584119890,
>>>> 0.07505490943478997090, 0.077660578370038813
>>>> ), 12, 12))
>>>>
>>>>
>>>> R-2.5.0 gives
>>>> > eigen(z, symmetric = FALSE, only.values = TRUE)$values
>>>> [1] 0.8465266+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>>>> [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>>>> [7] 0.6177419+0.0000000i -0.5604582+0.1958709i -0.5604582-0.1958709i
>>>> [10] 0.1458799+0.4909300i 0.1458799-0.4909300i 0.3378356+0.0000000i
>>>>
>>>> R-2.4.1 and many, many previous versions gave
>>>> > eigen(z, symmetric = FALSE, only.values = TRUE)$values
>>>> [1] 0.8794798+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>>>> [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>>>> [7] -0.5604582+0.1958709i -0.5604582-0.1958709i 0.5847887+0.0000000i
>>>> [10] 0.1458799+0.4909300i 0.1458799-0.4909300i 0.3378356+0.0000000i
>>>>
>>>> Sys.info()
>>>> sysname
>>>> release
>>>> "Linux"
>>>> "2.4.21-40.ELsmp"
>>>> version
>>>> nodename
>>>> "#1 SMP Thu Feb 2 22:13:55 EST 2006"
>>>> "mfa04559"
>>>> machine
>>>> "x86_64"
>>>>
>>>> Paul Gilbert
>>> Hmm, I don't get that
>>>
>>>> version$version.string
>>> [1] "R version 2.5.0 beta (2007-04-10 r41105)"
>>>> eigen(z, symmetric = FALSE, only.values = TRUE)$values
>>> [1] 0.8794798+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>>> [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>>> [7] -0.5604582+0.1958709i -0.5604582-0.1958709i 0.5847887+0.0000000i
>>> [10] 0.1458799+0.4909300i 0.1458799-0.4909300i 0.3378356+0.0000000i
>>>> Sys.info()
>>> sysname
>>> release
>>> "Linux"
>>> "2.6.20-1.2933.fc6"
>>> version
>>> nodename
>>> "#1 SMP Mon Mar 19 11:38:26 EDT 2007"
>>> "titmouse2.kubism.ku.dk"
>>> machine
>>> login
>>> "i686"
>>> "pd"
>>> user
>>> "pd"
>>>
>>>
>>>
>>> And
>>>
>>>> version$version.string
>>> [1] "R version 2.5.0 beta (2007-04-09 r41098)"
>>>> eigen(z, symmetric = FALSE, only.values = TRUE)$values
>>> [1] 0.8794798+0.0000000i -0.0280087+0.6244992i -0.0280087-0.6244992i
>>> [4] -0.2908409+0.5522274i -0.2908409-0.5522274i -0.6228929+0.0000000i
>>> [7] -0.5604582+0.1958709i -0.5604582-0.1958709i 0.5847887+0.0000000i
>>> [10] 0.1458799+0.4909300i 0.1458799-0.4909300i 0.3378356+0.0000000i
>>>> Sys.info()
>>> sysname release
>>> "Linux" "2.6.18.8-0.1-default"
>>> version nodename
>>> "#1 SMP Fri Mar 2 13:51:59 UTC 2007"
>>> "viggo"
>>> machine login
>>> "x86_64" "pd"
>>> user
>>> "pd"
>>>
>>>
>>>
>>> The latter should be the actual build used in the current beta tarball
>>> (which is what you used, right?).
>>>
>>> I would suspect one of the following:
>>>
>>> - RHEL compilers
>>> - over-optimizing compiler settings
>>> - system blas/lapack libraries
>>> - system glibc
>>>
>>> but Brian probably has more concrete information.
>>>
>>>> Prof Brian Ripley wrote:
>>>>
>>>>
>>>>> We are only aware of better behaviour from LAPACK 3.1 (which is what I
>>>>> suppose you are talking about, that is R compiled with its internal
>>>>> LAPACK).
>>>>>
>>>>> But in at least one case that means finding a complex set of
>>>>> eigenvalues where previously a real one was found.
>>>>>
>>>>> On Tue, 10 Apr 2007, Paul Gilbert wrote:
>>>>>
>>>>>
>>>>> I am having some trouble with a case where eigen in R-beta gives a
>>>>> different largest value than in previous versions of R. Other values
>>>>> seem to be the same. Before I spend too much time, is anyone aware
>>>>> of a
>>>>> problem (symmetric = FALSE, only.values = TRUE).
>>>>>>
>>>>> Paul Gilbert
>>>>> ====================================================================================
>>>>>>
>>>>>>
>>>>>>
>>>>> La version française suit le texte anglais.
>>>>>>
>>>>> ------------------------------------------------------------------------------------
>>>>>>
>>>>>>
>>>>>>
>>>>> This email may contain privileged and/or confidential
>>>>> inform...{{dropped}}
>>>>>>
>>>>> ______________________________________________
>>>>> R-devel at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>>>>
>>>>>>
>>>> ====================================================================================
>>>>
>>>>
>>>> La version française suit le texte anglais.
>>>>
>>>> ------------------------------------------------------------------------------------
>>>>
>>>>
>>>> This email may contain privileged and/or confidential
>>>> inform...{{dropped}}
>>>>
>>>> ______________________________________________
>>>> R-devel at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>>
>>>
>>> ______________________________________________
>>> R-devel at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>
>>
PaulG> ====================================================================================
PaulG> La version française suit le texte anglais.
PaulG> ------------------------------------------------------------------------------------
PaulG> This email may contain privileged and/or confidential information, and the Bank of
PaulG> Canada does not waive any related rights. Any distribution, use, or copying of this
PaulG> email or the information it contains by other than the intended recipient is
PaulG> unauthorized. If you received this email in error please delete it immediately from
PaulG> your system and notify the sender promptly by email that you have done so.
PaulG> ------------------------------------------------------------------------------------
PaulG> Le présent courriel peut contenir de l'information privilégiée ou confidentielle.
PaulG> La Banque du Canada ne renonce pas aux droits qui s'y rapportent. Toute diffusion,
PaulG> utilisation ou copie de ce courriel ou des renseignements qu'il contient par une
PaulG> personne autre que le ou les destinataires désignés est interdite. Si vous recevez
PaulG> ce courriel par erreur, veuillez le supprimer immédiatement et envoyer sans délai à
PaulG> l'expéditeur un message électronique pour l'aviser que vous avez éliminé de votre
PaulG> ordinateur toute copie du courriel reçu.
More information about the R-devel
mailing list