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:
Simulating Big Models in Julia with ModelingToolkit @ JuliaCon 2021 Workshop.
Symbolics.jl Github repo and its docs.
Caveats about Symbolics.jl#
Symbolics.jlcan 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
@registerto make it a primitive function.
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#
latexifyderivativegradientjacobiansubstitutesimplify
using Symbolics
using Latexify
@variables x y
x^2 + y^2
\[ \begin{equation}
x^{2} + y^{2}
\end{equation}
\]
You can use 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 expresseion 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}
\]
Symblic 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 evaulated 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{Num, 1}}:
xs[1:18]
xs[1]
\[ \begin{equation}
xs_1
\end{equation}
\]
Explicit vector form
collect(xs)
\[\begin{split} \begin{equation}
\left[
\begin{array}{c}
xs_1 \\
xs_2 \\
xs_3 \\
xs_4 \\
xs_5 \\
xs_6 \\
xs_7 \\
xs_8 \\
xs_9 \\
xs_{1 0} \\
xs_{1 1} \\
xs_{1 2} \\
xs_{1 3} \\
xs_{1 4} \\
xs_{1 5} \\
xs_{1 6} \\
xs_{1 7} \\
xs_{1 8} \\
\end{array}
\right]
\end{equation}
\end{split}\]
Operations on arrays are supported
sum(collect(xs))
\[ \begin{equation}
xs_1 + xs_{1 0} + xs_{1 1} + xs_{1 2} + xs_{1 3} + xs_{1 4} + xs_{1 5} + xs_{1 6} + xs_{1 7} + xs_{1 8} + xs_2 + xs_3 + xs_4 + xs_5 + xs_6 + xs_7 + xs_8 + 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{Num, 1}}:
xs[1:20]
A full list of vector components
xs = collect(xs)
\[\begin{split} \begin{equation}
\left[
\begin{array}{c}
xs_1 \\
xs_2 \\
xs_3 \\
xs_4 \\
xs_5 \\
xs_6 \\
xs_7 \\
xs_8 \\
xs_9 \\
xs_{1 0} \\
xs_{1 1} \\
xs_{1 2} \\
xs_{1 3} \\
xs_{1 4} \\
xs_{1 5} \\
xs_{1 6} \\
xs_{1 7} \\
xs_{1 8} \\
xs_{1 9} \\
xs_{2 0} \\
\end{array}
\right]
\end{equation}
\end{split}\]
rxs = rosenbrock(xs)
\[ \begin{equation}
\left( 1 - xs_1 \right)^{2} + \left( 1 - xs_{1 0} \right)^{2} + \left( 1 - xs_{1 1} \right)^{2} + \left( 1 - xs_{1 2} \right)^{2} + \left( 1 - xs_{1 3} \right)^{2} + \left( 1 - xs_{1 4} \right)^{2} + \left( 1 - xs_{1 5} \right)^{2} + \left( 1 - xs_{1 6} \right)^{2} + \left( 1 - xs_{1 7} \right)^{2} + \left( 1 - xs_{1 8} \right)^{2} + \left( 1 - xs_{1 9} \right)^{2} + \left( 1 - xs_2 \right)^{2} + \left( 1 - xs_3 \right)^{2} + \left( 1 - xs_4 \right)^{2} + \left( 1 - xs_5 \right)^{2} + \left( 1 - xs_6 \right)^{2} + \left( 1 - xs_7 \right)^{2} + \left( 1 - xs_8 \right)^{2} + \left( 1 - xs_9 \right)^{2} + 100 \left( xs_2 - xs_1^{2} \right)^{2} + 100 \left( xs_{1 1} - xs_{1 0}^{2} \right)^{2} + 100 \left( xs_{1 2} - xs_{1 1}^{2} \right)^{2} + 100 \left( xs_{1 3} - xs_{1 2}^{2} \right)^{2} + 100 \left( xs_{1 4} - xs_{1 3}^{2} \right)^{2} + 100 \left( xs_{1 5} - xs_{1 4}^{2} \right)^{2} + 100 \left( xs_{1 6} - xs_{1 5}^{2} \right)^{2} + 100 \left( xs_{1 7} - xs_{1 6}^{2} \right)^{2} + 100 \left( xs_{1 8} - xs_{1 7}^{2} \right)^{2} + 100 \left( xs_{1 9} - xs_{1 8}^{2} \right)^{2} + 100 \left( xs_{2 0} - xs_{1 9}^{2} \right)^{2} + 100 \left( xs_3 - xs_2^{2} \right)^{2} + 100 \left( xs_4 - xs_3^{2} \right)^{2} + 100 \left( xs_5 - xs_4^{2} \right)^{2} + 100 \left( xs_6 - xs_5^{2} \right)^{2} + 100 \left( xs_7 - xs_6^{2} \right)^{2} + 100 \left( xs_8 - xs_7^{2} \right)^{2} + 100 \left( xs_9 - xs_8^{2} \right)^{2} + 100 \left( xs_{1 0} - xs_9^{2} \right)^{2}
\end{equation}
\]
Gradient
grad = Symbolics.gradient(rxs, xs)
\[\begin{split} \begin{equation}
\left[
\begin{array}{c}
- 2 \left( 1 - xs_1 \right) - 400 xs_1 \left( xs_2 - xs_1^{2} \right) \\
- 2 \left( 1 - xs_2 \right) + 200 \left( xs_2 - xs_1^{2} \right) - 400 xs_2 \left( xs_3 - xs_2^{2} \right) \\
- 2 \left( 1 - xs_3 \right) + 200 \left( xs_3 - xs_2^{2} \right) - 400 xs_3 \left( xs_4 - xs_3^{2} \right) \\
- 2 \left( 1 - xs_4 \right) + 200 \left( xs_4 - xs_3^{2} \right) - 400 xs_4 \left( xs_5 - xs_4^{2} \right) \\
- 2 \left( 1 - xs_5 \right) + 200 \left( xs_5 - xs_4^{2} \right) - 400 xs_5 \left( xs_6 - xs_5^{2} \right) \\
- 2 \left( 1 - xs_6 \right) + 200 \left( xs_6 - xs_5^{2} \right) - 400 xs_6 \left( xs_7 - xs_6^{2} \right) \\
- 2 \left( 1 - xs_7 \right) + 200 \left( xs_7 - xs_6^{2} \right) - 400 xs_7 \left( xs_8 - xs_7^{2} \right) \\
- 2 \left( 1 - xs_8 \right) + 200 \left( xs_8 - xs_7^{2} \right) - 400 xs_8 \left( xs_9 - xs_8^{2} \right) \\
- 2 \left( 1 - xs_9 \right) + 200 \left( xs_9 - xs_8^{2} \right) - 400 xs_9 \left( xs_{1 0} - xs_9^{2} \right) \\
- 2 \left( 1 - xs_{1 0} \right) + 200 \left( xs_{1 0} - xs_9^{2} \right) - 400 xs_{1 0} \left( xs_{1 1} - xs_{1 0}^{2} \right) \\
- 2 \left( 1 - xs_{1 1} \right) + 200 \left( xs_{1 1} - xs_{1 0}^{2} \right) - 400 xs_{1 1} \left( xs_{1 2} - xs_{1 1}^{2} \right) \\
- 2 \left( 1 - xs_{1 2} \right) + 200 \left( xs_{1 2} - xs_{1 1}^{2} \right) - 400 xs_{1 2} \left( xs_{1 3} - xs_{1 2}^{2} \right) \\
- 2 \left( 1 - xs_{1 3} \right) + 200 \left( xs_{1 3} - xs_{1 2}^{2} \right) - 400 xs_{1 3} \left( xs_{1 4} - xs_{1 3}^{2} \right) \\
- 2 \left( 1 - xs_{1 4} \right) + 200 \left( xs_{1 4} - xs_{1 3}^{2} \right) - 400 xs_{1 4} \left( xs_{1 5} - xs_{1 4}^{2} \right) \\
- 2 \left( 1 - xs_{1 5} \right) + 200 \left( xs_{1 5} - xs_{1 4}^{2} \right) - 400 xs_{1 5} \left( xs_{1 6} - xs_{1 5}^{2} \right) \\
- 2 \left( 1 - xs_{1 6} \right) + 200 \left( xs_{1 6} - xs_{1 5}^{2} \right) - 400 xs_{1 6} \left( xs_{1 7} - xs_{1 6}^{2} \right) \\
- 2 \left( 1 - xs_{1 7} \right) + 200 \left( xs_{1 7} - xs_{1 6}^{2} \right) - 400 xs_{1 7} \left( xs_{1 8} - xs_{1 7}^{2} \right) \\
- 2 \left( 1 - xs_{1 8} \right) + 200 \left( xs_{1 8} - xs_{1 7}^{2} \right) - 400 xs_{1 8} \left( xs_{1 9} - xs_{1 8}^{2} \right) \\
- 2 \left( 1 - xs_{1 9} \right) + 200 \left( xs_{1 9} - xs_{1 8}^{2} \right) - 400 xs_{1 9} \left( xs_{2 0} - xs_{1 9}^{2} \right) \\
200 \left( xs_{2 0} - xs_{1 9}^{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 xs_1^{2} - 400 \left( xs_2 - xs_1^{2} \right) & - 400 xs_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
- 400 xs_1 & 202 - 400 \left( xs_3 - xs_2^{2} \right) + 800 xs_2^{2} & - 400 xs_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & - 400 xs_2 & 202 - 400 \left( xs_4 - xs_3^{2} \right) + 800 xs_3^{2} & - 400 xs_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & - 400 xs_3 & 202 - 400 \left( xs_5 - xs_4^{2} \right) + 800 xs_4^{2} & - 400 xs_4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & - 400 xs_4 & 202 + 800 xs_5^{2} - 400 \left( xs_6 - xs_5^{2} \right) & - 400 xs_5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & - 400 xs_5 & 202 - 400 \left( xs_7 - xs_6^{2} \right) + 800 xs_6^{2} & - 400 xs_6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & - 400 xs_6 & 202 + 800 xs_7^{2} - 400 \left( xs_8 - xs_7^{2} \right) & - 400 xs_7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_7 & 202 + 800 xs_8^{2} - 400 \left( xs_9 - xs_8^{2} \right) & - 400 xs_8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_8 & 202 + 800 xs_9^{2} - 400 \left( xs_{1 0} - xs_9^{2} \right) & - 400 xs_9 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_9 & 202 - 400 \left( xs_{1 1} - xs_{1 0}^{2} \right) + 800 xs_{1 0}^{2} & - 400 xs_{1 0} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 0} & 202 + 800 xs_{1 1}^{2} - 400 \left( xs_{1 2} - xs_{1 1}^{2} \right) & - 400 xs_{1 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 1} & 202 - 400 \left( xs_{1 3} - xs_{1 2}^{2} \right) + 800 xs_{1 2}^{2} & - 400 xs_{1 2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 2} & 202 - 400 \left( xs_{1 4} - xs_{1 3}^{2} \right) + 800 xs_{1 3}^{2} & - 400 xs_{1 3} & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 3} & 202 - 400 \left( xs_{1 5} - xs_{1 4}^{2} \right) + 800 xs_{1 4}^{2} & - 400 xs_{1 4} & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 4} & 202 - 400 \left( xs_{1 6} - xs_{1 5}^{2} \right) + 800 xs_{1 5}^{2} & - 400 xs_{1 5} & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 5} & 202 - 400 \left( xs_{1 7} - xs_{1 6}^{2} \right) + 800 xs_{1 6}^{2} & - 400 xs_{1 6} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 6} & 202 - 400 \left( xs_{1 8} - xs_{1 7}^{2} \right) + 800 xs_{1 7}^{2} & - 400 xs_{1 7} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 7} & 202 + 800 xs_{1 8}^{2} - 400 \left( xs_{1 9} - xs_{1 8}^{2} \right) & - 400 xs_{1 8} & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 8} & 202 - 400 \left( xs_{2 0} - xs_{1 9}^{2} \right) + 800 xs_{1 9}^{2} & - 400 xs_{1 9} \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 9} & 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 xs_1^{2} - 400 \left( xs_2 - xs_1^{2} \right) & - 400 xs_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
- 400 xs_1 & 202 - 400 \left( xs_3 - xs_2^{2} \right) + 800 xs_2^{2} & - 400 xs_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & - 400 xs_2 & 202 - 400 \left( xs_4 - xs_3^{2} \right) + 800 xs_3^{2} & - 400 xs_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & - 400 xs_3 & 202 - 400 \left( xs_5 - xs_4^{2} \right) + 800 xs_4^{2} & - 400 xs_4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & - 400 xs_4 & 202 + 800 xs_5^{2} - 400 \left( xs_6 - xs_5^{2} \right) & - 400 xs_5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & - 400 xs_5 & 202 - 400 \left( xs_7 - xs_6^{2} \right) + 800 xs_6^{2} & - 400 xs_6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & - 400 xs_6 & 202 + 800 xs_7^{2} - 400 \left( xs_8 - xs_7^{2} \right) & - 400 xs_7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_7 & 202 + 800 xs_8^{2} - 400 \left( xs_9 - xs_8^{2} \right) & - 400 xs_8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_8 & 202 + 800 xs_9^{2} - 400 \left( xs_{1 0} - xs_9^{2} \right) & - 400 xs_9 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_9 & 202 - 400 \left( xs_{1 1} - xs_{1 0}^{2} \right) + 800 xs_{1 0}^{2} & - 400 xs_{1 0} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 0} & 202 + 800 xs_{1 1}^{2} - 400 \left( xs_{1 2} - xs_{1 1}^{2} \right) & - 400 xs_{1 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 1} & 202 - 400 \left( xs_{1 3} - xs_{1 2}^{2} \right) + 800 xs_{1 2}^{2} & - 400 xs_{1 2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 2} & 202 - 400 \left( xs_{1 4} - xs_{1 3}^{2} \right) + 800 xs_{1 3}^{2} & - 400 xs_{1 3} & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 3} & 202 - 400 \left( xs_{1 5} - xs_{1 4}^{2} \right) + 800 xs_{1 4}^{2} & - 400 xs_{1 4} & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 4} & 202 - 400 \left( xs_{1 6} - xs_{1 5}^{2} \right) + 800 xs_{1 5}^{2} & - 400 xs_{1 5} & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 5} & 202 - 400 \left( xs_{1 7} - xs_{1 6}^{2} \right) + 800 xs_{1 6}^{2} & - 400 xs_{1 6} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 6} & 202 - 400 \left( xs_{1 8} - xs_{1 7}^{2} \right) + 800 xs_{1 7}^{2} & - 400 xs_{1 7} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 7} & 202 + 800 xs_{1 8}^{2} - 400 \left( xs_{1 9} - xs_{1 8}^{2} \right) & - 400 xs_{1 8} & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 8} & 202 - 400 \left( xs_{2 0} - xs_{1 9}^{2} \right) + 800 xs_{1 9}^{2} & - 400 xs_{1 9} \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 400 xs_{1 9} & 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)
[ Info: Precompiling IJuliaExt [2f4121a4-3b3a-5ce6-9c5e-1f2673ce168a]
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
ForwardDiff.gradient(rosenbrock, inxs) ≈ out
true
Sparce 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.