# Continuous Regularized Least Squares Polynomial Approximation on the Sphere

Authors: Yang Zhou; Yanan Kong
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we consider the problem of polynomial reconstruction of smooth functions on the sphere from their noisy values at discrete nodes on the two-sphere. The method considered in this paper is a weighted least squares form with a continuous regularization. Preliminary error bounds in terms of regularization parameter, noise scale, and smoothness are proposed under two assumptions: the mesh norm of the data point set and the perturbation bound of the weight. Condition numbers of the linear systems derived by the problem are discussed. We also show that sphericaltϵ-designs, which can be seen as a generalization of spherical t-designs, are well applied to this model. Numerical results show that the method has good performance in view of both the computation time and the approximation quality.


---

## Body

## 1. Introduction

In recent decades, the problem of finding a polynomial approximationg:S2⟶ℝ with degree of g up to an integer L of a continuous function f∈CS2 on the sphere S2:=x∈ℝ3:x2=1 has received much attention [1–6]. Given noisy data fδxi at points xi∈S2, i=1,…,N, regularized least squares schemes are usually considered, whereas the regularization varies for different tasks. In this paper, we consider approximating the function f with a continuous regularized least squares scheme (CRLS) as(1)ming∈ℙL∑i=1Nwigxi−fδxi2+λgℍs2,where wi, i=1,…,N are positive scalars as weights of the data fitting term. For normalization, we can always assume that ∑i=1Nwi=N. Moreover, λ is the regularized parameter to balance the data fitting term and the regularized term and ℙL denotes the polynomial space:(2)ℙL=spherical polynomials with degree≤L=spanYℓ,k,ℓ=0,…,L,k=1,…,2ℓ+1,where Yℓ,k is a fixed L2-orthonormal real spherical harmonic of degree ℓ and order k, which means(3)∫S2Yℓ,kxYℓ′,k′xdωx=δℓ,ℓ′δk,k′,ℓ,ℓ′=0,…,t;k,k′=1,…,2ℓ+1,where dωx denotes the surface measure on S2 and δℓ,ℓ′ is the Kronecker delta. It is well known that the dimension of ℙL is(4)dL:=dimℙL=∑ℓ=0L2ℓ+1=L+12.The Sobolev spaceℍs:=ℍsS2 is defined for s≥0 as the set of all functions f∈L2S2 whose Laplace–Fourier coefficients(5)f^ℓ,k=f,Yℓ,kL2S2=∫S2fxYℓ,kxdωxsatisfy(6)∑ℓ=0∞∑k=12ℓ+11+λℓsf^ℓ,k2<∞,where λℓ=ℓℓ+1, and the norm of ℍsS2 can be defined as(7)fℍs=∑ℓ=0∞∑k=12ℓ+11ϕℓsf^ℓ,k21/2,where the positive parameters ϕℓs are decreasing in ℓ and satisfy(8)ϕℓs≍1+λℓ−s≍ℓ+1−2s.Here, we say thatan≍bn if there exist positive constants c1 and c2 independent of n such that c1an≤bn≤c2an for all n. Obviously, by letting s=0, we can obtain ℍ0S2=L2S2. CRLS is to obtain the approximation of f with fitting the data values fδxi, wi=1…,N, and at the same time, we keep the approximation smoother. The motivation is often from the geophysical quantity problems such as satellite gravity gradiometry problem ([7], pp. 120, 262), in which a smoother polynomial representation of the gravitational potential than the Earth’s surface is to be found. On the other hand, CRLS can be seen as a generalization or approximation of some smooth-oriented regularized least squares models (see the appendix for discussions on some polynomial approximation models relative to this paper).Polynomial approximation problems arise in a variety of engineering areas, for which different models are used. To recover smooth functions onS2, a discrete regularized least squares method (DRLS) is proposed in [1, 5] as(9)ℒf=argming∈ℙL∑j=1Nwjgxj−fδxj2+λ∑j=1NwjRLgxj2,where RL:ℙL⟶ℙL is a linear “penalization” operator, which can be chosen in different ways. The points used are positive-weight cubature formula that is exact for all spherical polynomials of degree up to 2L, which will guarantee the uniqueness of the solution and furthermore derive its explicit form. A reconstruction error bound is obtained for the case that the nodes and weights establish a cubature formula that is exact for all spherical polynomials of degree up to 2L. To our best knowledge, the number of points needed to construct a 2L-degree polynomial-exact cubature on S2 is around 2L+12/3 [8], which is bigger than the dimension of the polynomial space ℙL in most practical cases. However, in real application, it is usually costly to acquire as much data points as expected, which motivates us to study the case that the number of data points is not large. In this paper, we do not assume that N≥L+12, that is, the number of data values is not necessary to reach or exceed the dimension of the polynomial space ℙL. The regularization term could still guarantee the uniqueness of the minimizer of CRLS when N<L+12 and avoid the overfitting.The rest of this paper is organized as follows. In Section2, we discuss the solution and approximation error of CRLS and the condition number of the linear system it derives. Based on the results, we further give a discussion on spherical tϵ-designs applied to this model. Numerical experiments are proposed in Section 3 in which the model is tested for its efficiency and accuracy.


