The paper analyses a spectral approach to reconstructing a scalar field on the sphere given only information about a masked version of the field together with precise information about the smooth mask. The theory is developed for a general mask and later specialised to the case of an axially symmetric mask. Numerical experiments are given for the case of an axial mask motivated by the cosmic microwave background assuming that the underlying field is a realisation of a Gaussian random field with an artificial angular power spectrum of moderate degree ell le 100. The recovery is highly satisfactory in the absence of noise and even in the presence of moderate noise.
['Mathematics - Numerical Analysis', 'General Relativity and Quantum Cosmology', '15A23', '15A29', '60G60']
Removing the mask  reconstructing a scalar field on the sphere from a masked field
{'REMOVING THE MASK - RECONSTRUCTING A SCALAR FIELD ON THE SPHERE FROM A MASKED FIELD': 'JAN HAMANN, QUOC T. LE GIA, IAN H. SLOAN, AND ROBERT S. WOMERSLEY  \nABSTRACT. The paper analyses a spectral approach to reconstructing a scalar field on the sphere, given only information about a masked version of the field together with precise information about the (smooth) mask. The theory is developed for a general mask, and later specialised to the case of an axially symmetric mask. Numerical experiments are given for the case of an axial mask motivated by the cosmic microwave background, assuming that the underlying field is a realisation of a Gaussian random field with an artificial angular power spectrum of moderate degree ( ℓ ≤ 100 ). The recovery is highly satisfactory in the absence of noise and even in the presence of moderate noise.', '1. INTRODUCTION': "In this paper, we study the reconstruction of a scalar field (for example temperature or pressure) on the unit sphere, given (possibly noisy) data on a masked version of the field, together with precise knowledge of the mask. The underlying motivation is the cosmic microwave background (CMB) for which the temperature observations, so important for the modern understanding of the early universe, are obscured over substantial portions of the sky by our own Milky Way, creating the need for masking some portions before attempting reconstruction.  \nThe paper aims at a proof-of-concept for a new spectral approach to such problems. While the theory is general, the numerical experiments are restricted to the case of an axially symmetric mask, and limit the field's angular power spectrum to polynomial degree ℓ ≤ 100 , corresponding to an angular resolution of approximately 2 · . Within these limitations, the recovery is shown to be highly satisfactory in the no-noise case, and also in the case of moderate noise.  \nThere is a rich literature on the inpainting problem for the particular case of the CMB [1, 20], with techniques based on harmonic methods [4, 9], iterative methods [15, 8], constrained Gaussian realisations [10, 5], group sparse optimisation methods [11] or neural networks [25, 17, 18, 13, 24].  \nClosest to the present approach is the work of Alonso et al. [2], which however differs in aim (which was to recover the angular power spectrum from a knowledge of the 'pseudo C ℓ ', through pooling together Fourier coefficients of different degrees).  \nThe problem is formulated in the next section, by reducing the problem to that of solving a large ill-posed linear system. Section 3 establishes properties of the matrix in that linear system. Section 4 outlines our stochastic approach to the solution of the linear system. Because the full linear system is currently beyond our resources, in Section 5 we obtain a large reduction in  \ndifficulty by specialising to the case of an axially symmetric mask. The final section is devoted to numerical experiments in both the no-noise and noisy cases.", '2. THE PROBLEM SETTING': "Taking r to be any point in the unit sphere S 2 in R 3 (i.e. r is a unit vector in R 3 ), it can be represented in spherical coordinates as  \nr = (sin θ cos ϕ, sin θ sin ϕ, cos θ ) , θ ∈ [0 , π ] , ϕ ∈ [0 , 2 π ) .  \nThe underlying real scalar field a ( r ) is assumed to be partially obscured by a known mask v = v ( r ) to give a masked field a v = a v ( r ) := a ( r ) v ( r ) .  \nIt is well-known [12] that the space of square-integrable functions L 2 ( S 2 ) admits an orthonormal basis formed by spherical harmonics  \n{ Y ℓ,m : ℓ = 0 , . . . ; m = -ℓ, . . . , ℓ } ,  \nwhere the Y ℓ,m is given explicitly by  \n(2.1) Y ℓ,m ( θ, ϕ ) = √ 2 ℓ +1 4 π ( ℓ -m )! ( ℓ + m )! P m ℓ (cos θ ) e imϕ , m ≥ 0 , Y ℓ,m ( θ, ϕ ) = ( -1) m Y ℓ, -m ( θ, ϕ ) , m< 0 ,  \nwhere { P m ℓ } denote the associated Legendre function, which is defined in terms of Legendre polynomials { P ℓ : ℓ -0 , 1 , . . . } by  \nP m ℓ ( t ) = ( -1) m (1 -t 2 ) m/ 2 d m d t m P ℓ ( t ) , m = 0 , 1 , . . . , ℓ ; ℓ = 0 , 1 , 2 , . . .  \nWe assume that the incompletely known field a is a spherical polynomial of degree at most L , and hence expressible as an expansion in terms of orthonormal (complex) spherical harmonics Y ℓ,m ( r ) of degree ℓ ≤ L ,  \n(2.2) a ( r ) = L ∑ ℓ =0 ℓ ∑ m = -ℓ a ℓ,m Y ℓ,m ( r ) , r ∈ S 2 ,  \nwhere  \na ℓ,m := ∫ S 2 a ( r ) Y ℓ,m ( r )d S ( r ) ,  \nwhere d S ( r ) = sin θ d θ d ϕ , following from the orthonormality of the spherical harmonics,  \n∫ S 2 Y ℓ,m ( r ) Y ℓ ' ,m ' ( r )d S ( r ) = δ ℓ,ℓ ' δ m,m ' , ℓ, ℓ ' ≥ 0 , m = -ℓ, . . . , ℓ, m ' = -ℓ ' , . . . , ℓ ' .  \nWe assume that the mask can be adequately approximated by a partial sum of spherical harmonics series as 1  \n(2.3) v ( r ) = K ∑ k =0 k ∑ ν = -k v k,ν Y k,ν ( r ) ,  \nwhere  \nv k,ν := ∫ S 2 v ( r ) Y k,ν ( r )d S ( r ) .  \nThus the masked signal a v ( r ) = a ( r ) v ( r ) is expressible as a spherical harmonic expansion of degree L + K ,  \n(2.4) a v ( r ) = a ( r ) v ( r ) = L + K ∑ j =0 j ∑ µ = -j a v j,µ Y j,µ ( r ) , r ∈ S 2 ,  \nwhere  \na v j,µ := ∫ S 2 a v J ( r ) Y j,µ ( r )d S ( r ) = ∫ S 2 a ( r ) v ( r ) Y j,µ ( r )d S ( r ) (2.5) = ∫ S 2 ( L ∑ ℓ =0 ℓ ∑ m = -ℓ a ℓ,m Y ℓ,m ( r ) )( K ∑ k =0 k ∑ ν = -k v k,ν Y k,ν ( r ) ) Y j,µ ( r )d S ( r ) = L ∑ ℓ =0 ℓ ∑ m = -ℓ K ∑ k =0 k ∑ ν = -k a ℓ,m v k,ν ∫ S 2 Y ℓ,m ( r ) Y k,ν ( r ) Y j,µ ( r )d S ( r ) = L ∑ ℓ =0 ℓ ∑ m = -ℓ E j,µ ; ℓ,m a ℓ,m ,  \nwhere, for j = 0 , . . . , L + K , µ = -j, . . . , j and ℓ = 0 , . . . , L , m = -ℓ, . . . , ℓ ,  \n(2.6) E j,µ ; ℓ,m = K ∑ k =0 k ∑ ν = -k ∫ S 2 Y ℓ,m ( r ) Y k,ν ( r ) Y j,µ ( r )d S ( r ) v k,ν .  \nThus the essential task in the reconstruction is to solve as accurately as possible the large linear system  \n(2.7) L ∑ ℓ =0 ℓ ∑ m = -ℓ E j,µ ; ℓ,m a ℓ,m = a v j,µ ,  \nwhere E is defined by (2.6).  \nEquation (2.7) can be considered as an overdetermined (but possibly not full-rank) set of linear equations for the a ℓ,m . We will find it convenient to replace the upper limit L + K in (2.7) by a more flexible upper limit J with L ≤ J ≤ L + K . Then the equation can be written as  \n(2.8) E a = a v ,  \nwhere E is a ( J +1) 2 × ( L +1) 2 matrix,  \n(2.9) a = ( a ℓ,m , ℓ = 0 , . . . , L, m = -ℓ, . . . , ℓ ) ∈ C ( L +1) 2 .  \nand  \n(2.10) a v = ( a v j,µ , ℓ = 0 , . . . , J, m = -ℓ, . . . , ℓ ) ∈ C ( J +1) 2 .  \nIn Section 4 we come to the most challenging part of the paper, which is the approximate solution of the ill-posed linear system, and before it the computation of the matrix E . Before then, however, it is useful to establish properties of the matrix E .", '3. PROPERTIES OF E': 'This section summarizes useful properties of the matrix E defined in (2.6). We first note that the product of three spherical harmonics (known as a Gaunt coefficient) can be evaluated in terms of Wigner 3j symbols, see eg. [6, Eq. 34.4.22] or [2]. Explicitly,  \nD ℓ,m ; k,ν ; j,µ := ∫ S 2 Y ℓ,m ( r ) Y k,ν ( r ) Y j,µ ( r )d S ( r ) (3.1) = ( -1) µ √ (2 ℓ +1)(2 k +1)(2 j +1) 4 π × ( ℓ k j 0 0 0 )( ℓ k j m ν -µ ) .  \nAs important special cases,  \n(3.2) D ℓ,m ; k,ν ; j,µ = 0 if          j + ℓ + k is odd, or k < | j -ℓ | , or k > j + ℓ, or m + ν = µ.  \n̸  \nThe following lemma gives several elementary properties of the matrix E , beginning with an explicit integral expression in terms of the mask function v .', 'Lemma 3.1. The elements of the matrix E satisfy': "(3.3) E j,µ ; ℓ,m = ∫ S 2 Y j,µ ( r ) Y ℓ,m ( r ) v ( r ) d S ( r ) ,  \nand, for a real mask v ,  \nE ℓ,m ; j,µ = E j,µ ; ℓ,m , (3.4)  \nE j,µ ; ℓ,m = K ∑ k =0 [ D ℓ,m ; k, 0; j,µ v k, 0 +2 k ∑ ν =1 ℜ ( D ℓ,m ; k,ν ; j,µ v k,ν ) ] , (3.5)  \nE j,µ ; ℓ, -m = ( -1) m -µ E j, -µ ; ℓ,m . (3.6)  \nProof. Firstly, from the definition of E in (2.6) we have  \nE j,µ ; ℓ,m = K ∑ k =0 k ∑ ν = -k ∫ S 2 Y ℓ,m ( r ) Y k,ν ( r ) Y j,µ ( r ) d S ( r ) v k,ν = ∫ S 2 Y j,µ ( r ) Y ℓ,m ( r ) K ∑ k =0 k ∑ ν = -k v k,ν Y k,ν ( r ) d S ( r ) = ∫ S 2 Y j,µ ( r ) Y ℓ,m ( r ) v ( r ) d S ( r ) ,  \nestablishing (3.3). From the definition of D ℓ,m ; k,ν ; j,µ in (3.1) as an integral, together with the spherical harmonic property  \n(3.7)  \nit follows that  \nD j,µ ; k,ν ; ℓ,m ; = ( -1) ν D ℓ,m ; k, -ν ; j,µ .  \nBecause both the mask v and the field a are real we have  \n(3.8) a ℓ,m = ( -1) m a ℓ, -m , v k,µ = ( -1) µ v k, -µ ,  \nfor all relevant values of ℓ, m, k and µ , and (3.4) then follows from (2.6). Also, from (3.1) and (3.7),  \nD ℓ, -m ; k,ν ; j,µ = ( -1) m D ℓ,m ; k,ν ; j,µ , D ℓ,m ; k, -ν ; j,µ = ( -1) ν D ℓ,m ; k,ν ; j,µ ,  \nso (3.8) gives  \nD ℓ,m ; k, -ν ; j,µ v k, -ν = D ℓ,m ; k,ν ; j,µ v k,ν ,  \nand (2.6) then yields (3.5). Finally, for a real mask (3.6) follows easily from (3.3).  \n□  \n- 3.1. Singular values of E . This subsection gives upper bounds on the singular values of the rectangular matrix E in terms of the real mask v . We use E ∗ to denote the complex conjugate transpose of the matrix E .  \nTheorem 3.2. Assume that the mask v ( r ) is approximated by the partial sum of its spherical Fourier series of degree K ≥ 1 and that L ≤ J ≤ L + K is the degree of the approximation a v J ( r ) to the masked field a v ( r ) , so that E is a ( J +1) 2 by ( L +1) 2 matrix. The singular values σ of E satisfy  \n(3.9) 0 ≤ σ ≤ v max  \n,  \nY ℓ, -m ( r ) = ( -1) m Y ℓ,m ( r ) ,  \nwhere  \nv max := max r ∈ S 2 | v ( r ) | .  \n̸  \nProof. Let u = 0 be an eigenvector of the positive semi-definite Hermitian matrix E ∗ E corresponding to the non-negative real eigenvalue σ 2 , so E ∗ E u = σ 2 u . Then, using (3.3), and writing the elements of u as u ℓ,m , ℓ = 0 , . . . , L, m = -ℓ, . . . , ℓ , we have  \n( E u ) j,µ = L ∑ ℓ =0 ℓ ∑ m = -ℓ u ℓ,m ∫ S 2 Y j,µ ( r ) Y ℓ,m ( r ) v ( r )d S ( r ) = ∫ S 2 Y j,µ ( r ) u ( r ) v ( r )d S ( r ) ,  \nfor j = 0 , . . . , J , µ = -j, . . . , j , where  \nu ( r ) := L ∑ ℓ =0 ℓ ∑ m = -ℓ u ℓ,m Y ℓ,m ( r ) ,  \ngiving  \nu ∗ E ∗ E u = ∥ E u ∥ 2 ℓ 2 = J ∑ j =0 j ∑ µ = -j ∣ ∣ ∣ ∣ ∫ S 2 u ( r ) v ( r ) Y j,µ ( r )d S ( r ) ∣ ∣ ∣ ∣ 2 .  \nThus, using Parseval's identity for uv and then u ,  \nσ 2 ∥ u ∥ 2 ℓ 2 = u ∗ ( σ 2 u ) = u ∗ ( E ∗ E u ) = J ∑ j =0 j ∑ µ = -j ∣ ∣ ∣ ∣ ∫ S 2 u ( r ) v ( r ) Y j,µ ( r )d S ( r ) ∣ ∣ ∣ ∣ 2 (3.10) = ∥ u v ∥ 2 L 2 ≤ ( v max ) 2 ∥ u ∥ 2 L 2 = ( v max ) 2 ∥ u ∥ 2 ℓ 2 .  \nThis gives the upper bound (3.9) on the singular values σ of E .  \n□  \n- 3.2. Eigenvalues of E . In this subsection we take J = L , making the matrix E square. From the second statement in Lemma 3.1, E is Hermitian thus its eigenvalues are real, and eigenvectors belonging to distinct eigenvalues are orthogonal.  \n̸  \nLet λ ∈ R be an eigenvalue of E and let q = 0 be a corresponding eigenvector, thus  \nE q = λ q .  \nThe following result provides both lower and upper bounds on λ , in terms of the minimum and maximum values of the mask v .  \nTheorem 3.3. Assume that the mask v is approximated by the partial sum of its spherical Fourier series of degree K ≥ 1 . Assume also that J = L , so that the matrix E is square. Then the eigenvalues of E lie in the interval ( v min , v max ) , where  \nv min := min r ∈ S 2 v ( r ) , v max := max r ∈ S 2 v ( r ) .  \nProof. Let λ ∈ R be an eigenvalue of E , with corresponding eigenvector q ∈ C ( L +1) 2 . Then from (3.3),  \nq ∗ E q = L ∑ j =0 j ∑ µ = -j L ∑ ℓ =0 ℓ ∑ m = -ℓ q j,µ E j,µ ; ℓ,m q ℓ,m (3.11) = L ∑ j =0 j ∑ µ = -j L ∑ ℓ =0 ℓ ∑ m = -ℓ q j,µ q ℓ,m ∫ S 2 Y j,µ ( r ) Y ℓ,m ( r ) v ( r ) d S ( r ) = ∫ S 2 q ( r ) q ( r ) v ( r )d S ( r ) ,  \nwhere  \nq ( r ) := L ∑ ℓ ∑ q ℓ,m Y ℓ,m ( r )  \n(3.12) ℓ =0 m = -ℓ .  \nIt follows that  \n∫ (3.13)  \n∥ q ∥ 2 L 2 v max = J ∑ j =0 j ∑ µ = -j | q j,µ | 2 v max = ∥ q ∥ 2 ℓ 2 v max ,  \nλ ∥ q ∥ 2 ℓ 2 = q ∗ ( λ q ) = q ∗ E q = S 2 | q ( r ) | 2 v ( r ) d S ( r ) <  \nwhere the inequality is strict because v , being a spherical polynomial of non-zero degree, cannot be identically equal to either its maximum or minimum value. Similarly, we have a lower bound  \nλ ∥ q ∥ 2 ℓ 2 > ∥ q ∥ 2 ℓ 2 v min ,  \ntogether proving λ ∈ ( v min , v max ) .  \n□  \nNote that even if the true (non-polynomial) mask lies in [0 , 1] for all r ∈ S 2 , the Gibbs phenomenon will typically produce oscillations in v , making v min < 0 and v max > 1 .  \nWe also note in passing that q is an eigenvector belonging to λ for the integral equation  \n∫ S 2 K L ( r , r ' ) q ( r ' ) d r ' = λq ( r ) ,  \nwhere K L ( r , r ' ) is the integral kernel given by  \nK L ( r , r ' ) = L ∑ ℓ =0 ℓ ∑ m = -ℓ Y ℓ,m ( r ) Y ℓ,m ( r ' ) v ( r ' ) (3.14) = L ∑ ℓ =0 2 ℓ +1 4 π P ℓ ( r · r ' ) v ( r ' ) ,  \nand in the last step we used the addition theorem for spherical harmonics. Here P ℓ is the Legendre polynomial of degree ℓ , normalised so that P ℓ (1) = 1 .", '4. SOLVING E a = a v': "As with any ill-posed system, it is essential to build in a priori knowledge of the solution. Neumaier [14, Section 8], knowing that the true solution of an ill-posed problem is generally smooth, controls the smoothness through a smoothing operator S . However, in this problem a smoothing operator would not be appropriate because the solution is the opposite of smooth, since for each ℓ, m the unknown quantity a ℓ,m is a realisation of an independent random variable. That is a property we must build into the solution. Accordingly, we assume, in accordance with the usual assumptions for the CMB, that the a ℓ,m are mean-zero uncorrelated random variables with covariance ( C ℓ ) L ℓ =0 , where C ℓ is real. Details on using Gaussian random fields to model the CMB can be found in the book by Marinucci and Peccati [12].  \nWe allow general J in the range L ≤ J ≤ L + K , giving a linear system with ( J + 1) 2 equations, so typically an over-determined linear system with more equations than unknowns, implying that an exact solution does not in general exist.  \nMoreover, we assume that the original field coefficients a ℓ,m are corrupted by noise, so the actual model is  \n(4.1) E a ε = a v + ε v = ( a + ε ) v ,  \nwhere ε is a vector of independent mean-zero random variables ε ℓ,m with a diagonal covariance matrix Υ , and a ε is an approximation to a . We also assume that the ε ℓ,m and the a ℓ,m are all statistically independent, so that in terms of expected values we have  \n⟨ ε ℓ,m ⟩ = 0 , ⟨ a ℓ,m ⟩ = 0 , ⟨ ε ℓ,m a ℓ ' ,m ' ⟩ = 0 , ⟨ a ℓ,m a ℓ ' ,m ' ⟩ = C ℓ δ ℓ,ℓ ' δ m,m ' , (4.2) ⟨ ε ℓ,m ε ℓ ' ,m ' ⟩ = Υ ℓ,m δ ℓ,ℓ ' δ m,m ' .  \nWe deduce the following expectations of quadratic forms:  \n⟨ aa ∗ ⟩ = Ω ⟨ a v a ∗ ⟩ = ⟨ ( E a ) a ∗ ⟩ = E Ω ⟨ a v ( a v ) ∗ ⟩ = ⟨ ( E a )( a ∗ E ∗ ) ⟩ = E Ω E ∗ (4.3) ⟨ εa ∗ ⟩ = 0 ⟨ εε ∗ ⟩ = Υ , ⟨ ε v ( ε v ) ∗ ⟩ = E Υ E ∗ ,  \nwhere  \n(4.4) Ω ℓ,m ; ℓ ' ,m ' = C ℓ δ ℓ,ℓ ' δ m,m ' .  \nLet Λ ∈ C ( L +1) 2 × ( L +1) 2 be a real symmetric-positive definite matrix, with associated norm ∥ a ∥ Λ = ( a ∗ Λ a ) 1 2 defined by  \n(4.5) ∥ a ∥ 2 Λ = a ∗ Λ a = tr [ aa ∗ Λ] = tr [Λ aa ∗ ] ,  \nwhere we used the matrix property tr( AB ) = tr( BA ) . The following theorem gives a condition for minimising the expected squared Λ -norm error of an approximate solution of (4.1). It is an extension/specialisation of [14, Theorem 8], which that author attributes to [3].  \nTheorem 4.1. Consider the over-determined linear system E a ε = a v + ε v , where a v = E a , ε v = E ε and a and ε have the stochastic properties in (4.3) . Assume that the ( J +1) 2 × ( L +1) 2 matrix E has full rank ( L +1) 2 . Among all approximations of the form a ε ≈ Q ( a v + ε v ) , where Q is a non-random ( L +1) 2 × ( J +1) 2 matrix, the expected squared error ⟨∥ a -Q ( a v + ε v ) ∥ 2 Λ ⟩ is minimized by any solution ̂ Q of the equation  \n(4.6) ̂ QE (Ω + Υ) = Ω .  \nThe resulting minimum expected squared error is  \n(4.7) 〈 ∥ a -̂ a ∥ 2 Λ 〉 = tr [ Λ ( Ω -Ω(Ω + Υ) -1 Ω )] ,  \nwhere ̂ a := ̂ Q ( a v + ε v ) .  \nRemark 1 . The minimizer ̂ a is in general not unique.  \nProof. Writing y ε := a v + ε v , a general linear approximation can be written as Q y ε = ̂ Q y ε + R y ε , where ̂ Q is an as yet unknown minimizer, and R = Q -̂ Q is a matrix in C ( L +1) 2 × ( J +1) 2 . The mean square error can now be written as  \n∥ a -( ̂ Q + R ) y ε ∥ 2 Λ = ∥ a -̂ Q y ε ∥ 2 Λ + ∥ R y ε ∥ 2 Λ -2 ℜ [( a -̂ Q y ε ) ∗ ) Λ R y ε ] = ∥ a -̂ Q y ε ∥ 2 Λ + ∥ R y ε ∥ 2 Λ -2 ℜ [ tr ( Λ R y ε ( a -̂ Q y ε ) ∗ )] .  \nOn taking expected values and using (4.3) we have  \n〈 tr(Λ R y ε ( a -̂ Q y ε ) ∗ 〉 = tr ( Λ R 〈 y ε ( a -̂ Q y ε ) ∗ 〉) = tr ( Λ R 〈 ( a v + ε v )( a -̂ Q ( a v + ε v )) ∗ 〉) = tr ( Λ R ( E ⟨ aa ∗ ⟩ -E ⟨ aa ∗ + εε ∗ ⟩ E ∗ ̂ Q ∗ )) = tr ( Λ R ( E Ω -E (Ω + Υ) E ∗ ̂ Q ∗ )) .  \nSo  \n〈 ∥ a -( ̂ Q + R ) y ε ∥ 2 Λ 〉 = 〈 ∥ a -̂ Q y ε ∥ 2 Λ 〉 + 〈 ∥ R y ε ∥ 2 Λ 〉 -2 ℜ tr ( Λ R ( E Ω -E (Ω + Υ) E ∗ ̂ Q ∗ )) .  \nBy definition, ̂ Q is a minimizer of ⟨∥ a -Q y ε ∥ 2 Λ ⟩ , so the linear term must vanish for all R . More precisely, we must have  \n(4.8) E Ω = E (Ω + Υ) E ∗ ̂ Q ∗ , or equivalently ̂ QE (Ω + Υ) E ∗ = Ω E ∗ ,  \nsince otherwise by taking R to be ( E Ω -E (Ω + Υ) E ∗ ̂ Q ∗ ) ∗ we obtain a contradiction. It is easily seen that the second equality in (4.8) is equivalent to (4.6): starting with (4.6) by right  \nmultiplying by E ∗ , starting with (4.8) by right multiplying by E and using the invertibility of E ∗ E , noting that E ∗ E is a square matrix of full rank ( L +1) 2 . If (4.6) holds then we have  \n〈 ∥ a -( ̂ Q + R ) y ε ∥ 2 Λ 〉 = 〈 ∥ a -̂ Q y ε ∥ 2 Λ 〉 + 〈 ∥ R y ε ∥ 2 Λ 〉 ≥ 〈 ∥ a -̂ Q y ε ∥ 2 Λ 〉  \nwith equality for R = 0 , corresponding to a = Q y .  \nThe expected squared error is  \n̂ ̂ ε  \n〈 ∥ a -̂ a ∥ 2 Λ 〉 = 〈 ∥ a -̂ Q y ε ∥ 2 Λ 〉 = 〈 ( a -̂ Q y ε ) ∗ Λ ( a -̂ Q y ε ) 〉 = tr [ Λ 〈 ( a -̂ Q y ε )( a -̂ Q y ε ) ∗ 〉] = tr [ Λ ( ⟨ aa ∗ ⟩ + 〈 ̂ Q y ε y ∗ ε ̂ Q ∗ 〉 -2 ℜ ( ̂ Q ⟨ y ε a ∗ ⟩ ) )] = tr [ Λ ( Ω+ ̂ Q ( E Ω E ∗ + E Υ E ∗ ) ̂ Q ∗ -2 ℜ ( ̂ QE Ω ) )] ,  \nNow by (4.8)  \n̂ QE Ω = ̂ Q ( E Ω E ∗ + E Υ E ∗ ) ̂ Q ∗ ,  \nand hence  \ntr[Λ ̂ QE Ω] = tr[Λ ̂ QE (Ω + Υ) E ∗ ̂ Q ∗ ] = tr[(Ω + Υ) 1 / 2 E ∗ ̂ Q ∗ Λ ̂ QE (Ω + Υ) 1 / 2 ] = ∥ ̂ QE (Ω + Υ) 1 / 2 ] ∥ 2 Λ ,  \nwhich is real, implying  \n〈 ∥ a -̂ a ∥ 2 Λ 〉 = tr [ Λ ( Ω -̂ QE Ω )] = tr [ Λ ( Ω -Ω(Ω + Υ) -1 Ω )] .  \n□  \nRemark 2 . Note that the equation (4.8) determining the minimizer ̂ Q y ε of the expected meansquare error does not depend on the matrix Λ , i.e. on the choice of quadratic norm. For example, using Λ = I or Λ = Ω does not change ̂ Q .  \nCorollary 4.2. Under the conditions of Theorem 4.1, let Γ be an arbitrary positive definite matrix of size ( J +1) 2 × ( J +1) 2 . Then a vector ̂ a ∈ R ( L +1) 2 that achieves the minimal error given in (4.7) is  \n(4.9) ̂ a := Ω(Ω + Υ) -1 α ,  \nwhere α ∈ R ( L +1) 2 is the unique solution of  \n(4.10) E ∗ Γ E α = E ∗ Γ( a v + ε v ) .  \nProof. The matrix ̂ Q defined by  \n(4.11) ̂ Q := Ω(Ω + Υ) -1 ( E ∗ Γ E ) -1 E ∗ Γ  \nis easily seen to satisfy the condition (4.6) in Theorem 4.1. Equally, it is easily seen that the corresponding minimizer  \n̂ a := ̂ Q ( a v + ε v )  \ncan be written exactly as stated in the corollary.  \n□  \nRemark 3 . The corollary gives our prescription for computing the coefficient vector ̂ a . Note that the postprocessing step in (4.9) is easily carried out given that the matrices Ω and Υ are diagonal, since each element of α is by this step merely reduced by a known factor. Note also that equation (4.10) is just the normal equation for the linear system if, as we shall assume in practice, Γ is the identity matrix. Formation of the normal equations can greatly increase the condition number of an already ill-conditioned system. In practice we shall address the ill-conditioning either by QR factorisation of the matrix E , or (less desirably) by adding a regularising term to the right-hand side, to obtain  \n(4.12) ( E ∗ Γ E +Σ) α = E ∗ Γ( a v + ε v ) ,  \nwhere Σ is an empirically chosen positive definite ( L +1) 2 × ( L +1) 2 matrix.  \nThe following proposition shows that if the elements of a v and ε v have the correct symmetry for real-valued fields a v = av and ϵ v = ϵv , then the computed values of ̂ a also follow the same symmetry. The practical importance of this result is that the symmetry property, since it occurs naturally, does not need to be enforced.  \nProposition 4.3. Assume that the components of a v and ε v satisfy  \na v j,µ = ( -1) µ a v j, -µ and ε v j,µ = ( -1) µ ε v j, -µ , µ = -j, . . . j, j ≥ 0 .  \nAssume also that the positive definite matrices Ω , Υ and Γ , and also Σ if present, are all diagonal, and that their diagonal elements are positive numbers independent of the second label µ or m . Then ̂ a given by (4.10) and (4.9) satisfies  \n̂ a ℓ,m = ( -1) m ̂ a ℓ, -m , m = -ℓ, . . . , ℓ, ℓ ≥ 0 .  \nProof. We first show that the components of b := E ∗ Γ( a v + ε v ) satisfy  \nb ℓ,m = ( -1) m b ℓ, -m , m = -ℓ, . . . , ℓ, ℓ ≥ 0 .  \nWe have, using (3.6),  \nb ℓ, -m = ∑ j ∑ µ ( E ∗ ) ℓ, -m ; j,µ Γ j ( a v j,µ + ε j,µ ) = ∑ j ∑ µ ( -1) m -µ ( E ∗ ) ℓ,m ; j, -µ Γ j ( a v j,µ + ε j,µ ) = ( -1) m ∑ j ∑ µ ( E ∗ ) ℓ,m ; j, -µ Γ j ( a v j, -µ + ε j, -µ ) = ( -1) m b ℓ,m ,  \nas required. A similar argument shows that  \n( E ∗ Γ E +Σ) ℓ, -m ; ℓ ' ,m ' = ( -1) m -m ' ( E ∗ Γ E +Σ) ℓ,m ; ℓ ' , -m ' .  \nSince α is the unique solution of  \n( E ∗ Γ E +Σ) α = b ,  \nby taking the ( ℓ, -m ) component of this equation we obtain  \n∑ ℓ ' ∑ m ' ( E ∗ Γ E +Σ) ℓ, -m ; ℓ ' m ' α ℓ ' ,m ' = b ℓ, -m ,  \nwhich with the above symmetry properties leads to  \n∑ ℓ ' ∑ m ' ( -1) m -m ' ( E ∗ Γ E +Σ) ℓ,m ; ℓ ' -m ' α ℓ ' ,m ' = ( -1) m b ℓ,m ;  \nOn taking the complex conjugate and dividing by ( -1) m this gives us  \n∗ ,  \n(4.13) ( E Γ E +Σ) c = b  \nwhere  \n(4.14)  \nc ℓ ' , -m ' := ( -1) m ' α ℓ ' ,m ' , m ' = -ℓ ' , . . . , ℓ ' , ℓ ' ≥ 0 .  \nWe see by uniqueness of the solution of (4.13) that c = α , thus by (4.14) the vector α has the desired symmetry. Multiplication by Ω(Ω+Υ) -1 clearly preserves the symmetry, thus the proof is complete. □", '5. AXIALLY SYMMETRIC MASKS': 'Ageneral mask v ( r ) leads to a large dense matrix E , of size ( J +1) 2 ( L +1) 2 × ( J +1) 2 ( L +1) 2 , see (2.6), a size beyond present resources if J and L are in the hundreds. In this section we consider the more tractable special case in which the mask v is axially symmetric, i.e.  \nv ( r ) = v ( θ, ϕ ) = v ( θ ) ,  \nwith v being a function of the polar angle θ and independent of the azimuthal angle ϕ . For this case we have  \nv k,ν = w k δ ν, 0 ,  \nwhere  \n(5.1) w k := v k, 0 = ∫ S 2 v ( r ) Y k, 0 ( r )d S ( r ) = √ π (2 k +1) ∫ π 0 v ( θ ) P k (cos( θ )) sin θ d θ.  \nNote that w k is real, and that  \nw k = 0 if k is odd and also v ( -r ) = v ( r ) .  \nIn this case it follows from (3.2) that E j,µ ; ℓ,m = 0 unless µ = m . Thus it is convenient to introduce a new notation,  \n(5.2) .  \nE ( m ) j,ℓ := E j,m ; ℓ,m  \nEquation (2.7) now becomes  \n(5.3) L ∑ ℓ =0 E ( m ) j,ℓ a ℓ,m = a v j,m  \n,  \nin which the coefficients belonging to different values of m are completely decoupled. This can be seen as just a special case of (2.8), albeit with uncoupled values of m , thus all of the analysis in Sections 3 and 4 remains applicable.  \nNote that from (3.5) of Lemma 3.1, E ( m ) j,ℓ is real and symmetric, E ( m ) j,ℓ = E ( m ) ℓ,j . Moreover  \nE ( m ) j,ℓ = 0 if j + ℓ is odd and v ( -r ) = v ( r ) , or if ℓ < | m | or j < | m | .  \nWe can treat E ( m ) j,ℓ as an ( J -| m | +1) × ( L -| m | +1) matrix, but how should we choose J ? The choice J = L inevitably leads to a poorly conditioned linear system. There would seem to be considerable benefit, at least in theory, in taking the largest value J = L + K , to ensure that the resulting overdetermined linear system makes use of all available information.  \nThe equation to be solved in practice is, instead of (4.10), now  \n(5.4) ( ( E ( m ) ) ∗ Γ a E ( m ) ) α = ( E ( m ) ) ∗ Γ a ( a v + ε v );  \nOr if regularisation is desired, then instead of (4.12) the equation to be solved becomes  \n(5.5) ( ( E ( m ) ) ∗ Γ a E ( m ) +Σ a ) α = ( E ( m ) ) ∗ Γ a ( a v + ε v ) .  \nHere Γ a and Σ a have the same diagonal values as Γ and Σ , but the second label on rows and columns has now disappeared, and the new matrices are of size J × J and L × L respectively.', '6. NUMERICAL EXPERIMENTS': "Recall that our goal is to reconstruct a scalar random field on the sphere given only a masked and noisy version of the field. We have seen in previous sections that the problem can be reduced to the solution of the overdetermined linear system  \n(6.1) E a = b v ,  \nwhere the given data b v = ( a + ε ) v = a v + ε v are the spherical harmonic coefficients of the masked noisy map with the coefficients a corrupted by independent Gaussian noise ε with mean ⟨ ε ⟩ = 0 and variance ⟨ εε ∗ ⟩ = Υ .  \nTo illustrate the potential of the method we consider numerical experiments where we know the 'true' solution a , so we can calculate errors to test performance and the effect of model parameters, including taking ε = 0 .  \nFIGURE 1. Original Gaussian random field  \n<!-- image -->  \nIn the experiments a specified angular power spectrum C ℓ , ℓ = 2 , . . . , L is used to generate an instance of a Gaussian random field with known spherical harmonic coefficients a at the N pix = 50 , 331 , 648 HEALPix 2 [7] points ( N side = 2048 ), using the HealPy 3 package [26]. Noise is then added as described in Subsection 6.3. The mask is then applied pointwise to the noisy map, and the masked noisy map used to calculate the spherical harmonic coefficients b v , again using the HealPy package. The next step is to estimate the original Fourier coefficients a ℓ,m using (4.9) together with (4.10), or alternatively using (4.9) with the regularised equation (4.12). The final step is to reconstruct the target field from its Fourier coefficients.  \nWe consider a Gaussian random field with the artificial angular power spectrum  \nC ℓ = g ( ℓ L +1 ) ,  \nwhere  \ng ( x ) = { 1 for 0 ≤ x ≤ 1 / 2 -2 x +2 for 1 / 2 ≤ x ≤ 1 .  \nIn the experiments, we assume that L = 100 and K = 900 . A realisation of the random field is shown in Figure 1. We shall use this realisation as the target field a ( r ) in all the following experiments.  \n- 6.1. An axially symmetric mask. When the mask applied to the noisy data is axially symmetric, the problem decomposes into independent problems for each value of m , as in Section 5. To construct the mask we first define the following non-decreasing function p ∈ C 3 ( R ) :  \n(6.2) p ( x ) =      0 for x ≤ 0 , x 4 (35 -84 x +70 x 2 -20 x 3 ) for 0 < x < 1 , 1 for x ≥ 1 .  \nThen as a function of the Cartesian coordinate z ∈ [ -1 , 1] of a point on the sphere, our mask is, for 0 < a z < b z < 1  \n(6.3) v ( z ) = p ( | z | -a z b z -a z ) .  \nIn Figure 2, an axially symmetric mask with a z = π 2 -10 π 180 and b z = π 2 -20 π 180 is plotted on the unit sphere. This mask has the value 1 (and hence has no masking effect) for points on the sphere  \nFIGURE 2. An axially symmetric C 3 mask with a z = π 2 -10 π 180 and b z = π 2 -20 π 180 .  \n<!-- image -->  \nmore than 20 · from the equator, the value 0 (complete masking) within 10 · from the equator, and smooth variation in between, through the function p , see (6.2), with an argument expressed as a function of the z coordinate of a point on the sphere. (We do not use a discontinuous mask because a discontinuous function has slow convergence of its Fourier series, leading to Gibbs' phenomenon for the truncated Fourier series.)  \nFIGURE 3. The C 3 mask v ( θ, ϕ ) with a = 10 π 180 , b = 20 π 180 .  \n<!-- image -->  \nThe transformation z := cos( θ ) = cos( π 2 -φ ) , with latitude φ = π 2 -θ ∈ [ -π 2 , π 2 ] , gives  \nv ( φ, ϕ ) = p ( | cos( π 2 -φ ) | -a z b z -a z ) .  \nIn terms of latitude the transition region is | φ | ∈ [ a, b ] where a = 10 π 180 and b = 20 π 180 , as illustrated in Figure 3.  \nThe masked field, including the addition of noise as described in Section 6.3, is illustrated in Figure 4.  \nUsing this axially symmetric mask all the rectangular matrices E = E ( m ) described in Section 6.1 for m = 0 , . . . , L of sizes ( J +1 -m ) × ( L +1 -m ) with J = L + K are pre-computed in parallel. Here K is the maximum multipole in the spherical harmonics approximation of v as in (2.3). We used the sympy package [21] to compute the entries for the matrix E . Fast quadrature methods on the unit sphere [16] could be used in a future implementation.  \n6.2. Numerical condition of the problem. Figures 5a and 5b illustrate the condition of the matrices E ( m ) for L = 100 , K = 900 and the mask in Figure 3. The largest singular values in Figure 5a are consistent with the upper bound on the singular values in Theorem 3.2. Even though the values of the mask lie in [0 , 1] , the polynomial approximation of the mask may have values slightly outside this interval due to the Gibbs phenomenon. The ill-conditioning of the matrices E ( m ) , as illustrated in Figure 5b, which can be severe especially for small m , arises from the smallest singular value.", 'REMOVING THE MASK': 'Reconstructed field with Gaussian noise with power spectrum 10-2C1  \nFIGURE 8. QR method for Gaussian noise with angular power spectrum Υ ℓ = 10 -2 C ℓ  \n<!-- image -->', '7. CONCLUSION': 'In this paper we have analysed a spectral method for recovering a scalar field from a masked and possibly noisy version of that field. The quality of the recovery might be considered acceptable even in the presence of noise. However, it is acknowledged that the quality will deteriorate as the noise level is increased and as the cutoff polynomial degree is increased from 100 .  \nFIGURE 7. QR method for Gaussian noise with angular power spectrum Υ ℓ = 10 -4 C ℓ  \n<!-- image -->', 'ACKNOWLEDGEMENTS': 'The assistance of Yu Guang Wang in the early stages of the project and constructive comments from anonymous referees are gratefully acknowledged. This research includes computations using the computational cluster Katana supported by Research Technology Services at', 'HAMANN AND LE GIA, SLOAN, WOMERSLEY': 'UNSW Sydney [19]. Some of the results in this paper have been derived using the healpy and HEALPix packages.', 'REFERENCES': "- [1] P. Abrial, Y. Moudden, J. L. Starck, J. Fadili, J. Delabrouille, and M. K. Nguyen. CMB data analysis and sparsity. Stat. Meth. , 5:289, 2008.\n- [2] David Alonso, Javier Sanchez, and Anˇze Slosar. A unified pseudoC ℓ framework. Mon. Not. Roy. Astron. Soc. , 484(3):4127-4151, 2019.\n- [3] M. Bertero, C. De Mol, and G. A. Viano. The stability of inverse problems. In Inverse scattering problems in optics , volume 20 of Topics Current Phys. , pages 161-214. Springer, Berlin-New York, 1980.\n- [4] Pawel Bielewicz, K. M. Gorski, and A. J. Banday. Low order multipole maps of CMB anisotropy derived from WMAP. Mon. Not. Roy. Astron. Soc. , 355:1283, 2004.\n- [5] Martin Bucher and Thibaut Louis. Filling in CMB map missing data using constrained Gaussian realizations. Mon. Not. Roy. Astron. Soc. , 424:1694, 2012.\n- [6] NIST Digital Library of Mathematical Functions . http://dlmf.nist.gov/, Release 1.1.4 of 2022-01-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.\n- [7] K. M. G'orski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. Astrophysical Journal , 622:759-771, April 2005.\n- [8] H. F. Gruetjen, J. R. Fergusson, M. Liguori, and E. P. S. Shellard. Using inpainting to construct accurate cut-sky CMB estimators. Phys. Rev. D , 95(4):043532, 2017.\n- [9] Kaiki Taro Inoue, Paolo Cabella, and Eiichiro Komatsu. Harmonic Inpainting of the Cosmic Microwave Background Sky I:Formulation and Error Estimate. Phys. Rev. D , 77:123539, 2008.\n- [10] Jaiseung Kim, Pavel Naselsky, and Nazzareno Mandolesi. Harmonic in-painting of CMB sky by constrained Gaussian realization. Astrophys. J. Lett. , 750:L9, 2012.\n- [11] Chao Li and Xiaojun Chen. Group sparse optimization for inpainting of random fields on the sphere. IMA Journal of Numerical Analysis , page drad071, 09 2023.\n- [12] Domenico Marinucci and Giovanni Peccati. Random fields on the sphere , volume 389 of London Mathematical Society Lecture Note Series . Cambridge University Press, Cambridge, 2011. Representation, limit theorems and cosmological applications.\n- [13] Gabriele Montefalcone, Maximilian H. Abitbol, Darsh Kodwani, and R. D. P. Grumitt. Inpainting CMB maps using Partial Convolutional Neural Networks. JCAP , 03:055, 2021.\n- [14] Arnold Neumaier. Solving ill-conditioned and singular linear systems: A tutorial on regularization. SIAM Review , 40(3):636-666, 1998.\n- [15] Atsushi J. Nishizawa and Kaiki Taro Inoue. Reconstruction of Missing Data using Iterative Harmonic Expansion. Mon. Not. Roy. Astron. Soc. , 462(1):588-600, 2016.\n- [16] Gerlind Plonka, Daniel Potts, Gabriele Steidl, and Manfred Tasche. Numerical Fourier analysis . Applied and Numerical Harmonic Analysis. Birkhauser/Springer, Cham, second edition, [2023] ©2023.\n- [17] Giuseppe Puglisi and Xiran Bai. Inpainting Galactic Foreground Intensity and Polarization Maps Using Convolutional Neural Networks. Astrophys. J. , 905(2):143, 2020.\n- [18] Alireza Vafaei Sadr and Farida Farsian. Filling in Cosmic Microwave Background map missing regions via Generative Adversarial Networks. JCAP , 03:012, 2021.\n- [19] D. Smith and L. Betbeder-Matibet. Katana, 2010.\n- [20] J. L. Starck, M. J. Fadili, and A. Rassat. Low-l CMB Analysis and Inpainting. Astron. Astrophys. , 550:A15, 2013.\n- [21] SymPy Development Team. Sympy.\n- [22] E. van den Berg and M. P. Friedlander. Probing the pareto frontier for basis pursuit solutions. SIAM Journal on Scientific Computing , 31(2):890-912, 2008.  \n- [23] E. van den Berg and M. P. Friedlander. SPGL1: A solver for large-scale sparse reconstruction, December 2019. https://friedlander.io/spgl1.\n- [24] Guo-Jian Wang, Hong-Liang Shi, Ye-Peng Yan, Jun-Qing Xia, Yan-Yun Zhao, Si-Yu Li, and Jun-Feng Li. Recovering the CMB Signal with Machine Learning. Astrophys. J. Supp. , 260(1):13, 2022.\n- [25] Kai Yi, Yi Guo, Yanan Fan Fan, Jan Hamann, and Yu Guang Wang. Cosmo vae: Variational autoencoder for cmb image inpainting. In 2020 International Joint Conference on Neural Networks (IJCNN) . IEEE, jul 2020.\n- [26] Andrea Zonca, Leo Singer, Daniel Lenz, Martin Reinecke, Cyrille Rosset, Eric Hivon, and Krzysztof Gorski. healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in python. Journal of Open Source Software , 4(35):1298, March 2019.  \nSCHOOL OF PHYSICS, UNSW SYDNEY, AUSTRALIA.  \nEmail address : jan.hamann@unsw.edu.au  \nSCHOOL OF MATHEMATICS AND STATISTICS, UNSW SYDNEY, AUSTRALIA.  \nEmail address : qlegia@unsw.edu.au, i.sloan@unsw.edu.au, r.womersley@unsw.edu.au"}