Model specification efficiency

[rileyjmurray]

All line numbers are with respect to dc6735c

Lines 38 to 53 can be reduced to three CVXPY constraint objects by using appropriate logical / block indexing.

Lines 56 to 65 can be reduced to a single CVXPY constraint object. The trick would be to generate a list of CVXPY Expression objects, where each element of the list was the expression defined from lines 59 to 64 (for appropriate k). Once that list is built you can turn it into one big CVXPY Expression object via "vstack"; this will give you the RHS of the constraint. You should then be able to write X_[range(K-1), :] == RHS.

The general pattern is

• Build two lists of CVXPY Expression objects "left_hand_sides", "right_hand_sides" by looping over your problem data.
• Define LHS = vstack(left_hand_sides), RHS = vstack( right_hand_sides ).
• Append "LHS <= RHS" or "LHS == RHS" to your list of CVXPY Constraint objects.

The function "vstack" can be found at cvxpy.atoms.affine.vstack.vstack( [list of expressions] ). It requires all arguments to have the same number of columns, so you may have to use the above pattern several times for different groups of constraints, depending on their dimensions.

[jonnyhyman]

Just curious, does this constitute a performance increase or a syntax improvement or both?

[rileyjmurray]

CVXPY's compile time will dramatically improve by vectorizing the constraints with the method described above. The time spent by the underlying solver (e.g. ECOS, SCS, Gurobi, MOSEK, etc...) will not be impacted.

I should clarify: when you call prob.solve() where "prob" is a CVXPY Problem, three things happen. First, CVXPY compiles the problem. Second, the compiled problem is passed off to a low-level solver (which subsequently solves the problem). Third, CVXPY transforms the low-level solver output back into the format given by the user. It is very hard to speed up this second step, but oftentimes the first and third steps are the most computationally expensive.

[EmbersArc]

Thanks for taking the time to help out! Greatly appreciated.
I've made the changes you proposed, but there doesn't seem to be a notable difference in the problem compile time. Can you quickly check whether I did it correctly?

Also, in its current form there's an error in line 79. Is the behavior of sum_squares different than that of sum in this case? It's not a big deal, I might be able to reformulate it but just curious.

[rileyjmurray]

Hmm. It's surprising that the compile time isn't meaningfully improved.

One thing you can try is to only concatenate expressions that are of the same curvature. That is- do not stack linear expressions with non-linear convex expressions such as "norm" or "sum_squares". By separating things in this way CVXPY can easily tell when a massive constraint block is linear (and hence requires little processing) as opposed to convex but non-linear.

To further improve performance, rather than vectorizing "norm" expressions via vstack, use the "axis" argument to norm. The axis argument isn't in the tutorial from what I can tell, but it does appear in the API (https://cvxgrp.github.io/cvxpy/api_reference/cvxpy.atoms.other_atoms.html#norm). Moreover, the axis argument for this cvxpy function behaves in an identical way to the axis argument in numpy.linalg.norm.

Regardling the value error on line 79, I'm not sure what's going on there. Sometimes funny things can happen when concatenating Expressions of shape (,) vs shape (1,) vs shape (1,1). For sanity's sake it's probably best to handle these constraints as you were doing before. Yes, there should be some performance benefit to rewriting this constraint, but you have to pick you battles when it comes to optimizing your code.

[EmbersArc]

Unfortunately using the norm with axis also led to some confusing error. Maybe it's the way my matrices are set up. But anyways thanks a lot! I'll try again at another point but for now we're working on the C++ implementation.

[EmbersArc]

@rileyjmurray The norm axis notation works just fine after switching around the dimensions. Before, it was X[k, i], now it's X[i, k] where k is the trajectory index and i is the state value index. Seems to be a 2x 5x speedup!