search

cvxgrp / scs

Splitting Conic Solver
MIT License
553 stars 136 forks source link

Solution not found on sparse quadratic program #234

Closed stephane-caron closed 1 year ago

stephane-caron commented 2 years ago

Specifications

Description

SCS does not find a solution to the attached problem: scs_dump.zip.

This problem should be feasible, with an objective value at the optimum around 562.443e9.

How to reproduce

The problem was generated by the following script:

import numpy as np
import scipy.sparse as spa
from qpsolvers import solve_ls

def generate_problem(n: int):
    # minimize || x - s ||^2
    R = spa.eye(n, format="csc")
    s = np.array(range(n), dtype=float)
    # such that G * x <= h
    G = spa.diags(
        diagonals=[
            [1.0 if i % 2 == 0 else 0.0 for i in range(n)],
            [1.0 if i % 3 == 0 else 0.0 for i in range(n - 1)],
            [1.0 if i % 5 == 0 else 0.0 for i in range(n - 1)],
        ],
        offsets=[0, 1, -1],
        format="csc",
    )
    h = np.ones(G.shape[0])
    # such that sum(x) == 42
    A = spa.csc_matrix(np.ones((1, n)))
    b = np.array([42.0]).reshape((1,))
    # such that x >= 0
    lb = np.zeros(n)
    return R, s, G, h, A, b, lb

if __name__ == "__main__":
    R, s, G, h, A, b, lb = generate_problem(n=15000)
    x = solve_ls(R, s, G, h, A, b, lb, solver="scs", verbose=False, write_data_filename="scs_dump")
    assert x is not None

Additional information

I tried varying the problem size n and comparing to other solvers:

Solver n Objective at returned solution
CVXOPT 1,500 561874994.0
HiGHS 1,500 561874994.0
OSQP 1,500 561875870.1
ProxQP 1,500 561874993.6
SCS 1,500 561874991.4
CVXOPT 15,000 562443121619.2
HiGHS 15,000 562443121619.0
OSQP 15,000 562443434746.2
ProxQP 15,000 562443121598.1
SCS 15,000 not found (this issue)

(All solvers use their default settings. The objective is $1/2 \Vert x - s \Vert^2$ at the returned solution $x$.)

The conversion from (R, s, G, h, A, b, lb) to SCS cone, data works for smaller values of n, but this doesn't mean it is clear of all suspicion. Hoping scs_dump.zip helps investigate.

Output

------------------------------------------------------------------
           SCS v3.2.2 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 15000, constraints m: 30002
cones:    z: primal zero / dual free vars: 1
      l: linear vars: 15000
      b: box cone vars: 15001
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
      alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
      max_iters: 100000, normalize: 1, rho_x: 1.00e-06
      acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
      nnz(A): 45500, nnz(P): 15000
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.52e+04  1.39e+04  2.97e+11 -2.30e+11  1.00e-01  5.38e-02 
------------------------------------------------------------------
status:  unbounded
timings: total: 5.42e-02s = setup: 4.87e-02s + solve: 5.42e-03s
     lin-sys: 5.40e-04s, cones: 4.99e-04s, accel: 4.40e-08s
------------------------------------------------------------------
objective = -inf
------------------------------------------------------------------
bodono commented 2 years ago

I think this is a case of the infeasible tolerances being too high, I reduced it to eps_infeas = 1.0e-12 and was able to solve the problem:

------------------------------------------------------------------
           SCS v3.2.2 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 15000, constraints m: 30002
