软件学报  2019, Vol. 30 Issue (7): 1903-1915   PDF    
基于SVM的多项式循环程序秩函数生成
李轶1, 蔡天训2, 樊建峰1,3, 吴文渊1, 冯勇1     
1. 中国科学院 重庆绿色智能技术研究院 自动推理与认知重庆市重点实验室, 重庆 400714;
2. 萨基姆通讯(深圳)有限公司, 广东 深圳 518000;
3. 中国科学院大学 计算机科学与技术学院, 北京 100093
摘要: 程序终止性问题是自动程序验证领域中的一个研究热点.秩函数探测是进行终止性分析的主要方法.针对单重无条件分支的多项式循环程序,将其秩函数计算问题归结为二分类问题,从而可利用支持向量机(SVM)算法来计算程序的秩函数.与基于量词消去技术的秩函数计算方法不同,该方法能在可接受的时间范围内探测到更为复杂的秩函数.
关键词: 程序终止性     SVM     机器学习     秩函数    
SVM-based Method for Detecting Ranking Functions in Polynomial Loop Programs
LI Yi1, CAI Tian-Xun2, FAN Jian-Feng1,3, WU Wen-Yuan1, FENG Yong1     
1. Chongqing Key Laboratory of Automated Reasoning and Cognition, Chongqing Institute of Green and Intelligent Technology, Chinese Academy of Sciences, Chongqing 400714, China;
2. Sagemcom Technologies Inc., Shenzhen 518000, China;
3. School of Computer and Control Engineering, University of Chinese Academy of Sciences, Beijing 100093, China
Foundation item: National Natural Science Foundation of China (61572024, 61103110, 11471307)
Abstract: Synthesizing ranking functions of polynomial loop programs is the dominant method for checking their termination. In this study, the synthesis of ranking functions of a class of polynomial loop program is reduced to the binary problem. The support vector machine (SVM) technique then is applied to solve such the binary problem. This naturally relates detection of ranking functions to SVM. Different from the CAD-based method for synthesizing ranking functions, the proposed method can get more expressive polynomial ranking functions in an acceptable time.
Key words: program termination     SVM     machine learning     ranking functions    

循环程序的终止性分析是程序验证的一个重要研究分支.在现代程序设计中, 几乎所有的程序都会含有循环.然而, 即使是简单的循环程序, 也很容易出错.循环中的很多错误往往需要执行多次或在某些特定情况下才能被发现.因此, 确保循环程序是可终止的是保障软件能够可靠运行的一个必要条件.尽管程序终止性问题早已被证明是一个不可判定问题[1, 2].但是人们发现, 具有某些特征的循环程序的终止性是可判定的.如:Tiwari在2004年证明了一类线性循环程序在实数域上的终止性是可判定[3].这类程序的终止性在文献[4-7]中被重新考虑.此外, Zhan等人在文献[8]中考虑了循环条件为等式的多项式程序的终止性问题, 并给出了该类程序可终止和不可终止的充分判准.在文献[9]中, Zhang等人建立了针对程序终止性验证的高级自动机算法.在循环非终止研究方面, Zhang给出了一种能够检测出简单循环中出现死循环的充分条件[10].针对确定型线性赋值循环程序, Leike等人在文献[11]中给出了GNTA条件去探测这类循环的不可终止性.

秩函数法是证明循环程序可终止性的主要方法.给定一个循环程序, 倘若它的秩函数被找到, 则表明该循环是终止的.为计算秩函数方便起见, 人们往往计算多项式型的秩函数.根据多项式秩函数的次数可以将多项式秩函数分为线性秩函数和非线性多项式秩函数.

在线性秩函数研究方面, 已有大量工作.譬如, 2001年, Colón和Sipma通过凸多面体理论以及Farkas’ Lemma给出了一种计算线性秩函数的方法[12, 13].2004年, Podelski和Rybalcheko提出了一种完备的方法来构造线性约束的单分支循环程序的线性秩函数[14].2012年, Bagnara等人[15]研究了线性程序的线性秩函数的计算复杂度问题.在2013年, 他们提出了最终的线性秩函数(eventual ranking function)这一概念[16], 并通过合成最终的线性秩函数来证明单分支线性约束循环程序的终止性.该方法首先利用一个单调增函数对循环程序的状态空间进行划分, 然后利用Farkas’ Lemma来合成最终的秩函数.文献[17]将最终的线性秩函数进行了推广, 提出了深度为L的最终线性秩函数.2014年, Leike等人针对线性循环程序, 给出了包括多阶段秩函数在内的多个秩函数定义[18]. 2017年, Ben-Amram等人[19]对多阶段秩函数进行了进一步的研究, 给出了一个能够在多项式时间内合成多阶段秩函数的方法.同时, 他们还指出了多阶段秩函数和线性字典序秩函数(lexicographic linear ranking function)、嵌套秩函数(nested ranking function)之间的相互联系.相较于线性秩函数的研究, 非线性多项式秩函数合成方面的工作较少, 主要有以下几篇文献讨论了非线性多项式秩函数的计算.Chen等人[20]利用柱形代数分解(CAD)算法给出了一种合成非线性多项式秩函数的方法.该方法首先设定循环程序的非线性秩函数模板, 然后将秩函数合成问题转化为半代数系统求解问题, 从而可利用符号计算工具DISCOVERER和QEPCAD来求解半代数系统, 最终得出非线性秩函数模板中参数的解.此外, 文献[21, 22]等人通过利用半正定规划SDP工具来计算非线性多项式秩函数.

但程序的秩函数未必是多项式的.也即, 可以构造一个循环程序, 它并没有多项式秩函数但具有其他形式的秩函数[23].因此, 当给定程序没有多项式秩函数时, 我们还不能说它是不可终止的, 因为它有可能有其他形式的秩函数.在本文中, 我们主要研究下面给出的单重无条件分支的循环程序:

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{while }}\;\mathit{\boldsymbol{\widetilde p}}(x) \triangleright 0\;{\rm{ do}}\\ \widetilde P\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\{ x{\rm{ :}} = {\rm{ }}\widetilde F(x)\} {\rm{ }}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{endwhile}} \end{array} $ (1)

其中, $x = {({x_1}, {x_2}, L, {x_n})^T} \in {R^n}$是循环程序的变量, $\mathit{\boldsymbol{\widetilde p}}(x)$为关于x的多项式向量, 即$\mathit{\boldsymbol{\widetilde p}}(x) = ({\widetilde p_1}(x), ..., {\widetilde p_s}(x)), $其中, $ \triangleright = ({ \triangleright _1}, ..., { \triangleright _m}), { \triangleright _i} \in \{ > , \ge \} , \mathit{\boldsymbol{\widetilde p}}(x) \triangleright 0$为循环程序$\widetilde p$的循环条件, $x: = \widetilde F(x)$为赋值语句, 其中, $\widetilde F(x) = (\widetilde f(x), \widetilde f(x), ..., \widetilde f(x)), $${f_j}(x) \in R[x].$既然程序中的所有表达式均是多项式的, 因此, 程序P被称为多项式循环程序.令$S = \{ x:\mathit{\boldsymbol{\widetilde p}}(x) \triangleright 0\} , $我们称程序$\mathit{\boldsymbol{\widetilde p}}$是不可终止的, 如果存在一点${x^*} \in {R^n}, $使得对任意的非负整数k, 有$\mathit{\boldsymbol{\widetilde p}}(x)({x^*}) \in S, $如果这样的点不存在, 则称程序$\widetilde p$在实数域上是可终止的.这里, ${\widetilde F^k} = \widetilde F \circ \widetilde F \circ ... \circ \widetilde F.$一般地, 程序$\widetilde p$中的多项式表达式${\widetilde p_i}$, ${\widetilde f_i}$并非是关于x的齐次多项式.但是, 通过增加新的变量z, 我们可将非齐次多项式${\widetilde p_i}$, ${\widetilde f_i}$变为齐次多项式.令${H_{n, {d_i}}}$为所有n元齐di次实系数多项式的集合, 则易知对任意给定的非齐次多项式程序$\widetilde P$, 其总可被转换为齐次多项式程序P.

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{while}}\;\mathit{\boldsymbol{p}}\left( {x, z} \right) \triangleright 0\;\;{\rm{and}}\;\;z > 0\;\;{\rm{do}}\\ P\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{\{ }}x: = F(x, z), z: = {z^d}\} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{endwhile}} \end{array} $

这里, $\mathit{\boldsymbol{p}}(x, z) = ({p_1}(x, z), ..., {p_s}(x, z)), F(x, z) = ({f_1}(x, z), ..., {f_n}(x, z)), $且对任意的$i \in [1, s], {p_i}(x, z) \in {H_{n + 1, {d_i}}}, $对任意的$j \in [1, n], {f_i}(x, z) \in {H_{n + 1, d}}.$这里, 对任意两个整数$a, b$$a \leqslant b, $$[a, b] = \{ a, a + 1, \ldots , b - 1, b\} , $由文献[3, 22]中的定理可知, 程序$\widetilde P$与程序P在终止性上是等价的.

因此, 不失一般性, 本文主要分析下列齐次多项式循环程序终止性.

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{while}}\;\mathit{\boldsymbol{c}}\left( \mathit{\boldsymbol{x}} \right) \triangleright 0\;\;{\rm{do}}\\ Q\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{\{ }}x: = M(x)\} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{endwhile}} \end{array} $

这里, $x \in {R^n}, C(x) = ({c_1}(x), ..., {c_s}(x)), M(x) = ({m_1}(x), ..., {m_n}(x)), $且对任意的$i \in [1, s], {c_j}(x) \in {H_{n, {d_i}}}, $对任意$j \in [1, n], $${m_j}(x) \in {H_{n, d}}, \triangleright = ({ \triangleright _1}, ..., { \triangleright _s}), { \triangleright _i} \in \{ > , \geqslant \} , $且必须存在某个j, 使得${ \triangleright _j} \triangleq > $.也即, 在程序Q中, 其循环条件中每个多项式都是齐次的且次数未必都相同.此外, 至少有一个不等式是严格的不等式.同时, 程序Q的迭代映射M(x)中的各个分量均是齐次多项式且具有相同次数.既然C(x)的分量均是齐次多项式, 故有:若$\mathit{\boldsymbol{c}}\left( \mathit{\boldsymbol{x}} \right) \triangleright 0$, 则对任意正数$\delta > 0$$\mathit{\boldsymbol{C}}(\delta \mathit{\boldsymbol{x}}) \triangleright 0$.程序Q能被其对应的变迁系统刻画:${\Delta _{\rm{Q}}} \buildrel \Delta \over = \mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \wedge x' = M(x).$其中, x表示当前状态, x'表示下一个状态.不失一般性, 下文中, 我们提到的程序Q均是指其对应的变迁系统.令$\Omega =\{x\in {{\mathbb{R}}^{n}}:\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}})\triangleright 0\}.$根据上述程序Q的定义, 既然$\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0$中至少存在一个严格不等式, 那么原点$O \notin \Omega .$对程序Q而言, 秩函数的存在表明程序Q必定是终止的.

在下文中, 我们将建立新的方法去探测程序Q的秩函数.该方法的主要思路是:将秩函数计算问题归结为二分类问题, 然后利用机器学习中的支持向量机(SVM)算法[24, 25]来合成所需秩函数, 最终验证合成的秩函数是否满足秩函数定义.相较于当前的基于CAD或SDP的多项式秩函数计算方法, 本文的主要贡献在于提出的基于SVM的方法用于计算非多项式型秩函数——代数分式型秩函数, 从而将秩函数的研究从多项式型秩函数扩展到代数分式型秩函数, 扩大了可计算秩函数类型的范围.

本文第1节介绍秩函数、支持向量机等相关概念.第2节建立基于SVM的秩函数探测算法, 并通过实例详细阐述该算法.第3节给出更多实验.第4节总结全文.

1 预备知识 1.1 秩函数

定义1.  记号同上.给定齐次多项式循环程序Q.令ρ(x)为一个n元函数.ρ(x)被称为程序Q的秩函数, 如果下列条件被满足:

(a) $\forall x.\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \Rightarrow \rho (x) \ge d$(d为常数);

