基于全变差和欧拉弹性杆的有监督模式识别方法

摘要

本发明公布了一种基于全变差和欧拉弹性杆的有监督模式识别方法，首先在所述最小二乘正则化框架下构造出基于全变差和欧拉弹性杆的能量泛函；然后利用变分法将该能量泛函最小化问题转换为求解对应的欧拉-拉格朗日偏微分方程；对所述偏微分方程求解，进而得到最终的分类器；运用该分类器对数据进行模式识别。本发明对有监督模式识别问题提出了新的方法，可以应用于一般情况下的分类问题，例如手写体数字识别；在绝大多数数据集上，本发明提出的方法能达到与现有的流行方法相媲美的效果。

说明书

技术领域

本发明属于模式识别技术领域，具体涉及将图像处理中基于全变差和欧拉弹性杆的模型 应用于能量最小化的最小二乘正则化(RLS)框架，用于有监督的模式识别。

背景技术

有监督学习技术是模式识别领域中的一大热门问题，在计算机视觉等多个领域有着广泛 的应用。模式识别究其本质研究的是分类技术，在不产生歧义的情况下，本说明书中用分类 这个术语来指代模式识别。有监督分类指的是利用有标签的训练样本集训练出一种规则(分 类器)来对新的样本进行预测，其一般形式描述为：

给定训练集{(x₁，y₁)，...(x_n，y_n)}，其中x_i∈R^d，R^d为d维实向量空间，对于二类分类问题 y_i∈{+1，-1}。其实质是要学习一个从向量x到标签y的映射，因此可以把分类问题归到一个 函数学习的框架下来考虑。一般的有监督分类的正则化框架可以表示为：

$\underset{u}{\min }{\Sigma }_{i=1}^{n}L(u\left(x_i\right),y_i)+\mathrm{\lambda S}\left(u\right)$

其中u(x)是要学习的映射函数，n是训练样本个数，L代表损失函数，S(u)是正则项用 来刻画映射u(x)的某些性质，λ是参数，用来调节损失项和正则项的比重。当L取预测值u(x) 和真实值y之差的平方时，就是最小二乘正则化框架，本发明就是基于这个最小二乘正则化 的连续情况下的框架。

现有的一些方法如SVM，采用的hinge损失函数，正则项惩罚的是分界面的间隙宽度， 但其仍是基于统计学习理论。常用的贝叶斯推理等统计学方法虽然取得了丰硕的成果，但也 逐渐暴露出某些缺点，比如参数过多从而导致维数灾难，经典的统计分布函数与真实数据的 实际分布不一致，难以有效估计数据的概率密度函数等。目前基于泛函能量极小化的函数学 习方法越来越受到重视，R.M.Rifkin在2003年提出基于Tikhonov正则化的最小二乘分类(RLSC)； 2010年Kush R等利用水平集的方法做有监督分类的研究，借鉴图像分割中的水平集方法， 将决策边界(水平集函数的零水平集)的表面积作为正则化能量泛函，都取得了不错的效果。

本发明从基础数学的角度出发，在一个连续的框架下建立能量泛函正则化模型，并利用 泛函分析变分法等工具，采用合理的求解方法来实现有监督分类。

在图像处理中，全变差图像去噪模型和欧拉弹性杆图像修补模型都取得了很大的成功， 均是通过建立泛函模型，利用变分法等数学手段来求解。可以把图像去噪、图像修补理解为 在已知的数据点上已有标签值(像素值)，要求得一个函数(要求的图像)来预测未知的点或 有污染的点上的值，将这种二维的情形上升到高维，就类似于分类问题。本发明将这两种模 型结合正则化的能量最小化框架来做有监督分类，这两种模型给出的正则项在不同程度上刻 画了所求函数的几何性质，使得从另一种全新的角度来解决分类问题。

发明内容

本发明的目的在于提出一种新的基于正则化框架的有监督模式识别方法。

本发明的技术方案如下：

一种基于全变差和欧拉弹性杆的有监督模式识别方法，该方法基于连续的最小二乘正则 化框架：其中x是用于模式识别的一个样本的特征向量表示，u(x) 为要求的函数，y为训练数据集中样本x的类标签值，Ω为训练样本所处的定义域，λ是参 数，S(u)为正则项，

所述正则项S(u)采用全变差和欧拉弹性杆的形式，从而得到新的有监督模式识别方法。

所述的有监督模式识别方法，包括如下步骤(参图1)：

1)首先在所述最小二乘正则化框架下构造出基于全变差和欧拉弹性杆的能量泛函；

2)然后利用变分法将该能量泛函最小化问题转换为求解对应的欧拉-拉格朗日偏微分方 程；

3)对所述偏微分方程求解，进而得到最终的分类器；

4)运用该分类器对数据进行模式识别。

所述的有监督模式识别方法，步骤3)中，采用梯度下降的时间搜索算法对所述偏微分 方程进行求解，即增加一个虚拟的时间变量，采用梯度下降法对所述偏微分方程进行迭代求 解。

所述的有监督模式识别方法，步骤3)中，采用滞后的线性方程组迭代法对所述偏微分 方程进行求解，即固定住优化代价函数中的非线性部分，使其变成一个线性方程组，对线性 方程组求解后再更新之前滞后的非线性部分，进行反复迭代求解。

所述的有监督模式识别方法，步骤1)中所述全变差能量泛函为：

$J_{\mathrm{TV}}\left[u\right]={\int }_{\Omega }{(u-y)}^2\mathrm{dx}+\lambda {\int }_{\Omega }|▿u|\mathrm{dx}$

其中，x是用于模式识别的一个样本的特征向量表示，u(x)为要求的函数，y为训练数据 集中一个样本的类标签值，J_TV[u]表示关于u(x)的基于全变差的能量泛函，第二项为在u(x)连续的情况下的全变差，Ω是自变量所在的定义域，是梯度算子，代表了 函数u(x)的梯度的L2模范数，λ是可预先设定的正则化参数。

所述的有监督模式识别方法，步骤1)中所述欧拉弹性杆的能量泛函为：

$J_{\mathrm{EE}}\left[u\right]={\int }_{\Omega }{(u-y)}^2\mathrm{dx}+\lambda {\int }_{\Omega }(a+b{\kappa }^2)|▿u|\mathrm{dx},$$\kappa =▿·\left(\frac{▿u}{|▿u|}\right)$

其中，x是用于模式识别的一个样本的特征向量表示，u(x)为要求的函数，y为训练数 据集中一个样本的类标签值，J_EE[u]代表关于u(x)的基于欧拉弹性能量的能量泛函，κ为函 数u的等值超曲面(水平集)的平均曲率，Ω是自变量所在的定义域，是函数u的梯度， 代表了函数u(x)的梯度向量的模范数，是散度算子的形式表达，λ是可预先设定的 一个正则化参数，a与b是控制全变差与欧拉弹性能量贡献度的两个参数。

本发明的有益效果：本发明对有监督模式识别问题提出了新的方法，可以应用于一般情 况下的分类问题，例如手写体数字识别；在绝大多数数据集上，本发明提出的方法能达到与 现有的流行方法相媲美的效果。

附图说明

图1是本发明的流程图。

图2是USPS手写体数字集中的样本示意图及分类结果示意。

具体实施方式

本发明的实施方式如下：

实施例一：USPS手写体数字识别

USPS(U.S.Postal Service)数据集使用美国邮政信封上扫描的手写体数字图片，每张图片 是16*16的灰度图片，包含一个数字，如图2(a)给出了其部分数据展示，从数字0到9， 共10类。本实施例从中随机抽取了1000个样本来做实验，由于原始数据维数较高，用主成 分分析方法(PCA)将其降至30维，并归一化到[0，1]区间。

步骤1：在正则化框架下构造能量泛函

(a)全变差能量泛函

数学上，函数f(x)的全变差(Total Variation)在一维情况下并且当函数连续时有形 式：表示f(x)的全变差，a和b为区间的端点，f′(x)为函数 的导数。可以看到全变差是对于函数的所有变化的总和的度量。

将关于函数u(x)的全变差作为泛函能量的正则项，得到的有监督分类的能量泛函为：

$J_{\mathrm{TV}}\left[u\right]={\int }_{\Omega }{(u-y)}^2\mathrm{dx}+\lambda {\int }_{\Omega }|▿u|\mathrm{dx}$

其中J_TV[u]表示关于u(x)的基于全变差的能量泛函(TV是Total Variation的缩写)，第 二项为在u(x)连续的情况下的全变差，Ω是自变量所在的区域，是梯度算子(例 如代表了函数u(x)的梯度的L2模范数，λ是参数。在图像去噪 中，全变差正则项的作用在于保持图像中物体边缘处的信息。一般的去噪方法往往会使得图 像边缘处也被模糊了，而全变差模型可以抑制这种模糊。这里在监督分类的能量泛函中，全 变差项不仅限制函数变化总和较小，使得在分类中可得到一个大部分区域平缓的函数，而且 在类与类的边界处允许有较大的变化，来使分类的边界更明显。最小化这个关于u的能量泛 函得到的函数即可给出最终的分类器，函数的零水平集为分类边界。

(b)欧拉弹性杆能量泛函

将欧拉弹性杆的模型推广到高维的有监督分类中，不像全变差那样直接，图像修补中的 模型为：

$\underset{I}{\min }\lambda {\int }_{\Omega }(a+{\mathrm{b\kappa }}^2)|▿I|\mathrm{dxdy}+{\int }_{\Omega D}{(I-I_0)}^2\mathrm{dxdy}$

其中I为所求图像，I₀为原始要修补图像，Ω为整个图像区域，D为要修补区域，a、b 为参数，κ为图像中等高线的曲率，形式上为：

$\kappa =▿·n,n=\frac{▿I}{|▿I|}$

为散度算子，n为等高线的单位外法向量，为图像的梯度，为图像的梯度模， 此处是L2模范数。欧拉弹性杆(Euler’s Elastica)是指达到如下弹性能量平衡态的曲线γ： E₂[γ]＝∫_γ(a+bκ²)ds，ds为弧长，κ为γ的曲率。修补模型中的第一项是对图像中所有等 高线的欧拉弹性能量的积分。将欧拉弹性杆运用于图像修补，其优势在于弹性杆的插值能力， 这种非线性样条是填补缺失区域的有力工具。把这种弹性能量推广到高维时，需要考虑曲率κ 推广到高维的情形，在二维情况下u(x，y)＝0(此处x，y为二维自变量)，决定一条等高线 y＝u(x)，而其曲率则具有形式在三维情况下u(x，y，z)＝0(x，y，z为三维自变 量)决定的是一个等值曲面，而该等值曲面的平均曲率具有和一维情况下曲线曲率相同的形 式。曲面的平均曲率是曲面的外在度量，它刻画了曲面嵌入到周围空间的弯曲程度。更高维 的情况，由u＝0决定的超曲面是函数u的等值超曲面，欧拉弹性能量迁移到高维就是刻画这 些超曲面的弹性能量。表1给出了从u的表达式通过隐函数定理决定的曲线或(超)曲面， 以及他们的曲率形式，其(平均)曲率除了一个和维数有关的常数因子外具有相同的形式。

表1由隐函数定理决定的曲线(超)曲面及其曲率

由此本发明给出了欧拉弹性杆正则项的函数学习能量泛函：

$J_{\mathrm{EE}}\left[u\right]={\int }_{\Omega }{(u-y)}^2\mathrm{dx}+\lambda {\int }_{\Omega }(a+b{\kappa }^2)|▿u|\mathrm{dx},$$\kappa =▿·\left(\frac{▿u}{|▿u|}\right)$

其中，J_EE[u]代表关于u(x)的基于欧拉弹性能量的能量泛函(EE是Euler’s Elastica的缩 写)，κ为函数u的等值超曲面的平均曲率。为了形式的简化，在这里仍然保留原始模型中κ 的形式，省去了前面的系数因子在计算时将这个因子添加到权重b上。其余记号含义 与之前TV的能量泛函中类似。能量泛函中第一项平方损失项限制函数尽量接近训练数据的真 实值，正则项限制函数u等值超曲面的弹性能量使其尽量的小，这就意味着最终得到的分类 边界有最低的弹性能量。如果令EE能量泛函中的a＝1，b＝0，这个能量泛函就退化成为没有 考虑由曲率项带来的能量的TV能量泛函。

步骤2：利用变分法将能量最小化转为偏微分方程

利用变分法可以将TV和EE的能量泛函最小化问题转变成求解对应的欧拉-拉格朗日偏微 分方程(PDE)。全变差能量泛函对应的PDE为：

$\lambda ▿·\left(\frac{▿u}{|▿u|}\right)-2(u-y)=0---\left(1\right)$

注意此处为散度算子，是函数u的梯度，代表了函数u(x)的梯度向量的模范数。 这是一个非线性椭圆方程。对于欧拉弹性杆能量泛函，对应的PDE为：

$\left(\begin{array}{c}\lambda ▿·V-2(u-y)=0\\ V=\phi \left(\kappa \right)\frac{▿u}{|▿u|}-\frac{1}{|▿u|}▿\left({\phi }^{\prime }\left(\kappa \right)\right|▿u\left|\right)+\frac{1}{{|▿u|}^3}▿u(▿u·▿\left({\phi }^{\prime }\left(\kappa \right)\right|▿u\left|\right))\end{array}\right)---\left(2\right)$

其中采用记号φ(κ)：＝a+bκ²，φ′(κ)表示关于曲率κ的导数。

步骤3：采用两种方法来求解，得到最终的分类器

偏微分方程在高维的情况下若采用常规的数值方法(例如有限差分、有限元)会遇到维 数灾难，而且上述的偏微分方程(1)和(2)形式十分复杂，数值上暂时没有适用的方法。 本发明采用函数近似的思想，将函数泛函最小化转化为关于线性系数的最小化问题，并采用 梯度下降的时间搜索算法(GD)和滞后的线性方程组迭代法(IagLE)两种方法来求解。

利用函数近似的思想，将函数u(x)表示成一些径向基函数的线性组合：

其中m为基函数的个数，w_i为线性系数，为高斯径向基函数，c和z_i均为参数， ||x-z_i||为向量的L2范数。令w＝(w₁，w₂，...，w_m)^T，于是求解函数u(x)的问题变为求解线性系 数w。在实施时基函数中的z_i取所有的训练样本点。

两种方法求解系数w：

(a)梯度下降的时间搜索算法(GD)

由于全变差方法是欧拉弹性能量方法的一个特例，从欧拉弹性能量方法的求解出发，全 变差的求解是类似的。引入时间变量t，从变分法得到的偏微分方程(2)中，得到函数u的 梯度为：

$\frac{\partial u(x,t)}{\partial t}=\lambda ▿·V(x,t)-2(u(x,t)-y)$

将u表示成矩阵乘以向量的形式：u＝Φw，其中Φ为矩阵，其元素为将Φ 固定，把u的梯度转换为系数向量w的梯度，于是梯度下降算法的迭代公式变为：

$W^{(k+1)}=w^{\left(k\right)}-\tau \frac{{\partial W}^{\left(k\right)}}{\partial t}=W^{\left(k\right)}-\tau {\Phi }^{-1}\left(\begin{array}{c}\frac{{\partial u}^{\left(k\right)}}{\partial t}|x=x_1\\ ·\\ ·\\ ·\\ \frac{{\partial u}^{\left(k\right)}}{\partial t}|x=x_n\end{array}\right)---\left(4\right)$

其中τ为时间步长，公式中的第一项由如下形式计算(经过一系列推导并略去u 的三阶及以上导数)：

$▿·V=\kappa -\frac{4\mathrm{b\kappa \Delta u}}{{|▿u|}^3}(▿\mathrm{uH}\left(u\right))·\frac{▿u}{|▿u|}{+\mathrm{b\kappa }}^3+2b(\frac{2\mathrm{\Delta u}}{{|▿u|}^5}+\frac{\kappa }{{|▿4|}^4})(▿\mathrm{uH}\left(u\right))·(▿\mathrm{uH}\left(u\right))$

$+2b\{-\frac{3}{{|▿u|}^5}(▿\mathrm{uH}\left(u\right))·▿u+\frac{▿u}{{|▿u|}^3}\}(-\frac{2\mathrm{\Delta u}}{{|▿u|}^2}+\frac{\kappa }{|▿u|})▿u·(▿\mathrm{uH}\left(u\right)).$

其中Δ是拉普拉斯算子，H(u)是u的Hessian阵。当b取零 时，就可得到全变差的下降公式。按照上述迭代公式(4)，初值选为w⁽⁰⁾＝(Φ^TΦ+ηI)^-1Φ^Ty， η为参数，y＝(y₁，y₂，...，y_n)^T，每步更新w^(k)直至其收敛。

(b)滞后的线性方程组迭代法(IagLE)

由于TV的欧拉拉格朗日偏微分方程(1)具有相对较简单的形式，首先从TV的求解来推 导。将PDE公式(1)中的非线性项的散度算子展开得到：

$-\frac{\lambda }{|▿u|}(\mathrm{\Delta u}-\frac{▿u^{T}H\left(u\right)▿u}{▿u^T▿u})+2(u-y)=0$

将u的高斯径向基函数组合的形式公式(3)代入上式，并将参数z_i用所有训练样本点x_i代替，此时基函数的个数m就等于训练样本的个数n，经过一系列推导后得到如下的方程组：

$\left(\begin{array}{c}\underset{i-1}{\overset{n}{\Sigma }}{w}_i{\phi }_i\{d-1-c{|x-x_i|}^2+c\frac{g^T(x-x_i){(x-x_i)}^{T}g}{g^{T}g}+\frac{\left|g\right|}{\lambda }\}=\left|g\right|y\\ g=\underset{i-1}{\overset{n}{\Sigma }}{w}_i{\phi }_i(x-x_i)\end{array}\right)---\left(5\right)$

其中d为数据的维数，|·|为向量的L2范数，把g固定可以得到关于w的线性方程组。于 是利用滞后的思想，本发明提供了另外一种迭代解法，称之为滞后的线性方程组迭代法(IagLE)：

1)给w一个随机的初值w⁽⁰⁾然后计算出g⁽⁰⁾。

