Error in multiply matrices with values (Cannot multiply with an MLinExpr from the left)
AnsweredHello friends,
I would like to request help with the following problem, please.
I need to multiply three matrices, as shown in the following model:
$$\min_{\textbf{x}} W\textrm{XC} + S\textrm{XR}$$
$$ \textrm{s.t.} \hspace{0.2cm} x_{ij} \in \{0, 1\}, \hspace{0.5cm} \sum_{j=1}^{n+1}=1, i=1, 2, ..., m. $$
Matrix X is a 0-1 matrix that must be optimized.
I'm considering the following dimensions:
W, S = (2, 2)
X = (2, 4)
C, R = (4, 4)
I created the following code (and already tried several possibilities), however I get the message "Cannot multiply with an MLinExpr from the left":
#!/usr/bin/python
import numpy as np
import scipy.sparse as sp
import gurobipy as gp
from gurobipy import GRB
import math
m = 2
n = 4
try:
model = gp.Model("matrix2")
W = np.matrix([[np.sqrt(3000000000.0), np.sqrt(1500000000.0)],
[0.0, 0.0]])
C = np.matrix([[(1.0/np.sqrt(2000000000.0)), 0.0, 0.0, 0.0],
[0.0, (1.0/np.sqrt(5000000000.0)), 0.0, 0.0],
[0.0, 0.0, (1.0/np.sqrt(1500000000.0)), 0.0],
[0.0, 0.0, 0.0, (1.0/np.sqrt(1000000000.0))]])
X = {}
X = model.addMVar((m, n), vtype=GRB.BINARY, name="X")
model.update()
model.setObjective((sum((W[i, :] @ X[:, j]) @ C[:, j] for i in range(m) for j in range(n))), GRB.MINIMIZE)
for i in range(m):
model.addConstr(X[i, :].sum() <= 1, name="row"+str(i))
model.optimize()
except gp.GurobiError as e:
print('Error code ' + str(e.errno) + ": " + str(e))
except AttributeError:
print('Encountered an attribute error')
How can I perform this multiplication?
Thanks in advance.
-
Official comment
-
Gurobi Staff
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?. -
-
-
The objective function should be a scalar function. That is, after plugging in values for the variables, the result should be a single number. This is not the case here - for \( W \in \mathbb{R}^{2 \times 2} \), \( X \in \{0,1\}^{2 \times 4} \), and \( C \in \mathbb{R}^{4 \times 4} \), \( WXC \) is a \( 2 \times 4 \) matrix. The same holds for \( SXR \). Could you describe mathematically what the objective function is trying to do?
0 -
-
Dear Eli,
Thanks for the feedback.
The matrix X is a positional matrix: in each row, value 1 means allocated task and value 0 is unallocated task, in one of the available resources (in each column).
The scalar must be the minimum possible time to process a task, among the various possible combinations of matrix X. Each resource has distinct cycles (W) and processing (C).
The scalar is obtained through \( \|W\textrm{XC}^T\|_2^2 \) .
So I would like to get the optimal scalar and matrix X whose combination 0-1 leads to the minimum value.
Ex.:
All combinations of X, following restrictions:
[1 0 0 0, 1 0 0 0]; [1 0 0 0, 0 1 0 0]; [0 1 0 0, 0 1 0 0]; [0 1 0 0, 1 0 0 0]; [0 0 1 0, 1 0 0 0]; [0 0 1 0, 0 1 0 0]; [1 0 0 0, 0 0 0 1]; [0 1 0 0, 0 0 0 1]; [0 0 1 0, 0 0 0 1]
Best choice: [0 1 0 0, 1 0 0 0]\( \|W\textrm{XC}^T\|_2^2 = \)
\(
\begin{bmatrix}
\sqrt{3 \times 10^9} & \sqrt{1,5 \times 10^9}\\
0 & 0
\end{bmatrix}
\begin{bmatrix}
0 & 1 & 0 & 0\\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\frac{1}{\sqrt{2 \times 10^9}} & 0 & 0 & 0\\
0 & \frac{1}{\sqrt{5 \times 10^9}} & 0 & 0\\
0 & 0 & \frac{1}{\sqrt{1,5 \times 10^9}} & 0\\
0 & 0 & 0 & \frac{1}{\sqrt{1 \times 10^9}}
\end{bmatrix} \) = 1,35 sec. (Miminum possible value)P.S.: I got a working version using numpy and itertools that chooses the matrix X and gets the scalar, but I need a version using solver without needing to use "brute force".
0 -
-
-
I see, thanks for the clarification. Currently, many operations in the Python matrix API involving MVar objects are limited to 1-D MVar objects. So I would introduce auxiliary variables to represent the elements of \( WX \) and \( WXC \). This makes it much easier to write all of the constraints:
import gurobipy as gp
from gurobipy import GRB
import numpy as np
m = 2
n = 4
model = gp.Model()
W = np.array([[np.sqrt(3000000000.0), np.sqrt(1500000000.0)],
[0.0, 0.0]])
C = np.array([[(1.0/np.sqrt(2000000000.0)), 0.0, 0.0, 0.0],
[0.0, (1.0/np.sqrt(5000000000.0)), 0.0, 0.0],
[0.0, 0.0, (1.0/np.sqrt(1500000000.0)), 0.0],
[0.0, 0.0, 0.0, (1.0/np.sqrt(1000000000.0))]])
X = model.addMVar((m, n), vtype=GRB.BINARY, name="X")
# Auxiliary variables
WX = model.addMVar((m, n), name='WX')
WXC = model.addMVar((m, n), name='WXC')
# Set WX variables equal to W @ X
model.addConstrs((WX[i, j] == W[i, :] @ X[:, j] for i in range(m) for j in range(n)), name='set_wx')
# Set WXC variables equal to WX @ C
model.addConstrs((WXC[i, j] == WX[i, :] @ C[:, j].reshape((n,)) for i in range(m) for j in range(n)), name='set_wxc')
# At least one nonzero per row
model.addConstrs((X[i, :].sum() >= 1 for i in range(m)), name='row')
# Minimize squared Frobenius norm of WXC
model.setObjective(sum(WXC[i, :] @ WXC[i, :] for i in range(m)), GRB.MINIMIZE)
model.optimize()Based on your description, I think the \( \texttt{row} \) constraints should be \( \geq 1 \) instead of \( \leq 1 \). Otherwise, the optimal solution is to set every \( X \) variable to \( 0 \).
The above code gives us an optimal solution of \( 1.35 \):
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 18 rows, 24 columns and 40 nonzeros
Model fingerprint: 0xca28091a
Model has 8 quadratic objective terms
Variable types: 16 continuous, 8 integer (8 binary)
Coefficient statistics:
Matrix range [1e-05, 5e+04]
Objective range [0e+00, 0e+00]
QObjective range [2e+00, 2e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Presolve removed 12 rows and 12 columns
Presolve time: 0.00s
Presolved: 6 rows, 12 columns, 20 nonzeros
Presolved model has 4 quadratic objective terms
Variable types: 4 continuous, 8 integer (8 binary)
Found heuristic solution: objective 4.3713203
Root relaxation: objective 9.202780e-01, 20 iterations, 0.00 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 0.92028 0 4 4.37132 0.92028 78.9% - 0s
H 0 0 2.3000000 0.92028 60.0% - 0s
H 0 0 1.8000000 0.92028 48.9% - 0s
0 0 0.92028 0 5 1.80000 0.92028 48.9% - 0s
0 0 1.02855 0 4 1.80000 1.02855 42.9% - 0s
0 0 1.02855 0 4 1.80000 1.02855 42.9% - 0s
H 0 0 1.3500000 1.02855 23.8% - 0s
0 0 infeasible 0 1.35000 1.35000 0.00% - 0s
Explored 1 nodes (31 simplex iterations) in 0.01 seconds
Thread count was 4 (of 4 available processors)
Solution count 4: 1.35 1.8 2.3 4.37132
Optimal solution found (tolerance 1.00e-04)
Best objective 1.350000000000e+00, best bound 1.350000000000e+00, gap 0.0000%0 -
-
Dear Eli,
Great explanation!
I have one last question to make the final calculation.
\( S\textrm{XR} \) is a scalar obtained with a simple diagonal trace, there is no need for the Forbenius norm calculation.
The end result should be:\( \|(W\textrm{XC})^T\|_2^2 + S\textrm{XR} = 1,35 + 0,09 = 1,44 \hspace{0.2cm} \textrm{sec.}\)
For:
W, S = (2, 2); X = (2, 4); C = (4, 4); R = (4, 2)I tried to follow your proposed lines, but with issues:
#!/usr/bin/env python3.7
import numpy as np
import scipy.sparse as sp
import gurobipy as gp
from gurobipy import GRB
import math
m = 2
n = 4
model = gp.Model()
W = np.array([[np.sqrt(3000000000.0), np.sqrt(1500000000.0)],
[0.0, 0.0]])
C = np.array([[(1.0/np.sqrt(2000000000.0)), 0.0, 0.0, 0.0],
[0.0, (1.0/np.sqrt(5000000000.0)), 0.0, 0.0],
[0.0, 0.0, (1.0/np.sqrt(1500000000.0)), 0.0],
[0.0, 0.0, 0.0, (1.0/np.sqrt(1000000000.0))]])
S = np.array([[60000000.0, 30000000.0],
[0.0, 0.0]])
R = np.array([[(1.0/1000000000.0), (1.0/1000000000.0)],
[(1.0/1000000000.0), (1.0/1000000000.0)],
[0.0, 0.0],
[0.0, 0.0]])
X = model.addMVar((m, n), vtype=GRB.BINARY, name="X")
# Auxiliary variables
WX = model.addMVar((m, n), name='WX')
WXC = model.addMVar((m, n), name='WXC')
SX = model.addMVar((m, n), name='SX')
SXR = model.addMVar((m, m), name='SXR')
SXR_trace = model.addVar(lb=0, vtype="C", name="SXR_trace")
model.addConstrs((WX[i, j] == W[i, :] @ X[:, j] for i in range(m) for j in range(n)), name='set_wx')
model.addConstrs((SX[i, j] == S[i, :] @ X[:, j] for i in range(m) for j in range(n)), name='set_sx')
model.addConstrs((WXC[i, j] == WX[i, :] @ C[:, j].reshape((n,)) for i in range(m) for j in range(n)), name='set_wxc')
model.addConstrs((SXR[i, i] == SX[i, :] @ R[:, i] for i in range(m)), name='set_sxr')
model.addConstrs((SXR_trace == np.trace(SXR[i, i], axis1=1, axis2=2) for i in range(m)), name='set_sxr_trace')
model.addConstrs((X[i, :].sum() >= 1 for i in range(m)), name='row')
model.setObjective(sum(WXC[i, :] @ WXC[i, :] + SXR_trace for i in range(m)), GRB.MINIMIZE)
model.optimize()Errata:
In the first post, I put R = (4, 4), but R = (4, 2).
In the first answer, I put \( \|W\textrm{XC}^T\|_2^2 \), but \( \|(W\textrm{XC})^T\|_2^2 \); I saw the use of reshape, then saw that the parentheses were missing.0 -
-
-
You could add a new matrix variable representing exactly the \( 2 \) diagonal elements of \( SXR \). Then, calling MVar.sum() on this object gives you \( \textrm{tr}(SXR) \).
SX = model.addMVar((m, n), name='SX')
SXR_diag = model.addMVar((m,), name='SXR_diag')
model.addConstrs((SX[i, j] == S[i, :] @ X[:, j] for i in range(m) for j in range(n)), name='set_sx')
model.addConstrs((SXR_diag[i] == SX[i, :] @ R[:, i] for i in range(m)), name='set_sxr_diag')
# Objective: min ||WXC||_F^2 + tr(SXR)
model.setObjective(sum(WXC[i, :] @ WXC[i, :] for i in range(m)) + SXR_diag.sum(), GRB.MINIMIZE)0 -
-
Works fine!
Thank you very much for help, and I congratulate you on excellent work.0 -
-
-
This is answered in a separate forum post.
0 -
Post is closed for comments.
Comments
8 comments