•    My current short term solution is to call this prune function occasionally but its pretty sloppy. Each array will go up and down in ram usage by 10x every function call. Also it is clearly pretty shape specific, but that would be solvable if I wanted. The other solution is to break my model into hundreds of intermediate steps. With each intermediate step having a constraint to equal the variable used in the next step. Best case though this stops the 10^50 factor of ram waste and caps it at likely a factor of only 10. (numbers being how much ram is wasted per useful piece of ram)

•  Gurobi Staff

Hello Griffin,

the bloat you are observing comes from the exponentially growing size of the symbolic representations of the linear expressions. Looking at your the example code (modified for brevity):

import gurobipy as gp
import numpy as np
import time

dims = 10
m = gp.Model("mip")
m.update()
matrix = np.eye(dims)

for i in range(7):
start = time.time()
u = np.dot(matrix, u)
nterms = sum(expr.size() for expr in u)
print("Iteration {:2d}, time={:5.2f}sec, nterms={:10d}".format(
i, time.time() - start, nterms))

On my machine this gives the following output:

Iteration 0, time= 0.00sec, nterms= 100
Iteration 1, time= 0.00sec, nterms= 1000
Iteration 2, time= 0.00sec, nterms= 10000
Iteration 3, time= 0.01sec, nterms= 100000
Iteration 4, time= 0.15sec, nterms= 1000000
Iteration 5, time= 1.53sec, nterms= 10000000
Iteration 6, time=26.03sec, nterms= 100000000

At each iteration u is an ndarray of (ten) LinExpr objects, but the multiplication with 'matrix' increases the representation size of each such linear expression by a factor of ten, leading eventually to an exponentially growing run time. (Obviously the LinExpr class could be implemented smarter to collect terms, remove zero coefficients, etc. -- maybe in the future).

This can easily be avoided if you change the order of evaluation from, say, A*(A*(A*x)) to ((A*A)*A)*x, so that all multiplications except for the last one is among ndarrays.  For the above example, this could be done as follows:

import gurobipy as gp
import numpy as np
import time

dims = 10
m = gp.Model("mip")
m.update()
matrix = np.eye(dims)

# Compute product of matrices
for i in range(7):
start = time.time()
matrix = np.dot(matrix, matrix)
print("Iteration {:2d}, time={:5.2f}sec".format(
i, time.time() - start))

# Now multiply the product with the variables
start = time.time()
u = np.dot(matrix, u)
nterms = sum(expr.size() for expr in u)
print("Final result: time={:5.2f}sec, nterms={:10d}".format(
time.time() - start, nterms))

which gives on my machine

Iteration 0, time= 0.00sec
Iteration 1, time= 0.00sec
Iteration 2, time= 0.00sec
Iteration 3, time= 0.00sec
Iteration 4, time= 0.00sec
Iteration 5, time= 0.00sec
Iteration 6, time= 0.00sec
Final result: time= 0.00sec, nterms= 100

Hope that helps!

Robert

•    We are totally on the same page on what is going on. Thank you for explaining it so much more succinctly than I did. The example case A*A*A*A*A*A*A*x is pretty contrived and is easily solved by precomputing the non symbolic components. I am a less confident in the ability to do that on my real model, at least in all cases.

I have been reorganizing my model such that I no longer have the exponential explosion. My task is a motion planner with lots of integer components. In my original formulation I was trying to minimize X at time step 50. Where for each timestep I declare a set of gurobi decision variables and pass the state and the decision through a simulator. Thus X at time step 50 had been passed through the simulator 50 times. In my new formulation I just say minimize X(50) with constraints for X(t+1) = simulator(X(t),decision). This should be an equivalent formulation however it means a single expression only ever passes through my simulator once. I am now able to actually run the optimization instead of crashing trying to formulate it.

I would however like to request this update being added to the roadmap (Obviously the LinExpr class could be implemented smarter to collect terms, remove zero coefficients, etc. -- maybe in the future).  We can still see from your example

Iteration 0, time= 0.00sec, nterms= 100

that 10 terms turned into 100 even when we precompute to the best of our ability. I would hope the Presolve step does this simplification somewhere? Or there is some other reason to expect these added terms don't also increase the optimization time?

Lastly I just would like to justify why I think this numpy approach is advantageous. I have existing functions designed without symbolic math in mind. For example a simulator that I use when evaluating the decision variables provided by my optimization. By loading them into numpy arrays I am able to pass the numpy arrays into the same function and instead of getting a float output, I get a symbolic one. If I switched to using special gurobi functions for handling matrix operations I would need to maintain 2 copies of the same function with all the obvious downsides that provides. I will need to look into the gurobi 9 matrix API to see if I can get the same benefit. (PS, thanks for modifying code for brevity, you clearly know more about your API than I do)

•  Gurobi Staff

Hello Griffin,

I have created a feature request to deal with term collection and dropping zero terms in LinExpr objects.  For sure such zero terms are disposed of some time before the actual optimization starts, but not while LinExpr objects are created.

For now you could work around this by explicitly simplifying the exploding LinExpr objects after each multiplication youself.  Here is an example how to do that.  Note that this technique is not very efficient, but if the goal is only to avoid exponential run time, even slow code shines.

import gurobipy as gp
import numpy as np
import time
from collections import defaultdict

def simplify_linexpr(le):
# Collect terms and drop zeros in the given gp.LinExpr object.
coeflist = [le.getCoeff(i) for i in range(le.size())]
varlist = [le.getVar(i) for i in range(le.size())]
constant = le.getConstant()
le.clear() # Remove all terms, and constant

# Rebuild le without zeros and collected terms
merged_pairs = defaultdict(lambda: 0.0)
for (coef, var) in zip(coeflist, varlist):
if coef != 0.0:
merged_pairs[var] += coef

dims = 10
m = gp.Model("mip")
m.update()
matrix = np.eye(dims)

for i in range(7):
start = time.time()
u = np.dot(matrix, u)

for expr in u:
simplify_linexpr(expr)

nterms = sum(expr.size() for expr in u)
print("Iteration {:2d}, time={:5.2f}sec, nterms={:10d}".format(
i, time.time() - start, nterms))

On my machine I now see

Iteration  0, time= 0.00sec, nterms=        10
Iteration 1, time= 0.00sec, nterms= 10
Iteration 2, time= 0.00sec, nterms= 10
Iteration 3, time= 0.00sec, nterms= 10
Iteration 4, time= 0.00sec, nterms= 10
Iteration 5, time= 0.00sec, nterms= 10
Iteration 6, time= 0.00sec, nterms= 10

I understand that this will be cumbersome to integrate if you want to maintain a single code path that eats both ndarrays of Var objects, and normal ndarrays, but this is as close as you can get by simple means.

Robert