Docs: Complex numbers and Autodiff Jacobians

[rakhusur]

Hi, I couldn't figure out from the docs how to use derivative-based methods for complex problems. Some things are mentioned in the news blog, others I still don't know. It would be nice if the following things were mentioned:

• By default, derivative-based methods fail because they try to use auto differentiation. Also, a descriptive error message would be very helpful here.
• With autodiff = false or a custom Jacobian, it works. The question here is, what is the Jacobian supposed to be for a non-holomorphic function? I would have expected that I need to pass two functions for derivatices w.r.t real and imaginary part (or the two Wirtinger derivatives). Is it at all possible to use non-holomorphic functions? If not, will I get nonsense if I use them with autodiff = false? Cheers.

[ChrisRackauckas]

Yeah, this is more of an upstream issue with the differentiation libraries. No autodiff library currently does this, and we need a good check on the numerical diff libraries to make sure we have this right (I don't think they do: I think they currently do the directional derivative along the real, which assumes analytic).

That Jacobian is just used in the Newton iterations, so we'd have to work out what is correct for non-holomorphic. For implicit methods it's just a line search so it's okay if it's only approximate, but this will be an issue with Rosenbrock methods if it's not exactly correct. It should correspond to the same thing as if you redefined the system as the vector of real and imaginary parts, making complex numbers just a nice abstraction for the user. Given this, it might just be best to add hook ins to implement them this way, though then it would ignore any optimizations that can be had on analytic functions.

[daviehh]

It should correspond to the same thing as if you redefined the system as the vector of real and imaginary parts

I used to create a new array twice the size of the complex-valued one, storing the real/imaginary parts separately. Assuming the initial condition/solution is stored as a vector for now, would this be a good way when working with e.g. Rosenbrock methods?

Define

to_reals(x) = [real(vec(x)); imag(vec(x))]

function to_complexes(xri)
    rleng = Int(length(xri)//2)
    [xri[i] + xri[i+rleng]*im for i in 1:rleng]
end

then, convert the complex-valued initial condition with to_reals(), and use to_complexes() inside the ode function on u like uc = to_complexes(u), carry out the complex-valued equation on uc to get duc and return du = to_reals(duc)? Thanks!

[ChrisRackauckas]

Yes. That adds a few allocations to the solving though, so it might slow it down.

[rakhusur]

This is also what I do for now. I'm accessing the content of the real vector through these functions to avoid converting it:

getreal(v::Vector{<: Real}, i::Integer) = v[i]
getimag(v::Vector{<: Real}, i::Integer) = v[i + end÷2]
getcomplex(v::Vector{<: Real}, i::Integer) = complex(getreal(v, i), getimag(v, i))
setreal!(v::Vector{<: Real}, val::Real, i::Integer) = v[i] = val
setimag!(v::Vector{<: Real}, val::Real, i::Integer) = v[i + end÷2] = val
function setcomplex!(v::Vector{<: Real}, val::Number, i::Integer)
    setreal!(v, real(val), i)
    setimag!(v, imag(val), i)
    val
end

But perhaps it's still slow because real and imaginary part are not adjacent in memory (didn't really test this).