\documentclass{siamltex}
%\documentclass{article}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{color}
\usepackage{cancel}
\usepackage{multirow, booktabs}

\hfuzz 50pt

 \title{Compressed cubature over polygonal domains \\
 {\MakeLowercase{\it \uppercase{v}ersion \uppercase{D}ecember 2 , 2018}}}
\author{
B. Bauman \thanks{Lawrence Livermore National Laboratory, Livermore, United States, ({\tt E-mail: bauman3@llnl.gov})} \and
A. Sommariva \thanks{University of Padua, Department of Mathematics, via Trieste 63, Padua, Italy, ({\tt E-mail: alvise@math.unipd.it})} \and
M. Vianello \thanks{University of Padua, Department of Mathematics, via Trieste 63, Padua, Italy, ({\tt E-mail: marcov@math.unipd.it})}  
}

\begin{document}
\maketitle
\begin{abstract}
In this paper we propose an algorithm to determine cubature rules of algebraic degree of precision $\delta$, on general polygons $\cal{P}$ by means of Matlab {\tt{polyshape}} objects and almost minimal rules on triangles, finally obtaining by Lawson and Hanson algorithm, a formula with positive weights, nodes in $\cal{P}$ and cardinality at most $N_{\delta}=(\delta+1)(\delta+2)/2$.
First, we test our algorithm on domains with different geometries and then on polygonal obscured and vignetted pupils, as required in some problems in optical astronomy.
\end{abstract}
\begin{keywords}
Numerical cubature, polygonal domains, minimal triangulation, minimal formulae, Caratheodory-Tchakaloff theorem.
\end{keywords}

\section{Introduction}

In this paper we propose how to compute cubature rules with algebraic degree of precision $ADE=\delta$ on polygons $\cal{P}$, that are in general not simple neither simply connected or connected.

These rules will have positive weights and nodes inside $\cal{P}$, i.e. of $PI$-type.

As general strategy we first triangulate the polygon using the Matlab {\tt{polyshape}} routines, using a number of triangles that is minimal or almost minimal, and then use an almost minimal rule of $PI$-type on each triangle, with $ADE=\delta$, we determine the desired rule of $PI$-type on $\cal{P}$.

Similar techniques are not new, but this one has the advantage of treating very general polygons as well as of introducing cubature formulae with nodes $\{x_i\}$ whose cardinality $M$ is in usually inferior than other general purpose codes.

Next, by solving a certain NNLS problem, by means of the Lawson-Hanson algorithm, we compute a rule of PI-type with nodes $\{t_i\} \subseteq \{x_i\}$ whose cardinality $m$ is at most $N_{\delta}=(\delta+1)(\delta+2)/2$. Such set exists in view of the well-known Caratheodory-Tchakaloff theorem.

We compare these new sets, with previous ones on some polygons, convex, concave or even not simply connected, showing the advantages of our technique.

Finally, we test our approach on polygonal obscured and vignetted pupils, useful to solve certain problems 
in optical astronomy.

\section{About the subdivision}
One of the key points to determine a rule with {\sl{few}} cubature nodes over a polygon ${\cal{P}}$, consists in its partitioning in few simpler regions $\Omega_1,\ldots,\Omega_N$ and the application of a known formula on each of them.

In \cite{SV07}, the authors considered the case of {\sl{simple}} polygons ${\cal{P}}$, i.e. without self intersections, having the so called {\sl{axis-property}} for which there exists a {\sl base-line} (say $l$), whose intersection with the polygon is connected, and such that in addition each line orthogonal to it (say $q$) has a connected intersection (if any) with the polygon ${\cal{P}}$, containing the point $l \cap q$. 
This class includes all convex polygons, for example by taking as $l$ the line connecting a pair of vertices with maximal distance but also certain nonconvex polygons. Once that this baseline is at hand, without the need of particular {\sl{triangulators}}, a formula with nodes in ${\cal{P}}$ and positive weights can be obtained via Gauss-Green theorem (with a well-established notation, we will name these rules of {\sl{PI-type}}, where $P$ stands for {\sl{positive weights}} and $I$ for {\sl{nodes in the domain}}). 

If such a baseline is not available, however a cubature rule can still be computed but without the warranty of being of PI-type (since the weights are still positive these formulae turn be of PO-type, i.e. some nodes are outside the given domain ${\cal{P}}$, though the weights are positive).


In order to enlarge the class of domains for which we can determine a rule of of PI-type, in {\cite{B09}} we implemented in Matlab a naive version of the so called {\sl{minimal triangulation}} (see also {\cite{BSA}} or {\cite{SU}} for a survey on the subject).
In few words, a triangulation of a simply connected and simple polygon ${\cal{P}}$ with more than 3 sides, is based on the two ears theorem (see \cite{ME75}), i.e on the fact that ${\cal{P}}$ has at least two {\sl{ears}}, which are triangles with two sides being the edges of ${\cal{P}}$ and the third one completely in its interior. The algorithm first computes such an ear, say $\Omega_i$, then removes it from ${\cal{P}}$, and repeats iteratively the procedure on ${\cal{P}} \backslash \Omega_i$ until only one triangle is left. The union of all the ears and of the final triangle provides a triangulation of the domain.

About the complexity of the algorithm, after several attempts by different scientists, it has been shown in  {\cite{CH91}} that the determination of such a triangulation on a simple polygon ${\cal{P}}$ with $n$ vertices costs $O(n)$ operations.

We point out that this decomposition is minimal, in the sense that for a simple and connected polygon that is also simply connected and with $n$ sides then ${\cal{P}}=\cup_{i=1}^{n-2} \Omega_i$, where $\Omega_i$ are triangles that only overlap on the sides, $n-2$ being the lower cardinality between all the possible triangulations. 

Alternatively, in \cite{G11}, it has been observed that any simple and simply connected polygon with $n$ vertices can be partitioned into $\mu$ quadrangles (some of which possibly degenerating into triangles) with $(n - 2)/2 \leq \mu \leq n -2$, where $\mu$ is often close to the lower bound. For example, any convex polygon is trivially partitioned into $\mu=(n - 2)/2$ quadrangles for $\mu$ even, or into $(\mu - 3)/2$ quadrangles plus one triangle for $\mu$ odd, simply by taking quadruples of consecutive vertices. Then by rules that are almost minimal and of PI type on the simplex or the square or alternatively of tensorial type, one can easily achieve rules on each triangle or quadrilateral subdomain of the partition, and more generally a rule of PI-type on ${\cal{P}}$.

Recently Matlab introduced the {\tt{polyshape}} environment, that manages a polygon defined by 2-D vertices via a polyshape object with properties describing its vertices, solid regions, and holes. This new class includes several facilities. For instance, it allows boolean operations, as intersection, difference, union, and symmetrical difference between polygons. Furthermore it operates rotation, scaling and translation of the given sets. 

In this framework, one can triangulate each polygon ${\cal{P}}$, that can be of very general nature, even not simply connected, not simple or disconnected. We stress that Matlab does not provide any reference about this partitioning. We tested the quality of the triangulation $\{\Omega_i\}$ on very different polygonal domains ${\cal{P}}$ with $n$ sides, achieving ${\cal{P}}=\cup_{k=1}^N \Omega_i$ with $N \sim n$ and very often $N=n-2$, as for the minimal triangulation in the case of simple, as well as simply connected polygons. The procedure turns out to be rather fast. In order to give a glimpse of the performance, we considered several polygons with a number of vertices $n$ ranging from 100 to 20000. 

