I have a quadratic program (see below) that causes an error due to my quad_form matrix being not symmetric or PSD. Why is this happening? (Note that the issue is with my matrix P below so feel free to ignore the first block of code in the function.)
def run_QP(n, p, L, Y, X, sigma, B_prop):
# Configuring our inputs to suit LP
Y_LP = Y.flatten('F')
X_LP = np.zeros((n*L,p*L))
for l in range(L):
X_LP[n*l:n*(l+1),p*l:p*(l+1)] = X
C_LP = np.eye(p)
for l in range(L-1):
C_LP = np.concatenate((C_LP, np.eye(p)), axis=1)
one_p = np.ones(p)
# Setting the objective matrix and vector
constant = 1/(p*cp.power(sigma,2))
P = constant * (X_LP.T @ X_LP)
I_pL = np.eye(p*L)
temp_term = -1*cp.log(B_prop[0]) * (one_p.T @ I_pL[0:p,:])
for l in range(1,L):
I_pL_trun = I_pL[p*l:p*(l+1),:]
temp_term -= cp.log(B_prop[l]) * (one_p.T @ I_pL_trun)
q = -1*constant*(Y_LP.T@X_LP) + temp_term
# Setting the equality constraints matrix
A = np.concatenate((X_LP, C_LP), axis=0)
# Setting the equality constraints vector
b = np.concatenate((Y_LP, one_p), axis=0)
# Define and solve the CVXPY problem.
x = cp.Variable(p*L)
prob = cp.Problem(cp.Minimize((1/2)*cp.quad_form(x, P) + q.T @ x),
[A @ x == b,])
prob.solve()
# Reconfiguring our outputs to suit pooled data
B_QP_est = prob.value
B_est = np.zeros((p,L))
for l in range(L):
B_col = B_QP_est[p*l:p*(l+1)]
B_est[:,l] = B_col
return B_est
B_bar_prob = [0.5,0.5]
L = 2
alpha = 0.5
p = 100
n = 20
sigma = 0.1
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) / n
B_hat = run_QP(n, p, L, Y, X, sigma, B_prop)
When I check for PSD using the following code (after swapping my cvxpy atomic functions out for numpy functions, I get true):
def isPSD(A, tol=1e-8):
E = np.linalg.eigvalsh(A)
return np.all(E > -tol)
Hence, I am very confused as to what is going on...
Edit: Thanks to the answer below from Michal Adamaszek, I have modified my program to work.
def run_QP(n, p, L, Y, X, sigma, B_prop):
# Configuring our inputs to suit LP
Y_LP = Y.flatten('F')
X_LP = np.zeros((n*L,p*L))
for l in range(L):
X_LP[n*l:n*(l+1),p*l:p*(l+1)] = X
C_LP = np.eye(p)
for l in range(L-1):
C_LP = np.concatenate((C_LP, np.eye(p)), axis=1)
one_p = np.ones(p)
# Setting the objective matrix and vector
constant = 1/(p*(sigma**2))
temp_term = 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),:]
temp_term -= np.log(B_prop[l]) * np.dot(one_p.T,I_pL_trun)
q = -1*constant*(Y_LP.T @ X_LP) + temp_term
# Setting the equality constraints matrix
A = np.concatenate((X_LP, C_LP), axis=0)
# Setting the equality constraints vector
b = np.concatenate((Y_LP, one_p), axis=0)
# Define and solve the CVXPY problem.
x = cp.Variable(p*L)
objective = cp.Minimize((1/2)*constant*cp.sum_squares(X_LP @ x) + (q.T @ x))
constraints = []
constraints.append(A @ x == b)
problem = cp.Problem(objective, constraints)
problem.solve()
print('optimal obj value:', problem.value)
# Reconfiguring our outputs to suit pooled data
B_QP_est = x.value
B_est = np.zeros((p,L))
for l in range(L):
B_col = B_QP_est[p*l:p*(l+1)]
B_est[:,l] = B_col
return B_est
However, this always gives inf
for the objective function, even when the conic part is zero. See the following output:
===============================================================================
CVXPY
v1.3.1
===============================================================================
(CVXPY) Jun 05 07:13:28 AM: Your problem has 200 variables, 1 constraints, and 0 parameters.
(CVXPY) Jun 05 07:13:28 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Jun 05 07:13:28 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Jun 05 07:13:28 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Jun 05 07:13:28 AM: Compiling problem (target solver=OSQP).
(CVXPY) Jun 05 07:13:28 AM: Reduction chain: CvxAttr2Constr -> Qp2SymbolicQp -> QpMatrixStuffing -> OSQP
(CVXPY) Jun 05 07:13:28 AM: Applying reduction CvxAttr2Constr
(CVXPY) Jun 05 07:13:28 AM: Applying reduction Qp2SymbolicQp
(CVXPY) Jun 05 07:13:28 AM: Applying reduction QpMatrixStuffing
(CVXPY) Jun 05 07:13:28 AM: Applying reduction OSQP
(CVXPY) Jun 05 07:13:28 AM: Finished problem compilation (took 8.957e-02 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Jun 05 07:13:28 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 = 240, constraints m = 180
nnz(P) + nnz(A) = 4156
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 -9.5106e-01 3.09e+01 1.59e+07 1.00e-01 6.26e-03s
50 1.0000e+30 2.49e-02 2.14e-04 1.00e-01 8.22e-03s
status: primal infeasible
number of iterations: 50
run time: 8.35e-03s
optimal rho estimate: 3.87e+00
So when the conic part is zero, a simple linear program below gets me the optimal solution:
from scipy.optimize import linprog
def run_LP(n, p, L, Y, X, B_prop):
# Configuring our inputs to suit LP
Y_LP = Y.flatten('F')
X_LP = np.zeros((n*L,p*L))
for l in range(L):
X_LP[n*l:n*(l+1),p*l:p*(l+1)] = X
C_LP = np.eye(p)
for l in range(L-1):
C_LP = np.concatenate((C_LP, np.eye(p)), axis=1)
one_p = np.ones(p)
# Setting the objective vector
c = 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),:]
c -= np.log(B_prop[l]) * np.dot(one_p.T,I_pL_trun)
# Setting the equality constraints matrix
A = np.concatenate((X_LP, C_LP), axis=0)
# Setting the equality constraints vector
b = np.concatenate((Y_LP, one_p), axis=0)
# Solving linear programming problem
res = linprog(c, A_eq=A, b_eq=b)
print("Linear program:", res.message)
# Reconfiguring our outputs to suit pooled data
B_LP_est = res.x
B_est = np.zeros((p,L))
for l in range(L):
B_col = B_LP_est[p*l:p*(l+1)]
B_est[:,l] = B_col
return B_est
Why is this happening?