\documentclass{beamer}

\input{hal_prefs.tex}
\usepackage{multirow}

\newcommand{\minm}{\min\limits_{p \in \mathbb{P}^d_n(\Omega)}}
\newcommand{\leastapprox}{\sum\limits_{j=1}^{M_{2n}}w_{j}(p({x}_j) - f({x}_j))^2}
\newcommand{\filteredhyper}{\cal{F}_{L,N}}

\newcommand{\lassohyper}{\mathbb{L}_{L}^{\lambda}}
% \newcommand{\filteredhyper}{\mathbb{F}_{L,N}}
\newcommand{\hyper}{\mathbb{L}_{L}}
\newcommand{\hybridhyper}{\mathbb{H}_{L}^{\lambda}}
\newcommand{\ltwohyper}{\mathbb{P}_{L}^{\lambda_2}}
\newcommand{\discreteinner}[1]{\left\langle#1\right\rangle_N}
\newcommand{\inner}[1]{\left\langle#1\right\rangle}
\newcommand{\brackets}[1]{\left(#1\right)}
\newcommand{\squarebrackets}[1]{\left[#1\right]}
\newcommand{\braces}[1]{\left\{#1\right\}}
\newcommand{\absolute}[1]{\left|#1\right|}
\newcommand{\softoperatordef}[1]{\mathbb{S}_{k}\left({#1}\right)}
\newcommand{\softoperatornew}[1]{\mathbb{S}_{\lambda\mu_{\ell}}\left(#1\right)}
\newcommand{\softoperator}[1]{\mathbb{S}_{\lambda\mu_{\ell}}\left(\discreteinner{#1}\right)}
\newcommand{\hybridoperator}[1]{\frac{\mathbb{S}_{\lambda\mu_{\ell}}\left(\discreteinner{#1}\right)}{1+\lambda_2b_{\ell}^2}}
\newcommand{\ltwooperator}[1]{\mathbb{R}_{L}{#1}}
\newcommand{\ltwooperatornew}{\mathbb{R}_{L}}
\newcommand{\ltworegmatrix}{\normtwo{\mathbf{W}^{1/2}\mathbf{R}_2\bm{\beta}}^2}
\newcommand{\loneregmatrix}{\normone{\mathbf{R}_1\bm{\beta}}}

\newcommand{\leastapproxmatrix}{{\|\mathbf{W}^{1/2} ( \mathbf{A}\bm{\beta} - \mathbf{f})\|_2^2}}
\newcommand{\measure}{\text{d}\omega}
\newcommand{\prob}{ \Phi_{\lambda_1, \lambda_2} (\bm{\beta})}
\newcommand{\probn}{ \Phi_{\lambda_{1,n}, \lambda_{2,n}} (\bm{\beta}^{n})}

\newcommand{\minmm}{\min\limits_{\bm{\beta}\in\mathbb{R}^{d}}}
\newcommand{\summ}[2]{\sum\limits_{#1 =1}^{#2}}
\newcommand{\norminfty}[1]{\|#1\|_{\infty}}
\newcommand{\normtwo}[1]{\|#1\|_2}
\newcommand{\normone}[1]{\|#1\|_1}
\newcommand{\normzero}[1]{\|#1\|_0}
\newcommand{\acronym}[1]{\text{#1}}
\newcommand{\thyper}{\text{hyperinterpolation}}
\newcommand{\Thyper}{\text{Hyperinterpolation}}


\title[On hyperinterpolation and its variants.\hspace{2em}\insertframenumber/
\inserttotalframenumber]{On hyperinterpolation and its variants}

\author{Alvise Sommariva}
\institute{Universit\`a degli Studi di Padova}
\date{{\footnotesize{A two days on Approximation Theory and Applications, \\ Conference in honor of Paolo Emilio Ricci on his 80th Birthday, \\ Rome (I), May 23-24, 2024}}}

\begin{document}

\frame{\titlepage}

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

The purpose of this talk is to

\begin{itemize}
\item recall the concept of the {\purples{classical hyperinterpolation}} and its basic properties;
\item outline the {\purples{main variants}} and their features;
\item show some {\purples{numerical examples}} and discuss the advantages of the techniques mentioned above;
\item conclude by mentioning some {\purples{new frontiers}} on this topic.
\end{itemize} 

\end{frame}





\begin{frame}[fragile]
\frametitle{Classical hyperinterpolation}

The concept of hyperinterpolation has been {\purples{introduced by I. Sloan}} in the paper 
\begin{center}
{\blues{Interpolation and hyperinterpolation over general regions}}
\end{center}
published in 1995 on {\it{Journal of Approx. Theory}}.

{\vspace{0.3cm}}

\begin{itemize}
\item It 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. 
\item It is called {\it{hyperinterpolation}}, since {\purples{under special assumptions it corresponds to interpolation}}.
\end{itemize}
\end{frame}





\begin{frame}[fragile]
\frametitle{Classical hyperinterpolation}
Let ${\mathbb{P}}^d_n(\Omega)$ be the vector space {\purples{$d$-variate polynomials of total-degree not exceeding $n$}}, restricted to a compact set or manifold $\Omega \subset R^d$.

\smsp

Given
\begin{itemize}
\item an {\purples{orthonormal basis}} $\{p_j\}$, $1 \leq j \leq N_n = {\mbox{ dim}} ({\mathbb{P}}^d_n (\Omega))$ of ${\mathbb{P}}^d_n(\Omega)$ w.r.t. an absolutely continuous measure $d \mu$ on $\Omega$,
\item a {\purples{quadrature formula}} exact for ${\mathbb{P}}^d_{\reds{2n}} (\Omega)$ w.r.t. $\mu$, with nodes $X = \{x_i\} \subset \Omega$ and positive weights $w_i$, $1 \leq i \leq M_{2n}$ with $M_{2n} \geq  N_n$,
\end{itemize}
the discretized orthogonal projection (hyperinterpolation) of $f \in C(\Omega)$ of degree $n$ is
\begin{eqnarray}
({\cal{L}}_n f)(x) &=& \sum_{j=1}^{M_{2n}} \langle f,p_j \rangle_{l_{2,{\bf{w}}}(X)}  p_j(x) =  {\purples{\sum_{i=1}^{M_{2n}}  w_i f(x_i) \sum_{j=1}^{N_n} p_j(x_i) p_j(x)}}. \nonumber
\end{eqnarray}
where 
$$
\langle f,g \rangle_{l_{2,{\bf{w}}}(X)} := \sum_{i=1}^{M_{2n}} w_i f(x_i) g(x_i).
$$
\end{frame}




\begin{frame}[fragile]
\frametitle{Classical hyperinterpolation}
I. H. Sloan proved the following results:

\smsp

\begin{itemize}
\item [1.] {\color{blue}{(LS)}} The polynomial ${\cal{L}}_n f \in {\mathbb{P}}^d_n(\Omega)$ is the solution to the {\purples{LS problem}}
\begin{equation}\label{HyperLeastApprox}
\minm {\frac{1}{2}} \langle f,p \rangle_{l_{2,{\bf{w}}}}  =\minm {\frac{1}{2}}  \leastapprox ,
\end{equation}

\item [2.]  {\color{blue}{(Operator norm)}}  Given $f \in C(\Omega)$,  $\|{\cal{L}}_n f\|_2 \leq {\sqrt{\mu(\Omega)}} \| f\|_{\infty}$,

\item [3.] {\color{blue}{(Error estimate)}}  If $f \in C(\Omega)$ then $$\|{\cal{L}}_n f-f\|_{{\reds{2}}} \leq 2 \sqrt{\mu(\Omega) }E_n(f,\Omega) $$ where $E_n(f,\Omega)=\min_{p \in P^d_n(\Omega)} \|f-p\|_{\reds{\infty}} $.
\end{itemize}

\end{frame}

\begin{frame}[fragile]
\frametitle{Classical hyperinterpolation}
\smsp

\begin{Definition}
An algebraic quadrature rule 
$$
\int_{\Omega} f(x) d\Omega \approx \sum_{i=1}^{M_{2n}} w_i f(x_i)
$$

\noindent is {\purples{minimal}} if it has algebraic degree of precision $2n$ and $M_{2n}={\mbox{dim}} ({\mathbb{P}}^d_{n}):=N_n$ 
\end{Definition}

\begin{Examples}
\begin{itemize}
\item [A.] Gaussian rules in the interval;
\item [B.] Angular equispaced rules in the circle;
\item [C.] No such rules in general on the sphere ($n \geq 3$).
\end{itemize}
\end{Examples}

\end{frame}

\begin{frame}[fragile]
\frametitle{Classical hyperinterpolation}
\smsp

\begin{itemize}
\item [4.] {\color{blue}{(Interpolation)}} The classical interpolation formula
$$
{\cal{L}}_n f(x_j)=f(x_j), \,\, j=1,\ldots,N_n
$$
holds for arbitrary $f$ if and only if the quadrature rule is minimal.
\item [5.]  {\color{blue}{(Projection)}}  If $p \in {\mathbb{P}}^d_{n}$ then ${\cal{L}}_n p \equiv p$.
\end{itemize}

\end{frame}


\begin{frame}[fragile]
\frametitle{Filtered hyperinterpolation: introduction}
\smsp
{\purples{Filtered hyperinterpolation}} was introduced by I.H.Sloan and R.S.Womersley on the unit sphere, in the paper published in 2012
\begin{center}
{\blues{Filtered hyperinterpolation: A constructive polynomial approximation on the sphere}}
\end{center}

\smsp

In the introduction they wrote:

{\scriptsize{
\begin{quote}
While there are many other ways (such as radial
basis functions, spline functions on a triangular mesh, wavelets) of representing scalar
physical quantities on a sphere, {\purples{polynomials}} continue to play an important role. 

{\vspace{0.2cm}}

Of course their resolving power is limited by the degree of the polynomial, but the principal
complaint about polynomials is often that {\purples{irregularities or discontinuities in the
quantity being approximated lead to fringes or oscillations}} (Gibbs’ phenomenon for
the case of discontinuities). 

{\vspace{0.2cm}}

A recognised way of {\purples{reducing the effect of fringes
is by filtering an initial polynomial approximation}}, that is by modifying the amplitude
of selected Fourier modes
\end{quote}
}}



\end{frame}


\begin{frame}[fragile]
\frametitle{Filtered hyperinterpolation: filters}

\begin{itemize}
\item They introduce a {\purples{filter function}} $h \in {C}([0,+\infty))$ that satisfies
\begin{equation*}
h(x)=\begin{cases}
1, &\text{ for }x\in[0,1/2],\\
0, &\text{ for }x\in[1,\infty).
\end{cases}
\end{equation*}
\item It is immediate to observe that depending on the behaviour in $[1/2,1]$, one can define {\purples{many filters}}. 
\item Example:  choose as filter $h \in {C}([0,+\infty))$ the function
\begin{equation}\label{equ:filtered}
h(x) = \left\{\begin{array}{cc}
1, & x \in [0,\frac{1}{2}] ,\\
\sin^2(\pi x), & x \in [\frac{1}{2},1] \\
0, &\text{ for }x\in[1,\infty).
\end{array}
\right.
\end{equation}

\end{itemize}
\end{frame} 

\begin{frame}[fragile]
\frametitle{Filtered hyperinterpolation: filters}

 \begin{figure}[!ht]\label{filter}
 \centering
   {\includegraphics[scale=0.40,clip]{filter.eps}}
\smnsp      
 \caption{The filter $h$ that is $\sin^2(\pi x)$,  for $x \in [\frac{1}{2},1]$.}
 \end{figure}

\end{frame} 


\begin{frame}[fragile]
\frametitle{Filtered hyperinterpolation: definition}

\smnsp

{\small{
\begin{definition}
Consider the quadrature rule 

\smnsp

$$
\int_{\Omega} f(x) d\Omega \approx \sum_{j=1}^{M_{\delta}} w_j f(x_j)
$$

\smnsp

\noindent and suppose that
\begin{itemize}
\item  $w_{j} > 0$, $x_j \in \Omega$, $j=1,\ldots,M_{\delta}$ (i.e. it is PI-type rule);
\item it has algebraic degree of exactness  {\reds{$\delta=n-1+\lfloor n/2 \rfloor$}}.
\end{itemize}

Introduce the discrete scalar product determined by such formula

\smnsp
\smnsp

$$
\langle f,g \rangle_{M_{\delta}}:= \sum_{j=1}^{M_{\delta}} w_j f(x_j) g(x_j)
$$

\smnsp

The {\purples{filtered hyperinterpolant}} ${\cal{F}}_n f \in \mathbb{P}^d_{\reds{n-1}}(\Omega)$ of $f \in C(\Omega)$ is defined as
\begin{equation}\label{equ:filteredhyper}
{\cal{F}}_n f:=\sum_{k=1}^{N_{n}} {\reds {h \left( \frac{deg(p_k)}{n}  \right) }} \langle f,p_k \rangle_{M_{\delta}} \, p_k
\end{equation}
where $N_n={\mbox{dim}}({\mathbb{P}}^d_n(\Omega))$, $\{p_k\}$ family of orthonormal polynomials*.
\end{definition}
}}
\end{frame} 


\begin{frame}[fragile]
\frametitle{Filtered hyperinterpolation: remark}

\smnsp

\begin{remark}
Since 

\begin{equation}
{\cal{F}}_n f:=\sum_{k=1}^{N_n} {\reds {h \left( \frac{deg(p_k)}{n}  \right) }} \langle f,p_k \rangle_{M_{\delta}} \, p_k
\end{equation}

and 

\begin{equation*}
h(x)=\begin{cases}
1, &\text{ for }x\in[0,1/2],\\
0, &\text{ for }x\in[1,\infty).
\end{cases}
\end{equation*}

depending on $\overline{\text{supp} (h)} \subseteq [0,1]$ one can {\reds{achieve more {\it{sparsity}}}} in the polynomial coefficients w.r.t. classical hyperinterpolation, i.e. some less relevant discrete Fourier coefficients are dismissed. 

\end{remark}

\begin{remark}[Projection]

Observe that ${\cal{F}}_n p =p$ for all $p \in {\mathbb{P}}^d_{\lfloor n/2 \rfloor}$.

\end{remark}
\end{frame} 




\begin{frame}[fragile]
\frametitle{Lasso hyperinterpolation: soft thresholding operator}

An alternative to filtered hyperinterpolation is the so called {\purples{Lasso hyperinterpolation}} based on the concept of 
{\purples{Soft thresholding operator}}.

\begin{definition}[Soft thresholding operator]
The {\purples{soft thresholding operator}} is defined as
\[{\cal{S}}_k (a):=\max(0,a-k)+\min(0,a+k),\]
where $k \geq 0$.
\end{definition}
Alternatively, we can define ${\cal{S}}_k (a)$ as follows
\[
{\cal{S}}_k (a)=\left\{
\begin{array}{cl}
a+k, & {\mbox{ if }} a < -k ,\\
0,  & {\mbox{ if }} -k \leq a \leq k, \\
a-k, & {\mbox{ if }}  a > k. \\
\end{array}
\right.
\]

\end{frame} 





\begin{frame}[fragile]
\frametitle{Lasso hyperinterpolation: definition}

\smnsp

{\footnotesize{
\begin{definition}
Consider the quadrature rule 

\smnsp

$$
\int_{\Omega} f(x) d\Omega \approx \sum_{j=1}^{M_{2n}} w_j f(x_j)
$$

\smnsp

\noindent and suppose that
\begin{itemize}
\item  $w_{j} > 0$, $x_j \in \Omega$, $j=1,\ldots,M_{2n}$ (i.e. it is PI-type rule);
\item it has algebraic degree of exactness  {\reds{$2n$}}.
\end{itemize}

Introduce the discrete scalar product determined by such formula

\smnsp
\smnsp

$$
\langle f,g \rangle_{N_{\delta}}:= \sum_{j=1}^{M_{2n}} w_j f(x_j) g(x_j)
$$

Then the {\purples{Lasso hyperinterpolation}} of $f$ is defined as
\begin{equation}\label{equ:lassohyper}
{\cal{L}}^{\lambda}_n{f}:=\sum_{j=1}^{N_n} {\purples{{\cal{S}}_{\lambda \mu_{j}} (\langle f, p_j \rangle_{M_{2n}})}} p_j,
\end{equation}
where ${\cal{S}}_{\lambda \mu_{j}}$, $j=1,\ldots,N_n$ are soft thresholding operators,
\begin{itemize}
\item $\lambda>0$ is the {\reds{regularization parameter}},
\item $\{\mu_{j}\}_{j=1,\ldots,N_n}$ is a set of {\reds{positive penalty parameters}}.
\end{itemize}
\end{definition}
}}
\end{frame} 


\begin{frame}[fragile]
\frametitle{Lasso hyperinterpolation: properties I}

\smnsp

Being
\begin{equation}\label{equ:lassohyper2}
{\cal{L}}^{\lambda}_n{f}:=\sum_{j=1}^{N_n} {\purples{{\cal{S}}_{\lambda \mu_{j}} (\langle f, p_j \rangle_{N_{\delta}})}} p_j,
\end{equation}
the effect of the soft threshold operator 
\[
{\cal{S}}_k (a)=\left\{
\begin{array}{cl}
a+k, & {\mbox{ if }} a < -k ,\\
0,  & {\mbox{ if }} -k \leq a \leq k, \\
a-k, & {\mbox{ if }}  a > k. \\
\end{array}
\right.
\]
is such that 
$$|\langle f, p_j \rangle_{M_{2n}} | \leq \lambda \mu_j \Rightarrow {\purples{{\cal{S}}_{\lambda \mu_{j}} (\langle f, p_j \rangle_{M_{2n}})}}=0$$
that is {\purples{more sparsity in the polynomial coefficients}} is achieved (w.r.t. the classical hyperinterpolation).

\end{frame} 



\begin{frame}[fragile]
\frametitle{Lasso hyperinterpolation: properties II}

\smnsp
\begin{theorem}
Let $f \in C(\Omega)$. Then 
\begin{equation}\label{L1HyperLeastApprox}
\min\limits_{p \in \mathbb{P}_n(\Omega)}  \braces{ \frac{1}{2}\summ{j}{M_{2n}} w_j( p ({x}_j) - f({x}_j))^2  + \lambda\summ{\ell}{N_n}\mu_{\ell}|\gamma_{\ell}|}
\end{equation}
with $$p({x})=\sum_{k=1}^{N_n} \gamma_{k} p_{k}({x}) \in \mathbb{P}^d_n(\Omega).$$
\end{theorem}

\begin{remark}
Notice that
\begin{itemize}
\item the first term in (\ref{L1HyperLeastApprox}) is exactly that of classical hyperinterpolation, 
\item the second term in (\ref{L1HyperLeastApprox})  is a penalization,
\item issue: how to choose suitable $\lambda$, $\mu_1, \ldots, \mu_{N_n}$.
\end{itemize}
\end{remark}

\end{frame} 



\begin{frame}[fragile]
\frametitle{Hybrid hyperinterpolation}

Now we introduce hybrid interpolation that mixes the features of filtered and Lasso hyperinterpolation.

\begin{definition}[Hybrid hyperinterpolation]
Suppose that
\begin{itemize}
\item $\Omega \subset {\mathbb{R}}^d$ is a compact domain and $f \in \mathcal{C}(\Omega)$,
\item $\langle f,g \rangle_{M_{2n}}$ is determined by an $M_{2n}$-point quadrature rule of PI-type in $\Omega$ with algebraic degree of exactness $2n$,
\item $h$ is a filter function,
\item $\{\mathcal{S}_{\lambda \mu_{\ell}}\}_{\ell=1}^{N_n}$ are soft thresholding operators with $\lambda>0$ and $\mu_1,\ldots,\mu_{N_n} >0$
\end{itemize}
The \emph{hybrid hyperinterpolation} of $f$ onto $\mathbb{P}^d_{n}(\Omega)$ is defined as
\begin{equation}\label{equ:hybridhyper}
{\cal{H}}^{\lambda}_n{f}:=\sum_{k=1}^{N_n} {\purples{h\left(\frac{\deg p_{k}}{n} \right)\mathcal{S}_{\lambda \mu_{k}} (\langle{f, p_{k}}\rangle_{M_{2n}})}} p_{k}.
\end{equation}
\end{definition}

\end{frame} 


\begin{frame}[fragile]
\frametitle{Hybrid hyperinterpolation: properties I}
 
Consider the $\ell_2^2 + \ell_1$-regularized LS problem
\begin{equation}\label{L22L1HyperLeastApprox}
\min\limits_{p \in \mathbb{P}^d_n(\Omega)}  \frac{1}{2}\summ{j}{N_n} w_j( p({x}_j) - f({x}_j))^2   + \frac{1}{2}\summ{j}{N} w_j({\mathbb{R}}_n{p({x}_j)})^2 + \lambda\summ{\ell}{d}\mu_{\ell}|\beta_{\ell}|,
\end{equation}
where
\begin{itemize}
\item $\lambda >0$, $\mu_1,\ldots,\mu_d >0$,
\item the operator ${\mathbb{R}}_n$ is defined as
\begin{equation}\label{R1}
{\mathbb{R}}_n p = \summ{\ell}{N_n} b_{\ell} \langle p_{\ell}, p \rangle_{M_n} p_{\ell} = \summ{\ell}{N_n} b_{\ell} \beta_{\ell} p_{\ell}
\end{equation}
with
\begin{equation}\label{coef:b}
b_{\ell} = \left\{\begin{array}{cc}
0, & \frac{\deg p_{\ell}}{n} \in [0,\frac{1}{2}] ,\\
\sqrt{\frac{1}{h_{\ell}} -1}, & \frac{\deg p_{\ell}}{n} \in [\frac{1}{2},1).
\end{array}
\right.
\end{equation}
\end{itemize}

\end{frame} 

\begin{frame}[fragile]
\frametitle{Hybrid hyperinterpolation: properties II}

\begin{theorem}\label{thm:HybridHyperSol}
Let $f \in C(\Omega)$. Then $\hybridhyper{f}$ is the solution to the regularized least squares approximation problem 
\begin{equation}\label{L22L1HyperLeastApproxA}
\min\limits_{p \in \mathbb{P}^d_n(\Omega)}  \frac{1}{2}\summ{j}{N_n} w_j( p({x}_j) - f({x}_j))^2   + \frac{1}{2}\summ{j}{N} w_j({\mathbb{R}}_n{p({x}_j)})^2 + \lambda\summ{\ell}{d}\mu_{\ell}|\beta_{\ell}|,
\end{equation}
with $$p(x)=\sum_{k=1}^{N_n} \gamma_{k} p_{k}(x) \in \mathbb{P}^d_n(\Omega).$$
\end{theorem}


 

\end{frame} 




\begin{frame}[fragile]
\frametitle{Technical difficulties of Hyperinterpolation}

There are technical difficulties.

\smsp

\begin{itemize}
\item {\purples{Low cardinality rules}}: for a fixed degree of precision $n$, minimal rules of PI-type exist in rare cases; alternatively one wants to determine those with the lowest number of nodes $M_n$ or not exceeding the dimension $N_n$ of the polynomial space ${\mathbb{P}}_n$.
\item {\purples{Computation of orthonormal basis}}: in some cases orthonormal basis w.r.t. a suitable weight function are available (e.g. interval, n-cube, sphere), but most of the times they are not available and one must determine them numerically.
\item {\purples{Computation of $\lambda$}}: in Lasso and Hybrid interpolation the performance depends on a parameter $\lambda$ and it is not trivial to determine a good one (argument of future research).
\end{itemize}



\end{frame}




\begin{frame}[fragile]
\frametitle{Numerical experiments: union of disks}

\smsp

The domain $\Omega$ is the {\purples{union of $M$ disks}} $B(C_k,r_k)$ with centers $C_k$ and radii $r_k$, i.e. 
$$\Omega=\cup_{k=1}^{M} B(C_k,r_k) \subset {\mathbb{R}}^2.$$

\begin{itemize}
\item low cardinality rules are available in which $N_n \leq M_n$;
\item triangular polynomial orthonormal basis are computed numerically by a QR based algorithm.
\end{itemize}

In the tests we set $\Omega=\Omega^{(r_1)} \cup \Omega^{(r_2)}$, where $\Omega_2^{(r_j)}$, $j=1,2$, is the union of $19$ disks, with
\begin{itemize}
\item centers $P_k^{(r_j)}=(r_j \cos(\theta_k),r_j \sin(\theta_k))$, with $r_1=2$ and $r_2=4$.
\item $\theta_k=2k\pi/19$, $k=0,\ldots,18$,
\item radii equal to $r_j/4$, $j=1,2$.
\end{itemize}



\end{frame}



\begin{frame}[fragile]
\frametitle{Some numerical experiments on union of disks}

\begin{figure}[!htbp]
  \centering
   %{\includegraphics[scale=0.15,clip]{expo_1A.eps}}
   {\includegraphics[scale=0.15,clip]{expo_1A-eps-converted-to.pdf}}
 \caption{The domain $\Omega=\Omega^{(r_1)} \cup \Omega^{(r_2)}$ in which we perform our tests. We represent in red the $N=496$ nodes of the cubature rule for algebraic degree of exactness=$30$, useful for classical hyperinterpolation at polynomial degree 15.}
 \label{fig_e1}
\end{figure}.

\end{frame}





\begin{frame}[fragile]
\frametitle{Some numerical experiments: union of disks}

We consider contaminated evaluations of
\begin{eqnarray*}
    f({\mathbf{x}})=(1-(x^2+y^2)) \exp ( x \, \cos(y)),
\end{eqnarray*}
contaminated by Gaussian noise ($\sigma=0.05$).

\smsp

Numerical tests show the favourable ratio between the sparsity and $L_2$ errors of the hybrid hyperinterpolants.

\begin{table}
\begin{center}
  \begin{tabular}{|c|c| c  |c|c| c  |c|}
    \hline
& $\lambda$ & {$L_2$ error} & {sparsity} & {$\lambda$ } & {$L_{\infty}$ error}& {sparsity}  \\
    \hline

{\em{Hyperint.}} & $ - $ & $ 0.045 $ & $ 136 $ & $ - $ & $0.1185$ & $136$ \\

{\em{Filtered}} & $ - $ & $ 0.0313 $ & $ 120 $ & $ - $ & $0.1043$ & $120$ \\

{\em{Lasso}} & $ 0.0059 $ & $ 0.0195 $ & $ 30 $ & $ 0.0072 $ & $0.0449$ & $16$ \\

{\em{Hybrid}} & $ 0.0048 $ & $ 0.0169 $ & $ 37 $ & $ 0.0068 $ & $0.0392$ & $21$ \\

    \hline
  \end{tabular}
  \caption{Tests at degree $15$ with $\lambda$ choosen to minimize $L_2$ and $L_{\infty}$ errors.}
  \label{union_disk_table}
  \end{center}
\end{table}

\end{frame}



\begin{frame}[fragile]
\frametitle{Some numerical experiments: union of disks}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=\textwidth]{udisks1.eps}
    \caption{Approximate  $f(x,y)=(1-(x^2+y^2)) \exp ( x \, \cos(y))$, perturbed by Gaussian noise ($\sigma=0.05$), over the union of disks, via
    hyperinterpolation $\mathcal{L}_L f^{\epsilon}$,
    filtered hyperinterpolation $\mathcal{F}_{L,N}f^{\epsilon}$, and
    hybrid hyperinterpolation $\mathcal{H}_{L}^{\lambda} f^{\epsilon}$ with $L=15$ (different notation w.r.t. the talk).}
    \label{figure:unionofdisks}
\end{figure}

\end{frame}




\begin{frame}[fragile]
\frametitle{Some numerical experiments: union of disks}
\begin{figure}[htbp]
    \centering
    \includegraphics[width=\textwidth]{udisks2.eps}
    \caption{The choices of regularization parameter $\lambda$   for hybrid hyperinterpolation $\mathcal{H}_n^{\lambda} f^{\epsilon}$ with $n=15$ for approximating
    $f(x,y)=(1-(x^2+y^2)) \exp ( x \, \cos(y))$, perturbed by Gaussian noise ($\sigma=0.05$)).}
    \label{figure:unionofdisks2}
\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, {\blues{Interpolation and hyperinterpolation over general regions}}, J. Approx. Theory 83 (1995), 238--254.

\item {\purples{Filtered hyperinterpolation:}} I.H. Sloan, R.S. Womersley  {\blues{Filtered hyperinterpolation: A constructive polynomial approximation on the sphere}}, GEM - International Journal on Geomathematics, 3 (2012), 95--117.

\item {\purples{Lasso hyperinterpolation:}} C. An, H.-N. Wu, {\blues{Lasso hyperinterpolation over general regions}}, SIAM J. Sci. Comput., 6 (2021), A3967--A3991.

\item {\purples{Hybrid hyperinterpolation:}} C. An, J. Ran, A. Sommariva, {\blues{Hybrid hyperinterpolation over general regions}}, submitted,

\item {\purples{Details on hyperinterpolation implementation:}} A. Sommariva and M. Vianello, {\blues{Numerical hyperinterpolation over nonstandard planar regions}}, Math. Comput. Simulation, 141 (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}




\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}



