eigenvalues of a circulant matrix
I think we need to make clear that eigen() by default relies on LAPACK routines and they in turn rely on BLAS routines. We have seen several instances in which LAPACK/BLAS return NaNs when they should not, but all that I am aware of are when (broken) external libraries were used. So the occurrence of NaNs should lead you to question other aspects of your computational environment. If you have such a broken environment, calling eigen(EISPACK=TRUE) may be a palliative, but it is better to track the problem down.
On Tue, 3 May 2005 Ted.Harding at nessie.mcc.ac.uk wrote:
I think the statement of the problem and the questions asked need clarifying. Some aspects puzzleme. See below. On 03-May-05 Globe Trotter wrote:
Looks like the files did not go through again. In any case, here is the kinv: please cut and paste and save to a file: [data snipped] --- Globe Trotter <itsme_410 at yahoo.com> wrote:
Date: Mon, 2 May 2005 19:51:24 -0700 (PDT)
From: Globe Trotter <itsme_410 at yahoo.com>
Subject: Re: [R] eigenvalues of a circulant matrix
To: r-help at stat.math.ethz.ch
OK, here we go:
I am submitting two attachments. The first is the datafile
called kinv used to create my circulant matrix, using the
following commands:
x<-scan("kinv")
y<-x[c(109:1,0:108)]
X=toeplitz(y)
eigen(X)
write(X,ncol=216,file="test.dat")
Having cut&pasted from the data placed in the body of the
message (omitted here) I get 216 numbers. Having put these
in a vector x (in my own way):
length(x)
##[1] 216
Question 1:
===========
Is this correct? Or has there been a problem with your
posting of the data?
If it is correct, given that you seem to only use x[1:109],
was there some point in giving the rest?
Question 2:
===========
Next, using your command:
y<-x[c(109:1,0:108)]
I now get
length(y)
##[1] 217
(as expected). The "0" in "0:108" seems to have been ignored
(again as expected), so this is equivalent to
y<-x[c(109:1,1:108)]
Is this as intended? If so, why use "0:108" instead of "1:108"?
Check:
y[1] ##[1] 19.4495
x[109] ##[1] 19.4495
y[109] ##[1] -0.00116801
x[1] ##[1] -0.00116801
y[110] ##[1] -0.00116801
x[1] ##[1] -0.00116801
y[217] ##[1] -6.28085
x[108] ##[1] -6.28085
Can you confirm that this is as intended?
Comment 3:
==========
You next command X=toeplitz(y): No apparent problems,
it gives a symmetric result:
which(X != t(X)) ## numeric(0)
with 217 rows and columns:
dim(X) ##[1] 217 217
and looks circulant:
X[(1:5),(1:5)]
[,1] [,2] [,3] [,4] [,5]
[1,] 19.449500 -6.280850 -0.486405 -0.826079 -0.167792
[2,] -6.280850 19.449500 -6.280850 -0.486405 -0.826079
[3,] -0.486405 -6.280850 19.449500 -6.280850 -0.486405
[4,] -0.826079 -0.486405 -6.280850 19.449500 -6.280850
[5,] -0.167792 -0.826079 -0.486405 -6.280850 19.449500
Question 4:
===========
Your next command, "eigen(X)", would simply output the results
to screen and does not assign to anything.
Your next command "write(X,ncol=216,file="test.dat")" as it
stands will write the toeplitz matrix X, constructed by
your command "X<-toeplitz(y)" to file, but with 216
columns instead of 217. However, the result consists
simply of numbers, and there is nothing like "NA" or "NaN"
in the file which I get.
Nor are there any NAs or NaNs in X itself, of course.
But, when you yourself did "write(X,ncol=216,file="test.dat")",
perhaps the "X" in this command was different from the "X"
which is the toeplitz matrix. So, was it the result of an
assignment from "eigen(X)" and, if so, which component or
components?
Question/Comment 5:
===================
So I have tried Z<-eigen(X). First of all, I get no problems
with NAs or NaNs:
which(is.na(Z$values)) ##numeric(0)
which(is.nan(Z$values)) ##numeric(0)
which(is.na(Z$vectors)) ##numeric(0)
which(is.nan(Z$vectors)) ##numeric(0)
Next, trying various options for wirting to file:
write(Z,ncol=216,file="test.dat")
simply does not work (not a writable structure), while
write(Z$values,ncol=216,file="test.dat")
produces simply a set of numbers, no NAs of NaNs, and likewise
write(Z$vectors,ncol=216,file="test.dat")
(the only occurrences of non-numeric characters are "e", as
in "e-05").
reports the following columns full of NaN's: 18, 58, 194, 200. (Note that eigen(X,symmetric=T) makes no difference and I get the same as above).
Therefore I find myself completely unable to reproduce your problem. However, for the various reasons stated in detail above, I am not at all sure that what you wrote as the statement of what you did in fact corresponds to what you really did! I even wonder whether Question 6: =========== Was the file "test.dat" the result of your "write" command? Or was it left over from a previous activity, the "write" from this session having failed to execute for some reason? (In which case the NaNs would have nothing to do with the results of "eigen(X)").
The second attachment contains only the eigenvectors obtained on calling a LAPACK routine directly (from C). The eigenvalues are essentially the same as that obtained using R. Here, I use the LAPACK-recommended double precision routine dspevd() routine for symmetric matrices in packed storage format. Note the absence of the NaN's....I would be happy to send my C programs to whoever is interested.
Well, I didn't get any NaNs in R either -- quite consistent with your C program! Please clarify according to the questions above. Best wishes, Ted. -------------------------------------------------------------------- E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 094 0861 Date: 03-May-05 Time: 12:06:57 ------------------------------ XFMail ------------------------------
______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595