text:
\subsection{Setup and Main Results}
\label{sec:Setup}

Let $L^2(\D,\mu)$ for $\D := \times_{j\in [D]} \D_j \subset \mathbbm{R}^D$ denote all functions $f: \D \rightarrow \mathbbm{C}$ that are square-integrable with respect to a given product of probability measures $\mu := \times_{j\in [D]} \mu_j$ over $\D$, and suppose that you are given a countable orthonormal (with respect to $\mu$) basis,
\begin{equation}
\B' := \left\{ T_{{\boldsymbol{n}}}: \D \rightarrow \C ~\big|~ {\boldsymbol{n}} \in \mathbbm{N}^D \right\},
\label{def:B'}
\end{equation}
of $L^2(\D,\mu)$ so that
\begin{equation*}
\langle T_{{\boldsymbol{k}}}, T_{{\boldsymbol{l}}}  \rangle_{(\D,\mu)} := \int_{\D} T_{{\boldsymbol{k}}}({\boldsymbol{\xi}}) \overline{T_{{\boldsymbol{l}}}({\boldsymbol{\xi}})} d\mu({\boldsymbol{\xi}}) = \delta_{{\boldsymbol{k}}, {\boldsymbol{l}}} =
\begin{cases}
1 & \textrm{if}~{\boldsymbol{k}} = {\boldsymbol{l}} \\
0 & \textrm{if}~{\boldsymbol{k}} \neq {\boldsymbol{l}}
\end{cases}.
\end{equation*}
Furthermore, suppose that $\B'$ is a tensor product basis so that
\begin{equation}
T_{{\boldsymbol{n}}}({\boldsymbol{\xi}}):=\prod _{j\in[D]} T_{j;n_j}(\xi_{j})
\label{def:T_n}
\end{equation}
holds for all ${\boldsymbol{n}} \in \mathbbm{N}^D$ and ${\boldsymbol{\xi}} \in \D$, where each set
$\B'_j := \left\{ T_{j;n_j}: \D_j \rightarrow \C ~\big|~ n_j \in \mathbbm{N} \right\}$
 with $j \in [D] := \{ 0, \dots, D-1 \}$ is itself an orthonormal (with respect to the probability measure $\mu_j$ over $\D_j \subset \mathbbm{R}$) basis of $L^2(\D_j,\mu_j)$.  We will call any such basis $\B'$ an Orthonormal Product Basis (OPB) with respect to $\mu$.

Our objective in this paper is to approximate smooth functions $f \in L^2(\D,\mu)$ as rapidly as possible using just a few point evaluations.  Toward this end we will take the traditional approach of considering only a finite subset $\B_{N,d}$ of $\B'$, and then approximating $f$ by approximating its projection $\tilde{f}$ onto the span of $\B_{N,d}$ (consider, e.g., hyperbolic cross/sparse grid methods for approximating functions of several variables \cite{shen2010sparse,dung2016hyperbolic,bungartz_griebel_2004}).  The potential improvement that the sublinear-time methods considered herein will then potentially provide over such standard methods in some cases will come from the fact that the finite basis $\B_{N,d}$ can be chosen to be extremely large herein (e.g., experiments were performed for Section~\ref{sec:Numerics} on a standard workstation using bases of cardinality $200^{100}$).  More specifically, herein we will consider two different types of bases $\B'$, each of which will allow us to demonstrate that the finite basis $\B_{N,d} \subset \B'$ we select below for approximation purposes also promote computational efficiency.

We will characterize $\B'$ below based on the behavior of its lowest order elements
\begin{equation*}
\B_N := \left\{ T_{{\boldsymbol{n}}} ~\big|~ \| {\boldsymbol{n}} \|_{\infty} < N \right\} \subset \B'
\end{equation*}
which we will assume throughout this paper is a finite Bounded Orthonormal System (BOS) with respect to the probability measure $\mu$ over $\D$ with a finite BOS constant
$$K' :=\max_{{\boldsymbol{n}} \in [N]^D} \|T_{{\boldsymbol{n}}}\|_\infty := \max_{{\boldsymbol{n}} \in [N]^D} {\rm sup}_{{\boldsymbol{\xi}} \in \D} \left| T_{{\boldsymbol{n}}} ({\boldsymbol{\xi}}) \right| \in [1,\infty).$$
Note that this implies that each set $\B_{j;N} := \left\{ T_{j;n_j}: \D_j \rightarrow \C ~\big|~ n_j \in [N] \right\} \subset \B'_j$ with $j \in [D]$ is itself also a BOS with respect to the probability measure $\mu_j$ over $\D_j \subset \mathbbm{R}$ with BOS constant
\begin{equation}
K_j :=\max_{n_j \in [N]} \|T_{j;n_j}\|_\infty \in [1,\infty).
\label{def:KjBOSC}
\end{equation}
Finally, we will further define $K^0_j$ to be
\begin{equation}
1 \leq K^0_j := \|T_{j;0}\|_\infty  \leq K_j
\label{def:K0jBOSC}
\end{equation}
for each $j \in [D]$.  Note that $K^0_j$ is strictly smaller than $K_j$ for many BOSs of interest (e.g., the cosine and Chebyshev polynomial bases as $\B'_j$ both have $K_j^0 = 1 < K_j = \sqrt{2}$).
From these two definitions we can also see, e.g., that $\prod _{j\in[D]}  {K}^0_j \leq K' = \prod _{j\in[D]}  {K}_j$ always holds.  Due to the boundedness of $K'$ assumed throughout the remainder of this paper we will always refer to $\B_N$ (as well as $\B'$ with slight abuse) as a Bounded Orthonormal Product Basis (BOPB) going forward.

We will approximate any given smooth $f \in L^2(\D,\mu)$ by approximating its projection $\tilde{f}$ onto the span of the finite BOS set
\begin{equation}
\B_{N,d} := \left \{ T_{{\boldsymbol{n}}}  ~\big|~ {\boldsymbol{n}} \in [N]^D ~{\rm and}~\| {\boldsymbol{n}} \|_0 \leq d \right \} \subseteq \B_N \subset \B'
\label{defBNd}
\end{equation}
for some $d \in [D+1]$, where $d$ is used to the help constrain the BOS constant.
The BOS constant $K \leq K'$ of $\B_{N,d}$ will be referred to as the {\it effective BOS constant} below.
As is usually the case in compressive sensing scenarios involving BOSs, its size will be a significant consideration with respect to sampling and computational efficiency.  In order to limit $K$'s size we will concentrate on the following two types of BOPBs going forward:
\begin{itemize}
\item \underline{\bf BOPBs of TYPE I:}  We will say a BOPB is of type I if the BOS constants $K_j$ are $1$ for all but at most $\tilde{d} \in \mathbbm{Z} \cap [0, D]$ BOS basis sets $\B_{j;N}$.  In this case we let $K_{\infty} := \max_{j\in [D]} K_j$ and note that $1 \leq K \leq K' \leq K_{\infty}^{\tilde{d}}$ so that $K$ will scale sub-exponentially in $D$ when $\tilde{d} \ll D$ independently of our choice of $d$ in \eqref{defBNd}.  We note that this type of BOPB includes several interesting examples of bases including the multidimensional Fourier basis (for which $\tilde{d} = 0$), and mixed BOPBs $\B'$ that have one-dimensional Fourier bases used for all but $\tilde{d}$ of their $\B'_j$ component bases.
\item \underline{\bf BOPBs of TYPE II:}  We will say a BOPB is of type II if $K_0 := \max_{j\in [D]} K^0_j = 1$.  This type of BOPB includes many bases where having a small number of interacting dimensions, $d$, helps to limit the effective BOS constant $K$ involved in the underlying sparse approximation problem.  Examples include the multivariate cosine, Chebyschev, and Legendre polynomial bases, as well as mixed polynomial bases where each one-dimensional component basis $\B'_j$ is, e.g., a potentially different Jacobi polynomial basis.
\end{itemize}

In either case above one can see that $K \leq K_{\infty}^d K^{D-d}_0$ will always hold.  In particular, $K = 2^{d/2}$ always holds if $\B'$ is either the multivariate cosine or Chebyshev basis in the type II case.  In the type I case we note that $1 \leq K \leq K_{\infty}^{\min(d,\tilde{d})} K^{\tilde{d}-\min(d,\tilde{d})}_0 \leq K_{\infty}^{\tilde{d}}$ will always hold so that $d$ can be set to $D$ without causing $K$ to become too large if, e.g., $\tilde{d} \ll D$.  This is certainly the case if $\B'$ is the multidimensional Fourier basis where $\tilde{d} = 0$.

Let $f \in L^2(\D,\mu)$ be smooth enough\footnote{Given that we will be recovering $f$ based on point samples we will require at least enough smoothness to guarantee that any particular point sample we might possibly utilize actually contains information about the given function's basis coefficients $\left\{ c_{{\boldsymbol{n}}} \right\}_{{\boldsymbol{n}} \in \mathbbm{N}^D}$.  Of course, the details regarding this smoothness requirement will vary with the choice of basis $\B'$.} that there exists a sequence $\left\{ c_{{\boldsymbol{n}}} \right\}_{{\boldsymbol{n}} \in \mathbbm{N}^D}$ such that
\begin{equation}
f({\boldsymbol{\xi}}) = \sum_{{\boldsymbol{n}} \in \mathbbm{N}^D} c_{{\boldsymbol{n}}} T_{{\boldsymbol{n}}}({\boldsymbol{\xi}})
\label{equ:FunctoApprox}
\end{equation}
holds pointwise for all ${\boldsymbol{\xi}} \in \D$.  Given such an $f: \D \rightarrow \mathbbm{C}$, we will denote its orthogonal projection onto the span of $\B_{N,d}$ by $\tilde{f}: \D \rightarrow \mathbbm{C}$.  Let
$$\mathcal{I}_{N,d} := \left\{ {\boldsymbol{n}} \in [N]^D ~\big|~ \| {\boldsymbol{n}} \|_0 \leq d \leq D \right\} \subset \mathbbm{N}^D$$
be the set of indices corresponding to the basis elements in $\B_{N,d}$.  We then have that
\begin{equation}
\tilde{f}({\boldsymbol{\xi}}) := \sum_{{\boldsymbol{n}} \in \mathcal{I}_{N,d}} \tilde{c}_{{\boldsymbol{n}}} T_{{\boldsymbol{n}}}({\boldsymbol{\xi}})
\label{equ:FinApproxf}
\end{equation}
for all ${\boldsymbol{\xi}} \in \D$, where ${\boldsymbol{\tilde{c}}}$ will be considered to be a vector in $\mathbbm{C}^{|\mathcal{I}_{N,d}|}$ indexed by $\mathcal{I}_{N,d}$.  Note further that the entries of ${\boldsymbol{\tilde{c}}}$ will satisfy $\tilde{c}_{{\boldsymbol{n}}} = c_{{\boldsymbol{n}}}$ for all ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$.

As mentioned above, we will ultimately approximate $f$ by producing a sparse approximation in $\B_{N,d}$ to $\tilde{f}$.  The best possible $s$-term approximation to $\tilde{f}$ in $\B_{N,d}$ will be denoted by $\tilde{f}_s^{\rm opt}: \D \rightarrow \mathbbm{C}$, and will be defined as follows:  Order the basis coefficients ${\boldsymbol{\tilde{c}}} \in \mathbbm{C}^{|\mathcal{I}_{N,d}|}$ of $\tilde{f}$ by their magnitudes so that
$$\left| \tilde{c}_{{\boldsymbol{n}}_1} \right| \geq \left| \tilde{c}_{{\boldsymbol{n}}_2} \right| \geq \left| \tilde{c}_{{\boldsymbol{n}}_3} \right| \geq \dots \geq  \left| \tilde{c}_{{\boldsymbol{n}}_{|\mathcal{I}_{N,d}|}} \right|,$$
where ties are broken lexicographically using the entries' indices in $\mathcal{I}_{N,d}$.  Then $\tilde{f}_s^{\rm opt}$ will be defined to be
$$\tilde{f}_s^{\rm opt} ({\boldsymbol{\xi}}) := \sum_{j = 1}^s \tilde{c}_{{\boldsymbol{n}}_j} T_{{\boldsymbol{n}}_j}({\boldsymbol{\xi}})$$
for all ${\boldsymbol{\xi}} \in \D$, and its (potentially) nonzero coefficients' indices will be denoted by
$${\Omega^{\rm opt}_{\tilde{f},s}} := \left\{ {\boldsymbol{n}}_1, \dots, {\boldsymbol{n}}_s \right\} \subset \mathcal{I}_{N,d}.$$

Note that $\tilde{f}_s^{\rm opt}$ will indeed have the property that
$$\left\|  \tilde{f} - \tilde{f}_s^{\rm opt} \right\|_{{L^2(\mathcal{D},\mu)}}~=~ \inf_{\small {\boldsymbol{z}} \in \mathbbm{C}^{|\mathcal{I}_{N,d}|},~\| {\boldsymbol{z}} \|_0 \leq s} \left\| \tilde{f} - \sum_{{\boldsymbol{n}} \in \mathcal{I}_{N,d}} z_{{\boldsymbol{n}}} T_{{\boldsymbol{n}}} \right\|_{{L^2(\mathcal{D},\mu)}}.$$
Furthermore, if we let ${\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} \in \mathbbm{C}^{|\mathcal{I}_{N,d}|}$ denote the $\B_{N,d}$ basis coefficients of $\tilde{f}_s^{\rm opt}$ then we can see that both $\left\| {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} \right\|_0 \leq s$ and $\left\|  {\boldsymbol{\tilde{c}}} - {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} \right\|_{2} ~=~\left\|  \tilde{f} - \tilde{f}_s^{\rm opt} \right\|_{{L^2(\mathcal{D},\mu)}}$ will hold.  As a result, norms involving the vector ${\boldsymbol{\tilde{c}}} - {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} \in \mathbbm{C}^{|\mathcal{I}_{N,d}|}$ can be interpreted as best $s$-term approximation errors of $\tilde{f}$ in a natural way.

Finally, to prove our main result below we will effectively be considering the point samples we take from $f$ to instead be point samples taken from $\tilde{f}$ that are contaminated with evaluation errors of size $\left(f - \tilde{f} \right) ({\boldsymbol{\xi}})$ at each evaluation point ${\boldsymbol{\xi}} \in \mathcal{G}$.  To bound all of these errors in a uniform fashion we will define
\begin{equation}
\gamma := \left\| f - \tilde{f} \right\|_\infty ~=~ {\rm sup}_{{\boldsymbol{\xi}} \in \D} \left|  \left(f - \tilde{f} \right) ({\boldsymbol{\xi}}) \right|.
\label{equ:gammaDef}
\end{equation}
The following theorem is proven in Section~\ref{sec:CoSaMPguarantee}.\\

\begin{mytheorem}{(Main Result).}
Let $\eta \in (0, \infty)$ and $s,d,N \in \mathbbm{N} \setminus \{1 \}$ with $d \leq D$ and $s < |\mathcal{I}_{N,d}| / 2$.  There exists a finite set of grid points $\mathcal{G} \subset \D$, an algorithm $\mathcal{H}: \mathbbm{C}^{\left| \mathcal{G} \right|} \rightarrow \left( \mathcal{I}_{N,d} \times \mathbbm{C} \right)^s$, and an absolute universal constant $C' \in \mathbbm{R}^+$ such that the function $a: \D \rightarrow \mathbbm{C}$ defined by $a({\boldsymbol{\xi}}) := \sum_{({\boldsymbol{n}},a_{{\boldsymbol{n}}}) \in \mathcal{H}(f(\mathcal{G}))} a_{\boldsymbol{n}} T_{\boldsymbol{n}}({\boldsymbol{\xi}})$ satisfies
\begin{equation}
\| f - a \|_{L^2(\mathcal{D},\mu)} \leq \left\| f - \tilde{f} \right\|_{L^2(\mathcal{D},\mu)} + C' \left( \sqrt{s} \left\| {\boldsymbol{\tilde{c}}} - {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} \right\|_2 + \left\| {\boldsymbol{\tilde{c}}} - {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}}  \right\|_1 + \gamma \sqrt{s} \right) + \eta
\label{eqn:MRintroerrorg}
\end{equation}
for all $f = \sum_{{\boldsymbol{n}} \in \mathbbm{N}^D} c_{{\boldsymbol{n}}} T_{{\boldsymbol{n}}} \in L^2(\mathcal{D},\mu)$ with $\gamma := \left\| f - \tilde{f} \right\|_\infty$ $=~ {\rm sup}_{{\boldsymbol{\xi}} \in \D} \left|  \left(f - \tilde{f} \right) ({\boldsymbol{\xi}}) \right| ~<~ \infty$, where $\tilde{f}: \D \rightarrow \mathbbm{C}$ is the finite dimensional approximation to $f$ defined as per \eqref{equ:FinApproxf}.

	If the BOPB $\B_{N,d}$ is of type I so that the BOS constants $K_j$ are $1$ for all but at most $\tilde{d} \in \mathbbm{Z} \cap [0, D]$ BOS basis sets $\B_{j;N}$, then
	$$\left| \mathcal{G} \right| = \mathcal{O} \left( s^3 D K^{4\tilde{d}}_\infty d^4 \cdot  \log^4 \left( \frac{DN}{d} \right) \log^2 (s) \log^2 (D)  \right),$$
	and the algorithm $\mathcal{H}$ will have runtime complexity
	$$\mathcal{O}\left( \left( s^5 + s^3 N \right) D^2 K^{4\tilde{d}}_\infty d^4 \cdot  \log^4 \left( \frac{DN}{d} \right) \log^2 (s) \log^2 (D) \log \left( \left\| {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} \right \|_2 / \eta\right) \right).$$
	If the BOPB $\B_{N,d}$ is of type II so that $K_0 = 1$, then
	$$\left| \mathcal{G} \right| = \mathcal{O} \left( s^3 D K^{4d}_\infty d^4 \cdot  \log^4 \left( \frac{DN}{d} \right) \log^2 (s) \log^2 (D) \right),$$
	and the algorithm $\mathcal{H}$ will have runtime complexity
	$$\mathcal{O} \left( \left( s^5 + s^3 N \right) D^2 K^{4d}_\infty d^4 \cdot  \log^4 \left( \frac{DN}{d} \right) \log^2 (s) \log^2 (D) \log \left( \left\| {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} \right \|_2 / \eta\right) \right).$$
	Here we have assumed that any desired basis function $T_{{\boldsymbol{n}}} \in \B_{N,d}$ can be evaluated at any desired point in $\D$ in $\mathcal{O}(N D)$-time (which will be the case, e.g., for polynomial product bases of degree $\leq ND$).
	\label{thm:MRintro}
\end{mytheorem}

\begin{proof} This is a restatement of Corollary~\ref{cor:MainRes}.
\end{proof}

Looking at Theorem~\ref{thm:MRintro} we can see that it effectively subsumes the theoretical recovery results of \cite{choi2018sparse}.  Consider, for example, the case where $f = \tilde{f} = \tilde{f}_s^{\rm opt}$ (so that $f$ is exactly $s$-sparse in $\B_{N,d}$).  In this setting we will have both $\left\| f - \tilde{f} \right\|_{L^2(\mathcal{D},\mu)} = \gamma = 0$ and ${\boldsymbol{\tilde{c}}} - {\boldsymbol{\tilde{c}}}_{\Omega^{\rm opt}_{\tilde{f},s}} = {\bf 0}$ hold true so that \eqref{eqn:MRintroerrorg} implies that $f$ is recovered exactly (up to any chosen tolerance $\eta$).  Unlike the results in \cite{choi2018sparse}, however, Theorem~\ref{thm:MRintro} also guarantees that the method $\mathcal{H}$ will work well for functions $f$ which have $\| f - \tilde{f}_s^{\rm opt} \|_{L^2(\mathcal{D},\mu)}$ relatively small, but nonzero.

The authors would also like to emphasize the generality of Theorem~\ref{thm:MRintro}, which is unique to the best of their knowledge in the literature related to sublinear-time sparse approximation methods.  If, for example, one chooses $\B_{N,d}$ to be the multidimensional Fourier basis with $d = D$ and $\tilde{d} = 0$ then one immediately obtains a new sparse Fourier transform result for functions of many variables whose sampling and runtime requirements scale only polylogarithmically in the total basis size $|\B_{N,d}| = N^D$.  Though this new Fourier result does not compare favorably to the best deterministic multidimensional Fourier results of this kind \cite{iwen2013improved,morotti2017explicit} with respect to achievable runtimes or error guarantees, it also does not use any of the specific algebraic structure of the Fourier basis.  This allows Theorem~\ref{thm:MRintro} to be significantly more flexible than these older Fourier results in that it can apply to situations which they don't cover.  For example, it can generate entirely discrete Fourier results where $\mathcal{G} \subset \{ \frac{j}{N}~\big|~ j \in [N] \}^D$ (unlike \cite{iwen2013improved}) by using a discrete and finite multidimensional Fourier BOBP with $f = \tilde{f}$ which will work for any choice of $N$ (unlike \cite{morotti2017explicit}, which requires $N$ to be prime).

Finally, the astute reader has likely noticed that Theorem~\ref{thm:MRintro} is phrased in the form of an existence result, which may be troubling to the practical numerical analyst who actually wants to know how to compute an accurate solution.  Let us allay any anxieties that this choice of theoretical statement may have birthed -- the algorithm $\mathcal{H}$ referred to above is a modified version of the well known CoSaMP algorithm \cite{needell2009cosamp} (see Algorithm~\ref{alg:main} in Section~\ref{sec:CoSaMPguarantee}).  It has been implemented and evaluated in Section~\ref{sec:Numerics}, and the code made publicly available.\footnote{See ``SHT II: Best s-Term Approximation Guarantees for Bounded Orthonormal Product Bases in Sublinear-Time'' on Mark Iwen's code page \url{https://www.math.msu.edu/~markiwen/Code.html}.}  In short, the result is entirely explicit and constructive with respect to the algorithm $\mathcal{H}$.  The grid $\mathcal{G} \subset \D$, which ultimately responsible for the form of Theorem~\ref{thm:MRintro} as an existence result, on the other hand, is a bit more nuanced with respect to its practical construction.

As we shall see below, the grid $\mathcal{G} \subset \D$ is constructed by randomly selecting points from $\D$ according to several prescribed probability distributions that are ultimately derived from the orthogonality measure $\mu$ (see Theorems~\ref{thm:iterInvariant}~and~\ref{thm:SuppIDWorks} and their proofs for details).  It is then proven that this randomly constructed grid will allow $\mathcal{H}$ to satisfy the error guarantee \eqref{eqn:MRintroerrorg} for all functions $f$ as per \eqref{equ:FunctoApprox} with high probability while simultaneously satisfying the stated upper bounds on its cardinality.  The runtime complexity of $\mathcal{H}$ follows from the boundedness of $|\mathcal{G}|$.  Hence, the existence result is proven by randomly constructing a grid $\mathcal{G}$ which is guaranteed to satisfy the conclusions of Theorem~\ref{thm:MRintro} with high probability.

In fact, this is completely analogous to the role of random sampling matrices in standard compressive sensing results involving the Restricted Isometry Property (RIP).  Many compressive sensing methods are guaranteed to be accurate if they are used in combination with a random sampling matrix that has the RIP, a condition which can only be achieved near-optimally with high probability.  Herein, the conclusions of Theorem~\ref{thm:MRintro} will hold for any grid $\mathcal{G}$ that can be used to form two associated random sampling matrices: one with the RIP, and another with a property known as the {\it Support Identification Property (SIP)} which will be defined in Section~\ref{subsec:SIP}.  The conclusions of Theorem~\ref{thm:MRintro} will hold whenever these two conditions are satisfied by $\mathcal{G}$, and it will be shown that a randomly constructed grid $\mathcal{G}$ will satisfy both conditions with high probability.

\subsection{An Outline of the Paper and of the Proof of Theorem~\ref{thm:MRintro}}

After reviewing some relevant compressive sensing results and establishing necessary notation in Section~\ref{sec:Notation}, we will begin proving Theorem~\ref{thm:MRintro} in Section~\ref{sec:CoSaMPguarantee}.  The first step in that process will be to prove a compressive sensing recovery guarantee for a generalized version of the well known CoSaMP method \cite{needell2009cosamp}.  This new theorem, Theorem~\ref{thm:iterInvariant}, will establish a best $s$-term recovery guarantee for the CoSaMP algorithm where the support identification step is performed by any algorithm $\mathcal{A}$ and grid $\mathcal{G}$ pair which has the SIP (see Section~\ref{subsec:SIP} below for details on the SIP).  With Theorem~\ref{thm:iterInvariant} in hand we will then turn our attention to constructing a sublinear-time algorithm $\mathcal{A}$ and grid $\mathcal{G}$ pair that have the SIP, an effort whose results are summarized by Theorem~\ref{thm:suppIdforsec3} (see also Proposition~\ref{coro:SIPrevealed}).  Combining Theorems~\ref{thm:iterInvariant}~and~\ref{thm:suppIdforsec3} then quickly establishes our main result above which appears in the form of Theorem~\ref{thm:NewCoSaMP} and Corollary~\ref{cor:MainRes} in Section~\ref{sec:CoSaMPguarantee}.

The vast majority of the effort in the paper will be focussed on proving Theorem~\ref{thm:suppIdforsec3} in Section~\ref{sec:SupportID}.  That is, to demonstrate that Algorithm~\ref{alg:suppid:impl} therein can be used together with a randomly constructed grid $\mathcal{G}$ in order to effectively achieve the SIP with high probability.  This is done by Theorem~\ref{thm:SuppIDWorks} (a specialized version of Theorem~\ref{thm:GenAlgSuppIDWorks} in Section~\ref{sec:SuppIDHeavyEls}) which formalizes the random sampling strategy one must use in order to construct $\mathcal{G}$ so that the SIP is achieved with high probability, and by Theorem~\ref{thm:suppId} which translates the conclusions of Theorem~\ref{thm:SuppIDWorks} into a SIP-type statement.  Theorem~\ref{thm:GenAlgSuppIDWorks}, in turn, follows from Theorem~\ref{thm:ExistFastFuncs} which is proven in Section~\ref{sec:proofBigSuppIDThm}.

Finally, the authors would like to note that the reader who is interested in seeing the proof of Theorem~\ref{thm:MRintro} unfold from basic compressive sensing principals in a more direct fashion (though without the benefit of waypoints explaining the relevance of each result to the final goal) might consider the following alternate reading order for the sections below:  Such readers can begin with Section~\ref{sec:proofBigSuppIDThm} after reviewing Section~\ref{sec:Notation}, followed by the first 4 paragraphs of Section~\ref{sec:SupportID}, then Section~\ref{sec:SuppIDHeavyEls}, and finally the remainder of Section~\ref{sec:SupportID} after which Theorem~\ref{thm:suppIdforsec3} will have been proven.  Reading Section~\ref{sec:SupportID} in this bottom up fashion first will then allow Section~\ref{sec:CoSaMPguarantee} to be read without having to temporarily take any of the theoretical statements therein for granted along the way.  For readers who are mostly interested in the numerical ramifications of the methods developed herein, we suggest skipping down to conduct a careful review of Algorithms~\ref{alg:main}~and~\ref{alg:suppid:impl} (together with the equations referred to therein) after reading Section~\ref{sec:Notation}, after which the careful numerical evaluation conducted in Section~\ref{sec:Numerics} should be understandable.

Before moving on to establish some additional required notation, however, we will first discuss the SIP in the next subsection.  This is crucial as the notion of the SIP will allow for easier sublinear-time methods to be developed in the future.  To emphasize this last point:  Any basis for which the SIP can be established via a sublinear-time algorithm can be combined with Theorem~\ref{thm:iterInvariant} below in order to produce a new sublinear-time compressive sensing method for that basis.  We expect that this new pathway for developing future sublinear-time algorithms will help to stimulate the further improvement and generalization of sparse Fourier transform techniques to other bases of interest going forward.

\subsection{The Support Identification Property (SIP)}
\label{subsec:SIP}

As above, let $[N] := \{ 0, \dots, N-1 \}$ for all $N \in \N$ and further define $\mathcal{P}([N])$ to be the power set of any such set $[N]$.  In Section~\ref{sec:CoSaMPguarantee} we will prove that CoSaMP will still produce accurate sparse approximations as long as its support identification step employs a triple with the {\it support identification property}.\\

\begin{mydef}[The Support Identification Property (SIP)]
Let $s \in [N]$, $\beta \in (0,1)$, $\Phi \in \mathbbm{C}^{m \times N}$, $\mathcal{A}: \mathbbm{C}^m \rightarrow \mathcal{P}([N])$, and $\Gamma: \mathbbm{C}^m \rightarrow [0,\infty)$ with $\Gamma({\boldsymbol{0}}) = 0$.  The triple $(\Phi,\mathcal{A},\Gamma)$ is said to have the Support Identification Property (SIP) of order $(s,\beta )$ if
$$\left\| {\boldsymbol{v}}_{\mathcal{A}(\Phi {\boldsymbol{v}} + {\boldsymbol{e}})^c} \right\|_2 \leq \beta \| {\boldsymbol{v}} \|_2$$
holds for all ${\boldsymbol{e}} \in \mathbbm{C}^m$ and ${\boldsymbol{v}} \in \mathbbm{C}^N$ with $\| {\boldsymbol{v}} \|_0 \leq s$ that also satisfy $\| {\boldsymbol{v}} \|_2 > \Gamma({\boldsymbol{e}})$.
\label{def:SIP}
\end{mydef}

Note that many triples with the SIP exist.  One completely trivial example is the triple consisting of the $N \times N$ identify matrix $I$, the function $\mathcal{A}$ which always outputs $[N]$, and the zero function $\Gamma$.  Of course this example is extremely unsatisfying -- generally for compressive sensing applications we prefer that the any SIP triple $(\Phi \in \mathbbm{C}^{m \times N},\mathcal{A}: \mathbbm{C}^m \rightarrow \mathcal{P}([N]),\Gamma)$ has $m \ll N$ and an efficient computational complexity for $\mathcal{A}$ (preferably sublinear-in-$N$ herein).  Thankfully these types of SIP triples also exist -- in fact it is easy to see that any fast and error-robust compressive sensing algorithm $\mathcal{A}$ must in fact be a member of such a triple.\\

\begin{lemma}
Let $\Gamma': \mathbbm{C}^m \rightarrow \mathbbm{R}^+$ be such that $\Gamma'({\bf 0}) = 0$, and let $\mathcal{A}: \mathbbm{C}^m \rightarrow \mathbbm{C}^N$ be a compressive sensing algorithm with an associated measurement matrix $\Phi \in \mathbbm{C}^{m \times N}$ that satisfies
$$\| {\boldsymbol{x}} - \mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right) \|_2 \leq \Gamma'({\boldsymbol{e}})$$
for all ${\boldsymbol{e}} \in \mathbbm{C}^m$ and ${\boldsymbol{x}} \in \mathbbm{C}^N$ with $\| {\boldsymbol{x}} \|_0 \leq s$.  Furthermore, let ${\rm supp}: \mathbbm{C}^N \rightarrow \mathcal{P}([N])$ output the indices of the nonzero entries of any given input vector.
Then, the triple $(\Phi,{\rm supp} \circ \mathcal{A},(1/ \beta') \Gamma')$ will have the SIP of order $(s,\beta')$ for all $\beta' < 1$.
\label{lem:EveryCSisSIP}
\end{lemma}

\begin{proof}
Let ${\boldsymbol{e}} \in \mathbbm{C}^m$ and note that
\begin{align*}
\left\|  {\boldsymbol{x}}_{{\rm supp}\left( \mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right) \right)^c} \right\|^2_2 ~&\leq~ \left\| {\boldsymbol{x}}_{{\rm supp}\left( \mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right) \right)^c} \right\|^2_2 + \left\| \left({\boldsymbol{x}} - \mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right) \right)_{{\rm supp}\left( \mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right) \right)}  \right\|^2_2\\
~&=~ \left\|  {\boldsymbol{x}} - \mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right) \right\|^2_2 ~\leq~ \left( \Gamma'({\boldsymbol{e}}) \right)^2.
\end{align*}
Thus, if $\| {\boldsymbol{x}} \|_2 > (1/\beta') \Gamma'({\boldsymbol{e}})$ then $\left\|  {\boldsymbol{x}}_{{\rm supp}\left( \mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right) \right)^c} \right\|_2 \leq \beta' \left((1/\beta')\Gamma'({\boldsymbol{e}}) \right) \leq \beta' \| {\boldsymbol{x}} \|_2.$
\end{proof}

