Fast evaluation of linear/bilinear expressions and Jacobians from strings

General Usage Performance
1

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:

  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 , 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.
  7. If only some of the variables change values at an iteration, can I exploit this for performance?

Thank you!

2

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).

1 Like