Greedy Recombination Interpolation Method (GRIM)

This article considers finding sparse approximations of functions with the aim of reducing computational complexity. Applications of sparse representations are wide ranging and include the theory of compressed sensing [DE03, CT05, Don06, CRT06], image processing [EMS08, BMPSZ08, BMPSZ09], facial recognition [GMSWY09], data assimilation [MM13], explainability [FMMT14, DF15], sensor placement in nuclear reactors [ABGMM16, ABCGMM18], reinforcement learning [BS18, HS19], DNA denoising [KK20], and inference acceleration within machine learning [ABDHP21, NPS22].

Loosely speaking, the commonly shared aim of sparse approximation problems is to approximate a complex system using only a few elementary features. In this generality, the problem is commonly tackled via Least Absolute Shrinkage and Selection Operator (LASSO) regression, which involves solving a minimisation problem under $l^{1}$-norm constraints. Constraining the $l^{1}$-norm is a convex relaxation of constraining the $l^{0}$-pseudo-norm, which counts the number of non-zero components of a vector. Since $p=1$ is the smallest positive real number for which the $l^{p}$-norm is convex, constraining the $l^{1}$-norm can be viewed as the best convex approximation of constraining the $l^{0}$-pseudo-norm (i.e. constraring the number of non-zero components). Implementations of LASSO techniques can lead to near optimally sparse solutions in certain contexts [DET06, Tro06, Ela10]. The terminology ”LASSO” originates in [Tib96], though imposing $l^{1}$-norm constraints is considered in the earlier works [CM73, SS86]. More recent LASSO-type techniques may be found in, for example, [LY07, DGOY08, GO09, XZ16, TW19].

Alternatives to LASSO include the Matching Pursuit (MP) sparse approximation algorithm introduced in [MZ93] and the CoSaMP algorithm introduced in [NT09]. The MP algorithm greedily grows a collection of non-zero weights one-by-one that are used to construct the approximation. The CoSaMP algorithm operates in a similar greedy manner but with the addition of two key properties. The first is that the non-zero weights are no longer found one-by-one; groups of several non-zero weights may be added at each step. Secondly, CoSaMP incorporates a naive prunning procedure at each step. After adding in the new weights at a step, the collection is then pruned down by retaining only the $m$-largest weights for a chosen integer $m\in{\mathbb{Z}}_{\geq 1}$. The ethos of CoSaMP is particularly notworthy for our purposes.

In this paper we assume that the system of interest is known to be a linear combination of some (large) number of features. Within this setting the goal is to identify a linear combination of a strict sub-collection of the features (i.e. not all the features) that gives, in some sense, a good approximation of the system (i.e. of the original given linear combination of all the features). Depending on the particular context considered, techniques for finding such sparse approximations include the pioneering Empirical Interpolation Method [BMNP04, GMNP07, MNPP09], its subsequent generalisation the Generalised Empirical Interpolation Method (GEIM) [MM13, MMT14], Pruning [Ree93, AK13, CHXZ20, GLSWZ21], Kernel Herding [Wel09a, Wel09b, CSW10, BLL15, BCGMO18, TT21, PTT22], Convex Kernel Quadrature [HLO21], and Kernel Thinning [DM21a, DM21b, DMS21]. The LASSO approach based on $l^{1}$-regularisation can still be used within this framework.

It is convenient for our purposes to consider the following loose categorisation of techniques for finding sparse approximations in this setting; those that are growth-based, and those that are thinning-based. Growth-based methods seek to inductively increase the size of a sub-collection of features, until the sub-collection is rich enough to well-approximate the entire collection of features. Thinning-based methods seek to inductively identify features that may be discarded without significantly affecting how well the remaining features can approximate the original entire collection. Of the techniques mentioned above, MP, EIM, GEIM and Kernel Herding are growth-based, whilst LASSO, Convex Kernel Quadrature and Kernel Thinning are thinning-based. An important observation is that the CoSaMP algorithm combines both growth-based and thinning-based techniques.

