EM梯度算法

EM梯度算法（精选三篇）

EM梯度算法 篇1

关键词：EM算法,缺失数据,参数估计

1、引言

在系统参数估计问题中极大似然估计的理论相对简单, 当模型中所有的变量都为可观测的时候, 该方法可以直接处理, 但模型中有隐含的未观测变量或者观测的数据不完整时候, 似然估计问题就变很复杂。为了解决不完全数据参数估计问题, A.P.Dempster, N.M.Laird和D.B.Rudin于1977年[1]提出EM (Expectation-Maximization) 算法, 并证明了该算法的收敛性, 但是当缺失的信息量很大时算法的收敛速度将很缓慢, 随后产生了各种改进的算法[6], 这些改进的算法都着眼于改善收敛速度和保证收敛性质。

本文讨论的是一种改进的A-ECM算法来提高EM算法的收敛速度, 该算法结合Aitken加速和ECM算法各自的优点, 既实现了算法加速也达到了稳定收敛的目的。

2、EM算法概述

EM算法是不完全数据问题中用于极大似然估计很常用的一种迭代算法, EM算法公式化了数据缺失机制这个概念。

2.1缺失数据极大似然估计

完全数据极大似然估计方法:设总体为X, 其分布律P{X=x}=f (x;θ) 形式为已知, θ为参数。设X1, X2, …, Xn是来自总体X的一个样本, 其对应的联合分布律为则定义如下函数:

L (θ) 称为样本的似然函数。极大似然估计就是求取θ使得L (θ) 在固定观测x1, x2, …, xn下取得最大值, 即

当导致缺失数据机制可以忽略时候, 根据这个似然Lign (θ|Y obs) 可作出关于θ的估计。

2.2EM算法

EM算法是不完全数据问题中计算极大估计的一种迭代算法。EM算法的每一步迭代中包括一个E步:期望步 (Expectation Step) 和一个M步:极大似然步 (Maximum Likelihood Step) 。E步在给定已观测到的数据和现有参数下, 求缺失数据的条件期望, 然后以这些期望值取代缺失数据。设θ (t) 是当前参数θ的估计, EM的E步在θ是θ (t) 下, 求完全数据似然的期望

M步比较简单, 即假设没有数据缺失的情况下进行θ的ML估计, EM的M步与一般似然估计中求l (θ|Y) 的θ估计一样。按极大化这个期望的完全数据似然确定θ (t+1) , 即

2.3EM算法的收敛性质

Dempster、Laird和Rubin[1977]证明了对于任意初始值EM算法都能收敛到局部极大。如果涉及的问题存在唯一的极大似然解, 那么EM算法最终迭代给出的结果就是极大似然的估计值。EM算法虽然具有稳定收敛的性质, 但收敛的速度缓慢。当缺失数据比例较大时候, 它的收敛速率也极其缓慢。

3、A-ECM算法与实现

3.1%Aitken算法

假设θ*为EM估计算法中θ的收敛值, 即

并设θk→θk+1的映射为θk+1=g (θk) , J=J (θ) 为g (θ) 在θ处的梯度。由θ在θ*附近的泰勒展开, 我们有如下近似公式

于是可以得到

其中θk是Aitken第k次迭代的估计值, θEMk+1是以θk为初始值通过EM算法得到的第k+1次的估计值。

Aitken算法的第k+1是按照如下步骤进行的:

1、θk是Aitken第k次迭代的估计值

2、以θk为标准EM算法的初始值, 进行迭代近计算第k+1步从估计值θEMk+1

3、然后将θEMk+1代入公式 (9) 进行计算得到第k+1步的Aitken估计值θk+1

4、重复以上的步骤直至收敛为止。

由于J (θ*) =I-i-1comiobs, 其中icom为以上已经定义过的完全数据信息量, iobs为观测数据信息量。Louis提出用i-1comiobs代替J (θ*) 进行计算, 得到

Jamshidian和Jenrish (1993) [3]指出实际上Louis采用i-1comiobs代替J (θ*) 的做法等价与数值计算中的Newton-Raphson.法。

1.1ECM算法

引理1[2]:EM算法在每次迭代后均能增加似然函数值。

引理2[2]:如果E (θ′, θ) 关于其两个变量都是连续的, 则EM能得到一个稳定的收敛值。

