>>>>> smallepsilon via R-help
>>>>> on Sun, 04 May 2025 18:11:57 +0000 writes:
> Peter, The eigenvalues are not identical(), but are
> all.equal(). When n is 20, the crossproduct is
> (numerically) a diagonal matrix with +-1 on the
> diagonal. When n is 50, this is not the case, but that
> could be an issue of nearly identical eigenvalues.
> Is there no way within R to require that the sequence of
> operations be the same for identical calls?
As Peter Dalgaard mentioned, you have the problem of underlying
system libraries that try to be fast and hence parallelize
computations, notably in this (and man similar calls), MKL uses
parallelized BLAS and/or LAPACK .. and that's what eigen() in R
(and Julia, python, matlab, ..) all use.
And parallelization is a killer of (strict / bit-level)
reproducibility, as you have just experienced.
If you know how to tell your OS / that you want only one
parallel "thread" / process / ...
you get back to reproducible (and slower) computations.
> The problem arose originally in a package test in which I wanted to
> verify that two ways of specifying something led to the
> execution of exactly the same calculations. The best proxy
> for this seemed to be to use identical() on the outputs,
> but if the same line of code (that should not involve an
> RNG) can lead to different results, this approach is
> doomed, yes? It is not absolutely critical that the
> outputs be identical(), but it would be much more
> reassuring than all.equal().
I understand and agree.
When I first became aware of the irreprodicibility of parallel
computations, only a few years ago, I was quite shocked and deillusionized..
> Thanks, Jesse
Martin Maechler
ETH Zurich and R Core team
> On Sunday, May 4th, 2025 at 12:27 PM, peter dalgaard
> <pda...@gmail.com> wrote:
>>
>>
>> Have you looked more closely at the differences?
>> Eigenvectors are only determined up to a sign change, so
>> different platforms often give results that differ by
>> sign. If you use a multitasking numerical library, the
>> same can happen within platform because the exact
>> sequence of computations differs between calls.
>>
>> You could check
>>
>> a) that e1$values and e2$values are the same b) that the
>> crossprod(e1$vectors, e2$vectors) is a diagonal matrix
>> with 1 or -1 in the diagonal. (This might fail if you
>> have eigenvalues that are almost identical, though.)
>>
>> -pd
>>
>> > On 4 May 2025, at 17.57, smallepsilon via R-help
>> r-help@r-project.org wrote:
>> >
>> > I am using MKL with R 4.5.0 on Linux, and eigen() is
>> producing different results with identical
>> calls. Specifically, when I run the code below, I get
>> "FALSE" from identical(). Just in case it had something
>> to do with a random number generator, I put identical
>> set.seed() calls immediately before each eigen() call,
>> but that did not help. Has anyone seen this behavior
>> before?
>> >
>> > (When n is 5, the identical() call almost always
>> returns "TRUE". As n increases, the proportion of FALSE
>> results increases, and it is nearly always FALSE when n
>> is 50.)
>> >
>> > Jesse
>> >
>> > ***
>> >
>> > n <- 50 > set.seed(20250504) > Sigma <- rWishart(1, df
>> = n, Sigma = diag(n))[,,1] > e1 <- eigen(Sigma) > e2 <-
>> eigen(Sigma) > identical(e1, e2)
>> >
>> > ______________________________________________ >
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and
>> more, see > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide
>> https://www.R-project.org/posting-guide.html > and
>> provide commented, minimal, self-contained, reproducible
>> code.
>>
>>
>> --
>> Peter Dalgaard, Professor, Center for Statistics,
>> Copenhagen Business SchoolSolbjerg Plads 3, 2000
>> Frederiksberg, Denmark Phone: (+45)38153501 Office: A
>> 4.23 Email: pd....@cbs.dk Priv: pda...@gmail.com
> ______________________________________________
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and
> more, see https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> https://www.R-project.org/posting-guide.html and provide
> commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide https://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
