<figure><img src="image_1.png" /></figure>

HAL
open science

# Rethinking multiple kernel learning under the lenses of Importance Weighted Monte Carlo Variational Inference

Davide Adamo, Manon Vuillien, Emmanuelle Vila, Marco Corneli

► To cite this version:

Davide Adamo, Manon Vuillien, Emmanuelle Vila, Marco Corneli. Rethinking multiple kernel learning under the lenses of Importance Weighted Monte Carlo Variational Inference. 2025. hal-05123496v1

HAL Id: hal-05123496

https://hal.science/hal-05123496v1

Preprint submitted on 20 Jun 2025 (v1), last revised 15 Jan 2026 (v2)

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rethinking multiple kernel learning under the lenses of Importance Weighted Monte Carlo Variational Inference

Davide Adamo¹,²*, Manon Vuillien¹, Emmanuelle Vila³, Marco Corneli¹,²

¹Université Côte d'Azur, UMR 7264 CEPAM, CNRS, Nice, France.

²Université Côte d'Azur, Inria, CNRS, Laboratoire J.A.Dieudonné, Maasai team, Nice, France.

³Université Lumière Lyon II, UMR 5133 Archéorient, CNRS, Lyon, France.

*Corresponding author(s). E-mail(s):
Davide.ADAMO@univ-cotedazur.fr;

## Abstract

Kernel methods have been widely used in machine learning as they are a powerful tool for implicitly mapping data into high-dimensional spaces, enabling the discovery of complex patterns that might be challenging to capture in the original feature space. Although some classification and regression problems can be successfully addressed with a single kernel, sometimes real-world scenarios exhibit complex structures, and it is desirable to employ several kernel types, one for each notion of similarity that one aims to take into account. This is where multiple kernel learning (MKL) comes into play. This paper revisits multi-kernel classification with a specific focus on kernel(s) selection in the light of recent developments in Monte Carlo (importance weighted) variational inference (MCVI). In the framework of kernelized logistic regression (KLR), we consider positive semidefinite linear combinations of kernels and treat the kernel weights as random variables. Proper choices of prior distributions coupled with the explicit derivation of the importance-weighted lower bound (IW-ELBO), generalization of the traditional variational lower bound (ELBO), allow us to both perform kernel selection via
shrinking and to perform posterior inference on the kernel weights, without need-
ing MCMC sampling. Unlike pure optimization-based approaches to MKL, our
optimization problem is fully unconstrained and can be optimized by standard
stochastic gradient descent.

**Keywords:** Multiple kernel learning, Monte Carlo variational inference, kernel(s) selection, importance-weighted lower bound

# 1 Introduction

## 1.1 Context

Kernel methods provide a powerful framework for capturing complex relationships in data by implicitly mapping input features into higher-dimensional spaces, the reproducing kernel Hilbert spaces (RKHS). In the RKHS, linear models for regression or classification can capture non-linear relations in the features space between the response variables and the covariates. The key ingredient in kernel methods is the kernel itself (Vapnik et al. 1998), a symmetric positive semi-definite function that determines the similarity between pairs of data points.

While a single kernel can capture specific types of relationships within the data,
in practice, data may exhibit multiple and complex structures that no single kernel
can fully describe. This limitation led to the development of Multiple Kernel Learning
(Gönen and Alpaydn 2011, MKL), whose goal is to exploit information from several
kernels and combining them into a unified model. Additionally, MKL approaches pro-
vide solutions to a number of standard challenges, such as hyper-parameters tuning
(even for a single kernel) and kernel selection, usually handled via cross-validation
(Stone 1974). Due to such flexibility, MKL methods have been effectively applied
in various domains. For instance, dealing with object recognition in images, Bucak
et al. (2013) represented each data using multiple feature sets and then solved an
MKL problem to demonstrate that combining multiple kernels outperformed single-kernel strategies. In bioinformatics, MKL has been used to integrate information from different sources, such as sequence data, expression data, and protein structures, in order to improve the prediction of biological outcomes like protein-protein interactions (Ben-Hur and Noble 2005).

In the last two decades, several approaches to MKL have been introduced in the literature. Here, we focus on two of them: optimization-based and fully Bayesian approaches.

**Optimization-based approaches.** One of the foundational works in multiple kernel learning is Simple MKL, introduced in Rakotomamonjy et al. (2008) and combining several kernels in a support vector machine (SVM). It employs a reduced gradient descent approach to solve a constrained optimization problem handling high-dimensional data and involves a strategic selection of a subset of pertinent kernels, leading to an improvement in predictive accuracy while preserving the advantages associated with integrating multiple kernels. This method tends to produce sparse kernel combinations ($l_1$ penalty on kernel weights). However, being a pure optimization approach, any form of inference on the kernel weights is prevented in Simple MKL. Bach (2008) proposed an MKL method based on group lasso. This method applies group lasso regularization to the kernel weights, promoting sparsity between groups of kernels rather than individual kernels. This technique has proven particularly effective in scenarios where kernels were naturally grouped, such as in hierarchical feature representations.

The methods mentioned above have been formulated for binary problems. Zien and Ong (2007) and Gehler and Nowozin (2009) provided extensions of MKL to multi-class problems.

Kloft et al. (2011) extended the idea of sparse MKL by proposing a general framework that combined different types of regularization (e.g., $l_1$, $l_2$) to balance sparsity and
smoothness in the learned kernel weights. This approach allowed for greater flexibility
in modeling and better performances using non-sparsity solutions.

**Fully Bayesian approaches.** While optimization-based methods play a central role in the landscape of MKL, fully Bayesian approaches emerged as a powerful alternative. These methods framed MKL within a probabilistic framework, treating the kernel weights as random variables with associated prior distributions. The Bayesian approach offered several advantages, including the ability to incorporate prior knowledge and quantify uncertainty in the learned kernel weights. Girolami and Rogers (2005) introduced a Bayesian MKL approach for binary classification. This method employs hierarchical models to represent the uncertainty in kernel weights, providing a probabilistic interpretation of the model. Damoulas and Girolami (2008) extended this framework to a multi-class formulation, enforcing a convex sum of input kernels using a Dirichlet prior on the kernel weights. Although powerful, these methods were computationally intensive and require sampling techniques for inference (e.g. MCMC). To address the computational cost challenges of training these algorithms, Gönen (2012) proposed the Bayesian Efficient MKL (BEMKL), a fully conjugate probabilistic model, and developed a deterministic variational approximation mechanism for inference. Thanks to variational inference this model is more suited for larger datasets than its Bayesian competitors. However, exact posterior inference is lost.

**Recent advances.** With the advent of deep learning, deep architectures offered a natural alternative to traditional (M)KL approaches. However, a number of works tried to extend MKL in the light of deep learning, sometimes by using kernel SVM as the last layer of deep architectures (e.g. Huang and LeCun 2006; Tang 2013), sometimes by combining kernels at different layers of deep architectures (convolutional kernel networks, Mairal et al. 2014; Mairal 2016; Mohammadnia-Qaraei et al. 2018), Other times kernels are built by non linearly transforming the feature vectors before applying
(Kingma et al. 2013, ELBO) that relies on importance sampling and allows us to sample from (a distribution converging to) the posterior distribution of the kernel weights easily and fast. The hybrid framework we introduce is referred to as IW-MKL.

Being suitable for large scale problems our approach aims to capture the best of two worlds, i.e. the simplicity of the optimization based methods in terms of formulation, possibly improved by the fact that our optimization problem is fully unconstrained; and the flexibility of Bayesian methods, allowing us to perform fast posterior inference, without needing MCMC sampling.

## 1.3 Organization of the article

The article is structured as follows. Section 2 presents the logistic regression model with multiple kernels and hybrid assumptions on the model parameters. In section 3 we detail our approach using Monte Carlo variational inference, extending the variational lower bound to the importance-weighted lower bound in two different distributional scenarios. Here, we also propose a straightforward and interpretable method for kernel selection. In section 4 we provide examples to assess our approach, including simulated data and comparisons with traditional MKL methods. In an archaeozoological context, section 5 presents a concrete application of our method for species classification. Finally, section 6 offers a discussion and some concluding remarks.

# 2 Multiple kernel logistic regression (MKLR)

Given a set of $N$ training data, $\{(x_i, y_i)\}_{i=1}^N$, $x_i$ represents a $P$-dimensional input (or feature) vector and $y_i \in \{0, 1\}$ its corresponding label. For simplicity we consider a binary classification problem, however, extension to the multi-class classification is straightforward. Let $D \in \mathbb{N}$ be the number of different kernel sources corresponding to different feature maps or, equivalently, kernel functions. We define a set of kernel matrices $K := \{K^{(1)}, K^{(2)}, \dots, K^{(D)}\}$, where each kernel matrix is of size $N \times N$, and
defined by

$$K^{(d)} := [k^{(d)}(x_i, x_j)]_{i,j=1}^{N},$$

with $k^{(d)}(\cdot, \cdot)$ being a positive semi-definite kernel function. Note that, these kernel matrices can either arise from different kernel functions (e.g. linear, polynomial or Gaussian kernels) or from the same function with different parameterization (for instance, different values of the scale parameter in the Gaussian or radial basis function kernel).

