張?bào)w琪，劉勇進(jìn)

(福州大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，福建 福州 350108)

0 引言

對于給定的噪聲信號y=(y1, …,yn)∈n和正整數(shù)k, 一般l1趨勢過濾問題具有以下的形式：

(1)

D(k, n)=D(1, n-k+1)D(k-1, n)(k≥2)

其中D(1, p)∈(p-1)×p為p上的一階差分矩陣.

非參數(shù)回歸是一個(gè)熱門的研究領(lǐng)域，學(xué)者們已經(jīng)提出了眾多的方法. 趨勢過濾是一種重要的非參數(shù)回歸模型，它是一種帶懲罰項(xiàng)的最小二乘擬合優(yōu)化方法. 趨勢過濾的應(yīng)用非常廣泛，包括宏觀經(jīng)濟(jì)學(xué)[1]、社會科學(xué)[2]、金融時(shí)間序列分析[3]、收益管理[4]、地球物理學(xué)[5]、噪聲處理[6-7]等領(lǐng)域. 為了解決一些實(shí)際的經(jīng)濟(jì)問題，許多具有線性時(shí)間復(fù)雜度的趨勢過濾方法已經(jīng)被提出，比如，指數(shù)平滑方法[8]、平均濾波方法[9]和Hodrick-Prescott(H-P)濾波方法[10]. Kim等[11]在2009年提出了l1趨勢過濾方法，它基于H-P濾波方法做了一些變化，即用絕對值的和(l1范數(shù))代替H-P濾波方法中用于懲罰估計(jì)趨勢變化的平方和(l2范數(shù)).但他們只解決了k=2的情況，因此，將這種思想推廣到更一般的情況，提出一種能夠解決一般l1趨勢過濾問題(k≥2)的原始對偶內(nèi)點(diǎn)法.

本研究針對一般l1趨勢過濾問題提出一種原始對偶內(nèi)點(diǎn)法.首先，介紹一些相關(guān)的預(yù)備知識，分析原始對偶內(nèi)點(diǎn)法中求解迭代方向的過程，進(jìn)而給出算法的迭代框架.然后，給出原始對偶內(nèi)點(diǎn)法的收斂性分析以及算法復(fù)雜度分析.最后，在合成數(shù)據(jù)集和真實(shí)數(shù)據(jù)集上進(jìn)行相關(guān)的實(shí)驗(yàn)，展示原始對偶內(nèi)點(diǎn)法與半光滑牛頓增廣拉格朗日方法、交替方向乘子法解決一般l1趨勢過濾問題的數(shù)值對比結(jié)果.實(shí)驗(yàn)結(jié)果表明, 對于不同的調(diào)節(jié)參數(shù)λ和k階差分矩陣，原始對偶內(nèi)點(diǎn)法都具有較好的性能.

1 原始對偶內(nèi)點(diǎn)法

原始對偶內(nèi)點(diǎn)法是目前求解二次規(guī)劃或其他凸優(yōu)化規(guī)劃問題最有效的求解方法之一. 該方法是Gonzalez-Lima等[12]求解線性規(guī)劃問題時(shí)，利用原始對偶方法產(chǎn)生的一個(gè)線性系統(tǒng)的簡化. 這種簡化的主要特征是它在解集上有定義并且保持了稀疏性，這些特性增加了算法的魯棒性、穩(wěn)定性和效率.

首先給出一些預(yù)備知識. 先將原問題(1)重構(gòu)為下面的等價(jià)形式

s.t.D(k, n)x-z=0

(2)

原問題(2)的拉格朗日函數(shù)為

其中，v∈n-k為對偶變量. 原問題(2)的KKT最優(yōu)化條件為

原問題(2)的對偶問題為

s.t.-λ1≤v≤λ1

(3)

D(k, n)D(k, n)Tv-D(k, n)y+μ1-μ2=0,v-λ1≤0, -v-λ1≤0

〈μ1,v-λ1〉=0, 〈μ2, -v-λ1〉=0 (μi≥0,i=1, 2, …,n)

對偶問題(3)是一個(gè)關(guān)于對偶變量v的二次規(guī)劃(QP)問題，如果有-λ1

1.1 迭代方向的求解

本節(jié)給出原始對偶內(nèi)點(diǎn)法中求解迭代方向的具體過程[11]，為了推導(dǎo)方便，本節(jié)將D(k, n)全部記為D.定義如下關(guān)于殘差的非線性系統(tǒng).即

(4)

