陳斯泳，安聰沛

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

對于 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

對于正則化算子RL，考慮以下兩類形式：

（1）第一類算子，考慮對正則化算子RL在 p∈PL上的作用進行定義，即

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

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

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

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

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

3 模型求解

先考慮求解模型（5），即RL=BL的情況?；谖墨I[10]中的定理2.1，有如下定理：

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

模型式（5）有唯一解

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

證明 式（6）證明見文獻[10]中的定理2.1。

下證式（7），對式（5）的目標函數(shù)中第一項進行展開，可得：

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

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

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

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

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

當(dāng)

綜合上面的分析可得：

定理得證。

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

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

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

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

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

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

對于子問題式（11），先令顯然，該問題是可分離的。那么，對于ζ的每個分量都有：

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

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

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

4 正則化算子

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

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

在該情況下，BL為對角元全為1的矩陣。此時，問題式（5）可以劃歸為Lasso問題[14]。

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

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

①三角多項式過濾函數(shù)[10]：

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

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

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

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

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

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

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

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

5 數(shù)值實驗

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

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

該函數(shù)為C∞(S2)中的函數(shù)。第二個函數(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ù)設(shè)定為

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

為檢驗逼近原始函數(shù)的效果，實驗中選取一致性誤差和L2誤差作為標準。

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3 誤差隨λ變化趨勢圖

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

6 結(jié)束語

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

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

[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.