一类新的无参数填充函数及其在最小二乘法中的应用

A New Type of Non-Parameter Filled Function and Its Application in Least Squares

摘要: 针对一般约束优化问题，本文提出了一类新的无参数填充函数，分析该函数的理论性质并设计算法，利用Python编程语言进行数值实验并且与前人结果进行比较，表明该填充函数算法不仅是有效的而且在处理某些函数时效果更好。进一步地，尝试将该算法与化学实验酶促反应的反应速度和反应浓度关系的处理结合，实验结果表明该算法在实际应用中也有较好的适应性。

Abstract: In this paper, a novel non-parameter filled function for solving general constrained global optimi-zation problem is proposed. Then the theoretical properties of the function are argued and corre-sponding algorithm is given in this paper. Numerical experiments using the Python programming language and comparison with previous results show that the proposed filled function algorithm is not only effective but also works better when dealing with certain function. Furthermore, the algorithm is tentatively combined with the treatment of the reaction rate and reactant concentra-tion of the enzymatic reaction of chemical experiments, indicating that the algorithm also has good adaptability in practical cases.

1. 引言

非线性规划，作为数学规划的重要组成部分在经营管理、工程设计、科学研究、军事指挥等方面得到广泛的应用。近年来有关全局优化的理论和算法层出不穷。全局优化方法大致可分为三类：1) 在局部极小点中搜寻全局极小点，一般来说就是调用辅助函数跳出当前局部极小点得到更优的局部极小点；2) 采用启发式或者随机搜索方法如模拟退火算法，遗传算法、人工神经网络法等；3) 针对特殊结构的问题采用特定的解决方法如D.C规划。而本文所研究的内容是属于第一类，构造填充函数来求取全局极小点。填充函数法首先由西安交通大学葛仁溥教授 [1] 在1984年提出，是一类较为有效的确定性全局优化方法。随后其在文献 [2] 中对最初提出的双参数填充函数进行改进，提出了一类性能更优的单参数填充函数。

上海大学的张连生教授 [3] 对填充函数的定义进行改进，并提出一系列性能更优的填充函数；而且在张教授科研团队的深入研究下，填充函数内容越来越丰富。如：针对非线性等式约束问题，杨永健教授在文献 [4] 中提出了一类单参数填充函数；针对一般约束性优化问题，王伟祥教授在文献 [5] 中提出了一类单参数填充函数；张莹在文献 [6] 中针对非光滑约束问题提出了一类单参数填充函数。

在分析现有填充函数时发现多数研究者所构造的填充函数或者带有指数项 [7] ，或者结构比较复杂 [8] ，本文提出一种新的不带指数项的无参数填充函数，具有结构简单的优点。通过分析该填充函数的性质知：若填充函数 $p\left(x,x^{*}\right)>0$ ，则说明在 $x$ 的邻域内一定存在比当前局部极小点的值更小的局部极小点。这个优良性质应用到算法设计中，可以提升算法的实用性。通过几个数值算例，表明该算法是有效的。最后通过对化学实验酶促反应中反应速度和反应浓度的关系数据的处理得知：该算法在实际应用中也有较好的适应性。

2. 新的无参数填充函数

考虑如下问题：

