关于ADMM的研究(一)

关于ADMM的研究(一)

最近在研究正则化框架如何应用在大数据平台上。找到了《Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers》这篇文章,感觉很适合现在的研究。下面转载的一篇博客,写的很细致,很有用。

业界一直在谈论大数据,对于统计而言,大数据其实意味着要不是样本量增加n,要不就是维度的增加p,亦或者两者同时增加,并且维度与样本量的增长速度呈线性或者指数型增长。在稀疏性的假设条件下,再加上一些正则性方法,统计学家可以证明各种加penalty的模型所给出的参数估计具有良好的统计性质,收敛速度也有保证,同时还会给出一些比较好的迭代算法,但是,他们并没有考虑真实环境下的所消耗的计算时间。虽然统计学家也希望尽量寻求迭代数目比较少的算法(比如one-step估计),但是面对真实的Gb级别以上的数据,很多时候我们还是无法直接用这些算法,原因是一般的硬件都无法支撑直接对所有数据进行运算的要求。如果想减少抽样误差,不想抽样,又想提高估计的精度,那么还是需要寻求其他思路,结合已有的模型思想来解决这些问题。在目前条件下,并行化、分布式计算是一种比较好的解决思路,利用多核和多机器的优势,这些好算法便可以大规模应用,处理大数据优势便体现出来了。对于统计而言,数据量越大当然信息越可能充分(假设冗余成分不是特别多),因为大样本性质本身就希望样本越多越好嘛。

本文是基于Stephen Boyd 2011年的文章《Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers》进行的翻译和总结。Boyd也给出了利用matlab的CVX包实现的多种优化问题的matlab示例

1. 优化的一些基本算法思想

ADMM算法并不是一个很新的算法,他只是整合许多不少经典优化思路,然后结合现代统计学习所遇到的问题,提出了一个比较一般的比较好实施的分布式计算框架。因此必须先要了解一些基本算法思想。

1.1 Dual Ascent

对于凸函数的优化问题,对偶上升法核心思想就是引入一个对偶变量,然后利用交替优化的思路,使得两者同时达到optimal。一个凸函数的对偶函数其实就是原凸函数的一个下界,因此可以证明一个较好的性质:在强对偶性假设下,即最小化原凸函数(primal)等价于最大化对偶函数(dual),两者会同时达到optimal。这种转化可以将原来很多的参数约束条件变得少了很多,以利于做优化。具体表述如下:

mins.t.f(x)Ax=bL(x,y)=f(x)+yT(Axb)g(y)=infxL(x,y)

在强对偶性的假设下,primal和dual问题同时达到最优。

x=argminxL(x,y)

因此,若对偶函数g(y)可导,便可以利用梯度上升法,交替更新参数,使得同时收敛到最优。迭代如下:

xk+1:yk+1:=argminxL(x,yk)(x-最小化步)=yk+αkg(y)=yk+αk(Axk+1b)(对偶变量更新,αk是步长)

g不可微的时候也可以将其转化下,成为一个所谓的subgradient的方法,虽然看起来不错,简单证明下即可知道xkyk同时可达到optimal,但是上述条件要求很苛刻:f(x)要求严格凸,并且要求α选择有比较合适。一般应用中都不会满足(比如f(x)是一个非零的仿射函数),因此dual ascent不会直接应用。

1.2 Dual Decomposition

虽然dual ascent方法有缺陷,要求有些严格,但是他有一个非常好的性质,当目标函数f是可分的(separable)时候(参数抑或feature可分),整个问题可以拆解成多个子参数问题,分块优化后汇集起来整体更新。这样非常有利于并行化处理。形式化阐述如下:

mins.t.f(x)=i=1Nfi(xi),xiRni,xRnAx=i=1NAixi=b,(对A矩阵按列切分开)L(x,y)=i=1NLi(xi,y)=i=1N(fi(xi)+yTAixi1NyTb)

因此可以看到其实下面在迭代优化时,x-minimization步即可以拆分为多个子问题的并行优化,对偶变量更新不变这对于feature特别多时还是很有用的。

xk+1i:yk+1:=argminxLi(xi,yk)(多个xi并行最小化步)=yk+αkg(y)=yk+αk(Axk+1b)(汇集整体的x,然后对偶变量更新)

对偶分解是非常经典的优化方法,可追溯到1960年代。但是这种想法对后面的分布式优化方法影响较大,比如近期的graph-structure优化问题。

1.3 Augmented Lagrangians and the Method of Multipliers

