Repository: JuliaDiff/FiniteDifferences.jl
Pull request #76 (2020-04-21 11:08:06)

Fixes handling of complex numbers
Addresses #75 .

This is a first attempt, and it's not clear to be whether this is precisely the correct API for this -- we treat complex numbers as a 2-vector of reals. This lets us re-cycle the `to_vec` interface and avoid having to modify any of the code that actually does differencing.

@sethaxen does this address your use case? Maybe take a look at the extra tests to get a flavour of the kind of results you could expect to see.

willtebbutt: Note that test failures are to be expected at the minute. The main thrust of the PR is independent of the failures.
oxinabox: I suggest the following test case.
Which I have verified is different on this branch than on master.

On this branch it is completely real output.
but on master it is complex

```julia
using LinearAlgebra, FiniteDifferences
B = Symmetric(randn(7,7))
seed = one(B^-0.5)
_fdm = central_fdm(5, 1; adapt=5)
FiniteDifferences.j′vp(_fdm, ^, seed, B, -0.5)
```

Not sure what is right, probably cross check against ForwardDiff

sethaxen: > @sethaxen does this address your use case? Maybe take a look at the extra tests to get a flavour of the kind of results you could expect to see.

Yes, this PR seems to resolve the issue raised in #75. I agree that the issue @oxinabox found does not look right.
willtebbutt: The reason that Lyndon is seeing this behaviour is subtle and requires a proper explanation as it's an instance of a more general property of `FiniteDifferences`.

### A more general case

Consider the following:

```julia
julia> FiniteDifferences.grad(central_fdm(5, 1), sum, Diagonal(randn(5)))[1]
5×5 Diagonal{Float64,Array{Float64,1}}:
 1.0   ⋅    ⋅    ⋅    ⋅
  ⋅   1.0   ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅
  ⋅    ⋅    ⋅   1.0   ⋅
  ⋅    ⋅    ⋅    ⋅   1.0
```
which is distinct from what you might expect:
```julia
julia> FiniteDifferences.grad(central_fdm(5, 1), sum, collect(Diagonal(randn(5))))[1]
5×5 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
```

It stems from the idea of viewing types as constraints on values. In this scenario, it manifests itself as the gradient of the off-diagonal elements of a `Diagonal` matrix being undefined, because there is no notion of being able to wiggle around the values of the off-diagonal elements of a `Diagonal` matrix. It is in this way that `Diagonal(x)` is distinct from `collect(Diagonal(x))` from `FiniteDifferences`' perspective.

Phrased differently, it is not correct to say that `Diagonal(x)` can be thought of as a `Matrix` that _happens_ to have off-diagonal elements.

### The Complex case

The the situation that Lyndon presents, `Complex` is analogous to `Matrix`, and `Real` to `Diagonal`. That is, the result is essentially constrained to only provide the derivative w.r.t. the real component of `B` because `B` is `Real` (or, rather, its elements are), and it is clear why the complex part of the result is "dropped".

### Why the other behaviour is bad

Consider again the `Diagonal` case, and what would happen if one were to try to obtain a non-zero gradient w.r.t. the off-diagonal elements (based on the idea of a `Diagonal` matrix being a `Matrix` whose off-diagonal elements _happen_ to be zero). In this case, the matrix would necessarily be perturbed outside of the `Diagonal` matrices, and would presumably be represented by a `Matrix`. It is straightforward to imagine situations in which this breaks your code eg. if the function of which the gradient is to be computed is only defined for `Diagonal` matrices. Moreover, even if your code does succeed it won't be checking the correct code path -- it will dispatch to the method of your function specialised to `Matrix`s rather than `Diagonal`s, which is clearly undesirable.

Why stop there? For many operations you could view a `Diagonal` matrix as a special case of an `Array{T, 3}` in which all elements beyond the first 2D slice are zero, and only the diagonal of the 1st 2D slice is non-zero. For example, `sum` gives the same answer under all interpretations.

