球面l1-正则化逼近模型及其应用研究

陈斯泳，安聪沛

暨南大学 信息科学技术学院 数学系，广州 510632

1 引言

球面上的逼近问题具有广泛的应用背景，例如测量地球表面的大气气流[1]、晶体结构[2]、静电学逼近[3]、球面几何的声波仿真[1]，热辐射传感[4]，宇宙中图像的恢复[5]和医疗图像重建[6]都是该问题的源泉。在数学上，上述实际问题可以抽象为球面上的最优化问题[3]、偏微分方程[5]、积分方程[7]、样条逼近[1]、鞍点问题[8]与最小二乘问题[9]等。

在二维球面S2上，文献[10]提出了在连续和离散情况下带有旋转不变性的正则化球面最小二乘多项式逼近，涵盖了球面多项式插值，超插值和过滤超插值等一系列最小二乘模型[10]。

本文主要研究一类在二维单位球面S2:={x=(x,y,z)T∈ℝ3|x2+y2+z2=1}上带l1-正则项的最小二乘逼近模型：

其中，f是给定的连续函数，且其在有N个点的点集XN={x1,x2,…,xN}⊂S2上的值是给定的（可能带有噪声）；PL:=PL(S2)是ℝ3中所有限定在球面上且次数小于等于L的多项式所组成的线性空间；正则化算子RL为线性算子，λ>0为正则化参数。该模型根据节点XN和正则化算子RL的选取，可以变换出多种不同的形式。

本文选取适当次数t的球面t-设计作为节点，其定义如下：

定义1[2]称球面S2上的点集XN={x1,x2,…, }xN⊂S2为球面t-设计，如果其满足于对球面上所有次数不超过t的多项式在XN上的值的算术平均值准确等于该多项式在球面上积分的几何平均值，即XN满足：

其中，dω(x)是单位球面上的表面测度。显然，得到一个球面t-设计，就相当于得到一个对次数不超过t的多项式准确成立的数值积分公式。在本文中，假定节点XN是t≥2L且节点数满足N=(t+1)2的球面t-设计。

2 基于球面t-设计的带l1-正则项最小二乘逼近模型

为了简化原模型式（1），选取球面调和函数

作为多项式空间PL的一组正交基[11]。对该正交基进行规范化，使得那么：

其维数为次数ℓ的球面调和函数Yℓ,k组成次数为ℓ的齐次调和多项式空间Hℓ，其维数为2ℓ+1。其正交性是关于L2内积正交：

其导出范数为于是，对任意函数 p∈PL，存在唯一的向量α=(αℓ,k)∈ℝ(L+1)2使得：

给定连续函数 f及点集XN={x1,x2,…,xN}，令f:=f(XN)为列向量 f=[f (x1),f(x2),…,f(xN)]T∈ℝN。

对于 L≥1，球面调和矩阵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]推导得到，即对于球面调和函数，有：

其中，x⋅y表示ℝ3中x和y的欧氏内积，Pℓ表示次数为ℓ的Legendre多项式，且经规范化后满足Pℓ(1)=1。此时：

其中，矩阵 BL∈ℝ(L+1)2×(L+1)2为任意非负数 β0,β1,…,βL组成的半正定对角矩阵：

那么，式（1）可以转化为如下最小二乘模型：

（2）第二类算子，考虑式（4）的特殊情况：正则化算子直接作用在系数α上，即=。此时，式（4）转化为如下最小二乘模型：

3 模型求解

先考虑求解模型（5），即RL=BL的情况。基于文献[10]中的定理2.1，有如下定理：

定理1给定L≥0，令XN={x1,x2,…,xN}⊂S2为球面t-设计且t≥2L。那么：

模型式（5）有唯一解

其中为软阈值算子，其定义为：

证明 式（6）证明见文献[10]中的定理2.1。

下证式（7），对式（5）的目标函数中第一项进行展开，可得：

由HL是非奇异的可知，问题式（5）是严格凸的，那么它存在唯一解。再由可微及式（5）的一阶最优性条件可得，它的唯一解满足：

其中，∂(⋅)表示次梯度[12]。由及 BL是对角阵可得，该问题是可分离的，故α是问题的解，当且仅当它的每个分量αℓ,k满足：

显然，记为问题的最优解，则

当2γℓ,k>λβℓ,k时那么

当 2γℓ,k<-λβℓ,k时 ，，那 么 ，此时

当

综合上面的分析可得：

定理得证。

下面考虑模型式（4），即RL=BLYL的情况。显然，这是一个l2-l1优化问题，本文使用交替方向法（Alternating Direction Method of Multipliers，ADMM）[13]求解该优化问题。首先，令ζ=RTLα，模型式（4）就转化为如下带约束优化问题：

式（9）的增广拉格朗日算子为：

其中，y为Lagrange乘子，ρ为罚参数。那么，ADMM包含如下迭代：

其中在每步迭代中，都需要求解两个子问题式（10）、（11），下面逐个进行分析。

对于子问题式（10），由其一阶最优性条件可知它的解αi+1满足如下线性方程组：

求解该线性方程组，可得：

对于子问题式（11），先令显然，该问题是可分离的。那么，对于ζ的每个分量都有：

其中，第一项不可微。依据次梯度的理论[12]与定理1的证明，可得该问题解的每个分量为：

其中，Sk(a)为软阈值算子式（8）。