其中:f1v=v-λ1；f2v=-v-λ1； 擾動(dòng)參數(shù)t>0； diag(x)表示對角元素為x的對角矩陣；rt(v,μ1,μ2)表示殘差;0表示零矩陣.當(dāng)t→∞時(shí)，rt(v,μ1,μ2)→0退化為對偶問題(3)的KKT最優(yōu)化條件.

rt(v+Δv,μ1+Δμ1,μ2+Δμ2)≈rt(v,μ1,μ2)+Drt(v,μ1,μ2)(Δv, Δμ1, Δμ2)T=0

其中:Drt為rt的雅可比矩陣.再進(jìn)一步將此牛頓方程改寫為以下的形式

Drt(v,μ1,μ2)(Δv, Δμ1, Δμ2)T=-rt(v,μ1,μ2)

(5)

方程中的Δμ1和Δμ2可以由如下的等式來計(jì)算

最后，將得到的非線性系統(tǒng)(4)的牛頓步Δw=(Δv, Δμ1, Δμ2)作為原始對偶內(nèi)點(diǎn)法所尋找的迭代方向.

1.2 原始對偶內(nèi)點(diǎn)法的算法框架

明確了牛頓迭代方向的求解過程后，下面給出原始對偶內(nèi)點(diǎn)法的算法迭代框架.

算法1: 原始對偶內(nèi)點(diǎn)法(PDIP)輸入: t>0, η∈(0, 0.5], ξ∈(0, 1), (v0, μ01, μ02)∈Ω, j=01. 計(jì)算非線性系統(tǒng) rt(vj, μj1, μj2)=0的牛頓步(Δvj, Δμj1, Δμj2)2. 令αj=ξmj, 其中mj為滿足下面不等式的最小非負(fù)整數(shù)m, 使得: rt(vj+ξmΔvj, μj1+ξmΔμj1, μj2+ξmΔμj2)≤(1-ηξm)rt(vj, μj1, μj2)(6)且(vj+ξmΔvj, μj1+ξmΔμj1, μj2+ξmΔμj2)∈Ω成立. 3. 更新 vj+1=vj+αjΔvj, μj+11=μj1+αjΔμj1, μj+12=μj2+αjΔμj24. 令j=j+1, 返回步驟1.

1.3 算法的收斂性和復(fù)雜度分析

原始對偶內(nèi)點(diǎn)法具有全局收斂性和局部收斂性且收斂速度為二階收斂，其算法復(fù)雜度為O(k2n). 首先給出一些結(jié)論：

根據(jù)文獻(xiàn)[13]的節(jié)10.2.4知道，結(jié)論2)成立等價(jià)于證明下面的性質(zhì).

性質(zhì)1對任意給定的正整數(shù)k,n, 矩陣D(k, n)D(k, n)T正定.

證明 矩陣D(k, n)D(k, n)T正定等價(jià)于證明D(k, n)Tx=0只有零解，使用數(shù)學(xué)歸納法來證明D(k, n)Tx=0只有零解.

當(dāng)k取值為2時(shí)，易知矩陣D(2, n)T為列滿秩矩陣，故齊次線性方程組D(2, n)Tx=0只有零解.

假設(shè)取值為k時(shí)結(jié)論成立，下面證明當(dāng)k+1時(shí)結(jié)論也成立，即證明齊次線性方程組

D(k+1, n)Tu=D(k, n)TD(1, n-k)Tu=0

(7)

只有零解. 由于取值為k時(shí)結(jié)論成立，即D(k, n)T列滿秩，根據(jù)式(7)得到D(1, n-k)Tu=0.由于D(1, n-k)T列滿秩，說明齊次方程組D(1, n-k)Tu=0只有零解，因此齊次方程組D(k+1, n)Tu=0只有零解.證畢.

性質(zhì)2對任意w∈Ω, 0≤α≤min{1,αmax}, Δw為在w處的牛頓方向，有下面的不等式成立.

(8)

此性質(zhì)的證明在文獻(xiàn)[13]中的節(jié)10.3.3有詳細(xì)的證明，故此處略去證明過程.

(9)

這說明α=1滿足回溯線搜索停止條件(6)，故算法最后的步長為1.

接下來給出算法的收斂性定理.

證明 當(dāng)j充分大時(shí)，αj=1, 有wj+1=wj+Δwj.對最優(yōu)解w*有rt(w*)=0, 可以得到

故算法二次收斂. 證畢.

2 數(shù)值實(shí)驗(yàn)與結(jié)果分析

