# Choosing a convention for complex numbers in DifferentiationInterface
**URL:** https://discourse.julialang.org/t/choosing-a-convention-for-complex-numbers-in-differentiationinterface/124433
**Category:** Machine Learning
**Tags:** autodiff, complex-numbers, di
**Created:** [January 4, 2025, 7:30pm UTC](https://discourse.julialang.org/t/choosing-a-convention-for-complex-numbers-in-differentiationinterface/124433 "2025-01-04T19:30:45Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [January 4, 2025, 7:30pm UTC](https://discourse.julialang.org/t/choosing-a-convention-for-complex-numbers-in-differentiationinterface/124433/1 "2025-01-04T19:30:45Z")
Hi everyone,
As DifferentiationInterface.jl gains more traction, I find myself forced to confront the issue of AD with complex numbers. Right now they are not officially supported, but they should be and I wonder what the community needs.
For instance, the derivative of a function \mathbb{C} \to \mathbb{C} can be represented as a matrix in \mathbb{R}^{2 \times 2}, or by a couple in \mathbb{C}^2, but not by a single number (except in the holomorphic case). I need to make such a decision for each operator, i.e. fill the following table and figure out for each cell:
- Should I error or return something?
- What is the correct convention for the result?
| operator | \mathbb{R} \to \mathbb{C} | \mathbb{C} \to \mathbb{R} | \mathbb{C} \to \mathbb{C} |
| --- | --- | --- | --- |
| `DI.derivative` (scalar to scalar) | ? | ? | ? |
| `DI.derivative` (scalar to array) | ? | ? | ? |
| `DI.gradient` (array to scalar) | ? | ? | ? |
| `DI.jacobian` (array to array) | ? | ? | ? |
| `DI.hessian` (array to scalar) | ? | ? | ? |
I also want to figure out how holomorphic behavior should be signaled to DI, if it can be exploited. Perhaps a function wrapper, or a keyword argument à la JAX?
Comments, wishes, pointers to relevant literature are all appreciated. Tagging a few opinionated people: @ChrisRackauckas @oxinabox @Mason @wsmoses @willtebbutt @avikpal
Related discussions:
>
>
> Find this while implementing https://github.com/SciML/BoundaryValueDiffEq.jl/pul…l/258
>
> This MWE is just proof of the idea and may not be meaningful, but I think it can still prove that the current DI lacks support for sparse Jacobian of complex numbers.
>
> \`\`\`julia
> using ForwardDiff, DifferentiationInterface, SparseMatrixColorings, ADTypes
> const DI = DifferentiationInterface
> \#backend = AutoForwardDiff() # this works
> backend = AutoSparse(
> AutoForwardDiff(),
> sparsity\_detector = ADTypes.KnownJacobianSparsityDetector(ones(3, 3)),
> coloring\_algorithm = ConstantColoringAlgorithm(ones(3, 3), ones(Int64, 3))
> )
> u0 = \[1.0; 0.0; 0.0\] .+1im
> jac\_cache = DI.prepare\_jacobian(nothing, similar(u0), backend, u0)
> \`\`\`
>
> stack trace
>
> \`\`\`julia
> ERROR: MethodError: no method matching DifferentiationInterfaceSparseMatrixColoringsExt.PushforwardSparseJacobianPrep(::DifferentiationInterface.BatchSizeSettings{…}, ::SparseMatrixColorings.ColumnColoringResult{…}, ::Matrix{…}, ::Vector{…}, ::Vector{…}, ::DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgPushforwardPrep{…})
>
> Closest candidates are:
> DifferentiationInterfaceSparseMatrixColoringsExt.PushforwardSparseJacobianPrep(::BS, ::C, ::M, ::S, ::R, ::E) where {BS\<:DifferentiationInterface.BatchSizeSettings, C\<:(AbstractColoringResult{:nonsymmetric, :column}), M\<:(AbstractMatrix{\<:Real}), S\<:(AbstractVector{\<:Tuple{Vararg{T, N}} where {N, T}}), R\<:(AbstractVector{\<:Tuple{Vararg{T, N}} where {N, T}}), E\<:DifferentiationInterface.PushforwardPrep}
> @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/bulUW/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:11
>
> Stacktrace:
> \[1\] \_prepare\_sparse\_jacobian\_aux\_aux(::DifferentiationInterface.BatchSizeSettings{…}, ::SparseMatrixColorings.ColumnColoringResult{…}, ::Vector{…}, ::Tuple{…}, ::AutoSparse{…}, ::Vector{…})
> @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/bulUW/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:107
> \[2\] \_prepare\_sparse\_jacobian\_aux(::DifferentiationInterface.PushforwardFast, ::Vector{…}, ::Tuple{…}, ::AutoSparse{…}, ::Vector{…})
> @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/bulUW/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:81
> \[3\] prepare\_jacobian(::Nothing, ::Vector{…}, ::AutoSparse{…}, ::Vector{…})
> @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/bulUW/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:49
> \[4\] top-level scope
> @ ~/Random/test2.jl:27
> Some type information was truncated. Use \`show(err)\` to see complete types.
> \`\`\`
>
> I believe the culprit is the annotation of the compressed matrix being too restricted in the SparseMatrixColorings extensions, I locally changed this into more generic ones and the errors are gone.
>
> https://github.com/JuliaDiff/DifferentiationInterface.jl/blob/cc8818a2bb0fb3dab2abf29ba213f89213a8613a/DifferentiationInterface/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl#L3-L33
> [@Taking Complex Autodiff Seriously in ChainRules](https://discourse.julialang.org/t/taking-complex-autodiff-seriously-in-chainrules/39317):
>
> So these conversations have been starting up recently again on the Slack #autodiff channel and in places like [https://github.com/JuliaDiff/ChainRulesCore.jl/issues/159](https://github.com/JuliaDiff/ChainRulesCore.jl/issues/159). I think there’s a lot of misconceptions out there about the derivatives of functions of complex numbers, caused in part by confusing notation and in part by many people’s education about complex numbers focusing too heavily on the holomorphic case. However, in an AD system where we want to deal with general code, we can not limi…
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [January 4, 2025, 8:02pm UTC](https://discourse.julialang.org/t/choosing-a-convention-for-complex-numbers-in-differentiationinterface/124433/2 "2025-01-04T20:02:16Z")
> [@gdalle](#):
>
> For instance, the derivative of a function \mathbb{C} \to \mathbb{C} can be represented as a matrix in \mathbb{R}^{2 \times 2}, or by a couple in \mathbb{C}^2, but not by a single number (except in the holomorphic case).
I’m partial to the CR calculus approach where you use a pair of complex numbers \partial f / \partial z and \partial f / \partial \bar{z}.
I find this easiest to use in the optimization context where you have functions \to \mathbb{R} in the end (in which case the two derivatives are conjugates), because it maps very easily onto a gradient that can be passed to optimization software. It also simplifies nicely in the holomorphic case (where the antilinear term \partial f / \partial \bar{z} = 0), and allows you to still think about things as complex numbers (rather than pairs of real and imaginary parts, which is often conceptually awkward).
But I’m not sure how this meshes with the internals of AD systems.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 4, 2025, 8:06pm UTC](https://discourse.julialang.org/t/choosing-a-convention-for-complex-numbers-in-differentiationinterface/124433/3 "2025-01-04T20:06:10Z")
I concour!
Also I am not sure if you should add in your table the JVP of a real function applied to complex vectors. I need this in bifurcationKit and use dispatch to handle this for now.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 7, 2025, 8:46am UTC](https://discourse.julialang.org/t/choosing-a-convention-for-complex-numbers-in-differentiationinterface/124433/4 "2025-01-07T08:46:17Z")
> [@gdalle](#):
>
> For instance, the derivative of a function \mathbb{C} \to \mathbb{C}C→C\mathbb{C} \to \mathbb{C} can be represented as a matrix in \mathbb{R}^{2 \times 2}R2×2\mathbb{R}^{2 \times 2}, or by a couple in \mathbb{C}^2C2\mathbb{C}^2, but not by a single number (except in the holomorphic case).
Actually using this would require a special Newton method and compiler analysis that do not currently exist. What we would need is:
1. Something in FunctionProperties.jl/DI for `isholomorphic`. This in theory can be done by constant propogation of the zeros through chain rules, where if the rules were writteen in the CR calculus approach then you would maybe be able to analyze whether you get a structural zero in the derivative of zbar.
2. Functionality in NonlinearSolve.jl, where if `isholomorphic` then create and use a Jacobian of a different size, and do the alternative Newton based on that property.
I’ve wanted such a thing for almost a decade now, but making sure we can properly declare something as holomorphic has been the hard part. The other thing that could be done is it could just be a user-set trait, that defaults to not holomorphic, and thus make the Newton default to the 2x2 form. The downside there is that a lot of “standard” ODE cases would then take a pretty big performance hit, and pretty much no one would know to actually set the trait.
So for now, we do the holomorphic Newton method simply because the holomorphic Newton requires 0 code changes, if you slap complex numbers into a standard Newton code it’s what you get, and if the user wants the other case they should just write their nonlinear solve using a real-valued u0 with the real and imaginary parts, and then inside of their `f` construct the complex numbers. This path of “try the simple thing, error often, but let the user solve this on their own if they want with a simple workaround” has been the okay equilibrium we’ve sat in.