[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@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.
Re: [R] complex conjugates roots from polyroot?
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: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Spencer Graves
Sent: Friday, November 23, 2007 12:08 PM
To: r-help@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@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@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.
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: [EMAIL PROTECTED]
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>
>
>
> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
> Behalf Of Spencer Graves
> Sent: Friday, November 23, 2007 12:08 PM
> To: r-help@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@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@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@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.
Re: [R] complex conjugates roots from polyroot?
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: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: Spencer Graves [mailto:[EMAIL PROTECTED]
Sent: Sunday, November 25, 2007 12:36 PM
To: Ravi Varadhan
Cc: r-help@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: [EMAIL PROTECTED]
>
Re: [R] complex conjugates roots from polyroot?
Spencer,
I just observed that the polynomial root calculation in the package,
"polynom", using the function solve() is more accurate than the polyroot()
function in the "base" package. Here is an example:
set.seed(1234)
p <- polynomial(sample(1:10, size=45, rep=T)) # degree 44
z <- solve(p)
findConjugates(z, complex.eps=.Machine$double.eps) # this identifies all 21
conjugate pairs
z1 <- polyroot(p)
findConjugates(z1, complex.eps=.Machine$double.eps) # this only identifies
only 3 conjugate pairs
As, I had mentioned earlier, I can't tell what tolerances are used by these
algorithms. However, it appears that "polynom" is more accurate.
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: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: Spencer Graves [mailto:[EMAIL PROTECTED]
Sent: Sunday, November 25, 2007 12:36 PM
To: Ravi Varadhan
Cc: r-help@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: [EMAIL PROTECTED]
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>
>
>
> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On
> Behalf Of Spencer Graves
> Sent: Friday, November 23, 2007 12:08 PM
> To: r-help@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 produce
Re: [R] complex conjugates roots from polyroot?
Hi, Ravi:
Thanks very much. Per your suggestion, I made
complex.eps=.Machine[["double.eps"]]), added your examples to the help
file. I also listed "author" as "Spencer Graves and Ravi Varadhan" in
the help file and uploaded the changes to R-forge, as I mentioned below.
Later, I will modify the code that computes the roots to use
roots(polynomial(...)), per your suggestion.
Best Wishes,
Spencer
Ravi Varadhan wrote:
> Spencer,
>
> I just observed that the polynomial root calculation in the package,
> "polynom", using the function solve() is more accurate than the polyroot()
> function in the "base" package. Here is an example:
>
> set.seed(1234)
> p <- polynomial(sample(1:10, size=45, rep=T)) # degree 44
> z <- solve(p)
> findConjugates(z, complex.eps=.Machine$double.eps) # this identifies all 21
> conjugate pairs
>
> z1 <- polyroot(p)
> findConjugates(z1, complex.eps=.Machine$double.eps) # this only identifies
> only 3 conjugate pairs
>
> As, I had mentioned earlier, I can't tell what tolerances are used by these
> algorithms. However, it appears that "polynom" is more accurate.
>
> 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: [EMAIL PROTECTED]
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> ----
>
>
> -Original Message-
> From: Spencer Graves [mailto:[EMAIL PROTECTED]
> Sent: Sunday, November 25, 2007 12:36 PM
> To: Ravi Varadhan
> Cc: r-help@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 go
