jump-dev / Convex.jl

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

Unbounded \ell_1-regularized least squares problem #181

Closed nan2ge1 closed 7 years ago

nan2ge1 commented 7 years ago
julia> using Convex

julia> m, n    = 50, 100;

julia> A       = complex( randn(m, n), randn(m, n) );

julia> b       = complex( randn(m), randn(m) );

julia> λ       = 1.;

julia> x       = ComplexVariable(n);

julia> problem = minimize( 0.5 * sumsquares(A*x - b) + λ * norm(x, 1) );

julia> solve!(problem)
WARN: m less than n, problem likely degenerate
----------------------------------------------------------------------------
        SCS v1.2.6 - Splitting Conic Solver
        (c) Brendan O'Donoghue, Stanford University, 2012-2016
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 20111
eps = 1.00e-04, alpha = 1.80, max_iters = 20000, normalize = 1, scale = 5.00
Variables n = 304, constraints m = 110
Cones:  primal zero / dual free vars: 1
        linear vars: 2
        soc vars: 107, soc blks: 3
Setup time: 2.21e-03s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0|      inf       inf      -nan      -inf       inf       inf  2.69e-04 
----------------------------------------------------------------------------
Status: Unbounded
Timing: Solve time: 2.75e-04s
        Lin-sys: nnz in L factor: 25475, avg solve time: 1.70e-04s
        Cones: avg projection time: 8.51e-07s
----------------------------------------------------------------------------
Certificate of dual infeasibility:
dist(s, K) = 2.0680e-25
|Ax + s|_2 * |c|_2 = 8.8319e-05
c'x = -1.0000
============================================================================
WARNING: Problem status Unbounded; solution may be inaccurate.

Since the original problem is feasible as long as λ > 0, I wonder what went wrong here? By the way, changing \ell_1 to \ell_q norm where q > 1 does not seem to share this issue.

julia> Pkg.status()
2 required packages:
 - Convex                        0.4.0+             master
 - SCS                           0.3.1
madeleineudell commented 7 years ago

Interesting! This is almost certainly due to a problem in the way we represent the l1 norm of complex variables. @Ayush-iitkgp , any ideas?

nan2ge1 commented 7 years ago

@madeleineudell Thank you very much for this analysis. For any x \in C^n, I would have just defined ||x||_1 := ||[x_re x_im]||_{2,1}, i.e.,

sum(sumabs2(hcat( real(x), imag(x) ), 2).^0.5),

if there is no other elegant way to do so.

Ayush-iitkgp commented 7 years ago

@nan2ge1 thanks for reporting this issue. @madeleineudell I think it has to do with the conic_form definition of abs atom as norm1 is defined as norm(x,1) = sum(abs(x))

https://github.com/JuliaOpt/Convex.jl/blob/master/src/atoms/lp_cone/abs.jl#L49 is where we have considered the special case of complex variables. Is the error something to do with it?

madeleineudell commented 7 years ago

I took a brief look at the abs atom, and it seem to do exactly the right thing: norm([real(x), complex(x)], 2). I'm afraid the way to debug this will be to look at the conic form of the expression itself, for a suitably small example. I agree with Nan that it's not clear why the reformulation would result in the solver declaring infeasibility.

On Fri, Mar 10, 2017 at 8:56 AM, Ayush Pandey notifications@github.com wrote:

@nan2ge1 https://github.com/nan2ge1 thanks for reporting this issue. @madeleineudell https://github.com/madeleineudell I think it has to do with the conic_form definition of abs atom as norm1 is defined as norm(x,1) = sum(abs(x))

https://github.com/JuliaOpt/Convex.jl/blob/master/src/ atoms/lp_cone/abs.jl#L49 is where we have considered the special case of complex variables. Is the error something to do with it?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaOpt/Convex.jl/issues/181#issuecomment-285674794, or mute the thread https://github.com/notifications/unsubscribe-auth/AAyp9DLr6RjxokDBGjly6cVHEfSUCbgIks5rkVaZgaJpZM4MXXG_ .

-- Madeleine Udell Assistant Professor, Operations Research and Information Engineering Cornell University https://people.orie.cornell.edu/mru8/ (415) 729-4115

Ayush-iitkgp commented 7 years ago

Thanks Nan. I must accept there was a slight mistake in the code. PR #182 fixes it. Let us know if the problem persist. Also, let us now about where you are using Convex.jl in optimization with complex variables, we are excited to know about it and spread the word. We also encourage you to make PR relating to the various applications in complex-domain.

nan2ge1 commented 7 years ago

@madeleineudell @Ayush-iitkgp Thank you so much for your efforts! I have tried the new version with several circular Gaussian random matrices of different sizes and it always converged to, up to some rounding errors, the same solution as e.g. ADMM. I will see how I could contribute to this project.