Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for Dual{Complex} #29

Merged
merged 1 commit into from
Dec 2, 2015
Merged

Add support for Dual{Complex} #29

merged 1 commit into from
Dec 2, 2015

Conversation

MikaelSlevinsky
Copy link
Contributor

Dear dualists (duellers?)

We would like to use DualNumbers in ApproxFun and dependent packages such as SingularIntegralEquations as a requirement. In ApproxFun, we create expansions of functions in numerous spectral bases (Chebyshev, Fourier, Legendre, Jacobi, etc...) and solve ordinary and partial differential equations in these bases as well.

One particular use of dual numbers will be in automatic differentiation to compute normal derivatives of incident waves to solve elliptic PDEs like the Helmholtz equation with Neumann boundary conditions in SingularIntegralEquations.

To make the transition possible, DualNumbers must be extended to work on Complex data types as well. This required a minor rewrite, and for the most part an improvement.

I would be happy to see this merged, though I can be patient given dependencies on DualNumbers that would break.

If merged, this would close #13, #23 (as a patch), and #24.


function squareroot(x)
it = x
while abs(it*it - x) > 1e-13
it = it - (it*it-x)/(2it)
it = (it+x/it)/2
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any reason for this change?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mlubin
Copy link
Contributor

mlubin commented Nov 20, 2015

This looks potentially very useful, thanks! We've been waiting for an application where this was needed, do you have an example or small write up that explains why dual numbers make sense with complex element types? What about Complex{Dual}?

Please deprecate the methods which are being removed, also use @Base.deprecate_binding for the renaming of du.

@mlubin
Copy link
Contributor

mlubin commented Nov 20, 2015

CC @goretkin

@MikaelSlevinsky
Copy link
Contributor Author

I guess this breaks with 0.3 as well.

Some toy examples are here, and with this PR, they carry over naturally to automatic differentiation of functions of a complex variable.

I haven't thought about Complex{Dual}, but since Julia's complex is actually Complex{T<:Real}, this would only be possible with Dual <: Real, which would be odd given Dual{Complex{Float64}}.

@mlubin
Copy link
Contributor

mlubin commented Nov 20, 2015

I'm okay with breaking 0.3, just be sure to update REQUIRE and the travis settings.

@goretkin
Copy link

I haven't followed along too carefully, but I imagine that ApproxFun needs to take some derivatives through, for example, a DFT.

You could define a DFT that doesn't use any complex arithmetic, but the way the DFT is usually expressed relies on complex arithmetic.

This is from a while back, but here's a sort-of clear example https://github.com/goretkin/ComplexDualNumbers.jl/blob/master/examples/dft_autodiff.jl

I defined two functions spec and spec_nocomplex which assigns some scalar to the spectrum of its inputs. They're the same, but internally one uses complex arithmetic. The point is that it's a R^n->R function, and you should be able to take its derivative.

@goretkin
Copy link

For me, the biggest question is still what to do about Complex{Dual} and Dual{Complex}. Current definitions aside, it does make sense to interpret a "Complex Dual" number both as a Complex number with Dual parts, and as a Dual number with Complex parts. The same goes for Dual Quaternions.

I don't think you necessarily have to answer that question in order to get this working. It wouldn't be so terrible to have Dual{Complex} be legal, but not Complex{Dual}, but I would love to find a way.
Perhaps defining Complex{T<:Number}<:Number and Dual{T<:Number}<:Number are the way to go. Contrary to what I might have said before, both Complex{Complex{Real}} and Dual{Dual{Real}} are meaningful.


immutable Dual{T<:ReComp} <: Number
value::T
dual::T
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First of all, +1 to this, thanks for writing it up. I'm really interested in seeing this in action.

Can we use epsilon instead of dual for the ɛ component? I feel like it would be clearer, and also less confusing for users transitioning to the new version of DualNumbers.jl (since before, dual was the name of the vectorized constructor).

@MikaelSlevinsky
Copy link
Contributor Author

Thanks for taking a look!

@mlubin, I deprecated the binding for du and modified the require and travis build file.

@jrevels, I tossed the typealias constructors, reverted dual to epsilon and re-added dual as a vectorized constructor.

Now, the type Dual contains a value and an epsilon, the display is an ɛ and the method epsilon() gets the epsilon component of a Dual. Very consistent. (The main issue was the use of real(::Dual) really.)

