Exploring Large Feature Spaces with Hierarchical Multiple Kernel Learning

Introduction

In the last two decades, kernel methods have been a prolific theoretical and algorithmic machine learning framework. By using appropriate regularization by Hilbertian norms, representer theorems enable to consider large and potentially infinite-dimensional feature spaces while working within an implicit feature space no larger than the number of observations. This has led to numerous works on kernel design adapted to specific data types and generic kernel-based algorithms for many learning tasks (see, e.g., ).

More precisely, we consider a positive definite kernel that can be expressed as a large sum of positive definite basis or local kernels. This exactly corresponds to the situation where a large feature space is the concatenation of smaller feature spaces, and we aim to do selection among these many kernels, which may be done through multiple kernel learning . One major difficulty however is that the number of these smaller kernels is usually exponential in the dimension of the input space and applying multiple kernel learning directly in this decomposition would be intractable.

Finally, we extend in Section 4 some of the known consistency results of the Lasso and multiple kernel learning , and give a partial answer to the model selection capabilities of our regularization framework by giving necessary and sufficient conditions for model consistency. In particular, we show that our framework is adapted to estimating consistently only the hull of the relevant variables. Hence, by restricting the statistical power of our method, we gain computational efficiency.

Hierarchical multiple kernel learning (HKL)

Our sum assumption corresponds to a situation where the feature map $\Phi(x)$ and feature space $\mathcal{F}$ for $k$ is the concatenation of the feature maps $\Phi_{v}(x)$ for each kernel $k_{v}$, i.e, $\mathcal{F}=\prod_{v\in V}\mathcal{F}_{v}$ and $\Phi(x)=(\Phi_{v}(x))_{v\in V}$. Thus, looking for a certain $\beta\in\mathcal{F}$ and a predictor function $f(x)=\langle\beta,\Phi(x)\rangle$ is equivalent to looking jointly for $\beta_{v}\in\mathcal{F}_{v}$, for all $v\in V$, and $f(x)=\sum_{v\in V}\langle\beta_{v},\Phi_{v}(x)\rangle$.

As mentioned earlier, we make the assumption that the set $V$ can be embedded into a directed acyclic graph. Directed acyclic graphs (referred to as DAGs) allow to naturally define the notions of parents, children, descendants and ancestors. Given a node $w\in V$, we denote by ${\rm A}(w)\subset V$ the set of its ancestors, and by ${\rm D}(w)\subset V$, the set of its descendants. We use the convention that any $w$ is a descendant and an ancestor of itself, i.e., $w\in{\rm A}(w)$ and $w\in{\rm D}(w)$. Moreover, for $W\subset V$, we let denote ${\rm sources}(W)$ the set of sources of the graph $G$ restricted to $W$ (i.e., nodes in $W$ with no parents belonging to $W$). Given a subset of nodes $W\subset V$, we can define the hull of $W$ as the union of all ancestors of $w\in W$, i.e., ${\rm hull}(W)=\bigcup_{w\in W}A(w)$. Given a set $W$, we define the set of extreme points of $W$ as the smallest subset $T\subset W$ such that ${\rm hull}(T)={\rm hull}(W)$ (note that it is always well defined, as $\bigcap_{T\subset V,\ {\rm hull}(T)={\rm hull}(W)}T$). See Figure 1 for examples of these notions.

The goal of this paper is to perform kernel selection among the kernels $k_{v}$, $v\in V$. We essentially use the graph to limit the search to specific subsets of $V$. Namely, instead of considering all possible subsets of active (relevant) vertices, we are only interested in estimating correctly the hull of these relevant vertices; in Section 2.2, we design a specific sparsity-inducing norms adapted to hulls.

In this paper, we primarily focus on kernels that can be expressed as “products of sums”, and on the associated $p$-dimensional directed grids, while noting that our framework is applicable to many other kernels. Namely, we assume that the input space $\mathcal{X}$ factorizes into $p$ components $\mathcal{X}=\mathcal{X}_{1}\times\cdots\times\mathcal{X}_{p}$ and that we are given $p$ sequences of length $q+1$ of kernels $k_{ij}(x_{i},x_{i}^{\prime})$, $i\in\{1,\dots,p\}$, $j\in\{0,\dots,q\}$, such that $k(x,x^{\prime})=\sum_{j_{1},\dots,j_{p}=0}^{q}\prod_{i=1}^{p}k_{ij_{i}}(x_{i},x_{i}^{\prime})=\prod_{i=1}^{p}\left(\sum_{j_{i}=0}^{q}k_{ij_{i}}(x_{i},x_{i}^{\prime})\right)$. We thus have a sum of $(q+1)^{p}$ kernels, that can be computed efficiently as a product of $p$ sums. A natural DAG on $V=\prod_{i=1}^{p}\{0,\dots,q\}$ is defined by connecting each $(j_{1},\dots,j_{p})$ to $(j_{1}\!+\!1,j_{2},\dots,j_{p})$, $\dots$, $(j_{1},\dots,j_{p-1},j_{p}\!+\!1)$. As shown in Section 2.2, this DAG will correspond to the constraint of selecting a given product of kernels only after all the subproducts are selected. Those DAGs are especially suited to nonlinear variable selection, in particular with the polynomial and Gaussian kernels. In this context, products of kernels correspond to interactions between certain variables, and our DAG implies that we select an interaction only after all sub-interactions were already selected.

