# Comparison of several approximation schemes on the Cubed Sphere

Jean-Pierre Croisille, Jean-Baptiste Bellet, Matthieu Brachet

► To cite this version:

Jean-Pierre Croisille, Jean-Baptiste Bellet, Matthieu Brachet. Comparison of several approximation schemes on the Cubed Sphere. EGU General Assembly 2024, Apr 2024, Vienna, Austria. 2024, 10.5194/egusphere-egu24-15777. hal-04620381

HAL Id: hal-04620381

https://hal.science/hal-04620381

Submitted on 21 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comparison of several approximation schemes on the Cubed Sphere

<figure><img src="image_1.png" /></figure>

## Goals

* Compare interpolation, weighted least squares and hyperinterpolation on the Cubed Sphere grid.

* Get an efficient scheme for spherical approximation using spherical harmonics.

## Cubed Sphere grid

Equiangular Cubed Sphere with step $\frac{\pi}{2N}$, $CS_N \subset S^2$

<figure><img src="image_2.png" /></figure>

Fig. CS<sub>N</sub> (white dots) is obtained by projection of a meshed circumscribed cube (N = 3).

## Data approximation

Find $f \in \mathcal{Y}_D$ such that $f(x) \approx y(x)$, $x \in CS_N$.

* Data: $y = [y(x)]_{x \in CS_N} \in \mathbb{R}^{6N^2+2}$

* Unknown: $f \in \mathcal{Y}_D$, with Fourier coefficients $\hat{f} \in \mathbb{R}^{(D+1)^2}$

* Matrix writing: find $\hat{f} \in \mathbb{R}^{(D+1)^2}$ such that $A_N^D \hat{f} \approx y$

* Application: $F: S^2 \to \mathbb{R}$ is approximated by $f$ from $y = F|_{CS_N}$

## Lexicographical interpolation

$$ \inf_{\substack{f \in \mathcal{Y}_{4N-1} \\ f|_{CS_N} = y}} (\|\hat{f}_{4N-1}\|, \|\hat{f}_{4N-2}\|, \dots, \|\hat{f}_0\|), \quad \text{with} \quad \hat{f}_n := [\hat{f}_n^m]_{m \le n} \in \mathbb{R}^{2n+1} \quad (\text{LI}) $$

For $D = 4N - 1$, there are spherical harmonics $f \in \mathcal{Y}_D$ which interpolate $y$; among those spherical harmonics, one selects the one which minimizes $(\|\hat{f}_D\|, \|\hat{f}_{D-1}\|, \dots, \|\hat{f}_0\|)$ for the lexicographical ordering in $\mathbb{R}^{D+1}$ [2, Chapter 6]. In practice, the numerical degree is $D = 3N$ and the computation is based on an orthogonal factorization of $A_N^{3N}$ under a suitable echelon form, as in [1].

## Weighted least squares

$$ \inf_{f \in \mathcal{Y}_{2N-1}} \sum_{x \in CS_N} \omega(x) |f(x) - y(x)|^2, \quad \text{with} \quad \omega: CS_N \to (0, +\infty) \quad (\text{WLS}) $$

We assume that $\omega$ is the interpolatory quadrature weight from [3]. The choice $D = 2N - 1$ is a numerical guarantee that the problem has a unique solution and a low condition number [4]. The normal equation associated to (WLS) is $(A_N^{2N-1})^\top \Omega A_N^{2N-1} \hat{f} = (A_N^{2N-1})^\top \Omega y$, with $\Omega = \text{diag}(\omega(x), x \in CS_N)$. ($\star$)

## Hyperinterpolation

$$ f \in \mathcal{Y}_{2N-1}, \quad \text{with} \quad \hat{f} = (A_N^{2N-1})^\top \Omega y \quad (\text{HI}) $$

In this way, $\hat{f}_n^m \approx \int_{S^2} F(x) Y_n^m(x) d\sigma$ if $y = F|_{CS_N}$ and the quadrature rule $\omega$ is accurate enough, as in [5]. The approximation (HI) can be understood as the zero-th order term of a Neumann series corresponding to ($\star$).

## Approximation of test functions

<figure><img src="image_3.png" /><figcaption>Fig. Four zonal functions, with axis $u = \frac{1}{\sqrt{14}}(1, 2, 3)$, and with various smoothness properties. (The grid $CS_4$ is superposed on the plot.)</figcaption></figure>

## Numerical results

<figure><img src="image_4.png" /><figcaption>Fig. Comparison of (LI), (WLS) and (HI) on the four test functions.</figcaption></figure>

## Properties of the normal matrix

$$ (A_N^{2N-1})^\top \Omega A_N^{2N-1} = I - E_N $$

* $E_N$ is block diagonal (12 blocks)

* $E_N$ is symmetric, sparse and “small”

<table><thead><tr><th>N</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>16</th><th>24</th><th>32</th></tr></thead><tbody><tr><th>Sparsity score (%)</th><td>0</td><td>0</td><td>0</td><td>0</td><td>6.5</td><td>5.1</td><td>6.9</td><td>5.9</td><td>7.1</td><td>6.3</td><td>7.2</td><td>6.6</td><td>6.8</td><td>6.8</td><td>6.6</td></tr><tr><th>Spectral radius</th><td>4e-16</td><td>7e-16</td><td>2e-15</td><td>2e-15</td><td>2e-2</td><td>9e-3</td><td>8e-3</td><td>9e-3</td><td>7e-3</td><td>7e-3</td><td>4e-3</td><td>4e-3</td><td>3e-3</td><td>1e-3</td><td>9e-4</td></tr></tbody></table>

Tab. Sparsity score (with threshold $10^{-14}$), and spectral radius of $E_N$.

## Convergence of the Neumann series for (WLS)

<figure><img src="image_5.png" /><figcaption>Fig. Approximation by a k-th order Neumann series: spectral radius of $[(A_N^{2N-1})^\top \Omega A_N^{2N-1}]^{-1} - (I + E_N + \dots + (E_N)^k)$ (thresholded at $10^{-17}$).</figcaption></figure>

## Theorem (Optimal quadrature rule on low resolution grids)

Let $1 \le N \le 4$. There is exactly one grid function $\omega: CS_N \to \mathbb{R}$ that defines an optimal quadrature rule with respect to the degree of accuracy. It has the octahedral symmetry, degree $4N - 1$, and is given as follows.

## Algebraic expression of optimal quadrature weights

<table><thead><tr><th>N</th><th>x<sub>1</sub></th><th>x<sub>2</sub></th><th>x<sub>3</sub></th><th>ω(x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub>)</th></tr></thead><tbody><tr><td rowspan="5">4</td><td>1/√3</td><td>1/√3</td><td>1/√3</td><td>729π/20020</td></tr><tr><td>1/(√2+s<sup>2</sup>)</td><td>1/(√2+s<sup>2</sup>)</td><td>s/(√2+s<sup>2</sup>)</td><td>2053π/51480 - 183π√2/80080</td></tr><tr><td>1/(√1+2s<sup>2</sup>)</td><td>1/(√1+2s<sup>2</sup>)</td><td>s/(√1+2s<sup>2</sup>)</td><td>2053π/51480 + 183π√2/80080</td></tr><tr><td>1/√2</td><td>1/√2</td><td>0</td><td>512π/15015</td></tr><tr><td>1/(√1+s<sup>2</sup>)</td><td>s/(√1+s<sup>2</sup>)</td><td>0</td><td>2048π/45045</td></tr><tr><td></td><td>1</td><td>0</td><td>0</td><td>736π/15015</td></tr></tbody></table>



## Conclusion and discussion

* (LI), (WLS) and (HI) look similar for the non-smooth functions $F_i$, $i \le 3$.

* (LI) and (WLS) look similar for the very smooth function $F_4$, and they outperform (HI) for $5 \le N \le 25$.

