inconsistent order of terms causes incorrect derivative calculation

[syitong]

using Symbolics
@variables f(..) h(..) A[1:8, 1:6] C[1:8, 1:20] Z(..)
foo = f(A[6, 2])*Differential(A[6, 2])(h(A[6, 2]))*C[2, 15] + Differential(A[6, 2])(h(A[6, 2]))*C[2, 7]*Z(A[6,2])
Symbolics.derivative(foo, Differential(A[6,2])(h(A[6,2])))

# ans: Differential(A[6, 2])(h(A[6, 2]))*C[2, 7] + C[2, 15]*f(A[6, 2])
# expected ans: f(A[6,2])*C[2,15]+C[2,7]*Z(A[6,2])

Julia version 1.8; Symbolics version 4.13.

It seems the bug is caused by the diff.jl line 273-274

        args = map(a->expand_derivatives(a, false), arguments(O))
        O1 = operation(O)(args...)

These lines will change the order of terms in the product, which becomes inconsistent to the order of occurrences in line 237.

The example code looks weird because I cannot reproduce the bug with other simpler case. I do not understand the ordering logic of operation action. It seems that the order of terms is not always messed up and thus the bug in the expand_derivatives does not always show up.

[bowenszhu]

Yes! You are right. The order of terms are inconsistent. Thanks for looking into the source code.
This is caused by

• <ₑ not fully reliable SymbolicUtils.jl#483

I've always been thinking JuliaSymbolics/SymbolicUtils.jl#483 would create implicit bugs but haven't found a good example. Thanks for providing one.

Symbolics.jl/src/diff.jl

Line 163 in 4682772

function expand_derivatives(O::Symbolic, simplify=false; occurances=nothing)

The problem is that the 3rd input argument occurances of expand_derivatives assumes some ordering of arguments. SymbolicUtils.arguments for Add and Mul uses SymbolicUtils.:<ₑ to sort arguments, but SymbolicUtils.:<ₑ has a problem as reported in JuliaSymbolics/SymbolicUtils.jl#483.

In the following, the versions are Julia v1.8.3, Symbolics v4.13.0, SymbolicUtils v0.19.11.

To bypass the problem of SymbolicUtils.:<ₑ, i.e. JuliaSymbolics/SymbolicUtils.jl#483, I firstly define an arbitrary ordering for the terms that are present in the test symbolic expression foo, and then redefine SymbolicUtils.arguments for Add and Mul.

using Symbolics
@variables f(..) h(..) A[1:8, 1:6] C[1:8, 1:20] Z(..)

# an arbitrary ordering
ord_vec = Symbolics.value.([
                               A[6, 2],
                               f(A[6, 2]),
                               h(A[6, 2]),
                               Differential(A[6, 2])(h(A[6, 2])),
                               C[2, 15],
                               C[2, 7],
                               Z(A[6, 2]),
                               f(A[6, 2]) * Differential(A[6, 2])(h(A[6, 2])) * C[2, 15],
                               Differential(A[6, 2])(h(A[6, 2])) * C[2, 7] * Z(A[6, 2]),
                           ])
ord_dict = Dict(t => i for (i, t) in enumerate(ord_vec))
function <ₐ(a, b) # new ordering
    aᵢ = get(ord_dict, a, 0)
    bᵢ = get(ord_dict, b, 0)
    if aᵢ > 0 && bᵢ > 0
        aᵢ < bᵢ
    else
        SymbolicUtils.:<ₑ(a, b) # call original ordering
    end
end
function SymbolicUtils.arguments(a::SymbolicUtils.Add)
    a.sorted_args_cache[] !== nothing && return a.sorted_args_cache[]
    args = sort!([v * k for (k, v) in a.dict], lt = <ₐ) # new ordering
    a.sorted_args_cache[] = iszero(a.coeff) ? args : vcat(a.coeff, args)
end
function SymbolicUtils.arguments(a::SymbolicUtils.Mul)
    a.sorted_args_cache[] !== nothing && return a.sorted_args_cache[]
    args = sort!([SymbolicUtils.unstable_pow(k, v) for (k, v) in a.dict], lt = <ₐ) # new ordering
    a.sorted_args_cache[] = isone(a.coeff) ? args : vcat(a.coeff, args)
end

foo = f(A[6, 2]) * Differential(A[6, 2])(h(A[6, 2])) * C[2, 15] +
      Differential(A[6, 2])(h(A[6, 2])) * C[2, 7] * Z(A[6, 2])
res = Symbolics.derivative(foo, Differential(A[6, 2])(h(A[6, 2])))

The result is f(A[6, 2])*C[2, 15] + C[2, 7]*Z(A[6, 2]) which is correct!

[karlwessel]

This seems to have been fixed at some point. For me with Symbolics 6.11 and SymbolicUtils 3.6 it gives the correct result:

julia> using Symbolics
julia> @variables f(..) h(..) A[1:8, 1:6] C[1:8, 1:20] Z(..)
5-element Vector{Any}:
 f⋆
 h⋆
 A[1:8,1:6]
 C[1:8,1:20]
 Z⋆
julia> foo = f(A[6, 2])*Differential(A[6, 2])(h(A[6, 2]))*C[2, 15] + Differential(A[6, 2])(h(A[6, 2]))*C[2, 7]*Z(A[6,2])
C[2, 15]*Differential(A[6, 2])(h(A[6, 2]))*f(A[6, 2]) + C[2, 7]*Differential(A[6, 2])(h(A[6, 2]))*Z(A[6, 2])
julia> Symbolics.derivative(foo, Differential(A[6,2])(h(A[6,2])))
C[2, 15]*f(A[6, 2]) + C[2, 7]*Z(A[6, 2])

[bowenszhu]

Thank you! @karlwessel