在A-ECM算法中E步采用和EM相同的计算过程, M步作如下的处理:把参数θ划分为θ=θ1, …, θs, gs (θ) 是关于θ=θ1, …, θs的函数。在M步的第k+1此迭代中, 要进行s次估计, 每次估计出θ中的一个参数。当估计θ1时, 假定其他s-1个估计的参数不变, 然后采用一般的似然估计的方法估计它, 接着用新的估计值代替原来的估计值进行对θ2参数的估计, 并保持其他参数不变。以此类推完成对θs的估计, 即估计第k+1次迭代中第i+1个参数时候

θi+1k+1=i进行对第i+2个参数进行估计, 直到完成θsk+1的估计, 然后对θk+2进行估计。

3.3Aitken-ECM算法

ECM算法:在迭代的前期能有很好的收敛速度, 其后期速度仍然不是很理想, 但ECM会在保证EM收敛的基本条件下收敛到一个稳定点[2]。Aitken法:当参数离ML较远的时候Aitken法显得不是很成功, 而末尾使用Aitken法, 收敛的速度很快。Aitken法基于数值近似的计算其缺点在于, 有时候会在收敛值附近振荡, 而且Aitken法迭代过程中不保证似然是稳定增长的。[3]

为能实现整个过程都有较好的收敛速度, 而且要稳定的收敛到一个极大值, ECM和Aitken法都不能做到, 但是从他们的特性来看他们是很好的互补, 以下提出一个混合ECM和Aitken法, 他们是交替计算的, 记为Aitken-ECM (A-ECM) 方法。首先在参数估计的初始阶段采用ECM算法, 当算法的收敛速度比较缓慢的时候我们切入Aitken算法来提速。考虑到Aitken在θ*附近无法稳定收敛, 则再转入ECM算法进行最后的稳定收敛过程。这里要注意的两个问题:

1、ECM算法收敛速度进入缓慢阶段, 切换为Aitken算法的判定可以根据以下准则:

2、可以根据似然函数L (θk) 的数值来判定Aitken算法是否出现小范围的振荡, L (θk) 出现振荡而非稳定增加的时候就可切换入ECM算法完成最后的收敛过程。

4实例分析

4.1系统仿真

以下通过一个实例从估计误差和收敛速度两方面来分别比较A-ECM和EM算法进行参数估计的效果, 这里采用的仿真系统如下, 其中H为单位阵列, N为计算步长。

系统在噪声影响下输出YN=[Y1, Y2, Y3, Y4, Y5], 如图1系统输出观测

将输出观测YN=[Y1, Y2, Y3, Y4, Y5]分别代入A-ECM和EM算法。由于我们更关心的是系统状态矩阵的参数估计, 这里就列出两种算法状态矩阵F估计的结果, 其他估计值在以下估计误差分析中列出他们的估计误差。

4.2估计误差分析

记估计误差为, 其中为EM的估计值, 其中用矩阵范数来度量估计值和原始值之间的误差。设矩阵M=mij, 则范数定义为

由以上估计误差表可以看出A-ECM算法的估计效果明显好于EM, 在300步的迭代过程中A-ECM算法已经很接近稳定值了, 而要想让EM能有这样的效果至少要2000步以上。

4.3算法收敛速度分析

有关A-ECM和EM算法的收敛情况可以通过下图看出

关于图2 A-ECM与EM算法似然函数收敛的说明:

1.从图中可见A-ECM算法在大概15步以内经过ECM算法很快提速到收敛值附近, 然后在切入Aitken进行加速, 大约在100步左右再切换回ECM进行最后的收敛过程, 而EM算法整个过程都是平稳缓慢的收敛。

2.可以看出ECM很快就对迭代的初始数值进行了修正, 所以在算法开始的阶段ECM就比EM有较好的初始值, EM在达到ECM的初始值时已经运算了15步左右, 而此刻ECM已经在收敛值附近了。

关于图3 A-ECM与EM算法收敛速度的说明:该图采用半对数有利于显示细微的变化过程。在开始阶段A-ECM的收敛速度大概是EM的10倍, 在20到100步之间采用Aitken算法, 使得速度有了提升并维持了100步左右。

5、小结

通过以上实例的仿真分析可以看出:在算法加速方面, Aitken算法弥补了ECM算法中期收敛慢的不足, 完成了加速;在算法收敛方面, ECM具有稳定单调收敛的性质, 弥补了Aitken数值算法在后期难于收敛的缺陷。A-ECM算法结合ECM收敛性和Aitken加速的特点, 实现了对EM算法的加速。

