張體琪，劉勇進

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

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

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

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

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

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)是一個關(guān)于對偶變量v的二次規(guī)劃(QP)問題，如果有-λ1

1.1 迭代方向的求解

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

(4)

其中:f1v=v-λ1；f2v=-v-λ1； 擾動參數(shù)t>0； diag(x)表示對角元素為x的對角矩陣；rt(v,μ1,μ2)表示殘差;0表示零矩陣.當t→∞時，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的雅可比矩陣.再進一步將此牛頓方程改寫為以下的形式

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

(5)

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

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

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

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

算法1: 原始對偶內(nèi)點法(PDIP)輸入: t>0, η∈(0, 0.5], ξ∈(0, 1), (v0, μ01, μ02)∈Ω, j=01. 計算非線性系統(tǒng) rt(vj, μj1, μj2)=0的牛頓步(Δvj, Δμj1, Δμj2)2. 令αj=ξmj, 其中mj為滿足下面不等式的最小非負整數(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 算法的收斂性和復雜度分析

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

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

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

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

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

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

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

(7)

只有零解. 由于取值為k時結(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ì)的證明在文獻[13]中的節(jié)10.3.3有詳細的證明，故此處略去證明過程.

(9)

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

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

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

故算法二次收斂. 證畢.

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

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

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

半光滑牛頓增廣拉格朗日方法是統(tǒng)計優(yōu)化中一類重要的方法，可以應(yīng)用于解決一些大規(guī)模的統(tǒng)計問題，包括Lasso問題[14]、SLBoxLSR問題[15]、OSCAR問題和SLOPE 問題[16]. 下面簡要介紹求解一般l1趨勢過濾問題的半光滑牛頓增廣拉格朗日方法. 對一個閉且正常凸的函數(shù)f:→(-∞, ∞]及一個常數(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. 計算: xj+1=argminx∈n{Φj(x):=infz Lσj(x, z; vj)}(10)2. 計算:zj+1=proxσ-1jλ·1(D(k, n)xj+1+σ-1jvj)3. 計算: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. 計算:xj+1∈argminx∈n Lσ(x, zj; vj)(12)2. 計算:zj+1∈argminz∈n-k Lσ(xj+1, z; vj)(13)3. 計算:vj+1=vj+ωσ(D(k, n)xj+1-zj+1)4. 令j=j+1, 返回步驟1.

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

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

(14)

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

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

2.2 停止條件

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)語

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