{'fraction_non_alphanumeric': 0.03229881669263486, 'fraction_numerical': 0.01250684296963827, 'mean_word_length': 4.812041116005874, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wavelets and wavelet-like transforms on the sphere and their application to geophysical data inversion

Many flexible parameterizations exist to represent data on the sphere. In addition to the venerable spherical harmonics, we have the Slepian basis, harmonic splines, wavelets and wavelet-like Slepian frames. In this paper we focus on the latter two: spherical wavelets developed for geophysical applications on the cubed sphere, and the Slepian"tree", a new construction that combines a quadratic concentration measure with wavelet-like multiresolution. We discuss the basic features of these mathematical tools, and illustrate their applicability in parameterizing large-scale global geophysical (inverse) problems.


INTRODUCTION
This paper is about parameterization, and its role in geophysical inverse problems: the analysis and representation of volumetric properties, with a particular emphasis on the three-dimensional ball and its surface, the twodimensional sphere. This is not the sole purview of geophysics (e.g. geodesy, geodynamics, seismology): there is a large literature on the subject in virtually every area of scientific inquiry (e.g. medical imaging, astronomy, cosmology, computer graphics, image processing, ...). For this reason we limit ourselves here to a few technical aspects, some of which have been published before (refs 1, 2) but are illustrated with different examples.
We begin with a new class of spherical wavelet basis on the ball designed for the analysis of geophysical models and for the tomographic inversion of global seismic data. 1,3,4 Its multiresolution character allows for modeling with an effective spatial resolution that varies with position within the Earth. We introduce two types of discrete wavelet transforms in the angular dimension of the "cubed sphere", a well-known Cartesian-to-spherical mapping. 5,6 These are applied to analyze the information of terrestrial topography data and in seismic wavespeed models of the Earth's mantle across scale space. The localization and sparsity properties of the wavelet bases allow finding a sparse solution to inverse problems by iterative minimization of a combination of the 2 norm of the data residuals and the 1 norm of the model wavelet coefficients. These have now been validated in realistic synthetic experiments to likely yield important gains in the inversion of seismic data for global seismic tomography, the procedure by which seismic waveforms are being used to image the three-dimensional distribution of elastic (compressional or P and shear or S ) wave speeds (V P and V S ) inside the Earth.
A second class of transforms is inspired by the type known as Slepian functions. These we understand to be families of orthogonal functions that are all defined on a common, e.g. geographical, domain, where they are either optimally quadratically concentrated or within which they are exactly limited, and which at the same time are exactly confined within a certain bandwidth, or maximally concentrated therein. Originally developed for the study of time series by Slepian,Landau and Pollak 7,8 in the 1960s, they have now been fully extended to the multidimensional Cartesian plane 9, 10 and the surface of the sphere, [11][12][13] where most of the geophysical applications lie. We have reported on some aspects of spherical Slepian functions and their role in the analysis and representation of geophysical data in previous papers in this series. 14,15 In this contribution we announce a hybrid construction that blends desirable aspects and properties of Slepian functions with ideas from wavelet and frame theory. We summarize the main ideas contained in Chapter 7, Multiscale Dictionaries of Slepian Functions on the Sphere, of the Ph. D. thesis by Eugene Brevdo, 2 deferring a full publication to a later date. arXiv:1109.1718v1 [physics.data-an] 8 Sep 2011

WAVELETS ON THE CUBED SPHERE
The use of wavelets is still no matter of routine in global geophysics, beyond applications in one and two Cartesian dimensions. This despite there being a wealth of available constructions on the sphere. [16][17][18][19][20][21][22][23][24][25][26][27][28] However, the cited studies are mostly concerned with surfaces, not volumes, and many of them require custom design and special algorithms. For the application to global seismic tomography that we envisage, we take advantage of the ease and flexibility of existing Cartesian constructions by implementing standard algorithms on a spherical-to-Cartesian map proposed by Ronchi et al. (1996), a grid that has long since proven its utility in the geosciences and beyond.
As in ref. 5, we define a coordinate 4-tuple (ξ, η, r, κ) for each of the κ = 1 → 6 "chunks" comprising this "cubed sphere" at radius r. In ref. 1 we list explicit formulas for the forward and inverse mapping of cubed-sphere to Cartesian coordinates. Here we suffice to say that a "master" surface chunk defined by the 2 × 2 N linearly spaced coordinate pairs −π/4 ≤ (ξ, η) ≤ π/4 maps, in a general sense, to the Cartesian coordinate vector which is rotated to occupy a total of six non-overlapping patches tiling the sphere, after which the entire construction is rotated over a standard set of Euler angles α = 0.0339, β = 1.1705, and γ = 1.1909. This results in the configuration shown in Figure 1 (left).
Alternatively, for reasons that will be made clear, we start from the coordinates −3π/8 ≤ (ξ, η) ≤ 3π/8 and apply the same rotations, thereby producing six overlapping "superchunks", as in Figure 1 (right). Throughout this paper we will quote N as the angular resolution level of our cubed sphere, which implies that it has 6 × 2 2N such elements, with typical seismic tomography grids having N = 7. The Euler angles used in our construction were chosen for geographical convenience, as can be seen from the continental outlines in Figure 2.
On these two grids we apply standard Cartesian transforms, e.g. orthogonal 29 and bi-orthogonal 30 ones, whereby for the case of the non-overlapping "chunks" construction we accommodate the presence of the "seams' by switching to special boundary filters at each of the edges, and applying preconditioners to the data prior to transformation in order to guarantee the usual polynomial cancellation throughout the closed rectangular interval, as is well established. 31 With this choice of bases, sparsity in the representation of many data types is to be generally expected. 32 In Figure 3 we explore coefficient statistics and the effects of thresholding on the reconstruction errors for a model of terrestrial topography. We focus on the fifth, or "North-American" chunk of our cubed sphere, and use the orthogonal D2 (Haar), D4 and D6 wavelet bases 29 on the interval, with preconditioning. 31 The top row uses the common conventions in plotting the wavelet and scaling coefficients in each of the bases after (hard) thresholding 32 such that only the coefficients larger than their value at the 85 th percentile level survive. The coefficients that have now effectively been zeroed out are left white in these top three panels. The middle series of panels of Figure 3 plots the spatial reconstruction after thresholding at this level; the root mean squared (rms) errors of these reconstructions are quoted as a percentage of the original root mean squared signal strengths. The thresholded wavelet transforms allow us to discard, as in these examples, 85% of the numbers required to make a map of North American topography in the cubed-sphere pixel basis: the percentage error committed is only 6.3%, 5.2% and 9.2% according to this energy criterion in the D2, D4 and D6 bases, respectively. From the map views it is clear that despite the relatively small error, the D2 basis leads to unsightly block artifacts in the reconstruction, which are largely avoided in the smoother and more oscillatory D4 and D6 bases. A view of the coefficient statistics is presented in the lowermost three panels of Figure 3. The coefficients are roughly lognormally distributed, which helps explain the success of the thresholded reconstruction approach. We conclude that the D4 basis is a good candidate for geophysical data representation. show the wavespeed anomalies from the average at that depth. The wavelet transform in the D4 basis (with special boundary filters and after preconditioning, and up until scale J = 3) was thresholded and the results re-expanded to the spatial grid, identically as we did for the topography in Figure 3. The results for specific values of the thresholding (quoted as the percentile of the original wavelet coefficients) are shown in Figs 4b and 4f for the 85 th , and Figs 4c and 4g for the 95 th , respectively. At each level of thresholding the number of nonzero wavelet/scaling expansion coefficients is quoted: at 0% thresholding this number is identical to the number of pixels in the surficial cubed sphere being plotted.
The reconstruction error can be visually assessed from the pictures; it is also quoted next to each panel as the percentage of the root mean squared error between the original and the reconstruction, normalized by the root mean squared value of the original in the original pixel representation, in percent. Specifically, we calculate and quote the ratio of 2 norms in the pixel-basis model vector m, which, in the lower-right annotations is called the "% error norm". We write A for any of the wavelet (analysis) transforms that can be used and S (synthesis) for their inverses, and T for the hard thresholding of the wavelet and scaling coefficients. In Figs 4e and 4j, the misfit quantity (2) is represented as a black line relevant to the left ordinate labeled " 2 error norm", which shows its behavior at 1% intervals of thresholding; the filled black circles correspond to the special cases shown in the map view. Only after about 80% of the coefficients have been thresholded does the error rise above single-digit percentage levels, but after that, the degradation is swift and inexorable. The blue curves in Figs 4e and 4j show another measure relevant in this context, namely the ratio of the 1 norms of the thresholded wavelet coefficients compared to the original ones, in percent, or The 2 ratios (2) in the black curves (and the left ordinate) evolve roughly symmetrically to the 1 ratios (3) in the blue curves (and the right ordinate), though evidently their range is different. The third measure that is being plotted as the red curve is the "total variation" norm ratio, in percent, namely

THE SPHERICAL SLEPIAN TREE TRANSFORM
Geophysical and cosmological signals are constrained by the physical processes that generate them and conditioned by the sensors that observe them. On the real line and in the plane, both physical and sampling constraints lead to assumptions of a bandlimit: that a signal contains no energy outside some supporting region in the frequency domain. Bandlimited signals on the sphere are zero save for the low-frequency sphericalharmonic components. Spherical harmonics are not orthogonal on arbitrary portions of the sphere, yet for the study of geophysical processes we may not have access to nor interest in signal originating from outside a specific geographic region of interest. By construction, Slepian functions satisfy the dual constraints of bandlimitation and spatial concentration by optimization of an energy criterion that concentrates as much energy as possible into the region of interest for a given bandlimit. 13 With Slepian function bases, however, it is "all or nothing": given a certain bandwidth and a spatial region of interest, the functions (variably) fill the entire bandwidth range and ultimately the entire spatial target. The wavelet constructions that we introduced in Section 2, on the other hand, were compactly supported and had multiresolution properties which the Slepian basis does not possess.
Here we develop an algorithm for the construction of Slepianesque dictionary elements that are bandlimited, localized, and multiscale. It is based on a subdivision scheme that constructs a binary tree from subdivisions of the region of interest. Such a dictionary has many nice properties: it closely overlaps with the most concentrated Slepian functions on the region of interest, and most element pairs have low coherence. Therefore, they, too, should be eminently suitable to solve ill-posed inverse problems in geophysics and cosmology. Though the new dictionary is no longer composed of purely orthogonal elements like the Slepian basis, it can also be combined with modern inversion techniques that promote sparsity in the solution. Thereby they provide significantly lower residual error after reconstruction when compared to "classically optimal" Slepian inversion techniques. 41

A Multiscale Dictionary of Slepian Functions
We now turn our focus to numerically constructing a dictionary D of functions that can be used to approximate bandlimited signals on the sphere and allows for the reconstruction of signals from their point samples. Let R ⊂ S 2 be a simply connected subset of the sphere. With L the bandwidth, the dictionary D will be composed of functions bandlimited to spherical harmonic degrees 0 ≤ l ≤ L. The construction is based on a binary tree. The node capacity, n, is a positive integer. Each node of the tree corresponds to the first n Slepian functions with bandlimit L and concentrated on a subset R ⊂ R. The top tree node is for the entire region R, and each node's children correspond to a division of R into two roughly equally sized subregions. As the child nodes will be concentrated in disjoint subsets of R , all of their corresponding functions and children are effectively incoherent. We now fix a height H of the tree: the number of times to subdivide R. The height is determined as the maximum number of binary subdivisions of R that can have n well concentrated functions. That is, we find the minimum integer H such that n ≥ N 2 −H |R|,L , with N the Shannon number, 13 which has the solution H = log 2 |R| 4π (L + 1) 2 n .
To complete the top-down construction, it remains to decide how to subdivide a region R into equally sized subregions. For roughly circular connected domains, the first Slepian function has no sign changes, and the second Slepian function has a single zero-level curve that subdivides the region into approximately equal areas; when R is a spherical cap, the subdivision is exact. 13 We thus subdivide a region R into the two nodal domains associated with the second Slepian function on that domain; see Figure 8 for a visualization of the subdivision scheme as applied to the African continent.

Concentration, Range, and Incoherence
The utility of the tree construction presented above depends on its ability to represent bandlimited functions in a region R, and its efficacy at reconstructing functions from point samples in R. These properties, in turn, reduce to questions of concentration, range, and incoherence. First, dictionary D is concentrated in R if its functions are concentrated in R. Second, the range of dictionary D is the subspace spanned by its elements. Ideally, the basis formed by the first few Slepian functions on R is a subspace of the range of D. Third, when D is incoherent, pairwise inner products of its elements have low amplitude: pairs of functions are approximately orthogonal. This, in turn, is a useful property when using D to estimate signals from point samples.
Unlike the concentration eigenvalues corresponding to the Slepian functions on R, not all of the eigenvalues of the elements of D R reflect their concentration within this top-level (parent) region. We thus define the modified concentration value Recalling that d (j,α) 2 = 1, the value ν is simply the percentage of energy of the (j, α) th element that is concentrated in R. This value is always larger than the element's eigenvalue, which relates its fractional energy within the smaller subset R (j) .
The size of dictionary D R,L,n is generally larger than the Shannon number N |R|,L = (L + 1) 2 / |R| /4/π for any node capacity n, and as a result it cannot form a proper basis: it has too many functions. Ideally, then, we require that elements of the dictionary span the space of the first N |R|,L Slepian functions. One possible answer to the question of the range of D is given by studying the angle between the subspaces 42 Here, P B is the orthogonal projection onto space B and all of the norms are with respect to the given subspace. The angle (A, B) is symmetric, non-negative, and zero iff A ⊂ B or B ⊂ A; furthermore it is invariant under unitary transforms applied to both on A and B, and admits a triangle inequality. It is thus a good indicator of distance between two subspaces; furthermore, it can be calculated accurately given two matrices whose columns span A and B. We can therefore identify the matrices A and B with the subspaces spanned by their columns.
Let G R,L 1:α denote the matrix containing the first α column vectors which are the spherical harmonic expansion coefficients of the traditional Slepian functions for a given region R and a bandwidth L. Let D denote the (L + 1) 2 × |D R,L,n | matrix containing the spherical harmonic representations of the elements of the dictionary D R,L,n . Figure 9 shows ( G 1:α , D) for R = Africa with L = 36. The Shannon number is N Africa,36 ≈ 79. It is clear that while the dictionaries D Africa,36,1 and D Africa,36,2 do not strictly span the space of functions bandlimited to L = 36 and optimally concentrated in Africa, they are a close approximation: the column span of G Africa,36 1:α is nearly linearly dependent with the spans of D Africa,36,1 and D Africa,36,2 , for α significantly larger than N Africa, 36 .
The requirement that the dictionary elements form an orthogonal basis is less important than the property of mutual incoherence. This we identify when the inner product between pairs of elements is almost always very low. Figure 10 shows that the two tree constructions on continental Africa have good incoherency properties: most dictionary element pairs are nearly orthogonal. As can be seen from Figure 10a, most pairwise inner products are nearly zero. More specifically, as expected, dictionary elements (j, 1), (2j, 1), (2j + 1, 1), (2(2j), 1), (2(2j) + 1, 1), (2(2j + 1), 1), (2(2j + 1) + 1, 1), . . ., tend to have large inner products, while those with non-overlapping borders do not. This exact property is also visible in the two diagonal submatrices of Figure 10b. In the off-diagonals, due to the orthonormality of the construction, elements of the form (j, 1) and (j, 2) are orthogonal. In contrast, due to the nature of the tree subdivision scheme, elements of the form (2j, 1) or (2j + 1, 1) and (j, 2) have a large inner product. However, the number of connections between nodes and their ancestors is O(n [2 H H]), while the total number of pairwise inner products is O([n2 H ] 2 ); and for reasonably sized values of L the ratio of ancestral connections to pairwise inner products gets to be small, see Figure 11.

CONCLUSIONS
We have introduced two classes of spherical parameterizations and discussed their properties in the context of geophysical model analysis and representation. The first involved the porting of traditional Cartesian wavelet transforms to the sphere via a cubed-sphere mapping, with and without special consideration for the boundaries between the six chunks constituting the entire spherical surface. The second was a novel elaboration of the classical ideas of signal concentration on the sphere. Starting from optimally spatially concentrated bandlimited spherical "Slepian" functions, we construct a dictionary of functions occupying a binary tree, where the elements of the tree are successive levels of Slepian functions calculated for an increasingly subdivided spatial domain.
Both the wavelet bases and the Slepian tree frames that we discuss in this paper are flexible and efficient ways of studying the information content and spatiospectral structure of geophysical and cosmological data. But in both cases their utility will also be derived from using them in the parameterization of geophysical inverse problems. In ref. 1 we discuss an iterative algorithm that solves an inverse problem in seismic tomography while promoting sparsity in the cubed-sphere wavelet basis in which the unknown model is expressed. Ref. 2 explores the inverse problem of approximating bandlimited, heterogeneously concentrated, essentially multiscale functions on the sphere from their samples.
In this paper and elsewhere, we have shown by example that many geophysical signals are sparse in the wavelet and Slepian (tree) "bases". To tackle large-scale and potentially ill-conditioned inverse problems, the sparsity of the solution should actively be encouraged. In either case the inverse problem comes down to minimizing a mixed 2 -1 functional of the form whereby d is a data set, K a linear operator, and S a certain synthesis map that transforms the set of unknown coefficients w into the domain of K. Thus, S could relate to the wavelet bases introduced in Section 2 or to the Slepian-tree dictionary introduced in Section 3. Numerical experiments conducted for a variety of experimental setups under geophysically realistic conditions have already indicated that both types of constructions will be amenable to solving inverse problems in geophysics and beyond. Further research will be reported elsewhere.

ACKNOWLEDGMENTS
We    . Wavelet and scaling coefficients (top), space-domain reconstructions after thresholding (middle), and "signed" histograms (bottom) of the wavelet and scaling coefficients of the "North American" face of the cubed-sphere version of the Earth's topography, to a 2 J dyadic subdivision with J = 3. We used preconditioned interval wavelet transforms. All coefficients were hard-thresholded, retaining only the 15% largest coefficients by absolute value. In the top row, the locations of zeroed coefficients are rendered white; those are also captured by the white bars in the histograms. The root mean squared (rms) error of the reconstruction after thresholding is indicated as a percentage of the signal rms. Tick marks on the color bars identify the 5 th , 50 th and 95 th percentile of the coefficients or the spatial reconstructions after thresholding, respectively. Interior ticks on the histograms roughly coincide with these same percentiles as applied to either the positive and negative coefficients when expressed on a logarithmic scale. Histograms for the positive coefficients point up and have ordinates in positive percentages, histograms for the negative coefficients point down and have ordinates in negative percentages; these percentages are with respect to the total number of positive and negative coefficients. The blue and red shaded areas of the histograms reflect the coefficients retained at the global 85 th thresholding level. , and to a 2 J dyadic subdivision (J = 3). As a function of the percentage of the coefficients that are being thresholded, and relatively to the original unthresholded values, the bottom panels quote the spatial 2 norms of the reconstruction error (in black), the total variation norms of the reconstructed images in the space domain (in red), and the 1 norms of the coefficients that remain (in blue). The reconstructions remain faithful to the originals even at elevated levels of thresholding.  d (2,·) d (3,·) d (7,·) d (6,·) d (5,·) d (4,·) Figure 6. The binary tree subdivision scheme and associated dictionary DR,L for the multiscale Slepian tree transform. We define the top-level region R as R (1) and the generic subsets R as R (j) .