参考文献

[1].A.P.Dempster, N.M.Laird, D.B.Rudin.Maximum likelihood from in-complete data via the EM algorithm.J.Roy.Statist.Soc.1977, B, 39, 1-38.

[2].X.L.Meng, , D.B.Rudin.Maximum likelihood estimation via the ECM algorithm:a general framework.Biometrika80, 267-278, 1993.

[3].M.Jamshidian, S.R.Jennrich, M.H.Chen.Conjugate gradient accelera-tion of the EM algorithm.J.Am.Statist.Assoc.1993, 88, 221-228.

[4].C.H.Lin, D.B.Rudin.The ECME algorithm:a simple extension of EM and ECM with fast monotone convergence.Biometrika, 1994, 81, 633-648.

[5].K.Lange.A gradient algorithm locally equivalent to the EM algorithm.J.Roy.Statist.Soc.1995, B57, 425-437.

斯坦福大学机器学习梯度算法总结 篇2

1基础概念和记号

线性代数对于线性方程组可以提供一种简便的表达和操作方式，例如对于如下的方程组：

4x1-5x2=13

-2x1+3x2=-9

可以简单的表示成下面的方式：

X也是一个矩阵，为(x1,x2)T，当然你可以看成一个列向量。

1.1基本记号

用A ∈表示一个矩阵A，有m行，n列，并且每一个矩阵元素都是实数。用x ∈ , 表示一个n维向量.通常是一个列向量.如果要表示一个行向量的话，通常是以列向量的转置（后面加T）来表示。

1.2向量的内积和外积

根据课内的定义，如果形式如xT y，或者yT x，则表示为内积，结果为一个实数，表示的是：，如果形式为xyT，则表示的为外积:。

1.3矩阵-向量的乘法

给定一个矩阵A ∈ Rm×n，以及一个向量x ∈ Rn，他们乘积为一个向量y = Ax ∈ Rm。也即如下的表示：

如果A为行表示的矩阵（即表示为），则y的表示为：

相对的，如果A为列表示的矩阵，则y的表示为：

即：y看成A的列的线性组合，每一列都乘以一个系数并相加，系数由x得到。

同理，yT=xT*A表示为：

yT是A的行的线性组合，每一行都乘以一个系数并相加，系数由x得到。

1.4矩阵-矩阵的乘法

同样有两种表示方式:

第一种：A表示为行，B表示为列

第二种，A表示为列，B表示为行：

本质上是一样的，只是表示方式不同罢了。

1.5矩阵的梯度运算（这是老师自定义的）

定义函数f，是从m x n矩阵到实数的一个映射，那么对于f在A上的梯度的定义如下：

这里我的理解是，f（A）=关于A中的元素的表达式，是一个实数，然后所谓的对于A的梯度即是和A同样规模的矩阵，矩阵中的每一个元素就是f(A)针对原来的元素的求导。

1.6其他概念

因为篇幅原因，所以不在这里继续赘述，其他需要的概念还有单位矩阵、对角线矩阵、矩阵转置、对称矩阵（AT=A）、反对称矩阵（A=-AT）、矩阵的迹、向量的模、线性无关、矩阵的秩、满秩矩阵、矩阵的逆（当且仅当矩阵满秩时可逆）、正交矩阵、矩阵的列空间(值域)、行列式、特征向量与特征值……

2用到的公式

在课程中用到了许多公式，罗列一下。嗯，部分公式的证明很简单，部分难的证明我也不会，也懒得去细想了，毕竟感觉上数学对于我来说更像是工具吧。

转置相关：

•(AT)T = A

•(AB)T = BT AT

•(A + B)T = AT + BT 迹相关： • For A ∈ Rn×n, trA = trAT.• For A, B ∈ Rn×n, tr(A + B)=trA + trB.• For A ∈ Rn×n, t ∈ R, tr(tA)= t trA.• For A, B such that AB issquare, trAB = trBA.• For A, B, C such that ABC issquare, trABC = trBCA = trCAB。当乘法变多时也一样，就是每次从末尾取一个矩阵放到前面去，这样的矩阵乘法所得矩阵的迹是一致的。秩相关

• For A ∈ Rm×n,rank(A)≤ min(m, n).If rank(A)= min(m, n), 则A称为满秩