cones:    z: primal zero / dual free vars: 1
      l: linear vars: 15000
      b: box cone vars: 15001
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-12
      alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
      max_iters: 100000, normalize: 1, rho_x: 1.00e-06
      acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
      nnz(A): 45500, nnz(P): 15000
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.52e+04  1.39e+04  2.97e+11 -2.30e+11  1.00e-01  2.87e-02
   250| 1.62e+00  2.34e+01  7.76e+05 -3.49e+07  1.02e+01  2.14e-01
   500| 1.29e+00  6.17e+01  1.01e+06 -1.57e+07  3.28e+01  3.94e-01
   750| 7.37e-01  3.70e+01  2.53e+04 -4.32e+06  3.28e+01  5.73e-01
  1000| 4.51e-01  2.02e+01  5.46e+04 -1.66e+06  3.28e+01  7.52e-01
  1250| 4.10e-01  1.16e+01  2.05e+04 -7.93e+05  3.28e+01  9.32e-01
  1500| 3.57e-01  1.87e+01  8.23e+04 -6.24e+05  1.05e+02  1.11e+00
  1750| 8.13e-02  6.47e+01  1.32e+05 -6.95e+05  1.19e+03  1.29e+00
  2000| 4.26e+00  5.08e+03  5.75e+03 -6.33e+05  1.19e+03  1.47e+00
  2250| 1.88e-04  9.20e+00  7.79e+03 -6.26e+05  1.19e+03  1.65e+00
  2500| 9.56e-05  1.41e+00  5.77e+02 -6.29e+05  9.23e+01  1.83e+00
  2600| 1.15e-04  1.07e-02  3.04e-01 -6.30e+05  9.23e+01  1.91e+00
------------------------------------------------------------------
status:  solved
timings: total: 1.91e+00s = setup: 2.65e-02s + solve: 1.88e+00s
     lin-sys: 1.36e+00s, cones: 2.27e-01s, accel: 4.21e-02s
------------------------------------------------------------------
objective = -629630.844620
------------------------------------------------------------------

I also tested decrease the solution tolerances to 1e-10 to get a more accurate solution:

------------------------------------------------------------------
           SCS v3.2.2 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 15000, constraints m: 30002
cones:    z: primal zero / dual free vars: 1
      l: linear vars: 15000
      b: box cone vars: 15001
settings: eps_abs: 1.0e-10, eps_rel: 1.0e-10, eps_infeas: 1.0e-12
      alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
      max_iters: 100000, normalize: 1, rho_x: 1.00e-06
      acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
      nnz(A): 45500, nnz(P): 15000
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.52e+04  1.39e+04  2.97e+11 -2.30e+11  1.00e-01  2.95e-02
   250| 1.62e+00  2.34e+01  7.76e+05 -3.49e+07  1.02e+01  2.16e-01
   500| 1.29e+00  6.17e+01  1.01e+06 -1.57e+07  3.28e+01  3.97e-01
   750| 7.37e-01  3.70e+01  2.53e+04 -4.32e+06  3.28e+01  5.78e-01
  1000| 4.51e-01  2.02e+01  5.46e+04 -1.66e+06  3.28e+01  7.60e-01
  1250| 4.10e-01  1.16e+01  2.05e+04 -7.93e+05  3.28e+01  9.41e-01
  1500| 3.57e-01  1.87e+01  8.23e+04 -6.24e+05  1.05e+02  1.12e+00
  1750| 8.13e-02  6.47e+01  1.32e+05 -6.95e+05  1.19e+03  1.30e+00
  2000| 4.26e+00  5.08e+03  5.75e+03 -6.33e+05  1.19e+03  1.48e+00
  2250| 1.88e-04  9.20e+00  7.79e+03 -6.26e+05  1.19e+03  1.66e+00
  2500| 9.56e-05  1.41e+00  5.77e+02 -6.29e+05  9.23e+01  1.84e+00
  2750| 1.70e-08  2.18e-06  9.26e-04 -6.30e+05  9.23e+01  2.02e+00
  3000| 2.36e-08  6.67e-10  2.59e-04 -6.30e+05  9.23e+01  2.20e+00
  3250| 8.95e-09  7.62e-07  1.77e-04 -6.30e+05  5.68e+03  2.38e+00
  3375| 2.01e-09  5.86e-07  5.98e-06 -6.30e+05  2.55e+04  2.47e+00
------------------------------------------------------------------
status:  solved
timings: total: 2.47e+00s = setup: 2.70e-02s + solve: 2.44e+00s
     lin-sys: 1.77e+00s, cones: 2.80e-01s, accel: 6.05e-02s
------------------------------------------------------------------
objective = -629631.000025
------------------------------------------------------------------
stephane-caron commented 1 year ago

Confirmed, I added that setting to the continuous-integration test :ok_hand: Thank you for taking a look :smiley: