Thank you @mikmoore and @gdalle for your great answers!
Regarding storing functions: The collections are indeed low-dimensional, e.g. <10. The reason I explored using SVectors is to be able to access the number derivative functions for compiler to reason about when iterating over, as discussed in this post, but I guess one can simply use bare Tuples as suggested in this post directly without the extra Static Array structure on top. I guess another way would be to use FunctionWrappers.jl for the target_function, and use the type of the output to inform compiler about the length of the derivatives Tuple (since if the codomain of the target_function is \mathbb{R}^n, then the number of derivatives, i.e. the length of
Regarding Symbolic D vs AD: I think the package would benefit from AD, since the target_functions are very nonlinear, and in most cases would lead to “expression swell” as you mentioned. As an example of a target_function (for some fixed parameters) h \colon \mathbb{R}^n \to \mathbb{R}^n:
h(y)=(h_1(y),...,h_n(y)), such that
h_i(y) = \int_X f_i(x,y)g(x) dx, where
f_i(x,y)=\frac{e^{1 - \lVert x- y_i \rVert}}{1+\sum_j{e^{1 - \lVert x- y_j \rVert}}} and g(x) is some pdf.
and the derivatives of interest are \frac{\partial^k h_i}{\partial y_i^k} for k \in {1,2}.
The use case: The practice is a multi-objective optimization: finding local Nash equilibria in a game theoretic setting. So the aim is to find the roots of \tilde{h}(y)=(\frac{\partial h_1}{\partial y_1},...,\frac{\partial h_n}{\partial y_n}) such that \frac{\partial^2 h_i}{\partial y_i^2}<0 \, \forall i (ignoring other cases for now). One approach that I have been trying is to use interval arithmetics (à la JuliaIntervals by @dpsanders et al.), especially using their IntervalConstraintProgramming.jl(*), however the current version seems to be not compatible with Dual type propagation needed to couple it with something like ForwardDiff.jl, since it internally converts constraints to ModelingToolkit.Operation types for which the Dual is not defined, and I couldn’t figure out how to make it work yet. (I guess I should open a new topic for this since it’s a different issue)
Long story short, that’s why I reverted to using symbolic differentiation using Symbolics.jl.
(*) The reason that I try to implement IntervalConstraintProgramming is that each y_i might be allowed to be vector-valued, e.g. points in 2D space, and thus the function for which the roots to be found becomes \tilde{h} \colon \mathbb{R}^{2n} \to \mathbb{R}^n, disallowing methods like Newton due to non-invertable Jacobian (as far as I understood).