WIP solving #174

[KnutAM]

Initial work on solving problem with ad for eigenvalues #174

[codecov-commenter]

Codecov Report

Merging #175 (bc6ec34) into master (c3b72b3) will decrease coverage by 0.07%.
The diff coverage is 66.66%.

@@            Coverage Diff             @@
##           master     #175      +/-   ##
==========================================
- Coverage   97.91%   97.84%   -0.08%     
==========================================
  Files          16       16              
  Lines        1296     1298       +2     
==========================================
+ Hits         1269     1270       +1     
- Misses         27       28       +1     

Impacted Files	Coverage Δ
src/Tensors.jl	81.81% <ø> (ø)
src/automatic_differentiation.jl	98.78% <ø> (ø)
src/eigen.jl	97.91% <66.66%> (-0.68%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c3b72b3...bc6ec34. Read the comment docs.

[KnutAM]

To start added internal get_max_iterations to make code run when called with dual number.
Does not resolve the correctness of the solution though. But

using Tensors

function M11f(A)
    l = eigvecs(A)
    return symmetric(l[:,1] ⊗ l[:,1])
end

A0 = SymmetricTensor{2,3}((1.2, 0., 0., 0.9, 0., 0.75))
B0 = rand(SymmetricTensor{2,3})

Δ = rand(SymmetricTensor{2,3})*1.e-6

println("A0: ", M11f(A0+Δ) ≈ M11f(A0) + gradient(M11f,A0)⊡Δ)
println("B0: ", M11f(B0+Δ) ≈ M11f(B0) + gradient(M11f,B0)⊡Δ)

produces

A0: false
B0: true

Seems like the wrong result is caused by a special case given by the A0 tensor in #174.
When debugging, no iterations are necessary for A0. My guess here is that we cannot expect AD to give the correct result for an iterative procedure if the input is such that no iterations occur? Perhaps derivatives must be specifically implemented?

Edit: The function could be discontinuous and AD would then not be properly defined.
@fredrikekre / @KristofferC is it then better to issue a warning/error. Or is it sufficient to add a warning about this in the docs for eigenvalue and/or automatic differentiation?

[KnutAM]

Based on Kristoffer's AD talk today, I think get_cos_sin() in eigen.jl could be to blame:

    function get_cos_sin(u::T,v::T) where {T}
        max_abs = max(abs(u), abs(v))
        if max_abs > 0
            u,v = (u,v) ./ max_abs
            len = sqrt(u^2 + v^2)
            cs, sn = (u,v) ./ len
            if cs > 0
                cs = -cs
                sn = -sn
            end
            T(cs), T(sn)
        else
            T(-1), T(0)
        end
    end

The problem is that len=0 (or at least real(len)) when we need to force it to go via the calculations with a Dual number, so we get division by zero. Using the following (I think equivalent code) doesn't work because differentiation produces NaN for cases when u=v=0

function get_cos_sin(u, v)
    α = atan(v, u) + π
    return cos(α), sin(α)
end

I don't have another solution on top of my head right now for this, so ideas are welcome!

[KnutAM]

function get_cos_sin(u, v)
    α = atan(v, u) + π
    return cos(α), sin(α)
end

Tested this with ForwardDiff's NaN-safe mode, but that did not work.

Note to possible future implementation: JuliaDiff/ForwardDiff.jl#111 (comment)

[KnutAM]

Solved by #182 and ForwardDiff.jl-PR575