$\left(P\right)\{\begin{array}{l}\min f\left(x\right),\\ \text{s}\text{.t}\text{.}x\in R^n.\end{array}$ (2.1)

其中 $f:R^n\to R$ 是连续可微函数。

假设2.1：函数 $f\left(x\right)$ 是强制性的，也就是当 $\Vert x\Vert \to +\infty$ 时，有 $f\left(x\right)\to +\infty$ 。

由假设2.1可以推出，一定存在一个紧集 $X\subset R^n$ ，它的内部包含了 $f\left(x\right)$ 的所有全局极小点和函数值较小的局部极小点，而且还可以认为目标函数在 $X$ 的边界上的函数值大于在其内部的函数值，因此，问题(2.1)可等价于如下问题：

$\{\begin{array}{l}\min f\left(x\right),\\ \text{s}\text{.t}\text{.}x\in X\end{array}$ (2.2)

假设2.2：问题(2.2)的不同局部极小点的个数是有限个。

下面给出葛仁浦教授在文献 [1] 中提出的填充函数的定义。

定义2.1：函数 $p\left(x,x^{*}\right)$ 称为 $f\left(x\right)$ 在局部极小点 $x^{*}$ 处的填充函数，如果满足：

1) $x^{*}$ 是 $p\left(x,x^{*}\right)$ 的一个严格局部极大点， $f\left(x\right)$ 在点 $x^{*}$ 处的盆谷 $B^{*}$ 成为 $p\left(x,x^{*}\right)$ 的峰的一部分。

2) $p\left(x,x^{*}\right)$ 在比 $B^{*}$ 高的盆谷里没有平稳点，即 $\nabla p\left(x,x^{*}\right)\ne 0$ 。

3) 如果存在比 $B^{*}$ 低的盆谷 $B_1^{*}$ ，则在 $x^{\prime }$ 和 $x^{*}$ 的连线上极小化 $p\left(x,x^{*}\right)$ 得到极小点 $x^{″}\in B^{*}$ ，其中， $N\left(x_1^{*},{\delta }_1\right)({\delta }_1>0)，x^{\prime }\in N\left(x_1^{*},{\delta }_1\right)$ 。

定义2.2：函数 $f\left(x\right)$ 在一局部极小点 $x^{*}$ 处的盆谷是指一连通域 $B^{*}$ ，具有下列性质：

1) $x^{*}\in B^{*}$ ；

2) 对于任意一点 $x\in B^{*}$ 使得 $x\ne x^{*}$ 及 $f\left(x\right)>f\left(x^{*}\right)$ ，存在一条从 $x$ 到c的下降路径。

针对问题(2.1)，提出了一类无参数填充函数如下：

$\begin{array}{r}\hfill p\left(x,x^{*}\right)=-\ln \left(1+{\Vert x-x^{*}\Vert }^2\right)\left(f\left(x\right)-f\left(x^{*}\right)\right).\end{array}$ (2.3)

其中 $x^{*}$ 是目标函数 $f\left(x\right)$ 的当前局部极小点。

下面来证明该函数是符合定义(2.1)的填充函数。

定理2.1： $x^{*}$ 是 $p\left(x,x^{*}\right)$ 的严格局部极大点。

证明：因为 $x^{*}$ 是 $f\left(x\right)$ 的一个局部极小点，则存在它的一个邻域 $N\left(x^{*},{\delta }^{*}\right)\left({\delta }^{*}>0\right)$ ，使得对于任意 $x\in N\left(x^{*},\delta \right)\backslash \left\{x^{*}\right\}$ ，有 $f\left(x\right)>f\left(x^{*}\right)$ 。则有，

$\begin{array}{c}p\left(x,x^{*}\right)=-\ln \left(1+\Vert x-x^{*}\Vert \right)\left(f\left(x\right)-f\left(x^{*}\right)\right)\\ <0=p\left(x^{*},x^{*}\right)\end{array}$

证毕。

定理2.2： $p\left(x,x^{*}\right)$ 在比 $B^{*}$ 高的盆谷 $B_1^{*}$ 里没有平稳点，即 $\nabla p\left(x,x^{*}\right)\ne 0$ 。

证明：由定理条件可知：对 $\forall x\in B_1^{*}$ ，有 $f\left(x\right)>f\left(x^{*}\right)$ ，又因为 $f\left(x\right)$ 在 $R^n$ 上连续可微，可以得到 $f\left(x\right)$ 在 $x^{*}$ 处的泰勒展开式为：

$f\left(x\right)=f\left(x^{*}\right)+{(x-x^{*})}^T\nabla f\left(x\right)+o\left(\Vert x-x^{*}\Vert \right)$