从上面可以看到dual ascent方法对于目标函数要求比较苛刻,为了放松假设条件,同时比较好优化,于是就有了Augmented Lagrangians方法,目的就是放松对于f(x)严格凸的假设和其他一些条件,同时还能使得算法更加稳健。

Lρ(x,y)=f(x)+yT(Axb)+ρ2Axb22mins.t.f(x)+ρ2Axb22Ax=b

从上面可以看到该问题等价于最初的问题,因为只要是可行解对目标函数就没有影响。但是加了后面的(ρ/2)Axb22惩罚项的好处是使得对偶函数gρ(y)=infxLρ(x,y)在更一般的条件下可导。计算过程与之前的dual ascent基本一样,除了最小化x时候加了扩增项。

xk+1yk+1=argminxLρ(x,yk)=yk+ρ(Axk+1b)

上述也称作method of multipliers,可能也是因为更新对偶变量y时步长由原来变化的αk转为固定的ρ了吧。该算法在即使f(x)不是严格凸或者取值为+情况都可以成立,适用面更广。同样可以简单证明primal变量x和对偶变量y可以同时达到最优。

虽然Augmented Lagrangians方法有优势,但也破坏了dual ascent方法的利用分解参数来并行的优势。当f是separable时,对于Augmented Lagrangians却是not separable的(因为平方项写成矩阵形式无法用之前那种分块形式),因此在xmin步时候无法并行优化多个参数xi。如何改进,继续下面的议题就可以慢慢发现改进思想的来源。

2. Alternating Direction Method of Multipliers(ADMM)

2.1 ADMM算法概述

为了整合dual ascent可分解性与method multiplers优秀的收敛性质,人们就又提出了改进形式的优化ADMM。目的就是想能分解原函数和扩增函数,以便于在对f更一般的假设条件下并行优化。ADMM从名字可以看到是在原来Method of Multipliers加了个Alternating Direction,可以大概猜想到应该是又想引入新变量,然后交叉换方向来交替优化。形式如下:

mins.t.f(x)+g(z)Ax+Bz=cLρ(x,z,y)=f(x)+g(z)+yT(Ax+Bzc)+(ρ/2)Ax+Bzc22

从上面形式确实可以看出,他的思想确实就是想把primal变量、目标函数拆分,但是不再像dual ascent方法那样,将拆分开的xi都看做是x的一部分,后面融合的时候还需要融合在一起,而是最先开始就将拆开的变量分别看做是不同的变量x和z,同时约束条件也如此处理,这样的好处就是后面不需要一起融合x和z,保证了前面优化过程的可分解性。于是ADMM的优化就变成了如下序贯型迭代(这正是被称作alternating direction的缘故):

xk+1zk+1yk+1=argminxLρ(x,zk,yk)=argminzLρ(xk+1,z,yk)=yk+ρ(Axk+1+Bzk+1c)

后面我们可以看到这种拆分思想非常适合统计学习中的1-norm等问题:loss + regulazition(注意:一定要保证z分解出来,ADMM借助的就是用一个z变量来简化问题,不管他是约束还是其他形式也罢,需要构造一个z出来,后面具体到细节问题我们会有更深的体会)。

为了简化形式,ADMM有一个scaled form形式,其实就是对对偶变量做了scaled处理。先定义每一步更新的残差为r=Ax+Bzc,于是稍加计算

yT(Ax+Bzc)+(ρ/2)Ax+Bzc22=yTr+(ρ/2)r22=(ρ/2)r+(1/ρ)y22(1/2ρ)y22=(ρ/2)r+u22(ρ/2)u22

此处u=(1/ρ)y称为scaled dual variable,并令每一步迭代的残差为rk=Axk+Bzkc,以及累计残差uk=u0+kj=1rj,于是ADMM形式就可以简化为如下形式

xk+1zk+1uk+1=argminxLρ(x,zk,yk)=argmin(f(x)+(ρ/2)Ax+Bzkc+uk22)=argminzLρ(xk+1,z,yk)=argmin(g(z)+(ρ/2)Axk+1+Bzc+uk)=uk+ρ(Axk+1+Bzk+1c)

写成这种形式有利于后面简化优化问题,当然可以不作任何处理。

2.2 ADMM算法性质和评价

(1)收敛性

关于收敛性,需要有两个假设条件:

  • f和g分别是扩展的实数函数Rn(Rm)R+,且是closed、proper和convex的;
  • 扩增的lagrangian函数L0有一个鞍点(saddle point);对于约束中的矩阵A,B都不需要满秩。

