Overdetermined linear programs

I have a linear program in conic form that I wish to solve in CVX:

cvx_begin;
    variable x(n) nonnegative;
    minimize (c'*x);
    A*x == b;
cvx_end

Here, \boldsymbol{A} \in \mathbb{R}^{m \times n} with m > n. I solved the problem with my own algorithm, and I verified the solution with the linear programming routine in MATLAB’s Optimization Toolbox:

z = linprog(c,[],[],A,b, zeros(size(c)), []);

Thus, I know that the feasible region is nonempty for my choice of \boldsymbol{A}, \boldsymbol{b}, and \boldsymbol{c}. But CVX claims to be unable to solve the problem:

Redundant equality constraints detected.
   CVX cannot solve this problem; but it is likely infeasible.
Status: Overdetermined
Optimal value (cvx_optval): NaN

Is there no way for me to convince CVX that the problem is actually solvable?

A small working example with m = 5 and n = 2 is

A = [-0.9439, -0.1537;  0.2576, -0.4602;  0.6777, 0.5686;  -0.3149, -1.5169;  -0.0058, -1.4912];
b = [-0.1302, -0.2308, 0.3366, -0.8209, -0.7913];
c = [-1.1548, -0.4573]';

I calculate the solution \hat{\boldsymbol{x}} = [0.0516, 0.5305]^T with optimum \boldsymbol{c}^T \hat{\boldsymbol{x}} = -0.3021176 using three methods: my algorithm, MATLAB’s linprog, and a short YALMIP program using the Mosek solver.

EDIT: Mark Stone pointed out that the problem with the aforementioned A, b, and c is infeasible. This was not the case when I originally initialized and ran the problem BUT it was true when I reran my code with A, b, and c generated from the code that I posted earlier.

However, I know that the original problem is feasible because I designed it that way:

m = 5;
n = 2;
A = randn(m,n);
v = rand(n,1);
c = randn(n,1);
b = A*v;

The feasible region must contain at least one point v, which is a least squares solution.

Since I was foolish and did not save the original A, b, and c from my working example, I had to initialize new problem variables to prove my point. However, I saved the matrices to 12 decimal places this time:

A = [4.147002930480000e-01    -5.148816329260000e-01;
     3.484411999520000e-01    -8.964461505020000e-01;
     3.492544166640000e-01    -1.203268186420000e+00;
    -7.292472676300000e-01     1.037815639490000e+00;
     3.268402487630000e-01    -8.459442123360000e-01];
b = [7.288209781380001e-02;
     3.630438960010000e-02;
     2.000855166540000e-02;
    -1.210455116160000e-01;
     3.378115026800000e-02];
c = [-2.971267999950000e-01
     -3.232037795940000e+00

My algorithm recovers \hat{\boldsymbol{x}} = [0.2424865589365387, 0.05375439221311855]^T. I verify it with the YALMIP program

y = sdpvar(n,1);
Constraints = [y >= 0; A*y == b];
Objective = c'*y;
options = sdpsettings('verbose',2,'solver','mosek');
sol = solvesdp(Constraints,Objective,options);
if sol.problem == 0
    solution = double(y);
end

and the Optimization Toolbox code from before:

z = linprog(c,[],[],A,b, zeros(size(c)), []);

Mosek gives me

Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL

and linprog gives me exitflag equal to 1, which I have interpreted (correctly, I hope) as “function converges to a solution”. When I find the difference between my vector x and those results, I get

[z solution] - [x x]

ans =

    -3.372579993055069e-13    -3.372579993055069e-13
    -4.644347406657090e-12    -4.644347406657090e-12

and when I check the feasibility of my solution, I get

A*x - b

ans =

    -2.251421271637355e-12
    -4.045888624126803e-12
    -5.547048931298093e-12
     4.582084711657330e-12
    -3.824808525454415e-12

My own algorithm uses a tolerance of 10^{-6} for constraint satisfaction, so the answer given above is fine for my purposes. But CVX still claims that the problem is infeasible even with cvx_precision set to low! I suspect that CVX demands stricter constraint satisfaction, which I do not need here. Is this true?

While we wait for mcg to weigh in, perhaps you could provide us with A and b (and you might as well include c as well).

A*x=b for   
A = [-0.9439, -0.1537;  0.2576, -0.4602;  0.6777, 0.5686;  -0.3149, -1.5169;  -0.0058, -1.4912];
b = [-0.1302, -0.2308, 0.3366, -0.8209, -0.7913]';

is infeasible, and was reported so when I tried solving it in CVX, YALMIP with several solvers, and with cplexlp (under MATLAB), even without the x >= 0 constraint.

Your solution, x=[0.051590059275147, 0.530452265520773]’ is the least squares solution (with or without the constraint x >= 0) to A*x=b
and can be obtained by several methods in MATLAB, including, ignoring for the moment the constraint x >=0,

A\b
pinv(A)*b

given that the solution satisfies x >= 0, and using CVX or YALMIP, say, to minimize norm(Ax-b) subject to x >= 0. The norm of the residual (norm(Ax-b) using the least squares solution x) is 4.485749216194889e-05, which is not close enough to zero to be considered zero, unless you have a very loose tolerance.

However, solving the LP using only the first 2 rows of A and b, results in a solution close to the least squares solution.

cplexlp(c',[],[],A(1:2,:),b(1:2),[0;0],[])-A\b
ans =
   1.0e-04 *
  -0.176481040301060
  -0.631887765759620

and indeed, it doesn’t matter what c is in this case, as there is only a single feasible solution.

In the interest of completeness, I ran it in linprog.

>> linprog(c',[],[],A,b,zeros(size(c)), [])
Exiting: The primal is infeasible; the equality constraints are dependent
 but not consistent.
ans =
   0.051572411171117
   0.530389076744197

and that last “answer” is the solution to

A(1:2,1:2,)*x=b(1:2)

which is not a feasible solution for A*x=b.

For the record, I would like to see your YALMIP program calling MOSEK, including the output produced. Use sdpsettings(‘solver’,‘mosek’,‘verbose’,2) to get maximum information reported. You claim MOSEK found the solution to the problem, but you misintepreted the output from linprog, so perhaps you misinterpreted MOSEK’s primal infeasibiliy reporting as well? Also, I am curious as to what your algorithm is. Based on my post, your algorithm, or at least your use of the algorithm or interpretation of the results, is erroneous.

Edit: The following applies to the second version of A, b, and c posted by the OP. A*x==b and x >= 0 is feasible. I successfully solved the LP in CVX (note that there are no degrees of freedom in the LP, so the optimal solution is the only feasible solution) when I used SeDuMi as the solver. Using sdtp3 as the solver, it appeared to solve the problem, then I got an error message (using build 1023 on 32 bit MATLAB R2013A under Windows 7 Professional Service Pack 1)

Calling SDPT3 4.0: 3 variables, 0 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints =  1
 dim. of linear var  =  3
 dim. of free   var  =  1 *** convert ublk to lblk
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|4.5e+00|1.4e+01|3.2e+02| 2.962410e+00  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|1.000|1.5e-06|1.0e-01|1.2e+01| 2.808319e+00 -8.044649e+00| 0:0:00| chol  1  1 
 2|1.000|0.989|6.1e-08|1.1e-02|2.8e-01| 1.762058e-01 -8.216665e-02| 0:0:00| chol  1  1 
 3|0.990|0.989|2.0e-08|1.1e-03|3.8e-03| 2.824753e-03 -3.510873e-04| 0:0:00| chol  1  1 
 4|0.989|0.989|3.0e-09|1.1e-04|4.2e-05| 3.124420e-05  5.135530e-05| 0:0:00| chol  1  1 
 5|0.989|0.989|4.6e-11|4.9e-06|4.7e-07| 3.433431e-07  5.643251e-07| 0:0:00| chol  1  1 
 6|0.993|0.989|7.6e-13|5.4e-08|6.8e-09| 4.448080e-09  5.186332e-09| 0:0:00| chol  1  1 
 7|0.994|0.989|6.7e-15|7.8e-10|9.6e-11| 5.514660e-11  4.226974e-11| 0:0:00| chol  1  1 
 8|0.996|0.989|0.0e+00|1.1e-11|1.4e-12| 6.531719e-13  2.601315e-13| 0:0:00| chol  1  1 
 9|1.000|0.981|5.6e-17|1.9e-13|3.1e-14| 8.822476e-15  7.175510e-16| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.82e-12
-------------------------------------------------------------------
 number of iterations   =  9
 primal objective value =  8.82247616e-15
 dual   objective value =  7.17551045e-16
 gap := trace(XZ)       = 3.14e-14
 relative gap           = 3.14e-14
 actual relative gap    = 8.10e-15
 rel. primal infeas     = 5.55e-17
 rel. dual   infeas     = 1.91e-13
 norm(X), norm(y), norm(Z) = 1.0e+00, 7.2e-16, 2.5e-01
 norm(A), norm(b), norm(C) = 2.0e+00, 2.0e+00, 1.2e+00
 Total CPU time (secs)  = 0.19  
 CPU time per iteration = 0.02  
 termination code       =  0
 DIMACS: 5.6e-17  0.0e+00  1.9e-13  0.0e+00  8.1e-15  3.1e-14
-------------------------------------------------------------------
Error using  * 
Inner matrix dimensions must agree.
Error in cvxprob/solve (line 467)
                bval = sgn * ( b' * y + d' );
Error in cvx_end (line 87)
        solve( prob ); 

I then ran CVX again, but this time omitted the redundant constraint x >= 0. The result was:

Trivial infeasibilities detected; solution determined analytically.
Status: Infeasible
Optimal value (cvx_optval): +Inf

So at this point, I will defer to mcg to assess the situation.

MOSEK says

Basic solution summary
Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal.  obj: 0.0000000000e+000   Viol.  con: 8e-014   var: 0e+000 
Dual.    obj: 0.0000000000e+000   Viol.  con: 0e+000   var: 0e+000 

for Ax=b x>=0.

This shows the problem has a solution. Maybe it is a tolerance issue because it might be hard to choose the right tolerances for the linear dependency check.

Avoiding linear dependencies should be avoided if at all possible.

There is an issue with formatting. I tried to fix it but with no luck.

I tried to get solution summary formatted nicely.

I fixed the formatting by using “preformatted text” (i.e., code) menu bar item.

Sorry folks, I’ve been swamped right now (with some CVX developments I am sure people are going to be happy with :-)) If someone would like to package this up as an official bug report to http://support.cvxr.com, I would appreciate it. (More in the next comment…)

What I’d need is a tidy example that reproduces the failure. Probably something could be constructed from the test here, I admit, but I’m more likely to actually get to it if it sits in my support queue.

Bug report received, thanks. I’ll look into this. I’m preparing for a version 2.1 rollout, so if I develop a fix soon it will likely end up in there. (Yes, I know, we never exited beta for 2.0, but a certain new feature of CVX merits a version update, so there you go.)