## 2. Solution, Approximation Error, and Condition Number

### 2.1. Solution of CRLS

For simplicity, denoteβℓ:=1/ϕℓs for ℓ=0,…,L. A matrix form of problem (1) could be written as(10)minc∈ℝL+12W1/2Yc−W1/2fδ22+λB1/2c22,where Y∈ℝN×L+12 is defined as Yi,ℓ2+k=Yℓ,kxi, W=diagw1,…,wN∈ℝN×N, fδ=fδx1,…,fδxNT∈ℝN, and(11)B=diagβ0,β1,β1,β1,…,βℓ,…,βℓ︸2ℓ+1,…,βL,…,βL︸2L+1∈ℝL+12×L+12.Problem (1) is convex, and thus by deriving its first-order optimality condition, we can obtain that solving problem (1) leads to find a solution of the linear system(12)Tc=YTWfδ,with(13)T=YTWY+λB.Remark 1.
Forλ>0, equation (12) will always have a unique solution by the fact that YTWY=W1/2YTW1/2Y is positive semidefinite and B is positive definite.
By the uniqueness of the solution of problem (1), we denote the minimizer of it by(14)fλδ=argming∈ℙL∑i=1Nwigxi−fδxi2+λgℍs2.
Furthermore, whenλ=0, problem (1) may have multiple solutions. In this case, fλδ represents the uniquely determined minimizer in ℙL of (1) with minimal norm.


### 2.2. Approximation Error

