I?m running R 3.2.2 on a Linux server (Redhat 4.4.7-16), and having the
following problem.
It works fine with the following:
require('MASS?)
var(mvrnorm(n = 1000, rep(0, 2), Sigma=matrix(c(10,3,3,2),2,2)))
However, when running the following in a loop with simulated data (Sigma):
# Sigma defined somewhere else
mvrnorm(n=1000, rep(0, 190), Sigma)
I get this opaque message:
*** caught illegal operation ***
address 0x7fe78f8693d2, cause 'illegal operand'
Traceback:
1: eigen(Sigma, symmetric = TRUE)
2: mvrnorm(n = nr, rep(0, NN), Sigma)
Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
I tried to do core dump (option 1), but it didn?t go anywhere (hanging
there forever). I also ran the same code on a Mac, and there was no problem
at all. What is causing the problem on the Linux server? In case the
variance-covariance matrix ?Sigma? is needed, I can provide its definition
later.
Thanks,
Gang
[[alternative HTML version deleted]]
I cobbled together a 190 x 190 positive definite matrix Sigma and ran
your example. I got a result instantaneously, with no error message.
(I'm running Linux; an ancient Fedora 17 system.)
So the problem is peculiar to your particular Sigma.
As the error message tells you, the problem comes from doing an
eigendecomposition of Sigma. So start your investigation by doing
E <- eigen(Sigma,symmetric=TRUE)
Presumably that will lead to the same error. How to get around this
error is beyond the scope of my capabilities.
You *might* get somewhere by using the singular value decomposition
(equivalent for a positive definite matrix) rather than the
eigendecomposition. I have the vague notion that the svd is more
numerically robust than eigendecomposition. However I might well have
that wrong.
Doing anything in 190 dimensions is bound to be fraught with numeric
peril.
cheers,
Rolf Turner
--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276
On 19/11/15 08:28, Gang Chen wrote:> I?m running R 3.2.2 on a Linux server (Redhat 4.4.7-16), and having the
> following problem.
>
> It works fine with the following:
>
> require('MASS?)
> var(mvrnorm(n = 1000, rep(0, 2), Sigma=matrix(c(10,3,3,2),2,2)))
>
> However, when running the following in a loop with simulated data (Sigma):
>
> # Sigma defined somewhere else
> mvrnorm(n=1000, rep(0, 190), Sigma)
>
> I get this opaque message:
>
> *** caught illegal operation ***
> address 0x7fe78f8693d2, cause 'illegal operand'
>
> Traceback:
> 1: eigen(Sigma, symmetric = TRUE)
> 2: mvrnorm(n = nr, rep(0, NN), Sigma)
>
> Possible actions:
> 1: abort (with core dump, if enabled)
> 2: normal R exit
> 3: exit R without saving workspace
> 4: exit R saving workspace
>
> I tried to do core dump (option 1), but it didn?t go anywhere (hanging
> there forever). I also ran the same code on a Mac, and there was no problem
> at all. What is causing the problem on the Linux server? In case the
> variance-covariance matrix ?Sigma? is needed, I can provide its definition
> later.
Thanks a lot for the pointer, Rolf!
You're correct that
E <- eigen(Sigma,symmetric=TRUE)
does lead to the same error on the RedHat machine. However, the same
computation on my Mac works fine:
E <- eigen(Sigma,symmetric=TRUE)
E$values
[1] 4.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6
2.6
[19] 2.6 2.6 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[37] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[55] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[73] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[91] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[109] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[127] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[145] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[163] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
0.8
[181] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
As you can see, all the eigenvalues are truly positive. Is it possible that
some numerical library or package is required by eigen() but missing on the
Linux server? In case the real Sigma is useful, here is how the matrix
Sigma is defined:
constrSigma <- function(n, sig) {
N <- n*(n-1)/2
Sigma <- diag(N)
bgn <- 1
for(ii in (n-1):1) {
end <- bgn+(ii-1)
Sigma[bgn:end,bgn:end][lower.tri(Sigma[bgn:end,bgn:end])] <- rep(sig,
ii*(ii-1)/2)
if(ii<(n-1)) {
xx <- cbind(rep(sig,ii), diag(ii)*sig)
yy <- xx
for(jj in 1:(n-1-ii)) {
if(jj>1) {
xx <- cbind(rep(0, ii), xx)
yy <- cbind(xx, yy)
}
}
Sigma[bgn:end,1:(bgn-1)] <- yy
}
bgn <- end+1
}
Sigma[upper.tri(Sigma)] <- t(Sigma)[upper.tri(t(Sigma))]
return(Sigma)
}
Sigma <- constrSigma(20, 0.1)
mvrnorm(n=1000, rep(0, 190), Sigma)
On Wed, Nov 18, 2015 at 3:56 PM, Rolf Turner <r.turner at auckland.ac.nz>
wrote:>
>
> I cobbled together a 190 x 190 positive definite matrix Sigma and ran
your example. I got a result instantaneously, with no error message. (I'm
running Linux; an ancient Fedora 17 system.)>
> So the problem is peculiar to your particular Sigma.
>
> As the error message tells you, the problem comes from doing an
eigendecomposition of Sigma. So start your investigation by
doing>
> E <- eigen(Sigma,symmetric=TRUE)
>
> Presumably that will lead to the same error. How to get around this
error is beyond the scope of my capabilities.>
> You *might* get somewhere by using the singular value decomposition
> (equivalent for a positive definite matrix) rather than the
eigendecomposition. I have the vague notion that the svd is more
numerically robust than eigendecomposition. However I might well have that
wrong.>
> Doing anything in 190 dimensions is bound to be fraught with numeric
> peril.
>
> cheers,
>
> Rolf Turner
>
> --
> Technical Editor ANZJS
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>
>
> On 19/11/15 08:28, Gang Chen wrote:
>>
>> I?m running R 3.2.2 on a Linux server (Redhat 4.4.7-16), and having
the
>> following problem.
>>
>> It works fine with the following:
>>
>> require('MASS?)
>> var(mvrnorm(n = 1000, rep(0, 2), Sigma=matrix(c(10,3,3,2),2,2)))
>>
>> However, when running the following in a loop with simulated data
(Sigma):>>
>> # Sigma defined somewhere else
>> mvrnorm(n=1000, rep(0, 190), Sigma)
>>
>> I get this opaque message:
>>
>> *** caught illegal operation ***
>> address 0x7fe78f8693d2, cause 'illegal operand'
>>
>> Traceback:
>> 1: eigen(Sigma, symmetric = TRUE)
>> 2: mvrnorm(n = nr, rep(0, NN), Sigma)
>>
>> Possible actions:
>> 1: abort (with core dump, if enabled)
>> 2: normal R exit
>> 3: exit R without saving workspace
>> 4: exit R saving workspace
>>
>> I tried to do core dump (option 1), but it didn?t go anywhere (hanging
>> there forever). I also ran the same code on a Mac, and there was no
problem>> at all. What is causing the problem on the Linux server? In case the
>> variance-covariance matrix ?Sigma? is needed, I can provide its
definition>> later.
>
>
>
[[alternative HTML version deleted]]
I have a sense of a tiny bell ringing faintly from the distant past... There may be an issue with some versions of BLAS/LAPACK on some systems showing up in eigendecompositions. Checking... hmm, there have been several issues over time but I don't see anything that would be directly related to the present issue. There was a really stupid bug in R's reference BLAS, but that was found and fixed two years ago, and besides, it involved complex eigenvalues which I don't see creeping into mvrnorm. Also PR#15211 located an issue that could cause a hang (not segfault) which ended in a WONTFIX situation because it really needed to be fixed in LAPACK. Anyways, some details about versions, hardware (you may be using a library optimized for the wrong CPU), would be needed. Even better if you can set up to run under Valgrind or a debugger to pinpoint the cause of the trouble. Also, check that you aren't simply running out of memory. Pragmatically speaking, couldn't you get away with using a Choleski decomposition of Sigma? -pd> On 18 Nov 2015, at 21:56 , Rolf Turner <r.turner at auckland.ac.nz> wrote: > > > I cobbled together a 190 x 190 positive definite matrix Sigma and ran your example. I got a result instantaneously, with no error message. (I'm running Linux; an ancient Fedora 17 system.) > > So the problem is peculiar to your particular Sigma. > > As the error message tells you, the problem comes from doing an eigendecomposition of Sigma. So start your investigation by doing > > E <- eigen(Sigma,symmetric=TRUE) > > Presumably that will lead to the same error. How to get around this error is beyond the scope of my capabilities. > > You *might* get somewhere by using the singular value decomposition > (equivalent for a positive definite matrix) rather than the eigendecomposition. I have the vague notion that the svd is more numerically robust than eigendecomposition. However I might well have that wrong. > > Doing anything in 190 dimensions is bound to be fraught with numeric > peril. > > cheers, > > Rolf Turner > > -- > Technical Editor ANZJS > Department of Statistics > University of Auckland > Phone: +64-9-373-7599 ext. 88276 > > On 19/11/15 08:28, Gang Chen wrote: >> I?m running R 3.2.2 on a Linux server (Redhat 4.4.7-16), and having the >> following problem. >> >> It works fine with the following: >> >> require('MASS?) >> var(mvrnorm(n = 1000, rep(0, 2), Sigma=matrix(c(10,3,3,2),2,2))) >> >> However, when running the following in a loop with simulated data (Sigma): >> >> # Sigma defined somewhere else >> mvrnorm(n=1000, rep(0, 190), Sigma) >> >> I get this opaque message: >> >> *** caught illegal operation *** >> address 0x7fe78f8693d2, cause 'illegal operand' >> >> Traceback: >> 1: eigen(Sigma, symmetric = TRUE) >> 2: mvrnorm(n = nr, rep(0, NN), Sigma) >> >> Possible actions: >> 1: abort (with core dump, if enabled) >> 2: normal R exit >> 3: exit R without saving workspace >> 4: exit R saving workspace >> >> I tried to do core dump (option 1), but it didn?t go anywhere (hanging >> there forever). I also ran the same code on a Mac, and there was no problem >> at all. What is causing the problem on the Linux server? In case the >> variance-covariance matrix ?Sigma? is needed, I can provide its definition >> later. > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.-- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com