• For A ∈ Rm×n,rank(A)= rank(AT).• For A ∈ Rm×n, B ∈ Rn×p,rank(AB)≤ min(rank(A), rank(B)).• For A, B ∈ Rm×n,rank(A + B)≤ rank(A)+rank(B).逆相关：

•(A−1)−1 = A

• If Ax = b, 左右都乘以A−1 得到 x = A−1b.•(AB)−1 = B−1A−1

•(A−1)T =(AT)−1.F通常表示为A−T.行列式相关：

• For A ∈ Rn×n, |A| = |AT |.• For A, B ∈ Rn×n, |AB| = |A||B|.• For A ∈ Rn×n, |A| = 0，表示矩阵A是奇异矩阵，不可逆矩阵

• For A ∈ Rn×n and A 可逆, |A|−1 = 1/|A|.梯度相关： • ∇x(f(x)+ g(x))= ∇xf(x)+ ∇xg(x).• For t ∈ R, ∇x(t f(x))= t∇xf(x).• ∇xbT x = b

• ∇xxT Ax = 2Ax(if A 对称)

• ∇2xxT Ax = 2A(if A 对称)

• ∇A|A| =(adj(A))T = |A|A−T.adj=adjoint

3梯度下降算法和正规方程组实例应用

例子用的是上节课的房价的例子，有一组数据，有房子面积和房子价格，输入格式举例：

老师定义的变量如下：

m:训练样本的数目

x：输入的变量（输入的特征，在这个例子中为房子面积，后来又加了一个房子的卧室数目）

y :输出变量（目标变量，这个例子中就是房价）

(x,y)：表示的是一个样本

：表示的第i个样本，表示为。

3.1监督学习概念

所谓的监督学习即为告诉算法每个样本的正确答案，学习后的算法对新的输入也能输入正确的答案。监督指的是在训练样本答案的监督下，h即为监督学习函数。

此例中我们假设输出目标变量是输入变量的线性组合，也就是说，我们的假设是存下如下的h（x）：

Theta表示是特征前面的参数（也称作特征权重）。也就是经过h(x)之后得到的就是预测的结果了。

如果假设x0=1，那么原来的h(x)就可以简单的表示为如下形式：，其中n为特征数目，我们为了表达简便，把theta和x都写成向量的形式。下面就是如何求出θ（向量）使得h(x)尽可能接近实际结果的,至少在训练集内接近训练集中的正确答案。

我们定义一个花费函数(costfunction)，针对每一组θ，计算出h(x)与实际值的差值。定义如下：

这也是用的最小二乘法的思想，但是之所以乘以1/2是为了简化后面的计算。针对训练集中的每一组数据。剩下的问题就是求得minJ（θ）时的θ取值，因为J(θ)是随着θ变化而变化，所以我们要求得minJ（θ）时的θ就是我们想要的θ（这个min也叫做最小花费函数)，怎么样求出这组theta呢？采用的方法就是梯度下降算法和正规方程组。我们首先来看梯度下降算法。

3.2梯度下降算法

梯度下降算法是一种搜索算法，基本思想可以这样理解：我们从山上的某一点出发，找一个最陡的坡走一步（也就是找梯度方向），到达一个点之后，再找最陡的坡，再走一步，直到我们不断的这么走，走到最“低”点（最小花费函数收敛点）。

如上图所示，x,y表示的是theta0和theta1，z方向表示的是花费函数，很明显出发点不同，最后到达的收敛点可能不一样。当然如果是碗状的，那么收敛点就应该是一样的。

算法的theta更新表示如下：

对每一个theta(j),都先求J(θ)对theta(j)的偏导(梯度方向)，然后减少α，然后将现在的theta(j)带入，求得新的theta(j)进行更新。其中α为步长，你可以理解为我们下山时走的步子的大小。步子太小了，收敛速度慢，步子太大了，可能会在收敛点附近来回摆动导致无法到达最低点。P.S.这个符号根据老师所说理解为程序中的赋值符号(=号)，如果是=号，则理解为值是相等的(编程里面的==号)。下面我们先理解下，假设现在训练集只有一组数据求关于theta(j)的偏导：

带入

可以得到关于一组数据的theta(j)的表达式，不妨，这组数据就是第i组，则表示为：

然后我们将这个更新theta(j)的方法扩充到m个训练样本中，就可以得到下面的式子：

P.S.最外面的那个xj(i)的理解为：第i组数据中的第j个特征(feature)值。

