Simplex法 randomizeしてワーストケース回避

hitonanode · hitonanode · commit eae8daddb002 · 2021-02-28T16:35:09.000+09:00
diff --git a/combinatorial_opt/simplex.hpp b/combinatorial_opt/simplex.hpp
@@ -1,15 +1,19 @@
 #pragma once
+#include <algorithm>
+#include <chrono>
 #include <cmath>
 #include <numeric>
+#include <random>
 #include <vector>
 
 // CUT begin
 // Maximize cx s.t. Ax <= b, x >= 0
 // Implementation idea: https://kopricky.github.io/code/Computation_Advanced/simplex.html
 // Refer to https://hitonanode.github.io/cplib-cpp/combinatorial_opt/simplex.hpp
-template <typename Float = double, int DEPS = 30> struct Simplex {
+template <typename Float = double, int DEPS = 30, bool Randomize = true> struct Simplex {
     const Float EPS = Float(1.0) / (1LL << DEPS);
     int N, M;
+    std::vector<int> shuffle_idx;
     std::vector<int> idx;
     std::vector<std::vector<Float>> mat;
     int i_ch, j_ch;
@@ -35,7 +39,7 @@ template <typename Float = double, int DEPS = 30> struct Simplex {
     inline Float abs_(Float x) noexcept { return x > -x ? x : -x; }
     void _solve() {
         std::vector<int> jupd;
-        for (j_ch = N;;) {
+        for (nb_iter = 0, j_ch = N;; nb_iter++) {
             if (i_ch < M) {
                 std::swap(idx[j_ch], idx[i_ch + N + 1]);
                 mat[i_ch][j_ch] = Float(1) / mat[i_ch][j_ch];
@@ -90,11 +94,38 @@ template <typename Float = double, int DEPS = 30> struct Simplex {
     }
 
 public:
-    Simplex(const std::vector<std::vector<Float>> &A, const std::vector<Float> &b, const std::vector<Float> &c) {
+    Simplex(std::vector<std::vector<Float>> A, std::vector<Float> b, std::vector<Float> c) {
         is_infty = infeasible = false;
+
+        if (Randomize) {
+            std::mt19937 rng(std::chrono::steady_clock::now().time_since_epoch().count());
+
+            std::vector<std::pair<std::vector<Float>, Float>> Abs;
+            for (unsigned i = 0; i < A.size(); i++) Abs.emplace_back(A[i], b[i]);
+            std::shuffle(Abs.begin(), Abs.end(), rng);
+            A.clear(), b.clear();
+            for (auto &&Ab : Abs) A.emplace_back(Ab.first), b.emplace_back(Ab.second);
+
+            shuffle_idx.resize(c.size());
+            std::iota(shuffle_idx.begin(), shuffle_idx.end(), 0);
+            std::shuffle(shuffle_idx.begin(), shuffle_idx.end(), rng);
+            auto Atmp = A;
+            auto ctmp = c;
+            for (unsigned i = 0; i < A.size(); i++) {
+                for (unsigned j = 0; j < A[i].size(); j++) A[i][j] = Atmp[i][shuffle_idx[j]];
+            }
+            for (unsigned j = 0; j < c.size(); j++) c[j] = ctmp[shuffle_idx[j]];
+        }
+
         _initialize(A, b, c);
         _solve();
+
+        if (Randomize) {
+            auto xtmp = x;
+            for (unsigned j = 0; j < c.size(); j++) x[shuffle_idx[j]] = xtmp[j];
+        }
     }
+    unsigned nb_iter;
     bool is_infty;
     bool infeasible;
     std::vector<Float> x;
diff --git a/combinatorial_opt/simplex.md b/combinatorial_opt/simplex.md
@@ -59,14 +59,19 @@ $
 
 - ひとまずその時点で $c_j$ が正となるものを選ぶ．
   - 複数存在する場合は添字番号最小（Bland's rule）
-- $\mathbf{c}$ の全要素が $0$ ならこれ以上解が改善しないので終了
+- $\mathbf{c}$ の全要素が $0$ ならこれ以上解が改善しないので終了．
 
 #### $y_i$ の決定
 
 - 最も不等式制約が厳しいものを選ぶ．
   - さもなくば原点が実行可能領域から飛び出るので
   - 複数存在する場合は添字番号最小（Bland's rule）
 
+#### ワーストケースの回避
+
+- 指数回のステップを要する恣意的なケースが知られている．
+- 特に [1] の p.43 (2.71) で挙げられているケースに関しては，変数・制約の順序をランダムに取り換えることで回避が可能であることを経験的に確認したので，デフォルトでこれを行うことにした．
+
 ## 問題例
 
 - [AtCoder jag2018summer-day3 H - Enlarge Circles](https://atcoder.jp/contests/jag2018summer-day3/tasks/jag2018summer_day3_h)