@goretkin, @dlfivefifty showed me that the DFT/DCT of Vector{Dual} can just be the dual(DCT(value),DCT(epsilon)), so Dual{Complex} satisfies our needs.

language: julia
os:
- linux
- osx
julia:
- 0.3
- 0.4
- release
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be 0.4

@MikaelSlevinsky
Copy link
Contributor Author

@mlubin, fixed!


typealias Dual64 Dual{Float32}
typealias Dual32 Dual{Float16}
typealias DualComplex256 Dual{Complex128}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you find these useful or is this just for consistency?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know if I'll use them more often than the generic constructor, but I guess their purpose is to mirror Float64 and Complex128 etc...

@mlubin
Copy link
Contributor

mlubin commented Nov 23, 2015

The name value seems too generic to export, it will lead to conflicts with other packages. Is there some other existing method in Base which we can extend?

@MikaelSlevinsky
Copy link
Contributor Author

The names I can think of are:

num
values # assumes there are many
var
base
read # this is a stretch

@dlfivefifty
Copy link
Collaborator

How about “diag”, since dual(a,b) is equivalent to

[a b; 0 a]

On 25 Nov 2015, at 10:41 PM, Richard Mikael Slevinsky [email protected] wrote:

The names I can think of are:

num
values # assumes there are many
var
base
read # this is a stretch

Reply to this email directly or view it on GitHub #29 (comment).

@mlubin
Copy link
Contributor

mlubin commented Nov 26, 2015

What about defining:

real{T<:Real}(x::Dual{T}) = x.value

In addition to avoiding needless deprecations in existing code, real just seems like the mathematically reasonable description of the first component of a dual number when the components themselves are real, by analogy to the complex numbers. We don't need to generalize this away... most users will be interacting with Dual{Float64}, so I think it's a useful concept to keep in place. For the special case of complex-valued components, maybe an unexported DualNumbers.value or diag or one of the other suggestions?

@mlubin
Copy link
Contributor

mlubin commented Nov 26, 2015

Another option is to define

getindex(x::Dual, idx) = (idx == 1) ? x.value : x.epsilon # plus some error checking

and use x[1] for the first component and x[2] for the second. There's a runtime penalty though, so this would be the cute way but not necessarily the fastest way to get the components.

@dlfivefifty
Copy link
Collaborator

How about just using “x.value” and “x.epsilon”?

On 26 Nov 2015, at 12:15 PM, Miles Lubin [email protected] wrote:

Another option is to define

getindex(x::Dual, idx) = (idx == 1) ? x.value : x.epsilon # plus some error checking
and use x[1] for the first component and x[2] for the second. There's a runtime penalty though, so this would be the cute way but not necessarily the fastest way to get the components.


Reply to this email directly or view it on GitHub #29 (comment).

@mlubin
Copy link
Contributor

mlubin commented Nov 26, 2015

I don't see a reason to remove the current accessors for the real case, but if the dot syntax helps for the complex case then I'm not opposed.

@jrevels
Copy link
Member

jrevels commented Nov 26, 2015

+1 for unexported DualNumbers.value, and just letting people use x.value/x.epsilon otherwise.

I'm also in favor of the real definition @mlubin proposed, with the caveat that there should probably be an explicit error message for real(::Dual{T<:Complex}).

@dlfivefifty
Copy link
Collaborator

real(x::Dual{T<:Complex}) should return dual(real(x.value),real(x.epsilon))

On 26 Nov 2015, at 4:58 PM, Jarrett Revels [email protected] wrote:

+1 for unexported DualNumbers.value, and just letting people use x.value/x.epsilon otherwise.

I'm also in favor of the real definition @mlubin https://github.com/mlubin proposed, with the caveat that there should probably be an explicit error message for real(::Dual{T<:Complex}).


Reply to this email directly or view it on GitHub #29 (comment).

@jrevels
Copy link
Member

jrevels commented Nov 26, 2015

real(x::Dual{T<:Complex}) should return dual(real(x.value),real(x.epsilon))

Isn't there an expectation that real(::Number) return a T<:Real? This expectation would be violated if real(::Dual) were to ever return a Dual.

@dlfivefifty
Copy link
Collaborator

I think “real” here is short for “realpart”, like “imag” is short for “imagpart”. So has nothing to do with the type Real.

I’d expect Real(::Number) to return a Real.

On 26 Nov 2015, at 5:07 PM, Jarrett Revels [email protected] wrote:

real(x::Dual{T<:Complex}) should return dual(real(x.value),real(x.epsilon))

Isn't there an expectation that real(::Number) to return a T<:Real? This expectation would be violated if real(::Dual) were to ever return a Dual.


Reply to this email directly or view it on GitHub #29 (comment).

@MikaelSlevinsky
Copy link
Contributor Author

👍 to unexported value. (For frequent use, it's a one-liner to import DualNumbers: value anyway.)

If we choose to be mathematically rigorous, shouldn't real{T<:Real}(z::Dual{T}) = z and real{T<:Complex}(z::Dual{T}) be undefined (or throw an error) since Re(z) is not complex differentiable? This makes deprecation complicated, however.

The same holds for conj and imag, but imag{T<:Real}(z::Dual{T}) should return zero(z).

One option is to add realdual and imagdual like conjdual for the other notions that are being discussed.

@dlfivefifty
Copy link
Collaborator

Doesn't dual tell the angle of differentiation?  So it’s not about complex differentiability, but differentiability at an angle?

On 26 Nov 2015, at 9:34 PM, Richard Mikael Slevinsky [email protected] wrote:

to unexported value. (For frequent use, it's a one-liner to import DualNumbers: value anyway.)

If we choose to be mathematically rigorous, shouldn't real{T<:Real}(z::Dual{T}) = z and real{T<:Complex}(z::Dual{T}) be undefined (or throw an error) since Re(z) is not complex differentiable? This makes deprecation complicated, however.

The same holds for conj and imag, but imag{T<:Real}(z::Dual{T}) should return zero(z).

One option is to add realdual and imagdual like conjdual for the other notions that are being discussed.


Reply to this email directly or view it on GitHub #29 (comment).

@MikaelSlevinsky
Copy link
Contributor Author

The dual part returns the epsilon component of z multiplied by f'(z.value), so:

julia> z = Dual(1.0,1.0+im)
1.0 + 0.0im  +  1.0 + 1.0imɛ

julia> exp(z)
2.718281828459045 + 0.0im  +  2.718281828459045 + 2.718281828459045imɛ

Doesn't this mean that the function must be complex differentiable?

@dlfivefifty
Copy link
Collaborator

No I think it's just

lim (f(z.value + h*z.epsilon) - f(z.value))/h

When f is complex differentiable, this reduces to f'(z.value)z.epsilon, but that's just a special case.

Sent from my iPhone

On 26 Nov 2015, at 9:55 PM, Richard Mikael Slevinsky [email protected] wrote:

The dual part returns the epsilon component of z multiplied by f'(z.value), so:

julia> z = Dual(1.0,1.0+im)
1.0 + 0.0im + 1.0 + 1.0imɛ

julia> exp(z)
2.718281828459045 + 0.0im + 2.718281828459045 + 2.718281828459045imɛ
Doesn't this mean that the function must be complex differentiable?


Reply to this email directly or view it on GitHub.

@dlfivefifty
Copy link
Collaborator

real returning a type that's not T<:Real feels like a jarring deviation from the expectations of the existing real definitions in Base.

Actually, it is not true that real always returns a real in Base:

julia> real([1+im 2; 3 4])
2x2 Array{Int64,2}:
 1  2
 3  4

@jrevels
Copy link
Member

jrevels commented Dec 1, 2015

Actually, it is not true that real always returns a real in Base:

I should've been more clear, I meant that real(::Number) always returns a Real (though, as I mentioned, real is vectorized, so obviously whatever element-wise behavior we define will carry over to the vectorized definition).

It seems that these definitions:

real{T<:Real}(x::Dual{T}) = x
real{T<:Complex}(x::Dual{T}) = Dual(real(value(x)), real(epsilon(x)))

...now have some pretty substantial arguments backing them, and I'm close to being won over. I especially feel like it would be very difficult to claim our Dual implementation supports complex autodiff without defining real properly, since real is sure to pop up everywhere in user code involving Complex. However, are we absolutely sure that the above definition of real delivers the correct results in the autodiff case? I basically believe it, but I'd need to test some edge cases to really convince myself.

I also do still share the same concerns as @mlubin - the standard nomenclature when discussing dual numbers is to refer to the non-perturbative component as the "real" component, and the function real as defined in Base means "give me the real component of this number."

I'm not sure we can resolve this without introducing inconsistency somewhere, but in the end, if inconsistency is introduced, I'd rather have that inconsistency be nominative than behavioral (i.e. I prefer distorting Base's assumption of real to breaking autodiff).

@MikaelSlevinsky
Copy link
Contributor Author

@jrevels, the dual part of the real function can be thought of as the dot product of the gradient of f(x,y) = x with reim(epsilon(z))...:

julia> using DualNumbers

julia> z = Dual(1.0+im,1.0)
1.0 + 1.0im + (1.0 + 0.0im)ɛ

julia> real(z)
1.0 + 1.0ɛ

julia> z = Dual(1.0+im,2.0+im)
1.0 + 1.0im + (2.0 + 1.0im)ɛ

julia> real(z)
1.0 + 2.0ɛ

julia> z = Dual(1.0+im,im)
1.0 + 1.0im + (0.0 + 1.0im)ɛ

julia> real(z)
1.0 + 0.0ɛ

As far as I can tell, this is correct.

Thoughts on either of these three ways forward?

  1. restart with
immutable Dual{T<:ReComp}
    primal::T
    dual::T
end

and rename the accessors value (EDIT: real) => primal and epsilon => dual. This would break dual, real and epsilon.

  1. Define realdual and imagdual in the same sense as absdual, abs2dual, conjdual. As far as I can see, the original purpose for these functions is computationally equivalent to the interchange im <=> ɛ. I would be willing to deprecate the current DualNumbers real for v0.1.6 and remove the deprecation in v0.1.7 EDIT: in favour of the limit definition real{T<:Real}(x::Dual{T}) = x.
  2. Keep this pull request on hold while I issue a pull request to Julia base for value. We have a strong case, but I fear that Base doesn't exist to resolve naming conflicts in external packages.

@mlubin
Copy link
Contributor

mlubin commented Dec 1, 2015

Who would've thought this would be so hard? :)

How about:

realpart{T<:Real}(x::Dual{T}) = x.value
value(x::Dual) = x.value
epsilon(x::Dual) = x.epsilon
export realpart, epsilon
# After a cycle of deprecations, define:
real{T<:Real}(x::Dual{T}) = x
real{T<:Complex}(x::Dual{T}) = Dual(real(value(x)), real(epsilon(x)))

I really don't want to give up the analogy of the first component being the real component in the case of Dual{Float64}. realpart seems less confusing than realdual.

@MikaelSlevinsky
Copy link
Contributor Author

Sounds good to me! I can make these changes if you're ok with me using your code.

@mlubin
Copy link
Contributor

mlubin commented Dec 1, 2015

Sure, go ahead

@MikaelSlevinsky
Copy link
Contributor Author

yet dual numbers can be generalized over any commutative ring R (real or complex numbers for instance) as the quotient ring R[X]/X^2 (so long as ɛ^2 = 0). https://en.wikipedia.org/wiki/Dual_number#Generalization

This page makes me want to have dualpart as well as realpart http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/other/dualNumbers/index.htm

@goretkin
Copy link

goretkin commented Dec 1, 2015

@MikaelSlevinsky I didn't understand whether that was in support of removing the derivative definition or not.

@MikaelSlevinsky
Copy link
Contributor Author

Sorry, I misread completely. That part is in the section on overriding Base functions, whereas the paragraphs above describe dual numbers more generally.

Usually, evaluating a function at a dual number returns a dual number, so what's wrong with the definition (in the context of that part of the readme)? Suggest any changes?

@jrevels
Copy link
Member

jrevels commented Dec 1, 2015

I don't want to hold up this PR with too much bikeshedding, but if we're leaning towards defining realpart/dualpart instead of value/epsilon, maybe it'd be better to say valuepart/dualpart (to make it less confusing when the x.value is Complex)? Otherwise, the proposed realpart/dualpartdefinitions seem like the right way to resolve this.

@mlubin
Copy link
Contributor

mlubin commented Dec 1, 2015

I'd much prefer realpart over valuepart.

@jrevels
Copy link
Member

jrevels commented Dec 1, 2015

Fair enough, you have more experience demoing this stuff then I do.

Either way, +1 to both component accessors having names of the form ___part. I like the consistency there.


value(z::Dual) = z.value
epsilon(z::Dual) = z.epsilon
value(z::Number) = z
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When would this definition be used?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can't tell exactly but if you mean line 32: value(z::Number) = z: if you know z's a Number but not necessarily a Dual it's good form to have a no-op.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd be in favor of removing the no-op until someone requests it for a good reason, but it's just a minor point.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with @mlubin. I'm not aware of other definitions of value, but I am thinking of it as an accessor of the non-dual part of a dual number, and not as a "mathematical function" so to speak, so I don't think we need to (yet) worry about finding a meaning for other types of numbers.

@MikaelSlevinsky
Copy link
Contributor Author

Ok I removed value(z::Number) and epsilon(z::Number). In the future, it may be useful to think of them as more than accessors - Val(z) and Eps(z) could be thought of as projection operators to the purely real/concrete and purely dual subspaces the same way as Re(z) and Im(z) act on complex z (EDIT: or real z).

I also added the constant imɛ and special show methods, so now

julia> using DualNumbers

julia> z = 1.0+2.0im-3.0ɛ+4.0imɛ
1.0 + 2.0im - 3.0ɛ + 4.0imɛ

julia> ɛ
ɛ

julia> im*ɛ
imɛ

@mlubin
Copy link
Contributor

mlubin commented Dec 2, 2015

Looks basically ready to merge. Could you squash the commits?

- remove dual() constructor (float() is deprecated in 0.4 anyways)

- replace real & epsilon with value and dual

- change du to ɛ addressing #24

- add Base.sinpi and Base.cospi as a patch addressing #23

- add complex function tests

add 2-norm(Vector{Dual})

add cis (cos+im*sin)

update readme to reflect more general support

fixes

dual => epsilon

re-add vectorized constructor dual

deprecate_binding du

release => 0.4 and deprecate real(::Dual)

dual => epsilon in readme

un-export value

incorporate some changes of #28:

__precompile()__

function => constructor in readme

import value in test

add definitions for non-complex differentiable functions

and tests

all generic norms work with isinf() and float()

fix print

so that you can copy paste and it will make a new Dual

add limit definition in readme

fix typo:

float is a no-op on Union{Dual{T},Dual{Complex{T}}}, not
Union{Dual{T},Complex{T}}

add no-op real{T<:Real}

add angle, realpart, and dualpart

deprecate real{T<:Real}(z::Dual{T}).

also, because sign := z/|z| in base, it works automatically and so #22
is no longer necessary.

remove value(z::Number) and epsilon(z::Number)

- add the constant imɛ
- add special dual_show methods for Dual{Bool}’s and
Dual{Complex{Bool}}’s

- replace real with DualNumbers.value

- change du to ɛ addressing #24

- add imɛ

- add Base.sinpi and Base.cospi as a patch addressing #23

- add cis (cos+im*sin)

- deprecate_binding du

- remove 0.3 from Travis build (and therefore support)

- deprecate real(::Dual)

- (real,epsilon) => (realpart,dualpart) in readme

- incorporate some changes of #28:

-- __precompile()__

- function => constructor in readme

- add isinf() and float() so that all generic norms work

- add show_dual{T<:Complex} so that you can copy paste and it will make a new Dual

- add show_dual{T<:Bool} and show_dual{T<:Complex{Bool}}

- add limit definition in readme

- add angle, realpart, and dualpart

- deprecate real{T<:Real}(z::Dual{T})

- sign := z/|z| in base, it works automatically and so #22
is no longer necessary

- add tests
@MikaelSlevinsky
Copy link
Contributor Author

Done! :)

mlubin added a commit that referenced this pull request Dec 2, 2015
Add support for Dual{Complex}
@mlubin mlubin merged commit 7d4eedf into JuliaDiff:master Dec 2, 2015
@mlubin
Copy link
Contributor

mlubin commented Dec 2, 2015

Thanks!
I'll let this sit on master for a week or so before tagging. That'll give me some time to check for unintended breakage. @goretkin, we're open for documentation tweaks.

@jrevels
Copy link
Member

jrevels commented Dec 2, 2015

🎉

@MikaelSlevinsky
Copy link
Contributor Author

Thanks again!

@MikaelSlevinsky
Copy link
Contributor Author

Is it time for a Pkg.tag? I'd like to push some code to SingularIntegralEquations but I want it to "require" the right version number of DualNumbers so users are unencumbered.

@mlubin
Copy link
Contributor

mlubin commented Dec 30, 2015

I'll tag tomorrow. I've been running on master for a while with no apparent issues.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Dual numbers and Complex numbers
5 participants