3.2.1批量梯度下降算法（batch gxxxxdxxxx algorithm）

重复执行上面的这个更新步骤，直到收敛，就可以得到这组θ的值了。就是这个过程：。

这个算法就是批量梯度下降算法，为什么叫批量梯度下降？因为注意到上式中每更新一个θj都需要计算所有的样本取值，所以当样本数目非常大的时候(例如上万条甚至数十万条的时候)，这样的更新非常慢，找θ也非常慢，所以就有了另外一种改进的梯度下降算法。

3.2.2随机梯度下降算法/增量梯度下降算法

做一个小小的改进，用一个样本做一个theta的更新，比如用样本1做theta(1)的更新，用样本2做theta(2)的更新，以此类推。这种方法的好处是速度上肯定比批量梯度下降算法快，而且样本数据越多，体现应该就越明显。劣势是得到的收敛点的值和批量梯度算法比起来也许不是最优的值。

3.2.3梯度下降算法总结

不管是批量梯度算法还是随机梯度下降算法，他们的共同点有以下：

1.时间复杂度都是O(mn)(m为样本数目，n为特征值/影响因子数目)

2.都有梯度下降性质：接近收敛时，每次“步子”（指实际减去的数，而不是前面定义的α，α是手动设置参数，人为改变才会变）会越来越小。其原因是每次减去α乘以梯度，但是随着收敛的进行，梯度会越来越小，所以减去的值会。

3.判定收敛的方法都是如下两种:

1)两次迭代值改变量极小极小

2)J(θ)的值改变量极小极小

3.3正规方程组

写在前面：这种方法是另一种方法了，和梯度下降算法就没啥联系了！！！

首先回顾下前面定义的矩阵梯度运算：

例如:

则：

定义一个矩阵，称作设计矩阵，表示的是所有的样本的输入：

因为前面得到的结论：

（θT*x(i)和x(i)的转置*θ结果是一样），可以得到

可以写成如下的形式：

又因为对于任意向量，所以可以得到：

运用下面介绍的一系列性质：

（5）是由(2)和(3)得到的，进行下面的推导

中间加tr不变的原因是因为是一个实数(看成1x1矩阵)加迹等于他本身。

将上式设为0，得到正规方程组

EM梯度算法 篇3

关键词：多个体系统；分布式优化；随机无梯度；通信时延

中图分类号：TP13 文献标志码：A文章编号：1672-1098（2016）01-0034-06

Abstract：Considering the delay of information communication among agents， which affect the convergence speed of the algorithm， the randomized gradient-free method for multi-agent optimization with communication delay was proposed， where its assumed that every agent only knows its own local objective function. The optimization goal is to minimize a sum of local objective functions through the interaction of delay information among agents in the system. Firstly， the optimization problem with delay was converted into the optimization problem without delay through augmenting delay nodes. Because the local objective function of agent is likely to be nonconvex， its subgradient does not exist or it is hard to be calculated， the distributed randomized gradient-free method was used. The theoretical analysis showed that the proposed algorithm is still convergent if the communication delays are upper bounded.

Key words：multi-agent system； distributed optimization； randomized gradient-free， communication delay

近年来，多个体分布式凸优化问题及其优化算法引起了人们的广泛关注，而多个体分布式优化算法是在集总式算法的基础上发展起来的。所谓集总式算法是指在多个体系统中，不是所有个体都发挥同样的作用，只有某个个体处于中心地位，负责处理其他个体的数据，并将数据反馈给其他个体。和集总式算法不同，分布式算法则是指多个体系统中的每个个体都对应着一个局部凸目标函数，并且个体之间进行信息交流，最终求得凸目标函数的最小值。和以往的集总式算法相比，分布式算法有很多优点，尤其在许多大规模的优化问题中占有很大优势，并且在生物工程、人工智能等许多领域有广泛应用，因此研究多个体的分布式优化算法有很大的意义。