For the case thatN≥t+12 and X is a spherical 2L-design or a cubature rule with 2L-degree polynomial exactness, the error bounds are discussed in [1, 5]. Model (1) can be seen as a variant of models considered in [1, 5], in which the difference is that we choose different regularization terms. Furthermore, to propose the error bounds, we do not need the assumption about the polynomial exactness the data point needs to have. The only assumption we need is on the mesh norm of the data point set X⊂Sd, which can be defined as follows.Definition 1.
(mesh norm). The mesh normhX of a point set X⊂Sd is(15)hX:=maxy∈S2minxi∈Xdisty,xi.
The mesh norm could be regarded as the largest geodesic distance of all points onSd to their nearest xi∈X or as the radius of the largest spherical cap that does not contain any xi∈X. The results in this section are based on a slight generalization of Theorem 5.1 in [9] which presents an L2 sampling inequality on the sphere Sd.Theorem 1.
Letd≥2 and let s>d/2. There exist positive constants c^s and c˜s (dependent on d and s) such that for any finite point set X=x1,…,xN of N distinct points on Sd with hX≤c^s and for any function g∈ℍs,(16)gL2≤c˜shXsgℍs+hXd/2∑i=1Ngxi21/2.For the setting of weighted least squares, we slightly revise the above theorem as the following.Theorem 2.
Letd≥2, s>d/2, and w∈1−ϵ,1−ϵ−1N with 0≤ϵ<1. There exist positive constants c^s and c˜s (dependent on d and s) such that for any finite point set X=x1,⋯,xN of N distinct points on Sd with hX≤c^s and for any function g∈ℍs,(17)gL2≤c˜shXsgℍs+hXd/21−ϵ∑i=1Nwigxi21/2.Proof.
The gap between Theorem 5.1 in [9] and the above one is the involvement of the weight vector w=w1,…,wN. By the fact that wi≥1−ϵ for ∀i=1,…,N and inequality (16), we can obtain the conclusion.Theorem 3.
Lets>1. For a point set X=x1,…,xN⊂S2 of N distinct points, a weight vector w=w1,…,wNT∈1−ϵ,1−ϵ−1N with 0≤ϵ<1 and a function f∈ℍs with values f=fx1,…,fxNT on X, and given a vector of deterministic noise δ=δ1,…,δNT and given λ≥0, there exist positive constants c^s and cs′ (dependent on s) such that if the point set X satisfies hX≤c^s, the following holds:(18)fλδ−fL2≤cs′hXs12λ1−ϵδ2+fℍs+hXd/21−ϵ2δ2+λ2fℍs.Proof.
For simplicity, we define the norm⋅w as yTWy for a vector y∈ℝN. Note that the norm is well defined since W is positive definite. Define the quadratic functional μg as(19)μg:=μλ,sg=∑i=1Nwigxi−fδxi2+λgℍs2,g∈ℙL.
Letfλδ, f∈ℝN be the vector of function value of fλδ, f at the point set X, respectively. By sampling inequality (17) in Theorem 2, we can obtain that there exists a positive constant cs″ such that(20)fλδ−fL2≤cs″hXsfλδ−fℍs+hXd/21−ϵfλδ−fW.
Then, by the fact thatfλδ is the minimizer of μg, the following holds:(21)fλδ−fW2=fλδ−fδ+δW2≤2fλδ−fW2+2δW2=2μfλδ−2λfλδℍs2+2δW2≤2μf+2δW2=2λfℍs2+4δW2.
Moreover, we have(22)fλδℍs2=1λμfλδ−fλδ−fδW2≤1λμf−fλδ−fδW2=1λδW2+λfℍs2−fλδ−fδW2≤1λδW2+fℍs2.
Thus, the termfλδ−fℍs can be bounded as(23)fλδ−fℍs≤fλδℍs+fℍs≤1λδW2+fℍs21/2+fℍs≤1λδW+2fℍs.
Substitute inequalities (21) and (23) into (20) with cs′=2cs″, and we have(24)fλδ−fL2≤12cs′hXs1λδW+2fℍs+hXd/21−ϵ2λfℍs2+4δW21/2≤cs′hXs12λδW+fℍs+hXd/21−ϵλ2fℍs+δW≤cs′hXs12λ1−ϵδ2+fℍs+hXd/21−ϵ2δ2+λ2fℍs,where the last inequality is derived by the fact that δw≤1/1−ϵδ2.
From the theorem, it is easy to see that the data point setX=x1,x2,…,xNand the weight w=w1,w2,…,wNT∈ℝN need to be chosen properly to improve the approximation quality of the model. The data set X should have its mesh norm as small as possible and w should be near to its mean value. Furthermore, one can assume that X and w establish a cubature rule with polynomial exactness so that it can give a fast way to find the solution of the problem. By this view, spherical tϵ-designs which are proposed and studied in [6, 10–14] are a suitable choice in this paper.Definition 2.
(sphericaltϵ-design). A spherical tϵ-design with 0≤ϵ<1 on Sd is a set of points X:=x1,…,xN⊂Sd, such that the quadrature rule(25)1N∑i=1Nwipxi=1Sd∫Sdpxdωdxis exact for all spherical polynomials p of degree at most t, with the weight vector w=w1,…,wNT satisfying(26)1−ϵ≤wi≤1−ϵ−1,i=1,…,N.
Naturally, we can choose a sphericaltϵ-design as the data set and its corresponding weight w as the weight in problem (1). Moreover, solving problem (1) can become trivial when t is large enough.Theorem 4.
LetX be a spherical tϵ-design with t≥2L and w be the vector of weights satisfying (25) and (26) with respect to X. Then,(27)YTWY=N4πIL+12,and (1) has the unique solution g=∑ℓ=0L∑k=12ℓ+1αℓ,kYℓ,kx with(28)αℓ,k=4πxiN+4πλβℓ∑i=1NwiYℓ,kxifδxi.Proof.
Note that whenX is a spherical tϵ-design,(29)YTWYℓ2+k,ℓ′2+k′=∑i=1NwiYℓ,kxiYℓ′,k′xi=N4π∫S2Yℓ,kxYℓ′,k′xdωx=N4πδℓℓ′δkk′,where the third equality is established by the orthonormality of spherical harmonics and the fact S2=4π. By Remark 1, problem (1) has a unique optimal solution satisfying(30)N4πIL+12+λBc=YTWfδ.
Note thatN/4πIL+12+λB is diagonal, and thus we can directly obtain the solution of the above linear system as equation (28). The proof is completed.


