Re: [Help-glpk] GLPK primal simplex never terminates but dual simplex solves it in seconds

[Andrew Makhorin]

> I have an LP model which seems to give troubles to GLPK's primal
> simplex algorithm. Here's the terminal output:

> Preprocessing...
> 1403 rows, 1884 columns, 2383489 non-zeros
> Scaling...
>  A: min|aij| = 1.000e-003  max|aij| = 1.000e+000  ratio = 1.000e+003
> GM: min|aij| = 3.164e-002  max|aij| = 3.161e+001  ratio = 9.990e+002
> EQ: min|aij| = 1.001e-003  max|aij| = 1.000e+000  ratio = 9.990e+002
> Constructing initial basis...
> Size of triangular part = 1403
> *     0: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *   500: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  1000: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  1500: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  2000: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  2500: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  3000: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  3500: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  4000: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  4500: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> *  5000: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
> .... and so on

> However, if I switch it to dual simplex then it cracks it in seconds:

> Preprocessing...
> 1403 rows, 1884 columns, 2383489 non-zeros
> Scaling...
>  A: min|aij| = 1.000e-003  max|aij| = 1.000e+000  ratio = 1.000e+003
> GM: min|aij| = 3.164e-002  max|aij| = 3.161e+001  ratio = 9.990e+002
> EQ: min|aij| = 1.001e-003  max|aij| = 1.000e+000  ratio = 9.990e+002
> Constructing initial basis...
> Size of triangular part = 1403
>       0:                          infeas = 4.780e+005 (0)
>       1:                          infeas = -2.274e-013 (0)
> |   500: obj =  1.000000000e+000  infeas = 2.045e-029 (0)
> |   718: obj =  1.000000000e+000  infeas = 8.318e-012 (0)
> OPTIMAL SOLUTION FOUND
> Time used:   4.4 secs
> Memory used: 220.4 Mb (231146952 bytes)

> I suspect that the problem could be due to potential numerical
> instability

No, I don't think so. Looks like your instance is highly degenerate
in the primal space (for example, if the feasible set is
 {x : Ax = b, x >= 0}, the vector b either is zero or has only few
non-zero elements).

>  because I have some small variable bounds (~ 4.0 10^-8)

Such bounds are equivalent to zero bounds, because by default
tol_bnd = 1e-7, i.e. if some basic variable x >= 0 takes on, say,
the value -0.7e-8, the simplex solver does not consider the zero
bound of x as violated.

> but it's not reported by simplex. Can there be another reason or I
> should just do some manual scaling to prevent it?

You could try to replace zero bounds by some non-zero value, say, by
1e-5 or even 1e-3, solve the perturbed instance with the primal simplex
to optimality, restore zero bounds (that will make the instance primal
infeasible), and then use the dual simplex to find optimal solution of
the original non-perturbed instance.

> Let me know if you want to see the model (it's >1M so I'm not
> attaching it to this email).

Please gzip it and post to me. Thanks.


Andrew Makhorin