所以有：

$f\left(x\right)-f\left(x^{*}\right)={(x-x^{*})}^T\nabla f\left(x\right)+o\left(\Vert x-x^{*}\Vert \right)>0$

则有： ${\left(x-x^{*}\right)}^T\nabla f\left(x\right)>0$ 。 $p\left(x,x^{*}\right)$ 梯度如下：

$\nabla p\left(x,x^{*}\right)=-\left[\frac{2\left(x-x^{*}\right)}{1+{\Vert x-x^{*}\Vert }^2}\left(f\left(x\right)-f\left(x^{*}\right)\right)+\ln \left(1+{\Vert x-x^{*}\Vert }^2\right)\nabla f\left(x\right)\right]$

于是

$\begin{array}{c}\nabla p{\left(x,x^{*}\right)}^T\left(x-x^{*}\right)=\left[\frac{2{\Vert x-x^{*}\Vert }^2}{1+{\Vert x-x^{*}\Vert }^2}\left(f\left(x\right)-f\left(x^{*}\right)\right)+\ln \left(1+{\Vert x-x^{*}\Vert }^2\right)\nabla f{\left(x\right)}^T\left(x-x^{*}\right)\right]\\ <0\end{array}$

所以 $\nabla P\left(x,x^{*}\right)\ne 0$ 。证毕。

定理2.3：如果 $x^{*}$ 不是全局极小点，并且 $x_1^{*}$ 是离 $x^{*}$ 最近的一个局部极小点并且满足 $f\left(x_1^{*}\right)<f\left(x^{*}\right)$ ，则对 $\forall x^{\prime }\in N\left(x_1^{*},{\delta }_1^{*}\right)\left({\delta }_1^{*}>0\right)$ 则在 $x^{\prime }$ 和 $x^{*}$ 的连线上极小化 $p\left(x,x^{*}\right)$ ，得到 $p\left(x,x^{*}\right)$ 局部极小点 $x^{″}\in B^{*}$ 。

证明：因为 $x^{*}$ 是 $f\left(x\right)$ 的局部极小点，则存在 $x^{*}$ 的小邻域 $N\left(x^{*},{\delta }^{*}\right)\left({\delta }^{*}>0\right)$ ，对任意的 $x\in N\left(x^{*},{\delta }^{*}\right)$ ，有 $f\left(x\right)\ge f\left(x^{*}\right)$ ，则有

$p\left(x,x^{*}\right)=-\ln \left(1+{\Vert x-x^{*}\Vert }^2\right)\left(f\left(x\right)-f\left(x^{*}\right)\right)<0$

其中 $x\in N\left(x^{*},{\delta }^{*}\right)$ 且 $x\ne x^{*}$ 。

相似地，对于 $x_1^{*}$ ，则存在它的一个小邻域 $N\left(x_1^{*},{\delta }_1^{*}\right)$ ( ${\delta }_1^{*}>0$ )，对任意的 $x^{\prime }\in N\left(x_1^{*},{\delta }_1^{*}\right)$ ，有 $f\left(x^{*}\right)>f\left(x^{\prime }\right)\ge f\left(x_1^{*}\right)$ ，显然有：

$p\left(x^{\prime },x^{*}\right)=-\ln \left(1+{\Vert x^{\prime }-x^{*}\Vert }^2\right)\left(f\left(x^{\prime }\right)-f\left(x^{*}\right)\right)>0$

