numerical (in)stability

[Ali Baharev]

Dear Andrew,

Thank you for your response. I decided to move our conversation

http://lists.gnu.org/archive/html/bug-glpk/2009-10/msg00018.html

here from the bug-glpk because we are not dealing with a bug anymore.

>> I just added a limit on the number of restarts in the simplex
>> algorithm. In this way i do avoid the infinite loop. Could it be
>> officially supported?
>
> In principle, yes. However, there is a technical problem. The point
> is that the numerical instability detected by the simplex solver,
> as a rule, does not lead to cycling, and it would be difficult to
> distinguish between cycling and non-cycling cases due to possible
> degeneracy.

Yes, i do understand that. My original idea was to offer an API
function that can set a bound on the restarts if the user wishes to do
so. By default, the number of restarts could be unlimited.

By the way, could you please check what i did? (The attached C
source.) It seems to work correctly but i do not known the internal
buildup of GLPK.

In other words: what functions should be called before jumping out
from the main loop? Am i allowed to jump out like this?

>> It would be nice to have a function that makes a hexadecimal snapshot
>> of the internal state of the solver and writes it to the disk so it
>> can be reproduced.
>
> Hexadecimal floats is a relatively new language feature, and I am
> not sure that it is supported by all C/C++ compilers used to-day.
> (Though this could be detected by configure script.)

Yes, that is true, i had a lot of pain with the compilers not
supporting this... How about making it optional? I mean that one could
pass a flag to the configure script and then the configure would check
if it works properly.

http://publib.boulder.ibm.com/infocenter/zos/v1r9/topic/com.ibm.zos.r9.cbcmg01/hexflpnot.htm

By the way, GCC with -std=c99 supports both reading and writing. I had
no problem writing with VS2005 but somehow reading failed, probably
requires some obscure flag.

>> Could you please help me how my LP problems should be change to
>> improve their numerical properties?
>
> Then you should explain me your source problem and how you model it
> with lp.

Unfortunately i cannot do anything with the source, it is external.
The LPs are automatically generated by linearizing arbitrary nonlinear
functions originally written in AMPL.

However the structure of the LP has some unique properties.

All columns have -1, 1 lower and upper bounds, resp. All rows are
either double bounded or fixed. The objective contains at most one
nonzero coefficient. If there is a nonzero coefficient in the
objective, it is either -1 or 1.

The number of columns equals the number of rows (it may change in the future).

>> I already removed the tiny coefficients as you proposed and i suspect
>> now the tiny row bounds are causing the trouble.
>
> Perhaps, yes, because row bounds are used as right-hand sides on
> ftran and btran, and if the basis matrix is ill-conditioned, non-zero
> rhs's may produce large round-off errors in the basic solution.

Let's try this one first.

I did realize that when i get into an infinite loop the otherwise
"normal" row bounds are tiny. It may or may not be a coincident
though. The problem might be that some row bounds are huge compared to
the others.

All rows are either double bounded or fixed. What shall i do with the
row bounds (tiny or huge)?

>> An example is at the link i sent yesterday. The main.cpp takes the
>> dump.txt as the input, builds the problem, writes it as dump.lp then
>> starts the iteration and produces the log.txt console output. I can
>> only reproduce this in the Windows environment.
>
> I need a time to check your code.

Thank you for your kind help,

Ali

[glpspx01_restart_limit.c]

---

/* glpspx01.c (primal simplex method) */

...

static int restart_counter = 0;
static int restart_limit = 0;

void set_restart_limit(const int limit) {
        restart_limit = limit;
}

int spx_primal(glp_prob *lp, const glp_smcp *parm)
{
          ...

          restart_counter = 0;
loop: /* main loop starts here */
      /* compute factorization of the basis matrix */
         
          ...

         /* make sure that the current basic solution remains primal
            feasible (or pseudo feasible on phase I) */
         if (check_stab(csa, parm->tol_bnd))
         {  /* there are excessive bound violations due to round-off
               errors */
            if (parm->msg_lev >= GLP_MSG_ERR)
               xprintf("Warning: numerical instability (primal simplex,"
                  " phase %s)\n", csa->phase == 1 ? "I" : "II");


                        ++restart_counter;

                        if (restart_counter > restart_limit) {
                           store_sol(csa, lp, GLP_UNDEF, GLP_UNDEF, 0);
               ret = GLP_EFAIL;
               goto done;
                        }

            /* restart the search */
            csa->phase = 0;
            binv_st = 0;
            rigorous = 5;
            goto loop;
         }

          ...
         
done: /* deallocate the common storage area */
      free_csa(csa);
      /* return to the calling program */
      return ret;
}

/* eof */

---

_______________________________________________
Help-glpk mailing list
Help-glpk@...
http://lists.gnu.org/mailman/listinfo/help-glpk