### 2.3. Condition Number

For a symmetric positive definite matrixM∈ℝdL×dL, let σ1M and σLM denote the largest and smallest eigenvalues of M and let κM=σ1M/σLM denote the condition number of M in the ⋅2 norm. The following theorem discusses the condition number and perturbation proposition of problem (1).Theorem 5.
Letfδ denote a perturbation of f. Then, the following holds:(31)κT≤σ1YTWY+λβLσdLYTWY+λβ0,(32)α−αδ2α2≤σ1YTWY+λβLσdLYTWY+λβ0⋅YTWf−fδ2YTWf2.Proof.
The proof of the theorem is similar to the proof of Theorem 4.1 in [1], and we provide quick illustration to keep the paper self-contained. The condition number of T could be obtained by the eigenvalue Weyl’s inequalities ([15], p. 181), which leads to(33)σ1T≤σ1YTWY+λσ1B=σ1YTWY+λβL,σLT≥σLYTWY+λσLB=σLYTWY+λβ0.
Note that problem (1) can be written as a ℓ2 least square problem:(34)minα∈ℝLdAα−b22,with(35)A=W1/2YλB1/2∈ℝN+dL×dL,b=W1/2f0∈ℝN+dL.
The perturbation bound could then be achieved by applying the standard least squares perturbation bound ([16], Theorem 1.4.6) as(36)α−αδ2α2≤ATb−bδ2ATb2κATA=YTWf−fδ2YWf2κT.
WhenX is a spherical tϵ‐ design with corresponding weight w=w1,…,wNT with t≥2L, by (27), we have that YTWY is a scalar multiple of the identity matrix, and from (13), we find(37)κT=N+4πλβLN+4πλβ0.
Therefore, estimate (31) of the condition number κT is sharp. It can be directly obtained by the above theorem that κT≍ℓ2s, which induces that s should be chosen properly as a balance of smoothing and solvability of the problem. However, the choice of s, which may depend on the nature of target function f∈ℂS2 and the scale of noise, will be an interesting question to be solved in the future.


## 2.1. Solution of CRLS

For simplicity, denoteβℓ:=1/ϕℓs for ℓ=0,…,L. A matrix form of problem (1) could be written as(10)minc∈ℝL+12W1/2Yc−W1/2fδ22+λB1/2c22,where Y∈ℝN×L+12 is defined as Yi,ℓ2+k=Yℓ,kxi, W=diagw1,…,wN∈ℝN×N, fδ=fδx1,…,fδxNT∈ℝN, and(11)B=diagβ0,β1,β1,β1,…,βℓ,…,βℓ︸2ℓ+1,…,βL,…,βL︸2L+1∈ℝL+12×L+12.Problem (1) is convex, and thus by deriving its first-order optimality condition, we can obtain that solving problem (1) leads to find a solution of the linear system(12)Tc=YTWfδ,with(13)T=YTWY+λB.Remark 1.
Forλ>0, equation (12) will always have a unique solution by the fact that YTWY=W1/2YTW1/2Y is positive semidefinite and B is positive definite.
By the uniqueness of the solution of problem (1), we denote the minimizer of it by(14)fλδ=argming∈ℙL∑i=1Nwigxi−fδxi2+λgℍs2.
Furthermore, whenλ=0, problem (1) may have multiple solutions. In this case, fλδ represents the uniquely determined minimizer in ℙL of (1) with minimal norm.


