`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).


#8

Agree, we also do not have Integer <: Rational and Real <: Complex, even though that is the case in math. So ship of consistency with math have sailed long before :slight_smile:


#10

Please avoid necroposting unless the new comment adds significant value. Otherwise you are effectively pinging everyone who has participated in an old conversation just to put in your two cents.

I’ve gone into this before (somewhere), but the way rationals and complexes are defined and constructed in Julia is completely consistent with how they are constructed in mathematics. The integers are not a subset of the rationals, the rationals are defined in terms of equivalence classes of pairs of integers; the integers are then identified with a subset of the rationals. The reals are then defined in terms of equivalence classes of rationals (Cauchy sequences or Dedekind cuts are two of the common constructions). The rationals are not a subset of the reals: they are identified with a subset of the reals. Similarly, the complex numbers are defined in terms of the reals and the reals are then identified with a subset of the complexes.


#11

Ok.

This is just wrong. There are axiomatic definitions of each of these numbers (such as, for example, Tarski’s axiomatization of the reals), which imply that ℕ ⊂ ℤ ⊂ ℚ ⊂ ℝ ⊂ ℂ. The fact that you choose to define these numbers constructively cannot contradict the implications of axiomatic definitions because two definitions are equivalent.

And even at the very beginning, with you can define it via Peano axioms or provide set-theoretic construction for them via the successor function. But the latter definition satisfies Peano axioms and hence cannot contradict it’s implications. So yes, it is true that integers are rationals and reals are complex numbers.

Look I understand why it was necessary to set the Julia’s type hierarchy the way it is right now. And indeed it probably is the best alternative of all other options (each of which has it’s own drawbacks). But we shouldn’t now pretend that it is consistent with hierarchy given by math and bend the mathematical definitions for the sake of it.


#12

Phrases such as “this is just wrong” and “we shouldn’t now pretend…” seem somewhat strong given that both Azamat’s and Stefan’s view points find support within the community of mathematicians. For example Stefan’s view point is also described by mathematics professor Ron Freiwald in “constructing the rationals”.

I’d like to observe that <: (is subtype of) does not always imply (is subset of). Or to quote Bill Clinton out of context: “It depends upon what the meaning of the word ‘is’ is”.


#13

I don’t think anyone was suggesting that number types used by computers match the mathematical definitions exactly. If you want to be very picky about this, then of course has nothing to do with Int, since the former is infinite while the latter isn’t, and surely isn’t related to <:AbstractFloat at all whatsoever. The very idea!

We could make this clear in all discussions, but that would become cumbersome. Hopefully most computer programmers understand these distinctions implicitly, and being pedantic about them is just a distraction.


#14

Sure, direct axiomatizations of the reals and complex numbers exist. You could also just construct ℂ and then forget all the constructions you used to get there and define the “real integers” as being the appropriate subset of ℂ. However, that doesn’t change the situation: in computers, unlike axiomatic mathematics, we have to construct the numbers that we use, not just give some axioms that they follow. After all, how would go from axiomatization of ℂ to an implementation of ℂ? Therefore the relationship between Julian integers, rationals, reals and complexes resemble the traditional mathematical construction rather than an axiomatic definition.