随着计算机的广泛应用，人们进入了大数据云计算的时代，因此对多个体分布式优化的研究也越来越深入。但这些方法主要是标准次梯度和一致性算法的结合。标准次梯度算法是将总的最优化任务分解，同时每个个体需要将自身的信息与周围邻居个体的信息进行加权组合，再根据自身的次梯度信息进行最优化，经过一系列的迭代运算，使得所有个体的状态都达到一致。事实上，一致性算法也是广泛研究的一个课题，即个体间通过信息交流使所有个体的状态最终达到一致并使结果达到最优。文献[1]最早给出了标准次梯度方法并分析了其收敛性。文献[2]922介绍了约束一致性和优化算法。在此基础上，文献[3]则介绍了基于随机投影的次梯度算法，文献[4]1715给出一种基于一致性算法的原始对偶次梯度方法。在文献[4]1720的启发下，文献[5-6]提出了一种分布式对偶平均算法（DDA）以及基于Push-sum的DDA算法。上述研究都是适用于多个体系统中的每个个体对应的局部目标函数存在凸函数且次梯度的情况。而文献[7-8]研究的是个体的局部目标函数是非凸的、次梯度不存在或很难计算的情况，因此文献[9]提出了一种分布式随机无梯度优化算法。

然而，以上都是假设在任何时刻每个个体都可以即时的和周围个体之间进行信息交流。而在文献[10-12]1108中，由于有限的带宽，随机的传输延迟以及不确定的连接拓扑，个体之间的信息交流可能出现延时的情况。为此，文献[12]1139给出了时延情形下的分布式次梯度优化算法以及具有通信时延的二阶系统一致性。

本文主要研究时延情形下的分布式随机无梯度优化算法。

1符号说明

网络中个体间的信息通信通常可以建模成一个有向图G（k）=（V，E（k），P（k）），其中V=（1，2，…，n）表示个体集合，n表示个体数目，E=（e1，e2，…，en）表示边集，P（k）表示网络拓扑图G在k时刻的权重邻接矩阵，并且G是一个有向图[13]。令Rn为一个n维向量空间，‖x‖为向量x的欧几里得范数，ΠX[x]表示向量x到集合X上的欧几里得投影，[x]i表示向量x的第i个分量，xT表示向量x的转置，[P]ij代表矩阵P的第i行j列的元素，E[x]代表向量x的期望值。f（x）则表示函数f（x）在x处的梯度。

2问题描述

考虑如下由n个个体组成的多个体系统分布式凸优化问题：

minx∈Rmh（x）=∑ni=1hi（x） （1）

这里h（x）为目标函数，x∈Rm是一个全局决策向量，hi∶Rm→R是个体i（i∈V）的自身目标函数并仅为个体i知道，同时假设其是Lipschitz连续的且Lipschitz常数为G0（hi），XRm是非空闭凸集。式（1）表明整个系统的目标函数可以分解成系统中各个个体自身目标函数之和。

考虑每个hi不一定是光滑的且可能是非凸函数。因此其次梯度不存在或很难计算。故已有相关基于次梯度的分布式优化算法对本文不再适用。为此，用一个高斯函数hiμi 来近似代替hi，此时式（1）的光滑版本如下：

minx∈Xhu（x）=∑ni=1hiui （x） （2）

这里hiui （x）是hi（x）的高斯近似[7]1116且满足：

hiui （x）=1ω∫Rmhi（x+μiξ）e-12‖ξ‖2dξ （3）

这里ω=∫Rme-12‖ξ‖2dξ=（2π）m2，μi代表函数fiμi 的光滑系数，是一个非负标量，ξ是随机向量。

3算法介绍

31DRGF法

在实际情况中，由于多个体在信息交流的过程中存在通信时延，为此提出如下的时变时延情形下的分布式随机无梯度优化算法（DRGF法）。

1） 初始化

① 选择一个随机向量xi（0）∈X，x∈X；

② 对每一个i，用高斯分布的方式生成一个局部随机序列{ξki}k≥0。

2） 迭代

① 计算平均权重θi（k+1）=∑nj=1[Φ（k）ij]xj（k-cji（k））=∑n（q+1）j=1[Φ（k）ij]xj；

② 计算xi（k+1）=ΠX[θi（k+1）-ηkgui（xi（k））]，这里ΠX[x]代表向量x到集合X上的投影，通过投影将约束集中的个体与其他个体分开，只讨论在约束集中的个体。

这里的ηk>0是迭代步长，且步长满足：

ηk>0，∑∞k=0ηk=∞，∑∞k=0η2k<∞。

随机无梯度的预测值可表示为：

gui（xi（k））=

hi（xi（k）+uiξi（k））-hi（xi（k））uiξi（k）

32系统扩维