## 2.2. Approximation Error

For the case thatN≥t+12 and X is a spherical 2L-design or a cubature rule with 2L-degree polynomial exactness, the error bounds are discussed in [1, 5]. Model (1) can be seen as a variant of models considered in [1, 5], in which the difference is that we choose different regularization terms. Furthermore, to propose the error bounds, we do not need the assumption about the polynomial exactness the data point needs to have. The only assumption we need is on the mesh norm of the data point set X⊂Sd, which can be defined as follows.Definition 1.
(mesh norm). The mesh normhX of a point set X⊂Sd is(15)hX:=maxy∈S2minxi∈Xdisty,xi.
The mesh norm could be regarded as the largest geodesic distance of all points onSd to their nearest xi∈X or as the radius of the largest spherical cap that does not contain any xi∈X. The results in this section are based on a slight generalization of Theorem 5.1 in [9] which presents an L2 sampling inequality on the sphere Sd.Theorem 1.
Letd≥2 and let s>d/2. There exist positive constants c^s and c˜s (dependent on d and s) such that for any finite point set X=x1,…,xN of N distinct points on Sd with hX≤c^s and for any function g∈ℍs,(16)gL2≤c˜shXsgℍs+hXd/2∑i=1Ngxi21/2.For the setting of weighted least squares, we slightly revise the above theorem as the following.Theorem 2.
Letd≥2, s>d/2, and w∈1−ϵ,1−ϵ−1N with 0≤ϵ<1. There exist positive constants c^s and c˜s (dependent on d and s) such that for any finite point set X=x1,⋯,xN of N distinct points on Sd with hX≤c^s and for any function g∈ℍs,(17)gL2≤c˜shXsgℍs+hXd/21−ϵ∑i=1Nwigxi21/2.Proof.
The gap between Theorem 5.1 in [9] and the above one is the involvement of the weight vector w=w1,…,wN. By the fact that wi≥1−ϵ for ∀i=1,…,N and inequality (16), we can obtain the conclusion.Theorem 3.
Lets>1. For a point set X=x1,…,xN⊂S2 of N distinct points, a weight vector w=w1,…,wNT∈1−ϵ,1−ϵ−1N with 0≤ϵ<1 and a function f∈ℍs with values f=fx1,…,fxNT on X, and given a vector of deterministic noise δ=δ1,…,δNT and given λ≥0, there exist positive constants c^s and cs′ (dependent on s) such that if the point set X satisfies hX≤c^s, the following holds:(18)fλδ−fL2≤cs′hXs12λ1−ϵδ2+fℍs+hXd/21−ϵ2δ2+λ2fℍs.Proof.
For simplicity, we define the norm⋅w as yTWy for a vector y∈ℝN. Note that the norm is well defined since W is positive definite. Define the quadratic functional μg as(19)μg:=μλ,sg=∑i=1Nwigxi−fδxi2+λgℍs2,g∈ℙL.
Letfλδ, f∈ℝN be the vector of function value of fλδ, f at the point set X, respectively. By sampling inequality (17) in Theorem 2, we can obtain that there exists a positive constant cs″ such that(20)fλδ−fL2≤cs″hXsfλδ−fℍs+hXd/21−ϵfλδ−fW.
Then, by the fact thatfλδ is the minimizer of μg, the following holds:(21)fλδ−fW2=fλδ−fδ+δW2≤2fλδ−fW2+2δW2=2μfλδ−2λfλδℍs2+2δW2≤2μf+2δW2=2λfℍs2+4δW2.
Moreover, we have(22)fλδℍs2=1λμfλδ−fλδ−fδW2≤1λμf−fλδ−fδW2=1λδW2+λfℍs2−fλδ−fδW2≤1λδW2+fℍs2.
Thus, the termfλδ−fℍs can be bounded as(23)fλδ−fℍs≤fλδℍs+fℍs≤1λδW2+fℍs21/2+fℍs≤1λδW+2fℍs.
Substitute inequalities (21) and (23) into (20) with cs′=2cs″, and we have(24)fλδ−fL2≤12cs′hXs1λδW+2fℍs+hXd/21−ϵ2λfℍs2+4δW21/2≤cs′hXs12λδW+fℍs+hXd/21−ϵλ2fℍs+δW≤cs′hXs12λ1−ϵδ2+fℍs+hXd/21−ϵ2δ2+λ2fℍs,where the last inequality is derived by the fact that δw≤1/1−ϵδ2.
From the theorem, it is easy to see that the data point setX=x1,x2,…,xNand the weight w=w1,w2,…,wNT∈ℝN need to be chosen properly to improve the approximation quality of the model. The data set X should have its mesh norm as small as possible and w should be near to its mean value. Furthermore, one can assume that X and w establish a cubature rule with polynomial exactness so that it can give a fast way to find the solution of the problem. By this view, spherical tϵ-designs which are proposed and studied in [6, 10–14] are a suitable choice in this paper.Definition 2.
(sphericaltϵ-design). A spherical tϵ-design with 0≤ϵ<1 on Sd is a set of points X:=x1,…,xN⊂Sd, such that the quadrature rule(25)1N∑i=1Nwipxi=1Sd∫Sdpxdωdxis exact for all spherical polynomials p of degree at most t, with the weight vector w=w1,…,wNT satisfying(26)1−ϵ≤wi≤1−ϵ−1,i=1,…,N.
Naturally, we can choose a sphericaltϵ-design as the data set and its corresponding weight w as the weight in problem (1). Moreover, solving problem (1) can become trivial when t is large enough.Theorem 4.
LetX be a spherical tϵ-design with t≥2L and w be the vector of weights satisfying (25) and (26) with respect to X. Then,(27)YTWY=N4πIL+12,and (1) has the unique solution g=∑ℓ=0L∑k=12ℓ+1αℓ,kYℓ,kx with(28)αℓ,k=4πxiN+4πλβℓ∑i=1NwiYℓ,kxifδxi.Proof.
Note that whenX is a spherical tϵ-design,(29)YTWYℓ2+k,ℓ′2+k′=∑i=1NwiYℓ,kxiYℓ′,k′xi=N4π∫S2Yℓ,kxYℓ′,k′xdωx=N4πδℓℓ′δkk′,where the third equality is established by the orthonormality of spherical harmonics and the fact S2=4π. By Remark 1, problem (1) has a unique optimal solution satisfying(30)N4πIL+12+λBc=YTWfδ.
Note thatN/4πIL+12+λB is diagonal, and thus we can directly obtain the solution of the above linear system as equation (28). The proof is completed.


