In this SciML update we have plenty of new features to check out including new symbolic solvers, direct Enzyme support on OrdinaryDiffEq.jl, and major loading improvements across the ecosystem. Time for the fun details!
A long-time requested feature for Symbolics.jl has been to add symbolic solvers. This is functionality to say give a system like x^2 - y = 5, sin(y) + cos(x) = 2, and let it spit out a symbolic exact solution. Well with this release it's finally here! Thanks to a Google Summer of Code project by Yassin ElBedwihy, the core pull-request adding the symbolic solver has finally landed. With this PR, you can do things like:
using Symbolics, Groebner
@variables x y z;
eqs = [x^2 + y + z - 1, x + y^2 + z - 1, x + y + z^2 - 1]
Symbolics.symbolic_solve(eqs, [x,y,z])and get the symbolic solution to the set of polynomial equations. Want to know the answer to the above? Well install it and find it!
It comes complete with extensions that make use of the Nemo Computer Algebra System (CAS), which is an abstract algebra system developed by other smart Julia developers that specializes in "abstract algebra" like computational group theory, computational ring theory, and more. It mixes some of these techniques in with rule-based techniques to give a solver that is "best of both worlds" and is easily extendable with other submodules. For Groebner basis calculations, it uses the Groebner.jl package which is demonstrated to be one of the most efficient implementations and thus serves as a solid building block.
In more detail, The symbolic_solve function uses 4 hidden solvers in order to solve the user's input. Its base, solve_univar, uses analytic solutions up to polynomials of degree 4 and factoring as its method for solving univariate polynomials. The function's solve_multipoly uses GCD on the input polynomials then throws passes the result to solve_univar. The function's solve_multivar uses Groebner basis and a separating form in order to create linear equations in the input variables and a single high degree equation in the separating variable. Each equation resulting from the basis is then passed to solve_univar. We can see that essentially, solve_univar is the building block of symbolic_solve. If the input is not a valid polynomial and can not be solved by the algorithm above, symbolic_solve passes it to ia_solve, which attempts solving by attraction and isolation. This only works when the input is a single expression and the user wants the answer in terms of a single variable. Say log(x) - a == 0 gives us [e^a]. This attraction isolation is then extendable via rules, so down the line we can add extensions to handle cases like LambertW.jl and SpecialFunctions.jl detection.
Its current feature completeness can be summarized as:
[x] Linear and polynomial equations
[x] Systems of linear and polynomial equations
[x] Some transcendental functions
[x] Systems of linear equations with parameters (via symbolic_linear_solve)
[ ] Systems of polynomial equations with parameters
[ ] Inequalities
[ ] Differential Equations (ODEs)
[ ] Integrals
With plans to continue developing and handle the next cases soon after.
A major was released on the Symbolics stack this month, signifying a breaking change. The major breaking change here is the re-adoption of TermInterface.jl, which is a core interface for symbolic terms. By having all symbolic libraries extend term interface, be it Metatheory.jl, Symbolics.jl, SymPy.jl, and more, the core interface gives a common specification for building and translating terms, making all of them interopable. This means that by re-adopting TermInterface, we now have bidirectional translation to and from SymPy as a well-maintained part of the interfaces, and moving between rule-based approaches and E-graphs based approaches is a standard part of the interfaces.
At the same time that we tacked TermInterface v2, we also tackled some of the long-standing problems in the Symbolics ecosystem. In particular, for a symbolic term like ex = f(x,y), while arguments(ex) == [x,y], if you build a symbolic expression of ex = x + y + z the internal data structures can be optimized by assuming no ordering in such a commutative operation. However before we guarenteed argument order on arguments(ex), so then arguments(ex) = [x,y,z] was enforced to always be sorted lexicographically. However, it turns out that after building out the symbolic stack that this choice is one of the most costly in the entire ecosystem! Thus we have changed arguments(ex) to not guarantee a sorting order, allowing symbolic manipulations which specialize on commutativity to skip spending 99% of their time calculating lexicographic sorts. In order to allow for printing in a stable manner, we created the new interface functions such as sorted_arguments(ex) which guarantee a sorting and thus take the performance hit, and this is used so that things like Latex outputs and displays are more stable.
Another major change was a breaking change to the maketerm syntax to remove the symtype argument. While symbolic terms like those in Symbolics.jl can still be typed, i.e. @variables x::Complex which changes their behavior in things like simplification rules, this information is now captured in the metadata instead of the term itself. This unifies more of the implementation between Symbolics.jl and Metatheory.jl to better allow usage of E-graphs on symbolic terms.
SymbolicUtils.jl v3 and Symboilcs v6 thus take these interface changes as their breaking bits. Most code should actually not be broken by these changes, we only had to update the code of approximately 10% of the upstream libraries to allow for this change, and many of those only because we the symbolics developers are using some deeper features. So it's generally a small break but we get some major performance improvements and nice new features that unify the ecosystem.
We had a few major Symbolics.jl updates, with this year having Symbolics v4, v5, and now v6, a few with SymbolicUtils, and new a few with TermInterface. We are happy to report that we believe TermInterface.jl is now finally stable, so these updates should have reached their conclusion. There is a common core change to SymbolicUtils.jl which should be a major performance improvement by changing the structure of the BasicSymbolic, however we believe this is likely to be non-breaking. Thus the 2023-2024 stream of majors on Symbolics seems to have come to its end and Symbolics is now in more of a feature building phase.
There are two kinds of adjoints, one is the continuous adjoint approaches which define a new ODE problem to solve, and another which uses automatic differentiation directly through the solver. There are pros and cons of each, as described in our recent review on differentiation of ODEs in exquisite detail.
However, improving support for both is always on the menu. With discrete adjoints, AD support directly through the solver has been supported by ForwardDiff, ReverseDiff, and Tracker since almost the dawn of DifferentialEquations.jl. Limited support for Zygote through specific solvers, such as those in SimpleDiffEq.jl, has also always existed as well. But Enzyme is a much more powerful system. We have used it rather extensively in the continuous adjoint infrastructure for almost half a decade now, but it always lacked the capability to directly differentiate the complexity of the ODE solvers...
Until now. As part of the JuliaCon 2024 Hackathon, the remaining issues were worked through and now the explicit methods in OrdinaryDiffEq.jl can be directly differentiated with Enzyme. This can be seen by using the AD passthrough setting to avoid the SciMLSensitivity adjoint catching as follows:
using Enzyme, OrdinaryDiffEq, StaticArrays
function lorenz!(du, u, p, t)
du[1] = 10.0(u[2] - u[1])
du[2] = u[1] * (28.0 - u[3]) - u[2]
du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
const _saveat = SA[0.0,0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.0,2.25,2.5,2.75,3.0]
function f(y::Array{Float64}, u0::Array{Float64})
tspan = (0.0, 3.0)
prob = ODEProblem{true, SciMLBase.FullSpecialize}(lorenz!, u0, tspan)
sol = DiffEqBase.solve(prob, Tsit5(), saveat = _saveat, sensealg = DiffEqBase.SensitivityADPassThrough())
y .= sol[1,:]
return nothing
end;
u0 = [1.0; 0.0; 0.0]
d_u0 = zeros(3)
y = zeros(13)
dy = zeros(13)
Enzyme.autodiff(Reverse, f, Duplicated(y, dy), Duplicated(u0, d_u0));This is a major leap forward in Enzyme support. SciMLSensitivity will soon update to include a EnzymeAdjoint option which makes use of this direct differentiation mode with automation of the process (such as unwrapping of function pointers for SciMLBase.FullSpecialize, so more standard definitions work). There are still a few details to work through in order to get compatability with all features, in particular support for ranges is the leading issue, but these are actively being worked on.
Support for implicit methods is currently lacking in this form because Enzyme is incompatible with PreallocationTools.jl structures due to being unable to being able to prove the lack of aliasing. Thus Enzyme adjoint sensitivity support for implicit methods will directly come with the OrdinaryDiffEq v7 changes to the autodiff API which is planning to change from the lagacy ForwardDiff-based Jacobian interface to using DifferentiationInterface.jl for the Jacobian specification, which would then allow for Enzyme-based Jacobians and will work due to Enzyme-over-Enzyme support. This is expected over the next month as mentioned at JuliaCon.
For many years SciMLSensitivity.jl only supported AbstractArray parameter types for its continuous adjoints because it required being able to solve differential equations based on the object type. This was relaxed a bit with the introduction of GaussAdjoint as the new standard adjoint method in 2023, but there were still limitations. Now thanks to a massive effort by Dhariya, this limitation has been lifted. Now any type which defines the SciMLStructures interface for tunable canonicalization can automatically be supported.
One major result of this is that the MTKParameters object from ModelingToolkit v9 is supported, which means that general adjoint differentiation is now compatible with MTK models. This includes compatability with the SymbolicIndexingInterface, meaning that lazy observed quantities can be differented with respect to, allowing all of the symbolic simplifications of MTK to be used within the context of AD.