Loading docs/math/linear-programming.md +0 −184 Original line number Diff line number Diff line Loading @@ -106,187 +106,3 @@ 显而易见的，打了蓝色斜线的区域为三块区域的交集，这就是这个线性规划的所有可行解。因为题目中说明，需要最小化 x 和 y，观察图像可知，点（2,3）为可行解中 x 和 y 最小的一个。因此， $x_{min}=2,y_{min}=3$ 把求解线性规划的图解法总结起来，就是先在坐标系中作出所有的约束条件，然后作出需要求极值的线性函数的定义域。定义域与约束条件的交集就是这个线性规划的解集，而所需求的极值由观察可以轻易得出。 ### 单纯形法 > $$ > \forall k,\sum_{i}a_{k,i}x_i\le b_{k}\\ > \forall i,x_i\ge 0 > $$ > > 求 > > $$ > \max(\sum_{i}c_{i}x_i) > $$ 我们考虑一个下标从 $0$ 开始的矩阵 $A_{0\sim m,0\sim n+m}$ : $$ A= \begin{bmatrix} 0 & c_1 & \ldots & c_n & 0 & 0 & \ldots & 0 \\ b_1 & a_{1,1} & \ldots & a_{1,n} & 1 & 0 & \dots & 0 \\ b_2 & a_{2,1} & \ldots & a_{2,n} & 0 & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ b_m & a_{m,1} & \ldots & a_{m,n} & 0 & 0 & \dots & 1 \end{bmatrix} $$ 其中 $A_{0,0}$ 为这个状态和初始时答案之差，新增的 $m$ 行相当于新建 $m$ 个变量，编号为 $x_{n+1\sim n+m}$ ，要求它们不小于 $0$ 于是 $$ \forall i\in[1,m]\\ \sum_{1\le j\le n+m}a_{i,j}x_j=b_i $$ 显然，要先求出一套可行的方案 如果 $$ x_i=\begin{cases}0 & i\in[1,n]\\b_{i-n} & i\in[n+1,n+m]\end{cases} $$ 那么显然满足，但是由于 $b_i$ 不一定不小于 $0$ ，所以这样可能导致 $x_{i}$ 小于 $0$ 然后，介绍一个魔法，通过矩阵变化改变 $b$ 的值来求合法解/最优解 #### 转轴 ~~之所以叫这个名称，是因为可行域是一个超立方体，这个魔法是经过一条轴，但是我不讲~~ 这个魔法是指： 选择两个数 $u,v$ ，考虑将 $x_v$ 变为非 $0$ 数而 $x_{u+n}$ 变为 $0$ ，也就是让 $x_v$ 来代替 $x_{u+n}$ 作为等于 $b_u$ 的位置。 如果要实现，就要运用矩阵变换。 我们可以先让第 $u$ 行全部除一个 $A_{u,v}$ ，然后把其它所有行 $i$ 的位置 $v$ 变为 $0$ （即使 $A_{i,v}=0$ ），可以直接把第 $i$ 行减去 $A_{i,v}$ 倍第 $u$ 行 然后就可以了！ 因为（不算 $u+n,v$ 列的话）后面 $m$ 列只有 $A_{i,i+n}=1$ ，前面 $n$ 列都是 $0$ ，所以矩阵变化不会影响其它的 $x$ 。 当然，最好把第 $u+n,v$ 行换一下，这样会方便很多 （还有，因为矩阵后 $m$ 列值的特殊性，一般不存这 $m$ 列） 顺便一提，这对答案的贡献为： $$ \dfrac{A_{0,v}A_{u,0}}{A_{u,v}} $$ ```cpp int id[N << 1]; double a[N][N]; void pivot(int u, int v) { double k = a[u][v]; swap(id[u + n], id[v]); a[u][v] = 1.; //已经知道第v列要被消掉了，就可以直接变为u+n列了（尽管要用一遍，就临时存一下） for (int i = 0; i <= n; i++) a[u][i] /= k; for (int i = 0; i <= m; i++) if (i != u && fabs(a[i][v]) > eps) { k = a[i][v]; a[i][v] = 0; //同上 for (int j = 0; j <= n; j++) a[i][j] -= k * a[u][i]; } } ``` 知道了基本的变（魔）化（法），就来看怎么求一个合法解 #### 初始化 引用一下 Candy? 的博客 > 算法导论上有一个辅助线性规划的做法 > > 但我发现好多人都用了 **随机初始化** 的黑科技 > > 在所有 $a_{i,0}<0$ 的约束中随机选一个作为 $u$ ，再随机选一个 $a_{u,v}<0$ 作为 $v$ ，然后 `pivot(u,v)` 后 $a_{i,0}$ 就变正了…… 但这可能出现循环，然后不断损精度，然后 WA/TLE，总之很玄学就对了 但是有「优先前面/后面」的选择倾向再加以小随机效果似乎不错 ```cpp int init() { while (1) { int u = 0, v = 0; for (int i = 1; i <= m; i++) if (a[i][0] < -eps && (!u || rand() % 3 == 0)) u = i; if (!u) //说明这个解是可行解 return 1; for (int i = 1; i <= n; i++) if (a[u][i] < -eps && (!v || rand() % 3 == 0)) v = i; if (!v) { //说明无解 printf("Infeasible\n"); return 0; } pivot(u, v); } } ``` #### 求解 每次找到一个 $v$ 使得 $A_{0,v}>0$ ，再找到一个 $u$ 使得 $A_{u,v}>0$ ，然后 `pivot(u,v)` 即可。 ##### 退化 在旋转的过程中，可能会存在保持答案不变的情况，这种现象称为退化。 ##### Bland 规则 在选择替入变量和替出变量的时候，我们总是选择满足条件的下标最小值。 - 替入变量 $v$ ：目标条件中，系数为正数（如果求最小值就是负数）的第一个作为替入变量 - 替出变量 $u$ ：对所有的约束条件中，选择对 $u$ 约束最紧的第一个（也就是 $\dfrac{A_{u,0}}{A_{u,v}}$ 最小的那个） ```cpp int simplex() { while (1) { int u = 0, v = 0; double mi = inf; for (int i = 1; i <= n; i++) if (a[0][i] > eps) { v = i; return; } if (!v) return 1; for (int i = 1; i <= m; i++) if (a[i][v] > eps && a[i][0] / a[i][v] < mi) { // bland法则 mi = a[i][0] / a[i][v]; u = i; } if (!u) { //没有约束 printf("Unbounded\n"); return 0; } pivot(u, v); } } ``` #### 其它说明 ##### 为什么这能对 每次矩阵变化后的子问题是等价于初始问题的（易证） 并且这类问题有一个特殊性：它的可行域是凸的 意味着不会有除最优解外的局部最优解出现 ##### 这能跑多大 ~~心有多大，梦想有多大，它就能跑多大~~ 对于这种玄学算法，鬼才知道它能跑多大。 不过一般来说， $nm\le10^6$ 是~~不~~可行的 Loading
docs/math/linear-programming.md +0 −184 Original line number Diff line number Diff line Loading @@ -106,187 +106,3 @@ 显而易见的，打了蓝色斜线的区域为三块区域的交集，这就是这个线性规划的所有可行解。因为题目中说明，需要最小化 x 和 y，观察图像可知，点（2,3）为可行解中 x 和 y 最小的一个。因此， $x_{min}=2,y_{min}=3$ 把求解线性规划的图解法总结起来，就是先在坐标系中作出所有的约束条件，然后作出需要求极值的线性函数的定义域。定义域与约束条件的交集就是这个线性规划的解集，而所需求的极值由观察可以轻易得出。 ### 单纯形法 > $$ > \forall k,\sum_{i}a_{k,i}x_i\le b_{k}\\ > \forall i,x_i\ge 0 > $$ > > 求 > > $$ > \max(\sum_{i}c_{i}x_i) > $$ 我们考虑一个下标从 $0$ 开始的矩阵 $A_{0\sim m,0\sim n+m}$ : $$ A= \begin{bmatrix} 0 & c_1 & \ldots & c_n & 0 & 0 & \ldots & 0 \\ b_1 & a_{1,1} & \ldots & a_{1,n} & 1 & 0 & \dots & 0 \\ b_2 & a_{2,1} & \ldots & a_{2,n} & 0 & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ b_m & a_{m,1} & \ldots & a_{m,n} & 0 & 0 & \dots & 1 \end{bmatrix} $$ 其中 $A_{0,0}$ 为这个状态和初始时答案之差，新增的 $m$ 行相当于新建 $m$ 个变量，编号为 $x_{n+1\sim n+m}$ ，要求它们不小于 $0$ 于是 $$ \forall i\in[1,m]\\ \sum_{1\le j\le n+m}a_{i,j}x_j=b_i $$ 显然，要先求出一套可行的方案 如果 $$ x_i=\begin{cases}0 & i\in[1,n]\\b_{i-n} & i\in[n+1,n+m]\end{cases} $$ 那么显然满足，但是由于 $b_i$ 不一定不小于 $0$ ，所以这样可能导致 $x_{i}$ 小于 $0$ 然后，介绍一个魔法，通过矩阵变化改变 $b$ 的值来求合法解/最优解 #### 转轴 ~~之所以叫这个名称，是因为可行域是一个超立方体，这个魔法是经过一条轴，但是我不讲~~ 这个魔法是指： 选择两个数 $u,v$ ，考虑将 $x_v$ 变为非 $0$ 数而 $x_{u+n}$ 变为 $0$ ，也就是让 $x_v$ 来代替 $x_{u+n}$ 作为等于 $b_u$ 的位置。 如果要实现，就要运用矩阵变换。 我们可以先让第 $u$ 行全部除一个 $A_{u,v}$ ，然后把其它所有行 $i$ 的位置 $v$ 变为 $0$ （即使 $A_{i,v}=0$ ），可以直接把第 $i$ 行减去 $A_{i,v}$ 倍第 $u$ 行 然后就可以了！ 因为（不算 $u+n,v$ 列的话）后面 $m$ 列只有 $A_{i,i+n}=1$ ，前面 $n$ 列都是 $0$ ，所以矩阵变化不会影响其它的 $x$ 。 当然，最好把第 $u+n,v$ 行换一下，这样会方便很多 （还有，因为矩阵后 $m$ 列值的特殊性，一般不存这 $m$ 列） 顺便一提，这对答案的贡献为： $$ \dfrac{A_{0,v}A_{u,0}}{A_{u,v}} $$ ```cpp int id[N << 1]; double a[N][N]; void pivot(int u, int v) { double k = a[u][v]; swap(id[u + n], id[v]); a[u][v] = 1.; //已经知道第v列要被消掉了，就可以直接变为u+n列了（尽管要用一遍，就临时存一下） for (int i = 0; i <= n; i++) a[u][i] /= k; for (int i = 0; i <= m; i++) if (i != u && fabs(a[i][v]) > eps) { k = a[i][v]; a[i][v] = 0; //同上 for (int j = 0; j <= n; j++) a[i][j] -= k * a[u][i]; } } ``` 知道了基本的变（魔）化（法），就来看怎么求一个合法解 #### 初始化 引用一下 Candy? 的博客 > 算法导论上有一个辅助线性规划的做法 > > 但我发现好多人都用了 **随机初始化** 的黑科技 > > 在所有 $a_{i,0}<0$ 的约束中随机选一个作为 $u$ ，再随机选一个 $a_{u,v}<0$ 作为 $v$ ，然后 `pivot(u,v)` 后 $a_{i,0}$ 就变正了…… 但这可能出现循环，然后不断损精度，然后 WA/TLE，总之很玄学就对了 但是有「优先前面/后面」的选择倾向再加以小随机效果似乎不错 ```cpp int init() { while (1) { int u = 0, v = 0; for (int i = 1; i <= m; i++) if (a[i][0] < -eps && (!u || rand() % 3 == 0)) u = i; if (!u) //说明这个解是可行解 return 1; for (int i = 1; i <= n; i++) if (a[u][i] < -eps && (!v || rand() % 3 == 0)) v = i; if (!v) { //说明无解 printf("Infeasible\n"); return 0; } pivot(u, v); } } ``` #### 求解 每次找到一个 $v$ 使得 $A_{0,v}>0$ ，再找到一个 $u$ 使得 $A_{u,v}>0$ ，然后 `pivot(u,v)` 即可。 ##### 退化 在旋转的过程中，可能会存在保持答案不变的情况，这种现象称为退化。 ##### Bland 规则 在选择替入变量和替出变量的时候，我们总是选择满足条件的下标最小值。 - 替入变量 $v$ ：目标条件中，系数为正数（如果求最小值就是负数）的第一个作为替入变量 - 替出变量 $u$ ：对所有的约束条件中，选择对 $u$ 约束最紧的第一个（也就是 $\dfrac{A_{u,0}}{A_{u,v}}$ 最小的那个） ```cpp int simplex() { while (1) { int u = 0, v = 0; double mi = inf; for (int i = 1; i <= n; i++) if (a[0][i] > eps) { v = i; return; } if (!v) return 1; for (int i = 1; i <= m; i++) if (a[i][v] > eps && a[i][0] / a[i][v] < mi) { // bland法则 mi = a[i][0] / a[i][v]; u = i; } if (!u) { //没有约束 printf("Unbounded\n"); return 0; } pivot(u, v); } } ``` #### 其它说明 ##### 为什么这能对 每次矩阵变化后的子问题是等价于初始问题的（易证） 并且这类问题有一个特殊性：它的可行域是凸的 意味着不会有除最优解外的局部最优解出现 ##### 这能跑多大 ~~心有多大，梦想有多大，它就能跑多大~~ 对于这种玄学算法，鬼才知道它能跑多大。 不过一般来说， $nm\le10^6$ 是~~不~~可行的