综合上面两点分析，运用ADMM算法求解模型，就是先将式（4）转化为带约束优化问题式（9），再将其分解为两个子问题，然后分别通过求解线性方程组与软阈值算子直接得到子问题的解，最后迭代求解两个子问题得到式（4）的解，即其迭代形式如下：

4 正则化算子

正则化算子RL由对角矩阵 BL的对角元 βℓ决定。不同的正则化算子，往往有不同的特性。选择合适的算子，往往能在特定的方面，如去噪、恢复等方面取得优异的效果。下面介绍三个具有旋转不变性[6]的正则化算子。

（1）单位算子，其定义如下：

在该情况下，BL为对角元全为1的矩阵。此时，问题式（5）可以划归为Lasso问题[14]。

（2）过滤算子，该算子对应的对角元βℓ的定义为：

其中，本文考虑如下过滤函数h(x)：

①三角多项式过滤函数[10]：

在式（14）中，由于当ℓ=L时，βL=∞，进而αL,k=0，故排除ℓ=L的情况。运用该类算子可以产生过滤多项式逼近。

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

球面调和函数是Laplace-Beltrami算子的特征函数，具有如下性质：

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

又由-Δ*为一个半正定算子[11]，对于任意s>0，定义(-Δ*)s为：

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

本文采用(-Δ*)s作为正则化算子。此时，对应的矩阵BL为：

运用该算子能够对带噪声的函数图像进行恢复[1]。

5 数值实验

本章主要展示相关数值实验的结果，以展现l2-l1模型在多项式逼近和函数图像恢复中的作用。

本文选取好条件球面t-设计作为节点，并选取如下两个测试函数。第一个函数是Franke函数[16]：

该函数为C∞(S2)中的函数。第二个函数为Franke函数 f1与球冠函数 fcap[17]之和：

f2=f1+fcap

其中

C(xc,r):={x ∈S2|arccos(x,xc)≤r}是以xc为中心，r为半径的球冠，η为一正常数。函数 f2在球面上连续但在球冠C(xc,r)边缘处不可微。在数值实验中，参数设定为

图1 误差随多项式次数L变化趋势图

为检验逼近原始函数的效果，实验中选取一致性误差和L2误差作为标准。

（1）一致性误差的定义由下式给定：

其中，XN为球面上一大规模的分布良好的点集。这里选取N=50 000的等区域节点[18]作为XN。对于 f2，增加在球冠边缘处布点。

（2）L2误差的定义由下式给定：

其中{x1,x2,…,xm}为次数t=160，m=25 921的好条件球面t-设计点集[19]。

5.1 过滤超插值算子应用于精确数据

本节展示应用过滤超插值算子于精确数据的数值实验，并比较不同的过滤函数。测试函数为 f1和 f2。对于给定多项式次数L，考虑节点次数t=2L及节点数N=(t+1)2。 βℓ由式（14）给定，过滤函数为 h1(x)和h2(x)，正则化参数设定为λ=1。

图1给出了多项式次数L=1,2,…,40时运用模型式（4）、（5）逼近函数 f1和 f2的误差趋势图。从图中可以看出，与h1(x)相比，过滤函数为h2(x)的过滤超插值算子具有较小的误差，模型（5）与模型（4）误差较为接近。

图2 微分算子应用于非光滑函数图像

5.2 微分算子应用于带噪声数据

本节展示重建带噪声的非光滑函数的数值实验。该实验使用模型（5），并选取微分算子(s=2)作为正则化算子和不同取值的λ作为正则化参数。

图2（a）展示了函数 f2的图像，图2（b）展示了带噪声的函数fδ2(x)=f2+δ(x)的图像，其中，对于每一个x，δ(x)为服从均值μ=0、标准差σ=0.2的正态分布随机变量的一个样本。在恢复带噪声的数据时，模型中正则化参数λ的选取是至关重要的。故先考察λ对误差的影响，以选取最优的λ。

实验选取给定多项式次数L=25及对应的节点次数t=2L=50，λ从10-15,10-14.5,…,104.5,105变化。图3给出一致性误差和L2误差变化的曲线。从图3中可以看出，随着λ的增大，一致性误差和L2误差皆呈“不变—减小—增大—不变”的趋势，显然存在最优的λ，使得误差最小。在后续的实验中，选取使得一致性误差最小的λ作为最优的λ。

接着，使用次数t=50，N=2 601的好条件球面t-设计恢复带噪声的数据。这里选取带l2-正则项的最小二乘模型[10]作为比较：

图3 误差随λ变化趋势图

图2（c）～（f）展示了两个模型的恢复效果和误差。从图2（d）、（f）可以看出，l2-l1模型和 l2-l2模型均能较好地恢复出原函数的图像，l2-l2模型整体恢复效果较好，但球冠函数和Franke函数交界处（不光滑）的误差较大，而l2-l1模型则能较好地处理交界处不光滑的现象，即交界处的误差较小，但整个球面上存在多处误差较大的区域。最后，图4给出误差随多项式次数变化的趋势图。

6 结束语

图4 误差随多项式次数变化图

本文研究球面上的多项式逼近问题，针对不同的正则化算子的选取，建立了一类球面上带l1-正则项最小二乘逼近模型。通过选取好条件球面t-设计点作为采样点，文中运用了直接法和交替方向法求解此逼近问题。最后，将此逼近模型应用到球面上函数逼近—精确数据和噪声污染的情形。测试结果表明，l2-l1模型能较好地逼近两种情形下的光滑和非光滑球面函数，特别是在边缘处有较好优势。

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