开放科学(资源服务)标志码(OSID)：

0 概述

支持向量机（Support Vector Machine，SVM）是20世纪90年代由VAPNIK^[1-2]提出的一种机器学习算法。与传统的以降低经验风险为目标的神经网络相比，SVM的主要思想是最小化结构风险，这使得SVM具有良好的泛化性能^[3-4]。目前，SVM在入侵检测^[5]、图像处理^[6]、故障诊断^[7]、干扰识别^[8]等领域得到广泛应用。然而，SVM在训练大规模数据集时，存在时间复杂度高导致训练速度慢的问题^[9]。KHEMCHANDANI等^[10]提出孪生支持向量机（Twin Support Vector Machine，TSVM）。与SVM相比，TSVM仅需求解两个较小规模的二次规划问题，因此训练时间仅为SVM的1/4左右。CHEN等^[11]提出投影孪生支持向量机（Projection Twin Support Vector Machine，PTSVM）。PTSVM通过为每个类寻找一个最优投影轴使投影点类内方差最小化，构建分类模型。

PENG^[12]把TSVM的思想用于回归领域，提出孪生支持向量回归（Twin Support Vector Regression，TSVR）算法。TSVR通过求解两个较小规模的二次规划问题寻求回归函数，与传统SVR相比，TSVR具有更好的泛化性能和更快的训练速度。CHEN等^[13]引入Sigmoid光滑函数，对TSVR的目标函数进行光滑处理，提出光滑孪生支持向量回归（Smooth Twin Support Vector Regression，STSVR）算法。STSVR通过求解一对无约束优化问题，获得了比TSVR更快的训练速度。PENG等^[14]基于PTSVM的思想，提出投影孪生支持向量回归（Projection Twin Support Vector Regression，PTSVR）算法，包括双边移位投影孪生支持向量回归（Pair-shifted PTSVR，PPTSVR）算法和单边移位投影孪生支持向量回归（Single-shifted PTSVR，SPTSVR）算法。PTSVR将训练集中的每个训练点进行上下移位得到两个新的移位集，从而利用PTSVM的思想求解两个最优超平面的法向量。PPTSVR和SPTSVR的主要区别在于，PPTSVR在两个移位集上构建回归函数，而SPTSVR则是在原始集和一个移位集上构建回归函数。实验结果表明，PTSVR有着比TSVR更好的预测性能。

然而，PPTSVR在寻找最优超平面的过程中，将所有训练样本对超平面的作用视为相同，没有反映数据在空间中的分布情况，当训练样本中存在异常点时，会削弱算法的拟合性能。本文采用孤立森林法赋予每个训练样本不同的权值，通过赋予潜在的异常点很小的权值，削弱异常点对超平面构造的影响，从而提高算法的拟合性能。同时，为了避免在对偶空间中求解二次规划问题，引入正号函数，将有约束优化问题转化为不光滑的无约束优化问题，并采用Sigmoid光滑函数进行光滑处理，提出一种加权光滑投影孪生支持向量回归（Weighted Smooth Projection Twin Support Vector Regression，WSPTSVR）算法，以证明其任意阶光滑且全局收敛的特性，并采用牛顿迭代法进行求解。

1 WSPTSVR算法

本节首先采用孤立森林法判断样本中的潜在异常点，并通过异常分数机制赋予每个样本不同的权值；其次引入Sigmoid光滑函数，提出WSPTSVR算法，并证明其具有全局唯一最优解；最后在原空间中使用牛顿迭代法进行求解。

1.1 基于孤立森林法的加权系数确定

孤立森林法能够有效地检测出样本中的潜在异常点，即使样本中没有异常点，也能够根据异常分数赋予每个样本相应的权值，从而提高算法的拟合性能^[15-17]。孤立森林法首先递归地分割给定样本集，直到每个样本都被单独分离出来，然后根据每个样本分离的路径长度计算其异常分数，根据异常分数大小判断样本是否为异常点并赋予每个样本点相应的权值。通过孤立森林法确定加权系数共分为以下3个步骤：

1）通过训练样本的子集构建$ t $棵孤立树，建立孤立森林。

2）在训练后的孤立森林中，根据每个样本点分离所需的路径长度赋予其相应的异常分数$ {S}_{i} $，$ {S}_{i} $的取值范围为0~1，且取值越大，该点为异常点的可能性越大，异常分数$ {S}_{i} $计算如下：

$ {S}_{i}={2}^{-\frac{E\left(h\right(i\left)\right)}{c\left(n\right)}} $

(1)

其中：$ h\left(i\right) $为样本$ i $的路径长度，即从根节点到叶子节点的边数，而异常点更容易被孤立，因此$ h\left(i\right) $较小；$ E\left(h\right(i\left)\right) $为样本$ i $在一组孤立森林中路径长度的均值；$ c\left(n\right) $是样本数为$ n $时路径长度的均值，用来标准化样本$ i $的路径长度$ h\left(i\right) $。$ c\left(n\right) $计算如下：

$ c\left(n\right)=2H(n-1)-2(n-1)/n $

(2)

其中：$ H(n-1)=\mathrm{l}\mathrm{n}(n-1)+0.577\mathrm{ \;}215\mathrm{\;}664\mathrm{\;}9 $为调和数。

3）通过孤立森林法中的异常分数机制为每个样本点赋予相应的权值：

$ {\rho }_{i}=1-{S}_{i}, i=\mathrm{1, 2}, \cdots , n $

(3)

其中：$ {\rho }_{i} $表示第$ i $个样本的权值。

文献[15]的研究表明，当样本点的异常分数$ {S}_{i} > 0.60 $时，即视该样本点为潜在异常点，应赋予很小的权值，则样本点的权值矩阵可表示为对角矩阵：

$ {\boldsymbol{W}}_{ij}=\left\{\begin{array}{l}{\rho }_{i},  \mathrm{当}i=j\mathrm{且}{S}_{i}\le 0.6\\ {10}^{-5}, \mathrm{当}i=j\mathrm{且}{S}_{i} > 0.6\\  0, \mathrm{其}\mathrm{他}\end{array}\right. $

(4)

1.2 线性情况

假设训练集用$ D=\left\{({\boldsymbol{x}}_{i};{\boldsymbol{y}}_{i}), i\in {\rm I}=\{\mathrm{1, 2}, \cdots , n\}\right\} $表示，其中，$ {\boldsymbol{x}}_{i}\in {R}^{m} $为输入，$ {\boldsymbol{y}}_{i}\in R $为输出，$ I $为指标集，$ n $为样本个数，$ m $为样本特征数。为了后续表示简单起见，分别用矩阵$ \boldsymbol{X}=[{\boldsymbol{x}}_{1}, {\boldsymbol{x}}_{2}, \cdots , {\boldsymbol{x}}_{n}{]}^{\mathrm{T}}\in {R}^{n\times m} $和向量$ \boldsymbol{Y}=[{\boldsymbol{y}}_{1}, {\boldsymbol{y}}_{2}, \cdots , {\boldsymbol{y}}_{n}{]}^{\mathrm{T}}\in {R}^{n} $表示输入和输出，用矩阵$ \boldsymbol{J}=[\boldsymbol{X}, \boldsymbol{Y}]\in {R}^{n\times (m+1)} $表示训练样本，分别用矩阵$ {\boldsymbol{A}}_{\epsilon }=[\boldsymbol{X}, \boldsymbol{Y}+\epsilon \boldsymbol{e}] $和$ {\boldsymbol{B}}_{\epsilon }=[\boldsymbol{X}, \boldsymbol{Y}-\epsilon \boldsymbol{e}]\in {R}^{n\times (m+1)} $表示上移和下移训练样本，其中$ \epsilon> 0 $为不敏感因子，$ \boldsymbol{e} $为全1的列向量。

在双边移位投影孪生支持向量回归算法的目标函数中引入权值矩阵$ {\boldsymbol{W}}_{ij} $，构造如下最优化问题：

$ \underset{{\boldsymbol{\omega }}_{1}, {\delta }_{1}, \boldsymbol{\xi }}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{{C}_{3}}{2}\left(\right|\left|{\boldsymbol{\omega }}_{1}\right|{|}_{2}^{2}+{\delta }_{1}^{2})+\frac{1}{2}(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}{)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}\right) +{C}_{1}{\boldsymbol{e}}^{\mathrm{T}}\boldsymbol{\xi }, \mathrm{s}.\mathrm{t}.\boldsymbol{F}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}\le -\boldsymbol{e}+\boldsymbol{\xi },     \boldsymbol{\xi }\ge 0 $

(5)

$ \underset{{\boldsymbol{\omega }}_{2}, {\delta }_{2}, \boldsymbol{\eta }}{\mathrm{m}\mathrm{i}\mathrm{n}}  \frac{{C}_{4}}{2}\left(\right|\left|{\boldsymbol{\omega }}_{2}\right|{|}_{2}^{2}+{\delta }_{2}^{2})+\frac{1}{2}(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}{)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}\right)+{C}_{2}{\boldsymbol{e}}^{\mathrm{T}}\boldsymbol{\eta }, \mathrm{s}.\mathrm{t}.\boldsymbol{G}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}\ge \boldsymbol{e}-\boldsymbol{\eta }, \;\;\;\boldsymbol{\eta }\ge 0 $

(6)

其中：$ \left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right] $和$ \left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]\in {R}^{m+1} $为两个超平面的法向量；$ {C}_{1}, {C}_{2} > 0 $为惩罚参数；$ {C}_{3}, {C}_{4} > 0 $为正则化参数；$ \boldsymbol{\xi }, \boldsymbol{\eta }\ge 0 $为不敏感损失参数；$ \boldsymbol{E}=\boldsymbol{J}-\boldsymbol{e}{\boldsymbol{e}}^{\mathrm{T}}\boldsymbol{J}/n $；$ \boldsymbol{F}={\boldsymbol{B}}_{\epsilon }-\boldsymbol{e}{\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{A}}_{\epsilon }/n $；$ \boldsymbol{G}={\boldsymbol{A}}_{\epsilon }-\boldsymbol{e}{\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{B}}_{\epsilon }/n $；$ \left|\right|  \cdot   |{|}_{2} $表示L₂范数。

然而，式（5）和式（6）需要在对偶空间中求解二次规划问题，训练速度较慢。为了加快训练速度，采用正号函数将其转化为无约束优化问题，在原空间中直接进行求解。

由Karush-Kuhn-Tucker（KKT）条件可得：

$ \boldsymbol{\xi }=\mathrm{m}\mathrm{a}\mathrm{x}\left\{0, \boldsymbol{e}+\boldsymbol{F}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}\right\} $

(7)

$ \boldsymbol{\eta }=\mathrm{m}\mathrm{a}\mathrm{x}\left\{0, \boldsymbol{e}-\boldsymbol{G}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}\right\} $

(8)

式（5）和式（6）可改写为如下无约束优化问题：

$ \underset{{\boldsymbol{\omega }}_{1}, {\delta }_{1}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{{C}_{3}}{2}\left(\right|\left|{\boldsymbol{\omega }}_{1}\right|{|}_{2}^{2}+{\delta }_{1}^{2})+\frac{1}{2}(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}{)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}\right)+{C}_{1}{\boldsymbol{e}}^{\mathrm{T}}\mathrm{m}\mathrm{a}\mathrm{x}\left\{0, \boldsymbol{e}+\boldsymbol{F}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}\right\} $

(9)

$ \underset{{\boldsymbol{\omega }}_{2}, {\delta }_{2}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{{C}_{4}}{2}\left(\right|\left|{\boldsymbol{\omega }}_{2}\right|{|}_{2}^{2}+{\delta }_{2}^{2})+\frac{1}{2}(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}{)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}\right)+{C}_{2}{\boldsymbol{e}}^{\mathrm{T}}\mathrm{m}\mathrm{a}\mathrm{x}\left\{0, \boldsymbol{e}-\boldsymbol{G}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}\right\} $

(10)

定理1  ^[18] 无约束优化问题式（9）和式（10）连续但不光滑。

证明  分别令正号函数$ ({\boldsymbol{x}}_{1}{)}_{+}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{0, \boldsymbol{e}+\boldsymbol{F}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}\right\} $和$ ({\boldsymbol{x}}_{2}{)}_{+}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{0, \boldsymbol{e}-\boldsymbol{G}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}\right\} $。显然，式（9）和式（10）的可微性和光滑性取决于正号函数$ \left({x}_{i}{)}_{+}\right(i=\mathrm{1, 2}) $。由于$ ({\boldsymbol{x}}_{i}{)}_{+}=\left\{\begin{array}{l}{\boldsymbol{x}}_{i},     {\boldsymbol{x}}_{i}\ge 0\\ 0    ,     {\boldsymbol{x}}_{i} < 0\end{array}\right. $，可得$ ({\boldsymbol{x}}_{i}{)}_{+} $在$ {\boldsymbol{x}}_{i}=0 $处的左右导数分别为$ 0 $和$ 1 $，因此$ ({\boldsymbol{x}}_{i}{)}_{+} $在$ {\boldsymbol{x}}_{i}=0 $处不可微。又由于$ ({\boldsymbol{x}}_{i}{)}_{+} $还可以表示为$ ({\boldsymbol{x}}_{i}{)}_{+}=\frac{\left|\right({\boldsymbol{x}}_{i}{)}_{+}|+({\boldsymbol{x}}_{i}{)}_{+}}{2} $，因此$ ({\boldsymbol{x}}_{i}{)}_{+} $是连续函数。所以，正号函数$ ({\boldsymbol{x}}_{i}{)}_{+} $连续但不光滑。

定理1表明式（9）和式（10）不光滑，因此无法使用牛顿迭代法等梯度优化方法求解。为此，采用Sigmoid光滑函数逼近正号函数$ ({\boldsymbol{x}}_{1}{)}_{+} $和$ ({\boldsymbol{x}}_{2}{)}_{+} $。

Sigmoid光滑函数$ p(\boldsymbol{x}, \alpha ) $定义如下：

$ p(\boldsymbol{x}, \alpha )=\boldsymbol{x}+\frac{1}{\alpha }\mathrm{l}\mathrm{o}{\mathrm{g}}_{a}(1+{\mathrm{e}}^{-\alpha x}) $

(11)

其中：$ \alpha $为正常数。

据此可得正号函数$ ({\boldsymbol{x}}_{1}{)}_{+} $和$ ({\boldsymbol{x}}_{2}{)}_{+} $的Sigmoid光滑函数分别如下：

$ ({\boldsymbol{x}}_{1}{)}_{+}=p(\boldsymbol{e}+\boldsymbol{F}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}, \alpha ) $

(12)

$ ({\boldsymbol{x}}_{2}{)}_{+}=p(\boldsymbol{e}-\boldsymbol{G}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}, \alpha ) $

(13)

将式（12）和式（13）分别代入式（9）和式（10），得到线性情况下加权光滑投影孪生支持向量回归算法的两个无约束优化问题，如式（14）、式（15）所示：

$ \underset{{\boldsymbol{\omega }}_{1}, {\delta }_{1}}{\mathrm{m}\mathrm{i}\mathrm{n}}   {P}_{1}=\frac{{C}_{3}}{2}\left(\right|\left|{\boldsymbol{\omega }}_{1}\right|{|}^{2}+{\delta }_{1}^{2})+\frac{1}{2}(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}{)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}\right)+{C}_{1}{\boldsymbol{e}}^{\mathrm{T}}p(\boldsymbol{e}+\boldsymbol{F}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}}, \alpha ) $

(14)

$ \underset{{\boldsymbol{\omega }}_{2}, {\delta }_{2}}{\mathrm{m}\mathrm{i}\mathrm{n}}   {P}_{2}=\frac{{C}_{4}}{2}\left(\right|\left|{\boldsymbol{\omega }}_{2}\right|{|}^{2}+{\delta }_{2}^{2})+\frac{1}{2}(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}{)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(\boldsymbol{E}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}\right) +{C}_{2}{\boldsymbol{e}}^{\mathrm{T}}p(\boldsymbol{e}-\boldsymbol{G}{\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}}, \alpha ) $

