Symbolics build_function issues on complicated expressions

Specific Domains Modelling & Simulations
1

I am trying to use Symbolics to define and manipulate my problem before moving to numerical methods to actually solve my problem.

After initial manipulations, I want to replace all symbolic variables with their actual values. This worked well for me as long as I used matrices of scalar symbolic variables. I am now trying to move to symbolic arrays, and substitute no longer substitutes all symbolic variables. I tried build_function which managed better, but which also seems to only work on sufficiently simple expressions.

The behaviour I expected is that all symbolic variables are replaced by numerical values in all cases.

I fail to solve my problem using either of substitute or build_function, with build_function managing more complicated expressions, but not the ones I need it for.

I would appreciate any help.

The following is the smallest test case I have showing the issue:

import Pkg
Pkg.activate(".")
Pkg.add("Symbolics")
Pkg.status(outdated=true)

import Symbolics
import Symbolics: @variables, get_variables, substitute, build_function
@variables K[1:2, 1:2]
exprs = [ -K[1,1] -(K[1,1]) (-K)[1,1] ]
subs = Dict(K[1,1] => 1)
for expr in exprs
    println("------------------------")
    @show expr get_variables(expr) 
    try 
        @show substitute(expr, subs)
    catch e
        @show e
    end
    f = build_function(expr, K, expression = false)
    @show f([1 2; 3 4])
end

@variables A[1:2, 1:2] B[1:2, 1:1] P[1:2, 1:2]

Pre1 = A*[1.0 0.0 0.0; 0.0 1.0 0.0]
Pre2 = broadcast(+, A*[1.0 0.0 0.0; 0.0 1.0 0.0], B*[0.0 0.0 1.0])
Post = Symbolics.Num[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 (broadcast(-, K))[1, 1] (broadcast(-, K))[1, 2]]

for expr in [ Pre1, Pre2, Post, Pre1*Post, Pre2*Post ]
    println("---------------------------")
    @show expr
    fs = build_function(expr, A, B, K, P, expression=false)
    for f in fs
        try 
            value = f(randn(2,2), randn(2,1), randn(1,2), randn(2,2))
            @show value unique(vcat(get_variables.(value)...))
        catch e
        end
    end        
end

The output I get is

Status `~/Workspace/Project.toml`
⌅ [0c5d862f] Symbolics v5.14.1 (<v5.16.1): julia
------------------------
expr = -K[1, 1]
get_variables(expr) = Any[K[1, 1]]
substitute(expr, subs) = -1
f([1 2; 3 4]) = -1
------------------------
expr = -K[1, 1]
get_variables(expr) = Any[K[1, 1]]
substitute(expr, subs) = -1
f([1 2; 3 4]) = -1
------------------------
expr = (broadcast(-, K))[1, 1]
get_variables(expr) = Any[(broadcast(-, K))[1, 1]]
e = MethodError(Core.kwcall, ((metadata = nothing,), SymbolicUtils.similarterm, broadcast(-, K), broadcast, Any[-, K], AbstractMatrix{Real}), 0x0000000000008337)
f([1 2; 3 4]) = -1
---------------------------
expr = (A*[1.0 0.0 0.0; 0.0 1.0 0.0])[Base.OneTo(2),Base.OneTo(3)]
value = [0.8426921154384049 1.3855828533066428 0.0; 2.245360817652516 0.45372334877928383 0.0]
unique(vcat(get_variables.(value)...)) = Any[]
---------------------------
expr = (broadcast(+, A*[1.0 0.0 0.0; 0.0 1.0 0.0], B*[0.0 0.0 1.0]))[Base.OneTo(2),Base.OneTo(3)]
value = [1.3634561500672868 -1.1694430577347337 1.6690534965623611; 1.581800165321597 1.1989213984343103 0.8634154223268747]
unique(vcat(get_variables.(value)...)) = Any[]
---------------------------
expr = Symbolics.Num[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 (broadcast(-, K))[1, 1] (broadcast(-, K))[1, 2]]
value = [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 -0.4432466306968838 0.3592732607253264]
unique(vcat(get_variables.(value)...)) = Any[]
---------------------------
expr = ((A*[1.0 0.0 0.0; 0.0 1.0 0.0])*Symbolics.Num[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 (broadcast(-, K))[1, 1] (broadcast(-, K))[1, 2]])[Base.OneTo(2),Base.OneTo(4)]
value = Symbolics.Num[-1.873960325836164 1.416262070648973 0.0 0.0; -1.4423420994229887 -0.31731200550770267 0.0 0.0]
unique(vcat(get_variables.(value)...)) = Any[]
---------------------------
expr = (broadcast(+, A*[1.0 0.0 0.0; 0.0 1.0 0.0], B*[0.0 0.0 1.0])*Symbolics.Num[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 (broadcast(-, K))[1, 1] (broadcast(-, K))[1, 2]])[Base.OneTo(2),Base.OneTo(4)]
value = Symbolics.Num[0.45315331007912363 -1.367295212832985 0.20383264427229145(broadcast(-, K))[1, 1] 0.20383264427229145(broadcast(-, K))[1, 2]; -1.3151182967142874 -1.1121032853446047 -0.25979763416747165(broadcast(-, K))[1, 1] -0.25979763416747165(broadcast(-, K))[1, 2]]
unique(vcat(get_variables.(value)...)) = Any[(broadcast(-, K))[1, 1], (broadcast(-, K))[1, 2]]

The first part of the output shows that substitute works for me in some cases, but fails as soon as broadcast is introduced into the expressions.

The second part of the output shows that build_function manages to replace all symbolic variables in the expressions Pre1, Pre2, Post and Pre1*Post, but fails in Pre2*Post.

2

Tjena Johan, I think the unfortunate answer here is that array variables are very buggy and mostly do not work as expected. The list of known issues is rather long, but there are likely several issues not reported here as well.

3

Tjena, och tack.

That is unfortunate, for me. I was hoping to use build_function to quickly repeat the same analysis when parameter values change, but if I have a large number of elements in my matrices, then I end up building functions that take a large number of scalar arguments, instead of a small number of array arguments.

(As you will already know from my name and the variable names A,B,K,P appearing, I am analysing stability properties of control systems.)

The problems “it works in simple expressions but not on complex ones” and “it fails with broadcast” appears in several bug reports and it seems to me like it should be relatively simple to solve. It also seems as if the, relatively simple, solution, should enable quite a bit of functionality to become reliable…

Would you agree?

4

You should stil be able to build function that take arrays as arguments. If you call array_of_syms = collect(symbolic_array) you can use array_of_syms as input argument to build_funciton,, which builds a function that expects a numeric array of the same dimensions.

In general, most things “should work”, if you scalarize the symbolic array expressions using collect. The problem you may get instead is that scalar expressions may be come very large, and the generated function may not perform very well and/or take a long time to compile if the scalar expression is very large.

5

My interpretation of your comment is that I should use the following:

    fs = build_function(Symbolics.scalarize(expr), collect(A), collect(B), collect(K), collect(P), expression=false)

That fails for me, because scalarize fails on the same subexpression of the same testcase where build_function failed.

The call to scalarize in

expr = (broadcast(+, A*[1.0 0.0 0.0; 0.0 1.0 0.0], B*[0.0 0.0 1.0])*Symbolics.Num[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 (broadcast(-, K))[1, 1] (broadcast(-, K))[1, 2]])[Base.OneTo(2),Base.OneTo(4)]
Symbolics.scalarize(expr)

returns

Symbolics.Num[A[1, 1] A[1, 2] B[1, 1]*(broadcast(-, K))[1, 1] B[1, 1]*(broadcast(-, K))[1, 2]; A[2, 1] A[2, 2] B[2, 1]*(broadcast(-, K))[1, 1] B[2, 1]*(broadcast(-, K))[1, 2]]

which did not scalarize the appearances of K. The function built by build_function throws an UndefVarError(:K), probably because it expects scalars like K[1,1], but not the K that appears as an argument to broadcast.

6

Try calling collect on K immediately after creating it so you never work with the symbolic arrays at all

import Symbolics
import Symbolics: @variables, get_variables, substitute, build_function
@variables K[1:2, 1:2]
K = collect(K)
exprs = [ -K[1,1] -(K[1,1]) (-K)[1,1] ]
subs = Dict(K[1,1] => 1)
for expr in exprs
    println("------------------------")
    @show expr get_variables(expr) 
    try 
        @show substitute(expr, subs)
    catch e
        @show e
    end
    f = build_function(expr, K, expression = false)
    @show f([1 2; 3 4])
end

@variables A[1:2, 1:2] B[1:2, 1:1] P[1:2, 1:2]
A,B,P = collect.((A,B,P))

Pre1 = A*[1.0 0.0 0.0; 0.0 1.0 0.0]
Pre2 = broadcast(+, A*[1.0 0.0 0.0; 0.0 1.0 0.0], B*[0.0 0.0 1.0])
Post = Symbolics.Num[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 (broadcast(-, K))[1, 1] (broadcast(-, K))[1, 2]]

for expr in [ Pre1, Pre2, Post, Pre1*Post, Pre2*Post ]
    println("---------------------------")
    @show expr
    fs = build_function(expr, A, B, K, P, expression=false)
    for f in fs
        try 
            value = f(randn(2,2), randn(2,1), randn(1,2), randn(2,2))
            @show value unique(vcat(get_variables.(value)...))
        catch e
        end
    end        
end
7

Excellent!

Thanks!