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

Add rules for evalpoly #190

Merged
merged 33 commits into from
Jun 30, 2020
Merged

Add rules for evalpoly #190

merged 33 commits into from
Jun 30, 2020

Conversation

sethaxen
Copy link
Member

@sethaxen sethaxen commented May 8, 2020

Most of the pushforward and pullback for evalpoly can be written as more calls to evalpoly, which is very efficient.

Edit: This is true in the scalar case but not the matrix. In the current implementation, both the matrix and scalar polynomials have forward and reverse rules on the same order as evalpoly itself.

Currently all but one of the tests fail for various reasons. I think the testers in ChainRulesTestUtils don't like when one of the inputs is a Tuple. Also, I don't expect the complex tests to pass until this package uses FiniteDifferences v0.10.0.

Edit: Oh, and I see that ChainRules doesn't test against 1.4 yet at all.

@oxinabox I also tested against FluxML/Zygote.jl#366 and may have hit a bug where Zygote doesn't like that the adjoint of p is a Tuple, but I'm not sure.

Edit: Now that I'm using Composite, it works fine in Zygote.

julia> using Zygote

julia> x = randn();

julia> pv = randn(10);

julia> p = Tuple(pv);

julia> evalpoly(x, pv)
126.48423475727726

julia> evalpoly(x, p)
126.48423475727726

julia> Zygote.gradient(evalpoly, x, pv) # seems to work fine
(1722.5531441128583, [1.0, 1.7894268524990125, 3.2020484604445225, 5.72983149812255, 10.253114343035136, 18.347198127169843, 32.83096899687731, 58.748617516574825, 105.12635373135284, 188.1159202721925])

julia> Zygote.gradient(evalpoly, x, p) # I guess not happy about the adjoint of p being a tuple?
ERROR: MethodError: no method matching conj(::NTuple{10,Float64})
Closest candidates are:
  conj(::Missing) at missing.jl:100
  conj(::ForwardDiff.Dual) at /Users/saxen/.julia/packages/ForwardDiff/cXTw0/src/dual.jl:367
  conj(::Real) at number.jl:167
  ...
Stacktrace:
 [1] wrap_chainrules_output at /Users/saxen/projects/Zygote.jl/src/compiler/chainrules.jl:49 [inlined]
 [2] map(::typeof(Zygote.wrap_chainrules_output), ::Tuple{ChainRulesCore.Zero,ChainRulesCore.Thunk{ChainRules.var"#51#56"{Float64,Float64,NTuple{10,Float64}}},ChainRulesCore.Thunk{ChainRules.var"#52#57"{Float64,Float64,NTuple{10,Float64}}}}) at ./tuple.jl:159
 [3] wrap_chainrules_output at /Users/saxen/projects/Zygote.jl/src/compiler/chainrules.jl:50 [inlined]
 [4] wrap_chainrules_pullback at /Users/saxen/projects/Zygote.jl/src/compiler/chainrules.jl:86 [inlined]
 [5] ZBack at /Users/saxen/projects/Zygote.jl/src/compiler/chainrules.jl:99 [inlined]
 [6] (::Zygote.var"#41#42"{Zygote.ZBack{ChainRules.var"#evalpoly_pullback#55"{Float64,NTuple{10,Float64}}}})(::Float64) at /Users/saxen/projects/Zygote.jl/src/compiler/interface.jl:45
 [7] gradient(::Function, ::Float64, ::Vararg{Any,N} where N) at /Users/saxen/projects/Zygote.jl/src/compiler/interface.jl:54
 [8] top-level scope at REPL[8]:1

@sethaxen
Copy link
Member Author

sethaxen commented May 10, 2020

This PR is basically feature complete, and tests pass. It's still marked as WIP because:

  • Pending clarification in Complex numbers ChainRulesCore.jl#159, the adjoints may need to be made transposes.
  • Complex tests are currently disabled until FD v0.10.0 is supported. They'll also depend on Unify the two rrule_test functions ChainRulesTestUtils.jl#33, since we use Horner's method for the polynomial, not the method that base uses, so they'll be approximately equal but not exactly. I don't think this should hold up this PR.
  • Matrix polynomials where either x or elements of p are UniformScaling are not supported. There are currently no rules for UniformScaling, so I don't think that should hold up this PR.

@sethaxen sethaxen changed the title [WIP] Add rules for evalpoly Add rules for evalpoly Jun 30, 2020
Copy link
Member

@oxinabox oxinabox left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should now be good to go, the tests look reasonable so I am happy to see this merged once they pass

@sethaxen
Copy link
Member Author

Yeah, the weird muladd in a macro giving an approximately equal but not equal result as muladd not in a macro on Travis issue came up again. No idea why. But tests pass now. Thanks!

@sethaxen sethaxen merged commit 224e553 into JuliaDiff:master Jun 30, 2020
@sethaxen sethaxen deleted the evalpoly branch June 30, 2020 23:56
@manuelbb-upb
Copy link

I just found a question on discourse asking why Zygote won't differentiate a polynomial created with Polynomials.jl.
I hope this is the right place to comment: I guess it could be due to the custom rules for evalpoly.
Take the example code from the original question:

using Zygote
using Polynomials
p=Polynomial([1.0,0.0,3.0,4.0])
@show derivative(p)(5.0)
@show gradient(p,5.0)

This errors because we need a similar method for p in _evalpoly_intermediate(x, p). No problem:

Base.similar( :: Polynomial, etype = eltype(p), len = length(p.coeffs)) = similar( p.coeffs, etype, len )

Now Zygote.gradient can calculate the gradient, but gives the wrong numerical result.
However,

gradient( evalpoly, 5.0, Tuple( p.coeffs ) )

gives the right result.

This might be due to the fact that the indexing for p starts at 0 and _evalpoly_intermediate(x,p), where p is not a Tuple, assumes indexing starts at 1.

p_tup = Tuple(p.coeffs)
p[0] == p_tup[1] # true
p[1] == p_tup[2] # true
# and so forth

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 this pull request may close these issues.

3 participants