As regions, we have taken 
\begin{enumerate}
\item a regular polygon ${\cal{P}}^{(1)}$ whose $n$ vertices are $$P_k=(\cos(t_k), \sin(t_k))$$ where $t_k=\frac{2\pi k}{n}$, with $k=1,\ldots,n$;
\item a polygonal cardioid ${\cal{P}}^{(2)}$ whose $n$ vertices are $$P_k=(\cos(t_k)\cdot(1-\cos(t_k)), \sin(t_k)\cdot(1-\cos(t_k)))$$ where $t_k=\frac{2\pi k}{n}$, with $k=1,\ldots,n$;
\item a polygonal Bernoulli lemniscate ${\cal{P}}^{(3)}$ whose $n$ vertices are 
$$P_k=(\sqrt{2} \cos(t_k)/(1+\sin^2(t_k),  \sqrt{2} \cos(t_k) \sin(t_k)/(1+\sin^2(t_k) ), 
$$ where $t_k=\frac{2\pi k}{n}$, with $k=1,\ldots,n$.
\end{enumerate}

\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.125,clip]{tri_ex1.png}}
   {\includegraphics[scale=0.125,clip]{tri_ex2.png}}
    {\includegraphics[scale=0.125,clip]{tri_ex3.png}}
 \caption{From left to the right, the three polygons ${\cal{P}}^{(1)}$, ${\cal{P}}^{(2}$, ${\cal{P}}^{(3)}$ for $n=32$, and their triangulation.}
 \label{fig_1}
 \end{figure}

Observe that the polygonal Bernoulli Lemniscate is a non simple domain, due to the self intersection in the origin, but can be still correctly treated until $n=1000$. From $n=2000$ on, the number of triangles becomes inferior to $n-4$ that corresponds to the cardinality of the minimal triangulation. 

We made 10 tests for each subcase, computing in each one first a polyshape object and then a triangulation. The average cputimes for each polygon ${\cal{P}}^{(k)}$, $k=1,\ldots,3$ are listed in Table {\ref{tab-00}}, respectively as $t_P$ and $t_T$, an asterisk meaning that the triangulation process presents warnings.

\begin{table}[!htbp]
\caption{Cputime of triangulation of several polygonal domains. Asterisk means that the triangulation process presents warnings. }
\begin{center}\footnotesize
\begin{tabular}{|c| c c |  c c | c c | }
\hline
\toprule
  vertices   &  \multicolumn{2}{c}{${\cal{P}}^{(1)}$}   & \multicolumn{2}{c}{${\cal{P}}^{(2)}$}   & \multicolumn{2}{c}{${\cal{P}}^{(3)}$}  \\
  \hline
     &  $t_P$   & $t_T$  & $t_P$  & $t_T$  & $t_P$  & $t_T$   \\
   \midrule 
$100$  & $4e-03$ & $2e-03$       & $3e-03$ & $1e-03$           & $3e-03$ & $1e-03$ \\
$500$  & $4e-03$ & $3e-03$      & $4e-03$ & $4e-03$         & $4e-03$ & $6e-03$ \\
$1000$  & $1e-02$ & $6e-03$     & $1e-02$ & $1e-02$          &  $1e-02$ & $2e-02$ \\
$2000$  & $4e-02$ & $7e-03$    & $5e-02^*$ & $4e-02^*$  &  $2e-02^*$ & $5e-02^*$ \\
$3000$  & $8e-02$ & $1e-02$    & $8e-02^*$ & $6e-02^*$  &  $4e-02^*$ & $1e-01^*$ \\
$4000$  & $1e-01$ & $1e-02$     & $2e-01^*$ & $9e-02^*$  &  $7e-02^*$ & $2e-01^*$ \\
$5000$  & $2e-01$ & $2e-02$    & $2e-01^*$ & $1e-01^*$   &  $1e-01^*$ & $3e-01^*$ \\
$10000$  & $8e-01$ & $4e-02$  & $1e+00^*$ & $5e-01^*$  &  $5e-01^*$ & $1e+00^*$ \\
$20000$  & $4e+00$ & $6e-02$ & $4e+00^*$ & $1e+00^*$ &  $2e+00^*$ & $3e+00^*$ \\
   \bottomrule
    
\end{tabular}
\end{center}
\label{tab-00}
\end{table}

All tests have been performed on a 2,7 GHz Intel Core i5 with 16 GB 1867 MHz DDR3 memory.

The results say that for polygons with less than $n \leq 3000$ vertices the triangulation cputime is negligible, while for $n > 3000$ the time $t_P$ needed for the polyshape construction  can be even dominant w.r.t. $t_T$.

\section{On the cubature nodes in the subdivision}

Once a subdivision ${\cal P}=\cup_{i=1}^n \Omega_i$ is at hand, in order to have a cubature rule of PI-type on $\cal P$, by the additivity of the integral operator, it is sufficient to define a rule of PI-type on each subdomain $\Omega_i$, that usually are triangles or quadrilaterals. 

 In the first instance, relevant in this work, knowing a formula on a reference triangle $\cal{T}^{*}$, after the conversion of the nodes in barycentric coordinates, it is straightforward to achieve one on each possible triangle $\cal{T}$ by varying the weights proportionally to their area. In other words, suppose that $\cal{T}^{*}$ is the unit-triangle with vertices $(0,0)$, $(0,1)$, $(1,0)$, $\{\phi_k\}_{k=1,\ldots,N_{\delta}}$ is a basis of the space of polynomials ${\mathbb P}_{\delta}$ of total degree $\delta$ and cardinality $N_{\delta}=(\delta+1)(\delta+2)/2$, e.g. the Proriol-Dubiner basis \cite{BP}, and that
\begin{equation}\label{VP105}
\sum_{i=1}^m  w^*_i \phi_k(\xi^*_i) = \int_{\cal{T}^{*}} \phi_k d{\cal{T}^{*}}, \, \, \xi_i \in {\cal{T}^{*}}, {\mbox{   }}k=1,\ldots,N_{\delta}.
\end{equation}
has algebraic degree of exactness $ADE=\delta$, i.e. able to integrate exactly each algebraic polynomial of total degree $\delta$ over the region ${\cal{T}^{*}}$.

Furthermore let $\xi^*_i=(x^*_i,y^*_i)$ and denote by ${\hat{\xi^*_i}}=(x^*_i,y^*_i,1-x^*_i-y^*_i)$ the barycentric coordinates of $\xi^*_i$, by $\mu(\cal{T})$ the area of $\cal{T}$, and by $V_1$, $V_2$, $V_3$ its vertices. Being $\mu({\cal{T}}^*)=1/2$, the cubature rule on $\cal{T}$ with nodes
$$
{\tilde{\xi}}_i=x^*_i \cdot V_1+y^*_i \cdot V_2+(1-x^*_i-y^*_i) \cdot V_3, \,\, {\mbox{   }} i=1,\ldots,m
$$
and weights
$$
w_i=\frac{\mu(\cal{T})}{\mu(\cal{T^*})} w^*_i =2{\mu(\cal{T})} w^*_i,  \,\,  {\mbox{   }} i=1,\ldots,m
$$
has $ADE=\delta$ on $\cal{T}$, i.e. integrates exactly each polynomial in ${\mathbb{P}}_{\delta}$ on $\cal T$.

In view of this well-known result, it is necessary to compute formulae with a certain ADE only on the reference triangle ${\cal{T}}^*$.

In this context, the more appealing ones to achieve a low cardinality rule are the so called {\sl{minimal rules}} that have the lowest number of nodes between those whose algebraic degree of exactness $ADE$ is $\delta$. Differently from the univariate case, where they are the so called {\sl{Gaussian rules}}, they are known only in few cases (see, e.g. {\cite{CO97}}).

In spite of this, the bibliography about rules with low cardinality $\cal{T}$,  with a certain algebraic degree of exactness (i.e. {\sl{almost minimal}}, since in general it is quite difficult to prove that they are minimal), is rather wide.
A common technique for their determination is the following. If $\{\phi_k\}_{k=1,\ldots,N_{\delta}}$ is a polynomial basis on $\Omega^*$ and all the moments $\gamma_k=\int_{\Omega^*} \phi_k d\Omega^*$, $k=1,\ldots,N_{\delta}$ are known, one computes via numerical optimisation the solutions $(\xi_k,w_k)$, with $k=1,\ldots,M$ of the problem (\ref{VP105})
with $M \leq N_{\delta}$ as low as possible. This task is not an easy one, especially when the number of equations $N_{\delta}$ becomes large (see, e.g., \cite{TWB07}). In order to lower the number of equations, symmetries of the nodes allow to search the optimal set in a less general family, as well as to solve smaller nonlinear systems. Variants of this successful strategy provided many almost minimal pointsets, even for rather large $\delta$ (see, e.g. \cite{D85}, \cite{LJ75}, \cite{PP15}, \cite{PP16}, \cite{XG10}, \cite{WX03}).

In Table {\ref{tab-01}} we have listed to the authors' knowledge the best of these almost minimal rules between those of PI-type. that are used to implement a formula with algebraic degree of precision $ADE=\delta$, not only on ${\cal{T^*}}$ but also any ${\cal{T}}$ and more generally, by the additivity of the integral operator, on any polygon ${\cal{P}}$.

\begin{table}[!htbp]
\caption{Cardinality $N^*_{\delta}$ of (almost) minimal rules on triangles with $ADE=\delta$.}
\begin{center}\footnotesize
\begin{tabular}{|| c | c | | c  | c | | c  | c | | c | c | | c | c ||}
\hline

  \hline
    $\delta$ &  $N^*_{\delta}$   &  $\delta$ &  $N^*_{\delta}$ &  $\delta$ &  $N^*_{\delta}$ &  $\delta$ &  $N^*_{\delta}$  &  $\delta$ &  $N^*_{\delta}$  \\
   \hline 
  	 $ 1$ & $  1$ & $11$ & $ 27$ & $21$ & $ 85$ & $31$ & $181$ & $41$ & $309$ \\
 	 $ 2$ & $  3$ & $12$ & $ 32$ & $22$ & $ 93$ & $32$ & $193$ & $42$ & $324$ \\
 	 $ 3$ & $  4$ & $13$ & $ 36$ & $23$ & $100$ & $33$ & $204$ & $43$ & $339$ \\
 	 $ 4$ & $  6$ & $14$ & $ 42$ & $24$ & $109$ & $34$ & $214$ & $44$ & $354$ \\
 	 $ 5$ & $  7$ & $15$ & $ 46$ & $25$ & $117$ & $35$ & $228$ & $45$ & $370$ \\
 	 $ 6$ & $ 11$ & $16$ & $ 52$ & $26$ & $130$ & $36$ & $243$ & $46$ & $385$ \\
 	 $ 7$ & $ 12$ & $17$ & $ 57$ & $27$ & $141$ & $37$ & $252$ & $47$ & $399$ \\
 	 $ 8$ & $ 16$ & $18$ & $ 66$ & $28$ & $150$ & $38$ & $267$ & $48$ & $423$ \\
 	 $ 9$ & $ 19$ & $19$ & $ 70$ & $29$ & $159$ & $39$ & $282$ & $49$ & $435$ \\
 	 $10$ & $ 24$ & $20$ & $ 78$ & $30$ & $171$ & $40$ & $295$ & $50$ & $453$ \\
	 \hline
\end{tabular}
\end{center}
\label{tab-01}
\end{table}

All the rules are stored as a Matlab file at \cite{Ssoft}, where we have corrected the formulae whenever necessary in order to match (\ref{VP105}) close to machine precision.

For higher degrees one can use the well-known Stroud conical rules  \cite{LC95}, whose cardinality is $(\delta+1)^2/4$ for odd $\delta$.

{\vspace{0.5cm}}

{\bf{Note:}}
We point out that though it is not the purpose of the paper, many authors determined adaptive algorithms for numerical integration of functions defined on triangles (see, e.g. \cite{BE92}, \cite{DR84}, \cite{G90}, \cite{K92}, \cite{KR87}, \cite{L82}, \cite{T09} and references therein). In particular one can take advantage of the vectorizing the procedure in Matlab as explained in {\cite{SH08}}, for integration on the square to implement a variant for adaptive cubature over polygons. It is our intention to investigate later this subject.


{\bf{Note:}} 
Observe that one can use any triangulation of ${\cal P}$ to determine a formula of PI-type on such domain, but the minimal one has the advantage of keeping as low as possible the cardinality of the cubature formula, since in any triangle a fixed number of points is used.

\section{Caratheodory-Tchakaloff subsampling}

The purpose of this section is to show how from a rule with high cardinality, positive weights, and $ADE=\delta$ we can extract another with the same degree of precision, positive weights, but with cardinality equal at most to the dimension $N_{\delta}$ of the vector space ${\mathbb P}_{\delta}$.

In some recent papers, the authors have applied a mathematical tool named CATCH (acronym of (Caratheodory-Tchakaloff) subsampling), for such compression of discrete measures, proposing its application to discrete polynomial Least squares by sparse moment matching.  In this framework, the method selects from a large discretization of a given region a much smaller number of (weighted) sampling points, even on a complex shape as can be the polygons investigated in this paper, keeping numerically invariant the Least Squares approximation estimates.

The key point is the following discrete version of Tchakaloff theorem \cite{T57} proved by Caratheodory theorem on finite-dimensional conic/convex combinations \cite{C07}, \cite{PSV17},

\vspace{0.25cm}

\begin{theorem}\label{TC167}
Let $\mu$ be a multivariate measure whose support is a ${\mathbb{P}}_{\delta}$-determining finite set $X = \{\xi_i\} \subset {\mathbb{R}}^d$ (i.e., $n$-degree polynomials vanishing there vanish ev- erywhere), with correspondent positive weights (masses) ${\bf{\lambda}} = \{\lambda_i\}$, $i = 1,\ldots,M$, $M = card(X) > N_{\delta}$, $N_{\delta}$ being the dimension of ${\mathbb{P}}_{\delta}$.
Then, there exist a cubature formula for the discrete measure $\mu$, with nodes $T_n = \{t_j\} \subset X$ and positive weights ${\bf{w}} = \{w_i\}$, $1 \leq j \leq m$, with $m \leq N_{\delta}$, such that
$$
\int_X p(x) d\mu = \sum_{i=1}^M \lambda_i p(x_i) = \sum_{j=1}^m w_j p(t_j), {\mbox{  }} \forall p \in {\mathbb{P}}_{\delta}.
$$
\end{theorem}

From the numerical side, given any polynomial basis $\{\phi_k\} $ of $ {\mathbb{P}}_{\delta}$,  define the Vandermonde-like matrix $V=\{V_n(X)\}_{i,j}=\phi_j(\xi_i)$  and let $\gamma$ the vector of moments of the polynomial basis $\{\phi_j\}$ with respect to the original discrete measure and consider the underdetermined moment system $V^Tu=\gamma$.

Caratheodory/Tchakaloff theorem \ref{TC167} asserts that there exists a sparse nonnegative solution $u$ to the system above, whose nonvanishing components (i.e., the weights $\{w_j \}$) are at most $N_{\delta}$ and determines the corresponding reduced sampling points $T_\delta = \{t_j\} \subseteq X$, that we may term the Caratheodory-Tchakaloff (CATCH) points of X.

The computation of these nodes has been considered in several papers (see e.g. \cite{SV15},  \cite{T15}). To our knowledge, essentially two approaches have been developed to get these rules, i.e. via Linear Programming (LP) and Quadratic Programming (QP).  

About the LP approach, it consists in solving via simplex-method
\begin{equation}
\left\{
\begin{array}{l}
\min_{u\geq 0} c^T u \\
V^Tu=b, {\hspace{0.5cm}} u \geq 0
\end{array}
\right.
\end{equation}

\noindent
where the constraints identify a polytope (the feasible region) in ${\mathbb{R}}^M$ and the vector $c$ is chosen to be linearly independent from the rows of $V^T$, so that the objective functional is not constant on the polytope \cite{PSV17}, \cite{T15}.

The QP based algorithm requires instead the solution of the NonNegative Least Squares (NNLS) problem
\begin{equation}\label{4.2}
{\mbox {compute }} u^* : \|Au^* - \gamma\|_2 =\min_u \|Au - \gamma\|_2 , u \geq 0, 
\end{equation}
in which $u^*$ can be obtained by the well-known {\sl{Lawson-Hanson}} active set optimization method {\cite{LH95}}, which determines a sparse solution to ({\ref{4.2}}). Its application gives a residual $\epsilon = \|Au^*- \gamma\|_2$ that is typically very small, say $< 10^{-14}$ for $\delta \leq 30$. 

Our numerical experience with now available Matlab software has shown that NNLS usually performs better than LP
in computing the CATCH weights, at least for moderate degrees $\delta$ (namely, $N_{\delta}$
for mild degrees) \cite{PSV17}. Consequently all our codes are based on this the application of Lawson-Hanson method to compute the Caratheodory/Tchakaloff sets.

We point out that there are several versions of NNLS codes available in Matlab. One is the built-in function {\tt{lsqnonneg}}, based on the Lawson-Hanson algorithm while an open-source version present in the package NNLSlab in \cite{SH}. Other alternatives are often obtained by MEX files and for generality purposed will not be used here.

{\vspace{0.25cm}}

{\bf{Note.}}  As alternative, one can compress the rule via the QR algorithm with column pivoting proposed in 1965 by Businger
and Golub {\cite{BG65}}, which is implemented for example by the Matlab {\sl{backslash}}
operator. This approach is equivalent to a greedy selection of the columns of $A$ in order to
maximize the successive volumes, and eventually the (absolute value of the)
determinant (the column selection problem being NP-hard). If the
matrix has full rank, the final result is a weight vector $w^*$ where only $N \leq N_{\delta}=(\delta+1)(\delta+2)/2$ components are nonzero, so that we can extract the cubature nodes $T_\delta = \{t_j\}$ from $X$ by
the column indexes corresponding to such components. A drawback is that the
resulting weights $w^*=\{w^*_k\}$ may not be all positive, but typically the
negative ones are few and of small size, so that the relevant stability parameter
$$
\sigma=\frac{\sum_{k=1}^N |w^*_k|}{|\sum_{k=1}^N w^*_k|} \geq 1
$$
is not far from $1$ (see also {\cite{SV15}}).

{\vspace{0.25cm}}

{\bf{Note.}} In the cubature framework an algorithm termed {\sl{Recursive
Halving Forest}}, based on a hierarchical SVD, has been proposed in \cite{T15}. Performances are reported for large scale problems (say that the order of $N_{\delta}$ is  $10^3$, $10^4$).  Unfortunately the software is not available and thus cannot be applied here as comparison.

\section{Numerical experiments}

As for the numerical experiments, we intend first to compare the results that we can obtain via this new algorithm with those in {\cite{SV07}} and next to show the effect of Caratheodory-Tchakaloff subsampling in the compression of a cubature rule. Finally we propose a method to obtain embedded rules on complicated polygonal domains.

In  {\cite{SV07}}, the numerical experiments where made respectively on a convex and on certain non-convex polygon, say ${{\Omega}}^{(C)}$, ${{\Omega}}^{(NC)}$, achieving rules of PI-type (see {\ref{fig_1a}} for the description of the polygons, baselines and pointset distribution for $\delta=10$). 

\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.125,clip]{2007C.png}}
   {\includegraphics[scale=0.125,clip]{2007NC.png}}
 \caption{From left to the right, the domains ${{\Omega}}^{(C)}$, ${{\Omega}}^{(NC)}$. The red line is the {\sl{axis}} used by the algorithm in {\cite{SV07}}. The pointset is relative to $ADE=\delta=10$.}
 \label{fig_1a}
 \end{figure}

