ADMM算法
https://joegaotao.github.io/2014/02/11/admm-stat-compute/
分布式計算、統(tǒng)計學(xué)習(xí)與ADMM算法
?
在整理舊電腦時，才發(fā)現(xiàn)13年下半年電腦里有不少殘文。老師說，東西擱下了再拿起來花費的時間和之前可能差不多。我一眼看過去這篇關(guān)于分布式計算的文章，貌似還真的沒有了當(dāng)時理解的深度和感覺。當(dāng)時還想利用ADMM算法，把統(tǒng)計中常見的帶懲罰的高維問題在此框架下用R重寫一下，但是中途多種事情一耽擱，就早已拋之腦后。看來任何事情，真的還是需要堅持，哪怕?lián)茳c時間都是好的。先把一篇殘文扔出來祭奠下過去的13年吧。公式多文字長，慎入！
業(yè)界一直在談?wù)摯髷?shù)據(jù)，對于統(tǒng)計而言，大數(shù)據(jù)其實意味著要不是樣本量增加n→∞n→∞，要不就是維度的增加p→∞p→∞，亦或者兩者同時增加，并且維度與樣本量的增長速度呈線性或者指數(shù)型增長。在稀疏性的假設(shè)條件下，再加上一些正則性方法，統(tǒng)計學(xué)家可以證明各種加penalty的模型所給出的參數(shù)估計具有良好的統(tǒng)計性質(zhì)，收斂速度也有保證，同時還會給出一些比較好的迭代算法，但是，他們并沒有考慮真實環(huán)境下的所消耗的計算時間。雖然統(tǒng)計學(xué)家也希望盡量尋求迭代數(shù)目比較少的算法（比如one-step估計），但是面對真實的Gb級別以上的數(shù)據(jù)，很多時候我們還是無法直接用這些算法，原因是一般的硬件都無法支撐直接對所有數(shù)據(jù)進行運算的要求。如果想減少抽樣誤差，不想抽樣，又想提高估計的精度，那么還是需要尋求其他思路，結(jié)合已有的模型思想來解決這些問題。在目前條件下，并行化、分布式計算是一種比較好的解決思路，利用多核和多機器的優(yōu)勢，這些好算法便可以大規(guī)模應(yīng)用，處理大數(shù)據(jù)優(yōu)勢便體現(xiàn)出來了。對于統(tǒng)計而言，數(shù)據(jù)量越大當(dāng)然信息越可能充分（假設(shè)冗余成分不是特別多），因為大樣本性質(zhì)本身就希望樣本越多越好嘛。
本文是基于Stephen Boyd 2011年的文章《Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers》進行的翻譯和總結(jié)。Boyd也給出了利用matlab的CVX包實現(xiàn)的多種優(yōu)化問題的matlab示例。
1. 優(yōu)化的一些基本算法思想
ADMM算法并不是一個很新的算法，他只是整合許多不少經(jīng)典優(yōu)化思路，然后結(jié)合現(xiàn)代統(tǒng)計學(xué)習(xí)所遇到的問題，提出了一個比較一般的比較好實施的分布式計算框架。因此必須先要了解一些基本算法思想。
1.1 Dual Ascent
對于凸函數(shù)的優(yōu)化問題，對偶上升法核心思想就是引入一個對偶變量，然后利用交替優(yōu)化的思路，使得兩者同時達(dá)到optimal。一個凸函數(shù)的對偶函數(shù)其實就是原凸函數(shù)的一個下界，因此可以證明一個較好的性質(zhì)：在強對偶性假設(shè)下，即最小化原凸函數(shù)（primal）等價于最大化對偶函數(shù)（dual），兩者會同時達(dá)到optimal。這種轉(zhuǎn)化可以將原來很多的參數(shù)約束條件變得少了很多，以利于做優(yōu)化。具體表述如下：
minf(x)s.t.Ax=b?L(x,y)=f(x)+yT(Ax?b)對偶函數(shù)（下界）?g(y)=infxL(x,y)minf(x)s.t.Ax=b?L(x,y)=f(x)+yT(Ax?b)?對偶函數(shù)（下界）g(y)=infxL(x,y)
在強對偶性的假設(shè)下，primal和dual問題同時達(dá)到最優(yōu)。
x?=argminxL(x,y?)x?=arg?minxL(x,y?)
因此，若對偶函數(shù)g(y)g(y)可導(dǎo)，便可以利用梯度上升法，交替更新參數(shù)，使得同時收斂到最優(yōu)。迭代如下：
xk+1:=argminxL(x,yk)(x-最小化步)yk+1:=yk+αk?g(y)=yk+αk(Axk+1?b)(對偶變量更新，αk是步長)xk+1:=arg?minxL(x,yk)(x-最小化步)yk+1:=yk+αk?g(y)=yk+αk(Axk+1?b)(對偶變量更新，αk是步長)
當(dāng)gg不可微的時候也可以將其轉(zhuǎn)化下，成為一個所謂的subgradient的方法，雖然看起來不錯，簡單證明下即可知道xkxk和ykyk同時可達(dá)到optimal，但是上述條件要求很苛刻：f(x)f(x)要求嚴(yán)格凸，并且要求αα選擇有比較合適。一般應(yīng)用中都不會滿足（比如f(x)f(x)是一個非零的仿射函數(shù)），因此dual ascent不會直接應(yīng)用。
1.2 Dual Decomposition
雖然dual ascent方法有缺陷，要求有些嚴(yán)格，但是他有一個非常好的性質(zhì)，當(dāng)目標(biāo)函數(shù)ff是可分的（separable）時候（參數(shù)抑或feature可分），整個問題可以拆解成多個子參數(shù)問題，分塊優(yōu)化后匯集起來整體更新。這樣非常有利于并行化處理。形式化闡述如下：
minf(x)=∑Ni=1fi(xi),xi∈Rni,x∈Rns.t.Ax=∑Ni=1Aixi=b,(對A矩陣按列切分開)?L(x,y)=N∑i=1Li(xi,y)=N∑i=1(fi(xi)+yTAixi?1NyTb)minf(x)=∑i=1Nfi(xi),xi∈Rni,x∈Rns.t.Ax=∑i=1NAixi=b,(對A矩陣按列切分開)?L(x,y)=∑i=1NLi(xi,y)=∑i=1N(fi(xi)+yTAixi?1NyTb)
因此可以看到其實下面在迭代優(yōu)化時，xx-minimization步即可以拆分為多個子問題的并行優(yōu)化，對偶變量更新不變這對于feature特別多時還是很有用的。
xk+1i:=argminxLi(xi,yk)(多個xi并行最小化步)yk+1:=yk+αk?g(y)=yk+αk(Axk+1?b)(匯集整體的x，然后對偶變量更新)xik+1:=arg?minxLi(xi,yk)(多個xi并行最小化步)yk+1:=yk+αk?g(y)=yk+αk(Axk+1?b)(匯集整體的x，然后對偶變量更新)
對偶分解是非常經(jīng)典的優(yōu)化方法，可追溯到1960年代。但是這種想法對后面的分布式優(yōu)化方法影響較大，比如近期的graph-structure優(yōu)化問題。
1.3 Augmented Lagrangians and the Method of Multipliers
從上面可以看到dual ascent方法對于目標(biāo)函數(shù)要求比較苛刻，為了放松假設(shè)條件，同時比較好優(yōu)化，于是就有了Augmented Lagrangians方法，目的就是放松對于f(x)f(x)嚴(yán)格凸的假設(shè)和其他一些條件，同時還能使得算法更加穩(wěn)健。
Lρ(x,y)=f(x)+yT(Ax?b)+ρ2∥Ax?b∥22?minf(x)+ρ2∥Ax?b∥22s.t.Ax=bLρ(x,y)=f(x)+yT(Ax?b)+ρ2‖Ax?b‖22?minf(x)+ρ2‖Ax?b‖22s.t.Ax=b
從上面可以看到該問題等價于最初的問題，因為只要是可行解對目標(biāo)函數(shù)就沒有影響。但是加了后面的(ρ/2)∥Ax?b∥22(ρ/2)‖Ax?b‖22懲罰項的好處是使得對偶函數(shù)gρ(y)=infxLρ(x,y)gρ(y)=infxLρ(x,y)在更一般的條件下可導(dǎo)。計算過程與之前的dual ascent基本一樣，除了最小化xx時候加了擴增項。
xk+1=argminxLρ(x,yk)yk+1=yk+ρ(Axk+1?b)xk+1=arg?minxLρ(x,yk)yk+1=yk+ρ(Axk+1?b)
上述也稱作method of multipliers，可能也是因為更新對偶變量yy時步長由原來變化的αkαk轉(zhuǎn)為固定的ρρ了吧。該算法在即使f(x)f(x)不是嚴(yán)格凸或者取值為+∞+∞情況都可以成立，適用面更廣。同樣可以簡單證明primal變量xx和對偶變量yy可以同時達(dá)到最優(yōu)。
雖然Augmented Lagrangians方法有優(yōu)勢，但也破壞了dual ascent方法的利用分解參數(shù)來并行的優(yōu)勢。當(dāng)ff是separable時，對于Augmented Lagrangians卻是not separable的（因為平方項寫成矩陣形式無法用之前那種分塊形式），因此在x?minx?min步時候無法并行優(yōu)化多個參數(shù)xixi。如何改進，繼續(xù)下面的議題就可以慢慢發(fā)現(xiàn)改進思想的來源。
2. Alternating Direction Method of Multipliers(ADMM)
2.1 ADMM算法概述
為了整合dual ascent可分解性與method multiplers優(yōu)秀的收斂性質(zhì)，人們就又提出了改進形式的優(yōu)化ADMM。目的就是想能分解原函數(shù)和擴增函數(shù)，以便于在對ff更一般的假設(shè)條件下并行優(yōu)化。ADMM從名字可以看到是在原來Method of Multipliers加了個Alternating Direction，可以大概猜想到應(yīng)該是又想引入新變量，然后交叉換方向來交替優(yōu)化。形式如下：
minf(x)+g(z)s.t.Ax+Bz=c?Lρ(x,z,y)=f(x)+g(z)+yT(Ax+Bz?c)+(ρ/2)∥Ax+Bz?c∥22minf(x)+g(z)s.t.Ax+Bz=c?Lρ(x,z,y)=f(x)+g(z)+yT(Ax+Bz?c)+(ρ/2)‖Ax+Bz?c‖22
從上面形式確實可以看出，他的思想確實就是想把primal變量、目標(biāo)函數(shù)拆分，但是不再像dual ascent方法那樣，將拆分開的xixi都看做是xx的一部分，后面融合的時候還需要融合在一起，而是最先開始就將拆開的變量分別看做是不同的變量xx和zz，同時約束條件也如此處理，這樣的好處就是后面不需要一起融合xx和zz，保證了前面優(yōu)化過程的可分解性。于是ADMM的優(yōu)化就變成了如下序貫型迭代（這正是被稱作alternating direction的緣故）：
xk+1=argminxLρ(x,zk,yk)zk+1=argminzLρ(xk+1,z,yk)yk+1=yk+ρ(Axk+1+Bzk+1?c)xk+1=arg?minxLρ(x,zk,yk)zk+1=arg?minzLρ(xk+1,z,yk)yk+1=yk+ρ(Axk+1+Bzk+1?c)
后面我們可以看到這種拆分思想非常適合統(tǒng)計學(xué)習(xí)中的?1?1-norm等問題：loss + regulazition（注意：一定要保證zz分解出來，ADMM借助的就是用一個zz變量來簡化問題，不管他是約束還是其他形式也罷，需要構(gòu)造一個zz出來，后面具體到細(xì)節(jié)問題我們會有更深的體會）。
為了簡化形式，ADMM有一個scaled form形式，其實就是對對偶變量做了scaled處理。先定義每一步更新的殘差為r=Ax+Bz?cr=Ax+Bz?c，于是稍加計算
yT(Ax+Bz?c)+(ρ/2)∥Ax+Bz?c∥22=yTr+(ρ/2)∥r∥22=(ρ/2)∥r+(1/ρ)y∥22?(1/2ρ)∥y∥22=(ρ/2)∥r+u∥22?(ρ/2)∥u∥22yT(Ax+Bz?c)+(ρ/2)‖Ax+Bz?c‖22=yTr+(ρ/2)‖r‖22=(ρ/2)‖r+(1/ρ)y‖22?(1/2ρ)‖y‖22=(ρ/2)‖r+u‖22?(ρ/2)‖u‖22
此處u=(1/ρ)yu=(1/ρ)y稱為scaled dual variable，并令每一步迭代的殘差為rk=Axk+Bzk?crk=Axk+Bzk?c，以及累計殘差uk=u0+∑kj=1rjuk=u0+∑j=1krj，于是ADMM形式就可以簡化為如下形式
xk+1=argminxLρ(x,zk,yk)=argmin(f(x)+(ρ/2)∥Ax+Bzk?c+uk∥22)zk+1=argminzLρ(xk+1,z,yk)=argmin(g(z)+(ρ/2)∥Axk+1+Bz?c+uk∥)uk+1=uk+Axk+1+Bzk+1?cxk+1=arg?minxLρ(x,zk,yk)=arg?min(f(x)+(ρ/2)‖Ax+Bzk?c+uk‖22)zk+1=arg?minzLρ(xk+1,z,yk)=arg?min(g(z)+(ρ/2)‖Axk+1+Bz?c+uk‖)uk+1=uk+Axk+1+Bzk+1?c
寫成這種形式有利于后面簡化優(yōu)化問題，當(dāng)然可以不作任何處理。
2.2 ADMM算法性質(zhì)和評價
（1）收斂性
關(guān)于收斂性，需要有兩個假設(shè)條件：
- ff和gg分別是擴展的實數(shù)函數(shù)Rn(Rm)→R?{+∞}Rn(Rm)→R?{+∞}，且是closed、proper和convex的；
- 擴增的lagrangian函數(shù)L0L0有一個鞍點（saddle point）；對于約束中的矩陣A,BA,B都不需要滿秩。
在此兩個假設(shè)下，可以保證殘差、目標(biāo)函數(shù)、對偶變量的收斂性。
Note：實際應(yīng)用而言，ADMM收斂速度是很慢的，類似于共軛梯度方法。迭代數(shù)十次后只可以得到一個acceptable的結(jié)果，與快速的高精度算法（Newton法，內(nèi)點法等）相比收斂就慢很多了。因此實際應(yīng)用的時候，其實會將ADMM與其他高精度算法結(jié)合起來，這樣從一個acceptable的結(jié)果變得在預(yù)期時間內(nèi)可以達(dá)到較高收斂精度。不過一般在大規(guī)模應(yīng)用問題中，高精度的參數(shù)解對于預(yù)測效果沒有很大的提高，因此實際應(yīng)用中，短時間內(nèi)一個acceptable的結(jié)果基本就可以直接應(yīng)用預(yù)測了。
（2）停止準(zhǔn)則
對于ADMM的能到到optimal的條件此處就不做贅述了，與基本的primal和dual feasibility 的條件差不多，即各primal variable的偏導(dǎo)和約束條件為0，從最優(yōu)條件中可以得到所謂的對偶?xì)埐?#xff08;dual residuals）和初始?xì)埐?#xff08;primal residuals）形式：
sk+1=ρATB(zk+1?zk)(dualresiduals)rk+1=Axk+1+Bzk+1?c(primalresiduals)sk+1=ρATB(zk+1?zk)(dualresiduals)rk+1=Axk+1+Bzk+1?c(primalresiduals)
相對而言，此處更難把握的其實是停止準(zhǔn)則，因為收斂速度問題，要想獲得一個還過得去可以拿來用的參數(shù)解，那么判斷迭代停止還是比較重要的。實際應(yīng)用中，一般都根據(jù)primal residuals和dual residuals足夠小來停止迭代，閾值包含了絕對容忍度（absolute tolerance）和相對容忍度（relative tolerance），設(shè)置還是非常靈活和難把握的（貌似網(wǎng)上有不少人吐槽這個停止準(zhǔn)則的不靠譜- -！），具體形式如下：
∥sk∥2≤?dual=√n?abs+?rel∥ATyk∥2∥rk∥2≤?pri=√p?abs+?relmax{∥Axk∥2,∥Bzk∥,∥c∥2}‖sk‖2≤?dual=n?abs+?rel‖ATyk‖2‖rk‖2≤?pri=p?abs+?relmax{‖Axk‖2,‖Bzk‖,‖c‖2}
上面的√pp和√nn分別是維度和樣本量。一般而言，相對停止閾值?rel=10?3?rel=10?3或者10?410?4，絕對閾值的選取要根據(jù)變量取值范圍來選取（咋選的呢？沒說額，具體比例都不給說- -！）
另外一些細(xì)節(jié)問題，比如原來懲罰參數(shù)ρρ是不變的，一些文獻(xiàn)也做了一些可變的懲罰參數(shù)，目的是為了降低對于懲罰參數(shù)初始值的依賴性。不過變動的ρρ會導(dǎo)致ADMM的收斂性證明比較困難，因此實際中假設(shè)經(jīng)過一系列迭代后ρρ也穩(wěn)定，邊可直接用固定的懲罰參數(shù)ρρ了。還有其他問題，諸如xx與zz迭代順序問題，實際操作下有所有不同，這些不是特別重要之處，可以忽略。其他與ADMM比較相關(guān)算法的有dual ADMM算法，distributed ADMM算法，還有整合了ADMM與proximal method of multiplier的算法
2.3 ADMM一般形式與部分具體應(yīng)用
當(dāng)構(gòu)造了ADMM算法中的f,g,A,Bf,g,A,B后，便可直接應(yīng)用該算法了。我們會經(jīng)常遇到如下三種一般形式的問題
- 二次目標(biāo)優(yōu)化項（quadratic objective terms）；
- 可分的目標(biāo)函數(shù)和約束（separable objective and constraints）；
- 光滑目標(biāo)函數(shù)項（smooth objective terms）。
為下面討論的方便，下面僅寫出xx-update的形式，根據(jù)ADMM簡化形式，zz-update對稱更新即可：
x+=argminx(f(x)+(ρ/2)∥Ax?v∥22),v=?Bz+c?ux+=arg?minx(f(x)+(ρ/2)‖Ax?v‖22),v=?Bz+c?u
上述更新xx時候zz和uu都定下來，是個常數(shù)，zz更新時后相同。
Proximity Operator（近鄰算子）
上述形式有種特殊情況：當(dāng)A=IA=I時，即約束條件沒有xx的線性組合形式，只是對于xx的可行區(qū)域進行限制。這種問題相當(dāng)常見，目前統(tǒng)計學(xué)習(xí)也有不少類似的高維優(yōu)化問題。此時xx-update如下
x+=argminx(f(x)+(ρ/2)∥x?v∥22),v=?Bz+c?ux+=arg?minx(f(x)+(ρ/2)‖x?v‖22),v=?Bz+c?u
上述右邊可以寫成vv的函數(shù)proxf,ρ(v)proxf,ρ(v)被稱作帶懲罰ρρ的ff的proximity operator（通常稱作proximal minimization，近鄰最小化），在變分分析中，還被稱作ff的Moreau-Yosida正則化。如果ff形式很簡單，可以寫出xx-update的解析解，比如ff是非空的凸包CC上的示性函數(shù)，那么xx-update就可以直接寫成投影形式
x+=argminx(f(x)+(ρ/2)∥x?v∥22)=ΠC(v)x+=arg?minx(f(x)+(ρ/2)‖x?v‖22)=ΠC(v)
投影與懲罰參數(shù)ρρ無關(guān)。若ff是非負(fù)象限Rn+R+n的投影，則直接有x+=(v)+x+=(v)+。
下面再談?wù)勆鲜鎏岬降娜N一般形式的優(yōu)化問題。
（1）Quadratic Objective Terms
假設(shè)ff是如下（凸）的二次函數(shù)
f(x)=12xTPx+qTx+rf(x)=12xTPx+qTx+r
PP是對稱的半正定矩陣P∈Sn+P∈S+n。這種形式問題也包含了ff是線性或者常數(shù)的特殊情況。若P+ρATAP+ρATA可逆，那么xx-update步求個導(dǎo)即有如下的顯示解，是vv的仿射函數(shù)
x+=(P+ρATA)?1(ρATv?q)x+=(P+ρATA)?1(ρATv?q)
因此在xx-minnimiztion步只需要做兩個矩陣運算即可，求逆與乘積，選用合適的線性運算庫即可以得到不錯的計算性能。當(dāng)然還可以利用一些矩陣分解技巧，這個要看矩陣大小和稀疏程度。因為對于Fx=gFx=g，可以將F=F1F2?FkF=F1F2?Fk，然后Fizi=zi?1,z1=F?11g,x=zkFizi=zi?1,z1=F1?1g,x=zk，這樣會更節(jié)省計算時間。其他矩陣計算技巧，基本都是如何對矩陣大規(guī)模求解，利用矩陣的稀疏性、緩存分解等來提高性能。此處不贅述，有個很重要的求逆的定理很有用：
(P+ρATA)?1=P?1?ρP?1AT(I+ρAP?1AT)?1AP?1(P+ρATA)?1=P?1?ρP?1AT(I+ρAP?1AT)?1AP?1
如果對于上述二次函數(shù)受限于某仿射集xx-update步就更復(fù)雜些，如
f(x)=12xTPx+qTx+rdormf={x∥Fx=g}f(x)=12xTPx+qTx+rdormf={x‖F(xiàn)x=g}
xx-update還有個重要的KKT方程可用：
(P+ρIFTF0)(xk+1v)+(q?ρ(zk?uk)?g)=0(P+ρIFTF0)(xk+1v)+(q?ρ(zk?uk)?g)=0
（2）Smooth Objective Terms
當(dāng)ff光滑時，那么求導(dǎo)即成為可能了。對于一些非線性優(yōu)化問題，包含梯度算法等方法的L-BFGS算法可以用。對于該算法有些小技巧如下：
- 早終止（early termination）：當(dāng)f(x)+(ρ/2)∥Ax?v∥22f(x)+(ρ/2)‖Ax?v‖22梯度很小時，早點終止迭代，否則后面就很慢了。
- 熱啟動（warm start）：即啟動迭代時，利用之前迭代過的值帶入即可。
（3）Separable objective and constraints?可分函數(shù)和約束對于并行計算和分布式計算來說是一個好消息。如果ATAATA是分塊的對角陣，那么約束中∥Ax∥22‖Ax‖22也是可分的，則擴增的拉格朗日函數(shù)LρLρ也是可分的。（注意，此處是指函數(shù)中的參數(shù)可分成小子塊，而不是說數(shù)據(jù)可分。）下面有一個很重要的例子，即soft thresholding（針對l1+l2l1+l2問題）:
當(dāng)f(x)=λ∥x∥1,λ>0f(x)=λ‖x‖1,λ>0，并且A=IA=I時，那么xx-update就變成了
x+=argminx(λ∥xi∥+(ρ/2)∥x?v∥22)x+=arg?minx(λ‖xi‖+(ρ/2)‖x?v‖22)
這種形式很常見在目前的高維統(tǒng)計中，雖然第一項在0處不可導(dǎo)，但是也有解析解，被稱作軟閾值（soft thresholding），也被稱作壓縮算子（shrinkage operator）。
x+i=Sλ/ρ(vi),→Sk(a)=?????a?k,a>k0,|a|≤ka+ka<?k→Sk(a)=(1?k|a|)+axi+=Sλ/ρ(vi),→Sk(a)={a?k,a>k0,|a|≤ka+ka<?k→Sk(a)=(1?k|a|)+a
在優(yōu)化領(lǐng)域，軟閾值被稱作是?1?1-norm問題的近鄰算子（proximity operator）。
3. 一些具體優(yōu)化應(yīng)用
3.1受約束的凸優(yōu)化問題
一般的受約束的凸優(yōu)化問題可以寫成如下形式
minf(x)s.tx∈Cminf(x)s.tx∈C
此類問題可以寫成ADMM形式
minf(x)+g(z)s.tx?z=0?Lρ(x,z,u)=f(x)+g(z)+(ρ/2)∥x?z+u∥22minf(x)+g(z)s.tx?z=0?Lρ(x,z,u)=f(x)+g(z)+(ρ/2)‖x?z+u‖22
其中的gg函數(shù)即CC的示性函數(shù)，上述是scaled形式，那么具體算法就是
xk+1=argmin(f(x)+(ρ/2)∥x?zk+uk∥22)zk+1=ΠC(xk+1+uk)uk+1=uk+xk+1?zk+1xk+1=arg?min(f(x)+(ρ/2)‖x?zk+uk‖22)zk+1=ΠC(xk+1+uk)uk+1=uk+xk+1?zk+1
則上述x?minx?min就變成了一個具體的受約束的優(yōu)化問題。比如對于經(jīng)典的二次規(guī)劃問題(QP)
min12xTPx+qTxs.tAx=b,x≥0min12xTPx+qTxs.tAx=b,x≥0
寫成ADMM形式
minf(x)+g(z)s.tx?z=0?f(x)=12xTPx+qTx,dormf={x|Ax=b}g(z)=I(ΠRn+(z))minf(x)+g(z)s.tx?z=0?f(x)=12xTPx+qTx,dormf={x|Ax=b}g(z)=I(ΠR+n(z))
即受約束的區(qū)域就是{x∣x≥0}{x∣x≥0}，gg是向非負(fù)象限投影的示性函數(shù)。而xx-update就變成了之前在Quadratic Objective Terms中談到的f(x)f(x)有仿射集定義域的優(yōu)化問題，根據(jù)KKT條件即可寫出來xx-update更新的形式，參見2.3節(jié)。
如果上述對xx限制不是限制x≥0x≥0上，而是一個錐約束（conic constraint）x∈Kx∈K，那么xx-update不變，繼續(xù)上述KKT方程，而只需要變一下zz-update，將向Rn+R+n投影改成向KK投影。比如將上述約束改成{Ax=b,x∈§n+}{Ax=b,x∈§+n}，即xx屬于半正定空間，那么向Sn+S+n投影就變成了一個半正定問題，利用特征值分解可以完成。這種受約束的凸優(yōu)化問題的形式化對后續(xù)許多問題，特別是我們很關(guān)注的?1?1-norm問題很重要，基本上都是轉(zhuǎn)化成這種形式來直接應(yīng)用ADMM算法，所以這里要好好把握其核心思想和形式。
雖然我對優(yōu)化不在行，但是感覺優(yōu)化問題還是挺有意思的，下面是一個經(jīng)典問題，即找到兩個非空凸包的交集中的一點。該算法都可以追溯到1930年代的Neumann交替投影算法（alternating projections algorithm）：
xk+1=ΠC(zk)zk+1=ΠD(xk+1)xk+1=ΠC(zk)zk+1=ΠD(xk+1)
ΠC,ΠDΠC,ΠD分別是兩個集合的歐式空間投影。寫成ADMM形式就是
xk+1=ΠC(zk?uk)zk+1=ΠD(xk+1+uk)uk+1=uk+xk+1?zk+1xk+1=ΠC(zk?uk)zk+1=ΠD(xk+1+uk)uk+1=uk+xk+1?zk+1
上述問題還可推廣至找到NN個非空凸包交集中一個點的問題，這樣其實在xx步是可以并行來做的，于是就有
xk+1i=ΠAi(zk?uki)zk+1=1NN∑i=1(xk+1i+uki)uk+1i=uki+xk+1i?zk+1?ui收斂均趨向于0,zk+1=ˉxk+1xk+1i=ΠAi(ˉxk?uki)uk+1i=uki+(xk+1i?ˉxk+1)xik+1=ΠAi(zk?uik)zk+1=1N∑i=1N(xik+1+uik)uik+1=uik+xik+1?zk+1?ui收斂均趨向于0,zk+1=xˉk+1xik+1=ΠAi(xˉk?uik)uik+1=uik+(xik+1?xˉk+1)
3.2??1?1-norm問題
高維統(tǒng)計理論的發(fā)展，如果要追溯起來我覺得可以從Lasso解法算起，類似的思想在往前追可能是Huber相關(guān)的工作。是對于lasso問題，由于當(dāng)年大家還沒搞清楚lasso和boosting之間關(guān)系，對于sparsity性質(zhì)不了解，誰也不知道如何很好地解決這個問題。直到后面Efron提出了LARS算法，對兩者的路徑解相似性做了很好的闡述，于是后面關(guān)于變量選擇，關(guān)于basis-pursuit，compressed sensing，sparse graphical models等各種新問題的產(chǎn)生，隨后各種優(yōu)化算法也隨之涌現(xiàn)出來，諸如Gradient Projection， Proximal methods，ADMM (Alternating Direction Method of Multipliers)， (Split) Bregman methods，Nesterov’s method。不過要能夠大規(guī)模部署?1?1-norm的解決方案，那么這些算法中ADMM可能是首選。此處?1?1-norm問題并不僅僅指Lasso問題，包含了多種?1?1-norm類型問題。下面均介紹下。
之所以說ADMM適合機器學(xué)習(xí)和統(tǒng)計學(xué)習(xí)的優(yōu)化問題，因為大部分機器學(xué)習(xí)問題基本都是“損失函數(shù)+正則項”形式，這種分法恰好可以套用到ADMM的框架f(x)+g(z)f(x)+g(z)。因此結(jié)合ADMM框架基本可以解決很多已有的問題，以及利用?1?1-norm構(gòu)造的新的優(yōu)化問題。下面將先介紹非分布式計算的版本，后面會單開一節(jié)來介紹如何分布式計算。
（1）Least Absolute Deviations
先從一個簡單的問題開始。在穩(wěn)健估計中，LAD是一個應(yīng)用很廣的模型，相對于直接優(yōu)化平方和損失∥Ax?b∥22‖Ax?b‖22，優(yōu)化絕對損失∥Ax?b∥1‖Ax?b‖1，它的抗噪性能更好。在ADMM框架下，往之前的受約束的凸優(yōu)化問題靠攏，這個問題有簡單的迭代算法
min∥z∥1s.t.Ax?b=z?letf(x)=0,g(z)=∥z∥1?xk+1=(ATA)?1AT(b+zk?uk)zk+1=S1/ρ(Axk+1?b+uk)uk+1=uk+Axk+1?zk+1?bmin‖z‖1s.t.Ax?b=z?letf(x)=0,g(z)=‖z‖1?xk+1=(ATA)?1AT(b+zk?uk)zk+1=S1/ρ(Axk+1?b+uk)uk+1=uk+Axk+1?zk+1?b
（2）Huber fitting
Huber問題與上面的其實差不多，只是損失函數(shù)形式不同，換成了Huber懲罰函數(shù)
minghub(z)s.t.Ax?b=z,ghub(z)={z2/2,|z|≤1|z|?12|z|>1minghub(z)s.t.Ax?b=z,ghub(z)={z2/2,|z|≤1|z|?12|z|>1
因此與LAD除了zz-update不在是proximity operator（或稱作軟閾值）之外，其余均是相同的
zk+1=ρ1+ρ(Axk+1?b+uk)+11+ρS1+1/ρ(Axk+1?b+uk)zk+1=ρ1+ρ(Axk+1?b+uk)+11+ρS1+1/ρ(Axk+1?b+uk)
看著像是proximity operator與一個殘差的加權(quán)。
LAD和Huber fitting這種問題只是一些傳統(tǒng)損失不加正則項的ADMM化，注意一定要構(gòu)造個zz出來即可，xx可以基本不用管，總是需要解的，下面的帶有正則項的優(yōu)化問題，ADMM形式就會更明顯。
（3）Basis Pursuit
基追蹤法師系數(shù)信號處理的一種重要方法。目的是想找到一組稀疏基可以完美恢復(fù)信號，換套話說就是為一個線性方程系統(tǒng)找到一個稀疏解。原始形式如下，與lasso有些像：
min∥x∥1s.t.Ax=bmin‖x‖1s.t.Ax=b
修改成ADMM形式，注意往之前受約束的凸優(yōu)化問題的那種形式回套，將?1?1看做約束，然后構(gòu)造帶定義域的f(x)f(x)，于是就有解
minf(x)+∥z∥1s.t.x?z=0f(x)=I({x∈Rn|Ax=b})indicator function?xk+1=Π(zk?uk)zk+1=S1/ρ(Axk+1+uk)uk+1=uk+xk+1?zk+1minf(x)+‖z‖1s.t.x?z=0f(x)=I({x∈Rn|Ax=b})indicator function?xk+1=Π(zk?uk)zk+1=S1/ρ(Axk+1+uk)uk+1=uk+xk+1?zk+1
其中Π(zk?uk)Π(zk?uk)是向一個線性約束的歐式空間中投影{x∈Rn∣Ax=b}{x∈Rn∣Ax=b}，這也是有直接的顯示解的
xk+1=(I?AT(ATA)?1A)(z?uk)+AT(AAT)?1bxk+1=(I?AT(ATA)?1A)(z?uk)+AT(AAT)?1b
對于矩陣求逆、分解等用之前矩陣那些小技巧即可加快計算，節(jié)省計算資源。
最近還有一類算法來解決?1?1問題，被稱作Bregman iteration methods，對于基追蹤相關(guān)問題，加正則項的Bregman iteration就是method of multiplier，而所謂的split Bregman iteration就等同于 ADMM。我沒有繼續(xù)深究，應(yīng)該就是類似于并行化的ADMM算法來解決基追蹤問題。
（4）一般化的損失函數(shù) +??1?1正則項問題
這類問題在高維統(tǒng)計開始時便是一個非常重要的問題，而即使到了現(xiàn)在也是一個非常重要的問題，比如group lasso，generalized lasso，高斯圖模型，Tensor型圖模型，與圖相關(guān)的?1?1問題等算法的開發(fā)，都可以在此框架上直接應(yīng)用和實施，這正是ADMM一個優(yōu)勢所在，便于快速實施，也便于可能的大規(guī)模分布式部署。
minl(x)+λ∥x∥1,?minl(x)+g(z)=l(x)+λ∥z∥1s.t.x?z=0?xk+1=argminx(l(x)+(ρ/2)∥x?zk+uk∥22)zk+1=S1/ρ(xk+1+uk)uk+1=uk+xk+1?zk+1minl(x)+λ‖x‖1,?minl(x)+g(z)=l(x)+λ‖z‖1s.t.x?z=0?xk+1=arg?minx(l(x)+(ρ/2)‖x?zk+uk‖22)zk+1=S1/ρ(xk+1+uk)uk+1=uk+xk+1?zk+1
可以看到與Basis Pursuit解法只是在xx-update上有區(qū)別：Basis Pursuit是構(gòu)造出來一個投影函數(shù)f(x)f(x)，而一般化的損失函數(shù)f(x)f(x)+?1?1正則項問題，用ADMM就更為自然。所以很適合作為框架來解決這一類問題：廣義線性模型（普通線性、logistic回歸、possion回歸、softmax回歸）+正則項；廣義可加模型+正則項；似然函數(shù)（高斯圖方向）+正則項。
- Lasso：f(x)=12∥Ax?b∥22f(x)=12‖Ax?b‖22，于是利用ADMM算法，xx-update的解析解就是xk+1=(ATA+ρI)?1(ATb+ρ(zk?uk))xk+1=(ATA+ρI)?1(ATb+ρ(zk?uk))；于是xx-update看起來是個嶺回歸了，因此ADMM對于lasso可以看做迭代的使用嶺回歸。至于矩陣求逆那些，利用之前的矩陣小技巧解決。
- Generalized lasso：這個問題可能不是那么為眾人所熟悉，他是Tibs的兒子搞出來的框羅類似fused lasso這種事先定義好的線性變化的懲罰項的模型，損失函數(shù)是平方損失，而懲罰變成了一個特殊的參數(shù)線性組合
min12∥Ax?b∥22+λ∥Fx∥1min12‖Ax?b‖22+λ‖F(xiàn)x‖1?1d fused lasso,A=IFij=?????1j=i+1?1j=i0otherwise?1d fused lasso,A=IFij={1j=i+1?1j=i0otherwise
?min12∥x?b∥22+λn?1∑i=1|xi+1?xi|?A=I,F二階差分矩陣，則被稱作L1 trend filtering?min12‖x?b‖22+λ∑i=1n?1|xi+1?xi|?A=I,F二階差分矩陣，則被稱作L1 trend filtering
若將上述這種寫成ADMM形式，同樣可以放到ADMM算法框架中解決
min12∥Ax?b∥22+λ∥z∥1s.t.Fx?z=0?xk+1=(ATA+ρFTF)?1(ATb+ρFT(zk?uk))zk+1=S1/ρ(Axk+1?b+uk)uk+1=uk+Fxk+1?zk+1?bmin12‖Ax?b‖22+λ‖z‖1s.t.Fx?z=0?xk+1=(ATA+ρFTF)?1(ATb+ρFT(zk?uk))zk+1=S1/ρ(Axk+1?b+uk)uk+1=uk+Fxk+1?zk+1?b
- Group lasso：graph lasso問題應(yīng)用比較廣，對不同組的參數(shù)同時進行懲罰，進行一組組參數(shù)的挑選，故曰group lasso。不同于lasso，其正則項變成了∑Ni=1∥xi∥2,xi∈Rni∑i=1N‖xi‖2,xi∈Rni，lasso其實是group lasso的一種特殊形式。正則項并不是完全可分的。此時只是zz-update變成了block的軟閾值形式
zk+1i=Sλ/rho(xk+1i+uk),i=1,…,N?Sk(a)=(1?k∥a∥2)+a,S(0)=0zik+1=Sλ/rho(xik+1+uk),i=1,…,N?Sk(a)=(1?k‖a‖2)+a,S(0)=0
這種形式還可以擴展到group間有重合的情況，即化成NN可能存在重合的組Gi?{1,…,n}Gi?{1,…,n}。一般來說這種問題會非常難解決，但是對于ADMM算法只需要換下形式就很直接（x,zx,z互換，會變成后面非常重要的一致性優(yōu)化問題（consensus optimization），局部xixi與全局真解zz子集^ziz^i的對應(yīng)。）
min12∥Az?b∥22+λ∑Ni=1∥xi∥2,xi∈R|Gi|s.t.xi?^zi=0,i=1,…,Nmin12‖Az?b‖22+λ∑i=1N‖xi‖2,xi∈R|Gi|s.t.xi?z^i=0,i=1,…,N
- Sparse Gaussian graph model：對于稀疏高斯圖，熟悉該問題的人知道這其實是lasso的圖上的推廣，損失函數(shù)寫成似然函數(shù)的負(fù)數(shù)即可l(x)=tr(SX)?logdetX,X∈Sn++l(x)=tr(SX)?log?detX,X∈S++n。于是原來向量的操作就變成了矩陣操作，ADMM算法也有點變化：
Xk+1=argminX(tr(SX)?logdetX+ρ2∥X?Zk+Uk∥F)Zk+1=argminZ(λ∥Z∥1+ρ2∥Xk+1?Z+Uk∥F)Uk+1=Uk+Xk+1?Zk+1Xk+1=arg?minX(tr(SX)?log?detX+ρ2‖X?Zk+Uk‖F(xiàn))Zk+1=arg?minZ(λ‖Z‖1+ρ2‖Xk+1?Z+Uk‖F(xiàn))Uk+1=Uk+Xk+1?Zk+1
上述算法繼續(xù)化簡，對于zz-update做逐個元素軟閾值操作即可Zk+1ij=Sλ/ρ(XK+1ij+Ukij)Zijk+1=Sλ/ρ(XijK+1+Uijk)。對于xx-update也類似操作，直接求導(dǎo)一階導(dǎo)為0，移項后對對稱矩陣做特征值分解即可
ρX?X?1=ρ(Zk?Uk)?S=QΛQT,QQT=I,Λ=diag(λ1,…,λn)ρX?X?1=ρ(Zk?Uk)?S=QΛQT,QQT=I,Λ=diag(λ1,…,λn)→ρ^X?^X?1=Λ,^X=QTXQ→ρX^?X^?1=Λ,X^=QTXQ
由于ΛΛ是對角陣，對于每個對角元素來說，上述問題就是解一個二次方程，解方程后，再將^XX^變化成XX即可
^Xii=λi+√λ2i+4ρ2ρ?X=Q^XQTX^ii=λi+λi2+4ρ2ρ?X=QX^QT
總之，上述跟?1?1相關(guān)的問題，基本都可以納入ADMM框架，并且可以快速求解。
4. Consensus and Sharing
本節(jié)講述的兩個優(yōu)化問題，是非常常見的優(yōu)化問題，也非常重要，我認(rèn)為是ADMM算法通往并行和分布式計算的一個途徑：consensus和sharing，即一致性優(yōu)化問題與共享優(yōu)化問題。
Consensus
4.1 全局變量一致性優(yōu)化（Global variable consensus optimization）（切割數(shù)據(jù)，參數(shù)（變量）維數(shù)相同）
所謂全局變量一致性優(yōu)化問題，即目標(biāo)函數(shù)根據(jù)數(shù)據(jù)分解成NN子目標(biāo)函數(shù)（子系統(tǒng)），每個子系統(tǒng)和子數(shù)據(jù)都可以獲得一個參數(shù)解xixi，但是全局解只有一個zz，于是就可以寫成如下優(yōu)化命題：
min∑Ni=1fi(xi),xi∈Rns.t.xi?z=0min∑i=1Nfi(xi),xi∈Rns.t.xi?z=0
注意，此時fi:Rn→R?{+∞}fi:Rn→R?{+∞}仍是凸函數(shù)，而xixi并不是對參數(shù)空間進行劃分，這里是對數(shù)據(jù)而言，所以xixi維度一樣xi,z∈Rnxi,z∈Rn，與之前的問題并不太一樣。這種問題其實就是所謂的并行化處理，或分布式處理，希望從多個分塊的數(shù)據(jù)集中獲取相同的全局參數(shù)解。
在ADMM算法框架下（先返回最初從擴增lagrangian導(dǎo)出的ADMM），這種問題解法相當(dāng)明確：
Lρ(x1,…,xN,z,y)=∑Ni=1(fi(xi)+yTi(xi?z)+(ρ/2)∥xi?z∥22)s.t.C={(x1,…,xN)|x1=…=xN}Lρ(x1,…,xN,z,y)=∑i=1N(fi(xi)+yiT(xi?z)+(ρ/2)‖xi?z‖22)s.t.C={(x1,…,xN)|x1=…=xN}
?xk+1i=argminx(fi(xi)+(yki)T(xi?zk)+(ρ/2)∥xi?z∥22))zk+1=1NN∑i=1(xk+1i+(1ρyki))yk+1i=yki+ρ(xk+1i?zk+1)?xik+1=arg?minx(fi(xi)+(yik)T(xi?zk)+(ρ/2)‖xi?z‖22))zk+1=1N∑i=1N(xik+1+(1ρyik))yik+1=yik+ρ(xik+1?zk+1)
對yy-update和zz-update的yk+1iyik+1和zk+1izik+1分別求個平均，易得ˉyk+1=0yˉk+1=0，于是可以知道zz-update步其實可以簡化為zk+1=ˉxk+1zk+1=xˉk+1，于是上述ADMM其實可以進一步化簡為如下形式：
xk+1i=argminx(fi(xi)+(yki)T(xi?ˉxk)+(ρ/2)∥xi?ˉxk∥22))yk+1i=yki+ρ(xk+1i?ˉxk+1)xik+1=arg?minx(fi(xi)+(yik)T(xi?xˉk)+(ρ/2)‖xi?xˉk‖22))yik+1=yik+ρ(xik+1?xˉk+1)
這種迭代算法寫出來了，并行化那么就是輕而易舉了，各個子數(shù)據(jù)分別并行求最小化，然后將各個子數(shù)據(jù)的解匯集起來求均值，整體更新對偶變量ykyk，然后再繼續(xù)回帶求最小值至收斂。當(dāng)然也可以分布式部署（hadoop化），但是說起來容易，真正工程實施起來又是另外一回事，各個子節(jié)點機器間的通信更新是一個需要細(xì)細(xì)揣摩的問題。
另外，對于全局一致性優(yōu)化，也需要給出相應(yīng)的終止迭代準(zhǔn)則，與一般的ADMM類似，看primal和dual的residuals即可
∥rk∥22=N∑i=1∥xki?ˉxk∥22,∥sk∥22=Nρ∥ˉxki?ˉxk?1∥22‖rk‖22=∑i=1N‖xik?xˉk‖22,‖sk‖22=Nρ‖xˉik?xˉk?1‖22
4.2 帶正則項的全局一致性問題
下面就是要將之前所談到的經(jīng)典的機器學(xué)習(xí)算法并行化起來。想法很簡單，就是對全局變量加上正則項即可，因此ADMM算法只需要改變下zz-update步即可
min∑Ni=1fi(xi)+g(z),xi∈Rns.t.xi?z=0?xk+1i=argminx+i(fi(xi)+(yki)T(xi?zk)(ρ/2)∥xi?z∥22))zk+1=argminz(g(z)+N∑i=1(?(yki)Tz+(ρ/2)∥xk+1i?z∥22))yk+1i=yki+ρ(xk+1i?zk+1)min∑i=1Nfi(xi)+g(z),xi∈Rns.t.xi?z=0?xik+1=arg?minx+i(fi(xi)+(yik)T(xi?zk)(ρ/2)‖xi?z‖22))zk+1=arg?minz(g(z)+∑i=1N(?(yik)Tz+(ρ/2)‖xik+1?z‖22))yik+1=yik+ρ(xik+1?zk+1)
同樣的，我們?nèi)詫z做一個平均處理，于是就有
zk+1=argminz(g(z)+(Nρ/2)∥z?ˉxk+1?(1/ρ)ˉyk∥22)zk+1=arg?minz(g(z)+(Nρ/2)‖z?xˉk+1?(1/ρ)yˉk‖22)
上述形式都取得是最原始的ADMM形式，簡化處理，寫成scaled形式即有
xk+1i=argminx(fi(xi)+(ρ/2)∥xi?zk+uki∥22))zk+1=argminz(g(z)+(Nρ/2)∥z?xk+1i?ˉuk∥22)uk+1i=uki+xk+1i?zk+1xik+1=arg?minx(fi(xi)+(ρ/2)‖xi?zk+uik‖22))zk+1=arg?minz(g(z)+(Nρ/2)‖z?xik+1?uˉk‖22)uik+1=uik+xik+1?zk+1
這樣對于后續(xù)處理問題就清晰明了多了。可以看到如果g(z)=λ∥z∥1g(z)=λ‖z‖1，即lasso問題，那么zz-update步就用軟閾值operator即可。因此，對于大規(guī)模數(shù)據(jù)，要想用lasso等算法，只需要對數(shù)據(jù)做切塊（切塊也最好切均勻點），納入到全局變量一致性的ADMM框架中，即可并行化處理。下面給出一些實例。
切割大樣本數(shù)據(jù)，并行化計算
在經(jīng)典的統(tǒng)計估計中，我們處理的多半是大樣本低維度的數(shù)據(jù)，現(xiàn)在則多是是大樣本高維度的數(shù)據(jù)。對于經(jīng)典的大樣本低維度數(shù)據(jù)，如果機器不夠好，那么就抽樣部分?jǐn)?shù)據(jù)亦可以實現(xiàn)較好估計，不過如果沒有很好的信息，就是想要對大樣本進行處理，那么切割數(shù)據(jù)，并行計算是一個好的選擇。現(xiàn)在的社交網(wǎng)絡(luò)、網(wǎng)絡(luò)日志、無線感應(yīng)網(wǎng)絡(luò)等都可以這么實施。下面的具體模型都在受約束的凸優(yōu)化問題中以及?1?1-norm問題中提過，此處只不過切割數(shù)據(jù)，做成分布式模型，思想很簡單，與帶正則項的global consensus問題一樣的處理。經(jīng)典問題lasso、sparse logistic lasso、SVM都可以納入如下框架處理。
有觀測陣A∈Rm×nA∈Rm×n和響應(yīng)值b∈Rmb∈Rm，可以對應(yīng)切分，即對矩陣AA和向量bb橫著切，
A=?? ??A1?AN?? ??b=?? ??b1?bN?? ??A=(A1?AN)b=(b1?bN)
于是原來帶正則項的優(yōu)化問題就可以按照數(shù)據(jù)分解到多個子系統(tǒng)上去分別優(yōu)化，然后匯集起來，形成一個global consensus問題。
min∑Ni=1li(Aixi?bi)+r(z)s.t.xi?z=0,i=1,…,Nxi,z∈Rnmin∑i=1Nli(Aixi?bi)+r(z)s.t.xi?z=0,i=1,…,Nxi,z∈Rn
結(jié)合受約束的凸優(yōu)化問題時所給出來的具體的ADMM算法解的形式，下面直接給出這些問題的ADMM迭代算法公式
（1）Lasso
xk+1i=(ATiAi+ρI)?1(ATibi+ρ(zk?uki))zk+1=S1/ρN(ˉxk+1?b+ˉuk)uk+1i=uki+xk+1i?zk+1xik+1=(AiTAi+ρI)?1(AiTbi+ρ(zk?uik))zk+1=S1/ρN(xˉk+1?b+uˉk)uik+1=uik+xik+1?zk+1
如果切割的數(shù)據(jù)量小于維數(shù)mi<nmi<n，那么求解時分解小的矩陣AiATi+ρIAiAiT+ρI即可；其他求逆采用矩陣加速技巧即可。
（2）Sparse Logistic Regression
xk+1i=argminxi(li(Aixi?bi)+(ρ/2)∥xi?zk+uki∥22zk+1=S1/ρN(ˉxk+1?ˉb+ˉuk)uk+1i=uki+xk+1i?zk+1xik+1=arg?minxi(li(Aixi?bi)+(ρ/2)‖xi?zk+uik‖22zk+1=S1/ρN(xˉk+1?bˉ+uˉk)uik+1=uik+xik+1?zk+1
在xx-update步是需要用一些有效的算法來解決?2?2正則的logistic回歸，比如L-BFGS，其他的優(yōu)化算法應(yīng)該問題不大吧。
（3）SVM
注意分類問題和回歸問題的損失函數(shù)不同，一般都是用l(sign(t)y)l(sign(t)y)形式來尋求最優(yōu)的分類權(quán)重使得分類正確。SVM使用Hinge Loss：?(y)=max(0,1?t?y)?(y)=max(0,1?t?y)，即將預(yù)測類別與實際分類符號相反的損失給凸顯出來。分布式的ADMM形式
xk+1i=argminxi(1T(Aixi+1)++(ρ/2)∥xi?zk+uki∥22zk+1=ρ(1/λ)+Nρ(ˉxk+1+ˉuk)uk+1i=uki+xk+1i?zk+1xik+1=arg?minxi(1T(Aixi+1)++(ρ/2)‖xi?zk+uik‖22zk+1=ρ(1/λ)+Nρ(xˉk+1+uˉk)uik+1=uik+xik+1?zk+1
4.3 一般形式的一致性優(yōu)化問題（切割參數(shù)到各子系統(tǒng)，但各子系統(tǒng)目標(biāo)函數(shù)參數(shù)維度不同，可能部分重合）
上述全局一致性優(yōu)化問題，我們可以看到，所做的處理不過是對數(shù)據(jù)分塊，然后并行化處理。但是更一般的優(yōu)化問題是，參數(shù)空間也是分塊的，即每個子目標(biāo)函數(shù)fi(xi)fi(xi)的參數(shù)維度不同xi,∈Rnixi,∈Rni，我們稱之為局部變量。而局部變量所對應(yīng)的的也將不再是全局變量zz，而是全局變量中的一部分zgzg，并且不是像之前的順序?qū)?yīng)，而可能是隨便對應(yīng)到zz的某個位置。可令g=G(i,?)g=G(i,?)，即將xixi映射到zz的某部位
(xi)j=zG(i,j)=^zi(xi)j=zG(i,j)=z^i
如果對所有ii有G(i,j)=jG(i,j)=j，那么xixi與zz就是順序映射，也就是全局一致性優(yōu)化問題，否則就不是。結(jié)合下圖就比較好理解
consensus
雖然如果用其他方法來做感覺會復(fù)雜，但是納入到上述ADMM框架，其實只不過是全局一致性優(yōu)化問題的一個局部化變形，不過此時不是對數(shù)據(jù)進行分塊，是對參數(shù)空間進行分塊
min∑Ni=1fi(xi)+g(z),xi∈Rnis.t.xi?^zi=0,i=1,…N?xk+1i=argminx(fi(xi)+(yki)Txi(ρ/2)∥xi?^zki∥22))zk+1=argminz(N∑i=1(?(yki)T^zi+(ρ/2)∥xk+1i?^zi∥22)))yk+1i=yki+ρ(xk+1i?^zk+1i)min∑i=1Nfi(xi)+g(z),xi∈Rnis.t.xi?z^i=0,i=1,…N?xik+1=arg?minx(fi(xi)+(yik)Txi(ρ/2)‖xi?z^ik‖22))zk+1=arg?minz(∑i=1N(?(yik)Tz^i+(ρ/2)‖xik+1?z^i‖22)))yik+1=yik+ρ(xik+1?z^ik+1)
后續(xù)想做平均化處理，即中間會發(fā)生重合的參數(shù)zizi取值一樣的，那么zz-update將只能找他對應(yīng)的那些xx進行平均化，也就是變成局部了，因為不是所有值都是要全局保持一致的。比如上面那個圖中的z1,z2,z3,z4z1,z2,z3,z4都分別只要求在部分xixi發(fā)生了共享需要保持一樣，而不是像之前全局要求每個xixi對應(yīng)的都是zz。即
zk+1g=∑G(i,j)=g((xk+1i)j+(1/ρ)(yki)j)∑G(x,y)=g1zgk+1=∑G(i,j)=g((xik+1)j+(1/ρ)(yik)j)∑G(x,y)=g1
該式子表示就是zz的第gg個變量的平均值來源于所有映射到該變量的xx與yy的平均值。與之前的global類似，此時對于yy的取均值會為0，因此zz-update就變成了更簡單的形式
zk+1g=1kg∑G(i,j)=g(xk+1i)zgk+1=1kg∑G(i,j)=g(xik+1)
同全局一致性優(yōu)化問題一樣，我們可以加上正則項，然后也可以變成帶正則項的一般形式的一致性優(yōu)化問題。此處不贅述，與全局基本類似。
Sharing
4.4 共享問題（sharing）（橫向切割數(shù)據(jù)，也可縱向切變量）
與之前的全局變量一致性優(yōu)化問題類似，共享問題也是一個非常一般而且常見的問題。他的形式如下：
minN∑i=1fi(xi)+g(N∑i=1xi)min∑i=1Nfi(xi)+g(∑i=1Nxi)
這里的第一部分局部損失fi(xi)fi(xi)與全局一致性優(yōu)化是一樣的，即所有的xi∈Rn,i=1,…,Nxi∈Rn,i=1,…,N同維度，而對于一個共享的目標(biāo)函數(shù)gg則是新加入的。在實際中，我們常常需要優(yōu)化每個子數(shù)據(jù)集上的損失函數(shù)，同時還要加上全局?jǐn)?shù)據(jù)所帶來的損失；或者需要優(yōu)化每個子系統(tǒng)的部分變量，同時還要優(yōu)化整個變量。共享問題是一個非常重要而靈活的問題，它也可以納入到ADMM框架中，形式如下：
min∑Ni=1fi(xi)+g(∑Ni=1zi)s.t.xi?zi=0,zi∈Rn,i=1,…,N,?xk+1i=argminxi(fi(xi)+(ρ/2)∥xi?zki+uki∥22))zk+1=argminz(g(N∑i=1zi)+ρ/2N∑i=1∥zi?xk+1i?uki∥22)uk+1i=uki+xk+1i?zk+1imin∑i=1Nfi(xi)+g(∑i=1Nzi)s.t.xi?zi=0,zi∈Rn,i=1,…,N,?xik+1=arg?minxi(fi(xi)+(ρ/2)‖xi?zik+uik‖22))zk+1=arg?minz(g(∑i=1Nzi)+ρ/2∑i=1N‖zi?xik+1?uik‖22)uik+1=uik+xik+1?zik+1
上述形式當(dāng)然還不夠簡潔，需要進一步化簡。因為xx-update可以不用擔(dān)心，分機并行處理優(yōu)化求解即可，而對于zz-update這里面需要對NnNn個變量求解，想加快速度，就減少變量個數(shù)。于是想辦法通過和之前那種平均方式一樣來簡化形式解。
對于zz-update步，令ai=uki+xk+1iai=uik+xik+1，于是zz-update步優(yōu)化問題轉(zhuǎn)化為
ming(Nˉz)+(ρ/2)∑Ni=1∥zi?ai∥22s.t.ˉz=1N∑Ni=1ziming(Nzˉ)+(ρ/2)∑i=1N‖zi?ai‖22s.t.zˉ=1N∑i=1Nzi
當(dāng)ˉzzˉ固定時，那么后面的最優(yōu)解（類似回歸）為zi=ai+ˉz?ˉazi=ai+zˉ?aˉ，帶入上式后于是后續(xù)優(yōu)化就開始整體更新（均值化）
xk+1i=argminxi(fi(xi)+(ρ/2)∥xi?xki+ˉxk?ˉzk+uk∥22))zk+1=argminz(g(Nˉz)+Nρ/2∥ˉz?ˉxk+1?uk∥22)uk+1=uki+ˉxk+1?ˉzk+1xik+1=arg?minxi(fi(xi)+(ρ/2)‖xi?xik+xˉk?zˉk+uk‖22))zk+1=arg?minz(g(Nzˉ)+Nρ/2‖zˉ?xˉk+1?uk‖22)uk+1=uik+xˉk+1?zˉk+1
另外，有證明如果強對偶性存在，那么global consensus問題與sharing問題是可以相互轉(zhuǎn)化的，可以同時達(dá)到最優(yōu)，兩者存在著很緊密的對偶關(guān)系。
本節(jié)開頭提過，sharing問題用來切分?jǐn)?shù)據(jù)做并行化，也可以切分參數(shù)空間做并行化。這對于高維、超高維問題是非常有好處的。因為高維統(tǒng)計中，大樣本是一方面問題，而高維度才是重中之重，如果能切分特征到低緯度中去求解，然后在合并起來，那么這將是一個很美妙的事情。上面利用regularized global consensus問題解決了切分大樣本數(shù)據(jù)的并行化問題，下面利用sharing思想解決常見的高維數(shù)據(jù)并行化問題
切割變量（特征）空間，并行化處理
同樣假設(shè)面對還是一個觀測陣A∈Rm×nA∈Rm×n和響應(yīng)觀測b∈Rnb∈Rn，此時有n>>mn>>m，那么要么就降維處理，要么就切分維度去處理，或者對于超高維矩陣，切分維度后再降維。此時AA矩陣就不是像之前橫著切分，而是豎著切分，這樣對應(yīng)著參數(shù)空間的切分：
A=[A1,…,AN],Ai∈Rm×ni,x=(x1,…,xN),x∈Rni,→Ax=N∑i=1AixiA=[A1,…,AN],Ai∈Rm×ni,x=(x1,…,xN),x∈Rni,→Ax=∑i=1NAixi
于是正則項也可以切分為r(x)=∑Ni=1ri(xi)r(x)=∑i=1Nri(xi)。那么最初的minl(Ax?b)+r(x)minl(Ax?b)+r(x)形式就變成了
minl(N∑i=1Aixi?b)+N∑i=1ri(xi)minl(∑i=1NAixi?b)+∑i=1Nri(xi)
這個與sharing問題非常接近了，做點變化那就是sharing問題了
minl(∑Ni=1zi?b)+∑Ni=1ri(xi)s.t.Aixi?zi=0,i=1,…,N?xk+1i=argminxi(ri(xi)+(ρ/2)∥Aixi?Aixki+ˉˉˉˉˉˉˉAxk?ˉzk+uk∥22))zk+1=argminz(l(Nˉz?b)+Nρ/2∥ˉz?ˉˉˉˉˉˉˉAxk+1?uk∥22)uk+1=uki+ˉˉˉˉˉˉˉAxk+1?ˉzk+1minl(∑i=1Nzi?b)+∑i=1Nri(xi)s.t.Aixi?zi=0,i=1,…,N?xik+1=arg?minxi(ri(xi)+(ρ/2)‖Aixi?Aixik+Axˉk?zˉk+uk‖22))zk+1=arg?minz(l(Nzˉ?b)+Nρ/2‖zˉ?Axˉk+1?uk‖22)uk+1=uik+Axˉk+1?zˉk+1
與之前的global consensus問題相比，ADMM框架xx-update與zz-update似乎是反過來了。于是將此形式直接套到Lasso等高維問題即有很具體的形式解了。
（1）Lasso
xk+1i=argminxi(λ∥xi∥1+(ρ/2)∥Aixi?Aixki+ˉˉˉˉˉˉˉAxk?ˉzk+uk∥22))ˉzk+1=1N+ρ(b+ρˉˉˉˉˉˉˉAxk+1+ρuk)uk+1=uk+ˉˉˉˉˉˉˉAxk+1?ˉzk+1xik+1=arg?minxi(λ‖xi‖1+(ρ/2)‖Aixi?Aixik+Axˉk?zˉk+uk‖22))zˉk+1=1N+ρ(b+ρAxˉk+1+ρuk)uk+1=uk+Axˉk+1?zˉk+1
當(dāng)∥ATi(Aixki+ˉzk?ˉˉˉˉˉˉˉAxk?uk)∥2≤λ/ρ‖AiT(Aixik+zˉk?Axˉk?uk)‖2≤λ/ρ時xk+1i=0xik+1=0（第ii塊特征不需要用），這樣加快了xx-update速度,不過這個對串行更有效，對并行起來也沒有多大用..
（2）Group Lasso?與lasso基本一樣，只是在xx-update上有一個正則項的不同，有?1?1-norm變成了?2?2-norm
xk+1i=argminxi(λ∥xi∥2+(ρ/2)∥Aixi?Aixki+ˉˉˉˉˉˉˉAxk?ˉzk+uk∥22)xik+1=arg?minxi(λ‖xi‖2+(ρ/2)‖Aixi?Aixik+Axˉk?zˉk+uk‖22)
該問題其實就是按組最小化(ρ/2)∥Aixi?v∥22+λ∥xi∥2(ρ/2)‖Aixi?v‖22+λ‖xi‖2，解為
if∥ATiv∥2≤λ/ρ,thenxi=0otherwisexi=(ATiAi+vI)?1ATivif‖AiTv‖2≤λ/ρ,thenxi=0otherwisexi=(AiTAi+vI)?1AiTv
涉及矩陣長短計算時，再看矩陣小技巧。
（3）Sparse Logstic Regression?也與lasso區(qū)別不大，只是zz-update的損失函數(shù)不同，其余相同于是
ˉzk+1=argminˉz(l(Nˉz)+(ρ/2)∥ˉz?ˉˉˉˉˉˉˉAxk+1?uk∥22)zˉk+1=arg?minzˉ(l(Nzˉ)+(ρ/2)‖zˉ?Axˉk+1?uk‖22)
（4）SVM
SVM與之前的global consensus時候優(yōu)化順序反了過來，與logistic rgression只是在zz-update步不同（損失函數(shù)不同）：
xk+1i=argminxi(λ∥xi∥22+(ρ/2)∥Aixi?Aixki+ˉˉˉˉˉˉˉAxk?ˉzk+uk∥22))ˉzk+1=argminˉz(1T(Nˉz+1)++(ρ/2)∥ˉz?ˉˉˉˉˉˉˉAxk+1?uk+1∥)uk+1=uk+ˉˉˉˉˉˉˉAxk+1?ˉzk+1xik+1=arg?minxi(λ‖xi‖22+(ρ/2)‖Aixi?Aixik+Axˉk?zˉk+uk‖22))zˉk+1=arg?minzˉ(1T(Nzˉ+1)++(ρ/2)‖zˉ?Axˉk+1?uk+1‖)uk+1=uk+Axˉk+1?zˉk+1
zz-update解析解可以寫成軟閾值算子
(ˉzk+1)i=?????vi?N/ρ,vi>?1/N+N/ρ?1/N,vi∈[?1/N,?1/N+N/ρ]vi,vi<?1/Nvi=(ˉˉˉˉˉˉˉAxk+1+ˉuk)i(zˉk+1)i={vi?N/ρ,vi>?1/N+N/ρ?1/N,vi∈[?1/N,?1/N+N/ρ]vi,vi<?1/Nvi=(Axˉk+1+uˉk)i
（5）Generalized Additive Models
廣義可加模型是一個很適合sharing框架的問題。它本身就是對各個各個特征做了變化后（非參方法），重新表示觀測的方式
b≈n∑j=1fj(xj)b≈∑j=1nfj(xj)
當(dāng)fifi是線性變化時，則退化成普通線性回歸。此時我們目標(biāo)優(yōu)化的問題是
minm∑i=1li(n∑j=1fj(xij)?bi)+n∑j=1rj(fj)min∑i=1mli(∑j=1nfj(xij)?bi)+∑j=1nrj(fj)
其中有mm個觀測，nn維特征（變量）。rjrj此時是對一個functional的正則，此時這個問題看起來似乎既可以對數(shù)據(jù)切分，也可以對特征切分，不過此時仍用sharing問題來做，相當(dāng)于對特征切分為一個特征為一個子系統(tǒng)，于是有
$$ \begin{split} f_j^{k+1} & =?{f_i?j} (r_j(f_j)+ (/2)^m{i=1}(f_j(x{ij}) - f^k_j(x_{ij}) + {z}^k_i + {f}^k_i) + uk_i\ {z}{k+1} & =?{{z}}(^m{i=1}l_i(N{z} - b_i) + /2^n_{j=1}|{z} - {f}^{k+1} - u^k| ,, ,{f}^k =?n_{j=1}fk_j(x_{ij})\ u^{k+1} & = u^k + {f}^{k+1} - {z}^{k+1} \ \end{split}
$$
fjfj是一個?2?2正則的損失，有直接求解的算法求解，zz可以一塊一塊的求解？
最后再說一個經(jīng)濟學(xué)中很重要的sharing問題的特例，即交換問題（exchange problem）：
min∑Ni=1fi(xi)s.t.∑Ni=1xi=0,xi∈Rn,i=1,…Nmin∑i=1Nfi(xi)s.t.∑i=1Nxi=0,xi∈Rn,i=1,…N
此時共享目標(biāo)函數(shù)g=0g=0。xixi可以表示不同物品在NN個系統(tǒng)上的交換數(shù)量，(xi)j(xi)j可以表示物品jj從子系統(tǒng)ii上收到的交換數(shù)目，約束條件就可以看做在這些系統(tǒng)中物品交換是保持均衡穩(wěn)定的。于是轉(zhuǎn)化為sharing問題，就有很簡單的ADMM解法（或者當(dāng)做之前講過的受約束的凸優(yōu)化問題來解，做投影）：
xk+1i=argminxi(fi(xi)+(ρ/2)∥xi?xki+ˉxk+uk∥22))uk+1=uki+ˉxk+1xik+1=arg?minxi(fi(xi)+(ρ/2)‖xi?xik+xˉk+uk‖22))uk+1=uik+xˉk+1
4.4 應(yīng)用小總結(jié)
感覺上通過consensus problem和general consensus problem，我們可以看到并行和分布式部署優(yōu)化方案的可行性。我們可以切分?jǐn)?shù)據(jù)以及相應(yīng)的目標(biāo)函數(shù)，也可以切分變量到各個子系統(tǒng)上去，分別作優(yōu)化，甚至我們可以大膽想象對不同類型數(shù)據(jù)塊用不同的優(yōu)化算法，結(jié)合consensus問題和ADMM算法，達(dá)到同一個global variable的優(yōu)化目的；或者對不同變量在不同類型數(shù)據(jù)塊上優(yōu)化，即使有重疊，也可以結(jié)合general consensus思想和ADMM算法來解決這個問題。當(dāng)然前提是能夠定義好需要估計的參數(shù)和優(yōu)化的目標(biāo)函數(shù)！大規(guī)模部署的前景還是很不錯的。下面具體分布式統(tǒng)計模型的構(gòu)建便是ADMM算法非常好的應(yīng)用。切分?jǐn)?shù)據(jù)、切分變量（不過每個子系統(tǒng)的目標(biāo)函數(shù)基本都是一樣的，其實應(yīng)該可以不同）
5. Nonconvex問題
5.1 變量選擇（Regressor Selection）
5.2 因子模型（Factor Model Fitting）
5.3 雙凸優(yōu)化（Bi-convex Problem）
非負(fù)矩陣分解（Nonnegative Matrix Factorization）
6. 具體實施與實際計算結(jié)果
這塊真的很實際，需要明白MPI的機理和Mapreduce、Graphlab等通信運作的機理，這樣才好部署ADMM算法，因為中間有很多迭代，需要做好子節(jié)點間參數(shù)與整體參數(shù)的通信，保持迭代時能同步更新參數(shù)。看實際運作，MPI和GraphLab可能更適合這種框架，Hadoop也是可以的，不過畢竟不是為迭代算法所生，要做好需要進行一些優(yōu)化。Boyd提到Hadoop其中的Hbase更適合這種框架，因為Hbase是一種大表格，帶有時間戳，適合記錄迭代的記錄，這樣就不容易導(dǎo)致分布計算時候搞不清是哪一步的迭代結(jié)果了，導(dǎo)致通信調(diào)整比較復(fù)雜。不過在MapReduce框架下實施ADMM算法是沒有什么問題的，只要熟稔Hadoop的一些細(xì)節(jié)部分，基本沒有太大問題。
8. 總結(jié)
一個好的一般性算法，我個人覺得是易實施，并可大規(guī)模應(yīng)用許多問題。可以讓統(tǒng)計學(xué)家卡在搞算法的瓶頸中解放出來，使得他們能快速用模擬，驗證自己構(gòu)建可能較為復(fù)雜的模型。只有當(dāng)看到一個令人感到欣慰的結(jié)果時，那些模型的統(tǒng)計性質(zhì)的證明才可能是有意義的，如果事先連希望都看不到，那證明起來都可能底氣不足，讓人難以信服，更難以大規(guī)模應(yīng)用統(tǒng)計學(xué)家所構(gòu)建的模型。現(xiàn)在是一個高維數(shù)據(jù)、海量數(shù)據(jù)的年代，算法的重要性更會凸顯出來，一個好的模型如果沒有一個有效的算法支撐，那么他將可能什么都不是，Lasso頭幾年所遭遇的冷遇也充分證明了這一點，再比如在沒有計算機年代，Pearson的矩估計應(yīng)用反而遠(yuǎn)多于Fisher的MLE估計方法也是一個道理。好的一般性的解決方案，我想這不管是優(yōu)化理論，還是統(tǒng)計等其他應(yīng)用學(xué)科，雖然知道沒有最牛最終極的方法，但是能涌現(xiàn)一些大范圍適用的方法，那就是再好不過了。一招鮮吃遍天，人還都是喜歡簡單、安逸愛偷懶的嘛..
總結(jié)
- 上一篇： java索引ref_java – 如何
- 下一篇： Comet OJ - Contest #