为了方便研究时延情形下的分布式随机无梯度优化算法的收敛性，下面将介绍如何利用系统扩维把有时延的优化问题转化为无时延的优化问题。在网络图G（k）中增加一些虚拟节点（见图1），这些节点的作用只是为了消耗信息在系统传递中的时延，自身没有自环。和增加的虚拟节点相比，网络中真实存在的点不仅可以传递信息，自己也可以处理信息。因此可以看出，虚拟节点只起到信息传递的作用，而不参与算法的迭代。

由于每个时刻通信网络每条边上的时延cij（k）都是不同的，所以令q=max ij。原来系统中的n个节点增加nq个节点即扩维后的多个体系统，共有n+nq个节点。为了不影响优化式（1），假设新增的时延个体所对应的局部目标函数值均为零。系统扩维后的权重矩阵为Φ（k）=[Φ1（k），Φ2（k）]T，其中[Φ1（k）]i，nl+j=（Φ（k））ij，l=cij（k）

0，otherwise，Φ2（k）=[Inq，0nq×n]，对于任意的l=1，2…q，i，j∈I，可以看出Φ（k）仍为行随机矩阵，但其主对角元素已不再是全部为正，因此无时延情形下的文献[7]，[8]，[13]的方法对本文都不再适用。

4收敛性分析

首先，做出如下假设：

假设1对于有向图G以及任意有向边上的时变时延cij（k），有：

1）有向图G是强连通图，其所对应的邻接矩阵Φ（k）是行随机的，而不一定是双随机的；

2）假设0≤cij（k）≤ij（k）<∞，即时变时延有界，令q=max（vi，vj）∈Eij

为了方便分析DRGF算法的收敛性，给出下列引理：

引理1[7]1 122对每一个i∈V有以下结论成立：

1） hiui （x）是可微凸集并且满足：

hi（x）≤hiui （x）≤muiG0（hi）

2） 梯度hiui （xi（k））满足：

E[giuixi（k）|Hk]=hiui （xi（k））

3）随机无梯度的预测值gui（xi（k））满足：

E[giui xi（k）|Hk]≤（m+4）2G20 （hi）

引理2[2]777令Ψ（k·s）=[Φ（k）Φ（k-1），…Φ（s）]，则可得以下结论：

1） 当k→+∞时，Ψ（k，s）的每一列[Ψ（k，s）]j收敛到一个正向量j（s）1，即：limk→+∞（[Ψ（k，s）]j-j（s）1）=0对所有j∈{1，2，…，n（q+1）}都成立。这里j（s）>0且∑n（q+1）j=1j（s）=1。

2） 存在一个正整数ν和一个非负数0≤α<1，使得：‖[Ψ（k，s）]ij-j（s）‖≤（q+1）4n4αk-s-3M-MvMv‖Ψ（s，s）‖j对所有i，j∈{1，…，q}和所有k≥s成立。

定理1在假设1，引理1成立的情况下，用上述DRGF法生成一个序列{xi（k）}k≥0，这里假设每一个hi都是Lipschitz连续的， 则对每一个t≥1和 j∈V，有：

定理1给出了目标函数的期望值E和最优值h（x*）之间的误差，不难看出，该误差和个体数n，系统扩维后增加的节点数nq，迭代步长η和Lipschitz常数有关。由定理1看出只要通信时延cij有上界，算法依然收敛。当cij=0时，则转化为无时延的情况。此外，还可以看出，时延情形下的迭代误差比无时延的迭代误差大，这表明若个体在信息交流的过程中存在通信时延，则会对个体间协同解决优化问题造成一定的影响。

推论1当假设1成立时，用上述方法生成一个序列{xi（k）}k≥0，步长为{ηk}k≥0且limk→∞ηk=。则有：

由推论1可以看出个体i在k时刻的状态xi（k）收敛到k时刻所有个体状态的平均值k。即此时多个体状态达到平均一致。此外多个体优化问题此时也取得最优值。

定理2在上述引理1，引理2以及推论1成立的前提下，又假设：若limk→∞ηk=且=0，则∑∞k=0ηk=∞。则有如下结论成立：

lim supt→∞E[f（j（t））]≤m 0∑n（q+1）i=1μi+D0 （16）

其中D=（n（q+1））（m+5）2（92+2（n（（q+1））/B2（1-）））

证明首先易得：

∑∞k=0ηk=∞，limt→∞1∑t-1k=0ηk12∑n（q+1）i=1E[‖xi（0）-x*‖2]=0 （17）

