形式化10：线性算术理论

算术理论"/>

形式化10：线性算术理论

线性算术理论

语法

A-原子要素 E-表达式 R-关系式 P-命题

A : : = x ∣ c ∣ c × x ( 变量集合 常量集合 常量 × 变量 ) E : : = A ∣ A + E ( 减法可以转换为加法 ) R : : = E = E ∣ E < = E ∣ E < E P : : = R ∣ P ∧ P \begin{align} A::=&x|c|c\times x\ (变量集合\quad 常量集合 \quad常量\times变量)\\ E::=&A|A+E\ (减法可以转换为加法) \\ R::=&E=E|E<=E|E<E \\ P::=&R \\ |&P \land P \end{align} A::=E::=R::=P::=∣​x∣c∣c×x (变量集合常量集合常量×变量)A∣A+E (减法可以转换为加法)E=E∣E<=E∣E<ERP∧P​​

例如（示例中的{相当于/\）

{ x + y = 0.8 x − y = 0.2 ⟹ { x + y = 0.8 x + ( − 1 ) ∗ y = 0.2 \begin{cases} x+y=0.8 \\ x-y=0.2 \end{cases} \Longrightarrow \begin{cases} x+y=0.8 \\ x+(-1)*y=0.2 \end{cases} {x+y=0.8x−y=0.2​⟹{x+y=0.8x+(−1)∗y=0.2​

# Z3 code: 
x, y = Reals(‘x y’)
solve(x+y == 0.8, x-y == 0.2)

解决线性算术理论的算法

傅里叶-莫茨金消元法（Fourier-Motzkin variable elimination）

1. 将线性算术进行标准化，寻找某一个变量 x i x_i xi​，使其仅出现正负值，且不等式右边都为 ≥ 0 \ge 0 ≥0的形式

{ x i + P 1 ( x ) ≥ 0 ( p 1 ) . . . . . . x i + P m ( x ) ≥ 0 ( p m ) − x i + Q 1 ( x ) ≥ 0 ( q 1 ) . . . . . . − x i + Q n ( x ) ≥ 0 ( q n ) R ( x ) ≥ 0 ( r 1 ) \begin{cases} x_i+P_1(x) &\ge 0 \quad &(p_1)\\ ......\\ x_i+P_m(x) &\ge 0 \quad &(p_m)\\ -x_i+Q_1(x) &\ge 0 \quad &(q_1)\\ ......\\ -x_i+Q_n(x) &\ge 0 \quad &(q_n)\\ R(x) &\ge 0 \quad &(r_1) \end{cases} ⎩ ⎨ ⎧​xi​+P1​(x)......xi​+Pm​(x)−xi​+Q1​(x)......−xi​+Qn​(x)R(x)​≥0≥0≥0≥0≥0​(p1​)(pm​)(q1​)(qn​)(r1​)​

其中 P ? ( x ) , Q ? ( x ) , R ( x ) P_?(x),Q_?(x),R(x) P?​(x),Q?​(x),R(x)中都不含有变量 x i x_i xi​。

2. 进行消元，将公式 p 1 p_1 p1​与 q 1 q_1 q1​到 q n q_n qn​分别相加，消去 x i x_i xi​，可得到n个不等式
3. 将 p 2 p_2 p2​到 p m p_m pm​执行和第二步相同的操作，共可得到 m × n m\times n m×n个不等式

{ P 1 ( x ) + Q 1 ( x ) ≥ 0 . . . . . . P 1 ( x ) + Q n ( x ) ≥ 0 P 2 ( x ) + Q 1 ( x ) ≥ 0 . . . . . . P 2 ( x ) + Q n ( x ) ≥ 0 . . . . . . P m ( x ) + Q 1 ( x ) ≥ 0 . . . . . . P m ( x ) + Q n ( x ) ≥ 0 R ( x ) ≥ 0 \begin{cases} P_1(x) + Q_1(x) &\ge 0\\ ......\\ P_1(x) + Q_n(x) &\ge 0\\ P_2(x) + Q_1(x) &\ge 0\\ ......\\ P_2(x) + Q_n(x) &\ge 0\\ ......\\ P_m(x) + Q_1(x) &\ge 0\\ ......\\ P_m(x) + Q_n(x) &\ge 0\\ R(x) &\ge 0 \end{cases} ⎩ ⎨ ⎧​P1​(x)+Q1​(x)......P1​(x)+Qn​(x)P2​(x)+Q1​(x)......P2​(x)+Qn​(x)......Pm​(x)+Q1​(x)......Pm​(x)+Qn​(x)R(x)​≥0≥0≥0≥0≥0≥0≥0​

4. 回到第一步，继续消除下一个变量
5. 如果线性不等式中的所有变量都被消除，那么我们会得到一个常不等式。如果常不等式为真，那么原不等式有解，否则无解

例如

此算法时间复杂度很高

单纯形法（Simplex）

1. 将不等式转化为标准范式
   { P 1 ( x ) ≥ c 1 . . . . . . P n ( x ) ≥ c n ⇒ 令 P i ( x ) = s i { P 1 ( x ) − s 1 = 0 . . . . . . P n ( x ) − s n = 0 s 1 ≥ c 1 . . . . . . s n ≥ c n \begin{cases} P_1(x) \ge c_1 \\ ...... \\ P_n(x) \ge c_n \end{cases} \xRightarrow{令P_i(x)=s_i} \begin{cases} P_1(x) - s_1 = 0\\ ...... \\ P_n(x) - s_n = 0 \\ s_1 \ge c_1 \\ ......\\ s_n \ge c_n \end{cases} ⎩ ⎨ ⎧​P1​(x)≥c1​......Pn​(x)≥cn​​令Pi​(x)=si​ ​⎩ ⎨ ⎧​P1​(x)−s1​=0......Pn​(x)−sn​=0s1​≥c1​......sn​≥cn​​

其中x称为基本变量，s称为额外变量。
例如：
{ x + y ≥ 2 2 x − y ≥ 0 − x + 2 y ≥ 1 → { x + y − s 1 = 0 2 x − y − s 2 = 0 − x + 2 y − s 3 = 0 s 1 ≥ 2 s 2 ≥ 0 s 3 ≥ 1 \begin{cases} x+y \ge 2 \\ 2x-y \ge 0 \\ -x+2y \ge 1 \end{cases} \to \begin{cases} x+y-s1=0 \\ 2x-y-s2=0 \\ -x+2y-s3=0 \\ s1 \ge 2 \\ s2 \ge 0 \\ s3 \ge 1 \end{cases} ⎩ ⎨ ⎧​x+y≥22x−y≥0−x+2y≥1​→⎩ ⎨ ⎧​x+y−s1=02x−y−s2=0−x+2y−s3=0s1≥2s2≥0s3≥1​

2. 构建单纯形表

	x	y
s1	1	1
s2	2	-1
s1	-1	2

其中需满足：

{ s 1 ≥ 2 s 2 ≥ 0 s 3 ≥ 1 \begin{cases} s1 \ge 2 \\ s2 \ge 0 \\ s3 \ge 1 \end{cases} ⎩ ⎨ ⎧​s1≥2s2≥0s3≥1​

3. 尝试并修正（Trial and Fix）
   对x和y选取任意值进行尝试，然后通过单纯形表的变换进行修正，具体如下：

   • 令 x = y = 0 x=y=0 x=y=0，可得：
     x = 0 , y = 0 → s 1 = 0 , s 2 = 0 , s 3 = 0 x=0,y=0 \to \textcolor{red}{s1=0},s2=0,\textcolor{red}{s3=0} x=0,y=0→s1=0,s2=0,s3=0
     其中明显的红色部分不满足要求，需要对其进行修正
   • 首先修正 s 1 s1 s1：通过移项将某个基础变量由 s 1 s1 s1进行表示
     s 1 = x + y → x = s 1 − y s1=x+y \to x=s1-y s1=x+y→x=s1−y
     然后将 s 2 , s 3 s2,s3 s2,s3中的 x x x替换为 s 1 s1 s1

   s 2 = 2 x − y = 2 ( s 1 − y ) − y = 2 s 1 − 3 y s 3 = − x + 2 y = − ( s 1 − y ) + 2 y = − s 1 + 3 y \begin{align} &s2=2x-y=2(s1-y)-y=2s1-3y\\ &s3=-x+2y=-(s1-y)+2y=-s1+3y \end{align} ​s2=2x−y=2(s1−y)−y=2s1−3ys3=−x+2y=−(s1−y)+2y=−s1+3y​​

   • 更新单纯形表

   	s1	y
   x	1	-1
   s2	2	-3
   s3	-1	3

   • 继续尝试，由于此次 s 1 s1 s1有明确范围，可以将 s 1 s1 s1恰好取在边界处（注意，不再是假设 x , y x,y x,y的值，而是 s 1 , y s1,y s1,y）
     s 1 = 2 , y = 0 → x = 2 , s 2 = 4 , s 3 = − 2 s1=2,y=0 \to x=2,s2=4,\textcolor{red}{s3=-2} s1=2,y=0→x=2,s2=4,s3=−2
     其中明显的红色部分不满足要求，需要对其进行修正（ x x x取任意值都可以）
   • 修正 s 3 s3 s3：通过移项将单纯形表上方的某个基础变量由 s 3 s3 s3进行表示（显然只有 y y y）
     s 3 = − s 1 + 3 y → y = s 1 + s 3 3 s3=-s1+3y \to y=\frac{s1+s3}{3} s3=−s1+3y→y=3s1+s3​
     然后将 x , s 2 x,s2 x,s2中的 y y y替换为 s 3 s3 s3

   x = s 1 − y = s 1 − s 1 + s 3 3 = 2 s 1 − s 3 3 s 2 = 2 s 1 − 3 y = 2 s 1 − 3 ( s 1 + s 3 3 ) = s 1 − s 3 \begin{align} &x=s1-y=s1-\frac{s1+s3}{3}=\frac{2s1-s3}{3} \\ &s2=2s1-3y=2s1-3(\frac{s1+s3}{3})=s1-s3 \end{align} ​x=s1−y=s1−3s1+s3​=32s1−s3​s2=2s1−3y=2s1−3(3s1+s3​)=s1−s3​​

   • 更新单纯形表

   	s1	s2
   x	2 3 \frac{2}{3} 32​	− 1 3 -\frac{1}{3} −31​
   s2	1 1 1	− 1 -1 −1
   s3	1 3 \frac{1}{3} 31​	1 3 \frac{1}{3} 31​

   • 继续尝试，由于 s 1 , s 3 s1,s3 s1,s3都有明确范围，因此直接取边界值
     s 1 = 2 , s 3 = 1 → x = 1 , s 2 = 1 , y = 1 s1=2,s3=1 \to x=1,s2=1,y=1 s1=2,s3=1→x=1,s2=1,y=1
     条件全部满足，这样就得到了一组可行的模型（model） [ x = 1 , y = 1 ] [x=1,y=1] [x=1,y=1]

分支定界法（Branch and Bound）

分支定界法应用在ILP问题中，类似于LP（Linear Problem）但是将问题限制在了整数域中。例如

S = { x + y ≥ 2 2 x − y ≥ 0 − x + 2 y ≥ 1 ( x , y ∈ Z ) S= \begin{cases} x+y &\ge 2 \\ 2x-y &\ge 0 \\ -x+2y &\ge 1 \end{cases} (x,y \in \mathbb{Z}) S=⎩ ⎨ ⎧​x+y2x−y−x+2y​≥2≥0≥1​(x,y∈Z)

问题的求解可以利用分治法：

1. 首先解决 x , y ∈ R x,y \in \mathbb{R} x,y∈R上的解(利用单纯形法)

2. 如果无解，那么此问题UNSAT

3. 如果有解 [ x = r 0 , y = r 1 ] [x=r0,y=r1] [x=r0,y=r1]，

   • 如果 r 0 , r 1 ∈ Z r0,r1 \in \mathbb{Z} r0,r1∈Z，则SAT
   • 否则建立分支： ( 1 ) S ∪ [ x ≥ ⌈ c 0 ⌉ ] ( 2 ) S ∪ [ x ≤ ⌊ c 0 ⌋ ] (1)\ S\cup [x \ge \lceil c0 \rceil] \quad (2)S\cup [x \le \lfloor c0 \rfloor] (1) S∪[x≥⌈c0⌉](2)S∪[x≤⌊c0⌋]
     例如 [ x = 1.7 , y = 3.5 ] [x=1.7, y=3.5] [x=1.7,y=3.5]，那么分支后变为

   { x + y ≥ 2 2 x − y ≥ 0 − x + 2 y ≥ 1 x ≥ 2 和 { x + y ≥ 2 2 x − y ≥ 0 − x + 2 y ≥ 1 x ≤ 1 \begin{cases} x+y &\ge 2 \\ 2x-y &\ge 0 \\ -x+2y &\ge 1 \\ x &\ge 2 \end{cases} \quad 和 \quad \begin{cases} x+y &\ge 2 \\ 2x-y &\ge 0 \\ -x+2y &\ge 1 \\ x &\le 1 \end{cases} ⎩ ⎨ ⎧​x+y2x−y−x+2yx​≥2≥0≥1≥2​和⎩ ⎨ ⎧​x+y2x−y−x+2yx​≥2≥0≥1≤1​

4. 对于以上两个分支递归求解，直到SAT或UNSAT

线性运算的应用

编译器优化

for (int i = 1; i < 10; ++i) {A[j + i] = A[j];
}

在上面的代码中可以看到A[j]是不变的，因此可以采用循环变量外提（Loop invariant hoisting）的操作优化代码的执行效率

int t = A[j];
for (int i = 1; i < 10; ++i) {A[j + i] = t;
}

编译器通常会执行这样的优化，但是在此之前，编译器必须检查 i > = 1 ∧ i < 10 ∧ j + i = j i>=1 \land i<10 \land j+i=j i>=1∧i<10∧j+i=j这个命题是UNSAT的，即保证内存没有交叠。

N皇后问题

1. 对系统建立数学模型，设棋盘为N阶矩阵：int board[n][n]，
   • b o a r d [ i ] [ j ] = 1 board[i][j]=1 board[i][j]=1：这个格子有皇后
   • b o a r d [ i ] [ j ] = 0 board[i][j]=0 board[i][j]=0：这个格子没有皇后
2. 每行存在一个皇后：遍历所有的 i i i（即遍历每一行）需满足 ∑ j b o a r d [ i ] [ j ] = 1 \sum_{\substack{j}}board[i][j]=1 ∑j​​board[i][j]=1
3. 每列存在一个皇后：遍历所有的 j j j（即遍历每一列）需满足 ∑ i b o a r d [ i ] [ j ] = 1 \sum_{\substack{i}}board[i][j]=1 ∑i​​board[i][j]=1
4. 每条对角线最多有一个皇后
5. 每条反对角线最多有一个皇后

用z3实现的代码如下

def n_queen_la(board_size: int):solver = Solver()n = board_size# Each position of the board is represented by a 0-1 integer variable:#   ...    ...    ...    ...#   x_2_0  x_2_1  x_2_2  ...#   x_1_0  x_1_1  x_1_2  ...#   x_0_0  x_0_1  x_0_2  ...#board = [[Int(f"x_{row}_{col}") for col in range(n)] for row in range(n)]# only be 0 or 1 in boardfor row in board:for pos in row:solver.add(Or(pos == 0, pos == 1))# @exercise 11: please fill in the missing code to add# the following constraint into the solver:#   each row has just 1 queen,#   each column has just 1 queen,#   each diagonal has at most 1 queen,#   each anti-diagonal has at most 1 queen.# raise Todo("exercise 11: please fill in the missing code.")for i in range(n):arr = []for j in range(n):arr.append(board[i][j])solver.add(sum(arr) == 1)for j in range(n):arr = []for i in range(n):arr.append(board[i][j])solver.add(sum(arr) == 1)# 主对角线# 从(0,i)出发遍历主对角线for i in range(n):arr = []for j in range(n):if j + i < n:print(f"[{j},{j + i}]")arr.append(board[j][j + i])print()solver.add(sum(arr) <= 1)# 从(i,0)出发遍历主对角线for i in range(1, n):arr = []for j in range(n):if j + i < n:print(f"[{j + i},{j}]")arr.append(board[j + i][j])solver.add(sum(arr) <= 1)# 副对角线# 从(i,0)开始遍历副对角线# 则副对角线的所有点满足：横纵坐标和为ifor i in range(n):arr = []for j in range(i + 1):arr.append(board[i - j][j])solver.add(sum(arr) <= 1)# 从(i,n-1)开始遍历副对角线# 则副对角线的所有点满足：横纵坐标和为i+n-1for i in range(1, n):arr = []for j in range(n):row = i + n - 1 - jif row < n:arr.append(board[row][j])solver.add(sum(arr) <= 1)return solver.check()

子集求和问题

在一个集合 S = { x 1 , x 2 , . . . , x n } S=\{x_1,x_2,...,x_n\} S={x1​,x2​,...,xn​}中，是否存在一个子集 T ⊆ S T \subseteq S T⊆S，使得子集中所有的和满足 ∑ T = 0 \sum T=0 ∑T=0

1. 建立数学模型，设int select[n]
   • select[i]=1： x i + 1 x_{i+1} xi+1​被选中
   • select[i]=0： x i + 1 x_{i+1} xi+1​未被选中
2. 描述约束： ∑ i s e l e c t [ i ] × S [ i ] = 0 \sum_{\substack{i}}select[i]\times S[i]=0 ∑i​​select[i]×S[i]=0

任务调度问题

有n个任务，给定其启动时间和结束时间，如何安排这些任务可以获得最多的任务执行数
S = { s 1 , . . . , s n } F = { f 1 , . . . , f n } \begin{align} S=\{s1,...,sn\} \\ F=\{f1,...,fn\} \end{align} S={s1,...,sn}F={f1,...,fn}​​

1. 数学建模，设int select[n]表示是否选定某个任务, times[i][j]表示任务i是否在时刻j中需要执行
2. 设定约束条件
   • t i m e s [ i ] [ j ] = s e l e c t [ i ] × t i m e s [ i ] [ j ] times[i][j]=select[i] \times times[i][j] times[i][j]=select[i]×times[i][j], 例如如果没有选择任务1, 那么第1行就会被清零
   • ∑ i b o a r d [ i ] [ j ] ≤ 1 \sum_{\substack{i}}board[i][j]\le 1 ∑i​​board[i][j]≤1,即同一时刻最多只能执行一项任务(每一列相加不超过1)
3. 设定目标： m a x ( ∑ i s e l e c t [ i ] ) max(\sum_{\substack{i}}select[i]) max(∑i​​select[i])。

任务分配问题

有n个人p1-pn，有n个任务t1-tn，特定的人pi做特定的任务ti有特定的收益cij，如下表所示

	t1	…	tn
p1	c11	…	c1n
…
pn	cn1	…	cnn

如何分配任务可以取得最大收益？

1. 建模：int assign[n][n]，assign[i][j]表示任务j是否分配给人i
2. 设定约束：显然每一行只有一个1，每一列只有一个1（一个人完成一个任务）
3. 设定目标： m a x ( ∑ i , j p r o f i t [ i ] [ j ] ) max(\sum_{\substack{i,j}}profit[i][j]) max(∑i,j​​profit[i][j])

0-1背包问题

给定 n n n个物品，其重量为 W i Wi Wi，价值为 V i Vi Vi，现有一个最大承重 W W W的背包, 如何进行选择才能使包中的物品价值总和最大?

1. 建模: 设置int select[n]表示某个物品是否被选取
2. 约束: ∑ i W [ i ] ∗ s e l e c t [ i ] ≤ W \sum_{\substack{i}}W[i]*select[i]\le W ∑i​​W[i]∗select[i]≤W
3. 目标: m a x ( ∑ i V [ i ] × s e l e c t [ i ] ) max(\sum_{\substack{i}}V[i] \times select[i]) max(∑i​​V[i]×select[i])

总结

通过先行理论和现成的工具，我们无需设计算法，仅需客观的描述问题的约束就可以得到需要的结果。实际上我们并没有降低算法的复杂性，只是将其丢给了数学运算工具（如z3）去完成复杂的工作！