\documentclass[aspectratio=149, sansserif, 10pt]{beamer}

\input{hal_prefs.tex}

\title{On some theorems by Tchakaloff, Davis and Wilhelmsen and their applications}

\author{{\purples{\bf{Alvise Sommariva}}}}
\date{{\blues{Probabilistic methods, Signatures, Cubature and Geometry.}} \\ University of York, January 9-11, 2023 \\ {\vspace{0.4cm}} {\reds{Thanks}}: Conference Organizers, co-workers, RITA, UMI group TAA, GNCS project {\sl{Methods and software for multivariate integral models}}}


\begin{document}

\thispagestyle{empty}
\frame{\titlepage}


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


In this talk, 
\begin{itemize}
\item we start introducing the well-known {\reds{Tchakaloff theorem}} (and one of its variants) \\ {\scriptsize{existence theorem of certain algebraic and {\it{low}} cardinality cubature rules with positive weights on a multivariate compact domain $\Omega \subset {\mathbb{R}}^d$,;}} 
\item  we compress cubature rules, \\ {\scriptsize{we show how, from cubature rules on $\Omega$ with positive weights and interior nodes (i.e. of PI-type), whose algebraic degree of precision $ADE$ is equal to $m$ and the number of nodes higher than the dimension \dlacc $r$" of the polynomial space ${\mathbb{P}}_m(\Omega)$ of total degree $m$, we can {\reds{extract rules of PI-type}} but with at most \dlacc$r$" nodes, by means of Lawson-Hanson algorithm}};
\item we recall {\reds{some results by Davis and Wilhelmsen}}  \\ {\scriptsize{show how they imply that we can determine rules of PI-type, with ADE equal to $m$, if the moments w.r.t. to a basis of ${\mathbb{P}}_m(\Omega)$ are available as well as if we can evaluate numerically the characteristic function ${\cal{X}}_{\Omega}$}};
\item {\reds{examples of this novel meshless approach}}  \\ {\scriptsize{time permitting}}
\end{itemize}

\begin{vuoto}
{\blues{Important}}: all the Matlab routines used in this talk are available at the author's homepage.
\end{vuoto}
\end{frame}



\begin{frame}[fragile]
\frametitle{Example: quadrature over polygons (meshfree approach)}

