search

jump-dev / Convex.jl

A Julia package for disciplined convex programming
https://jump.dev/Convex.jl/stable/
Other
564 stars 120 forks source link

SCS fails on 2-variable least-squares problem (on simplex) #108

Closed mpf closed 8 years ago

mpf commented 8 years ago

Below is the output of running Convex on a 2-variable least-squares problem over the simplex, i.e.,

prob = minimize( sum_squares(A*x-b), sum(x)==τ, x≥0)

I'm not familiar enough with the internals of Convex.jl to know if the failure has something to do with the way the problem is reformulated, or if it's a numerical issue with SCS. In either case, the problem seems simple enough that it should be solvable.

Here is the code:

using Convex, SCS, JLD

data = load("/tmp/data.jld")
A, b, τ = data["A"], data["b"], data["τ"]
x = Variable( size(A,2) )
prob = minimize( sum_squares(A*x-b), sum(x)==τ, x≥0)
solve!(prob, SCS.SCSSolver(verbose=1))
@show x

The data can be downloaded here.

julia> include("test/prob.jl")
----------------------------------------------------------------------------
    SCS v1.1.5 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 5386
eps = 1.00e-04, alpha = 1.80, max_iters = 20000, normalize = 1, scale = 5.00
Variables n = 5, constraints m = 2697
Cones:  primal zero / dual free vars: 2
    linear vars: 3
    soc vars: 2692, soc blks: 2
Setup time: 6.88e-04s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0|      inf       inf       nan      -inf       inf       inf  2.46e-04
   100|      inf       inf       nan       inf       inf       inf  8.15e-03
<clip>...</clip>
  3280|      inf       inf       nan       inf       inf       inf  1.96e-01
----------------------------------------------------------------------------
Status: Infeasible
Timing: Total solve time: 1.96e-01s
    Lin-sys: nnz in L factor: 8089, avg solve time: 3.08e-05s
    Cones: avg projection time: 4.05e-06s
----------------------------------------------------------------------------
Certificate of primal infeasibility:
dist(y, K*) = 2.7105e-20
|A'y|_2 * |b|_2 = 9.9951e-05
b'y = -1.0000
============================================================================
WARNING: Problem status Infeasible; solution may be inaccurate.
x = Variable of
size: (2, 1)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()
value: [NaN
 NaN]
Variable of
size: (2, 1)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()
value: 2x1 Array{Float64,2}:
 NaN
 NaN

julia>
madeleineudell commented 8 years ago

That's very strange. You can take a look at the conic form of the problem by calling c,A,b,cones,_ = conic_problem(prob). The problem passed to the solvers is 

minimize c^T x
subject to b - Ax in cones

(see the MPB docs).

I tried this problem with Mosek and ECOS and they also failed. I'm not sure why this particular problem is giving the solvers such trouble. With random data of the same size (and using the same formulation) the solvers succeed.

mlubin commented 8 years ago

The optimal objective value is pretty large (~1e8), but not too bad. Ipopt doesn't have any issues:

using JuMP, Ipopt, JLD

data = load("data.jld")
A, b, τ = data["A"], data["b"], data["τ"]
n = size(A,2)

m = Model(solver=IpoptSolver())
@defVar(m, x[1:n] ≥ 0)
@setObjective(m, Min, sum((A*x-b).^2))
@addConstraint(m, sum(x) == τ)

status = solve(m)
This is Ipopt version 3.12.1, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        2
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:    10752

Total number of variables............................:        2
                     variables with only lower bounds:        2
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.4973601e+08 0.00e+00 1.12e+04   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.3703469e+08 2.84e-14 4.32e-07  -1.4 1.07e+02    -  1.00e+00 1.00e+00f  1
   2  1.3703469e+08 0.00e+00 8.68e-02 -11.0 6.98e-04    -  9.90e-01 1.00e+00f  1
   3  1.3703469e+08 0.00e+00 8.68e-04 -11.0 1.99e-07    -  9.90e-01 1.00e+00h  1
   4  1.3703469e+08 2.84e-14 8.68e-06 -11.0 3.28e-13    -  9.90e-01 1.00e+00f  1
   5  1.3703469e+08 0.00e+00 1.65e-08 -11.0 1.13e-13    -  9.98e-01 1.00e+00h  1
   6  1.3703469e+08 2.84e-14 2.23e-12 -11.0 4.51e-14    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   4.7095083704558457e+05    1.3703469152281412e+08
Dual infeasibility......:   2.2252975357988175e-12    6.4750487179865272e-10
Constraint violation....:   2.8421709430404007e-14    2.8421709430404007e-14
Complementarity.........:   1.3202291009656215e-11    3.8415302269173957e-09
Overall NLP error.......:   1.3202291009656215e-11    3.8415302269173957e-09