In this paper we work in the setting considered by Maday et al. when introducing GEIM [MM13, MMT14, MMPY15]. Namely, we let $X$ be a real Banach space and ${\cal N}\in{\mathbb{Z}}_{\geq 1}$ be a (large) positive integer. Assume that ${\cal F}=\{f_{1},\ldots,f_{{\cal N}}\}\subset X$ is a collection of non-zero elements, and that $a_{1},\ldots,a_{{\cal N}}\in{\mathbb{R}}\setminus\{0\}$. Consider the element $\varphi\in X$ defined by $\varphi:=\sum_{i=1}^{{\cal N}}a_{i}f_{i}$. Let $X^{\ast}$ denote the dual of $X$ and suppose that $\Sigma\subset X^{\ast}$ is a finite subset with cardinality $\Lambda\in{\mathbb{Z}}_{\geq 1}$. We use the terminology that the set ${\cal F}$ consists of the features whilst the set $\Sigma$ consists of data. Then we consider the following sparse approximation problem. Given $\varepsilon>0$, find an element $u=\sum_{i=1}^{{\cal N}}b_{i}f_{i}\in\operatorname{Span}({\cal F})\subset X$ such that the cardinality of the set $\{i\in\{1,\ldots,{\cal N}\}:b_{i}\neq 0\}$ is less than ${\cal N}$ and that $u$ is close to $\varphi$ throughout $\Sigma$ in the sense that, for every $\sigma\in\Sigma$, we have $|\sigma(\varphi-u)|\leq\varepsilon$.

The tasks of approximating sums of continuous functions on finite subsets of Euclidean space, cubature for empirical measures, and kernel quadrature for kernels defined on finite sets are included as particular examples within this general mathematical framework (see Section 2) for full details). The inclusion of approximating a linear combination of continuous functions within this general framework ensures that several machine learning related tasks are included. For example, each layer of a neural network is typically given by a linear combination of continuous functions. Hence the task of approximating a layer within a neural network by a sparse layer is covered within this general framework. Observe that there is no requirement that the original layer is fully connected; in particular, it may itself already be a sparse layer. Consequently, any approach to the approximation problem covered in the general framework above could be used to carry out the reduction step (in which a fixed proportion of a models weights are reduced to zero) in the recently proposed sparse training algorithms [GLMNS18, EJOPR20, LMPY21, JLPST21, FHLLMMPSWY23].

The GEIM approach developed by Maday et al. [MM13, MMT14, MMPY15] involves dynamically growing a subset $F\subset{\cal F}$ of the features and a subset $L\subset\Sigma$ of the date. At each step a new feature from ${\cal F}$ is added to $F$, and a new piece of data from $\Sigma$ is added to $L$. An approximation of $\varphi$ is then constructed via the following interpolation problem. Find a linear combination of the elements in $F$ that coincides with $\varphi$ throughout the subset $L\subset\Sigma$. This dynamic growth is feature-driven in the following sense. The new feature to be added to $F$ is determined, and this new feature is subsequently used to determine how to extend the collection $L$ of linear functionals on which an approximation will be required to match the target.

The growth is greedy in the sense that the element chosen to be added to $F$ is, in some sense, the one worst approximated by the current collection $F$. If we let $f\in{\cal F}$ be the newly selected feature and $J[f]$ be the linear combination of the previously selected features in $F$ that coincides with $f$ on $L$, then the element to be added to $L$ is the one achieving the maximum absolute value when acting on $f-J[f]$. A more detailed overview of GEIM can be found in Section 2 of this paper; full details of GEIM may be found in [MM13, MMT14, MMPY15].

Momentarily restricting to the particular task of kernel quadrature, the recent thinning-based Convex Kernel Quadrature approach proposed by Satoshi Hayakawa, the first author, and Harald Oberhauser in [HLO21] achieves out performs existing techniques such as Monte Carlo, Kernel Herding [Wel09a, Wel09b, CSW10, BLL15, BCGMO18, PTT22], Kernel Thinning [DM21a, DM21b, DMS21] to obtain new state-of-the-art results. Central to this approach is the recombination algorithm [LL12, Tch15, ACO20]. Originating in [LL12] as a technique for reducing the support of a probability measure whilst preserving a specified list of moments, at its core recombination is a method of reducing the number of non-zero components in a solution of a system of linear equations whilst preserving convexity.

An improved implementation of recombination is given in [Tch15] that is significantly more efficient than the original implementation in [LL12]. A novel implementation of the method from [Tch15] was provided in [ACO20]. A modified variant of the implementation from [ACO20] is used in the convex kernel quadrature approach developed in [HLO21]. A method with the same complexity as the original implementation in [LL12] was introduced in the works [FJM19, FJM22].

Returning to our general framework, we develop the Greedy Recombination Interpolation Method (GRIM) which is a novel hybrid combination of the dynamic growth of a greedy selection algorithm, in a similar spirit to GEIM [MM13, MMT14, MMPY15], with the thinning reduction of recombination that underpins the successful convex kernel quadrature approach of [HLO21]. GRIM dynamically grows a collection of linear functionals $L\subset\Sigma$. After each extension of $L$, we apply recombination to find an approximation of $\varphi$ that coincides with $\varphi$ throughout $L$ (cf. the recombination thinning Lemma 3.1). Subsequently, for a chosen integer $m\in{\mathbb{Z}}_{\geq 1}$, we extend $L$ by adding the $m$ linear functionals from $\Sigma$ achieving the $m$ largest absolute values when applied to the difference between $\varphi$ and the current approximation of $\varphi$ (cf. Section 4). We inductively repeat these steps a set number of times. Evidently GRIM is a hybrid of growth-based and thinning-based techniques in a similar spirit to the CoSaMP algorithm [NT09].

