I have a set of *m* equality constraints with linear and/or bilinear terms, defined as strings. For example:

[ “`P0_I_m`

− `F1_m`

− `F2_m`

”,

“`P0_E_m`

− `F6_m`

− `F7_m`

”,

“`P0_dS_m`

+ `P0_E_m`

− `P0_I_m`

”,

“`F1_m`

+ `F2_m`

+ `F4_m`

− `F3_m`

”,

“`F3_m`

− `F4_m`

− `F5_m`

”,

“`F5_m`

− `F6_m`

− `F7_m`

”,

“`F4_m`

− `P2_T_F3_F4 * F3_m`

”,

“`F5_m`

− `P2_T_F3_F5 * F3_m`

”,

“`SUM_TC`

− `P2_T_F3_F4`

− `P2_T_F3_F5`

”]

Here, each constraint must equal 0. Altogether *n* atomic terms, e.g. F4_m, F6_m, F7_m, P2_T_F3_F5 and P0_I_m, are the variables, while others, e.g. F2_m and SUM_TC, are constants. I also need a *m by n* Jacobian matrix of the constraints w.r.t. the variables. Most its entries are constant, eg 0 or -1, while others are linear expressions in the variables, eg `-F3_m`

or `F3_m - F5_m + F6_m`

. Some of these are repeated across the matrix.

The constraint set is generated by visual modelling on a TypeScript frontend, meaning users don’t write actual equations. A problem can have up to 3000 variables, with a constraint having up to 500 terms, and there can be up to 15,000 constraints.

**Task:** Efficiently compute the constraint values and the Jacobian at each of 1 – 100 iterations of an algorithm.

This requires an initiation stage, where I parse the constraint strings into expressions that can be evaluated by Julia. Afterwards, I evaluate the expressions several times. Both stages can be time consuming, even with multi-processing. I’m wondering what the fastest approach could be (eg, by avoiding Meta.parse and eval.)

I’m currently using DynamicExpressions.jl by @MilesCranmer. The code is quite long, so I summarise the steps:

- Define
`OrderedDict{Symbol,Node{Float64}}`

for the variables (as features). - Parse the linear/bilinear terms in a constraint and define the corresponding Node expression with the “Noded” variables.
- Create a Jacobian matrix with Float64 constant entries and with NaNs for variable-dependent (“symbolic”) entries.
- Define an OrderedDict with keys being those “symbolic” entries as Nodes, and values eg [11, 45, 80], being the linear indices of the corresponding Jacobian matrix entries.
- After each iteration, evaluate the constraints and Jacobian terms, and insert the latter at the corresponding Jacobian matrix indices.

**Questions:**

- What alternative approaches could be faster (eg, using a different package)?
- Should I keep known values as nodes in expressions, eg
`Node(Float64, val = 1.0) + Node(Float64, feature = 2)`

or take the constants outside and add separately? - What can I gain performance-wise by eliminating the “-” binary and/or unary operator from
`OperatorEnum(binary_operators=[+, -, *], unary_operators=[-])`

. - When does a Node expression become too long, so that splitting it would improve performance?
- Is the individual evaluation of each of the many simple Jacobian terms, like
`-F3_m`

, efficient? Or is it better to somehow combine them into an array and evaluate jointly? - Alternatively, I could evaluate the symbolic Jacobian entries without DynamicExpressions.jl , but how do I handle multi-term entries like
`F3_m - F5_m + F6_m`

? I could use DynamicExpressions.jl for those terms only, but that’s more complicated. - If only some of the variables change values at an iteration, can I exploit this for performance?

Thank you!