`Complex{Dual}` or `Dual{Complex}`


#1

TLDR: I think Dual should be <:Real and that Complex{Dual} leads to more generic code than Dual{Complex}.

Hi all,

I am toying with mixed complex and dual numbers.

My use case is numerical integration of a nonlinear PDE (2D Navier-Stokes) and the variational equations, i.e. equations that describe the evolution of infinitesimal around the reference trajectory. The two equations can be marched simultaneously time by letting the solution be described using dual numbers, where the \epsilon component the perturbation. With julia type generic code, there is not much to be done to get the equations solved when the code for the main PDE is available. Quite neat!

I solve the PDE using a Fourier spectral technique, and thus the state variables are complex coefficients. FFTs are used to transform the solution to the physical domain, evaluate nonlinearity there, and then backtransform to the spectral domain.

When the solution is represented using dual numbers, what you really want to do is to perform transforms of the real and perturbation parts separately. This can be achieved using guru FFTW interface, which allows specifying the stride pattern and the number of to perform (it requires hacking into FFTW.jl, at the moment, but it can be done).

The main issue (an annoyance, really) I have is that, currently [1], dual numbers are defined like this

const ReComp = Union{Real,Complex}
immutable Dual{T<:ReComp} <: Number
    value::T
    epsilon::T
end

This definition does not play nice with the signature of methods that you would define in type-generic code, and requires (AFAICS), specialised signatures for dual types. For instance, what would be written naturally as

# this method signature would work for `T` begin a Float64 or a Dual{Float64} if Dual <: Real
function rfft!(U::AbstractVector{Complex{T}}, u::AbstractVector{T}) where {T<:Real}
    # transform u to U
    # if T < Dual, transform both value and epsilon parts separately, using 
    # FFTW guru interface with custom strides pattern
end

needs custom Dual number methods such as

function rfft!(U::AbstractVector{Dual{Complex{T}}}, u::AbstractVector{Dual{T}}) where {T<:Real}
    # do stuff
end

This, and many other use cases where you have method signature such as

foo(a::Complex{T}, b::T) where {T} = # ...

as well as in type definitions, e.g.

struct Foo{T<:Real, A::AbstraArray{T}}
    data::A
end

struct CFoo{T<:Complex, A::AbstraArray{T}}
    data::A
end

would be solved gracefully if dual numbers where defined as

struct Dual{T<:Real} <: Real
    value::T
    epsilon::T
end

i.e. they where subtypes of real. One would then use Complex{Dual{T}} for automatic differentiation of functions with complex numbers, or other purposes.

I would be interested in knowing if there are cases where the current Dual{Complex} is better than Complex{Dual}.

Best,

Davide

[1] https://github.com/JuliaDiff/DualNumbers.jl/pull/29


#2

See also https://github.com/JuliaDiff/DualNumbers.jl/issues/45

I think there is broadly agreement for this. The only argument against is that dual numbers are not “real numbers” in the mathematical sense. But I think in Julia Real correlates more with a real interface than a mathematical real number. (Also, NaN and Inf are not real numbers either.)


#3

Thanks, did not go that far in the Dual repo issues.

Dual <: Real might not be mathematically correct, but would make sense in many other practical situations.

Maybe also worth linking this FFTW.jl issue https://github.com/JuliaMath/FFTW.jl/issues/31, for exposing the full FFTW guru functionality.


#4

Obligatory link to ForwardDiff issue on this:

https://github.com/JuliaDiff/ForwardDiff.jl/issues/157

The OP was actually trying to solve a spectral discretization of a PDE as well, but we ran into the issue that the stiff solvers wouldn’t work because the differentiation tools (Calculus.jl and ForwardDiff.jl) don’t support complex numbers. In a few days we’ll have DiffEqDiffTools.jl up and running to all of the stiff solvers will be able to run on complex by using autodiff=false, but yes I would really like to see dual numbers compatible with complex sooner rather than later.

Edit: I see you know about that thread :stuck_out_tongue:.


#5

Slightly off topic. I have coded up low storage 3rd/4th order IMEX Runge-Kutta schemes as stiff solvers, specifically developed for Navier-Stokes applications. You can find them here https://github.com/OpenRotatingPlaneCouetteFlow/IMEXRKCB.jl.


#6

I agree, I want autodiff for complex numbers quite badly as well :slight_smile:

See also https://github.com/JuliaDiff/ForwardDiff.jl/issues/43, which is not hard in principle and would do what you want with a nice API and without modifications to your code (also provided FFTs are taught about dual numbers). That would be pretty awesome.


#7

My point was that the Real interface should be written down, in a similar way to the AbstractArray interface. Then it is safe for a type to be a subtype provided it follows the rules of the interface.

I think worrying about “mathetically correct” is impractical, unless the restriction on Complex{T<:Real} is changed to Complex{T} (as well as similar restrictions).