[R-pkg-devel] mvrnorm, eigen, tests, and R CMD check

[Kevin Coombes]

Hi,

I wrote and maintain the Thresher package. It includes code to do
simulations. In the "tests" directory of the package, I do some simple
simulations and run the main algorithm, then write out summaries of the
results

The initial submission of the package to CRAN was delayed because the
"Rout.save" files matched the "Rout" files on 64-bit R but *not*  on 32-bit
R on Windows. After investigating, I realized that when my simulation code
called "MASS::mvrnorm", I got different results from 64-bit and 32-bit
versions of R on the same machine.  Pushing further, I determined that this
was happening because mvrnorm used "eigen" to compute the eigenvalues and
eigenvectors, and "eigen" itself gave different answers in the two R
versions..

The underlying issue (mathematically) is that the correlation/covariance
matrix I was using had repeated eigenvalues, and so there is no unique
choice of basis for the associated eigenspace. This observation suggests
that the issue is potentially more general than 32-bit versus 64-bit; the
results will depend on the implementation of the eigen-decomposition in
whatever linear algebra module is compiled along with R, so it can change
from machine to machine.

I "solved" (well, worked around) the immediate problem with package
submission by changing the test code to not write out anything that might
differ between versions.

With all of that as background, here are my main questions:

[1] Is there any way to put something into the "tests" directory that would
allow me to use these simulations for what computer scientists call
regression testing? (That is, to make sure my changes to the code haven't
changed results in an unexpected way.)

[2] Should there be a flag or instruction to R CMD check that says to only
run or interpret this particular test on a specific version or machine? (Or
is there already such a flag that I don't know about?)

[3] Should the documentation (man page) for "eigen" or "mvrnorm" include a
warning that the results can change from machine to machine (or between
things like 32-bit and 64-bit R on the same machine) because of difference
in linear algebra modules? (Possibly including the statement that
"set.seed" won't save you.)

You can reproduce my example by running this code in different R
versions/machines:

sig <- matrix(0, 16, 16)
sig[1:10, 1:10] <- 0.5
sig[11:16,11:16] <- 0.5
diag(sig) <- 1
eigen(sig)

Best,
  Kevin

	[[alternative HTML version deleted]]