Fwd: Re: [R] eigenvalues of a circulant matrix
Prof Brian Ripley
ripley at stats.ox.ac.uk
Tue May 3 13:51:42 CEST 2005
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
More information about the R-help
mailing list