In Table {\ref{tab-02}}, we compare the two approaches, as well as the cardinality of each cubature pointset, with cardinalities $M_{\delta}^{(old)}$, $M_{\delta}^{(new)}$, ${\tilde{M}}_{\delta}^{(new)}$ respectively applying the algorithm proposed in {\cite{SV07}}, the one presented here before and after the compression. 

\begin{table}[!htbp]
\caption{Cardinality $M_{\delta}^{(old)}$, $M_{\delta}^{(new)}$,  ${\tilde{M}}_{\delta}^{(new)}$, respectively applying the algorithm proposed in {\cite{SV07}}, the one presented here before and after the compression, on two polygons ${{\Omega}}^{(C)}$, ${{\Omega}}^{(NC)}$, with $ADE=\delta=5,10,\ldots,40$.}
\begin{center}\footnotesize
\begin{tabular}{|c| c c c | c c c |}
\hline
\toprule
    &  \multicolumn{3}{c}{${{\Omega}}^{(C)}$}   & \multicolumn{3}{c}{${{\Omega}}^{(NC)}$}   \\
  \hline
  $\delta$   & $M_{\delta}^{(old)}$  &  $M_{\delta}^{(new)}$ &  ${\tilde{M}}_{\delta}^{(new)}$  & $M_{\delta}^{(old)}$ &  $M_{\delta}^{(new)}$ &  ${\tilde{M}}_{\delta}^{(new)}$      \\
   \midrule 
         $ 5$  & $ 180$ &  $  28$   & $21$         & $235$     & $49$ & $ 21$  \\
 	 $10$& $ 660$  & $  96$   & $66$          & $870$  & $168$ & $ 66$  \\
 	 $15$ & $ 1440$  & $ 184$  & $136$    & $1905$ & $322$ & $136$   \\
 	 $20$ & $ 2520$ & $  312$ & $231$    & $3340$ &  $546$ & $231$   \\
 	 $25$ & $ 3900$ & $  468$ & $351$    & $5175$    & $819$ & $351$ \\
 	 $30$ & $ 5580$ & $ 684$  & $496$    & $7410$   & $1197$ & $496$  \\
 	 $35$ & $ 7560$ & $ 912$  & $666$     & $10045$  & $1596$ & $666$  \\
 	 $40$ & $ 9840$ & $ 1180$ & $861$      & $13080$ & $2065$ & $861$  \\
	 \hline
   \bottomrule
    