在合成數(shù)據(jù)集和真實(shí)數(shù)據(jù)集上分別進(jìn)行數(shù)值實(shí)驗(yàn). 將原始對偶內(nèi)點(diǎn)法與半光滑牛頓增廣拉格朗日方法、交替方向乘子法進(jìn)行對比，以進(jìn)一步說明原始對偶內(nèi)點(diǎn)法在解決一般l1趨勢過濾問題的較好性能. 所有的數(shù)值測試都在MATLAB-R2019a中進(jìn)行，運(yùn)行環(huán)境為Dell臺式電腦(Intel(R) Core(TM) i5-8500 CPU @3.0 GHZ，8 GB RAM).

2.1 半光滑牛頓增廣拉格朗日方法與交替方向乘子法

半光滑牛頓增廣拉格朗日方法是統(tǒng)計(jì)優(yōu)化中一類重要的方法，可以應(yīng)用于解決一些大規(guī)模的統(tǒng)計(jì)問題，包括Lasso問題[14]、SLBoxLSR問題[15]、OSCAR問題和SLOPE 問題[16]. 下面簡要介紹求解一般l1趨勢過濾問題的半光滑牛頓增廣拉格朗日方法. 對一個(gè)閉且正常凸的函數(shù)f:→(-∞, ∞]及一個(gè)常數(shù)γ>0，γf的臨近映射[17]具有如下的形式:

給定κ>0，1范數(shù)的臨近映射，也被稱為軟閾值算子，具有以下的形式:

其中sign表示符號函數(shù)，1n表示分量全為1的n維向量.給定σ>0，定義原問題(2)的增廣拉格朗日函數(shù)為

接下來給出半光滑牛頓增廣拉格朗日法的算法框架.

算法2: 半光滑牛頓增廣拉格朗日方法(SSNAL)輸入: σ0>0, λ≥0, (x0, z0; v0)∈n×n-k×n-k, j=0.1. 計(jì)算: xj+1=argminx∈n{Φj(x):=infz Lσj(x, z; vj)}(10)2. 計(jì)算:zj+1=proxσ-1jλ·1(D(k, n)xj+1+σ-1jvj)3. 計(jì)算:vj+1=vj+σj(D(k, n)xj+1-zj+1)4. 更新σj+1↑σ∞≤+∞, 令j=j+1, 返回步驟1.

(11)

交替方向乘子法也是求解凸優(yōu)化問題的一種重要的方法，由Gabay等[18]在1976年首次提出. 交替方向乘子法求解一般l1趨勢過濾問題的算法框架如下:

算法3: 交替方向乘子法(ADMM)輸入: ω∈0, 1+52(), σ>0, (x0, z0;v0)∈n×n-k×n-k, j=0.1. 計(jì)算:xj+1∈argminx∈n Lσ(x, zj; vj)(12)2. 計(jì)算:zj+1∈argminz∈n-k Lσ(xj+1, z; vj)(13)3. 計(jì)算:vj+1=vj+ωσ(D(k, n)xj+1-zj+1)4. 令j=j+1, 返回步驟1.

求解子問題(12)可以轉(zhuǎn)化為求解下面的等價(jià)形式

(In+σD(k, n)D(k, n)T)xj+1=y+D(k, n)T(σzj-vj)

(14)

可以進(jìn)一步應(yīng)用共軛梯度法或者Cholesky分解法來得到線性系統(tǒng)(14)的解. 同時(shí)注意到子問題(13)具有閉式解，其形式如下:

由于求解一般l1趨勢過濾問題的交替方向乘子法的關(guān)鍵也在于子問題(12)的求解，根據(jù)其等價(jià)形式(14)，D(k, n)D(k, n)T的計(jì)算成本為O(n(n-k)2)，故求解子問題(12)的等價(jià)形式(14)所需的總時(shí)間復(fù)雜度為O(n3).

2.2 停止條件

使用基于原問題(2)的KKT條件的相對殘差來測量所有算法獲得的近似最優(yōu)解的精度.

2.3 在合成數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果

在合成數(shù)據(jù)集上對比原始對偶內(nèi)點(diǎn)法和半光滑牛頓增廣拉格朗日法、交替方向乘子法解決一般l1趨勢過濾問題的數(shù)值結(jié)果，下面簡要介紹合成數(shù)據(jù)集的構(gòu)造過程.

yt=xt+zt,t=1, …，n；xt+1=xt+vt,t=1, …,n-1