在此两个假设下,可以保证残差、目标函数、对偶变量的收敛性。

Note:实际应用而言,ADMM收敛速度是很慢的,类似于共轭梯度方法。迭代数十次后只可以得到一个acceptable的结果,与快速的高精度算法(Newton法,内点法等)相比收敛就慢很多了。因此实际应用的时候,其实会将ADMM与其他高精度算法结合起来,这样从一个acceptable的结果变得在预期时间内可以达到较高收敛精度。不过一般在大规模应用问题中,高精度的参数解对于预测效果没有很大的提高,因此实际应用中,短时间内一个acceptable的结果基本就可以直接应用预测了。

(2)停止准则

对于ADMM的能到到optimal的条件此处就不做赘述了,与基本的primal和dual feasibility 的条件差不多,即各primal variable的偏导和约束条件为0,从最优条件中可以得到所谓的对偶残差(dual residuals)和初始残差(primal residuals)形式:

sk+1rk+1=ρATB(zk+1zk)(dualresiduals)=Axk+1+Bzk+1c(primalresiduals)

相对而言,此处更难把握的其实是停止准则,因为收敛速度问题,要想获得一个还过得去可以拿来用的参数解,那么判断迭代停止还是比较重要的。实际应用中,一般都根据primal residuals和dual residuals足够小来停止迭代,阈值包含了绝对容忍度(absolute tolerance)和相对容忍度(relative tolerance),设置还是非常灵活和难把握的(貌似网上有不少人吐槽这个停止准则的不靠谱- -!),具体形式如下:

sk2ϵdualrk2ϵpri=n√ϵabs+ϵrelATyk2=p√ϵabs+ϵrelmax{Axk2,Bzk,c2}

上面的p√n√分别是维度和样本量。一般而言,相对停止阈值ϵrel=103或者104,绝对阈值的选取要根据变量取值范围来选取(咋选的呢?没说额,具体比例都不给说- -!)

另外一些细节问题,比如原来惩罚参数ρ是不变的,一些文献也做了一些可变的惩罚参数,目的是为了降低对于惩罚参数初始值的依赖性。不过变动的ρ会导致ADMM的收敛性证明比较困难,因此实际中假设经过一系列迭代后ρ也稳定,边可直接用固定的惩罚参数ρ了。还有其他问题,诸如x与z迭代顺序问题,实际操作下有所有不同,这些不是特别重要之处,可以忽略。其他与ADMM比较相关算法的有dual ADMM算法,distributed ADMM算法,还有整合了ADMM与proximal method of multiplier的算法

2.3 ADMM一般形式与部分具体应用

当构造了ADMM算法中的f,g,A,B后,便可直接应用该算法了。我们会经常遇到如下三种一般形式的问题

  • 二次目标优化项(quadratic objective terms);
  • 可分的目标函数和约束(separable objective and constraints);
  • 光滑目标函数项(smooth objective terms)。

为下面讨论的方便,下面仅写出x-update的形式,根据ADMM简化形式,z-update对称更新即可:

x+=argminx(f(x)+(ρ/2)Axv22),v=Bz+cu

上述更新x时候z和u都定下来,是个常数,z更新时后相同。

Proximity Operator(近邻算子)

上述形式有种特殊情况:当A=I时,即约束条件没有x的线性组合形式,只是对于x的可行区域进行限制。这种问题相当常见,目前统计学习也有不少类似的高维优化问题。此时x-update如下

x+=argminx(f(x)+(ρ/2)xv22),v=Bz+cu

上述右边可以写成v的函数proxf,ρ(v)被称作带惩罚ρ的f的proximity operator(通常称作proximal minimization,近邻最小化),在变分分析中,还被称作f的Moreau-Yosida正则化。如果f形式很简单,可以写出x-update的解析解,比如f是非空的凸包C上的示性函数,那么x-update就可以直接写成投影形式

x+=argminx(f(x)+(ρ/2)xv22)=ΠC(v)

投影与惩罚参数ρ无关。若f是非负象限Rn+的投影,则直接有x+=(v)+

下面再谈谈上述提到的三种一般形式的优化问题。

(1)Quadratic Objective Terms

假设f是如下(凸)的二次函数

f(x)=12xTPx+qTx+r

P是对称的半正定矩阵PSn+。这种形式问题也包含了f是线性或者常数的特殊情况。若P+ρATA可逆,那么x-update步求个导即有如下的显示解,是v的仿射函数

x+=(P+ρATA)1(ρATvq)

