Proper `epsilon` for complex?

[ChrisRackauckas]

It seems that the finite difference routines would "just work" if epsilon was changed to the "correct value". Unless there's some kind of norm argument, I would propose changing each part using the current epsilon. For example, if epsilon1 is from real(x) and epsilon2 is from imag(x), then epsilon = epsilon1 + im*epsilon2.

Would this be fine? This definition would make it safe for the case where a user is just "packing" floats into complex numbers, which can be common in physical applications. If this is acceptable I'll make a PR.

[johnmyleswhite]

To be clear, you want to do finite differencing on ℂ → ℂ functions?

[ChrisRackauckas]

Yes.

[ChrisRackauckas]

To clarify, ℂ → ℂ functions which sometimes are just R^2n -> R^2n packed in ℂ^N @AshtonSBradley

[antoine-levitt]

Not in favour of this. For instance, https://github.com/JuliaNLSolvers/Optim.jl/pull/435 would silently fail, because (unless I'm mistaken) the finite-differencing of a C^n to R function would be well-defined, but wrong (it would be in R^n). As I mentioned in https://github.com/JuliaDiff/ForwardDiff.jl/issues/157, I think the way to go is to define only those derivatives which make mathematical sense, ie functions C^n -> R.

What do you mean by C to C functions that are R^2n packed in C^n ? In this case the jacobian should be 2n x 2n and can not be done with a n x n complex matrix, right?

edit: to be clear, what I'm in favour of is that finite_difference assumes (and checks) a C^n to R function, computes the finite-differences as a R^2n to R function (i.e. 2n funcalls), and reassembles as a C^n vector. This is the only way that I can see that makes sense: the other possibility is to return a doubled-size jacobian (ie converting input and output to R^2n), but this conversion to reals should rather be done in the caller.

[ChrisRackauckas]

Not in favour of this. For instance, JuliaNLSolvers/Optim.jl#435 would silently fail, because (unless I'm mistaken) the finite-differencing of a C^n to R function would be well-defined, but wrong (it would be in R^n). As I mentioned in JuliaDiff/ForwardDiff.jl#157, I think the way to go is to define only those derivatives which make mathematical sense, ie functions C^n -> R.

It's the same result as doing this: https://github.com/JuliaDiff/DualNumbers.jl/pull/29#issuecomment-159878148 . ℂ → ℂ is the most standard case, so that should have first class treatment.

edit: to be clear, what I'm in favour of is that finite_difference assumes (and checks) a C^n to R function, computes the finite-differences as a R^2n to R function (i.e. 2n funcalls), and reassembles as a C^n vector. This is the only way that I can see that makes sense: the other possibility is to return a doubled-size jacobian (ie converting input and output to R^2n), but this conversion to reals should rather be done in the caller.

What about a keyword argument to enable this? That's double the function calls and very unnecessary in most cases, but is something that should be treated correctly.

[ChrisRackauckas]

Thinking about this more, probably default to safe, but keyword arg to opt into performance (holomorphic defaults to false, set it to true). I think that makes every application (including @antoine-levitt 's and @AshtonSBradley) safe when the derivative is well-defined in R^2n, but lets users take the fast route if complex differentiable.

[ChrisRackauckas]

Note that this, and many performance issues, have been solved by DiffEqDiffTools.jl, which is more of a dev library than a user-facing finite differencing package, but can be incorporated into packages similarly but give these features and a hefty speedup.