[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