You can imagine the same situation for `Real` vs `Complex` numbers. To correctly compute the gradient w.r.t. the "complex bit" of a `Real` number, you would have to perturb it into a `Complex` number. This in turn could well yield issues for problems in code that isn't define to accept complex numbers, and won't test the execution paths that the user specifies.

Moreover, why stop at `Complex` numbers, when both `Real` and `Complex` numbers are in fact special cases of quaternions? You can see where I'm going here...

### Why can AD potentially do this differently?

You'll notice that AD has the option to do this differently, in that it can return eg. a complex gradient and _will_ execute the intended code. This is because the way that AD and finite differencing compute / estimate gradients are fundamentally different. AD doesn't have to modify the input, whereas input modification is fundamental to finite differencing.

(Whether or not returning a complex gradient for a real number in an AD system is a good idea is another matter)

### Conclusion

It seems to be the case that you have to truncate derivative information to get sensible results. Consequently, I'll fix up this PR and get it merged.

Technically the change is breaking, so a new minor version is in order.

#### nb.

the execution path argument above is why in the next minor release of FiniteDifferences we need to prevent users trying to differentiate w.r.t. integer inputs.
oxinabox: > Technically the change is breaking, so a new minor version is in order.

Why is this breaking?
Bug-fixes are not breaking
willtebbutt: Hmm maybe you're right? When I said technically, I guess what I really meant was that this feels pretty breaky. If we're happy to view it purely as a bug fix, then I'm happy to bump patch.
oxinabox: bug fixes are often breaking to if one assumes it is ok to depend on bugs.
but that path leads to madness