\end{tabular}
  	

\end{center}
\label{tab-02}
\end{table}

From this table, it is evident compression via NNLS provides a cubature pointset with cardinality $(\delta+1)(\delta+2)/2$, i.e. the dimension of the vector space ${\mathbb{P}}_{\delta}$ and that the compression ratio $M_{\delta}^{(new)}/{\tilde{M}}_{\delta}^{(new)}$ depends on the number of sides of the polygon.

As for the cputime, the computation of any of these rules is not an issue since their determination takes less than $10^{-2}$ seconds.  Compression using Matlab built-in routine {\tt{lsqnonneg}} or the open-source NNLSlab is in general not too time consuming for low degrees, while it becomes relevant when $\delta \geq 20$.

In Table \ref{tab-03}, we approximate the moments $\{\gamma_k\}$ of the basis $\{\phi_k\}$ consisting of the product of scaled Chebyshev polynomials on the smallest rectangle with sides parallel to the cartesian axes, containing the domain in analysis and whose total degree is inferior or equal to $\delta$.
In particular, in such polynomial basis, first we evaluate the moments $\{\gamma^{(T)}_k\}$ by the rule via triangulation and then compare them with those of the formulae obtained by compression (say $\{ {\tilde{\gamma}}^{(T)}_k \}$) and Gauss-Green approach (say $\{\gamma^{(GG)}_k\}$), reporting $E_{T,TC}:=\|\gamma^{(T)}_k-{\tilde{\gamma}}^{(T)}_k\|_2$, $E_{T,GG}=\|\gamma^{(T)}_k-\gamma^{(GG)}_k\|_2$.

