hi everyone,
I've been using the solver for a month or two now and there are many things
I really like about it, in particular the cleanliness of the code, it is
very easy to see what is going on and make experimental changes. The main
reason I'm writing is to ask if there is any built-in provision for avoiding
degeneracy cycling and what people's plans/thoughts are on the issue. I
don't mean the problem of losing feasibility, re-entering phase 1 with
consequent worsening of the objective function and thereby repeating a
basis, I'm talking about when you perform a number of consecutive degenerate
pivots that don't improve the objective and eventually return to the same
basis.
The manual and the code don't make any reference to this, so I assume there
isn't any anti-cycling, for a while I was thinking that the head[m+1..m+n]
array of non-basic variables was maintaining an implicit ordering somewhat
like Bland's rule, but then I realized that there is no ordering on the
leaving variables, see chuzr in glpspx01.c for example, in case of a tie
between leaving variables it tries to choose the one with the largest
coefficient in that column of the tableau. (Note: No tolerance is used
when making the test for a tie, I suggest adding one could enable choosing a
larger coefficient and thus a better conditioned matrix, does everyone
agree?). So I guess cycling can occur, has anybody tried entering one of
the textbook examples to see whether it really happens?
So anyway, I tried implementing Bland's rule, it was easy to do, I just
check if the previous pivot was degenerate, and if so, use Bland's pivot
selection instead of the normal one. Unfortunately this unreasonably hurt
performance on the normal (non cycling) case, i.e. it led to lengthy stalls,
as my problems contain a lot of degeneracy and the Bland pivot rule doesn't
attempt to choose a row with a good reduced cost for entry. So I also tried
the Cacioppi & Schrage rule described by Richard Kipp Martin (Large Scale
Linear & Integer Optimization, p171-172, you can find this in Google
Books). I don't know if I implemented it correctly though, because it's
supposed to address this problem by forcing very few restrictions on the
entering/leaving variables, but I observed much stalling still.
Eventually what solved the problem for me was a perturbation of the
right-hand sides, note that I hadn't actually observed cycling in practice
but I did observe lengthy stalls with the standard code, the perturbation
fixed this and led to a dramatic improvement. I did it by adding, to each
lb[1..m] and ub[1..m] pair, a random number in the range -1e-5..1e-5. I
didn't bother removing the perturbation afterwards, as the accuracy of the
results was suitable for my purpose, but I would suggest that if doing e.g.
primal simplex, then after reaching optimality we have a primal/dual
feasible solution, and removing the perturbation amounts to changing the
dual objective, so we wouldn't lose dual feasibility and we could continue
to optimize with dual simplex until the true optimum is reached.
Another way of doing the perturbation is to note that the solver doesn't
normally work with a right-hand side, e.g. an equation like x+y+z<=5 is
changed to s-x-y-z=0 with s<=5, rather than x+y+z-s=5 with s>=0 as in the
textbook way of doing things. We could introduce a right-hand side,
containing only the small perturbations, so e.g. the calculation of x_B
(bbar) would be done as
x_B = A_B^{-1} (b - A_N x_N)
where b is the perturbation, rather than
x_B = -A_B^{-1} A_N x_N
as we are currently doing. This kind of approach might be a bit more
elegant and avoid having to save the lb/ub arrays for when the perturbation
is removed.
Some other suggestions: We maintain a copy of the row-representation of N
(as suggested by Bixby) for the cbar and gamma recurrences, but the code to
delete a column of N as it's pivoted into the basis doesn't look
particularly efficient, I made a version where the N is constructed once at
the start, really just a row-representation of the full coefficient matrix
(I | -A), and extra columns are simply ignored when updating cbar and
gamma. I also made the cbar and gamma arrays a bit bigger so there is a
space for every variable rather than only the non-basic variables, so that
they correspond with the columns of (I | -A) and make the addressing a bit
easier. For the gamma array I put an inner-loop test to avoid updating the
basic variables, but for cbar I didn't bother with such a test.
I also think more diagnostic output would be good, I added a counter for how
many degenerate pivots have occurred (next to the iteration counter), I'm
also planning to output a count of how many nonzeros occur in the basis
factorization, as I'm keen to know the effect of different pivoting
strategies (including adding the tolerance in chuzr that I mentioned
earlier) on the sparseness of the basis inverse. Does anyone have a good
reference for a devex-like pivoting strategy that attempts to keep the basis
as sparse as possible? Because I think this would be really helpful for my
(very large) problems. Another change I was thinking of making, is that
when we construct the pivot row, we should apply a tolerance in the test for
zero, improving the sparseness, but I haven't tried this yet.
Another thing I did was to implement a GLP-native format for storing the LP,
I did this rather reluctantly as the world doesn't really need yet another
competing format, but the problem with MPS is that it doesn't save the
minimization/maximization flag and with CPLEX LP format it doesn't save the
objective function coefficient. Also, at least with CPLEX LP, the rows are
loaded in a different order (they are entered as they are encountered in the
file), which wasn't acceptable for my purpose, in fact it was very
misleading and I wasted a fair bit of time wondering why my post-solution
analyzer wasn't working properly. The new format is very similar to the
solution files already output, I did it by modifying glp_write_sol and
glp_read_sol to create glp_write_glp and glp_read_glp. I also added a new
--basis option for glpsol, so that it can read (non-optimal) a solution file
written by glp_write_sol and use it as the starting basis.
One last comment, the current arrangements aren't very good for a caller
that wishes to solve, make some changes in the model, resolve etc. In
particular I'm talking about the init_csa routine e.g. in glpspx01.c, which
insists on making a copy of the entire problem every time it is entered,
could we make this incremental and only merge in the changes to the problem
to its private data structures?
Sorry, this turned into a bit of an essay, but what does everyone think?
Should I upload some of my code? It's not particularly clean at the moment,
but it would probably be a starting point.
cheers, Nick
_______________________________________________
Bug-glpk mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-glpk