## 2.1 Frequentist model statement

In order to define a linear combination of kernels still being a kernel, let us define a new matrix $K \in \mathbb{R}^{N \times N}$

$$K := \sum_{d=1}^{D} \beta_d K^{(d)} \quad (1)$$

where $\beta_1, \dots, \beta_D$ are assumed to be non-negative weights that control the contribution of each kernel matrix. Although it is common in the literature to constrain $\beta_1, \dots, \beta_D$ to lie within the $D-1$ simplex (i.e., $\sum_{d=1}^{D} \beta_d = 1$ and $\beta_d \ge 0$), this assumption is not needed in order for $K$ to be a kernel. In this paper we treat both the scenarios where $\beta$ belongs to the simplex or not.

Let $K_i^{(d)}$ denote the $i$-th row of the matrix $K^{(d)}$. The $i$-th observation $y_i$ is then described by the row vectors $K_i^{(d)}$ for $d \in \{1, \dots, D\}$. We see $y_i$ as the outcome of a Bernoulli random variable $Y_i$, whose probability of success is denoted by $p_i$. Assuming the independence of the $N$ random variables $Y_1, \dots, Y_N$, the model is specified via the following equation

$$\log \left( \frac{p_i}{1 - p_i} \right) = K_i \alpha = \left[ \sum_{d=1}^{D} \beta_d K_i^{(d)} \right] \alpha,$$
where $\alpha$ is an unknown parameter in $\mathbb{R}^N$. The right-hand side of the above equation can be reformulated as follows:

$$ \left[ \sum_{d=1}^{D} \beta_d K_i^{(d)} \right] \alpha = \sum_{n=1}^{N} \sum_{d=1}^{D} \beta_d K_{in}^{(d)} \alpha_n = \beta^T L_i \alpha, $$

where $\beta := [\beta_1, \dots, \beta_D]^T \in \mathbb{R}_+^D$ and $L_i := [K_i^{(1)}, \dots, K_i^{(D)}]^T \in \mathbb{R}^{D \times N}$ is a matrix with $K_i^{(d)}$ as its rows. The likelihood of the observed data is

$$ p(y|\mathcal{K}, \alpha, \beta) = \prod_{i=1}^{N} p_i^{y_i} (1 - p_i)^{1-y_i}, $$

where $y := (y_1, \dots, y_N)$ and the probabilities $p_i$ are given by:

$$ p_i = \frac{e^{\beta^T L_i \alpha}}{1 + e^{\beta^T L_i \alpha}}. $$

The log-likelihood function

$$ \begin{aligned} \log p(y|\mathcal{K}, \alpha, \beta) &= \sum_{i=1}^{N} \left[ y_i \log \left( \frac{p_i}{1-p_i} \right) + \log(1-p_i) \right] = \\ &= \sum_{i=1}^{N} \left[ y_i \beta^T L_i \alpha - \log \left( 1 + e^{\beta^T L_i \alpha} \right) \right] \end{aligned} \quad (2) $$

could be directly maximized with respect to $(\alpha, \beta)$. However, this approach leads to a somehow challenging constrained optimization problem, due to the requirement that the kernel weights $\beta$ are either non-negative or in the simplex. More important, a direct maximization would only provide us with point estimates of the parameters, offering limited insights into their variability or uncertainty. Fortunately, a Bayesian view of the model allows us to overcome both this issues.
## 2.2 Bayesian model statement

We now see the logistic regression parameter $\alpha$ as a random variable drawn from a Gaussian distribution $N(0, (\lambda K)^{-1})$, where $\lambda$ is a user defined parameter controlling the strength of the prior. Thus:

$$ \log p(\alpha | \lambda, K) = -\frac{\lambda}{2} \alpha^T K \alpha + \frac{1}{2} \log (|\lambda K|) - \frac{N}{2} \log 2\pi. \quad (3) $$

Consequently, the log-likelihood of the complete data becomes

$$ \begin{aligned} l_R(\alpha, \beta) &:= \log p(y, \alpha | \mathcal{K}, \beta) = \log p(y | \mathcal{K}, \alpha, \beta) + \log p(\alpha | \lambda, K) \\ &= \sum_{i=1}^{N} \left[ y_i \beta^T L_i \alpha - \log \left( 1 + e^{\beta^T L_i \alpha} \right) \right] - \frac{\lambda}{2} \alpha^T K \alpha + \frac{1}{2} \log (|\lambda K|) - \frac{N}{2} \log 2\pi. \end{aligned} \quad (4) $$

Knowing that under our assumptions $K$ in Eq. (1) is a positive semi-definite kernel matrix, the log-likelihood in Eq. (4) corresponds to a kernelized LR via the Representer Theorem (Schölkopf et al. 2001), with the $l_2$ penalty in the feature space translating into the quadratic form $\alpha^T K \alpha$. The last two terms in the above equation allow us to have a penalized log-likelihood $l_R(\alpha, \beta)$ which is a well defined joint probability and we will need it in Section 3. As for the term $\frac{1}{2} \log (|\lambda K|)$, as usual in the literature (see e.g. Bishop 1991, Section 6.4.5), we introduce a diagonal perturbation of the matrix $K$ matrix in order to make it non-singular.

We now introduce a prior distribution over the kernel weights $\beta$. Although several probability distributions might be taken into account, here we only inspect two simple solutions: an Exponential prior, in order to preserve the positivity of each kernel weight and a symmetric Dirichlet prior in order put the kernel weights onto the $D-1$ simplex.
### Exponential prior.

Let us assume that each $\beta_d$ is independently and identically distributed according to an Exponential distribution with parameter $\delta$:

$$p(\beta|\delta) = \prod_{d=1}^{D} p(\beta_d|\delta) = \delta^D e^{-\delta \sum_{d=1}^{D} \beta_d}.$$

This leads to the following expression for the log-prior:

$$\log p(\beta|\delta) = D \log \delta - \delta \sum_{d=1}^{D} \beta_d = D \log \delta - \delta \|\beta\|_1,$$

where the last equality comes from the fact that each $\beta_d$ is non-negative. Thus the log-likelihood of the complete data $(y, \alpha, \beta)$ reads:

$$\log p(y, \alpha, \beta|K, \delta) = l_R(\alpha, \beta) + \log p(\beta|\delta) = l_R(\alpha, \beta) - \delta \|\beta\|_1 + D \log \delta. \quad (5)$$

As it can be seen a Lasso penalty naturally arises on $\beta$, making this objective function interesting to optimize since it allows to select the relevant kernels via shrinking. Similar optimization problems have been explored in the context of SVMs (Rakotomamonjy et al. 2008) and LR (He et al. 2021) from a purely optimization perspective.

### Dirichlet prior.

As an alternative to the Exponential prior and in order to constrain $\beta$ to lie within the $D-1$ simplex, let us assume that $\beta$ is sampled from a symmetric Dirichlet distribution with positive parameter $\delta$:

$$p(\beta|\delta) = \prod_{d=1}^{D} p(\beta_d|\delta) = \frac{1}{B(\delta)} \prod_{d=1}^{D} (\beta_d)^{\delta-1},$$
where $B(\cdot)$ denotes the beta function. In the above equation, with an abuse of notation $\delta$ represents both the user-defined scalar and the D-dimensional vector of $(\delta, \dots, \delta)$ when evaluated in the beta function. Thus:

$$ \log p(\beta|\delta) = (\delta - 1) \sum_{d=1}^{D} \log(\beta_d) - \log(B(\delta)). $$

The log-likelihood of the complete data $(y, \alpha, \beta)$ reads:

$$ \begin{aligned} \log p(y, \alpha, \beta|\mathcal{K}, \delta) &= l_R(\alpha, \beta) + \log p(\beta|\delta) \\ &= l_R(\alpha, \beta) + (\delta - 1) \sum_{d=1}^{D} \log(\beta_d) - \log(B(\delta)). \end{aligned} \qquad (6) $$

## 3 Monte Carlo variational inference (MCVI)

The joint optimization of $\alpha$ and $\beta$ in Eq. (4) is not trivial and no closed-form solution exists, not even when one fixes $\alpha$ (respectively $\beta$) and optimizes with respect to $\beta$ ($\alpha$) although in that case our problem reduces to a standard (kernelized) LR. Although the probabilistic framework introduced in Section 2.2 and leading to the log-likelihood in Eqs. (5) and (6) would allow us to adopt a fully Bayesian approach (no point estimates), we do adopt a hybrid perspective and infer the posterior density of $\beta$ and the MAP estimate of $\alpha$. This choice is driven by the desire to keep the inference as simple as possible while still benefiting from the flexibility and robustness offered by the Bayesian framework.

That said, the proposed probabilistic formulation relies onto variational inference in order to maximize a lower bound of the following integrated log-likelihood of the complete data