2)利用最小二乘方法求解关于w的线性方程组(5)，给出一个新的w^(k)并计算g^(k)。

3)迭代执行第2)步直至w收敛。

对于EE方法具有复杂的欧拉拉格朗日方程(2)，将含有曲率的项固定，记为K＝a+bκ²， 于是可以得到类似的方程组：

$\left(\begin{array}{c}\underset{i=1}{\overset{n}{\Sigma }}{w}_i{\phi }_i\{d-1-c{|x-x_i|}^2+c\frac{g^T(x-x_i){(x-x_i)}^{T}g}{g^{T}g}+\frac{\left|g\right|}{\mathrm{\lambda K}}\}=\left|g\right|y\\ g=\Sigma w_i{\phi }_i(x-x_i)\\ K=a+b{(\frac{1}{\left|g\right|}\Sigma w_i{\phi }_i(1-d+c{\left|\right|x-x_i\left|\right|}^2-c\frac{g^T(x-x_i){(x-x_i)}^{T}g}{g^{T}g}))}^2\end{array}\right)$

将g和K都固定即可得到关于w的线性方程组。同样给w一个初值，用与TV相同的迭 代方法，每步更新g和K直至w收敛。

在求解线性方程组的时候，假设该线性方程组形式上为Ψw＝y，Ψ为形式上简化的矩 阵。由于w的模有时会非常大使得线性方程组变成坏条件的，一个比较好的求解策略是求解 如下的正则化的最小二乘问题：