For the compression we used in particular the software available from the package NNLSlab in \cite{SH}.

\begin{table}[!htbp]
\caption{Comparison of the moments of cubature rules obtained resp. via triangulation (with/without compression) and Gauss-Green approach on two polygons ${{\Omega}}^{(C)}$, ${{\Omega}}^{(NC)}$, with $ADE=\delta$.}
\begin{center}\footnotesize
\begin{tabular}{|c| c c | c c |}
\hline
\toprule
    &  \multicolumn{2}{c}{${{\Omega}}^{(C)}$}   & \multicolumn{2}{c}{${{\Omega}}^{(NC)}$}   \\
  \hline
  $\delta$   &  $E_{T,TC}$ &  $E_{T,GG}$   & $E_{T,TC}$  &  $E_{T,GG}$    \\
   \midrule 
$ 5$ & $3e-16$ & $3e-16$    & $2e-16$ & $2e-16$ \\
$10$ & $3e-16$ & $3e-16$   & $3e-16$ & $3e-16$  \\
$15$ & $4e-16$ & $4e-16$    & $5e-16$ & $5e-16$ \\
$20$ &$8e-16$ & $8e-16$    & $6e-16$ & $6e-16$ \\
$25$ & $6e-16$ & $6e-16$   & $1e-15$ & $1e-15$ \\
$30$ & $2e-15$ & $2e-15$   & $1e-15$ & $1e-15$ \\
$35$ & $1e-15$ & $1e-15$   & $2e-15$ & $2e-15$ \\
$40$ & $2e-15$ & $2e-15$   & $2e-15$ & $2e-15$  \\
	   
	    \hline
   \bottomrule
    
\end{tabular}
  	
\end{center}
\label{tab-03}
\end{table}
	
The estimates in Table {\ref{tab-03}} are relevant, since in \cite{SV15}  we proved that if $f \in C(\Omega)$ and $\lambda$ is a measure over $\Omega$ (e.g. the scaled Chebyshev measure mentioned above or even a fully discrete one), $\nu$ the discrete measure such that $\int_{\Omega} f(P) d\nu=\sum_k w_k f(P_k)$, ${\mathbb{P}}_{\delta}={\mbox{span}}(\{\phi_k\})$, then denoting by $E_{\delta}(f):=\min_{p \in {\mathbb{P}}_{\delta}} \|f-p\|_{\infty}$, $\epsilon_{mom}:=\|{\bf{m}}-{\bf{m^{\ast}}}\|_2, $
where ${\bf{m}}_k:=\int_{\Omega} \phi_k(P) d\nu$, ${\bf{m_k^{\ast}}}:= \sum_j w_j\phi_k(Q_j)$ we have
$$
\left|\int_{\Omega} f(P) d\nu- \sum_k w_k f(Q_k)\right| \leq C E_{\delta}(f) +\|f\|_{L^2_{\lambda(\Omega)}} \epsilon_{mom}, {\mbox{ }} \forall f\in C(\Omega)
$$
in which 
$$C \leq 2 (\nu(\Omega)+\sqrt{\lambda(\Omega)} \epsilon_{mom})$$
if the cubature rule has nonnegative weights, as in our instances.

As consequence, we can expect that the behaviours of the three rules mentioned above is numerically {\sl{equivalent}} when integrating a function $f \in C(\Omega)$, since all the quantities $\epsilon_{mom}$ are very small.


To show the advantages of this new approach via polyshape triangulation, we consider as examples two polygons ${{\Omega}}^{(QF_1)}$, ${{\Omega}}^{(QF_2)}$, that we can think as discretization of a quatrefoil, that could not be treated by the previous algorithms since the domains are not simple and have not the axis property.

More precisely,
$$
{{\Omega}}^{(QF_i)}=(\cos(t_k) \cdot \sin(2t_k), \sin(t_k) \cdot \sin(2t_k)))
$$
where $t_k=\frac{2k\pi}{M}$, with $k=1,\ldots,M$. In particular, we set $M=129$, and $M=513$ respectively for ${{\Omega}}^{(QF_1)}$, ${{\Omega}}^{(QF_2)}$. The polygon ${{\Omega}}^{(QF_1)}$ can be partioned in $N^{(QF_1)}_{tri}=120$ triangles, while ${{\Omega}}^{(QF_2)}$ in $N^{(QF_2)}_{tri}=504$ triangles.



\begin{table}[!htbp]
\caption{Comparison of the moments of cubature rules obtained via triangulation (with/without compression) as well as their cardinalities $M_{\delta}^{(new)}$,  ${\tilde{M}}_{\delta}^{(new)}$, on two not simple polygons ${{\Omega}}^{(QF_1)}$, ${{\Omega}}^{(QF_2)}$, with $ADE=\delta=5,10,\ldots,35$.}
\begin{center}\footnotesize
\begin{tabular}{|c| c c c c | c c c c |}
\hline
\toprule
    &  \multicolumn{4}{c}{${{\Omega}}^{(QF_1)}$}   & \multicolumn{4}{c}{${{\Omega}}^{(QF_2)}$}   \\
  \hline
  $\delta$   &  $E_{T,TC}$ &  $M_{\delta}^{(new)}$  &  ${\tilde{M}}_{\delta}^{(new)}$ &  ${\mbox{Ratio}}_{\delta}$   &  $E_{T,TC}$ &  $M_{\delta}^{(new)}$  &  ${\tilde{M}}_{\delta}^{(new)}$ &  ${\mbox{Ratio}}_{\delta}$     \\
   \midrule 
$  5$ & $1e-15$ & $    840$ & $     21$ & $   40.0$        & $1e-15$ & $   3528$ & $     21$ & $  168.0$ \\
$ 10$ & $1e-15$ & $   2880$ & $     66$ & $   43.6$      & $2e-15$ & $  12096$ & $     66$ & $  183.3$ \\
$ 15$ & $2e-15$ & $   5520$ & $    136$ & $   40.6$     & $2e-15$ & $  23184$ & $    136$ & $  170.5$\\
$ 20$ & $3e-15$ & $   9360$ & $    231$ & $   40.5$     & $4e-15$ & $  39312$ & $    231$ & $  170.2$ \\
$ 25$ & $1e-14$ & $  14040$ & $    351$ & $   40.0$    &$6e-15$ & $  58968$ & $    351$ & $  168.0$ \\
$ 30$ & $9e-15$ & $  20520$ & $    496$ & $   41.4$    & $1e-14$ & $  86184$ & $    496$ & $  173.8$ \\
$ 35$ & $9e-15$ & $  27360$ & $    666$ & $   41.1$    & $1e-14$ & $ 114912$ & $    666$ & $  172.5$ \\
% $ 40$ & $1e-14$ & $  35400$ & $    861$ & $   41.1$    & $5e-02^*$ & $ 148680$ & $    534^*$ & $  278.4^*$ \\
	   
	    \hline
   \bottomrule
    
