鐘軼君， 李崇君

(1. 浙江理工大學(xué)理學(xué)院數(shù)學(xué)科學(xué)系，浙江 杭州 310018；2. 大連理工大學(xué)數(shù)學(xué)與科學(xué)學(xué)院，遼寧 大連 116024)

散亂數(shù)據(jù)擬合(函數(shù)擬合/逼近)是一個在信號處理、計算機圖形學(xué)中被廣泛研究的問題。其通過給定的散亂的、帶噪音的點集來擬合一個較好的曲線或曲面，需要解決2個基本問題：①如何選擇一個簡單的函數(shù)空間作為描述對象函數(shù)的逼近空間；②其空間要有很好的逼近度。為此，隨著散亂點曲面擬合工具多樣性的不斷探索，有學(xué)者開始利用探索多層結(jié)構(gòu)性質(zhì)來改進(jìn)逼近方法，如LEE等[1]提出的基于B樣條的多層逼近方法；CASTANO等[2]提出了基于小波的多層正則化方法等。JOHNSON等[3]進(jìn)一步探討了在平移不變空間(principal shift invariant，PSI)中基于 B樣條和小波的散亂點曲面擬合問題。以上方法均力求在所構(gòu)建的函數(shù)類中與研究對象函數(shù)本身誤差最小。對于第②個問題，學(xué)者們希望在其空間中得到散亂數(shù)據(jù)的稀疏表示解，因此推動了稀疏逼近和優(yōu)化方法在函數(shù)擬合領(lǐng)域的廣泛應(yīng)用。隨著深層挖掘曲面中的稀疏性，幾何建模中很多的啟發(fā)式方法得到了進(jìn)一步發(fā)展和推廣。近些年，信號的稀疏性研究也取得了豐碩的成果，其中細(xì)化研究信號的稀疏分布(非零元素的分布)成為了熱點，如HUANG[4]提出了結(jié)構(gòu)化稀疏的概念，即一類非零元具有一定結(jié)構(gòu)性的稀疏向量，如廣泛研究的非零元是成塊出現(xiàn)的塊稀疏[5-9]，在統(tǒng)計中也稱為Group Lasso[10-17]，除此之外還有樹稀疏[18]、圖稀疏[19]。劉翼鵬[19]討論了向量服從全局稀疏、局部分布不同的凸優(yōu)化模型，討論了全局稀疏局部稠密，即非零元素呈簇狀分布的稀疏向量。提取出了可以描述整個信號的結(jié)構(gòu)性信息，有利于使用更少的樣本和運算量來恢復(fù)信號，且能夠決定信號是否可以有效地得到恢復(fù)。對于不同的實際場景，信號的稀疏表示向量具有不同的結(jié)構(gòu)，更好地研究這些結(jié)構(gòu)有利于理解和解決實際問題。隨著對稀疏性的細(xì)化探索，學(xué)者們嘗試在曲面重構(gòu)，尤其是在PSI空間中的曲面重構(gòu)問題所具有的結(jié)構(gòu)型稀疏表示進(jìn)行研究。HAO 等[20]提出了利用凸優(yōu)化模型來考慮PSI空間中基于B樣條的散亂點曲面擬合問題，利用多塊LASSO (MLASSO)得到散亂點曲面擬合問題的稀疏解；文獻(xiàn)[21]得到了曲面在PSI空間中每一層 B樣條上的稀疏表示解?；贛LASSO的散亂點曲面擬合問題不僅可以得到稀疏解(稀疏表示系數(shù))而且擬合效果較好，但其多個正則化參數(shù)的選取是一個待解決的問題。OSHER等[22]提出 Bregman逆尺度空間算法(inverse scale space，ISS)，即

此方法是貪婪算法和凸優(yōu)化算法的一種結(jié)合，不僅規(guī)避了凸優(yōu)化算法中參數(shù)的選取問題而且運行速率較快。之后LI和ZHONG[23]提出了一種帶有刪除機制的分片逆尺度空間算法，核心是將分片稀疏性引入到ISS算法中。在稀疏恢復(fù)問題上，相較于ISS算法，該算法能夠更好地保護(hù)信號的小尺度非零元素，在應(yīng)用于散亂數(shù)據(jù)擬合問題中得到了較好的擬合效果。分片稀疏性是一類描述向量中非零元素散亂分布現(xiàn)象的稀疏性。DONOHO和HUO[24]針對片稀疏信號對應(yīng)觀測矩陣具有的分塊結(jié)構(gòu)，引入分塊矩陣互相干(mutual coherence)條件，可將分塊正交矩陣的精確恢復(fù)信號的稀疏度理論上界推廣到一般分塊矩陣之中, 改進(jìn)了稀疏信號可精確恢復(fù)的充分性條件。其結(jié)果不僅給出了分片稀疏恢復(fù)的理論保證，也提升了L0模型和壓縮感知恢復(fù)整體稀疏信號的可信度。向量的分片稀疏結(jié)構(gòu)所對應(yīng)的觀測矩陣A的分塊結(jié)構(gòu)恰好符合PSI空間基于 B樣條所生成基函數(shù)的層數(shù)，因此引入分片稀疏性到稀疏優(yōu)化算法(模型)中，并分析散亂點曲面擬合問題是有意義的。

然而分片ISS算法(P_ISS)[23]也有其缺陷，算法中的刪除機制增加了新的參數(shù)選取問題。目前考慮基于散亂點的曲面擬合的稀疏表示的文獻(xiàn)較少?；诖耍疚奶岢鲆环N自適應(yīng)的分片逆尺度空間算法(adaptive piecewise ISS，aP_ISS)，由此得到散亂點曲面擬合問題的稀疏表示。此算法不僅規(guī)避了帶有刪除機制的分片逆尺度算法中需要人工選取參數(shù)的問題，而且通過對“分片符號一致性”的分析，使得算法具有自適應(yīng)性，可更好地應(yīng)用到散亂點曲面逼近問題中。

本文主要貢獻(xiàn)如下：

(1) 通過將分片稀疏的理念引入到散亂數(shù)據(jù)擬合問題中，拓廣了稀疏技術(shù)在函數(shù)逼近中的應(yīng)用。

(2) 改進(jìn)了文獻(xiàn)[22]中的分片逆尺度空間算法，規(guī)避了算法中的刪除機制。利用“分片符號一致性”使得算法更具自適應(yīng)性。

(3) 稀疏算法作為高效算法已經(jīng)成功地應(yīng)用到了曲面擬合問題中，本文建立了一個稀疏表示系數(shù)的先驗信息，從而可以更好擬合曲面的細(xì)節(jié)信息，達(dá)到更好地逼近效果。與已有的稀疏重建算法相比，不僅有更松弛的理論保證并且擬合效果更好。

1 預(yù)備知識

1.1 散亂數(shù)據(jù)曲面擬合問題

典型的散亂點曲面擬合問題是指：存在一個散亂點集合X={x,x,… ,x} ?Ω?Rd和相對應(yīng)的函數(shù)值f|X= {f1,f2, …,fn}，且需要找到一個屬于空間H的函數(shù)g能夠擬合給定的數(shù)據(jù) { (xk,fk) }nk=1。經(jīng)典的光滑樣條散亂點重構(gòu)模型為

由于g屬于Beppo-Levi空間，所以當(dāng)數(shù)據(jù)點規(guī)模變大時計算量也會變大。為了解決該問題，文獻(xiàn)[3]考慮在一個由緊致支撐函數(shù)移位擴張生成的空間(稱為 PSI空間)中求解以上正則化問題。利用PSI空間不僅避免了計算量增大的問題，也由于其和小波、B樣條的聯(lián)系使得稀疏逼近數(shù)據(jù)函數(shù)成為可能。

利用分塊矩陣A= [A1, … ,AN]表示由定義的觀測矩陣

向量bσ表示被噪音污染的數(shù)據(jù){fi}，則曲面擬合問題g(xi) ≈fi可以記為bσ=Ax+e，其中x為曲面表示系數(shù)；e為噪音。

1.2 符號說明

定義1[24].觀測矩陣A的列互相干為

其中，ai為矩陣A的第i列。