Whilst the dynamic growth of GRIM is in a similar spirit to that of GEIM, there are several important distinctions between GRIM and GEIM. The growth in GRIM is data-driven rather than feature-driven. The extension of the data to be interpolated with respect to in GRIM does not involve making any choices of features from ${\cal F}$. The new information to be matched is determined by examining where in $\Sigma$ the current approximation is furthest from the target $\varphi$ (cf. Section 4).

Expanding on this point, we only dynamically grow a subset $L\subset\Sigma$ of the data and do not grow a subset $F\subset{\cal F}$ of the features. In particular, we do not pre-determine the features that an approximation will be a linear combination of. Instead, the features to be used are determined by recombination (cf. the recombination thinning Lemma 3.1). Indeed, besides an upper bound on the number of features being used to construct an approximation (cf. Section 4), we have no control over the features used. Allowing recombination and the data $\Sigma$ to determine which features are used means there is no requirement to use the features from one step at subsequent steps. There is no requirement that any of the features used at one specific step must be used in any of the subsequent steps.

GRIM is not limited to extending the subset $L\subset\Sigma$ by a single linear functional at each step. GRIM is capable of extending $L$ by $m$ linear functionals for any given integer $m\in{\mathbb{Z}}_{\geq 1}$ (modulo restrictions to avoid adding more linear functionals than there are in the original subset $\Sigma\subset X^{\ast}$ itself). The number of new linear functionals to be added at a step is a hyperparameter that can be optimised during numerical implementations.

Unlike [HLO21], our use of recombination is not restricted to the setting of reducing the support of a discrete measure. After each extension of $L\subset\Sigma$, we use recombination [LL12, Tch15, ACO20] to find an element $u\in\operatorname{Span}({\cal F})$ satisfying that $u\equiv\varphi$ throughout $L$. Recombination is applied to the linear system determined by the combination of the set $\{\sigma(\varphi):\sigma\in L\}$, for a given $\sigma\in L$ we get the equation $\sum_{i=1}^{{\cal N}}a_{i}\sigma(f_{i})=\sigma(\varphi)$, and the sum of the coefficients $a_{1},\ldots,a_{{\cal N}}$ (i.e. the trivial equation $a_{1}+\ldots+a_{{\cal N}}=\sum_{i=1}^{{\cal N}}a_{i}$) (cf. the recombination thinning Lemma 3.1).

Our use of recombination means, in particular, that GRIM can be used for cubature and kernel quadrature (cf. Sections 2 and 4). Since recombination preserves convexity [LL12], the benefits of convex weights enjoyed by the convex kernel quadrature approach in [HLO21] are inherited by GRIM (cf. Section 2).

Moreover, at each step we optimise our use of recombination over multiple permutations of the orderings of the equations determining the linear system to which recombination is applied (cf. the Banach Recombination Step in Section 4). The number of permutations to be considered at each step gives a parameter that may be optimised during applications of GRIM.

Whilst we analyse the complexity cost of the Banach GRIM algorithm (cf. Section 5), computational efficiency is not our top priorities. GRIM is designed to be a one-time tool; it is applied a single time to try and find a sparse approximation of the target system $\varphi\in X$. The cost of its implementation is then recouped through the repeated use of the resulting approximation for inference via new inputs. Thus GRIM is ideally suited for use in cases where the model will be repeatedly computed on new inputs for the purpose of inference/prediction. Examples of such models include, in particular, those trained for character recognition, medical diagnosis, and action recognition.

With this in mind, the primary aim of our complexity cost considerations is to verify that implementing GRIM is feasible. We verify this by proving that, at worst, the complexity cost of GRIM is ${\cal O}(\Lambda^{2}{\cal N}+s\Lambda^{4}\log({\cal N}))$ where ${\cal N}$ is the number of features in ${\cal F}$, $\Lambda$ is the number of linear functionals forming the data $\Sigma$, and $s$ is the maximum number of shuffles considered during each application of recombination (cf. Lemma 5.3).

The remainder of the paper is structured as follows. In Section 2 we fix the mathematical framework that will be used throughout the article and motivate its consideration. In particular, we illustrate some specific examples covered by our framework. Additionally, we summarise the GEIM approach of Maday et al. [MM13, MMT14, MMPY15] and highlight the fundamental differences in the philosophies underlying GEIM and GRIM.

