Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Open
sethaxen opened this issue Apr 29, 2020 · 6 comments
Open

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

sethaxen opened this issue Apr 29, 2020 · 6 comments

Comments

@sethaxen
Copy link
Member

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],)
@oxinabox oxinabox added the bug Something isn't working label Apr 29, 2020
@sethaxen
Copy link
Member Author

sethaxen commented Apr 29, 2020 •
edited
Loading

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
Copy link
Member Author

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
Copy link
Member

willtebbutt commented Apr 29, 2020 •
edited
Loading

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
Copy link
Contributor

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

@oxinabox oxinabox removed the bug Something isn't working label Apr 29, 2020
@oxinabox
Copy link
Member

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

@willtebbutt
Copy link
Member

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.

@oxinabox oxinabox mentioned this issue Nov 26, 2024
Open
@willtebbutt willtebbutt mentioned this issue Apr 30, 2020
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@willtebbutt @oxinabox @sethaxen @nickrobinson251