slow quadratic constraint creation in Pyomo

Question

trying to construct a large scale quadratic constraint in Pyomo as follows:

import pyomo as pyo
from pyomo.environ import *

scale   = 5000
pyo.n   = Set(initialize=range(scale))
pyo.x   = Var(pyo.n, bounds=(-1.0,1.0))
# Q is a n-by-n matrix in numpy array format, where n equals <scale>
Q_values = dict(zip(list(itertools.product(range(0,scale), range(0,scale))), Q.flatten()))
pyo.Q   = Param(pyo.n, pyo.n, initialize=Q_values)

pyo.xQx = Constraint( expr=sum( pyo.x[i]*pyo.Q[i,j]*pyo.x[j] for i in pyo.n for j in pyo.n ) <= 1.0 )

turns out the last line is unbearably slow given the problem scale. tried several things mentioned in PyPSA, Performance of creating Pyomo constraints and pyomo seems very slow to write models. but no luck.

any suggestion (once the model was constructed, Ipopt solving was also slow. but that's independent from Pyomo i guess)?

ps: construct such quadratic constraint directly as follows didnt help either (also unbearably slow)

pyo.xQx = Constraint( expr=sum( pyo.x[i]*Q[i,j]*pyo.x[j] for i in pyo.n for j in pyo.n ) <= 1.0 )

@ErlingMOSEK it is. but converting it to SOCP doesnt quite helpful since the fundamental issue here is the slow second order constraint generation in Pyomo (such large scale second order constraint requires loop/list comprehension anyway) — leo, Dec 21 '20 at 12:10
@ErlingMOSEK or u probably meant to use Schur's trick. however I'm not sure how to construct such constraint in Pyomo. Maybe u can help illustrate? — leo, Dec 21 '20 at 13:22
but unfortunately mosek is neither free nor GPL...and here it is supposed to a technically issue — leo, Dec 21 '20 at 13:55

score 1 · Answer 1 · answered Dec 22 '20 at 10:51

You can get a small speed-up by using quicksum in place of sum. To measure the performance of the last line, I modified your code a little bit, as shown:

import itertools
from pyomo.environ import *
import time
import numpy as np

scale = 5000
m = ConcreteModel()
m.n = Set(initialize=range(scale))
m.x = Var(m.n, bounds=(-1.0, 1.0))
# Q is a n-by-n matrix in numpy array format, where n equals <scale>
Q = np.ones([scale, scale])
Q_values = dict(
    zip(list(itertools.product(range(scale), range(scale))), Q.flatten()))
m.Q = Param(m.n, m.n, initialize=Q_values)

t = time.time()
m.xQx = Constraint(expr=sum(m.x[i]*m.Q[i, j]*m.x[j]
                            for i in m.n for j in m.n) <= 1.0)
print("Time to make QuadCon = {}".format(time.time() - t))

The time I measured with sum was around 174.4 s. With quicksum I got 163.3 seconds.

Not satisfied with such a modest gain, I tried to re-formulate as a SOCP. If you can factorize Q like so: Q= (F^T F), then you could easily express your constraint as a quadratic cone, as shown below:

import itertools
import time
import pyomo.kernel as pmo
from pyomo.environ import *
import numpy as np

scale = 5000
m = pmo.block()
m.n = np.arange(scale)
m.x = pmo.variable_list()
for j in m.n:
    m.x.append(pmo.variable(lb=-1.0, ub=1.0))
# Q = (F^T)F factors (eg.: Cholesky factor)
_F = np.ones([scale, scale])
t = time.time()
F = pmo.parameter_list()
for f in _F:
    _row = pmo.parameter_list(pmo.parameter(_e) for _e in f)
    F.append(_row)
print("Time taken to make parameter F = {}".format(time.time() - t))
t1 = time.time()
x_expr = pmo.expression_tuple(pmo.expression(
    expr=sum_product(f, m.x, index=m.n)) for f in F)
print("Time for evaluating Fx = {}".format(time.time() - t1))
t2 = time.time()
m.xQx = pmo.conic.quadratic.as_domain(1, x_expr)
print("Time for quad constr = {}".format(time.time() - t2))

Running on the same machine, I observed a time of around 112 seconds in the preparation of the expression that gets passed to the cone. Actually preparing the cone takes very little time (0.031 s).

Naturally, the only solver that can handle Conic constraints in pyomo is MOSEK. A recent update to the Pyomo-MOSEK interface has also shown promising speed-ups.

You can try MOSEK for free by getting yourself a MOSEK trial license. If you want to read more about Conic reformulations, a quick and thorough guide can be found in the MOSEK modelling cookbook. Lastly, if you are affiliated with an academic institution, then we can offer you a personal/institutional academic license. Hope you find this useful.

Btw using cvxpy.org instead of Pyomo with a free optimizer is likely to be much faster regarding input. cvxpy.org with Mosek will be even faster of course :-). — ErlingMOSEK, Dec 23 '20 at 06:39

slow quadratic constraint creation in Pyomo

1 Answers1