Sparse multiresolution representations with adaptive kernels

TL;DR

This paper introduces a novel sparse multiresolution approach with adaptive kernels for RKHS-based methods, enabling flexible, efficient, and locally adaptive function estimation that overcomes prior limitations in complexity, smoothness heterogeneity, and kernel selection.

Contribution

It proposes a new integral representation and sparse functional programming framework for adaptive, multiresolution kernel methods with theoretical guarantees and practical efficiency.

Findings

Exact and efficient solutions via duality for non-convex sparse problems

Improved handling of heterogeneous smoothness in functions

Reduced computational complexity in kernel-based estimation

Abstract

Reproducing kernel Hilbert spaces (RKHSs) are key elements of many non-parametric tools successfully used in signal processing, statistics, and machine learning. In this work, we aim to address three issues of the classical RKHS based techniques. First, they require the RKHS to be known a priori, which is unrealistic in many applications. Furthermore, the choice of RKHS affects the shape and smoothness of the solution, thus impacting its performance. Second, RKHSs are ill-equipped to deal with heterogeneous degrees of smoothness, i.e., with functions that are smooth in some parts of their domain but vary rapidly in others. Finally, the computational complexity of evaluating the solution of these methods grows with the number of data points, rendering these techniques infeasible for many applications. Though kernel learning, local kernel adaptation, and sparsity have been used to address these issues, many of these approaches are computationally intensive or forgo optimality guarantees. We tackle these problems by leveraging a novel integral representation of functions in RKHSs that allows for arbitrary centers and different kernels at each center. To address the complexity issues, we then write the function estimation problem as a sparse functional program that explicitly minimizes the support of the representation leading to low complexity solutions. Despite their non-convexity and infinite dimensionality, we show these problems can be solved exactly and efficiently by leveraging duality, and we illustrate this new approach in simulated and real data.

Tables

Table 1. Table I : Classification results for ( PI ′ superscript PI ′ \text{PI}^{\prime} ) using the training samples as kernel centers and using centers selected from k-means and ( PII ′′ superscript PII ′′ \text{PII}^{\prime\prime} ) using the centers selected from k-means

Method	Number of Kernels per Classifier	Accuracy
(${\text{PI}}^{\prime }$)	12000	98.83 $\%$
(${\text{PI}}^{\prime }$) with k-means	800	98.16 $\%$
(${\text{PII}}^{\prime \prime }$)	788	98.12 $\%$
(${\text{PII}}^{\prime \prime }$)	731	96.71 $\%$
(${\text{PII}}^{\prime \prime }$)	53	85.66 $\%$

Equations

$$f(\mathbf{x})=\underset{n\to \infty }{\lim }\sum _{j=1}^n{a}_{j}k(\mathbf{x},{\mathbf{z}}_j;w_0)\text{,}$$

$$\underset{f\in {\mathcal{H}}_0}{\mathrm{minimize}}$$

$$subject\; to$$

$$f^{\star }(\cdot )=\sum _{i=1}^N{a}_i^{\star }k(\cdot ,{\mathbf{x}}_i;w_0)\text{.}$$

$$\underset{\{a_i\}\in \mathbb{R}}{\mathrm{minimize}}$$

$$subject\; to$$

$${\hat{y}}_i=f({\mathbf{x}}_i)=\sum _{j=1}^N{a}_{j}k({\mathbf{x}}_i,{\mathbf{x}}_j;w_0)\text{.}$$

$$y_i=\exp \left[-\frac{(x_i-2.5{)}^2}{2}\right]\text{,}i=1,\dots ,N\text{,}$$

$$f(\mathbf{x})=\underset{n\to \infty }{\lim }\sum _{j=1}^n{a}_{j}k(\mathbf{x},{\mathbf{z}}_j;w_j)\text{.}$$

$$h(\cdot )={\int }_{\mathcal{X}\times \mathcal{W}}\alpha (\mathbf{z},w)k(\cdot ,\mathbf{z};w)d\mathbf{z}dw\text{,}$$

$${\alpha }_m(\mathbf{z},w)=\sum _{j=1}^n{a}_j{r}_m(w-w_j)\prod _{k=1}^p{r}_m([\mathbf{z}{]}_k-[{\mathbf{z}}_j{]}_k)\text{,}$$

$$\underset{\alpha \in L_2(\mathcal{X}\times \mathcal{W})}{\mathrm{minimize}}$$

$$subject\; to$$

$${\hat{y}}_i=f({\mathbf{x}}_i)={\int }_{\mathcal{X}\times \mathcal{W}}\alpha (\mathbf{z},w)k({\mathbf{x}}_i,\mathbf{z};w)d\mathbf{z}dw\text{,}$$

$$\Vert \alpha {\Vert }_{L_0}={\int }_{\mathcal{X}\times \mathcal{W}}I\left[\alpha (\mathbf{z},w)\ne 0\right]d\mathbf{z}dw\text{,}$$

$$\underset{\alpha \in L_2(\mathcal{X})}{\mathrm{minimize}}$$

$$subject\; to$$

$${\hat{y}}_i=f({\mathbf{x}}_i)={\int }_{\mathcal{X}}\alpha (\mathbf{z})k({\mathbf{x}}_i,\mathbf{z};w_0)d\mathbf{z}\text{.}$$

$$\underset{\alpha \in L_2(\mathcal{W})}{\mathrm{minimize}}$$

$$subject\; to$$

$${\hat{y}}_i=f({\mathbf{x}}_i)=\sum _{j=1}^M{\int }_{\mathcal{W}}{\alpha }_j(w)k({\mathbf{x}}_i,{\mathbf{z}}_j;w)dw\text{,}$$

$$\mathcal{L}(\alpha ,\hat{\mathbf{y}},\mathbf{\lambda },\mathbf{\mu })$$

$$-\sum _{i=1}^N{\mathbf{\lambda }}_i\int \alpha (\mathbf{z},w)k({\mathbf{x}}_i,\mathbf{z};w)d\mathbf{z}dw$$

$$+\sum _{i=1}^N{\mathbf{\lambda }}_i{\hat{y}}_i+\sum _{i=1}^N{\mathbf{\mu }}_{i}c({\hat{y}}_i,y_i)\text{.}$$

$$g(\mathbf{\lambda },\mathbf{\mu })=\underset{\begin{array}{c}\alpha \in L_2\\ {\hat{y}}_i\in \mathbb{R}\end{array}}{\min }\mathcal{L}(\alpha ,\hat{\mathbf{y}},\mathbf{\lambda },\mathbf{\mu })\text{.}$$

$$\underset{\mathbf{\lambda }\in {\mathbb{R}}^N,\mathbf{\mu }\in {\mathbb{R}}_{+}^N}{\mathrm{maximize}}$$

$$g(\mathbf{\lambda },\mathbf{\mu })$$

$${\mathcal{L}}_{\hat{y}}(\hat{\mathbf{y}},\mathbf{\lambda },\mathbf{\mu })=\sum _{i=1}^N{\mathbf{\mu }}_{i}c({\hat{y}}_i,y_i)+\sum _{i=1}^N{\mathbf{\lambda }}_i{\hat{y}}_i$$

$$\begin{gathered} {\mathcal{L}}_{\alpha }(\alpha ,\mathbf{\lambda })=\int [\frac{1}{2}{\alpha }^2(\mathbf{z},w)+\gamma I\left[\alpha (\mathbf{z},w)\ne 0\right] \\ -\sum _{i=1}^N{\mathbf{\lambda }}_i\alpha (\mathbf{z},w)k({\mathbf{x}}_i,\mathbf{z};w)]d\mathbf{z}dw\text{,} \end{gathered}$$

