F 2 \mathbb{F}_2 F2上MQ方程的Wu消元法实现
由于Groebner基方法求解MQ问题的困难性,我们先后又尝试了MQ转SAT方法,MQ转SMT方法来求解,但是求解效果仍然有限.在进一步调研后,我们又设计了基于Wu消元法求解算法;更重要的是,我们明确了接下来的道路:只有坚持走实验算法学的道路才能有希望改进 F 2 \mathbb{F}_2 F2上MQ方程求解这样的NP-问题;
本文源代码:MQ方程的Wu消元法实现
F 2 \mathbb{F}_2 F2上MQ方程解的搜索困难性
我们首先来看看 F 2 \mathbb{F}_2 F2上MQ方程的解分布特点:在图中,我们并没有真实地展示解的相对位置(因为那是在 n n n-维空间里的,无法可视化),我们把 n n n-维空间里的点(在 F 2 n \mathbb{F}^n_2 F2n上)嵌入到了 x x x轴上,因此 n n n-维空间里的相邻点并不会在线上体现出相邻,但并不妨碍我们对解的可满足方程数量的随机分布特点的直观认识,可见要找到 M = 16 M=16 M=16的点(也就是方程的解)是一个困难的搜索问题;
Wu消元法求解 F 2 \mathbb{F}_2 F2上MQ问题
Wu消元法的核心算法是"拟除法",我们假设此时有多项式 f , g ∈ k [ x 1 , . . . , x n ] f, g \in k[x_1,...,x_n] f,g∈k[x1,...,xn],它们都包含变元 x i x_{i} xi,也假定变元 x i x_{i} xi的项序最靠前,是下一个待消元的变元,那么我们通过拟除法 REM ( f , g , x i ) \operatorname{REM}(f,g,x_i) REM(f,g,xi)可以得到拟除结果的商和余多项式 q , r q, r q,r,其中 r r r是不含变元 x i x_i xi的,这样我们就得到了消元结果(当然我们不希望看见的是 r = 0 r=0 r=0);注意,Wu消元法解方程在实际应用中往往比Groebner基方法的效率更高 [1];
我们假设现有多项式方程组 g 1 , . . . , g m ∈ k [ x 1 , . . . , x n ] g_{1},..., g_{m} \in k[x_1,...,x_n] g1,...,gm∈k[x1,...,xn],我们希望通过Wu消元法,得到如下形式的特征列:
f 1 = f 1 ( x 1 ) f 2 = f 2 ( x 1 , x 2 ) ⋮ f n = f n ( x 1 , … , x n ) \begin{aligned} f_{1} &=f_{1}(x_{1}) \\ f_{2} &=f_{2}(x_{1}, x_{2}) \\ & \vdots \\ f_{n} &=f_{n}(x_{1}, \ldots, x_{n}) \end{aligned} f1f2fn=f1(x1)=f2(x1,x2)⋮=fn(x1,…,xn)
我们期望得到 f 1 ≠ 0 , . . . , f n ≠ 0 f_{1} \neq 0,..., f_{n} \neq 0 f1=0,...,fn=0,这样我们就可以通过这个严格的三角列逐个解出 x 1 . . . x n x_{1}...x_{n} x1...xn,但是在现实中,多项式方程组 g 1 , . . . , g m g_{1},..., g_{m} g1,...,gm具有多个特征列,我们大概率不能得到这样的严格三角列(我们的实验也证明了这一点),在大多数情况下,我们只能得到形如如下形式的三角列:
f 1 = f 1 ( x 1 ) = 0 f 2 = f 2 ( x 1 , x 2 ) = 0 ⋮ f n = f n ( x 1 , … , x i ) = 0 ⋮ f n = f n ( x 1 , … , x n − 1 ) ≠ 0 f n = f n ( x 1 , … , x n − 1 , x n ) ≠ 0 \begin{aligned} f_{1} &=f_{1}(x_{1}) = 0 \\ f_{2} &=f_{2}(x_{1}, x_{2}) = 0 \\ & \vdots \\ f_{n} &=f_{n}(x_{1}, \ldots, x_{i}) = 0 \\ & \vdots \\ f_{n} &=f_{n}(x_{1}, \ldots, x_{n-1}) \neq 0 \\ f_{n} &=f_{n}(x_{1}, \ldots, x_{n-1}, x_{n}) \neq 0 \end{aligned} f1f2fnfnfn=f1(x1)=0=f2(x1,x2)=0⋮=fn(x1,…,xi)=0⋮=fn(x1,…,xn−1)=0=fn(x1,…,xn−1,xn)=0
这样的三角列也是原方程组的特征列,但是我们就没法依靠它解出原方程了;我们将在后文里讨论如果去求一个严格的三角列,它的本质是一个搜索问题(我们本来想通过消元法来逃避原来的解空间搜索问题,却陷入了另一个搜索问题,所以直观感觉就是MQ问题的本质就是一个搜索问题);
例.1.(Wu消元法求特征列):考虑 k [ x 1 . . . x 4 ] k[x_1...x_4] k[x1...x4]上的多项式方程组:命 P = { P 1 , P 2 , P 3 } \mathcal{P}=\{P_{1}, P_{2}, P_{3}\} P={P1,P2,P3},其中
P 1 = x 1 x 4 2 + x 4 2 − x 1 x 2 x 4 − x 2 x 4 + x 1 x 2 + 3 x 2 P 2 = x 1 x 4 + x 3 − x 1 x 2 P 3 = x 3 x 4 − 2 x 2 2 − x 1 x 2 − 1 \begin{aligned} &P_{1}=x_{1} x_{4}^{2}+x_{4}^{2}-x_{1} x_{2} x_{4}-x_{2} x_{4}+x_{1} x_{2}+3 x_{2} \\ &P_{2}=x_{1} x_{4}+x_{3}-x_{1} x_{2} \\ &P_{3}=x_{3} x_{4}-2 x_{2}^{2}-x_{1} x_{2}-1 \end{aligned} P1=x1x42+x42−x1x2x4−x2x4+x1x2+3x2P2=x1x4+x3−x1x2P3=x3x4−2x22−x1x2−1
消元过程如下(输入多项式集合 P \mathcal{P} P,输出特征列 C \mathcal{C} C):
P = F 1 = { P 1 , P 2 , P 3 } ⊂ F 2 = { P 1 , ⋯ , P 5 } ⊂ F 3 = { P 1 , ⋯ , P 6 } B 1 = [ P 2 ] B 2 = [ P 4 , P 2 ] B 3 = [ P 6 , P 4 , P 2 ] = C R 1 = { P 4 , P 5 } R 2 = { P 6 } R 3 = ∅ , \begin{aligned} \mathcal{P}=&\mathcal{F}_{1}=\{P_{1}, P_{2}, P_{3}\} \quad \subset \mathcal{F}_{2}=\{P_{1}, \cdots, P_{5}\} \quad \subset \mathcal{F}_{{3}}=\{P_{1}, \cdots, P_{{6}}\} \\ &\mathcal{B}_{1}=[P_{2}] \quad \quad \quad \mathcal{B}_{2}=[P_{4}, P_{2}] \quad \quad \quad \mathcal{B}_{3}=[P_{6}, P_{4}, P_{2}]=\mathcal{C}\\ &\mathcal{R}_{1}=\{P_{4}, P_{5}\} \quad \quad \quad \mathcal{R}_{2}=\{P_{6}\} \quad \quad \quad \mathcal{R}_{3}=\varnothing, \end{aligned} P=F1={P1,P2,P3}⊂F2={P1,⋯,P5}⊂F3={P1,⋯,P6}B1=[P2]B2=[P4,P2]B3=[P6,P4,P2]=CR1={P4,P5}R2={P6}R3=∅,
其中, P 4 = REM ( P 1 , P 2 , x 4 ) P_4 = \operatorname{REM}(P_1,P_2,x_4) P4=REM(P1,P2,x4), P 4 = REM ( P 3 , P 2 , x 4 ) P_4 = \operatorname{REM}(P_3,P_2,x_4) P4=REM(P3,P2,x4),而 P 6 = REM ( P 4 , P 5 , x 3 ) P_6 = \operatorname{REM}(P_4,P_5,x_3) P6=REM(P4,P5,x3);最终我们可以得到三角特征列:
C = [ C 1 , C 2 , C 3 ] = [ x 1 ( 2 x 1 x 2 2 + 2 x 2 2 − 2 x 1 x 2 + x 1 + 1 ) x 1 x 3 2 + x 3 2 − x 1 2 x 2 x 3 − x 1 x 2 x 3 + x 1 3 x 2 + 3 x 1 2 x 2 ] x 1 x 4 + x 3 − x 1 x 2 . \begin{aligned} \mathbb{C} &=[C_{1}, C_{2}, C_{3}] \\ &=[\begin{array}{l} x_{1}(2 x_{1} x_{2}^{2}+2 x_{2}^{2}-2 x_{1} x_{2}+x_{1}+1) \\ x_{1} x_{3}^{2}+x_{3}^{2}-x_{1}^{2} x_{2} x_{3}-x_{1} x_{2} x_{3}+x_{1}^{3} x_{2}+3 x_{1}^{2} x_{2}] \\ x_{1} x_{4}+x_{3}-x_{1} x_{2} \end{array}. \end{aligned} C=[C1,C2,C3]=[x1(2x1x22+2x22−2x1x2+x1+1)x1x32+x32−x12x2x3−x1x2x3+x13x2+3x12x2]x1x4+x3−x1x2.
上例中我们用到的求特征列的算法形式化如下,其中 BasSet ( F , o r d ) \operatorname{BasSet}(\mathcal{F}, ord) BasSet(F,ord)是按照项序 o r d ord ord求方程组 F \mathcal{F} F的基列,其实质就是找出下一个用来做拟除的除式.
正如我们前面所说,如上算法求出的特征列并不一定是严格的三角列,事实上,找到一个严格的三角列是很困难的(虽然它在MQ问题的求解中往往是存在的);那么是否存在一种一定找到严格的三角列的算法呢?显然可以,只需要按照项序
x
1
,
.
.
.
,
x
n
x_1,...,x_n
x1,...,xn,逐个变元做两两拟除消元即可,但是这样的计算复杂度特别高,我们现在假定
F
\mathcal{F}
F是变元数
n
n
n,方程数
m
m
m的方程组,第一轮消元
x
1
x_1
x1后我们得到
C
m
2
C^{2}_m
Cm2个方程,第二轮消元
x
2
x_2
x2后我们得到
C
m
(
m
−
1
)
2
C^{2}_{m(m-1)}
Cm(m−1)2个方程(当然其中会有0方程,但是即使是0方程我们也是完成拟除运算才得到的)…以此类推,对于
n
n
n个变元,消元算法的计算复杂度是
O
(
m
n
)
\mathcal{O}(m^n)
O(mn)的,这简直比暴力搜索解的复杂度
O
(
m
2
n
)
\mathcal{O}(m2^n)
O(m2n)还高;
基于以上分析,我们可以看出,计算复杂度是 O ( m n ) \mathcal{O}(m^n) O(mn)的消元算法的搜索策略是广度优先的,也就是每个变元的特征消元多项式它都要找出来,这其实是没有必要的,因为其实只需要找到对应严格三角列的那个"消元路径"就行了,因此我们设计了深度优先的消元算法;
Wu消元算法运行过程中,按照项序排列的变元对应的消元多项式数量变化示意图,我们从 x 1 x_1 x1开始消元,期待得到非零的 f n f_n fn,这个过程中,一开始多项式数量会急剧膨胀,然后随着变元数量变少,又会很快下降,以至于我们的消元路径常常会在中途某个 x i x_i xi时得到0多项式而终止(如图中路径 w 2 w_2 w2所示)我们期望通过深度优先搜索的方式找到路径 w 1 w_1 w1,得到严格的非零三角特征列:
深度优先的Wu消元算法的实现如下:
while VAR_ELIMINATED != ITEMS_VARS_USED[-1]:
print(' ============ '+ str(VAR_ELIMINATED) +' ================= ');
INDEX_VAR = ITEMS_VARS_USED.index(VAR_ELIMINATED);
print('R_m = '+str(len(CHARASTIC_SETS[VAR_ELIMINATED])));
print('R_m+1 = '+str(len(CHARASTIC_SETS[ITEMS_VARS_USED[INDEX_VAR+1]])));
if RUNNING_TIMES>50000:break;
# --- checking stopping ---
VAR_ELIMINATED_NEXT = ITEMS_VARS_USED[ ITEMS_VARS_USED.index(VAR_ELIMINATED_START)+1 ];
if len(CHARSET_NUM_RECORD[VAR_ELIMINATED_NEXT])>2000:
if CHARSET_NUM_RECORD[VAR_ELIMINATED_NEXT][-2000] == len(CHARASTIC_SETS[VAR_ELIMINATED_NEXT]):
VAR_ELIMINATED_START = ITEMS_VARS_USED[ ITEMS_VARS_USED.index(VAR_ELIMINATED_START) +1];
VAR_ELIMINATED = VAR_ELIMINATED_START;continue;
if len(CHARASTIC_SETS[VAR_ELIMINATED])<2:
VAR_ELIMINATED = VAR_ELIMINATED_START;continue; # back-tracing;
if len(CHARASTIC_SETS[ITEMS_VARS_USED[-4]])>200 and (ITEMS_VARS_USED.index(VAR_ELIMINATED_START) < ITEMS_VARS_USED.index(ITEMS_VARS_USED[-5])):
VAR_ELIMINATED_START = ITEMS_VARS_USED[-5];
VAR_ELIMINATED = VAR_ELIMINATED_START;continue;
# --- choosing polynomials ---
if VAR_ELIMINATED==VAR_ELIMINATED_START:POLYNOMIAL_WUBAS = CHARASTIC_SETS[VAR_ELIMINATED][ np.random.randint( len(CHARASTIC_SETS[VAR_ELIMINATED]) ) ];
POLYNOMIAL_DIV = CHARASTIC_SETS[VAR_ELIMINATED][ np.random.randint( len(CHARASTIC_SETS[VAR_ELIMINATED]) ) ];
while POLYNOMIAL_DIV == POLYNOMIAL_WUBAS:POLYNOMIAL_DIV = CHARASTIC_SETS[VAR_ELIMINATED][ np.random.randint( len(CHARASTIC_SETS[VAR_ELIMINATED]) ) ];
INDEX_i = CHARASTIC_SETS[VAR_ELIMINATED].index(POLYNOMIAL_WUBAS);INDEX_j = CHARASTIC_SETS[VAR_ELIMINATED].index(POLYNOMIAL_DIV);
while (VAR_ELIMINATED,(INDEX_i,INDEX_j)) in VISITED_POINTS:
if VAR_ELIMINATED==VAR_ELIMINATED_START:POLYNOMIAL_WUBAS = CHARASTIC_SETS[VAR_ELIMINATED][ np.random.randint( len(CHARASTIC_SETS[VAR_ELIMINATED]) ) ];
POLYNOMIAL_DIV = CHARASTIC_SETS[VAR_ELIMINATED][ np.random.randint( len(CHARASTIC_SETS[VAR_ELIMINATED]) ) ];
while POLYNOMIAL_DIV == POLYNOMIAL_WUBAS:POLYNOMIAL_DIV = CHARASTIC_SETS[VAR_ELIMINATED][ np.random.randint( len(CHARASTIC_SETS[VAR_ELIMINATED]) ) ];
INDEX_i = CHARASTIC_SETS[VAR_ELIMINATED].index(POLYNOMIAL_WUBAS);INDEX_j = CHARASTIC_SETS[VAR_ELIMINATED].index(POLYNOMIAL_DIV);
RUNNING_TIMES_VIS+=1;
if RUNNING_TIMES_VIS>2000:
NEED_TO_BACKT = True;RUNNING_TIMES_VIS=0;break;
if NEED_TO_BACKT:
NEED_TO_BACKT = False;
VAR_ELIMINATED_START = ITEMS_VARS_USED[ ITEMS_VARS_USED.index(VAR_ELIMINATED_START) +1];
VAR_ELIMINATED = VAR_ELIMINATED_START;continue;
CHARSET_NUM_RECORD[ITEMS_VARS_USED[INDEX_VAR+1]].append( len(CHARASTIC_SETS[ITEMS_VARS_USED[INDEX_VAR+1]]) );
VISITED_POINTS.append( (VAR_ELIMINATED,(INDEX_i,INDEX_j)) );
# --- Doing reduction on x_n ---
POLYNOMIAL_QUO,POLYNOMIAL_REM = pseudo_division(POLYNOMIAL_WUBAS,POLYNOMIAL_DIV,VAR_ELIMINATED);
if POLYNOMIAL_REM==0 or (POLYNOMIAL_REM in CHARASTIC_SETS[ITEMS_VARS_USED[INDEX_VAR+1]]):
VAR_ELIMINATED = VAR_ELIMINATED_START;continue; # back-tracing;
CHARASTIC_SETS[ITEMS_VARS_USED[INDEX_VAR+1]].append(POLYNOMIAL_REM);
VAR_ELIMINATED = ITEMS_VARS_USED[INDEX_VAR+1];
POLYNOMIAL_WUBAS = POLYNOMIAL_REM;
RUNNING_TIMES+=1;print('Running --- '+str(RUNNING_TIMES));
if VAR_ELIMINATED==ITEMS_VARS_USED[-1]:print(POLYNOMIAL_REM);
其实我们如上代码采用的是以深度优先为主,实时探测消元进行位置的思路,当消元过程已经迈过了消元特征方程数量的峰值区域后,方程数量急剧减少,这时就更容易得到 0 0 0多项式,这时我们就需要调整搜索起始位置,适当地进行广度优先的策略,这样更容易得到严格特征列,以下是对 m = 16 , n = 8 m=16,n=8 m=16,n=8的MQ方程组运行的结果,成功解出了方程:
============ t1 =================
R_m = 16
R_m+1 = 0
Running --- 1
============ t2 =================
R_m = 1
R_m+1 = 0
============ t1 =================
R_m = 16
R_m+1 = 1
Running --- 2
============ t2 =================
R_m = 2
R_m+1 = 0
Running --- 3
============ t3 =================
R_m = 1
R_m+1 = 0
============ t1 =================
R_m = 16
R_m+1 = 2
Running --- 4
============ t2 =================
R_m = 3
R_m+1 = 1
Running --- 5
============ t3 =================
R_m = 2
R_m+1 = 0
Running --- 6
============ t4 =================
R_m = 1
R_m+1 = 0
============ t1 =================
R_m = 16
R_m+1 = 3
... ...
... ...
============ t5 =================
R_m = 533
R_m+1 = 41
... ...
============ t5 =================
R_m = 533
R_m+1 = 41
Running --- 13272
============ t6 =================
R_m = 42
R_m+1 = 5
Running --- 13273
============ t7 =================
R_m = 6
R_m+1 = 0
Running --- 13274
f_n = t8 + 1;
其中变元消元项序为 t 1 , . . . , t 8 t_1,...,t_8 t1,...,t8, R m \mathcal{R}_m Rm代表当前变元对应特征多项式 f m f_m fm(非0的)的数量, R m + 1 \mathcal{R}_{m+1} Rm+1代表当前变元对应特征多项式 f m + 1 f_{m+1} fm+1(非0的)的数量;最终得到 f n = t 8 + 1 = 0 f_n = t_8+1=0 fn=t8+1=0,也就解出了 t 8 = 1 t_8 = 1 t8=1,再逐次代入之前的方程就可以解出全部解(根据零点定理这就是原方程组(问题里它是 0 0 0维理想)对应的簇);
总结
Wu消元法实现MQ方程求解的过程中,我们在用拟除这个操作按某个项序进行变元消去时,不总是能得到非 0 0 0的多项式,这时消元就会中断,这样就不能得到一个非 0 0 0的三角列以反向解出方程,因此我们不得不采用搜索问题的处理策略来解决这个问题;
参考文献
[1] Ideals, Varieties, and Algorithms,2004,DA Cox and Little, J. and D O’Shea;
[2] 可满足性问题的算法设计与分析,1997,贺思敏;
更多推荐
【Applied Algebra】GF(2)上的MQ方程求解:Wu消元法实现
发布评论