Lemma~\ref{lem:EveryCSisSIP} demonstrates that many nontrivial SIP triples of the type we are interested in exist.  Of course, using a compressive sensing method in order to create a SIP triple seems slightly nonsensical given that one would generally want to create a SIP triple in order to develop a new compressive sensing method in the first place.  This immediately raises the question of whether nontrivial SIP triples exist which do not in themselves already effectively serve as a compressive sensing method.  The answer to that question is ``yes'', and the easiest example is the SIP triple which the original CoSaMP method is itself already implicitly utilizes.  Given ${\boldsymbol{x}} \in \mathbbm{C}^N$ and $s \in  [N]$ let ${\boldsymbol{x}} \big|_{s} \in \mathbbm{C}^N$ be the vector obtained from ${\boldsymbol{x}}$ by setting all but its $s$-largest magnitude entries to $0$.  The following lemma explicitly demonstrates the SIP triple on which the original CoSaMP algorithm \cite{needell2009cosamp} is implicitly based.\\

\begin{lemma}[The CoSaMP SIP Triple]
Let $\Phi \in \mathbbm{C}^{m \times N}$ have the RIP of order $(2s, 0.1)$ (so that its Restricted Isometry Constants (RIC)s satisfy $\delta_{s} \leq \delta_{2s} \leq 0.1$), $s \in [N]$, $\beta \in (0.2223, 1)$, and define $\mathcal{A}: \mathbbm{C}^m \rightarrow \mathcal{P}([N])$ by $\mathcal{A}({\boldsymbol{y}}) := {\rm supp} \left( (\Phi^* {\boldsymbol{y}})\big|_{s}~ \right)$, and $\Gamma: \mathbbm{C}^m \rightarrow [0,\infty)$ by $\Gamma({\boldsymbol{e}}) := \left( \frac{2.34}{\beta - 0.2223} \right) \| {\boldsymbol{e}} \|_2$.  Then, the triple $\left( \Phi,\mathcal{A}, \Gamma \right)$ has the SIP of order $(s,\beta )$.
\end{lemma}

\begin{proof}
Lemma 4.2 of \cite{needell2009cosamp} implies that $\left\|  {\boldsymbol{x}}_{\mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right)^c} \right\|_2 \leq 0.2223 \| {\boldsymbol{x}} \|_2 + 2.34 \| {\boldsymbol{e}} \|_2$ holds for all ${\boldsymbol{e}} \in \mathbbm{C}^m$ and ${\boldsymbol{x}} \in \mathbbm{C}^N$ with $\| {\boldsymbol{x}} \|_0 \leq s$.  Suppose, furthermore, that $\| {\boldsymbol{x}} \|_2 > \Gamma({\boldsymbol{e}}) = \left( \frac{2.34}{\beta - 0.2223} \right) \| {\boldsymbol{e}} \|_2$. Then,
\begin{align*}
\left\|  {\boldsymbol{x}}_{\mathcal{A} \left( \Phi {\boldsymbol{x}} + {\boldsymbol{e}} \right)^c} \right\|_2 &~\leq~0.2223 \| {\boldsymbol{x}} \|_2 + 2.34 \| {\boldsymbol{e}} \|_2 ~<~ 0.2223 \| {\boldsymbol{x}} \|_2 + 2.34 \left( \frac{\beta - 0.2223}{2.34} \right) \| {\boldsymbol{x}} \|_2\\
&~= \beta \| {\boldsymbol{x}} \|_2
\end{align*}
holds.  Also, $\Gamma({\bf 0}) = 0$.
\end{proof}

In Section~\ref{sec:CoSaMPguarantee} we will demonstrate that the original SIP triple implicitly used by the  CoSaMP algorithm can be replaced with any other SIP triple of similar quality without substantively changing the performance of the resulting CoSaMP variant as a compressive sensing algorithm.  Before we can do this, however, we will require some additional notation and preliminary infrastructural results.

\section{Introduction} \label{sec:Introduction}
\setcounter{equation}{0}
\input{Introduction.tex}

\section{Notation and Preliminaries} \label{sec:Notation}
\input{Notation.tex}

\section{Robust Sublinear-Time Sparse Approximation via CoSaMP with Fast Support Identification}
\label{sec:CoSaMPguarantee}
\input{CoSaMP_w_Supp_ID.tex}

\section{Sublinear-Time Support Identification}
\label{sec:SupportID}
\input{Support_ID.tex}

\section{Empirical Evaluation}
\label{sec:Numerics}
\input{Numerics.tex}

\section*{Acknowledgements}
Mark Iwen was supported in part by NSF DMS-1912706, and would like to dedicate this paper to his ever bright, hard working, and spirited wife Tsveta, and to the prosperity of their newborn daughter Evgenia.  Evgenia -- I am anxious to know you are healthy, eager to see you are happy, and already sad at the distant prospect your moving out.  May you be more like your mother than like me for your own sake!

\bibliographystyle{abbrv}

\subsection{Restrictions and Partial Evaluations}

The following notation will be utilized heavily during the analysis of the proposed support identification procedure.  Let $\S \subset [D]$, ${\boldsymbol{w}} \in \times_{j\in\S} \D_j$ with $w_j \in \D_j$, and ${\boldsymbol{n}} \in [N]^D$.  The function $T_{\S;{\boldsymbol{n}}}: \times_{j\in\S} \D_j \rightarrow \C$ is defined to be
\begin{equation}
T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}) := \prod_{j \in \S} T_{j;n_j} (w_j).
\label{equ:RestictedBOS}
\end{equation}
Then, the set
\begin{equation}
\B_{\S} := \left\{ T_{\S;{\boldsymbol{n}}} ~\big|~{\boldsymbol{n}} \in [N]^D ~{\rm with}~\| {\boldsymbol{n}} \|_0 \leq d \right\}
\label{def:B_S}
\end{equation}
is a BOS with respect to the probability measure $\mu_{\S} := \otimes_{j \in \S} \mu_j$ over $\displaystyle \D_{\S} := \times_{j\in\S} \D_j \subset \mathbbm{R}^{|\S|}$ with BOS constant $K_{\S} \leq \min \left\{ \prod_{j \in \S} K_j, K_{\infty}^d K^{\max\{|\S|-d,0\}}_0 \right\}$.
For any set $E$, let $\mathcal{P}(E)$ denote the power set of $E$ containing all possible subsets of $E$.
Given any vector ${\boldsymbol{v}} \in \C^p$ and $\T \subset [p]$ we will let ${\boldsymbol{v}}_{\T} \in \C^p$ have entries
\begin{equation*}
\left(v_{\T} \right)_j =
\begin{cases}
v_j & \textrm{if}~j \in \T\\
0 & \textrm{if}~j \notin \T
\end{cases}.
\end{equation*}
For $t \in [p]$, we let ${\boldsymbol{v}}_t \in \C^p$ a vector restricting ${\boldsymbol{v}}$ to its $t$ largest-magnitude entries.
Let $\S^c := [D] \setminus \S$ for all $\S \subset [D]$.  We will then construct $f_{\S;{\boldsymbol{w}}}:  \D_{\S^c} \rightarrow \C$ from $f: \D \rightarrow \C$ by defining
\begin{equation*}
f_{\S;{\boldsymbol{w}}}({\boldsymbol{z}}) := f({\boldsymbol{\xi}})
\end{equation*}
where ${\boldsymbol{\xi}} \in \D \subset \mathbbm{R}^D$ is the unique vector with ${\boldsymbol{\xi}}_{\S} = {\boldsymbol{w}}$ and ${\boldsymbol{\xi}}_{\S^c} = {\boldsymbol{z}}$.
In this context, we define the permutation function
$\varrho_\S\colon \D_{\S}\times\D_{\S^c}\rightarrow\D$ given by
\begin{equation}\label{def:rho}
\varrho_\S({\boldsymbol{w}}, {\boldsymbol{z}})={\boldsymbol{\xi}} \text{ such that } {\boldsymbol{\xi}}_\S={\boldsymbol{w}} \text{ and } {\boldsymbol{\xi}}_{\S^c}={\boldsymbol{z}}.
\end{equation}
This yields the alternative characterization $f_{\S;{\boldsymbol{w}}}({\boldsymbol{z}}) = f\left(\varrho_\S({\boldsymbol{w}}, {\boldsymbol{z}})\right)$.

The restricted vectors of the input vector ${\boldsymbol{\xi}}$ such as ${\boldsymbol{\xi}}_\S$ and ${\boldsymbol{\xi}}_{\S^c}$ have the reduced dimensions. However, the coefficient vectors such as ${\boldsymbol{\tilde{c}}}$ and ${\boldsymbol{\tilde{r}}}$ will maintain the full dimension even though they are restricted to some subset of indices.

If ${\boldsymbol{n}}, {\boldsymbol{m}} \in \mathcal{I}_{N,d}$ and $\S \subset [D]$ then we will define $({\boldsymbol{n}},{\boldsymbol{m}})_\S \in \mathcal{I}_{N,d}$ to be the vector ${\boldsymbol{n}}_\S + {\boldsymbol{m}}_{\S^c}$.  Furthermore, for a given $\vect{v}\in \C^{|\mathcal{I}_{N,d}|}$, ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$, and $\S \subset [D]$, we will let the vector $\vect{v}_{\S;{\boldsymbol{n}}} \in \C^{|\mathcal{I}_{N,d}|}$ indexed by $\vect{k} \in \mathcal{I}_{N,d}$ have entries given by
\begin{align}
\left({v}_{\S;{\boldsymbol{n}}}\right)_{\vect{k}}=
   \begin{cases}
       v_{\vect{k}},&\text{ if } {\boldsymbol{k}}_{\S} = {\boldsymbol{n}}_{\S} \\
       0 &\text{ otherwise}\\
   \end{cases}. \label{def:vecprefixfixed}
\end{align}
Note that $\vect{v}_{\S;{\boldsymbol{n}}} \in \C^{|\mathcal{I}_{N,d}|}$ will only have at most $N^{|\S^c|}$ nonzero entries corresponding to the entries of ${\boldsymbol{v}} \in \C^{|\mathcal{I}_{N,d}|}$, $v_{{\boldsymbol{m}}} \in \C$, whose indices ${\boldsymbol{m}}$ match those of ${\boldsymbol{n}}$ on $\S$ (i.e., so that ${\boldsymbol{m}}_{\S} = {\boldsymbol{n}}_{\S}$).

The following calculation will be repeated sufficiently often that it merits being referred to as a lemma.  It concerns the partial sum approximation to $f$ from \eqref{equ:FunctoApprox} in $\B$ given by \eqref{equ:FinApproxf}.
Recall that ${\boldsymbol{\tilde{c}}}$ contains only the entries of the sequence ${\boldsymbol{c}}$ corresponding to the indices in $\mathcal{I}_{N,d}$.

Moreover, the next lemma also demonstrates the usage of the newly introduced notation. Its statement will be used later in the proofs of Lemmas ~\ref{lem:EntryID_ExactInnerProd} and~\ref{lem:SupportLem2}.
\begin{lemma}
Let $\S \subset [D]$, ${\boldsymbol{w}} \in \D_{\S} = \times_{j\in\S} \D_j$ with $w_j \in \D_j$, and ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$.  Then
$$\left \langle \tilde{f}_{\S;{\boldsymbol{w}}}, T_{\S^c;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S^c}, \mu_{\S^c}\right)} = \left \langle~ {\boldsymbol{\tilde{c}}}_{\S^c;{\boldsymbol{n}}},  \overline{\Phi_{\S;{\boldsymbol{n}};{\boldsymbol{w}}} } ~\right \rangle$$
where $\tilde{f}$ is as in \eqref{equ:FinApproxf}, and $\Phi_{\S;{\boldsymbol{n}};{\boldsymbol{w}}} \in \C^{|\mathcal{I}_{N,d}|}$ is a vector indexed by $\vect{k} \in \mathcal{I}_{N,d}$ with entries
\begin{equation}
\left(\Phi_{\S;{\boldsymbol{n}};{\boldsymbol{w}}}\right)_{\boldsymbol{k}} :=
\begin{cases}
T_{\S;{\boldsymbol{k}}}({\boldsymbol{w}}) & \textrm{if}~{\boldsymbol{k}}_{\S^c} = {\boldsymbol{n}}_{\S^c}\\
0 & \textrm{otherwise}
\end{cases}.
\label{equ:RowRIPlem1}
\end{equation}
\label{lem:PartEvalinnerP}
\end{lemma}