(b) $\forall x\forall x'.\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \wedge x' = M(x) \Rightarrow \rho (x) - \rho (x') \ge 1.$

上述条件中, (b)能被简化为$\forall x.\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \Rightarrow \rho (x) - \rho (M(x)) \ge 1.$因此, 等价地, 如果ρ(x)满足下列条件:

(a) $\forall x.\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \Rightarrow \rho (x) \ge d;$

(c) $\forall x.\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \Rightarrow \rho (x) - \rho (M(x)) \ge 1, $

ρ(x)为程序Q的秩函数.

1.2 支持向量机

支持向量机是一种性能良好的有监督的学习算法, 其核心思想是通过构造超平面将数据进行分离.给定一组训练样本集$D = \{ ({\mathit{\boldsymbol{x}}_{\bf{1}}}, {y_{\bf{1}}}), ..., ({\mathit{\boldsymbol{x}}_m}, {y_m})\} , {y_i} \in \{ - 1, 1\} , $其中, x1是样本的属性, y1是对应样本的标签.如果存在某个超平面将不同标签的样本点分开, 我们就称该训练样本集是线性可分的.在训练样本集是线性可分的情况下, 能够将训练样本集分开的超平面有许多, 而我们的目的是寻找到泛化能力最强的那一个超平面.如图 1所示, 存在着多条超平面将图中的样本点分开, 但我们需要找到的是中间那条线, 也就是最“居中”的那一条.

Fig. 1 Support vector and margin 图 1 支持向量与间隔

在样本空间中, 划分超平面可用线性方程aTx+k=0来描述, 其中, a=(a1, …, ad)是超平面的法向量, k是截距项.支持向量机可用于找到最“居中”的那一个超平面.如图 1所示, 支持向量机可以找到划分超平面aTx+k=0, 该超平面将图中的两类样本点划分开, 其中, 落在aTx+k=1和aTx+k=–1上的点叫作支持向量.支持向量机的求解通常使用SMO算法[24]来实现.LIBSVM[25]是一个关于SVM的集成软件包, 由台湾大学的Chih-Wei Hsu、Chih-Chung Chang以及Chih-Jen Lin编写完成.软件的主要部分采用C语言编写完成, 可在MATLAB上调用.本文的主要计算部分即采用该软件包完成.

2 主要结果

本节阐述如何将型如程序Q的秩函数计算问题归结为二分类问题, 以便可以利用支持向量机(SVM)工具来计算得到一个候选秩函数, 最后再对该候选秩函数进行验证, 以检验其是否满足秩函数的两个条件(a)和(c).

一般地, 计算程序Q的秩函数, 需要首先预设一个秩函数模板.秩函数模板形式可以是多种多样的.但为了方便计算, 人们常将秩函数模板设定为多项式模板.但是, 由于程序Q是齐次多项式程序, 即其中的所有表达式都是齐次多项式的, 因此下面的定理告诉我们, 程序Q没有多项式秩函数.

定理2.  记号同上.程序Q没有多项式秩函数.

证明:假设程序Q具有多项式秩函数.那么, 不失一般性, 设多项式$\rho (x) = {\sum _{{C_\alpha }}}{C_\alpha }{x^\alpha }$为程序Q的秩函数模板.这里, α=(α1, …, αn), αi为非负整数, ${x^\alpha } = x_1^{{\alpha _1}} \cdot x_2^{{\alpha _2}} \cdot ... \cdot x_n^{{\alpha _n}}, $Cα为单项式xα前的系数.既然ρ(x)为程序Q秩函数, 根据秩函数定义1, ρ(x)必然满足条件(a)和条件(c):

$ \forall x\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \Rightarrow \rho (x) = {\sum _{{C_\alpha }}}{C_\alpha }{x^\alpha } \ge d $ (2)
$ \forall x.\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \Rightarrow \rho (x) - \rho (M(x)) = {\sum _{{C_\alpha }}}{C_\alpha }{x^\alpha } - {\sum _{{C_\alpha }}}{C_\alpha }M{(x)^\alpha } = {\sum _{{C_\alpha }}}{C_\alpha }({x^\alpha } - M{(x)^\alpha }) \ge 1 $ (3)

这里, $M{(x)^\alpha } = {m_1}{(x)^{{\alpha _1}}} \cdot {m_2}{(x)^{{\alpha _2}}} \cdot ... \cdot {m_n}{(x)^{{\alpha _n}}}.$因为C(x)的所有多项式都是齐次的, 则有

$ \forall x\forall \delta .\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0 \wedge \delta {\rm{ > 0}} \Rightarrow \mathit{\boldsymbol{C}}(\delta \mathit{\boldsymbol{x}}) \triangleright 0 $ (4)

任取x*满足$\mathit{\boldsymbol{C}}(x^*) \triangleright 0$, 由式(4)可知, 对任意的正数d, 均有$\mathit{\boldsymbol{C}}(\delta x^*) \triangleright 0, $则有, 当$\delta \to 0$时,

$ {\sum _{{C_\alpha }}}{C_\alpha }(\delta x{^*\alpha } - M{(\delta x^*)^\alpha }) \to 0 \ngeqslant 1. $

这显然与式(3)矛盾.故程序Q没有多项式秩函数.

上述定理表明, 程序Q不可能存在多项式秩函数, 因此, 在计算其秩函数时, 不能将其秩函数模板设置为多项式形式.下文中, 我们将表明程序Q的秩函数模板可以设置为代数分式型模板.这里, 代数分式是两个代数表达式的商, 如$\frac{{x - y}}{{x{y^4} + {x^2} - 1}}, \frac{{\sqrt {xy} }}{{x + 3y + 2}}$均是代数分式.

2.1 程序Q的秩函数计算归结为二分类问题

既然程序Q已被证明不可能有多项式型的秩函数, 那么本节将讨论如何将程序Q的秩函数计算问题归结为二分类问题, 从而便于利用支持向量机(SVM)算法去探测程序Q的代数分式型秩函数.首先给出一些必要的记号.记$\Upsilon (x) = ({x^{{\alpha _1}}}, {x^{{\alpha _2}}}, ..., {x^{{\alpha _k}}}).$这里, $\Upsilon (x)$中的每个分量xkα为关于x的单项式, 且αk为非负整数向量.并记$|{\alpha _k}|$为单项式xkα的次数.若x=(x1, x2), αk=(1, 2), 则${x^{{\alpha _k}}} = x_1^1x_2^2.$$\Upsilon $可被看作是${\mathbb{R}^n} \to {\mathbb{R}^k}$单项式P映射.不失一般性, 我们可以按次数将Υ中的单项式分类排列, 即

$ \Upsilon (x) = ({x^{{\alpha _{{i_1}}}}}, {x^{{\alpha _{{i_2}}}}}, ..., {x^{{\alpha _{{i_{{j_1}}}}}}}, {x^{{\alpha _{{i_{{j_1} + 1}}}}}}, ..., {x^{{\alpha _{{i_{{j_2}}}}}}}, {x^{{\alpha _{{i_{{j_2} + 1}}}}}}, ..., {x^{{\alpha _{{i_{{j_3}}}}}}}, ..., {x^{{\alpha _{{i_{{j_v} + 1}}}}}}, ..., {x^{{\alpha _{{i_k}}}}}). $

这里, 向量Υ被分割成v+1段.其中, 第1段${\Upsilon _1}(x) = ({x^{{\alpha _{{i_1}}}}}, {x^{{\alpha _{{i_2}}}}}, ..., {x^{{\alpha _{{i_{{j_1}}}}}}})$具有相同的次数, 即$|{\alpha _{{i_1}}}| = ... = |{\alpha _{{i_{{j_1}}}}}|;$第2段${\Upsilon _2}(x) = ({x^{{\alpha _{{i_{{j_1} + 1}}}}}}, ...., {x^{{\alpha _{{i_{{j_2}}}}}}})$具有相同的次数, 即$|{\alpha _{{i_{{j_1} + 1}}}}| = ... = |{\alpha _{{i_{{j_2}}}}}|$, …, 第v+1段${\Upsilon _{v + 1}}(x) = ({x^{{\alpha _{{i_{{j_v} + 1}}}}}}, ..., {x^{{\alpha _{{i_{{j_k}}}}}}})$具有相同的次数, 即$|{\alpha _{{i_{{j_v} + 1}}}}| = ... = |{\alpha _{{i_{{j_k}}}}}|.$$\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} $表示将第v${\Upsilon _v}(x)$中的所有分量各自平方后相加再开根号.令

$ T(x) = \left( {\frac{{{\Upsilon _1}(x)}}{{\sqrt {\sum ({\Upsilon _1}{{(x)}^2})} }}, \frac{{{\Upsilon _2}(x)}}{{\sqrt {\sum ({\Upsilon _2}{{(x)}^2})} }}, ..., \frac{{{\Upsilon _{v + 1}}(x)}}{{\sqrt {\sum ({\Upsilon _{v + 1}}{{(x)}^2})} }}} \right), $

显然, T(x)可被看作是${\mathbb{R}^n} \to {\mathbb{R}^k}$映射.并且, T(x)各个分量的分子分母齐次且次数相同, 故对任意的λ > 0, 有T(x)= T(λx).令

$ U(x) = T(x) - T(M(x)) = \left( \begin{array}{l} \frac{{{\Upsilon _1}(x)}}{{\sqrt {\sum ({\Upsilon _1}{{(x)}^2})} }} - \frac{{{\Upsilon _1}(M(x))}}{{\sqrt {\sum ({\Upsilon _1}{{(M(x))}^2})} }}, \frac{{{\Upsilon _2}(x)}}{{\sqrt {\sum ({\Upsilon _2}{{(x)}^2})} }} - \frac{{{\Upsilon _2}(M(x))}}{{\sqrt {\sum ({\Upsilon _2}{{(M(x))}^2})} }}, ..., \\ \frac{{{\Upsilon _{v + 1}}(x)}}{{\sqrt {\sum ({\Upsilon _{v + 1}}{{(x)}^2})} }} - \frac{{{\Upsilon _{{\rm{v + 1}}}}(M(x))}}{{\sqrt {\sum ({\Upsilon _{v + 1}}{{(M(x))}^2})} }} \end{array} \right). $

同样, U(x)可被看作是${\mathbb{R}^n} \to {\mathbb{R}^k}$映射, 并且由于U(x)的特殊构造型式, 使得对任意的λ > 0, 有U(x)=U(λx), 即在同一条射线上的任意两个非零点x, λx, 它们在函数U(x)的映射下的像是相同的.

下面的命题表明, 若T(x)在Ω上有定义, 那么, 函数T(x)为有界函数.也即, 存在正数δ, 使得对任意的x∈Ω有|T(x)|≤δ.

命题3.  记号同上.若T(x)在Ω上有定义, 那么, T(x)在Ω上有界.

证明:既然T(x)在Ω上有定义, 则对Ω中任意一点x, T(x)中各个分量的分母均不会为0.考虑T(x)中第v段, 即

$ \frac{{{\Upsilon _v}(x)}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }} = \left( {\frac{{{x^{{\alpha _{{i_{{j_{v - 1}} + 1}}}}}}}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}, ..., \frac{{{x^{{\alpha _{{i_{{j_v}}}}}}}}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}} \right). $

这里,

$ \sqrt {\sum ({\Upsilon _v}{{(x)}^2})} = \sqrt {{{({x^{{\alpha _{{i_{{j_{v - 1}} + 1}}}}}})}^2} + ... + {{({x^{{\alpha _{{i_{{j_v}}}}}}})}^2}} . $

显然有,

$ \left| {\frac{{{\Upsilon _v}(x)}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}} \right| = \left| {\left( {\frac{{{x^{{\alpha _{{i_{{j_{v - 1}} + 1}}}}}}}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}, ..., \frac{{{x^{{\alpha _{{i_{{j_v}}}}}}}}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}} \right)} \right| \le \left| {\frac{{{x^{{\alpha _{{i_{{j_{v - 1}} + 1}}}}}}}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}} \right| \\ + ... + \left| {\frac{{{x^{{\alpha _{{i_{{j_v}}}}}}}}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}} \right| \le 1 + ... + 1 \le {j_v} - {j_{v - 1}} < \infty . $

因此, $\frac{{{\Upsilon _v}(x)}}{{\sqrt {\sum ({\Upsilon _v}{{(x)}^2})} }}$在Ω上有界.既然T(x)中的每一段都是有界的, 故T(x)在Ω上有界.

既然T(x), U(x)可被看作是${\mathbb{R}^n} \to {\mathbb{R}^k}$映射, 若它们在Ω上有定义, 则记

$ T(\Omega ) = \{ t \in {\mathbb{R}^k}:t = T(x), x \in \Omega \} , U(\Omega ) = \{ u \in {\mathbb{R}^k}:u = U(x), x \in \Omega \} . $

根据命题3, 若T(x)在Ω上有定义, 那么T(Ω)是${\mathbb{R}^k}$空间中的有界域.倘若存在某$a \in {\mathbb{R}^k}$, 使得连续函数$L(u) = {a^T} \cdot u$U(Ω)上有下界1, 即

$ \forall u \in U(\Omega ) \Rightarrow L(u) = {a^T} \cdot u \ge 1 $ (5)

那么, 根据U(Ω)的定义有,

$ \forall x \in \Omega \Rightarrow {a^T} \cdot (T(x) - T(M(x))) \ge 1 $ (6)

根据上述向量a, 构造T(Ω)上的连续函数$L'(t) = {a^T} \cdot t.$既然T(Ω)有界, 故L'(t)在T(Ω)上必有下确界, 即存在c, 使得$\forall t \in T(\Omega ) \Rightarrow L'(t) = {a^T} \cdot t \geqslant c.$根据T(Ω)的定义, 上式等价于

$ \forall x \in \Omega \Rightarrow {a^T} \cdot T(x) \ge c $ (7)

$\rho (x) = {a^T} \cdot T(x) - c.$不难看出, 上述式(6)对应于秩函数定义中的条件(c), 式(7)对应于条件(a).因此, 根据式(6)和式(7)以及秩函数定义, 易知$\rho (x) = {a^T} \cdot T(x) - c$恰好是程序Q的秩函数, 且是代数分式型的秩函数.从上述分析可以看出, 若使得式(5)或者式(6)成立的向量a被找到, 那么式(7)自然成立, 从而可以构造得到程序Q的秩函数$\rho (x) = {a^T} \cdot T(x) - c.$因此, 计算程序Q的秩函数, 我们可以通过计算可满足式(5)或式(6)的向量a来得到.由此可得下列命题.

命题4.  记号同上.若存在aT使得$\forall u \in U(\Omega ) \Rightarrow L(u) = {a^T} \cdot u \geqslant 1, $那么程序Q有代数分式秩函数.要寻找使得式(5)成立的向量aT, 我们可先对Ω的像集U(Ω)进行分析.

情形(Ⅰ):OU(Ω).不难看出, 倘若${\mathbb{R}^k}$空间中的原点O落在U(Ω)中, 即OU(Ω), 那么, 对任意固定的向量aT, 式(5)都不可能成立.这是因为${a^T} \cdot O \equiv 0 \ngeqslant 1.$这等价于, 倘若存在点$x \in \Omega \subseteq {\mathbb{R}^n}, $使得$T(x) - T(M(x)) = 0.$那么, 对任意固定的向量aT, ${a^T} \cdot (T(x) - T(M(x))) \equiv 0 \ngeqslant 1.$而程序Q被给定时, 像集U(Ω)将依赖于T(x)的选择.秩函数的计算依赖于T(x)的选择.因此, 在选择T(x)时, 需要首先满足下面的限制条件:

$ \forall x \in \Omega \subseteq {\mathbb{R}^n} \Rightarrow T(x) - T(M(x)) \ne 0 $ (C1)

限制条件(C1)成立等价于$O \notin U(\Omega )$.

命题5.  记号同上.假设U(x)在Ω上有定义.倘若${\mathbb{R}^k}$空间中的原点$O \notin U(\Omega ), $那么, 计算使得式(5)成立的向量aT的问题就等价于:在空间${\mathbb{R}^k}$中寻找可将原点$O \in {\mathbb{R}^k}$与像集$U(\Omega ) \in {\mathbb{R}^k}$严格分离的超平面.

证明:若存在aT使得式(5)成立, 则$\forall u \in U(\Omega ) \Rightarrow {a^T} \cdot u - 1 \geqslant 0.$$\pi (u) = {a^T} \cdot u - \frac{1}{2}, $则有

$ \forall u \in U(\Omega ) \Rightarrow {a^T} \cdot u - 1 = \left( {{a^T} \cdot u - \frac{1}{2}} \right) - \frac{1}{2} = \pi (u) - \frac{1}{2} \geqslant 0. $

$\forall u \in U(\Omega ) \Rightarrow \pi (u) \geqslant \frac{1}{2} > 0.$同时, 因为选择的T(x)满足限制条件(C1), 故$O \notin U(\Omega ), $$\pi (O) = 0 - \frac{1}{2} = - \frac{1}{2} < 0.$故超平面$\pi (u) = {a^T} \cdot u - \frac{1}{2}$严格分离原点O与像集U(Ω).反之, 若存在超平面$\pi (u) = {a^T} \cdot u + k$严格分离原点O与像集U(Ω), 则有两种情况:

(ⅰ) $\forall u \in U(\Omega ) \Rightarrow \pi (u) = {a^T} \cdot u + k > 0$$\pi (O) < 0;$

(ⅱ) $\forall u \in U(\Omega ) \Rightarrow \pi (u) = {a^T} \cdot u + k < 0$$\pi (O) > 0.$

下面将证明无论哪种情况发生, 满足式(5)的向量aT都必然存在.不失一般性, 这里仅证明(ⅰ), (ⅱ)的证明是类似的.考虑情形(ⅰ), 即$\forall u \in U(\Omega ) \Rightarrow \pi (u) = {a^T} \cdot u + k > 0$$\pi (O) < 0, $则有$\pi (O) = 0 + k = k < 0.$因此, ${a^T} \cdot u - ( - k) > $$0 \Rightarrow \frac{{{a^T}}}{{ - k}} \cdot u - 1 > 0.$显然, 向量$\frac{{{a^T}}}{{ - k}}$满足式(5).对情形(ⅱ), 同理可得, 向量$\frac{{{a^T}}}{{ - k}}$满足式(5).