where $c^{2}=a^{2}+2ab$, $A=a+b+c$, and $H_{k}$ is the $k$-th Hermite polynomial. By appropriately truncating the sum, i.e, by considering that the first $q$ basis kernels are obtained from the first $q$ single Hermite polynomials, and the $(q+1)$-th kernel is summing over all other kernels, we obtain a decomposition of a uni-dimensional Gaussian kernel into $q+1$ components ($q$ of them are one-dimensional, the last one is infinite-dimensional, but can be computed by differencing). The decomposition ends up being close to a polynomial kernel of infinite degree, modulated by an exponential . One may also use an adaptive decomposition using kernel PCA (see, e.g., ), which is equivalent to using the eigenvectors of the empirical covariance operator associated with the data (and not the population one associated with the Gaussian distribution with same variance). In simulations, we tried both with no significant differences.

Finally, by taking product over all variables, we obtain a decomposition of the $p$-dimensional Gaussian kernel into $(q+1)^{p}$ components, that are adapted to nonlinear variable selection. Note that for $q=1$, we obtain ANOVA-like decompositions .

2 Graph-based structured regularization

Given $\beta\in\prod_{v\in V}\mathcal{F}_{v}$, the natural Hilbertian norm $\|\beta\|$ is defined through $\|\beta\|^{2}=\sum_{v\in V}\|\beta_{v}\|^{2}$. Penalizing with this norm is efficient because summing all kernels $k_{v}$ is assumed feasible in polynomial time and we can bring to bear the usual kernel machinery; however, it does not lead to sparse solutions, where many $\beta_{v}$ will be exactly equal to zero.

As said earlier, we are only interested in the hull of the selected elements $\beta_{v}\in\mathcal{F}_{v}$, $v\in V$; the hull of a set $I$ is characterized by the set of $v$, such that ${\rm D}(v)\subset I^{c}$, i.e., such that all descendants of $v$ are in the complement $I^{c}$: ${\rm hull}(I)=\{v\in V,{\rm D}(v)\subset I^{c}\}^{c}$. Thus, if we try to estimate ${\rm hull}(I)$, we need to determine which $v\in V$ are such that ${\rm D}(v)\subset I^{c}$. In our context, we are hence looking at selecting vertices $v\in V$ for which $\beta_{{\rm D}(v)}=(\beta_{w})_{w\in{\rm D}(v)}=0$.

where $(d_{v})_{v\in V}$ are positive weights. Penalizing by such a norm will indeed impose that some of the vectors $\beta_{{\rm D}(v)}\in\prod_{w\in{\rm D}(v)}\mathcal{F}_{w}$ are exactly zero. We thus consider the following minimization problemFollowing , we consider the square of the norm, which does not change the regularization properties, but allow simple links with multiple kernel learning.:

Finally, note that in certain settings (finite dimensional Hilbert spaces and distributions with absolutely continuous densities), these norms have the effect of selecting a given kernel only after all of its ancestors . This is another explanation why hulls end up being selected, since to include a given vertex in the models, the entire set of ancestors must also be selected.

Optimization problem

In this section, we give optimality conditions for the problems in Eq. (1), as well as optimization algorithms with polynomial time complexity in the number of selected kernels. In simulations we consider total numbers of kernels larger than $10^{30}$, and thus such efficient algorithms are essential to the success of hierarchical multiple kernel learning (HKL).

The problem in Eq. (1) is thus equivalent to

The pair $(\alpha,\eta)$ is optimal for Eq. (1), with $\forall w,\beta_{w}\!=\!\zeta_{w}\sum_{i=1}^{n}\alpha_{i}\Phi_{w}(x_{i})$, if and only if $(a)$ given $\eta$, $\alpha$ is optimal for the single kernel learning problem with kernel matrix $K=\sum_{w\in V}\zeta_{w}(\eta)K_{w}$, and $(b)$ given $\alpha$, $\eta\in H$ maximizes

