[R] complex conjugates roots from polyroot?
Ravi Varadhan
rvaradhan at jhmi.edu
Mon Nov 26 16:31:22 CET 2007
Hi Spencer,
The default tolerance in your function might be a bit too conservative (i.e.
too small). Here is an example:
> set.seed(123)
> z <- polyroot(sample(1:5, size=40, rep=T)) # polynomial of degree 39
> zmat <- which(outer(z, z, FUN="is.conj"), arr.ind=T)
> zmat[zmat[,1] < zmat[,2], ]
row col
[1,] 3 5
[2,] 2 7
[3,] 6 10
[4,] 9 11
[5,] 13 15
[6,] 16 18
[7,] 1 22
[8,] 21 24
[9,] 17 25
[10,] 19 26
[11,] 23 29
[12,] 12 30
[13,] 20 32
[14,] 31 33
[15,] 8 34
[16,] 4 35
[17,] 28 36
[18,] 27 37
[19,] 38 39
> findConjugates(z)
[1] -0.0168724+0.8593702i -0.4707594-0.7991353i -0.8955664-0.2000173i
[4] 0.1495329-0.9612410i 0.3762447-0.8971735i 0.5654040-0.7332919i
[7] 0.6648998-0.6811790i -0.9087387+0.3750093i -0.5663986-0.8890497i
[10] 0.7937910-0.6069981i 0.3565431-0.9877840i
>
You can see that your function finds only 11 conjugate pairs as opposed to
the correct 19 pairs (since there is only 1 real root, the rest must be
conjugate pairs). Increasing the tolerance to 1.e-10, gives all the 19
pairs. Therefore, I would suggest a default of 1.e-10 or even larger, say,
sqrt(.Machine$double.eps).
I wouldn't consider multiple real roots to be conjugates, since they are not
distinct points on the complex plane, as well as for the reason that you
have given.
Best,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: Spencer Graves [mailto:spencer.graves at pdf.com]
Sent: Sunday, November 25, 2007 12:36 PM
To: Ravi Varadhan
Cc: r-help at r-project.org
Subject: Re: [R] complex conjugates roots from polyroot?
Hi, Ravi:
Question: Are duplicate real numbers complex conjugates? For
some purposes, perhaps. However, for identifying (damped) periodicity
on a difference or differential equation, we must exclude repeated real
roots.
In any event, your suggestion inspired me to create a function
"findConjugates" (see below), which I've added to the "FinTS" package
currently under development on R-Forge. The source code is currently
available via "svn checkout
svn://svn.r-forge.r-project.org/svnroot/fints". In another day or so,
it should be available (with documentation) via
'install.packages("FinTS",repos="http://r-forge.r-project.org")'. In
another couple of months, it should appear on CRAN.
Thanks again for your suggestion.
Best Wishes,
Spencer
######################################
findConjugates <- function(x,
complex.eps=1000*.Machine[["double.neg.eps"]]){
##
## 1. compute normalization
##
if(length(x)<1)return(complex(0))
ax <- abs(x)
m2 <- outer(ax, ax, pmax)
##
## 2. Compute complex differences
##
c2 <- (abs(outer(x, Conj(x), "-") / m2) < complex.eps)
c2[m2==0] <- FALSE
c2 <- (c2 & lower.tri(c2))
##
## 3. Any differences exceed complex.eps?
##
if(any(c2)){
# check standard differences
d2 <- (abs(outer(x, x, "-") / m2) > complex.eps)
d2[m2==0] <- FALSE
#
cd2 <- (c2 & d2)
if(any(cd2)){
ic <- sort(unique(row(cd2)[cd2]))
return(x[ic])
}
}
complex(0)
}
######################################
Ravi Varadhan wrote:
> Hi Spencer,
>
> Here is a simple approach to detect conjugate pairs:
>
> is.conj <- function(z1, z2, tol=1.e-10) {
> # determine if two complex numbers are conjugates
> cond1 <- abs(Re(z1) - Re(z2)) < tol
> cond2 <- abs(Im(z1) + Im(z2)) < tol
> cond1 & cond2
> }
>
> set.seed(123)
> z <- polyroot(sample(1:5, size=8, rep=T))
> zmat <- which(outer(z, z, FUN="is.conj"), arr.ind=T)
> zmat[zmat[,1] < zmat[,2], ]
>
> # result
> row col
> [1,] 1 3
> [2,] 5 6
> [3,] 4 7
>
>
> We see that (1,3), (4,7), and (5,6) are the conjugate pairs.
>
> This doesn't address the issue of numerical round-off (there is no
argument
> in polyroot that governs the accuracy of the roots).
>
> Best,
> Ravi.
>
>
----------------------------------------------------------------------------
> -------
>
> Ravi Varadhan, Ph.D.
>
> Assistant Professor, The Center on Aging and Health
>
> Division of Geriatric Medicine and Gerontology
>
> Johns Hopkins University
>
> Ph: (410) 502-2619
>
> Fax: (410) 614-9625
>
> Email: rvaradhan at jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>
----------------------------------------------------------------------------
> --------
>
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
> Behalf Of Spencer Graves
> Sent: Friday, November 23, 2007 12:08 PM
> To: r-help at r-project.org
> Subject: [R] complex conjugates roots from polyroot?
>
> Hi, All:
>
> Is there a simple way to detect complex conjugates in the roots
> returned by 'polyroot'? The obvious comparison of each root with the
> complex conjugate of the next sometimes produces roundoff error, and I
> don't know how to bound its magnitude:
>
> (tst <- polyroot(c(1, -.6, .4)))
> tst[-1]-Conj(tst[-2])
> [1] 3.108624e-15+2.22045e-16i
> abs(tst[-1]-Conj(tst[-2]))/abs(tst[-1])
> 1.971076e-15
> .Machine$double.neg.eps
> 1.110223e-16
>
> Testing (abs(tst[-1]-Conj(tst[-2]))/abs(tst[-1]) <
> (20*.Machine$double.neg.eps)) would catch this example, but it might not
> catch others.
>
> The 'polyroot' help page says, "See Also ... the 'zero' example in
> the demos directory." Unfortunately, I've so far been unable to find
> that.
>
> Any suggestions?
> Thanks in advance.
> Spencer Graves
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
>
More information about the R-help
mailing list