MISDP/MISOCP 在 cvxpy 中
MISDP/MISOCP in cvxpy
我正在尝试解决 CVXPY 中的以下问题。
由于我们正在求解的 PSD 矩阵,问题是混合整数 SDP。然而,根据 this list 看来 none 的求解器可以处理这样的问题。
我们能否利用 A
是 2x2 矩阵这一事实以某种方式将其转换为混合整数 SOCP 问题?
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(271828)
m = 2; n = 50
x = np.random.randn(m,n)
off = cp.Variable(boolean=True)
A = cp.Variable((2,2), PSD=True)
b = cp.Variable(2)
obj = cp.Maximize(cp.log_det(A))
constraints = [ cp.norm(A@x[:,i] + b) <= 1 + 20*off for i in range(n) ]
constraints += [cp.sum(off) <= 20]
prob = cp.Problem(obj, constraints)
optval = prob.solve(solver='XPRESS', verbose=False) # seems to work, although it's not super accurate
print(f"Optimum value: {optval}")
# plot the ellipse and data
angles = np.linspace(0, 2*np.pi, 200)
rhs = np.row_stack((np.cos(angles) - b.value[0], np.sin(angles) - b.value[1]))
ellipse = np.linalg.solve(A.value, rhs)
plt.scatter(x[0,:], x[1,:])
plt.plot(ellipse[0,:].T, ellipse[1,:].T)
plt.xlabel('Dimension 1'); plt.ylabel('Dimension 2')
plt.title('Minimum Volume Ellipsoid')
plt.show()
比方说A=[[x,z], [z,y]]
,那么你可以最大化sqrt(det(A))
(相当于你的objective)。注意
det(A) = xy-z^2
因此最大化 sqrt(det(A))
与最大化 u
相同
xy - z^2 >= u^2
相当于
xy >= z^2 + u^2
这(几乎)是 https://docs.mosek.com/modeling-cookbook/cqo.html#rotated-quadratic-cones
意义上的旋转二阶锥
我觉得像
x >= quad_over_lin([z,u], y)
(没测试语法)可能是用cvxpy最方便的表达方式。
请注意,圆锥体和 quad_over_lin 的定义还强加了 x,y>=0
,因此您不需要单独使用它,并且圆锥约束会自动保证 A
的 PSDness。
我正在尝试解决 CVXPY 中的以下问题。
由于我们正在求解的 PSD 矩阵,问题是混合整数 SDP。然而,根据 this list 看来 none 的求解器可以处理这样的问题。
我们能否利用 A
是 2x2 矩阵这一事实以某种方式将其转换为混合整数 SOCP 问题?
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(271828)
m = 2; n = 50
x = np.random.randn(m,n)
off = cp.Variable(boolean=True)
A = cp.Variable((2,2), PSD=True)
b = cp.Variable(2)
obj = cp.Maximize(cp.log_det(A))
constraints = [ cp.norm(A@x[:,i] + b) <= 1 + 20*off for i in range(n) ]
constraints += [cp.sum(off) <= 20]
prob = cp.Problem(obj, constraints)
optval = prob.solve(solver='XPRESS', verbose=False) # seems to work, although it's not super accurate
print(f"Optimum value: {optval}")
# plot the ellipse and data
angles = np.linspace(0, 2*np.pi, 200)
rhs = np.row_stack((np.cos(angles) - b.value[0], np.sin(angles) - b.value[1]))
ellipse = np.linalg.solve(A.value, rhs)
plt.scatter(x[0,:], x[1,:])
plt.plot(ellipse[0,:].T, ellipse[1,:].T)
plt.xlabel('Dimension 1'); plt.ylabel('Dimension 2')
plt.title('Minimum Volume Ellipsoid')
plt.show()
比方说A=[[x,z], [z,y]]
,那么你可以最大化sqrt(det(A))
(相当于你的objective)。注意
det(A) = xy-z^2
因此最大化 sqrt(det(A))
与最大化 u
相同
xy - z^2 >= u^2
相当于
xy >= z^2 + u^2
这(几乎)是 https://docs.mosek.com/modeling-cookbook/cqo.html#rotated-quadratic-cones
意义上的旋转二阶锥我觉得像
x >= quad_over_lin([z,u], y)
(没测试语法)可能是用cvxpy最方便的表达方式。
请注意,圆锥体和 quad_over_lin 的定义还强加了 x,y>=0
,因此您不需要单独使用它,并且圆锥约束会自动保证 A
的 PSDness。