\begin{proof}
Computing the inner product one quickly sees that
\begin{align*}
\left \langle \tilde{f}_{\S;{\boldsymbol{w}}}, T_{\S^c;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S^c}, \mu_{\S^c}\right)} &= \int_{\D_{\S^c}} \tilde{f}_{\S;{\boldsymbol{w}}}({\boldsymbol{z}}) ~\overline{T_{\S^c;{\boldsymbol{n}}}({\boldsymbol{z}})}~d \mu_{\S^c}({\boldsymbol{z}})\\
&= \int_{\D_{\S^c}} \left( \sum_{{\boldsymbol{k}} \in \mathcal{I}_{N,d}} \tilde{c}_{{\boldsymbol{k}}} ~T_{\S;{\boldsymbol{k}}}({\boldsymbol{w}})~T_{\S^c;{\boldsymbol{k}}}({\boldsymbol{z}}) \right) \overline{T_{\S^c;{\boldsymbol{n}}}({\boldsymbol{z}})}~d \mu_{\S^c}({\boldsymbol{z}})\\
&= \sum_{{\boldsymbol{k}} \in \mathcal{I}_{N,d}} \left(\tilde{c}_{\S^c;{\boldsymbol{n}}}\right)_{\vect{k}}~T_{\S;{\boldsymbol{k}}}({\boldsymbol{w}}).
\end{align*}
The stated result follows.
\end{proof}

Let $m\in \mathbbm{N}$ and ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$. For any matrix $A \in \mathbbm{C}^{m \times |\mathcal{I}_{N,d}|}$, we define $(A)_{{\boldsymbol{n}}}$ be the column of $A$ corresponding to the index ${\boldsymbol{n}}$. Also, we can choose multiple columns, e.g., for ${\boldsymbol{n}}_1, {\boldsymbol{n}}_2 \in \mathcal{I}_{N,d}$, $(A)_{\{{\boldsymbol{n}}_1,{\boldsymbol{n}}_2\}}$ refers the columns of $A$ corresponding to the indices ${\boldsymbol{n}}_1$ and ${\boldsymbol{n}}_2$.  More generally, for any $\S \subset \mathcal{I}_{N,d}$ the matrix $(A)_S = A_S \in \mathbbm{C}^{m \times |\S|}$ will consist of the columns of $A$ indexed by $\S$.

\subsection{Sampling Matrices associated to a BOS and Restricted Isometry Constants}

Given a BOS as in \eqref{def:BOS_B}, let $\{{\boldsymbol{\xi}}_\ell\}_{\ell \in [m]} \subset \D$ be sampling points drawn independently at random according to $\mu$ with corresponding samples $\{\tilde{{\boldsymbol{y}}}_{\ell} :=\tilde{f}({\boldsymbol{\xi}}_\ell)\}_{\ell \in [m]}$ from $\tilde{f}$ in \eqref{equ:FinApproxf}.  The {\it random sampling matrix} $\Phi \in\C^{m\times |\mathcal{I}_{N,d}|}$ associated with the points $\{{\boldsymbol{\xi}}_\ell\}_{\ell \in [m]}$ and the BOS has entries given by
\begin{equation}
\Phi_{\ell,{\boldsymbol{n}}}=T_{\vect{n}}({\boldsymbol{\xi}}_\ell)
\label{def:A}
\end{equation}
with indices $\ell \in [m]$ and ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$.  One can see that, e.g., $\tilde{{\boldsymbol{y}}} = \Phi {\boldsymbol{\tilde{c}}}$ will hold in this case.  Furthermore, results from the compressive sensing literature guarantee that $\frac{1}{\sqrt{m}} \Phi$ will also have well-behaved {\it restricted isometry constants} as soon as $m$ is sufficiently large.\\

\begin{mydef}[See Definition 6.1 in \cite{foucart2013mathematical}]
The $s$-th restricted isometry constant $\delta_s$ of a matrix $A \in \C^{m \times N}$ is the smallest $\delta \geq 0$ such that
\begin{equation*}
(1-\delta)\|{\boldsymbol{x}}\|_2^2\leq \|A{\boldsymbol{x}}\|_2^2 \leq (1+\delta)\|{\boldsymbol{x}}\|_2^2
\end{equation*}
holds for all $s$-sparse vectors ${\boldsymbol{x}} \in \C^N$.  The matrix $A$ is said to satisfy the restricted isometry property (RIP) of order $(s,\delta)$ if $1 > \delta \geq \delta_s \geq 0$.
\end{mydef}

\begin{mytheorem}[See Theorem 12.32 and Remark 12.33 in \cite{foucart2013mathematical}]
\label{thm:BOS_RIP}
Let $\Phi\in \C^{m\times \left|\mathcal{I}_{N,d} \right|}$ be the random sampling matrix associated to a {BOS} with constant $K\geq 1$ for $m,\left|\mathcal{I}_{N,d} \right|,s \in \mathbbm{Z}^+ \setminus \{1\}$. If, for $\delta,p \in (0,1)$,
\begin{equation*}
m\geq a K^2 \delta^{-2}s \cdot \max \left\{ \ln^2(s) \ln \left( \left|\mathcal{I}_{N,d} \right| \right) \ln(m), \ln(p^{-1}) \right\},
\end{equation*}
then with probability at least $1-p$ the restricted isometry constant $\delta_s$ of $\frac{1}{\sqrt{m}}\Phi$ satisfies $\delta_s \leq \delta$ so that $\Phi$ has the RIP of order $(s,\delta)$. Here the constant $a \in \mathbbm{R}^+$ is universal.
\end{mytheorem}

In addition, one can also, e.g., bound the $\ell_2$ operator norm of the random sampling matrix $\Phi$ in the event that it has the RIP.  We have the following consequence of Theorem~\ref{thm:BOS_RIP}.\\

\begin{lemma}[See Proposition 3.5 in \cite{needell2009cosamp}]
\label{lem:OpBoundRSM}
Suppose $A \in \C^{m \times \left|\mathcal{I}_{N,d} \right| }$ has the restricted isometry property (RIP) of order $(s,\delta)$.  Then,
$$\| A {\boldsymbol{x}} \|_2 \leq \sqrt{1 + \delta} \left( \frac{\| {\boldsymbol{x}} \|_1}{\sqrt{s}} + \| {\boldsymbol{x}} \|_2 \right)$$
holds for all ${\boldsymbol{x}} \in \C^{\left|\mathcal{I}_{N,d} \right|}$.
\end{lemma}

We are now prepared to develop the new CoSaMP variants on which our new sublinear-time algorithms will be based.

\subsection{Exactly sparse case and noisy measurements}\label{sec:Numerics:exactly_sparse}

  We start with the exactly sparse case, and we consider tensor product basis functions with different bases in $d=D$ spatial dimensions, where we choose $d$ up to 100. We set $N=200$ and use $\mathcal{I}_{N,d}=\mathcal{I}_{200,d}$ as search space of possible basis indices, where e.g. $|\mathcal{I}_{200,50}|\approx 10^{115}$ and $|\mathcal{I}_{200,100}|\approx 10^{230}$. We set the maximum number of iterations $\kappa:=20$, and we always use $m_{\rm CE}:=50\,s$ samples for the coefficient estimation where $s=|\S|$.
  For every data point in every plot below, we use 100 different randomly generated trial signals
  \begin{equation} \label{equ:trial_signal_f_sparse}
   f({\boldsymbol{\xi}}) = \sum_{{\boldsymbol{n}}\in\S} c_{{\boldsymbol{n}}} \, T_{{\boldsymbol{n}}}({\boldsymbol{\xi}}),
  \end{equation}
  where we draw the function's support set $\S\subset\mathcal{I}_{N,d}$ uniformly at random without repetition and the coefficients $c_{{\boldsymbol{n}}}\in\{-1,1\}$ uniformly at random.

  Below, a trial will always refer to the execution of Algorithm~\ref{alg:main} on a particular randomly generated trial function $f$ as defined in~\eqref{equ:trial_signal_f_sparse}. A failed trial will refer to any trial where Algorithm~\ref{alg:main} failed to recover the correct support set $\S$ for $f$.

  We assume that the function evaluations of $f$ are contaminated with (white) Gaussian noise, i.e., we provide Algorithm~\ref{alg:main} with noisy samples
  $$
   \boldsymbol{y'} = {\boldsymbol{y}} + \boldsymbol{g'} = {\boldsymbol{y}} + \sigma \frac{\|{\boldsymbol{y}}\|_2}{\|\boldsymbol{g}\|_2} \boldsymbol{g},
  $$
  where ${\boldsymbol{y}}$ contains noiseless samples from $f$, $\boldsymbol{g}\sim \mathcal{N}({\boldsymbol{0}},I)$, and $\sigma\in\mathbbm{R}^+$ is used to control the Signal to Noise Ratio (SNR) defined herein by
  $$
   \mathrm{SNR}_{\mathrm{db}} := 10 \, \log_{10} \left(\frac{\|{\boldsymbol{y}}\|_2^2}{\|\boldsymbol{g'}\|_2^2}\right) = -10 \, \log_{10} (\sigma^2).
 $$

In the following subsections, we consider different types of basis functions. First, in Section~\ref{sec:Numerics:exactly_sparse:mixed1}, we use mixed bases in up to 100 spatial dimensions, which consist of Fourier, Chebyshev, and Legendre bases. Afterwards, we use bases which only consist of Fourier type in Section~\ref{sec:Numerics:exactly_sparse:fourier}, Chebyshev type in Section~\ref{sec:Numerics:exactly_sparse:chebyshev}, and Legendre type in Section~\ref{sec:Numerics:exactly_sparse:legendre}.

\subsubsection{Mixed bases}\label{sec:Numerics:exactly_sparse:mixed1}

First, we consider basis functions $T_{n_j}$ of mixed type: $T_{n_0}$, $T_{n_{\lfloor d/2 \rfloor}}$, and $T_{n_{d-2}}$ are of Chebyshev type; $T_{n_1}$, $T_{n_{\lfloor d/2 \rfloor - 1}}$, and $T_{n_{d-1}}$ are of (preconditioned) Legendre type; and the remaining $d-6$ basis functions $T_{n_j}$, $j\in [d] \setminus\{0,1,\lfloor d/2 \rfloor - 1,\lfloor d/2 \rfloor,d-2,d-1\}$ are of Fourier type.
Preconditioned Legendre type means that instead of using standard Legendre polynomials $L_n(x)=(2n-1)/n \, x \, L_{n-1}(x) - (n-1)/n \, L_{n-2}(x)$, $L_1:=x$, $L_0:=1$, with BOS constant $K=\sqrt{2n+1}$, we apply the preconditioning method from~\cite{R12}, i.e., we use the preconditioned Legendre polynomials $Q_n(x) := \sqrt{\pi/2}\,(1-x^2)^{1/4}\, L_n(x)$ with BOS constant $K=\sqrt{3}$ and choose the sampling nodes randomly with respect to the Chebyshev measure for the basis functions $T_{n_1}$, $T_{n_{\lfloor d/2 \rfloor - 1}}$, and $T_{n_{d-1}}$. Consequently, the overall BOS constant is $\sqrt{2}^3\cdot\sqrt{3}^3\cdot 1^{d-6}=\sqrt{6}^3$ independent of the spatial dimension~$d=D$.
For the entry identification and pairing steps, we set the parameter $m_2=\#{\boldsymbol{z}}_{j,k}:=4\,s$ for different sparsities $s=|\S|$.
The parameter $m_1=\#{\boldsymbol{w}}_{j,\ell}$ is chosen as $c \, s$, where the constant $c\geq 1$ does not depend on the sparsity~$s$, which is distinctly smaller than the theoretical results of $m_1 \sim s^2$ in Theorem~\ref{thm:SuppIDWorks} and Lemma~\ref{lem:SamplingBounds}.

\begin{figure}[!h]
\subfloat[number of samples vs.\ spatial dimension $d$]{
	\label{fig:mixed1:d100:samples_runtime_iter_success_vs_d:samples}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.3\textwidth,
			grid=major,
			width=0.45\textwidth,
			xtick={6,10,25,50,100},
            ymin=0.6e4,
            xmode=log,
            ymode=log,
      xticklabel={
        \pgfkeys{/pgf/fpu=true}
        \pgfmathparse{exp(\tick)}%
        \pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
        \pgfkeys{/pgf/fpu=false}
      },
   			xlabel={$d$},
			ylabel={\#samples},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 2,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,3.570e+04) (10,6.130e+04) (25,1.573e+05) (50,3.173e+05) (100,6.373e+05)
};
\addlegendentry{$s=10$}
\addplot [black,domain=6:100, samples=100, dashed]{6.1e3*x};
\addlegendentry{$6100\,d$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,2.212e+05) (10,3.812e+05) (25,9.812e+05) (50,1.981e+06) (100,3.981e+06)
};
\addlegendentry{$s=25$}
\addplot [blue,domain=6:100, samples=100, dashed]{4e4*x};
\addlegendentry{$4\cdot 10^4 d$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[average runtime vs.\ spatial dimension $d$]{
	\label{fig:mixed1:d100:samples_runtime_iter_success_vs_d:runtime}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.30\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,10,25,50,100},
            ymin=0.2,
            xmode=log,
            ymode=log,
            log ticks with fixed point,
			xlabel={$d$},
			ylabel={avg. runtime in seconds},
		y label style={xshift=-0.8em},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 2,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,2.428e+00) (10,4.067e+00) (25,1.140e+01) (50,2.922e+01) (100,8.350e+01)
};
\addlegendentry{$s=10$}
\addplot [black,domain=6:100, samples=100, dashed]{5 + 0.008*x*x};
\addlegendentry{$0.008d^2\!+\!5$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,1.136e+01) (10,1.797e+01) (25,5.903e+01) (50,1.561e+02) (100,4.912e+02)
};
\addlegendentry{$s=25$}
\addplot [blue,domain=6:100, samples=100, dashed]{26 + 0.045*x*x};
\addlegendentry{$0.045d^2\!+\!26$}
		\end{axis}
		\end{tikzpicture}
}
\\
\subfloat[average iteration vs.\ spatial dimension $d$]{
	\label{fig:mixed1:d100:samples_runtime_iter_success_vs_d:iterations}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.25\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,10,25,50,100},
			ytick={1,2,3,4,5,6},
			ymin=1,ymax=5.7,
            xmode=log,
            log ticks with fixed point,
			xlabel={$d$},
			ylabel={avg. iteration},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = -1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(6,4.190e+00) (10,4.070e+00) (25,4.230e+00) (50,4.730e+00) (100,5.380e+00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(6,4.280e+00) (10,3.630e+00) (25,3.720e+00) (50,3.650e+00) (100,3.730e+00)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[success rate vs.\ spatial dimension~$d$]{
	\label{fig:mixed1:d100:samples_runtime_iter_success_vs_d:success_rate}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.25\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,10,25,50,100},
			ytick={0,0.5,0.8,1},
			ymin=0.0,ymax=1,
            xmode=log,
            log ticks with fixed point,
			xlabel={$d$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = -1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(6,1.000e+00) (10,1.000e+00) (25,9.900e-01) (50,9.900e-01) (100,9.900e-01)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(6,1.000e+00) (10,1.000e+00) (25,1.000e+00) (50,1.000e+00) (100,1.000e+00)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\caption{Number of samples, runtime, number of iterations, success rates
 vs.\ spatial dimension $d=D\in\{6,10,25,50,100\}$ for mixed bases ($3$~Chebyshev, $3$~Legendre, $d-6$ Fourier), $N=200$, sparsity $s\in\{10,25\}$, $\mathrm{SNR}_{\mathrm{db}}=10$, $m_1=8s$, $m_2=4s$.}
\label{fig:mixed1:d100:samples_runtime_iter_success_vs_d}
\end{figure}

In Figure~\ref{fig:mixed1:d100:samples_runtime_iter_success_vs_d}, we visualize the obtained results in dependence of the spatial dimensions $d=D\in\{6,10,25,50,100\}$ for sparsity $s\in\{10,25\}$ and signal to noise ratio $\mathrm{SNR}_{\mathrm{db}}=10$.
In Figure~\ref{fig:mixed1:d100:samples_runtime_iter_success_vs_d:samples}, we plot the number of samples with respect to the spatial dimension $d$. We observe that the number of samples grows nearly linearly in~$d$. Additionally, we plot the average runtime of the 100 test runs with respect to $d$ in Figure~\ref{fig:mixed1:d100:samples_runtime_iter_success_vs_d:runtime}, and we observe that it grows approximately like~$\sim d^2$. When having a look at the average number of iterations in Figure~\ref{fig:mixed1:d100:samples_runtime_iter_success_vs_d:iterations}, we observe that 4.1 to 5.4 iterations were required for sparsity $s=10$ and around 4 iterations for sparsity $s=25$. For the considered test setting, the observed success rate was 100\% for sparsity $s=25$ and at least 99\% for sparsity $s=10$, cf.\ Figure~\ref{fig:mixed1:d100:samples_runtime_iter_success_vs_d:success_rate}.

\begin{figure}[htbp]
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=0$]{
	\label{fig:mixed1:d50:success_vs_m1s:snr0}
\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,8,10},
			ytick={0,0.5,0.8,1},
			xmin=2,xmax=11.5,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(5.0,0.00) (5.5,0.05) (6.0,0.12) (6.5,0.28) (7.0,0.52) (7.5,0.63) (8.0,0.78) (8.5,0.89) (9.0,0.91) (9.5,0.98) (10.0,0.99) (10.5,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(4.0,0.00) (4.5,0.14) (5.0,0.62) (5.5,0.88) (6.0,0.94) (6.5,0.99) (7.0,0.99) (7.5,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(3.0,0.00) (3.5,0.19) (4.0,0.77) (4.5,0.95) (5.0,0.98) (5.5,0.99) (6.0,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=3$]{
	\label{fig:mixed1:d50:success_vs_m1s:snr3}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,8,10},
			ytick={0,0.5,0.8,1},
			xmin=2,xmax=11.5,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(4.5,0.00) (5.0,0.15) (5.5,0.20) (6.0,0.35) (6.5,0.56) (7.0,0.87) (7.5,0.88) (8.0,0.96) (8.5,0.96) (9.0,0.96) (9.5,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.5,0.00) (4.0,0.25) (4.5,0.84) (5.0,0.95) (5.5,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(2.5,0.00) (3.0,0.23) (3.5,0.93) (4.0,0.99) (4.5,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\\
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=6$]{
	\label{fig:mixed1:d50:success_vs_m1s:snr6}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,8,10},
			ytick={0,0.5,0.8,1},
			xmin=2,xmax=11.5,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(3.5,0.00) (4.0,0.03) (4.5,0.03) (5.0,0.23) (5.5,0.31) (6.0,0.58) (6.5,0.76) (7.0,0.91) (7.5,0.95) (8.0,0.99) (8.5,0.99) (9.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.0,0.00) (3.5,0.12) (4.0,0.62) (4.5,0.99) (5.0,0.99) (5.5,1.00)

};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(2.5,0.01) (3.0,0.68) (3.5,0.98) (4.0,0.99) (4.5,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=10$]{
	\label{fig:mixed1:d50:success_vs_m1s:snr10}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,8,10},
			ytick={0,0.5,0.8,1},
			xmin=2,xmax=11.5,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(4.0,0.01) (4.5,0.07) (5.0,0.32) (5.5,0.38) (6.0,0.61) (6.5,0.85) (7.0,0.94) (7.5,0.97) (8.0,0.99) (8.5,0.99) (9.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.0,0.00) (3.5,0.29) (4.0,0.89) (4.5,0.99) (5.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(2.0,0.00) (2.5,0.06) (3.0,0.86) (3.5,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\caption{Success rate vs.\ $m_1/s$ for mixed bases, $d=D=50$, $N=200$, sparsity $s\in\{10,25,50\}$, $m_2=\#{\boldsymbol{z}}_{j,k}=4\,s$, and $\mathrm{SNR}_{\mathrm{db}}\in\{0,3,6,10\}$.}
\label{fig:mixed1:d50:success_vs_m1s}
\end{figure}

For different choices of the parameter $m_1=\#{\boldsymbol{w}}_{j,\ell}\in\{2.5s,3s,3.5s,\ldots,10s\}$, we investigate the success rate for spatial dimension $d=D=50$ and sparsities $s\in\{10,25,50\}$ in Figure~\ref{fig:mixed1:d50:success_vs_m1s}, where we set the signal to noise ratio $\mathrm{SNR}_{\mathrm{db}}$ to 0, 3, 6, and 10 in Figure~\ref{fig:mixed1:d50:success_vs_m1s:snr0}, \ref{fig:mixed1:d50:success_vs_m1s:snr3}, \ref{fig:mixed1:d50:success_vs_m1s:snr6}, and~\ref{fig:mixed1:d50:success_vs_m1s:snr10}, respectively. We observe that the success rates increase for growing parameter~$m_1$. Moreover, the transition between 0\% success rate and 99\%--100\% success rate occurs relatively fast. Additionally, the value $m_1/s$ where the success rate reaches 99\% seems to decrease for increasing sparsity~$s$ and for increasing signal to noise ratio.

\begin{figure}[htbp]
\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.33\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={10,25,50},
            ymin=6e4,ymax=6e6,
            xmode=log,
            ymode=log,
      xticklabel={
        \pgfkeys{/pgf/fpu=true}
        \pgfmathparse{exp(\tick)}%
        \pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
        \pgfkeys{/pgf/fpu=false}
      },
			xlabel={$s$},
			ylabel={\#samples},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 2,
		]
		\addplot[black,mark=*,mark size=2pt,mark options={solid}] coordinates {
(10,396500)  (25,1614950)  (50,5447500)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=0$}
		\addplot[blue,mark=*,mark size=2pt,mark options={solid}] coordinates {
(10,376700)  (25,1367450)  (50,3962500)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=3$}
		\addplot[red,mark=*,mark size=2pt,mark options={solid}] coordinates {
(10,317300)  (25,1119950)  (50,3962500)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=6$}
		\addplot[cyan,mark=*,mark size=2pt,mark options={solid}] coordinates {
(10,317300)  (25,1119950)  (50,3467500)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=10$}
		\end{axis}
		\end{tikzpicture}
\hfill
\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.33\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={10,25,50},
			ymin=3,ymax=1.5e3,
            xmode=log,
            ymode=log,
      xticklabel={
        \pgfkeys{/pgf/fpu=true}
        \pgfmathparse{exp(\tick)}%
        \pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
        \pgfkeys{/pgf/fpu=false}
      },
			xlabel={$s$},
			ylabel={avg. runtime in seconds},
		y label style={xshift=-0.8em},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 2,
		]
		\addplot[black,mark=*,mark size=2pt,mark options={solid}] coordinates {
 (10,32.78) (25,169.43) (50,1.305e3)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=0$}
		\addplot[blue,mark=*,mark size=2pt,mark options={solid}] coordinates {
 (10,30.46) (25,156.52) (50,750.97)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=3$}
		\addplot[red,mark=*,mark size=2pt,mark options={solid}] coordinates {
 (10,30.56) (25,166.59) (50,680.28)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=6$}
		\addplot[cyan,mark=*,mark size=2pt,mark options={solid}] coordinates {
 (10,29.20) (25,153.32) (50,659.03)
};
\addlegendentry{$\mathrm{SNR}_{\mathrm{db}}=10$}
		\end{axis}
		\end{tikzpicture}
\caption{Number of samples and average runtime vs.\ sparsity $s$ for mixed bases, $\geq 99\%$ success rate, $d=D=50$, $N=200$.}
\label{fig:mixed1:d50:samples_runtime_vs_s}
\end{figure}

In Figure~\ref{fig:mixed1:d50:samples_runtime_vs_s}, we plot the used number of samples and average runtime as a function of the sparsity $s\in\{10,25,50\}$ for spatial dimensoin $d=D=50$ and for each signal to noise ratio $\mathrm{SNR}_{\mathrm{db}}\in\{0,3,6,10\}$. We observe that the plots only differ slightly for the different signal to noise ratios $\mathrm{SNR}_{\mathrm{db}}\in\{3,6,10\}$, i.e.\ the numbers of samples and runtimes seem to depend only mildly on the signal to noise ratios for $\geq 99\%$ success rate. In the case $\mathrm{SNR}_{\mathrm{db}}=0$, i.e., when the energy of the signal and of the noise match, the runtimes for sparsities $s\in\{10,25\}$ are similar to the ones of $\mathrm{SNR}_{\mathrm{db}}\in\{3,6,10\}$ and approximately double for $s=50$.

\begin{figure}[htb]
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=0$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={4,5,6,7,8,11},
			ytick={0,0.5,0.8,1},
			xmin=3,xmax=11,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(5.0,0.00) (5.5,0.05) (6.0,0.05) (6.5,0.16) (7.0,0.30) (7.5,0.48) (8.0,0.63) (8.5,0.75) (9.0,0.94) (9.5,0.93) (10.0,0.98) (10.5,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(4.0,0.00) (4.5,0.05) (5.0,0.26) (5.5,0.84) (6.0,0.98) (6.5,1.00)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=3$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,9,10},
			ytick={0,0.5,0.8,1},
			xmin=3,xmax=11,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(4.0,0.00) (4.5,0.01) (5.0,0.07) (5.5,0.16) (6.0,0.27) (6.5,0.41) (7.0,0.62) (7.5,0.81) (8.0,0.88) (8.5,0.93) (9.0,0.99) (9.5,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.5,0.00) (4.0,0.01) (4.5,0.62) (5.0,0.89) (5.5,0.99) (6.0,1.00)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\\
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=6$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,8,9},
			ytick={0,0.5,0.8,1},
			xmin=3,xmax=11,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(4.0,0.00) (4.5,0.06) (5.0,0.13) (5.5,0.27) (6.0,0.31) (6.5,0.54) (7.0,0.73) (7.5,0.88) (8.0,0.96) (8.5,0.97) (9.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.5,0.00) (4.0,0.28) (4.5,0.89) (5.0,0.99) (5.5,1.00)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$\mathrm{SNR}_{\mathrm{db}}=10$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,8,9},
			ytick={0,0.5,0.8,1},
			xmin=3,xmax=11,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(3.5,0.01) (4.0,0.02) (4.5,0.11) (5.0,0.14) (5.5,0.29) (6.0,0.50) (6.5,0.61) (7.0,0.80) (7.5,0.91) (8.0,0.99) (8.5,0.99) (9.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.5,0.01) (4.0,0.56) (4.5,0.99) (5.0,1.00)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\caption{Success rate vs.\ $m_1/s$ for mixed bases, $d=D=100$, $N=200$, sparsity $s\in\{10,25\}$, $m_2=\#{\boldsymbol{z}}_{j,k}=4\,s$, and $\mathrm{SNR}_{\mathrm{db}}\in\{0,3,6,10\}$.}
\label{fig:mixed1:d100:success_vs_m1s}
\end{figure}

Additionally, we repeat the tests of Figure~\ref{fig:mixed1:d50:success_vs_m1s} for spatial dimension $d=D=100$ and sparsities $s\in\{10,25\}$, and we visualize the corresponding results in Figure~\ref{fig:mixed1:d100:success_vs_m1s}. We obtain results analogously to the previous ones.

\FloatBarrier

\newpage

\subsubsection{Fourier bases}\label{sec:Numerics:exactly_sparse:fourier}

As for the case of mixed bases in Section~\ref{sec:Numerics:exactly_sparse:mixed1}, we now perform the numerical tests for tensor products of Fourier bases and show the results in Figure~\ref{fig:fourier:d100:samples_runtime_iter_success_vs_d}. Here, the overall BOS constant $K$ is $1$ independent of the spatial dimension~$d$. Due to the smaller BOS constant, we can reduce the parameters $m_1$ to $5s$ and $m_2$ to $s$ while still obtaining a success rate of 100\%, cf.\ Figure~\ref{fig:fourier:d100:samples_runtime_iter_success_vs_d:success}. As in Section~\ref{sec:Numerics:exactly_sparse:mixed1}, the numbers of samples in Figure~\ref{fig:fourier:d100:samples_runtime_iter_success_vs_d:samples} grow nearly linearly in~$d$ and the average runtimes in Figure~\ref{fig:fourier:d100:samples_runtime_iter_success_vs_d:runtimes} approximately like $\sim d^2$. The average number of iterations in Figure~\ref{fig:fourier:d100:samples_runtime_iter_success_vs_d:iterations} is between $3$ and $4$.

\begin{figure}[!h]
\subfloat[number of samples vs.\ spatial dimension $d$]{
\label{fig:fourier:d100:samples_runtime_iter_success_vs_d:samples}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.3\textwidth,
			grid=major,
			width=0.45\textwidth,
			xtick={6,10,25,50,100},
            ymin=2e3,
            xmode=log,
            ymode=log,
      xticklabel={
        \pgfkeys{/pgf/fpu=true}
        \pgfmathparse{exp(\tick)}%
        \pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
        \pgfkeys{/pgf/fpu=false}
      },
			xlabel={$d$},
			ylabel={\#samples},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = -1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,6.000e+03) (10,1.000e+04) (25,2.500e+04) (50,5.000e+04) (75,7.500e+04) (100,1.000e+05)
};
\addlegendentry{$s=10$}
\addplot [black,domain=6:100, samples=100, dashed]{1e3*x};
\addlegendentry{$\sim d$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,3.562e+04) (10,6.062e+04) (25,1.544e+05) (50,3.106e+05) (75,4.669e+05) (100,6.231e+05)
};
\addlegendentry{$s=25$}
\addplot [forget plot,blue,domain=6:100, samples=100, dashed]{6.2e3*x};
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[average runtime vs.\ spatial dimension $d$]{
\label{fig:fourier:d100:samples_runtime_iter_success_vs_d:runtimes}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.3\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,10,25,50,100},
			ymin=0.006,ymax=100,
            xmode=log,
            ymode=log,
            log ticks with fixed point,
			xlabel={$d$},
			ylabel={avg. runtime in seconds},
		y label style={xshift=-0.8em},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = -1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,7.290e-02) (10,1.499e-01) (25,6.645e-01) (50,2.167e+00) (75,4.625e+00) (100,7.923e+00)
};
\addlegendentry{$s=10$}
\addplot [black,domain=6:100, samples=100, dashed]{8.8e-04*x*x};
\addlegendentry{$\sim d^2$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
(6,4.564e-01) (10,1.078e+00) (25,5.712e+00) (50,2.043e+01) (75,4.718e+01) (100,8.028e+01)
};
\addlegendentry{$s=25$}
\addplot [forget plot,blue,domain=6:100, samples=100, dashed]{0.0082*x*x};
		\end{axis}
		\end{tikzpicture}
}
\\
\subfloat[average iteration vs.\ spatial dimension $d$]{
\label{fig:fourier:d100:samples_runtime_iter_success_vs_d:iterations}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.25\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,10,25,50,100},
			ytick={1,2,3,4,5},
			ymin=1,ymax=4,
            xmode=log,
            log ticks with fixed point,
			xlabel={$d$},
			ylabel={avg. iteration},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = -1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(6,3.210e+00) (10,3.320e+00) (25,3.430e+00) (50,3.370e+00) (75,3.420e+00) (100,3.480e+00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(6,3.540e+00) (10,3.620e+00) (25,3.670e+00) (50,3.500e+00) (75,3.670e+00) (100,3.600e+00)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[success rate vs.\ spatial dimension~$d$]{
	\label{fig:fourier:d100:samples_runtime_iter_success_vs_d:success}
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.25\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,10,25,50,100},
			ytick={0,0.5,0.8,1},
			ymin=0.0,ymax=1,
            xmode=log,
            log ticks with fixed point,
			xlabel={$d$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = -1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(6,1.000) (10,1.000) (25,1.000) (50,1.000) (75,1.000) (100,1.000)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(6,1.000) (10,1.000) (25,1.000) (50,1.000) (75,1.000) (100,1.000)
};
\addlegendentry{$s=25$}
		\end{axis}
		\end{tikzpicture}
}
\caption{Number of samples, runtime, number of iterations, success rate vs.\ spatial dimension $d=D\in\{6,10,25,50,100\}$ for Fourier bases, $N=200$, sparsity $s\in\{10,25\}$, $\mathrm{SNR}_{\mathrm{db}}=10$, $m_1=5s$, $m_2=s$.}
\label{fig:fourier:d100:samples_runtime_iter_success_vs_d}
\end{figure}

\begin{figure}[!h]
\subfloat[$d=6$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={1,1.5,2,2.5,3},
			ytick={0,0.5,0.8,1},
			xmax=3.6,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.09) (2.0,0.85) (2.5,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.01) (1.5,0.97) (2.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(1.0,0.00) (1.5,0.70) (2.0,1.00)
};
\addlegendentry{$s=50$}
		\addplot[brown,mark=*,mark size=1pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.01) (2.0,0.99) (2.5,1.00)
};
\addlegendentry{$s=100$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$d=8$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={1,1.5,2,2.5,3},
			ytick={0,0.5,0.8,1},
			xmax=3.6,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.02) (2.0,0.62) (2.5,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.92) (2.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(1.0,0.00) (1.5,0.69) (2.0,1.00)
};
\addlegendentry{$s=50$}
		\addplot[brown,mark=*,mark size=1pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.00) (2.0,0.96) (2.5,1.00)
};
\addlegendentry{$s=100$}
		\end{axis}
		\end{tikzpicture}
}
\\
\subfloat[$d=10$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={1,1.5,2,2.5,3},
			ytick={0,0.5,0.8,1},
			xmax=3.6,
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.00) (2.0,0.42) (2.5,0.98) (3.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.96) (2.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(1.0,0.00) (1.5,0.68) (2.0,1.00)
};
\addlegendentry{$s=50$}
		\addplot[brown,mark=*,mark size=1pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.00) (2.0,0.96) (2.5,1.00)
};
\addlegendentry{$s=100$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$d=12$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={1,1.5,2,2.5,3},
			ytick={0,0.5,0.8,1},
			ymin=0,ymax=1,
			xmax=3.6,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.00) (2.0,0.34) (2.5,0.94) (3.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.91) (2.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(1.0,0.00) (1.5,0.70) (2.0,1.00)
};
\addlegendentry{$s=50$}
		\addplot[brown,mark=*,mark size=1pt,mark options={solid}] coordinates {
(1.0,0.00) (1.5,0.02) (2.0,0.96) (2.5,1.00)
};
\addlegendentry{$s=100$}
		\end{axis}
		\end{tikzpicture}
}
\caption{Success rate vs.\ $m_1/s$ for Fourier bases, spatial dimension $d=D\in\{6,8,10,12\}$, $N=200$, sparsity $s\in\{10,25,50,100\}$, $\mathrm{SNR}_{\mathrm{db}}=10$, $m_2=s$.}
\label{fig:fourier:success_vs_m1s}
\end{figure}