通过分析 $x^{\prime }$ 和 $x^{*}$ 的轨迹可知：在最初 $p\left(x^{*},x^{*}\right)=0$ ，x从 $x^{*}$ 处开始迭代搜索，有 $p\left(x,x^{*}\right)<0$ ，而当 $x$ 接近 $x^{\prime }$ 时，可知 $p\left(x,x^{*}\right)>0$ ，所以可以知道在 $f\left(x\right)\ge f\left(x^{*}\right)$ 时， $p\left(x,x^{*}\right)$ 递减，在 $f\left(x\right)<f\left(x^{*}\right)$ 时， $p\left(x,x^{*}\right)$ 递增。由此可以知道在 $x^{\prime }$ 和 $x^{*}$ 的连线上存在 $x^{″}$ 极小化 $p\left(x,x^{*}\right)$ 。证毕。

定理2.1~定理2.3表明 $p\left(x,x^{*}\right)$ 是符合定义2.1的填充函数。下面的定理2.4表明该填充函数具有一个很好的性质，能在算法执行中快速地找到更优的局部极小点。

定理2.4：若 $x^{*}$ 是 $f\left(x\right)$ 的全局极小点，则对任意的 $x\in X$ 都有 $p\left(x,x^{*}\right)\le 0$ 成立；否则存在点 $x^{**}$ 以及它的邻域 $N\left(x^{**},\delta \right)\left(\delta >0\right)$ ，使得对任意的 $x\in N\left(x^{**},\delta \right)$ ，成立 $p\left(x,x^{*}\right)>0$ 。

证明：若 $x^{*}$ 是 $f\left(x\right)$ 的全局极小点，则对任意 $x\in X$ ，都有 $f\left(x\right)\ge f\left(x^{*}\right)$ ，所以 $p\left(x,x^{*}\right)=-\ln \left(1+{\Vert x-x^{*}\Vert }^2\right)\left(f\left(x\right)-f\left(x^{*}\right)\right)\le 0$ ；若 $x^{*}$ 不是 $f\left(x\right)$ 的全局极小点，则存在点 $x^{**}$ 使得 $f\left(x^{**}\right)<f\left(x^{*}\right)$ ，由 $f\left(x\right)$ 的连续性可知存在 $x^{**}$ 的邻域 $N\left(x^{**},{\delta }^{**}\right)\left({\delta }^{**}>0\right)$ ，则 $\forall x\in N\left(x^{**},{\delta }^{**}\right)$ 有 $f\left(x\right)<f\left(x^{*}\right)$ 。所以有 $p\left(x,x^{*}\right)>0$ 。证毕。

注：由定理2.2的证明可以发现 $d=x-x^{*}$ 是填充函数 $p\left(x,x^{*}\right)$ 的一个下降方向；通过定理2.4可以发现，当遇到 $p\left(x,x^{*}\right)>0$ 时，可以知道 $f\left(x\right)<f\left(x^{*}\right)$ ，则在算法实施过程中说明找到一个比当前局部极小点的值更小的点。

3. 填充函数算法

填充函数法是求解全局优化问题的有效方法，该方法的基本思想是：在可行域 $X$ 中选取一点 $x_0$ ，利用已有的局部优化算法如最速下降法、共轭梯度法和拟牛顿法等来求出目标函数 $f\left(x\right)$ 的一个局部极小点 $x^{*}$ ，然后在 $x^{*}$ 处构造一个填充函数，在 $x^{*}$ 的某个邻域内选取一点作为初始点，极小化填充函数，得到一点 $x_1^{*}$ ，使得 $f\left(x_1^{*}\right)<f\left(x^{*}\right)$ ；然后用 $x_1^{*}$ 代替 $x_0$ 并重复上述过程，直到找到全局极小点。根据这种思想，结合本文所提出填充函数的特点以及文献的算法 [9] ，提出相应改进算法。步骤如下：

步骤0：初始化数据，设定搜索步长 $\delta >0,k=1$ ，选择方向 $d_i\left(i=1,2,\cdots ,m\right),m\ge 2n$ ，其中 $n$ 为变量的维度，选取一个初始点 $x_0\in X$ 。

步骤1：以 $x_0$ 为初始点，应用局部优化算法如拟牛顿法求出 $f\left(x\right)$ 的一个局部极小点 $x^{*}$ ，并求出局部极小值 $f\left(x^{*}\right)$ 。