Moreover, the total duality gap can be upperbounded as the sum of the two separate duality gaps for the two optimization problems, which will be useful in Section 3.2 (see Appendix for more details). Note that in the case of “flat” regular multiple kernel learning, where the DAG has no edges, we obtain back usual optimality conditions .

Following a common practice for convex sparsity problems , we will try to solve a small problem where we assume we know the set of $v$ such that $\|\beta_{{\rm D}(v)}\|$ is equal to zero (Section 3.3). We then “simply” need to check that variables in that set may indeed be left out of the solution. In the next section, we show that this can be done in polynomial time although the number of kernels to consider leaving out is exponential (Section 3.2).

2 Conditions for global optimality of reduced problem

We let denote $J$ the complement of the set of norms which are set to zero. We thus consider the optimal solution $\beta$ of the reduced problem (on $J$), namely,

with optimal primal variables $\beta_{J}$, dual variables $\alpha$ and optimal pair $(\eta_{J},\zeta_{J})$. We now consider necessary conditions and sufficient conditions for this solution (augmented with zeros for non active variables, i.e., variables in $J^{c}$) to be optimal with respect to the full problem in Eq. (1). We denote by $\delta=\sum_{v\in J}d_{v}\|\beta_{D(v)\cap J}\|$ the optimal value of the norm for the reduced problem.

If the reduced solution is optimal for the full problem in Eq. (1) and all kernels in the extreme points of $J$ are active, then we have

If $\max_{t\in{\rm sources}(J^{c})}\textstyle\sum_{w\in{\rm D}(t)}\alpha^{\top}K_{w}\alpha/(\sum_{v\in{\rm A}(w)\cap{\rm D}(t)}d_{v})^{2}\leqslant\delta^{2}+\varepsilon/\lambda$, then the total duality gap is less than $\varepsilon$.

The proof is fairly technical and can be found in the Appendix; this result constitutes the main technical contribution of the paper: it essentially allows to solve a very large optimization problem over exponentially many dimensions in polynomial time.

The necessary condition $(N_{J})$ does not cause any computational problems. However, the sufficient condition $(S_{J,\varepsilon})$ requires to sum over all descendants of the active kernels, which is impossible in practice (as shown in Section 5, we consider $V$ of cardinal often greater than $10^{30}$). Here, we need to bring to bear the specific structure of the kernel $k$. In the context of directed grids we consider in this paper, if $d_{v}$ can also be decomposed as a product, then $\sum_{v\in{\rm A}(w)\cap{\rm D}(t)}d_{v}$ is also factorized, and we can compute the sum over all $v\in{\rm D}(t)$ in linear time in $p$. Moreover we can cache the sums $\sum_{w\in{\rm D}(t)}K_{w}/(\sum_{v\in{\rm A}(w)\cap{\rm D}(t)}d_{v})^{2}$ in order to save running time.

3 Dual optimization for reduced or small problems

4 Kernel search algorithm

We are now ready to present the detailed algorithm which extends the feature search algorithm of . Note that the kernel matrices are never all needed explicitly, i.e., we only need them (a) explicitly to solve the small problems (but we need only a few of those) and (b) implicitly to compute the sufficient condition $(S_{J,\varepsilon})$, which requires to sum over all kernels, as shown in Section 3.2.

Initialization: set $J={\rm sources}(V)$, compute $(\alpha,\eta)$ solutions of Eq. (3), obtained using Section 3.3

while $(N_{J})$ and $(S_{J,\varepsilon})$ are not satisfied and $\#(V)\leqslant Q$

If $(N_{J})$ is not satisfied, add violating variables in ${\rm sources}(J^{c})$ to $J$ else, add violating variables in ${\rm sources}(J^{c})$ of $(S_{J,\varepsilon})$ to $J$

Recompute $(\alpha,\eta)$ optimal solutions of Eq. (3)

The previous algorithm will stop either when the duality gap is less than $\varepsilon$ or when the maximal number of kernels $Q$ has been reached. In practice, when the weights $d_{v}$ increase with the depth of $v$ in the DAG (which we use in simulations), the small duality gap generally occurs before we reach a problem larger than $Q$. Note that some of the iterations only increase the size of the active sets to check the sufficient condition for optimality; forgetting those does not change the solution, only the fact that we may actually know that we have an $\varepsilon$-optimal solution.