命题5将程序Q的秩函数计算问题与两个点集(单点集{O}与像集U(Ω)的隔离问题建立了等价关系.下面的情形(Ⅱ)表明, 即便原点$O \notin U(\Omega ), $能严格分离单点集{O}与像集U(Ω)的超平面也未必存在.

情形(Ⅱ):虽然原点$O \notin U(\Omega ), $但原点O却落在U(Ω)的边界上, 即OU(Ω)之间没有间隙.譬如:$U(\Omega ) = $ $\{ ({u_1}, {u_2}), u_2^2 - {u_1} > 0\} .$显然, 原点$O = (0, 0) \notin U(\Omega )$, 但在边界上.由命题5可知, 既然寻找使得式(5)成立的向量aT等价于寻找可将原点O与像集U(Ω)严格分离的超平面, 那么, 当情形(Ⅱ)发生时, 那样的严格分离超平面显然不可能存在.但是, 倘若U(Ω)是个闭的, 那么, 当$O \notin U(\Omega )$时, 原点OU(Ω)具有间隙, 因而有可能存在超平面将OU(Ω)严格分开.不幸的是, 根据程序Q定义及Ω的定义不等式系统可知, 既然不等式系统$\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0$中至少含有一个严格不等式, 故Ω不是闭的, 因而其在U映射下的像U(Ω)不是闭的.尽管U(Ω)不是闭的, 但O与像集U(Ω)之间的严格分离超平面的探测问题可以放松为去探测$\mathit{\boldsymbol{C}}(\mathit{\boldsymbol{x}}) \triangleright 0$空间中可严格分离原点O与某个包含U(Ω)的闭集的严格分离超平面π(u)的问题.为此, 首先给出一些记号.

$\overline{\Omega }=\{x\in {{\mathbb{R}}^{n}}:\mathrm{C}(\mathrm{x})\ge 0\}, \Omega^*=\overline{\Omega }\backslash \{O\}, $这里, {O}为包含原点O的单点集.显然, Ω是包含Ω的闭集.假设U(x)在Ω*上有定义, 即U(x)中各分量的分母在Ω*上均不为0, 既然$\Omega \subseteq \Omega^*, $那么, U(x)在Ω上也有定义.不难看出, $U(\Omega ) \subseteq U(\Omega^*).$下面我们将证明U*)是有界闭集.

命题6.  记号同上.U*)是${\mathbb{R}^k}$中有界闭集.

证明:既然Ω*在几何上是一个由从原点出发的射线构成的锥, 根据此前分析可知, 对任意的λ > 0, 有U(x)= U(λx).这表明, 同一条射线上的任意两个非零点在映射U(x)下的像是相同的.因此, 既然同一条射线上的任意两个非零点在映射U(x)下的像是相同的, 故可在Ω*中的每条射线上各自寻找代表元.既然$O \notin \Omega^*, $则那样的代表元可以设置为$\frac{x}{{|x|}}, x \in \Omega^*.$该代表元刚好是Ω*中的每条从原点出发的射线与单位球Θ相交的交点, 即$\frac{x}{{|x|}} \in \Omega^* \cap \Theta .$因此, $U(x) = U\left( {\frac{x}{{|x|}}} \right), $故而有$U(\Omega^*) = U(\Omega^* \cap \Theta ).$此外, 因为$\overline \Omega \cap \Theta = (\Omega^* \cup \{ O\} ) \cap \Theta = \Omega^* \cap \Theta .$显然, 因为$\overline \Omega $是闭的, Θ是有界闭的, 故$\overline \Omega \cap \Theta $也是有界闭的.由上式可知$\Omega^* \cap \Theta $也是有界闭的.既然前面已经假设U(x)在Ω*上有定义, 因此U(x)在Ω*上连续.由连续函数性质可知, $\Omega^* \cap \Theta $U(x)的像$U(\Omega^* \cap \Theta )$是有界闭的.

既然$U(\Omega^*) = U(\Omega^* \cap \Theta ), $故有U*)是有界闭的.

如果$O \notin U(\Omega^*), $根据命题6, 既然U*)是有界闭集, 那么在O与像集U*)之间存在间隙, 因而有可能在二者之间存在严格分离超平面.既然$U(\Omega ) \subseteq U(\Omega^*), $$O \notin U(\Omega^*)$时, 如果存在超平面π(u)能够严格分离原点O与像集U*), 那么π(u)也能严格分离原点OU(Ω).类似于命题5的证明, 我们有下列结论.

定理7.  记号同上.假设U(x)在Ω*上有定义且原点$O \notin U(\Omega^*) \subseteq {\mathbb{R}^k}.$如果存在超平面π(u)可将原点O与像集U*)严格分离, 那么一定存在向量aT使得式(5)成立.

证明:如果存在超平面π(u)=aT·u+k可将原点O与像集U*)严格分离, 那么, 既然$U(\Omega ) \subseteq U(\Omega^*), $则超平面π(u)=aT·u+k也将严格地分离原点O与像集U*).根据命题5及其证明可知, 一定存在向量${a^T}: = \frac{{{a^T}}}{{ - k}}$满足式(5), 即$\forall u \in U(\Omega ) \Rightarrow \frac{{{a^T}}}{{ - k}} \cdot u \geqslant 1.$

定理7表明, 当定理的前提假设被满足时, 程序Q的秩函数计算问题可以最终归结为原点OU*)的超平面分离问题.倘若那样的超平面被找到, 那么$\rho (x) = \frac{{{a^T}}}{{ - k}} \cdot T(x)$就恰好是程序Q的代数分式秩函数.而后者显然是一个二分类问题.通常, 二分类问题可以利用支持向量机(SVM)算法来解决.这里, 令A={O}为一单点集, B=U*).我们接下来的问题就是如何利用SVM算法求得一个可将集合AB严格分离的超平面.一般地, 在通常的二分类问题中, 涉及到的两个数据集均是有限点集.而这里的二分类问题中, 集合B显然是一个无穷点集.因此, 为了便于利用SVM算法计算AB的超平面, 我们将采取下面的步骤.

S1:从B集中取出有限个点记为${B_0} \subseteq B.$转S2;

S2:调用SVM算法计算单点集A与有限点集B0的分离超平面, 记为π(u)=aT·u+k.转S3;

S3:验证式(5):$\forall u \in U(\Omega ) \Rightarrow {\pi _0}(u) = \frac{{{a^T}}}{{ - k}} \cdot u \geqslant 1$是否成立.若成立, 则表明程序Q有代数分式型秩函数, 即$\rho (x) = \frac{{{a^T}}}{{ - k}} \cdot T(x).$反之, 从B集中取出新的有限点集${B_1}({B_0} \subset {B_1}).$${B_0}: = {B_1}, $转S2.

上述步骤S1中, 采用打格子的方式进行B0的构造.这一点将在下一小节中加以具体介绍.在S2中, 当调用SVM算法计算单点集A与有限点集B0的分离超平面π0(u)时, 理论上会有失败的可能, 即那样的超平面可能不存在, 但在具体的大量实验中我们发现, 那样的超平面总存在.在S3中, 式(5)的验证可归结为验证:$\forall x \in \Omega \Rightarrow $$\frac{{{a^T}}}{{ - k}} \cdot (T(x) - T(M(x))) \geqslant 1.$同时, 既然T(x)含有根式表达式, 故我们利用强有力不等式验证器BOTTEMA来验证上式是否成立.在S3中, 若划分AB0的超平面经过验证不可能对应一个分式秩函数时, 我们将通过更加细分格子的思路来得到含有更多点的有限点集B1.下面是基于上述3个步骤建立的基于SVM的秩函数探测算法.

2.2 有限点集B0的构造

在上述的S1中, 需要从B=U*)中取出有限个点构成B0.由命题8的证明可知, $U(\Omega^*) = U(\Omega^* \cap \Theta ).$这里, $\Omega^* \cap \Theta $是一个有界闭集, Θ为单位球.因此, 要构造B0, 我们可以先从有界闭集$\Omega^* \cap \Theta $中取出有限点集C, 然后通过U(x)映射得到点集U(C), 从而得到B0=U(C).显然, 直接在球面Θ上通过打格点的方式取点并不方便.但单位球Θ总是可以被包含在某个超立方体△中.故, 为了构造$\Omega^* \cap \Theta $上的点集C, 我们采取如下思路:先在区域$\Omega^* \cap \Delta $上构造点集S.这里, $\Delta = {[ - m, m]^n}$为一超立方体.然后对点集S中的点进行归一化, 得到归一化后的点集S0.显然, ${S^0} \subseteq \Omega^* \cap \Theta .$最后令$C: = {S^o}.$而在超立方体△上打格子是相对容易的.因此, 在区域$\Omega^* \cap \Delta $上构造点集S时我们可以通过对超立方体D打格子, 然后将落入Ω*的格点收集起来构成点集S.图 2所示为一个算法化的取点过程.

