0

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?

Resu
  • 187
  • 1
  • 12

1 Answers1

1

Internally CVXPY does something else to factorize the matrix, so if a check with eigvalsh reveals you are just borderline good, then the internal method used in CVXPY can tip the balance the other way.

More importantly, there is no need to formulate your problem as a QP if you already have it on conic form to start with. This is a standard SOCP reformulation. Namely,

(1/2)*cp.quad_form(x, P)

is in your case equivalent to

(1/2) * constant * cp.sum_squares(X_LP @ x)

because

x.T@P@x = constant * norm(X_LP@x)^2

This way you avoid the roundabout way where you first construct P=X_LP.T@X_LP and then CVXPY has to work to factorize it back again. You also get a problem convex by construction, so no numerical issues will come in the way of the convexity check. More about this in https://docs.mosek.com/modeling-cookbook/cqo.html#conic-quadratic-modeling

  • Thank you so much, the program runs now. However, I am getting -inf for the objective function with no optimal values. This occurs even if I set the conic part "(1/2) * constant * cp.sum_squares(X_LP @ x)" to zero. This shouldn't be happening since a simple LP from scipy is able to solve it in that case. Do you know why that happens? I will add more to main question to explain what I am getting at. – Resu Jun 05 '23 at 07:06
  • 1
    @Resu You need as a minimum post the log output (use verbose=True in solve()) or fully reproducible data. Otherwise it is impossible to know what kind of numerical issues you are experiencing. – Michal Adamaszek Jun 05 '23 at 07:11
  • ok i will do it and add the log output to the question above. Thanks. – Resu Jun 05 '23 at 07:13
  • Done. You may have a look. – Resu Jun 05 '23 at 07:14
  • 1
    Well, you see yourself. You must have some weird data scaling and the solver internally explodes with objective values like 1e+30. Try different solvers or work on the scaling of your data. Maybe https://docs.mosek.com/modeling-cookbook/practical.html#scaling and the following can help – Michal Adamaszek Jun 05 '23 at 07:33