步骤2：设置 $i=1$ 。

步骤3：如果 $i\le m$ ，转步骤4；否则 $x^{*}$ 作为全局极小点，停止计算。

步骤4：令 $x_1=x^{*}+\delta d_i$ ，构造 $f\left(x\right)$ 在 $x^{*}$ 处的填充函数 $P\left(x,x^{*}\right)$ 。然后从 $x_1$ 开始使用局部优化算法得到 $P\left(x,x^{*}\right)$ 的局部极小点 $x_1^{*}$ 。

步骤5：令 $x=x_1^{*}$ 。

步骤6：若 $p\left(x,x^{*}\right)>0$ ，则令 $k=k+1$ 且 $x_0=x$ ，转步骤1，否则令 $i=i+1$ ，转步骤2，若 $x\notin X$ ，则令 $i=i+1$ ，转到步骤3。

在步骤4中的 $x_1$ 一般是这样选取的： $m=2n$ ，令 $x_1$ 的选取是对称的。例如对 $n=2$ ， $x_1$ 为：

$x^{*}+\delta \left(1,0\right),x^{*}+\delta \left(0,1\right),x^{*}-\delta \left(1,0\right),x^{*}-\delta \left(0,1\right)$

4. 数值实验

为了验证本文所提出的无参数填充函数算法是有效的，应用上述算法，使用常用的几个测试函数，使用Python3.5.6进行编程计算。

4.1. 测试函数

1) 6-Hump Back Camel函数 [10]

$\{\begin{array}{l}\min f\left(x\right)=4x_1^2-2.1x_1^4+\frac{1}{3}{x}_1^6-x_1{x}_2-4x_2^2+4x_2^4\\ \text{s}\text{.t}\text{.}-3\le x_1,x_2\le 3\end{array}$

2) Restrign函数 [10]

$\{\begin{array}{l}\min f\left(x\right)=x_1^2+x_2^2-\cos \left(18x_1\right)-\cos \left(18x_2\right)\\ \text{s}\text{.t}\text{.}-1\le x_1,x_2\le 1\end{array}$

3) 文献 [11] 中的一个二维函数(c = 0.05)