An explanation of the recombination algorithm and its use within our setting is provided in subsection 3. In particular, we prove the Recombination Thinning Lemma 3.1 detailing our use of recombination to find $u\in\operatorname{Span}({\cal F})$ coinciding with $\varphi$ on a given finite subset $L\subset X^{\ast}$.

The Banach GRIM algorithm is both presented and discussed in Section 4. We formulate the Banach Extension Step governing how we extend an existing collection of linear functionals $L\subset\Sigma$, the Banach Recombination Step specifying how we use the recombination thinning Lemma 3.1 in the Banach GRIM algorithm, and provide an upper bounds for the number of features used to construct the approximation at each step.

The complexity cost of the Banach GRIM algorithm is considered in Section 5. We prove Lemma 5.3 establishing the complexity cost of any implementation of the Banach GRIM algorithm, and subsequently establish an upper bound complexity cost for the most expensive implementation.

The theoretical performance of the Banach GRIM algorithm is considered in Section 6. Theoretical guarantees in terms of a specific geometric property of $\Sigma$ are established for the Banach GRIM algorithm in which a single new linear functional is chosen at each step (cf. the Banach GRIM Convergence Theorem 6.2). The specific geometric property is related to a particular packing number of $\Sigma$ in $X^{\ast}$ (cf. Subsection 6.2). The packing number of a subset of a Banach space is closely related to the covering number of the subset. Covering and packing numbers, first studied by Kolmogorov [Kol56], arise in a variety of contexts including eigenvalue estimation [Car81, CS90, ET96], Gaussian Processes [LL99, LP04], and machine learning [EPP00, SSW01, Zho02, Ste03, SS07, Kuh11, MRT12, FS21]

In Section 7 we compare the performance of GRIM against other reduction techniques on three tasks. The first is an $L^{2}$-approximation task motivated by an example in [MMPY15]. The second is a kernel quadrature task using machine learning datasets as considered in [HLO21]. In particular, we illustrate that GRIM matches the performance of the tailor-made convex kernel quadrature technique developed in [HLO21]. The third task is approximating the action recognition model from [JLNSY17] for the purpose of inference acceleration. In particular, this task involves approximating a function in a pointwise sense that is outside the Hilbert space framework of the proceeding two examples. Acknowledgements: This work was supported by the DataSıg Program under the EPSRC grant ES/S026347/1, the Alan Turing Institute under the EPSRC grant EP/N510129/1, the Data Centric Engineering Programme (under Lloyd’s Register Foundation grant G0095), the Defence and Security Programme (funded by the UK Government) and the Hong Kong Innovation and Technology Commission (InnoHK Project CIMDA). This work was funded by the Defence and Security Programme (funded by the UK Government).

In this section we rigorously formulate the sparse approximation problem that GRIM will be designed to tackle. Further, we briefly summarise the Generalised Empirical Interpolation Method (GEIM) of Maday et al. [MM13, MMT14, MMPY15] to both highlight some of the ideas we utilise in GRIM, and to additionally highlight the key novel properties not satisfied by GEIM that will be satisfied by GRIM. We first fix the mathematical framework in which we will work for the remainder of this paper.

Let $X$ be a Banach space and ${\cal N}\in{\mathbb{Z}}_{>0}$ be a (large) positive integer. Assume that ${\cal F}=\{f_{1},\ldots,f_{{\cal N}}\}\subset X$ is a collection of non-zero elements, and that $a_{1},\ldots,a_{{\cal N}}\in{\mathbb{R}}\setminus\{0\}$. Consider the element $\varphi\in X$ defined by

Let $X^{\ast}$ denote the dual of $X$ and suppose that $\Sigma\subset X^{\ast}$ is a finite subset of cardinality $\Lambda\in{\mathbb{Z}}_{\geq 1}$. We use the terminology that the set ${\cal F}$ consists of features whilst the set $\Sigma$ consists of data. We consider the following sparse approximation problem. Given $\varepsilon>0$, find an element $u=\sum_{i=1}^{{\cal N}}b_{i}f_{i}\in\operatorname{Span}({\cal F})$ such that the cardinality of the set $\left\{i\in\{1,\ldots,{\cal N}\}~{}:~{}b_{i}\neq 0\right\}$ is less than ${\cal N}$ and that $u$ is close to $\varphi$ throughout $\Sigma$ in the sense that, for every $\sigma\in\Sigma$, we have $|\sigma(\varphi-u)|\leq\varepsilon$.