(15)

引理1^[19]   若矩阵$ \boldsymbol{A}\in {R}^{n\times n} $是实对称矩阵，则$ \boldsymbol{A} $是正定矩阵当且仅当存在矩阵$ \boldsymbol{B}\in {R}^{m\times n} $使得$ \boldsymbol{A}={\boldsymbol{B}}^{\mathrm{T}}\boldsymbol{B} $。

定理2   式（14）中的$ {P}_{1} $和式（15）中的$ {P}_{2} $是严格凸函数。

证明   对任意的正常数$ \alpha $，$ {P}_{1} $显然是连续可微的，以下证明其是严格凸函数：

令$ {\boldsymbol{u}}_{1}={\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}} $，则由式（14）可得：

$ \nabla {P}_{1}\left({\boldsymbol{u}}_{1}\right)=({C}_{3}\boldsymbol{I}+{\boldsymbol{E}}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\boldsymbol{E}){\boldsymbol{u}}_{1}+\\                       {C}_{1}{\boldsymbol{F}}^{\mathrm{T}}\frac{1}{1+\mathrm{e}\mathrm{x}\mathrm{p}(-\alpha (\boldsymbol{F}{\boldsymbol{u}}_{1}+\boldsymbol{e}\left)\right)} $

(16)

$ \begin{array}{l}{\nabla }^{2}{P}_{1}\left({\boldsymbol{u}}_{1}\right)={C}_{3}\boldsymbol{I}+{\boldsymbol{E}}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\boldsymbol{E}+\\                      \alpha {C}_{1}{\boldsymbol{F}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\frac{\mathrm{e}\mathrm{x}\mathrm{p}(-\alpha (\boldsymbol{F}{\boldsymbol{u}}_{1}+\boldsymbol{e}\left)\right)}{{\left(1+\mathrm{e}\mathrm{x}\mathrm{p}(-\alpha (\boldsymbol{F}{\boldsymbol{u}}_{1}+\boldsymbol{e}\left)\right)\right)}^{2}}\right)\boldsymbol{F}\end{array} $

(17)

考虑到权值矩阵和对角矩阵都是正定阵，两者可以分别分解为$ {\boldsymbol{P}}^{\mathrm{T}}\boldsymbol{P} $和$ {\boldsymbol{Q}}^{\mathrm{T}}\boldsymbol{Q} $的形式，其中$ \boldsymbol{P} $和$ \boldsymbol{Q} $分别为两者的平方根矩阵。因此，$ {\nabla }^{2}{P}_{1}\left({\boldsymbol{u}}_{1}\right) $可以分解如下：

$ \begin{array}{l}{\nabla }^{2}{P}_{1}\left({\boldsymbol{u}}_{1}\right)={C}_{3}\boldsymbol{I}+{\boldsymbol{E}}^{\mathrm{T}}{\boldsymbol{P}}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{E}+\alpha {C}_{1}{\boldsymbol{F}}^{\mathrm{T}}{\boldsymbol{Q}}^{\mathrm{T}}\boldsymbol{Q}\boldsymbol{F}=\\                         {C}_{3}\boldsymbol{I}+{\left(\boldsymbol{P}\boldsymbol{E}\right)}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{E}+\alpha {C}_{1}{\left(\boldsymbol{Q}\boldsymbol{F}\right)}^{\mathrm{T}}\boldsymbol{Q}\boldsymbol{F}\end{array} $

(18)

其中：$ {C}_{1}, {C}_{3} $和$ \alpha $都是正标量；$ \boldsymbol{I} $为正定矩阵。由引理1可得$ {\nabla }^{2}{P}_{1}\left({\boldsymbol{u}}_{1}\right) $是正定的，因此$ {P}_{1} $是严格凸函数。同理可证$ {P}_{2} $也是连续可微的严格凸函数。

式（14）和式（15）二阶可微，因此可以使用具有快速收敛能力的牛顿迭代法进行求解^[20-21]。由定理2可知式（14）和式（15）是严格凸的，因此使用牛顿迭代法求解可以全局收敛，并得到唯一最优解。牛顿迭代法可表示如下：

$ {\boldsymbol{u}}_{1}^{{k}+1}={\boldsymbol{u}}_{1}^{{k}}-\left({\nabla }^{2}{P}_{1}\right({\boldsymbol{u}}_{1}{\left)\right)}^{-1}(\nabla {P}_{1}({\boldsymbol{u}}_{1}\left)\right) $

(19)

$ {\boldsymbol{u}}_{2}^{{k}+1}={\boldsymbol{u}}_{2}^{{k}}-\left({\nabla }^{2}{P}_{2}\right({\boldsymbol{u}}_{2}{\left)\right)}^{-1}(\nabla {P}_{2}({\boldsymbol{u}}_{2}\left)\right) $

(20)

其中：$ {\boldsymbol{u}}_{2}={\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}} $。$ \nabla {P}_{2}\left({\boldsymbol{u}}_{2}\right) $和$ {\nabla }^{2}{P}_{2}\left({\boldsymbol{u}}_{2}\right) $分别如下：