In Figure~\ref{fig:fourier:success_vs_m1s}, we depict the success rate as a function of $m_1/s\in\{1,1.5,2,2.5,3\}$ for sparsities $s\in\{10,25,50,100\}$ and signal to noise ratio $\mathrm{SNR}_{\mathrm{db}}=10$ in spatial dimensions $d\in\{6,8,10,12\}$. We observe a very small dependence on the spatial dimension~$d$. For $m_1=3s$, the success rate is 100\% in each considered case. Furthermore, there is a rapid transition between full and zero success rate, i.e., the success rate is 0\% for $m_1=s$ each time.

\FloatBarrier

\subsubsection{Chebyshev bases}\label{sec:Numerics:exactly_sparse:chebyshev}

Next, we consider the tensor products of Chebyshev basis functions. Here we expect larger numbers of samples and runtimes compared to the Fourier case in Section~\ref{sec:Numerics:exactly_sparse:fourier} due to the BOS constant $K=\sqrt{2}^d$ for Chebyshev and $K=1$ for Fourier. In particular, for fixed sparsity~$s$ and fixed success rate, the numbers of samples and runtimes might grow for increasing spatial dimension~$d$.

\begin{figure}[!h]
\subfloat[$d=6$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={2,3,4,5,6,7},
			ytick={0,0.5,0.8,1},
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(2.5,0.00) (3.0,0.01) (3.5,0.12) (4.0,0.32) (4.5,0.47) (5.0,0.81) (5.5,0.92) (6.0,0.96) (6.5,0.97) (7.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(2.0,0.00) (2.5,0.01) (3.0,0.29) (3.5,0.84) (4.0,0.98) (4.5,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(2.0,0.00) (2.5,0.10) (3.0,0.74) (3.5,0.95) (4.0,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$d=8$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={3,4,5,6,7,8,10,12,14},
			ytick={0,0.5,0.8,1},
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(4.0,0.00) (5.0,0.03) (6.0,0.14) (7.0,0.34) (8.0,0.65) (9.0,0.78) (10.0,0.91) (11.0,0.92) (12.0,0.97) (13.0,0.98) (14.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.0,0.00) (4.0,0.04) (5.0,0.34) (6.0,0.76) (7.0,0.97) (8.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(3.0,0.00) (3.5,0.02) (4.0,0.23) (4.5,0.59) (5.0,0.87) (5.5,0.98) (6.0,0.98) (6.5,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\\
\subfloat[$d=10$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,8,10,14,18,26,30},
			ytick={0,0.5,0.8,1},
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(8.0,0.00) (10.0,0.06) (12.0,0.25) (14.0,0.53) (16.0,0.70) (18.0,0.81) (20.0,0.92) (22.0,0.92) (24.0,0.96) (26.0,0.99) (28.0,0.99) (30.0,0.99) (32.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(6.0,0.00) (8.0,0.09) (10.0,0.55) (12.0,0.81) (14.0,0.94) (16.0,0.98) (18.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(5.0,0.00) (6.0,0.02) (7.0,0.19) (8.0,0.46) (9.0,0.76) (10.0,0.92) (11.0,0.94) (12.0,0.98) (13.0,0.96) (14.0,0.99) (15.0,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\hfill
\subfloat[$d=12$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.26\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={12,20,28,36,44,52,62},
			ytick={0,0.5,0.8,1},
			ymin=0,ymax=1,
			xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
			ylabel={success rate},
			legend style={legend cell align=left}, legend pos=south east,
			legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(12.0,0.00) (16.0,0.09) (20.0,0.14) (24.0,0.37) (28.0,0.52) (32.0,0.60) (36.0,0.78) (40.0,0.83) (44.0,0.88) (48.0,0.93) (52.0,0.97) (56.0,0.96) (60.0,0.98) (62.0,1.00)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(12.0,0.00) (14.0,0.04) (16.0,0.19) (18.0,0.20) (20.0,0.42) (22.0,0.59) (24.0,0.69) (26.0,0.76) (28.0,0.85) (30.0,0.90) (32.0,0.95) (34.0,0.94) (36.0,0.94) (38.0,0.97) (40.0,0.97) (42.0,0.98) (44.0,0.98) (46.0,1.00)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(10.0,0.00) (12.0,0.02) (14.0,0.26) (16.0,0.41) (18.0,0.73) (20.0,0.82) (22.0,0.83) (24.0,0.90) (26.0,0.98) (28.0,0.97) (30.0,0.97) (32.0,0.96) (34.0,0.97) (36.0,0.97) (38.0,1.00)
};
\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
}
\caption{Success rate vs.\ $m_1/s$ for Chebyshev bases, spatial dimension $d=D\in\{6,8,10,12\}$, $N=200$, sparsity $s\in\{10,25,50\}$, $\mathrm{SNR}_{\mathrm{db}}=10$, $m_2=4s$.}
\label{fig:chebyshev:success_vs_m1s}
\end{figure}

In Figure~\ref{fig:chebyshev:success_vs_m1s}, we depict the success rate as a function of $m_1/s$ for sparsities $s\in\{10,25,50\}$ and signal to noise ratio $\mathrm{SNR}_{\mathrm{db}}=10$ in spatial dimensions $d\in\{6,8,10,12\}$. As predicted, we observe that we have to increase $m_1=\#{\boldsymbol{w}}_{j,\ell}$ distinctly for growing spatial dimension~$d$ and fixed sparsity~$s$. For instance, for $s=25$, we observe a 98\% success rate for $m_1=4s$ and $d=6$, but obtain a success rate of only 4\% for $d=8$. For $d=12$, we had to choose $m_1=42s$ to achieve a success rate of 98\%.

\begin{figure}[htbp]
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.33\textwidth,
			grid=major,
			width=0.48\textwidth,
			xtick={6,8,10,12},
			ytick={1e4,1e5,1e6,1e7},
            xmode=log,
            ymode=log,
      xticklabel={
        \pgfkeys{/pgf/fpu=true}
        \pgfmathparse{exp(\tick)}%
        \pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
        \pgfkeys{/pgf/fpu=false}
      },
			ymin=1e4,ymax=1e8,
			xlabel={$d$},
			ylabel={\#samples},
			legend style={legend cell align=left,at={(0.97,1.1)}},
			legend columns = 3,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
 (6,3.130e+04) (8,8.450e+04) (10,1.981e+05) (12,570900)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
 (6,125550) (8,3.012e+05) (10,8.562e+05) (12,2646250)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180},only marks] coordinates {
 (6,4.425e+05) (8,977500) (10,2.662e+06) (12,8742500)
};
\addlegendentry{$s=50$}
\addplot [black,domain=6:12, samples=100, dashed]{2^10.7965 * 2^(0.69*x)};
\addlegendentry{$\sim 2^{0.69d}$}
\addplot [blue,domain=6:12, samples=100, dashed]{2^12.220807844665455 * 2^(0.75*x)};
\addlegendentry{$\sim 2^{0.75d}$}
\addplot [red,domain=6:12, samples=100, dashed]{2^14.3032 * 2^(0.72*x)};
\addlegendentry{$\sim 2^{0.72d}$}
		\end{axis}
		\end{tikzpicture}
\hfill
		\begin{tikzpicture}[baseline]
		\begin{axis}[
			font=\footnotesize,
			enlarge x limits=true,
			enlarge y limits=true,
			height=0.33\textwidth,
			grid=major,
			width=0.48\textwidth,
			ymode=log,
			xmode=log,
			log ticks with fixed point,
			xtick={4,6,8,10,12},
			ytick={0.1,1,10,100},
			ymin=1e-1,ymax=2e3,
			xlabel={$d$},
			ylabel={avg. runtime in seconds},
			legend style={legend cell align=left,at={(0.97,1.1)}},
			legend columns = 3,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
		 (6,0.1606) (8,0.3843) (10,0.8007) (12,1.9203)
};
\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
		 (6,0.9689) (8,2.2925) (10,5.7969) (12,18.3088)
};
\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180},only marks] coordinates {
		 (6,7.0094) (8,15.8067) (10,39.7017) (12,209.2183)
};
\addlegendentry{$s=50$}
\addplot [black,domain=6:12, samples=100, dashed]{2^(-6.191117657342546) * 2^(0.58*x)};
\addlegendentry{$\sim 2^{0.58d}$}
\addplot [blue,domain=6:12, samples=100, dashed]{2^(-4.230602954903033) * 2^(0.69*x)};
\addlegendentry{$\sim 2^{0.69d}$}
\addplot [red,domain=6:12, samples=100, dashed]{2^(-2.2594) * 2^(0.80*x)};
\addlegendentry{$\sim 2^{0.80d}$}
		\end{axis}
		\end{tikzpicture}
\caption{Number of samples and average runtime vs.\ spatial dimension $d\in\{6,8,10,12\}$ for Chebyshev bases, $N=200$, sparsity $s\in\{10,25,50\}$, $\mathrm{SNR}_{\mathrm{db}}=10$, $m_2=4s$, success rate $\geq 99\%$.}
\label{fig:chebyshev:samples_runtime_vs_d}
\end{figure}

In Figure~\ref{fig:chebyshev:samples_runtime_vs_d}, we investigate the dependence of the spatial dimension~$d$ on the number of samples and average runtime for the case of $\geq 99\%$ success rate in more detail. For our test cases, we observe that the numbers of samples grow approximately like between $\sim 2^{0.69d}$ and $\sim 2^{0.75d}$ as well as the runtimes approximately like between $\sim 2^{0.58d}$ and $\sim 2^{0.80d}$. In each case, this is distinctly less than the worst case upper bounds in Theorem~\ref{thm:NewCoSaMP} suggest.

\FloatBarrier

\subsubsection{Preconditioned Legendre bases}\label{sec:Numerics:exactly_sparse:legendre}

Here, we consider the tensor products of preconditioned Legendre basis functions $Q_n$ with BOS constant $K=\sqrt{3}^d$, cf.\ Section~\ref{sec:Numerics:exactly_sparse:mixed1}.
In Figure~\ref{fig:legendre:success_vs_m1s}, we show the success rates as a function of $m_1/s$ for sparsities $s\in\{10,25,50\}$ and signal to noise ratio $\mathrm{SNR}_{\mathrm{db}}=10$ in spatial dimensions $d\in\{6,8,10,12\}$.
As in the case of Chebyshev bases, we observe that we have to increase $m_1$ distinctly for growing spatial dimension~$d$ and fixed sparsity~$s$. For instance, for $s=25$, we observe a 97\% success rate for $m_1=4s$ and $d=6$ as well as 100\% for $m_1=4.5s$ and $d=6$, but obtained a success rate of only 2\% for $m_1=4s$ and $d=8$. Moreover, we had to choose $m_1=18s$ to have a success rate of 100\% for $d=10$ and $m_1=56s$ for $d=12$.

When comparing the obtained results with the ones for the Chebyshev bases, we do not numerically observe the higher BOS constant $K=\sqrt{3}^d$ here for $d=6,8,10$. The plots in Figure~\ref{fig:legendre:success_vs_m1s} look very similar to the ones in Figure~\ref{fig:chebyshev:success_vs_m1s}. For $d=12$, the values of $m_1$ where a success rate of $\geq 99\%$ is reached are slightly larger than the ones in the Chebyshev case.

\begin{figure}[!h]
	\subfloat[$d=6$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={2,3,4,5,6,7},
		ytick={0,0.5,0.8,1},
		ymin=0,ymax=1,
		xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
		ylabel={success rate},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(3.0,0.00) (3.5,0.09) (4.0,0.24) (4.5,0.64) (5.0,0.85) (5.5,0.88) (6.0,0.98) (6.5,0.99) (7.0,1.00)
		};
		\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(2.0,0.00) (2.5,0.06) (3.0,0.39) (3.5,0.82) (4.0,0.97) (4.5,1.00)
		};
		\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(2.0,0.00) (2.5,0.05) (3.0,0.69) (3.5,0.91) (4.0,1.00)
		};
		\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[$d=8$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={3,4,5,6,7,8,10,12,14},
		ytick={0,0.5,0.8,1},
		ymin=0,ymax=1,
		xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
		ylabel={success rate},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(4.0,0.00) (5.0,0.01) (6.0,0.10) (7.0,0.33) (8.0,0.71) (9.0,0.80) (10.0,0.95) (11.0,0.99) (12.0,1.00)
		};
		\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(3.0,0.00) (4.0,0.02) (5.0,0.36) (6.0,0.80) (7.0,0.97) (8.0,0.99) (9.0,1.00)
		};
		\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(3.0,0.00) (4.0,0.25) (5.0,0.90) (6.0,0.98) (7.0,1.00)
		};
		\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\\
	%
	\subfloat[$d=10$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={6,8,10,15,18,30,34},
		ytick={0,0.5,0.8,1},
		ymin=0,ymax=1,
		xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
		ylabel={success rate},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(6.0,0.00) (8.0,0.01) (10.0,0.07) (12.0,0.28) (14.0,0.52) (16.0,0.80) (18.0,0.87) (20.0,0.98) (22.0,0.96) (24.0,0.97) (26.0,1.00) (28.0,0.98) (30.0,0.99) (32.0,0.99) (34.0,1.00)
		};
		\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(6.0,0.00) (8.0,0.09) (10.0,0.59) (12.0,0.88) (14.0,0.95) (16.0,0.96) (18.0,1.00)
		};
		\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(5.0,0.00) (6.0,0.01) (7.0,0.17) (8.0,0.52) (9.0,0.80) (10.0,0.91) (11.0,0.96) (12.0,0.98) (13.0,0.99) (14.0,0.99) (15.0,1.00)
		};
		\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[$d=12$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={12,28,36,48,56,64,72,80},
		ytick={0,0.5,0.8,1},
		ymin=0,ymax=1,
		xlabel={$m_1/s\;=\;$(\#${\boldsymbol{w}}_{j,\ell})/s$},
		ylabel={success rate},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(12.0,0.00) (16.0,0.07) (20.0,0.18) (24.0,0.47) (28.0,0.58) (32.0,0.74) (36.0,0.84) (40.0,0.88) (44.0,0.95) (48.0,0.92) (52.0,0.91) (56.0,0.97) (60.0,0.98) (64.0,0.99) (68.0,0.99) (72.0,0.98) (76.0,0.99) (80.0,1.00)
		};
		\addlegendentry{$s=10$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(12.0,0.00) (14.0,0.04) (16.0,0.11) (18.0,0.38) (20.0,0.44) (22.0,0.59) (24.0,0.71) (26.0,0.85) (28.0,0.89) (30.0,0.93) (32.0,0.90) (34.0,0.96) (36.0,0.96) (38.0,0.97) (40.0,0.98) (42.0,0.97) (44.0,0.98) (46.0,0.98) (48.0,0.99) (50.0,0.99) (52.0,0.99) (54.0,0.99) (56.0,1.00)
		};
		\addlegendentry{$s=25$}
		\addplot[red,mark=pentagon,mark size=2.5pt,mark options={solid,rotate=180}] coordinates {
(10.0,0.00) (12.0,0.10) (14.0,0.24) (16.0,0.44) (18.0,0.68) (20.0,0.79) (22.0,0.88) (24.0,0.95) (26.0,0.96) (28.0,0.98) (30.0,0.96) (32.0,0.97) (34.0,0.97) (36.0,0.99) (38.0,1.00)
		};
		\addlegendentry{$s=50$}
		\end{axis}
		\end{tikzpicture}
	}
	\caption{Success rate vs.\ $m_1/s$ for preconditioned Legendre bases, spatial dimension $d=D\in\{6,8,10,12\}$, $N=200$, sparsity $s\in\{10,25,50\}$, $\mathrm{SNR}_{\mathrm{db}}=10$, $m_2=4s$.}
	\label{fig:legendre:success_vs_m1s}
\end{figure}

\FloatBarrier

\subsection{Approximately sparse case}
After considering exactly sparse test functions~$f$ in Section~\ref{sec:Numerics:exactly_sparse}, we continue with examples for the approximately sparse case, i.e., our test functions under consideration will have infinitely many non-zero basis coefficients~$c_{\boldsymbol{n}}$.

\subsubsection{Fourier type with $D=10$}
We use the 10-variate periodic test function $f\colon\mathbbm{T}^{10}\rightarrow\mathbbm{R}$,
\begin{equation} \label{equ:f:10}
f({\boldsymbol{\xi}}):=\prod_{t\in\{0,2,7\}}N_2(\xi_t) + \prod_{t\in\{1,4,5,9\}}N_4(\xi_t) + \prod_{t\in\{3,6,8\}}N_6(\xi_t),
\end{equation}
from \cite[Section~3.3]{potts2016sparse} and \cite[Section~5.3]{kammerer2017high} with infinitely many non-zero Fourier coefficients~$c_{\boldsymbol{n}}$,
where $\mathbbm{T}\simeq[0,1)$ is the torus and $N_m:\mathbbm{T}\rightarrow\mathbbm{R}$ is the B-Spline of order $m\in\mathbbm{N}$,
$$N_m(x) := C_m \sum_{n\in\mathbbm{Z}} \operatorname{sinc}\left(\frac{\pi}{m}n\right)^m (-1)^n \,\mathrm{e}^{2\pi\mathrm{i}nx},$$
with a constant $C_m>0$ such that $\Vert N_m \Vert_{L^2(\mathbbm{T})}=1$.
We remark that each B-Spline $N_m$ of order $m\in\mathbbm{N}$ is a piece-wise polynomial of degree $m-1$.
We approximate the function $f$ by multivariate trigonometric polynomials $a$ using Algorithm~\ref{alg:main}.
The obtained basis index sets $\tilde{\Omega}$ should ``consist of'' the union of three lower dimensional manifolds,
a three-dimensional hyperbolic cross in the dimensions $1,3,8$;
a four-dimensional hyperbolic cross in the dimensions $2,5,6,10$;
and a three-dimensional hyperbolic cross in the dimensions $4,7,9$.
All tests are performed 10 times and the relative $L^2(\mathbbm{T}^{10})$  approximation error
$$
\frac{ \Vert f-a\Vert_{L^2(\mathbbm{T}^{10})} }{ \Vert f\Vert_{L^2(\mathbbm{T}^{10})} }
=
\frac{
\sqrt{\Vert f\Vert_{L^2(\mathbbm{T}^{10})}^2 - \sum_{\boldsymbol{n}\in \tilde{\Omega}}\vert\hat{f}_{\boldsymbol{n}}\vert^2 + \sum_{\boldsymbol{n}\in \tilde{\Omega}}\vert a_{\boldsymbol{n}}-\hat{f}_{\boldsymbol{n}}\vert^2} }{ \Vert f\Vert_{L^2(\mathbbm{T}^{10})} }
$$
is computed each time,
where the approximant $a:=\sum_{\boldsymbol{n}\in \tilde{\Omega}} a_{\boldsymbol{n}} \,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{n}\cdot\circ}$.

We set the parameters $N=64$, $d=D=10$, $m_2=s$, and we always use $m_{\rm CE}:=50\,s$ samples for the coefficient estimation where $s=|\S|$.
 %
For our tests, we consider two different parameter combinations: $m_1=3s$ and $\kappa=20$, as well as $m_1=8s$ and $\kappa=10$. The obtained results, i.e., the numbers of samples, average runtimes, average iterations, and relative $L^2(\mathbbm{T}^{10})$ errors are plotted as a function of the sparsity~$s\in\{100,200,500,1000\}$ in Figure~\ref{fig:approx_sparse:fourier_d10:samples_runtime_iter_success_vs_s}. Due to the parameter choices for $m_1$ and $m_2$, we observe that the numbers of samples grow quadratically for increasing sparsity~$s$. The average runtimes grow approximately like $\sim s^3\min\{s,N\}$ and this means $\sim s^3$ for fixed $N$.
Moreover, the average numbers of iterations are much smaller than its imposed maximum $\kappa$ in most cases. The relative $L^2(\mathbbm{T}^{10})$ errors decrease for increasing sparsity $s$ having a value of approximately $10^{-2}$ for sparsity $s=1000$. Again, we emphasize the extremely high power of Algorithm~\ref{alg:main}, which is able to determine the $s=1000$ approximately largest bases coefficients and the corresponding indices for our test function out of $| \mathcal{I}_{N,d}| = N^d = 64^{10} \approx 10^{18}$ allowed indices.

\begin{figure}[!h]
	\subfloat[numbers of samples vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.3\textwidth,
		grid=major,

		width=0.47\textwidth,
		xtick={100,200,500,1000},
		xmode=log,
		ymode=log,
		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={\#samples},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(100,5.750e+05) (200,2.290e+06) (500,1.428e+07) (1000,5.705e+07)
		};
	    \addlegendentry{$m_1=3s$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(100,1.525e+06) (200,6.090e+06) (500,3.802e+07) (1000,1.520e+08)
        };
        \addlegendentry{$m_1=8s$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[average runtimes vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.3\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={100,200,500,1000},
		ytick={1e2,1e3,1e4},

        ymax=8e4,
		xmode=log,
		ymode=log,
		ytick={85,430,5400,42000},
		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. runtime in seconds},
		y label style={xshift=-0.8em},
		legend style={legend cell align=left,at={(1,1.3)}},
		legend columns = 2,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
(100,1.102e+02) (200,3.762e+02) (500,3.934e+03) (1000,3.260e+04)
		};
	    \addlegendentry{$m_1=3s$, $\kappa=20$}
		\addplot [black,domain=100:1000, samples=100, dashed]{29+3e-5*x*x*x};
		\addlegendentry{$3\cdot 10^{-5}s^3+29$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
(100,1.139e+02) (200,3.930e+02) (500,5.946e+03) (1000,5.023e+04)
		};
	    \addlegendentry{$m_1=8s$, $\kappa=10$}
		\addplot [black,domain=100:1000, samples=100, dashed]{5e-5*x*x*x};
\addlegendentry{$5\cdot 10^{-5}s^3$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\\
	%
	\subfloat[average iteration vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.3\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={100,200,500,1000},
		ytick={3,5,7,18},

		ymin=3,ymax=20,
		ytick={3,5,8,10,18},
		xmode=log,

		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. iteration},
		legend style={legend cell align=left}, legend pos=north east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(100,1.770e+01) (200,9.200e+00) (500,7.500e+00) (1000,8.200e+00)
		};
	    \addlegendentry{$m_1=3s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(100,4.600e+00) (200,4.100e+00) (500,4.600e+00) (1000,5.000e+00)
		};
	    \addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[$L^2(\mathbbm{T}^{10})$ error vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.3\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={100,200,500,1000},

		ymin=1e-2,ymax=1,
		xmode=log,
		ymode=log,

		xticklabel={
	\pgfkeys{/pgf/fpu=true}
	\pgfmathparse{exp(\tick)}%
	\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
	\pgfkeys{/pgf/fpu=false}
},
		xlabel={$s$},
		ylabel={relative error},
		legend style={legend cell align=left}, legend pos=south west,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(100,3.969e-01) (200,1.870e-01) (500,5.690e-02) (1000,1.240e-02)
		};
	    \addlegendentry{$m_1=3s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(100,3.504e-01) (200,1.856e-01) (500,5.683e-02) (1000,1.238e-02)
		};
	    \addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	\caption{Number of samples, runtime, number of iterations, $L^2(\mathbbm{T}^{10})$ error
		vs.\ sparsity $s\in\{100,200,500,1000\}$ for Fourier basis and test function~\eqref{equ:f:10}.}
	\label{fig:approx_sparse:fourier_d10:samples_runtime_iter_success_vs_s}
\end{figure}

\FloatBarrier

\subsubsection{Chebyshev and Legendre type with $D=7$}

\begin{figure}[!h]
	\centering
	\begin{tikzpicture}[font=\scriptsize]
	\begin{axis}[width=0.4\textwidth,height=0.22\textwidth,enlarge x limits=false,enlarge y limits=true,xmin=-5,xmax=2,ymin=0,ymax=0.8,xtick={-1,0,1}
	]
	\addplot[samples=100,black!70,domain=-5:-4.5] {0};
	\addplot[samples=100,black!70,domain=-4.5:-2.5] {x*x/8+9*x/8+81/32};
	\addplot[samples=100,black!70,domain=-2.5:-1] {-x*x/4-3*x/4+3/16};
	\addplot[samples=100,blue,domain=-1:-0.5,ultra thick] {-x*x/4-3*x/4+3/16};
	\addplot[samples=100,blue,domain=-0.5:1,ultra thick] {x*x/8-3*x/8+9/32};
	\addplot[samples=100,black!70,domain=1:1.5] {x*x/8-3*x/8+9/32};
	\addplot[samples=100,black!70,domain=1.5:10] {0};
	\node [right] at (axis cs:-4.2,0.7) {$B_3$};
	\end{axis}
	\end{tikzpicture}
	\hspace{1em}
	\begin{tikzpicture}[font=\scriptsize]
	\begin{axis}[width=0.6\textwidth,height=0.22\textwidth,enlarge x limits=false,enlarge y limits=true,xmin=-8,xmax=3,ymin=0,ymax=0.8,xtick={-1,0,1}
	]
	\addplot[samples=100,black!70,domain=-10:-7.5] {0};
	\addplot[samples=100,black!70,domain=-7.5:-5.5] {(2*x+15)^4/6144};
	\addplot[samples=100,black!70,domain=-5.5:-3.5] {-5645/1536-205/48*x-95/64*x*x-5/24*x^3-1/96*x^4};
	\addplot[samples=100,black!70,domain=-3.5:-1.5] {715/3072+25/128*x+55/128*x*x+5/32*x^3+1/64*x^4};
	\addplot[samples=100,black!70,domain=-1.5:-1] {155/1536-5/32*x+5/64*x*x-1/96*x^4};
	\addplot[samples=100,blue,domain=-1:0.5,ultra thick] {155/1536-5/32*x+5/64*x*x-1/96*x^4};
	\addplot[samples=100,blue,domain=0.5:1,ultra thick] {(2*x-5)^4/6144};
	\addplot[samples=100,black!70,domain=1:2.5] {(2*x-5)^4/6144};
	\addplot[samples=100,black!70,domain=2.5:10] {0};
	\node [right] at (axis cs:-5,0.7) {$B_5$};
	\end{axis}
	\end{tikzpicture}
	\vspace{-0.6em}
	\caption{B-splines $B_3$ and $B_5$ considered in interval $[-1,1]$.}
	\label{fig:Bspline:nonper}
\end{figure}

Next, we apply Algorithm~\ref{alg:main} on the 7-variate test function $f\colon [-1,1]^{7}\rightarrow\mathbbm{R}$,
\begin{equation} \label{equ:f_cheb:7}
f({\boldsymbol{\xi}}):=\prod_{t\in\{0,2,5\}}B_3(\xi_t) + \prod_{t\in\{1,3,4,6\}}B_5(\xi_t)
\vspace{-0.5em}
\end{equation}
similar as in \cite{potts2017multivariate},
where
$B_3\colon\mathbbm{R}\rightarrow\mathbbm{R}$ is a shifted, scaled and dilated B-spline of order~3
and
$B_5\colon\mathbbm{R}\rightarrow\mathbbm{R}$ is a shifted, scaled and dilated B-spline of order~5, see Figure~\ref{fig:Bspline:nonper} for illustration.
We remark that the absolute values of the Chebyshev coefficients $c_n$, $n\in\mathbbm{N}_0$,
of $B_3$ and $B_5$ decay like $\sim n^{-3}$ and $\sim n^{-5}$, respectively.
The obtained basis index sets $\tilde{\Omega}$ should ``consist of'' the union of two lower dimensional manifolds,
a three-dimensional hyperbolic cross in the dimensions $0,2,5$; and
a four-dimensional hyperbolic cross in the dimensions $1,3,4,6$.
All tests are performed 10 times and the relative $L^2([-1,1]^7,\mu_\mathrm{C})$ approximation error
$\Vert f-a\Vert_{L^2([-1,1]^7,\mu_\mathrm{C})} / \Vert f\Vert_{L^2([-1,1]^7,\mu_\mathrm{C})}$
is computed each time,
where the approximant $a:=\sum_{{\boldsymbol{n}}\in \tilde{\Omega}} a_{{\boldsymbol{n}}} \,T_{\boldsymbol{n}}$, $T_{\boldsymbol{n}}$ is the Chebyshev product basis, and $\mu_\mathrm{C}({\boldsymbol{\xi}}):=\pi^{-D} \prod_{j\in [D]} (1-\xi_j^2)^{-1/2}$ is the Chebyshev product measure.

We set the parameters $N=64$, $d=D=7$, $m_2=4s$, and we always use $m_{\rm CE}:=50\,s$ samples for the coefficient estimation where $s=|\S|$.
 %
We consider two different parameter combinations: $m_1=4s$ and $\kappa=20$, as well as $m_1=8s$ and $\kappa=10$. The obtained results, i.e., the numbers of samples, average runtimes, average iterations, and relative $L^2([-1,1]^7,\mu_\mathrm{C})$ errors are plotted as a function of the sparsity~$s\in\{25,50,100,200,500\}$ in Figure~\ref{fig:approx_sparse:chebyshev_d7:samples_runtime_iter_success_vs_s}. Due to the parameter choices for $m_1$ and $m_2$, we observe that the numbers of samples grow quadratically for increasing sparsity~$s$. The average runtimes grow approximately like $\sim s^3\min\{s,N\}$ and this means $\sim s^3$ for fixed $N$.
Moreover, the average numbers of iterations are well below its imposed maximum $\kappa$ for $m_1=4s$ and $\kappa=20$ as well as close to $\kappa$ for $m_1=8s$ and $\kappa=10$. The relative $L^2([-1,1]^7,\mu_\mathrm{C})$ errors decrease for increasing sparsity $s$ having a value of approximately $2.3\cdot 10^{-4}$ for sparsity $s=500$.
We emphasize that Algorithm~\ref{alg:main} is able to easily determine the $s=500$ approximately largest basis coefficients and the corresponding basis indices for our test function out of $| \mathcal{I}_{N,d}| = N^d = 64^{7} \approx 4.4\cdot 10^{12}$ allowed indices.

\begin{figure}[!h]
	\subfloat[number of samples vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.47\textwidth,
		xtick={25,50,100,200,500},
		xmode=log,
		ymode=log,
		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={\#samples},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
        	(25,131250) (50,522500) (100,2085000) (200,8330000) (500,52025000)
        };
        \addlegendentry{$m_1=4s$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
			(25,261250) (50,1042500) (100,4165000) (200,16650000) (500,104025000)
		};
		\addlegendentry{$m_1=8s$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[average runtime vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={1,100,20000},

		ymin=1,
		ymax=2e5,
		ytick={1,10,190,1000,15000},
		xmode=log,
		ymode=log,
		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. runtime in sec.},
		y label style={xshift=-0.8em},
		legend style={legend cell align=left,at={(1,1.1)}},
		legend columns = -1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid},only marks] coordinates {
(25,1.428e+00) (50,1.301e+01) (100,1.772e+02) (200,7.521e+02) (500,1.327e+04)
        };
        \addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot [black,domain=25:500, samples=100, dashed]{1.06e-4*x*x*x};
\addlegendentry{$\sim s^3$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid},only marks] coordinates {
(25,1.422e+00) (50,1.523e+01) (100,1.867e+02) (200,1.072e+03) (500,1.935e+04)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\\
	%
	\subfloat[average iteration vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={5,7,10,13},

		ymin=1,ymax=13,
		xmode=log,

		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. iteration},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,9.900e+00) (50,1.200e+01) (100,1.320e+01) (200,1.100e+01) (500,1.260e+01)
        };
        \addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,5.400e+00) (50,6.900e+00) (100,7.200e+00) (200,8.000e+00) (500,9.300e+00)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[relative {$L^2([-1,1]^7,\mu_\mathrm{C})$} error vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},

		ymin=6e-6,ymax=1,
		xmode=log,
		ymode=log,

		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={relative error},
		legend style={legend cell align=left}, legend pos=south west,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,4.726e-02) (50,1.556e-02) (100,5.592e-03) (200,1.495e-03) (500,2.287e-04)
        };
        \addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,4.242e-02) (50,1.549e-02) (100,5.514e-03) (200,1.493e-03) (500,2.279e-04)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	\caption{Number of samples, runtime, number of iterations, relative $L^2([-1,1]^7,\mu_\mathrm{C})$ error
		vs.\ sparsity $s\in\{25,50,100,200,500\}$ for Chebyshev basis and test function~\eqref{equ:f_cheb:7}.}
	\label{fig:approx_sparse:chebyshev_d7:samples_runtime_iter_success_vs_s}
\end{figure}

In addition, we use the preconditioned Legendre polynomials $Q_n$ from Section~\ref{sec:Numerics:exactly_sparse:mixed1} as basis functions, i.e.\ $T_{\boldsymbol{n}}$ is now the Legendre product basis in the approximant $a:=\sum_{{\boldsymbol{n}}\in \tilde{\Omega}} a_{{\boldsymbol{n}}} \,T_{\boldsymbol{n}}$. Besides that, we keep all parameters identical but determine the relative $L^2([-1,1]^7,\mu_\mathrm{L})$ approximation error
$\Vert f-a\Vert_{L^2([-1,1]^7,\mu_\mathrm{L})} / \Vert f\Vert_{L^2([-1,1]^7,\mu_\mathrm{L})}$, which corresponds to the Legendre basis and uses the probability measure $\mu_\mathrm{L} \equiv 2^{-D}$. The results are shown in Figure~\ref{fig:approx_sparse:legendre_d7:samples_runtime_iter_success_vs_s}. Here, we observe that the numbers of iterations are higher by up to $\approx$ 50\% compared to the Chebyshev case in Figure~\ref{fig:approx_sparse:chebyshev_d7:samples_runtime_iter_success_vs_s}, and that they reach the imposed maximum of $\kappa:=20$ for $m_1=4s$ and $\kappa:=10$ for $m_2=8s$ in several cases. Correspondingly, the runtimes are also higher by up to $\approx$ 50\%. The obtained relative $L^2([-1,1]^7,\mu_\mathrm{L})$ errors are similar, but we also remark that we cannot compare these errors directly to the relative $L^2([-1,1]^7,\mu_\mathrm{C})$ errors of the Chebyshev basis since they are measured with respect to different probability measures, $\mu_\mathrm{C}({\boldsymbol{\xi}}):= \pi^{-D} \prod_{j\in [D]} (1-\xi_j^2)^{-1/2}$ for Chebyshev and $\mu_\mathrm{L} \equiv 2^{-D}$ for Legendre.