The task of finding a sparse approximation of a sum of continuous functions, defined on a finite set of Euclidean space, is within this framework. More precisely, let $N,d,E\in{\mathbb{Z}}_{\geq 1}$, $a_{1},\ldots,a_{N}\in{\mathbb{R}}$, $\Omega\subset{\mathbb{R}}^{d}$ be a finite subset, and, for each $i\in\{1,\ldots,N\}$, $f_{i}\in C^{0}\left(\Omega;{\mathbb{R}}^{E}\right)$ be a continuous function $\Omega\to{\mathbb{R}}^{E}$. Then finding a sparse approximation of the continuous function $F:=\sum_{i=1}^{N}a_{i}f_{i}$ is within this framework. To see this, first let $e_{1},\ldots,e_{E}\in{\mathbb{R}}^{E}$ be the standard basis of ${\mathbb{R}}^{E}$ that is orthonormal with respect to the Euclidean dot product $\left<\cdot,\cdot\right>_{{\mathbb{R}}^{E}}$ on ${\mathbb{R}}^{E}$. Then note that for each point $p\in\Omega$ and every $j\in\{1,\ldots,E\}$ the mapping $f\mapsto\left<f(p),e_{j}\right>_{{\mathbb{R}}^{E}}$ determines a linear functional $C^{0}(\Omega;{\mathbb{R}}^{E})\to{\mathbb{R}}$ that is in the dual space $C^{0}(\Omega;{\mathbb{R}}^{E})^{\ast}$. Denote this linear functional by $\delta_{p,j}:C^{0}(\Omega;{\mathbb{R}}^{E})\to{\mathbb{R}}$. Therefore by choosing $X:=C^{0}(\Omega;{\mathbb{R}}^{E})$, ${\cal N}:=N$, and $\Sigma:=\left\{\delta_{p,j}:p\in\Omega~{}\text{and}~{}j\in\{1,\ldots,E\}\right\}\subset X^{\ast}$, we see that this problem is within our framework. Here we are also using the observation that if $f,h\in C^{0}(\Omega;{\mathbb{R}}^{E})$ satisfy, for every $p\in\Omega$ and every $j\in\{1,\ldots,E\}$, that $|\delta_{p,j}[f-h]|\leq\varepsilon$, then we have $||f-h||_{C^{0}(\Omega;{\mathbb{R}}^{E})}\leq C\varepsilon$ for a constant $C\geq 1$ depending on the particular norm chosen on ${\mathbb{R}}^{E}$. Thus finding an approximation $u$ of $F$ that satisfies, for every $\sigma\in\Sigma$, that $|\sigma(F-u)|\leq\varepsilon/C$ allows us to conclude that $||F-u||_{C^{0}(\Omega;{\mathbb{R}}^{E})}\leq\varepsilon$.

The cubature problem [Str71] for empirical measures, which may be combined with sampling to offer an approach to the cubature problem for general signed measures, is within this framework. More precisely, let $N\in{\mathbb{Z}}_{\geq 1}$, $a_{1},\ldots,a_{N}>0$, $\Omega=\{x_{1},\ldots,x_{N}\}\subset{\mathbb{R}}^{d}$, and ${\cal M}[\Omega]$ denote the collection of finite signed measures on $\Omega$. Recall that ${\cal M}[\Omega]$ can be viewed as a subset of $C^{0}(\Omega)^{\ast}$ by defining, for $\nu\in{\cal M}[\Omega]$ and $\psi\in C^{0}(\Omega)$, $\nu[\psi]:=\int_{\Omega}\psi(x)d\nu(x)$. Consider the empirical measure $\mu:=\sum_{i=1}^{N}a_{i}\delta_{x_{i}}$ and, for $e=(e_{1},\ldots,e_{d})\in{\mathbb{Z}}_{\geq 0}^{d}$, define $p_{e}:{\mathbb{R}}^{d}\to{\mathbb{R}}$ by $p_{e}(x_{1},\ldots,x_{d}):=x_{1}^{e_{1}}\ldots x_{d}^{e_{d}}$. Then the choices that $X:={\cal M}[\Omega]$, ${\cal N}:=N$, and, for a given $K\in{\mathbb{Z}}_{\geq 0}$, $\Sigma:=\left\{p_{e}:e=(e_{1},\ldots,e_{d})\text{ with }e_{1}+\ldots+e_{d}\leq K\right\}\subset C^{0}(\Omega)\subset{\cal M}[\Omega]^{\ast}$ illustrate that the cubature problem of reducing the support of $\mu$ whilst preserving its moments of order no greater than $K$ is within our framework.

