Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Special support for roots when the coefficients are real? #19

Closed
nsajko opened this issue Oct 8, 2023 · 2 comments · Fixed by #20
Closed

Special support for roots when the coefficients are real? #19

nsajko opened this issue Oct 8, 2023 · 2 comments · Fixed by #20

Comments

@nsajko
Copy link
Contributor

nsajko commented Oct 8, 2023

Empirically, it seems that roots(::Vector{Float64}) returns the roots such that they exactly conform to the form that follows from the complex conjugate root theorem. That is, a real root has its imaginary part as exactly zero, and complex roots come in complex conjugate pairs, where the complex conjugate operation is exact. Example:

julia> import AMRVW

julia> const A = AMRVW
AMRVW

julia> A.roots(rand(20))
19-element Vector{ComplexF64}:
  -1.2606264533481746 + 0.0im
  -0.9310479653026866 - 0.48344571715462364im
  -0.9310479653026866 + 0.48344571715462364im
  -0.7452679583988157 + 0.0im
 -0.48547663249542244 - 0.8197940002283542im
 -0.48547663249542244 + 0.8197940002283542im
 -0.30132852172591407 - 1.010361296238549im
 -0.30132852172591407 + 1.010361296238549im
  -0.2338481178786882 - 1.3556065317036028im
  -0.2338481178786882 + 1.3556065317036028im
  -0.1383958798698237 + 0.0im
  0.34088495469868396 - 1.0901157173978522im
  0.34088495469868396 + 1.0901157173978522im
  0.34435589931079486 - 0.8281753378613422im
  0.34435589931079486 + 0.8281753378613422im
    0.754774828425891 - 0.5684903527552175im
    0.754774828425891 + 0.5684903527552175im
   0.9427060520193593 - 0.3181694681365702im
   0.9427060520193593 + 0.3181694681365702im

I wonder if (or to what degree) is this guaranteed or (alternatively) an accident?

If it is guaranteed, I think it would make sense if this package exposed a new interface that would return the roots in a less redundant and more convenient form, as expected from the complex conjugate root theorem. So, instead of returning a Vector{Complex{<:AbstractFloat}}, this new interface would return an instance of a type like this, where the complex conjugate root pairs are stored as only the root whose imaginary part is positive:

struct ComplexConjugateRootTheoremRoots{F<:AbstractFloat}
  r::Vector{F}

  # The imaginary part of each element is guaranteed to be positive.
  c::Vector{Complex{F}}
end

How much sense am I making?

@jverzani
Copy link
Owner

jverzani commented Oct 8, 2023

It is guaranteed, the algorithm reduces the matrix until there is a linear or quadratic to evaluate. I'm not sure what is more convenient here though. Without thinking too much, this seems a bit complicated, but you could play around with it.

The key lines in the https://github.com/jverzani/AMRVW.jl/blob/master/src/AMRVW_algorithm.jl file are 50, 79, 102, and 112. Unless I missed something, here is where EIGS is written to.

@nsajko
Copy link
Contributor Author

nsajko commented Oct 8, 2023

Without thinking too much, this seems a bit complicated, but you could play around with it.

The key lines in the https://github.com/jverzani/AMRVW.jl/blob/master/src/AMRVW_algorithm.jl file are 50, 79, 102, and 112. Unless I missed something, here is where EIGS is written to.

The linked PR implements the feature, and it's actually trivial. I suppose your comment discusses a hypothetical lower-level implementation that would be more efficient? Still, I suppose my proof-of-concept PR is good for a start.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants