鐘軼君， 李崇君

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

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

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

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

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

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

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

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

1 預(yù)備知識(shí)

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

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

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

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

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

1.2 符號(hào)說明

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

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

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

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

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

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

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

1.3 逆尺度空間算法

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

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

輸出：x(tk)。

步驟1.初始化

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

在Sk=Ik上的解。

步驟3.更新時(shí)間變量為

步驟4.更新次梯度為

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

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

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

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

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

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

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

1.4 分片稀疏

對(duì)于給定的向量

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

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

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

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

分片逆尺度空間系統(tǒng)(piecewise inverse scale space，P_ISS)[22]旨在恢復(fù)分片稀疏向量中每一個(gè)子向量∈Rni的支集，即P_ISS方法分別且同時(shí)求解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.并行時(shí)間路徑更新。

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

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

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

滿足不等式

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

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

2.1 P_ISS算法中刪除機(jī)制的缺陷

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

2.2 分片符號(hào)一致性

本文考慮對(duì)算法進(jìn)行分片符號(hào)一致性的改進(jìn)，即當(dāng)某一“片”到達(dá)符號(hào)一致性時(shí), 停止此“片”的支集更新。從而保證在最小二乘計(jì)算中，該“片”不會(huì)隨著晚到達(dá)符號(hào)一致性的“片”的更新而增加冗余的元素。該算法自行地對(duì)每一“片”實(shí)行分片符號(hào)一致性的“監(jiān)督”，使得算法自適應(yīng)性更強(qiá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.并行時(shí)間路徑更新。

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

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

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

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

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

其中，

證明：在每個(gè)子矩陣Ai中，選擇指標(biāo)集為Si的列向量組成的子矩陣ASi?？紤]Grassmannian矩陣

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

證畢

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

其中，

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

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

3 數(shù)值實(shí)驗(yàn)

3.1 實(shí)驗(yàn)設(shè)計(jì)

本文的散亂數(shù)據(jù)曲面擬合問題，對(duì)應(yīng)于分片稀疏模型，矩陣A= [A1, … ,AN]是由定義的觀測矩陣帶噪音的觀測數(shù)據(jù){fi}是用向量bσ表示。此部分通過給定的散亂點(diǎn)集{xk,yk}和其相對(duì)應(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)榱艘粋€(gè)求解分片稀疏向量(曲面的系數(shù)向量)，即

的分片稀疏逼近問題。

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

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

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

(2) 帶有刪除機(jī)制的分片逆尺度空間算法(P_ISS)，其中刪除機(jī)制中參數(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)、分片稀疏度及算法運(yùn)行時(shí)間(CPU)進(jìn)行4種算法比較，RMS為

其中，

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

3.2 實(shí)驗(yàn)結(jié)果分析

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

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

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

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

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

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

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

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

4 結(jié) 束 語

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

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