$\{\begin{array}{l}\min f\left(x\right)={\left[1-2x_2+c*\sin \left(4\text{π}{x}_2\right)-x_1\right]}^2+{\left[x_2-0.5\sin \left(2\text{π}{x}_1\right)\right]}^2\\ \text{s}\text{.t}\text{.}-10\le x_1,x_2\le 10\end{array}$

4.2. 数值结果

以表格形式给出数值结果，表1~3分别对应测试函数(1)~(3)的数值结果，表4是文献 [11] 算例与本文算例的数值对照表。下面是表格中符号的解释。

$k$ ：求解第 $k$ 次局部极小点的迭代步数；

$x_k$ ：满足 $x_k\in X$ 的第 $k$ 次初始点；

$x_k^{*}$ ：第 $k$ 个局部极小点；

$f\left(x_k^{*}\right)$ ：第 $k$ 个局部极小点的目标函数值；

$\text{iter}$ ：算法迭代步数；

$x^{*}$ ：算法最终得到的全局极小点；

$f\left(x^{*}\right)$ ：算法最终得到的全局极小值。

Table 1. Original value x = ( 1 , 2 ) T

表1. 初始点为 $x={\left(1,2\right)}^T$

Table 2. Original value x = ( 1 , 1 ) T

表2. 初始点为 $x={\left(1,1\right)}^T$

Table 3. c = 0.05 ，Original value x = ( 6 , − 2 ) T

表3. $c=0.05$ ，初始点为 $x={\left(6,-2\right)}^T$

Table 4. Compared with paper [11]

表4. 本文算例与文献 [11] 算例对照表

本文算例与文献 [11] 算例对比表明，应用本文算法所得结果精度更高，迭代步数更少，特别对算例3而言，更好的迭代步数得到更优的全局极小点。

5. 填充函数在实验数据处理中的应用

5.1. 背景和问题

酶是一种具有特异性的高效生物催化剂，绝大多数的酶是活细胞产生的蛋白质。酶的催化条件温和，在常温、常压下即可进行。酶催化的反应称为酶促反应，要比相应的非催化反应要快 ${10}^3\sim {10}^{17}$ 倍。酶促反应动力学简称酶促动力学，主要研究酶促反应的速度和底物(即反应物)浓度以及其他因素的关系。在底物浓度很低时酶促反应是一级反应；当底物浓度处于中间范围时，是混合级反应；当底物浓度增加时，向零级反应过度 [12] 。

某生化系学生为了研究嘌呤霉素在某项酶促反应中对反应速度与底物浓度之间关系的影响，设计了两个实验，一个实验中所使用的酶是经过嘌呤霉素处理的，而另一个实验所用的酶是未经嘌呤霍素处理过的，所得的实验数据见表5。试根据问题的背景和这些数据建立一个合适的数学模型，来反映这项酶促反应的速度与底物浓度以及嘌呤霉素处理与否之间的关系。

Table 5. Reaction rate and substrate concentration data in puromycin experiments

表5. 嘌呤霉素实验中的反应速度与底物浓度数据

注：ppm = 0.001%。

本文仅对嘌呤毒素处理过的数据进行拟合，为了方便处理，对数据取平均值，数据如表6：

Table 6. Average the data after treatment with scorpion toxin

表6. 对嘌呤毒素处理后数据取平均值处理

5.2. 分析与假设

记酶促反应速度为 $y$ ，底物浓度为 $x$ ，二者之间的关系写作 $y=f\left(x,\beta \right)$ ，其中 $\beta$ 为参数。由酶促反应的基本性质可知，当底物浓度较小时，反应速度大致与浓度成正比(即一级反应)；当底物浓度很大，渐进饱和时，反应速度将趋于一个固定值–最终反应速度(即零级反应)。下面是针对该种性质提出的模型：

指数增长模型

$y=f\left(x,\beta \right)={\beta }_1\left(1-e^{-{\beta }_{2}x}\right)$

本文将求解系数 $\beta$ 转换为下面无约束优化问题：

$\begin{array}{c}\min f\left(x\right)={{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(y_i-{\hat{y}}_i\right)}}^2\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left({\beta }_1\left(1-e^{-{\beta }_2{x}_{\text{}}}\right)-{\hat{y}}_i\right)}^2}\end{array}$

给定初始值 ${\beta }_0={\left(155,14\right)}^T$ ，经过2次迭代得到 ${\beta }^{*}={\left(181.54213939,13.10716834\right)}^T$ ，而文献 [12] 中所给出参数估计 ${\beta }_1={\left(192.0969,11.3846\right)}^T$ 。本文目标函数值为 $918.2109$ ，剩余标准差 $\text{RMSE}=17.4949$ ，而

Figure 1. Fitting curve

图1. 拟合曲线

文献 [12] 所得的目标函数值为 $1041.9798$ ，剩余标准差 $\text{RMSE}=18.6367$ 。由此可知，使用本文算法所得的全局极小值更小，拟合效果更好。所得拟合曲线如图1所示。

6. 结语

本文给出了一个形式简单，性能更优的填充函数，分析其性质并根据该性质设计算法，数值实验表明该算法是有效的，与文献 [11] 数值结果的比较可知，该算法在处理某些函数时，有着较好的适应性。同时，本文也尝试使用该算法求解最小二乘问题，实验结果表明该算法在实际应用中也有较好的适应性。

基金项目

浙江省科技厅资助项目(LGN19C040001)；浙江省教育厅科研项目(Y201636028)；国家社科基金资助项目(17BTJ028)。