Symbolic calculations in Julia

Symbolic calculations in Julia#

Symbolics.jl is a computer Algebra System (CAS) for Julia. The symbols are number-like and follow Julia semantics so we can put them into a regular function to get a symbolic counterpart. Symbolics.jl is the backbone of ModelingToolkit.jl.

Source:

Caveats about Symbolics.jl#

  1. Symbolics.jl can only handle traceble, quasi-static expressions. However, some expressions are not quasi-static e.g. factorial. The number of operations depends on the input value.

Use @register to make it a primitive function.

  1. Some code paths is untraceable, such as conditional statements: if…else…end.

You can use ifelse(cond, ex1, ex2) to make it traceable.

Basic operations#

  • latexify

  • derivative

  • gradient

  • jacobian

  • substitute

  • simplify

using Symbolics
using Latexify

@variables x y
x^2 + y^2
\[ \begin{equation} x^{2} + y^{2} \end{equation} \]

You can use Latexify.latexify() to see the LaTeX code.

A = [
    x^2+y 0 2x
    0 0 2y
    y^2+x 0 0
]

latexify(A)
\[\begin{split}\begin{equation} \left[ \begin{array}{ccc} y + x^{2} & 0 & 2 x \\ 0 & 0 & 2 y \\ x + y^{2} & 0 & 0 \\ \end{array} \right] \end{equation}\end{split}\]

Derivative: Symbolics.derivative(expr, variable)

Symbolics.derivative(x^2 + y^2, x)
\[ \begin{equation} 2 x \end{equation} \]

Gradient: Symbolics.gradient(expr, [variables])

Symbolics.gradient(x^2 + y^2, [x, y])
\[\begin{split} \begin{equation} \left[ \begin{array}{c} 2 x \\ 2 y \\ \end{array} \right] \end{equation} \end{split}\]

Jacobian: Symbolics.jacobian([exprs], [variables])

Symbolics.jacobian([x^2 + y^2; y^2], [x, y])
\[\begin{split} \begin{equation} \left[ \begin{array}{cc} 2 x & 2 y \\ 0 & 2 y \\ \end{array} \right] \end{equation} \end{split}\]

Substitute: Symbolics.substitute(expr, mapping)

Symbolics.substitute(sin(x)^2 + 2 + cos(x)^2, Dict(x => y^2))
\[ \begin{equation} 2 + \sin^{2}\left( y^{2} \right) + \cos^{2}\left( y^{2} \right) \end{equation} \]
Symbolics.substitute(sin(x)^2 + 2 + cos(x)^2, Dict(x => 1.0))
\[ \begin{equation} 3 \end{equation} \]

Simplify: Symbolics.simplify(expr)

Symbolics.simplify(sin(x)^2 + 2 + cos(x)^2)
\[ \begin{equation} 3 \end{equation} \]

This expression gets automatically simplified because it’s always true

2x - x
\[ \begin{equation} x \end{equation} \]
ex = x^2 + y^2 + sin(x)
isequal(2ex, ex + ex)

false

isequal(simplify(2ex), simplify(ex + ex))

true

ex / ex
\[ \begin{equation} 1 \end{equation} \]

Symbolic integration: use SymbolicNumericIntegration.jl. The Youtube video by doggo dot jl gives a concise example.

Custom functions used in Symbolics#

With @register, the rand() in foo() will be evaluated as-is and will not be expanded by Symbolics.jl.

By default, new functions are traced to the primitives and the symbolic expressions are written on the primitives. However, we can expand the allowed primitives by registering new functions.

More number types#

Complex number

@variables z::Complex
\[\begin{split} \begin{equation} \left[ \begin{array}{c} \mathrm{real}\left( z \right) + \mathrm{imag}\left( z \right) \mathit{i} \\ \end{array} \right] \end{equation} \end{split}\]

Array types with subscript

@variables xs[1:18]

1-element Vector{Symbolics.Arr{Symbolics.Num, 1}}:
 xs[1:18]