Fig. 2 The algorithm of getting the sample points 图 2 取样本点算法

在上述算法中, 需要设置m的值.在我们的实验中, m被取为100.此外, 取点的间距h需要根据计算机的计算能力进行适当的调整, 在计算机的实际计算能力范围内, 间距越小, 取得的样本点越多, 最终计算出的代数分式函数经过验证成为真正秩函数的可信度就越高.

2.3 例子

本节将通过一个例子来阐述上文提出的方法.

例1:考虑下列循环:

$ {Q_1}\;\;{\rm{while}}\;\;4x_1^2 + x_2^2 + 16 \le 16{x_1} + 6{x_2} \wedge {x_1} + 4{x_2} \ge x_2^2 + 5\;\;\;\;{\rm{do}} $
$ {\rm{\{ }}{x'_1}{\rm{ = }}x_1^2 - x_2^2{\rm{;}}{x'_2}{\rm{ = }}{x_1} + {x_2}{\rm{ + 1;\} }} $

该循环若使用RegularChains[26]或者redlog[27]来求解线性或者非线性秩函数, 均会因复杂度太高, 而无法在有效的时间内得出结论.接下来将演示如何通过上文提出的基于SVM的方法来探测循环程序P1的代数分式秩函数.

(a) 首先将循环程序Q1进行齐次化处理, 变成终止性等价的循环程序${Q'_1}{\rm{ , }}$如下所示:

$ {{Q'}_1}\;\;{\rm{while}}\;\;4x_1^2 + x_2^2 + 16x_3^2 \le 16{x_1}{x_3} + 6{x_2}{x_3} \wedge {x_1}{x_3} + 4{x_2}{x_3} \ge x_2^2 + 5x_3^2 \wedge {x_3} > {\rm{0}}\;\;\;{\rm{do}} $
$ {\rm{\{ }}{x'_1}{\rm{ = }}x_1^2 - x_2^2{\rm{;}}{x'_2}{\rm{ = }}{x_1}{x_3} + {x_2}{x_3}{\rm{ + }}x_3^2{\rm{;}}{x'_3}{\rm{ = }}x_3^2;{\rm{\} }} $

其中, $\Omega {\rm{ = \{ }}x:4x_1^2 + x_2^2 + 16x_3^2 \leqslant 16{x_1}{x_3} + 6{x_2}{x_3} \wedge {x_1}{x_3} + 4{x_2}{x_3} \geqslant x_2^2 + 5x_3^2 \wedge {x_3}{\rm{ > 0\} }}$是齐次循环程序${Q'_{\rm{1}}}{\rm{ }}$的循环条件所围成的空间, 它是一个锥, 也即, 针对任意的$x \in \Omega , $都有$\lambda x \in \Omega , \lambda > 0.$

$\bar \Omega {\rm{ = \{ }}x:4x_1^2 + x_2^2 + 16x_3^2 \leqslant 16{x_1}{x_3} + 6{x_2}{x_3} \wedge {x_1}{x_3} + 4{x_2}{x_3} \geqslant x_2^2 + 5x_3^2 \wedge {x_3}{\rm{ > 0\} }}{\rm{.}}$$\Omega^* {\rm{ = }}\bar \Omega \backslash \{ O\} .$这里, O为原点(0, 0, 0).对该例, 我们可令$T(x) = \left( {\frac{{{x_1}}}{{|x|}}, \frac{{{x_2}}}{{|x|}}, \frac{{{x_3}}}{{|x|}}} \right).$T(x)的每一个分量的奇点均为原点O, 而Ω*中并不包含原点, 故T(x)在Ω*上连续.这里, $x = ({x_1}, {x_2}, {x_3}, ), |x| = \sqrt {x_1^2 + x_2^2 + x_3^2} .$$T(x') = \left( {\frac{{{{x'}_1}}}{{|x'|}}, \frac{{{{x'}_2}}}{{|x'|}}, \frac{{{{x'}_3}}}{{|x'|}}} \right), $其中, $x' = ({x'_1}, {x'_2}, {x'_3}).$既然$x' = M(x) = $ $(x_1^2 - x_2^2, {x_1}{x_3} + {x_2}{x_3} + x_3^2, x_3^2), $因此,

$ U(x) = T(x) - T(x') = T(x) - T(M(x)) = \left( {\frac{{{x_1}}}{{|x|}} - \frac{{{{x'}_1}}}{{|x'|}}, \frac{{{x_2}}}{{|x|}} - \frac{{{{x'}_2}}}{{|x'|}}, \frac{{{x_3}}}{{|x|}} - \frac{{{{x'}_3}}}{{|x'|}}} \right). $

这里, $|x'| = \sqrt {{{({{x'}_1})}^2} + {{({{x'}_2})}^2} + {{({{x'}_3})}^2}} .$因此, U(x)将原三维空间(原像空间)中的点映射到了另外一个三维空间(像空间).根据定理7, 我们可通过支持向量机SVM的方法在像空间中去探测可将原点O与像点U*)严格分隔开的超平面π(u)=aT·u+k.若那样的超平面被找到, 则$\rho (x) = \frac{{{a^T}}}{{ - k}} \cdot T(x)$就恰好是程序${Q'_{\rm{1}}}{\rm{ }}$的代数分式秩函数.但是, 我们需要事先验证该例是否满足定理7的前提条件, 即判断:U(x)在Ω*上是否有定义且是否原点$O \notin U(\Omega^*).$

(ⅰ)判断U(x)在Ω*上是否有定义.这等价于判断是否存在使得$|x| = 0 \vee |x'| = 0$成立的点会落入到Ω*中.利用Maple中的命令solve不难验证, 使得$|x| = 0 \vee |x'| = 0$成立的点仅有原点O=(0, 0, 0).而Ω*并不包含原点.故U(x)在Ω*上有定义.

(ⅱ)判断是否原点$O \notin U(\Omega^*).$这等价于判断下列系统是否有解:

$ \left. \begin{array}{l} \frac{{{x_1}}}{{|x|}} - \frac{{{{x'}_1}}}{{|x'|}} = 0 \wedge \frac{{{x_2}}}{{|x|}} - \frac{{{{x'}_2}}}{{|x'|}} = 0 \wedge \frac{{{x_3}}}{{|x|}} - \frac{{{{x'}_3}}}{{|x'|}} = 0 \wedge {{x'}_1}{\rm{ = }}x_1^2 - x_2^2 \wedge {{x'}_2} = {x_1}{x_3} + {x_2}{x_3}{\rm{ + }}x_3^2 \wedge {{x'}_3}\\ {\rm{ = }}x_3^2 \wedge 4x_1^2 + x_2^2 + 16x_3^2\\ \le 16{x_1}{x_3} + 6{x_2}{x_3} \wedge {x_1}{x_3} + 4{x_2}{x_3}\\ \ge x_2^2 + 5x_3^2 \wedge {x_3}\\ \ge 0 \wedge ({x_1} \ne 0 \vee {x_2} \ne 0 \vee {x_3} \ne 0) \wedge |x'|\\ {\rm{ = }}\sqrt {{{({{x'}_1})}^2} + {{({{x'}_2})}^2} + {{({{x'}_3})}^2}} \end{array} \right\} $ (8)

若式(8)无解, 则表明$O \notin U(\Omega^*).$由于代数分式以及根式的出现使得我们直接计算上述系统是否有解较为困难, 因此, 我们采取了一些化简措施.比如, 将上述系统中的代数分式化为多项式, 利用平方手段消除根号等.由此我们得到下列多项式系统:

$ \left. \begin{array}{l} x_1^2({({{x'}_1})^2} + {({{x'}_2})^2} + {({{x'}_3})^2}) = {({{x'}_1})^2}(x_1^2 + x_2^2 + x_3^2) \wedge x_2^2({({{x'}_1})^2} + {({{x'}_2})^2} + {({{x'}_3})^2})\\ {\rm{ }} = {({{x'}_2})^2}(x_1^2 + x_2^2 + x_3^2) \wedge x_3^2({({{x'}_1})^2} + {({{x'}_2})^2} + {({{x'}_3})^2})\\ {\rm{ }} = {({{x'}_3})^2}(x_1^2 + x_2^2 + x_3^2) \wedge {{x'}_1} = x_1^2 - x_2^2 \wedge {{x'}_2} = {x_1}{x_3} + {x_2}{x_3} + x_3^2 \wedge {{x'}_3}\\ {\rm{ }} = x_3^2 \wedge 4x_1^2 + x_2^2 + 16x_3^2 \le 16{x_1}{x_3} + 6{x_2}{x_3} \wedge {x_1}{x_3} + 4{x_2}{x_3} \ge x_2^2 + 5x_3^2 \wedge {x_3} \ge 0 \end{array} \right\} $ (9)

显然, 式(8)的任何一个解也必是式(9)的解.调用Maple中的命令solve计算可得式(9)仅有零解(0, 0, 0).显然, $(0, 0, 0) \notin \Omega^*.$因为零解不可能是式(8)的解, 故式(8)无解.这表明原点$O \notin U(\Omega^*).$

综上所述, 该循环程序满足定理7的前提条件.因此, 根据定理7, 秩函数计算问题可归结为OU*)之间严格分离超平面的计算问题.

(b) 取样本点.

采用第2.2节介绍的取点算法, 并将该算法中的参数m=100.在本例中我们令取点间距h=1.由于在Ω约束条件中有x3≥0的条件, 故x3不必从-100开始取点.通过取点算法, 我们得到U*)中的有限点集B0=U(C).将B0中的点的标签置为1, 而令单点集A=(O)的标签为-1.接下来调用SVM算法计算AB0的严格分离超平面π0(u)=aT·u+k.

(c) 使用LIBSVM软件包计算出像空间中的划分超平面.

本文采用MATLAB调用LIBSVM软件包的方式来计算一个“硬分隔”的超平面.令$\overline D = A \cup {B_0}.$在求得划分超平面之后, 需要对D中的点进行一次完整的测试, 测试结果必须是100%完全精确才能进行下一步操作.对该例通过计算得到划分超平面中aTk的值,

$ a = (2.798739925884132, - 6.058645127373750, 6.167554957795089), k = - 1.000000000061875. $

由上述计算得到的aTk的值可确定一个严格分离超平面π0(u)=aT·u+k, 它将A={O}与${B_0}( \subseteq U(\Omega^*))$严格分离.即:

(ⅰ) $\forall u \in {B_0} \Rightarrow {\pi _{\rm{0}}}(u) = {a^T} \cdot u + k > 0$${\pi _{\rm{0}}}(O) < 0;$

(ⅱ) $\forall u \in {B_0} \Rightarrow {\pi _{\rm{0}}}(u) = {a^T} \cdot u + k < 0$${\pi _{\rm{0}}}(O) > 0.$

由类似于命题5的证明可知, 无论哪种情形发生, 都将有

$ \forall u \in {B_0} \Rightarrow \frac{{{a^T}}}{{ - k}} \cdot u - 1 > 0. $

接下来将验证超平面π0(u)是否能将A={O}与B=U*)严格分离.为此, 就需要验证下列两种情形之一是否成立.

(ⅰ) $\forall u \in U(\Omega^*) \Rightarrow {\pi _{\rm{0}}}(u) = {a^T} \cdot u + k > 0$${\pi _{\rm{0}}}(O) < 0;$

(ⅱ) $\forall u \in U(\Omega^*) \Rightarrow {\pi _{\rm{0}}}(u) = {a^T} \cdot u + k < 0$${\pi _{\rm{0}}}(O) > 0.$

既然像点uU*)且由U*)的定义可知$U(\Omega^*) = \{ u \in {\mathbb{R}^k}:u = U(x), x \in \Omega^*\} , $那么验证上述(ⅰ), (ⅱ)两式就等

价于验证:

(ⅰ*) $\forall x \in \Omega^* \Rightarrow \beta (x) = {a^T}[T(x) - T(M(x))] + k > 0$k < 0;

(ⅱ*) $\forall x \in \Omega^* \Rightarrow \beta (x) = {a^T}[T(x) - T(M(x))] + k < 0$k > 0.

这是因为, 显然地, 若(ⅰ)(或(ⅱ))成立, 那么(ⅰ*)(或(ⅱ*))也必成立.反之亦然.

(e) 使用BOTTEMA软件包进行验证.

既然k=–1.000000000061875 < 0, 那么我们需要验证:

$ \forall x \in \Omega^* \Rightarrow \beta (x) = {a^T}[T(x) - T(M(x))] + k > 0 $

是否成立.由于最终需要验证式(5)或式(6)是否成立, 所以, 在本例中, 我们直接验证:

$ \forall x \in \Omega \Rightarrow \beta (x) = {a^T}[T(x) - T(M(x))] + k > 0 $

是否成立.BOTTEMA是一款强有力的不等式证明软件.它可以验证带复杂根式的代数表达式的不等式是否成立.我们可以在Maple上调用BOTTEMA中的各类函数.通过验证发现,

$ \forall x \in \Omega \Rightarrow \beta (x) = \frac{{{a^T}}}{{ - k}}[T(x) - T(M(x))] - 1 > 0 $

的确成立.故,

$ \rho (x) = \frac{{{a^T}}}{{ - k}} \cdot T(x) = \frac{{{{(2.798\;739\;925\;884\;132, - 6.058\;645\;127\;373\;750, 6.167\;554\;957\;795\;089)}^T}}}{{ - {\rm{1}}{\rm{.000}}\;{\rm{000}}\;{\rm{000}}\;{\rm{061}}\;{\rm{875}}}}\left( {\frac{{{x_1}}}{{|x|}}, \frac{{{x_2}}}{{|x|}}, \frac{{{x_3}}}{{|x|}}} \right) $

是该程序的代数分式型秩函数.因此循环程序是可终止的.

3 实验

下面给出基于SVM的多项式循环程序秩函数探测的具体算法.

实验是在一台装有7.86GB的可用内存、处理器频率为2.59GHz、操作系统64位的Windows操作系统的计算机上进行的.下面我们选取了4个循环程序用于对比.每个循环程序具体的计算方法可以参照前文中的例1.

针对表 1中的6个循环程序, 我们先后用量词消去工具Redlog和Maple中的RegularChains软件包来计算各自的线性秩函数, 从表 2可以看出, 除了循环程序P4P5外, 其他5个循环程序都不能通过量词消去方法来证明其是终止的.原因在于, 量词消去算法的复杂度太高, 当循环程序的赋值语句或循环条件中含有非线性项时, 无法在有效的时间内给出计算结果.在表 2中, unknown表示计算时间超过10个小时仍未得出结果.Terminating表示相应的线性秩函数或代数分式秩函数被成功找到, 故表明循环是终止的.特别地, 在Terminating后面的括号中给出了计算秩函数所花费的时间.显然, 当循环中的表达式是线性时, 基于量词消去的工具即能有效地计算它们的线性秩函数.但, 当表达式中含有非线性项时, 实验结果表明, 量词消去工具较难在可接受的时间内计算线性秩函数.此外, 本文提出的基于SVM的方法所涉及的时间开销包括:样本点集的构造所需时间、SVM计算时间以及用BOTTEMA验证代数分式是否为秩函数所需时间.在实验中我们发现, 本文方法最主要的时间开销在于样本点集的构造.譬如, 循环P3, 样本点集的构造花费了76分钟.当然, 我们也可以尝试使用量词消去的方法来合成非线性秩函数, 但是, 非线性项的引入同样会带来高复杂度的问题.从表 2中可以看出, 本文提出的基于SVM的秩函数探测方法能够探测出上述6个循环程序的确具有代数分式型秩函数, 从而证明了这6个程序均是可终止的.上述6个循环程序在使用SVM的方法时, 各自的取点间距h以及所产生的划分超平面aT·u+k中的参数aT和参数k的计算值见表 3.

Fig. 3 Ranking function detection of polynomial loop programs based on SVM 图 3 基于SVM的多项式循环程序秩函数探测

Table 1 More examples 表 1 更多实例

Table 2 Comparison of several tools 表 2 几种工具的比较

Table 3 The parameters generated using the SVM method 表 3 使用SVM的方法所产生的参数

表 3可以看出, 循环程序P3P4所采用的取点间距为0.5, 而循环程序P2采用的取点间距为0.2.从理论上来说, 取点间距越小, 所选取的样本点越多, 最终计算出的代数分式函数被成功验证为秩函数的可能性就越大.但是考虑到计算机的实际计算能力, 针对不同的循环程序, 取点间距要作适当的调整.

4 总结

本文提出了一种基于支持向量机(SVM)理论的求解单重无条件分支多项式循环程序秩函数的方法.该方法将秩函数计算问题归结为二分类问题, 然后利用SVM工具计算出一个候选的代数分式型函数, 最后借助不等式自动验证工具BOTTEMA对其进行验证, 以确定该代数分式型函数是否为循环程序的秩函数.由于SVM方法内部采用数值计算, 因此相较于现有的基于量词消去的秩函数计算方法, 本文方法能在可接受时间内探测形式较为复杂的秩函数.实验结果表明, 针对某些循环程序, 传统的量词消去方法不能在合理时间内判断出它们是否有线性秩函数, 但通过本文提出的方法却可以找到其代数分式型秩函数.

参考文献
[1]
Turing A. On computable numbers, with all application to the entscheidungsproblem. London Mathematical Socity, 1936, 42(2): 230-265. https://www.onacademic.com/detail/journal_1000035257000410_1f50.html
[2]
Yang L, Zhou CC, Zhan NJ, Xia BC. Recent advances in program verification through computer algebra. Frontiers of Computer Science, 2010, 4(1): 1-16. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=zggdxxxswz-jsjkx201001001
[3]
Tiwari A. Termination of linear programs. Computer Aided Verification, 2004, 3114: 70-82. [doi:10.1007/b98490]
[4]
Braverman M. Termination of integer linear programs. In: Ball T, Jones RB, eds. Proc. of the Computer Aided Verification (CAV 2006). Springer-Verlag, 2006. 372-385.
[5]
Xia BC, Yang L, Zhan NJ, Zhang ZH. Symbolic decision procedure for termination of linear programs. Formal Aspects of Computing, 2009, 23(2): 171-190. http://d.old.wanfangdata.com.cn/NSTLQK/10.1007-s00165-009-0144-5/
[6]
Xia BC, Zhang ZH. Termination of linear programs with nonlinear constraints. Journal of Symbolic Computation, 2010, 45(11): 1234-1249. [doi:10.1016/j.jsc.2010.06.006]
[7]
Li Y. Witness to non-termination of linear programs. Theoretical Computer Science, 2017, 681: 75-100. [doi:10.1016/j.tcs.2017.03.036]
[8]
Liu J, Xu M, Zhan NJ, Zhao HJ. Discovering non-terminating inputs for multi-path polynomial programs. Journal of Systems Science and Complexity, 2014, 27: 1284-1304. http://d.old.wanfangdata.com.cn/Periodical/xtkxysx201406013
[9]
Chen YF, Heizmann M, Lengal O, Li Y, Tsai MH, Turrini A, Zhang LJ. Advanced automate-based algorithms for program termination checking. In: Proc. of the ACM SIGPLAN Conf. on Programming Language Design and Implementation (PLDI 2018). ACM, 2018. 135-150.
[10]
Zhang J. A path-based approach to the detection of infinite looping. In: Proc. of the 2nd Asia-Pacific Conf. on Quality Software (APAQS 2001). Hong Kong: IEEE Computer Society, 2001. 88-96.
[11]
Leike J, Heizmann M. Geometric nontermination arguments. In: Proc. of the 24th Int'l Conf. on Tools and Algorithms for the Construction and Analysis of Systems (TACAS 2018). Berlin, Heidelberg: Springer-Verlag, 2018, 10806: 266-283.
[12]
Colón M, Sipma H. Synthesis of linear ranking functions. In: Margaria T, Yi W, eds. In: Proc. of the 7th Int'l Conf. on Tools and Algorithms for Construction and Analysis of Systems (TACAS 2001). Springer-Verlag, 2001. 67-81.
[13]
Podelski A, Rybalchenko A. A complete method for the synthesis of linear ranking functions. In: Proc. of the Verification, Model Checking, and Abstract Interpretation (VMCAI). Berlin: Springer-Verlag, 2004. 239-251.
[14]
Bagnara R, Mesnard F, Pescetti A, Zaffanella E. A new look at the automatic synthesis of linear ranking functions. Information and Computation, 2012, 215: 47-67. [doi:10.1016/j.ic.2012.03.003]
[15]
Bagnara R, Mesnard F. Eventual linear ranking functions. In: Proc. of the 15th Symp. on Principles and Practice of Decla-rative Programming (PPDP). Madrid: ACM Press, 2013. 229-238.
[16]
Li Y, Zhu G, Feng Y. The L-depth eventual linear ranking functions for single-path linear constraint loops. In: Proc. of the 10th Int'l Symp. on Theoretical Aspects of Software Engineering (TASE 2016). IEEE, 2016. 30-37.
[17]
Leike J, Heizmann M. Ranking templates for linear loops. In: Proc. of the 20th Int'l Conf. on Tools and Algorithms for Construction and Analysis of Systems (TACAS 2014). Berlin, Heidelberg: Springer-Verlag, 2014, 8413: 172-186.
[18]
Ben-Amram AM, Genaim S. On multiphase-linear ranking functions. In: Computer Aided Verification. Cham: Springer-Verlag, 2017, 10427: 601-620.
[19]
Chen YH, Xia BC, Yang Lu, Zhan NJ, Zhou CC. Discovering non-linear ranking functions by solving semi-algebraic systems. In: Proc. of the Int'l Colloquium on Theoretical Aspects of Computing (ICTAC 2007). Berlin: Springer-Verlag, 2007. 34-49.
[20]
Cousot P. Proving program invariance and termination by parametric abstraction, langrangian relaxation and semidefinite programming. In: Proc. of the 6th Int'l Conf. on Verification, Model Checking, and Abstract Interpretation (VMCAI 2005). Springer-Verlag, 2005. 1-24.
[21]
Shen LY, Wu M, Yang ZF, Zeng ZB. Generating exact nonlinear ranking functions by symbolic-numeric hybrid method. Journal of Systems Science and Complexity, 2013, 26(2): 291-301. [doi:10.1007/s11424-013-1004-1]
[22]
Li Y. Termination of semi-algebraic loop programs. In: Proc. of the Int'l Symp. on Dependable Software Engineering: Theories, Tools and Applications (SETTA 2017). Springer-Verlag, 2017. 131-146.
[23]
Platt J. Sequential minimal optimization: A fast algorithm for training support vector machines. Technical Report, MSR-TR-98-14, Microsoft Research, 1998.
[24]
Hsu CW, Chang CC, Lin CJ. A practical guide to support vector classification. Technical Report, Department of Computer Science, Taiwan University, 2003.
[25]
Zhou ZH. Machine Lerning. Beijing: Tsinghua University Press, 2016.
[26]
The RegularChains Library. http://regularchains.org/index.html
[27]
Redlog-Computing with Logic. http://redlog.eu/
[25]
周志华. 机器学习. 北京: 清华大学出版社, 2016.