Number of objective function evaluations             = 7
Number of objective gradient evaluations             = 7
Number of equality constraint evaluations            = 8
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      0.416
Total CPU secs in NLP function evaluations           =      0.128
mpf commented 8 years ago

Thanks for the input, folks! SCS (via Convex) behaves better if it's given a hint:

# Conic formulation
x = Variable(n)
r = Variable(m)
t = Variable(1)
prob = minimize( t, norm(r) ≤ t, A*x+r==b, sum(x)==τ, x≥0)
solve!(prob, SCS.SCSSolver(verbose=1))

Here's the output:

----------------------------------------------------------------------------
    SCS v1.1.5 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 10761
eps = 1.00e-04, alpha = 1.80, max_iters = 20000, normalize = 1, scale = 5.00
Variables n = 2693, constraints m = 5382
Cones:  primal zero / dual free vars: 2690
    linear vars: 3
    soc vars: 2689, soc blks: 1
Setup time: 2.01e-03s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0|      inf       inf       nan      -inf       inf       inf  5.44e-04
   100|      inf       inf       nan       inf       inf       inf  2.12e-02
   200|      inf       inf       nan       inf       inf       inf  3.80e-02
   300|      inf       inf       nan       inf       inf       inf  5.84e-02
   400|      inf       inf       nan       inf       inf       inf  7.50e-02
   500| 1.22e-02  1.69e-02  1.65e-08  1.17e+04  1.17e+04  0.00e+00  9.32e-02
   600| 9.81e-03  8.99e-03  5.34e-09  1.17e+04  1.17e+04  0.00e+00  1.10e-01
   700| 6.54e-03  7.23e-03  2.95e-09  1.17e+04  1.17e+04  0.00e+00  1.26e-01
   800| 4.44e-03  4.80e-03  1.04e-09  1.17e+04  1.17e+04  0.00e+00  1.43e-01
   900| 3.59e-03  3.09e-03  2.61e-09  1.17e+04  1.17e+04  0.00e+00  1.60e-01
  1000| 2.52e-03  2.45e-03  1.24e-09  1.17e+04  1.17e+04  0.00e+00  1.75e-01
  1100| 1.69e-03  1.67e-03  1.06e-09  1.17e+04  1.17e+04  0.00e+00  1.90e-01
  1200| 1.49e-03  1.18e-03  3.23e-10  1.17e+04  1.17e+04  0.00e+00  2.07e-01
  1300| 1.04e-03  9.57e-04  4.27e-10  1.17e+04  1.17e+04  0.00e+00  2.22e-01
  1400| 6.28e-04  7.16e-04  7.82e-11  1.17e+04  1.17e+04  0.00e+00  2.38e-01
  1500| 8.12e-04  4.52e-04  3.11e-10  1.17e+04  1.17e+04  0.00e+00  2.56e-01
  1600| 3.83e-04  4.25e-04  7.66e-11  1.17e+04  1.17e+04  0.00e+00  2.71e-01
  1700| 3.94e-04  3.07e-04  1.33e-10  1.17e+04  1.17e+04  0.00e+00  2.86e-01
  1800| 4.28e-04  2.02e-04  4.11e-11  1.17e+04  1.17e+04  0.00e+00  3.05e-01
  1900| 1.67e-04  1.96e-04  7.63e-11  1.17e+04  1.17e+04  0.00e+00  3.21e-01
  2000| 2.68e-04  1.27e-04  2.96e-11  1.17e+04  1.17e+04  0.00e+00  3.38e-01
  2100| 2.04e-04  1.08e-04  2.65e-11  1.17e+04  1.17e+04  0.00e+00  3.57e-01
  2200| 1.00e-04  9.33e-05  3.22e-12  1.17e+04  1.17e+04  0.00e+00  3.77e-01
  2240| 8.15e-05  8.79e-05  1.91e-12  1.17e+04  1.17e+04  0.00e+00  3.83e-01
----------------------------------------------------------------------------
Status: Solved
Timing: Total solve time: 3.83e-01s
    Lin-sys: nnz in L factor: 24213, avg solve time: 1.01e-04s
    Cones: avg projection time: 3.89e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 1.3217e-14, dist(y, K*) = 0.0000e+00, s'y/m = 3.0304e-15
|Ax + s - b|_2 / (1 + |b|_2) = 8.1460e-05
|A'y + c|_2 / (1 + |c|_2) = 8.7927e-05
|c'x + b'y| / (1 + |c'x| + |b'y|) = 1.9135e-12
----------------------------------------------------------------------------
c'x = 11707.1375, -b'y = 11707.1375
============================================================================