xs[1]
\[ \begin{equation} \mathtt{xs}_{1} \end{equation} \]

Explicit vector form

collect(xs)
\[\begin{split} \begin{equation} \left[ \begin{array}{c} \mathtt{xs}_{1} \\ \mathtt{xs}_{2} \\ \mathtt{xs}_{3} \\ \mathtt{xs}_{4} \\ \mathtt{xs}_{5} \\ \mathtt{xs}_{6} \\ \mathtt{xs}_{7} \\ \mathtt{xs}_{8} \\ \mathtt{xs}_{9} \\ \mathtt{xs}_{10} \\ \mathtt{xs}_{11} \\ \mathtt{xs}_{12} \\ \mathtt{xs}_{13} \\ \mathtt{xs}_{14} \\ \mathtt{xs}_{15} \\ \mathtt{xs}_{16} \\ \mathtt{xs}_{17} \\ \mathtt{xs}_{18} \\ \end{array} \right] \end{equation} \end{split}\]

Operations on arrays are supported

sum(collect(xs))
\[ \begin{equation} \mathtt{xs}_{1} + \mathtt{xs}_{10} + \mathtt{xs}_{11} + \mathtt{xs}_{12} + \mathtt{xs}_{13} + \mathtt{xs}_{14} + \mathtt{xs}_{15} + \mathtt{xs}_{16} + \mathtt{xs}_{17} + \mathtt{xs}_{18} + \mathtt{xs}_{2} + \mathtt{xs}_{3} + \mathtt{xs}_{4} + \mathtt{xs}_{5} + \mathtt{xs}_{6} + \mathtt{xs}_{7} + \mathtt{xs}_{8} + \mathtt{xs}_{9} \end{equation} \]

Example: Rosenbrock function#

Wikipedia: https://en.wikipedia.org/wiki/Rosenbrock_function We use the vector form of Rosenbrock function.

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

rosenbrock (generic function with 1 method)

The function is at minimum when xs are all one’s

rosenbrock(ones(100))

0.0

N = 20
@variables xs[1:N]

1-element Vector{Symbolics.Arr{Symbolics.Num, 1}}:
 xs[1:20]

A full list of vector components

