DualNumbers.jl
DualNumbers.jl copied to clipboard
Complex Differentiation
I'd really like to be able to differentiate complex functions with ForwardDiff. For almost all analytic functions, the derivative can be taken in the same way as in the real case, so Julia's dynamic type system should make complex diferentiation work out of the box. It seems to me that the only reason it doesn't is because many functions in ForwardDiff
accept Real
arguments, but if we simply changed these to accept Number
s then we should be most of the way there. Functions that are not analytic such as real
, imag
, abs
, conj
, abs2
, reim
should all return errors when their derivatives are taken at a complex input. For example, consider these two functions:
f(z) = z^2
g(z) = Complex(real(z)^2 - imag(z)^2, 2 * real(z) * imag(z))
These functions are equal and both analytic, but it is not obvious by looking at g
that it is analytic. In my opinion, it would be acceptable if derivative(f, im)
worked while derivative(g, im)
threw an error.
I already started experimenting with this. I copied the definitions of several functions in ForwardDiff
and just replaced Real
with Complex
, like so:
ForwardDiff.can_dual(::Type{Complex{T}}) where {T<:Real} = true
@inline function ForwardDiff.derivative(f::F, x::R) where {F,R<:Complex}
T = typeof(Tag(f, R))
return extract_derivative(T, f(Dual{T}(x, one(x))))
end
@inline function Base.:*(partials::Partials, x::Complex)
return Partials(scale_tuple(partials.values, x))
end
@inline function Base.:/(partials::Partials, x::Complex)
return Partials(div_tuple_by_scalar(partials.values, x))
end
@inline Base.:*(x::Complex, partials::Partials) = partials * x
@inline Base.:*(partials::Partials{0,V}, x::Complex) where {V} = Partials{0,promote_type(V, typeof(x))}(tuple())
@inline Base.:*(x::Complex, partials::Partials{0,V}) where {V} = Partials{0,promote_type(V, typeof(x))}(tuple())
@inline Base.:/(partials::Partials{0,V}, x::Complex) where {V} = Partials{0,promote_type(V, typeof(x))}(tuple())
This was enough for me to correctly differentiate sin
. I expect that if we just do this throughout ForwardDiff
(and replace the Real
definitions with Number
, instead of adding Complex
definitions like I did here), then we'll be most of the way there. Here's my plan:
- Find and replace
Real
withNumber
throughout the repository - Tweak things until all the test cases pass
- Add test cases for complex functions
I'm tempted to actually do this and make a pull request as soon as I have some time to do so.
Just use a Complex{<:Dual}
.
Note DualNumbers.jl takes the other approach you advocate but is planned to be redesigned as ForwardDiff.Dual
What could be added is support for derivative(f, im)
but that doesn't require a redesign of Dual
.
Thanks for the response. I tried that solution and there were some problems.
Complex{Dual}
worked for some functions like polynomials, but didn't work for functions like sin
. If we look at the source code for sin(::Complex)
, which gets called when we try sin(::Complex{Dual})
:
function sin(z::Complex{T}) where T
F = float(T)
zr, zi = reim(z)
if zr == 0
Complex(F(zr), sinh(zi))
elseif !isfinite(zr)
if zi == 0 || isinf(zi)
Complex(F(NaN), F(zi))
else
Complex(F(NaN), F(NaN))
end
else
s, c = sincos(zr)
Complex(s * cosh(zi), c * sinh(zi))
end
end
So if we call derivative(sin, im)
(after overloading derivative
and extract_derivative
to use Complex{Dual}
) we end up in the first code path because zr === Dual(0, 1) == 0
. However, this gives the wrong answer because zr != sin(zr) * cosh(zi)
. These expressions have the same real part, but different infinitesimal parts.
Simply passing a Dual
into the implementation of sin(::Real)
from Base
wouldn't work for similarish reasons, which is why ForwardDiff
provides a method for sin(::Dual)
function Base.sin(d::Dual{T}) where T
s, c = sincos(value(d))
return Dual{T}(s, c * partials(d))
end
If we wanted to use Complex{Dual}
, we would have to make this method and methods like it take both Dual
and Complex{Dual}
. Alternatively, we could use Dual{Complex}
. It seems to me that the trade off is between replacing Dual
with Union{Dual, Complex{Dual}}
throughout the codebase or replacing Real
with Number
. The second option seems preferable, but the first might work too with some trickery.
As for DualNumbers.jl
, that is an option, but it doesn't have a derivative
interface to hide the dual numbers, and it isn't actively maintained. Maybe I could look into using TaylorSeries.jl
. In any case I think ForwardDiff
would benefit from this feature.
What could be added is support for
derivative(f, im)
but that doesn't require a redesign ofDual
.
To be clear, I tried this when I was experimenting with Complex{Dual}
, and ran into the problems described above.
Changing Dual
s supertype to Number
will be incredibly disruptive and break a lot of code. I think there is no chance a PR doing this will be merged so I would recommend not putting time into that.
I don't think you would have to add Union
s everywhere, only special cases where the default breaks (such as sin
). For example, exp
already works:
julia> exp(Complex(Dual(0,1),Dual(1,0)))
Dual{Nothing}(0.5403023058681398,0.5403023058681398) + Dual{Nothing}(0.8414709848078965,0.8414709848078965)*im
Changing
Dual
s supertype toNumber
will be incredibly disruptive and break a lot of code. I think there is no chance a PR doing this will be merged so I would recommend not putting time into that.
Okay, fair enough. Also, I just noticed that I posted this issue on DualNumbers
and not ForwardDiff
, my bad 😬
Taking a closer look at ForwardDiff
, we wouldn't actually have to change that many places to take Complex{<:Dual}
in addition to Dual
. Most of the special function definitions for Dual
s are done with unary_dual_definition
and binary_dual_definition
, so we just have to change those and a couple other places. I might give that a go.