$ \nabla {P}_{2}\left({\boldsymbol{u}}_{2}\right)=({C}_{4}\boldsymbol{I}+{\boldsymbol{E}}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\boldsymbol{E}){\boldsymbol{u}}_{2}-\\ {C}_{2}{\boldsymbol{G}}^{\mathrm{T}}\frac{1}{1+\mathrm{e}\mathrm{x}\mathrm{p}(-\alpha (-\boldsymbol{G}{\boldsymbol{u}}_{2}+\boldsymbol{e}\left)\right)} $

(21)

$ \begin{array}{l}{\nabla }^{2}{P}_{2}\left({\boldsymbol{u}}_{2}\right)={C}_{4}\boldsymbol{I}+{\boldsymbol{E}}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\boldsymbol{E}+\\ \alpha {C}_{2}{\boldsymbol{G}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\frac{\mathrm{e}\mathrm{x}\mathrm{p}(-\alpha (-\boldsymbol{G}{\boldsymbol{u}}_{2}+\boldsymbol{e}\left)\right)}{(1+\mathrm{e}\mathrm{x}\mathrm{p}(-\alpha (-\boldsymbol{G}{\boldsymbol{u}}_{2}{+\boldsymbol{e}\left)\right))}^{2}}\right)\boldsymbol{G}\end{array} $

(22)

首先，通过上述牛顿迭代法求解式（14）和式（15）可得两个最优超平面的法向量$ {\boldsymbol{u}}_{1}={\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}} $和$ {\boldsymbol{u}}_{2}={\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}} $。

然后，通过求解式（23）和式（24）两个无约束最小化问题，可得两个最优超平面的偏置分别如式（25）、式（26）所示：

$ \mathrm{m}\mathrm{i}\mathrm{n}\left\{\frac{1}{2}\sum\limits _{i=1}^{n}({\boldsymbol{\omega }}_{1}^{\mathrm{T}}{\boldsymbol{x}}_{i}+{\delta }_{1}({\boldsymbol{y}}_{i}+\epsilon )+{b}_{1}{)}^{2}\right\} $

(23)

$ \mathrm{m}\mathrm{i}\mathrm{n}\left\{\frac{1}{2}\sum\limits _{i=1}^{n}({\boldsymbol{\omega }}_{2}^{\mathrm{T}}{\boldsymbol{x}}_{i}+{\delta }_{2}({\boldsymbol{y}}_{i}-\epsilon )+{b}_{2}{)}^{2}\right\} $

(24)

$ {b}_{1}=-{\boldsymbol{\omega }}_{1}^{\mathrm{T}}\stackrel{-}{\boldsymbol{x}}-{\delta }_{1}(\stackrel{-}{\boldsymbol{y}}+\epsilon ) $

(25)

$ {b}_{2}=-{\boldsymbol{\omega }}_{2}^{\mathrm{T}}\stackrel{-}{\boldsymbol{x}}-{\delta }_{2}(\stackrel{-}{\boldsymbol{y}}-\epsilon ) $

(26)

在通常情况下，$ {\delta }_{1}, {\delta }_{2}\ne 0 $^{[14, 21]}，可得移位函数$ {f}_{1}\left(\boldsymbol{x}\right)=-({\boldsymbol{\omega }}_{1}^{\mathrm{T}}\boldsymbol{x}+{b}_{1})/{\delta }_{1} $和$ {f}_{2}\left(\boldsymbol{x}\right)=-({\boldsymbol{\omega }}_{2}^{\mathrm{T}}\boldsymbol{x}+{b}_{2})/{\delta }_{2} $。

最终，可得回归函数如下：

$ \begin{array}{l}f\left(\boldsymbol{x}\right)=\frac{1}{2}\left({f}_{1}\left(\boldsymbol{x}\right)+{f}_{2}\left(\boldsymbol{x}\right)\right)=\\ \;\;\;\;\;\;\;\;\;\;\;\;\;-\frac{1}{2}{\left(\frac{{\boldsymbol{\omega }}_{1}}{{\delta }_{1}}+\frac{{\boldsymbol{\omega }}_{2}}{{\delta }_{2}}\right)}^{\mathrm{T}}\boldsymbol{x}-\frac{1}{2}\left(\frac{{b}_{1}}{{\delta }_{1}}+\frac{{b}_{2}}{{\delta }_{2}}\right)\end{array} $

(27)

1.3 非线性情况

在非线性情况下，利用核技巧，构造加权光滑投影孪生支持向量回归算法的两个无约束优化问题如下：

$ \begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}{P}_{1}=\frac{1}{2}\left(φ \right(\boldsymbol{E}\left)\right[{\boldsymbol{\omega }}_{1}^{\mathrm{T}}{\delta }_{1}{\left]\right)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(φ \right(\boldsymbol{E}\left)\right[{\boldsymbol{\omega }}_{1}^{\mathrm{T}}{\delta }_{1}{]}^{\mathrm{T}})+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{C}_{3}}{2}\left(\left|\right|{\boldsymbol{\omega }}_{1}^{}|{|}^{2}+{\delta }_{1}^{2}\right)+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{C}_{1}{\boldsymbol{e}}^{\mathrm{T}}p(\boldsymbol{e}+φ(\boldsymbol{F}\left)\right[{\boldsymbol{\omega }}_{1}^{\mathrm{T}}{\delta }_{1}{]}^{\mathrm{T}}, \alpha ) \end{array}$

(28)