$$ \log p(y, \alpha|\mathcal{K}, \delta) = \log \int_{\beta} p(y, \alpha, \beta|\mathcal{K}, \delta)d\beta. $$
Our main contribution is to improve the variational lower bound (ELBO) of the above log-likelihood making use of the importance weighted lower bound (IW-ELBO), thus reinterpreting the approach introduced in Burda et al. (2015) in the context of Importance Weighted Autoencoders (IWAEs). Let us describe our strategy in more detail.

## 3.1 Variational lower bound

Focusing on the integrated log-likelihood with respect to $\beta$, we have that for any distribution $q(\cdot)$ on $\beta$, with the same support of the prior $p(\cdot|\delta)$, it holds that

$$ 
\begin{aligned}
\log p(y, \alpha | \mathcal{K}) &= \log \int_{\beta} p(y, \alpha, \beta | \mathcal{K}) d\beta \\
&= \log \int_{\beta} \frac{p(y, \alpha, \beta | \mathcal{K})}{q(\beta)} q(\beta) d\beta \qquad (7) \\
&= \log \mathbb{E}_{\beta \sim q} \left[ \frac{p(y, \alpha, \beta | \mathcal{K})}{q(\beta)} \right] \geq \mathbb{E}_{\beta \sim q} \left[ \log \frac{p(y, \alpha, \beta | \mathcal{K})}{q(\beta)} \right],
\end{aligned}
$$

where the inequality follows from Jensen's inequality. It is also well known that the above inequality turns into an equality when $q(\beta) = p(\beta|y, K, \alpha)$ which is however not analytically tractable here. This lower bound on the log-likelihood is commonly referred to as the evidence lower bound (ELBO) in the context of variational auto-encoding inference (Kingma et al. 2013). We will denote it by $\mathcal{L}(q, \alpha)$.

### Exponential ELBO.

With respect to the Exponential prior on $\beta$, leading to Eq. (5), one option is to set the variational (approximate) posterior distribution as

$$ q(\beta) = \prod_{d=1}^{D} q(\beta_d) $$
with $q(\beta_d) := q(\beta_d|\gamma_d)$ being the p.d.f. of an Exponential distribution of parameter $\gamma_d > 0$. By denoting $\gamma := (\gamma_1, \dots, \gamma_D)$, the lower-bound $\mathcal{L}(q, \alpha)$ can be developed as:

$$ 
\begin{aligned}
\mathcal{L}(q, \alpha) &= \mathbb{E}_{\beta \sim q_{\gamma}}[\log p(y, \alpha, \beta | \mathcal{K})] - \mathbb{E}_{\beta \sim q_{\gamma}}[\log q(\beta | \gamma)] = \\
&= \mathbb{E}_{\beta \sim q_{\gamma}}[l_R(\alpha, \beta)] + \mathbb{E}_{\beta \sim q_{\gamma}}\left[\log \frac{p(\beta | \delta)}{q(\beta | \gamma)}\right] = \\
&= \mathbb{E}_{\beta \sim q_{\gamma}}[l_R(\alpha, \beta)] - \sum_{d=1}^{D} \left(\frac{\delta}{\gamma_d} + \log \gamma_d\right) + D(\log \delta + 1)
\end{aligned}
\quad (8) $$

where the last two terms come from the explicit calculation of the negative Kullback-Leibler (K-L) divergence between the approximate posterior $q(\cdot|\gamma)$ and the prior distribution $p(\cdot|\delta)$. The expectation cannot be computed explicitly and this is where Monte Carlo (MC) comes into play. First, we approximate it as

$$ \mathbb{E}_{\beta \sim q_{\gamma}}[l_R(\alpha, \beta)] \approx \frac{1}{S} \sum_{s=1}^{S} l_R(\alpha, \beta^{(s)}), $$

where $S$ is the number of MC i.i.d. samples from the variational distribution $q_{\gamma}$. Then, for the Exponential distribution we can adopt the reparametrization trick, i.e. $\beta_d = \frac{1}{\gamma_d} \eta_d$, with $\eta_d \sim \text{Exponential}(1)$, i.i.d. for all $d$. Thus, following Kingma et al. (2013) we can reparameterize the above approximation by replacing $\beta^s$ with $(\frac{1}{\gamma_1} \eta_1^{(s)}, \dots, \frac{1}{\gamma_D} \eta_D^{(s)})$. This requires sampling $SD$ independent Exponential random variables of parameter one. After these steps, the ELBO is approximated by a quantity that is differentiable with respect to the model parameters $\alpha$ and the variational parameters $\gamma$. That approximation is an unbiased estimate of ELBO and the model and the variational parameters can be optimized by gradient descent (GD).
## Dirichlet ELBO.

With respect to the Dirichlet prior on $\beta$, leading to Eq. (6) and by denoting $\gamma := (\gamma_1, \dots, \gamma_D)$, the ELBO can be developed as:

$$ \mathcal{L}(q, \alpha) = \mathbb{E}_{\beta \sim q_{\gamma}} [l_R(\alpha, \beta)] + \sum_{d=1}^{D} (\delta - \gamma_d)(\psi(\gamma_d) - \psi(\gamma_0)) + \log \left( \frac{B(\gamma)}{B(\delta)} \right), \quad (9) $$

where the last two terms come from the explicit calculation of the negative K-L divergence between the variational posterior $q(\cdot|\gamma)$ and the prior $p(\cdot|\delta)$ in which $\psi(\cdot)$ denotes the digamma function and $\gamma_0 := \sum_{d=1}^{D} \gamma_d$.

The Monte Carlo approximation of the expectation is similar to the one described for the Exponential prior. However, unlike the previous scenario, reparametrizing samples from a Dirichlet distribution is more involved due to the lack of a straightforward differentiable transformation. In this case we leverage PyTorch's `.rsample()` function<fnref n="2" /> to draw differentiable samples, which relies on an implicit reparametrization of the Dirichlet distribution through Gamma-distributed variables (see Figurnov et al. 2018; Jankowiak and Obermeyer 2018, for more details). This allows gradients to propagate through the sampling process and enabling the optimization of the variational parameters. Explicit derivation of both the K-L divergences are reported in the supplementary materials.

## 3.2 Importance weighted lower bound

In this subsection, we extend the standard variational lower bound (ELBO) following the approach introduced in Burda et al. (2015). The main intuition is: instead of relying on a single random variable $\beta$ from the approximate posterior $q(\beta)$, we can draw multiple samples and average them. Indeed, the integrated log-likelihood can be

<footnote n="2">²https://docs.pytorch.org/docs/stable/distributions.html</footnote>
re-expressed as

$$
\begin{align}
\log p(y, \alpha | \mathcal{K}) &= \log \int_{\beta} p(y, \alpha, \beta | \mathcal{K}) d\beta = \log \int_{\beta} \frac{p(y, \alpha, \beta | \mathcal{K})}{q(\beta)} q(\beta) d\beta \nonumber \\
&= \log \mathbb{E}_{\beta \sim q_{\gamma}} \left[ \frac{p(y, \alpha, \beta | \mathcal{K})}{q(\beta)} \right] = \log \mathbb{E}_{\beta^{(1)}, \dots, \beta^{(R)} \stackrel{\text{i.i.d.}}{\sim} q_{\gamma}} \left[ \frac{1}{R} \sum_{j=1}^{R} \frac{p(y, \alpha, \beta^{(j)} | \mathcal{K})}{q(\beta^{(j)})} \right] \nonumber \\
&\geq \mathbb{E}_{\beta^{(1)}, \dots, \beta^{(R)} \stackrel{\text{i.i.d.}}{\sim} q_{\gamma}} \left[ \log \left( \frac{1}{R} \sum_{j=1}^{R} \frac{p(y, \alpha, \beta^{(j)} | \mathcal{K})}{q(\beta^{(j)})} \right) \right] =: \mathcal{L}_{(\alpha, \gamma)}^{\text{IW}}, \tag{10}
\end{align}
$$

where we used again the Jensen inequality<fnref n="3" />. As the number of samples *R* increases, the importance weighted lower bound (here denoted by $\mathcal{L}_{(\alpha, \gamma)}^{\text{IW}}$ for our particular case) monotonically increase up to convergence to the integrated log-likelihood it bounds. See Theorem 1 in Burda et al. (2015) for a formal statement of this result and relative assumptions. In particular, since $\mathcal{L}_{(\alpha, \gamma)}^{\text{IW}} = \mathcal{L}_{(q, \alpha)}$ for $R = 1$, it bounds the log-likelihood $\log p(y, \alpha | \mathcal{K})$ better than the ELBO for all $R \ge 2$.

For a fixed the sample $j \in \{1, \dots, R\}$, the term $p(y, \alpha, \beta^{(j)} | \mathcal{K})$ factorises as

$$p(y, \alpha, \beta^{(j)} | \mathcal{K}) = p(y | \mathcal{K}, \alpha, \beta^{(j)}) p(\alpha | \lambda, \mathcal{K}) p(\beta^{(j)} | \delta),$$

where the first and the second term on the r.h.s. of the equality are the exponential of Eqs. (2) and (3), respectively, and the third depends on the prior distribution we attach to $\beta$ (Exponential or Dirichlet). By making the second term explicit (not the first in order to keep the exposition readable), the IW-ELBO reads

$$
\mathcal{L}_{(\alpha, \gamma)}^{\text{IW}} = \mathbb{E}_{\beta^{(1)}, \dots, \beta^{(R)} \stackrel{\text{i.i.d.}}{\sim} q_{\gamma}} \left[ \log \left( \frac{1}{R} \sum_{j=1}^{R} \frac{p(y|\mathcal{K}, \alpha, \beta^{(j)}) \cdot p(\beta^{(j)}|\delta) \cdot |\lambda K|^{\frac{1}{2}}}{\exp(\frac{\lambda}{2}\alpha^T K \alpha) \cdot q(\beta^{(j)}|\gamma)} \right) \right] - \frac{N}{2} \log(2\pi), \quad (11)
$$

<footnote n="3">³The reader should not confuse the sampling scheme described in the above equation (R samples) with the one previously introduced to detail the MC approximation to approximated the ELBO (S samples). In particular, in the above Eq. (10) no MC approximation is adopted.</footnote>
where, in order to facilitate optimisation, non-sampling-dependent terms are taken out of the sum. As we did for the ELBO, here too $q_γ$ is assumed to be either an Exponential or a Dirichlet density function and the expectation is approximated via Monte Carlo and reparameterizations are used in order to update the variational parameters via GD.

### 3.3 Estimating $\beta$ and kernel(s) selection

Maximizing the lower bounds in Eqs. (10)-(11) provides us with estimates of $\alpha$ and the variational parameters $\gamma$. However, in order to compute the predictive accuracy of our model on a test data set and perform kernel selection, we also need an estimate of $\beta$. In the ELBO case, variational MAP estimates could be used by computing the mode of the variational distribution $q_γ$, where $\gamma$ denotes here the optimised variational parameter. However, the Exponential and Dirichlet distributions are unfit for such a strategy, since the former has mode zero (independently on $\gamma$) and the latter does not have a mode as soon as one $\gamma_d$ is smaller equal than one. As an alternative we chose to adopt the variational posterior mean in our experiments.

In the IW context, after training, the estimated kernel weights are sampled from an approximate posterior distribution as detailed in Algorithm 1. This distribution is referred to the *expected importance weighted* distribution and it was shown (Cremer et al. 2017) that i) the expected importance weighted distribution better approximates the posterior than $q_γ$ and ii) as soon as $R$ increases it converges toward the true posterior.

A simple way to obtain point estimates of the kernel weights $\beta$ is through the posterior mean (i.e. by averaging the samples obtained from several repetitions of Algorithm 1). This is computationally efficient and (as shown in later sections) leads to accurate test predictions. However, posterior means come with at least one important limitation: it prevents any effective shrinkage between kernels, assigning non-zero
<caption>**Algorithm 1** Sampling from the approximate posterior</caption>

1: $R \leftarrow$ number of importance samples

2: **for** $i \in (1, \dots, R)$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp; $\beta^{(i)} \sim q(\beta|\gamma)$

4: &nbsp;&nbsp;&nbsp;&nbsp; $w_i \leftarrow \frac{p(y, \alpha, \beta^{(i)}|K)}{q(\beta^{(i)}|\gamma)}$

5: **end for**

6: $\tilde{w}_i \leftarrow w_i / \sum_{i=1}^{R} w_i$

7: $j \sim \text{Categorical}(\tilde{\mathbf{w}})$

8: **Return** $\beta^* \leftarrow \beta^{(j)}$

weights to all kernel weights. This drawback can impair the ability to identify and discard irrelevant kernels, as we will further observe in our experiments.

To address this issue, we propose a straightforward and interpretable method for kernel selection based on the posterior histograms of the kernel weights. For each kernel $d \in \{1, \dots, D\}$, we analyze the mode of the histogram, i.e. the bin with the highest density. Let $P$ denote the number of bins used to compute the histogram. If the highest density for a kernel's posterior samples falls in the first bin (i.e., the left-most bin closest to zero), we consider that kernel irrelevant and set its weight to zero. For kernels whose highest density falls in a bin not being the left-most one we fix their weights to the left edge of the interval corresponding to the mode bin. This heuristic provides a natural way to discard redundant or uninformative kernels while preserving those with meaningful contributions. If the prior and approximate posterior are modeled using Exponential distributions, the proposed method can be used directly without modification. However, if the kernel weights are constrained to lie in the $(D-1)$-simplex, which ensures they sum to 1, as in the case of a Dirichlet distributions, an additional projection step is required (e.g. softmax transform)

The proposed method offers several advantages. First, it leverages the uncertainty inherent in the variational samples to assess the relevance of each kernel, providing a principled way to remove uninformative kernels. Additionally, the use of histograms makes the procedure intuitive, easy to implement, and interpretable. Moreover, the
hyper-parameter *P*, denoting the number of bins, controls the granularity of the selection process, allowing for adjustments based on the desired sensitivity. As a thumb rule, in the experiments we set *P* to be the minimum bin width between the Sturges (Sturges 1926) and the Freedman Diaconis (Freedman and Diaconis 1981) estimators. However, using more sophisticated approaches such as Kernel Density Estimation (KDE) could offer a more refined view of the posterior distributions and will be explored in future work.

# 4 Numerical examples

In this section we present experimental results that highlight the effectiveness of our MKL approach in a binary classification context.

## 4.1 Motivational example

We begin by presenting a toy example, in order to highlight how, for different distributional choices of priors for the kernel weights, our MKL approach performs kernel(s) selection in simple binary classification scenarios.

We consider two synthetic datasets, `moons` and `circles` (from the `sklearn`<fnref n="4" /> library), with *N* = 100 data points. We split the entire dataset into 60 points for training 40 for testing. We aim at solving the classification task by selecting kernels from a predefined set of kernel functions. Specifically, we test five kernels<fnref n="5" />: linear kernel, two polynomial kernels (of degrees 2 and 3), sigmoid kernel and Gaussian kernel.

We first consider an Exponential prior on the kernel weights and we optimize the variational lower bound, referring to it as Exp ELBO approach. We then move to the optimization of the importance-weighted lower bound, where two variants are presented. Hence, Exp IW-ELBO-f stands for the full model which uses the posterior

<footnote n="4">⁴https://scikit-learn.org/stable/</footnote>

<footnote n="5">⁵In this example, we adopted the default hyper-parameters of the sklearn library to compute kernels. In the case of the Gaussian kernel, therefore, $\sigma = 1/2$. A natural MKL application is to let the model learn the most appropriate hyper-parameters for a given kernel. In order to present a simple and intuitive example, we decided to skip this step here, while an application with several kernels can be found in Section 4.2.</footnote>
<table><caption>**Table 1:** Test accuracy and optimal kernel weights ($\beta_1^*, \dots, \beta_5^*$) associated to the linear, polynomials (degrees 2 and 3), sigmoid and Gaussian kernel respectively. Fixed hyper-parameters for the Exponential distributions: $\delta = 7$, $\lambda = 10^{-3}$. Fixed hyper-parameters for the Dirichlet distributions: $\delta = 1$, $\lambda = 10^{-3}$.</caption><thead><tr><th rowspan="2">Algorithm</th><th rowspan="2">Test Acc. (%)</th><th colspan="5">Kernel weights prediction via posterior mean</th></tr><tr><th>β<sub>1</sub><sup>*</sup></th><th>β<sub>2</sub><sup>*</sup></th><th>β<sub>3</sub><sup>*</sup></th><th>β<sub>4</sub><sup>*</sup></th><th>β<sub>5</sub><sup>*</sup></th></tr></thead><tbody><tr><td colspan="7" style="text-align:center;">moons</td></tr><tr><td>ℓ<sub>2</sub> GD</td><td>91.37 ± 4.14</td><td>0.2530 ± 0.1333</td><td>0.0359 ± 0.0438</td><td>0.5145 ± 0.0987</td><td>0.4635 ± 0.1326</td><td><b>0.8473 ± 0.0989</b></td></tr><tr><td>ℓ<sub>1</sub> GD</td><td>91.62 ± 3.73</td><td>0.0960 ± 0.0688</td><td>1.08e-09 ± 1e-10</td><td>5.01e-09 ± 1e-09</td><td>3.44e-09 ± 1e-09</td><td><b>1.1585 ± 0.1726</b></td></tr><tr><td>Exp ELBO</td><td>90.25 ± 3.25</td><td>0.1090 ± 0.0156</td><td>0.1531 ± 0.0081</td><td>0.2706 ± 0.0202</td><td>0.1232 ± 0.0149</td><td><b>0.5347 ± 0.1133</b></td></tr><tr><td>Exp IW-ELBO-f</td><td>91.12 ± 3.11</td><td>0.1250 ± 0.0138</td><td>0.1323 ± 0.0079</td><td>0.1998 ± 0.0113</td><td>0.1302 ± 0.0093</td><td><b>0.5434 ± 0.0908</b></td></tr><tr><td>Exp IW-ELBO-ks</td><td><b>92.37 ± 3.58</b></td><td></td><td></td><td colspan="3" style="text-align:center;">#Kernels 3.1 ± 0.64</td></tr><tr><td>Dir ELBO</td><td>91.25 ± 3.11</td><td>0.0978 ± 0.0073</td><td>0.1077 ± 0.0116</td><td>0.1624 ± 0.0133</td><td>0.1176 ± 0.0084</td><td><b>0.5145 ± 0.0342</b></td></tr><tr><td>Dir IW-ELBO-f</td><td>91.37 ± 3.01</td><td>0.1171 ± 0.0099</td><td>0.1241 ± 0.0104</td><td>0.1732 ± 0.0139</td><td>0.1258 ± 0.0085</td><td><b>0.4598 ± 0.0324</b></td></tr><tr><td>Dir IW-ELBO-ks</td><td>92.12 ± 3.66</td><td></td><td></td><td colspan="3" style="text-align:center;">#Kernels 3.6 ± 0.68</td></tr><tr><td colspan="7" style="text-align:center;">circles</td></tr><tr><td>ℓ<sub>2</sub> GD</td><td>87.00 ± 5.16</td><td>0.3187 ± 0.2229</td><td>0.3092 ± 0.0989</td><td>0.6220 ± 0.1923</td><td>0.1225 ± 0.0851</td><td><b>0.7547 ± 0.2460</b></td></tr><tr><td>ℓ<sub>1</sub> GD</td><td>87.12 ± 4.42</td><td>0.0897 ± 0.0834</td><td>3.70e-09 ± 1e-10</td><td>0.2216 ± 0.4691</td><td>1.55e-09 ± 1e-10</td><td><b>0.8333 ± 0.4558</b></td></tr><tr><td>Exp ELBO</td><td>87.63 ± 4.14</td><td>0.1339 ± 0.0063</td><td>0.2293 ± 0.0256</td><td>0.5223 ± 0.0615</td><td>0.0278 ± 0.0025</td><td><b>0.5380 ± 0.0298</b></td></tr><tr><td>Exp IW-ELBO-f</td><td>88.75 ± 4.22</td><td>0.1395 ± 0.0040</td><td>0.1972 ± 0.0143</td><td><b>0.4258 ± 0.0311</b></td><td>0.0524 ± 0.0055</td><td>0.3966 ± 0.0233</td></tr><tr><td>Exp IW-ELBO-ks</td><td>87.50 ± 4.47</td><td></td><td></td><td colspan="3" style="text-align:center;">#Kernels 2.4 ± 0.50</td></tr><tr><td>Dir ELBO</td><td>88.37 ± 4.13</td><td>0.0727 ± 0.0040</td><td>0.1143 ± 0.0062</td><td><b>0.5030 ± 0.0231</b></td><td>0.0354 ± 0.0025</td><td>0.2746 ± 0.0149</td></tr><tr><td>Dir IW-ELBO-f</td><td><b>90.00 ± 4.21</b></td><td>0.1154 ± 0.0045</td><td>0.1791 ± 0.0118</td><td><b>0.3447 ± 0.0147</b></td><td>0.0493 ± 0.0056</td><td>0.3115 ± 0.0051</td></tr><tr><td>Dir IW-ELBO-ks</td><td>87.87 ± 4.07</td><td></td><td></td><td colspan="3" style="text-align:center;">#Kernels 3.5 ± 0.51</td></tr></tbody></table>

mean of the kernel weights as a point estimate of $\beta$ and Exp IW-ELBO-ks stands for a sparse model. The sparsity is handled via the hyperparameter $P$, the number of bins in the histogram of the approximate posterior distribution of $\beta$, as presented in Section 3.3. We then constrain the kernel weights to lie within the simplex by considering Dirichlet distributions, leading to Dir ELBO, Dir IW-ELBO-f and Dir IW-ELBO-ks.

Additionally, we include comparisons with two gradient-based methods (GD) which optimize in $(\alpha, b)$ the following objective

$$f(\alpha, \beta) := \sum_{i=1}^{N} \left[ y_i \beta^T L_i \alpha - \log \left( 1 + e^{\beta^T L_i \alpha} \right) \right] - \frac{\lambda}{2} \alpha^T K \alpha - \delta \|\beta\|_p,$$

where $p \in \{1, 2\}$ and we set $\beta_j := \exp(b_j)$, for all $j$ in order to enforce the constraints $\beta_j > 0$. Hence $b = (b_1, \dots, b_5)$ in an unknown real vector and maximization of the objective with respect to it is performed via standard gradient descent. Sparsity promoting is achieved by setting $p = 1$ ($\ell_1$ GD) while a non-sparse solution is achieved by a Ridge penalization, e.g. set $p = 2$ ($\ell_2$ GD). Finally, in order to account for the data variability, 10 independent train/test split were realized.
<figure><img src="image_3.png" /><figcaption>Fig. 1: Frequency of selection of the five tested kernels, up to the 20 random splits.</figcaption></figure>

Table 1 reports the average test accuracy and the point estimates of the kernel weights for all methods. Although accuracy is not the primary objective of this example, all the proposed approaches achieve solid classification performance overall. Looking at the results on the *moons* dataset, we observe that effective kernel selection generally leads to better model generalization. However, for the *circles* dataset, the best accuracy is obtained by the IW-ELBO-f models, indicating that the underlying structure of the data is especially well captured by certain combinations of kernels. This highlights how specific data geometries may favor individual kernels, while in other cases, leveraging combinations, including apparently less influential kernels, can prove advantageous.

In both scenarios, we observe that less informative kernels, such as the linear one, still receive non-zero weight estimates. As discussed in Section 3.3, this highlights a limitation of using the posterior mean as estimator, which tends to vanish sparsity. However, in the IW-ELBO settings, a noticeable shrinkage effect emerges and can be better appreciated by inspecting the estimated posterior densities shown in Figures 2 and 3. Moreover, IW-ELBO-ks allow for an explicit inspection of kernel selection through histograms that represent how often each kernel was selected across different
runs (Figure 1). This provides a more interpretable and direct understanding of which
kernels are useful versus those that are redundant or irrelevant.

<figure><img src="image_4.png" /><figcaption>Fig. 2: Kernel weights densities for a random split of the moons dataset.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3: Kernel weights densities for a random split of the circles dataset.</figcaption></figure>

4.2 Benchmark datasets

In this section we present a comparative evaluation of our approach against two state-of-the-art MKL methods: SimpleMKL (Rakotomamonjy et al. 2008) and BEMKL (Gönen 2012).
The comparison focuses on number of selected kernels and test accuracy. We employ three benchmark data sets from the UCI<sup>6</sup> repository: ionosphere, sonar and wpbc.

We adopted the very same setting as described in Gönen (2012) that we now recall. Each dataset undergoes 20 independent randomized train/test splits, where 70% of the data being allocated for training and 30% for test. The training data are normalized in order to have zero mean and unit variance, while the test data are normalized using the corresponding training statistics. For each data set, we consider 10 Gaussian kernels with weights {2<sup>-3</sup>, 2<sup>-2</sup>, ..., 2<sup>6</sup>} and 3 polynomial kernels, respectively of degree 1, 2, and 3. Kernels are applied both globally (across all features) and individually on each feature. Thus, denoting *d* the number of features, the total number of input kernel matrices is *D* = 13(1 + *d*). All kernel matrices are then normalized to have unit diagonal. The results reported in Table 2 demonstrate that IW-ELBO-f consistently achieves superior accuracy across the three datasets when compared to classical MKL approaches. This performance gain might be attributed to the tighter variational lower bound provided by importance-weighted inference (IW-ELBO), which, as shown also in Example 1, enables more expressive posterior approximations and improves generalization. For instance, on the sonar dataset, IW-ELBO-f variants (Exp IW-ELBO-f and Dir IW-ELBO-f) significantly outperform both SimpleMKL and BEMKL, highlighting their ability to capture complex, distributed signal structures that sparse solutions may overlook. However, the sparse counterpart IW-ELBO-ks, while highly effective at kernel selection, often reducing the model to fewer than 10% of the original kernels, tends to underperform when informative signals are spread across many kernels. This is evident, for example, in the sonar dataset, where Exp IW-ELBO-ks achieves an average test accuracy of only 72%. A similar property can be observed also in the ionosphere and wpbc datasets, where the performance gap between IW-ELBO-f and IW-ELBO-ks is again notable. These results support the idea that aggressively

<footnote n="6"><sup>6</sup>https://archive.ics.uci.edu</footnote>
sparse solutions can impair generalization by discarding complementary information,
underscoring the need to carefully balance model expressiveness and interpretability.
BEMKL remains a strong competitor in terms of computational efficiency and sim-
plicity, offering closed form variational updates with reasonable predictive accuracy.
However, it relies on conjugate priors that impose strong assumptions on the form
of the approximate posterior and the type of uncertainty that can be modeled. As a
result, it is limited to a narrow family of posterior distributions, which may fail to
capture the complexity or multimodality of the true posterior over kernel weights, and
thus restrict its flexibility in modeling nuanced dependencies between kernels. Overall,
while IW-ELBO-f emerges as the most accurate and expressive option, IW-ELBO-ks
remains an attractive choice when interpretability and feature selection are priori-
tized, as long as one remains cautious about its impact on generalization, particularly
in noisy or highly redundant kernel spaces.

5 Application to an archaeological dataset

The aim of this section is to assess the potential of our approach in a archaeozoological
use case. Being the observed data a collection of bones, the classification problem that
we attack is the taxonomic identification of two different taxa: the roe deer (Capreolus
capreolus) and the gazelle (Gazella sp.). In archaeozoology, a similar task can be far
from trivial due to the morphological similarity between the two species, with only a
few distinguishing anatomical features.

In more detail, the observed data consists into 60 3D point clouds, resulting
from the 3D scan of modern and archaeological bones (astragali) from the two taxa
(Figure 4a). Each bone corresponds to a specimen and the sixty specimens are part
of a larger dataset recently published in Vuillien et al. (2025). The dataset is divided
into two groups which correspond to the taxa: 30 roe deer and 30 gazelle. Part of the
modern bones come from the reference collections of the National Museum of Natural
<table><caption>**Table 2:** Average results on benchmark MKL datasets. *N* and *D* denote the number of training instances and input kernels, respectively. Fixed hyper-parameters for Exponential distributions: *δ* = 10, *λ* = 10⁻³. Fixed hyper-parameters for Dirichlet distributions: *δ* = 1, *λ* = 10⁻³.</caption><thead><tr><th>Algorithm</th><th>Time (s)</th><th>#Kernels</th><th>Test Acc. (%)</th></tr></thead><tbody><tr><td colspan="4"><strong>ionosphere</strong></td></tr><tr><td>SimpleMKL (Rakotomamonjy et al. 2008)</td><td>123.0 ± 46.0</td><td>23.60 ± 2.60</td><td>91.50 ± 2.50</td></tr><tr><td>BEMKL (Gönen 2012)</td><td>22.72 ± 0.05</td><td>41.90 ± 4.42</td><td>92.03 ± 1.72</td></tr><tr><td>Exp ELBO</td><td>64.58 ± 3.14</td><td>-</td><td>92.41 ± 1.87</td></tr><tr><td>Exp IW-ELBO-f</td><td>110.18 ± 5.63</td><td>-</td><td><strong>92.88 ± 1.75</strong></td></tr><tr><td>Exp IW-ELBO-ks</td><td>-</td><td>7.20 ± 2.29</td><td>87.83 ± 4.34</td></tr><tr><td>Dir ELBO</td><td>70.81 ± 0.97</td><td>-</td><td>91.79 ± 1.86</td></tr><tr><td>Dir IW-ELBO-f</td><td>123.29 ± 2.15</td><td>-</td><td>92.26 ± 1.92</td></tr><tr><td>Dir IW-ELBO-ks</td><td>-</td><td>10.05 ± 1.77</td><td>88.87 ± 2.73</td></tr><tr><td colspan="4"><strong>sonar</strong></td></tr><tr><td>SimpleMKL (Rakotomamonjy et al. 2008)</td><td>163.0 ± 93.0</td><td>36.70 ± 5.10</td><td>80.60 ± 5.10</td></tr><tr><td>BEMKL (Gönen 2012)</td><td>51.34 ± 0.19</td><td>15.30 ± 3.40</td><td>76.88 ± 3.95</td></tr><tr><td>Exp ELBO</td><td>46.31 ± 4.09</td><td>-</td><td>84.44 ± 3.68</td></tr><tr><td>Exp IW-ELBO-f</td><td>217.5 ± 12.45</td><td>-</td><td>84.61 ± 2.75</td></tr><tr><td>Exp IW-ELBO-ks</td><td>-</td><td>12.40 ± 2.76</td><td>71.75 ± 4.25</td></tr><tr><td>Dir ELBO</td><td>22.25 ± 1.19</td><td>-</td><td>82.22 ± 3.39</td></tr><tr><td>Dir IW-ELBO-f</td><td>197.8 ± 16.23</td><td>-</td><td><strong>85.08 ± 3.84</strong></td></tr><tr><td>Dir IW-ELBO-ks</td><td>-</td><td>64.60 ± 6.92</td><td>80.31 ± 4.87</td></tr><tr><td colspan="4"><strong>wpbc</strong></td></tr><tr><td>SimpleMKL (Rakotomamonjy et al. 2008)</td><td>20.60 ± 6.20</td><td>15.80 ± 2.40</td><td>76.70 ± 1.20</td></tr><tr><td>BEMKL (Gönen 2012)</td><td>12.72 ± 0.11</td><td>28.25 ± 2.75</td><td>75.25 ± 1.01</td></tr><tr><td>Exp ELBO</td><td>6.38 ± 4.57</td><td>-</td><td>77.62 ± 4.95</td></tr><tr><td>Exp IW-ELBO-f</td><td>127.0 ± 6.53</td><td>-</td><td>77.62 ± 5.96</td></tr><tr><td>Exp IW-ELBO-ks</td><td>-</td><td>6.3 ± 2.1</td><td>70.51 ± 4.75</td></tr><tr><td>Dir ELBO</td><td>7.91 ± 0.47</td><td>-</td><td>77.97 ± 3.22</td></tr><tr><td>Dir IW-ELBO-f</td><td>100.7 ± 10.9</td><td>-</td><td><strong>78.98 ± 5.98</strong></td></tr><tr><td>Dir IW-ELBO-ks</td><td>-</td><td>35.90 ± 5.52</td><td>74.23 ± 5.40</td></tr></tbody></table>

History (Mammalian and Birds collection), in Paris, and the osteological collections from Archéorient (UMR 5133) laboratory. The archaeological bones come from the archaeological site of “Tell Sheikh Hassan” located in the actual Syria and dated from the IVth millenium BC.

## 5.1 Context

Roe deer (*Capreolus capreolus*) and gazelle (*Gazella* sp.) are taxa frequently encountered in archaeological contexts dating to the Holocene. They are mostly found in the form of more or less well-preserved bones. These herbivorous ungulates possess a relatively similar morphology, which can often be difficult to distinguish depending on
the degree of preservation of the bones and the specific bone under study (astragalus). Taxonomic identification is typically performed using anatomical and biometric characteristics found in their modern counterparts, which are available in physical or digital reference collections. However, the observation of these anatomical criteria depends largely on the conservation of the bone and the presence of diagnostic criteria. Furthermore, the morphological comparison between modern and ancient specimens is not always possible due to the evolution and adaptation of animals over time. This poses a significant challenge in identifying these taxa in an archaeological context.

More recently, molecular techniques such as palaeogenetics and palaeoproteomics have provided new approaches for identifying these species (e.g. Le Meillour et al. 2020, 2023; Rodrigues Soares et al. 2022). Despite their success, these methods have some limitations. Poor preservation of DNA or proteins, alterations to bone surfaces and gaps in reference data can hinder their applicability. As a result, alternative approaches based on the digital acquisition of 3D data and their analysis are becoming more and more popular. They include geometric morphometric methods for other wild ungulates (e.g. Curran 2012; Pelletier et al. 2022) and, more recently, machine and deep learning approaches (Vuillien et al. 2025).

## 5.2 Data pre-processing and kernel construction

After that the point clouds have been centered and scaled in such a way to be enclosed in the unit ball of $\mathbb{R}^3$, centered at the origin, they still have a different pose, number of points and there is no matching between points. Thus, defining a notion of similarity between them is far from being trivial. In order to avoid additional pre-processing steps such as the registration of the entire data set or its down-sampling$^7$, we relied on features extraction via Topological Data Analysis (TDA, Chazal and Michel 2021), a branch of mathematics designed to extract meaningful topological patterns from data

<footnote>$^7$Our point clouds are in the order of one million points. As such, the idea of using deep architectures for feature extraction would require down-sampling and, potentially, information loss.</footnote>
sets. Although an in depth presentation of TDA is beyond the scope of this paper, we sketch the main ideas in a very informal way. Consider a single point cloud. The points forming the cloud are initially considered as balls of null radius. Then, we let the radius increase progressively and, as far as the balls intersect with each other a simplicial complex (i.e. the multidimensional generalization of a neighbor graph) forms and becomes more and more dense. Global and local geometric features such as connected components (referred as topological features of dimension 0), loops or cycles (dimension 1), and three-dimensional holes or voids (dimension 2) are collected from the simplicial complex via persistent homology. The topological features are summarized in a persistence diagram (PD) (Figure 4c), which represents the lifetime of each feature across multiple scales.

We use the similarity between two PDs as a proxy of the similarity between the corresponding 3D point clouds. In more detail, fix a topological dimension (0,1 or 2), take two PDs and remove from each one all the points not corresponding to the chosen dimension. We then compare the two PDs via the 2-Wasserstein distance (Villani 2009; Peyré et al. 2019) between them and convert it into a similarity score via

$$K_{ij} := \exp\left(-\frac{W_2^2(PD_i, PD_j)}{\eta}\right), \qquad (12)$$

where $\eta$ is a positive temperature hyper-parameter and $W_2(\cdot, \cdot)$ denotes the 2-Wasserstein distance. Thus the above quantity is the $(i, j)$-th entry of a kernel matrix$^8 K \in \mathbb{R}_{+}^{N \times N}$. We might compute such a matrix for any topological dimension, however connected components (feature of dimension 0) are of no interest since they mostly account for the empty interior of the 3D cloud which has not a counterpart in the real astragalus.

<footnote>$^8$Properly speaking K is not a kernel matrix since it is not positive definite. However some simple manipulations related with its spectral decomposition can be made in order to make it positive definite (De Plaen et al. 2020).</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4: 3D data analysis pipeline via topological data analysis (TDA).</figcaption></figure>

Implementation.

For TDA we have used the gudhi toolbox (Maria et al. 2014) and POT toolbox
(Flamary et al. 2021) for the OT Wasserstein computations.

5.3 Species identification

Following the TDA framework illustrated above, we computed two Wasserstein dis-
tance matrices, where the first is based on topological cycles (features of dimension
1) and the second is based on topological voids (dimension 2). As input kernel matri-
ces we tested these two distance matrices, each of them translated into its kernelized
version via Eq. (12), with hyper-parameters η = { σ/₁₀, σ/₅, σ, 5σ, 10σ } and where σ is
the standard deviation computed from the respective distance matrix. We thus have
D = 10 input kernel matrices.

The entire data set is randomly split with stratification into 80% for training (48 bones) and 20% for testing (12 bones). 20 independent splits are repeated in order to assess the uncertainty around an average predictive accuracy. In Table 3 are reported the results in terms of average test accuracy and kernel weights prediction.

From the table, we observe that the average test accuracy is approximately 92-93%,
which aligns with human-level accuracy. More interesting, however, are the results
related to kernel selection performed by the IW-ELBO-ks models. In Figure 5, we
first notice that the highest weight is assigned to the matrix based on 1-dimensional
<table><caption>**Table 3:** Binary classification results between roe deer and gazelle. Fixed hyper-parameters for the Exponential distributions: $\delta = 10$, $\lambda = 10^{-3}$. Fixed hyper-parameters for Dirichlet distributions: $\delta = 1$, $\lambda = 10^{-3}$.</caption><thead><tr><th>Algorithm</th><th>#Kernels</th><th>Test Acc. (%)</th></tr></thead><tbody><tr><td>Exp ELBO</td><td>-</td><td>91.25 ± 8.11</td></tr><tr><td>Exp IW-ELBO-f</td><td>-</td><td><b>93.75 ± 5.81</b></td></tr><tr><td>Exp IW-ELBO-ks</td><td>1.4 ± 0.75</td><td>91.25 ± 9.67</td></tr><tr><td>Dir ELBO</td><td>-</td><td>91.25 ± 8.11</td></tr><tr><td>Dir IW-ELBO-f</td><td>-</td><td>92.50 ± 7.40</td></tr><tr><td>Dir IW-ELBO-ks</td><td>2.5 ± 1.10</td><td>90.42 ± 8.02</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5:** Frequency of selection of the ten tested kernels, up to the 20 random splits. The first half of the kernels ($\beta_1, \dots, \beta_5$) is associated to the first topological dimension while the second half ($\beta_6, \dots, \beta_{10}$) represents the second dimension.</figcaption></figure>

cycles. Specifically, $\beta_3$ corresponds to the weight of the Gaussian kernel matrix where $\eta$ is equal to $\sigma$. This matrix, shown in Figure 6a, highlights two clusters corresponding to the roe deer and gazelle species. Upon closer inspection of the matrix, we can also observe some intraspecific differences, particularly linked to certain gazelle specimens. On the other hand, the kernel matrices based on topological voids (dimension
<figure><img src="image_8.png" /><figcaption>**Fig. 6:** Most informative kernel matrices corresponding to weights $\beta_3$ and $\beta_8$, respectively.</figcaption></figure>

2) seem to be less informative. This observation is confirmed by the matrix shown in Figure 6b, which does not reveal significant properties related to structural or morphological differences among specimens. Nevertheless, their combined use with the matrix in Figure 6a proves crucial. In all tested methods, the model consistently learns a kernel matrix that heavily relies on the 1-dimensional cycles matrix while also incorporating some voids-based matrices (dimension 2). This combination helps to densify the two clusters corresponding to roe deer and gazelle species and to provide better generalization (Figure 7).

# 6 Conclusion

We revisited a multiple kernel logistic regression problem through the lenses of modern Monte Carlo importance-weighted variational inference, proposing a hybrid framework that combines the simplicity of optimization based methods with the expressiveness and uncertainty modeling of Bayesian approaches. Our formulation relies on the use of the IW-ELBO, e.g. the importance-weighted lower bound, which allows for tighter variational approximations of the marginal log-likelihood and efficient posterior
<figure><img src="image_9.png" /><figcaption>**Fig. 7:** Estimated linear combination for a single random split of the data set. It has been computed within the Exp IW-ELBO-f approach and using the posterior mean as point estimate for $\beta$.</figcaption></figure>

sampling over the kernel weights. The two resulting variants IW-ELBO-f and IW-ELBO-ks, adopting different estimates of $\beta$, offer complementary strengths. Although both methods achieve state-of-the-art accuracy, the former consistently provides better accuracy by retaining the full richness of the posterior, while the latter provides interpretable and efficient kernel(s) selection via posterior, histogram-based, thresholding. An inspection on standard benchmark datasets confirms that IW-ELBO is especially efficient in scenarios where relevant information is distributed across multiple kernels. Furthermore, the application to 3D point clouds from an archaeological context illustrates the model's ability to learn meaningful combinations of kernels, here based on persistence diagrams, and to identify the most informative among them, such as those associated with topological features that better capture variation between
shapes. Despite these promising results, the current approach presents some limita-
tions. As already mentioned, the use of posterior means as kernel weights estimates
in IW-ELBO-f tends to dilute sparsity, making it less suitable when explicit feature
selection is desired. On the other hand, the heuristic shrinking strategy adopted in
IW-ELBO-ks could benefit from more principled alternatives, such as Kernel Density
Estimation (KDE) or automatic kernel hyperparameter learning.

Funding

This paper constitutes the doctoral thesis in applied mathematics of the first
author. This work has received financial support from the CNRS through the MITI
interdisciplinary programs.

Conflict of Interest

The authors have no competing interests to declare.

Author Contribution

All authors participated in writing, revising, and approving the final manuscript.
The main text, including the theoretical and experimental sections and the method
implementations were done by D.A. and M.C. The archaeological context and
zooarchaeological data preparation were done by M.V. and E.V.

Data Availability

The data used in this study consist of publicly available benchmark datasets (Ionosphere, Sonar, WPBC from the UCI repository), which can be found at https://archive.ics.uci.edu/. The archaeological dataset comprises of 3D scans of modern and archaeological bones, as described in Vuillien et al. (2025). The 3D models of the
modern gazelle and roe deer astragali are available on the 3D platform of the National Museum of Natural History of Paris, which can be found in the “2024_Astragales” file (https://3dtheque.mnhn.fr/work/1039).

# Appendix A K-L divergence derivations

## Exponential distributions

Proof

$$
\begin{align*}
\mathbb{E}_q \left[ \log \frac{q(\beta|\gamma)}{p(\beta|\delta)} \right] &= \mathbb{E}_q \left[ \log \frac{\prod_{d=1}^D \gamma_d e^{-\gamma_d \beta_d}}{\prod_{d=1}^D \delta e^{-\delta \beta_d}} \right] = \mathbb{E}_q \left[ \log \left( \prod_{d=1}^D \frac{\gamma_d}{\delta} e^{\beta_d (\delta - \gamma_d)} \right) \right] \\
&= \mathbb{E}_q \left[ \sum_{d=1}^D \left( \log(\gamma_d) - \log(\delta) + \beta_d (\delta - \gamma_d) \right) \right] \\
&= \sum_{d=1}^D \log(\gamma_d) - D \log(\delta) + \sum_{d=1}^D \mathbb{E}_q[\beta_d](\delta - \gamma_d) \\
&= \sum_{d=1}^D \log(\gamma_d) - D \log(\delta) + \sum_{d=1}^D \frac{1}{\gamma_d} (\delta - \gamma_d) \\
&= \sum_{d=1}^D \log(\gamma_d) - D \log(\delta) + \sum_{d=1}^D \frac{\delta}{\gamma_d} - D \\
&= \sum_{d=1}^D \left( \frac{\delta}{\gamma_d} + \log(\gamma_d) \right) - D(\log(\delta) + 1)
\end{align*}
$$

□
Dirichlet distributions

Proof

$$
\begin{align*}
\mathbb{E}_q \left[ \log \frac{q(\beta|\gamma)}{p(\beta|\delta)} \right] &= \mathbb{E}_q \left[ \log \frac{\frac{1}{B(\gamma)} \prod_{d=1}^{D} (\beta_d)^{\gamma_d-1}}{\frac{1}{B(\delta)} \prod_{d=1}^{D} (\beta_d)^{\delta-1}} \right] = \mathbb{E}_q \left[ \log \left( \frac{B(\delta)}{B(\gamma)} \prod_{d=1}^{D} (\beta_d)^{\gamma_d-\delta} \right) \right] \\
&= \mathbb{E}_q \left[ \log(B(\delta)) - \log(B(\gamma)) + \sum_{d=1}^{D} (\gamma_d - \delta) \log(\beta_d) \right] \\
&= \log(B(\delta)) - \log(B(\gamma)) + \sum_{d=1}^{D} (\gamma_d - \delta) \mathbb{E}_q [\log(\beta_d)] \\
&= \log(B(\delta)) - \log(B(\gamma)) + \sum_{d=1}^{D} (\gamma_d - \delta) (\psi(\gamma_d) - \psi(\gamma_0)).
\end{align*}
$$

□

References

Bach, F.R.: Consistency of the group lasso and multiple kernel learning. Journal of Machine Learning Research 9(6) (2008)

Burda, Y., Grosse, R., Salakhutdinov, R.: Importance weighted autoencoders. arXiv preprint arXiv:1509.00519 (2015)

Ben-Hur, A., Noble, W.S.: Kernel methods for predicting protein-protein interactions. Bioinformatics 21(suppl.1), 38–46 (2005)

Bishop, C.M.: Pattern Recognition and Machine Learning vol. 4. Springer, ??? (1991)

Bucak, S.S., Jin, R., Jain, A.K.: Multiple kernel learning for visual object recognition: A review. IEEE Transactions on Pattern Analysis and Machine Intelligence 36(7), 1354–1369 (2013)

Blei, D.M., Kucukelbir, A., McAuliffe, J.D.: Variational inference: A review for statisticians. Journal of the American statistical Association 112(518), 859–877 (2017)
Chazal, F., Michel, B.: An introduction to topological data analysis: fundamental and practical aspects for data scientists. *Frontiers in artificial intelligence* **4**, 667963 (2021)

Cremer, C., Morris, Q., Duvenaud, D.: Reinterpreting importance-weighted autoencoders. *arXiv preprint arXiv:1704.02916* (2017)

Curran, S.C.: Expanding ecomorphological methods: geometric morphometric analysis of cervidae post-crania. *Journal of Archaeological Science* **39**(4), 1172–1182 (2012)

Damoulas, T., Girolami, M.A.: Probabilistic multi-class multi-kernel learning: on protein fold recognition and remote homology detection. *Bioinformatics* **24**(10), 1264–1270 (2008)

De Plaen, H., Fanuel, M., Suykens, J.A.: Wasserstein exponential kernels. In: 2020 International Joint Conference on Neural Networks (IJCNN), pp. 1–6 (2020). IEEE

Flamary, R., Courty, N., Gramfort, A., Alaya, M.Z., Boisbunon, A., Chambon, S., Chapel, L., Corenflos, A., Fatras, K., Fournier, N., Gautheron, L., Gayraud, N.T.H., Janati, H., Rakotomamonjy, A., Redko, I., Rolet, A., Schutz, A., Seguy, V., Sutherland, D.J., Tavenard, R., Tong, A., Vayer, T.: Pot: Python optimal transport. *Journal of Machine Learning Research* **22**(78), 1–8 (2021)

Freedman, D., Diaconis, P.: On the histogram as a density estimator: L 2 theory. *Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete* **57**(4), 453–476 (1981)

Figurnov, M., Mohamed, S., Mnih, A.: Implicit reparameterization gradients. *Advances in neural information processing systems* **31** (2018)

Gönen, M., Alpaydın, E.: Multiple kernel learning algorithms. *The Journal of Machine*
Gehler, P., Nowozin, S.: On feature combination for multiclass object classification. In: 2009 IEEE 12th International Conference on Computer Vision, pp. 221–228 (2009). IEEE

Gönen, M.: Bayesian efficient multiple kernel learning. In: Proceedings of the 29th International Conference on International Conference on Machine Learning, pp. 91–98 (2012)

Girolami, M., Rogers, S.: Hierarchicsbayesian models for kernel learning. In: Proceedings of the 22nd International Conference on Machine Learning, pp. 241–248 (2005)

He, X., Huang, J., Zeng, Z.: Logistic regression based multi-task, multi-kernel learning for emotion recognition. In: 2021 6th IEEE International Conference on Advanced Robotics and Mechatronics (ICARM), pp. 572–577 (2021). IEEE

Huang, F.J., LeCun, Y.: Large-scale learning with svm and convolutional for generic object categorization. In: 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06), vol. 1, pp. 284–291 (2006). IEEE

Jankowiak, M., Obermeyer, F.: Pathwise derivatives beyond the reparameterization trick. In: International Conference on Machine Learning, pp. 2235–2244 (2018). PMLR

Kloft, M., Brefeld, U., Sonnenburg, S., Zien, A.: Lp-norm multiple kernel learning. The Journal of Machine Learning Research **12**, 953–997 (2011)

Kingma, D.P., Welling, M., et al.: Auto-encoding variational bayes. Banff, Canada (2013)
Le Meillour, L., Zirah, S., Zazzo, A., Cersoy, S., Détroit, F., Imalwa, E., Lebon, M., Nankela, A., Tombret, O., Pleurdeau, D., et al.: Palaeoproteomics gives new insight into early southern african pastoralism. Scientific Reports **10**(1), 14427 (2020)

Le Meillour, L., Zazzo, A., Zirah, S., Tombret, O., Barriel, V., Arthur, K.W., Arthur, J.W., Cauliez, J., Chaix, L., Curtis, M.C., et al.: The name of the game: palaeoproteomics and radiocarbon dates further refine the presence and dispersal of caprines in eastern and southern Africa. Royal Society Open Science **10**(11), 231002 (2023)

Mairal, J.: End-to-end kernel learning with supervised convolutional kernel networks. Advances in neural information processing systems **29** (2016)

Maria, C., Boissonnat, J.-D., Glisse, M., Yvinec, M.: The gudhi library: Simplicial complexes and persistent homology. In: Mathematical Software–ICMS 2014: 4th International Congress, Seoul, South Korea, August 5-9, 2014. Proceedings 4, pp. 167–174 (2014). Springer

Mairal, J., Koniusz, P., Harchaoui, Z., Schmid, C.: Convolutional kernel networks. Advances in neural information processing systems **27** (2014)

Mohammadnia-Qaraei, M.R., Monsefi, R., Ghiasi-Shirazi, K.: Convolutional kernel networks based on a convex combination of cosine kernels. Pattern Recognition Letters **116**, 127–134 (2018)

Peyré, G., Cuturi, M., et al.: Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning **11**(5-6), 355–607 (2019)

Pelletier, M., Kotiaho, A., Niinimäki, S., Salmi, A.-K.: Impact of selection and domestication on hindlimb bones of modern reindeer populations: Archaeological implications for early reindeer management by sámi in fennoscandia. Historical
Biology **34**(5), 802-820 (2022)

Rakotomamonjy, A., Bach, F., Canu, S., Grandvalet, Y.: Simplemkl. Journal of Machine Learning Research **9**, 2491–2521 (2008)

Rodrigues Soares, A.E., Boroffka, N., Schröder, O., Sverchkov, L., Benecke, N., Günther, T.: Ancient dna from a 2700-year-old goitered gazelle (gazella subgutturosa) supports gazelle hunting in iron age central asia. Royal Society Open Science **9**(6), 220104 (2022)

Schölkopf, B., Herbrich, R., Smola, A.J.: A generalized representer theorem. In: International Conference on Computational Learning Theory, pp. 416–426 (2001). Springer

Stone, M.: Cross-validatory choice and assessment of statistical predictions. Journal of the royal statistical society: Series B (Methodological) **36**(2), 111–133 (1974)

Sturges, H.A.: The choice of a class interval. Journal of the american statistical association **21**(153), 65–66 (1926)

Tang, Y.: Deep learning using linear support vector machines. arXiv preprint arXiv:1306.0239 (2013)

Vapnik, V., et al.: Statistical learning theory. John Wiley & Sons, New York (1998)

Vuillien, M., Adamo, D., Vila, E., Agraw, A., Helmer, D., Mashkour, M., Moussous, A., Notter, O., Rossoni-Notter, E., Théry, A., Corneli, M.: Topological data analysis and multiple kernel learning for species identification of modern and archaeological small ruminants. Journal of Computer Applications in Archaeology **8**(1), 170–187 (2025)

Villani, C.: The wasserstein distances. Optimal transport: old and new, 93–111 (2009)
Wilson, A.G., Hu, Z., Salakhutdinov, R.R., Xing, E.P.: Stochastic variational deep kernel learning. Advances in neural information processing systems 29 (2016)

Wilson, A.G., Hu, Z., Salakhutdinov, R., Xing, E.P.: Deep kernel learning. In: Artificial Intelligence and Statistics, pp. 370–378 (2016). PMLR

Wang, T., Zhang, L., Hu, W.: Bridging deep and multiple kernel learning: A review. Information Fusion 67, 3–13 (2021)

Zien, A., Ong, C.S.: Multiclass multiple kernel learning. In: Proceedings of the 24th International Conference on Machine Learning, pp. 1191–1198 (2007)