* For $1 \le N \le 4$, (WLS) reduces to (HI), because the chosen weight has the degree of accuracy $4N - 1$ (it is the optimal one).

* More generally, (WLS) can be solved efficiently by a few iterations of suitable iterative solvers, due to the structure of the normal matrix.

* For instance, compromises between (HI) and (WLS) have been achieved using truncated Neumann series.

* Perspectives include the use of fast spherical transforms [6] to reduce the computational cost.

## References

[1] J.-B. Bellet, M. Brachet, and J.-P. Croisille, Interpolation on the Cubed Sphere with Spherical Harmonics, Numerische Mathematik, 153 (2023), pp. 249-278.

[2] J.-B. Bellet. Mathematical and numerical methods for three-dimensional reflective tomography and for approximation on the sphere. Habilitation thesis, Université de Lorraine, 2023.

[3] J.-B. Bellet, M. Brachet, and J.-P. Croisille, Quadrature and symmetry on the Cubed Sphere, Journal of Computational and Applied Mathematics, 409 (2022).

[4] J.-B. Bellet and J.-P. Croisille, Least Squares Spherical Harmonics Approximation on the Cubed Sphere, Journal of Computational and Applied Mathematics, 429 (2023).

[5] C. An and H.-N. Wu. Bypassing the quadrature exactness assumption of hyperinterpolation on the sphere. Journal of Complexity, 80 (2024), 101789.

[6] S. Kunis and D. Potts. Fast spherical Fourier algorithms. Journal of Computational and Applied Mathematics, 161 (2003), pp. 75-98.

## Acknowledgements