定義 2[23].設(shè)矩陣A= [A1, … ,AN]是N個矩陣的拼接，定義第i塊的子陣互相干為

命題1[23].第i塊的子陣互相干μi,i滿足

其中，αi∈[0,1]為衡量第i個子矩陣塊內(nèi)互相干與矩陣整體互相干比例關(guān)系的參數(shù)。

Oracle估計[21]是最小二乘解在真實支撐集S上的部分，oracle的非零元素可表達(dá)為

分片oracle估計[22]是oracle估計限制在子空間Si上的部分

1.3 逆尺度空間算法

文獻(xiàn)[25]中提出了一個自適應(yīng)的逆尺度空間算法(adaptive inverse scale space，aISS)，實現(xiàn)算法步驟如下：

輸入：A，bσ,threshold ≥ 0 。

輸出：x(tk)。

步驟1.初始化

步驟2. 計算x(tk)為問題

在Sk=Ik上的解。

步驟3.更新時間變量為

步驟4.更新次梯度為

其中，t=tk+1，ejn∈R為第j個分量為1其余分量為0的向量。

步驟5.更新指標(biāo)集為

文獻(xiàn)[21]中給出了如下ISS系統(tǒng)(1)的理論保證條件。

假設(shè) 1[21].限制強凸性(restricted strong convexity)：存在γ∈ ( 0,1]，使得

假設(shè) 2[21].不可表示條件或精確恢復(fù)條件(irre-presentable condition or exact recovery condition)：存 在η∈ ( 0,1)， 使 得其中

定理1[21]. 設(shè)假設(shè)1和2均成立，那么存在一個時刻τ使得Bregman ISS路徑(1)在時刻τ之前沒有錯誤選擇，即?t≤τ，supp(x(t)) ? supp(x)以大于的概率成立。若

的概率有sign(x(τ))=sign(x)，即Bregman ISS求解的最優(yōu)解達(dá)到了符號一致性。

1.4 分片稀疏

對于給定的向量

(1) 全局稀疏。x稱為s-稀疏的若x至多含有個非零元素；

(2) 塊稀疏[4-9]。x稱為塊s-稀疏的若x至多含有s個含有非零元素的塊，即塊l0范數(shù)2,0||x||=不超過s；

(3) 分片稀疏[23]。向量可以分成N個部分，且每部分xiT∈Rni都是一個稀疏向量，若已知si=||xi||0(i=1,…,N)，則稱x是分片(s1,… ,sN)-稀疏向量。

1.5 分片逆尺度空間算法(P_ISS)

分片逆尺度空間系統(tǒng)(piecewise inverse scale space，P_ISS)[22]旨在恢復(fù)分片稀疏向量中每一個子向量∈Rni的支集，即P_ISS方法分別且同時求解xi，再由所有的xi共同構(gòu)成x，即

算法1[22].分片逆尺度空間算法

輸入：觀測矩陣A= [A1, … ,AN]，觀測向量bσ，噪音程度參數(shù)σ。

輸出：分片稀疏向量x=(x1T, … ,xTN)T∈Rn。

步驟1.初始化。

步驟2.最小二乘。更新為

的解。其中I為指標(biāo)集；IP為在指標(biāo)集I上的投影算子。

步驟3.并行時間路徑更新。

步驟4.并行的次梯度向量更新。

步驟5.并行的指標(biāo)集更新。

步驟6.刪除機制。設(shè)置向量

滿足不等式

的分量為0。引入分片稀疏性后提出的帶刪除機制的分片Bregman逆尺度空間(P_ISS)算法。因該算法結(jié)構(gòu)類似于正交匹配算法，所以其運行速率較快，又因求解的是分片凸優(yōu)化問題，故其收斂到l1稀疏解。數(shù)值實驗表明相比于 ISS算法，分片Bregman逆尺度空間算法能夠更好地保護(hù)小尺度元素不被噪音混淆。

2 自適應(yīng)分片逆尺度空間算法(aP_ISS)

2.1 P_ISS算法中刪除機制的缺陷