In the directed $p$-grid case, the total running time complexity is a function of the number of observations $n$, and the number $R$ of selected kernels; with proper caching, we obtain the following complexity, assuming $O(n^{3})$ for the single kernel learning problem, which is conservative: $O(n^{3}R+n^{2}Rp^{2}+n^{2}R^{2}p)$, which decomposes into solving $O(R)$ single kernel learning problems, caching $O(Rp)$ kernels, and computing $O(R^{2}p)$ quadratic forms for the sufficient conditions. Note that the kernel search algorithm is also an efficient algorithm for unstructured MKL.

Consistency conditions

As said earlier, the sparsity pattern of the solution of Eq. (1) will be equal to its hull, and thus we can only hope to obtain consistency of the hull of the pattern, which we consider in this section.

Following , we make the following assumptions on the underlying joint distribution of $(X,Y)$: (a) the joint covariance matrix $\boldsymbol{\Sigma}$ of $(\Phi(x_{v}))_{v\in V}$ (defined with appropriate blocks of size $f_{v}\times f_{w}$) is invertible, (b) $E(Y|X)=\sum_{w\in\boldsymbol{W}}\langle\boldsymbol{\beta}_{w},\Phi_{w}(x)\rangle$ with $\boldsymbol{W}\subset V$ and $\mathop{\rm var}(Y|X)=\boldsymbol{\sigma}^{2}>0$ almost surely. With these simple assumptions, we obtain (see proof in the Appendix):

then $\boldsymbol{\beta}$ and the hull of $\boldsymbol{W}$ are consistently estimated when $\lambda_{n}n^{1/2}\to\infty$ and $\lambda_{n}\to 0$.

If the $\boldsymbol{\beta}$ and the hull of $\boldsymbol{W}$ are consistently estimated for some sequence $\lambda_{n}$, then

Note that the last two propositions are not consequences of the similar results for flat MKL , because the groups that we consider are overlapping. Moreover, the last propositions show that we indeed can estimate the correct hull of the sparsity pattern if the sufficient condition is satisfied. In particular, if we can make the groups such that the between-group correlation is as small as possible, we can ensure correct hull selection. Finally, it is worth noting that if the ratios $d_{w}/\max_{v\in{\rm A}(w)}d_{v}$ tend to infinity slowly with $n$, then we always consistently estimate the depth of the hull, i.e., the optimal interaction complexity. We are currently investigating extensions to the non parametric case , in terms of pattern selection and universal consistency.

Simulations

We can see from Table 1, that HKL outperforms other methods, in particular for the datasets bank-32nm, bank-32nh, pumadyn-32nm, pumadyn-32nh, which are datasets dedicated to non linear regression. Note also, that we efficiently explore DAGs with very large numbers of vertices $\#(V)$.

For binary classification datasets, we compare HKL (with the logistic loss) to two other methods (L2, greedy) in Table 2. For some datasets (e.g., spambase), HKL works better, but for some others, in particular when the generating problem is known to be non sparse (ringnorm, twonorm), it performs slightly worse than other approaches.

Conclusion

Appendix A Optimization results

with equality if and only if $\eta_{v}=d_{v}^{-1}\|\beta_{{\rm D}(v)}\|(\sum_{v\in V}d_{v}\|\beta_{{\rm D}(v)}\|)^{-1}$.

When the DAG is a tree (i.e., when each vertex has at most one parent), then, without loss of generality we may consider that only one vertex has no parent (the root $r$) while all others $w$ have exactly one parent $\pi(w)$. In this situation, we have for all $v\neq r$, $\zeta_{\pi(v)}^{-1}-\zeta_{v}^{-1}=-\eta_{\pi(v)}^{-1}$. Moreover, for all leaves $v$, $\zeta_{v}=\eta_{v}$. This implies that the constraint $\eta\geqslant 0$ is equivalent to $\zeta\geqslant 0$ and for all $v\neq r$, $\zeta_{\pi(v)}\geqslant\zeta_{v}$. The final constraint $\sum_{v\in V}\eta_{v}d_{v}^{2}\leqslant 1$, may then be written as:

which is clearly convex . When the DAG is not a tree, we conjecture that the set $Z$ is not convex.

A.2 Fenchel conjugates

In particular, we have for the following standard examples:

for least-squares regression, we have $\varphi_{i}(a)=\frac{1}{2}(y_{i}-a)^{2}$ and $\psi_{i}(b)=\frac{1}{2}b^{2}+by_{i}$,