\end{tabular}
  	
\end{center}
\label{tab-04}
\end{table}


 \begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.15,clip]{convex_polygon.png}}
   {\includegraphics[scale=0.15,clip]{nonconvex_polygon.png}}
   {\includegraphics[scale=0.15,clip]{quaterfoil.png}}
 \caption{Compressed set at degree $\delta=10$, for the polygons ${{\Omega}}^{(C)}$, ${{\Omega}}^{(NC)}$, ${{\Omega}}^{(QF_1)}$.}
 \label{3434}
 \end{figure}

By a careful look at Table {\ref{tab-04}}, we observe that the compression ratio $${\mbox{Ratio}}_{\delta}={{M}}_{\delta}^{(new)}/{\tilde{M}}_{\delta}^{(new)}$$ is almost constant. The reason is that at degree $\delta$ the cubature rule has approximatively $\delta^2/6$ points for triangle, thus $M_{\delta}^{(new)}\approx N^{(QF_i)}_{tri} \delta^2/6$ and since the compressed set has cardinality ${\tilde{M}}_{\delta}^{(new)}\approx \delta^2/2$, necessarily
$$
{\mbox{Ratio}}_{\delta}=\frac{M_{\delta}^{(new)}} {{\tilde{M}}_{\delta}^{(new)}} \approx \frac{N^{(QF_i)}_{tri} \delta^2/6}{\delta^2/2}=\frac{N^{(QF_i)}_{tri}}{3}.
$$

At a first glance the value $E_{T,TC}$  seems to deteriorate increasing $\delta$, but one has to consider that the dimension $N_{\delta}=(\delta+1)(\delta+2)/2$ of the vector of the moments increases, so the {\sl{relative}} norm two $E_{T,TC}/\sqrt{N_{\delta}}$ is not too far from to machine precision. 


% In general for rather high degrees, say $\delta \geq 30$, compression may have problems, as in the case of the second domain ${{\Omega}}^{(Q_2)}$ for $\delta=40$, where the suggested pointset was too small and not % matching the wanted moments. We report that using {\tt{lsqnonneg}} the results were not improving i.e. $E_{T,TC} \approx 3e-02$ using $601$ points. Applying instead the QR based method, the error is $E_{T,TC} 
% \approx 1e-14$, determining correctly $861$ points, of which $201$ possess negative weights, but whose sum is $14$ per cent of the measure of the domain, more precisely $-2.2e-01$ of $1.57e+00$, i.e. still providing % a stable NI-type cubature rule since $\sigma=\frac{\sum_{k=1}^N |w_k|}{|\sum_{k=1}^N w_k|} \approx 1.28$ where $w_k$ are the weights of the compressed rule by QR method.
  
We finally observe that by Caratheodory-Tchakaloff subsampling we can determine embedded rules over a polygon $\cal P$, i.e. cubature formulae $Q^{(L)}$, $Q^{(H)}$,
$$
Q^{(L)}(f)=\sum_{k=1}^{M_L} w^{(L)}_k f(P^{(L)}_k), {\mbox{ }} Q^{(H)}(f)=\sum_{k=1}^{M_H} w^{(H)}_k f(P^{(H)}_k)
$$ 
with different degree of precision, say $\delta_{Q^{(L)}} < \delta_{Q^{(H)}} $, where $\{P^{(L)}_k\} \subseteq \{P^{(H)}_k\}$ (see also  \cite{SV18}).

A typical application is that they allow stopping criteria, based on rules with different degrees of precision $\delta_1$, $\delta_2$, with $\delta_1 < \delta_2$, with the property that the the nodes of the rule with $ADE=\delta_1$ are between those of the rule with $ADE=\delta_2$. Their purpose is to minimize the number of function evaluations, for providing an estimate of the cubature error at degree $\delta_1$.

More generally, given a sequence of degrees $\delta_1<\delta_2<\dots <\delta_k$, 
we can compute the nested CATCH sequence $\{T_{\delta_j}\}$ 
\begin{equation} \label{nested}
{\cal P} \supset T_{\delta_k}\supset T_{\delta_{k-1}}\supset \dots\supset T_{\delta_2}
\supset T_{\delta_1}
\end{equation}
together with the corresponding sequence of positive weight vectors, 
say $\{{\bf{w}}_{\delta_j}\}$, 
by solving backward the sequence of NLLS problems
\begin{equation} \label{nestedNNLS}
\mbox{compute}\;{\bf{u}}_{j-1}^\ast\,:\;\|A_j{\bf{u}}_{j-1}^\ast-{\bf{b}}_j\|_2
=\min\|A_j{\bf{u}}_{j-1}-{\bf{b}}_j\|_2\;,\;\;{\bf{u}}_{j-1}\geq 
{\bf{0}}\;,
\end{equation}
for $j=k+1,k,\dots,2$, 
where $A_j=(V_{\delta_{j-1}}(T_{\delta_j}))^t$, 
${\bf{b}}_j=A_j{\bf{u}}_j^\ast$, and we set 
$T_{\delta_{k+1}}=X$, ${\bf{u}}_{k+1}^\ast={{\bf{\lambda}}}$.




   
\section{Application to polygonal obscured and vignetted pupils}
As application of these new cubature formulae, we consider a problem arising in optical design. 

In \cite{F88}, Forbes suggested to use product Gaussian quadrature over bivariate domains $\Omega$ that were filled, circular or
elliptical apertures (pupils), to compute root-mean-square (rms) spot size for an optical design.

Later, Bauman and Xiao \cite{BX10} studied cubature methods based on prolate spheroidal wave functions \cite{XRY03}, to treat situations where the pupil is obscured and vignetted (a feature that occurs, for example, in optical astronomy; see figure {\ref{fig_1b}}). Indeed, Large Synoptic Survey Telescope (LSST) has a large central obscuration (about 62 per cent obscuration by diameter) as well as a considerable vignetting (of up to 10 per cent by area) making the techniques developed by Forbes not appropriate for the problem.


\begin{figure}[!htbp]
  \centering
% {\includegraphics[scale=0.175,clip]{disks_intersection_lens.jpg}}
% {\includegraphics[scale=0.20,clip]{disks_intersection_lens_10.jpg}}
{\includegraphics[scale=0.4,clip]{lens_dec18.eps}}
 \caption{We have drawn the circles determining the domain ${{\Omega}}^{(M_1)}$, i.e. a polygonal obscured and vignetted pupil, the compressed pointset for $\delta=10$ is determined by the green circles, while the initial one by small black dots.}
 \label{fig_1b}
 \end{figure}

 In \cite{DFSV13} the authors took into account an example of circular pupil (the unit disk) which is obscured by a central smaller disk and clipped by a circular arc of larger radius, similar to that appearing in \cite{BX10}, providing algebraic rules of PI-type with $ADE=\delta$.
 
In all those examples the pupils where circular, but in real applications they can also have a different geometry, e.g. regular polygons with many sides. To this purpose, we introduce here a novel approach to determine rules of PI-type having $ADE=\delta$ on domains that are based on {\sl polygonal} pupils.
 
As example, we take into account systems where the polygon ${\cal{P}}$ is the set subtraction of two polygonal regions ${\cal{P}}^{(O)}$, ${\cal{P}}^{(I)}$ where