在改良的貪婪算法中，比如壓縮采樣匹配追蹤算法 (compressive sampling MP，CoSaMP)[26]，利用刪除機制對算法當(dāng)前恢復(fù)的解進(jìn)行自適應(yīng)的修正是較為常用且有效的方法。CoSaMP在算法最后一步保留最大的s(s=||x||0)個非零元素，其余的元素均設(shè)置為0以克服類似StOMP算法[27]在支集添加冗余元素的缺陷。在算法P_ISS迭代過程中，由于分片稀疏向量x整體到達(dá)符號一致性[22]的時刻其中iτ為第i片到達(dá)符號一致性的時刻。即所有的“片”到達(dá)了符號一致性算法才會停止更新支集。在所有的N個子向量都到達(dá)符號一致性之前，早到達(dá)符號一致性的“片”會由于不斷的迭代產(chǎn)生錯誤或冗余的非零元素，所以在P_ISS算法中設(shè)置刪除機制可以有效地避免冗余元素的產(chǎn)生。但是，刪除機制中引入的參數(shù)C使得算法增加了選取問題，即參數(shù)C與數(shù)據(jù)密切相關(guān)，需要經(jīng)驗性選取，消耗大量時間。

2.2 分片符號一致性

本文考慮對算法進(jìn)行分片符號一致性的改進(jìn)，即當(dāng)某一“片”到達(dá)符號一致性時, 停止此“片”的支集更新。從而保證在最小二乘計算中，該“片”不會隨著晚到達(dá)符號一致性的“片”的更新而增加冗余的元素。該算法自行地對每一“片”實行分片符號一致性的“監(jiān)督”，使得算法自適應(yīng)性更強。

2.3 自適應(yīng)分片逆尺度空間算法(aP_ISS)

算法2.自適應(yīng)分片逆尺度空間算法

輸入：觀測矩陣A= [A1, … ,AN],觀測向量bσ，噪音程度參數(shù)σ。

輸出：

步驟1.初始化。

步驟2.最小二乘。更新為

的解。其中，I為指標(biāo)集；PI為在指標(biāo)集I上的投影算子。

步驟3.并行時間路徑更新。

步驟4.并行的次梯度向量更新。

步驟5.并行的自適應(yīng)的指標(biāo)集更新。若第i片未達(dá)到符號一致性，則

實際算法操作中，通過算法當(dāng)前迭代向量和oracle向量的符號一致性，即考慮第i片的向量xi和分片oracle估計向量是否達(dá)到了符號一致性

并衡量是否更新指標(biāo)集。下面給出 aP_ISS的性能保證定理。首先給出一些相關(guān)的引理。

引理1.矩陣AS= [AS1,…,ASN]的奇異值Λ滿足

其中，

證明：在每個子矩陣Ai中，選擇指標(biāo)集為Si的列向量組成的子矩陣ASi。考慮Grassmannian矩陣

其中，Φij為ΦS中的元素，且ΦS中對角線元素滿足非對角線元素滿足|Φij| ≤μ。因此

證畢

引理 2.設(shè)矩陣A= [A1, … ,AN]的每一個子矩陣Ai的互相干參數(shù)是αi，假設(shè)1成立，則

其中，

定理2.設(shè)假設(shè)條件1和2均成立，若分片稀疏信號滿足

其中，smax為si=|Si|的最大值；

3 數(shù)值實驗

3.1 實驗設(shè)計

本文的散亂數(shù)據(jù)曲面擬合問題，對應(yīng)于分片稀疏模型，矩陣A= [A1, … ,AN]是由定義的觀測矩陣帶噪音的觀測數(shù)據(jù){fi}是用向量bσ表示。此部分通過給定的散亂點集{xk,yk}和其相對應(yīng)的觀測數(shù)據(jù){f(xk,yk)}擬合曲面并得到其稀疏表示g(x,y)。其逼近函數(shù)是由張量積 B樣條基函數(shù)(x,y)(i= 1 ,…,N;j= 1,…,ni)表示，其中N為多層B樣條的層數(shù)，基函數(shù)的總數(shù)為n=n1+…+nN，于是待擬合曲面為

考慮到N層B樣條基函數(shù)之間是相關(guān)的，于是此問題轉(zhuǎn)變?yōu)榱艘粋€求解分片稀疏向量(曲面的系數(shù)向量)，即