\begin{paper}
 {\teals{\href{https://link.springer.com/article/10.1007/s10543-007-0131-2}{Product Gauss cubature over polygons based on Green's integration formula}}} (2007)
\end{paper}

\begin{itemize}
\item {\purples{purpose}}:  {{cubature formula over convex, nonconvex or even multiply connected polygons}} $\Omega$. 
\item  {\purples{ADE}}:  the formula with nodes $\{ {\mathbf{x}}_k \}_k \subset {\mathbb{R}}^2$ and weights $\{ w_k \}_k \subset {\mathbb{R}}$ has {\reds{algebraic degree of exactness ADE=$2n-1$}}, i.e. denoting by ${\mathbb{P}}_{m}$ the space of polynomials of total degree $m$, as well as ${\mathbf{x}}_k=(x_k,y_k)$,
$$
\int_{\Omega} p(x,y) dx dy=\sum_{k=1}^N w_k p(x_k,y_k), {\mbox{ for all }} p \in {\mathbb{P}}_{2n-1},
$$
with $N \approx m \, n^2$,  where ``{\em{m}}"  is the number of sides that are not orthogonal to a given line, and not lying on it. 
\item  {\purples{preprocessing}}: {\reds{it does not need any preprocessing like triangulation of the domain}}, but relies directly on univariate Gauss-Legendre quadrature via Green’s integral formula. 
\end{itemize}
\end{frame}



\begin{frame}[fragile]
\frametitle{Example: quadrature over polygons (meshfree approach)}

\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.2425,clip]{fig_01.eps}}
   {\includegraphics[scale=0.2425,clip]{fig_03.eps}}
 \caption{Examples of polygonal domains (ADE=$9$). 
 
 {\vspace{0.1cm}} 
 
 {\purples{Left:}} a convex domain with 6 sides (180 nodes), {\purples{Right:}} a non-convex domain with 9 sides (255 nodes).
{\vspace{0.2cm}}


 {\vspace{0.1cm}} 


{\reds{Red dots}}: nodes with negative weights. {\blues{Blue dots}}: nodes with positive weights.}
 \label{fig_2CRI}
 \end{figure}

\end{frame}




\begin{frame}[fragile]
\frametitle{Example: quadrature over polygons (meshfree approach)}

\begin{itemize}
\item {\blues{\bf Pros}}: {\purples{meshfree approach}}.
\item {\reds{\bf Cons}}: In general these rules may have {\purples{external nodes}} as well as {\purples{negative weights}}.
\item {\blues{\bf Pros}}: With some trick based on roto-translations we can compute rules over {\purples{convex domains}} of PI-type (positive weights and internal nodes).
\item {\reds{\bf Cons}}: With some trick based on roto-translations we can compute rules over some {\purples{non-convex domains}} but in general this is not always possible.
\item {\reds{\bf Cons}}: Rule may have {\purples{high cardinality}} if the polygon has many sides.
\end{itemize}


{\vspace{0.3cm}}

We prefer rules with

\begin{itemize}
\item {\purples{internal nodes}} (integrand may not be defined outside the domain),
\item {\purples{positive weights}} (more numerical stability and application to hyperinterpolation),
\end{itemize}
i.e. of PI-type. 


{\vspace{0.3cm}}
Furthermore, we look for rules with {\purples{low cardinality}} (few nodes, hence few samples of the function).

\end{frame}




\begin{frame}[fragile]
\frametitle{Example: quadrature over polygons (triang. based approach)}

\begin{paper}
 {\teals{\href{https://www.sciencedirect.com/science/article/abs/pii/S0377042719306636}{Compressed cubature over polygons with applications to optical design}}} (2020)
\end{paper}

\begin{itemize}
\item {\purples{purpose}}:  {{cubature formula over convex, nonconvex or even multiply connected polygons}} $\Omega$. 
\item  {\purples{strategy}}: once a minimal triangulation is available (see Matlab {\tt{polyshape}} toolbox), we obtain the rule by applying an almost-minimal rule of PI-type on each triangle with the wanted ADE, summing the contributions.
\end{itemize}

Some remarks
\begin{itemize}
\item {\purples{minimal triangulation}}: one can triangulate a general $M$-sides polygon via $M-2$ triangles (easy task in a convex polygon, not trivial for a general one),
\item {\purples{almost-minimal rule}} the number of its nodes is almost minimal between those having a certain degree of precision and requirement on the nodes (e.g. internal) and weights (e.g. positive).
\end{itemize}

\end{frame}



\begin{frame}[fragile]
\frametitle{Example: quadrature over polygons (triang. based approach)}

 {\vspace{-0.1cm}} 

\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.29,clip]{fig_05.eps}}
   {\includegraphics[scale=0.29,clip]{fig_06.eps}}
 \caption{Examples of polygonal domains (ADE=$9$). 
 
 {\vspace{0.1cm}} 
 
 {\purples{Left:}} a convex domain with 6 sides (77 nodes, the previous rule had 180 nodes), 
 {\purples{Right:}} a non-convex domain with 9 sides (133 nodes, the previous rule had 255 nodes).

 {\vspace{0.1cm}} 

{\reds{Note}}: all the weights are positive.}
 \label{fig_2CRI}
 \end{figure}

\end{frame}


\begin{frame}[fragile]
\frametitle{Example: quadrature over polygons (triang. based approach)}

\begin{itemize}
\item {\reds{\bf Cons}}: {\purples{requires triangulation}}.
\item {\blues{\bf Pros}}: In general these rules always have {\purples{internal nodes}} as well as {\purples{positive weights}}.
\item {\reds{\bf Cons}}: Rule still may have {\purples{high cardinality}} if the polygon has many sides.
\end{itemize}

{\vspace{0.3cm}}

Observe that in the examples above for ADE$=9$

\begin{itemize}
\item {\purples{convex domain}}: the rule has 77 nodes,
\item {\purples{not convex domain}}:  the rule has 133 nodes.
\end{itemize}

{\vspace{0.3cm}}

\begin{remark}
In both cases the number of nodes is higher than the {\purples{dimension of the polynomial space}} ${\mathbb{P}}_9$ that is equal to $(9+1)(9+2)/2=55$.

{\vspace{0.2cm}}

Our project is to quickly {\purples{extract}} from the previous one, another rule of PI-type with the same degree of precision but with a number of nodes at most equal to ${\mbox{dim}}({\mathbb{P}}_9)$ (i.e. a {\purples{rule compression}}).

\end{remark}

\end{frame}



\begin{frame}[fragile]
\frametitle{Carath\'eodory-Tchakaloff Subsampling}

Tchakaloff theorem, a cornerstone of quadrature theory, substantially asserts that:

\begin{vuoto}
{\purples{For every compactly supported measure there exists a positive algebraic quadrature formula of ADE=$m$ with cardinality not exceeding the dimension of ${\mathbb{P}}_m$ (restricted to the measure support)}}
\end{vuoto}

i.e. the goal of the previous section can be achieved.

{\vspace{0.2cm}}

\begin{itemize}
\item Originally proved by V. Tchakaloff in 1957 for absolutely continuous measures, it has then be {\purples{extended to any measure with finite polynomial moments}}, and to arbitrary finite dimensional spaces of integrable functions.

\item We begin by stating a {\purples{discrete version of Tchakaloff theorem}} whose proof is based on a theorem by Carath\'eodory. See also 

\end{itemize}

\begin{paper}
{{
M. Putinar, {\teals{\href{https://www.jstor.org/stable/2162136}{A Note on Tchakaloff's Theorem}}}, Proceedings of the American Mathematical Society
Vol. 125, No. 8 (Aug., 1997), pp. 2409--2414.
}}
\end{paper}

\end{frame}



\begin{frame}[fragile]
\frametitle{Carath\'eodory-Tchakaloff Subsampling}


\begin{theorem}[Carath\'eodory-Tchakaloff, see more general Putinar theorem]
Let
\begin{enumerate}
\item $\mu$ be a {\purples{multivariate discrete measure}} supported at a finite set $X=\{ {\mathbf{x}}_k \}_{k=1,\dots,N} \subset {\mathbb{R}}^d$, with correspondent positive weights $\{w_k\}_{k=1,\dots,N}$, 
\item $\Phi = {\mbox{span}}(\phi_1,\ldots,\phi_r)$ a finite dimensional space of {\purples{d-variate functions}} defined on $\Omega \supseteq X$, with ${\mbox{dim}}(\Phi |_X ) \leq r$.
\end{enumerate}
Then there exist a quadrature formula with nodes $T = \{{\mathbf{t}}_k\}_{k=1,\dots,N_c} \subseteq X$ and positive weights $\{{u}_k\}_{k=1,\dots,N_c}$, such that {\reds{$N_c \leq {\mbox{dim}}(\Phi |_X )$}} and
$$
\int_{\Omega} f({\bf{x}}) d\mu:=\sum_{k=1}^N w_k f({\bf{x}}_k)=\sum_{i=1}^{N_c} {u}_i f({\bf{t}}_i), {\mbox{ for all }} f \in \Phi |_X.
$$
\end{theorem}

Roughly speaking, adapting the theorem to our needs,

\begin{vuoto}
If we have a rule of PI-type with ADE=$m$ and cardinality $N$ higher than $r={\mbox{dim}}({\mathbb{P}}_m(\Omega))$ then we can extract one of PI-type with ADE=$m$ and cardinality {\reds{$N_c \leq r$}}.
\end{vuoto}
\end{frame}



\begin{frame}[fragile]
\frametitle{Carath\'eodory-Tchakaloff Subsampling}

Given 
\begin{itemize}
\item a formula of PI-type with ADE=$m$, nodes $X=\{ {\mathbf{x}}_k \}_{k=1,\dots,N} \subset {\mathbb{R}}^d$ and positive weights $\{{{w}}_k\}_{k=1,\dots,N}$,
\item a basis $\{\phi_1,\ldots,\phi_r\}$ of ${\mathbb{P}}_m(\Omega)$,
\end{itemize}

\noindent let 

\begin{itemize}
\item $V_{i,j}=(\phi_j ({\mathbf{x}}_i))$ the Vandermonde matrix at the nodes,
\item ${\mathbf{b}}=(b_j)_{j=1,\ldots,r}$ where $b_j=\int_{\Omega} \phi_j d\mu=\sum_{i=1}^N w_i \phi_j({\mathbf{x}}_i)$, the vector of the $\mu$ moments.
\end{itemize}
The problem mentioned above resorts into {\purples{computing a nonnegative solution with at most \dlacc$r$" nonvanishing}} components to the underdetermined linear system
$$
V^T {\mathbf{u}}={\mathbf{b}}.
$$
\end{frame}




\begin{frame}[fragile]
\frametitle{Carath\'eodory-Tchakaloff Subsampling}

The computation of a nonnegative solution with at most $r={\mbox{dim}}({\mathbb{P}}_{m}(\Omega))$ nonvanishing components to the underdetermined linear system
$V^T {\mathbf{u}}={\mathbf{b}}$ can be perfomed finding a sparse solution to the quadratic minimum problem
$$
{\mbox{NNLS:}} 
\left\{ 
\begin{array}{l}
\min_u \|V^T {{\mathbf{u}}}-{\mathbf{b}}\|_2 \\
{\mathbf{u}} \geq 0
\end{array}
\right.
$$
via {\purples{Lawson-Hanson active set method}} for NonNegative Least Squares (NNLS).

{\vspace{0.2cm}}

In Matlab this can be done by means of the Matlab built-in routine {\purples{\tt{lsqnonneg}}} as well as by the more recent {\purples{\tt{LHDM}}} by Dessole, Marcuzzi and Vianello.

\begin{remark}
\begin{itemize}
\item The approach mentioned above is effective for {\reds{mild ADE}}, say on the order of ADE=20 for bivariate domains and ADE=10 for trivariate domains.
\item There are also {\reds{other approaches}}, e.g. by based on linear programming or by a different combinatorial algorithm (recursive Halving Forest), based on SVD.
\end{itemize}
\end{remark}
\end{frame}



\begin{frame}[fragile]
\frametitle{Carath\'eodory-Tchakaloff Subsampling}
As example, we can consider the application of the technique mentioned above to extract a rule of PI-type, for computing a similar one on the polygonal domains treated above.

{\vspace{0.1cm}}

\begin{algorithm}
{\tt{input}}: the nodes $\{ {\mathbf{x}}_k \}_{k=1,\ldots,N}$, the weights $\{ w_k \}_{k=1,\ldots,N}$ of a PI-type rule with $N > r$ ($r$ is the dimension of the polynomial space ${\mathbb{P}}_m$ when ADE=$m$) and a polynomial basis $\{\psi_j\}_{j=1,\ldots,r}$;
\begin{enumerate}
\item {\purples{Vandermonde matrix}}: compute $U=(\phi_k({\mathbf{x}}_i))$;
\item {\purples{fight ill-conditioning}}: compute the QR factorization with column pivoting ${\sqrt{W}} U(:,\pi)=QR$, where ${\sqrt{W}}=diag(\{w_k\})$ and $\pi$ is a permutation vector; this corresponds to a {\purples{change of basis}} $(\phi_1,\ldots,\phi_r)=(\psi_1,\ldots,\psi_r)R^{-1}$, so obtaining an orthonormal basis w.r.t. the discrete measure defined by the nodes $\{ {\mathbf{x}}_k \}_{k=1,\ldots,N}$ and the weights $\{ w_k \}_{k=1,\ldots,N}$;
\item {\purples{moments}}: evaluate the vector ${\mathbf{b}}=Q^T {\mathbf{w}}$ where ${\mathbf{w}}_k=w_k$; 
\item {\purples{compute a positive sparse solution}}: solve $Q^T {\mathbf{u}}={\mathbf{b}}$ by Lawson-Hanson algorithm (or its alternatives).
\end{enumerate}
\end{algorithm}

\end{frame}

\begin{frame}[fragile]
\frametitle{Carath\'eodory-Tchakaloff Subsampling}

By this algorithm, 
\begin{itemize}
\item adopting as {\purples{basis}} $\{\psi_j\}$ the total-degree product Chebyshev basis of the smallest Cartesian rectangle $[a_1,b_1] \times [a_2,b_2]$ containing $\Omega$, with the graded lexicographical ordering, 
\item {\purples{from the PI-rules}} with ADE=9 obtained via triangulation,
\end{itemize}
 we get the PI-rules below with cardinality $55={\mbox{dim}}({\mathbb{P}}_9)$.
\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.18,clip]{fig_05C.eps}}
   {\includegraphics[scale=0.18,clip]{fig_06C.eps}}
 \caption{Examples of polygonal domains (ADE=$9$). 
 
 {\vspace{0.05cm}} 
 
 {\purples{Left:}} a convex domain with 6 sides (55 nodes, the previous rule had 77 nodes), 
 {\purples{Right:}} a non-convex domain with 9 sides (55 nodes, the previous rule had 133 nodes).
 }
 \label{fig_2CRIcom}
 \end{figure}

\end{frame}


\begin{frame}[fragile]
\frametitle{Carath\'eodory-Tchakaloff Subsampling}

{\vspace{-0.15cm}}

\begin{remark}[When do not apply this technique]
We observe that this approach is useful only when the initial rule of PI-type with ADE=m has cardinality higher the dimension $L$ of ${\mathbb{P}}_m(\Omega)$.

{\smsp}

Thus it is {\purples{worthless}} in the case of classical domains as the interval, the disk, simplex, cube, sphere, where there are explicit rules of PI-type with cardinality inferior to $L$.
\end{remark}

{\vspace{-0.1cm}}

\begin{remark}[Cputimes on the previous domains]
For mild ADE the {\purples{computation of these compressed rules is fast}}. Running Matlab R2022a, on a computer with an Apple M1 processor and 16GB of RAM, we had average cputimes as in the table below:

 \begin{table}
\begin{center}
  \begin{tabular}{| c | c | c|}
    \hline
    {\em{Domain}} & {\em{tri. rule}} &  {\em{compress.}} \\
    \hline
    {\em{Convex domain}} & 1.1-3s & 1.7-2s \\
    {\em{Non-convex domain}}  & 5.2-3s & 1.1-2s \\
    \hline
  \end{tabular}
  \caption{Average cputime for computing rules with ADE=9 in the polygonal domains used in the tests.}
  \label{cube_table}
  \end{center}
\end{table}

{\vspace{-0.33cm}}
\end{remark}

\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem}
In the paper

\begin{vuoto}
D. R. Wilhelmsen, {\teals{\href{https://www.jstor.org/stable/2005429}{A Nearest Point Algorithm for Convex Polyhedral 
Cones and Applications to Positive Linear approximation}}}, Math. Comp., (30) 1976, pp. 48--57,
\end{vuoto}

\smsp
the author extended a result from

\begin{vuoto}
P. J. Davis, {\teals{\href{https://www.ams.org/journals/mcom/1967-21-100/S0025-5718-1967-0222534-4/S0025-5718-1967-0222534-4.pdf}{A construction of nonnegative approximate quadratures}}}, Math. Comp., (21) 1967, pp. 578--582.
\end{vuoto}
\smsp
In his own words

\begin{vuoto}
A constructive proof of the Tchakaloff theorem was given by Davis. Although his paper deals only with the integration functional, his results are easily adapted to more general functionals.
\end{vuoto}


In view of both the contributions we will refer to the next result as {\reds{Davis-Wilhelmsen theorem}} and show its impact to the numerical construction of algebraic cubature rules.
\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem}

\smnsp

\begin{definition}[Tchakaloff set]
\begin{enumerate}
\item Consider a finite-dimensional function space 

\smnsp
$${\Phi}=span\{\phi_1,\dots,\phi_r\}, {\hspace{0.3cm}}r={\mbox{dim}}(\Phi),$$ on a compact domain 
$\Omega \subset{R^d}$, 
satisfying the {\reds{Krein condition}}, i.e. there exists a function in
$\Phi$ not vanishing on $\Omega$ (e.g., a polynomial space).
\item Let ${\mathcal{L}}$ be a {\reds{strictly positive linear functional}} on ${\Phi}$, i.e. ${\mathcal{L}}(\phi)>0$ 
for every $\phi\in \Phi$, $\phi\geq 0$ not identically vanishing on $\Omega$.  
\end{enumerate}
A Tchakaloff set for ${\mathcal{L}}$ on $\Phi$ is a subset, say ${T}\subset 
\Omega$, that 
contains the support 
of a Tchakaloff-like representation for ${\mathcal{L}}$, 
i.e. 

\smnsp
{\reds{
$${\mathcal{L}}(\phi)=\sum_{j=1}^{N_c}
{w_j\,\phi(t_j)}, \smhsp \forall \phi\in \Phi,$$ 
}}
where {\reds{$N_c\leq r={\mbox{dim}}(\Phi)$, $\{t_j\}\subset 
{T}$}} and {\reds{$w_j>0, \,\, j=1,\ldots,\nu$}}.   

\end{definition} 

\smnsp

\begin{example}[Wilhelmsen, p.49]
Set ${\Phi}={\mathbb{P}}_L$ and ${\mathcal{L}}(\phi)=\int_{\Omega} \phi({\mathbf{x}}) w({\mathbf{x}}) \, d{\mathbf{x}}$, $w \geq 0$ weight function.
\end{example}
\end{frame}






\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem}
\begin{theorem}  
Let
\begin{enumerate}
\item $\Phi$ be the linear {\reds{span of continuous, 
real-valued, linearly independent functions}} $\{\phi_k\}_{k=1,\ldots,r}$ 
defined on a compact set $\Omega \subset {\mathbb{R}}^d$. 
\item Assume that $\Phi$ satisfies the {\reds{Krein condition}} (i.e. there is at least 
one $f \in \Phi$ which does not vanish on $\Omega$) and that ${\mathcal{L}}$ is a {\reds{strictly positive 
linear functional}} on ${\Phi}$, i.e. ${\mathcal{L}}(\phi)>0$ 
for every $\phi\in \Phi$, $\phi\geq 0$ not identically vanishing on $\Omega$.  
\item $\{P_i\}_{i=1}^{+\infty}$ is an everywhere dense subset of $\Omega$
\end{enumerate}
Then {\purples{for sufficiently large $\mathcal{I}$, the set $X=\{P_i\}_{i=1}^
{\mathcal{I}}$ is 
a Tchakaloff set}}, i.e.
\begin{equation} \label{tchcub}
{\mathcal{L}}(f)=\sum_{j=1}^{N_c} w_j f(Q_j), {\mbox{   }} \forall f \in \Phi
\end{equation}
where $w_j > 0, \, j=1,\ldots,N_c$ and $\{Q_j\}_{j=1}^{N_c} 
\subset X \subset \Omega$, 
with $N_c={\mbox{card}}(\{Q_j\}) \leq r$.
\end{theorem}
\end{frame}







\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem}
In our numerical framework this theorem is important because it says that:

\begin{tcolorbox}
{\purples{a sufficiently dense set of the compact domain $\Omega$ contains the nodes of a PI-type rule with algebraic degree of precision ADE=$m$}}
\end{tcolorbox}

\begin{algorithm}[sketch]
\begin{itemize}
\item {\purples{Moment computation}}: determine $\{\gamma_k\}_{k=1,\ldots,r}$ of a certain polynomial basis $\{\phi_k\}_{k=1,\ldots,r}$ of ${\mathbb{P}_m(\Omega)}$; 
\item {\purples{Pointset}}: using an {\tt{in-domain}} routine on a mesh in a domain ${\cal{R}}$ containing $\Omega$, determine {\reds{a sufficient number of points}} $\{ {\tilde{P}}_l\}_{l=1,\ldots,N}$ inside $\Omega$ so that the overdetermined linear system $V^T w=\gamma$, with $V_{l,k}=(\phi_k({\tilde{P}}_l))$, has a nonnegative solution $w$ with at most $N_c \leq r$ positive components.
\item {\purples{Rule extraction}}: solve $V^T w=\gamma$ via fast Lawson-Hanson algorithm (or alternatives).
\end{itemize}
\end{algorithm}
In spite of the simplicity of this approach there are many aspects that deserve explanations, on the implementation side as well as on the theoretical one.
\end{frame}



\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: polynomial basis}
The choice of the polynomial basis $\{\phi_k\}_{k=1,\ldots,r}$ of ${\mathbb{P}}_m$ is not trivial. It must be
\begin{itemize}
\item not too badly conditioned,
\item allow a fast computation of the moments $$\gamma_k=\int_{\Omega} \phi_k({\mathbf{x}}) \, d\mu, \smhsp k=1,\ldots,r.$$
\end{itemize}

\begin{remark}[Choice of the polynomial basis]
{\small{
For many multivariate regions, it is enough to use a specific basis based on {\purples{total-degree product Chebyshev basis}} of the smallest Cartesian hyper-rectangle $\prod_{k=1}^d [a_k,b_k]$ containing $\Omega$, with the graded lexicographical ordering.

\smsp

However the choice may be much more difficult for regions over manifolds. 

\smsp

{\blues{Example.}} If $\Omega$ is {\purples{compact domain on the unit-sphere ${\mathbb{S}}_2 \subset {\mathbb{R}}^3$}}, e.g. a spherical polygon, 
\begin{itemize}
\item the classical {\reds{spherical harmonics}} restricted on $\Omega$ may be severely ill-conditioned, 
\item alternatively, using the {\reds{tensorial basis of ${\mathbb{R}}^3$}}, restricted to the sphere, some elements of the trivariate basis mentioned above must be somehow discarded, since the dimension of the polynomials ${\mathbb{P}}_m({\mathbb{S}}_2)$ is $(m+1)^2$ while ${\mathbb{P}}_m({\mathbb{R}}^3)$ is $(m+1)(m+2)(m+3)/6$ (i.e. for $m=9$ we have ${\mbox{dim}} ({\mathbb{P}}_9({\mathbb{S}}_2))=100$ while ${\mbox{dim}} ({\mathbb{P}}_9({\mathbb{R}}^3))=220$).
\end{itemize}
}}
\end{remark}
\end{frame}





\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: in-domain routines}
A second issue is the availability of the {\purples{in-domain routine}}.

\begin{itemize}
\item {\purples{Verification of certain inequalities}}: for example, the unit-ball $B({\mathbf{0}},1)$ is defined as the set
$$
B({\mathbf{0}},1):=\{ {\mathbf{x}}=(x,y,z) \in {\mathbb{R}}^3 {\mbox{ such that }} x^2+y^2+z^2 \leq 1 \}.
$$
\item {\purples{Specific codes}}: Matlab built-in {\tt{inpolygon}} (polygonal domains), {\tt{inpolyhedra}} (polyhedral domains).
\end{itemize}

At the same time, it is important the choice of a {\purples{sufficient dense set}}. 

\begin{itemize}
\item Too low cardinality: many iterations of the algorithm.
\item Too high cardinality: expensive cost of the indomain routine.
\end{itemize}

\smsp

Thus, to determine quickly a sufficiently dense set $T$, 

\begin{itemize}
\item its {\purples{cardinality and distribution must be well-thought a priori}} 
\item the {\purples{in-domain routine must be fast}}.
\end{itemize}
\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: moment computation}
A third issue is the {\purples{moment computation}}.

\smsp
\smsp

This can be done in various ways, for instance resorting to one of the following:

\begin{enumerate}
\item they are {\purples{explicitly known}};
\item a {\purples{Gauss-Green approach}} in bivariate domains, or a {\purples{divergence theorem based}} in trivariate domains;
\item {\purples{specific rules}} even with negative weights or external nodes,
\item {\purples{adaptive rules}}.
\end{enumerate}

\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: a bivariate example}

As first example we consider a bivariate domain $\Omega$ whose boundary $\partial \Omega$ is parametrically defined by {\purples{piecewise rational functions}}  i.e. $$\partial \Omega=\{{\mathbf{x}}=(x,y), x=x(t), y=y(t), t \in [a,b] \} $$ where $x,y$ are certain rational functions defined on $[a,b]$.

\smsp

\begin{example}
This class of bivariate domains $\Omega$ includes those such that $\partial \Omega$  is defined by
\begin{itemize}
\item piecewise by {\purples{NURBS}} curves,
\item {\purples{composite Bezier}} curves,
\item parametric {\purples{splines}}.
\end{itemize}
\end{example}

{\purples{Applications}}:
\begin{itemize}
\item NEFEM with NURBS-shaped curvilinear elements,
\item VEM with NURBS-shaped curvilinear elements.
\end{itemize}

\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: a bivariate example}
\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.175,clip]{ex1.eps}}
   {\includegraphics[scale=0.175,clip]{ex2.eps}}
      {\includegraphics[scale=0.175 ,clip]{ex3.eps}}
 
 \caption{Examples of integration domains.}
 \label{fig_2CRI}
 \end{figure}
 
\end{frame}



\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: a bivariate example}

{\vspace{-0.05cm}}

The pointset in $\Omega$ is extracted from a fine mesh on the minimal rectangle ${\cal{R}}^* \supseteq \Omega$.

The {\tt{in-domain}} routine can be based on {\purples{Jordan curve theorem}}: 

\begin{tcolorbox}

{\vspace{-0.1cm}}

{\reds{a point $P$ belongs to a Jordan domain $\Omega$ if and only if, 
having taken a point $P^* \notin \Omega$ then the segment 
${\overline{P^*P}}$ crosses $\partial \Omega$ an odd number $c(P)$ of times}}. 

{\vspace{-0.1cm}}

\end{tcolorbox}

\smnsp

\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.25,clip]{FIG1.eps}}
   {\includegraphics[scale=0.25,clip]{FIG2.eps}}
   
   \smnsp
      \smnsp
   
 \caption{Points and boundary intersections. On the left $c(P)=1$ and the point $P$ is in the domain. On the right $c(P)=2$ and the point $P$ is outside the domain.}
 \label{fig_2CRI}
 \end{figure} 
 
 \smnsp
 \smnsp
 
 \begin{vuoto}
 For details see:
{\teals{\href{https://www.sciencedirect.com/science/article/abs/pii/S0893965922000696}{inRS: implementing the indicator function for NURBS-shaped planar domains}}}, Applied Mathematics Letters, Volume 130, August 2022.
 \end{vuoto}


\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: a bivariate example, moments computation}
Having in mind to compute a rule with algebraic degree of precision $ADE=m$ by moments equations, we 

\begin{itemize}
\item define a {\purples{suitable basis}} $\{\phi_j\}$ of the polynomial space ${\mathbb{P}}_m$ ({\it{tensorial Chebyshev}} on the bounding box ${\cal R}^*$ of $\Omega$), 
\item compute the {\purples{moments}} $\gamma_1,\ldots,\gamma_r$, where
$$
{\reds{
\gamma_j:=\int_{\cal S} \phi_j(x,y) dx dy.
}}
$$
\end{itemize}
\smsp

To this purpose:

\begin{enumerate}
\item By applying the {\purples{Gauss-Green}} theorem, each $\gamma_j$ is the sum of some line integrals, that after some computation are shown to require the {\purples{integration in $[-1,1]$ of continuous rational functions}}.
\item We compute these integrals in $[-1,1]$ by {\purples{high-order Gauss-Legendre  rule}} (other techniques may be used).
\end{enumerate}
\end{frame}






\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: a bivariate example, test case}
 \begin{figure}[!htbp]
  \centering
    {\includegraphics[scale=0.35,clip]{domain05_06.eps}}
 \caption{{A curvilinear domain $\Omega$, the mesh points $P$ outside the domain or on its boundary
(in red), those inside the domain (in green) and the nodes of
a cubature formula of PI-type for $n=6$ (28 magenta dots). The control 
points of the NURBS curve are represented as cyan squares, joined 
to represent the so called {\it{control points polygon}}.}}
 \label{fig_1N2}
 \end{figure}

\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: a bivariate example, test case}

{{
\begin{table}[ht]
\begin{center}\small
{\renewcommand{\arraystretch}{1.1}
\begin{tabular}{| c | c | c | c | c | c | c |  }
\hline
 {\mbox{ADE}} & $\#$ & {\mbox{\#}} trial pts & {\mbox{cond}} & {\mbox{moment res}} & {\mbox{cpu}}  \\
\hline
 $ 2 $  & $6$ & $50$ ($121$) & $1$ &  5$e-$16 
&  5.3$e-$3 \\
 $4 $   & $15$ & $50$ ($121$) & $1$ 
&   7$e-$16 &  6.1$e-$3\\
 $6 $    & $28$ & $89$ ($196$) & $1$ 
&   1$e-$15 &  7.6$e-$3\\
 $8 $   & $45$ & $239$ ($484$) & $1$ 
&   2$e-$15 &  1.1$e-$2\\
 $10 $   & $66$ & $491$ ($961$) & $1$ 
&   4$e-$15 &  1.6$e-$2\\		 
\hline
\end{tabular}
        }
\caption{{Degree of precision $ADE=2,4,6,8,10$ of the rule, cardinality $\#$ 
of the extracted nodes, 
cubature conditioning and moment residual of the rule 
on curvilinear domain $\Omega$, number of trial points used 
in the extraction, cubature condition number {\tt{cond}}, moment residual 
of the rule 
and median of the cputime over 50 tests.}}
\label{tab01}
\end{center}
\end{table}
}}

\begin{vuoto}
For details see: {\teals{\href{https://www.math.unipd.it/~marcov/pdf/nurbscatch.pdf}{Low cardinality Positive Interior cubature on NURBS-shaped domains}}} (to appear on BIT Numer. Math.).
\end{vuoto}
\end{frame}



\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: other examples}

The same approach can be used to compute {\blues{meshless cubature rules over polyhedra}}:

\begin{itemize}
\item {\purples{polynomial basis}}: tensorial-type Chebyshev basis on the bounding box ${\cal{R}}$ of the domain $\Omega$;
\item {\purples{mesh points}}: {\it{sufficiently dense}} Halton points in the bounding box ${\cal{R}}$;
\item {\purples{in-domain routine}}: {\tt{inpolyhedron.m}};
\item {\purples{moment computation}}: via divergence theorem on the facets (routine requested: cubature formula over a polygon).
\end{itemize}

\end{frame}



\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: other examples}
Furthermore it can be used to compress {\blues{Quasi-Montecarlo cubature rules on $\Omega \subset {\mathbb{R}}^d$}} obtained by set operations of compact domains $\Omega_1,\ldots,\Omega_{\nu} \subset {\mathbb{R}}^d$:

\begin{itemize}
\item {\purples{polynomial basis}}: tensorial-type Chebyshev basis ${\mathbb{P}}_m$ on the bounding box ${\cal{R}}$ of the domain $\Omega$;
\item {\purples{mesh points}}:  {\it{sufficiently dense}} Halton points in the bounding box ${\cal{R}}$;
\item {\purples{in-domain routine}}: in domain routine on each domain $\Omega_k$, $k=1,\ldots,\nu$ followed by suitable set operations;
\item {\purples{moment computation}}: via Quasi-Montecarlo cubature.
\end{itemize}


This allows to achieve a rule with few nodes that equals the results of the QMC rule applied to polynomials in ${\mathbb{P}}_m$.

\smsp
\smsp

{\purples{Purpose}}: {\reds{retaining the approximation power of the original QMC formula}} (up to a quantity proportional to the best polynomial approximation error of degree $m$ to $f$, in the uniform norm on $\Omega$),  {\reds{using much fewer nodes}}.

\end{frame}

\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: other examples}

\begin{figure}[!ht]
  \centering
   {\includegraphics[scale=0.265,clip]{Nurbs_U.eps}
\includegraphics[scale=0.265,clip]{Nurbs_I.eps} 
\includegraphics[scale=0.265,clip]{Nurbs_SD.eps} 
}
\caption{231 compressed QMC nodes with exactness degree $n=20$, 
on complex shapes arising 
from union (top-left), 
intersection (top-right) and symmetric difference (bottom) of two 
NURBS-shaped domains (extraction from a million Halton points 
of domain bounding boxes, basis $\Phi$ obtained by orthonormalization of a tensorial type basis in the bounding box $\cal R$ of the domain $\Omega$).}
 \end{figure}

\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: other examples}


\begin{table}[ht]

\begin{center}

{\footnotesize

\begin{tabular}{| c | c | c | c | c |}

\hline

deg & 5 & 10 & 15 & 20\\

\hline \hline

card. $Q_c$ & 21 & 66 & 136 & 231\\

compr. ratio & 1.2e+04 & 3.9e+03 & 1.9e+03 & 1.1e+03\\

\hline \hline 

cpu $Q_c$ & 4.0e-02 & 1.2e-01 & 2.8e-01 & 5.8e+00 \\

speed-up & 6.0 & 9.0 & 11.0 & 1.4 \\

\hline \hline  

mom. resid. $Q_c^{new}$ & 5.8e-16 & 1.4e-15 & 2.4e-15  & 7.0e-15 \\

\hline

\end{tabular}
}
\caption{\small{Compression parameters of QMC cubature
with $N=255923$ Halton points on the intersection of two NURBS-shaped domains
as in Figure above top-right.}}

\label{nurbs-I}

\end{center}

\end{table}




\begin{table}[ht]
\begin{center}
{\footnotesize
\begin{tabular}{| c | c | c | c | c |}
\hline
deg & 5 & 10 & 15 & 20\\
\hline \hline
$E(f_1)$ & 2.7e-04 & 1.4e-08 & 3.0e-13 & 4.5e-16\\
\hline
$E(f_2)$ & 2.3e-04 & 2.4e-05 & 1.1e-05 & 5.6e-06\\
\hline
\end{tabular}
}
\caption{\small{Relative QMC compression errors $E(f_k)$, $k=1,2$
for the two test functions 
$f_1(P)=\exp(-|P-P_0|^2)$, $f_2(P)=|P-P_0|^5$
on the intersection of Fig. 1 top-right.}}
\label{errs-U}
\end{center}
\end{table}

\end{frame}




\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: other examples}

 \begin{figure}[!htbp]
  \centering
    {\includegraphics[scale=0.20,clip]{tetrasphere-eps-converted-to.pdf}}
 \caption{84 compressed QMC nodes with exactness degree n = 6, on intersection (red bullets) and difference (green bullets) of a tetrahedral element with a ball (extraction from a million Halton points of domain bounding boxes, cputime: $\approx 5 \cdot 10^{-2}s$, basis $\Phi$ obtained by orthonormalization of a tensorial type basis in the bounding box $\cal R$ of the domain $\Omega$).}
 \label{tetraball}
 \end{figure}

\end{frame}


\begin{frame}[fragile]
\frametitle{The Davis-Wilhelmsen theorem: other examples}

\begin{table}[ht]
\begin{center}
{\footnotesize
\begin{tabular}{| c | c | c | c |}
\hline
deg & 2 & 4 & 6\\
\hline \hline
card. $Q_c$ & 10 & 35 & 84\\
compr. ratio &2.2e+04 & 6.2e+03 & 2.6e+03\\
\hline \hline
cpu $Q_c$ & 4.1e-02 & 4.1e-02 & 1.8e-01\\
\hline
mom. resid. $Q_c$  & 1.7e-16 & 6.0e-16 & 1.2e-15\\
\hline
\end{tabular}
}
\caption{\small{Compression parameters of QMC cubature
with $N=216217$ Halton points on the intersection of a tetrahedral element 
with a ball as in the last figure.}}
\label{tetraball-int}
\end{center}
\end{table}


\begin{table}[ht]
\begin{center}
{\footnotesize
\begin{tabular}{| c | c | c | c |}
\hline
deg & 2 & 4 & 6\\
\hline \hline
card. $Q_c$ & 10 & 35 & 84\\
compr. ratio &5.9e+03 & 1.7e+03 & 7.0e+02\\
\hline \hline
cpu $Q_c$ & 3.3e-02 & 2.1e-02 & 5.5e-02\\
speed-up & 3.4 & 6.9 & 8.8\\
\hline
mom. resid. $Q_c$ & 5.0e-16 & 6.1e-16 & 1.2e-15\\
\hline
\end{tabular}
}
\caption{\small{Compression parameters of QMC cubature
with $N=58561$ Halton points on the difference of a tetrahedral element
with a ball as in the last figure.}}
\label{tetraball-diff}
\end{center}
\end{table}

\end{frame}



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

{\footnotesize{

\begin{enumerate}



\item M. Dessole, F. Marcuzzi, M. Vianello,  {\teals{\href{https://drna.padovauniversitypress.it/2020/1/3}{Accelerating the Lawson-Hanson NNLS solver for large-scale Tchakaloff regression designs}}}, DRNA 13, 20--29  (2020). \\
{\blues{\scriptsize{Fast Lawson-Hanson algorithm and Matlab codes.}}}

\item M. Tchernychova, Carath\'eodory cubature measures, Ph.D. dissertation in Mathematics (supervisor: T. Lyons), University of Oxford, 2015. 

\item G. Elefante, A. Sommariva, and M. Vianello, {\teals{\href{https://drna.padovauniversitypress.it/system/files/papers/08_DRNA_SA2022.pdf}{CQMC: an improved code for low-dimensional Compressed Quasi-MonteCarlo cubature}}}, DRNA 15 (2), 92--1000  (2022). \\
{\blues{\scriptsize{Compression of Quasi-Montecarlo rules.}}}

 \item A. Sommariva, {\teals{\href{https://www.math.unipd.it/\string~alvise/software.html}{Matlab codes used in the numerical experiments}}} 

\item A. Sommariva, M. Vianello, 
  {\teals{\href{https://www.tandfonline.com/doi/abs/10.1080/01630563.2015.1062394}{Compression of multivariate discrete measures and applications}}}, 
  {\em Numer. Funct. Anal. Optim.}, 36, 1198--1223 (2015). \\
{\blues{\scriptsize{(Details on cubature compression)}}}

 \item A. Sommariva, M. Vianello, {\teals{\href{https://arxiv.org/abs/2211.16620}{TetraFreeQ: tetrahedra-free quadrature on polyhedral elements}}}, submitted \\
{\blues{\scriptsize{(Meshless cubature rule over polyhedra)}}}

\item M.W. Wilson, {\teals{\href{https://www.ams.org/journals/mcom/1969-23-106/S0025-5718-1969-0242374-1/S0025-5718-1969-0242374-1.pdf}{A General Algorithm for Nonnegative Quadrature Formulas}}}, Mathematics of Computation, Vol. 23, No. 106, pp. 253--258.(1969). \\
{\blues{\scriptsize{(Pioneering paper on using Davis-Wilhelmsen theorem for determining cubature rules)}}}

\end{enumerate}

}}


\end{frame}

\end{document}


