Error computing j′vp of a scalar -> matrix function

[sethaxen]

I can't compute the j′vp for the function p -> x^p on FiniteDifferences v0.10.0, where x is a matrix and p is a scalar:

julia> using FiniteDifferences

julia> fdm = central_fdm(5, 1);

julia> x = [1.0 2.0; 3.0 4.0]; # real input

julia> x^3.0 # real output 
2×2 Array{Float64,2}:
 37.0   54.0
 81.0  118.0

julia> seed = [1.0 1.0; 1.0 1.0]; # real seed

julia> j′vp(fdm, p -> x^p, seed, 3.0)  # but an error is raised that one of the vectors has length 8?
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(8),), b has dims (Base.OneTo(4),), mismatch at 1")
Stacktrace:
 [1] promote_shape at ./indices.jl:178 [inlined]
 [2] promote_shape at ./indices.jl:169 [inlined]
 [3] +(::Array{Float64,1}, ::Array{Float64,1}) at ./arraymath.jl:45
 [4] add_sum(::Array{Float64,1}, ::Array{Float64,1}) at ./reduce.jl:21
 [5] _mapreduce(::FiniteDifferences.var"#22#24"{FiniteDifferences.var"#49#51"{Int64,Base.var"#64#65"{Base.var"#64#65"{Base.var"#64#65"{typeof(first),typeof(to_vec)},var"#7#8"},FiniteDifferences.var"#Real_from_vec#30"},Array{Float64,1}},Float64,UnitRange{Int64},Array{Float64,1},Float64}, ::typeof(Base.add_sum), ::IndexLinear, ::Base.OneTo{Int64}) at ./reduce.jl:403
 [6] _mapreduce_dim(::Function, ::Function, ::NamedTuple{(),Tuple{}}, ::Base.OneTo{Int64}, ::Colon) at ./reducedim.jl:312
 [7] #mapreduce#580 at ./reducedim.jl:307 [inlined]
 [8] mapreduce at ./reducedim.jl:307 [inlined]
 [9] _sum at ./reducedim.jl:657 [inlined]
 [10] #sum#584 at ./reducedim.jl:653 [inlined]
 [11] sum at ./reducedim.jl:653 [inlined]
 [12] fdm(::FiniteDifferences.Central{UnitRange{Int64},Array{Float64,1}}, ::FiniteDifferences.var"#49#51"{Int64,Base.var"#64#65"{Base.var"#64#65"{Base.var"#64#65"{typeof(first),typeof(to_vec)},var"#7#8"},FiniteDifferences.var"#Real_from_vec#30"},Array{Float64,1}}, ::Float64, ::Val{true}; condition::Int64, bound::Float64, eps::Float64, adapt::Int64, max_step::Float64) at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/methods.jl:270
 [13] fdm(::FiniteDifferences.Central{UnitRange{Int64},Array{Float64,1}}, ::FiniteDifferences.var"#49#51"{Int64,Base.var"#64#65"{Base.var"#64#65"{Base.var"#64#65"{typeof(first),typeof(to_vec)},var"#7#8"},FiniteDifferences.var"#Real_from_vec#30"},Array{Float64,1}}, ::Float64, ::Val{true}) at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/methods.jl:229
 [14] #fdm#25 at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/methods.jl:281 [inlined]
 [15] fdm at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/methods.jl:281 [inlined] (repeats 2 times)
 [16] #_#7 at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/methods.jl:93 [inlined]
 [17] Central at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/methods.jl:93 [inlined]
 [18] #48 at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/grad.jl:16 [inlined]
 [19] iterate at ./generator.jl:47 [inlined]
 [20] _collect(::Base.OneTo{Int64}, ::Base.Generator{Base.OneTo{Int64},FiniteDifferences.var"#48#50"{FiniteDifferences.Central{UnitRange{Int64},Array{Float64,1}},Base.var"#64#65"{Base.var"#64#65"{Base.var"#64#65"{typeof(first),typeof(to_vec)},var"#7#8"},FiniteDifferences.var"#Real_from_vec#30"},Array{Float64,1}}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:678
 [21] collect_similar(::Base.OneTo{Int64}, ::Base.Generator{Base.OneTo{Int64},FiniteDifferences.var"#48#50"{FiniteDifferences.Central{UnitRange{Int64},Array{Float64,1}},Base.var"#64#65"{Base.var"#64#65"{Base.var"#64#65"{typeof(first),typeof(to_vec)},var"#7#8"},FiniteDifferences.var"#Real_from_vec#30"},Array{Float64,1}}}) at ./array.jl:607
 [22] map(::Function, ::Base.OneTo{Int64}) at ./abstractarray.jl:2072
 [23] #jacobian#47 at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/grad.jl:15 [inlined]
 [24] jacobian at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/grad.jl:10 [inlined]
 [25] _j′vp(::FiniteDifferences.Central{UnitRange{Int64},Array{Float64,1}}, ::Function, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/grad.jl:79
 [26] j′vp(::FiniteDifferences.Central{UnitRange{Int64},Array{Float64,1}}, ::Function, ::Array{Float64,2}, ::Float64) at /Users/saxen/.julia/packages/FiniteDifferences/VHfqf/src/grad.jl:73
 [27] top-level scope at REPL[16]:1

 julia> j′vp(fdm, p -> p*x, seed, 3.0) # this scalar -> matrix function is fine
 (10.000000000000309,)

julia> j′vp(fdm, x->x^3.0, seed, x) # this matrix -> matrix function is also fine
([51.00000000000287 86.99999999999915; 66.9999999999964 110.9999999999985],)

[sethaxen]

Maybe part of what is going on here is that the output type of ^ is tightly linked to how close to an integer the power is. e.g.

julia> x^3.0
2×2 Array{Float64,2}:
 37.0   54.0
 81.0  118.0

julia> x^(3.0+sqrt(eps()))
2×2 Array{Complex{Float64},2}:
 37.0-1.83838e-9im   54.0+8.40924e-10im
 81.0+1.26139e-9im  118.0-5.76992e-10im

[sethaxen]

Actually, I'm not convinced anymore that this is a bug. This works

julia> j′vp(fdm, p -> complex(x)^real(p), complex(seed), complex(3.0))
(487.58111835736804 - 4.440892098500626e-14im,)

julia> j′vp(fdm, p -> complex(x)^imag(p), complex(seed), complex(3.0))
(8.326672684688674e-16 + 3.0165251948230645im,)

where I've made sure the x is complex to break the dependency of the output type on p's exact value, p is complex because its derivative will be complex, and seed is complex because the output is complex.

[willtebbutt]

Hi Seth. Yes, this is almost certainly what's going on. The doubling of the dimensionality is being caused by the output type being complex for non-integer exponent.

I'm not entirely sure what to do about this if I'm completely honest. From FiniteDifferences' perspective this function's output dimensionality depends on the value of its input, which is something it's definitely not designed to handle, so this probably isn't going to be a quick fix if you're really interested in computing the jpvp at exactly 3.0. If you're content to settle for computing it at 3.0 + sqrt(eps()), then you'll almost certainly be completely fine -- you'll just have to provide a complex-valued seed.

Just to be up-front, unless I'm missing something and it's actually really easy to fix this corner case, it's probably going to take a long time for us to implement a solution to this particular problem.

edit: I wrote this before seeing your last comment. I suspect that we came to similar conclusions about what's going on.

[nickrobinson251]

even if not making changes, may still warrant a mention in the docs as a "gotcha"?

[oxinabox]

Can we detect this case, or a subset of this case, and throw a clearer error?

[willtebbutt]

even if not making changes, may still warrant a mention in the docs as a "gotcha"?

Definitely.

Can we detect this case, or a subset of this case, and throw a clearer error?

We can definitely do something here. A simple option would be to highlight that

• the size of to_vec(output_of_function_being_differenced) has changed.
• this is probably due to either the type of the output of the function changing for different inputs, or the output size changing for different inputs.

Additionally, it's currently not completely straightforward to figure out exactly what inputs FiniteDifferences is testing your functions at. As such, it would probably make sense to add a debug mode to make it easier to figure out what inputs actually get evaluated, and what outputs they produced. You could probably already do this quite easily with the grid if you know what you're doing, but it's a bit of a faff.