因此在x-minnimiztion步只需要做两个矩阵运算即可,求逆与乘积,选用合适的线性运算库即可以得到不错的计算性能。当然还可以利用一些矩阵分解技巧,这个要看矩阵大小和稀疏程度。因为对于Fx=g,可以将F=F1F2Fk,然后Fizi=zi1,z1=F11g,x=zk,这样会更节省计算时间。其他矩阵计算技巧,基本都是如何对矩阵大规模求解,利用矩阵的稀疏性、缓存分解等来提高性能。此处不赘述,有个很重要的求逆的定理很有用:

(P+ρATA)1=P1ρP1AT(I+ρAP1AT)1AP1

如果对于上述二次函数受限于某仿射集x-update步就更复杂些,如

f(x)=12xTPx+qTx+rdormf={xFx=g}

x-update还有个重要的KKT方程可用:

(P+ρIFFT0)(xk+1v)+(qρ(zkuk)g)=0

(2)Smooth Objective Terms

f光滑时,那么求导即成为可能了。对于一些非线性优化问题,包含梯度算法等方法的L-BFGS算法可以用。对于该算法有些小技巧如下:

  • 早终止(early termination):当f(x)+(ρ/2)|Axv|22梯度很小时,早点终止迭代,否则后面就很慢了。
  • 热启动(warm start):即启动迭代时,利用之前迭代过的值带入即可。

(3)Separable objective and constraints 可分函数和约束对于并行计算和分布式计算来说是一个好消息。如果ATA是分块的对角阵,那么约束中|Ax|22也是可分的,则扩增的拉格朗日函数Lρ也是可分的。(注意,此处是指函数中的参数可分成小子块,而不是说数据可分。)下面有一个很重要的例子,即soft thresholding(针对l1+l2问题):

f(x)=λ|x|1,λ>0,并且A=I时,那么x-update就变成了

x+=argminx(λxi+(ρ/2)xv22)

这种形式很常见在目前的高维统计中,虽然第一项在0处不可导,但是也有解析解,被称作软阈值(soft thresholding),也被称作压缩算子(shrinkage operator)。

x+i=Sλ/ρ(vi),Sk(a)=⎧⎩⎨⎪⎪ak0,a+k,a>k|a|ka<kSk(a)=(1k|a|)+a

在优化领域,软阈值被称作是1-norm问题的近邻算子(proximity operator)。

3. 一些具体优化应用

3.1受约束的凸优化问题

一般的受约束的凸优化问题可以写成如下形式

mins.tf(x)xC

此类问题可以写成ADMM形式

mins.tf(x)+g(z)xz=0Lρ(x,z,u)=f(x)+g(z)+(ρ/2)xz+u22

其中的g函数即C的示性函数,上述是scaled形式,那么具体算法就是

xk+1zk+1uk+1=argmin(f(x)+(ρ/2)xzk+uk22)=ΠC(xk+1+uk)=uk+xk+1zk+1

则上述xmin就变成了一个具体的受约束的优化问题。比如对于经典的二次规划问题(QP)

mins.t12xTPx+qTxAx=b,x0

写成ADMM形式

mins.tf(x)+g(z)xz=0f(x)g(z)=12xTPx+qTx,dormf={x|Ax=b}=I(ΠRn+(z))

即受约束的区域就是{xx0}g是向非负象限投影的示性函数。而x-update就变成了之前在Quadratic Objective Terms中谈到的f(x)有仿射集定义域的优化问题,根据KKT条件即可写出来x-update更新的形式,参见2.3节。

如果上述对x限制不是限制x0上,而是一个锥约束(conic constraint)xK,那么x-update不变,继续上述KKT方程,而只需要变一下z-update,将向Rn+投影改成向K投影。比如将上述约束改成{Ax=b,x§n+},即x属于半正定空间,那么向(S^n_{+})投影就变成了一个半正定问题,利用特征值分解可以完成。这种受约束的凸优化问题的形式化对后续许多问题,特别是我们很关注的1-norm问题很重要,基本上都是转化成这种形式来直接应用ADMM算法,所以这里要好好把握其核心思想和形式。

虽然我对优化不在行,但是感觉优化问题还是挺有意思的,下面是一个经典问题,即找到两个非空凸包的交集中的一点。该算法都可以追溯到1930年代的Neumann交替投影算法(alternating projections algorithm):

xk+1zk+1=ΠC(zk)=ΠD(xk+1)

ΠC,ΠD分别是两个集合的欧式空间投影。写成ADMM形式就是

xk+1zk+1uk+1=ΠC(zkuk)=ΠD(xk+1+uk)=uk+xk+1zk+1

上述问题还可推广至找到N个非空凸包交集中一个点的问题,这样其实在x步是可以并行来做的,于是就有

xk+1izk+1uk+1i=ΠAi(zkuki)=1Ni=1N(xk+1i+uki)=uki+xk+1izk+1ui收敛均趋向于0,zk+1=x¯k+1xk+1iuk+1i=ΠAi(x¯kuki)=uki+(xk+1ix¯k+1)

3.2 1-norm问题

高维统计理论的发展,如果要追溯起来我觉得可以从Lasso解法算起,类似的思想在往前追可能是Huber相关的工作。是对于lasso问题,由于当年大家还没搞清楚lasso和boosting之间关系,对于sparsity性质不了解,谁也不知道如何很好地解决这个问题。直到后面Efron提出了LARS算法,对两者的路径解相似性做了很好的阐述,于是后面关于变量选择,关于basis-pursuit,compressed sensing,sparse graphical models等各种新问题的产生,随后各种优化算法也随之涌现出来,诸如Gradient Projection, Proximal methods,ADMM (Alternating Direction Method of Multipliers), (Split) Bregman methods,Nesterov’s method。不过要能够大规模部署1-norm的解决方案,那么这些算法中ADMM可能是首选。此处1-norm问题并不仅仅指Lasso问题,包含了多种1-norm类型问题。下面均介绍下。

之所以说ADMM适合机器学习和统计学习的优化问题,因为大部分机器学习问题基本都是“损失函数+正则项”形式,这种分法恰好可以套用到ADMM的框架f(x)+g(z)。因此结合ADMM框架基本可以解决很多已有的问题,以及利用1-norm构造的新的优化问题。下面将先介绍非分布式计算的版本,后面会单开一节来介绍如何分布式计算。

(1)Least Absolute Deviations

先从一个简单的问题开始。在稳健估计中,LAD是一个应用很广的模型,相对于直接优化平方和损失|Axb|22,优化绝对损失|Axb|1,它的抗噪性能更好。在ADMM框架下,往之前的受约束的凸优化问题靠拢,这个问题有简单的迭代算法

mins.t.z1Axb=zletf(x)=0,g(z)=z1xk+1zk+1uk+1=(ATA)1AT(b+zkuk)=S1/ρ(Axk+1b+uk)=uk+Axk+1zk+1b

(2)Huber fitting

Huber问题与上面的其实差不多,只是损失函数形式不同,换成了Huber惩罚函数

mins.t.ghub(z)Axb=z,ghub(z)=⎧⎩⎨z2/2,|z|12|z|1|z|>1

因此与LAD除了z-update不在是proximity operator(或称作软阈值)之外,其余均是相同的

zk+1=ρ1+ρ(Axk+1b+uk)+11+ρS1+1/ρ(Axk+1b+uk)

看着像是proximity operator与一个残差的加权。

LAD和Huber fitting这种问题只是一些传统损失不加正则项的ADMM化,注意一定要构造个z出来即可,x可以基本不用管,总是需要解的,下面的带有正则项的优化问题,ADMM形式就会更明显。

(3)Basis Pursuit

基追踪法师系数信号处理的一种重要方法。目的是想找到一组稀疏基可以完美恢复信号,换套话说就是为一个线性方程系统找到一个稀疏解。原始形式如下,与lasso有些像:

mins.t.x1Ax=b

修改成ADMM形式,注意往之前受约束的凸优化问题的那种形式回套,将1看做约束,然后构造带定义域的f(x),于是就有解

mins.t.f(x)+z1xz=0f(x)=I({xRn|Ax=b})indicator functionxk+1zk+1uk+1=Π(zkuk)=S1/ρ(Axk+1+uk)=uk+xk+1zk+1

其中Π(zkuk)是向一个线性约束的欧式空间中投影xRnAx=b,这也是有直接的显示解的

xk+1=(IAT(ATA)1A)(zuk)+AT(AAT)1b

对于矩阵求逆、分解等用之前矩阵那些小技巧即可加快计算,节省计算资源。

最近还有一类算法来解决1问题,被称作Bregman iteration methods,对于基追踪相关问题,加正则项的Bregman iteration就是method of multiplier,而所谓的split Bregman iteration就等同于 ADMM。我没有继续深究,应该就是类似于并行化的ADMM算法来解决基追踪问题。

(4)一般化的损失函数 + 1正则项问题

