vancouver-slides.md

Extensible Symbolic Programming in Julia

       Vancouver Julia Meetup

            04/20/2021
           Shashi Gowda

         PhD candidate in
Computational Science and Engineering
         Julia Lab, MIT

---

Outline

• Tour of Symbolics.jl
• Code generation
• Rule-based rewriting of expressions
• Interface and Interoperability

---

Symbolics.jl/SymbolicUtils.jl - Goals

• Easy to use (like SymPy/Mathematica)
• Fast!: use the fastest implementation
• Extensible: allow packages with term-like types to use the features
• Code-generation: turn symbolic results into executable, parallel Julia code.

---

Symbolics.jl - Tour

Expressions

using Symbolics
@variables t x y u(..)

x + y*x + x*y, x^2 + y

u(x,y,t)

---

Symbolics.jl - Tour

Expressions

julia> typeof(x+x*y)
Num # <- we will come back to this later

julia> typeof(x+x*y) <: Number
true

---

Symbolics.jl - Tour

Arrays of Expressions

julia> M = [t^2 + t + t^2  2t + 4t
            x + y + y + 2t  x^2 - x^2 + y^2]
2×2 Matrix{Num}:
  t + 2(t^2)   6t
 x + 2t + 2y  y^2

---

Symbolics.jl - Tour

Arrays of Expressions

julia> using LinearAlgebra

julia> det(M)
(t + 2(t^2))*(y^2) - (6t*(x + 2t + 2y))

julia> M \ M
2×2 Matrix{Num}:
 1  0
 0  1

julia> M * M

---

Symbolics.jl - Tour

Codegen!

using LinearAlgebra, SparseArrays

julia> f, f! = build_function(M, [x, y, t])

julia> f, f! = build_function(sparse(Diagonal(M)), [x, y, t])

---

Symbolics.jl - Tour

Differentiation

julia> Symbolics.jacobian([x + x*y, x^2 + y],
                          [x, y])
2×2 Matrix{Num}:
 1 + y  x
    2x  1

rosenbrock(X) = sum(1:length(X)-1) do i
    100 * (X[i+1] - X[i]^2)^2 + (1 - X[i])^2
end

@variables X[1:100]

rr = rosenbrock(X)

Symbolics.hessian(rr, X)

---

Symbolics.jl - Tour

Sparsity

using BenchmarkTools

julia> subsets = getindex.((X,), UnitRange.((1:5:100),(5:5:100)));

julia> exprs = rosenbrock.(subsets);

julia> @btime Symbolics.jacobian(exprs, X)
  1.078 s (4935567 allocations: 198.67 MiB)

julia> @btime Symbolics.sparsejacobian(exprs, X);
  23.134 ms (105622 allocations: 5.24 MiB)

julia> Symbolics.jacobian_sparsity(exprs, X)

---

Symbolics.jl - Tour

Sparsity

using BenchmarkTools

julia> subsets = getindex.((X,), UnitRange.((1:5:100),(5:5:100)))

julia> exprs = rosenbrock.(subsets);

julia> @btime Symbolics.jacobian(exprs, X)
  1.078 s (4935567 allocations: 198.67 MiB)

julia> @btime Symbolics.sparsejacobian(exprs, X);
  23.134 ms (105622 allocations: 5.24 MiB)

julia> Symbolics.jacobian_sparsity(exprs, X)

---

Symbolics.jl - Missing Parts

• Integration (definite integral)
• Polynomial division (easy)
• Groebner basis
• Polynomial equation solving
• Simplification under assumptions (Z3)

---

Under the hood - What is an Expression?

julia> @variables x y

julia> sin(cos(x))

julia> typeof(sin(cos(y)))

julia> using Symbolics: value

julia> value(sin(cos(x)))

---

Under the hood - What is an Expression?

julia> @variables x y

julia> value(sin(cos(x)))

julia> using SymbolicUtils

julia> SymbolicUtils.print_tree(stdout,
                value(sin(cos(x))))

---

Under the hood - What is an Expression?

Add -- Asymptotically better for storage

using SymbolicUtils

@syms x y

x + y

z1 = 10
for i=1:1000
    z1 += sin(rand([x, y]))
end

---

Under the hood - What is an Expression?

Mul -- Asymptotically better for storage

z2 = 10
for i=1:1000
    z2 *= sin(rand([x, y]))
end

---

Rule-based rewriting

using SymbolicUtils

@syms w z α::Real β::Real

r1 = @rule sin(2(~x)) => 2sin(~x)*cos(~x)

r1(sin(2z))

r1(sin(3z))

r1(sin(2*(w-z)))

r2 = @rule sin(~x + ~y) => sin(~x)*cos(~y) + cos(~x)*sin(~y);

r2(sin(α+β))

---

How do rules look at a term?

• istree
• operation
• arguments

function show_structure(x)
    @show typeof(x)
    @show istree(x)
    if istree(x)
        @show operation(x)
        @show arguments(x)
    end
end

show_structure(hypot(sin(x), cos(x)))
show_structure(z1)
show_structure(z2)

---

Benefits of the Interface

• Using the term interface allows term-rewriting of any type.
• It also allows packages like Metatheory to be an alternative term-rewriting system

---

Thank you!

(Thanks to Yingbo Ma for helping me prepare this talk!)

Work presented here is by:

• me
• Yingbo Ma
• Chris Rackauckas
• Mason Protter
• many other contributors