球面l1-正則化逼近模型及其應(yīng)用研究

陳斯泳，安聰沛

暨南大學(xué) 信息科學(xué)技術(shù)學(xué)院 數(shù)學(xué)系，廣州 510632

1 引言

球面上的逼近問(wèn)題具有廣泛的應(yīng)用背景，例如測(cè)量地球表面的大氣氣流[1]、晶體結(jié)構(gòu)[2]、靜電學(xué)逼近[3]、球面幾何的聲波仿真[1]，熱輻射傳感[4]，宇宙中圖像的恢復(fù)[5]和醫(yī)療圖像重建[6]都是該問(wèn)題的源泉。在數(shù)學(xué)上，上述實(shí)際問(wèn)題可以抽象為球面上的最優(yōu)化問(wèn)題[3]、偏微分方程[5]、積分方程[7]、樣條逼近[1]、鞍點(diǎn)問(wèn)題[8]與最小二乘問(wèn)題[9]等。

在二維球面S2上，文獻(xiàn)[10]提出了在連續(xù)和離散情況下帶有旋轉(zhuǎn)不變性的正則化球面最小二乘多項(xiàng)式逼近，涵蓋了球面多項(xiàng)式插值，超插值和過(guò)濾超插值等一系列最小二乘模型[10]。

本文主要研究一類(lèi)在二維單位球面S2:={x=(x,y,z)T∈?3|x2+y2+z2=1}上帶l1-正則項(xiàng)的最小二乘逼近模型：

其中，f是給定的連續(xù)函數(shù)，且其在有N個(gè)點(diǎn)的點(diǎn)集XN={x1,x2,…,xN}?S2上的值是給定的（可能帶有噪聲）；PL:=PL(S2)是?3中所有限定在球面上且次數(shù)小于等于L的多項(xiàng)式所組成的線性空間；正則化算子RL為線性算子，λ>0為正則化參數(shù)。該模型根據(jù)節(jié)點(diǎn)XN和正則化算子RL的選取，可以變換出多種不同的形式。

本文選取適當(dāng)次數(shù)t的球面t-設(shè)計(jì)作為節(jié)點(diǎn)，其定義如下：

定義1[2]稱(chēng)球面S2上的點(diǎn)集XN={x1,x2,…, }xN?S2為球面t-設(shè)計(jì)，如果其滿足于對(duì)球面上所有次數(shù)不超過(guò)t的多項(xiàng)式在XN上的值的算術(shù)平均值準(zhǔn)確等于該多項(xiàng)式在球面上積分的幾何平均值，即XN滿足：

其中，dω(x)是單位球面上的表面測(cè)度。顯然，得到一個(gè)球面t-設(shè)計(jì)，就相當(dāng)于得到一個(gè)對(duì)次數(shù)不超過(guò)t的多項(xiàng)式準(zhǔn)確成立的數(shù)值積分公式。在本文中，假定節(jié)點(diǎn)XN是t≥2L且節(jié)點(diǎn)數(shù)滿足N=(t+1)2的球面t-設(shè)計(jì)。

2 基于球面t-設(shè)計(jì)的帶l1-正則項(xiàng)最小二乘逼近模型

為了簡(jiǎn)化原模型式（1），選取球面調(diào)和函數(shù)

作為多項(xiàng)式空間PL的一組正交基[11]。對(duì)該正交基進(jìn)行規(guī)范化，使得那么：

其維數(shù)為次數(shù)?的球面調(diào)和函數(shù)Y?,k組成次數(shù)為?的齊次調(diào)和多項(xiàng)式空間H?，其維數(shù)為2?+1。其正交性是關(guān)于L2內(nèi)積正交：

其導(dǎo)出范數(shù)為于是，對(duì)任意函數(shù) p∈PL，存在唯一的向量α=(α?,k)∈?(L+1)2使得：

