Can't take a gradient using Symbolic vector

Cade_Armstrong · November 28, 2023, 10:02pm

I’m running into some trouble when taking a gradient of a symbolic function using a symbolic vector.

using Symbolics

@variables x[1:2]

xₒ = [1;1]

f = (x[1] + 2*x[2])^3

f_fun = build_function(f, x)

vals = eval(f_fun)

vals(xₒ)

df = Symbolics.gradient(f, x)

This code works just fine up until I try to use Symbolics.gradient. At this point, I get this error:

ERROR: MethodError: no method matching similar(::Type{Vector{Num}}, ::Tuple{UnitRange{Int64}})

Does anyone have any experience with why this might be happening and how to fix it? Any help would be greatly appreciated.

EDIT:

So I tried this:

df = Symbolics.gradient(f,[x[1],x[2]])

which worked just fine. However, when I run isequal(x, [x[1],x[2]]), it returns true. So now I’m extra confused. Any advice on why one way works and the other doesn’t would be great!

contradict · November 28, 2023, 10:16pm

What version of Symbolics are you using?

julia> using Symbolics

julia> xₒ = [1;1]
2-element Vector{Int64}:
 1
 1

julia> @variables x[1:2]
1-element Vector{Symbolics.Arr{Num, 1}}:
 x[1:2]

julia> xₒ = [1;1]
2-element Vector{Int64}:
 1
 1

julia> f = (x[1] + 2*x[2])^3
(x[1] + 2x[2])^3

julia> f_fun = build_function(f, x)
:(function (x,)
      #= /home/russel/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:373 =#
      #= /home/russel/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:374 =#
      #= /home/russel/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:375 =#
      (^)((+)((getindex)(x, 1), (*)(2, (getindex)(x, 2))), 3)
  end)

julia> vals = eval(f_fun)
#7 (generic function with 1 method)

julia> vals(xₒ)
27

julia> df = Symbolics.gradient(f, x)
2-element OffsetArray(::Vector{Num}, 1:2) with eltype Num with indices 1:2:
 3((x[1] + 2x[2])^2)
 6((x[1] + 2x[2])^2)

(MTK) pkg> st
Project MTK v0.1.0
Status `~/Desktop/Julia/MTK/Project.toml`
  [b0b7db55] ComponentArrays v0.15.5
  [82cc6244] DataInterpolations v4.5.0
  [0c46a032] DifferentialEquations v7.11.0
  [961ee093] ModelingToolkit v8.73.1
  [16a59e39] ModelingToolkitStandardLibrary v2.3.4
  [8913a72c] NonlinearSolve v2.8.2
  [1dea7af3] OrdinaryDiffEq v6.59.3
  [d96e819e] Parameters v0.12.3
  [91a5bcdd] Plots v1.39.0
  [731186ca] RecursiveArrayTools v2.38.10
  [0c5d862f] Symbolics v5.10.0

Cade_Armstrong · November 28, 2023, 10:32pm

I have the same version that you have. Interesting that it is working for you. Do you think that another package that you have installed could be making it work for you in some weird way?

contradict · November 28, 2023, 10:50pm

Works in a temporary environment with just Symbolics as well. How about Julia versions?

(jl_BnaUHa) pkg> st
Status `/tmp/jl_BnaUHa/Project.toml`
  [0c5d862f] Symbolics v5.10.0