这类问题在高维统计开始时便是一个非常重要的问题,而即使到了现在也是一个非常重要的问题,比如group lasso,generalized lasso,高斯图模型,Tensor型图模型,与图相关的1问题等算法的开发,都可以在此框架上直接应用和实施,这正是ADMM一个优势所在,便于快速实施,也便于可能的大规模分布式部署。

minl(x)+λx1,mins.t.l(x)+g(z)=l(x)+λz1xz=0xk+1zk+1uk+1=argminx(l(x)+(ρ/2)xzk+uk22)=S1/ρ(Axk+1+uk)=uk+xk+1zk+1

可以看到与Basis Pursuit解法只是在x-update上有区别:Basis Pursuit是构造出来一个投影函数f(x),而一般化的损失函数f(x)+1正则项问题,用ADMM就更为自然。所以很适合作为框架来解决这一类问题:广义线性模型(普通线性、logistic回归、possion回归、softmax回归)+正则项;广义可加模型+正则项;似然函数(高斯图方向)+正则项。

  • Lasso:f(x)=12|Axb|22,于是利用ADMM算法,x-update的解析解就是xk+1=(ATA+ρI)1(ATb+ρ(zkuk));于是x-update看起来是个岭回归了,因此ADMM对于lasso可以看做迭代的使用岭回归。至于矩阵求逆那些,利用之前的矩阵小技巧解决。
  • Generalized lasso:这个问题可能不是那么为众人所熟悉,他是Tibs的儿子搞出来的框罗类似fused lasso这种事先定义好的线性变化的惩罚项的模型,损失函数是平方损失,而惩罚变成了一个特殊的参数线性组合

min12Axb22+λFx1 1d fused lasso,A=IFij=⎧⎩⎨⎪⎪110j=i+1j=iotherwise

min12xb22+λi=1n1|xi+1xi|A=I,F二阶差分矩阵,则被称作L1 trend filtering

若将上述这种写成ADMM形式,同样可以放到ADMM算法框架中解决

mins.t.12Axb22+λz1Fxz=0xk+1zk+1uk+1=(ATA+ρFTF)1(ATb+ρFT(zkuk))=S1/ρ(Axk+1b+uk)=uk+Fxk+1zk+1b
  • Group lasso:graph lasso问题应用比较广,对不同组的参数同时进行惩罚,进行一组组参数的挑选,故曰group lasso。不同于lasso,其正则项变成了Ni=1|xi|2,xiRni,lasso其实是group lasso的一种特殊形式。正则项并不是完全可分的。此时只是z-update变成了block的软阈值形式
zk+1i=Sλ/rho(xk+1i+uk),i=1,,NSk(a)=(1ka2)+a,S(0)=0

这种形式还可以扩展到group间有重合的情况,即化成N可能存在重合的组Gi1,,n。一般来说这种问题会非常难解决,但是对于ADMM算法只需要换下形式就很直接(x,z互换,会变成后面非常重要的一致性优化问题(consensus optimization),局部xi与全局真解z子集z^i的对应。)

mins.t.12Azb22+λi=1Nxi2,xiR|Gi|xiz^i=0,i=1,,N
  • Sparse Gaussian graph model:对于稀疏高斯图,熟悉该问题的人知道这其实是lasso的图上的推广,损失函数写成似然函数的负数即可l(x)=tr(SX)logdetX,XSn++。于是原来向量的操作就变成了矩阵操作,ADMM算法也有点变化:
Xk+1Zk+1Uk+1=argminX(tr(SX)logdetX+ρ2XZk+UkF)=argminZ(λZ1+ρ2Xk+1Z+UkF)=Uk+Xk+1Zk+1

上述算法继续化简,对于z-update做逐个元素软阈值操作即可Zk+1ij=Sλ/ρ(XK+1ij+Ukij)。对于x-update也类似操作,直接求导一阶导为0,移项后对对称矩阵做特征值分解即可

ρXX1=ρ(ZkUk)S=QΛQT,QQT=I,Λ=diag(λ1,,λn) ρX^X^1=Λ,X^=QTXQ

由于Λ是对角阵,对于每个对角元素来说,上述问题就是解一个二次方程,解方程后,再将X^变化成X即可

X^ii=λi+λ2i+4ρ−−−−−−√2ρX=QX^QT

总之,上述跟1相关的问题,基本都可以纳入ADMM框架,并且可以快速求解。

posted @ 2015-07-09 20:00  菜鸡一枚  阅读(15969)  评论(0编辑  收藏  举报