給定連續(xù)函數(shù) f及點(diǎn)集XN={x1,x2,…,xN}，令f:=f(XN)為列向量 f=[f (x1),f(x2),…,f(xN)]T∈?N。

對(duì)于 L≥1，球面調(diào)和矩陣YL:=YL(XN)∈?(L+1)2×N的元素為：

Y?,k(xj),?=0,1,…,L,k=1,2,…,2?+1;j=1,2,…,N

對(duì)于正則化算子RL，考慮以下兩類(lèi)形式：

（1）第一類(lèi)算子，考慮對(duì)正則化算子RL在 p∈PL上的作用進(jìn)行定義，即

其中，第二行公式由加法定理[11]推導(dǎo)得到，即對(duì)于球面調(diào)和函數(shù)，有：

其中，x?y表示?3中x和y的歐氏內(nèi)積，P?表示次數(shù)為?的Legendre多項(xiàng)式，且經(jīng)規(guī)范化后滿足P?(1)=1。此時(shí)：

其中，矩陣 BL∈?(L+1)2×(L+1)2為任意非負(fù)數(shù) β0,β1,…,βL組成的半正定對(duì)角矩陣：

那么，式（1）可以轉(zhuǎn)化為如下最小二乘模型：

（2）第二類(lèi)算子，考慮式（4）的特殊情況：正則化算子直接作用在系數(shù)α上，即=。此時(shí)，式（4）轉(zhuǎn)化為如下最小二乘模型：

3 模型求解

先考慮求解模型（5），即RL=BL的情況。基于文獻(xiàn)[10]中的定理2.1，有如下定理：

定理1給定L≥0，令XN={x1,x2,…,xN}?S2為球面t-設(shè)計(jì)且t≥2L。那么：

模型式（5）有唯一解

其中為軟閾值算子，其定義為：

證明 式（6）證明見(jiàn)文獻(xiàn)[10]中的定理2.1。

下證式（7），對(duì)式（5）的目標(biāo)函數(shù)中第一項(xiàng)進(jìn)行展開(kāi)，可得：

由HL是非奇異的可知，問(wèn)題式（5）是嚴(yán)格凸的，那么它存在唯一解。再由可微及式（5）的一階最優(yōu)性條件可得，它的唯一解滿足：

其中，?(?)表示次梯度[12]。由及 BL是對(duì)角陣可得，該問(wèn)題是可分離的，故α是問(wèn)題的解，當(dāng)且僅當(dāng)它的每個(gè)分量α?,k滿足：

顯然，記為問(wèn)題的最優(yōu)解，則

當(dāng)2γ?,k>λβ?,k時(shí)那么

當(dāng) 2γ?,k<-λβ?,k時(shí) ，，那 么 ，此時(shí)

當(dāng)

綜合上面的分析可得：

定理得證。

下面考慮模型式（4），即RL=BLYL的情況。顯然，這是一個(gè)l2-l1優(yōu)化問(wèn)題，本文使用交替方向法（Alternating Direction Method of Multipliers，ADMM）[13]求解該優(yōu)化問(wèn)題。首先，令ζ=RTLα，模型式（4）就轉(zhuǎn)化為如下帶約束優(yōu)化問(wèn)題：

式（9）的增廣拉格朗日算子為：

其中，y為L(zhǎng)agrange乘子，ρ為罰參數(shù)。那么，ADMM包含如下迭代：

其中在每步迭代中，都需要求解兩個(gè)子問(wèn)題式（10）、（11），下面逐個(gè)進(jìn)行分析。

對(duì)于子問(wèn)題式（10），由其一階最優(yōu)性條件可知它的解αi+1滿足如下線性方程組：

求解該線性方程組，可得：

對(duì)于子問(wèn)題式（11），先令顯然，該問(wèn)題是可分離的。那么，對(duì)于ζ的每個(gè)分量都有：

其中，第一項(xiàng)不可微。依據(jù)次梯度的理論[12]與定理1的證明，可得該問(wèn)題解的每個(gè)分量為：