julia> versioninfo()
Julia Version 1.9.4
Commit 8e5136fa297 (2023-11-14 08:46 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
  Threads: 1 on 8 virtual cores
Environment:
  JULIA_REVISE_INCLUDE = 1

trung · November 28, 2023, 10:57pm

I am reproducing the same issue as @Cade_Armstrong and I am currently running Julia 1.9.4 and Symbolics 5.10.0. My platform is on Mac OS X 13.5.1 and the hardware is a M2 MacBook Pro.

Dan · November 28, 2023, 11:00pm

Error replicated also.

julia> df = Symbolics.gradient(f, x[1:2])
2-element Vector{Num}:
 0
 0

works but it gets a wrong answer.

julia> Symbolics.gradient(f, [x[1],x[2]])
2-element Vector{Num}:
 3((x[1] + 2x[2])^2)
 6((x[1] + 2x[2])^2)

is more like it.

contradict · November 28, 2023, 11:18pm

Wow, that’s quite a diversity of results! Can someone who gets the exception post a full backtrace?

nsajko · November 28, 2023, 11:43pm

julia> df = Symbolics.gradient(f, x)
ERROR: MethodError: no method matching similar(::Type{Vector{Num}}, ::Tuple{UnitRange{Int64}})

Closest candidates are:
  similar(::VSCodeServer.JuliaInterpreter.Compiled, ::Any)
   @ VSCodeServer ~/.vscode-oss/extensions/julialang.language-julia-1.60.2-universal/scripts/packages/JuliaInterpreter/src/types.jl:7
  similar(::Type{A}, ::Type{T}, ::StaticArraysCore.Size{S}) where {A<:Array, T, S}
   @ StaticArrays ~/.julia/packages/StaticArrays/yXGNL/src/abstractarray.jl:136
  similar(::Type{A}, ::StaticArraysCore.Size{S}) where {A<:AbstractArray, S}
   @ StaticArrays ~/.julia/packages/StaticArrays/yXGNL/src/abstractarray.jl:124
  ...

Stacktrace:
 [1] _array_for(::Type{Num}, ::Base.HasShape{1}, axs::Tuple{UnitRange{Int64}})
   @ Base ./array.jl:662
 [2] _array_for(::Type{Num}, itr::Symbolics.Arr{Num, 1}, isz::Base.HasShape{1})
   @ Base ./array.jl:665
 [3] gradient(O::Num, vars::Symbolics.Arr{Num, 1}; simplify::Bool)
   @ Symbolics ~/.julia/packages/Symbolics/gBKZv/src/diff.jl:435
 [4] gradient(O::Num, vars::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/gBKZv/src/diff.jl:434
 [5] top-level scope
   @ REPL[12]:1

(@v1.11) pkg> st
Status `~/.julia/environments/v1.11/Project.toml`
  [233ec0c9] AMRVW v1.2.1
  [9b6a8646] AllocCheck v0.1.0
  [6e4b80f9] BenchmarkTools v1.3.2
  [bcf9a6e7] BugReporting v0.3.4
  [31c24e10] Distributions v0.25.103
  [7c1d4256] DynamicPolynomials v0.5.3 `~/tmp/jl/DynamicPolynomials.jl`
  [af72fe11] FindMinimaxPolynomial v0.4.0 `~/src/gitlab.com/nsajko/FindMinimaxPolynomial`
  [c91e804a] Gadfly v1.4.0
  [c3a54625] JET v0.8.20
  [b8f27783] MathOptInterface v1.22.0
  [102ac46a] MultivariatePolynomials v0.5.3 `~/tmp/jl/MultivariatePolynomials.jl`
  [d8a4904e] MutableArithmetics v1.3.3 `~/tmp/jl/MutableArithmetics.jl`
  [cae6f51d] OptimalSortingNetworks v1.0.0 `~/src/gitlab.com/nsajko/OptimalSortingNetworks`
  [e4faabce] PProf v3.1.0
  [14b8a8f1] PkgTemplates v0.7.46
  [f27b6e38] Polynomials v4.0.6
  [0c5d862f] Symbolics v5.10.0
  [6dd1b50a] Tulip v0.9.5
  [ce3070e3] TupleShuffling v1.0.0 `~/src/gitlab.com/nsajko/TupleShuffling`
  [4fd249ec] TupleSorting v1.0.0 `~/src/gitlab.com/nsajko/TupleSorting`
  [afbbf031] TypedPolynomials v0.4.1
  [99140e25] UnsafeAssume v1.0.1 `~/src/gitlab.com/nsajko/UnsafeAssume`

nsajko · November 28, 2023, 11:53pm

github.com

JuliaSymbolics/Symbolics.jl/blob/8729f1ccf8978fb6df2480e9ae4c87f05c3d66e4/src/diff.jl#L434-L436


      
          function gradient(O, vars::AbstractVector; simplify=false)
              Num[Num(expand_derivatives(Differential(v)(value(O)),simplify)) for v in vars]
          end

The code hasn’t been changed in two years, so I guess this is a Symbolics.Arr bug?

Cade_Armstrong · November 29, 2023, 12:07am

Yep, it looks like its a bug that was reported at the beginning of October. Here’s the link:

github.com/JuliaSymbolics/Symbolics.jl

Symbolic array incompatible with Symbolics.gradient and Symbolics.hessian

opened 09:27AM - 06 Oct 23 UTC

Krebbekx

Symbolic arrays seem to be incompatible with the gradient in the following sense…. ```julia using Symbolics @variables x[1:3] fun = x[1] gradfunx = Symbolics.gradient(fun, x) # does not work ``` The above code raises the error: ``` MethodError: no method matching similar(::Type{Vector{Num}}, ::Tuple{UnitRange{Int64}}) ``` A workaround I use is to unpack `x` and put it into an array by hand: ```julia y = [x...] gradfuny = Symbolics.gradient(fun, y) # works ``` I think it is desirable that `Symbolics` supports the usage of `gradient` as in `gradfunx` above, because the same *does* work for the Jacobian: ```julia funarray = [x[1]] # make an array of functions jacfunx = Symbolics.jacobian(funarray, x) # works jacfunx2 = Symbolics.jacobian(funarray, [x...]) # works, and gives the same answer ``` Furthermore, when using the `hessian`, there is another issue. ```julia fun = x[1]^3 hessfunx = Symbolics.hessian(fun, x) ``` outputs ``` 3×3 Matrix{Num}: 0 0 0 0 0 0 0 0 0 ``` whereas ```julia fun = x[1]^3 hessfunx = Symbolics.hessian(fun, [x...]) ``` outputs ``` 3×3 Matrix{Num}: 6x[1] 0 0 0 0 0 0 0 0 ``` This means that you get a different answer, depending on whether you use a symbolic array, or an array of symbolic expressions.

In the bug report, the author gives a pretty simple work around.

using Symbolics

@variables x[1:2]

f = (x[1] + 2*x[2])^3

y = [x...]

df = Symbolics.gradient(f, y)

seems to work fine.