求解一般l1趨勢過濾問題的原始對偶內點法

張體琪，劉勇進

(福州大學數學與統計學院，福建 福州 350108)

0 引言

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

(1)

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

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

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

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

1 原始對偶內點法

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

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

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

(2)

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

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

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

1.1 迭代方向的求解

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

(4)

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

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可以由如下的等式來計算

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

1.2 原始對偶內點法的算法框架

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

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

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

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

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

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

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

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

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

(7)

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

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

(8)

此性質的證明在文獻[13]中的節10.3.3有詳細的證明，故此處略去證明過程.

(9)

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

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

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

故算法二次收斂. 證畢.

2 數值實驗與結果分析

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

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

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

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

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

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

算法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)

交替方向乘子法也是求解凸優化問題的一種重要的方法，由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)可以轉化為求解下面的等價形式

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

(14)

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

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

2.2 停止條件

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

2.3 在合成數據集上的實驗結果

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

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

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

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

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

2.4 在真實數據集上的實驗結果

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

表2 測試數據集介紹

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

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

3 結語

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