$ \begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}{P}_{2}=\frac{1}{2}\left(φ \right(\boldsymbol{E}\left)\right[{\boldsymbol{\omega }}_{2}^{\mathrm{T}}{\delta }_{2}{\left]\right)}^{\mathrm{T}}{\boldsymbol{W}}_{ij}\left(φ \right(\boldsymbol{E}\left)\right[{\boldsymbol{\omega }}_{2}^{\mathrm{T}}{\delta }_{2}{]}^{\mathrm{T}})+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{C}_{4}}{2}\left(\left|\right|{\boldsymbol{\omega }}_{2}^{}|{|}^{2}+{\delta }_{2}^{2}\right)+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{C}_{2}{\boldsymbol{e}}^{\mathrm{T}}p(\boldsymbol{e}-φ (\boldsymbol{G}\left)\right[{\boldsymbol{\omega }}_{2}^{\mathrm{T}}{\delta }_{2}{]}^{\mathrm{T}}, \alpha )\end{array} $

(29)

其中：$ φ (\cdot ) $表示从实空间$ R $到核特征空间$ H $的映射。

定理2同样适用于非线性加权光滑投影孪生支持向量回归算法，即式（28）和式（29）是连续可微的严格凸函数，亦可采用具有快速收敛能力的牛顿迭代法进行求解。

通过求解式（28）和式（29）可得两个最优超平面的法向量$ {\boldsymbol{u}}_{1}={\left[\begin{array}{cc}{\boldsymbol{\omega }}_{1}^{\mathrm{T}}& {\delta }_{1}\end{array}\right]}^{\mathrm{T}} $和$ {\boldsymbol{u}}_{2}={\left[\begin{array}{cc}{\boldsymbol{\omega }}_{2}^{\mathrm{T}}& {\delta }_{2}\end{array}\right]}^{\mathrm{T}} $。令$ \stackrel{-}{\boldsymbol{x}}=\sum\limits _{i=1}^{n}φ \left({\boldsymbol{x}}_{i}\right)/n $，则非线性加权光滑投影孪生支持向量回归算法的偏置$ {b}_{1} $和$ {b}_{2} $也可以分别用式（25）和式（26）求得，相应的回归函数转换如下：

$ \begin{array}{l}f\left(\boldsymbol{x}\right)=\frac{1}{2}\left({f}_{1}\left(\boldsymbol{x}\right)+{f}_{2}\left(\boldsymbol{x}\right)\right)=\\ \;\;\;\;\;\;\;\;\;\;\;\;\;-\frac{1}{2}{\left(\frac{{\boldsymbol{\omega }}_{1}}{{\delta }_{1}}+\frac{{\boldsymbol{\omega }}_{2}}{{\delta }_{2}}\right)}^{\mathrm{T}}φ \left(\boldsymbol{x}\right)-\frac{1}{2}\left(\frac{{b}_{1}}{{\delta }_{1}}+\frac{{b}_{2}}{{\delta }_{2}}\right)\end{array} $

(30)

1.4 算法步骤

在线性情况下，采用牛顿迭代法训练加权光滑投影孪生支持向量回归算法的过程如下：

输入  惩罚参数$ {C}_{1}\mathrm{和}{C}_{2} $，正则化参数$ {C}_{3}\mathrm{和}{C}_{4} $，不敏感损失参数$ \epsilon $，精度$ {t}_{\mathrm{t}\mathrm{o}\mathrm{l}} $，最大迭代次数$ {i}_{\mathrm{m}\mathrm{a}\mathrm{x}} $和训练数据集矩阵$ \boldsymbol{J} $

输出  回归函数$ f\left(\boldsymbol{x}\right) $

步骤1  通过式（4）计算权值矩阵$ {\boldsymbol{W}}_{ij} $。

步骤2  给定任意初始点$ {\boldsymbol{u}}_{1}^{0} $和$ {\boldsymbol{u}}_{2}^{0} $，并令迭代步骤$ i=0 $。

步骤3  基于$ {\boldsymbol{u}}_{1}^{i} $和$ {\boldsymbol{u}}_{2}^{i} $，通过式（19）和式（20）分别计算$ {\boldsymbol{u}}_{1}^{i+1} $和$ {\boldsymbol{u}}_{2}^{i+1} $直到$ \left|\right|{\boldsymbol{u}}_{1}^{i+1}-{\boldsymbol{u}}_{1}^{i}\left|\right|, \left|\right|{\boldsymbol{u}}_{2}^{i+1}-{\boldsymbol{u}}_{2}^{i}\left|\right| < {t}_{\mathrm{t}\mathrm{o}\mathrm{l}} $或$ i > {i}_{\mathrm{m}\mathrm{a}\mathrm{x}} $，并通过式（25）和式（26）分别计算偏置$ {b}_{1} $和$ {b}_{2} $。

步骤4  通过式（27）计算$ f\left(\boldsymbol{x}\right) $。

非线性情况跟线性情况类似，此处省略。

1.5 时间复杂度分析

本节将WSPTSVR算法的时间复杂度与TSVR^[12]、STSVR^[13]、PTSVR^[14]以及快速聚类加权孪生支持向量回归（Fast Clustering-based Weighted Twin Support Vector Regression，FCWTSVR）算法^[9]进行对比。由于加法的时间复杂度远小于乘法的时间复杂度，因此本文只考虑矩阵乘法运算的次数。另外，核矩阵的求解是所有算法都不可避免的，因此也不作考虑。

