Symmetric/Hermitian matrix function rules

[sethaxen]

This PR adapts Zygote's adjoint rules for power series functions (analytic functions) matrix functions of real Symmetric and Hermitian inputs, originally introduced in FluxML/Zygote.jl#355 (derivation here)
It additionally adds frules for the same functions. Hermitian below also includes Symmetric{<:Real}.

A few notes (Hermitian below also includes Symmetric{<:Real}):

• It also adds rules for eigen and eigvals, since they share a number of utility functions with the power series rules. The eigen rule should probably be a fallback for [WIP] rrule for eigen decomposition  #179. Moved to Add rules for Symmetric/Hermitian eigen and svd #323
• The rules assume that the (co)tangents are of the same type as their corresponding primal. It also assumes they are Hermitian, etc, not Composite{<:Hermitian} (relates Representing (co)tangents of structured matrices #191).
• For the moment, it only adapts Zygote's rules for the type-stable functions, i.e. those whose scalar codomain is always real when their domain is real. (relates Matrix Powers are type-unstable: sometime Real sometimes Complex JuliaLang/julia#35667).
• The above combined with Hermitian matrices having real eigenvalues means that the rrules are not dependent on the conjugation convention of the cotangent.

To-do:

• Write tests (this will include resolving Power Series adjoints are subtly wrong in a way gradcheck can't catch FluxML/Zygote.jl#608, since it relates these same rules). This gist tests most of the cases for Zygote's versions.

Because only a few utility functions differ between the Hermitian{<:Complex} case and the Symmetric{<:Real} and Hermitian{<:Real} cases, I don't think not having FD v0.10.0 support to test the rules with complex inputs should hold up this PR.

[sethaxen]

These rules will probably not be usable with Zygote without some modification even after FluxML/Zygote.jl#366 is merged, due to Zygote not enforcing type constraints. e.g.

julia> using Zygote, ChainRules

julia> A = Symmetric(randn(3, 3));

julia> Zygote.pullback(sum ∘ exp, A)[2](1.0)[1]
3×3 Array{Float64,2}:
 16.1831    4.00247     2.91251
  4.00247  -0.0429689  -0.331864
  2.91251  -0.331864   -0.546545

julia> rrule(exp, A)[2](ones(size(A)))[2] |> unthunk # this is identical, and hey, we have type constraints!
3×3 Symmetric{Float64,Array{Float64,2}}:
 16.1831    4.00247     2.91251
  4.00247  -0.0429689  -0.331864
  2.91251  -0.331864   -0.546545

julia> Zygote.pullback(sum ∘ exp, A)[2](im)[1]
3×3 Array{Complex{Float64},2}:
 0.0+16.1831im  0.0+4.00247im    0.0+2.91251im
 0.0+4.00247im  0.0-0.0429689im  0.0-0.331864im
 0.0+2.91251im  0.0-0.331864im   0.0-0.546545im

julia> rrule(exp, A)[2](im*ones(size(A)))[2] |> unthunk # that's...not even close
3×3 Symmetric{Complex{Float64},Array{Complex{Float64},2}}:
 0.0-6.6388im   0.0+5.58432im   0.0+5.86321im
 0.0+5.58432im  0.0-0.608647im  0.0-0.610593im
 0.0+5.86321im  0.0-0.610593im  0.0-1.26019im

I'm not certain if there's an easy fix to make the two compatible.

[nickrobinson251]

Cool! And thanks!

I am not well placed to review the maths here. On a quick look thing seem sensible.
A few minor comments :) Hopefully we can find someone who can who's able review more thoroughly :)

[antoine-levitt]

oh that's more clever than what I do and gets you eps^2/3 (I only get eps^1/2), nice! Of course this is all assuming that all quantities are order 1 (or else you need to be careful about relative errors rather than absolute ones), but it's fine.

[sethaxen]

That's right, yeah. TBH, error analysis is not my thing. I plan to come back and try to improve this in the future. For the moment, this seems to be fine. I haven't been able to construct a random almost-degenerate matrix for which the pushforward/pullback disagrees with finite differences so far.

[antoine-levitt]

well it'll agree to eps^2/3 if what you wrote is correct. I can check the math if you want

[sethaxen]

That'd be great if you have the time!

[antoine-levitt]

Yup, checks out. The approximation (f(l1) - f(l2))/(l1-l2) ~= (f'(l1) + f'(l2))/2 is accurate to order dl^2. The roundoff error when computing (f(l1) - f(l2))/(l1-l2) is order eps/dl. So it's advantageous to switch to the approximation when eps/dl ~= dl^2 => dl = cbrt(eps). The worst case accuracy is indeed eps^2/3.

[sethaxen]

Awesome! Thanks for checking that!

[antoine-levitt]

Regarding reference, I'd just cite Higham, pretty sure he discusses DK theorem somewhere.

[sethaxen]

I made a few changes to ensure that the rules do the right thing when the wrapped array is immutable. This isn't tested though.

[sethaxen]

Blocked by JuliaDiff/ChainRulesCore.jl#279 and need to stabilize output type for jvp.

[oxinabox]

I don't love that we are using an explictly nonexported function.
Do we even need to be testing @inferred here, or can we rely on frule_test to do that for us later

[sethaxen]

Unfortunately we cannot use frule_test or rrule_test here because we need to work around the outputs of some functions being type-unstable. (and rrule's won't agree with FiniteDifferences for Symmetric/Hermitian, only with them composed with the constructor). I copied over a simplified is_inferrable.