线性演算(一)

线性演算(Linear Arithmetic)

线性演算文法定义: \[ \begin{aligned} formula&:formula\land formula\ |\ (formula)\ |\ atom\\ atom&:sum\ op\ sum\\ op&:=\ |\ \le\ |\ < \\ sum&:term\ |\ sum+term\\ term&:identifier\ |\ constant\ |\ constant\ identifier \end{aligned} \] 根据文法构造良定义的线性演算公式, 如下: \[ 3x_1+2x_2\le5x_3\ \land\ 2x_1-2x_2=0. \] 从上面的公式可以看出, 等式逻辑(equality logic)是作为线性演算的一部分的.

1 线性演算的求解器

Simplex 方法是针对数值优化的一个算法. 给定实数变量的线性约束的合取公式, Simplex 算法用以找到目标函数的最佳值. 目标函数和约束合称为线性程序(Linear program). 我们只关注判定问题而不是优化问题, 所以, 讨论对于 Simplex 方法的一种变体, 称作通用 Simplex 方法(general Simplex). 通用 Simplex 方法接收在实数集上的线性约束的合取公式, 并且没有目标函数, 并判定集合的可满足性质.另外, 整型线性规划(Integer linear programming, ILP)也是一类基于整型的约束问题.

2 Simplex 算法

Simplex 算法, 判定弱线性不等式合取公式(weak linear inequalities)的可满足行问题.就公式的变量而言, 约束集常伴有一个线性目标函数. 如果约束集合是可满足的, Simplex 算法就提供一个可满足的赋值, 该赋值使得目标函数的值最大化. Simplex 算法是最坏指数复杂度的算法. 也有多项式时间的算法解决此类问题(ellipsoid method, 1979). Simplex 算法在实践中仍然是一个值得考虑的非常高效并广泛使用的方法, 显然是因为现实中指数复杂度的问题很少.

2.1 范式

因为只考虑判定问题而不是优化问题, 所以涉及到的只是 Simplex 算法的变型, 通用 Simplex(general Simplex). 后者并不涉及到目标函数求最值问题. 通用 Simplex 算法接收两种约束类型作为输入:

  1. 等式形式的公式 \[ a_1x_1+...+a_nx_n=0. \]

  2. 变量的上下界 \[ l_i\le x_i\le u_i. \]

以上输入公式的形式称为一般形式(general form). 一般形式的公式并入限制弱线性约束的建模能力, 因为我们可以将任意的弱线性约束转换成上述的形式. 所谓弱线性约束即, \(L\Join R\), 其中的 \(\Join\in\{=,\le,\ge\}\). \(m\) 表示约束的个数. 对于第 \(i\) 个约束, \(1\le i\le m\):

  1. \(R\) 中的所有加数移到左侧以得到 \(L^\prime\Join b\), 其中 \(b\) 是一个常数.

  2. 引入一个新的变量 \(s_i\). 添加约束 \[ L^\prime-s_i=0\ {\rm and}\ s_i\Join b. \] 如果 \(\Join\) 是等式运算符, 重写 \(s_i = b\)\(s_i\ge b\ \)\(\ s_i\le b\).

转换过后, 原公式和现公式是等价的.

比如, 考虑如下约束合取式子: \[ \begin{aligned} x+y\ge 2 &\ \land \\ 2x-y \ge 0 &\ \land \\ -x+ 2y \ge 1 &\ . \end{aligned} \] 该问题可以改写为如下一般形式: \[ \begin{aligned} x+y-s_1=0&\ \land\\ 2x-y-s_2=0&\ \land\\ -x+2y-s_3=0&\ \land\\ s_1\ge 2&\ \land\\ s_2\ge 0&\ \land\\ s_3\ge 1&\ . \end{aligned} \] 新引入的变量 \(s_1,...s_m\) 称为附加变量(additional variables). 原约束中的变量 \(x_1,...,x_n\) 称为问题变量(problem variables). 因此, 转换之后的约束由 \(n\) 个问题变量和 \(m\) 个附加变量组成.