PPTSVR、SPTSVR与TSVR一样，在对偶空间中求解二次规划问题，其时间复杂度为$ O\left({n}^{3}\right) $。而WSPTSVR与STSVR一样，在原空间中采用牛顿迭代法进行求解，由于迭代过程中均为对角矩阵相乘和对角矩阵求逆，因此一次迭代的时间复杂度为$ O\left(n\right) $，而牛顿迭代法具有快速收敛能力，因此STSVR与WSPTSVR算法的时间复杂度要小于TSVR及PPTSVR。另外，由于WSPTSVR算法在目标函数中加入了加权项，因此时间复杂度比STSVR稍大。对于FCWTSVR，其时间复杂度虽然也为$ O\left({n}^{3}\right) $，但是还采用了快速聚类的方法剔除异常点，因此时间复杂度要大于PPTSVR、SPTSVR和TSVR。

2 实验结果与分析 2.1 实验设计与参数选择

为了验证WSPTSVR算法的有效性，将其与TSVR、STSVR、PPTSVR、SPTSVR以及FCWTSVR分别在基准测试数据集和人工测试函数上进行比较。首先对数据集进行归一化处理，然后采用标准的五折交叉验证法将数据集划分为训练集和测试集。采用RMSE、MAE、ET这3个性能指标来综合评价回归算法的性能：

$ {R}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt{\frac{1}{n}\sum\limits _{i=1}^{n}({\widehat{\boldsymbol{y}}}_{i}-{\boldsymbol{y}}_{i}{)}^{2}} $

(31)

$ {M}_{\mathrm{M}\mathrm{A}\mathrm{E}}=\frac{1}{n}\sum\limits _{i=1}^{n}|{\widehat{\boldsymbol{y}}}_{i}-{\boldsymbol{y}}_{i}| $

(32)

$ {E}_{\mathrm{E}\mathrm{T}}=\frac{\sum\limits _{i=1}^{n}({\widehat{\boldsymbol{y}}}_{i}-{\boldsymbol{y}}_{i}{)}^{2}}{\sum\limits _{i=1}^{n}({\boldsymbol{y}}_{i}{-\stackrel{-}{\boldsymbol{y}})}^{2}} $

(33)

其中：$ n $为样本个数；$ {\boldsymbol{y}}_{i} $为实际输出；$ {\widehat{\boldsymbol{y}}}_{i} $为$ {\boldsymbol{y}}_{i} $的预测值；$ \stackrel{-}{\boldsymbol{y}}=\sum\limits _{i=1}^{n}{\boldsymbol{y}}_{i}/n $为$ {\boldsymbol{y}}_{i} $的平均值。

通常地，RMSE的值越小，算法的预测性能越好；MSE的值越小，算法的预测误差越小；ET的值越小，表明预测值与真实值的一致性越好。一个好的回归算法应该综合考量RMSE、MAE和ET的值。

此外，还分别统计了各算法的CPU运行时间。所有实验都是在MATLAB 2019a平台上用Intel i5-9400F@2.90 GHz处理器、16 GB内存的计算机完成的，并且所有的实验结果均取20次独立运行结果的平均值。

参数选择与算法性能密切相关，实验采用网格搜索法选取最优参数。由于$ {ε}_{1} $和$ {ε }_{2} $的设置对预测性能不会产生很大影响^{[13, 22-23]}，因此将$ {ε}_{1} $和$ {ε }_{2} $的值都设置为0.01。在TSVR中，令$ {C}_{1}={C}_{2} $并且在集合$ \left\{{2}^{i}|i=-8, \cdots , 0, \cdots , 8\right\} $中进行寻优。为了保证算法对比的公平性^{[9, 24]}，令STSVR、PPTSVR、SPTSVR和WSPTSVR算法中$ {C}_{1}={C}_{2} $、$ {C}_{3}={C}_{4} $，且在集合$ \left\{{2}^{i}|i=-8, \cdots , 0, \cdots , 8\right\} $中进行寻优。此外，在STSVR和WSPTSVR算法中，设置精度t_tol=10^-6，最大迭代次数$ {i}_{\mathrm{m}\mathrm{a}\mathrm{x}}=50 $。在非线性实验中，核函数统一采用高斯径向基函数$ K({\boldsymbol{x}}_{i}, {\boldsymbol{x}}_{j})=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{‖{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}‖}^{2}/2{\sigma }^{2}\right) $，参数$ \sigma $在集合$ \left\{{2}^{i}\right|i=-5, \cdots , 0, \cdots , 5\} $中进行寻优。

2.2 基准测试数据集

为了对比各个算法的综合性能，将其在7个基准测试数据集上进行测试，使用的基准测试数据集如表 1所示，其中样本数最小为103，最大为9 568，它们可以从UCI机器学习库下载。

表 2统计了各个算法在7个基准测试数据集上的实验结果，为了清楚起见，最优结果加粗表示。

由表 2可以看出：1）与PPTSVR相比，WSPTSVR算法在大多数数据集上有着相似或者更小的RMSE、MAE和ET，这是由于WSPTSVR算法采用孤立森林法赋予样本中的异常点很小的权值，并通过异常分数赋予每个样本点不同的权值，能够有效地降低潜在噪声或异常点对回归超平面的影响，从而提升算法的预测性能；2）与采用聚类方法去除异常点的FCWTSVR相比，WSPTSVR算法在4个数据集上的RMSE值最优，在3个数据集上的MAE值最优，而FCWTSVR的RMSE、MAE和ET值均在3个数据集上最优，表明WSPTSVR算法的预测性能与当前提出的回归算法相比具有可比性。