$${\alpha }_d(\mathbf{z},w;\mathbf{\lambda })=\{\begin{array}{ll}{\bar{\alpha }}_d(\mathbf{z},w;\mathbf{\lambda })\text{,}& \left|{\bar{\alpha }}_d(\mathbf{z},w;\mathbf{\lambda })\right|>\sqrt{2\gamma }\\ 0\text{,}& \text{otherwise}\end{array}$$

$$\underset{\alpha \in L_2}{\inf }\int F(\alpha (x),x)dx=\int \underset{\bar{\alpha }\in \mathbb{R}}{\inf }F(\bar{\alpha },x)dx\text{.}$$

$$F(\bar{\alpha },\mathbf{z},w)=\frac{{\bar{\alpha }}^2}{2}+\gamma I\left[\bar{\alpha }\ne 0\right]-\sum _{i=1}^N{\mathbf{\lambda }}_{i}k({\mathbf{x}}_i,\mathbf{z};w)\bar{\alpha }\text{,}$$

$$F^{\prime }(\bar{\alpha },\mathbf{z},w)=\frac{{\bar{\alpha }}^2}{2}-\sum _{i=1}^N{\mathbf{\lambda }}_{i}k({\mathbf{x}}_i,\mathbf{z};w)\bar{\alpha }+\gamma \text{,}$$

$${\bar{\alpha }}^{\star }(\mathbf{z},w)=\underset{\bar{\alpha }\in \mathbb{R}}{\mathrm{argmin}}{F}^{\prime }(\bar{\alpha },\mathbf{z},w)=\sum _{i=1}^N{\mathbf{\lambda }}_{i}k({\mathbf{x}}_i,\mathbf{z};w)\text{,}$$

$${\alpha }^{\star }(\cdot ,w)={\bar{\alpha }}_w^{\star }(\cdot )=\sum _{i=1}^N{a}_{i}k({\mathbf{x}}_i,\cdot ;w)\text{.}$$

$$h^{\star }(\cdot )={\int }_{\overline{\mathcal{W}}}\left[{\int }_{\overline{\mathcal{X}}}{\bar{\alpha }}_w(\mathbf{z})k(\cdot ,\mathbf{z};w)d\mathbf{z}\right]dw\text{,}$$

Full text

I Introduction

Reproducing kernel Hilbert spaces (RKHSs) are at the core of non-parametric techniques in signal processing, statistics, and machine learning [3, 4, 5, 6, 7, 8, 9]. Their success stems from the fact that they combine the versatility of functional spaces with the tractability of parametric methods. Indeed, despite the richness of functions found in RKHSs, they can be represented as a (possibly infinite) linear combination of simple basis functions called reproducing kernels (or Mercer kernels) [7, 4]. For smooth functions, however, a celebrated variational result known as the representer theorem states that the number of basis is finite and that they are given by the kernels evaluated at the data points [10, 11]. In other words, learning smooth functions in RKHSs is effectively a parametric, finite dimensional problem. Despite their success, RKHS learning methods have two major practical drawbacks: (i) the RKHS must be known a priori and (ii) evaluating the learned function can be computationally prohibitive.

Indeed, the appropriate functional space of the solution is seldom known in practice. Moreover, since the RKHS dictates the shape and smoothness properties of its functions, its choice is application-specific and ultimately affects the learning performance [12, 5, 13, 14, 15]. Although there exist classes of kernels (and thus RKHSs) that can approximate continuous function arbitrarily well [16], they may need a large number of data points to do so. In fact, it is well-known that RKHS methods are not sample efficient for learning functions with varying degrees of smoothness, i.e., functions that are smooth in some parts of their domain but vary rapidly in others [17].

Kernel learning approaches have been proposed to address this issue by fitting a combination of kernels from a predefined set [12, 16, 18] or by using spectral representations of positive-definite functions [19, 20, 21]. Even when the general form of the kernel is known, the choice of smoothness parameter remains (e.g., selecting the bandwidth of Gaussian kernels). In this case, common approaches include grid search with cross-validation [13, 22] or application-specific heuristic such as maximizing the margin of the support-vectors machine (SVM) [14, 15]. These methods, however, quickly become impractical as they often must search over fine grids, use a large set of kernels, or require additional data. They also fail to address the issue of estimating functions with heterogeneous degrees of smoothness. Although this can be done by choosing different RKHSs for different regions of the domain by means of plug-in rules [23], binary optimization [24], hypothesis testing [25], or gradient descent and alternating minimization [26, 27], these solutions come with no optimality guarantee due to the non-convexity of these locally adapted smoothness formulations.

The issue of computational complexity stems directly from the use of the representer theorem [10, 11, 28]. Despite the reduction from the infinite dimensional problem of estimating smooth functions in RKHSs to that of estimating a finite number of parameters, solving of the problem can become computationally expensive because the number of parameters is proportional to the number of data points. Thus, evaluating the function at any point requires as many kernel evaluations as training samples, which in many applications is prohibitively high. This issue is often addressed by imposing a sparsity penalty on the coefficients to reduce the number of kernel evaluations. Greedy heuristics [29, 30] and $\ell_{1}$-norm relaxations [31, 32, 33, 34] are then often used to cope with the combinatorial nature of this problem, which is known to be NP-hard in general [35, 36]. These methods, however, often implicitly rely on the classical representer theorems [10, 11], despite the fact that they no longer hold in the presence of sparsity penalties (see Remark 1). Hence, even if the sparse optimization program could be solved exactly, the solution would remain suboptimal with respect to the original function estimation problem.

In this work, we propose to tackle these issues by simultaneously (i) adapting the kernels (RKHSs) locally and (ii) allowing arbitrary centers instead of constraining the kernels to be evaluated at the training samples. To do so, we put forward an integral representation of functions that lie in the sum space of an uncountable number of RKHSs taken from a parametrized family. This representation does not assume that the kernels are evaluated at any particular points of the domain. Enforcing sparsity on the coefficients of this representation allows us to determine the optimal parameter and center of each kernel locally. Moreover, sparsity has the added benefit of minimizing the number of kernels used. Despite the non-convex (sparsity) and infinite dimensional (functional) nature of the resulting optimization problem, we can leverage the results from [37] to solve it exactly and efficiently using duality. These results also give rise to a new sparse representer theorem.

The paper is structured as follows: Section II gives an overview of reproducing kernel Hilbert spaces and defines the problems associated with learning in these spaces. In section III we formulate the problem using an integral representation of a function in an RKHS. In Section IV we show the solution to our problem in the dual domain, prove that the duality gap between the primal and the dual problem is zero, and present our algorithm. Numerical examples are then used to illustrate the effectiveness of this method at locally identifying the correct kernel parameters as well as the kernel centers in different applications (Section V).

II Learning in RKHSs: The Classical Problem

Given a training set of data pairs $({\mathbf{x}}_{i},y_{i})$, $i=1,\dots,N$, where ${\mathbf{x}}_{i}\in{\mathcal{X}}$ are the observations or independent variables, with ${\mathcal{X}}\subset{\mathbb{R}}^{p}$ compact, and $y_{i}\in{\mathbb{R}}$ is the label or dependent variable, we seek a function $f:{\mathcal{X}}\to{\mathbb{R}}$ in the RKHS ${\mathcal{H}}_{0}$ that fits these data, i.e., such that $c(f({\mathbf{x}}_{i}),y_{i})$ is small for some convex figure of merit $c$, e.g., quadratic loss, hinge loss, or logistic log-likelihood. Formally, an RKHS is a complete, linear function space endowed with a unique reproducing kernel. A reproducing kernel $k:{\mathbb{R}}^{p}\times{\mathbb{R}}^{p}\to{\mathbb{R}}$ is a positive-semidefinite function with the property $\langle f(\cdot),k(\cdot,{\mathbf{x}})\rangle_{{\mathcal{H}}_{0}}=f({\mathbf{x}})$ for any function $f\in{\mathcal{H}}_{0}$ and point ${\mathbf{x}}\in{\mathbb{R}}^{p}$, where $\langle\cdot,\cdot\rangle_{{\mathcal{H}}_{0}}$ denotes the inner product of the Hilbert space ${\mathcal{H}}_{0}$ [7]. Reproducing kernels are often parametrized by a constant that characterizes the smoothness/richness of the RKHS, such as the bandwidth of sincs, the order of polynomial kernels, or the scale/variance of Gaussian (radial basis function, RBF) kernels [3, 4]. Just as each RKHS has a unique kernel, each kernel defines a unique RKHS. Explicitly, every function $f\in{\mathcal{H}}_{0}$ can be written as the pointwise limit of a linear combination of kernels, i.e.,

[TABLE]

where $w_{0}$ denotes the kernel parameter and the ${\mathbf{z}}_{j}\in{\mathcal{X}}$ are called the kernel centers. In other words, the RKHS ${\mathcal{H}}_{0}$ is the completion of the space $\operatorname{span}\{k(\cdot,{\mathbf{z}}\,;\,w_{0})\mid{\mathbf{z}}\in{\mathcal{X}}\}$ [38, Sec. 2]. Note that different parameters $w_{0}$ result in different RKHSs.

There are infinitely many representations of the form (1) that can interpolate a finite set of points. To avoid overfitting the data and obtain a unique solution, a smoothness prior is often used by minimizing the RKHS norm of the solution [3, 4] as in

[TABLE]

Despite its infinite dimensional nature, this optimization problem admits a solution written as a finite linear combination of kernels centered at the sample points, i.e., there exists a solution $f^{\star}$ of (PI) of the form [11, 10]

[TABLE]

This celebrated variational result is known as the representer theorem and lies at the core of the success of RKHS methods by reducing the functional (PI) to the finite dimensional

[TABLE]

Their success notwithstanding, RKHS-based methods have limitations that can hinder their use in practice. Firstly, the kernel function $k$ (i.e., the RKHS ${\mathcal{H}}_{0}$) must be chosen a priori and given that it determines the shape of the solution [see (2)], its choice directly affects the method performance. Even when the form of $k$ is known, selecting the smoothness parameter $w_{0}$ in (1) can be challenging without application-specific prior knowledge. Secondly, functions RKHSs have limited capability to fit functions with heterogeneous degrees of smoothness, i.e., functions that vary slowly in some regions of their domains and rapidly in others [17]. Although adapting $k$ or the parameter $w_{0}$ in (2) has been proposed to address this issue, ($\text{PI}^{\prime}$) then becomes a non-convex program, foregoing efficient solutions and/or optimality guarantees. Finally, although the representer theorem allows us to use the finite dimensional ($\text{PI}^{\prime}$) to solve the functional (PI), note that evaluating the solution (2) requires as many kernel evaluations as observations. Hence, the complexity of the representation grows with the sample size, which is infeasible for large data sets. Moreover, even if finding a sparse set of coefficients $a_{j}$ in ($\text{PI}^{\prime}$) were tractable, which it is not [35, 36], the representation in (2) does not hold in the presence of a sparsity regularization. In fact, there often exist more parsimonious representations that fit the samples as well as any solution of ($\text{PI}^{\prime}$) (see Remark 1).

To overcome these issues, the next section puts forward an integral parametrization of functions in RKHSs. This representation can be used to pose problems that can locally adapt not only the kernel centers, but also the kernel itself (i.e., the RKHS). Moreover, it allows for regularizations beyond smoothness, most notably sparsity. Although the resulting optimization programs are non-convex and infinite dimensional, we show they have zero duality gap and can therefore be solved efficiently and exactly using conventional methods, such as (stochastic) (sub)gradient descent. This also allows us to derive a new sparse representer theorem.

Remark 1**.**

When seeking parsimonious representations of functions in RKHSs, i.e., representations as in (1) for which the $a_{j}$ are sparse, the classical representer theorems leading to (2) do not apply. To see this is the case, take the counter-example illustrated in Figure 1. Let the sample points be taken from an underlying function composed of a single Gaussian kernel, namely

[TABLE]

where $x_{i}\neq 2.5$ for all $i$. What is more, assume that we know the correct RKHS ${\mathcal{H}}_{0}$, i.e., that the kernel function in (1) and (2) is $k(x,z)=\exp\left[-(x-z)^{2}/2\right]$. Then, it is clear that the most parsimonious function in ${\mathcal{H}}_{0}$ that fits the data is $f^{\prime}(\cdot)=\exp\left[-\frac{(\cdot-2.5)^{2}}{2}\right]$. However, $f^{\prime}$ is not in the feasible set of ($\text{PI}^{\prime}$). Hence, though it can find functions with the same approximation error, they will not be the simplest representation, as illustrated in Figure 1.

III Learning an Integral Representation Instead

In the previous section, we have argued that classical RKHS-based methods suffer from three main drawbacks: (i) the RKHS must be fixed a priori, (ii) they are suboptimal for functions with heterogeneous degrees of smoothness, and (iii) the computational complexity of evaluating solutions is proportional to the sample size. Our first step towards addressing these issues is to extend (1) by allowing different RKHSs at different centers. Throughout this work, we assume that the reproducing kernel family $k$ is fixed (e.g., Gaussian kernels) and adapt its smoothness parameter (e.g., bandwidth). Nevertheless, all results hold for any parametrized dictionary of kernels, i.e., they need not belong to the same family. Explicitly, we write

[TABLE]

The function $f$ in (4) now lies in the sum space of a countable number of RKHSs ${\mathcal{H}}=\bigoplus_{j=1}^{\infty}{\mathcal{H}}_{j}$, where the ${\mathcal{H}}_{j}$ are defined by the parameters $w_{j}$. Note that taking $w_{j}=w_{0}$ for all $j$ in (4) recovers the representation (1) of a function in ${\mathcal{H}}_{0}$.

Similar formulations have been proposed to deal with the aforementioned issues, although optimally selecting $w_{j}$ [issue (i)–(ii)] and ${\mathbf{z}}_{j}$ [issue (iii)] remains an open problem [29, 30, 31, 32, 33, 34]. Our second step is therefore to tackle this hurdle by introducing an integral counterpart of (4). Explicitly, we define the function

[TABLE]

where ${\mathcal{W}}$ is a compact subset of ${\mathbb{R}}$ and $\alpha:{\mathcal{X}}\times{\mathcal{W}}\rightarrow{\mathbb{R}}$ is in $L_{2}({\mathcal{X}}\times{\mathcal{W}})$. Notice that, in contrast to (1) and (4), the expression in (5) does not depend on a choice of centers or kernel parameters, thus addressing issues (i) and (ii). Before proceeding, we show that (4) and (5) are essentially equivalent, i.e., that (5) can essentially represent the functions in ${\mathcal{H}}$.

Proposition 1**.**

Let $k$ be a continuous reproducing kernel, i.e., $k(\cdot,{\mathbf{z}};\cdot)$ is continuous over the compact set ${\mathcal{X}}\times{\mathcal{W}}$ for each ${\mathbf{z}}\in{\mathcal{X}}$. Then, for each $f\in{\mathcal{H}}_{p}$ there exists a sequence $\{h_{m}\}$ of functions as in (5) such that $h_{m}\to f$ pointwise.

Proof.

Consider the approximation of the identity $r_{m}(x)=m\operatorname{\mathbb{I}}\left[|x|<1/m\right]$ and note that $r_{m}(x)\to\delta(x)$ weakly in the vague topology, i.e., $\int_{\mathcal{D}}r_{m}(x)\varphi(x)\to\varphi(0)$ for all $\varphi$ continuous and ${\mathcal{D}}$ compact [39]. Now let $f\in{\mathcal{H}}_{p}$ be written as $f(\cdot)=\sum_{j=1}^{n}a_{j}k(\cdot,{\mathbf{z}}_{j}\,;\,w_{j})$ and take $h_{m}(\cdot)=\int_{{\mathcal{X}}\times{\mathcal{W}}}\alpha_{m}({\mathbf{z}},w)k(\cdot,{\mathbf{z}}\,;\,w)d{\mathbf{z}}dw$ with

[TABLE]

where $[{\mathbf{z}}]_{k}$ indicates the $k$-th element of the vector ${\mathbf{z}}$. Note that $\alpha_{m}\in L_{2}$, so that $h_{m}$ is indeed of the form (5). Since the reproducing kernel is continuous, it readily holds from (6) that $g_{m}\to f$ pointwise for all ${\mathbf{x}}\in{\mathcal{X}}$. ∎

Proposition 1 shows that there is no loss in using (5) since it can essentially represent all functions in the sum space ${\mathcal{H}}$ of interest. Though this result is straightforward when (5) defines the inner product in ${\mathcal{H}}$ or if $\alpha$ allowed distributions (Dirac deltas), the former is typically not the case (except for sinc kernels) and the latter does not hold since $\alpha\in L_{2}$. Also, it would appear that this integral representation did nothing but aggravate the computational complexity problems [(iii)], now that the coefficients $\alpha$ are functions. What is more, estimating $\alpha$ is once again a functional problem. In the sequel, we address the computational complexity issue by explicitly minimizing the support of $\alpha$ using sparse functional optimization programs (SFPs). In Section IV, we then show that these optimization problems can be solved efficiently and exactly using duality, providing an explicit algorithm to compute $\alpha$.

III-A A Sparse Functional Formulation

Proposition 1 suggests that the centers and kernel parameters can be obtained by leveraging sparsity. Indeed, notice from (6) that functions in ${\mathcal{H}}$ admit an integral representation (5) in which $\alpha$ is a superposition of bump functions centered around $({\mathbf{z}}_{j},w_{j})$, i.e., a function that vanishes over most of its domain. In other words, functions in ${\mathcal{H}}$ admit a sparse integral parametrization. This observation motivates estimating $f$, equivalently $\alpha$, using the following SFP:

[TABLE]

where $\gamma\geq 0$ is a regularization factor that trades-off smoothness and sparsity; $\|\cdot\|_{L_{2}}$ denotes the $L_{2}$-norm, which induces smoothness, enhances robustness to noise, and improves the numerical properties of the optimization problem; and $\|\cdot\|_{L_{0}}$ denotes the “$L_{0}$-norm,” defined as the measure of the support of a function, i.e.,

[TABLE]

where $\operatorname{\mathbb{I}}[x\neq 0]=1$ if $x\neq 0$ and zero otherwise. Notice that (7) is the functional counterpart of the discrete “$\ell_{0}$-norm”.

Problem (PII) seeks the function $f\in{\mathcal{H}}$ with the sparsest integral representation that fits the data according to the convex loss $c$, which is both non-convex and infinite dimensional. As it is, it therefore appears to be intractable. Before addressing this matter, however, we argue that its solutions would indeed cope with the statistical and computational issues of classical RKHS methods. In Section IV, we then derive efficient algorithms to obtain these solutions.

To be sure, (PII) precludes the choice of a specific RKHS [issue (i)] or kernel centers by leveraging the integral representation (5). Simultaneously, it enables the solution to locally adapt the RKHS over the domain to account for functions with heterogeneous degrees of smoothness [issue (ii)]. Finally, the sparse solutions of (PII) can be used to obtain low complexity representations of functions in the sum space ${\mathcal{H}}$, i.e., representations with a small number of kernels [issue (iii)]. Intuitively, we can obtain a finite series from the integral representation (5) by using the bumps in $\alpha$ to determine the pair $({\mathbf{z}}_{j},w_{j})$ that define (4), as suggested by (6). The $a_{j}$ can then be obtained by directly minimizing $c$.

Although the remainder of this work studies the general (PII), two particular cases are of marked interest. First, the functional space of the solution is sometimes dictated by the application. For instance, one may seek bandlimited functions of a specific bandwidth. In this case, the reproducing kernel $k$ and its parameter $w_{0}$ are fixed and (PII) reduces to

[TABLE]

Problem ($\text{PII}^{\prime}$) sets to find solutions that use as few kernels as possible and its often tackled using greedy methods such as KOMP [30]. Despite its success, its applications either implicitly relies on the representer theorem [9], which does not hold for sparse problems (see Remark 1), or use a grid search over the space, which can quickly become prohibitive. In contrast, the solution of ($\text{PII}^{\prime}$) is guaranteed to provide the sparsest integral representation. These can then be approximated using the aforementioned peak finding method to yield low complexity, discrete solutions. Problem ($\text{PII}^{\prime}$) can therefore obtain solutions with the same cost and lower complexity, as illustrated in Section V.

Second, a set of candidate kernel centers might be available a priori, given by domain experts or unsupervised learning techniques such as clustering. Problem ($\text{PII}^{\prime\prime}$) can then be used to optimally select a subset of these centers as well as determine the appropriate RKHS. Explicitly,

[TABLE]

where ${\mathbf{z}}_{j}$, for $j=1,\dots,M$, are the predefined candidate centers. Observe that ($\text{PII}^{\prime\prime}$) promotes the sparsity of each $\alpha_{j}$, so that the coefficient of those centers that do not contribute to satisfy the fit constraint will vanish. Hence, the solution of ($\text{PII}^{\prime\prime}$) effectively selects the smallest subset of candidate centers. Furthermore, by locally adapting the RKHS, it can further reduce the number of kernels (centers) in the final solution by using less kernels to fit smoother portions of the data. The formulation ($\text{PII}^{\prime\prime}$) is particularly attractive for high dimensional problems for which evaluating the integrals in (PII) may be challenging.

In the next section, we derive a method for solving (PII) exactly and efficiently. Solutions to problems ($\text{PII}^{\prime}$) and ($\text{PII}^{\prime\prime}$) were presented in [2] and [1] respectively. We do so by formulating the dual problem of (PII) and showing that its solution can be used to obtain an optimal $\alpha^{\star}$ (strong duality). This then allows us to propose efficient solutions for (PII) by means of its dual problem. This strong duality result is also exploited to derive a new integral representer theorem (Corollary 2) that accounts for sparsity.

IV Learning in the Dual Domain

Having argued that (PII) [or ($\text{PII}^{\prime}$)–($\text{PII}^{\prime\prime}$)] is worth solving, we now return to the issue of how. To understand the challenge of directly tackling (PII), observe that it is a non-convex, infinite dimensional optimization program. Moreover, its discrete version is in general NP-hard to solve [35]. In the sequel, we address these hurdles using duality. It is worth noting that duality is an established approach to solve semi-infinite convex programs [40, 41, 42]. Indeed, dual problems are convex and their dimension is equal to the number of constraints. Thus, they can be solved efficiently. Moreover, it is well-known that when the original problem is convex, its solutions can be obtained from solutions of the dual problem under mild conditions [43]. Nevertheless, this is not the case of (PII).

In the sequel, we derive a method to solve (PII) by first obtaining its dual problem in closed-form (Section IV-A). Then, we show that strong duality holds (Section IV-B). Consequently, solutions of (PII) can be obtained from solutions of its dual problem. We then conclude by showing how to efficiently solve the latter (Section IV-C).

IV-A The dual problem of (PII)

To derive the dual problem of (PII), start by introducing the Lagrange multipliers $\bm{\lambda}\in{\mathbb{R}}^{N}$, associated with its equality constraints and $\bm{\mu}\in{\mathbb{R}}^{N}_{+}$, associated with its inequality constraints. Its Lagrangian is then defined as

[TABLE]

For conciseness, we omit the set ${\mathcal{X}}\times{\mathcal{W}}$ over which the integrals are computed. To proceed, obtain the dual function by minimizing the Lagrangian as in

[TABLE]

The dual function is the objective of the dual problem, defined explicitly as

[TABLE]

Notice that the dual function (9) is concave regardless of the convexity of the original problem since it is the minimum over a set of affine functions in $(\bm{\lambda},\bm{\mu})$. What is more, despite the infinite dimensionality of the original problem, it is defined over $2N$ variables. Hence, (DII) is a finite dimensional convex program, which can be solved efficiently as long as its objective can be evaluated efficiently.

Yet, computing $g$ involves a functional, non-convex problem. Indeed, notice that (9) can be separated as

[TABLE]

where

[TABLE]

is the value a convex optimization problem, since $c$ is convex and $\bm{\mu}_{i}\geq 0$, and

[TABLE]

where we used the integral definitions of the “$L_{0}$-norm” in (7) and the $L_{2}$-norm. Hence, evaluating $g$ involves solving a non-convex functional optimization problem similar to the original (PII). Here, however, we can exploit separability to obtain the closed-form thresholding solution presented in the following proposition.

Proposition 2**.**

A minimizer $\alpha_{d}$ of (12) is given by

[TABLE]

where $\bar{\alpha}_{d}({\mathbf{z}},w;\bm{\lambda})=\sum_{i=1}^{N}\bm{\lambda}_{i}k({\mathbf{x}}_{i},{\mathbf{z}};w)$.

Proof.

To obtain (13), we start by separating the objective of (12) across ${\mathbf{z}}$ and $w$. To do so, we leverage the following lemma:

Lemma 1**.**

Let $F(\alpha,x)$ be a normal integrand, i.e., continuous in $\alpha$ for fixed $x$ and measurable in $x$ for fixed $\alpha$. Then,

[TABLE]

Proof.

See [44, Thm. 3A]. ∎

Taking

[TABLE]

in Lemma 1, yields that minimizing (12) is equivalent to minimizing $F$ individually for each $({\mathbf{z}},w)$. Although still non-convex, this problem is now scalar and admits a simple solution.

Indeed, notice that the indicator function in (15) can only take two values depending on $\bar{\alpha}$. Hence, its optimal value is the minimum of two cases: (i) if $\bar{\alpha}=0$, then $F(0,{\mathbf{z}},w)=0$ for all $({\mathbf{z}},w)$; alternatively, (ii) if $\bar{\alpha}\neq 0$, then (15) becomes

[TABLE]

whose minimization is a quadratic problem with closed-form solution

[TABLE]

so that $\min_{\bar{\alpha}\neq 0}F(\bar{\alpha}^{\star},{\mathbf{z}},w)=\gamma-\bar{\alpha}^{\star}({\mathbf{z}},w)^{2}/2$. Immediately, $\alpha_{d}({\mathbf{z}},w)=\bar{\alpha}^{\star}({\mathbf{z}},w)$ if $\gamma-\bar{\alpha}_{d}({\mathbf{z}},w)^{2}/2<0$ or $\alpha({\mathbf{z}},w)$ vanishes, which yields (13). ∎

Hence, despite the non-convexity and infinite dimensionality nature of (12), it admits an explicit solution in the form of the thresholded function (13). We can thus evaluate the convex objective of (DII), that can then be solved using classical convex optimization tools such as (stochastic) (sub)gradient ascent (see Section IV-C). Still, the question remains of whether this is a worthwhile endeavor. Indeed, (9) can be interpreted as a relaxation of the hard constraints in (PII), so that (DII) provides a lower bound on its optimal value. Though we can evaluate (13) on the solution $(\bm{\lambda}^{\star},\bm{\mu}^{\star})$ of (DII), there are no guarantees that it is a solution of (PII). We tackle this challenge in the sequel.

IV-B Strong duality and the integral representer theorem

Due to the non-convexity of the original problem, the only immediate guarantee we can give about the optimal value of (DII) is that it is a lower bound on the optimal value of (PII). The central technical result of this section shows that strong duality holds for (PII), i.e., it has null duality gap (Theorem 1). This result has deep implications for the problems we posed in Section II. First, it implies that we can obtain a solution of (PII) by solving its (DII) (Corollary 1). Second, it allows us to write a representer theorem similar to the original ones stating that the solution $\alpha^{\star}$ of (PII) is a linear combination of kernels evaluated at the data points.

Let us start with the main theorem:

Theorem 1**.**

Strong duality holds for (PII) if the kernel $k(\cdot,{\mathbf{z}}\,;\,\omega)$ has no point masses and Slater’s condition is met. In other words, if $P$ is the optimal solution of (PII) and $D$ is the optimal solution of (DII), then the duality gap $P-D$ vanishes.

Proof.

This result can be found in [37]. For ease of reference, a proof is provided in Appendix A. ∎

Theorem (1) states that, though it is not a convex program, (PII) has null duality gap. Immediately, we obtain the following corollary:

Corollary 1**.**

Let $(\bm{\lambda}^{\star},\bm{\mu}^{\star})$ be a solution of (DII) and assume $k\in L_{2}$ and analytic. Then, $\alpha_{d}^{\star}(\cdot,\cdot)=\alpha_{d}(\cdot,\cdot;\bm{\lambda}^{\star})$ is a solution of (PII) for $\alpha_{d}$ as in (13).

Proof.

See Appendix B. ∎

Thus, despite the apparent challenges, (PII) is tractable and can be solved efficiently and exactly using duality (as detailed in Section IV-C). It is worth noting that the technical hypotheses of Corollary 1 are mild and hold for all commonly used reproducing kernels, such as Gaussian, polynomial, and sinc kernels. A fundamental feature of this approach is that it solves the functional problem without relying on discretizations. This is of utmost importance since discretizing (PII) can lead to NP-hard, large dimensional, and potentially ill-conditioned problems. Moreover, it tackles the sparse problem directly instead of using convex relaxations.

Another fundamental implication of Theorem 1 is the following integral representer theorem.

Corollary 2** (Integral representer theorem).**

A solution $\alpha^{\star}$ of (PII) can be obtained by thresholding a parametrized family of functions $\bar{\alpha}^{\star}_{w}\in{\mathcal{H}}_{w}$, where ${\mathcal{H}}_{w}$ is the RKHS induced by the kernel $k(\cdot,\cdot;w)$. In fact, $\bar{\alpha}^{\star}_{w}$ lives in a finite dimensional subspace of ${\mathcal{H}}_{w}$ spanned by the kernels evaluated at the data points. Explicitly, there exist $a_{i}\in{\mathbb{R}}$ such that

[TABLE]

Proof.

From Corollary 1, $\alpha^{\star}=\alpha_{d}^{\star}$ almost everywhere with $\alpha_{d}^{\star}({\mathbf{z}},w)=\alpha_{d}({\mathbf{z}},w;\bm{\lambda}_{i}^{\star})$. Thus, the corollary stems from (13) for $\alpha_{w}(\cdot)=\sum_{i=1}^{N}\bm{\lambda}_{i}^{\star}k({\mathbf{x}}_{i},\cdot;w)$, so that $\alpha_{w}\in{\mathcal{H}}_{w}$ by definition and $\alpha^{\star}({\mathbf{z}},w)=\alpha_{w}({\mathbf{z}})\operatorname{\mathbb{I}}(|\alpha_{w}({\mathbf{z}})|>\sqrt{2\gamma})$. ∎

Recall that the classical representer theorem [10, 11] reduces the functional problem (PI) to the finite dimensional problem ($\text{PI}^{\prime}$) by showing that it has a solution of the form (2) that lies in the span of the kernels evaluated at the data points. Likewise, Theorem 2 states that a solution of (PII) can be obtained from a family of functions with the similar finite representation (18). Note, however, that although classical representer theorems do not account for sparse solution (as argued in Remark 1), Theorem 2 holds in the presence of the “$L_{0}$-norm” regularizer. The cost of doing so is adding a layer of indirection with respect to the original functional problem: whereas the classical representer theorem yields a function in the RKHS directly from the parameter $a_{i}$ through (2), Corollary 2 states that the functional parameters of the integral representation (5) of this function can be obtained from the $a_{i}$. Still, the problem it solves is considerably more complex and comprehensive than (PI).

It is worth noting that although (PII) searches for $\alpha\in L_{2}$, Corollary 2 shows its solution depends only on a family of functions belonging to the RKHSs considered. Explicitly, the solution $\alpha^{\star}(\cdot,w)$ of (PII) is a thresholded version of a function in the RKHS ${\mathcal{H}}_{w}$ with reproducing kernel $k(\cdot,\cdot;w)$. In particular, when considering the problem without sparsity ($\gamma=0$), we have that $\alpha^{\star}(\cdot,w)\in{\mathcal{H}}_{w}$. In the case of ($\text{PII}^{\prime}$), this further simplifies to $\alpha^{\star}\in{\mathcal{H}}_{0}$.

This observation allows us to think of (5) as building the function $h^{\star}$ point-by-point by integrating the value of partial $L_{2}$-inner products between the reproducing kernel of ${\mathcal{H}}_{w}$ and a function $\bar{\alpha}_{w}\in{\mathcal{H}}_{w}$. To be more specific, notice that (5) can be written as the iterated integrals

[TABLE]

where $\overline{{\mathcal{X}}}\subseteq{\mathcal{X}}$, $\overline{{\mathcal{W}}}\subseteq{\mathcal{W}}$, and $\overline{{\mathcal{X}}}\times\overline{{\mathcal{W}}}$ is the set induced by the support of $\alpha^{\star}$, i.e., $\{({\mathbf{z}},w)\in{\mathcal{X}}\times{\mathcal{W}}\mid|\alpha({\mathbf{z}},w)|>\sqrt{2\gamma}\}$. The innermost integral in (19) can be interpreted as an inner product in $L_{2}$ between $\bar{\alpha}_{w}$ and $k(\cdot,{\mathbf{z}};w)$ computed only where the magnitude of $\bar{\alpha}_{w}$ is large enough, defined by the regularization parameter $\gamma$. This sort of trimmed inner product is linked to robust projections found in different statistical methods [45, 46]. The outer integral then accumulates the projections of $\bar{\alpha}_{w}$ over the relevant subset $\overline{{\mathcal{W}}}$ of RKHSs considered to form the functional solution.

Before proceeding, it is worth noting that Theorem 1 holds under very mild conditions. Indeed, the reproducing kernels typically used in applications, such as polynomial or Gaussian kernels, do not have point masses. In fact, if the function of interest is in $L_{2}$, i.e., ${\mathcal{H}}\subseteq L_{2}$, then we need not consider kernels containing Dirac deltas since they are not square integrable. As for Slater’s condition [43], the infinite dimensionality of $\alpha$ makes it so we can always find one that perfectly interpolates the data, though it may neither be smooth nor have a sparse representation. Hence, finding a strictly feasible solution of (PII) is straightforward for most $c$.

In the next section, we leverage the closed-form of the dual function from Proposition 2 and the strong duality result from Theorem 1 to obtain an explicit algorithm to solve (PII).

IV-C Dual gradient ascent

Theorem 1 shows that we can obtain a solution of (PII) through (13). Still, although (13) can be evaluated using the closed form expression from Proposition 2, it requires the optimal dual variables $(\bm{\lambda}^{\star},\bm{\mu}^{\star})$. In this section, we propose a projected supergradient ascent method to solve (DII) (Algorithm 1). Alternatively, other standard convex optimization algorithms can be used to exploit structure in the solution of (DII). For instance, efficient solvers based on coordinate ascent can be leveraged to solve large-scale instances [47].

Start by recalling that a supergradient of a function $f:{\mathcal{D}}\to{\mathbb{R}}$ at ${\mathbf{x}}\in{\mathcal{D}}\subseteq{\mathbb{R}}^{n}$ is any vector ${\mathbf{d}}$ such that $f({\mathbf{y}})\leq f({\mathbf{x}})+{\mathbf{d}}^{T}({\mathbf{y}}-{\mathbf{x}})$ for all ${\mathbf{y}}\in{\mathcal{D}}$. Though supergradients may not be an ascent direction at ${\mathbf{x}}$, taking small steps in its direction decreases the distance to any maximizer of a convex function $f$ [43]. Thus, we can solve (DII) by repeating, for $t=0,1,\dots$,

[TABLE]

where ${\mathbf{d}}_{\bm{\lambda}}^{(t)},{\mathbf{d}}_{\bm{\mu}}^{(t)}$ are the supergradients of $\bm{\lambda}$ and $\bm{\mu}$, respectively, $\eta_{\lambda},\eta_{\mu}>0$ are step sizes, and $[x]_{+}=\max(0,x)$. The projection of $\bm{\mu}_{i}$ on the non-negative numbers guarantees that the constraints of (DII) are satisfied. The supergradient in (20) are readily obtained from the constraint violation of the dual minimizers [43]. Explicitly, let $\hat{y}_{d,i}(\bm{\lambda},\bm{\mu})$ and $\alpha_{d}({\mathbf{z}},w;\bm{\lambda})$ be minimizers of (11) and (12) respectively. Then,

[TABLE]

Since $\hat{y}_{d,i}$ is the solution of the convex optimization problem (11), the update for the dual variables $\bm{\mu}_{i}$ in (20b) can be efficiently evaluated using (21b). The update expression for $\bm{\lambda}_{i}$ in (20a), however, requires that the integral in (21a) be evaluated. To do so, we can either use numerical integration methods, since $\alpha_{d}$ is available in closed-form from (13), or rely on Monte Carlo methods. The latter approach is especially interesting because it can be integrated with the optimization iterations in (20) to obtain a stochastic supergradient ascent algorithm summarized in Algorithm 1. Since Monte Carlo gives an unbiased estimate of $d_{\bm{\lambda}_{i}}$, typical convergence guarantees for stochastic optimization apply [48, 49, 50].

V Numerical Experiments

In the previous sections we have claimed that our algorithm can estimate (i) kernel widths ($\text{PII}^{\prime\prime}$), (ii) kernel centers ($\text{PII}^{\prime}$), and (iii) kernels of varying centers and widths (PII). In this section we show through a sample signal, how we can achieve claim (i). Then we show how moving from (i) to (iii) reduces complexity. In our discussion about the complexity of the representation in section V-A, we show how we can achieve (ii) on random signals of fixed width. In section V-B, we solve (PII) for a signal of varying degrees of smoothness and show how we can reduce complexity regardless of sample size. Lastly, in sections V-C and V-D we apply our algorithm to solve (PII) and ($\text{PII}^{\prime\prime}$) on two examples of real applications: a user localization problem and a digit classification problem.

For the estimation, we search over functions in the family of RKHSs, which have Gaussian functions as kernels

[TABLE]

where width of the kernel is directly proportional to the hyper-parameter $w$.

To start, the effect of the choice of RKHS on the performance of a learning algorithm is examined. To this end, a signal, which lies in the RKHS with a Gaussian kernel of width $w_{0}=0.453$ is constructed. The classical problem in ($\text{PI}^{\prime}$) is compared to the problem presented in ($\text{PII}^{\prime\prime}$). A grid search is used to examine the performance of ($\text{PI}^{\prime}$) for different values of $w$. The value of $w_{0}$ was chosen such that it would not be directly on the grid, since in practice it is unlikely to include the value of the width of the originating signal. We generate $S$ signals of the form

[TABLE]

with $j=1\ldots S$. For each $f_{j}$ a training set of $N=50$ samples was generated with $m=10$. The amplitude $a_{i}$ of each function is selected at random from a uniform distribution $\mathcal{U}(1,2)$. The $\tilde{{\mathbf{x}}}_{i}$ are i.i.d random variables drawn from the uniform distribution $\mathcal{U}(1,2)$ and the $\xi_{j}$ are i.i.d. random variables drawn from ${\mathcal{N}}(0,10^{-3})$, which represent the noise.

It should be noted that, given sufficient iterations, well chosen step sizes, and a large $\gamma$, our method can approximate point masses. However, smoother approximations of the point masses can be obtained by using only few iterations. Additionally, these smooth approximations are more robust to the choice of the tuning parameters. Kernel centers and widths can subsequently be obtained by selecting the extreme points of the function $\alpha({\mathbf{z}},w)$, since the optimal $\alpha({\mathbf{z}},w)$ is a function of $w$, the kernel width, and ${\mathbf{z}}$ the kernel centers. Kernels using the widths and centers approximated from the extreme points of the $\alpha({\mathbf{z}},w)$ are used to train a least squared estimator.

A grid search is performed for problem ($\text{PI}^{\prime}$) by uniformly sampling $w$ over the interval [0,1] at 0.1 increments. The problem ($\text{PII}^{\prime\prime}$) is solved using $\gamma=4000$, $\eta_{\lambda}=0.001$, $\eta_{\mu}=0.1$ and $T=5000$. The performance of the two algorithms is compared over $1000$ repetitions of the sampled signal, each with a training set of size $N=100$ and a test set of size $N_{test}=1000$. The MSE of ($\text{PI}^{\prime}$) decreases as the value of $w$ increases—see Figure 2. Due to the non-uniform sampling of the signal, smoother kernels on average have a better performance. In areas, in which the sampling is sparse, the thinner kernels cannot represent the signal between the samples. Additionally, the thinner kernels are more likely to overfit to the noise than the smoother kernels. However, the smoother kernels cannot model the faster variation in the signal well. In contrast ($\text{PII}^{\prime\prime}$) finds a sparse solution, which uses 14 kernels on average, of varying smoothness, with an average MSE of $0.0457$, which can both take advantage of the ability of smoother kernels to avoid overfitting as well as thinner kernels to model fast variation. Indeed, we observe in Figure 3 that our algorithm chooses a mixture of kernels of width around $0.453$ and kernels of width $1$.

Smoother kernels perform better because of the random sampling combined with the restriction of only using kernels centered at the sample points. Therefore, we investigate the effect of solving problem (PII) which finds both kernel centers and kernel widths. Problem (PII) is solved using $\gamma=1000$, $\eta_{\lambda}=0.01$, $\eta_{\mu}=1$ and $T=1000$ over 1000 randomly sampled training sets, and results in an MSE of $0.0588$. Although the MSE of (PII) is similar to that of ($\text{PII}^{\prime\prime}$), it is important to note that by placing kernels arbitrarily we are able to better estimate the width of the kernel: by comparing Figure 4 to Figure 3 it can be seen that (PII) uses only 1 to 2 kernels per representation of width 1 whereas ($\text{PII}^{\prime\prime}$) uses on average 6 kernels of width 1. Moreover, we consistently obtain representations of lower complexity when solving (PII)—see Figure 5.

V-A Examining the Complexity of the Solution

So far we have shown that the complexity of the formulation can be reduced by moving centers in addition to moving the width. To further explore the effect of kernel centers on the complexity of the solution, we compare the performance of ($\text{PII}^{\prime}$) to that of kernel orthogonal matching pursuit (KOMP) with pre-fitting (see [30, 9]), for a simulated signal as in (23). KOMP takes an initial function and a set of sample points and tries to estimate it by a parsimonious function of a lower complexity. As a backwards feature selection metheod, the algorithm starts by including all samples and then reduces the complexity of the function by reducing one feature at a time. The KOMP algorithm in [30, 9] was modified by changing the stopping criteria to be the estimation error, rather than the distance to the original function. This stopping criteria allows us to compare the sparsity needed to obtain similar estimation error.

The signal was sampled from the function in (23) using $w_{0}=0.5$ and $m=[5,10,20]$ by generating $N=[2m,4m,6m]$ samples for each function, thus creating $9$ different sample size and signal pairs. The problem in ($\text{PII}^{\prime}$) was solved using $\gamma=30$, $\eta_{\lambda}=0.05$, $\eta_{\mu}=0.1$ and $T=1000$. Subsequently, a least squares algorithm was trained using kernels at the location found by our algorithm. Both our method and KOMP used $w=0.5$ as the kernel hyper-parameter.

The number of kernels needed to obtain the same MSE is compared over $1000$ realizations of each signal between ($\text{PII}^{\prime}$) and KOMP. When the number of samples is at least $30$, our method is able to find a sparser representation $100\%$ of the times. In the cases with fewer samples the problem is likely undersampled, such that the estimation of the function is more difficult. Figure 6 shows two cases in which $20$ samples are simulated, where $m=5$ and $m=10$. In both cases, our method finds sparser representations in $99\%$ of the realizations. When $10$ kernels and $5$ kernels are superimposed, our method finds a representation which is less sparse in only $0.4\%$ and $0.3\%$ of realization respectively. Lastly, when the signal is a weighted sum of 5 functions and only 10 samples are generated, our method cannot find a sparser solution for $3.7\%$ of the realizations.

The generalization MSE was compared between the two methods for different levels of sparsity. Figure 7 shows the changes in generalization MSE as the number of kernels used in the representation increases. $1000$ realizations of a signal with m = $10$ and a training set of size N = $100$ were used. The ability of our method to place kernels at any location, beyond the training set, allows it to achieve significantly lower errors compared to KOMP at any sparsity level. As the number of kernels used increases, the difference in performance between the two methods decreases. At approximately 25 kernels the performance of our method plateaus. Comparatively, KOMP achieves a plateau when the representation holds 50 kernels.

V-B Varying Degrees of Smoothness

In the previous sections we have only considered signals from functions belonging to an RKHS in the family of RKHSs with Gaussian kernels. In this section we explore the effect of sample size on the complexity of the representation and the MSE on a signal of varying degrees of smoothness. To this end, a signal of varying smoothness is simulated using the following equation:

[TABLE]

where $\xi_{i}\in\mathcal{N}(0,10^{-3})$ represents the noise.

The solution of problem (PII) was compared to destructive KOMP, with the stopping criteria set to be the desired number of kernels rather than the distance from the original function. This stopping criteria allows us to have a fair comparison between our method and KOMP by using equally sparse functions. The problem in (PII) was solved using $\gamma=2$, $\eta_{\lambda}=0.001$, $\eta_{\mu}=30$ and $T=1000$. Sample sizes of 51, 101, 201, 301, 401, and 501 were created by uniformly sampling in the interval $[-5,5]$. Test sets of 1000 samples randomly selected on the interval $[-5,5]$ were created. Using the method of selecting kernel centers and widths by selecting the peaks of the function $\alpha({\mathbf{z}},w)$, our method finds a representation with 26 kernels regardless of the sample size. It can be seen in Figure 8 that in addition to the number of kernels being consistent across all sample sizes the MSE is also consistent for our method. The MSE of the estimation using KOMP, however, increases as the sample size grows.

The problem of reducing features is a combinatorial problem which grows exponentially with the sample size. The backwards approach used by KOMP is a greedy approach which removes only one kernel at a time. As the sample size increases, there are more misleading paths of removal it can take. Additionally, it is only using kernels placed at the sample points, which means it will need more kernels when the true kernel is centered between two sample points.

V-C User Localization Problem

In the remainder of this section we will apply our method to real world application for which the class of functions the signal belongs to is unknown. We consider the problem of using RF signals to identify the location of a receiver. Specifically, given the Wi-Fi signal strength from seven routers we wish to identify the room in which our receiver is located [51]. The signal strength varies depending on the location of the router. The signal was received from 7 routers spread throughout an office building. The data was collected using an Android device. At each location the signal strength from each router was observed at $1s$ intervals. The data was then categorized into 4 groups, each representing the room in which the signal strength was observed. All the rooms are on the same floor, with the rooms representing the conference room, the kitchen, the indoor sports room, and work areas [51]. The goal is to be able to accurately detect the location of the android device given the measured signal strength.

We use 10-fold cross-validation in order to estimate the generalization error of our algorithm as was used in [51]. The dataset was split into 10 sets of equal size with equal distribution of each label. At each turn one of the sets was used for testing while the others were concatenated and used to train the algorithm. This multiclass classification problem was solved using the one-vs-one strategy, which required 6 comparisons. The final class assignment is made through voting. Each comparison makes a prediction on the class of a sample and thus casts a vote for a particular class. The class with the majority of votes is assigned to the sample. The cost function was for this classification problem is

[TABLE]

Solving problem (PII) we obtain an average accuracy of $98\%$, similar to the performance observed in [51], in which a fuzzy decision tree algorithm with 50 rules was used to obtain an accuracy of $96.65\%$. This result has been observed to be consistent over increasing values of the sparsity parameter $\gamma$.

V-D Mnist Digits Classification

We use data of handwritten digits from the MNIST data set [52], which consists of a training set of $60,000$ sample-label pairs and a testing set of $10,000$ images and labels. Each sample is a 28-pixel by 28-pixel grayscale image, which was vectorized to form 784 dimensional features. The labels are between 0 and 9 and correspond to the digit written. There are a total of 10 classes.

The number of features is too large to estimate the value of $\alpha({\mathbf{z}},w)$ at every ${\mathbf{z}}\in{\mathbb{R}}^{784}$. In order to find a set ${\mathcal{X}}$ over which $\alpha({\mathbf{z}},w)$ is defined, we use k-means with $400$ clusters for each digit. Then ${\mathcal{X}}^{\prime}$ in ($\text{PII}^{\prime\prime}$) is defined as the set of all cluster centers and the cost function in (25) is used. We then run our algorithm using a one-to-one strategy for multi-class classification and achieve an accuracy of $98.12\%$ for an average of 788 features per classification which is comparable to the accuracy found using ($\text{PI}^{\prime}$) using the training set as kernel centers and ($\text{PI}^{\prime}$) using the centers found through k-means. The complexity of the representation can be further be reduced, however it comes at the cost of the classification accuracy.

Although the dimensionality of the features in the original data makes the use of (PII) impractical, we can solve that problem, by projecting the data into a lower dimensional space by using principal component analysis (PCA). The formulation in (PII) has the advantage that the found kernel centers can give some intuition about the distribution of the signal. Particularly, in the case of digits they can describe digits which are representative of written digits. To illustrate that we have performed the classification of the digits ’0’ and ’1’ using the first 3 principal components. The low dimensional feature set allows us to find the ${\mathbf{x}}$ which result in the highest value for $\alpha(w,{\mathbf{x}})$. From these points we can reconstruct the corresponding digits. Figure 9 shows the resulting images. These are not part of the initial written digit data set but rather represent an image that is closest to all written digit. The accuracy of the classification is $99.62\%$.

VI Conclusions

In this paper, we have introduced a method for function estimation in RKHS which can model signals of varying degrees of smoothness. The algorithm finds a sparse representation of a function in the sum space of a family of RKHSs, and determines the kernel parameter for each kernel. Additionally, due to the sparsity, traditional representer theorems no longer hold, so our algorithm placed kernels arbitrarily. While the problem was not convex, a change in the representation of the function to the integral over the product of the coefficient function and the kernel function, allowed us to solve the problem in the dual domain. By leveraging the results on strong duality, we were able to solve the problem in the dual domain.

The theoretical results were validated though simulated signals. We showed that our algorithm finds kernel centers and widths which can represent the function. Furthermore, we showed that our algorithm is able to find sparser representations than KOMP, with the same error. The sparseness of the representation was shown to be independent of sample size, which is not true for greedy kernel reduction methods. We also validated our method on a localization dataset, for which we were able to reduce the complexity by $86\%$ while maintaining the high accuracy. Similarly, a sparse kernel representation was obtained for classifying the digits from the MNIST dataset.

Appendix A Proof of Theorem 1

Proof.

In order to show strong duality, it is sufficient to show that the perturbed function $P(\xi)$ in (26) is convex [53, 54]. Consider the perturbed version of the optimization function

[TABLE]

In equation (26) $f_{0}(\alpha)$ represents the objective function of our original problem: $f_{0}(\alpha)=\gamma\int\mathbb{I}(\alpha({\mathbf{z}},w)\neq 0)+0.5\alpha^{2}(w,{\mathbf{x}})\,dw\,d{\mathbf{z}}$. The optimal solution for $P(0)$ is the solution to the primal problem.

Convexity of the perturbed problem can be shown, by proving that given an arbitrary pair of perturbations $\xi_{1}$ and $\xi_{2}$ and the corresponding optimal values $P(\xi_{1})$ and $P(\xi_{2})$, for any $\beta\in[0,1]$ the solution $P(\xi_{\beta})$ has the following property, where $\xi_{\beta}$ is defined by $\xi_{\beta}=\beta\xi_{1}+(1-\beta)\xi_{2}$:

[TABLE]

In order to prove the convexity of the perturbed problem we need to introduce the following lemma.

Lemma 2**.**

The set of constraints given by

[TABLE]

is convex.

Proof.

Given an arbitrary pair $b_{1},b_{2}\in\mathcal{B}$, there exists a corresponding $\alpha_{1},\alpha_{2}\in L_{2}$ such that $b_{1}=f_{0}(\alpha_{1})$ and $b_{2}=f_{0}(\alpha_{2})$. In order to prove the convexity of the set, we will show that there exists a feasible $\alpha_{\beta}\in L_{2}$ such that for any $\beta\in[0,1]$

[TABLE]

Let $\mathbb{B}$ be the Borel field of all possible subsets of $\mathcal{U}$, where $\mathcal{U}=\{\mathcal{X}\times\mathcal{W}\}$ is the set of all possible kernel centers and kernel widths and. Let us construct a measure over $\mathbb{B}$, where $\mathcal{V}\subset\mathbb{B}$.

[TABLE]

The first $2N$ elements of the measure represent the estimated function of the signal ${\mathbf{y}}$ using $\alpha_{1}$ and $\alpha_{2}$ and a subset of the kernels, where ${\mathbf{k}}({\mathbf{v}})_{i}=k({\mathbf{X}},{\mathbf{z}}_{v};w_{v})$ and ${\mathbf{v}}=[{\mathbf{z}}_{v}^{T},w_{v}]^{T}$. The last two elements of $\bm{m}(\mathcal{V})$ measure the sparsity of functions $\alpha_{1}$ and $\alpha_{2}$ over the set $\mathcal{V}$ respectively. Two sets are of interest, the empty set and $\mathcal{U}$. The measure of the former is $\bm{m}(\emptyset)=0$ and the measure for $\mathcal{U}$ can be inferred from our optimization problem.

[TABLE]

Lyapunov’s convexity theorem [55] states that a non-atomic measure vector on a Borel field is convex. Note that the representation in (5) allows us to construct the measure with non-atomic masses and is essential to the proof of strong duality. Since $\alpha$ does not contain any point masses our measure $\bm{m}$ is convex . Therefore, for any $\beta\in[0,1]$, there exists a set $\mathcal{V}_{\beta}\subset\mathbb{B}$ such that:

[TABLE]

The measure of the complement of the set $\mathcal{V}_{\beta}$, as defined by $\mathcal{V}_{\beta}^{c}\cup\mathcal{V}_{\beta}=\mathcal{U}$ and $\mathcal{V}_{\beta}^{c}\cap\mathcal{V}_{\beta}=\emptyset$, can be computed, due to the additivity property of measures:

[TABLE]

We can define the function $\alpha_{\beta}$ from (32) and (33):

[TABLE]

From this construction of $\alpha_{\beta}({\mathbf{v}})$ it can be easily seen that $f_{0}(\alpha_{\beta})=\beta f_{0}(\alpha_{1})+(1-\beta)f_{0}(\alpha_{2})$. Next we will show that $\alpha_{\beta}$ is feasible. Define $\hat{{\mathbf{y}}}_{\beta}$ as:

[TABLE]

Since $c(z_{i},y_{i})$ is convex it follows that for any $i\in[1,N]$, $c(\beta z_{i,1}+(1-\beta)z_{i,2},y_{i})\leq\beta c(z_{i,1},y_{i})+(1-\beta)c(z_{i,2},y_{i})$. We can use this property to show that $c(z_{i,\beta},y_{i})\leq\xi_{\beta}$.

[TABLE]

Thus it was proven that $\alpha_{\beta}$ is also feasible and therefore the set of constraints is convex. ∎

Let $(\alpha_{1},{\mathbf{z}}_{1},\xi_{1})$ and $(\alpha_{2},{\mathbf{z}}_{2},\xi_{2})$ be the pair of optimal solutions to the two perturbed problems $P(\xi_{1})$ and $P(\xi_{2})$. We have shown that there exists a feasible point $\alpha_{\beta}$ for the problem perturbed by $\xi_{\beta}=\beta\xi+(1-\beta)\xi^{\prime}$, which satisfies $f_{0}(\alpha_{\beta})=P(\xi_{1})+(1-\beta)P(\xi_{2})$. Given that it is a feasible point the objective function is greater or equal to the solution of the problem

[TABLE]

Since the perturbed problem is convex, the original problem has zero duality gap.

∎

Appendix B Proof of Corollary 1

Proof.

Theorem 1 implies that any solution $(\alpha^{\star},\bm{\hat{y}}^{\star})$ of (PII) is such that $(\alpha^{\star},\bm{\hat{y}}^{\star})\in\operatorname*{argmin}_{\alpha,\,\hat{y}_{i}}{\mathcal{L}}(\alpha,\bm{\hat{y}},\bm{\lambda}^{\star},\bm{\mu}^{\star})$ [43]. Since ${\mathcal{L}}$ in (8) separates across $\alpha$ and $\bm{\hat{y}}$, we can consider the minimizations individually to obtain

[TABLE]

for ${\mathcal{L}}_{\alpha}$ from (12). We also know from Proposition 2 that $\alpha_{d}^{\star}$ is in the $\operatorname*{argmin}$ set of (38). In the sequel, we show that it is (essentially) its only element.

To do so, we construct $\alpha^{\star}$ piece by piece by partitioning the integral in ${\mathcal{L}}_{\alpha}$ into three disjoint sets depending on the value of $\alpha_{d}^{\star}$. On ${\mathcal{A}}_{>}=\{({\mathbf{z}},w)\in{\mathcal{X}}\times{\mathcal{W}}\mid|\alpha_{d}^{\star}({\mathbf{z}},w)|>\sqrt{2\gamma}\}$, we know that $\alpha_{d}^{\star}$ takes values from (17), the unique minimizer of ${\mathcal{L}}_{\alpha}$ since it stems from the minimization of the strongly convex function (16). Moreover, our assumptions on the reproducing kernel together with (17) imply that $\alpha_{d}^{\star}\in L_{2}$ when restricted to ${\mathcal{A}}_{>}$. Hence,

[TABLE]

Over the set ${\mathcal{A}}_{<}=\{({\mathbf{z}},w)\in{\mathcal{X}}\times{\mathcal{W}}\mid|\alpha_{d}^{\star}({\mathbf{z}},w)|<\sqrt{2\gamma}\}$, notice from (12) that the integrand of ${\mathcal{L}}_{\alpha}$ is always non-negative. What is more, it is always strictly positive unless $\alpha\equiv 0$. This is ready by applying Lemma (1) and the fact that the minimum of (16) is positive. Thus, from the monotonicity of the integral operator, $\alpha^{\star}$ is again unique and equal to zero on ${\mathcal{A}}_{<}$. From (13), so is $\alpha_{d}^{\star}$ and we obtain

[TABLE]

Immediately, we have that $\alpha^{\star}\in L_{2}$ over ${\mathcal{A}}_{>}\cup{\mathcal{A}}_{<}$.

To conclude the proof, observe that $\mathfrak{m}\left[{\mathcal{A}}_{>}\cup{\mathcal{A}}_{<}\right]=\mathfrak{m}\left[{\mathcal{X}}\times{\mathcal{W}}\right]$, where $\mathfrak{m}$ denotes the Lebesgue measure. Indeed, the complement of ${\mathcal{A}}_{>}\cup{\mathcal{A}}_{<}$ is the set ${\mathcal{A}}_{=}=\{({\mathbf{z}},w)\in{\mathcal{X}}\times{\mathcal{W}}\mid|\alpha_{d}^{\star}({\mathbf{z}},w)|=\sqrt{2\gamma}\}$. From our assumption on the reproducing kernel, ${\mathcal{A}}_{=}$ is the set of zeros of a real analytic function, which are isolated and therefore countable [56]. In other words, $\alpha^{\star}$ and $\alpha_{d}^{\star}$ are in the same equivalence class in $L_{2}$ since they are equal except perhaps on a set of measure zero. ∎