其中，Sk(a)為軟閾值算子式（8）。

綜合上面兩點(diǎn)分析，運(yùn)用ADMM算法求解模型，就是先將式（4）轉(zhuǎn)化為帶約束優(yōu)化問(wèn)題式（9），再將其分解為兩個(gè)子問(wèn)題，然后分別通過(guò)求解線性方程組與軟閾值算子直接得到子問(wèn)題的解，最后迭代求解兩個(gè)子問(wèn)題得到式（4）的解，即其迭代形式如下：

4 正則化算子

正則化算子RL由對(duì)角矩陣 BL的對(duì)角元 β?決定。不同的正則化算子，往往有不同的特性。選擇合適的算子，往往能在特定的方面，如去噪、恢復(fù)等方面取得優(yōu)異的效果。下面介紹三個(gè)具有旋轉(zhuǎn)不變性[6]的正則化算子。

（1）單位算子，其定義如下：

在該情況下，BL為對(duì)角元全為1的矩陣。此時(shí)，問(wèn)題式（5）可以劃歸為L(zhǎng)asso問(wèn)題[14]。

（2）過(guò)濾算子，該算子對(duì)應(yīng)的對(duì)角元β?的定義為：

其中，本文考慮如下過(guò)濾函數(shù)h(x)：

①三角多項(xiàng)式過(guò)濾函數(shù)[10]：

在式（14）中，由于當(dāng)?=L時(shí)，βL=∞，進(jìn)而αL,k=0，故排除?=L的情況。運(yùn)用該類(lèi)算子可以產(chǎn)生過(guò)濾多項(xiàng)式逼近。

（3）微分算子，Laplace-Beltrami算子 Δ*[11]的定義為：

球面調(diào)和函數(shù)是Laplace-Beltrami算子的特征函數(shù)，具有如下性質(zhì)：

Δ*Y?,k(x)=-?(?+1)Y?,k(x)

又由-Δ*為一個(gè)半正定算子[11]，對(duì)于任意s>0，定義(-Δ*)s為：

(-Δ*)sY?,k(x)=[?(?+1)]sY?,k(x)

本文采用(-Δ*)s作為正則化算子。此時(shí)，對(duì)應(yīng)的矩陣BL為：

運(yùn)用該算子能夠?qū)г肼暤暮瘮?shù)圖像進(jìn)行恢復(fù)[1]。

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

本章主要展示相關(guān)數(shù)值實(shí)驗(yàn)的結(jié)果，以展現(xiàn)l2-l1模型在多項(xiàng)式逼近和函數(shù)圖像恢復(fù)中的作用。

本文選取好條件球面t-設(shè)計(jì)作為節(jié)點(diǎn)，并選取如下兩個(gè)測(cè)試函數(shù)。第一個(gè)函數(shù)是Franke函數(shù)[16]：

該函數(shù)為C∞(S2)中的函數(shù)。第二個(gè)函數(shù)為Franke函數(shù) f1與球冠函數(shù) fcap[17]之和：

f2=f1+fcap

其中

C(xc,r):={x ∈S2|arccos(x,xc)≤r}是以xc為中心，r為半徑的球冠，η為一正常數(shù)。函數(shù) f2在球面上連續(xù)但在球冠C(xc,r)邊緣處不可微。在數(shù)值實(shí)驗(yàn)中，參數(shù)設(shè)定為

圖1 誤差隨多項(xiàng)式次數(shù)L變化趨勢(shì)圖

為檢驗(yàn)逼近原始函數(shù)的效果，實(shí)驗(yàn)中選取一致性誤差和L2誤差作為標(biāo)準(zhǔn)。

（1）一致性誤差的定義由下式給定：

其中，XN為球面上一大規(guī)模的分布良好的點(diǎn)集。這里選取N=50 000的等區(qū)域節(jié)點(diǎn)[18]作為XN。對(duì)于 f2，增加在球冠邊緣處布點(diǎn)。

