[R] complex conjugates roots from polyroot?

[Spencer Graves]

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.
>