$\underset{w}{\min }\left|\right|\mathrm{\Psi w}-y{\left|\right|}_2^2+\eta {\left|\right|w\left|\right|}_2^2$

其中η为正则化参数，其解为w＝(Ψ^TΨ+ηI)^-1Ψ^Ty。于是滞后的线性方程组迭代法中有三 个参数：高斯径向基函数中的参数c，能量泛函中的正则化因子λ和求解线性方程组的参数η。

将得到的w代入基函数的线性组合表达式(3)，即可得到函数u的近似解析表达式。该 近似解析表达式的零水平集为二类分类边界，最初得到的是一个二类分类器，即在预测时， 将样本代入u中，若函数值大于零则预测标签为1，小于等于零则标签为-1。而此实施例中 的手写体数字共10类，是个多类分类问题，采用一种常用的一对多的策略，即训练10个二 类分类器，再将这10个二类分类器组合得到最终的分类器。图2(b)表示的是对于图2(a) 的识别结果。

在求解过程中对于EE方法固定b＝0.01，于是对于TV和EE的GD解法均只有两个参数c 和λ，lagLE解法有三个参数c，λ和η，此处采用对数据集五折交叉验证来选择最优参数， 并与SVM和BPNN(反向传播神经网络算法)进行分类效果的对比。表2给出的是各种方法 的准确率对比，此处的准确率是各方法五折交叉验证中最优参数对应的最优准确率。从表2 中可以看出本发明提出的基于全变差和欧拉弹性能量的两种方法在USPS数据集上效果超过 了SVM和BPNN。

表2各种方法在USPS数据集上准确率的对比结果(％)

实施例二：应用于其他一些常用分类数据

本实施例基于Iibsvm网站以及UCI机器学习库中的8个分类的数据集，包括两类及多类 数据。对于所有数据集都归一化到[0，1]区间，操作步骤与实施例一中相同。对于所有方法都 仍然用五折交叉验证来选取最优参数，各方法的最优准确率对比如表3所示：

表3各方法的分类准确率对比(％)

从表3中可以看出本发明提出的TV和EE模式识别方法的IagLE解法在总体上比GD的要 好，EE的效果略胜于TV的效果，他们在所有数据集上都比BPNN的准确率要高。EE方法在 IagLE解法下在六个数据集上超过了SVM，在其余两个数据上的效果与SVM也相差不大，可 以看出本发明提出的基于全变差和欧拉弹性能量的模式识别方法可以与现在很成熟的SVM算 法相媲美。