Incorrect solution for LP

[blegat]

Consider the problem

min -x2
- x1 - 2x2 = 0
  x1 - 3x2 = 0
x1 free
x2 >= 0

with dual

max 0
-y1 +  y2 = 0
2y1 + 3y2 >= 1

The only feasible solution of the primal is (0, 0) and for the dual, the feasible solution is the ray (1, 1) starting at (1/5, 1/5).
SeDuMi returns the correct x (0, 0) but an incorrect ỳ (-1.2, 0.6).

The problem is the 13th test problem of MOI and it is used to verify the correctness of the Julia SeDuMi wrapper.
Here is what I get when using sedumi directly:

julia> A = [-1.0 -2.0; 1.0 -3.0]
2×2 Array{Float64,2}:
 -1.0  -2.0
  1.0  -3.0

julia> b = [0.0, 0.0]
2-element Array{Float64,1}:
 0.0
 0.0

julia> c = [0.0; -1.0]
2-element Array{Float64,1}:
  0.0
 -1.0

julia> using SeDuMi

julia> sedumi(A, b, c, SeDuMi.Cone(1, 1))
>> SeDuMi 1.32 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Put 1 free variables in a quadratic cone
eqs m = 2, order n = 4, dim = 4, blocks = 2
nnz(A) = 4 + 0, nnz(ADA) = 4, nnz(L) = 3
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            1.25E+01 0.000
  1 :   0.00E+00 1.44E-06 0.000 0.0000 1.0000 1.0000   1.00  1  1  1.0E-06
  2 :   0.00E+00 5.75E-10 0.000 0.0004 0.9999 0.9999   1.00  1  2  3.5E-10

iter seconds digits       c*x               b*y
  2      0.4   Inf -3.2329888893e-11  0.0000000000e+00
|Ax-b| =   9.8e-10, [Ay-c]_+ =   1.6E-16, |x|=  1.4e+00, |y|=  1.3e+00

Detailed timing (sec)
   Pre          IPM          Post
3.601E-01    2.107E-01    2.660E-02    
Max-norms: ||b||=0, ||c|| = 1,
Cholesky |add|=0, |skip| = 1, ||L.L|| = 2.
([-4.67511e-10, 3.23299e-11], [-1.2, 0.6], Dict{String,Any}("timing"=>[0.360066 0.210701 0.0266042],"dinf"=>0.0,"pinf"=>0.0,"cpusec"=>1.66,"wallsec"=>0.597372,"feasratio"=>1.0,"r0"=>9.75092e-10,"iter"=>2.0,"numerr"=>0.0))

Any idea what's the issue ?

[siko1056]

Does your input data match your theoretical computation? Your first constraint is - x1 - x2 = 0 whereas your first data row is A = [-1.0 -2.0; and not A = [-1.0 -1.0;

When looking at the 13th test problem of MOI there your dual version is given to check feasibility. Thus what is the purpose of your derived primal LP?

Maybe you not dualize it. Just reformulate the given problem for SeDuMi by introducing a slack variable for the inequality constraint:

min              x3
s.t. -x1 +  x2       = 0
     2x1 + 3x2 + x3  = 1
                 x3 >= 0

For Octave 4.4.1 and SeDuMi, I get a result much closer to your desired values (note that x and y are swapped compared to your example):

c = [0 0 1];
b = [0 1];
A = [1 -1 0; 2 3 1];
K.f = 2;
K.l = 1;

[x, y] = sedumi (A, b, c, K)

SeDuMi 1.32 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Put 2 free variables in a quadratic cone
eqs m = 2, order n = 4, dim = 5, blocks = 2
nnz(A) = 5 + 0, nnz(ADA) = 4, nnz(L) = 3
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            6.04E+00 0.000
  1 :   3.81E-02 4.27E-01 0.000 0.0707 0.9900 0.9900   1.40  1  1  9.0E-01
  2 :  -9.07E-06 1.68E-04 0.000 0.0004 0.9999 0.9999   1.03  1  1  4.8E-02
  3 :  -1.02E-12 1.90E-11 0.000 0.0000 1.0000 1.0000   1.00  1  1  5.4E-09

iter seconds digits       c*x               b*y
  3      0.1   6.3  5.9391019751e-12 -1.0242923753e-12
|Ax-b| =   2.4e-12, [Ay-c]_+ =   2.6E-12, |x|=  2.3e+00, |y|=  1.1e-12

Detailed timing (sec)
   Pre          IPM          Post
2.242E-02    7.471E-02    5.411E-03
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.
x =
   2.0000e-01
   2.0000e-01
   5.9391e-12

y =
  -5.1214e-13
  -1.0243e-12

[blegat]

Does your input data match your theoretical computation? Your first constraint is - x1 - x2 = 0 whereas your first data row is A = [-1.0 -2.0; and not A = [-1.0 -1.0;

Indeed, sorry for the mistake, I have corrected it in the first post.

When looking at the 13th test problem of MOI there your dual version is given to check feasibility. Thus what is the purpose of your derived primal LP?

I have decided to always use the dual form:
https://github.com/blegat/SeDuMi.jl/blob/f1df5c28a9ab0138184d3fcb7d8ece1806837c46/src/MOI_wrapper.jl#L11-L20
Maybe an option could be added to use the primal form (or automatically choose one based on some heuristic).

Maybe you not dualize it. Just reformulate the given problem for SeDuMi by introducing a slack variable for the inequality constraint:

For Octave 4.4.1 and SeDuMi, I get a result much closer to your desired values (note that x and y are swapped compared to your example):

Indeed, thanks for the tip.
However, if the user entered the dual of the 13th MOI test problem to SeDuMi, e.g.,

min -x2
- x1 - 2x2 = 0
  x1 - 3x2 = 0
x1 free
x2 >= 0

he would still get an incorrect dual. Isn't it a bug in SeDuMi ?

[siko1056]

However, if the user entered the dual of the 13th MOI test problem to SeDuMi, e.g.,

min -x2
- x1 - 2x2 = 0
  x1 - 3x2 = 0
x1 free
x2 >= 0

he would still get an incorrect dual. Isn't it a bug in SeDuMi ?

Not really a bug. Your matrix A is perfectly square. This makes it hard for SeDuMi to decide what is the constraint dimension m and the cone dimension N. Your problem data is unfortunately not transposed (internally SeDuMi always works with At instead of A). See

sedumi/sedumi.m

Lines 10 to 11 in b4a462b

	% If size(A,2)==length(b), then it solves the linear program
	% MINIMIZE c'*x SUCH THAT A'*x = b, x >= 0

and

sedumi/pretransfo.m

Lines 165 to 167 in b4a462b

	if size(At,2) ~= m
	if m == size(At,1)
	At = At'; %user gave A instead of At.

Just try to input the transposed matrix A.

[blegat]

Indeed, that resolves the issue, thanks @siko1056 !