## 2.3. Condition Number

For a symmetric positive definite matrixM∈ℝdL×dL, let σ1M and σLM denote the largest and smallest eigenvalues of M and let κM=σ1M/σLM denote the condition number of M in the ⋅2 norm. The following theorem discusses the condition number and perturbation proposition of problem (1).Theorem 5.
Letfδ denote a perturbation of f. Then, the following holds:(31)κT≤σ1YTWY+λβLσdLYTWY+λβ0,(32)α−αδ2α2≤σ1YTWY+λβLσdLYTWY+λβ0⋅YTWf−fδ2YTWf2.Proof.
The proof of the theorem is similar to the proof of Theorem 4.1 in [1], and we provide quick illustration to keep the paper self-contained. The condition number of T could be obtained by the eigenvalue Weyl’s inequalities ([15], p. 181), which leads to(33)σ1T≤σ1YTWY+λσ1B=σ1YTWY+λβL,σLT≥σLYTWY+λσLB=σLYTWY+λβ0.
Note that problem (1) can be written as a ℓ2 least square problem:(34)minα∈ℝLdAα−b22,with(35)A=W1/2YλB1/2∈ℝN+dL×dL,b=W1/2f0∈ℝN+dL.
The perturbation bound could then be achieved by applying the standard least squares perturbation bound ([16], Theorem 1.4.6) as(36)α−αδ2α2≤ATb−bδ2ATb2κATA=YTWf−fδ2YWf2κT.
WhenX is a spherical tϵ‐ design with corresponding weight w=w1,…,wNT with t≥2L, by (27), we have that YTWY is a scalar multiple of the identity matrix, and from (13), we find(37)κT=N+4πλβLN+4πλβ0.
Therefore, estimate (31) of the condition number κT is sharp. It can be directly obtained by the above theorem that κT≍ℓ2s, which induces that s should be chosen properly as a balance of smoothing and solvability of the problem. However, the choice of s, which may depend on the nature of target function f∈ℂS2 and the scale of noise, will be an interesting question to be solved in the future.


