Jacobian shape error with FillArrays.Zeros vector

[matbesancon]

I think this repo is where the issue should be posted, otherwise I'll transfer it to FillArrays

Versions:

  [1a297f60] FillArrays v0.11.7
  [26cc04aa] FiniteDifferences v0.12.12

the issue:

julia> v = rand(3);
julia> FiniteDifferences.jacobian(ffdm, x -> zeros(3), v)[1]
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> FiniteDifferences.jacobian(ffdm, x -> FillArrays.Zeros(3), v)[1]
1×3 Matrix{Float64}:
 0.0  0.0  0.0

[matbesancon]

The bug seems absent using 0.12.0, I'm trying to bisect now

[matbesancon]

the bug seems to have been introduced in FiniteDifferences 0.12.3, .2 was not failing similarly

[nickrobinson251]

suggests it might have been introduced by #156

[mzgubic]

What's ffdm?

[matbesancon]

const ffdm = FiniteDifferences.forward_fdm(5, 1)

[mzgubic]

The issue is that zeros and Zeros have a different meaning.

zeros(3) is just a convenience function to create a vector which values happen to be zero, and which has 3 "degrees of freedom", i.e. it can be perturbed in any of the three entries, and it is still a vector.

Zeros(3) struct on the other hand represents a structure where all the entries are identically zero. It is not possible to perturb the first element by epsilon and still represent the underlying vector by a Zeros struct.

What is happening under the hood is that Zeros is being treated as a struct with one field (3, representing the size of the vector), indeed introduced by #156 as @nickrobinson251 suggests. That is why there is only row in the jacobian, and not three as with zeros.

Note that even with 0.12.3 the to_vec does not really work, i.e. you can't call back on the vector and get the Zeros back:

julia> v, back = to_vec(Zeros(3))
([0.0, 0.0, 0.0], FiniteDifferences.var"#Vector_from_vec#32"{Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}, Fill{FiniteDifferences.var"#Real_from_vec#22", 1, Tuple{Base.OneTo{Int64}}}, Fill{Vector{Float64}, 1, Tuple{Base.OneTo{Int64}}}}(3-element Zeros{Float64}, 3-element Fill{FiniteDifferences.var"#Real_from_vec#22"}: entries equal to Real_from_vec, 3-element Fill{Vector{Float64}}: entries equal to [0.0]))

julia> back(v)
ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}
Closest candidates are:
  convert(::Type{T}, ::LinearAlgebra.Factorization) where T<:AbstractArray at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:58
  convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:14
  convert(::Type{T}, ::T) where T at essentials.jl:205
  ...
Stacktrace:
 [1] oftype(x::Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}, y::Vector{Float64})
   @ Base ./essentials.jl:375
 [2] (::FiniteDifferences.var"#Vector_from_vec#32"{Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}, Fill{FiniteDifferences.var"#Real_from_vec#22", 1, Tuple{Base.OneTo{Int64}}}, Fill{Vector{Float64}, 1, Tuple{Base.OneTo{Int64}}}})(x_vec::Vector{Float64})
   @ FiniteDifferences ~/JuliaEnvs/test/bug/dev/FiniteDifferences/src/to_vec.jl:49
 [3] top-level scope
   @ REPL[18]:1

Could you tell me a bit more about where you are using this? I assume that it is not easy to swap using Zeros for zeros?

[mzgubic]

If it is just for testing what you have in https://github.com/matbesancon/MathOptSetDistances.jl/pull/46/checks?check_run_id=2865353507, a quick and dirty fix would something like

function to_vec(x::FillArrays.Zeros)
    v = collect(x)
    return v, _ -> error("can't create `Zeros` from a vector")
end

which will work as much as it worked before 0.12.3 (i.e. turning Zeros to a vector is fine, but going back errors)

[matbesancon]

thanks that will work. This is indeed the use case. However, this still means we produce a different behavior from what would be expected, except if the to_vec implementation is pirated in the package source

[mzgubic]

I may have misread the tests, but my impression was that in the test

dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), v, s)
grad_fdm1 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'
grad_fdm2 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(MOD.DefaultDistance(), x, s), v)[1]'

dΠ is the gradient that is computed in the package source, and you are testing it with finite differences. Computing dΠ does not depend on to_vec, right?

If you do want the dΠ to be a 3x3 matrix then just testing it with a hack seems fine?

[matbesancon]

indeed I confused myself yesterday evening, the type piracy can live in the tests and works there

[mzgubic]

I guess we can close this, unless anyone considers this to be unwanted behaviour?