Commons Math
  1. Commons Math
  2. MATH-842

Investigate why Bland's rule in Simplex Solver still creates cycles

    Details

    • Type: Improvement Improvement
    • Status: Closed
    • Priority: Major Major
    • Resolution: Fixed
    • Affects Version/s: 3.1
    • Fix Version/s: 3.3
    • Labels:
      None

      Description

      As a consequence of MATH-828, Bland's rule has been introduced to prevent cycling. Now there are cases where cycles can still occur (see testMath828Cycle in SimplexSolverTest). These cases have for now been solved with a simple heuristic:

      • after maxIterations / 2 -> ignore Bland's rule

      This issue has been created to further investigate the problem and come up with a cleaner solution.

      1. MATH-842.patch
        7 kB
        Thomas Neidhart

        Issue Links

          Activity

          Hide
          Thomas Neidhart added a comment - - edited

          In r1523495 applied the following changes to the old, deprecated SimplexSolver in the optimization.linear package:

          • apply Bland's rule also for pivot column selection as described in http://www.stanford.edu/class/msande310/blandrule.pdf
          • add getter/setter for cyclePrevention, default is true
          • increase default cut-off threshold to 1e-10
          • added test to demonstrate a cycle of the simplex algorithm

          With these changes, the solver works very stable, also on the problem attached to MATH-828 which resulted in the addition of Bland's rule.

          The change has also to be applied to the new SimplexSolver in the optim.linear package. When copying the code there has been an additional mistake: Bland's rule was never applied as the criteria was based on getEvaluations(), which is never incremented by the simplex solver and not supported.

          Show
          Thomas Neidhart added a comment - - edited In r1523495 applied the following changes to the old, deprecated SimplexSolver in the optimization.linear package: apply Bland's rule also for pivot column selection as described in http://www.stanford.edu/class/msande310/blandrule.pdf add getter/setter for cyclePrevention, default is true increase default cut-off threshold to 1e-10 added test to demonstrate a cycle of the simplex algorithm With these changes, the solver works very stable, also on the problem attached to MATH-828 which resulted in the addition of Bland's rule. The change has also to be applied to the new SimplexSolver in the optim.linear package. When copying the code there has been an additional mistake: Bland's rule was never applied as the criteria was based on getEvaluations(), which is never incremented by the simplex solver and not supported.
          Hide
          Thomas Neidhart added a comment -

          I am unsure how to make the change to the new SimplexSolver:

          • create a class which implements OptimizationData, e.g. CyclePrevention and pass it as part of the call to optimize
          • make it a constructor argument
          • with a getter/setter
          Show
          Thomas Neidhart added a comment - I am unsure how to make the change to the new SimplexSolver: create a class which implements OptimizationData, e.g. CyclePrevention and pass it as part of the call to optimize make it a constructor argument with a getter/setter
          Hide
          Gilles added a comment -

          I do not understand what is done exactly; from the above, it seems that the behaviour is related to "maxIterations"; hence a subclass of "MaxIter" might be appropriate.

          Show
          Gilles added a comment - I do not understand what is done exactly; from the above, it seems that the behaviour is related to "maxIterations"; hence a subclass of "MaxIter" might be appropriate.
          Hide
          Thomas Neidhart added a comment -

          After a second thought, I think I will just remove the getter/setters again and make the prevention mechanism the default one for now.
          The heuristic with half of the max iterations was just a temporary solution, but with a proper implementation of Bland's rule it should be obsolete.

          Show
          Thomas Neidhart added a comment - After a second thought, I think I will just remove the getter/setters again and make the prevention mechanism the default one for now. The heuristic with half of the max iterations was just a temporary solution, but with a proper implementation of Bland's rule it should be obsolete.
          Hide
          Thomas Neidhart added a comment -

          The reason why there is still cycling is identified: Bland's rule was not correctly implemented.
          I will revert my changes so far, as I will not be able to complete them in the next 2 weeks.

          The full fix should be as follows (similar to other lp solvers):

          • support selection of different pivot selection rules
            • Dantzig (the standard selection rule, which is implemented now)
            • Bland (first index selection)
            • Steepest edge
            • Devex
          • support adaptive mode: revert to Bland only when cycling is detected

          There are many more things other solvers like lp_solve or glpk support, and we should consider creating an Options object that bundles all solver related settings so that a user can easily adapt the behavior for his/her purpose.

          Show
          Thomas Neidhart added a comment - The reason why there is still cycling is identified: Bland's rule was not correctly implemented. I will revert my changes so far, as I will not be able to complete them in the next 2 weeks. The full fix should be as follows (similar to other lp solvers): support selection of different pivot selection rules Dantzig (the standard selection rule, which is implemented now) Bland (first index selection) Steepest edge Devex support adaptive mode: revert to Bland only when cycling is detected There are many more things other solvers like lp_solve or glpk support, and we should consider creating an Options object that bundles all solver related settings so that a user can easily adapt the behavior for his/her purpose.
          Hide
          Thomas Neidhart added a comment -

          Patch that contains the first, reverted changes.

          Show
          Thomas Neidhart added a comment - Patch that contains the first, reverted changes.
          Hide
          Gilles added a comment -

          [...] we should consider creating an Options object that bundles all solver related settings [...]

          That furiously looks like the "OptimizationData" which you didn't like.

          Show
          Gilles added a comment - [...] we should consider creating an Options object that bundles all solver related settings [...] That furiously looks like the "OptimizationData" which you didn't like.
          Hide
          Thomas Neidhart added a comment -

          No, not really.

          The options object would be just a container for lots of int/double/boolean/enum fields/flags that would be provided to the solver, either before or during the optimize call. These options are also very specific to the solver, so there is no need to generalize them in one way or another. Implementing OptimizationData would allow to specify it during the call to optimize, but I would also be fine to make it part of the solver interface, e.g. SimplexSolver.withSolverOptions(SimplexSolverOptions).

          What I did not like about the OptimizationData is that we attempted to create a common interface for things that are different, and that it was not clear which optimizer uses or requires which input. It would make sense to unify things if you could replace one instance with another, but this is clearly not the case: you can not use e.g. a BrentOptimizer to solve a linear problem, yet the API allows it (at least theoretically).

          We also did not gain anything wrt thread-safety, as each optimizer stores all the passed OptimizationData in its instance (not that I think thread-safety for algorithm classes should be the main design goal), and I did not see the benefit of providing all the required inputs as part of the call to a generic optimize method.

          What we should have done: exclude the optimize method from the BaseOptimizer class and provide specific ones for classes of optimizers that support the Liskov substitution principle, e.g. the LinearOptimizer defines a optimize method with a signature that all LinearOptimizer instances support. Then you could create e.g. a DualSimplex or PrimalSimplex optimizer both implementing the same interface, and you can exchange one implementation with another.

          Show
          Thomas Neidhart added a comment - No, not really. The options object would be just a container for lots of int/double/boolean/enum fields/flags that would be provided to the solver, either before or during the optimize call. These options are also very specific to the solver, so there is no need to generalize them in one way or another. Implementing OptimizationData would allow to specify it during the call to optimize, but I would also be fine to make it part of the solver interface, e.g. SimplexSolver.withSolverOptions(SimplexSolverOptions). What I did not like about the OptimizationData is that we attempted to create a common interface for things that are different, and that it was not clear which optimizer uses or requires which input. It would make sense to unify things if you could replace one instance with another, but this is clearly not the case: you can not use e.g. a BrentOptimizer to solve a linear problem, yet the API allows it (at least theoretically). We also did not gain anything wrt thread-safety, as each optimizer stores all the passed OptimizationData in its instance (not that I think thread-safety for algorithm classes should be the main design goal), and I did not see the benefit of providing all the required inputs as part of the call to a generic optimize method. What we should have done: exclude the optimize method from the BaseOptimizer class and provide specific ones for classes of optimizers that support the Liskov substitution principle, e.g. the LinearOptimizer defines a optimize method with a signature that all LinearOptimizer instances support. Then you could create e.g. a DualSimplex or PrimalSimplex optimizer both implementing the same interface, and you can exchange one implementation with another.
          Hide
          Gilles added a comment -

          No, not really.

          [...] e.g. SimplexSolver.withSolverOptions(SimplexSolverOptions).

          If withSolverOptions is defined in an interface, then SimplexSolverOptions will most probably be an implementation of an interface (say, Options). Options is thus equivalent to OptimizationData: a marker interface (since you indicated that a bunch of options are implementation-specific).

          What I did not like about the OptimizationData is that we attempted to create a common interface for things that are different, and that it was not clear which optimizer uses or requires which input.

          "Not clear" could refer to a documentation issue, to better highlight required "options" (which is why "OptimizationData" was not named "OptionData").

          The purpose was to allow the passing of the same list of options to different optimizers, and each would pick the data it needs. I, at least, had a use for that; but I can agree that it is not a common pattern, especially in strongly typed languages.

          you can not use e.g. a BrentOptimizer to solve a linear problem, yet the API allows it [...]

          This was a trade-off in order to inherit code which all optimizers need (e.g. the counters). Sometimes multiple inheritance is useful.
          If we completely depart from the original design (based on a single hierarchy), many improvements become possible (as with Evan's proposal, and what Luc proposes for ODE, IIUC).

          What we should have done [...]

          No, we could not have done that; because nobody proposed it.

          Show
          Gilles added a comment - No, not really. [...] e.g. SimplexSolver.withSolverOptions(SimplexSolverOptions). If withSolverOptions is defined in an interface, then SimplexSolverOptions will most probably be an implementation of an interface (say, Options ). Options is thus equivalent to OptimizationData : a marker interface (since you indicated that a bunch of options are implementation-specific). What I did not like about the OptimizationData is that we attempted to create a common interface for things that are different, and that it was not clear which optimizer uses or requires which input. "Not clear" could refer to a documentation issue, to better highlight required "options" (which is why "OptimizationData" was not named "OptionData"). The purpose was to allow the passing of the same list of options to different optimizers, and each would pick the data it needs. I, at least, had a use for that; but I can agree that it is not a common pattern, especially in strongly typed languages. you can not use e.g. a BrentOptimizer to solve a linear problem, yet the API allows it [...] This was a trade-off in order to inherit code which all optimizers need (e.g. the counters). Sometimes multiple inheritance is useful. If we completely depart from the original design (based on a single hierarchy), many improvements become possible (as with Evan's proposal, and what Luc proposes for ODE, IIUC). What we should have done [...] No, we could not have done that; because nobody proposed it.
          Hide
          Thomas Neidhart added a comment -

          I dont see the need to make it generic, withOptions(SimplexSolverOptions) should be specific to the SimplexSolver. Other solvers might have other options specific to them, that alter the behavior of the solver.

          Without specifying such options the solver would be initialized with default options.

          I understand that there are "generic" options almost each solver supports, like maxIterations, but thats almost everything the solver have in common, afaict.

          Show
          Thomas Neidhart added a comment - I dont see the need to make it generic, withOptions(SimplexSolverOptions) should be specific to the SimplexSolver. Other solvers might have other options specific to them, that alter the behavior of the solver. Without specifying such options the solver would be initialized with default options. I understand that there are "generic" options almost each solver supports, like maxIterations, but thats almost everything the solver have in common, afaict.
          Hide
          Thomas Neidhart added a comment -

          Added PivotSelectionRule in r1550975.

          The SimplexSolver now supports two rules:

          • Dantzig, default
          • Bland: prevent cycles but may be slower to find a optimal solution
          Show
          Thomas Neidhart added a comment - Added PivotSelectionRule in r1550975. The SimplexSolver now supports two rules: Dantzig, default Bland: prevent cycles but may be slower to find a optimal solution
          Hide
          Luc Maisonobe added a comment -

          Closing all resolved issue now available in released 3.3 version.

          Show
          Luc Maisonobe added a comment - Closing all resolved issue now available in released 3.3 version.

            People

            • Assignee:
              Thomas Neidhart
              Reporter:
              Thomas Neidhart
            • Votes:
              0 Vote for this issue
              Watchers:
              4 Start watching this issue

              Dates

              • Created:
                Updated:
                Resolved:

                Development