This work was supported by the French National program LEFE (Les Enveloppes Fluides et l'Environnement).
# Comparison of several approximation schemes on the Cubed Sphere

<figure><img src="image_2.png" /></figure>

## Goals

* Compare interpolation, weighted least squares and hyperinterpolation on the Cubed Sphere grid.

* Get an efficient scheme for spherical approximation using spherical harmonics.

## Cubed Sphere grid

Equiangular Cubed Sphere with step $\frac{\pi}{2N}$, $CS_N \subset S^2$

<figure><img src="image_3.png" /></figure>

Fig. CS<sub>N</sub> (white dots) is obtained by projection of a meshed circumscribed cube (N = 3).

## Data approximation

Find $f \in \mathcal{Y}_D$ such that $f(x) \approx y(x)$, $x \in CS_N$.

* Data: $y = [y(x)]_{x \in CS_N} \in \mathbb{R}^{6N^2+2}$

* Unknown: $f \in \mathcal{Y}_D$, with Fourier coefficients $\hat{f} \in \mathbb{R}^{(D+1)^2}$

* Matrix writing: find $\hat{f} \in \mathbb{R}^{(D+1)^2}$ such that $A_N^D \hat{f} \approx y$

* Application: $F: S^2 \to \mathbb{R}$ is approximated by $f$ from $y = F|_{CS_N}$

## Lexicographical interpolation

$\inf_{\substack{f \in \mathcal{Y}_{4N-1} \\ f|_{CS_N} = y}} (\|\hat{f}_{4N-1}\|, \|\hat{f}_{4N-2}\|, \dots, \|\hat{f}_0\|)$, with $\hat{f}_n := [\hat{f}_n^m]_{m \le n} \in \mathbb{R}^{2n+1}$ (LI)

For $D = 4N - 1$, there are spherical harmonics $f \in \mathcal{Y}_D$ which interpolate $y$; among those spherical harmonics, one selects the one which minimizes $(\|\hat{f}_D\|, \|\hat{f}_{D-1}\|, \dots, \|\hat{f}_0\|)$ for the lexicographical ordering in $\mathbb{R}^{D+1}$ [2, Chapter 6]. In practice, the numerical degree is $D = 3N$ and the computation is based on an orthogonal factorization of $A_N^{3N}$ under a suitable echelon form, as in [1].

## Weighted least squares

$\inf_{f \in \mathcal{Y}_{2N-1}} \sum_{x \in CS_N} \omega(x) |f(x) - y(x)|^2$, with $\omega: CS_N \to (0, +\infty)$ (WLS)

We assume that $\omega$ is the interpolatory quadrature weight from [3]. The choice $D = 2N - 1$ is a numerical guarantee that the problem has a unique solution and a low condition number [4]. The normal equation associated to (WLS) is $(A_N^{2N-1})^\top \Omega A_N^{2N-1} \hat{f} = (A_N^{2N-1})^\top \Omega y$, with $\Omega = \text{diag}(\omega(x), x \in CS_N)$. ($\star$)

## Hyperinterpolation

$f \in \mathcal{Y}_{2N-1}$, with $\hat{f} = (A_N^{2N-1})^\top \Omega y$ (HI)

In this way, $\hat{f}_n^m \approx \int_{S^2} F(x) Y_n^m(x) d\sigma$ if $y = F|_{CS_N}$ and the quadrature rule $\omega$ is accurate enough, as in [5]. The approximation (HI) can be understood as the zero-th order term of a Neumann series corresponding to ($\star$).

## Approximation of test functions

<figure><img src="image_4.png" /><figcaption>Fig. Four zonal functions, with axis $u = \frac{1}{\sqrt{14}}(1, 2, 3)$, and with various smoothness properties. (The grid $CS_4$ is superposed on the plot.)</figcaption></figure>

## Numerical results

<figure><img src="image_5.png" /><figcaption>Fig. Comparison of (LI), (WLS) and (HI) on the four test functions.</figcaption></figure>

## Properties of the normal matrix

$$(A_N^{2N-1})^\top \Omega A_N^{2N-1} = I - E_N$$

* $E_N$ is block diagonal (12 blocks)

* $E_N$ is symmetric, sparse and “small”

<table><thead><tr><th>N</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>16</th><th>24</th><th>32</th></tr></thead><tbody><tr><th>Sparsity score (%)</th><td>0</td><td>0</td><td>0</td><td>0</td><td>6.5</td><td>5.1</td><td>6.9</td><td>5.9</td><td>7.1</td><td>6.3</td><td>7.2</td><td>6.6</td><td>6.8</td><td>6.8</td><td>6.6</td></tr><tr><th>Spectral radius</th><td>4e-16</td><td>7e-16</td><td>2e-15</td><td>2e-15</td><td>2e-2</td><td>9e-3</td><td>8e-3</td><td>9e-3</td><td>7e-3</td><td>7e-3</td><td>4e-3</td><td>4e-3</td><td>3e-3</td><td>1e-3</td><td>9e-4</td></tr></tbody></table>

Tab. Sparsity score (with threshold $10^{-14}$), and spectral radius of $E_N$.

## Convergence of the Neumann series for (WLS)

<figure><img src="image_6.png" /><figcaption>Fig. Approximation by a k-th order Neumann series: spectral radius of $[(A_N^{2N-1})^\top \Omega A_N^{2N-1}]^{-1} - (I + E_N + \dots + (E_N)^k)$ (thresholded at $10^{-17}$).</figcaption></figure>

## Theorem (Optimal quadrature rule on low resolution grids)

Let $1 \le N \le 4$. There is exactly one grid function $\omega: CS_N \to \mathbb{R}$ that defines an optimal quadrature rule with respect to the degree of accuracy. It has the octahedral symmetry, degree $4N - 1$, and is given as follows.

## Algebraic expression of optimal quadrature weights

<table><thead><tr><th>N</th><th>x<sub>1</sub></th><th>x<sub>2</sub></th><th>x<sub>3</sub></th><th>ω(x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub>)</th></tr></thead><tbody><tr><td>1</td><td>1/√3</td><td>1/√3</td><td>1/√3</td><td>π/2</td></tr><tr><td>2</td><td>1/√3</td><td>1/√3</td><td>1/√3</td><td>9π/70</td></tr><tr><td>3</td><td>1/√2</td><td>1/√2</td><td>0</td><td>16π/105</td></tr><tr><td>4</td><td>1/√2</td><td>1/√2</td><td>0</td><td>4π/21</td></tr><tr><td>5</td><td>1/√3</td><td>1/√3</td><td>1/√3</td><td>9π/140</td></tr><tr><td>6</td><td>1/√2</td><td>1/√2</td><td>1/√2</td><td>61π/840 - 3π√3/560</td></tr><tr><td>7</td><td>1/√2</td><td>1/√2</td><td>1/√2</td><td>61π/840 + 3π√3/560</td></tr></tbody></table>

<table><thead><tr><th>N</th><th>x<sub>1</sub></th><th>x<sub>2</sub></th><th>x<sub>3</sub></th><th>ω(x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub>)</th></tr></thead><tbody><tr><td>1</td><td>1/√3</td><td>1/√3</td><td>1/√3</td><td>729π/20020</td></tr><tr><td>2</td><td>1/√2 + s<sup>2</sup></td><td>1/√2 + s<sup>2</sup></td><td>s/√2 + s<sup>2</sup></td><td>2053π/51480 - 183π√2/80080</td></tr><tr><td>3</td><td>1/√2 + s<sup>2</sup></td><td>1/√2 + s<sup>2</sup></td><td>s/√2 + s<sup>2</sup></td><td>2053π/51480 + 183π√2/80080</td></tr><tr><td>4</td><td>1/√2 + s<sup>2</sup></td><td>1/√2 + s<sup>2</sup></td><td>0</td><td>512π/15015</td></tr><tr><td>5</td><td>1/√2 + s<sup>2</sup></td><td>1/√2 + s<sup>2</sup></td><td>0</td><td>2048π/45045</td></tr><tr><td>6</td><td>1/√2 + s<sup>2</sup></td><td>1/√2 + s<sup>2</sup></td><td>0</td><td>736π/15015</td></tr></tbody></table>

with $s = \sqrt{2} - 1$

## Conclusion and discussion

* (LI), (WLS) and (HI) look similar for the non-smooth functions $F_i$, $i \le 3$.

* (LI) and (WLS) look similar for the very smooth function $F_4$, and they outperform (HI) for $5 \le N \le 25$.

* For $1 \le N \le 4$, (WLS) reduces to (HI), because the chosen weight has the degree of accuracy $4N - 1$ (it is the optimal one).

* More generally, (WLS) can be solved efficiently by a few iterations of suitable iterative solvers, due to the structure of the normal matrix.

* For instance, compromises between (HI) and (WLS) have been achieved using truncated Neumann series.

* Perspectives include the use of fast spherical transforms [6] to reduce the computational cost.

## References

[1] J.-B. Bellet, M. Brachet, and J.-P. Croisille, Interpolation on the Cubed Sphere with Spherical Harmonics, Numerische Mathematik, 153 (2023), pp. 249-278.

[2] J.-B. Bellet, M. Brachet, and J.-P. Croisille, Mathematical and numerical methods for three-dimensional reflective tomography and for approximation on the sphere. Habilitation thesis, Université de Lorraine, 2023.

[3] J.-B. Bellet, M. Brachet, and J.-P. Croisille, Quadrature and symmetry on the Cubed Sphere, Journal of Computational and Applied Mathematics, 409 (2022).

[4] J.-B. Bellet and J.-P. Croisille, Least Squares Spherical Harmonics Approximation on the Cubed Sphere, Journal of Computational and Applied Mathematics, 429 (2023).

[5] C. An and H.-N. Wu. Bypassing the quadrature exactness assumption of hyperinterpolation on the sphere. Journal of Complexity, 80 (2024), 101789.

[6] S. Kunis and D. Potts. Fast spherical Fourier algorithms. Journal of Computational and Applied Mathematics, 161 (2003), pp. 75-98.

## Acknowledgements

This work was supported by the French National program LEFE (Les Enveloppes Fluides et l'Environnement).