I have an optimization problem that is convex and bounded: \begin{align*} \text{Minimize}_{\text{wrt} B_{\text{opt}}}\qquad &\frac{1}{2p\sigma^2}\|Y_{\text{opt}}-X_{\text{opt}}B_{\text{opt}}\|_2^2-\sum_{l=1}^L(\log\pi_l)1_p^\top\{I_{pL}\}_{p(l-1):pl,:}B_{\text{opt}} \\ \text{subject to}\qquad & C_{\text{opt}}B_{\text{opt}}=1_p, \end{align*} where:
- $p$ and $\sigma$ are fixed constants;
- $\sum_l\pi_l=1$ are fixed probabilities;
- $Y_{\text{opt}}\in\mathbb{R}^{nL}$, $X_{\text{opt}}\in\{0,1\}^{nL\times pL}$, and $B_{\text{opt}}\in\{0,1\}^{pL}$, where $n$, $p$, $L$ are all fixed constants;
- $1_p$ is just a $p$-length vector of ones;
- $\{I_{pL}\}_{p(l-1):pl,:}$ is just the sub-matrix obtained by taking the $p(l-1)$ to $pl$ rows of the identity matrix $I_{pL}$;
- $C_{\text{opt}}=[I_p,\dots,I_p]\in\{0,1\}^{p\times pL}$, where $I_p$ is the $p\times p$ identity matrix.
This is convex and the solution can analytically be derived if there was no constraint. However, when I run CVXPY for it, I get an unbounded objective function. In fact CVXPY gives an unbounded objective for this optimization even when there is no constraint (what is going on?)! I'm quite confused as to why it is not working...
I have added a sample code below, along with the verbose outputs of CVXPY.
def run_QP(n, p, L, Y, X, sigma, B_prop):
# Configuring our inputs to suit LP
Y_opt = Y.flatten('F')
X_opt = np.zeros((n*L,p*L))
for l in range(L):
X_opt[n*l:n*(l+1),p*l:p*(l+1)] = X
C_opt = np.eye(p)
for l in range(L-1):
C_opt = np.concatenate((C_opt, np.eye(p)), axis=1)
one_p = np.ones(p)
# Setting the objective matrix and vector
q = np.zeros(p*L)
for l in range(L):
I_pL = np.eye(p*L)
I_pL_trun = I_pL[p*l:p*(l+1),:]
q -= np.log(B_prop[l]) * np.dot(one_p.T,I_pL_trun)
# Setting the equality constraints matrix
A = C_opt
# Setting the equality constraints vector
b = one_p
# Define and solve the CVXPY problem
constant = 1/(2*p*(sigma**2))
B_opt = cp.Variable(p*L)
objective = cp.Minimize(constant*cp.sum_squares(Y_opt - X_opt @ B_opt) + (q.T @ B_opt))
constraints = []
constraints.append(A @ B_opt == b)
problem = cp.Problem(objective, constraints)
problem.solve(solver=cp.OSQP,verbose=True)
print('optimal obj value:', problem.value)
return
Running the above function:
B_bar_prob = [0.5,0.5]
alpha = 0.5
p = 100
n = 50
sigma = 0.1
L = len(B_bar_prob)
indices = np.random.choice(np.array([0, 1]), size=p, p=B_bar_prob)
B = np.eye(np.max(indices)+1)[indices]
X = binomial(1, alpha, (n, p))
Psi = normal(0, sigma, (n, L))
Y = np.dot(X, B) + Psi
B_prop = np.sum(B, axis=0) / p
run_QP(n, p, L, Y, X, sigma, B_prop)
and the output is
===============================================================================
CVXPY
v1.3.1
===============================================================================
(CVXPY) Jun 06 03:18:27 AM: Your problem has 200 variables, 1 constraints, and 0 parameters.
(CVXPY) Jun 06 03:18:27 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Jun 06 03:18:27 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Jun 06 03:18:27 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Jun 06 03:18:27 AM: Compiling problem (target solver=OSQP).
(CVXPY) Jun 06 03:18:27 AM: Reduction chain: CvxAttr2Constr -> Qp2SymbolicQp -> QpMatrixStuffing -> OSQP
(CVXPY) Jun 06 03:18:27 AM: Applying reduction CvxAttr2Constr
(CVXPY) Jun 06 03:18:27 AM: Applying reduction Qp2SymbolicQp
(CVXPY) Jun 06 03:18:27 AM: Applying reduction QpMatrixStuffing
(CVXPY) Jun 06 03:18:27 AM: Applying reduction OSQP
(CVXPY) Jun 06 03:18:27 AM: Finished problem compilation (took 6.337e-02 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Jun 06 03:18:27 AM: Invoking solver OSQP to obtain a solution.
-----------------------------------------------------------------
OSQP v0.6.2 - Operator Splitting QP Solver
(c) Bartolomeo Stellato, Goran Banjac
University of Oxford - Stanford University 2021
-----------------------------------------------------------------
problem: variables n = 300, constraints m = 200
nnz(P) + nnz(A) = 5400
settings: linear system solver = qdldl,
eps_abs = 1.0e-05, eps_rel = 1.0e-05,
eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
rho = 1.00e-01 (adaptive),
sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: on, time_limit: off
iter objective pri res dua res rho time
1 -1.7228e+05 3.90e+01 9.80e+04 1.00e-01 9.51e-03s
25 -1.0000e+30 1.20e-02 4.28e-01 1.00e-01 1.22e-02s
status: dual infeasible
number of iterations: 25
run time: 1.24e-02s
optimal rho estimate: 3.00e-03
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Jun 06 03:18:27 AM: Problem status: unbounded
(CVXPY) Jun 06 03:18:27 AM: Optimal value: -inf
(CVXPY) Jun 06 03:18:27 AM: Compilation took 6.337e-02 seconds
(CVXPY) Jun 06 03:18:27 AM: Solver (including time spent in interface) took 2.242e-02 seconds
optimal obj value: -inf