How do I add a matrix constraint `Ax=b` to a Pyomo model efficiently?

Question

I want to add the constraints Ax=b to a Pyomo model with my numpy arrays A and b as efficient as possible. Unfortunately, the performance is very bad currently. For the following example

    import time
    import numpy as np
    import pyomo.environ as pyo

    start = time.time()
    rows = 287
    cols = 2765
    A = np.random.rand(rows, cols)
    b = np.random.rand(rows)
    mdl = pyo.ConcreteModel()
    mdl.rows = range(rows)
    mdl.cols = range(cols)
    mdl.A = A
    mdl.b = b
    mdl.x_var = pyo.Var(mdl.cols, bounds=(0.0, None))

    mdl.constraints = pyo.ConstraintList()
    [mdl.constraints.add(sum(mdl.A[row, col] * mdl.x_var[col] for col in mdl.cols) <= mdl.b[row]) for row in mdl.rows]

    mdl.obj = pyo.Objective(expr=sum(mdl.x_var[col] for col in mdl.cols), sense=pyo.minimize)
    end = time.time()
    print(end - start)
   

is takes almost 30 seconds because of the add statement and the huge amount of columns. Is it possible to pass A, x, and b directly and fast instead of adding it row by row?

Answer 1

The main thing that is slowing down your construction above is the fact that you are building the constraint list elements within a list comprehension, which is unnecessary and causes a lot of bloat.

This line:

[mdl.constraints.add(sum(mdl.A[row, col] * mdl.x_var[col] for col in mdl.cols) <= mdl.b[row]) for row in mdl.rows]

Constructs a list of the captured results of each ConstraintList.add() expression, which is a "rich" return. That list is an unnecessary byproduct of the loop you desire to do over the add() function. Just change your generation scheme to either a loop or a generator (by using parens) to avoid that capture, as such:

(mdl.constraints.add(sum(mdl.A[row, col] * mdl.x_var[col] for col in mdl.cols) <= mdl.b[row]) for row in mdl.rows)

And the model construction time drops to about 0.02 seconds.