且有：lim supk→∞∑t-1k=0η2k∑t-1k=0ηk≤lim supk→∞ηk=

则

lim supt→∞1∑t-1k=0ηk12（n（q+1））（m+4）2 20 ∑t-1k=0η2k≤

12（n（q+1））（m+4）220 （18）

将推论1中的结论带入式（18）得

lim supt→∞1∑t-1k=0ηk（n（q+1））（m+

5）0∑t-1k=0ηkΩk≤2（n（q+1））（m+4）（m+

5）（2+n（q+1）α2（1-））20 （19）

利用定理1，将式（18）～式（19）带入式（10），最终可得：

lim supt→∞E[h（j（t））]-h*≤m 0∑n（q+1））i=1ui+

12（n（q+1））（m+4）220 +2（n（q+1））（m+

4）（m+5）×（2+n（q+1）α2（1-））20 （20）

事实上，定理2表明误差值lim supt→∞

E[h（j（t））]-h*是由两个方面来决定的。一方面是式（7）中所呈现的误差值，由于原始局部目标函数hi（x）的次梯度不存在或很难计算，用高斯近似函数hiμi （x）来代替hi（x），从而导致了该误差。另一方面，迭代步长在迭代的过程中逐步减小而ηk>0始终成立， 事实上， 当ηk=0时， 式（17）右侧变为lim supt→∞E[h（j（t））]-h*≤m0∑n（q+1）i=1μi，可见此时的误差值和约束集的维数m，lipschitz常数以及系统扩维所增加的时延节点数nq有关。

6结束语

在文献[2，7]的基础上，本文研究了时延情形下的分布式随机无梯度优化算法。首先给出时延情形下的分布式随机无梯度方法，然后通过系统扩维将时延优化问题转化为无时延优化问题，并利用转移矩阵的相关结论分析并证明了其收敛性。当然，还有一些问题有待解决，比如在无向图网络拓扑下中有时延的随机无梯度算法以及切换网络下的优化算法。

参考文献：

[1]NEDIA，OZDAGLAR A. Distributed Subgradient Methods for Multi-Agent Optimization[J].IEEE Transactions on Automatic Control，2009， 54（1）：48-61.

[2]NEDI A， OZDAGLAR A， PARRILO P A. Constrained Consensus and Optimization in Multi-Agent Networks [J]. IEEE Transactions on Automatic Control， 2010， 55（4）：922-938.

[3]SUNDHAR RAM S， NEDI A， V VEERAVALLI V. Distributed Stochastic Subgradient Projection Algorithms for Convex Optimization [J]. Journal of Optimization Theory & Applications， 2010， 147（3）：516-545.

[4]D Y，S X，H Z.Distributed Primal—Dual Subgradient Method for Multiagent Optimization via Consensus Algorithms[J].IEEE Transactions on Systems Man & Cybernetics Part B Cybernetics， 2011，41（6）：1 715-1 724.

[5]DUCHI J， AGARWAL A， WAINWRIGHT M. Dual Averaging for Distributed Optimization： Convergence Analysis and Network Scaling [J]. IEEE Transactions on Automatic Control， 2012， 57（3）：592-606.

[6]TSIANOS K I， LAWLOR S， RABBAT M G. Push-Sum Distributed Dual Averaging for convex optimization [C]// IEEE Conference on Decision & Control， 2012.

[7]NESTEROV Y. Random gradient-free minimization of convex functions [J]. General Information， international association for research and teaching， 2011， 36（16）：1 112-1 142

[8]LI J， WU C， WU Z， et al. Gradient-free method for nonsmooth distributed optimization [J]. Journal of Global Optimization， 2015， 61：325-340.

[9]YUAN D， HO D W C. Randomized Gradient-Free Method for Multiagent Optimization Over Time-Varying Networks [J]. IEEE Transactions on Neural Networks & Learning Systems， 2015， 26：1 342-1 347.

[10]TSIANOS K I， RABBAT M G. Distributed consensus and optimization under communication delays[C]//49th Annual Allerton Conference on Communication， Control， and Computing， 2011：974-982.

[11]李德权，张晓倩. 时延情形下的分布式Push-sum 次梯度优化算法的研究[J].安徽理工大学学报（自然科学版），2015，35（2）：6-12.

[12]刘德进， 刘成林. 具有通信时延的离散时间二阶多个体网络的一致性问题[J]. 控制理论与应用， 2010， 27（8）：1 108-1 158.