matrix exp(A) returns real-valued adjoint for real-valued A

[sdewaele]

The current adjoint for the matrix exp(A) where A is real-valued can return a complex-valued matrix. The complex components are only small numerical errors, close to machine precision. However, the fact that it is complex-valued leads to issues in the adjoint computation, see below for an example. This PR converts the return value to a real-valued array for real-valued arguments.

BTW - This issue has some similarities to #402, in that

1. Having indexing in the test code revealed the issue. For that reason, it may be an idea to add a gradcheck_getindex, to complement the current gradcheck? Or a test of compatibility of the return value of the adjoint with the input type?

2. An alternative fix would be to allow the adjoint of the indexed array to be complex-valued. This is similar to allowing unconstrainted matrices as adjoints in Adjoints for functions of specialized matrices #402. This would also require changes to the getindex adjoint.

# Note: A randomly generated 8 × 8 also fails (almost?) always. Using this confirmed example to be sure.
A = [ 0.0    1.0    0.0
      0.0    0.0    1.0
     -4.34 -18.31  -0.43]
x = [1.0]

f(A,x) = exp(A*x[1])
Y,Bf = Zygote.pullback(f,A,x)
Ȳ = ones(3,3)
Ā,x̄ = Bf(Ȳ) # fails

Error:

ERROR: InexactError: Float64(14.999150171071268 + 9.678155693021492e-16im)
Stacktrace:
 [1] Type at .\complex.jl:37 [inlined]
 [2] convert at .\number.jl:7 [inlined]
 [3] setindex!(::Array{Float64,1}, ::Complex{Float64}, ::Int64) at .\array.jl:767
 [4] #980 at C:\Users\user\.julia\packages\Zygote\ycnjm\src\lib\array.jl:32 [inlined]
 [5] #2619#back at C:\Users\user\.julia\packages\ZygoteRules\6nssF\src\adjoint.jl:49 [inlined]
 [6] literal_getindex at C:\Users\user\.julia\packages\Zygote\ycnjm\src\lib\lib.jl:77 [inlined]
 [7] (::typeof(∂(Zygote.literal_getindex)))(::Complex{Float64}) at C:\Users\user\.julia\packages\Zygote\ycnjm\src\compiler\interface2.jl:0
 [8] f at .\REPL[7]:1 [inlined]
 [9] (::typeof(∂(f)))(::Array{Float64,2}) at C:\Users\user\.julia\packages\Zygote\ycnjm\src\compiler\interface2.jl:0
 [10] (::getfield(Zygote, Symbol("##28#29")){typeof(∂(f))})(::Array{Float64,2}) at C:\Users\user\.julia\packages\Zygote\ycnjm\src\compiler\interface.jl:38
 [11] top-level scope at none:0

[mcabbott]

Edit: reading failure, I just realised this is a PR not an issue! If the gradient of exp is returning complex numbers only due to numerical rounding, that seems like a problem worth solving.

Earlier: There is also the question of whether real parameters should by default be allowed to return complex gradients, as is the current behaviour, see #342 for a discussion. But while they are, the fix here seems too narrow, as it does not allow for real A but complex dA on the backward pass.

And finally, given this behaviour, the gradient of getindex should promote as necessary. This is #376, which is fixed by #256. After this fix indexing will not be a diagnostic in the manner suggested.

[sdewaele]

Thanks for your comments! I will update the PR accordingly, and fix the test issue.

[sdewaele]

The check failed for Julia 1.3 with error

No output has been received in the last 10m0s,
this potentially indicates a stalled build or something wrong with the build itself.

I am not sure why this is happening; when I test on my machine with Julia 1.3 running the tests takes 5 minutes. Versioninfo:

julia>  versioninfo()
Julia Version 1.3.0
Commit 46ce4d7933 (2019-11-26 06:09 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-7500U CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)

[MikeInnes]

Travis probably just has a slow CPU. Bors should be fine:

bors r+

[bors]

Build succeeded

• ci/gitlab/gitlab.com