\begin{figure}[!h]
	\subfloat[number of samples vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.47\textwidth,
		xtick={25,50,100,200,500},
		xmode=log,
		ymode=log,
		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={\#samples},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
	(25,131250) (50,522500) (100,2085000) (200,8330000) (500,52025000)
};
\addlegendentry{$m_1=4s$}
\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
	(25,261250) (50,1042500) (100,4165000) (200,16650000) (500,104025000)
};
\addlegendentry{$m_1=8s$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[average runtime vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={10,1000,30000},

		ymax=2e5,
		xmode=log,
		ymode=log,
		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. runtime in sec.},
		y label style={xshift=-0.8em},
		legend style={legend cell align=left}, legend pos=north west,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,7.210e+00) (50,5.395e+01) (100,3.555e+02) (200,1.931e+03) (500,2.612e+04)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,7.130e+00) (50,5.524e+01) (100,3.171e+02) (200,1.971e+03) (500,2.727e+04)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\\
	%
	\subfloat[average iteration vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={5,10,15,20},

		ymin=-1,ymax=20,
		xmode=log,

		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. iteration},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,1.020e+01) (50,1.610e+01) (100,1.840e+01) (200,1.820e+01) (500,1.940e+01)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,8.300e+00) (50,9.300e+00) (100,9.400e+00) (200,9.200e+00) (500,1.000e+01)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[relative {$L^2([-1,1]^7,\mu_\mathrm{L})$} error vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={1,1e-2,1e-4},

		ymin=6e-6,ymax=1,
		xmode=log,
		ymode=log,

		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={relative error},
		legend style={legend cell align=left}, legend pos=south west,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,7.240e-02) (50,2.124e-02) (100,5.603e-03) (200,1.460e-03) (500,2.178e-04)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,3.757e-02) (50,1.394e-02) (100,4.879e-03) (200,1.421e-03) (500,2.109e-04)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	\caption{
		Number of samples, runtime, number of iterations, relative $L^2([-1,1]^7,\mu_\mathrm{L})$ error
		vs.\ sparsity $s\in\{25,50,100,200,500\}$ for Legendre basis and test function~\eqref{equ:f_cheb:7}.

	}
	\label{fig:approx_sparse:legendre_d7:samples_runtime_iter_success_vs_s}
\end{figure}

\FloatBarrier

\subsubsection{Mixed type with $D=10$}

\begin{sloppypar}
Finally, we combine parts of the test functions from the previous two subsections. We consider the 10-variate test function
$f\colon \tilde{\D}\rightarrow\mathbbm{R}$, $\tilde{\D}:=[-1,1]\times\mathbbm{T}\times[-1,1]^3\times\mathbbm{T}^2\times[-1,1]\times\mathbbm{T}^2$,
\begin{equation} \label{equ:f_mixed:10}
f({\boldsymbol{\xi}}):=B_3(\xi_0) B_3(\xi_2) N_4(\xi_8) + B_5(\xi_3) B_5(\xi_4) N_2(\xi_1) N_2(\xi_6) + B_3(\xi_7) N_2(\xi_5) N_2(\xi_9)
\vspace{-0.5em}
\end{equation}
In spatial dimensions $j=0,2,3,4,7$, we use Chebyshev basis functions as well as Fourier basis functions in the remaining spatial dimensions $j=1,5,6,8,9$.
All tests are performed 10 times and the relative $L^2(\tilde{\D},\mu_\mathrm{F,C})$ approximation error
$\Vert f-a\Vert_{L^2(\tilde{\D},\mu_\mathrm{F,C})} / \Vert f\Vert_{L^2(\tilde{\D},\mu_\mathrm{F,C})}$
is computed each time,
where the approximant $a:=\sum_{{\boldsymbol{n}}\in \tilde{\Omega}} a_{{\boldsymbol{n}}} \,T_{\boldsymbol{n}}$, $$T_{\boldsymbol{n}}({\boldsymbol{\xi}}):=\left(\prod_{j\in\{0,2,3,4,7\}} \cos(n_j\arccos{\xi_j})\right)\,\left(\prod_{j\in\{1,5,6,8,9\}} \mathrm{e}^{2\pi\mathrm{i}n_j\xi_j}\right)$$ is the mixed product basis, and $\mu_\mathrm{F,C}({\boldsymbol{\xi}}):= \left(\frac{2}{\pi}\right)^5 \prod_{j\in \{0,2,3,4,7\}} (1-\xi_j^2)^{-1/2}$ is the corresponding probability measure.
\end{sloppypar}

Here we set the parameters $N=64$, $d=D=10$, $m_2=4s$, and we always use $m_{\rm CE}:=50\,s$ samples for the coefficient estimation where $s=|\S|$.
We consider two different parameter combinations from the previous subsection: $m_1=4s$ and $\kappa=20$, as well as $m_1=8s$ and $\kappa=10$. The obtained results, i.e., the numbers of samples, average runtimes, average iterations, and relative $L^2(\tilde{\D},\mu_\mathrm{F,C})$ errors are plotted as a function of the sparsity~$s\in\{25,50,100,200\}$ in Figure~\ref{fig:approx_sparse:mixed_d10:samples_runtime_iter_success_vs_s}. As before, the numbers of samples grow quadratically for increasing sparsity~$s$. The average runtimes grow approximately like $\sim s^3\min\{s,N\}$ and this means $\sim s^3$ for fixed $N$.
Moreover, the average numbers of iterations are well below its imposed maximum $\kappa$. The relative $L^2(\tilde{\D},\mu_\mathrm{F,C})$ errors decrease for increasing sparsity $s$ having a value of approximately $4.9\cdot 10^{-3}$ for sparsity $s=500$.
We emphasize that Algorithm~\ref{alg:main} is able to easily determine the $s=500$ approximately largest basis coefficients and the corresponding indices for our test function out of $| \mathcal{I}_{N,d}| = N^d = 64^{10} \approx 10^{18}$ possible indices.

\begin{figure}[!h]
	\subfloat[number of samples vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.47\textwidth,
		xtick={25,50,100,200,500},
		xmode=log,
		ymode=log,
		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={\#samples},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,1.912e+05) (50,7.625e+05) (100,3.045e+06) (200,1.217e+07) (500,7.602e+07)
		};
		\addlegendentry{$m_1=4s$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,3.812e+05) (50,1.522e+06) (100,6.085e+06) (200,2.433e+07) (500,1.520e+08)
		};
		\addlegendentry{$m_1=8s$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[average runtime vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={1,100,10000},

		ymin=1,
		ymax=9e5,
		ytick={2,25,400,2000,25000},
		xmode=log,
		ymode=log,
		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. runtime in sec.},
		y label style={xshift=-0.8em},
		legend style={legend cell align=left}, legend pos=north west,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,2.132e+00) (50,2.788e+01) (100,3.105e+02) (200,1.449e+03) (500,2.258e+04)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,3.173e+00) (50,4.525e+01) (100,4.079e+02) (200,1.874e+03) (500,3.065e+04)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\\
	%
	\subfloat[average iteration vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={3,5,7,9},

		ymin=2,ymax=9,
		xmode=log,

		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. iteration},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,5.800e+00) (50,8.100e+00) (100,8.200e+00) (200,8.400e+00) (500,9.200e+00)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,4.200e+00) (50,5.200e+00) (100,5.600e+00) (200,5.600e+00) (500,6.200e+00)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[relative $L^2(\tilde{\D},\mu_\mathrm{F,C})$ error vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},

		ymin=5e-4,ymax=1,
		xmode=log,
		ymode=log,

		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={relative error},
		legend style={legend cell align=left}, legend pos=south west,
		legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,1.871e-01) (50,8.152e-02) (100,4.288e-02) (200,1.689e-02) (500,4.923e-03)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,1.750e-01) (50,8.123e-02) (100,4.275e-02) (200,1.683e-02) (500,4.914e-03)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	\caption{Number of samples, runtime, number of iterations, relative $L^2(\tilde{\D},\mu_\mathrm{F,C})$ error
		vs.\ sparsity $s\in\{25,50,100,200\}$ for mixed Fourier+Chebyshev basis and test function~\eqref{equ:f_mixed:10}.}
	\label{fig:approx_sparse:mixed_d10:samples_runtime_iter_success_vs_s}
\end{figure}

Additionally, we use the preconditioned Legendre polynomials $Q_n$ from Section~\ref{sec:Numerics:exactly_sparse:mixed1} as basis functions in the spatial dimensions $j=0,2,3,4,7$ instead of the Chebyshev polynomials. Besides that, we keep all parameters identical but determine now the relative $L^2(\tilde{\D},\mu_\mathrm{F,L})$ approximation error with respect to the probability measure $\mu_{F,L}({\boldsymbol{\xi}})\equiv 2^{-5}$ which corresponds to the current choice of bases. The results are presented in Figure~\ref{fig:approx_sparse:mixed_l_d10:samples_runtime_iter_success_vs_s}. As before, we observe that the numbers of iterations are higher, now by up to $\approx$ 100\% compared to using Chebyshev polynomials in Figure~\ref{fig:approx_sparse:mixed_l_d10:samples_runtime_iter_success_vs_s}. Correspondingly, the runtimes also double in some cases. The obtained relative errors are similar, but we again remark that we cannot compare these errors directly since they are measured with respect to different probability measures.

\begin{figure}[!h]
	\subfloat[number of samples vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.47\textwidth,
		xtick={25,50,100,200,500},
		xmode=log,
		ymode=log,
		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={\#samples},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,1.912e+05) (50,7.625e+05) (100,3.045e+06) (200,1.217e+07) (500,7.602e+07)
		};
		\addlegendentry{$m_1=4s$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,3.812e+05) (50,1.522e+06) (100,6.085e+06) (200,2.433e+07) (500,1.520e+08)
		};
		\addlegendentry{$m_1=8s$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[average runtime vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},

		xmax=500,
		ymin=4,
		ymax=8e5,
		ytick={4,40,500,2500,4e4},
		xmode=log,
		ymode=log,
		log ticks with fixed point,

		xlabel={$s$},
		ylabel={avg. runtime in sec.},
		y label style={xshift=-0.8em},
		legend style={legend cell align=left}, legend pos=north west,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,4.505e+00) (50,3.624e+01) (100,5.236e+02) (200,2.404e+03) (500,3.490e+04)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,4.725e+00) (50,5.789e+01) (100,5.122e+02) (200,2.905e+03) (500,4.061e+04)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\\
	%
	\subfloat[average iteration vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},
		ytick={5,10,15,20},

		ymin=-1,ymax=15,
		xmode=log,

		log ticks with fixed point,
		xlabel={$s$},
		ylabel={avg. iteration},
		legend style={legend cell align=left}, legend pos=south east,
		legend columns = 1,

		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,8.700e+00) (50,1.010e+01) (100,1.450e+01) (200,1.320e+01) (500,1.380e+01)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
        \addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,5.300e+00) (50,6.700e+00) (100,7.400e+00) (200,8.200e+00) (500,8.100e+00)
        };
        \addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	%
	\hfill
	%
	\subfloat[relative $L^2(\tilde{\D},\mu_\mathrm{F,L})$ error vs.\ sparsity $s$]{
		\begin{tikzpicture}[baseline]
		\begin{axis}[
		font=\footnotesize,
		enlarge x limits=true,
		enlarge y limits=true,
		height=0.26\textwidth,
		grid=major,

		width=0.48\textwidth,
		xtick={25,50,100,200,500},

		ymin=5e-4,ymax=1,
		xmode=log,
		ymode=log,

		xticklabel={
			\pgfkeys{/pgf/fpu=true}
			\pgfmathparse{exp(\tick)}%
			\pgfmathprintnumber[fixed relative, precision=3]{\pgfmathresult}
			\pgfkeys{/pgf/fpu=false}
		},
		xlabel={$s$},
		ylabel={relative error},
		legend style={legend cell align=left}, legend pos=south west,
		legend columns = 1,
		]
		\addplot[black,mark=square,mark size=2.5pt,mark options={solid}] coordinates {
(25,1.639e-01) (50,8.973e-02) (100,4.164e-02) (200,1.619e-02) (500,4.675e-03)
		};
		\addlegendentry{$m_1=4s$, $\kappa=20$}
		\addplot[blue,mark=triangle,mark size=2.5pt,mark options={solid}] coordinates {
(25,1.608e-01) (50,7.551e-02) (100,3.873e-02) (200,1.599e-02) (500,4.613e-03)
		};
		\addlegendentry{$m_1=8s$, $\kappa=10$}
		\end{axis}
		\end{tikzpicture}
	}
	\caption{Number of samples, runtime, number of iterations, relative $L^2(\tilde{\D},\mu_\mathrm{F,L})$ error
		vs.\ sparsity $s\in\{25,50,100,200\}$ for mixed Fourier+Legendre basis and test function~\eqref{equ:f_mixed:10}.}
	\label{fig:approx_sparse:mixed_l_d10:samples_runtime_iter_success_vs_s}
\end{figure}

\FloatBarrier

\subsection{Proof of Theorem~\ref{thm:SuppIDWorks}:  Identifying $\Omega^{\alpha,2s}_{[D]}$ for $\tilde{h}^{\rm opt}_{2s}$ Using Samples from $h = \tilde{h}^{\rm opt}_{2s} + e_h$}
\label{sec:SuppIDHeavyEls}
\input{Support_ID_sec4_1.tex}

\subsection{Proof of Theorem~\ref{thm:ExistFastFuncs}: Generalized Entry Identification \& Pairing}
\label{sec:proofBigSuppIDThm}

In the vast majority of this subsection we will be considering an arbitrary function $\tilde{h}\colon \D \rightarrow \C$ of $D$ variables as per \eqref{def:htilde} whose coefficient vector ${\boldsymbol{\tilde{r}}} \in \C^{\mathcal{I}_{N,d}}$ is only nonzero for entries indexed by index vectors ${\boldsymbol{q}} \in \mathcal{I}_{N,d}$.  In particular, we will be focussing almost exclusively on the development of efficient strategies for learning about the support of the coefficient vector ${\boldsymbol{\tilde{r}}}$ of such $\tilde{h}$ in the special case where ${\boldsymbol{\tilde{r}}}$ is $s'$-sparse so that $\tilde{h} = \tilde{h}^{\rm opt}_{s'}$ and ${\boldsymbol{\tilde{r}}} = {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}}$.
Our first lemma does this by telling us how to estimate the $\ell^2$-norm of any ${\boldsymbol{r''}} = \left({\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}}\right)_{\S;{\boldsymbol{n}}} = {\boldsymbol{\tilde{r}}}_{\left\{ {\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'} ~\big| ~ {\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} \right\}} \in \C^{\mathcal{I}_{N,d}}$ in that case
(i.e., how to estimate of the energy of all the coefficients of $\tilde{h}^{\rm opt}_{s'}$ whose index vectors ${\boldsymbol{q}} \in \mathcal{I}_{N,d}$ match another fixed index vector ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$ in all index positions $\S \subset [D]$) by using just a few inner products with ``simpler'' functions of only $|\S| < D$ variables.  The idea is that these inner products will be easy to approximate numerically for $|\S|$ small.  As a result, one can hope to learn about the index vectors of the nonzero entries of any such ${\boldsymbol{r''}}$ by approximately computing just a few inner products involving functions of just a few variables in order to, e.g., discover values of ${\boldsymbol{n}}$ for which $\| {\boldsymbol{r''}} \|_2$ is large.\\

\begin{lemma}
Let $\delta \in (0, 3/4]$, $\S \subset [D]$, and $\{{\boldsymbol{z}}_{k}\}_{k \in [m_2]} \subset \D_{\S^c}$ be $m_2$ sampling points drawn independently at random according to $\mu_{\S^c}$ in order to form a zero-padded random sampling matrix $\Phi_{\S^c;{\boldsymbol{0}}} \in\C^{m_2\times |\mathcal{I}_{N,d}|}$ for the BOS $\B_{\S^c}$ as in \eqref{def:B_S}  with entries
\begin{equation}
\left( \Phi_{\S^c;{\boldsymbol{0}}} \right)_{k,{\boldsymbol{q}}} :=  \begin{cases}
       T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) & \textrm{if}~{\boldsymbol{q}}_{\S} = {\boldsymbol{0}} ~\&~ {\boldsymbol{q}} \in \mathcal{I}_{N,d}\\%\textrm{if}~{\boldsymbol{q}}_{\S^c} = {\boldsymbol{0}}_{\S^c} ~\&~ {\boldsymbol{q}} \in \mathcal{I}_{N,d}\\
       0 &\text{ otherwise}\\
   \end{cases}
\label{equ:EI_RIPmat}
\end{equation}
indexed by $k \in [m_2]$ and ${\boldsymbol{q}} \in \mathcal{I}_{N,d}$.   Suppose the nonzero columns of $\frac{1}{\sqrt{m_2}}\Phi_{\S^c;{\boldsymbol{0}}}$ have the restricted isometry property (RIP) of order $(s,\delta)$.  Then, for all ${\boldsymbol{n}} \in \mathcal{I}_{N,d} \subseteq [N]^D$, vectors of additive evaluation errors ${\boldsymbol{e}}^h \in \C^{m_2}$, and functions $\tilde{h}$ as per \eqref{def:htilde} one will have
\begin{equation}
\left| \sqrt{\sum_{k\in [m_2]} \frac{1}{m_2} \left| \left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)} + e^h_k \right|^2 } - \| {\boldsymbol{r''}} \|_2 \right|
\leq
\frac{2}{3}
\delta \| {\boldsymbol{r''}} \|_2 +\frac{\| {\boldsymbol{e}}^h \|_2}{\sqrt{m_2}}
\label{equ:EI_error_guar}
\end{equation}
where ${\boldsymbol{r''}} := {\boldsymbol{\tilde{r}}}_{\left\{ {\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'} ~\big| ~ {\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} \right\}}$.
\label{lem:EntryID_ExactInnerProd}
\end{lemma}

\begin{proof}
Consider the zero-padded random sampling matrix $\Phi_{\S^c;{\boldsymbol{n}}} \in\C^{m_2\times |\mathcal{I}_{N,d}|}$ for the BOS $\B_{\S^c}$ as in \eqref{def:B_S} with entries
\begin{equation}
\left( \Phi_{\S^c;{\boldsymbol{n}}} \right)_{k,{\boldsymbol{q}}} :=  \begin{cases}
       T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) & \textrm{if}~{\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} ~\&~ {\boldsymbol{q}} \in \mathcal{I}_{N,d}\\
       0 &\text{ otherwise}\\
   \end{cases}
\label{equ:EI_RIPmatn}
\end{equation}
indexed by $k \in [m_2]$ and ${\boldsymbol{q}} \in \mathcal{I}_{N,d}$.  Note that ${\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} ~\&~ {\boldsymbol{q}} \in \mathcal{I}_{N,d} \implies ({\boldsymbol{q}},{\boldsymbol{0}})_{\S^c} \in \mathcal{I}_{N,d}$ for all $d \in [D+1] \setminus \{0\}, \S^c \subset [D],$ and ${\boldsymbol{n}},{\boldsymbol{q}} \in \mathcal{I}_{N,d}$. As a result, the nonzero columns of $\Phi_{\S^c;{\boldsymbol{0}}}$ will contain the nonzero columns of $\Phi_{\S^c;{\boldsymbol{n}}}$ as a subset.\footnote{Note that the nonzero columns of $\Phi_{\S^c;{\boldsymbol{n}}}$ will be indexed by different ${\boldsymbol{q}}$ in $\Phi_{\S^c;{\boldsymbol{0}}}$.  However, this reindexing will ultimately just represent a permutation of the nonzero columns of $\Phi_{\S^c;{\boldsymbol{n}}}$ as a submatrix of $\Phi_{\S^c;{\boldsymbol{0}}}$.  And, permuting the columns of a matrix does not change its restricted isometry constants.}  This further implies that the matrix consisting of the nonzero columns of $\frac{1}{\sqrt{m_2}}\Phi_{\S^c;{\boldsymbol{n}}}$ will also have the restricted isometry property (RIP) of order $(s,\delta)$.

Applying Lemma~\ref{lem:PartEvalinnerP} together with the definition of $\tilde{h}_{s'}^{\rm opt}$ we now have that
$$\sum_{k \in [m_2]} \frac{1}{m_2} \left| \left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)} + e^h_k \right|^2 =\sum_{k \in [m_2]}  \left| \frac{1}{\sqrt{m_2}} \left \langle~ \left({\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}}\right)_{\S;{\boldsymbol{n}}},  \overline{\Phi_{\S^c;{\boldsymbol{n}};{\boldsymbol{z}}_k}} ~\right \rangle +   \frac{e^h_k}{\sqrt{m_2}} \right|^2 .$$
Noting now that each vector $\Phi_{\S^c;{\boldsymbol{n}};{\boldsymbol{z}}_k}$ as per \eqref{equ:RowRIPlem1} can be replaced by an equivalent row of $\Phi_{\S^c;{\boldsymbol{n}}}$ in \eqref{equ:EI_RIPmatn} we can further see that
\begin{align*}
\sqrt{\sum_{k \in [m_2]} \frac{1}{m_2} \left| \left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)} + e^h_k \right|^2 } &= \sqrt{\sum_{k \in [m_2]}  \left| \frac{1}{\sqrt{m_2}} \left(\Phi_{\S^c;{\boldsymbol{n}}} ({\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}})_{\S;{\boldsymbol{n}}} \right)_k +  \frac{e^h_k}{\sqrt{m_2}} \right|^2 }\\
&= \left \| \frac{1}{\sqrt{m_2}} \Phi_{\S^c;{\boldsymbol{n}}} ({\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}})_{\S;{\boldsymbol{n}}} + \frac{{\boldsymbol{e}}^h}{\sqrt{m_2}} \right\|_2\\
&= \left \| \frac{1}{\sqrt{m_2}} \Phi_{\S^c;{\boldsymbol{n}}} \, {\boldsymbol{r''}} + \frac{{\boldsymbol{e}}^h}{\sqrt{m_2}} \right\|_2.
\end{align*}
Using that $\Phi_{\S^c;{\boldsymbol{n}}}$ has the restricted isometry property (RIP) of order $(s,\delta)$ together with the (reverse) triangle inequality on this last line we now can see that
\begin{align*}
\sqrt{1 - \delta} \| {\boldsymbol{r''}} \|_2 - \frac{\| {\boldsymbol{e}}^h \|_2}{\sqrt{m_2}} &\leq \sqrt{\sum_{k \in [m_2]} \frac{1}{m_2} \left| \left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)} + e^h_k \right|^2 } \\ &\leq \sqrt{1 + \delta} \| {\boldsymbol{r''}} \|_2 +\frac{\| {\boldsymbol{e}}^h \|_2}{\sqrt{m_2}}.
\end{align*}
After subtracting $\| {\boldsymbol{r''}} \|_2$ from the quantities in the inequality above, we finally use the bounds
$(\sqrt{1 + \delta} - 1) \frac{\sqrt{1 + \delta} + 1}{\sqrt{1 + \delta} + 1} = \frac{\delta}{\sqrt{1 + \delta} + 1} < \frac{\delta}{2}$ and $(\sqrt{1 - \delta} - 1) \frac{\sqrt{1 - \delta} + 1}{\sqrt{1 - \delta} + 1} = \frac{-\delta}{\sqrt{1 - \delta} + 1} \geq -\frac{2}{3} \delta$ to finish the proof.
\end{proof}