（2）L2誤差的定義由下式給定：

其中{x1,x2,…,xm}為次數(shù)t=160，m=25 921的好條件球面t-設(shè)計(jì)點(diǎn)集[19]。

5.1 過(guò)濾超插值算子應(yīng)用于精確數(shù)據(jù)

本節(jié)展示應(yīng)用過(guò)濾超插值算子于精確數(shù)據(jù)的數(shù)值實(shí)驗(yàn)，并比較不同的過(guò)濾函數(shù)。測(cè)試函數(shù)為 f1和 f2。對(duì)于給定多項(xiàng)式次數(shù)L，考慮節(jié)點(diǎn)次數(shù)t=2L及節(jié)點(diǎn)數(shù)N=(t+1)2。 β?由式（14）給定，過(guò)濾函數(shù)為 h1(x)和h2(x)，正則化參數(shù)設(shè)定為λ=1。

圖1給出了多項(xiàng)式次數(shù)L=1,2,…,40時(shí)運(yùn)用模型式（4）、（5）逼近函數(shù) f1和 f2的誤差趨勢(shì)圖。從圖中可以看出，與h1(x)相比，過(guò)濾函數(shù)為h2(x)的過(guò)濾超插值算子具有較小的誤差，模型（5）與模型（4）誤差較為接近。

圖2 微分算子應(yīng)用于非光滑函數(shù)圖像

5.2 微分算子應(yīng)用于帶噪聲數(shù)據(jù)

本節(jié)展示重建帶噪聲的非光滑函數(shù)的數(shù)值實(shí)驗(yàn)。該實(shí)驗(yàn)使用模型（5），并選取微分算子(s=2)作為正則化算子和不同取值的λ作為正則化參數(shù)。

圖2（a）展示了函數(shù) f2的圖像，圖2（b）展示了帶噪聲的函數(shù)fδ2(x)=f2+δ(x)的圖像，其中，對(duì)于每一個(gè)x，δ(x)為服從均值μ=0、標(biāo)準(zhǔn)差σ=0.2的正態(tài)分布隨機(jī)變量的一個(gè)樣本。在恢復(fù)帶噪聲的數(shù)據(jù)時(shí)，模型中正則化參數(shù)λ的選取是至關(guān)重要的。故先考察λ對(duì)誤差的影響，以選取最優(yōu)的λ。

實(shí)驗(yàn)選取給定多項(xiàng)式次數(shù)L=25及對(duì)應(yīng)的節(jié)點(diǎn)次數(shù)t=2L=50，λ從10-15,10-14.5,…,104.5,105變化。圖3給出一致性誤差和L2誤差變化的曲線。從圖3中可以看出，隨著λ的增大，一致性誤差和L2誤差皆呈“不變—減小—增大—不變”的趨勢(shì)，顯然存在最優(yōu)的λ，使得誤差最小。在后續(xù)的實(shí)驗(yàn)中，選取使得一致性誤差最小的λ作為最優(yōu)的λ。

接著，使用次數(shù)t=50，N=2 601的好條件球面t-設(shè)計(jì)恢復(fù)帶噪聲的數(shù)據(jù)。這里選取帶l2-正則項(xiàng)的最小二乘模型[10]作為比較：

圖3 誤差隨λ變化趨勢(shì)圖

圖2（c）～（f）展示了兩個(gè)模型的恢復(fù)效果和誤差。從圖2（d）、（f）可以看出，l2-l1模型和 l2-l2模型均能較好地恢復(fù)出原函數(shù)的圖像，l2-l2模型整體恢復(fù)效果較好，但球冠函數(shù)和Franke函數(shù)交界處（不光滑）的誤差較大，而l2-l1模型則能較好地處理交界處不光滑的現(xiàn)象，即交界處的誤差較小，但整個(gè)球面上存在多處誤差較大的區(qū)域。最后，圖4給出誤差隨多項(xiàng)式次數(shù)變化的趨勢(shì)圖。

