# Fast evaluation of linear/bilinear expressions and Jacobians from strings

**URL:** https://discourse.julialang.org/t/fast-evaluation-of-linear-bilinear-expressions-and-jacobians-from-strings/120721
**Category:** Performance
**Tags:** question
**Created:** [September 30, 2024, 12:40pm UTC](https://discourse.julialang.org/t/fast-evaluation-of-linear-bilinear-expressions-and-jacobians-from-strings/120721 "2024-09-30T12:40:25Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![Olegg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olegg/32/51316_2.png) [@Olegg](https://discourse.julialang.org/u/Olegg)
#### Post date: [September 30, 2024, 12:40pm UTC](https://discourse.julialang.org/t/fast-evaluation-of-linear-bilinear-expressions-and-jacobians-from-strings/120721/1 "2024-09-30T12:40:25Z")

</div>

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](https://juliahub.com/ui/Packages/General/DynamicExpressions) by @MilesCranmer. The code is quite long, so I summarise the steps:

1. Define `OrderedDict{Symbol,Node{Float64}}` for the variables (as features).
2. Parse the linear/bilinear terms in a constraint and define the corresponding Node expression with the “Noded” variables.
3. Create a Jacobian matrix with Float64 constant entries and with NaNs for variable-dependent (“symbolic”) entries.
4. 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.
5. After each iteration, evaluate the constraints and Jacobian terms, and insert the latter at the corresponding Jacobian matrix indices.

**Questions:**

1. What alternative approaches could be faster (eg, using a different package)?
2. 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?
3. What can I gain performance-wise by eliminating the “-” binary and/or unary operator from `OperatorEnum(binary_operators=[+, -, *], unary_operators=[-])`.
4. When does a Node expression become too long, so that splitting it would improve performance?
5. 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?
6. Alternatively, I could evaluate the symbolic Jacobian entries without [DynamicExpressions.jl](https://juliahub.com/ui/Packages/General/DynamicExpressions) , but how do I handle multi-term entries like `F3_m - F5_m + F6_m`? I could use [DynamicExpressions.jl](https://juliahub.com/ui/Packages/General/DynamicExpressions) for those terms only, but that’s more complicated.
7. If only some of the variables change values at an iteration, can I exploit this for performance?

Thank you!

---

<div class="post-metadata">

### Author: ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [October 1, 2024, 5:33pm UTC](https://discourse.julialang.org/t/fast-evaluation-of-linear-bilinear-expressions-and-jacobians-from-strings/120721/2 "2024-10-01T17:33:23Z")

</div>

I think I would just use SparseArrays.jl for this. So you could parse each constraint into a SparseMatrix like C^i = v\_k M^i\_{kl} v\_l where

- C^i is the value of the i th constraint
- \vec v = \begin{pmatrix}1&v\_1&v\_2&\ldots\end{pmatrix}^T is the state vector, i.e. the vector of the values of your variables, augmented by a constant 1 in the first entry (to allow for modeling linear terms in M\_{ij})
- M^i the matrix defining the i th constraint. These you parse from the string input essentially.

Then all your desired quantities are just vector-(sparse)matrix multiplications:

- The constraints are shown above C^i = v\_k M^i\_{kl} v\_l
- The derivative of the i th constraint w.r.t to the k th variable is essentially just a matrix-vector product with the symmetrized M^i: \partial\_j C^i = v\_k M^i\_{kj}+M^i\_{jk}v\_k = (M^i\_{jk}+(M^i)^T\_{jk})v\_k . There is a possible optimization here if you take care that the M^i are symmetric or if you precompute the inner matrix (still sparse).

I don’t think you can be much faster except if you somehow have more structure in your system (e.g. if you only have very few bilinear constraint then it could make sense to treat them separately).
