亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        機(jī)器學(xué)習(xí)隨機(jī)優(yōu)化方法的個(gè)體收斂性研究綜述*

        2017-02-25 02:38:32張夢晗
        數(shù)據(jù)采集與處理 2017年1期
        關(guān)鍵詞:優(yōu)化方法

        陶 卿 馬 坡 張夢晗 陶 蔚

        (1.中國人民解放軍陸軍軍官學(xué)院十一系,合肥,230031; 2.解放軍理工大學(xué)指揮信息系統(tǒng)學(xué)院,南京,210007)

        機(jī)器學(xué)習(xí)隨機(jī)優(yōu)化方法的個(gè)體收斂性研究綜述*

        陶 卿1馬 坡1張夢晗1陶 蔚2

        (1.中國人民解放軍陸軍軍官學(xué)院十一系,合肥,230031; 2.解放軍理工大學(xué)指揮信息系統(tǒng)學(xué)院,南京,210007)

        隨機(jī)優(yōu)化方法是求解大規(guī)模機(jī)器學(xué)習(xí)問題的主流方法,其研究的焦點(diǎn)問題是算法是否達(dá)到最優(yōu)收斂速率與能否保證學(xué)習(xí)問題的結(jié)構(gòu)。目前,正則化損失函數(shù)問題已得到了眾多形式的隨機(jī)優(yōu)化算法,但絕大多數(shù)只是對迭代進(jìn)行平均的輸出方式討論了收斂速率,甚至無法保證最為典型的稀疏結(jié)構(gòu)。與之不同的是,個(gè)體解能很好保持稀疏性,其最優(yōu)收斂速率已經(jīng)作為open問題被廣泛探索。另外,隨機(jī)優(yōu)化普遍采用的梯度無偏假設(shè)往往不成立,加速方法收斂界中的偏差在有偏情形下會隨迭代累積,從而無法應(yīng)用。本文對一階隨機(jī)梯度方法的研究現(xiàn)狀及存在的問題進(jìn)行綜述,其中包括個(gè)體收斂速率、梯度有偏情形以及非凸優(yōu)化問題,并在此基礎(chǔ)上指出了一些值得研究的問題。

        機(jī)器學(xué)習(xí);隨機(jī)優(yōu)化;個(gè)體收斂性;有偏梯度估計(jì);非凸問題

        引 言

        通俗地說,機(jī)器學(xué)習(xí)的主要目的是讓計(jì)算機(jī)系統(tǒng)具有類似于人的學(xué)習(xí)能力。如何設(shè)計(jì)基于訓(xùn)練數(shù)據(jù)的學(xué)習(xí)算法從有限樣本集合得到分布意義下最優(yōu)的正確率,是統(tǒng)計(jì)機(jī)器學(xué)習(xí)研究的主要內(nèi)容[1]?;谟邢迾颖緦W(xué)習(xí)問題泛化界的統(tǒng)計(jì)分析,人們在機(jī)器學(xué)習(xí)算法設(shè)計(jì)中廣泛采用正則化損失函數(shù)的優(yōu)化準(zhǔn)則[2]為

        minr(w)+l(w)

        (1)

        隨著數(shù)據(jù)規(guī)模的急劇增加和各種應(yīng)用驅(qū)動學(xué)習(xí)范式的不斷涌現(xiàn),如何求解大規(guī)模不同學(xué)習(xí)范式導(dǎo)致的多樣性正則化損失函數(shù)優(yōu)化問題已經(jīng)成為機(jī)器學(xué)習(xí)領(lǐng)域亟需解決的關(guān)鍵性科學(xué)問題。實(shí)際上,自從支持向量機(jī)(Support vector machine,SVM)出現(xiàn)后,學(xué)習(xí)問題的優(yōu)化算法一直是機(jī)器學(xué)習(xí)領(lǐng)域眾多研究者普遍關(guān)注的一項(xiàng)議題。計(jì)算機(jī)領(lǐng)域內(nèi)的一些高水平國際會議經(jīng)常出現(xiàn)以優(yōu)化算法為主導(dǎo)的tutorial,機(jī)器學(xué)習(xí)頂級會議ICML(International conference on machine leaning),NIPS(Neural information processing systems)和COLT (Conference on learning theory)也曾多次舉辦workshop 進(jìn)行專門討論,機(jī)器學(xué)習(xí)頂級刊物JMLR(Journal of machine learning research)更是為此設(shè)立了special topic 和出版了special issue。一些在數(shù)學(xué)規(guī)劃領(lǐng)域有重要影響的團(tuán)隊(duì)也積極投身于大規(guī)模機(jī)器學(xué)習(xí)優(yōu)化問題研究的熱潮中,數(shù)學(xué)規(guī)劃領(lǐng)域的權(quán)威期刊SIOPT(SIAM journal on optimization)等頻頻出現(xiàn)求解機(jī)器學(xué)習(xí)優(yōu)化問題的論文。目前,依靠機(jī)器學(xué)習(xí)自身特點(diǎn)驅(qū)動而迅速發(fā)展起來的隨機(jī)優(yōu)化算法成為解決大規(guī)模問題的有效手段。如何在這些優(yōu)化算法中保持學(xué)習(xí)問題的正則化項(xiàng)結(jié)構(gòu)以及獲取最優(yōu)收斂速率是研究中的核心問題。在正則化機(jī)器學(xué)習(xí)問題及其優(yōu)化算法方面,已經(jīng)發(fā)表了很多綜述性論文,如文獻(xiàn)[3~8]。特別,在文獻(xiàn)[7]中以損失函數(shù)和優(yōu)化求解為主線,對統(tǒng)計(jì)機(jī)器學(xué)習(xí)的正則化損失函數(shù)框架進(jìn)行了綜述和分析,具體討論了幾何margin和損失函數(shù)兩種理論分析方法,重點(diǎn)講述了損失函數(shù)的漸近最優(yōu)性,闡述了數(shù)學(xué)優(yōu)化方法和機(jī)器學(xué)習(xí)優(yōu)化算法之間的區(qū)別和聯(lián)系,強(qiáng)調(diào)了任何關(guān)于損失函數(shù)和優(yōu)化求解方法的進(jìn)展都將會促進(jìn)機(jī)器學(xué)習(xí)的發(fā)展。文獻(xiàn)[8]首先指出大規(guī)模機(jī)器學(xué)習(xí)問題的訓(xùn)練樣本集合往往具有冗余和稀疏的特點(diǎn),機(jī)器學(xué)習(xí)優(yōu)化問題中的正則化項(xiàng)和損失函數(shù)也蘊(yùn)含著特殊的結(jié)構(gòu)含義,直接使用整個(gè)目標(biāo)函數(shù)梯度的批處理黑箱方法不僅難以處理大規(guī)模問題,也無法滿足機(jī)器學(xué)習(xí)對結(jié)構(gòu)的要求,然后針對L1正則化問題,介紹了依據(jù)機(jī)器學(xué)習(xí)自身特點(diǎn)驅(qū)動而迅速發(fā)展起來的坐標(biāo)優(yōu)化、在線和隨機(jī)優(yōu)化方法的一些研究進(jìn)展。目前,很多經(jīng)典的一階優(yōu)化方法經(jīng)過適當(dāng)?shù)母脑?,不僅可以保證正則化項(xiàng)的結(jié)構(gòu),同時(shí)還具有不依賴正則化項(xiàng)光滑性的最優(yōu)收斂速率,但由于最終解大多都采取了加權(quán)平均的輸出方式,仍然存在著一些弊端。特別是對L1 正則化問題,盡管每一步迭代產(chǎn)生的解具有稀疏性,但平均求和的最終輸出方式卻破壞了這種求解算法最初極力維護(hù)的稀疏性。為了更好地保持正則化項(xiàng)的結(jié)構(gòu),應(yīng)該研究個(gè)體收斂性問題,但即使是對數(shù)學(xué)規(guī)劃中討論的標(biāo)準(zhǔn)凸優(yōu)化問題,很多經(jīng)典文獻(xiàn)和書籍也都缺乏個(gè)體收斂速率的論述[9]。一般地說,時(shí)空復(fù)雜性是衡量一個(gè)算法優(yōu)劣的主要標(biāo)準(zhǔn)。結(jié)合機(jī)器學(xué)習(xí)問題的實(shí)際含義,人們評判學(xué)習(xí)問題(1)優(yōu)化算法的依據(jù)主要是收斂速率的界、正則化項(xiàng)結(jié)構(gòu)的保持以及存儲復(fù)雜性。因此,本文綜述也在首先介紹隨機(jī)結(jié)構(gòu)優(yōu)化算法的基礎(chǔ)上,按照個(gè)體最優(yōu)速率這一評判標(biāo)準(zhǔn)展開,這也是本文與文獻(xiàn)[3~8]的不同之處。另外,當(dāng)標(biāo)準(zhǔn)隨機(jī)優(yōu)化算法一些普遍使用的如樣本集合獨(dú)立同分布假設(shè)不成立時(shí),本文還介紹了一些重要的拓廣問題。

        1 隨機(jī)優(yōu)化算法

        目前,機(jī)器學(xué)習(xí)領(lǐng)域的研究者們主要借助數(shù)學(xué)規(guī)劃特別是凸優(yōu)化領(lǐng)域的一階梯度方法,考慮無約束的優(yōu)化問題

        minf(w)

        (2)

        式中:f(w)為定義在Rn上可微的凸函數(shù),記w*是問題(2)的一個(gè)最優(yōu)解。對于問題(2),批處理形式經(jīng)典的梯度下降方法的迭代步驟為

        (3)

        式中:at為迭代步長,f(wt)為f(w)在wt處的梯度。梯度下降方法主要在每一點(diǎn)處目標(biāo)函數(shù)梯度的相反方向會使目標(biāo)函數(shù)下降最快這一事實(shí)。在梯度下降方法的基礎(chǔ)上,針對約束或者非光滑問題,優(yōu)化方法出現(xiàn)了很多變形,其中包括投影次梯度算法[9]、鏡像下降[10]、對偶平均[11]和交替方向乘子方法[12]等。但對于正則化加損失函數(shù)這種具有明確機(jī)器學(xué)習(xí)含義的優(yōu)化問題(1)來說,仍然屬于黑箱方法。對于求解特定領(lǐng)域的問題, 優(yōu)化理論方面著名學(xué)者Nesterov曾經(jīng)指出“黑箱方法在凸優(yōu)化問題上的重要性將不可逆轉(zhuǎn)地消失,徹底地取而代之的是巧妙運(yùn)用問題結(jié)構(gòu)的新算法”[11]。因此,為了利用這些經(jīng)典算法有效處理大規(guī)模具有實(shí)際含義的機(jī)器學(xué)習(xí)問題,還需要在使用形式上進(jìn)行必要的改變。

        首先,這些批處理形式的優(yōu)化方法由于每次迭代都要涉及到損失函數(shù)梯度的計(jì)算,從而不可避免地需要遍歷訓(xùn)練樣本集合,該操作方式顯然無法適用于大規(guī)模機(jī)器學(xué)習(xí)問題的求解。對于優(yōu)化問題(1),梯度下降方法隨機(jī)形式可表示為

        (4)

        式中:i為從樣本集合中第t+1次隨機(jī)抽取樣本的序號。廣義上來說,隨機(jī)優(yōu)化算法的每步迭代僅需要知道目標(biāo)函數(shù)梯度的無偏估計(jì),特別對于有限樣本的機(jī)器學(xué)習(xí)問題來說,由于假設(shè)樣本是獨(dú)立同分布的,單個(gè)樣本對應(yīng)目標(biāo)函數(shù)的梯度就是整個(gè)訓(xùn)練集上目標(biāo)函數(shù)梯度的無偏估計(jì),從而隨機(jī)優(yōu)化方法只需要計(jì)算部分甚至單個(gè)樣本對應(yīng)目標(biāo)函數(shù)的梯度,這就克服了批處理算法每次迭代都需要遍歷訓(xùn)練樣本集合的固有缺陷。另一方面,由于大規(guī)模學(xué)習(xí)問題訓(xùn)練樣本往往存在著冗余現(xiàn)象,實(shí)際應(yīng)用中往往只需運(yùn)行隨機(jī)優(yōu)化算法少許迭代步驟后,學(xué)習(xí)精度就已經(jīng)呈現(xiàn)出穩(wěn)定的趨勢[13]。2007 年,Shalev-Shwartz 使用隨機(jī)投影次梯度對大規(guī)模SVM 進(jìn)行求解(稱為Pegasos)[14],取得了轟動一時(shí)的實(shí)際效果。即使是在10年后的今天,很多新算法的比較對象中仍然還會有Pegasos 的身影出現(xiàn)。

        其次,文獻(xiàn)[15]指出,當(dāng)將L1 正則化項(xiàng)和損失函數(shù)整體作為目標(biāo)函數(shù)使用一階梯度隨機(jī)優(yōu)化方法時(shí),卻無法獲得L1 正則化項(xiàng)應(yīng)該帶來的稀疏性,這表明黑箱優(yōu)化方法難以保證優(yōu)化問題中正則化項(xiàng)的結(jié)構(gòu)。為了解決這一問題,Xiao 等在2009 年提出了一種對保持正則化項(xiàng)結(jié)構(gòu)具有重要影響的算法稱之為RDA(Regularized dual average),成功地將對偶平均優(yōu)化算法推廣至正則化情形[16]。隨后,Duchi 等將鏡像下降算法推廣至正則化情形,稱之為COMID(Composite objective mirror descent)[17]。對于問題(2),批處理形式鏡像下降方法的迭代步驟為

        (5)

        (6)

        當(dāng)r(w)為L1范數(shù)時(shí),不難發(fā)覺COMID 在優(yōu)化過程中將正則化項(xiàng)和損失函數(shù)分別看待,只是對損失函數(shù)采用相應(yīng)的一階梯度方法,從而很好地保證了L1正則化問題每一步迭代解的稀疏性。這樣的研究思路目前已被廣泛采用,有時(shí)人們也將使用了這種技巧的優(yōu)化方法稱為算法的Proximal 形式[18-19]或者更加明確地稱為結(jié)構(gòu)優(yōu)化方法[16]。除了能夠保證正則化項(xiàng)的結(jié)構(gòu)之外,結(jié)構(gòu)優(yōu)化方法的重要理論意義還在于其收斂速率僅與損失函數(shù)的光滑性有關(guān),即盡管目標(biāo)函數(shù)非光滑,但只要損失函數(shù)滿足光滑性條件,結(jié)構(gòu)優(yōu)化方法仍然能夠獲得目標(biāo)函數(shù)光滑時(shí)才能達(dá)到的最優(yōu)收斂速率,這表明結(jié)構(gòu)優(yōu)化方法的最優(yōu)收斂速率界比黑箱方法有數(shù)量級的提升。因此,理論分析和實(shí)際應(yīng)用均表明經(jīng)典優(yōu)化方法的隨機(jī)結(jié)構(gòu)化形式在處理大規(guī)模學(xué)習(xí)問題方面具有獨(dú)特的優(yōu)勢。值得指出的是,在近期優(yōu)化算法的專題討論中,經(jīng)典優(yōu)化方法的隨機(jī)結(jié)構(gòu)化形式占了相當(dāng)大的比重,一些新的研究動向與進(jìn)展甚至對經(jīng)典數(shù)學(xué)規(guī)劃理論本身的發(fā)展產(chǎn)生了積極的影響。

        2 最優(yōu)收斂速率

        眾所周知,對于不同凸性和光滑性的優(yōu)化問題,確定型批處理形式優(yōu)化算法的最優(yōu)收斂速率已經(jīng)有了很多定論[20]。顯然,評判在機(jī)器學(xué)習(xí)中起重要作用的隨機(jī)結(jié)構(gòu)優(yōu)化算法的收斂性是理論分析首要面對的問題。對非光滑損失函數(shù)問題,鏡像下降方法和對偶平均方法均已推廣至隨機(jī)結(jié)構(gòu)化情形,RDA 和COMID 對于無強(qiáng)凸性假設(shè)的一般凸損失函數(shù)情形均能達(dá)到最優(yōu)的收斂速率[16-17]??傮w來說,非光滑一般凸問題的研究較為平淡,在獲取平均方式輸出的最優(yōu)收斂速率界方面,很多算法都沒有碰到特別棘手的困難。相比于非光滑情形,光滑損失函數(shù)的優(yōu)化研究比較豐富。這主要是因?yàn)楫?dāng)目標(biāo)函數(shù)光滑時(shí),存在著一種加速技巧,可以將梯度方法的收斂速率提升一個(gè)數(shù)量級。1983 年,Nesterov 首先提出了這種里程碑式的加速方法[21],獲得了目標(biāo)函數(shù)光滑時(shí)批處理優(yōu)化算法的最優(yōu)收斂速率。隨后,這種技巧出現(xiàn)了多種形式的變形[18,22]。因此,大量的隨機(jī)優(yōu)化研究都圍繞光滑損失函數(shù)優(yōu)化問題的加速上,Nesterov 的加速方法及其變形也已成功地推廣至隨機(jī)結(jié)構(gòu)化情形[16,23]。這些較早期的隨機(jī)優(yōu)化算法收斂性分析幾乎都是先建立相應(yīng)在線算法的regret界,隨后通過在線與隨機(jī)算法之間的切換技巧[24],得到隨機(jī)算法對所有迭代進(jìn)行平均作為輸出方式的收斂速率。但是,這種通過在線與隨機(jī)切換而獲得收斂速率的方法在處理強(qiáng)凸優(yōu)化時(shí)卻遭遇了一些困境。這是因?yàn)楫?dāng)目標(biāo)函數(shù)強(qiáng)凸時(shí),即使是采用對所有迭代進(jìn)行平均的輸出方式,目前標(biāo)準(zhǔn)的SGD(Stochastic gradient descent) 也未能證明能獲得最優(yōu)的收斂速率界。由此一些學(xué)者甚至對SGD 在經(jīng)典優(yōu)化理論中的統(tǒng)治地位產(chǎn)生了質(zhì)疑,強(qiáng)凸問題優(yōu)化算法最優(yōu)收斂速率問題也引起了研究者的普遍關(guān)注。相比于一般凸問題,強(qiáng)凸目標(biāo)函數(shù)的收斂速率研究可謂是精彩紛呈?;\統(tǒng)地說,存在著兩種主要的思路:(1)對算法的迭代過程進(jìn)行修改,如Hazan 等在SGD 中嵌入適當(dāng)數(shù)目的內(nèi)循環(huán),提出一種Epoch-GD 算法[24],獲得了平均輸出方式的最優(yōu)收斂速率。(2)僅對輸出方式進(jìn)行修改,如Rakhlin 等以后半部分迭代解平均的α-suffix 技巧來代替全部平均的輸出方式[25],Lacoste-Julien等提出了一種加權(quán)平均的輸出方式[26],這些輸出方式均使得SGD 獲得了最優(yōu)的收斂速率,其中加權(quán)平均克服了α-suffix 技巧輸出不能on-the-fly 計(jì)算的問題。

        上述討論的優(yōu)化方法不僅獲得了最優(yōu)收斂速率,而且計(jì)算復(fù)雜性和收斂速率的界還與樣本數(shù)目無關(guān),因此理論上不需要存儲樣本集合,可以解決任意規(guī)模的問題。此時(shí),也可以認(rèn)為在訓(xùn)練樣本獨(dú)立同分布的假設(shè)下,這些算法所優(yōu)化目標(biāo)函數(shù)中含有基于損失函數(shù)的期望風(fēng)險(xiǎn)[2,19]。但是在這種情形下,即使是黑箱方法,當(dāng)目標(biāo)函數(shù)強(qiáng)凸且光滑時(shí),標(biāo)準(zhǔn)的隨機(jī)梯度優(yōu)化方法只能獲得次線性的收斂速率界,與全梯度方法的線性收斂速率未能很好匹配,這種差異說明僅僅依靠目標(biāo)函數(shù)梯度的無偏估計(jì)就期望得到全梯度方法的效果不現(xiàn)實(shí)[19,27],還需要對優(yōu)化問題進(jìn)行一定的限制。

        在實(shí)際應(yīng)用中,通常給定樣本集合,此時(shí)的機(jī)器學(xué)習(xí)優(yōu)化問題表現(xiàn)為正則化經(jīng)驗(yàn)風(fēng)險(xiǎn)的形式。通過充分利用經(jīng)驗(yàn)風(fēng)險(xiǎn)這種“有限和”形式的結(jié)構(gòu),優(yōu)化算法往往可以得到更好的收斂速率界。自然,處理這種“有限和”形式目標(biāo)函數(shù)的優(yōu)化算法也可以采取每一步迭代僅對其中一個(gè)損失函數(shù)進(jìn)行操作的方式。為了區(qū)別起見,將這樣的優(yōu)化算法稱為增量算法[9,28],如果以隨機(jī)方式選取操作的損失函數(shù),此時(shí)的算法稱為隨機(jī)增量優(yōu)化算法。從形式上看,增量算法的迭代方式似乎與SGD完全相同,但關(guān)鍵的不同是,增量算法可以存儲部分梯度或迭代信息,這正是增量算法比標(biāo)準(zhǔn)隨機(jī)方法具有更好收斂速率界的主要原因。

        關(guān)于“有限和”形式增量算法的研究成果層出不窮,其中代表性的論文有SAG[27,29],SAGA[30],F(xiàn)inito[28]和MISO[31]等。不論是處理黑箱還是Proximal 形式的問題,當(dāng)目標(biāo)函數(shù)強(qiáng)凸且光滑(損失函數(shù))時(shí),這些方法均可獲得和批處理全梯度算法相同階的收斂速率界。更進(jìn)一步,為了避免增量算法在處理大規(guī)模問題時(shí)的存儲瓶頸,Johnson 等將減小方差策略嵌入到SGD中,稱為SVRG[32]。SVRG 的主要思路是使用全梯度對SGD 迭代中的梯度進(jìn)行修正,從而達(dá)到了減少方差的目的,同時(shí)也付出了不得不周期性地遍歷樣本集合的代價(jià)。但對于光滑強(qiáng)凸優(yōu)化問題,SVRG 在和標(biāo)準(zhǔn)SGD 具有幾乎相同存儲需求的條件下,實(shí)現(xiàn)了與SAG 一樣的線性收斂速率。隨后,SVRG 被成功拓廣至Proximal情形,對目標(biāo)函數(shù)強(qiáng)凸但損失函數(shù)光滑的問題也獲得了線性的收斂速率[19]。另外,一些研究者還對強(qiáng)凸與只有光滑條件的增量算法作了統(tǒng)一的處理[33],還有一些研究工作使用Nesterov 技巧對增量算法也進(jìn)行了加速,得到了最優(yōu)收斂速率界[34-35]。

        3 個(gè)體收斂速率

        使用迭代過程中的個(gè)體直接作為輸出在稀疏學(xué)習(xí)問題中具有更明確的實(shí)際意義,能夠比平均方式的解獲得更高的稀疏效果。這說明為了更好地保持正則化項(xiàng)的結(jié)構(gòu),應(yīng)該研究個(gè)體收斂性問題。但即使是對數(shù)學(xué)規(guī)劃中專門討論的投影次梯度批處理優(yōu)化方法,如果沒有目標(biāo)函數(shù)光滑性假設(shè)保證下的單調(diào)性[18,22],很多經(jīng)典文獻(xiàn)和書籍也都缺乏個(gè)體收斂速率的論述[9]。因此,個(gè)體收斂性在數(shù)學(xué)規(guī)劃理論研究中也有其重要的地位。下面對個(gè)體收斂速率的主要進(jìn)展和存在的問題進(jìn)行簡單的梳理。

        對于標(biāo)準(zhǔn)隨機(jī)優(yōu)化算法SGD,只有當(dāng)目標(biāo)函數(shù)不僅強(qiáng)凸而且光滑時(shí),才能相對比較容易地獲得個(gè)體收斂速率界[25]。對于單純的強(qiáng)凸問題,如果不改變算法本身,也必須對平均的輸出方式進(jìn)行修改,才能獲得最優(yōu)收斂速率[25-26],但研究者對這些結(jié)果似乎不滿意。實(shí)際上,人們最為期待的莫過于SGD 對于強(qiáng)凸問題能否達(dá)到最優(yōu)個(gè)體收斂速率,這個(gè)問題看似簡單,截至目前卻始終沒有答案。為了捍衛(wèi)SGD 的經(jīng)典與尊嚴(yán),Shamir 在2012 年的機(jī)器學(xué)習(xí)頂級會議COLT 上把強(qiáng)凸目標(biāo)函數(shù)下SGD 的個(gè)體最優(yōu)收斂速率作為open 問題提出[36]。2013 年,Shamir 等提出一種由平均輸出方式收斂速率得到個(gè)體收斂速率的一般技巧[37],盡管得到了眾盼所歸的SGD 個(gè)體收斂速率,但獲得的收斂速率界與平均輸出方式的收斂速率界相差一個(gè)對數(shù)因子,顯然未能達(dá)到最優(yōu)收斂速率界。在隨機(jī)優(yōu)化算法的個(gè)體收斂速率研究方面,Chen 等提出的Optimal RDA 獲得了比較全面而又非常理想的結(jié)果[38]。其主要的思路是對對偶平均方法進(jìn)行改進(jìn),在每一步迭代中增加一個(gè)不同形式的子優(yōu)化問題求解,對研究者通常獨(dú)立討論的一般凸、強(qiáng)凸或光滑等類型的問題,均獲得了個(gè)體最優(yōu)收斂速率。2015 年,Nesterov 等在對偶平均方法的迭代中巧妙地嵌入了一種線性插值操作,證明了該方法在一般凸情形下具有最優(yōu)的個(gè)體收斂速率,并且這種個(gè)體收斂呈現(xiàn)出與平均方式收斂同樣的穩(wěn)定性[39]。從理論角度來說,這種改動與標(biāo)準(zhǔn)的對偶平均方法區(qū)別極小,是對偶平均方法一種很好的擴(kuò)展,也是對一階梯度方法個(gè)體最優(yōu)收斂速率比較接近大家期待的一種回答。

        仔細(xì)分析可以發(fā)現(xiàn),文獻(xiàn)[38]中的Optimal RDA 和文獻(xiàn)[39]中的線性插值操作技巧獲得個(gè)體收斂速率的原理不同,以至于獲得的學(xué)習(xí)結(jié)果也不相同。一般地說,在處理L1 正則化問題時(shí),Optimal RDA 的個(gè)體解具有很好的稀疏性,但卻不具有收斂的穩(wěn)定性,這也是子優(yōu)化問題的解直接作為最終輸出的通用弊?。欢€性插值技巧嵌入在迭代過程的梯度運(yùn)算后,從而使最終的個(gè)體解具有很好的收斂穩(wěn)定性,但卻不具有稀疏性,這也是將插值累積作為最終輸出的通用弊病。另外,這兩種個(gè)體收斂速率分析的思路目前僅適用于步長策略靈活的對偶平均方法。值得指出,文獻(xiàn)[38]對對偶平均算法的改動很大,這實(shí)際上與標(biāo)準(zhǔn)SGD 個(gè)體收斂速率open問題的本意已有所偏離,而線性插值操作技巧[39]對于強(qiáng)凸或光滑等目標(biāo)函數(shù)情形能否得到個(gè)體最優(yōu)收斂速率也未討論。最近,文獻(xiàn)[40]提出了一種嵌入線性插值操作的投影次梯度方法,證明了其在一般凸情形下具有個(gè)體最優(yōu)收斂速率。更進(jìn)一步,將所獲結(jié)論推廣至隨機(jī)方法情形。

        4 隨機(jī)優(yōu)化算法的拓廣

        目前大多數(shù)的機(jī)器學(xué)習(xí)隨機(jī)優(yōu)化方法研究都假設(shè)隨機(jī)抽取的部分甚至單個(gè)樣本點(diǎn)對應(yīng)的目標(biāo)函數(shù)梯度是整個(gè)目標(biāo)函數(shù)梯度的無偏估計(jì)。在這種假設(shè)條件下,特別是對光滑損失函數(shù)的優(yōu)化問題,如前所述,人們得到了眾多形式的加速算法,這些算法具有最優(yōu)的收斂速率且收斂速率的界與樣本數(shù)目無關(guān)。

        但在實(shí)際問題中,梯度無偏估計(jì)的假設(shè)往往是不成立的。(1)當(dāng)在給定訓(xùn)練集合的條件下企圖優(yōu)化基于損失函數(shù)期望形式的目標(biāo)函數(shù)時(shí),根本無法知曉樣本集合是否滿足獨(dú)立同分布條件,此時(shí)當(dāng)然認(rèn)為梯度估計(jì)出現(xiàn)有偏更為合理。(2)L1 正則化問題是一種最為簡單典型的正則化損失函數(shù)優(yōu)化問題,其簡單之處在于在每一步迭代過程中涉及的優(yōu)化子問題可以解析求解,從而整個(gè)隨機(jī)算法具有非常理想的計(jì)算代價(jià),但對于 Total-variation,Nuclear-norm等類型的正則化項(xiàng)和Fused-Lasso 等問題,子問題只能近似求解,這種求解方式也可以視為是一種梯度估計(jì)有偏情形下的精確求解[41]。(3)在光滑損失函數(shù)的優(yōu)化算法中,需要知道損失函數(shù)的Lipschitz 常數(shù),這個(gè)決定步長的參數(shù)對算法性能會產(chǎn)生嚴(yán)重的影響,Lipschitz 常數(shù)估計(jì)中的誤差問題也可以歸結(jié)為一種廣義的梯度無偏估計(jì)問題[42]。優(yōu)化領(lǐng)域著名學(xué)者Nesterov 領(lǐng)導(dǎo)的研究小組對梯度有偏情形進(jìn)行了深入的研究,他們定義一般形式的梯度無偏估計(jì)問題,目前很多文獻(xiàn)中所關(guān)注的具體情形都是其特例[42]。一個(gè)令人驚訝的事實(shí)是所有研究者對所涉及的具體梯度估計(jì)有偏問題都得到了一致相同的結(jié)論,即不論是批處理還是隨機(jī)優(yōu)化算法,非加速算法收斂速率界保持梯度偏差為一常量,而加速算法平均輸出形式收斂速率的界卻會出現(xiàn)隨著迭代的增加而累加定性誤差的現(xiàn)象,這使得加速算法無任何理論保證, 處于十分尷尬的境地[41,43-44]。對于非加速算法,在假定了廣義梯度偏差以一定速率衰減的條件下,可以獲得通常意義下極限為0 的收斂速率界,且這個(gè)結(jié)論可以用來求解一種矩陣分解問題,其中的偏差由原問題和對偶問題之間的間隙控制[41]。

        對于正則化損失函數(shù)形式的優(yōu)化問題,有時(shí)人們對正則化項(xiàng)或損失函數(shù)項(xiàng)強(qiáng)加了更多的學(xué)習(xí)含義。例如,為了獲得更稀疏的支持向量,人們對hinge 損失采取了強(qiáng)行截?cái)嗟氖侄蝃45-46];同樣地,為了獲得更好的特征稀疏度,人們對L1 正則化項(xiàng)也采取了截?cái)嗟氖址╗47]。與標(biāo)準(zhǔn)的凸優(yōu)化模型相比,這些截?cái)嗄P鸵簿哂泻芎玫慕y(tǒng)計(jì)含義,并體現(xiàn)了正則化項(xiàng)和損失函數(shù)對學(xué)習(xí)問題的原本目標(biāo)更精確的逼近。但不可回避的是,這些模型都帶來一個(gè)難以處理的問題,即如何求解這些對學(xué)習(xí)問題有額外需求而導(dǎo)致的非凸優(yōu)化問題。非凸優(yōu)化問題一直是數(shù)學(xué)規(guī)劃領(lǐng)域的難解問題,只是在一些特定的假設(shè)下,才能得到一些局部收斂的算法。MM(Majorization-minimization)方法可以視為處理非凸優(yōu)化問題的一般框架,在迭代過程中,MM 算法通過優(yōu)化目標(biāo)函數(shù)局部凸上界的形式避開了目標(biāo)函數(shù)的非凸性[48]。著名算法CCCP(Concave convex procedure)[49]或DC (Difference convex)規(guī)劃[50]實(shí)際上是求解兩個(gè)凸函數(shù)之差形式非凸問題的特例,并且人們已經(jīng)使用CCCP 和DC 規(guī)劃處理截?cái)嗯幚硇问降膬?yōu)化問題,均獲得了比求解凸優(yōu)化問題更令人滿意的學(xué)習(xí)結(jié)果[45-47],尤其是關(guān)于非線性截?cái)郤VM 的求解[45],由于支持向量數(shù)目的減少,大大縮短了求解優(yōu)化問題的時(shí)間,這篇論文也獲得了ICML2006 的最佳論文獎(jiǎng)。自然,人們會想到能否在非凸批處理算法的基礎(chǔ)上進(jìn)一步拓廣隨機(jī)優(yōu)化方法求解大規(guī)模非凸優(yōu)化問題。2015 年,針對“有限和”形式的非凸優(yōu)化問題,受SAG 方法[27]的啟發(fā),Mairal 給出了針對光滑目標(biāo)函數(shù)增量形式的MM 算法[31]。梯度有偏情形的大規(guī)模優(yōu)化問題和非凸問題的研究同屬于凸優(yōu)化問題隨機(jī)優(yōu)化算法應(yīng)用范圍的一種拓廣。目前,也只是在比較嚴(yán)格的假設(shè)下,才獲得了與理想情形類似的收斂速率界,關(guān)于個(gè)體收斂速率界的結(jié)果還未見報(bào)道,非“有限和”形式目標(biāo)函數(shù)的探討更是極少。

        5 結(jié)束語

        綜上所述,隨機(jī)優(yōu)化算法是求解大規(guī)模機(jī)器學(xué)習(xí)問題的主流優(yōu)化方法之一,具有堅(jiān)實(shí)的理論基礎(chǔ)和良好的發(fā)展前景,但在大家普遍關(guān)注的問題中仍然存在著很多令人無法回避而又亟待需要克服的缺陷。根據(jù)上述對隨機(jī)優(yōu)化研究進(jìn)展的梳理,以下幾個(gè)方面值得研究:

        (1)對偶平均方法經(jīng)過一些改進(jìn),無論是對一般凸、強(qiáng)凸還是光滑的損失函數(shù),已經(jīng)具有了最優(yōu)的個(gè)體收斂速率,但在收斂穩(wěn)定性或保持正則化項(xiàng)結(jié)構(gòu)方面仍然不能令人完全滿意。因此,能否對目前主流的一階隨機(jī)梯度方法進(jìn)行適當(dāng)?shù)母倪M(jìn),獲得最優(yōu)個(gè)體收斂速率并高效保證正則化項(xiàng)結(jié)構(gòu),同時(shí)具有收斂穩(wěn)定性,這是經(jīng)典梯度方法在機(jī)器學(xué)習(xí)優(yōu)化問題應(yīng)用中必須解決的問題,也是對SGD 個(gè)體收斂速率open 問題比較接近的回答。

        (2)實(shí)際問題中往往會出現(xiàn)目標(biāo)函數(shù)梯度估計(jì)有偏情形,假設(shè)梯度估計(jì)有偏是從L1 正則化損失函數(shù)優(yōu)化問題過渡到一般正則優(yōu)化問題的一種有效途徑,也是處理參數(shù)不確定性的一種手段。遺憾的是,無偏情形下一些具有最優(yōu)收斂速率隨機(jī)算法的收斂速率界卻會隨迭代次數(shù)累計(jì)增加偏差項(xiàng)。因此,一些應(yīng)用中必須研究的問題,也是個(gè)體收斂性研究中值得研究的問題。

        (3)由于實(shí)際應(yīng)用的驅(qū)動,機(jī)器學(xué)習(xí)領(lǐng)域也出現(xiàn)了一些具有特定含義的非凸正則化問題,目前蓬勃發(fā)展并取得空前成功的深度神經(jīng)網(wǎng)絡(luò)也強(qiáng)烈地依賴于非凸優(yōu)化問題的求解[32]。非凸問題研究雖然有一些進(jìn)展,但還存在相當(dāng)多的問題沒有涉及。因此,研究非凸問題的個(gè)體最優(yōu)收斂速率界與穩(wěn)定性問題具有理論意義和很好的應(yīng)用前景。

        總之,隨機(jī)優(yōu)化方法的個(gè)體最優(yōu)收斂速率、收斂穩(wěn)定性及其拓廣的相關(guān)研究具有重要的理論意義和實(shí)用價(jià)值,在當(dāng)前機(jī)器學(xué)習(xí)優(yōu)化方法發(fā)展中具有一定的挑戰(zhàn)性問題,這項(xiàng)研究對算法的進(jìn)一步分布式實(shí)現(xiàn)也具有指導(dǎo)意義。

        [1] Vapnik V N. Statistical learning theory[M]. New York: Wiley-Interscience, 1998.

        [2] Zhang T. Statistical behavior and consistency of classification methods based on convex risk minimization[J]. Annals of Statistic, 2004,32:56-85.

        [3] Duchi J. Introductory lectures on stochastic convex optimization [EB/OL]. Park City Mathematics Institute, Graduate Summer School Lectures, http://web.stanford.edu/~jduchi/, 2016.

        [4] Bottou L, Curtis F, Nocedal J. Optimization methods for large-scale machine learning[EB/OL].https://arxiv.org/abs/1606.04838, 2016-06-15[2016-06-15].

        [5] 吳啟暉,邱俊飛,丁國如.面向頻譜大數(shù)據(jù)處理的機(jī)器學(xué)習(xí)方法[J].數(shù)據(jù)采集與處理,2015,30(4):703-731.

        Wu Qihui, Qiu Junfei, Ding Guoru. Maching learning methods for big spectrum data processing[J]. Journal of Data Acquistion and Processing,2015,30(4):703-731.

        [6] 潘志松,唐斯琪,邱俊洋,等.在線學(xué)習(xí)算法綜述[J].數(shù)據(jù)采集與處理,2016,31(6):1067-1082.

        Pan Zhisong,Tang Siqi,Qiu Junyang,et al. Survey on online learning algorithms[J]. Journal of Data Acquistion and Processing,2016,31(6):1067-1082.

        [7] 孫正雅,陶卿.統(tǒng)計(jì)機(jī)器學(xué)習(xí)綜述:損失函數(shù)與優(yōu)化求解[J].中國計(jì)算機(jī)學(xué)會通訊,2009,5(8):7-14.

        Sun Zhengya,Tao Qing.The Loss function and optimizer in statistical machine learning[J]. Communications of the Chinese Computer Federation,2009,5(8):7-14.

        [8] 陶卿,高乾坤,姜紀(jì)遠(yuǎn),等.稀疏學(xué)習(xí)優(yōu)化問題的求解綜述[J]. 軟件學(xué)報(bào),2013,24(11):2498-2507.

        Tao Qing, Gao Qiankun, Jiang Jiyuan, et al. Survey of solving the optimization problems for sparse learning[J]. Journal of Software,2013,24(11):2498-2507.

        [10]Beck A, Teboulle M. Mirror descent and nonlinear projected sub-gradient methods for convex optimization[J]. Oper Res Lett, 2003,31:167-175.

        [11]Nesterov Y. Primal-dual subgradient methods for convex problems[J]. Math Program, 2009,120:261-283.

        [12]Boyd S, Parikh N, Chu E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends in Machine Learning, 2011, 3(1):1-122.

        [13]Zhang T. Solving large scale linear prediction problems using stochastic gradient descent algorithms[C]//Proceedings of the International Conference on Machine Leaning. New York: ACM, 2004: 919-926.

        [14]Shalev-Shwartz S, Singer Y, Srebro N. Pegasos: Primal estimated sub-gradient solver for SVM[J]. Mathematical Programming, 2011, 127(1):3-30.

        [15]Langford J. Sparse online learning via truncated gradient[J]. Journal of Machine Learning Research, 2008, 10(2):777-801.

        [16]Xiao L. Dual averaging methods for regularized stochastic learning and online optimization[J]. Journal of Machine Learning Research, 2010,11:2543-2596.

        [17]Duchi J, Shalev-Shwartz S, Singer Y, et al. Composite objective mirror descent[C]// Proceedings of the Conference on Learning Theory. Haifa, Israel:[s.n.],2010:14-26.

        [18]Tseng P. Approximation accuracy, gradient methods, and error bound for structured convex optimization[J]. Mathematical Programming, 2010,125(2):263-295.

        [19]Xiao L, Zhang T. A proximal stochastic gradient method with progressive variance reduction[J]. SIAM Journal on Optimization, 2014, 24(4):2057-2075.

        [20]Nemirovsky A S, Yudin D B. Problem complexity and method efficiency in optimization[M]. New York: Wiley, 1983.

        [21]Nesterov Y. A method of solving a convex programming problem with convergence rateO(1/k2)[J]. Soviet Mathematics Doklady, 1983,27(2):372-376.

        [22]Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM Journal on Imaging Sciences, 2009, 2(1):183-202.

        [23]Hu C, Kwok J T, Pan W. Accelerated gradient methods for stochastic optimization and online learning[C]//Advances in Neural Information Processing Systems.Vancouver,British Columbia,Canada:[s.n.] ,2009:781-789.

        [24]Hazan E, Kale S. Beyond the regret minimization barrier: An optimal algorithm for stochastic strongly-convex optimization[J].Journal of Machine Learning Research, 2011,15(1):2489-2512.

        [25]Rakhlin A, Shamir O, Sridharan K. Making gradient descent optimal for strongly convex stochastic optimization[C]//Proceedings of the International Conference on Machine Learning.Edinburgh,Scotland:[s.n.], 2012:449-456.

        [26]Lacoste-julien S, Schmidt M, Bach F. A simpler approach to obtaining an O(1/t) convergence rate for the projected stochastic subgradient method[EB/ OL] . https://arxiv.org/abs/1212.2002, 2012-12-20[2012-12-20].

        [27]Schmidt M, Roux N L, Bach F. Minimizing finite sums with the stochastic average gradient[J]. Mathematical Programming, 2013,26(5):405-411.

        [28]Defazio A J, Caetano T S, Domke J. Finito: A faster, permutable incremental gradient method for big data problems[C]//Proceedings of the International Conference on Machine Leaning.Beijing:[s.n.], 2014:1125-1133.

        [29]Roux N L, Schmidt M, Bach F R. A stochastic gradient method with an exponential convergence _rate for finite training sets[C]//Advances in Neural Information Processing Systems. Lake Tahoe:[s.n.], 2012: 2663-2671.

        [30]Defazio A, Bach F, Lacoste-julien S. SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives[C]//Advances in Neural Information Processing Systems. Canada:[s.n.] ,2014:1646-1654.

        [31]Mairal J. Incremental majorization-minimization optimization with application to large-scale machine learning[J]. SIAM Journal on Optimization, 2015,25(2):829-855.

        [32]Johnson R, Zhang T. Accelerating stochastic gradient descent using predictive variance reduction[C]// Advances in Neural Information Processing Systems. Lake Tahoe: [s.n.], 2013:315-323.

        [33]Allen Zhu Z , Yuan Y. Univer: A universal variance reduction framework for proximal stochastic gradient method[EB/OL]. https://arxiv.org/abs/1506.01972, 2015-07-05[2015-07-05].

        [34]Atsushi N. Stochastic proximal gradient descent with acceleration techniques[C]//Advances in Neural Information Processing System.Montreal,Canada:[s.n.], 2014:1574-1582.

        [35]Lin H,Mairal J,Harchaoui Z. A universal catalyst for first-order optimization[C]//Advances in Neural Information Processing System.Montreal,Canada:[s.n.], 2015:1604-1616.

        [36]Shamir O. Open problem: Is averaging needed for strongly convex stochastic gradient descent[C]//Proceedings of the Conference on Learning Theory. USA:MIT Press,2012:471-475.

        [37]Shamir O, Zhang T. Stochastic gradient descent for non-smooth optimization: Convergence results and optimal averaging schemes[C]//Proceedings of the International Conference on Machine Learning. Atlanta,USA:[s.n.],2013:71-79.

        [38]Chen X, Lin Q, Pena J. Optimal regularized dual averaging methods for stochastic optimization[C]//Advances in Neural Information Processing Systems. Lake Tahoe: [s.n.], 2012, 404-412.

        [39]Nesterov Y, Shikhman V. Quasi-monotone subgradient methods for nonsmooth convex minimization[J]. J of Opti Theory and Appl, 2015, 165(3): 917-940.

        [40]陶蔚,潘志松,儲德軍,等.線性插值投影次梯度方法的最優(yōu)個(gè)體收斂速率[J].計(jì)算機(jī)研究與發(fā)展,2017,54(3):1-8.

        Tao Wei, Pan Zhisong, Chu Dejun, et al. The optimal individual convergence rate for the projected subgradient methods with linear interpolation operation[J]. Journal of Computer Research and Development, 2017,54(3):1-8.

        [41]Schmidt M, Le Roux N, Bach F. Convergence rates of inexact proximal-gradient methods for convex optimization[J]. Advances in Neural Information Processing Systems, 2011, 24:1458-1466.

        [42]Devolder O, Glineur F, Nesterov Y. First-order methods of smooth convex optimization with inexact oracle[J]. Mathematical Programming, 2014, 146(1/2):37-75.

        [43]Devolder O. Stochastic first order methods in smooth convex optimization[EB/OL]. CORE Discussion Papers, http://www.uclouvain.be/en-357992.html, 2011.

        [44]Honorio J. Convergence rates of biased stochastic optimization for learning sparse ising models[C] //Proceedings of the International Conference on Machine Learning. Edinburgh, Scotland:[s.n.], 2012:257-264.

        [45]Collobert R, Sinz F, Weston J, et al. Trading convexity for scalability[C]// Proceedings of the International Conference on Machine Learning. USA: MIT Press, 2006:201-208.

        [46]Wu Y, Lin Y. Robust truncated hinge loss support vector machines[J]. Journal of the American Statistical Association, 2007,102(479): 974-983.

        [47]Zhang T. Mutil-stage convex relaxation for learning with sparse regularization[C]// Advances in Neural Information Processing System. Vancouver, British Columbia, Canada:[s.n.],2008:1929-1936.

        [48]Lange K, Hunter D R, Yang I. Optimization transfer using surrogate objective functions[J]. J Comput Graph Statist, 2000(9):1-20.

        [49]Yuille A L, Rangarajan A. The concave-convex procedure[J]. Neural Computation, 2003,15(4):915-936.

        [50]Horst R, Thoai N V. DC programming: Overview[J]. J Optim Theory Appl, 1999,103:1-43.

        Individual Convergence of Stochastic Optimization Methods in Machine Learning

        Tao Qing1, Ma Po1, Zhang Menghan1, Tao Wei2

        (1.11st Department, Army Officer Academy of PLA, Hefei, 230031, China; 2.College of Command System, The PLA University of Science and Technology, Nanjing, 210007, China)

        The stochastic optimization algorithm is one of the state-of-the-art methods for solving large-scale machine learning problems, where the focus is on whether or not the optimal convergence rate is derived and the learning structure is ensured. So far, various kinds of stochastic optimization algorithms have been presented for solving the regularized loss problems. However, most of them only discuss the convergence in terms of the averaged output, and even the simplest sparsity cannot be preserved. In contrast to the averaged output, the individual solution can keep the sparsity very well, and its optimal convergence rate is extensively explored as an open problem. On the other hand, the commonly-used assumption about unbiased gradient in stochastic optimization often does not hold in practice. In such cases, an astonishing fact is that the bias in the convergence bound of accelerated algorithms will accumulate with the iteration, and this makes the accelerated algorithms inapplicable. In this paper, an overview of the state-of-the-art and existing problems about the stochastic first-order gradient methods is given, which includes the individual convergence rate, biased gradient and nonconvex problems. Based on it, some interesting problems for future research are indicated.

        machine learning; stochastic optimization; individual convergence; biased gradient estimation; non-convex problems

        國家自然科學(xué)基金(61673394,61273296)資助項(xiàng)目。

        2016-12-05;

        2017-01-10

        TP391

        A

        陶卿(1965-),男,教授,博士生導(dǎo)師,研究方向:機(jī)器學(xué)習(xí)、模式識別和應(yīng)用數(shù)學(xué),E-mail:qing.tao@ia.ac.cn, taoqing@gmail.com。

        馬坡(1993-),男,碩士研究生,研究方向:凸優(yōu)化及其在機(jī)器學(xué)習(xí)中的應(yīng)用。

        張夢晗(1994-),男,碩士研究生,研究方向:模式識別、人工智能。

        陶蔚(1991-),男,碩士研究生,研究方向:機(jī)器學(xué)習(xí)、網(wǎng)絡(luò)安全和模式識別。

        猜你喜歡
        優(yōu)化方法
        超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
        民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
        關(guān)于優(yōu)化消防安全告知承諾的一些思考
        一道優(yōu)化題的幾何解法
        由“形”啟“數(shù)”優(yōu)化運(yùn)算——以2021年解析幾何高考題為例
        學(xué)習(xí)方法
        可能是方法不對
        用對方法才能瘦
        Coco薇(2016年2期)2016-03-22 02:42:52
        四大方法 教你不再“坐以待病”!
        Coco薇(2015年1期)2015-08-13 02:47:34
        賺錢方法
        在线亚洲日本一区二区| 亚洲羞羞视频| 阿v视频在线| 国产在线观看黄片视频免费| 日韩精品免费在线视频| 久久老熟女一区二区三区福利 | 亚洲日韩中文字幕无码一区| 亚洲欧美日本| 国内精品久久久久国产盗摄| av一区二区三区观看| 天天做天天爱夜夜爽女人爽| 男人扒开女人下面狂躁小视频| 国产精品福利小视频| 吃下面吃胸在线看无码| 日本免费一区二区久久久| 亚洲人成网站色7799| 国产精品久久久久9999赢消| 美女视频一区| 国产一区二区三区白浆在线观看| 日本在线一区二区三区视频| 亚洲色偷偷偷综合网| 一本一道波多野结衣一区| 亚洲电影久久久久久久9999| 视频区一区二在线观看| 无码人妻aⅴ一区二区三区| 久久久久香蕉国产线看观看伊| 精品性高朝久久久久久久| 视频二区精品中文字幕| 护士人妻hd中文字幕| 97在线观看播放| 男人天堂免费视频| 国产三级在线观看不卡| 成年人一区二区三区在线观看视频 | 最新亚洲无码网站| 丰满老熟女性生活视频| 先锋影音人妻啪啪va资源网站| 色www永久免费视频| 98bb国产精品视频| 亚洲一区二区丝袜美腿| 色综合悠悠88久久久亚洲| 国产成人精品久久一区二区三区|