Moreover the Kernel Quadrature problem for empirical probability measures is within our framework. To be more precise, let ${\cal X}$ be a finite set and ${\cal H}_{k}$ is a Reproducing Kernel Hilbert Space (RKHS) associated to a positive semi-definite symmetric kernel function $k:{\cal X}\times{\cal X}\to{\mathbb{R}}$ (appropriate definitions can be found, for example, in [BT11]). In this case ${\cal H}_{k}=\operatorname{Span}\left(\left\{k_{x}:x\in{\cal X}\right\}\right)$ where, for $x\in{\cal X}$, the function $k_{x}:{\cal X}\to{\mathbb{R}}$ is defined by $k_{x}(z):=k(x,z)$ (see, for example, [BT11]).

In this context one can consider the Kernel Quadrature problem, for which the Kernel Herding [Wel09a, Wel09b, CSW10, BLL15, BCGMO18, TT21, PTT22], Convex Kernel Quadrature [HLO21] and Kernel Thinning [DM21a, DM21b, DMS21] methods have been developed. Given a probability measure $\mu\in{\mathbb{P}}[{\cal X}]$ the Kernel Quadrature problem involves finding, for some $n\in{\mathbb{Z}}_{\geq 1}$, points $z_{1},\ldots,z_{n}\in{\cal X}$ and weights $w_{1},\ldots,w_{n}\in{\mathbb{R}}$ and such that the measure $\mu_{n}:=\sum_{j=1}^{n}w_{j}\delta_{z_{j}}$ approximates $\mu$ in the sense that, for every $f\in{\cal H}_{k}$, we have $\mu_{n}(f)\approx\mu(f)$. Additionally requiring the weights $w_{1},\ldots,w_{n}$ to be convex in the sense that they are all positive and sum to one (i.e. $w_{1},\ldots,w_{n}>0$ and $w_{1}+\ldots+w_{n}=1$) ensures both better robustness properties for the approximation $\mu_{n}\approx\mu$ and better estimates when the $m$-fold product of quadrature formulas is used to approximate $\mu^{\otimes m}$ on ${\cal X}^{\otimes m}$; see [HLO21].

A consequence of ${\cal H}_{k}=\operatorname{Span}\left(\left\{k_{x}:x\in{\cal X}\right\}\right)$ is that linearity ensures that this approximate equality will be valid for all $f\in{\cal H}_{k}$ provided it is true for every $f\in\{k_{x}:x\in{\cal X}\}$. The inclusion ${\cal H}_{k}\subset C^{0}({\cal X})$ means that the subset $\{k_{x}:x\in{\cal X}\}\subset{\cal H}_{k}$ can be viewed as a finite subset of the dual space ${\cal M}[{\cal X}]^{\ast}$. The choice of $X:={\cal M}[{\cal X}]$ and $\Sigma:=\{k_{x}:x\in{\cal X}\}$ illustrates that the Kernel Quadrature problem for empirical probability distributions $\mu\in{\mathbb{P}}[{\cal X}]$ is within the framework we consider.

The framework we consider is the same as the setup for which GEIM was developed by Maday et al. [MM13, MMT14, MMPY15]. It is useful to briefly recall this method. GEIM

Firstly, the growth of the subset $L_{n-1}$ to $L_{n}$ is feature-driven. That is, a new feature from ${\cal F}$ to be used by the next approximation is selected, and then this new feature is used to determine the new functional from $\Sigma$ to be added to the collection on which we require the approximation to coincide with the target $\varphi$ (cf. GEIM (B)). Since we only seek to approximate $\varphi$ over the data $\Sigma$, data-driven growth would be preferable. That is, we would rather select the new information from $\Sigma$ to be matched by the next approximation before any consideration is given to determining which features from ${\cal F}$ will be used to construct the new approximation.

Secondly, related to the first aspect, the features from ${\cal F}$ to be used to construct ${\cal J}_{n}[\varphi]$ are predetermined. Further, the features used to construct ${\cal J}_{n}[\varphi]$ are forced to be included in the features used to construct ${\cal J}_{m}[\varphi]$ for any $m>n$. This restriction is not guaranteed to be sensible; it is conceivable that the features that work well at one stage are disjoint from the features that work well at another. We would prefer that the features used to construct an approximation be determined by the target $\varphi$ and the data selected from $\Sigma$ on which we require the approximation to match $\varphi$. This would avoid retaining features used at one step that become less effective at later steps, and could offer insight regarding which of the features in ${\cal F}$ are sufficient to capture the behaviour of $\varphi$ on $\Sigma$.

Thirdly, at each step GEIM can provide an approximation for any $w\in\operatorname{Span}({\cal F})$. That is, for $n\in{\mathbb{Z}}_{\geq 1}$, the element ${\cal J}_{n}[w]$ is a linear combination of the elements in $S_{n}$ that coincides with $w$ on $L_{n}$. Requiring the operator ${\cal J}_{n}:\operatorname{Span}({\cal F})\to\operatorname{Span}(S_{n})$ to provide such an approximation for every $w\in\operatorname{Span}({\cal F})$ stops the method from exploiting any advantages available for the particular choice of $\varphi\in\operatorname{Span}({\cal F})$. As we only aim to approximate $\varphi$ itself, we would prefer to remove this requirement and allow the method the opportunity to exploit advantages resulting from the specific choice of $\varphi\in\operatorname{Span}({\cal F})$.