6 結(jié)束語(yǔ)

圖4 誤差隨多項(xiàng)式次數(shù)變化圖

本文研究球面上的多項(xiàng)式逼近問(wèn)題，針對(duì)不同的正則化算子的選取，建立了一類(lèi)球面上帶l1-正則項(xiàng)最小二乘逼近模型。通過(guò)選取好條件球面t-設(shè)計(jì)點(diǎn)作為采樣點(diǎn)，文中運(yùn)用了直接法和交替方向法求解此逼近問(wèn)題。最后，將此逼近模型應(yīng)用到球面上函數(shù)逼近—精確數(shù)據(jù)和噪聲污染的情形。測(cè)試結(jié)果表明，l2-l1模型能較好地逼近兩種情形下的光滑和非光滑球面函數(shù)，特別是在邊緣處有較好優(yōu)勢(shì)。

[1]Freeden W，Gervens T，Schreiner M.Constructive approximation on the sphere with applications to geomathematics[M].[S.l.]：Oxford University Press on Demand，1998.

[2]Delsarte P，Goethals J M，Seidel J J.Spherical codes and designs[J].Geometriae Dedicata，1977，6：363-388.

[3]Saff E B，Kuijlaars A B J.Distributing many points on a sphere[J].The Mathematical Intelligencer，1997，19（1）：5-11.

[4]Grella K，Schwab C.Sparse discrete ordinates method in radiativetransfer[J].ComputationalMethodsin Applied Mathematics，2011，11（3）：305-326.

[5]Le Gia Q T，Sloan I，Wendland H.Multiscale analysis in sobolev spaces on the sphere[J].SIAM Journal on Numerical Analysis，2010，48（6）：2065-2090.

[6]Reimer M.Multivariate polynomial approximation[M].Basel：Birkh?user，2003：144.

[7]Atkinson K E.The numerical solution of integral equations of the second kind[M].[S.l.]：Cambridge University Press，1997.

[8]Le Gia Q T，Sloan I H，Wathen A J.Stability and preconditioning for a hybrid approximation on the sphere[J].Numerische Mathematik，2011，118（4）：695-711.

[9]Le Gia Q T，Narcowich F J，Ward J D.Continuous and discrete least-squares approximation by radial basis functions on spheres[J].Journal of Approximation Theory，2006，143（1）：124-133.

[10]An C，Chen X，Sloan I H.Regularized least squares approximations on the sphere using spherical designs[J].SIAM Journal on Numerical Analysis，2012，50：1513-1534.

[11]Müller C.Spherical harmonics[M].Berlin，Heidelberg：Springer，1966：17.

[12]Boyd S，Vandenberghe L.Convex optimization[M].[S.l.]：Cambridge University Press，2004.

[13]Boyd S，Parikh N，Chu E.Distributed optimization and statistical learning via the alternating direction method of multipliers[J].Found Trends Mach Learn，2011，3（1）：1-122.

[14]Tibshirani R.Regression shrinkage and selection via the lasso[J].Journal of the Royal Statistical Society Series B：Methodological，1996：267-288.

[15]Filbir F，Themistoclakis W.Polynomial approximation on the sphere using scattered data[J].Mathematische Nachrichten，2008，281（5）：650-668.

[16]Renka R J.Multivariate interpolation of large sets of scattered data[J].ACM Transactions on Mathematical Software（TOMS），1988，14：139-148.

[17]Williamson D L，Drake J B，Hack J J.A standard test set for numerical approximations to the shallow water equations in spherical geometry[J].Journal of Computational Physics，1992，102：211-224.

[18]Leopardi P.Diameter bounds for equal area partitions of the unit sphere[J].Electronic Transactions on Numerical Analysis，2009，35：1-16.

[19]An C，Chen X，Sloan I H.Well conditioned spherical designs for integration and interpolation on the twosphere[J].SIAM Journal on Numerical Analysis，2010，48：2135-2157.