在算法的训练时间方面，WSPTSVR算法的时间复杂度要小于TSVR、PPTSVR以及SPTSVR，实验结果也表明，WSPTSVR算法的训练速度快于TSVR、PPTSVR以及SPTSVR。但是由于WSPTSVR算法在目标函数中引入了权值矩阵，因此与STSVR相比，训练速度稍慢。而FCWTSVR由于采用了快速聚类算法剔除异常点，其时间复杂度在6种算法中是最大的，实验结果也表明，FCWTSVR的训练速度最慢。

2.3 人工测试函数

为了进一步验证WSPTSVR算法在非线性情况下的拟合性能，在$ \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\boldsymbol{x}\right) $函数上进行实验，并且人为添加两种不同类型的噪声。$ \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\boldsymbol{x}\right) $函数的具体形式如下：

$ {\boldsymbol{y}}_{i}=\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left({\boldsymbol{x}}_{i}\right)+{n}_{i}=\frac{\mathrm{s}\mathrm{i}\mathrm{n}\left({\boldsymbol{x}}_{i}\right)}{{\boldsymbol{x}}_{i}}+{n}_{i},   {\boldsymbol{x}}_{i}\in [-3\mathrm{\pi }, 3\mathrm{\pi }] $

(34)

其中：$ {\boldsymbol{x}}_{i} $表示输入；$ {\boldsymbol{y}}_{i} $表示对应于$ {\boldsymbol{x}}_{i} $的输出；$ {n}_{i} $表示噪声。两种不同类型噪声的具体形式如下：

$ {n}_{{i}_{1}}=\left(-\frac{\left|{\boldsymbol{x}}_{i}\right|}{8\mathrm{\pi }}+0.5\right)\times {\xi }_{i}, {\xi }_{i}\sim U[-\mathrm{0.1, 0.1}] $

(35)

$ {n}_{{i}_{2}}=\left(-\frac{\left|{\boldsymbol{x}}_{i}\right|}{8\mathrm{\pi }}+0.5\right)\times {\xi }_{i}, {\xi }_{i}\sim N[\mathrm{0, 0}.{1}^{2}] $

(36)

其中：$ U[-\mathrm{0.1, 0.1}] $表示在闭区间$ [-\mathrm{0.1, 0.1}] $内服从均匀分布；$ N[\mathrm{0, 0}.{1}^{2}] $表示服从均值为$ 0 $、方差为$ 0.{1}^{2} $的高斯分布。

随机生成47个混入噪声的独立样本作为训练样本和150个混入噪声的独立样本作为测试样本。另外，在训练样本中人为加入3个异常点。

表 3给出了2种不同类型噪声下6种回归算法的平均RMSE、MAE、ET以及CPU运行时间，最优结果加粗表示。

从表 3可以看出：当样本中没有异常点时，FCWTSVR和WSPTSVR的RMSE、MAE以及ET值要小于其他4种回归算法，表明FCWTSVR和WSPTSVR算法具备更好的预测性能和抗噪声能力；当样本中存在异常点时，WSPTSVR算法的RMSE、MAE和ET值均明显小于PPTSVR，以RMSE为例，WSPTSVR算法比PPTSVR平均约小44.2%，表明本文算法在样本中存在异常点时，有着更好的预测性能。这是由于本文算法采用孤立森林法赋予异常点很小的权值，因此受异常点的影响较小。另外，FCWTSVR采用聚类的方法剔除样本中的异常点，同样获得了较好的预测性能。在训练时间方面，由于WSPTSVR算法直接在原空间中采用牛顿迭代法进行求解，训练速度比TSVR、PPTSVR和SPTSVR更快，而与STSVR相比，由于目标函数中添加了权值矩阵，因此训练速度稍慢。

图 1给出了高斯分布噪声下不同回归算法在$ \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(x\right) $函数上的拟合曲线。从图 1可以看出，与PPTSVR相比，WSPTSVR算法在噪声和异常点的干扰下更接近真实曲线，拟合效果更好。这是由于WSPTSVR算法采用孤立森林法判断样本中的异常点，并利用样本的异常分数赋予噪声和异常点很小的权值，从而获得较好的拟合效果。反观TSVR、STSVR、PPTSVR和SPTSVR在异常点处的拟合曲线扭曲较为明显，受异常点的影响较大，拟合效果较差。原因是这4种回归算法赋予所有样本点相同的权值，而异常点的存在迫使拟合曲线偏向异常点，从而使得拟合效果变差。另外，FCWTSVR采用聚类的方法确定并剔除异常点，因此在异常点处也获得了较好的拟合效果。

4 结束语

本文提出一种加权光滑投影孪生支持向量回归算法。采用孤立森林法判断样本中潜在的异常点，并通过赋予潜在的异常点很小的权值，有效地削弱了其对超平面构造的影响。引入Sigmoid光滑函数，通过在原空间中采用牛顿迭代法求解无约束优化问题，获得了比双边移位投影孪生支持向量回归算法更快的训练速度。实验结果表明，与其他代表性回归算法相比，该算法受异常点的影响更小，拟合效果更佳。下一步将从寻找逼近能力更优的光滑函数入手，使WSPTSVR算法获得更好的拟合性能。