\begin{enumerate}
\item the polygon ${\cal{P}}^{(O)}$ is the intersection of other three regular polygons ${\cal{P}}^{(O)}_1$, ${\cal{P}}^{(O)}_2$, ${\cal{P}}^{(O)}_3$, sharing the same center but with different apothems; 
\item ${\cal{P}}^{(I)}$ consists in the union of two regular and polygons ${\cal{P}}^{(I)}_1$, ${\cal{P}}^{(I)}_2$ whose intersection is not empty;
\item ${\cal{P}}^{(I)} \subset {\cal{P}}^{(O)}$.
\end{enumerate}




In general these operations are not easy to perform, since it is not straighforward to compute the intersection, the union and the subtraction of polygons, but the new release of Matlab R2017b, allows the solution of this problem via the built-in routines {\tt{intersect}}, {\tt{union}}, {\tt{subtract}} applicable to polyshape objects. This allows to analyze complicated systems of lenses, even more general than ours and typical of optical design.

In the test, we consider regular polygons ${{\Omega}}^{(M_1)}$, ${{\Omega}}^{(M_2)}$
\begin{itemize}
\item ${\cal{P}}^{(O)}_1$ with center $C^{(O)}_1=(0,0)$ and apothem $r^{(O)}_1=1$,
\item ${\cal{P}}^{(O)}_2$ with center $C^{(O)}_2=(0,0)$ and apothem $r^{(O)}_2=0.6120$,
\item ${\cal{P}}^{(O)}_3$ with center $C^{(O)}_3=(0,-0.1184)$ and apothem $r^{(O)}_3=0.5663$,
\item ${\cal{P}}^{(I)}_1$ with center $C^{(I)}_1=(0,-0.1184)$ and apothem $r^{(I)}_1=1.0761$,
\item ${\cal{P}}^{(I)}_2$ with center $C^{(I)}_2=(0,-0.3761)$ and apothem $r^{(I)}_1=1.2810$,
\end{itemize}
each one with vertices 
$$
C^{(s)}_k+r^{(s)}_k (\cos(t_j),\sin(t_j)), {\mbox{ }} t_j=\frac{2j\pi}{M_i}, {\mbox{ }} j=1,\ldots,M_i
$$
being $s=O$ and $k=1,2,3$ or  $s=I$ and $k=1,2$ and $M_i$ depending on ${{\Omega}}^{(M_i)}$ (i.e. $100$ and $500$).
The polygons ${{\Omega}}^{(M_1)}$, ${{\Omega}}^{(M_2)}$, are partitioned respectively in $202$ and $998$ triangles.

 
\begin{table}[!htbp]
\caption{Comparison of the moments of cubature rules obtained via triangulation (with/without compression) as well as their cardinalities $N_{T}$,  $N_{TC}$, on two polygons ${{\Omega}}^{(M_1)}$, ${{\Omega}}^{(M_2)}$, with $ADE=\delta$, $M_1=100$, $M_2=500$.}
\begin{center}\footnotesize
\begin{tabular}{|c| c c c c | c c c c |}
\hline
\toprule
    &  \multicolumn{4}{c}{${{\Omega}}^{(M_1)}$}   & \multicolumn{4}{c}{${{\Omega}}^{(M_2)}$}   \\
  \hline
  $\delta$   &  $E_{T,TC}$ &  $N_{T}$  &  $N_{TC}$ &  ${\mbox{Ratio}}_{T,TC}$   &  $E_{T,TC}$ &  $N_{T}$  &  $N_{TC}$ &  ${\mbox{Ratio}}_{T,TC}$     \\
   \midrule 
 	 $  5$ & $1e-15$ & $   1414$ & $     21$ & $   67.3$ & $2e-15$ & $   6986$ & $     21$ & $  332.7$ \\
 	 $ 10$ & $3e-15$ & $   4848$ & $     66$ & $   73.5$ & $1e-14$ & $  23952$ & $     66$ & $  362.9$  \\ 
 	 $ 15$ & $5e-15$ & $   9292$ & $    136$ & $   68.3$ & $3e-15$ & $  45908$ & $    136$ & $  337.6$  \\
 	 $ 20$ & $5e-15$ & $  15756$ & $    231$ & $   68.2$ & $9e-15$ & $  77844$ & $    231$ & $  337.0$\\
 	 $ 25$ & $9e-15$ & $  23634$ & $    351$ & $   67.3$ & $9e-15$ & $ 116766$ & $    351$ & $  332.7$ \\
	    \hline
   \bottomrule
    
\end{tabular}
  	
\end{center}
\label{tab-05}
\end{table}

These computations are performed to approximate integrals of certain algebraic polynomials necessary in the determination of lens aberration, and in general the degrees of exactness of interest are $\delta \leq 25$. We report that we experimented moments matching deterioration with the compression for higher degrees with Matlab NNLS code {\tt{lsqnonneg}}, while these problems where not found with the open-source routine present in the package NNLSlab in \cite{SH}.

All the numerical Matlab software is available at {\cite{SH}}.

\section{Acknowledgements}

The work of A. Sommariva and M. Vianello is partially supported by the BIRD163015, PRAT-CPDA143275, DOR funds of the University of Padova, and by the GNCS-INdAM. Their research has been accomplished within the RITA {\sl{Research ITalian network on Approximation}}.
 
\bibliographystyle{plain}
\begin{thebibliography}{99}

 
\bibitem{B09} % CITED
 {\sc F.~Basaglia}, A new cubature method over polygons (Italian), Laurea Thesis in Mathematics, University of Padova, 2009 (advisors A. Sommariva and M. Vianello).

\bibitem{BX10} % CITED
 {\sc B.~Bauman and H.~Xiao}, Gaussian quadrature for optical design with noncircular pupils and fields, and broad wavelength range, Proc. SPIE, 7652(2) (2010), pp.1--12.
 
\bibitem{BSA} % CITED
 {\sc M.~Bern, J.R.~Shewchuk, and N.~Amenta}, Triangulations and mesh generation, Handbook of Discrete and Computational Geometry, Third Edition, CRC Press.

\bibitem{BE92} % CITED
 {\sc J.~Berntsen, T.O.~Espelid}, Algorithm 706: DCUTRI: an algorithm for adaptive cubature over a collection of triangles, ACM Transactions on Mathematical Software, 18 (1992), pp.329--342.

\bibitem{BG65} % CITED
 {\sc P.A.~Businger and G.H.~Golub}, Linear least-squares solutions by Householder
transformations, Numer. Math. 7 (1965), pp.269--276.

\bibitem{BP} % CITED
 {\sc M.G.~Blyth and C.~Pozrikidis}, A Lobatto integration grid over the triangle, IMA J. Appl. Math (2005), pp. 1--17.