xs = collect(xs)
\[\begin{split} \begin{equation} \left[ \begin{array}{c} \mathtt{xs}_{1} \\ \mathtt{xs}_{2} \\ \mathtt{xs}_{3} \\ \mathtt{xs}_{4} \\ \mathtt{xs}_{5} \\ \mathtt{xs}_{6} \\ \mathtt{xs}_{7} \\ \mathtt{xs}_{8} \\ \mathtt{xs}_{9} \\ \mathtt{xs}_{10} \\ \mathtt{xs}_{11} \\ \mathtt{xs}_{12} \\ \mathtt{xs}_{13} \\ \mathtt{xs}_{14} \\ \mathtt{xs}_{15} \\ \mathtt{xs}_{16} \\ \mathtt{xs}_{17} \\ \mathtt{xs}_{18} \\ \mathtt{xs}_{19} \\ \mathtt{xs}_{20} \\ \end{array} \right] \end{equation} \end{split}\]
rxs = rosenbrock(xs)
\[ \begin{equation} \left( 1 - \mathtt{xs}_{1} \right)^{2} + \left( 1 - \mathtt{xs}_{10} \right)^{2} + \left( 1 - \mathtt{xs}_{11} \right)^{2} + \left( 1 - \mathtt{xs}_{12} \right)^{2} + \left( 1 - \mathtt{xs}_{13} \right)^{2} + \left( 1 - \mathtt{xs}_{14} \right)^{2} + \left( 1 - \mathtt{xs}_{15} \right)^{2} + \left( 1 - \mathtt{xs}_{16} \right)^{2} + \left( 1 - \mathtt{xs}_{17} \right)^{2} + \left( 1 - \mathtt{xs}_{18} \right)^{2} + \left( 1 - \mathtt{xs}_{19} \right)^{2} + \left( 1 - \mathtt{xs}_{2} \right)^{2} + \left( 1 - \mathtt{xs}_{3} \right)^{2} + \left( 1 - \mathtt{xs}_{4} \right)^{2} + \left( 1 - \mathtt{xs}_{5} \right)^{2} + \left( 1 - \mathtt{xs}_{6} \right)^{2} + \left( 1 - \mathtt{xs}_{7} \right)^{2} + \left( 1 - \mathtt{xs}_{8} \right)^{2} + \left( 1 - \mathtt{xs}_{9} \right)^{2} + 100 \left( \mathtt{xs}_{2} - \left( \mathtt{xs}_{1} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{11} - \left( \mathtt{xs}_{10} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{12} - \left( \mathtt{xs}_{11} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{13} - \left( \mathtt{xs}_{12} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{14} - \left( \mathtt{xs}_{13} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{15} - \left( \mathtt{xs}_{14} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{16} - \left( \mathtt{xs}_{15} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{17} - \left( \mathtt{xs}_{16} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{18} - \left( \mathtt{xs}_{17} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{19} - \left( \mathtt{xs}_{18} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{20} - \left( \mathtt{xs}_{19} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{3} - \left( \mathtt{xs}_{2} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{4} - \left( \mathtt{xs}_{3} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{5} - \left( \mathtt{xs}_{4} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{6} - \left( \mathtt{xs}_{5} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{7} - \left( \mathtt{xs}_{6} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{8} - \left( \mathtt{xs}_{7} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{9} - \left( \mathtt{xs}_{8} \right)^{2} \right)^{2} + 100 \left( \mathtt{xs}_{10} - \left( \mathtt{xs}_{9} \right)^{2} \right)^{2} \end{equation} \]

Gradient

grad = Symbolics.gradient(rxs, xs)
\[\begin{split} \begin{equation} \left[ \begin{array}{c} - 2 \left( 1 - \mathtt{xs}_{1} \right) - 400 \mathtt{xs}_{1} \left( \mathtt{xs}_{2} - \left( \mathtt{xs}_{1} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{2} \right) + 200 \left( \mathtt{xs}_{2} - \left( \mathtt{xs}_{1} \right)^{2} \right) - 400 \mathtt{xs}_{2} \left( \mathtt{xs}_{3} - \left( \mathtt{xs}_{2} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{3} \right) + 200 \left( \mathtt{xs}_{3} - \left( \mathtt{xs}_{2} \right)^{2} \right) - 400 \mathtt{xs}_{3} \left( \mathtt{xs}_{4} - \left( \mathtt{xs}_{3} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{4} \right) + 200 \left( \mathtt{xs}_{4} - \left( \mathtt{xs}_{3} \right)^{2} \right) - 400 \mathtt{xs}_{4} \left( \mathtt{xs}_{5} - \left( \mathtt{xs}_{4} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{5} \right) + 200 \left( \mathtt{xs}_{5} - \left( \mathtt{xs}_{4} \right)^{2} \right) - 400 \mathtt{xs}_{5} \left( \mathtt{xs}_{6} - \left( \mathtt{xs}_{5} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{6} \right) + 200 \left( \mathtt{xs}_{6} - \left( \mathtt{xs}_{5} \right)^{2} \right) - 400 \mathtt{xs}_{6} \left( \mathtt{xs}_{7} - \left( \mathtt{xs}_{6} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{7} \right) + 200 \left( \mathtt{xs}_{7} - \left( \mathtt{xs}_{6} \right)^{2} \right) - 400 \mathtt{xs}_{7} \left( \mathtt{xs}_{8} - \left( \mathtt{xs}_{7} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{8} \right) + 200 \left( \mathtt{xs}_{8} - \left( \mathtt{xs}_{7} \right)^{2} \right) - 400 \mathtt{xs}_{8} \left( \mathtt{xs}_{9} - \left( \mathtt{xs}_{8} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{9} \right) + 200 \left( \mathtt{xs}_{9} - \left( \mathtt{xs}_{8} \right)^{2} \right) - 400 \mathtt{xs}_{9} \left( \mathtt{xs}_{10} - \left( \mathtt{xs}_{9} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{10} \right) + 200 \left( \mathtt{xs}_{10} - \left( \mathtt{xs}_{9} \right)^{2} \right) - 400 \mathtt{xs}_{10} \left( \mathtt{xs}_{11} - \left( \mathtt{xs}_{10} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{11} \right) + 200 \left( \mathtt{xs}_{11} - \left( \mathtt{xs}_{10} \right)^{2} \right) - 400 \mathtt{xs}_{11} \left( \mathtt{xs}_{12} - \left( \mathtt{xs}_{11} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{12} \right) + 200 \left( \mathtt{xs}_{12} - \left( \mathtt{xs}_{11} \right)^{2} \right) - 400 \mathtt{xs}_{12} \left( \mathtt{xs}_{13} - \left( \mathtt{xs}_{12} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{13} \right) + 200 \left( \mathtt{xs}_{13} - \left( \mathtt{xs}_{12} \right)^{2} \right) - 400 \mathtt{xs}_{13} \left( \mathtt{xs}_{14} - \left( \mathtt{xs}_{13} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{14} \right) + 200 \left( \mathtt{xs}_{14} - \left( \mathtt{xs}_{13} \right)^{2} \right) - 400 \mathtt{xs}_{14} \left( \mathtt{xs}_{15} - \left( \mathtt{xs}_{14} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{15} \right) + 200 \left( \mathtt{xs}_{15} - \left( \mathtt{xs}_{14} \right)^{2} \right) - 400 \mathtt{xs}_{15} \left( \mathtt{xs}_{16} - \left( \mathtt{xs}_{15} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{16} \right) + 200 \left( \mathtt{xs}_{16} - \left( \mathtt{xs}_{15} \right)^{2} \right) - 400 \mathtt{xs}_{16} \left( \mathtt{xs}_{17} - \left( \mathtt{xs}_{16} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{17} \right) + 200 \left( \mathtt{xs}_{17} - \left( \mathtt{xs}_{16} \right)^{2} \right) - 400 \mathtt{xs}_{17} \left( \mathtt{xs}_{18} - \left( \mathtt{xs}_{17} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{18} \right) + 200 \left( \mathtt{xs}_{18} - \left( \mathtt{xs}_{17} \right)^{2} \right) - 400 \mathtt{xs}_{18} \left( \mathtt{xs}_{19} - \left( \mathtt{xs}_{18} \right)^{2} \right) \\ - 2 \left( 1 - \mathtt{xs}_{19} \right) + 200 \left( \mathtt{xs}_{19} - \left( \mathtt{xs}_{18} \right)^{2} \right) - 400 \mathtt{xs}_{19} \left( \mathtt{xs}_{20} - \left( \mathtt{xs}_{19} \right)^{2} \right) \\ 200 \left( \mathtt{xs}_{20} - \left( \mathtt{xs}_{19} \right)^{2} \right) \\ \end{array} \right] \end{equation} \end{split}\]

Hessian = Jacobian of gradient

hes1 = Symbolics.jacobian(grad, xs)
\[\begin{split} \begin{equation} \left[ \begin{array}{cccccccccccccccccccc} 2 + 800 \left( \mathtt{xs}_{1} \right)^{2} - 400 \left( \mathtt{xs}_{2} - \left( \mathtt{xs}_{1} \right)^{2} \right) & - 400 \mathtt{xs}_{1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - 400 \mathtt{xs}_{1} & 202 - 400 \left( \mathtt{xs}_{3} - \left( \mathtt{xs}_{2} \right)^{2} \right) + 800 \left( \mathtt{xs}_{2} \right)^{2} & - 400 \mathtt{xs}_{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & - 400 \mathtt{xs}_{2} & 202 - 400 \left( \mathtt{xs}_{4} - \left( \mathtt{xs}_{3} \right)^{2} \right) + 800 \left( \mathtt{xs}_{3} \right)^{2} & - 400 \mathtt{xs}_{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 400 \mathtt{xs}_{3} & 202 - 400 \left( \mathtt{xs}_{5} - \left( \mathtt{xs}_{4} \right)^{2} \right) + 800 \left( \mathtt{xs}_{4} \right)^{2} & - 400 \mathtt{xs}_{4} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - 400 \mathtt{xs}_{4} & 202 + 800 \left( \mathtt{xs}_{5} \right)^{2} - 400 \left( \mathtt{xs}_{6} - \left( \mathtt{xs}_{5} \right)^{2} \right) & - 400 \mathtt{xs}_{5} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{5} & 202 - 400 \left( \mathtt{xs}_{7} - \left( \mathtt{xs}_{6} \right)^{2} \right) + 800 \left( \mathtt{xs}_{6} \right)^{2} & - 400 \mathtt{xs}_{6} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{6} & 202 + 800 \left( \mathtt{xs}_{7} \right)^{2} - 400 \left( \mathtt{xs}_{8} - \left( \mathtt{xs}_{7} \right)^{2} \right) & - 400 \mathtt{xs}_{7} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{7} & 202 + 800 \left( \mathtt{xs}_{8} \right)^{2} - 400 \left( \mathtt{xs}_{9} - \left( \mathtt{xs}_{8} \right)^{2} \right) & - 400 \mathtt{xs}_{8} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{8} & 202 + 800 \left( \mathtt{xs}_{9} \right)^{2} - 400 \left( \mathtt{xs}_{10} - \left( \mathtt{xs}_{9} \right)^{2} \right) & - 400 \mathtt{xs}_{9} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{9} & 202 - 400 \left( \mathtt{xs}_{11} - \left( \mathtt{xs}_{10} \right)^{2} \right) + 800 \left( \mathtt{xs}_{10} \right)^{2} & - 400 \mathtt{xs}_{10} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{10} & 202 + 800 \left( \mathtt{xs}_{11} \right)^{2} - 400 \left( \mathtt{xs}_{12} - \left( \mathtt{xs}_{11} \right)^{2} \right) & - 400 \mathtt{xs}_{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{11} & 202 - 400 \left( \mathtt{xs}_{13} - \left( \mathtt{xs}_{12} \right)^{2} \right) + 800 \left( \mathtt{xs}_{12} \right)^{2} & - 400 \mathtt{xs}_{12} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{12} & 202 - 400 \left( \mathtt{xs}_{14} - \left( \mathtt{xs}_{13} \right)^{2} \right) + 800 \left( \mathtt{xs}_{13} \right)^{2} & - 400 \mathtt{xs}_{13} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{13} & 202 - 400 \left( \mathtt{xs}_{15} - \left( \mathtt{xs}_{14} \right)^{2} \right) + 800 \left( \mathtt{xs}_{14} \right)^{2} & - 400 \mathtt{xs}_{14} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{14} & 202 - 400 \left( \mathtt{xs}_{16} - \left( \mathtt{xs}_{15} \right)^{2} \right) + 800 \left( \mathtt{xs}_{15} \right)^{2} & - 400 \mathtt{xs}_{15} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{15} & 202 - 400 \left( \mathtt{xs}_{17} - \left( \mathtt{xs}_{16} \right)^{2} \right) + 800 \left( \mathtt{xs}_{16} \right)^{2} & - 400 \mathtt{xs}_{16} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{16} & 202 - 400 \left( \mathtt{xs}_{18} - \left( \mathtt{xs}_{17} \right)^{2} \right) + 800 \left( \mathtt{xs}_{17} \right)^{2} & - 400 \mathtt{xs}_{17} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{17} & 202 + 800 \left( \mathtt{xs}_{18} \right)^{2} - 400 \left( \mathtt{xs}_{19} - \left( \mathtt{xs}_{18} \right)^{2} \right) & - 400 \mathtt{xs}_{18} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{18} & 202 - 400 \left( \mathtt{xs}_{20} - \left( \mathtt{xs}_{19} \right)^{2} \right) + 800 \left( \mathtt{xs}_{19} \right)^{2} & - 400 \mathtt{xs}_{19} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{19} & 200 \\ \end{array} \right] \end{equation} \end{split}\]

call hessian() directly

hes2 = Symbolics.hessian(rxs, xs)
\[\begin{split} \begin{equation} \left[ \begin{array}{cccccccccccccccccccc} 2 + 800 \left( \mathtt{xs}_{1} \right)^{2} - 400 \left( \mathtt{xs}_{2} - \left( \mathtt{xs}_{1} \right)^{2} \right) & - 400 \mathtt{xs}_{1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - 400 \mathtt{xs}_{1} & 202 - 400 \left( \mathtt{xs}_{3} - \left( \mathtt{xs}_{2} \right)^{2} \right) + 800 \left( \mathtt{xs}_{2} \right)^{2} & - 400 \mathtt{xs}_{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & - 400 \mathtt{xs}_{2} & 202 - 400 \left( \mathtt{xs}_{4} - \left( \mathtt{xs}_{3} \right)^{2} \right) + 800 \left( \mathtt{xs}_{3} \right)^{2} & - 400 \mathtt{xs}_{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 400 \mathtt{xs}_{3} & 202 - 400 \left( \mathtt{xs}_{5} - \left( \mathtt{xs}_{4} \right)^{2} \right) + 800 \left( \mathtt{xs}_{4} \right)^{2} & - 400 \mathtt{xs}_{4} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - 400 \mathtt{xs}_{4} & 202 + 800 \left( \mathtt{xs}_{5} \right)^{2} - 400 \left( \mathtt{xs}_{6} - \left( \mathtt{xs}_{5} \right)^{2} \right) & - 400 \mathtt{xs}_{5} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{5} & 202 - 400 \left( \mathtt{xs}_{7} - \left( \mathtt{xs}_{6} \right)^{2} \right) + 800 \left( \mathtt{xs}_{6} \right)^{2} & - 400 \mathtt{xs}_{6} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{6} & 202 + 800 \left( \mathtt{xs}_{7} \right)^{2} - 400 \left( \mathtt{xs}_{8} - \left( \mathtt{xs}_{7} \right)^{2} \right) & - 400 \mathtt{xs}_{7} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{7} & 202 + 800 \left( \mathtt{xs}_{8} \right)^{2} - 400 \left( \mathtt{xs}_{9} - \left( \mathtt{xs}_{8} \right)^{2} \right) & - 400 \mathtt{xs}_{8} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{8} & 202 + 800 \left( \mathtt{xs}_{9} \right)^{2} - 400 \left( \mathtt{xs}_{10} - \left( \mathtt{xs}_{9} \right)^{2} \right) & - 400 \mathtt{xs}_{9} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{9} & 202 - 400 \left( \mathtt{xs}_{11} - \left( \mathtt{xs}_{10} \right)^{2} \right) + 800 \left( \mathtt{xs}_{10} \right)^{2} & - 400 \mathtt{xs}_{10} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{10} & 202 + 800 \left( \mathtt{xs}_{11} \right)^{2} - 400 \left( \mathtt{xs}_{12} - \left( \mathtt{xs}_{11} \right)^{2} \right) & - 400 \mathtt{xs}_{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{11} & 202 - 400 \left( \mathtt{xs}_{13} - \left( \mathtt{xs}_{12} \right)^{2} \right) + 800 \left( \mathtt{xs}_{12} \right)^{2} & - 400 \mathtt{xs}_{12} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{12} & 202 - 400 \left( \mathtt{xs}_{14} - \left( \mathtt{xs}_{13} \right)^{2} \right) + 800 \left( \mathtt{xs}_{13} \right)^{2} & - 400 \mathtt{xs}_{13} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{13} & 202 - 400 \left( \mathtt{xs}_{15} - \left( \mathtt{xs}_{14} \right)^{2} \right) + 800 \left( \mathtt{xs}_{14} \right)^{2} & - 400 \mathtt{xs}_{14} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{14} & 202 - 400 \left( \mathtt{xs}_{16} - \left( \mathtt{xs}_{15} \right)^{2} \right) + 800 \left( \mathtt{xs}_{15} \right)^{2} & - 400 \mathtt{xs}_{15} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{15} & 202 - 400 \left( \mathtt{xs}_{17} - \left( \mathtt{xs}_{16} \right)^{2} \right) + 800 \left( \mathtt{xs}_{16} \right)^{2} & - 400 \mathtt{xs}_{16} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{16} & 202 - 400 \left( \mathtt{xs}_{18} - \left( \mathtt{xs}_{17} \right)^{2} \right) + 800 \left( \mathtt{xs}_{17} \right)^{2} & - 400 \mathtt{xs}_{17} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{17} & 202 + 800 \left( \mathtt{xs}_{18} \right)^{2} - 400 \left( \mathtt{xs}_{19} - \left( \mathtt{xs}_{18} \right)^{2} \right) & - 400 \mathtt{xs}_{18} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{18} & 202 - 400 \left( \mathtt{xs}_{20} - \left( \mathtt{xs}_{19} \right)^{2} \right) + 800 \left( \mathtt{xs}_{19} \right)^{2} & - 400 \mathtt{xs}_{19} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 \mathtt{xs}_{19} & 200 \\ \end{array} \right] \end{equation} \end{split}\]
isequal(hes1, hes2)

true

Sparse matrix#

Sparse Hessian matrix of the Hessian matrix of the Rosenbrock function w.r.t. to vector components.

hes_sp = Symbolics.hessian_sparsity(rosenbrock, xs)

20×20 SparseArrays.SparseMatrixCSC{Bool, Int64} with 58 stored entries:
⎡⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦

Visualize the sparse matrix with Plots.spy()

using Plots
spy(hes_sp)
_images/f52d4518e711d429b941280d703d4f75ada1866b2b0d90badbeb7edf6d99b41f.png

Generate functions from symbols#

  • build_function(ex, args...) generates out-of-place (oop) and in-place (ip) function expressions in a pair.

  • build_function(ex, args..., parallel=Symbolics.MultithreadedForm()) generates a parallel algorithm to evaluate the output. See the example in the official docs.

  • build_function(ex, args..., target=Symbolics.CTarget()) generates a C function from Julia.

fexprs = build_function(grad, xs);

Get the Out-of-place f(input) version

foop = eval(fexprs[1])

#3 (generic function with 1 method)

Get the In-place f!(out, in) version

fip = eval(fexprs[2])

#5 (generic function with 1 method)

inxs = rand(N)
out = similar(inxs)
fip(out, inxs)  ## The inplace version returns nothing. The results are stored in out parameter.

foop(inxs)
isapprox(foop(inxs), out)

true

You can save the generated function for later use: write("function.jl", string(fexprs[2])) Here, ForwardDiff.jl checks if our gradient generated from Symbolics.jl is correct.

using ForwardDiff: gradient
gradient(rosenbrock, inxs) ≈ out

true

Sparse Hessian matrix, only non-zero expressions are calculated.

hexprs = build_function(hes_sp, xs)
hoop = eval(hexprs[1])
hip = eval(hexprs[2])
hoop(rand(N))

20×20 SparseArrays.SparseMatrixCSC{Bool, Int64} with 58 stored entries:
⎡⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦

This notebook was generated using Literate.jl.