## 3. Numerical Experiments

In this section, we report the numerical results for the CRLS model proposed in this paper. To illustrate the effectiveness of the model, we first compare the approximation errors of our model and the DRLS model (9) which is proposed in [1, 5]. The computation is implemented in Matlab 2016a and performed on a HP notebook computer equipped with Intel Core i7-6500U 2.5G Hz CPU, 16 GB RAM running Windows 10. The solver for linear systems generated by the two models is selected as the left matrix divide command embedded in Matlab, which indeed performs fast in the experiment.For the choice of the regularized parameterλ and the penalized parameters βℓ,ℓ=0,…,L for the DRLS model, one can refer to [5], in which βℓ are chosen based on the theoretical error bounds and real applications and λ is chosen by the balancing principle. The parameter choice of CRLS is also of great importance for the approximation quality of our model, and this will be a crucial problem we will tackle in the future work. The task for the numerical experiment here is to compare the numerical performance between CRLS and DRLS. What should also be noticed in the numerical experiments is that for an approximation ℋf of f, both the uniform and L2 approximation errors as ℋf−fCS2 and ℋf−fL2 are difficult to obtain exactly. In such case, the errors will be approximately computed by(38)f−ℋfL2=∫S2f−ℋf2dωx1/2≈4πN∑i=1Nfzi−ℋfzi21/2,f−ℋfCS2=maxz∈S2f−ℋf≈maxzi∈Xtfzi−ℋfzi,where Xt=z1,…,zNt⊂S2 is a large-scaled and well-distributed point set and we use it to estimate the errors. In this paper, Xt is selected as a maximal determinant (MD) (see [6, 17]) point set with 104 points.Example 1.
In the first example, we will test the two models, CRLS and DRLS, and record their approximation errors for comparison. The data point setX is also selected as the MD node of size 900. The penalized parameters are selected as βℓ=1+ℓ+1s and L=29 so that both models will achieve unique minimizers. The target function supposed to approximate is the Franke function [18]:(39)f1x≔f1x,y,z=0.75exp−9x−224−9y−224−9z−224+0.75exp−9x+1249−9y+110−9z+110+0.5exp−9x−724−9y−342−9z−542−0.2exp−9x−42−9y−72−9z−52,x,y,zT∈S2,which is not a spherical polynomial but continuously differentiable on the whole sphere. The noises in the data fδ obey a uniform distribution in −0.5,0.5 i.i.d. Images of f1 and its noisy version are shown in Figures 1(a) and 1(b).
Figure2 shows uniform and L2 approximation errors of the both models for s=0.15j, j=1,…,20 and λ=10−3.5+0.15k, k=1,…,20. It should be noted that both models will only be effective to recover the target function by choosing the parameters s and λ properly. Therefore, thresholds as a success of the approximation are introduced into the figure, that is, L2 error ≤1.75 and uniform error ≤1. As can be seen from the figure, CRLS shows its effectiveness and robustness in this example, which seems to be achieving success with more choices of the parameters than DRLS. Also, the deep blue area of the first row is larger than that in the second row, which indicates that CRLS has more chance than DRLS to achieve small errors.Figure 1
Franke function (noised) defined in (39): (a) f1; (b)f1δ.
(a)(b)Figure 2
Approximation errors onf1 with MD nodes: N=29 and L=29 (X-axis: logλ; Y-axis: s). (a) CRLS L∞ error. (b) CRLS L2 error. (c) DRLS L∞ error. (d) DRLS L2 error.
(a)(b)(c)(d)Example 2.
In this example, we will focus on the solving speed of the linear systems induced by the models. The condition number and CPU time about the linear systems derived by CRLS and DRLS will be reported. Two types of data point set with different scales will be chosen in this example. One is a well-distributed type which is used in Example1, and the other type is ill-distributed, called Gauss–Legendre point set (GL) [19] with Gauss points in cos(polar angle) and equally spaced points in azimuthal angle.
In Example1, the two models cannot achieve each one’s smallest errors in the same setting λ and s, even when other settings are the same. To deal with this point, for a fixed data set, we search for the optimal parameter setting of λ and s for each model to achieve least uniform errors separately with the scheme shown in Example 1 and report their errors, condition numbers, and CPU time in Table 1. From Table 1, it is obvious to see that for both MD and GL cases, CRLS and DRLS can achieve similar uniform errors at their best situations. However, in most cases, CRLS is more well conditioned than DRLS and costs less CPU time.
It is an interesting question to explain that CRLS performs better than DRLS in the simulation. First, we can compare the theoretical results on the condition numbers of the linear systems induced by CRLS and DRLS and then get the reason that CRLS has smaller condition numbers and solving times. However, the theoretical underlying reason that CRLS achieves smaller errors is unknown to us until now. We can try to explain it in a heuristic way: we know that the difference between CRLS and DRLS is the choice of the regularization term. In CRLS, we use the⋅ℍsS2 norm of the polynomial approximation g while in DRLS, we use the ℓ2 norm of the value of g with some operator transform. Thus, it can be easy to obtain that the regularization term in DRLS is determined on the choice of the data point set XN but in CRLS, it is not. This may possibly be the reason for CRLS obtaining smaller errors.Table 1
Numerical results of approximatef1 with different data point sets and models.
N, LECEDκCκDTCTDMD900,290.21120.21261.84×1045.53×1040.0180.0371600, 390.16940.17045.96×1043.16×1050.0380.0982500, 490.15280.15559.32×1041.22×1060.110.293600, 590.12830.12851.06×1053.22×1050.310.674900, 690.11990.11989.14×1051.18×1060.761.52GL496, 150.2719 (0.2700)8.77×1022.83×1030.0008 (0.0013)1326, 250.1811 (0.1893)4.30×1031.34×1040.0059 (0.010)1891, 300.1577 (0.1661)1.57×1043.01×1050.0095 (0.035)2556, 350.1578 (0.1574)9.28×1056.27×1060.024 (0.048)3321, 400.1648 (0.1656)1.72×1052.97×1040.039 (0.097)5151, 500.1192 (0.1209)1.91×1051.24×1050.13 (0.28)E: uniform errors of the model; κ: condition numbers of the model; T: CPU time. ⋅C refers to results from CRLS and ⋅D refers to the ones from DRLS.


## 4. Conclusion and Future Work

In this paper, a continuous regularized least squares polynomial approximation model on the two-sphere is discussed. Condition number and error bound of the model are proposed in the situation that the data are noisy and their cardinality is not necessarily larger than the dimension of the polynomial space. To guarantee the error bound, we need two assumptions: the mesh norm of the data point set should be small enough and the weights should be close to their mean. Numerical results show that the novel model is effective to approximate smooth functions on the sphere.We conclude this paper by pointing out a few problems for future directions. The first is that it is of great importance to study the strategy of choosing the parameterss and λ, which is crucial to the approximation quality of the model. Secondly, it is an interesting topic to explain that CRLS could achieve smaller errors than DRLS in the numerical experiments, since the difference between them is only the choices of the regularization term.



---

*Source: 9172385-2020-08-20.xml*