From p.dalgaard at biostat.ku.dk  Wed Apr 11 18:54:59 2007
From: p.dalgaard at biostat.ku.dk (Peter Dalgaard)
Date: Wed, 11 Apr 2007 18:54:59 +0200
Subject: [Rd] eigen in beta
In-Reply-To: <461D032D.40304@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>	<Pine.LNX.4.64.0704111555430.28108@gannet.stats.ox.ac.uk>
	<461D032D.40304@bank-banque-canada.ca>
Message-ID: <461D12E3.5040206@biostat.ku.dk>

Paul Gilbert wrote:
> Prof Brian Ripley wrote:
>   
>> On Wed, 11 Apr 2007, Paul Gilbert wrote:
>>     
>>> Hmmm.  It is a bit disconcerting that make check passes and I can get 
>>> a fairly seriously wrong answer. Perhaps this test could be added to 
>>> make check.
>>>       
>> Well, you have learnt something new about software engineering!  'make 
>> check' is supposed to test the operation of R, not compilers.
>>     
> I'm not sure I would call this engineering. Engineers are suppose to 
> worry about the integrity of the whole system.  In any case, if it is 
> not make check,  then what is the system for checking that you have a 
> good R build?
>   

Well, there's make check-all... 

That doesn't check everything either, though. The point is that there is 
a limit to what we can check: We don't check whether the CPU gets 
floating point operations wrong in the 5th decimal place on rare 
occasions either. (Same thing applies to structural engineering: You 
can't test the chemical composition of every steel rod. At some point 
you have to shift reponsibility to your subcontractors.)

We do have regression checks though, i.e. we try to ensure that known 
errors do not reappear, and compiler issues have been worked around 
occasionally. E.g. (IIRC) the original qbeta() function was written so 
as to get the very last bit of accuracy squeezed out of the solution to 
a nonlinear equation, but modern optimizing compilers would reorder some 
instructions, losing a few bits of accuracy and sending  the code into 
an infinite loop in some cases.

    -pd


>>> Whether or not the one answer is correct may be questionable, but 
>>> there is no question that
>>>
>>> prod(eigen(z, symmetric = FALSE, only.values = TRUE)$values ) *
>>> prod(eigen(solve(z), symmetric = FALSE, only.values = TRUE)$values )
>>>  [1] 1.01677-0i
>>>
>>> is wrong. (The product of the determinants should equal the 
>>> determinant of the product, and the determinant of I is 1.)
>>>
>>> On this machine I am using GNU Fortran (GCC 3.2.3 20030502 (Red Hat 
>>> Linux 3.2.3-54)) 3.2.3 20030502 (Red Hat Linux 3.2.3-54),  which is a 
>>> bit old but not really old. (I am using gcc 4.1.1 on my home machine, 
>>> but the failing machine is suppose to be fairly new and "supported." 
>>> (There is an OS upgrade planned.)  Should I really be thinking of 
>>> this as an old compiler?  Is the
>>>       
>> Yes, a known bad compiler, as searching the list archives would have 
>> found out.  Don't use anything older that gcc 3.4.x on x86_64 Linux.
>>     
> Thanks, I'll bug the sys admin.
>
> Paul
>   
>>> compiler the most likely problem or is it possible I have a bad BLAS 
>>> configuration, or something else? Previous versions of R have 
>>> compiled without problems on this machine. (I am never very sure 
>>> where to find all the information to report for a problem like this. 
>>> Is there a simple way to get all the relevant information?)
>>>
>>> Paul Gilbert
>>>
>>> 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
>>>>>
>>>>>           
>>> ==================================================================================== 
>>>
>>>
>>> La version fran?aise suit le texte anglais.
>>>
>>> ------------------------------------------------------------------------------------ 
>>>
>>>
>>> This email may contain privileged and/or confidential information, 
>>> and the Bank of
>>> Canada does not waive any related rights. Any distribution, use, or 
>>> copying of this
>>> email or the information it contains by other than the intended 
>>> recipient is
>>> unauthorized. If you received this email in error please delete it 
>>> immediately from
>>> your system and notify the sender promptly by email that you have 
>>> done so. 
>>> ------------------------------------------------------------------------------------ 
>>>
>>>
>>> Le pr?sent courriel peut contenir de l'information privil?gi?e ou 
>>> confidentielle.
>>> La Banque du Canada ne renonce pas aux droits qui s'y rapportent. 
>>> Toute diffusion,
>>> utilisation ou copie de ce courriel ou des renseignements qu'il 
>>> contient par une
>>> personne autre que le ou les destinataires d?sign?s est interdite. Si 
>>> vous recevez
>>> ce courriel par erreur, veuillez le supprimer imm?diatement et 
>>> envoyer sans d?lai ?
>>> l'exp?diteur un message ?lectronique pour l'aviser que vous avez 
>>> ?limin? de votre
>>> ordinateur toute copie du courriel re?u.
>>>
>>>       
> ====================================================================================
>
> 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
>



[Rd] eigen in beta
Peter Dalgaard