其中：vt為趨勢斜率；xt為“真實(shí)”的潛在趨勢；zt表示不規(guī)則分量或噪聲.初始條件為x1=0, 噪聲zt服從獨(dú)立同分布N(0,β2).趨勢斜率vt從一個(gè)簡單的馬爾可夫過程中選取得到.當(dāng)概率為p的時(shí)候，有vt+1=vt, 即潛在趨勢沒有斜率變化.當(dāng)概率為1-p時(shí)，從均勻分布[-b,b]上選取vt+1.同樣地，從均勻分布[-b,b]上選取初始斜率v1.在實(shí)驗(yàn)中，設(shè)置p=0.01,β=1,b=0.5.選取參數(shù)λ=1, 100, 10 000;k=2, 4; 數(shù)據(jù)維數(shù)n=i×105,i=1, …, 5.設(shè)定算法的最長運(yùn)行時(shí)間為3 h，即10 800 s. 表格中黑色橫線表示算法在該參數(shù)下所用時(shí)間超出設(shè)定的最長運(yùn)行時(shí)間10 800 s. 由于篇幅有限，本文只展示k=4的結(jié)果.

表1展示了當(dāng)k=4、tol=1×10-7時(shí), PDIP、SSNAL、ADMM算法在合成數(shù)據(jù)集上解決一般l1趨勢過濾問題的數(shù)值對比結(jié)果，可以看到PDIP算法不管是在運(yùn)行時(shí)間還是在求解精度上都仍然保持著明顯的優(yōu)勢，并且參數(shù)λ取得越大，PDIP算法的優(yōu)勢更明顯. 隨著λ逐漸增大，SSNAL和ADMM算法已經(jīng)達(dá)不到求解的精度并且還需要更長的運(yùn)行時(shí)間. 故總結(jié)出：PDIP算法相比于SSNAL和ADMM算法在合成數(shù)據(jù)集上解決一般l1趨勢過濾問題時(shí)更加高效和穩(wěn)健.

表1 當(dāng)k=4、 tol=1×10-7時(shí)，PDIP、SSNAL、ADMM算法在合成數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果

2.4 在真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果

在真實(shí)數(shù)據(jù)集上比較PDIP和SSNAL、ADMM算法解決一般l1趨勢過濾問題的數(shù)值結(jié)果. 真實(shí)數(shù)據(jù)來自PJM數(shù)據(jù)庫. 在PJM數(shù)據(jù)庫中選取了4個(gè)數(shù)據(jù)集作為實(shí)驗(yàn)的測試集. 實(shí)驗(yàn)選取調(diào)節(jié)參數(shù)λ=1, 100, 10 000. 數(shù)據(jù)的維數(shù)依賴于所選取的數(shù)據(jù)集的大小. 表2給出測試數(shù)據(jù)集的介紹. 表3展示了當(dāng)k=4、tol=1×10-7時(shí)PDIP、SSNAL、ADMM算法在真實(shí)數(shù)據(jù)集上解決一般l1趨勢過濾問題的數(shù)值結(jié)果.

表2 測試數(shù)據(jù)集介紹

表3 當(dāng)k=4、 tol=1×10-7時(shí)，PDIP、SSNAL、ADMM算法在真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果

通過表3數(shù)據(jù)結(jié)果可以看出，PDIP算法的優(yōu)勢依舊明顯. 在λ=10 000,n=803 000時(shí), PDIP算法仍具有出色的性能表現(xiàn)，僅僅花費(fèi)約40.81 s就求解到了滿足精度的解； 但SSNAL和ADMM算法的性能表現(xiàn)比較差，不僅求解精度不高而且求解時(shí)間也很長. 因此可以得出結(jié)論：PDIP算法在真實(shí)數(shù)據(jù)集上解決一般l1趨勢過濾問題時(shí)的性能表現(xiàn)明顯優(yōu)于SSNAL和ADMM算法的性能表現(xiàn).

3 結(jié)語

為求解一般l1趨勢過濾問題，提出一種原始對偶內(nèi)點(diǎn)法，給出原始對偶內(nèi)點(diǎn)法的算法框架并對算法進(jìn)行收斂性分析和時(shí)間復(fù)雜度分析，最后在合成數(shù)據(jù)集和真實(shí)數(shù)據(jù)集上進(jìn)行數(shù)值實(shí)驗(yàn)，通過將原始對偶內(nèi)點(diǎn)法與半光滑牛頓增廣拉格朗日方法和交替方向乘子法比較，進(jìn)一步證明原始對偶內(nèi)點(diǎn)法的高效性和穩(wěn)健性.