All three aspects are addressed in GRIM. The greedy growth of a subset $L\subset\Sigma$, determining the functionals in $\Sigma$ at which an approximation is required to agree with $\varphi$ is data-driven. At each step the desired approximation is found using recombination [LL12, Tch15] so that the features used to construct the approximation are determined by $\varphi$ and the subset $L$ and, in particular, are not predetermined. This use of recombination ensures both that GRIM only produces approximations of $\varphi$, and that GRIM can exploit advantages resulting from the specific choice of $\varphi\in\operatorname{Span}({\cal F})$.

In this section we illustrate how, given a finite collection $L\subset X^{\ast}$ of linear functionals, recombination [LL12, Tch15] can be used to find an element $u\in\operatorname{Span}({\cal F})\subset X$ that coincides with $\varphi$ throughout $L$, provided we can compute the values $\sigma(f_{i})$ for every $i\in\{1,\ldots,{\cal N}\}$ and every linear functional $\sigma\in L$. The recombination algorithm was initially introduced by Christian Litterer and the first author in [LL12]; a substantially improved implementation was provided in the PhD thesis of Maria Tchernychova [Tch15]. A novel implementation of the method from [Tch15] was provided in [ACO20]. A method with the same complexity as the original implementation in [LL12] was introduced in the works [FJM19, FJM22]. Recombination has been applied in a number of contexts including particle filtering [LL16], kernel quadrature [HLO21], and mathematical finance [NS21].

For the readers convenience we briefly overview the ideas involved in the recombination algorithm. For this illustrative purpose consider a linear system of equations ${\bf A}{\bf x}={\bf y}$ where ${\bf A}\in{\mathbb{R}}^{m\times k}$, ${\bf x}\in{\mathbb{R}}^{m}$, ${\bf y}\in{\mathbb{R}}^{k}$, and $k,m\in{\mathbb{Z}}_{\geq 1}$ with $m\geq k$. We assume, for every $i\in\{1,\ldots,k\}$, that ${\bf x}_{i}>0$. Recombination relies on the simple observation that this linear system is preserved under translations of ${\bf x}$ by elements in the kernel of the matrix ${\bf A}$. To be more precise, if ${\bf x}$ satisfies that ${\bf A}{\bf x}={\bf y}$, and if ${\bf e}\in\ker({\bf A})$ and $\theta\in{\mathbb{R}}$, then ${\bf x}+\theta{\bf e}$ also satisfies that ${\bf A}\left({\bf x}+\theta{\bf e}\right)={\bf y}$. The recombination algorithm makes use of linearly independent elements in the $\ker({\bf A})$ to reduce the number of non-zero entries in the solution vector ${\bf x}$.

As outlined in [LL12], this could in principle be done as follows. First, we choose a basis for the kernel of ${\bf A}$. Computing the Singular Value Decomposition (SVD) of ${\bf A}$ gives a method of finding such a basis that is well-suited to dealing with numerical instabilities. Supposing that $\ker({\bf A})\neq\{0\}$, let ${\bf e}_{1},\ldots,{\bf e}_{l}$ be a basis of $\ker({\bf A})$ found via SVD. For each $j\in\{1,\ldots,l\}$ denote the coefficients of ${\bf e}_{j}$ by ${\bf e}_{j,1},\ldots,{\bf e}_{j,m}\in{\mathbb{R}}$.

Consider the element ${\bf e}_{1}\in{\mathbb{R}}^{m}$. Choose $i\in\{1,\ldots,m\}$ such that

Replace ${\bf x}$ by the vector ${\bf x}-\left({\bf x}_{i}/{\bf e}_{1,i}\right){\bf e}_{1}$. The new ${\bf x}$ remains a solution to ${\bf A}{\bf x}={\bf y}$, and now additionally satisfies that ${\bf x}_{i}=0$. Moreover, for every $j\in\{1,\ldots,m\}$ such that $j\neq i$, our choice of $i$ in (3.1) ensures that ${\bf x}_{j}\geq 0$. Finally, for $j\in\{2,\ldots,l\}$, replace ${\bf e}_{j}$ by ${\bf e}_{j}-({\bf e}_{j,i}/{\bf e}_{1,i}){\bf e}_{1}$ to ensure that ${\bf e}_{j,i}=0$. This final alteration allows us to repeat the process for the new ${\bf x}$ using the new ${\bf e}_{2}$ in place of ${\bf e}_{1}$ since the addition of any scalar multiple of ${\bf e}_{2}$ to ${\bf x}$ will not change the fact that ${\bf x}_{i}=0$.

