数学专业毕业论文-基于Wasserstein 距离的图像识别算法研究

目录致谢 I目录 III插图索引 V算法索引 VII摘要 IXABSTRACT XI第一章绪论 11.1引言 11.2相关的工作 21.3论文结构 3第二章Wasserstein度量 52.1 Wasserstein度量的定义 52.2 Wasserstein度量与最优传输问题之间的关系 52.3与KL散度的区别和联系 6第三章基于Wasserstein度量聚类的图像分类识别算法 73.1转换图片为待聚类的经验分布 73.2基于Wasserstein度量的K-means聚类算法 73.3目标优化问题模型的建立 9第四章WDP和WBP问题的求解和主要结果分析 114.1基于正则W距离求解WDP和WBP 114.2基于Bregman-ADMM算法求解WDP和WBP 134.3两种算法的结果比较和分析 164.4基于聚类的图像识别算法结果展示 20第五章总结和展望 235.1总结 235.2未来的展望 23参考文献 25III中国科学技术大学学士学位论文插图索引1.1 Google以图找图功能. . . . . . . . . . . . . . . . . . . . . . . . . . 12.1 Euclidean插值(左)与Wasserstein插值(右)(参考来源:1) . . . . . . 64.1测试用输入图片. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184.2同一WBP问题两种算法的比较. . . . . . . . . . . . . . . . . . . . 194.3测试用三维体素. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194.4同一三维WBP问题两种算法的结果比较. . . . . . . . . . . . . . . 194.5聚类程序输入. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204.6 K = 3聚类结果示意. . . . . . . . . . . . . . . . . . . . . . . . . . . 204.7聚类结果示意:数字2 . . . . . . . . . . . . . . . . . . . . . . . . . . 214.8聚类结果示意:数字4 . . . . . . . . . . . . . . . . . . . . . . . . . . 214.9聚类结果示意:数字6 . . . . . . . . . . . . . . . . . . . . . . . . . . 214.10聚类结果示意:数字2;3;4;5;6 . . . . . . . . . . . . . . . . . . . . . 214.11聚类结果和真实标签的可视化比较. . . . . . . . . . . . . . . . . . 21V中国科学技术大学学士学位论文算法索引3.1图像分类算法框架. . . . . . . . . . . . . . . . . . . . . . . . . . . . 73.2基于Wasserstein度量的K-means算法. . . . . . . . . . . . . . . . . 83.3基于Wasserstein度量的K-means+聚类初始化算法. . . . . . . . . 84.1 Sinkhorn迭代. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124.2 WBP问题求解算法一. . . . . . . . . . . . . . . . . . . . . . . . . . . 134.3 BADMM求解WBP问题. . . . . . . . . . . . . . . . . . . . . . . . . 17VII中国科学技术大学学士学位论文摘要Wasserstein度量是一种用于衡量概率测度间差异程度的度量，该度量相比于信息论中的KL散度，满足度量性质，且求解插值问题时可以保持分布函数的几何性质。随着近些年数据驱动算法的火热，模式识别问题被广泛研究，图像识别和分类课题迅速发展。传统的欧式度量不能保留图像几何性质的缺点在复杂和真实的识别环境中逐渐暴露出来，而Wasserstein度量既可以保留分布的几何属性，又可以作为一种度量应用于统计学习的模型中，从而开始受到关注。本文中，我们主要聚焦该度量在图像识别算法中应用可能性的探究。我们提出将K-means聚类算法中的欧式度量计算替换为Wasserstein度量的计算，由此使得对图像的直接聚类成为可能。而算法框架中引出的两个最优化问题:Wasserstein度量计算和Wasserstein重心的高效求解，是我们算法设计中的一大障碍：与其表述等价的最优运输问题很难又快又好地求出精确解。针对这一困难，我们使用了两种最优化的方法求解：第一种是通过给定一个正则项修正，使得该问题可以借由Sinkhorn迭代和次梯度迭代求得近似解；第二种则是利用交替乘子下降法(BADMM)分解问题，近似求解与原问题等价的最优化问题。最后，为了验证我们的算法的实际效果，我们实现了若干例子，并作了相关分析。关键词：Wasserstein度量Wasserstein重心最优传输问题图像识别和聚类K-means Sinkhorn迭代BADMMIX中国科学技术大学学士学位论文ABSTRACTGiventwoprobabilitymeasures,Wassersteinmetricisknownasamethodtojudgethe dissimilarity between these two measures. Traditional statistics are interested in itsproperty to remain geometry information while solving the interpolation-type problem.With the rapid development in Statistical Learning for Image Recognition and Cluste-ring Issues, the Euclidean metric and its related mutation cannot appropriately describethe information from image when it comes to a complex scene and more realistic sce-nario. However, Wasserstein metric can attain the geometric information in images.Therefore, researchers are focusing on Wasserstein metric and to see how it can becombined with traditional learning method, or apply in metric learning schema.In this thesis, we proposed the possibility to use Wasserstein metric directly inclustering method,i.e. K-means, by replacing Euclidean metric to Wasserstein Metric.This can be applied to the recognition problem, a hot issue in CV and CG. To solve outWasserstein Distance and Barycentric Problem, we follow the former research and givetwo algorithm to overcome the computational intractability. One for sinkhorn iterationwith Entropic Regularized Wasserstein Distance, and one for Bregman Alternative Di-rection Method of Multipliers(BADMM). Finally, we analyze our result and verifiedout method by various examples.Keywords: Wasserstein Metric, Wasserstein Barycentric Problem, Optimal Transport,Image Clustering, K-means, Sinkhorn Iteration, Bregman-ADMMXI中国科学技术大学学士学位论文第一章绪论1.1引言随着时代和社会的发展，越来越多的公司和学者开始关注大数据时代下，基于数据驱动的图形学和计算机视觉的算法。图像和三维模型的识别(Recognition)问题近几年逐渐火热起来。对于图像识别算法的研究也在推动当今社会的发展，在人脸识别、自动驾驶、医疗、艺术设计等诸多领域，都有识别算法的应用。例如下图就反映了生活中的图像识别的一个应用:以图找图检索功能。(a)检索图片(b)以图找图结果图1.1 Google以图找图功能在模式识别的问题中，无论是有监督地统计学习，还是图片与数据库的检索和匹配问题，其中都要涉及到一个图片与图片之间，或者二者所提取的特征量之间比较的问题。用数学的语言描述，若我们认为fxg与fyg是两簇不同分类的图像的集合，则f(x) - f(y)用于反映图像之间的相似程度或者差异。这里的映射f是一个抽象概念，表明对图像数据的处理。它既可以是对高维特征向量的降维(如PCA,t-SNE)，也可以是对图像的特征的描述(如利用卷积神经网络提取特征)。在很长一段时间内，科学家们在这个具体的工程问题上，主要关注如何设计更好的统计模型，对这样这样一个高维数据进行合理表达，如采用SIFT角点提取特征点，边缘检测提取轮廓，或者关注纹理和颜色直方图中的信息，并结合统计学和深度学习的方法加以对图像进行区分，取得了很多的进展，尤其是以近几年火热的卷积神经网络架构为首的一批深度学习模型，在固定数据集的学习问题上取得了其他方法难以匹敌的准确率。但是，这些方法在实际的应用仍然存在一定问题。如图1.1所示，虽然以图找图功能能找到完全匹配的图片(上图未放出)，但很明显可以看出，找图结果大部分依赖的是颜色，而不是图片中八万的字样。随着未来图像输入数据的复杂化(如室内场景的照片)，几何信息的重要性也会越来越高，这是因为图像和三维1中国科学技术大学学士学位论文模型作为具有高维特征的一种数据，除了具有一般数据中点的信息特征(如颜色，亮度)，还具有点与点，边与边之间的几何信息特征。但是，在过去的研究中，只要是基于欧式度量体系的特征提取或者学习算法，都不能在特征上有效地保留几何信息。即使是实际结果优秀的神经网络，其中间层的可视化结果也暗示了它对几何特征并不敏感。由此，我们自然提出一个问题：能否定义或寻找一种度量，使得它在计算过程中可以保持或者充分利用几何信息，且这样的度量可以直接计算图像和图像之间的差异性。Hausdorff距离(详见2 )虽然可以被应用于图像的比较和识别问题，但它的仍然不是保几何的。所以我们决定探究Wasserstein度量在相关问题中应用的可能性。本文将聚焦Wasserstein度量在图像识别算法中的应用。针对传统欧式度量的局限性，我们利用Wasserstein度量取代传统的基于欧式度量的l1,l2范数。该度量不再简单地把图像看作一个普通的高维特征向量，而是将其看为一个统计量，这样在求解过程中可以保留样本分布的几何信息；之后，在计算此度量的基础上，修改最简单的聚类算法，从而提出一种设计图像识别算法的可行性；为了实现聚类算法，我们需要求解对应的优化问题：最后我们将对算法加以测试和分析。我们遇到的主要挑战是:如何高效计算Wasserstein度量聚类算法和迭代过程中参数的选取我们为克服这些困难，尝试了如下的解决手段:调研并实现了快速计算wasserstein度量的方法，单次迭代复杂度降为线性对于没有显式表达的参数：反复实验，耐心求解；大胆假设，小心求证1.2相关的工作Wasserstein度量的首次提出位于上个世纪3 ，它的一个最知名的等价定义为离散情形下的最优传输问题(Optimal Transport)，因而该问题在几何和图像处理中很早就有应用(例如4 ，5 )。该问题在初期只能使用网络流和单纯形法求解，计算速度较慢，所以对于一阶情形，即EMD(EarthMoversDistance)的分析较多(如6 )。近些年，7 和8 分别提出了两种不同的求解Wasserstein度量的算法，二者都利用了KL散度计算具有线性时间的性质，将难以计算的Wasserstein距离转化为求解一个近似解，而且这两种算法都是具有良好的收敛性质，这为后来该度量的研究打下了坚实的基础。后来也有学者将其推广到更为一般的问题上9 ，总结了各种Wasserstein度量和其变种问题的求解;10 则提出了一个高效，可并行的计算Wasserstein重心问题的算法模型。2中国科学技术大学学士学位论文此后，越来越多的学者关注其在不同领域中的各种应用，11 总结了前人工作，并提出了套用卷积核，使用三角网格上测地曲率等诸多改进，将Wasser-stein度量应用到对应问题和BRDF等诸多经典图形学问题中;1 将颜色直方图和Wasserstein度量相结合，提出了其在颜色空间和形状推断上的相关应用;12 将Wasserstein度量应用到了图像分割问题。此外，十分热门的WassersteinGAN( 13 )中的训练过程也应用到了Wasserstein度量，网络经过训练后可以生成足以以假乱真的人脸图像。1.3论文结构本文主要探究Wasserstein度量在图像识别问题中的一种应用，论文结构安排如下：第一章阐述了我们研究的背景、目的和方法第二章将详细介绍Wasserstein度量，包括其定义和良好的性质第三章提出我们的主要识别算法，将Wasserstein度量和经典聚类算法K-means结合并应用至图像分类，并由此引出的两个需要解决的最优化问题。第四章为主要算法的介绍，将基于正则Wasserstein距离和交替乘子迭代法给出分别给出两种第三章最优化问题的算法，并加以比较。同时，我们将算法应用到第三章的聚类识别框架中，并展示若干实际的例子。第五章是对全文的总结和展望3中国科学技术大学学士学位论文第二章Wasserstein度量2.1 Wasserstein度量的定义对于定义在概率空间上的任意两个概率测度X ,Y 以及概率空间上的一个度量d(x;y),定义这两个概率测度之间的Wasserstein距离为:Wp( ; ) def= inf2 ( ; )(E (d(X;Y)p)1p= inf2 ( ; )( d(x;y)pd (x;y)1p:(2.1)我们可以很容易验证式2.1满足正定性和对称性，而由Minkovski不等式可推出式2.1满足三角不等式，因此这的确是一个定义良好的度量。2.2 Wasserstein度量与最优传输问题之间的关系依据上节定义，我们可以发现Wasserstein度量实际上刻画了两个概率测度，或者实际情况中两个样本分布之间的一个”距离”表征。如果我们把两个分布比作一堆沙子的两个状态，沙子总量都是1个单位，那么该度量大致反映了把这堆沙子从一个状态搬运到另一个状态的所做的”功”。而做功过程中，沙子的几何形状是可以保留的，图2.1解释了这一现象，插值过程中分布的几何性质不会被破坏，分布是真正地”平移”过去的，而传统的欧氏距离是直接“混合”了两个分布。下面我们考虑式2.1的离散情形，在实际生活和应用中，观测到的样本往往将会以经验分布的形式给出，如下式2.2 :=mi=1ai xi; =ni=1bi yi: (2.2)其中fxig和fyig为上的有限离散点集。而以这两个离散分布为边缘分布的联合分布实际上由一个二元函数转化为了一个传输矩阵。(即 ( ; )2Rm n)为方便下文描述，我们记离散情形下的分布为 (x;!),其中x表示位置，!表示对应的位置的频率且!需要满足归一化条件。因此，在离散情形下，Wasserstein度量表现了从若干位置到若干位置的一个最佳对应问题，其表达式与经典的最优传输问题(OT)相同，因而求解离散的Wasserstein度量的传输矩阵，就是求解一个最优运输的优化问题。这是我们下文设计算法的一个重要理论基础。5中国科学技术大学学士学位论文图2.1 Euclidean插值(左)与Wasserstein插值(右)(参考来源:1 )2.3与KL散度的区别和联系传统的信息论和统计理论中，KL散度(参见14 )经常用来刻画两个分布之间的差异。假定我们有同一个样本空间上的两个分布X ,Y ,则它们之间的KL散度定义为:KL( jj ) =x(x)ln (x) (x)dx: (2.3)因为KL散度不对称，所以它不是一个良好的度量，而且它不能刻画两个过于相近的分布，例如:N(0;3) N(;3)计算得KL( jj ) = 12。当!0时，KL散度是发散的。但是相较于Wasser-stein度量，KL散度的一个最大的优势是:离散情形它的计算是线性的。我们在第四章提出的求解Wasserstein度量的相关优化问题的算法，都是利用了KL散度的这一性质简化了算法的计算时间复杂度。6中国科学技术大学学士学位论文第三章基于Wasserstein度量聚类的图像分类识别算法3.1转换图片为待聚类的经验分布我们将主要探究基于Wasserstein度量的图像识别算法的研究。对于图像数据，由于其特征空间一般维数都比较高，所以必须要将特征维数降低，而使用Wasserstein度量计算距离并不需要降维；相反，它可以利用样本分布的几何性质。我们使用最为基础的无监督学习方法聚类算法，对图像进行分类和区分。为了使用聚类算法，我们需要将输入图像转换为经验分布，即直方图。在忽略颜色的情况下，我们直接利用图像的各个像素的灰度值的归一化结果表示一幅图片，把图像看作一个二维空间0;1 0;1上的分布，由此得到了分类算法的主体框架，如算法3.1所示:1读入聚类个数K，输入个数N;2 for i = 1 to N do3读入第i张图像Imgi并对其预处理;4 p Imgi的灰度矩阵;5以p中非0位置构建分布 i ;6 /构建方法，若p(x;y)= 0则Prob(x=width;y=height) = p(x;y)7 end8 labels;centers Clustering (f ig;N;K);算法3.1:图像分类算法框架3.2基于Wasserstein度量的K-means聚类算法K-means是最经典也是最常用的聚类算法之一。一个朴素的思想是把其中关于欧氏距离的计算替换成Wasserstein度量的计算，把重心计算替换为求一个分布作为一些直方图的“重心”，我们就得到了一个基于Wasserstein度量的最简单的聚类算法。详见算法3.2此外，若采用K-means+初始化，我们把其中利用欧式度量计算的部分替换成Wasserstein度量，就会得到如下的一个基于Wasserstein度量的初始化算法，详见算法3.3 :7中国科学技术大学学士学位论文input :一些分布f ig和个数N,聚类个数Koutput:分类标签label,重心C = fcenter 1:center Kg1初始化center 1:center K ;2 while重心的变化量不收敛do3 for i = 1 to N do4 for k = 1 to K do5 dist(i) Wasserstein_Distance ( i;center k);6 end7 mindist;labeli min(dist);8 end9 for k = 1 to K do10 center k Wasserstein_Barycenter (f i j labeli = kg);11 end12 end13 return label,C;算法3.2:基于Wasserstein度量的K-means算法input :一些分布f ig和个数N,聚类个数Koutput:初始化重心C = fcenter 1:center Kg1 center 1 random(1;N) ;2 for k = 2 to K do3 for i = 1 to N do4 for j = 1 to k - 1 do5 dist(k) Wasserstein_Distance ( i;center j);6 end7 mindist(i);labeli min(dist);8 end9 /mindist为每个分布到目前最近“重心”的距离10对mindist归一化;11以mindist为概率，随机选取下标s ;12 center k s;13 end14 return C;算法3.3:基于Wasserstein度量的K-means+聚类初始化算法8中国科学技术大学学士学位论文3.3目标优化问题模型的建立下面我们将给出在算法3.2中需要具体设计的两个函数:Wasserstein_Distance和Wasserstein_Barycenter所求解问题的数学形式。这是我们在接下来将要优化的目标函数。3.3.1 Wasserstein Distance Problem(WDP)我们要解决的第一个最优化问题Wasserstein距离的求解。采用与式2.2的相同记号，我们首先给出一般的离散情形的Wasserstein距离的定义。若我们记a = (a1;:;am)，b = (b1;:;bn)，D = d(x i;y j)pij，则Wasserstein度量的离散形式定义为：Wp( ; ) def= argmin2 ( ; )1p= argmin2 ( ; )tr(D T)1p: (3.1)该优化问题的可行域为 ( ; ) = fT 2R+m n j T1n = a;TT1m = bg:我们将求解该问题在p = 2时的特殊情形。取度量函数d(x;y)为空间上的欧式度量，并对上式平方，我们就得到了待求解的最优化问题(记作WDP):W( ; ) def= argmin2 ( ; )tr(D T): (3.2)3.3.2 Wasserstein Barycentric Problem(WBP)第二个目标问题为Wasserstein重心的求解，即求解分布的加权平均。为此，我们首先定义平权意义下的Wasserstein重心如下:定义3.3.1. 称为一族经验分布f igNi=1的Wasserstein重心，若它是下列最优化问题的解:argmin1NNi=1Wp( ; i)p: (3.3)我们仍然求解p = 2时刻的最优化问题(该问题被记作WBP)。同样地，取d(x;y)为空间上的欧式度量。将式3.2代入后，目标优化问题有如下形式:argmin; i2 ( ; i)Ni=1tr(D( ; i) Ti ): (3.4)9中国科学技术大学学士学位论文第四章WDP和WBP问题的求解和主要结果分析本章我们将解决上一节所提出的两个最优化问题。我们通过调研，结合自身推导，将给出两种不同的算法分别求解WBP和WDP问题，并对结果进行比较和分析。而后，将求解算法带回到上一章中提出的图像识别算法框架中，同时展示部分图像聚类识别的结果。4.1基于正则W距离求解WDP和WBP4.1.1正则W距离与Sinkhorn迭代考虑一个关于离散Wasserstein度量的正则化：Wpp; ( ; ) = argmin2 ( ; )-1 H( ): (4.1)其中H( ) = -i;j i;j ln( i;j)为相对熵。在 充分大时，这个结果可以看作原问题的一个近似解。如果我们记K = e- D，则式3.4的优化问题转化为求解以下最优化问题：W ( ; ) def= argmin2 ( ; )1ijij ln( ijKij)= argmin2 ( ; )KL( jjK): (4.2)这个问题被转化为了关于KL散度的最优化问题。这个问题的对偶问题具有如下表示:argmin2Rm; 2RnTa + Tb -mi=1nj=1e- (Dij- i- j): (4.3)上述式4.2和式4.3的解均由Sinkhorn迭代算法给出(详见Cuturi 15 )。若我们记 和 分别为原问题和对偶问题的最优解，则它有表示形式:=diag(! )Kdiag(! ); (4.4)= - ln(! ) + ln(! )T1nn 1n: (4.5)该迭代算法的解的存在唯一性的证明参考了15 ，下面我们给出该算法的简要流程，如下算法4.1所示：11中国科学技术大学学士学位论文input : D;a;b; output: Wasserstein正则距离dist1 K exp(- D) ;2 u (1;1;:;1);3 while u不收敛do4 v b:=(Ku);5 u a:=(Kv);6 end7 dist trace(D(diag(u)Kdiag(v)T);算法4.1: Sinkhorn迭代经过大量实验，我们发现取定50mean(D) 6 6 200mean(D)的时候，都可以得到较为不错的近似解。最后我们在实验中取定了 = 100mean(D)4.1.2基于次梯度和正则W距离的WBP问题求解在引入正则Wasserstein距离后，式3.4可被改写为:argmin1NNi=1W2 ( ; i) = argmin(x;!); i2 ( ; i)Ni=1KL( ijjexp(- Di): (4.6)考虑!和 i固定时，因为我们取定了p = 2，所以该优化问题为关于x的二次规划，可以直接求得x的迭代公式为:x (1 - )x + N( ki=1xi(T i; )T)diag( 1!): (4.7)算法实现中我们将取定 = 1。下面考虑在x固定的情况下，如何求解!，因为 i的迭代公式由Sinkhorn迭代法给出。如果记f(x;!) = 1NNi=1W2 ( (x;!); i)，则我们有= 1NNi=1i; :为映射!f(x;!)的次梯度(subgradient)，其中 i; 是每个正则W距离问题对偶问题的最优解。函数f(x;!)虽然是凸的，但是非光滑，因此我们需要先对该函数做一次光滑化，之后再进行迭代(理论推导详见Nesterov 16 )。x固定时，以此次梯度作迭代求解!即可。这里对!迭代求解的算法参考了Cuturi7 中的方法，它可以较为快速地求得!的一个近似最优解。综上，我们得到了第一个求解WBP问题的算法，整理如下算法4.2所示:12中国科学技术大学学士学位论文input : f ig;N;lambdaoutput: Wasserstein重心centroid1初始化猜测 (x;!) ;2 while x,!不收敛do3固定x，迭代更新!;4 for i = 1 to N do5 i Sinkhorn (!;! i);6 end7 x 1N(ki=1xi(T i; )T)diag( 1!) ;8 end9 centroid (x;!);算法4.2: WBP问题求解算法一4.2基于Bregman-ADMM算法求解WDP和WBP4.2.1利用BADMM算法求解WDP问题ADMM(Alternative Direction Method of Multipliers)算法是一种将优化问题中变元拆分，然后分别进行迭代的算法。一般地，对于如下的凸优化问题:argmin f(x) + g(y)s:t: Ax + By = c算法按如下方式迭代:记增广函数为L (x;y;z) = f(x) + g(y) + zT(Ax + By - c) + 2Ax + By - c22: (4.8)则有迭代公式:xk+1 = argminxL (x;yk;zk);yk+1 = argminyL (xk+1;y;zk);zk+1 = zk + (Axk+1 + Byk+1 - c):(4.9)而BADMM(Bregman ADMM)算法将式4.8中增光Lagrange函数L 中二次罚函数项替换成了Bregman散度的形式，迭代形式与ADMM算法相同。Bregman散度17 的数学表达式如下:BDf(x;y) = f(x) - f(y)- : (4.10)若我们取式4.10中凸可微函数为f(x) = xln(x)，则散度表达式与KL相同，下面的算法中我们将利用这一性质改写目标优化函数。为了在式3.2中利用BADMM算法，我们考虑与之等价的一个最优化问题:13中国科学技术大学学士学位论文argminA2A; B2Btr(D( ; ) TA);s:t: A = B:(4.11)其中A = fT 2R+m n j TT1m = bg ; B = fT 2R+m n j T1n = ag:这符合使用ADMM算法的条件，因此我们可以直接写出基于BADMM算法得到的求解迭代形式，如下所示:L ( A; B; ) = tr(D( ; ) TA) + tr( ( A - B)T)+ (KL( A; B) + KL( B; A);(4.12)n+1A = argminA2Atr(D( ; ) TA) + tr( n TA) + KL( A; nB); (4.13)n+1B = argminB2B- tr( n TB) + KL( k; n+1A ); (4.14)n+1 = n + ( n+1A - n+1B ): (4.15)4.2.2利用BADMM算法求解WBP问题为了本节算法叙述所需，在这里提出一些符号约定。假定所有观测样本是位于d维样本空间上的，样本个数为N，第k个样本中有mk个位置观测频率非0。记n =Ni=1mk假定 k = (xk;!k)为一个样本分布，其中xk = (xk1;:;xkmk)2Rd mk为位置，!k = (!k1;:;!kmk)为观测到的频率假定 = (x;!)为待求分布，m为概率密度非0位置的个数。(即x2Rd m)假定Di = D( ; i)为度量矩阵，即Di(i;j) = d( (xi); (xj)p假定 k = ( ki;j)2Rm mk+为 和 k的传输矩阵，并且记 = ( 1; 2;:; k)利用上述符号，那么式3.2的目标优化问题可被表述为:argmin(x;!); Nk=1tr(Dk( k)T);k2 k(!) = fA2Rm mk j Ai;j 0 A1mk = ! 1TmA = !kg:(4.16)首先，和4.1节方法相同，在p = 2且取d(x;y)为欧氏距离的前提下，若固定!; ，则求解x是一个无约束的关于x的二次规划的问题，它的迭代公式为:记X = (x1;:;xn)2Rd n ,则:14中国科学技术大学学士学位论文x = 1NX Tdiag(1!): (4.17)在固定x的前提下，同WDP问题求解类似，我们考虑与式4.16的目标相同的一个等价优化问题：argmin!; kA2k;A; kB2k;B(!)Nk=1tr(Dk( kA)T);s:t: kA = kB:(4.18)其中k;A = f ki;j 0jmi=1ki;j = !kj;j = 1;:;mkg;k;B(!) = f ki;j 0jmkj=1ki;j = !i;i = 1;:;mg:此问题若!固定，即可采用BADMM算法迭代求解，但是只固定x求解!和 是一个很困难的事情，因为这样该问题是一个非凸的优化问题(由于!是作用在 kB的可行域上，这个问题对于变量!是非凸的)，很难处理。为了求解该优化问题，参考10 ，首先先写出固定!情形下的BADMM迭代公式:n+1A = argminkA2k;Aki=1(tr(Dk( kA)T) + tr( k;n( kA)T) + KL( kA; k;nB ); (4.19)n+1B = argmin!; kB2k;B(!)ki=1(-tr( k;n( kB)T) + KL( kB; k;n+1A ); (4.20)n+1 = n + ( n+1A - n+1B ): (4.21)观察发现，求解(4.19 )式的问题只需要按列使用Lagranian乘数法求导计算即可，问题在于(4.20 )式，如果把!看作固定常量，(4.20 )式的求解方法和(4.19 )式类似，每个变量按行求解Lagranian乘数法，我们会得到如下迭代公式： k;n+1A = k;nB exp(-Dk + k); (4.22)k;n+1A (i;j) = k;n+1A (i;j)mi=1 k;n+1A (i;j)!kj; (4.23) k;n+1B = k;n+1A exp(k); (4.24)k;n+1B (i;j) = k;n+1B (i;j)mkj=1 k;n+1B (i;j)!i: (4.25)15中国科学技术大学学士学位论文其中为Hadammard乘积，表示矩阵对应位置元素相乘。如果对于一个固定的!做迭代，则该迭代过程和BADMM算法是完全相同的，其结果也是收敛的。下面考虑迭代!，该问题由于!的存在成为了一个非凸的优化问题，较难处理，我们使用基于KL散度的Gradient Decent Projection 9 的迭代方法求解，即求解以下优化问题:argmin!Ni=1KL(!jj k;n+1B 1m):该优化问题的解可以直接求得:!i = Nk=1(mkj=1 k;n+1B (i;j)1Nmi=1Nk=1(mkj=1 k;n+1B (i;j)1N: (4.26)为了保证迭代的稳定性，我们采用一个近似求解的手段:交换非对称的KL散度的两个变量。这样做虽然不能严格得到其收敛性的证明，但是我们将会得到形式上更简单，实际情况下更为稳定，效果尚可的一个近似公式。交换变量后，新的优化问题转变为：argmin!Ni=1KL( k;n+1B 1mjj!): (4.27)求解该问题，我们最终得到关于!的一个近似迭代公式为:!i = B1n1Tm B1n: (4.28)下面我们给出该算法的具体流程，如下算法4.3所示:在初始化该问题时，我们既可以根据给定的样本集生成一个服从多维正态分布的随机变量，对其采样作为初始值，也可以根据实际问题自行生成初始分布(比如在矩形上均匀撒点)。此外，罚函数因子是一个重要参数。由于算法代数敏感度较高，罚因子过大或过小都会对结果造成剧烈影响。经过实验，取= 2NNk=1mean(D( ; k):较佳。为了简化计算量，算法4.3中经实验可取 = 10，表示每迭代10次更新一次w4.3两种算法的结果比较和分析4.3.1 WDP问题求解结果我们取定8张图，在相同的迭代收敛条件下(残差6 0:0001或2000步停止)，计算得到的结果如表4.1和表4.2所示。在实际测试过程中，基于正则Wasserstein16中国科学技术大学学士学位论文input :一些分布f ig和个数N, output: Wasserstein重心centroid1初始化猜测 (x;!) ;2初始化 kB;3 while A - Bfro不收敛do4 for k = 1 to N do5按照式4.23迭代更新 kA;6按照式4.24求解 kB;7 end8按照式4.28迭代更新!;9 for k = 1 to N do10按照式4.25迭代更新 kB;11 end12每 次循环更新x;13 end14 centroid (x;!);算法4.3: BADMM求解WBP问题距离的Sinkhorn迭代方法迭代收敛速度要远远优于BADMM，运算速度也很快，表4.1的生成只用了0.3秒，而BADMM算法用了约70秒。不过，由于正则Wasserstein距离是近似解，所以表4.1的每个元素，包括对角元均有一定误差。BADMM算法迭代求解WDP问题是收敛的，但是对于偏差过大的图像(如汉字“我”的手写体和其他图像)，罚函数因子会变大，导致迭代收敛速度变慢，所以出现了在2000次迭代后矩阵仍不对称的现象。表4.1基于正则W距离的WDP结果0.0006 0.0073 0.0017 0.0100 0.0288 0.0328 1.4202 0.01240.0073 0.0005 0.0099 0.0130 0.0268 0.0336 1.4244 0.01180.0017 0.0099 0.0009 0.0114 0.0308 0.0351 1.3920 0.01410.0100 0.0130 0.0114 0.0011 0.0284 0.0358 1.5089 0.02040.0288 0.0268 0.0308 0.0284 0.0021 0.0102 1.3205 0.01230.0328 0.0336 0.0351 0.0358 0.0