我想使用cvxpy解决以下优化问题:
Maximize(y @ A @ x)
，其中A是一个(n,m)大小的正数值元素矩阵;X是长度为m的布尔向量;Y是长度为n的布尔向量。
约束条件是:sum(x)==1和sum(y)==1
例如:A = [[ 87, 96, 127, 46], [155, 166, 92, 11], [111, 163, 126, 112]]那么解决方案是:x = [0,1,0,0]y = [0,1,0]结果是:y @ A @ x = 166，这是给定约束条件下的最大值。
但是，cvvxpy抛出错误并说:
DCPError:问题不遵循DCP规则。具体地说:目标不是DCP。它的以下子表达式不是:Var27 @[[87]。96. 127. 46岁。)(155年。166. 92. 11。)(111年。163. 126. 112.]] @ var26
表达式y @ A @ x求值为表达式(UNKNOWN, non - negative，())，表明曲率是未知的，即使A @ x和y @ A都是仿射和non - negative。
使用cvvxpy解决这个看似简单的问题的方法是什么?
我的代码是:import cvxpy as cp
import numpy as np
rows = 3
cols = 4
A = np.random.randint(0,255,size=(rows,cols))
x = cp.Variable(cols, boolean=True)
y = cp.Variable(rows, boolean=True)
objective = cp.Maximize(y @ A @ x)
constraints = [cp.sum(x)==1,cp.sum(y)==1]
prob = cp.Problem(objective,constraints)
prob.solve()
我试着做以下事情来给随机矩阵一个属性:
a = np.random.randint(0,255,size=(rows,cols))
A = cp.Variable((rows,cols),nonneg=True)
A.value = A.project(a)
我试着让它与PSD=真方阵。试过了,pos = True。我还尝试设置probe .solve(qcp=True)。所有这些仍然导致上述关于DCP的错误信息。
非线性优化
如果您决定不进行线性化，并且您可以承担某些类别A的不精确解的风险，那么
import numpy as np
from scipy.optimize import minimize, LinearConstraint, Bounds
A = np.array([
[ 87, 96, 127, 46],
[155, 166, 92, 11],
[111, 163, 126, 112],
])
n, m = A.shape
def objective(xy: np.ndarray) -> float:
x = xy[:m]
y = xy[-n:]
return -y @ A @ x
Ac = np.zeros((2, m+n))
Ac[0, :m] = 1
Ac[1, -n:] = 1
x0 = np.empty(m + n)
x0[:m] = 1/m
x0[-n:] = 1/n
result = minimize(
fun=objective, x0=x0,
bounds=Bounds(lb=np.zeros_like(x0), ub=np.ones_like(x0)),
constraints=LinearConstraint(A=Ac, lb=(1,1), ub=(1,1)),
)
print(result)
print('x ~', result.x[:m].round(3))
print('y ~', result.x[-n:].round(3))
print('Constraint residual:', Ac.dot(result.x) - 1)
message: Optimization terminated successfully
success: True
status: 0
fun: -166.0
x: [ 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00
1.000e+00 0.000e+00]
nit: 7
jac: [-1.550e+02 -1.660e+02 -9.200e+01 -1.100e+01 -9.600e+01
-1.660e+02 -1.630e+02]
nfev: 48
njev: 6
x ~ [0. 1. 0. 0.]
y ~ [0. 1. 0.]
Constraint residual: [0. 0.]
匹配您演示的输出。
线性规划
作为线性程序，
import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
A = np.array([
( 87, 96, 127, 46),
(155, 166, 92, 11),
(111, 163, 126, 112),
])
n, m = A.shape
# Objective coefficients: maximize the n*m implied products in A
c = np.full(n + m + n*m, -1)
c[:n+m] = 0
# y, x are integral
integrality = np.zeros(n + m + n*m)
integrality[:n+m] = 1
# y, x are binary; Axy vary between 0 and each A
lb = np.zeros(n + m + n*m)
ub = np.ones(n + m + n*m)
ub[n+m:] = A.ravel()
# constraints: for each Axy, Axy <= A*y 0 <= A*y - Axy
# Axy <= A*x 0 <= A*x - Axy
Acy = np.repeat(np.eye(n), m, axis=0); Acy[Acy.nonzero()] = A.ravel()
Acx = np.tile(np.eye(m), (n, 1)); Acx[Acx.nonzero()] = A.ravel()
Auy = np.zeros_like(lb); Auy[:n] = 1
Aux = np.zeros_like(lb); Aux[n: n+m] = 1
Ac = np.vstack((
np.hstack((Acy, np.zeros((m*n, m)), -np.eye(m*n))),
np.hstack((np.zeros((m*n, n)), Acx, -np.eye(m*n))),
Auy, Aux,
))
lbc = np.ones(2*m*n + 2); lbc[:-2] = 0
ubc = np.ones_like(lbc); ubc[:-2] = np.inf
result = milp(
c=c, integrality=integrality,
bounds=Bounds(lb=lb, ub=ub),
constraints=LinearConstraint(A=Ac, lb=lbc, ub=ubc),
options={'disp': True},
)
print(result)
print('y =', result.x[:n])
print('x =', result.x[n: n+m])
print('Axy =')
print(result.x[-n*m:].reshape((n, m)))
Running HiGHS 1.2.0 [date: 2021-07-09, git hash: n/a]
Copyright (c) 2022 ERGO-Code under MIT licence terms
Presolving model
26 rows, 19 cols, 55 nonzeros
24 rows, 17 cols, 58 nonzeros
Objective function is integral with scale 1
Solving MIP model with:
24 rows
17 cols (5 binary, 0 integer, 12 implied int., 0 continuous)
58 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% -1292 inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% -374.3333333 -96 289.93% 0 0 0 17 0.0s
L 0 0 0 0.00% -206.6666667 -163 26.79% 19 4 7 39 0.0s
H 0 0 0 0.00% -186.2666667 -166 12.21% 21 6 17 46 0.0s
40.0% inactive integer columns, restarting
Model after restart has 11 rows, 8 cols (3 bin., 0 int., 5 impl., 0 cont.), and 24 nonzeros
0 0 0 0.00% -186.2666667 -166 12.21% 2 0 0 56 0.0s
Solving report
Status Optimal
Primal bound -166
Dual bound -166
Gap 0% (tolerance: 0.01%)
Solution status feasible
-166 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.03 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 62 (total)
0 (strong br.)
29 (separation)
10 (heuristics)
message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: -165.99999999999997
x: [ 0.000e+00 1.000e+00 ... 0.000e+00 -0.000e+00]
mip_node_count: 1
mip_dual_bound: -165.99999999999997
mip_gap: 0.0
y = [ 0. 1. -0.]
x = [ 0. 1. -0. 0.]
Axy =
[[ 0.00000000e+00 -4.73695157e-15 0.00000000e+00 -0.00000000e+00]
[-0.00000000e+00 1.66000000e+02 0.00000000e+00 -0.00000000e+00]
[-0.00000000e+00 -0.00000000e+00 0.00000000e+00 -0.00000000e+00]]
import numpy as np
A = np.array([
( 87, 96, 127, 46),
(155, 166, 92, 11),
(111, 163, 126, 112),
])
n, m = A.shape
i = A.argmax()
y = np.zeros(n, dtype=int)
x = np.zeros(m, dtype=int)
y[i // m] = 1
x[i % m] = 1
print(y)
print(x)
[0 1 0]
[0 1 0 0]
我可以通过求解一个编码x和y值的变量w来线性化这个问题。由于x和y的约束都是只有1个非零元素，y @ A @ x只从矩阵A中选择1个元素，那么我们可以解出A_flat @ w，其中A_flat是平坦矩阵，w的约束是只有一个非零元素。这样，可以使用cvvxpy来解决这个问题。解题中的例子:
A = np.asarray([[87, 96, 127, 46], [155, 166, 92, 11], [111, 163, 126, 112]])
Af = cp.vec(A) #Flattened into vector array
w = cp.Variable(A.size, boolean=True) #encodes x and y boolean array values
objective = cp.Maximize(Af @ w) #Equivalent to y @ A @ x
constraints = [cp.sum(w)==1] #constraint boolean array to only have only 1 non-zero value
prob = cp.Problem(objective,constraints)
prob.solve() #solve using cvxpy
xy = np.reshape(w.value,(A.shape[0],A.shape[1]),order="F")
index = np.argwhere(xy==1)
x = xy[index[0][0],:] #solution for the x boolean array
y = xy[:,index[0][1]] #solution for the y boolean array
那么结果是:
A @ w = 166.0
y @ A @ x = 166.0
新配方等于旧配方，问题解决了。
编辑:更短的方法
与上面类似，但我们可以通过使变量w与矩阵A形状相同并将目标函数的形式更改为cp.sum(cp.multiply(A,w)):
w = cp.Variable(A.shape, boolean=True)
objective = cp.Maximize(cp.sum(cp.multiply(A,w)))
constraints = [cp.sum(w)==1] #constraint boolean array to only have only 1 non-zero value
prob = cp.Problem(objective,constraints)
prob.solve() #solve using cvxpy
index = np.argwhere(w.value==1)
x = w.value[index[0][0],:] #solution for the x boolean array
y = w.value[:,index[0][1]] #solution for the y boolean array