Lemma~\ref{lem:EntryID_ExactInnerProd} yields an alternate entry identification technique to that provided in Section 4.1 of \cite{choi2018sparse}.  In particular, if $\S^c = [D] \setminus \{ j \}$ for some $j \in [D]$ the inner products $\left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)}$ are just one-dimensional integrals that can be computed to high accuracy for any desired ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$ using only $\mathcal{O}(N)$ function evaluations of $(\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}: \D_{\S} \rightarrow \C$ via, e.g., a quadrature rule whenever the the basis functions in the $j^{\rm th}$-dimension, $\B_j = \left \{ T_{\S;{\boldsymbol{n}}} ~|~ {\boldsymbol{n}} \in \mathcal{I}_{N,d} \right \}$, are polynomials of degree at most $N$.  If $\B_j$ is either the Fourier or Chebyshev basis and $N$ is very large then these one-dimensional integrals can also be computed for all ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$ in sublinear-in-$N$ time since $(\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}$ will be $\B_j$-sparse (see, e.g., \cite{gilbert2014recent,gilbert2005improved,iwen2007empirical,iwen2008deterministic,iwen2010combinatorial,bailey2012design,hassanieh2012simple,iwen2013improved,segal2013improved,merhi2017new,hu2017rapidly,bittens2019deterministic}).

When $|\S^c| \ll D-1$ the situation becomes more difficult.  However, to efficiently evaluate the higher-dimensional inner products $\left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)}$ that arise in these settings one can instead utilize non-adaptive random sampling techniques motivated by compressive sensing theory.  The following lemma does this by quantifying how well the estimator $$\frac{1}{m_1} \sum_{\ell \in [m_1]} (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k} ({\boldsymbol{w}}_\ell) ~\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)}$$ based on the $m_1$ randomly chosen grid points $\left\{ {\boldsymbol{w}}_\ell \right\}_{\ell \in [m_1]} \subset \D_{\S}$ approximates all such $\left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)}$.\\

\begin{lemma}
Let $\tilde{\delta} \in (0,1)$, $\S \subset [D]$, and $\{{\boldsymbol{w}}_{\ell}\}_{\ell \in [m_1]} \subset \D_{\S}$ be $m_1$ sampling points drawn independently at random according to $\mu_{\S}$ in order to form a zero-padded random sampling matrix $\Phi_{\S;{\boldsymbol{0}}} \in\C^{m_1 \times |\mathcal{I}_{N,d}|}$ for the BOS $\B_{\S}$  in \eqref{def:B_S} with entries
\begin{equation}
\left( \Phi_{\S;{\boldsymbol{0}}} \right)_{\ell,{\boldsymbol{q}}} :=  \begin{cases}
       T_{\S;{\boldsymbol{q}}}({\boldsymbol{w}}_\ell) & \textrm{if}~{\boldsymbol{q}}_{\S^c} = {\boldsymbol{0}} ~\&~ {\boldsymbol{q}} \in \mathcal{I}_{N,d}\\
       0 &\text{ otherwise}\\
   \end{cases}
\label{equ:EI_RIPmat2}
\end{equation}
indexed by $\ell \in [m_1]$ and ${\boldsymbol{q}} \in \mathcal{I}_{N,d}$.   Suppose the nonzero columns of $\frac{1}{\sqrt{m_1}}\Phi_{\S;{\boldsymbol{0}}}$ have the restricted isometry property (RIP) of order $(2,\tilde{\delta})$, and let
$$e^h_k := \frac{1}{m_1} \sum_{\ell \in [m_1]} (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k} ({\boldsymbol{w}}_\ell) ~\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)} - \left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)}$$
for any desired ${\boldsymbol{n}} \in \mathcal{I}_{N,d} \subseteq [N]^D$, function $\tilde{h}$ as per \eqref{def:htilde}, and point ${\boldsymbol{z}}_k \in \D_{\S^c}$.
Then,
\begin{equation}
\left| e^h_k \right| \leq \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2 \tilde{\delta} \left \| \left(\Phi_{\S^c;{\boldsymbol{0}};{\boldsymbol{z}}_k}\right)_{\Omega^{\rm opt}_{\tilde{h},s'}} \right \|_2
\end{equation}
where $\Phi_{\S^c;{\boldsymbol{0}};{\boldsymbol{z}}_k}$ is defined as in \eqref{equ:RowRIPlem1}.
\label{lem:SupportLem2}
\end{lemma}

\begin{proof}
We begin by noting that
\begin{align}
\frac{1}{m_1} \sum_{\ell \in [m_1]} (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k} ({\boldsymbol{w}}_\ell) ~\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)} ~=~ \sum_{{\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'}} \tilde{r}_{\boldsymbol{q}} \, T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) ~\nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right)
\label{equ:ProofNewLem7_1}
\end{align}
where
\begin{align*}
\nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right) &:= \sum_{\ell \in [m_1]} \frac{1}{\sqrt{m_1}}T_{\S;{\boldsymbol{q}}}({\boldsymbol{w}}_\ell) ~\frac{1}{\sqrt{m_1}}\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)}\\
&= \left \langle \frac{1}{\sqrt{m_1}} \left( \Phi_{\S;{\boldsymbol{0}}} \right)_{({\boldsymbol{q}},{\boldsymbol{0}})_{\S}}, \frac{1}{\sqrt{m_1}} \left( \Phi_{\S;{\boldsymbol{0}}} \right)_{({\boldsymbol{n}},{\boldsymbol{0}})_{\S}} \right \rangle.
\end{align*}
Appealing to standard results concerning coherence in, e.g., Chapter 6 of \cite{foucart2013mathematical} one can see that
$ \left| \nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right) - 1 \right | \leq \tilde{\delta}$ holds if ${\boldsymbol{n}}_{\S} = {\boldsymbol{q}}_{\S}$, and that $ \left| \nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right) \right | \leq \tilde{\delta}$ holds if ${\boldsymbol{n}}_{\S} \neq {\boldsymbol{q}}_{\S}$.

