ADMM算法
https://joegaotao.github.io/2014/02/11/admm-stat-compute/
分布式計算、統計學習與ADMM算法
?
在整理舊電腦時，才發現13年下半年電腦里有不少殘文。老師說，東西擱下了再拿起來花費的時間和之前可能差不多。我一眼看過去這篇關于分布式計算的文章，貌似還真的沒有了當時理解的深度和感覺。當時還想利用ADMM算法，把統計中常見的帶懲罰的高維問題在此框架下用R重寫一下，但是中途多種事情一耽擱，就早已拋之腦后。看來任何事情，真的還是需要堅持，哪怕撥點時間都是好的。先把一篇殘文扔出來祭奠下過去的13年吧。公式多文字長，慎入！
業界一直在談論大數據，對于統計而言，大數據其實意味著要不是樣本量增加n→∞n→∞，要不就是維度的增加p→∞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。這種轉化可以將原來很多的參數約束條件變得少了很多，以利于做優化。具體表述如下：
minf(x)s.t.Ax=b?L(x,y)=f(x)+yT(Ax?b)對偶函數（下界）?g(y)=infxL(x,y)minf(x)s.t.Ax=b?L(x,y)=f(x)+yT(Ax?b)?對偶函數（下界）g(y)=infxL(x,y)
在強對偶性的假設下，primal和dual問題同時達到最優。
x?=argminxL(x,y?)x?=arg?minxL(x,y?)
因此，若對偶函數g(y)g(y)可導，便可以利用梯度上升法，交替更新參數，使得同時收斂到最優。迭代如下：
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是步長)
當gg不可微的時候也可以將其轉化下，成為一個所謂的subgradient的方法，雖然看起來不錯，簡單證明下即可知道xkxk和ykyk同時可達到optimal，但是上述條件要求很苛刻：f(x)f(x)要求嚴格凸，并且要求αα選擇有比較合適。一般應用中都不會滿足（比如f(x)f(x)是一個非零的仿射函數），因此dual ascent不會直接應用。
1.2 Dual Decomposition
雖然dual ascent方法有缺陷，要求有些嚴格，但是他有一個非常好的性質，當目標函數ff是可分的（separable）時候（參數抑或feature可分），整個問題可以拆解成多個子參數問題，分塊優化后匯集起來整體更新。這樣非常有利于并行化處理。形式化闡述如下：
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)
因此可以看到其實下面在迭代優化時，xx-minimization步即可以拆分為多個子問題的并行優化，對偶變量更新不變這對于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，然后對偶變量更新)
對偶分解是非常經典的優化方法，可追溯到1960年代。但是這種想法對后面的分布式優化方法影響較大，比如近期的graph-structure優化問題。
1.3 Augmented Lagrangians and the Method of Multipliers
從上面可以看到dual ascent方法對于目標函數要求比較苛刻，為了放松假設條件，同時比較好優化，于是就有了Augmented Lagrangians方法，目的就是放松對于f(x)f(x)嚴格凸的假設和其他一些條件，同時還能使得算法更加穩健。
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
從上面可以看到該問題等價于最初的問題，因為只要是可行解對目標函數就沒有影響。但是加了后面的(ρ/2)∥Ax?b∥22(ρ/2)‖Ax?b‖22懲罰項的好處是使得對偶函數gρ(y)=infxLρ(x,y)gρ(y)=infxLρ(x,y)在更一般的條件下可導。計算過程與之前的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轉為固定的ρρ了吧。該算法在即使f(x)f(x)不是嚴格凸或者取值為+∞+∞情況都可以成立，適用面更廣。同樣可以簡單證明primal變量xx和對偶變量yy可以同時達到最優。
雖然Augmented Lagrangians方法有優勢，但也破壞了dual ascent方法的利用分解參數來并行的優勢。當ff是separable時，對于Augmented Lagrangians卻是not separable的（因為平方項寫成矩陣形式無法用之前那種分塊形式），因此在x?minx?min步時候無法并行優化多個參數xixi。如何改進，繼續下面的議題就可以慢慢發現改進思想的來源。
2. Alternating Direction Method of Multipliers(ADMM)
2.1 ADMM算法概述
為了整合dual ascent可分解性與method multiplers優秀的收斂性質，人們就又提出了改進形式的優化ADMM。目的就是想能分解原函數和擴增函數，以便于在對ff更一般的假設條件下并行優化。ADMM從名字可以看到是在原來Method of Multipliers加了個Alternating Direction，可以大概猜想到應該是又想引入新變量，然后交叉換方向來交替優化。形式如下：
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變量、目標函數拆分，但是不再像dual ascent方法那樣，將拆分開的xixi都看做是xx的一部分，后面融合的時候還需要融合在一起，而是最先開始就將拆開的變量分別看做是不同的變量xx和zz，同時約束條件也如此處理，這樣的好處就是后面不需要一起融合xx和zz，保證了前面優化過程的可分解性。于是ADMM的優化就變成了如下序貫型迭代（這正是被稱作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)
后面我們可以看到這種拆分思想非常適合統計學習中的?1?1-norm等問題：loss + regulazition（注意：一定要保證zz分解出來，ADMM借助的就是用一個zz變量來簡化問題，不管他是約束還是其他形式也罷，需要構造一個zz出來，后面具體到細節問題我們會有更深的體會）。
為了簡化形式，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
寫成這種形式有利于后面簡化優化問題，當然可以不作任何處理。
2.2 ADMM算法性質和評價
（1）收斂性
關于收斂性，需要有兩個假設條件：
- ff和gg分別是擴展的實數函數Rn(Rm)→R?{+∞}Rn(Rm)→R?{+∞}，且是closed、proper和convex的；
- 擴增的lagrangian函數L0L0有一個鞍點（saddle point）；對于約束中的矩陣A,BA,B都不需要滿秩。
在此兩個假設下，可以保證殘差、目標函數、對偶變量的收斂性。
Note：實際應用而言，ADMM收斂速度是很慢的，類似于共軛梯度方法。迭代數十次后只可以得到一個acceptable的結果，與快速的高精度算法（Newton法，內點法等）相比收斂就慢很多了。因此實際應用的時候，其實會將ADMM與其他高精度算法結合起來，這樣從一個acceptable的結果變得在預期時間內可以達到較高收斂精度。不過一般在大規模應用問題中，高精度的參數解對于預測效果沒有很大的提高，因此實際應用中，短時間內一個acceptable的結果基本就可以直接應用預測了。
（2）停止準則
對于ADMM的能到到optimal的條件此處就不做贅述了，與基本的primal和dual feasibility 的條件差不多，即各primal variable的偏導和約束條件為0，從最優條件中可以得到所謂的對偶殘差（dual residuals）和初始殘差（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)
相對而言，此處更難把握的其實是停止準則，因為收斂速度問題，要想獲得一個還過得去可以拿來用的參數解，那么判斷迭代停止還是比較重要的。實際應用中，一般都根據primal residuals和dual residuals足夠小來停止迭代，閾值包含了絕對容忍度（absolute tolerance）和相對容忍度（relative tolerance），設置還是非常靈活和難把握的（貌似網上有不少人吐槽這個停止準則的不靠譜- -！），具體形式如下：
∥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，絕對閾值的選取要根據變量取值范圍來選取（咋選的呢？沒說額，具體比例都不給說- -！）
另外一些細節問題，比如原來懲罰參數ρρ是不變的，一些文獻也做了一些可變的懲罰參數，目的是為了降低對于懲罰參數初始值的依賴性。不過變動的ρρ會導致ADMM的收斂性證明比較困難，因此實際中假設經過一系列迭代后ρρ也穩定，邊可直接用固定的懲罰參數ρρ了。還有其他問題，諸如xx與zz迭代順序問題，實際操作下有所有不同，這些不是特別重要之處，可以忽略。其他與ADMM比較相關算法的有dual ADMM算法，distributed ADMM算法，還有整合了ADMM與proximal method of multiplier的算法
2.3 ADMM一般形式與部分具體應用
當構造了ADMM算法中的f,g,A,Bf,g,A,B后，便可直接應用該算法了。我們會經常遇到如下三種一般形式的問題
- 二次目標優化項（quadratic objective terms）；
- 可分的目標函數和約束（separable objective and constraints）；
- 光滑目標函數項（smooth objective terms）。
為下面討論的方便，下面僅寫出xx-update的形式，根據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都定下來，是個常數，zz更新時后相同。
Proximity Operator（近鄰算子）
上述形式有種特殊情況：當A=IA=I時，即約束條件沒有xx的線性組合形式，只是對于xx的可行區域進行限制。這種問題相當常見，目前統計學習也有不少類似的高維優化問題。此時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的函數proxf,ρ(v)proxf,ρ(v)被稱作帶懲罰ρρ的ff的proximity operator（通常稱作proximal minimization，近鄰最小化），在變分分析中，還被稱作ff的Moreau-Yosida正則化。如果ff形式很簡單，可以寫出xx-update的解析解，比如ff是非空的凸包CC上的示性函數，那么xx-update就可以直接寫成投影形式
x+=argminx(f(x)+(ρ/2)∥x?v∥22)=ΠC(v)x+=arg?minx(f(x)+(ρ/2)‖x?v‖22)=ΠC(v)
投影與懲罰參數ρρ無關。若ff是非負象限Rn+R+n的投影，則直接有x+=(v)+x+=(v)+。
下面再談談上述提到的三種一般形式的優化問題。
（1）Quadratic Objective Terms
假設ff是如下（凸）的二次函數
f(x)=12xTPx+qTx+rf(x)=12xTPx+qTx+r
PP是對稱的半正定矩陣P∈Sn+P∈S+n。這種形式問題也包含了ff是線性或者常數的特殊情況。若P+ρATAP+ρATA可逆，那么xx-update步求個導即有如下的顯示解，是vv的仿射函數
x+=(P+ρATA)?1(ρATv?q)x+=(P+ρATA)?1(ρATv?q)
因此在xx-minnimiztion步只需要做兩個矩陣運算即可，求逆與乘積，選用合適的線性運算庫即可以得到不錯的計算性能。當然還可以利用一些矩陣分解技巧，這個要看矩陣大小和稀疏程度。因為對于Fx=gFx=g，可以將F=F1F2?FkF=F1F2?Fk，然后Fizi=zi?1,z1=F?11g,x=zkFizi=zi?1,z1=F1?1g,x=zk，這樣會更節省計算時間。其他矩陣計算技巧，基本都是如何對矩陣大規模求解，利用矩陣的稀疏性、緩存分解等來提高性能。此處不贅述，有個很重要的求逆的定理很有用：
(P+ρATA)?1=P?1?ρP?1AT(I+ρAP?1AT)?1AP?1(P+ρATA)?1=P?1?ρP?1AT(I+ρAP?1AT)?1AP?1
如果對于上述二次函數受限于某仿射集xx-update步就更復雜些，如
f(x)=12xTPx+qTx+rdormf={x∥Fx=g}f(x)=12xTPx+qTx+rdormf={x‖Fx=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
當ff光滑時，那么求導即成為可能了。對于一些非線性優化問題，包含梯度算法等方法的L-BFGS算法可以用。對于該算法有些小技巧如下：
- 早終止（early termination）：當f(x)+(ρ/2)∥Ax?v∥22f(x)+(ρ/2)‖Ax?v‖22梯度很小時，早點終止迭代，否則后面就很慢了。
- 熱啟動（warm start）：即啟動迭代時，利用之前迭代過的值帶入即可。
（3）Separable objective and constraints?可分函數和約束對于并行計算和分布式計算來說是一個好消息。如果ATAATA是分塊的對角陣，那么約束中∥Ax∥22‖Ax‖22也是可分的，則擴增的拉格朗日函數LρLρ也是可分的。（注意，此處是指函數中的參數可分成小子塊，而不是說數據可分。）下面有一個很重要的例子，即soft thresholding（針對l1+l2l1+l2問題）:
當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)
這種形式很常見在目前的高維統計中，雖然第一項在0處不可導，但是也有解析解，被稱作軟閾值（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
在優化領域，軟閾值被稱作是?1?1-norm問題的近鄰算子（proximity operator）。
3. 一些具體優化應用
3.1受約束的凸優化問題
一般的受約束的凸優化問題可以寫成如下形式
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函數即CC的示性函數，上述是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就變成了一個具體的受約束的優化問題。比如對于經典的二次規劃問題(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))
即受約束的區域就是{x∣x≥0}{x∣x≥0}，gg是向非負象限投影的示性函數。而xx-update就變成了之前在Quadratic Objective Terms中談到的f(x)f(x)有仿射集定義域的優化問題，根據KKT條件即可寫出來xx-update更新的形式，參見2.3節。
如果上述對xx限制不是限制x≥0x≥0上，而是一個錐約束（conic constraint）x∈Kx∈K，那么xx-update不變，繼續上述KKT方程，而只需要變一下zz-update，將向Rn+R+n投影改成向KK投影。比如將上述約束改成{Ax=b,x∈§n+}{Ax=b,x∈§+n}，即xx屬于半正定空間，那么向Sn+S+n投影就變成了一個半正定問題，利用特征值分解可以完成。這種受約束的凸優化問題的形式化對后續許多問題，特別是我們很關注的?1?1-norm問題很重要，基本上都是轉化成這種形式來直接應用ADMM算法，所以這里要好好把握其核心思想和形式。
雖然我對優化不在行，但是感覺優化問題還是挺有意思的，下面是一個經典問題，即找到兩個非空凸包的交集中的一點。該算法都可以追溯到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問題
高維統計理論的發展，如果要追溯起來我覺得可以從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?1-norm的解決方案，那么這些算法中ADMM可能是首選。此處?1?1-norm問題并不僅僅指Lasso問題，包含了多種?1?1-norm類型問題。下面均介紹下。
之所以說ADMM適合機器學習和統計學習的優化問題，因為大部分機器學習問題基本都是“損失函數+正則項”形式，這種分法恰好可以套用到ADMM的框架f(x)+g(z)f(x)+g(z)。因此結合ADMM框架基本可以解決很多已有的問題，以及利用?1?1-norm構造的新的優化問題。下面將先介紹非分布式計算的版本，后面會單開一節來介紹如何分布式計算。
（1）Least Absolute Deviations
先從一個簡單的問題開始。在穩健估計中，LAD是一個應用很廣的模型，相對于直接優化平方和損失∥Ax?b∥22‖Ax?b‖22，優化絕對損失∥Ax?b∥1‖Ax?b‖1，它的抗噪性能更好。在ADMM框架下，往之前的受約束的凸優化問題靠攏，這個問題有簡單的迭代算法
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問題與上面的其實差不多，只是損失函數形式不同，換成了Huber懲罰函數
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與一個殘差的加權。
LAD和Huber fitting這種問題只是一些傳統損失不加正則項的ADMM化，注意一定要構造個zz出來即可，xx可以基本不用管，總是需要解的，下面的帶有正則項的優化問題，ADMM形式就會更明顯。
（3）Basis Pursuit
基追蹤法師系數信號處理的一種重要方法。目的是想找到一組稀疏基可以完美恢復信號，換套話說就是為一個線性方程系統找到一個稀疏解。原始形式如下，與lasso有些像：
min∥x∥1s.t.Ax=bmin‖x‖1s.t.Ax=b
修改成ADMM形式，注意往之前受約束的凸優化問題的那種形式回套，將?1?1看做約束，然后構造帶定義域的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
對于矩陣求逆、分解等用之前矩陣那些小技巧即可加快計算，節省計算資源。
最近還有一類算法來解決?1?1問題，被稱作Bregman iteration methods，對于基追蹤相關問題，加正則項的Bregman iteration就是method of multiplier，而所謂的split Bregman iteration就等同于 ADMM。我沒有繼續深究，應該就是類似于并行化的ADMM算法來解決基追蹤問題。
（4）一般化的損失函數 +??1?1正則項問題
這類問題在高維統計開始時便是一個非常重要的問題，而即使到了現在也是一個非常重要的問題，比如group lasso，generalized lasso，高斯圖模型，Tensor型圖模型，與圖相關的?1?1問題等算法的開發，都可以在此框架上直接應用和實施，這正是ADMM一個優勢所在，便于快速實施，也便于可能的大規模分布式部署。
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上有區別：Basis Pursuit是構造出來一個投影函數f(x)f(x)，而一般化的損失函數f(x)f(x)+?1?1正則項問題，用ADMM就更為自然。所以很適合作為框架來解決這一類問題：廣義線性模型（普通線性、logistic回歸、possion回歸、softmax回歸）+正則項；廣義可加模型+正則項；似然函數（高斯圖方向）+正則項。
- 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這種事先定義好的線性變化的懲罰項的模型，損失函數是平方損失，而懲罰變成了一個特殊的參數線性組合
min12∥Ax?b∥22+λ∥Fx∥1min12‖Ax?b‖22+λ‖Fx‖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問題應用比較廣，對不同組的參數同時進行懲罰，進行一組組參數的挑選，故曰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互換，會變成后面非常重要的一致性優化問題（consensus optimization），局部xixi與全局真解zz子集^ziz^i的對應。）
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的圖上的推廣，損失函數寫成似然函數的負數即可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)Zk+1=arg?minZ(λ‖Z‖1+ρ2‖Xk+1?Z+Uk‖F)Uk+1=Uk+Xk+1?Zk+1
上述算法繼續化簡，對于zz-update做逐個元素軟閾值操作即可Zk+1ij=Sλ/ρ(XK+1ij+Ukij)Zijk+1=Sλ/ρ(XijK+1+Uijk)。對于xx-update也類似操作，直接求導一階導為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相關的問題，基本都可以納入ADMM框架，并且可以快速求解。
4. Consensus and Sharing
本節講述的兩個優化問題，是非常常見的優化問題，也非常重要，我認為是ADMM算法通往并行和分布式計算的一個途徑：consensus和sharing，即一致性優化問題與共享優化問題。
Consensus
4.1 全局變量一致性優化（Global variable consensus optimization）（切割數據，參數（變量）維數相同）
所謂全局變量一致性優化問題，即目標函數根據數據分解成NN子目標函數（子系統），每個子系統和子數據都可以獲得一個參數解xixi，但是全局解只有一個zz，于是就可以寫成如下優化命題：
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?{+∞}仍是凸函數，而xixi并不是對參數空間進行劃分，這里是對數據而言，所以xixi維度一樣xi,z∈Rnxi,z∈Rn，與之前的問題并不太一樣。這種問題其實就是所謂的并行化處理，或分布式處理，希望從多個分塊的數據集中獲取相同的全局參數解。
在ADMM算法框架下（先返回最初從擴增lagrangian導出的ADMM），這種問題解法相當明確：
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)
這種迭代算法寫出來了，并行化那么就是輕而易舉了，各個子數據分別并行求最小化，然后將各個子數據的解匯集起來求均值，整體更新對偶變量ykyk，然后再繼續回帶求最小值至收斂。當然也可以分布式部署（hadoop化），但是說起來容易，真正工程實施起來又是另外一回事，各個子節點機器間的通信更新是一個需要細細揣摩的問題。
另外，對于全局一致性優化，也需要給出相應的終止迭代準則，與一般的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 帶正則項的全局一致性問題
下面就是要將之前所談到的經典的機器學習算法并行化起來。想法很簡單，就是對全局變量加上正則項即可，因此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)
同樣的，我們仍對zz做一個平均處理，于是就有
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
這樣對于后續處理問題就清晰明了多了。可以看到如果g(z)=λ∥z∥1g(z)=λ‖z‖1，即lasso問題，那么zz-update步就用軟閾值operator即可。因此，對于大規模數據，要想用lasso等算法，只需要對數據做切塊（切塊也最好切均勻點），納入到全局變量一致性的ADMM框架中，即可并行化處理。下面給出一些實例。
切割大樣本數據，并行化計算
在經典的統計估計中，我們處理的多半是大樣本低維度的數據，現在則多是是大樣本高維度的數據。對于經典的大樣本低維度數據，如果機器不夠好，那么就抽樣部分數據亦可以實現較好估計，不過如果沒有很好的信息，就是想要對大樣本進行處理，那么切割數據，并行計算是一個好的選擇。現在的社交網絡、網絡日志、無線感應網絡等都可以這么實施。下面的具體模型都在受約束的凸優化問題中以及?1?1-norm問題中提過，此處只不過切割數據，做成分布式模型，思想很簡單，與帶正則項的global consensus問題一樣的處理。經典問題lasso、sparse logistic lasso、SVM都可以納入如下框架處理。
有觀測陣A∈Rm×nA∈Rm×n和響應值b∈Rmb∈Rm，可以對應切分，即對矩陣AA和向量bb橫著切，
A=?? ??A1?AN?? ??b=?? ??b1?bN?? ??A=(A1?AN)b=(b1?bN)
于是原來帶正則項的優化問題就可以按照數據分解到多個子系統上去分別優化，然后匯集起來，形成一個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
結合受約束的凸優化問題時所給出來的具體的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
如果切割的數據量小于維數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，其他的優化算法應該問題不大吧。
（3）SVM
注意分類問題和回歸問題的損失函數不同，一般都是用l(sign(t)y)l(sign(t)y)形式來尋求最優的分類權重使得分類正確。SVM使用Hinge Loss：?(y)=max(0,1?t?y)?(y)=max(0,1?t?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 一般形式的一致性優化問題（切割參數到各子系統，但各子系統目標函數參數維度不同，可能部分重合）
上述全局一致性優化問題，我們可以看到，所做的處理不過是對數據分塊，然后并行化處理。但是更一般的優化問題是，參數空間也是分塊的，即每個子目標函數fi(xi)fi(xi)的參數維度不同xi,∈Rnixi,∈Rni，我們稱之為局部變量。而局部變量所對應的的也將不再是全局變量zz，而是全局變量中的一部分zgzg，并且不是像之前的順序對應，而可能是隨便對應到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就是順序映射，也就是全局一致性優化問題，否則就不是。結合下圖就比較好理解
consensus
雖然如果用其他方法來做感覺會復雜，但是納入到上述ADMM框架，其實只不過是全局一致性優化問題的一個局部化變形，不過此時不是對數據進行分塊，是對參數空間進行分塊
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)
后續想做平均化處理，即中間會發生重合的參數zizi取值一樣的，那么zz-update將只能找他對應的那些xx進行平均化，也就是變成局部了，因為不是所有值都是要全局保持一致的。比如上面那個圖中的z1,z2,z3,z4z1,z2,z3,z4都分別只要求在部分xixi發生了共享需要保持一樣，而不是像之前全局要求每個xixi對應的都是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)
同全局一致性優化問題一樣，我們可以加上正則項，然后也可以變成帶正則項的一般形式的一致性優化問題。此處不贅述，與全局基本類似。
Sharing
4.4 共享問題（sharing）（橫向切割數據，也可縱向切變量）
與之前的全局變量一致性優化問題類似，共享問題也是一個非常一般而且常見的問題。他的形式如下：
minN∑i=1fi(xi)+g(N∑i=1xi)min∑i=1Nfi(xi)+g(∑i=1Nxi)
這里的第一部分局部損失fi(xi)fi(xi)與全局一致性優化是一樣的，即所有的xi∈Rn,i=1,…,Nxi∈Rn,i=1,…,N同維度，而對于一個共享的目標函數gg則是新加入的。在實際中，我們常常需要優化每個子數據集上的損失函數，同時還要加上全局數據所帶來的損失；或者需要優化每個子系統的部分變量，同時還要優化整個變量。共享問題是一個非常重要而靈活的問題，它也可以納入到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
上述形式當然還不夠簡潔，需要進一步化簡。因為xx-update可以不用擔心，分機并行處理優化求解即可，而對于zz-update這里面需要對NnNn個變量求解，想加快速度，就減少變量個數。于是想辦法通過和之前那種平均方式一樣來簡化形式解。
對于zz-update步，令ai=uki+xk+1iai=uik+xik+1，于是zz-update步優化問題轉化為
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
當ˉzzˉ固定時，那么后面的最優解（類似回歸）為zi=ai+ˉz?ˉazi=ai+zˉ?aˉ，帶入上式后于是后續優化就開始整體更新（均值化）
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問題是可以相互轉化的，可以同時達到最優，兩者存在著很緊密的對偶關系。
本節開頭提過，sharing問題用來切分數據做并行化，也可以切分參數空間做并行化。這對于高維、超高維問題是非常有好處的。因為高維統計中，大樣本是一方面問題，而高維度才是重中之重，如果能切分特征到低緯度中去求解，然后在合并起來，那么這將是一個很美妙的事情。上面利用regularized global consensus問題解決了切分大樣本數據的并行化問題，下面利用sharing思想解決常見的高維數據并行化問題
切割變量（特征）空間，并行化處理
同樣假設面對還是一個觀測陣A∈Rm×nA∈Rm×n和響應觀測b∈Rnb∈Rn，此時有n>>mn>>m，那么要么就降維處理，要么就切分維度去處理，或者對于超高維矩陣，切分維度后再降維。此時AA矩陣就不是像之前橫著切分，而是豎著切分，這樣對應著參數空間的切分：
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
當∥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區別不大，只是zz-update的損失函數不同，其余相同于是
ˉ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時候優化順序反了過來，與logistic rgression只是在zz-update步不同（損失函數不同）：
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)
當fifi是線性變化時，則退化成普通線性回歸。此時我們目標優化的問題是
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的正則，此時這個問題看起來似乎既可以對數據切分，也可以對特征切分，不過此時仍用sharing問題來做，相當于對特征切分為一個特征為一個子系統，于是有
$$ \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可以一塊一塊的求解？
最后再說一個經濟學中很重要的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
此時共享目標函數g=0g=0。xixi可以表示不同物品在NN個系統上的交換數量，(xi)j(xi)j可以表示物品jj從子系統ii上收到的交換數目，約束條件就可以看做在這些系統中物品交換是保持均衡穩定的。于是轉化為sharing問題，就有很簡單的ADMM解法（或者當做之前講過的受約束的凸優化問題來解，做投影）：
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 應用小總結
感覺上通過consensus problem和general consensus problem，我們可以看到并行和分布式部署優化方案的可行性。我們可以切分數據以及相應的目標函數，也可以切分變量到各個子系統上去，分別作優化，甚至我們可以大膽想象對不同類型數據塊用不同的優化算法，結合consensus問題和ADMM算法，達到同一個global variable的優化目的；或者對不同變量在不同類型數據塊上優化，即使有重疊，也可以結合general consensus思想和ADMM算法來解決這個問題。當然前提是能夠定義好需要估計的參數和優化的目標函數！大規模部署的前景還是很不錯的。下面具體分布式統計模型的構建便是ADMM算法非常好的應用。切分數據、切分變量（不過每個子系統的目標函數基本都是一樣的，其實應該可以不同）
5. Nonconvex問題
5.1 變量選擇（Regressor Selection）
5.2 因子模型（Factor Model Fitting）
5.3 雙凸優化（Bi-convex Problem）
非負矩陣分解（Nonnegative Matrix Factorization）
6. 具體實施與實際計算結果
這塊真的很實際，需要明白MPI的機理和Mapreduce、Graphlab等通信運作的機理，這樣才好部署ADMM算法，因為中間有很多迭代，需要做好子節點間參數與整體參數的通信，保持迭代時能同步更新參數。看實際運作，MPI和GraphLab可能更適合這種框架，Hadoop也是可以的，不過畢竟不是為迭代算法所生，要做好需要進行一些優化。Boyd提到Hadoop其中的Hbase更適合這種框架，因為Hbase是一種大表格，帶有時間戳，適合記錄迭代的記錄，這樣就不容易導致分布計算時候搞不清是哪一步的迭代結果了，導致通信調整比較復雜。不過在MapReduce框架下實施ADMM算法是沒有什么問題的，只要熟稔Hadoop的一些細節部分，基本沒有太大問題。
8. 總結
一個好的一般性算法，我個人覺得是易實施，并可大規模應用許多問題。可以讓統計學家卡在搞算法的瓶頸中解放出來，使得他們能快速用模擬，驗證自己構建可能較為復雜的模型。只有當看到一個令人感到欣慰的結果時，那些模型的統計性質的證明才可能是有意義的，如果事先連希望都看不到，那證明起來都可能底氣不足，讓人難以信服，更難以大規模應用統計學家所構建的模型。現在是一個高維數據、海量數據的年代，算法的重要性更會凸顯出來，一個好的模型如果沒有一個有效的算法支撐，那么他將可能什么都不是，Lasso頭幾年所遭遇的冷遇也充分證明了這一點，再比如在沒有計算機年代，Pearson的矩估計應用反而遠多于Fisher的MLE估計方法也是一個道理。好的一般性的解決方案，我想這不管是優化理論，還是統計等其他應用學科，雖然知道沒有最牛最終極的方法，但是能涌現一些大范圍適用的方法，那就是再好不過了。一招鮮吃遍天，人還都是喜歡簡單、安逸愛偷懶的嘛..
總結
- 上一篇： java索引ref_java – 如何
- 下一篇： Comet OJ - Contest #