\documentclass{beamer}

\input{hal_prefs.tex}

\title[Quadrature/Hyperinterp. over Spherical Polygons.\hspace{2em}\insertframenumber/
\inserttotalframenumber]{Numerical Quadrature and Hyperinterpolation over Spherical Triangles/Polygons by the dCATCH Package.}

\author{Alvise Sommariva (with M. Dessole and M. Vianello)}
\institute{Universit\`a degli Studi di Padova}
\date{June 23, 2021}

\begin{document}

\frame{\titlepage}

\begin{frame}[fragile]
\frametitle{Purpose}
The main purpose of the talk is to show a novel strategy for determining
by Caratheodory-Tchakaloff compression implemented in the d-CATCH package
\begin{itemize}
\item the {\purples{computation of nodes and weights}} of a low-cardinality positive quadrature formula, nearly exact for polynomials of a given degree,

\item {\purples{hyperinterpolation}} of mild degree,
\end{itemize}

\noindent over spherical polygons.

\smsp

The {\purples{Matlab software}} used in this talk is available at the authors' homepage.

\end{frame}







\begin{frame}[fragile]
\frametitle{Spherical polygon}

A {\reds{great circle}} is the intersection of the unit-sphere ${\mathbb{S}}^2$ with a plane passing through the origin.

\smsp

Let $P_1, \ldots,P_L$ be distinct points of ${\mathbb{S}}^2$ and set $P_0=P_L$. 

\smsp

A {\reds{spherical polygon}} is the region ${\cal{P}} \subset {\mathbb{S}}^2$ whose boundary $\delta {\cal{P}}$ 
\begin{itemize}
\item is determined by the geodesic arcs $\{\gamma_k\}_{k=0,\ldots,L}$, where each $\gamma_k$ is the portion of the great circles joining $P_k$ with $P_{k+1}$. 
\item is oriented counterclockwise.
\end{itemize}

\smsp

In this talk we suppose that ${\cal{P}}$ 

\begin{itemize}
\item is contained in a cap $\Omega$ whose {\sl{polar angle}} is strictly inferior than $\pi$;
\item $\cal{P}$ is {\purples{simple}}, i.e. it has no self intersections.
\end{itemize}


\end{frame}






\begin{frame}[fragile]
\frametitle{Spherical polygon triangulation}
Under these assumptions, 
\begin{itemize}
\item if $C$ is the center of such a cap $\Omega$ containing ${\cal{P}}$ then determine the {\purples{stereographic projection ${\cal{P}}^{\prime}$ of ${\cal{P}}$ on the plane $\pi_C$ tangent in $C$}} to the unit-sphere;
\item {\purples{compute a triangulation over the planar polygon ${\cal{P}}^{\prime}$}} (e.g. by Matlab built-in environment {\tt{polyshape}}), i.e. ${\cal{P}}^{\prime}=\cup_{i=1}^m T^{\prime}_i$ where $T^{\prime}_i \subset \pi_C$ are planar triangles whose interior do not overlap, i.e. if $j \neq k$ then ${\mbox{int}}(T^{\prime}_j) \cap {\mbox{int}}(T^{\prime}_k)=\emptyset$;
\item {\purples{map back to the sphere}}, by means of the inverse of the stereographic projection, each planar triangle $T^{\prime}_k$ into a spherical triangle $T_k$.
\end{itemize}
The spherical triangles $\{T_k\}_{k=1,\ldots,M}$ determine a spherical triangulation of ${\cal{P}}$.
\end{frame}






\begin{frame}[fragile]
\frametitle{Cubature on spherical polygons}
The purpose of this section is to determine a {\reds{cubature rule over the spherical polygon}} ${\cal{P}}$, that has {\purples{internal nodes}} $\{Q_k\}_{k=1,\ldots,N}$ and {\purples{positive weights}} $\{w_k\}_{k=1,\ldots,N}$ so that if $f$ is a continuous function then
$$
\int_{{\cal{P}}} f(x,y,z)d\sigma  \approx \sum_{k=1}^N w_k f(x_k,y_k,z_k).
$$

Since $\{T_k\}_{k=1,\ldots,M}$ is a spherical triangulation of ${\cal{P}}$, if we determine such a {\reds{cubature rule over each spherical triangle}} $T_k$ then by the addivity of integration we have such a rule on ${\cal{P}}$ either.
\end{frame}





\begin{frame}[fragile]
\frametitle{Cubature on spherical polygons: spherical triangles}

{\vspace{-0.3cm}}

\begin{figure}\label{1}
\centering
{\includegraphics[scale=0.4,clip]{umbrella.eps}}
 \caption{What we will get: quadrature nodes (and weights) on a spherical triangle lifted from the projected elliptical triangle.}
 \end{figure}

\end{frame}



\begin{frame}[fragile]
\frametitle{Cubature on spherical polygons: spherical triangles}
With no loss of generality, up to a suitable rotation, we concentrate on spherical triangles ${\cal{T}} = ABC$ with centroid $(A+B+C)/\|A+B+C\|_2$ at the north pole, that are completely contained in the northern-hemisphere, and do not touch the equator. 

\smsp

Then, if $f \in C({\cal{T}})$, $g(x,y)=\sqrt{1-x^2-y^2}$,
$$I_{{\cal{T}}}:=\int_{\cal{T}} f(x,y,z)d\sigma =  \int_{{\cal{T}}^{\perp}} f(x,y,g(x,y)) \frac{1}{g(x,y)}dxdy, $$
where ${\cal{T}}^{\perp}$ is the projection of ${\cal{T}}$ onto the xy-plane, that is the curvilinear triangle whose vertices, say $\hat{A}$, $\hat{B}$, $\hat{C}$, are the xy-coordinates of ${A}$, ${B}$, ${C}$, respectively.
\end{frame}




\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: spherical triangles}
{\small{
The sides of ${\cal{T}}^{\perp}$ are arcs of ellipses centered at the origin, being the projections of great circle arcs. Then we can split the planar integral into the sum of the integrals on three elliptical sectors $S_i$ with $i=1,2,3$, obtained by joining the origin with the vertices $\hat{A}$, $\hat{B}$, $\hat{C}$, namely
$$I_{{\cal{T}}}=\int_{{\cal{T}}^{\perp}} f(x,y,g(x,y)) \frac{1}{g(x,y)}dxdy=\sum_{i=1}^3 \int_{S_i} f(x,y,g(x,y)) \frac{1}{g(x,y)}dxdy
$$

\smsp

If the purpose is to compute an algebraic rule over $T$, with degree of precision $n$, positive weights and internal nodes (i.e. rules of PI-type) we face the problem that while $f$ is a polynomial of total degree $n$ then being $g(x,y)=\sqrt{1-x^2-y^2}$, we have that in general $f(x,y,g(x,y)) \frac{1}{g(x,y)}$ may not be a polynomial.
}}
\end{frame}


\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: spherical triangles}
To see this properly let 
\begin{itemize}
\item $
f(x,y,z)=x^{\alpha} y^{\beta} z^{\gamma}, \quad 0 \leq {\alpha} + {\beta} + {\gamma} \leq n, \,\, {\alpha}, {\beta}, {\gamma} \in \N
$;
\item $g(x,y)=\sqrt{1-x^2-y^2}$.
\end{itemize}
Thus
$$
f(x,y,g(x,y)) \frac{1}{g(x,y)}=x^{\alpha} y^{\beta} (1-x^2-y^2)^{(\gamma-1)/2}.
$$
Thus if $\gamma$ is 
\begin{itemize}
\item {\purples{odd}} then $f(x,y,g(x,y)) \frac{1}{g(x,y)}$ is a {\reds{polynomial}} of degree at most $n$,
\item {\purples{even}} then $f(x,y,g(x,y)) \frac{1}{g(x,y)}$ is {\reds{$1/g$ multiplied for a polynomial}} of degree at most $n$.
\end{itemize}
\end{frame}


\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: spherical triangles}
Key points: 
\begin{itemize}
\item {\purples{Approximate $1/g$ by a polynomial $p_{\epsilon}$ of degree $m=m_{\epsilon}$}} such that $|p_{\epsilon}-1/g| \leq \epsilon \cdot 1/|g|$.
\item Since $f/g \approx f \cdot p_{\epsilon} \in {\mathbb{P}}_{n+m}$, we {\purples{integrate $f \cdot p_{\epsilon}$}} instead of $f/g$ on the elliptical sectors $S_i$, $i=1,2,3$.
\end{itemize}

By {\purples{determining a rule of algebraic degree of precision $n+m$ over each elliptical sector $S_i$, $i=1,2,3$}}, with internal nodes, and positive weights then we have a rule on ${\cal{T}}^{\perp}:=\cup_{i=1}^3 S_i$ with nodes $(x_k,y_k)_{k=1,\ldots,N_{n+m}}$, weights $w_{k=1,\ldots,N_{n+m}}$ of PI-type.

\vspace{0.15cm}

Mapping back the nodes on the sphere, we have a rule over the spherical triangle that is near algebraic with ADE $n$ since
{\purples{
\begin{equation}
\int_{\cal T} f(x,y,z) d\sigma \approx \sum_{j=1}^{N_{n+m}} \frac{w_i}{\sqrt{1-x_j^2-y_j^2}} f(x_j,y_j,\sqrt{1-x_j^2-y_j^2}).
\end{equation}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: spherical triangles}
Some details: 
\begin{itemize}
\item {\purples{Rule over elliptical sector $S_i$}}: since each $S_i$ is an affine transformation of a circular sector of the unit-disk $S^*_i$, we determine a formula on $S^*_i$ and map it to $S_i$ (some care on the weights that must be multiplied by absolute value of the transformation matrix determinant);
\item {\purples{Computation of $m=m_{\epsilon}$}}: it is sufficient to find the degree of a (near) optimal univariate polynomial approximation (up to $\epsilon$) to  
$
{\reds{1/\sqrt{1-t}}}$ where $t \in [0,\max\{\|{\hat{A}}\|_2^2,\|{\hat{B}}\|_2^2,\|{\hat{C}}\|_2^2\}]
$. 

\vspace{0.15cm}

Thus $m=m_{\epsilon}$ can be estimated by Chebfun (even stored in tables!).
\item {\purples{Caratheodory-Tchakaloff rule compression}}:  from the nodes $\{P_k\}$ and weights $\{w_k\}$ of the PI-type rule on $\cal{T}$, we extract one with cardinality at most $(n+1)^2$, by Lawson-Hanson algorithm (see {\tt{dCATCH}} suite implementation).
\end{itemize}

\end{frame}









\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: example}

{\vspace{-0.3cm}}

\begin{figure}\label{2}
\centering
{\includegraphics[scale=0.4,clip]{america_cub_ADE8_sphF.eps}}
 \caption{Quadrature nodes on a spherical polygon of a rule of PI-type, with ADE=8, 380544 points, before Caratheodory-Tchakaloff compression. Cputime: about 10 seconds.}
 \end{figure}

\end{frame}

\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: example}

\smnsp
\smnsp

\begin{figure}\label{3}
\centering
{\includegraphics[scale=0.4,clip]{america_cub_ADE8_sphC.eps}}
 \caption{Triangulation of a spherical polygon (967 spherical triangles) and quadrature rule of PI-type with ADE=8, 81 points, after Caratheodory-Tchakaloff compression (magenta). Cputime: 3.5 seconds.}
 \end{figure}

\end{frame}








\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: hyperinterpolation}
As introduced by I.H. Sloan in 1995, given
\begin{itemize}
\item an {\purples{orthonormal basis of ${\mathbb{P}}^d_n(\Omega)$}} (the subspace of $d$-variate polynomials of total-degree not exceeding $n$, restricted to a compact set or manifold $\Omega \subset R^d$) {\purples{w.r.t. a given measure}} $d \mu$ on $\Omega$), say $\{p_j\}$, $1 \leq j \leq N_n = {\mbox{ dim}} ({\mathbb{P}}^n_d (\Omega))$, 
\item a {\purples{quadrature formula}} exact for ${\mathbb{P}}^{\reds{2n}}_d (\Omega)$ with nodes $X = \{x_i\} \subset \Omega$ and positive weights ${\bf{w}} = \{w_i\}$, $1 \leq i \leq M$ with $M \geq  N_n$, 
\end{itemize}
the discretized orthogonal projection (hyperinterpolation) of $f \in C(\Omega)$ is
\begin{eqnarray}
({\cal{L}}_n f)(x) &=& \sum_{j=1}^M \langle f,p_j \rangle_{l_{2,{\bf{w}}}(X)} p_j(x) =  {\purples{\sum_{i=1}^M  w_i f(x_i) \sum_{j=1}^{N_n} p_j(x_i) p_j(x)}}. \nonumber
\end{eqnarray}
\end{frame}


\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: hyperinterpolation}
Letting $\Omega$ be the spherical polygon,
\begin{itemize}
\item by means of the routines in our Matlab package {\tt{dCATCH}} we determine the required {\purples{orthonormal basis}} $\{p_j\}$,
\item we {\purples{compute the quadrature rule}} of degree {\reds{$2n$}} on $\Omega$,
\item we finally get the hyperinterpolant of degree {\reds{$n$}}.
\end{itemize}

\smsp

As previously mentioned, the Matlab software implementing this approach is available at the authors' homepage.


\end{frame}




\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: hyperinterpolation test}

\smnsp

As example, we take into account a course map of Australia (without taking into account its smaller islands).

 
\begin{figure}\label{4}
\centering
{\includegraphics[scale=0.3,clip]{australia_island_cub_ADE10_2D.eps}}

\smnsp
\smnsp

 \caption{Quadrature nodes of PI-type on a coarse approximation of Australia, useful for hyperinterpolation. ADE is 10, there are 167 sph. triangles, the full rule has 81711 points, the compressed one 121, with moments error of  $\approx 5\cdot 10^{-15}$. The whole process took about 3 seconds.}
 \end{figure}

\end{frame}


\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: hyperinterpolation test}

\smnsp
Setting $(x_0,y_0,z_0)\approx (-6.325e-01,6.668e-01,-3.908e-01)$ as an approximation of australian centroid, and $$h(x,y,z) = (x-x_0)^2+(y-y_0)^2+(z-z_0)^2$$ we define

\begin{enumerate}
\item {\dkgreens{$f_1(x,y,z)=1+x+y^2+x^2 \cdot y+x^4+y^5+x^2 \cdot y^2 \cdot z^2$}} (polynomial, total degree 6);
\item {\reds{$f_2(x,y,z)=\cos(10\cdot (x+y+z))$}} (oscillations);
\item {\dkgreens{$f_3(x,y,z)=\sin(-h(x,y,z))$}} (regular);
\item {\dkgreens{$f_4(x,y,z)=\exp(-h(x,y,z))$}} (regular);
\item {\reds{$f_5(x,y,z)=((x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2)^{3/2}$}} (hard);
\item {\purples{$f_6(x,y,z)=((x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2)^{5/2}$}} (medium);
\end{enumerate}
\end{frame}




\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: hyperinterpolation test}

\smnsp

 
\begin{figure}\label{hL2a}
\centering
{\includegraphics[scale=0.4,clip]{hyp_INF_err_aus}}
 \caption{Inf-Norm hyperinterpolation error.}
 \end{figure}

\end{frame}

\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: hyperinterpolation test}
\begin{figure}\label{hL2b}
\centering
{\includegraphics[scale=0.4,clip]{hyp_norm_aus.eps}}
 \caption{Hyperinterpolation operator Inf-Norm.}
 \end{figure}
\end{frame}



\begin{frame}[fragile]
\frametitle{Bibliography}

{\scriptsize{
\begin{itemize}
\item {\purples{On quadrature rules over spherical triangles:}} K. Atkinson, {\sl{Numerical integration on the sphere}}, J. Austral. Math. Soc. Ser. B 23 (1981/82), 332--347.

\item {\purples{On quadrature rules over spherical triangles:}} A. Sommariva and M. Vianello, {\sl{Near-algebraic Tchakaloff-like quadrature on spherical triangles}}, Appl. Math. Lett. (120) 2021.

\item {\purples{Lawson-Hanson seminal work:}}  C.L. Lawson, R.J. Hanson,  {\sl{Solving Least Squares Problems}}, Classics in Applied Mathematics 15, SIAM, Philadelphia, 1995.

\item {\purples{On Lawson-Hanson algorithm:}} M. Dessole, F. Marcuzzi and M. Vianello, {\sl{Accelerating the Lawson-Hanson NNLS solver for large-scale
Tchakaloff regression designs}}, Dolomites Res. Notes Approx. DRNA 13
(2020), 20--29.

\item {\purples{On spherical triangulation:}} G.Barequet, M.Dickerson, D.Eppstein, {\sl{On triangulating three-dimensional polygons}}, Comp. Geometry 10 (1998) 155--170.

\item {\purples{Hyperinterpolation seminal paper:}} I.H. Sloan, {\sl{Interpolation and hyperinterpolation over general regions}}, J. Approx. Theory 83 (1995), 238--254.

\item {\purples{Details on hyperinterpolation implementation:}} A. Sommariva and M. Vianello, {\sl{Numerical hyperinterpolation over nonstandard planar regions}}, Mathematics and Computers in Simulation, Volume 141, 10 (2017) 110--120.

\item {\purples{Matlab software:}} A. Sommariva homepage, {\url{https://www.math.unipd.it/~alvise/software.html}}.
\end{itemize}
}}




\end{frame}

\end{document}



\begin{frame}[fragile]
\frametitle{Quadrature on spherical polygons: hyperinterpolation}

As soon as an algebraic quadrature formula is at hand, it is possible to reconstruct a continuous function from its values at the quadrature nodes by {\purples{hyperinterpolation}}, a technique that avoids the inherent difficulties of multivariate polynomial interpolation.

\smsp

Introduced by I.H. Sloan in 1995, hyperinterpolation is essentially an {\purples{orthogonal (Fourier-like) projection on polynomial spaces}}, w.r.t. the discrete measure associated with a positive algebraic quadrature formula, or in other words a weighted least-squares polynomial approximation at the quadrature nodes. 

\smsp

In the last twenty years the subject has been developed and {\purples{extended}} from the sphere to several d-dimensional domains, such as cubes and balls but also less standard ones, from both the theoretical and the modelling/computational point of view.

\end{frame}