After iteratively repeating this procedure for $j=2,\ldots,l$, we obtain a vector ${\bf x}^{\prime}\in{\mathbb{R}}^{m}$ whose coefficients are all non-negative, with at most $m-l$ of the coefficients being non-zero, and still satisfying that ${\bf A}{\bf x}^{\prime}={\bf y}$. That is, the original solution ${\bf x}$ has been reduced to a new solution ${\bf x}^{\prime}$ with at least $l$ fewer non-zero coefficients than the original solution ${\bf x}$.

Observe that the positivity of the original coefficients is weakly preserved. That is, if we let $x^{\prime}_{1},\ldots,x^{\prime}_{m}\in{\mathbb{R}}$ denote the coefficients of the vector ${\bf x}^{\prime}\in{\mathbb{R}}^{m}$, then for each $i\in\{1,\ldots,m\}$ we have that $x^{\prime}_{i}\geq 0$. One consequence of this preservation is that recombination can be used to reduce the support of an empirical measure whilst preserving a given finite collection of moments [LL12]. Moreover, this property is essential to the ground-breaking state-of-the-art convex kernel quadrature method developed by Satoshi Hayakawa, the first author, and Harald Oberhauser in [HLO21].

The implementation of recombination proposed in [LL12] iteratively makes use of the method outlined above, for $m=2k$, applied to linear systems arising via sub-dividing the original system. The improved tree-based method developed by Maria Tchernychova in [Tch15] provides a significantly more efficient implementation. A geometrically greedy algorithm is proposed in [ACO20] to implement the algorithm of [Tch15]. In the setting of our example above, it follows from [Tch15] that the complexity of the recombination algorithm is $O(mk+k^{3}\log(m/k))$.

Having outlined the recombination method, we turn our attention to establishing the following result regarding the use of recombination in our setting.

In this section we detail the Banach GRIM algorithm. The following Banach Extension Step is used to grow a collection of linear functionals from $\Sigma$ at which we require our next approximation of $\varphi$ to agree with $\varphi$. Banach Extension Step

Recall (cf. Section 3) that if recombination is applied to a linear system corresponding to a matrix $A$, then a Singular Value Decomposition SVD of the matrix $A$ is used to find a basis for $\ker(A)$. Re-ordering the rows of the matrix $A$ (i.e. changing the order in which the equations are considered) can potentially result in a different basis for $\ker(A)$ being selected. Thus shuffling the order of the equations can affect the approximation returned by recombination via the recombination thinning Lemma 3.1. We exploit this by optimising the approximation returned by recombination over a chosen number of shuffles of the equations forming the linear system. This is made precise in the following Banach Recombination Step detailing how, for a given subset $L\subset\Sigma$, we use recombination to find an element $u\in\operatorname{Span}({\cal F})$ that is within $\varepsilon_{0}$ of $\varphi$ throughout $\Sigma$. Banach Recombination Step

GRIM uses data-driven growth. For each $m\in\{2,\ldots,M\}$, GRIM first determines the new linear functionals in $\Sigma$ to be added to $\Sigma_{m-1}$ to form $\Sigma_{m}\subset\Sigma$ before using recombination to find an approximation coinciding with $\varphi$ on $\Sigma_{m}$. That is, we first choose the new information that we want our approximation to match before using recombination to both select the elements from ${\cal F}$ and use them to construct our approximation.

Evidently we have the nesting property that $\Sigma_{t_{1}}\subset\Sigma_{t_{2}}$ for integers $t_{1}\leq t_{2}$, ensuring that at each step we are increasing the amount of information that we require our approximation to match. For each integer $m\in\{1,\ldots,M\}$ let $S_{m}\subset{\cal F}$ denote the sub-collection of elements from ${\cal F}$ used to form the approximation $u_{m}$. Recombination is applied to a system of $1+k_{1}+\ldots+k_{m}$ linear equations when finding $u_{m}$, hence we may conclude that $\#(S_{m})\leq\min\left\{1+k_{1}+\ldots+k_{m},{\cal N}\right\}$ (cf. the recombination thinning Lemma 3.1). Besides this upper bound for $\#(S_{m})$, we have no control on the sets $S_{m}$. We impose only that the linear functionals are greedily chosen; the selection of the elements from ${\cal F}$ to form the approximation $u_{m}$ is left up to recombination and determined by the data. In contrast to GEIM, there is no requirement that elements from ${\cal F}$ used to form $u_{m}$ must also be used for $u_{l}$ for $l>m$.