2.2 Simplex算法基础

通常将线性约束满足问题视为几何问题. 在几何体系下, 每一个变量对应一个维度, 每一个约束就定义了一个凸子空间(convex subspace): 特别地, 不等式定义半空间(half-spaces), 等式定义超平面(hyperplanes). 满足赋值的(闭合)子空间由半空间和超平面的交点定义, 并形成凸多面体. 凸子空间之间的交点也是凸的这一事实暗示了这一点. 同样地, 用一个 \(m\times(n+m)\) 的矩阵 \(A\) 来代表输入问题的系数. 变量 \(x_1,...,x_n,s_1,...,s_m\) 重写为向量 x. 据此, 上述问题的等价形式为存在一个向量 x 使得 \[ A{\rm x} = 0\ {\rm and}\ \bigwedge_{i=1}^{m}l_i\le s_i\le u_i, \] 其中, \(l_i\in\{-\infty\}\bigcup\mathbb{Q}\) 是变量 \(x_i\) 的下界, \(u_i\in\{+\infty\}\bigcup\mathbb{Q}\) 是变量 \(x_i\) 的上界. 继续上面的例子, 将合取公式按 \(x,y,s_1,s_2,s_3\) 写成系数矩阵的形式 \[ \begin{pmatrix} 1& 1& -1& 0& 0\\ 2& -1& 0& -1& 0\\ -1& 2& 0& 0& -1 \end{pmatrix}. \]\(m\) 个变量的集合对应的是列, 称为基本变量(basic variables), 记为 \(\mathcal{B}\), 也称为因变量, 因为这些变量的值由哪些非基本变量(nonbasic variables)决定. 非基本变量记为 \(\mathcal{N}\). 方便存储和操作矩阵 \(A\) 的一种表现形式称作 tableau, 这是简化了矩阵 \(A\), 使其没有了对角子矩阵. 因此, tableau 是一个 \(m\times n\) 的矩阵, 其中, 每一列对应的是非基本变量, 每一行与一个基本变量相关联. 相同的基本变量在矩阵 \(A\) 中对角子矩阵的该行上具有"-1"条目. 因此, 原来存储在对角矩阵中的信息现在由标记行的变量表示.

继续我们的例子, 转换后的 tableau 和范围约束如下

\(s\) \(y\)
\(s_1\) 1 1
\(s_2\) 2 -1
\(s_3\) -1 2

\[ \begin{aligned} 2 \le s_1\\ 0 \le s_2\\ 1 \le s_3 \end{aligned} \]

tableau 简化了矩阵 \(A\) 的表达形式, 因此 \(A\)x=0 可以重写成 \[ \bigwedge_{x_i\in\mathcal{B}}(x_i = \Sigma_{x_j\in\mathcal{N}}a_{ij}x_j). \]

2.3 带有上下界的 Simplex 算法

通用 Simplex 算法除了 tableau 之外, 还维护一组赋值 \(\alpha: \mathcal{B\bigcup Q}\rightarrow\mathbb{Q}\). 算法按如下步骤初始化数据结构:

  • 基本变量集合 \(\mathcal{B}\) 就是附加变量的集合.
  • 非基本变量结合 \(\mathcal{N}\) 就是问题变量的集合.
  • 对于任意的 \(x_i, i\in{1,...,n+m}, \alpha(x_i) = 0\).