\bibitem{C07} % CITED
 {\sc C.~Caratheodory}, {\"U}ber den Variabilittsbereich der Koeffizienten von Potenzreihen, die gegebene Werte nicht annehmen, Math. Ann. 64 (1907), pp.95--115.

\bibitem{CH91} % CITED
 {\sc B.~Chazelle}, Triangulating a Simple Polygon in Linear Time, Discrete and Computational Geometry, 6 (1991), pp.485--524.

\bibitem{CO97} % CITED
 {\sc R.~Cools}, Constructing cubature formulae: the science behind the art, Volume 6 January 1997, pp.1--54.
 
\bibitem{DFSV13} % CITED
 {\sc G.~Da Fies, A.~Sommariva and M..~Vianello}, Algebraic cubature by linear blending of elliptical arcs, Appl.Num.Math. 74 (2013), pp.49--61.
 
  \bibitem{DR84} % CITED
 {\sc E.~de Doncker and I.~Robinson}, An Algorithm for Automatic Integration over a Triangle Using Nonlinear Extrapolation, ACM Trans. Math. Soft., 10 (1984), pp.1--16.
 
 \bibitem{D85} % CITED
 {\sc D.A.~Dunavant}, High degree efficient symmetrical gaussian quadrature rules for triangles, Int. J. for Num. Meth. Eng., 21 (1985), pp.1129--1148.

%\bibitem{EGET93}
%{\sc H.~ElGindy, H.~Everett, and G.T.~Toussaint}, Slicing an ear using prune-and-search, Pattern Recognition Letters, (1993), pp.719--722.

\bibitem{FS12} % CITED
{\sc M.~Festa and A.~Sommariva}, Computing almost minimal formulas on the square, J. Comput. Appl. Math. 236 (2012), pp.4296--4302. 

\bibitem{F88} % CITED
{\sc G.W.~Forbes}, Optical system assessment for design: numerical ray tracing in the Gaussian pupil, J. Opt. Soc. Am. A 5 (1988), pp.1943--1956.

\bibitem{G11} % CITED
{\sc M.~Gentile, A.~Sommariva and M.~Vianello}, Polynomial interpolation and cubature over polygons, J. Comput. Appl. Math. 235 (2011), pp.5232--5239. 

\bibitem{Gsoft}
{\sc M.~Gentile, A.~Sommariva and M.~Vianello}, POLYGINT: a Matlab code for interpolation and cubature at Discrete Extremal Sets of polygons (available online at: {\tt{www.math.unipd.it/$\sim$marcov/CAAsoft.html}}).

\bibitem{G90} % CITED
 {\sc A.~Genz}, An adaptive numerical integration algorithm for simplices, Lecture Notes in Computer Science, 507, Springer-Verlag, (1991), pp.279--285.

\bibitem{KR87} % CITED
 {\sc D.K.~Kahaner, O.W.~Rechard}, TWODQD an adaptive routine for two-dimensional integration, J. Comput. Appl. Math. 17 (1987), pp.215--234. 
 
 
\bibitem{K92} % CITED
 {\sc P.~Keast}, Software for Integration Over Triangles and General Simplices, NATO ASI Series (Series C: Mathematical and Physical Sciences), vol 357 (1992), pp..283-294.
 
 \bibitem{L82} %CITED
 {\sc D.P.~Laurie}, CUBTRI: Automatic Integration over a Triangle, ACM Trans. Math. Soft., 8 (1982), pp.210--218.

 \bibitem{LH95} %CITED
 {\sc C.L.~Lawson and R.J.~Hanson}, Solving least squares problems. Revised reprint of the 1974 original, Classics in Applied Mathematics 15, SIAM, Philadelphia, 1995.


 \bibitem{LC95} %CITED
 {\sc J.N.~Lyness and R.~Cools}, A Survey of Numerical Cubature over Triangles, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill, 1994.

\bibitem{LJ75} % CITED
{\sc J.N.~Lyness, D.~Jespersen}, Moderate Degree Symmetric Quadrature Rules for the Triangle, J. Inst. Maths Applies, 15 (1975), pp.19-32.

\bibitem{MATpolyshape}
{\sc Mathworks}, Polyshape documentation, (available in Matlab R2017b and online at: {\tt{https://it.mathworks.com/help/matlab/ref/polyshape.html}}).

\bibitem{ME75} % CITED 
{\sc G.H.~Meisters}, Polygons have ears, American Mathematical Monthly, (1975), pp.648--651.

%\bibitem{MFLS17} % CITED
%{\sc J.M.~Myre, E.~Frahm, D.J.~Lilja, M.O.~Saar}, TNT-NN: A Fast Active Set Method for Solving %Large Non-Negative Least Squares Problems, Procedia Computer Science, 108 (2017), pp.%755--764. See {\tt{https://github.com/MyreResearchGroup/tnt-nn}} for the related software.

\bibitem{PP15} % CITED 
{\sc S-A.~Papanicolopulos}, Computation of moderate-degree fully symmetric cubature rules on the triangle using symmetric polynomials and algebraic solving, Comput. Math. Appl. 69 (2015), pp.650--666. 

\bibitem{PP16} % CITED 
{\sc S-A.~Papanicolopulos}, New fully symmetric and rotationally symmetryc cubature rules on the triangle using minimal orthonormal bases, J. Comput. Appl. Math. 294 (2016), pp.39--48. 

\bibitem{PSV17} % CITED 
{\sc F.~Piazzon, A. Sommariva and M. Vianello}, Caratheodory-Tchakaloff Subsampling, Dolomites Res. Notes Approx. DRNA 10 (2017), pp.5--14.

\bibitem{SH08} % CITED 
{\sc L.~Shampine}, Vectorized adaptive quadrature in Matlab, J. Comput. Appl. Math., 211 (2008), pp.131--140.

\bibitem{SH} % CITED 
{\sc M.~Slawski}, Non-negative least squares: comparison of algorithms, {\tt{https://sites.google.com/site/slawskimartin/code}}.

\bibitem{SH} % CITED
 {\sc A.~Sommariva}, {\tt{http://www.math.unipd.it/$\sim$alvise/software.html}}.
 
\bibitem{Ssoft} % CITED
 {\sc A.~Sommariva}, {\tt{http://www.math.unipd.it/$\sim$alvise/sets.html}}.


\bibitem{SV07} % CITED
 {\sc A.~Sommariva and M.~Vianello}, Product Gauss cubature over polygons based on Green's integration formula, BIT Numerical Mathematics 47 (2007), pp.441--453.

\bibitem{SV15} % CITED
 {\sc A.~Sommariva and M.~Vianello}, Compression of multivariate discrete measures and applications, Numer. Funct. Anal. Optim. 36 (2015), pp.1198--1223.
 
 \bibitem{SV18} % CITED
 {\sc A.~Sommariva and M.~Vianello}, Nearly optimal nested sensors location for polynomial regression on complex geometries, to appear on Sampling Theory in Signal and Image Processing.
 
 \bibitem{SU} % CITED
 {\sc S.~Suri}, Polygon Triangulation, (available online at: {\tt{https://www.cs.ucsb.edu/$\sim$suri/cs235/Triangulation.pdf}}).

\bibitem{TWB07} % CITED
 {\sc M.A.~Taylor, B.A.~Wingate, L.P.~Bos}, A cardinal function algorithm for computing multivariate quadrature points, SIAM Journal on Numerical Analysis, volume 45, no. 1 (2007), pp.193--205.

\bibitem{T57} % CITED
 {\sc V.~Tchakaloff}, Formules de cubatures m\'ecaniques \`a coefficients non n\'egatifs, (French) Bull. Sci. Math. 81 (1957), pp.123--134.
 
 \bibitem{T15} % CITED
 {\sc M. Tchernychova}, {\sl{Caratheodory}} cubature measures, Ph.D. dissertation in Mathematics (supervisor: T. Lyons), University of Oxford, 2015.

\bibitem{T09} % CITED
 {\sc R.T.~Tr{\^{i}}mbitas}, Adaptive Cubatures on Triangle, Results. Math.,, 53 (2009), pp.453?462,

\bibitem{XG10}
{\sc H.~Xiao and Z.~Gimbutas}, A numerical algorithm for the construction of efficient quadrature rules in two dimensions, Comput. Math. Appl. 59 (2010), pp.663--676.

\bibitem{XRY03} % CITED
{\sc H.~Xiao, V.~Rokhlin and N.~Yarvin}, Prolate SpheroidalWave Functions, Quadrature and Interpolation, Inverse Problems 17 (2001), pp.805--838.

\bibitem{WX03} % CITED
{\sc S.~Wandzura and H.~Xiao}, Symmetric Quadrature Rules on a Triangle, Comput. Math. Appl. 45 (2003), pp.1829--1840.




\end{thebibliography}

\end{document}