的分片稀疏逼近問題。

本文通過對一個函數(shù)f添加高斯噪音來人工生成一個帶噪音的數(shù)據(jù)集{((xk,yk) ,f(xk,yk)):k=1,… ,9 00}，即

對比以下算法求解散亂數(shù)據(jù)曲面擬合問題的稀疏表示的數(shù)值效果：

(1) 逆尺度空間算法(aISS)。

(2) 帶有刪除機制的分片逆尺度空間算法(P_ISS)，其中刪除機制中參數(shù)C=0.3。

(3) 多塊LASSO (MLASSO)利用Jacobi型的逼近交替方向乘子算法(alternating direction method of multipilers，ADMM)[28]求解MLASSO，其中正則化參數(shù)的值取為γ1=0.3，γ2= 0 .15，γ3=0.2，γ4=0.3。

(4) 本文提出的自適應(yīng)分片逆尺度空間算法(aP_ISS)。通過均方根誤差(root-mean-square，RMS)、分片稀疏度及算法運行時間(CPU)進(jìn)行4種算法比較，RMS為

其中，

3個噪音分別取值為σ1= 0 .25,σ2= 0 .1,σ3= 0 .01。

3.2 實驗結(jié)果分析

從表1～3的數(shù)值結(jié)果和圖1～3的擬合效果可以發(fā)現(xiàn)，若不考慮分片稀疏性的a_ISS算法，本文算法(aP_ISS)較其他 2種算法在近似散亂點集上優(yōu)勢較明顯，同時在保護(hù)分片稀疏度上效果更好，即每一層支集均不為空集且稀疏。相較于分片稀疏性的MLASSO算法，aP_ISS算法在運行效率上有明顯的優(yōu)勢，MLASSO模型存在事先人工輸入和調(diào)節(jié)多個正則化參數(shù)的缺陷，本文算法的自適應(yīng)性保證了算法可以更快地運行，且 MLASSO算法近似的曲面效果遠(yuǎn)不如本文算法。P_ISS算法的主要缺陷在于刪除機制中參數(shù)C的選取依賴數(shù)據(jù)、噪音等因素，需要花費大量時間，而本文算法不僅避免了參數(shù)的選取而且能夠得到和 P_ISS算法同等的近似效果。

另從表1～3的結(jié)果還可以觀察到，分片稀疏逼近的目的是為了在每一層上都有非零元素表示，且是稀疏表示，本文算法除了近一步改進(jìn)已有的分片算法 P_ISS，更重要的是為了通過“分片符號一致性”來體現(xiàn)分片稀疏的思想。實驗結(jié)果表明本文的算法不僅在擬合效果上更好，并且達(dá)到了分片稀疏表示的效果。

表1 觀測函數(shù)f1對應(yīng)曲面擬合的數(shù)值結(jié)果

表2 觀測函數(shù)f2對應(yīng)曲面擬合的數(shù)值結(jié)果

表3 觀測函數(shù)f3對應(yīng)曲面擬合的數(shù)值結(jié)果

圖1 觀測函數(shù)f1對應(yīng)4種算法擬合出的近似曲面的對比圖

圖2 觀測函數(shù)f2對應(yīng)4種算法擬合出的近似曲面對比圖

圖3 觀測函數(shù)f3對應(yīng)4種算法擬合出的近似曲面對比圖

4 結(jié) 束 語

本文利用 aP_ISS處理散亂點曲面擬合問題,通過引入分片稀疏性，可以有效地保護(hù)在PSI空間中多層支集的分片稀疏性；進(jìn)一步通過對分片逆尺度空間系統(tǒng)的分片符號一致性的考慮，避免了分片算法中參數(shù)的選取問題。實驗結(jié)果表明，該算法近似散亂點曲面效果較好，驗證了利用稀疏優(yōu)化算法處理散亂點曲面擬合問題的優(yōu)勢。

分片稀疏思想應(yīng)用到散亂數(shù)據(jù)擬合等問題仍然欠缺較為完善的研究，通過對稀疏優(yōu)化算法的自適應(yīng)分片處理來求解曲面擬合問題，將是未來研究的重點。