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.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} - \mathtt{xs}_{1}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{11} - \mathtt{xs}_{10}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{12} - \mathtt{xs}_{11}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{13} - \mathtt{xs}_{12}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{14} - \mathtt{xs}_{13}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{15} - \mathtt{xs}_{14}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{16} - \mathtt{xs}_{15}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{17} - \mathtt{xs}_{16}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{18} - \mathtt{xs}_{17}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{19} - \mathtt{xs}_{18}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{20} - \mathtt{xs}_{19}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{3} - \mathtt{xs}_{2}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{4} - \mathtt{xs}_{3}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{5} - \mathtt{xs}_{4}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{6} - \mathtt{xs}_{5}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{7} - \mathtt{xs}_{6}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{8} - \mathtt{xs}_{7}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{9} - \mathtt{xs}_{8}^{2} \right)^{2} + 100 \left( \mathtt{xs}_{10} - \mathtt{xs}_{9}^{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} - \mathtt{xs}_{1}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{2} \right) + 200 \left( \mathtt{xs}_{2} - \mathtt{xs}_{1}^{2} \right) - 400 \mathtt{xs}_{2} \left( \mathtt{xs}_{3} - \mathtt{xs}_{2}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{3} \right) + 200 \left( \mathtt{xs}_{3} - \mathtt{xs}_{2}^{2} \right) - 400 \mathtt{xs}_{3} \left( \mathtt{xs}_{4} - \mathtt{xs}_{3}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{4} \right) + 200 \left( \mathtt{xs}_{4} - \mathtt{xs}_{3}^{2} \right) - 400 \mathtt{xs}_{4} \left( \mathtt{xs}_{5} - \mathtt{xs}_{4}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{5} \right) + 200 \left( \mathtt{xs}_{5} - \mathtt{xs}_{4}^{2} \right) - 400 \mathtt{xs}_{5} \left( \mathtt{xs}_{6} - \mathtt{xs}_{5}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{6} \right) + 200 \left( \mathtt{xs}_{6} - \mathtt{xs}_{5}^{2} \right) - 400 \mathtt{xs}_{6} \left( \mathtt{xs}_{7} - \mathtt{xs}_{6}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{7} \right) + 200 \left( \mathtt{xs}_{7} - \mathtt{xs}_{6}^{2} \right) - 400 \mathtt{xs}_{7} \left( \mathtt{xs}_{8} - \mathtt{xs}_{7}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{8} \right) + 200 \left( \mathtt{xs}_{8} - \mathtt{xs}_{7}^{2} \right) - 400 \mathtt{xs}_{8} \left( \mathtt{xs}_{9} - \mathtt{xs}_{8}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{9} \right) + 200 \left( \mathtt{xs}_{9} - \mathtt{xs}_{8}^{2} \right) - 400 \mathtt{xs}_{9} \left( \mathtt{xs}_{10} - \mathtt{xs}_{9}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{10} \right) + 200 \left( \mathtt{xs}_{10} - \mathtt{xs}_{9}^{2} \right) - 400 \mathtt{xs}_{10} \left( \mathtt{xs}_{11} - \mathtt{xs}_{10}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{11} \right) + 200 \left( \mathtt{xs}_{11} - \mathtt{xs}_{10}^{2} \right) - 400 \mathtt{xs}_{11} \left( \mathtt{xs}_{12} - \mathtt{xs}_{11}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{12} \right) + 200 \left( \mathtt{xs}_{12} - \mathtt{xs}_{11}^{2} \right) - 400 \mathtt{xs}_{12} \left( \mathtt{xs}_{13} - \mathtt{xs}_{12}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{13} \right) + 200 \left( \mathtt{xs}_{13} - \mathtt{xs}_{12}^{2} \right) - 400 \mathtt{xs}_{13} \left( \mathtt{xs}_{14} - \mathtt{xs}_{13}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{14} \right) + 200 \left( \mathtt{xs}_{14} - \mathtt{xs}_{13}^{2} \right) - 400 \mathtt{xs}_{14} \left( \mathtt{xs}_{15} - \mathtt{xs}_{14}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{15} \right) + 200 \left( \mathtt{xs}_{15} - \mathtt{xs}_{14}^{2} \right) - 400 \mathtt{xs}_{15} \left( \mathtt{xs}_{16} - \mathtt{xs}_{15}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{16} \right) + 200 \left( \mathtt{xs}_{16} - \mathtt{xs}_{15}^{2} \right) - 400 \mathtt{xs}_{16} \left( \mathtt{xs}_{17} - \mathtt{xs}_{16}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{17} \right) + 200 \left( \mathtt{xs}_{17} - \mathtt{xs}_{16}^{2} \right) - 400 \mathtt{xs}_{17} \left( \mathtt{xs}_{18} - \mathtt{xs}_{17}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{18} \right) + 200 \left( \mathtt{xs}_{18} - \mathtt{xs}_{17}^{2} \right) - 400 \mathtt{xs}_{18} \left( \mathtt{xs}_{19} - \mathtt{xs}_{18}^{2} \right) \\
- 2 \left( 1 - \mathtt{xs}_{19} \right) + 200 \left( \mathtt{xs}_{19} - \mathtt{xs}_{18}^{2} \right) - 400 \mathtt{xs}_{19} \left( \mathtt{xs}_{20} - \mathtt{xs}_{19}^{2} \right) \\
200 \left( \mathtt{xs}_{20} - \mathtt{xs}_{19}^{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 \mathtt{xs}_{1}^{2} - 400 \left( \mathtt{xs}_{2} - \mathtt{xs}_{1}^{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} - \mathtt{xs}_{2}^{2} \right) + 800 \mathtt{xs}_{2}^{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 + 800 \mathtt{xs}_{3}^{2} - 400 \left( \mathtt{xs}_{4} - \mathtt{xs}_{3}^{2} \right) & - 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} - \mathtt{xs}_{4}^{2} \right) + 800 \mathtt{xs}_{4}^{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 \mathtt{xs}_{5}^{2} - 400 \left( \mathtt{xs}_{6} - \mathtt{xs}_{5}^{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} - \mathtt{xs}_{6}^{2} \right) + 800 \mathtt{xs}_{6}^{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 \mathtt{xs}_{7}^{2} - 400 \left( \mathtt{xs}_{8} - \mathtt{xs}_{7}^{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 \mathtt{xs}_{8}^{2} - 400 \left( \mathtt{xs}_{9} - \mathtt{xs}_{8}^{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 \mathtt{xs}_{9}^{2} - 400 \left( \mathtt{xs}_{10} - \mathtt{xs}_{9}^{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} - \mathtt{xs}_{10}^{2} \right) + 800 \mathtt{xs}_{10}^{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 \mathtt{xs}_{11}^{2} - 400 \left( \mathtt{xs}_{12} - \mathtt{xs}_{11}^{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} - \mathtt{xs}_{12}^{2} \right) + 800 \mathtt{xs}_{12}^{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} - \mathtt{xs}_{13}^{2} \right) + 800 \mathtt{xs}_{13}^{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} - \mathtt{xs}_{14}^{2} \right) + 800 \mathtt{xs}_{14}^{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} - \mathtt{xs}_{15}^{2} \right) + 800 \mathtt{xs}_{15}^{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} - \mathtt{xs}_{16}^{2} \right) + 800 \mathtt{xs}_{16}^{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} - \mathtt{xs}_{17}^{2} \right) + 800 \mathtt{xs}_{17}^{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 \mathtt{xs}_{18}^{2} - 400 \left( \mathtt{xs}_{19} - \mathtt{xs}_{18}^{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} - \mathtt{xs}_{19}^{2} \right) + 800 \mathtt{xs}_{19}^{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 \mathtt{xs}_{1}^{2} - 400 \left( \mathtt{xs}_{2} - \mathtt{xs}_{1}^{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} - \mathtt{xs}_{2}^{2} \right) + 800 \mathtt{xs}_{2}^{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 + 800 \mathtt{xs}_{3}^{2} - 400 \left( \mathtt{xs}_{4} - \mathtt{xs}_{3}^{2} \right) & - 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} - \mathtt{xs}_{4}^{2} \right) + 800 \mathtt{xs}_{4}^{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 \mathtt{xs}_{5}^{2} - 400 \left( \mathtt{xs}_{6} - \mathtt{xs}_{5}^{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} - \mathtt{xs}_{6}^{2} \right) + 800 \mathtt{xs}_{6}^{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 \mathtt{xs}_{7}^{2} - 400 \left( \mathtt{xs}_{8} - \mathtt{xs}_{7}^{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 \mathtt{xs}_{8}^{2} - 400 \left( \mathtt{xs}_{9} - \mathtt{xs}_{8}^{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 \mathtt{xs}_{9}^{2} - 400 \left( \mathtt{xs}_{10} - \mathtt{xs}_{9}^{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} - \mathtt{xs}_{10}^{2} \right) + 800 \mathtt{xs}_{10}^{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 \mathtt{xs}_{11}^{2} - 400 \left( \mathtt{xs}_{12} - \mathtt{xs}_{11}^{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} - \mathtt{xs}_{12}^{2} \right) + 800 \mathtt{xs}_{12}^{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} - \mathtt{xs}_{13}^{2} \right) + 800 \mathtt{xs}_{13}^{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} - \mathtt{xs}_{14}^{2} \right) + 800 \mathtt{xs}_{14}^{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} - \mathtt{xs}_{15}^{2} \right) + 800 \mathtt{xs}_{15}^{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} - \mathtt{xs}_{16}^{2} \right) + 800 \mathtt{xs}_{16}^{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} - \mathtt{xs}_{17}^{2} \right) + 800 \mathtt{xs}_{17}^{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 \mathtt{xs}_{18}^{2} - 400 \left( \mathtt{xs}_{19} - \mathtt{xs}_{18}^{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} - \mathtt{xs}_{19}^{2} \right) + 800 \mathtt{xs}_{19}^{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)
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
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.