https://xkcd.com/1172/
![](https://imgs.xkcd.com/comics/workflow.png )
willtebbutt: 😂 excellent xkcd reference -- point taken.
simeonschaub: I just wanted to clarify, because I believe I was the one who initially implemented this for complex numbers: My original version assumes these are analytical `C -> C` functions, whereas it looks like with this change only `C -> R` functions are supported, but they don't have to be analytical, only differentiable in the reals. I am fine with this change, as this seems more useful, but it's not that the behavior before was a bug, it just has a different definition of derivatives.
willtebbutt: > whereas it looks like with this change only C -> R functions are supported, but they don't have to be analytical

Ish. It's the case for `grad`, but not for `jacobian`, `jvp`, or `jpvp`. With the last three it's still perfectly possible to work with functions that map to objects involving numbers in `C`, just the dimensionality is doubled.
simeonschaub: What I mean is that this is now treating complex vectors just as `R^2N` vectors, which is different from "complex differentiability" in the mathematical sense.
willtebbutt: @oxinabox based on @simeonschaub 's comment, we probably _should_ consider this breaking after all. The behaviour has changed from "assume function is analytic" to "don't assume analytic, and change conventions accordingly"
willtebbutt: > What I mean is that this is now treating complex vectors just as R^2N vectors, which is different from "complex differentiability" in the mathematical sense.

Yes. Although on the bug vs not-a-bug discussion, it's worth noting that the implementation before wasn't working as intended in the sense that gradients were being dropped (as per Seth's original [issue](https://github.com/JuliaDiff/FiniteDifferences.jl/issues/75))
simeonschaub: Yes, it just always assumed analytical functions.
willtebbutt: @oxinabox ready for review.
willtebbutt: Note that the CI appears to be lying here -- tests have passed.

Patch:
diff --git a/Project.toml b/Project.toml
index 94b6770..4b333df 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FiniteDifferences"
 uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
-version = "0.9.6"
+version = "0.9.7"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/grad.jl b/src/grad.jl
index 2635c3f..461f0c1 100644
--- a/src/grad.jl
+++ b/src/grad.jl
@@ -2,11 +2,11 @@
     jacobian(fdm, f, x...)
 
 Approximate the Jacobian of `f` at `x` using `fdm`. Results will be returned as a
-`Matrix{<:Number}` of `size(length(y_vec), length(x_vec))` where `x_vec` is the flattened
+`Matrix{<:Real}` of `size(length(y_vec), length(x_vec))` where `x_vec` is the flattened
 version of `x`, and `y_vec` the flattened version of `f(x...)`. Flattening performed by
 [`to_vec`](@ref).
 """
-function jacobian(fdm, f, x::Vector{<:Number}; len=nothing)
+function jacobian(fdm, f, x::Vector{<:Real}; len=nothing)
     len !== nothing && Base.depwarn(
         "`len` keyword argument to `jacobian` is no longer required " *
         "and will not be permitted in the future.",
@@ -38,11 +38,11 @@ end
 replace_arg(x, xs::Tuple, k::Int) = ntuple(p -> p == k ? x : xs[p], length(xs))
 
 """
-    _jvp(fdm, f, x::Vector{<:Number}, ẋ::AbstractVector{<:Number})
+    _jvp(fdm, f, x::Vector{<:Real}, ẋ::AbstractVector{<:Real})
 
 Convenience function to compute `jacobian(f, x) * ẋ`.
 """
-function _jvp(fdm, f, x::Vector{<:Number}, ẋ::Vector{<:Number})
+function _jvp(fdm, f, x::Vector{<:Real}, ẋ::Vector{<:Real})
     return fdm(ε -> f(x .+ ε .* ẋ), zero(eltype(x)))
 end
 
@@ -75,7 +75,7 @@ end
 
 j′vp(fdm, f, ȳ, xs...) = j′vp(fdm, xs->f(xs...), ȳ, xs)[1]
 
-function _j′vp(fdm, f, ȳ::Vector{<:Number}, x::Vector{<:Number})
+function _j′vp(fdm, f, ȳ::Vector{<:Real}, x::Vector{<:Real})
     return transpose(first(jacobian(fdm, f, x))) * ȳ
 end
 
diff --git a/src/to_vec.jl b/src/to_vec.jl
index 002512d..b42feb3 100644
--- a/src/to_vec.jl
+++ b/src/to_vec.jl
@@ -4,15 +4,22 @@
 Transform `x` into a `Vector`, and return the vector, and a closure which inverts the
 transformation.
 """
-function to_vec(x::Number)
-    function Number_from_vec(x_vec)
+function to_vec(x::Real)
+    function Real_from_vec(x_vec)
         return first(x_vec)
     end
-    return [x], Number_from_vec
+    return [x], Real_from_vec
 end
 
-# Base case -- if x is already a Vector{<:Number} there's no conversion necessary.
-to_vec(x::Vector{<:Number}) = (x, identity)
+function to_vec(z::Complex)
+    function Complex_from_vec(z_vec)
+        return Complex(z_vec[1], z_vec[2])
+    end
+    return [real(z), imag(z)], Complex_from_vec
+end
+
+# Base case -- if x is already a Vector{<:Real} there's no conversion necessary.
+to_vec(x::Vector{<:Real}) = (x, identity)
 
 function to_vec(x::AbstractVector)
     x_vecs_and_backs = map(to_vec, x)


Test patch:
diff --git a/test/grad.jl b/test/grad.jl
index e47667a..44097ce 100644
--- a/test/grad.jl
+++ b/test/grad.jl
@@ -2,7 +2,7 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec
 
 @testset "grad" begin
 
-    @testset "jvp(::$T)" for T in (Float64, ComplexF64)
+    @testset "jvp(::$T)" for T in (Float64,)
         rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
         x, y = randn(rng, T, N), randn(rng, T, M)
         ẋ, ẏ = randn(rng, T, N), randn(rng, T, M)
@@ -14,7 +14,7 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec
         @test ż_manual ≈ ż_multi
     end
 
-    @testset "grad(::$T)" for T in (Float64, ComplexF64)
+    @testset "grad(::$T)" for T in (Float64,)
         rng, fdm = MersenneTwister(123456), central_fdm(5, 1)
         x = randn(rng, T, 2)
         xc = copy(x)
@@ -43,9 +43,12 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec
         @test xc == x
     end
 
-    @testset "jacobian / _jvp / _j′vp (::$T)" for T in (Float64, ComplexF64)
+    @testset "jacobian / _jvp / _j′vp (::$T)" for T in (Float64,)
         rng, P, Q, fdm = MersenneTwister(123456), 3, 2, central_fdm(5, 1)
-        ȳ, A, x, ẋ = randn(rng, T, P), randn(rng, T, P, Q), randn(rng, T, Q), randn(rng, T, Q)
+        ȳ = randn(rng, T, P)
+        A = randn(rng, T, P, Q)
+        x = randn(rng, T, Q)
+        ẋ = randn(rng, T, Q)
         Ac = copy(A)
 
         check_jac_and_jvp_and_j′vp(fdm, x->A * x, ȳ, x, ẋ, A)
@@ -54,7 +57,8 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec
         @test Ac == A
 
         # Prevent regression against https://github.com/JuliaDiff/FiniteDifferences.jl/issues/67
-        @test first(jacobian(fdm, identity, x)) ≈ one(Matrix{T}(undef, length(x), length(x)))
+        J = first(jacobian(fdm, identity, x))
+        @test J ≈ one(Matrix{T}(undef, size(J)))
     end
 
     @testset "multi vars jacobian/grad" begin
@@ -113,7 +117,7 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec
         end
     end
 
-    @testset "j′vp(::$T)" for T in (Float64, ComplexF64)
+    @testset "j′vp(::$T)" for T in (Float64,)
         rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
         x, y = randn(rng, T, N), randn(rng, T, M)
         z̄ = randn(rng, T, N + M)
@@ -124,4 +128,53 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec
         @test x̄ȳ_manual ≈ vcat(x̄ȳ_auto...)
         @test x̄ȳ_manual ≈ vcat(x̄ȳ_multi...)
     end
+
+    # Tests for complex numbers, to prevent regressions against
+    # https://github.com/JuliaDiff/FiniteDifferences.jl/pull/76 and
+    # https://github.com/JuliaDiff/FiniteDifferences.jl/issues/75
+    @testset "Complex correctness - $T" for T in (ComplexF64,)
+        rng = MersenneTwister(123456)
+        x = randn(rng, T)
+        y = randn(rng, T)
+        fdm = FiniteDifferences.Central(5, 1)
+
+        # Addition.
+        dx, dy = FiniteDifferences.jacobian(fdm, +, x, y)
+        @test dx ≈ [1 0; 0 1]
+        @test dy ≈ [1 0; 0 1]
+
+        # Negation.
+        dx, = FiniteDifferences.jacobian(fdm, -, x)
+        @test dx ≈ [-1 0; 0 -1]
+
+        # Multiplication.
+        dx, dy = FiniteDifferences.jacobian(fdm, *, x, y)
+        @test dx ≈ [real(y) -imag(y); imag(y) real(y)]
+        @test dy ≈ [real(x) -imag(x); imag(x) real(x)]
+
+        # Magnitude
+        dx, = FiniteDifferences.grad(fdm, abs2, x)
+        @test real(dx) ≈ 2 * real(x)
+        @test imag(dx) ≈ 2 * imag(x)
+
+        # Concatenate tangents / cotangents.
+        ẋ = randn(rng, T)
+        ẏ = randn(rng, T)
+        ẋẏ = vcat(first(to_vec(ẋ)), first(to_vec(ẏ)))
+        z̄ = randn(rng, T)
+        z̄_vec = first(to_vec(z̄))
+
+        # Check jvp against jacobian result.
+        f(x, y) = 5 * sin(x) * cos(y)
+        Jx, Jy = FiniteDifferences.jacobian(fdm, f, x, y)
+        ż = jvp(fdm, f, (x, ẋ), (y, ẏ))
+        ż_vec = [real(ż), imag(ż)]
+        ż_manual = hcat(Jx, Jy) * ẋẏ
+        @test ż_vec ≈ ż_manual
+
+        # Check j′vp against jacobian result.
+        x̄, ȳ = j′vp(fdm, f, z̄, x, y)
+        @test [real(x̄), imag(x̄)] ≈ Jx'z̄_vec
+        @test [real(ȳ), imag(ȳ)] ≈ Jy'z̄_vec
+    end
 end