MISDP/MISOCP em cvxpy

0

Pergunta

Eu estou tentando resolver o seguinte problema no CVXPY.

O problema é um mixed-integer SDP devido ao PSD matriz estamos resolvendo. No entanto, de acordo com essa lista , parece que nenhum dos solucionadores pode lidar com tal problema.

Podemos usar o fato de que A é uma matriz 2x2, de alguma forma, converter para um mixed-integer SOCP problema?

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()
1

Melhor resposta

1

Vamos dizer que A=[[x,z], [z,y]]e , em seguida, você pode maximizar sqrt(det(A)) (o que é equivalente ao seu objetivo). Note que

det(A) = xy-z^2

para maximizar sqrt(det(A)) é o mesmo que maximizar u sujeito a

xy - z^2 >= u^2

equivalentemente

xy >= z^2 + u^2

Isso é (quase) uma rodada de segunda ordem cone no sentido de https://docs.mosek.com/modeling-cookbook/cqo.html#rotated-quadratic-cones

Eu acho que algo como

x >= quad_over_lin([z,u], y)

(ainda não testado a sintaxe) pode ser a maneira mais conveniente de se expressar em cvxpy.

Note que a definição de cone e de quad_over_lin também impõe x,y>=0 assim você não precisa que, separadamente, e a cônica restrições garante automaticamente PSDness de A.

2021-11-24 09:08:37

E aqui está um pouco de luz de entretenimento relacionados com o seu problema: twitter.com/mosektw/status/1190259248250445824
Michal Adamaszek

Em outros idiomas

Esta página está em outros idiomas

Русский
..................................................................................................................
Italiano
..................................................................................................................
Polski
..................................................................................................................
Română
..................................................................................................................
한국어
..................................................................................................................
हिन्दी
..................................................................................................................
Français
..................................................................................................................
Türk
..................................................................................................................
Česk
..................................................................................................................
ไทย
..................................................................................................................
中文
..................................................................................................................
Español
..................................................................................................................
Slovenský
..................................................................................................................