Let $\Omega' := \left\{ {\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'} ~|~ {\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} \right \} \subset \Omega^{\rm opt}_{\tilde{h},s'}$, and $\Omega'' := \Omega^{\rm opt}_{\tilde{h},s'} \setminus \Omega'$.
Using \eqref{equ:ProofNewLem7_1} one has that
\begin{align*}
\Bigg| \frac{1}{m_1} \sum_{\ell \in [m_1]} &(\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k} ({\boldsymbol{w}}_\ell) ~\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)} - \sum_{{\boldsymbol{q}} \in \Omega'} \tilde{r}_{\boldsymbol{q}} T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) \Bigg| \\ &= \left| \sum_{{\boldsymbol{q}} \in \Omega'} \tilde{r}_{\boldsymbol{q}} T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) \left( \nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right) - 1 \right) + \sum_{{\boldsymbol{q}} \in \Omega''} \tilde{r}_{\boldsymbol{q}} T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) \nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right) \right|\\
&\leq \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2 \sqrt{ \sum_{{\boldsymbol{q}} \in \Omega'} \left | T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) \left( \nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right) - 1 \right) \right |^2 + \sum_{{\boldsymbol{q}} \in \Omega''} \left | T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) ~\nu_{\S}\left( {\boldsymbol{q}},{\boldsymbol{n}} \right) \right |^2}
\end{align*}
where the last inequality follows from Cauchy-Schwarz.  Continuing from this the last line we can further see that
\begin{align*}
\Bigg| \frac{1}{m_1} \sum_{\ell \in [m_1]} &(\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k} ({\boldsymbol{w}}_\ell) ~\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)} - \sum_{{\boldsymbol{q}} \in \Omega'} \tilde{r}_{\boldsymbol{q}} T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) \Bigg| \\
&\leq \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2 \tilde{\delta} \sqrt{ \sum_{{\boldsymbol{q}} \in \Omega'} \left | T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) \right |^2 + \sum_{{\boldsymbol{q}} \in \Omega''} \left | T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) \right |^2} \\
&= \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2 \tilde{\delta} \left \| \left(\Phi_{\S^c;{\boldsymbol{0}};{\boldsymbol{z}}_k}\right)_{\Omega^{\rm opt}_{\tilde{h},s'}} \right \|_2.
\end{align*}

To finish we note that
\begin{equation*}
\sum_{{\boldsymbol{q}} \in \Omega'} \tilde{r}_{\boldsymbol{q}} T_{\S^c;{\boldsymbol{q}}}({\boldsymbol{z}}_k) ~= ~\left \langle ~ \left({\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}}\right)_{\S;{\boldsymbol{n}}}, \overline{\Phi_{\S^c;{\boldsymbol{n}};{\boldsymbol{z}}_k}} ~\right \rangle ~=~ \left \langle \left(\tilde{h}_{s'}^{\rm opt}\right)_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)}
\end{equation*}
by Lemma~\ref{lem:PartEvalinnerP}.  The desired result follows.
\end{proof}

Choose any ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$ and $\S \subset [D]$ you like.  Using Lemma~\ref{lem:SupportLem2} to approximate the inner product appearing in Lemma~\ref{lem:EntryID_ExactInnerProd}'s \eqref{equ:EI_error_guar} then yields the following estimator for accurately approximating
the $\ell^2$-norm of ${\boldsymbol{r''}} := {\boldsymbol{\tilde{r}}}_{\left\{ {\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'} ~\big| ~ {\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} \right\}} \in \C^{\mathcal{I}_{N,d}}$ for the coefficient vector ${\boldsymbol{\tilde{r}}}$ of any function $\tilde{h}$ as in \eqref{def:htilde}.  The estimator is defined for any function $u: \D \rightarrow \C$, $\S \subset [D]$, and ${\boldsymbol{n}} \in \mathcal{I}_{N,d}$ to be
\begin{equation}
E^u_{\S; {\boldsymbol{n}}} := \frac{1}{m_2} \sum_{k\in [m_2]} \left| \frac{1}{m_1} \sum_{\ell \in [m_1]} u_{\S^c;{\boldsymbol{z}}_k} ({\boldsymbol{w}}_\ell) ~\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)}  \right|^2
\label{equ:EstimatorGen}
\end{equation}
for fixed nodes $\{ {\boldsymbol{w}}_\ell \}_{\ell \in [m_1]} \subset \D_{\S}$, and $\{ {\boldsymbol{z}}_k \}_{k \in [m_2]} \subset \D_{\S^c}$.  Note that \eqref{equ:EstimatorGen} is essentially identical to the pairing energy estimator defined in Section 4.2 of \cite{choi2018sparse}.  The following lemma provides an error guarantee for this estimator that matches the quality of those in \cite{choi2018sparse} despite having a simpler proof (see Lemma 7 in \cite{choi2018sparse}).\\

\begin{lemma}
Let $\S \subset [D]$, $\delta \in (0, 3/4]$, $\tilde{\delta} \in (0,1/s']$, and $m_1, m_2 \in \N$.  Furthermore, suppose that $\{ {\boldsymbol{w}}_\ell \}_{\ell \in [m_1]} \subset \D_{\S}$, and $\{ {\boldsymbol{z}}_k \}_{k \in [m_2]} \subset \D_{\S^c}$ satisfy the RIP assumptions concerning \eqref{equ:EI_RIPmat2} and \eqref{equ:EI_RIPmat} in Lemmas~\ref{lem:SupportLem2} and~\ref{lem:EntryID_ExactInnerProd}, respectively.
Then, for all ${\boldsymbol{n}} \in \mathcal{I}_{N,d} \subseteq [N]^D$ and functions $\tilde{h}$ as per \eqref{def:htilde} one will have
$$\left| \sqrt{E^{\tilde{h}_{s'}^{\rm opt}}_{\S; {\boldsymbol{n}}}} - \| {\boldsymbol{r''}} \|_2 \right|
\leq
\frac{2}{3}
\delta \| {\boldsymbol{r''}} \|_2 +
\sqrt{\frac{7}{4}} \sqrt{\tilde{\delta}} \cdot \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2, $$
where ${\boldsymbol{r''}} := {\boldsymbol{\tilde{r}}}_{\left\{ {\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'} ~\big| ~ {\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} \right\}}$.
\label{lem:EstErrorforEstimator}
\end{lemma}

\begin{proof}
Applying Lemma~\ref{lem:EntryID_ExactInnerProd} we can immediately see that
\begin{equation}
\left| \sqrt{E^{\tilde{h}_{s'}^{\rm opt}}_{\S; {\boldsymbol{n}}}} - \| {\boldsymbol{r''}} \|_2 \right| \leq
\frac{2}{3}
\delta
\| {\boldsymbol{r''}} \|_2 +\frac{\| {\boldsymbol{e}}^h \|_2}{\sqrt{m_2}}
\label{equ:FinalNewlem7b}
\end{equation}
where ${\boldsymbol{e}}^h \in \C^{m_2}$ has its entries given by
$$e^h_k := \frac{1}{m_1} \sum_{\ell \in [m_1]} (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k} ({\boldsymbol{w}}_\ell) ~\overline{T_{\S;{\boldsymbol{n}}}({\boldsymbol{w}}_\ell)} - \left \langle (\tilde{h}_{s'}^{\rm opt})_{\S^c;{\boldsymbol{z}}_k}, T_{\S;{\boldsymbol{n}}} \right \rangle_{\left( \D_{\S}, \mu_{\S}\right)}.$$
Thus, it suffices to bound $\| {\boldsymbol{e}}^h \|_2$ in order to obtain our final result.

Applying Lemma~\ref{lem:SupportLem2} we can see that
\begin{align*}
\| {\boldsymbol{e}}^h \|_2^2 &\leq \sum_{k\in [m_2]} \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2^2 \tilde{\delta}^2 \left \| \left(\Phi_{\S^c;{\boldsymbol{0}};{\boldsymbol{z}}_k}\right)_{\Omega^{\rm opt}_{\tilde{h},s'}} \right \|_2^2\\
&= \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2^2 \tilde{\delta}^2 \left \| \left(\Phi_{\S^c;{\boldsymbol{0}}}\right)_{\Omega^{\rm opt}_{\tilde{h},s'}} \right \|_{\rm F}^2,
\end{align*}
where we have used that $\left(\Phi_{\S^c;{\boldsymbol{0}};{\boldsymbol{z}}_k}\right)_{\Omega^{\rm opt}_{\tilde{h},s'}}$ are the rows of the submatrix $\left( \Phi_{\S^c;{\boldsymbol{0}}} \right)_{\Omega^{\rm opt}_{\tilde{h},s'}} \in \C^{m_2 \times s'}$ of $\Phi_{\S^c;{\boldsymbol{0}}}$ in \eqref{equ:EI_RIPmat}.  Using the RIP property of the nonzero columns of $\frac{1}{\sqrt{m_2}} \Phi_{\S^c;{\boldsymbol{0}}}$ we can now finish bounding $\| {\boldsymbol{e}}^h \|_2^2$ by noting that

\begin{align*}
\| {\boldsymbol{e}}^h \|_2^2 ~\leq~ \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2^2 \tilde{\delta}^2 \cdot s'\cdot m_2 (1 + \delta) ~\leq~
\frac{7}{4}
 m_2 \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2^2 \tilde{\delta}
\end{align*}
where we have used that $\delta \in (0, 3/4]$ and that $\tilde{\delta} \in (0,1/s']$.  Substituting this last bound into \eqref{equ:FinalNewlem7b} now finishes the proof.
\end{proof}

Though useful, Lemma~\ref{lem:EstErrorforEstimator} presupposes that one has access to noiseless samples from $\tilde{h}_{s'}^{\rm opt}$.  This will rarely be the case in practice.  The next lemma bounds the error of the estimator \eqref{equ:EstimatorGen} in the setting where one instead has noisy samples from $\tilde{h}_{s'}^{\rm opt}$.  Such noisy samples will be represented with the help of an arbitrary additive noise/error function, $e_h: \D \rightarrow \C$, below.\\

\begin{lemma}
Let $\S \subset [D]$, $\delta \in (0, 3/4]$, $\tilde{\delta} \in (0,1/s']$, and $m_1, m_2 \in \N$.  Furthermore, suppose that $\{ {\boldsymbol{w}}_\ell \}_{\ell \in [m_1]} \subset \D_{\S}$, and $\{ {\boldsymbol{z}}_k \}_{k \in [m_2]} \subset \D_{\S^c}$ satisfy the RIP assumptions concerning \eqref{equ:EI_RIPmat2} and \eqref{equ:EI_RIPmat} in Lemmas~\ref{lem:SupportLem2} and~\ref{lem:EntryID_ExactInnerProd}, respectively.
Then, for all ${\boldsymbol{n}} \in \mathcal{I}_{N,d} \subseteq [N]^D$, $\tilde{h}$ as per \eqref{def:htilde}, and additive error functions $e_h: \D \rightarrow \C$ one will have
$$\left| \sqrt{E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}}} - \| {\boldsymbol{r''}} \|_2 \right| \leq
\frac{2}{3}
\delta \| {\boldsymbol{r''}} \|_2 +
\sqrt{\frac{7}{4} \tilde{\delta}} \cdot \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2 + \frac{
\sqrt{2}\| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1m_2}}, $$
where ${\boldsymbol{r''}} := {\boldsymbol{\tilde{r}}}_{\left\{ {\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'} ~\big| ~ {\boldsymbol{q}}_{\S} = {\boldsymbol{n}}_{\S} \right\}}$ and $\mathcal{E}^h_{\S} \in \C^{m_1 \times m_2}$ has entries $(\mathcal{E}^h_{\S})_{\ell, k}=e_h(\varrho_{\mathcal{S}}({\boldsymbol{w}}_\ell, {\boldsymbol{z}}_k))$ with the permutation function
$\varrho_\S\colon \D_{\S}\times\D_{\S^c}\rightarrow\D$ defined in
\eqref{def:rho}.
\label{lem:EstErrorforEstimatorNoise}
\end{lemma}

\begin{proof}
Note that
\begin{equation*}
\sqrt{E^{\tilde{h}_{s'}^{\rm opt}}_{\S; {\boldsymbol{n}}}}-\sqrt{E^{e_h}_{\S; {\boldsymbol{n}}}} \leq \sqrt{E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}}} \leq \sqrt{E^{\tilde{h}_{s'}^{\rm opt}}_{\S; {\boldsymbol{n}}}} + \sqrt{E^{e_h}_{\S; {\boldsymbol{n}}}}
\end{equation*}
by the (reverse) triangle inequality.  As a result one can immediately see that
\begin{align*}
\left| \sqrt{E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}}} - \| {\boldsymbol{r''}} \|_2 \right| &\leq \left| \sqrt{E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}}} -  \sqrt{E^{\tilde{h}_{s'}^{\rm opt}}_{\S; {\boldsymbol{n}}}}\right| + \left| \sqrt{E^{\tilde{h}_{s'}^{\rm opt}}_{\S; {\boldsymbol{n}}}} - \| {\boldsymbol{r''}} \|_2 \right|\\
&\leq  \sqrt{E^{e_h}_{\S; {\boldsymbol{n}}}}  +
\frac{2}{3}
\delta \| {\boldsymbol{r''}} \|_2 +  \sqrt{ \frac{7}{4} \tilde{\delta}} \cdot \left\| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2
\end{align*}
where the bound on the second term above follows from Lemma~\ref{lem:EstErrorforEstimator}. It remains to show that $\sqrt{E^{e_h}_{\S; {\boldsymbol{n}}}} \leq \frac{\sqrt{2} \| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1m_2}}$.

Define ${\boldsymbol{v}} \in \C^{m_2}$ by $v_k:=\frac{1}{\sqrt{m_1}} \left\langle (\mathcal{E}^h_{\S})_k,  \frac{1}{\sqrt{m_1}} \left( \Phi_{\S;{\boldsymbol{0}}} \right)_{({\boldsymbol{n}},{\boldsymbol{0}})_{\S}} \right\rangle$ where $\Phi_{\S;{\boldsymbol{0}}}$ is defined in \eqref{equ:EI_RIPmat2}, and note that  $|v_k| \leq  \frac{1}{\sqrt{m_1}} \left\| (\mathcal{E}^h_{\S})_k \right\|_2 \left\|  \frac{1}{\sqrt{m_1}} \left( \Phi_{\S;{\boldsymbol{0}}} \right)_{({\boldsymbol{n}},{\boldsymbol{0}})_{\S}} \right\|_2$. Furthermore, $$\left\|  \frac{1}{\sqrt{m_1}} \left( \Phi_{\S;{\boldsymbol{0}}} \right)_{({\boldsymbol{n}},{\boldsymbol{0}})_{\S}} \right\|_2 = \left\|  \frac{1}{\sqrt{m_1}} \left( \Phi_{\S;{\boldsymbol{0}}} \right)_{\left\{({\boldsymbol{n}},{\boldsymbol{0}})_{\S}, ({\boldsymbol{\tilde{n}}},{\boldsymbol{0}})_{\S} \right\}} \begin{bmatrix} 1\\0 \end{bmatrix} \right\|_2 \leq \sqrt{1+\widetilde{\delta}}
\leq \sqrt{1+\frac{1}{s'}}
\leq \sqrt{2}$$ for any ${\boldsymbol{\tilde{n}}} \neq {\boldsymbol{n}}$. Thus, $|v_k| \leq \sqrt{ \frac{2}{m_1} } \left\| (\mathcal{E}^h_{\S})_k \right\|_2$. As a result, $\sqrt{E^{e_h}_{\S; {\boldsymbol{n}}}} =\frac{1}{\sqrt{m_2}} \| {\boldsymbol{v}} \|_2 \leq \frac{\sqrt{2} \| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1m_2}}$.
\end{proof}

For any given $\S \subseteq [D]$ we denote the power set of $N^{\S}$ by $\mathcal{P}\left(N^{\S} \right)$.  In the final theorem of this subsection we will prove that the energy estimator in \eqref{equ:EstimatorGen} can be used for an arbitrary $s'$-sparse function $\tilde{h}^{\rm opt}_{s'}$ to define a new set-valued function $\F^{s'}_{\S}: \mathcal{P}\left(N^{\S} \right) \rightarrow \mathcal{P}\left(N^{\S} \right)$ for each $\S$ which, when given any subset $\T \subset N^{\S}$ containing the heavy set $\Omega^{\alpha,{s'}}_{\S}$as per \eqref{def:HeavySet} as input, will output a smaller subset $\T' \subset \T$  which still contains $\Omega^{\alpha,{s'}}_{\S}$.  These set-valued functions were also called ``energetic-index sieve function'' in Section~\ref{sec:SuppIDHeavyEls} and will then be used to iteratively build up subsets $\T$ of controlled cardinality for larger and larger sets of indices $\S$ until we eventually have a set of full index vectors $\T'' \subset \mathcal{I}_{N,d}$ which contains all of $\Omega^{\alpha,{s'}}_{[D]}$.  This set of full index vectors $\T''$ will then be able to be used as an accurate estimate of $\Omega^{\rm opt}_{\tilde{h},s'}$, the support of $\tilde{h}^{\rm opt}_{s'}$.

Before we can state our final theorem we must define the set-valued functions $\F^{s'}_{\S}: \mathcal{P}\left(N^{\S} \right) \rightarrow \mathcal{P}\left(N^{\S} \right)$ in question. For a given $\S \subseteq [D]$, $\T \subseteq N^{\S}$, $\tilde{h}$ as per \eqref{def:htilde}, and additive error function $e_h: \D \rightarrow \C$, let an ordering of the elements of $\T$, ${\boldsymbol{n}}_1, {\boldsymbol{n}}_2, \cdots, {\boldsymbol{n}}_{|\T|} \in \T$, be defined by
\begin{equation}
E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}_1} ~\geq~ E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}_2} ~\geq~ E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}_3} ~\geq~ \dots ~\geq~ E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{n}}_{|\T|}}
\label{def:SizeOrder}
\end{equation}
with ties broken lexicographically.
We define $\F^{s'}_{\S}$ based on this ordering by
\begin{equation}
\F^{s'}_{\S} (\T) := \left \{ {\boldsymbol{n}}_1, {\boldsymbol{n}}_2, \cdots, {\boldsymbol{n}}_{\min({s'},|\T|)} \right\} \subseteq \T.
\label{def:SetValCardRedux}
\end{equation}
The following theorem proves that $\Omega^{\alpha,{s'}}_{\S} \cap \T \subseteq \F^{s'}_{\S} (\T)$ provided that the additive error~$e_h$ is sufficiently mild.\\

\begin{mytheorem}[{Entry Identification and Pairing}]
Let $\S \subseteq [D]$ with $|\S| > 0$, $\delta \in (0, 1/2]$, $\tilde{\delta} \in \left(0,\frac{1}{256\alpha^2s'} \right]$, and $m_1, m_2 \in \N$.  Furthermore, suppose that $\{ {\boldsymbol{w}}_\ell \}_{\ell \in [m_1]} \subset \D_{\S}$, and $\{ {\boldsymbol{z}}_k \}_{k \in [m_2]} \subset \D_{\S^c}$ satisfy the RIP assumptions concerning \eqref{equ:EI_RIPmat2} and \eqref{equ:EI_RIPmat} in Lemmas~\ref{lem:SupportLem2} and~\ref{lem:EntryID_ExactInnerProd}, respectively.  Then, $\Omega^{\alpha,s'}_{\S} \cap \T \subseteq \F^{s'}_{\S} (\T)$ for all $s'$-sparse $\tilde{h} = \tilde{h}^{\rm opt}_{s'}$, $\T \subseteq N^{\S}$, and additive error functions $e_h: \D \rightarrow \C$ provided that
\begin{equation}
\| {\boldsymbol{\tilde{r}}} \|_2 = \| {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \|_2 >  \frac{6 \alpha \sqrt{{s'}}}{\sqrt{m_1m_2}} \cdot \| \mathcal{E}^h_{\S} \|_{\rm F} \label{equ:threshold}
\end{equation}
holds, where
$\mathcal{E}^h_{\S} \in \C^{m_1 \times m_2}$ has entries $(\mathcal{E}^h_{\S})_{\ell, k}=e_h(\varrho_\S ({\boldsymbol{w}}_\ell, {\boldsymbol{z}}_k))$ with the permutation function
$\varrho_\S\colon \D_{\S}\times\D_{\S^c}\rightarrow\D$
defined in \eqref{def:rho}.
\label{thm:pairingFunc}
\end{mytheorem}

\begin{proof}
We will focus on the case where $|\T| > s'$ since the result holds trivially when $|\T| \leq s'$.  Suppose for the sake of contradiction that ${\boldsymbol{m}} \in \Omega^{\alpha,s'}_{\S} \cap \T$, but that ${\boldsymbol{m}} \notin \F^{s'}_{\S} (\T)$.  It must then be the case that $E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{k}}} ~\geq~ E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{m}}}$ for some ${\boldsymbol{k}} \in \T$ with ${\boldsymbol{k}}_{\S} \notin \Omega^{\rm opt}_{s',\S} := \left\{ {\boldsymbol{q}}_{\S}  ~\big| ~ {\boldsymbol{q}} \in \Omega^{\rm opt}_{\tilde{h},s'} \right\}$ since $\Omega^{\alpha,s'}_{\S} \subset \Omega^{\rm opt}_{s',\S}$ and $| \F^{s'}_{\S} (\T) | = s' \geq |\Omega^{\rm opt}_{s',\S}|$.  Thus, $\| {\boldsymbol{\tilde{r}}}_{\S;{\boldsymbol{k}}} \|_2 = 0$.  As a result, Lemma~\ref{lem:EstErrorforEstimatorNoise} implies that
\begin{equation}\label{equ:estE_upper_bound_epsilon}
\sqrt{E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{k}}}} ~\leq~ \sqrt{ \frac{7}{4} \tilde{\delta}} \cdot \left\| {\boldsymbol{\tilde{r}}} \right\|_2 + \frac{\sqrt{2} \| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1m_2}} ~\leq~ \frac{\sqrt{7}}{32} \frac{\left\| {\boldsymbol{\tilde{r}}} \right\|_2}{\alpha \sqrt{s'}} + \frac{\sqrt{2} \| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1m_2}} =: \epsilon.
\end{equation}
On the other hand, Lemma~\ref{lem:EstErrorforEstimatorNoise} also implies that $\sqrt{E^{\tilde{h}_{s'}^{\rm opt}+e_h}_{\S; {\boldsymbol{m}}}} \geq \| {\boldsymbol{\tilde{r}}}_{\S;{\boldsymbol{m}}} \|_2 (1 - \frac{2}{3}\delta) - \epsilon$.
Combining this with~\eqref{equ:estE_upper_bound_epsilon}, we have $\frac{2 \epsilon}{1 - \frac{2}{3}\delta} ~\geq~ \| {\boldsymbol{\tilde{r}}}_{\S;{\boldsymbol{m}}} \|_2$.
Since $\delta\leq\frac{1}{2}$, it must also be the case that
$$3 \epsilon ~\geq~ \frac{2 \epsilon}{1 - \frac{2}{3}\delta} ~\geq~ \| {\boldsymbol{\tilde{r}}}_{\S;{\boldsymbol{m}}} \|_2 ~\geq~ \frac{\| {\boldsymbol{\tilde{r}}} \|_2}{\alpha \sqrt{s'}}.$$
However, it is impossible that $4 \alpha \sqrt{s'} \epsilon \geq \| {\boldsymbol{\tilde{r}}} \|_2$ since by assumption
\begin{align*}
3 \alpha \sqrt{s'} \epsilon ~=~
3 \alpha \sqrt{s'} \left( \frac{\sqrt{7}}{32} \frac{\left\| {\boldsymbol{\tilde{r}}} \right\|_2}{\alpha \sqrt{s'}} + \frac{ \sqrt{2}\| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1m_2}} \right)
&< 3 \alpha \sqrt{s'} \left( \frac{\sqrt{7}}{32} \frac{\left\| {\boldsymbol{\tilde{r}}} \right\|_2}{\alpha \sqrt{s'}} +  \frac{ \sqrt{2}\left\| {\boldsymbol{\tilde{r}}} \right\|_2}{6 \alpha \sqrt{s'}} \right) < \left\| {\boldsymbol{\tilde{r}}} \right\|_2.
\end{align*}
Hence, ${\boldsymbol{m}} \in \Omega^{\alpha,s'}_{\S} \cap \T \implies {\boldsymbol{m}} \in \F^{s'}_{\S} (\T)$.
\end{proof}

Theorem~\ref{thm:pairingFunc} forms the basis of our support identification strategy.  As such, it behooves us to investigate its associate resource demands and error performance more closely.  We do this in the next subsection.

\subsubsection{Associated Runtime, Sampling, and Error Bounds}
\label{sec:runtimeerrorsamp}

The following lemmas provide evaluation complexity, sampling, and error bounds for the set valued functions $\F^{s'}_{\S}: \mathcal{P}\left(N^{\S} \right) \rightarrow \mathcal{P}\left(N^{\S} \right)$ defined in \eqref{def:SizeOrder} -- \eqref{def:SetValCardRedux}.  We will begin by providing more meaningful error bounds for the case where the function $\tilde{h}$ in question is not exactly BOPB-sparse.\\

\begin{lemma}
Let $\S \subseteq [D]$ with $|\S| > 0$, $\delta \in (0, 1/2]$, $\tilde{\delta} \in \left(0,\frac{1}{256\alpha^2s'} \right]$, $\gamma \in \mathbbm{R}^+$, and $m_1, m_2 \in \N$.  Furthermore, suppose that $\{ {\boldsymbol{w}}_\ell \}_{\ell \in [m_1]} \subset \D_{\S}$, and $\{ {\boldsymbol{z}}_k \}_{k \in [m_2]} \subset \D_{\S^c}$ satisfy the RIP assumptions concerning \eqref{equ:EI_RIPmat2} and \eqref{equ:EI_RIPmat} in Lemmas~\ref{lem:SupportLem2} and~\ref{lem:EntryID_ExactInnerProd}, respectively.  Finally, suppose also that $e_h = \tilde{h} - \tilde{h}^{\rm opt}_{s'} + e'$ for an arbitrary function $e': \D \rightarrow \C$ with $\displaystyle {\rm sup}_{{\boldsymbol{\xi}} \in \D} |e'({\boldsymbol{\xi}})| \leq \gamma$.  Then, the additive sampling error $\mathcal{E}^h_{\S} \in \C^{m_1 \times m_2}$ satisfies
$$\frac{\| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1 m_2}} ~\leq~
\sqrt{\frac{771}{512}} \left\| {\boldsymbol{\tilde{r}}} - {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2 +
\sqrt{\frac{771}{512\,s'}} \left\| {\boldsymbol{\tilde{r}}} - {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_1 + \gamma$$
where $\mathcal{E}^h_{\S}$ has entries $(\mathcal{E}^h_{\S})_{\ell, k}=e_h\left(\varrho_\S({\boldsymbol{w}}_\ell, {\boldsymbol{z}}_k)\right)$ as in Theorem~\ref{thm:pairingFunc}.
\label{lem:errorBound}
\end{lemma}

\begin{proof}
Note that
\begin{align}
\frac{\| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1 m_2}} &\leq \gamma + \frac{1}{\sqrt{m_1 m_2}} \sqrt{\sum_{\ell,k} \left| \left( \tilde{h} - \tilde{h}^{\rm opt}_{s'} \right) \left(\varrho_{\mathcal{S}}( {\boldsymbol{w}}_\ell, {\boldsymbol{z}}_k)\right) \right|^2} \nonumber \\
&= \gamma  + \frac{1}{\sqrt{m_1 m_2}} \sqrt{\sum_{\ell,k} \left| \left\langle {\boldsymbol{\tilde{r}}} - {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}}, {\Phi}^*_{(\ell,k)}\right\rangle \right|^2} \nonumber\\
&=\left\| \frac{1}{\sqrt{m_1 m_2}} {\Phi} \left( {\boldsymbol{\tilde{r}}} - {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right) \right\|_2 + \gamma \label{equ:energybound},
\end{align}
where ${\Phi} \in \C^{m_1 m_2 \times \left| \mathcal{I}_{N,d} \right|}$ has entries given by ${\Phi}_{(\ell,k),{\boldsymbol{n}}} = T_{{\boldsymbol{n}}}(\varrho_{\mathcal{S}}({\boldsymbol{w}}_\ell, {\boldsymbol{z}}_k))$.  Note also that $\frac{1}{\sqrt{m_1 m_2}} {\Phi}$ consists of a subset of the columns of the Kronecker product $\left( \frac{1}{\sqrt{m_1}} {\Phi}_{\S;{\boldsymbol{0}}} \right) \otimes \left( \frac{1}{\sqrt{m_2}} {\Phi}_{\S^c;{\boldsymbol{0}}} \right)$
where ${\Phi}_{\S;{\boldsymbol{0}}} \in \C^{m_1 \times \left| \mathcal{I}_{N,d} \right|}$ is defined in \eqref{equ:EI_RIPmat2}, and ${\Phi}_{\S^c;{\boldsymbol{0}}} \in \C^{m_2 \times \left|  \mathcal{I}_{N,d} \right|}$ is defined in \eqref{equ:EI_RIPmat}.  Furthermore, Proposition 6.6 of \cite{foucart2013mathematical} implies that the nonzero columns of $\frac{1}{\sqrt{m_1}} {\Phi}_{\S;{\boldsymbol{0}}}$ also has the RIP of order $(s',\frac{1}{256\alpha^2})$ since it has the RIP of order $(2,\frac{1}{256\alpha^2s'})$.  Hence, $\frac{1}{\sqrt{m_1 m_2}} {\Phi}$ has the RIP of order $\left(s',\left(1+\frac{1}{256\alpha^2}\right)\left(1+\frac{1}{2}\right)-1\right)$ by Lemma 2 of \cite{duarte2012kronecker}, and consequently of order $\left(s',\frac{259}{512}\right)$ for $\alpha\geq 1$.
Returning to \eqref{equ:energybound}, we can now use
Lemma~\ref{lem:OpBoundRSM} to see that
$$\frac{\| \mathcal{E}^h_{\S} \|_{\rm F}}{\sqrt{m_1 m_2}} ~\leq~ \sqrt{\frac{771}{512}} \left\| {\boldsymbol{\tilde{r}}} - {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_2 + \sqrt{\frac{771}{512s'}} \left\| {\boldsymbol{\tilde{r}}} - {\boldsymbol{\tilde{r}}}_{\Omega^{\rm opt}_{\tilde{h},s'}} \right\|_1 + \gamma$$
as desired.
\end{proof}

The next lemma tells us how many evaluation points we need to randomly generate in Lemmas~\ref{lem:SupportLem2} and~\ref{lem:EntryID_ExactInnerProd} before we can be sure to have the RIP properties required by both Theorem~\ref{thm:pairingFunc} and Lemma~\ref{lem:errorBound} above hold with high probability.\\

\begin{lemma}
Let $\S \subseteq [D]$, $\delta \in (0, 1/2]$, and $\tilde{\delta} = \left(0,\frac{1}{256\alpha^2s'}\right]$.
Furthermore, suppose that $m_1,m_2,s',N,D \in \mathbbm{Z}^+ \setminus \{1\}$, $d \in \mathbbm{Z} \cap [1, D]$, and $p \in (0,1)$ satisfy
$$m_1 \geq a_1 \alpha^4 K_{\S}^2  (s')^2 \cdot \max \left\{ d \ln \left( \frac{DN}{d} \right) \ln(m_1), \ln \left(p^{-1} \right) \right\},$$
and
$$m_2 \geq a_2 K_{\S^c}^2 \delta^{-2} s' \cdot \max \left\{ d \ln^2(s') \ln \left( \frac{DN}{d} \right) \ln(m_2), \ln \left(p^{-1} \right) \right\},$$
where $a_1,a_2 \in \mathbbm{R}^+$ are universal constants.
Then, the samples $\{ {\boldsymbol{w}}_\ell \}_{\ell \in [m_1]} \subset \D_{\S}$ and $\{ {\boldsymbol{z}}_k \}_{k \in [m_2]} \subset \D_{\S^c}$ will both simultaneously satisfy their respective RIP assumptions concerning \eqref{equ:EI_RIPmat2} and \eqref{equ:EI_RIPmat} in Lemmas~\ref{lem:SupportLem2} and~\ref{lem:EntryID_ExactInnerProd} above with probability at least $1 - p$.
\label{lem:SamplingBounds}
\end{lemma}

\begin{proof}

The bounds on both $m_1$ and $m_2$ follow from applications of Theorem~\ref{thm:BOS_RIP}.  To bound $m_1$ we note that the normalized nonzero columns of \eqref{equ:EI_RIPmat2} need to have the RIP of order $(2,\tilde{\delta})$, and have an associated BOS constant of $K_{\S}$.  Furthermore, there will never be more than
\begin{equation}
\left|\mathcal{I}_{N,d} \right|  = {D \choose d} N^d \leq \left( \frac{\mathbbm{e}D}{d} \right)^d N^d = \left( \frac{\mathbbm{e}DN}{d} \right)^d
\label{equ:IdNCardBound}
\end{equation}
nonzero columns of \eqref{equ:EI_RIPmat2} for any choice of $\S \subseteq [D]$.  As a consequence we can see that it suffices to have
\begin{align*}
m_1 \geq&~ a'_1 K_{\S}^2 \tilde{\delta}^{-2} \cdot \max \left\{ d \ln \left( \frac{DN}{d} \right) \ln(m_1), \ln \left(p^{-1} \right) \right\}\\
\geq&~ a_1 \alpha^4 K_{\S}^2  (s')^2 \cdot \max \left\{ d \ln \left( \frac{DN}{d} \right) \ln(m_1), \ln \left(p^{-1} \right) \right\}
\end{align*}
in order to satisfy the required RIP conditions for \eqref{equ:EI_RIPmat2} with probability at least $1 - p/2$.

To bound $m_2$ we note that the normalized nonzero columns of \eqref{equ:EI_RIPmat} need to have the RIP of order $(s',\delta)$, and have an associated BOS constant of $K_{\S^c}$.  As a result, \eqref{equ:IdNCardBound} together with Theorem~\ref{thm:BOS_RIP} implies that it suffices to have
$$m_2 \geq a_2 K_{\S^c}^2 \delta^{-2} s' \cdot \max \left\{ d \ln^2(s') \ln \left( \frac{DN}{d} \right) \ln(m_2), \ln \left(p^{-1} \right) \right\}$$
in order to satisfy the required RIP conditions for \eqref{equ:EI_RIPmat} with probability at least $1 - p/2$.
The final desired probability of success now results from the union bound.
\end{proof}

\begin{remark}
\label{rem:SampBound}
To simplify the appearance of our bounds from Lemma~\ref{lem:SamplingBounds} we will make use of the following additional facts and mild assumptions.  First, we will assume hereafter that both $m_1$ and $m_2$ are less than $\left|\mathcal{I}_{N,d} \right|$.  We consider this a reasonable assumption given that the techniques presented herein should only be used in situations where this is the case.  Furthermore, we will use $\delta = 1/2$ above as this is its largest valid parameter setting, and will also consider $\alpha$ to be a universal constant given that it is ultimately set to a fixed value. Finally, we will also replace our probability of failure parameter $p$ by $c / 2D$ for some small constant $c < 0.01$ (for example) in anticipation of wanting to survive a union bound involving $2D-1$ applications of Lemma~\ref{lem:SamplingBounds} for $2D-1$ different sets of random samples.  This will allow us to assert that any at most $2D-1$ different set valued functions $\F^{s'}_{\S}$ will all simultaneously satisfy both Theorem~\ref{thm:pairingFunc} and Lemma~\ref{lem:errorBound} with a ``high probability'' of at least $0.99$.  Utilizing these simplifications we obtain the simplified sufficient sampling conditions
\begin{align*}
m_1 &\geq~ c'_1 K_{\S}^2  (s')^2 \cdot d^2 \ln^2 \left( \frac{DN}{d} \right) \ln \left(D \right),\\
m_2 &\geq~ c'_2 K_{\S^c}^2 s' \cdot d^2 \ln^2(s') \ln^2 \left( \frac{DN}{d} \right) \ln \left(D \right)
\end{align*}
for new absolute constants $c'_1, c'_2 \in \mathbbm{R}^+$.

Finally, and perhaps most controversially, we will make the additional assumption above that either $(i)$ the BOS constants $K_j$ are $1$ for all but at most $\tilde{d} \in \mathbbm{Z} \cap [0, D]$ BOS basis sets $\B_j$ (note that $\tilde{d}$ can be independent of $d$), or else that $(ii)$ $K_0 = 1$.  In either case we will have that both $K_{\S}$ and $K_{\S^c}$ will be bounded above by a constant that depends on only $\tilde{d}$ or $d$, respectively, for all $\S \subseteq [D]$.  In particular, in case $(i)$ we will have that $K_{\S}$ and $K_{\S^c}$ are both at most $K^{\tilde{d}}_\infty$, and in case $(ii)$ that $K_{\S}$ and $K_{\S^c}$ are both at most $K^{d}_\infty$.  Utilizing this final assumption now allows us to bound the total number of samples we need in order to compute any $2D-1$ informative $\F^{s'}_{\S}$ sets for any given $2D-1$ sets $\S$ with high probability (w.h.p.) by either
$$m_1 m_2 \geq c_1 K^{4\tilde{d}}_\infty (s')^3 d^4 \cdot  \ln^4 \left( \frac{DN}{d} \right) \ln^2 (s') \ln^2 (D)$$
in case $(i)$ (note here that letting $d = D$ still avoids exponential dependence on $D$ in this setting), and by
$$m_1 m_2 \geq c_2 K^{4d}_\infty (s')^3 d^4 \cdot  \ln^4 \left( \frac{DN}{d} \right) \ln^2 (s') \ln^2 (D)$$
in case $(ii)$, where $c_1,c_2 \in \mathbbm{R}^+$ are absolute constants.
\end{remark}

We are now ready to demonstrate the numerical performance of our proposed method.