From maechler at stat.math.ethz.ch  Wed Apr 11 17:53:17 2007
From: maechler at stat.math.ethz.ch (Martin Maechler)
Date: Wed, 11 Apr 2007 17:53:17 +0200
Subject: [Rd] eigen in beta
In-Reply-To: <461CF5EA.4050001@bank-banque-canada.ca>
References: <461BCE69.2020409@bank-banque-canada.ca>
	<Pine.LNX.4.64.0704101854560.32417@gannet.stats.ox.ac.uk>
	<461BD7E9.1050200@bank-banque-canada.ca>
	<461BEDCB.7090805@biostat.ku.dk>
	<Pine.LNX.4.64.0704110734210.1007@gannet.stats.ox.ac.uk>
	<461CF5EA.4050001@bank-banque-canada.ca>
Message-ID: <17949.1133.668179.769813@stat.math.ethz.ch>

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



[Rd] eigen in beta
Martin Maechler