Quadratic optimization too slow
Answered
Hello everyone,
I am trying to solve the following optimization problem using python interface:
Modelling the problem as in the code below returns correct results. However, it takes very long time to do so.
Is it the way how I am formulating the problem causing the bottleneck?
Thanks in advance.
# 0, average nodes
# 1, min-nodes
# 2, max-nodes
# 3 and 4 are sink(0) and sink(1)
nodes, edges = [2, 1, 2, 0, 3, 4] [(1, 4), (3, 5), (4, 3), (0, 4)]
def run(self):
# Create a model
m = gp.Model("qp")
m.Params.LogToConsole = 0
m.setParam('NonConvex', 2)
# Create variables
for i in range(len(self.nodes)):
m.addVar(lb=0.0, ub=1.0, name='x' + str(i))
m.update() # model stops working without it.
# Add the constraints
m.addConstr(m.getVars()[-1] == 1.0, 'SINK(1)')
m.addConstr(m.getVars()[-2] == 0.0, 'SINK(0)')
# Create the constraints and objective function
obj = 0
for i in range(len(self.nodes) - 2):
variable = m.getVars()[i]
node_type = self.nodes[i]
n1_index = self.edges[i][0]
n2_index = self.edges[i][1]
n1_var = m.getVars()[n1_index]
n2_var = m.getVars()[n2_index]
if node_type == 0: # AVE
m.addConstr(variable == self.ave_probs[i][0] * n1_var + self.ave_probs[i][-1] * n2_var, name='AVE')
if node_type == 1:
m.addConstr(variable - n1_var <= 0, 'Min_1')
m.addConstr(variable - n2_var <= 0, 'Min_2')
obj += (variable ** 2 - (variable * n1_var) - (variable * n2_var) + (n1_var * n2_var))
if node_type == 2:
m.addConstr(-variable + n1_var <= 0, 'Max_1')
m.addConstr(-variable + n2_var <= 0, 'Max_2')
obj += (variable ** 2 - (variable * n1_var) - (variable * n2_var) + (n1_var * n2_var))
m.setObjective(obj)
m.optimize()
v = []
# Print the variables names and values
for var in m.getVars():
v.append(round(var.X, 3))
return v
-
I have followed the advices from https://support.gurobi.com/hc/en-us/articles/360013420111-How-do-I-improve-the-time-to-build-my-model-#:~:text=creating%20new%20expressions.-,Matrix%20API,-When%20used%20properly. So I have modeled the problem in matrix form, but it doesn't seem to make much difference. Could somebody help me improve the speed? I have tried every thing I could to remove any excessive computation steps and that's why the code now looks less readable (sorry for that).
Minimize v(x)² - v(x)v(y) - v(x)v(z) + v(y)v(z)
Constrains:
IF x from V_MAX
-v(x) <= v(y) <=> -v(x) + v(y) <= 0
-v(x) <= -v(z) <=> -v(x) + v(z) <= 0
IF x from V_MIN
v(x) <= v(y) <=> v(x) - v(y) <= 0
v(x) <= v(z) <=> v(x) - v(z) <= 0
IF x from V_AVE
v(x) = 0.5 * v(y) + 0.5 * v(z) <=> v(x) - 0.5 * v(y) - 0.5 * v(z) = 0
IF x is a SINK(0)
v(x) = 0
IF x is a SINK(1)
v(x) = 1
For all vertices
0 <= v(x) <= 1, <=> v(x) <= 1 & -v(x) <= 0
0 <= v(y) <= 1, <=> v(y) <= 1 & -v(y) <= 0
0 <= v(z) <= 1, <=> v(z) <= 1 & -v(z) <= 0
AVE = 0, MIN = 1, MAX = 2, SINK(0) = 3, SINK(1) = 4
"""
class QP2:
def run(self, nodes, edges, ave_probs):
# print(60 * '-' + '|Quadratic Programming|----' + 60 * '-')
# Create a model
m = gp.Model("qp")
m.Params.LogToConsole = 0
m.setParam('NonConvex', 2)
n = len(nodes)
x = m.addMVar(shape=n, ub=1.0)
# Matrices for Equality constraints
eq_count = nodes.count(0) + 2 # sinks
Aeq = [[0.] * n for _ in range(eq_count)]
beq = ([0.] * (eq_count-1))
beq.append(1.0)
# Objective Function = x² - xz - xy + yz where y, z are the neighbors of x.
H = [[0.] * n for _ in range(n)]
# Matrices for Equality constraints
uneq_count = (n - eq_count) * 2
A = [[0.] * n for _ in range(uneq_count)]
b = [0.] * uneq_count
# Populating Aeq, beq, A, b, H
c_row = 0
row = 0
for i in range(n):
if nodes[i] == 0: # AVE
y = edges[i][0]
z = edges[i][1]
Aeq[c_row][i] = 1
Aeq[c_row][y] -= ave_probs[i][0]
Aeq[c_row][z] -= ave_probs[i][1]
c_row += 1
elif nodes[i] == 2: # MAX
y = edges[i][0]
z = edges[i][1]
A[row][i] = -1
A[row][y] += 1
row += 1
A[row][i] = -1
A[row][z] += 1
row += 1
# Generate H for MAX
H[i][i] += 1 # x²
# Subtract the value of the neighbors.
# if wie had xy in the expression we have to make it 2 * xy because we will divide on 2.
H[i][y] -= 0.5
H[y][i] -= 0.5
H[i][z] -= 0.5
H[z][i] -= 0.5
# because of yz
H[y][z] += 0.5
H[z][y] += 0.5
elif nodes[i] == 1: # MIN
y = edges[i][0]
z = edges[i][1]
A[row][i] = 1
A[row][y] -= 1
row += 1
A[row][i] = 1
A[row][z] -= 1
row += 1
# Generate H for MIN
H[i][i] += 1 # x²
# Subtract the value of the neighbors.
# if wie had xy in the expression we have to make it 2 * xy because we will divide on 2.
H[i][y] -= 0.5
H[y][i] -= 0.5
H[i][z] -= 0.5
H[z][i] -= 0.5
# because of yz
H[y][z] += 0.5
H[z][y] += 0.5
elif nodes[i] == 3: # SINKS
Aeq[c_row][i] = 1
c_row += 1
elif nodes[i] == 4:
Aeq[c_row][i] = 1
A = np.array(A)
b = np.array(b)
Aeq = np.array(Aeq)
beq = np.array(beq)
H = np.array(H)
m.addConstr(A @ x <= b)
m.addConstr(Aeq @ x == beq)
m.setObjective(x @ H @ x)
m.optimize()
v = [round(x.X, 3) for x in m.getVars()]
# Write the constraints to a file to see if they are correct.
# m.write("model.lp") #todo compare with old model
return v
if __name__ == '__main__':
nodes, edges = [1, 2, 0, 1, 3, 4], [(5, 4), (2, 4), (5, 1), (5, 4)]
non_sinks = nodes[0:-2]
# Create average probabilities
ave_probs = [[0.0, 0.0]] * len(non_sinks)
for i in range(len(non_sinks)):
if nodes[i] == 0:
n1_prob = 0.5
n2_prob = 0.5
assert n1_prob + n2_prob == 1.0
ave_probs[i] = [n1_prob, n2_prob]
x = QP2().run(nodes, edges, ave_probs)
print(x)0 -
Hi Mohamad,
I ran your code, but found that the model is formulated and solved in a fraction of a second. So it is hard to say where the bottleneck is. I assume you are running this with a much larger dataset in practice?
Are you sure the model formulation time is the bottleneck? Non-convex quadratic models can take a long time to solve. I would suggest not suppressing the solver output, i.e. remove this line:
m.Params.LogToConsole = 0
so that you can see when your formulation code is completed and the solve run starts.
0 -
Hi Simon,
as you mentioned, I want to run the model on larger number of nodes. I understand that such problems might take a long time, but starting from 30 nodes (variables) it starts to take a long time and the PC becomes very hot!
I even tried to run the script outside PyCharm to avoid unnecessary overhead. As you suggested, I don't think the formulation is the problem, I am just trying to improve everything possible. Do you think that changing parameters would help increase the speed?Even though I increased the tolerance:
It seems to sometimes get stuck in this step:
Should I maybe make the model accept greater gap value?
Thanks in advance
0 -
The following set of parameters seems to make the optimization much faster:
Could you please explain to me what exactly the SolutionLimit parameter does? (The docs were not clear for me as a beginner)
Is there a parameter that I can set to tell the optimizer to give the current best solution after a specific amount of time?0 -
Hi Mohamad,
Since the Gurobi logging output has started, your formulation code is complete and the model is being solved (i.e. the model.optimize() statement has been reached). So I don't think you have an issue with the time taken the formulate the model here. I would suggest sticking with the simpler version of your code to avoid confusion.
Evidently, Gurobi is spending a non-trivial amount of time solving the model. This is not unusual; some models are computationally more challenging than others to solve to optimality. So, you need to think carefully about the termination criteria for this model:
- How much time do you want the solver to spend? If you set the TimeLimit parameter, Gurobi will stop solving once the given time has elapsed, and return the best known feasible solution at that time. However, there is no guarantee of the quality of solutions found after this time; it may even occur that no feasible solution was found in this time (this doesn't seem to be a problem in your case, the first incumbent is found quickly).
- Do you require solutions of a certain quality? If you set the MIPGap parameter, Gurobi will stop solving once this gap value has been reached. This guarantees solutions of the required quality, but there is no guarantee of the time taken to find such solutions.
SolutionLimit is not one I would recommend using at first. It will terminate the solve after the given number of unique solutions are found, regardless of their objective value. This setting applies no solve time or solution quality criteria, so is not a good one to use unless you understand very well how the solver performs on your specific models.
I suggest reading the following articles first to get a better understanding of this:
- Read What is the MIPGap? to understand MIPGap as a solution quality metric
- Read MIP Logging to understand the information Gurobi is reporting in the log
0
Please sign in to leave a comment.
Comments
5 comments