JuliaDiff / ForwardDiff.jl

Forward Mode Automatic Differentiation for Julia
Other
887 stars 141 forks source link

LinAlg.normalize fails on Vector{ForwardDiff.Dual} #175

Open tpapp opened 7 years ago

tpapp commented 7 years ago

example:

import ForwardDiff: Dual

v = ForwardDiff.Dual{6,Float64}[Dual(-0.0927781,-0.0679932,0.0650883,-0.0217832,0.00395549,0.0415781,0.128862),Dual(0.370172,0.0333549,0.0891681,0.0439229,0.0209036,-0.011965,0.0331594),Dual(-0.456713,0.0373191,0.0385198,0.0652645,0.00863201,-0.0826074,-0.20347),Dual(-0.46295,-0.00857001,-0.116858,-0.0657061,-0.0169481,0.0535432,0.0957022),Dual(0.656853,-0.00849306,-0.0966359,-0.0287605,-0.0171647,-0.00708444,-0.0745087)]

normalize(v, 1)

fails with

ERROR: MethodError: no method matching __normalize!(::Array{ForwardDiff.Dual{6,Float64},1}, ::ForwardDiff.Dual{6,Float64})

I isolated the problem and would be willing to make a PR, but I am unsure what an elegant fix would be. Possibly a method __normalize!{N, T <: AbstractFloat}(v::AbstractVector, nrm::ForwardDiff.Dual{N,T}) that does the same thing, branching on the size of nrm, but using typemax(T) instead here?

On a related note, I don't understand why Dual does not restrict T to AbstractFloat (I see little practical utility for AD using other types), but I am new to AD so perhaps I need to study the library more.

jrevels commented 7 years ago

Thanks for bringing this up!

I isolated the problem and would be willing to make a PR, but I am unsure what an elegant fix would be.

The most elegant fix, IMO, would be to implement Base.prevfloat(n::Dual) = prevfloat(value(n)) and open a PR in Base to change the type restriction here to Real. I haven't tried it yet, but I believe it should work after that.

Otherwise, ForwardDiff (or Base) could implement a generic normalize which does the naive thing (x ./ norm(x)), but then Dual numbers wouldn't get the benefit of Base's special implementation that tries to avoid overflow/underflow for floating point numbers.

On a related note, I don't understand why Dual does not restrict T to AbstractFloat (I see little practical utility for AD using other types), but I am new to AD so perhaps I need to study the library more.

One example is composing interval types with AD, which are also (I believe) usually subtyped to Real.