Dual numbers and Complex numbers

[goretkin]

Many functions of a complex argument are not differentiable, and there's some math here that I don't quite understand, but if I define the DFT as

function dft{T<:Number}(x::Vector{T})
    N = length(x)

    #T should not have any imaginary part 

    C = zeros(T,N)
    S = zeros(T,N)

    for k=0:(N-1)
        for n=0:(N-1)
            cs = exp(-im*2pi* k*n/N)
            c,s = real(cs), imag(cs)
            @inbounds C[k+1] += c*x[n+1]
            @inbounds S[k+1] += s*x[n+1]
        end
    end
    return C,S 
end

which just returns the real and imaginary parts of the transform separately,

and also define some function that maps, for example, a 1D time-series x to a Real

function spec(x)
    C,S = dft(x)
    h = floor(length(x) / 2)
    w = cat(1, [1.0:h ], zeros(length(x)-int(h)) )
    sum( (C.^2 + S.^2) .* w ) #some weighting on different frequencies
end

then I can take the derivative of spec with respect to the input timeseries

function grad{T<:Real}(fun::Function,x::Vector{T})
    N = length(x)
    G = zeros(T,N)
    for i = 1:N
        xd = [Dual(x[j],1.0*(j==i)) for j=1:N]
        G[i] = epsilon(fun(xd))
    end
    return G
end

julia> grad(spec,rand(10,))
10-element Array{Float64,1}:
 12.1547 
 19.9224 
 -3.22121
 18.8699 
 21.0265 
 -9.73494
  1.77734
 18.524  
 24.9879 
  9.14384

hopefully that example made sense. If I wanted to use the complex version of the DFT and suitably modify spec ((C.^2 + S.^2) would be taking the squared-magnitude of those complex numbers), then as far as Real numbers go, spec would return the same output for the same inputs. But I don't see how I would put a Dual number through that function.

[mlubin]

The library wasn't designed to mix dual numbers with complex numbers, and I don't think any thought has been put into it. If you can figure out a way to make it work, we can definitely look at it.

[papamarkou]

That's true, dual numbers was meant to be a neat short impelementation of dual numbers. Imo we shouldn't be mixing dual and imaginary numbers in the scope of the present package. Mathematically speaking, both dual and imaginary numbers are special cases of the more general algebra of hypercomplex numbers. If there is an interest on merging treatment of complex and dual numbers, it may be a good idea, yet it would take some thinking.

[goretkin]

Yes, that makes perfect sense. I think there are mathematically unsound things to do in my specific example, which would work in my case but seems like it could not work in general. If Complex{Number} or Dual{Number} were possible instead of just {Real}, you could construct Dual{Complex} or Complex{Dual} which would carry around four parts (1, eps, im, eps*im). But, even if that works out, you wouldn't want to be able to have a Dual{Dual} or Complex{Complex}. On Dec 9, 2014 5:38 PM, "Theodore Papamarkou" notifications@github.com wrote:

That's true, dual numbers was meant to be a neat short impelementation of dual numbers. Imo we shouldn't be mixing dual and imaginary numbers in the scope of the present package. Mathematically speaking, both dual and imaginary numbers are special cases of the more general algebra of hypercomplex numbers. If there is an interest on merging treatment of complex and dual numbers, it may be a good idea, yet it would take some thinking.

— Reply to this email directly or view it on GitHub https://github.com/JuliaDiff/DualNumbers.jl/issues/13#issuecomment-66371773 .

[goretkin]

I got a ComplexDual type working. With minor modifications to DualNumbers (due to method ambiguity), it coexists with and requires DualNumbers. Here is the modification I made to DualNumbers https://github.com/goretkin/DualNumbers/tree/forcomplexdual I imagine that these changes are not consequential (perhaps except for performance) to any usage.

The implementation is here https://github.com/goretkin/ComplexDualNumbers.jl I discussed some design struggles in the README (probably related to issues of multiple inheritance).

The examples contain the problem I opened the issue with (taking the derivative of a real-valued function of multiple real arguments), along with another example that takes the derivative of a complex-valued function of a complex-argument, assuming the derivative exists.

[mlubin]

Pretty cool! Do you have any examples? You might want to raise the issue of Complex{T <: Real} in julia base if removing that restriction would make this work better.

[goretkin]

Some examples here: https://github.com/goretkin/ComplexDualNumbers.jl/tree/master/examples

I think it's going to get really tricky to not do Complex{T <: Real}. Ideally the type system could capture the structure of mathematics (or at least hypercomplex numbers), whatever that means. I would like to define Complex{T} over any T that makes sense, including T=Dual. But Complex{Dual} should look exactly like Dual{Complex}, and I shouldn't be able to do Complex{Complex} if we're assuming the same imaginary unit in both complex numbers.

Now Complex{Complex} could make sense, as one friend pointed out, if the imaginary unit for the ouer most Complex is i, the unit for the innermost Complex is k and then you call j=ik=ki (see https://en.wikipedia.org/wiki/Bicomplex_number), but I don't know how to achieve that sort of thing.

It is tempting to define a HyperComplex type that is parameterized by the multiplication table.

[mlubin]

Related: http://en.wikipedia.org/wiki/Dual_quaternion Just found out that this is actually implemented in https://github.com/forio/Quaternions.jl

[goretkin]

Awesome to learn that Dual Quaternions exist.

I think that for my use of ComplexDual numbers, I should have real and imag return Dual numbers corresponding to the "real" and "imaginary" parts of Complex{Dual}. I think, possibly controversially, that DualNumbers should also have real return the whole dual number because then

1. (x+conj(x))/2 == real(x) for x Dual
2. real is a no-op for real numbers. Any function that is defensive against complex numbers by taking the real of its input, you can still take the derivative

Of course, reasonable people expect real and imag to return real Real numbers, so this might also cause problems. For ComplexDual, though, this seems to make the most sense. someone could take real(x)^2+imag(x)^2 instead of abs2(x), and then you couldn't take the AD of that function.

I think there's a problem that we want Dual numbers to behave exactly the same in any code that is written to work with Real numbers (and similarly ComplexDual numbers for any code that is written to work with Complex numbers), but it doesn't feel good to have Dual <: Real (and you can't even have ComplexDual <: Complex, at least not in v0.3)