for logistic regression, we have $\varphi_{i}(a)=\log(1+\exp(-y_{i}a_{i}))$, where $y_{i}\in\{-1,1\}$, and $\psi_{i}(b)=(1+by_{i})\log(1+by_{i})-by_{i}\log(-by_{i})$ if $by_{i}\in$, $+\infty$ otherwise.

for support vector machine classification, we have $\varphi_{i}(a)=\max(0,1-y_{i}a)$, where $y_{i}\in\{-1,1\}$, and $\psi_{i}(b)=y_{i}b$ if $by_{i}\in$, $+\infty$ otherwise.

A.3 Preliminary propositions

and the optimal $\beta$ can be found from an optimal $\alpha$ as $\beta_{w}=\sum_{i=1}^{n}\alpha_{i}\Phi_{w}(x_{i})$.

Proof We introduce auxiliary variables $u_{i}=\sum_{v\in V}\langle\beta_{v},\Phi_{v}(x_{i})\rangle$ and consider the Lagrangian:

Minimizing with respect to the primal variables $u,\beta$, we get the dual problem.

We will use the following simple result, which implies that each component $\zeta_{w}(\eta)$ is a concave function of $\eta$:

The minimum of $\sum_{j=1}^{m}a_{j}x_{j}^{2}$ subject to $\sum_{j=1}^{m}x_{j}=1$ is equal to $\left(\sum_{j=1}^{m}a_{i}^{-1}\right)^{-1}$ and is attained at $x_{i}=a_{i}^{-1}\left(\sum_{j=1}^{m}a_{i}^{-1}\right)^{-1}$.

The following proposition derives the dual of the problem in $\eta$:

which can be minimized in closed form with respect to $\delta^{2}$ and $\kappa\in L$, and leads to (using Lemma 1):

A.4 Duality gaps

This function is convex in $\eta$ (because of Lemma 1) and concave in $\alpha$, standard arguments (e.g., primal and dual strict feasibilities) show that there is no duality gap to the variational problems:

We can decompose the duality gap, given a pair $(\eta,\alpha)$ as

We thus get the desired upper bound from which proposition 1 (of the main paper) follows, as well as the upper bound on the duality gap.

A.5 Necessary and sufficient conditions - truncated problem

We assume that we know the optimal solution of a truncated problem where the entire set of decendants of some nodes have been removed. We let denote $J$ the hull of the set of active variables. We now consider necessary conditions and sufficient conditions for this solution to be optimal with respect to the full problem. This will lead to Proposition 2 and 3 of the main paper.

We first use Proposition 2 of the Appendix, to get a set of $\kappa_{vw}$ for $(v,w)\in J$ for the reduced problem; the goal here is to get necessary conditions by relaxing the dual problem defining $\kappa\in L$ and find an approximate solution, while for the sufficient condition, any candidate leads to a sufficient condition. It turns out that we will use the solution of the relaxed solution required for the necessary condition for the sufficient condition.

If we assume that all variables in $J$ are indeed active, then any optimal $\kappa\in L$ must be such that $\kappa_{vw}=0$ if $v\in J$ and $w\in J^{c}$. We then let free $\kappa_{vw}$ for $v,w$ in $J$. Our goal is to find good candidates for those free dual parameters.

We first derive necessary conditions by lowerbounding the sums by maxima:

which can be minimized in closed form with respect to $\kappa$ leading to

For sufficient conditions, we simply take the value obtained before for $\kappa$, which leads to

A.6 Optimality conditions for the primal formulation

We know derive optimality conditions for the problem in the paper, which we will need in Section B, i.e.:

and thus $\beta$ if optimal if and ony if, we have, with $\delta=\sum_{v\in J}d_{v}\|\beta_{{\rm D}(v)\cap J}\|$:

Note that when regularizing by ${\lambda}\sum_{v\in V}d_{v}\|\beta_{{\rm D}(v)}\|$ instead of $\frac{\lambda}{2}\left(\sum_{v\in V}d_{v}\|\beta_{{\rm D}(v)}\|\right)^{2}$, we have the same optimality condition with $\delta=1$.

Appendix B Consistency conditions

The consistency condition is then obtained by studying when the first order expansion indeed has the correct sparsity pattern (for more precise statements and arguments, see ). We let denote $\gamma_{\boldsymbol{W}}$ the solution of the previous problem, restricted to $\gamma_{{\boldsymbol{W}}^{c}}=0$. We have:

Following the previous section, it is optimal if and only for all $\Delta\in{\boldsymbol{W}}^{c}$,

The condition for good pattern selection is that for all $\Delta\in{\boldsymbol{W}}^{c}$,

References