如果给所有变量的初始 0 赋值满足所有基本变量的上下界, 那么可以说该公式是可满足的. 否则, 该算法开始修改赋值的处理程序. $$ \[\begin{aligned} &\mathbf{Algorithm\ 1.}\ \rm{GENERAL\mbox{-}SIMPLEX}\\ &\mathbf{Input:}\ \rm{A\ linear\ system\ of\ constraints}\ \mathcal{S}\\ &\mathbf{Output:}\ \rm{''Satisfiable''\ if\ the\ system\ is\ satisfiable\ and\ ''Unsatisfiable''\ otherwise}\\ &\\ &1.\hspace{2mm}\rm{Transform\ the\ system\ into\ the\ general\ form}\\ &\hspace{49mm}A\mathbf{x}=0\hspace{5mm} {\rm and}\hspace{5mm}\bigwedge^m_{i=1}l_i\le s_i\le u_i \\ &2.\hspace{2mm}{\rm Set\ }\mathcal{B}\ {\rm to\ be\ the\ set\ of\ additional\ variables\ }s_1,...,s_m.\\ &3.\hspace{2mm}{\rm Construct\ the\ tableau\ for\ }A.\\ &4.\hspace{2mm}{\rm Determine\ a\ fixed\ order\ on\ the\ variables}.\\ &5.\hspace{2mm}{\rm If\ there\ is\ no\ basic\ variable\ that\ violates\ its\ bounds,\ return\ ''Satisfiable''.\ }\\ &\hspace{5.5mm}{\rm Otherwise,\ let\ }x_i\ {\rm be the first\ basic\ variable\ in\ the\ order\ that\ violates\ its\ bounds}.\\ &6.\hspace{2mm}{\rm Search\ for\ the\ first\ suitable\ nonbasic\ variable\ }x_j\ {\rm in\ the\ order\ for\ pivoting\ with\ }x_i.\\ &\hspace{5.5mm}{\rm If\ there\ is\ no\ such\ variable,\ return\ ''Unsatisfiable''}.\\ &7.\hspace{2mm}{\rm Perform\ the\ pivot\ operation\ on\ }x_i\ {\rm and}\ x_j.\\ &8.\hspace{2mm}{\rm Go\ to\ step\ 5}. \end{aligned}\]

$$ 算法1 总结了通用 Simplex 算法的过程, 该算法包含两个不变式:

  • In-1. \(A\)x=0

  • In-2. 非基础变量的值总是在界定的范围内: \[ \forall x_j\in \mathcal{N}.l_j\le \alpha(x_j)\le u_j. \]

显然, 这些不变式最初是成立的, 因为向量 x 中的所有变量都设置为 0, 并且非基本变量没有界定上下界.

算法的主循环检查是否存在任一基本变量超过所规定的阈值. 如果没有这样一个变量, 那么基本变量和非基本变量同时满足其自身的界限. 由于不变式 In-1, 这意味着当前的赋值 \(\alpha\) 是满足的, 算法返回"Satisfiable".

否则, 让变量的 \(x_i\) 作为超出阈值的基本变量, 假定, 不失一般性地, \(\alpha(x_i)>u_i\), 也就是说 \(x_i\) 的赋值超过了规定的上界. 那我们如何改变对于 \(x_i\) 的赋值使其满足阈值的要求? 我们需要找到减少 \(x_i\) 值的方法. \(x_i\) 的值是这样指定的: \[ x_i = \Sigma_{x_j\in \mathcal{N}} a_{ij}x_j. \] 可以通过减小非基本变量 \(x_j\) 的值来减小 \(x_i\) 的值, 此时 \(a_{ij}> 0\), 并且其当前赋值大于其下界 \(l_j\), 或者通过增加非基本变量 \(x_j\) 的值来减小 \(x_i\) 的值, 此时 \(a_{ij}<0\), 并且其当前赋值小于其上界 \(u_j\). 任一变量 \(x_j\) 满足其中一个条件就称为合适(suitable). 如果不存在合适的变量, 那么该问题就是不可满足的并且算法终止.

