Use Leja ordering for fromroots
#586
Conversation
martinholters
force-pushed
the
mh/leja
branch
from
db1f724 to
c012b05
Compare
martinholters
force-pushed
the
mh/leja
branch
from
c012b05 to
4f4a80e
Compare
|
Thanks! This looks good. I guess a test would be good to catch any regressions. |
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## master #586 +/- ##
==========================================
- Coverage 75.69% 74.01% -1.69%
==========================================
Files 38 38
Lines 4040 3817 -223
==========================================
- Hits 3058 2825 -233
- Misses 982 992 +10 ☔ View full report in Codecov by Sentry. |
There is one? If this looks good otherwise, I'll look into how to (elegantly) work around the missing |
|
It looks good. Thanks again. |
...as the required `argmax` method has only been added in Julia v1.7
| if P !== FactoredPolynomial | ||
| @test fromroots(P, roots(p)) * p[end] ≈ p |
There was a problem hiding this comment.
For FactoredPolynomial, fromroots is trivially correct; the actual conversion happens when promoting to a common type with p. That goes through evalpoly which uses an arbitrary ordering of the roots. I wonder whether it would generally be preferable to use the Leja ordering in evalpoly here, in which case it would probably be worth it to replace the Dict in FactoredPolynomial with an OrderedDict and do the Leja ordering upon construction. No clue whether that's a good idea.
There was a problem hiding this comment.
Well that's a question for another day, but it seems like a reasonable thing to do. The extra package has negligible load time.
|
Quite some failing checks, but failures all look unrelated/pre-existing to me. |
| #ii = argmax(j -> abs(roots[j]), eachindex(roots)) | ||
| ii = firstindex(roots) | ||
| rmax = abs(roots[ii]) | ||
| for j in eachindex(roots) | ||
| rabs = abs(roots[j]) | ||
| if rabs > rmax | ||
| ii = j | ||
| rmax = rabs | ||
| end | ||
| end |
There was a problem hiding this comment.
Alternatively, Compat.jl provides that argmax method. Would depending on Compat be an option or should we keep this loop? (Same below.)
There was a problem hiding this comment.
I don't mind having some small copy of package code in this package for easier compatability. I'd keep the loop unless you have a strong preference.
Yes, seems like code coverage changed recently. |
|
I'm happy merging now, unless you wanted to adjust the FactoredPolynomial storage. Thanks again! |
martinholters
changed the title
fromrootsfromroots
I'm not, merging as is is fine with me. (The other loose end would be the handling of zeros with higher arity where I'm not sure the current implementation follows the paper as I find that part a bit ambiguous. However, even if wrong, that will just be sub-optimal, not a catastrophic failure, and likely still better than what we had before.)
Welcome! |
Improve accuracy of
fromrootsby utilizing Leja ordering, see https://doi.org/10.1023/A:1025555803588 and ref. JuliaDSP/DSP.jl#584.I'm not sure about code organization and whether this hooks up in all the right places.
Another issue I'm not quite sure about is the correct handling of roots with multiplicity larger than 1.
Finally, this adds substantially to the runtime, but IMHO, accuracy is usually more important than performance here.