\(\theta\) 表示增加或者减少 \(\alpha(x_j)\) 的大小, 以符合 \(x_i\) 的上界要求: \[ \theta \doteq \frac{u_i - \alpha(x_i)}{a_{ij}}. \]\(x_j\) 增加(或减少) \(\theta\), 将 \(x_i\) 置于其自身的值域内. 另一方面, \(x_j\) 不必再满足其值域, 因此可能违反不变式 In-2. 因此, 我们在 tableau 中交换 \(x_i\)\(x_j\), 使 \(x_i\) 为非基本变量, 而 \(x_j\) 为基本变量. 这需要对 tabeau 进行转换, 这称为变基操作. 重复执行变基操作, 直到找到可满足的赋值, 或者判定系统不可满足.

2.4 变基操作

在进行变量 \(x_i\)\(x_j\) 的变换之前, 需要引入如下定义:

给定两个变量 \(x_i\)\(x_j\), 系数 \(a_{ij}\) 称为轴元素(pivot element). 变量 \(x_j\) 所在的列称为轴列(pivot column). 行 \(i\) 称为轴行(pivot row).

变换两个变量的前提条件是它们的枢轴元素即系数不能是 0, 即 \(a_{ij}\ne0\). 变基操作如下:

  1. 首先解决变量 \(x_j\) 对应的行 \(i\).
  2. 对于所有的不等于 \(i\) 的行 \(l\), 使用从行 \(i\) 取得的关于 \(x_j\) 的等式来消除 \(x_j\).

这个过程其实就是高斯消元法(Gaussian variable elimination).

比如说, 对于上述的例子. 初始化赋值 \(\alpha(x_i)=0\). 针对上文的 tableau 和约束条件.显然, \(s_1\) 是大于等于 2 的, 而初始赋值为 0, 所以不符合约束条件. 按顺序来说, 最低的非基本变量是 \(x\). 变量 \(x\) 的系数是正值, 没有上限约束, 所以可以作为变基操作的变量. 需要对变量 \(s_1\) 进行加 2 操作, 使得其能够满足大于等于 2 的约束条件, 也就意味着变量 \(x\) 同时也要自增 2(\(\theta=2\)). 变基操作的第一步是对于 \(s_1\) 所在行对 \(x\) 进行变换操作: \[ s_1=x+y\Longleftrightarrow x=s_1-y. \] 该等式带入到其他行的等式中去: \[ \begin{aligned} s_2=2(s_1-y)-y&\Longleftrightarrow s_2=2s_1-3y\\ s_3=-(s_1-y)+2y&\Longleftrightarrow s_3=-s_1+3y \end{aligned} \] 写成表格形式:

\(s_1\) \(y\)
\(x\) 1 -1
\(s_2\) 2 -3
\(s_3\) -1 3

\[ \begin{aligned} \alpha(x)&=2\\ \alpha(y)&=0\\ \alpha(s_1)&=2\\ \alpha(s_2)&=4\\ \alpha(s_3)&=-2\\ \end{aligned} \]

第一步变基操作之后, 可以看到变量 \(s_3\) 也不满足其约束条件; 这是下一步要进行变基操作的变量. 结合上文, 适合参与变换的是变量 \(y\). 需要在 \(s_3\) 的基础上加上 3 才能满足其下界的要求. 也就是说, \[ \theta=\frac{1-(-2)}{3}=1. \] 在对变量 \(s_3\)\(y\) 进行变换操作之后, 最终的 tableau 为

\(s_1\) \(s_3\)
\(x\) 2/3 -1/3
\(s_2\) 1 -1
\(y\) 1/3 1/3

对应的赋值为 \[ \begin{aligned} \alpha(x)&=1\\ \alpha(y)&=1\\ \alpha(s_1)&=2\\ \alpha(s_2)&=1\\ \alpha(s_3)&=1\\ \end{aligned} \] 该组赋值 \(\alpha\) 满足约束条件, 因此 \(\{x\mapsto 1, y\mapsto1\}\) 是一组可满足的赋值.


线性演算(一)
https://socod.github.io/2022/02/39dd3a3dcdf6/
作者
socod
发布于
2022年2月13日
许可协议