Computing the eigenvectors and eigenvalues of a perturbed matrix can be remarkably difficult when the unperturbed matrix has repeated eigenvalues. In this work, we show how the limiting eigenvectors and eigenvalues of a symmetric matrix K(epsilon) as epsilon -> 0 can be obtained relatively easily from successive Schur complements provided that the entries scale in different orders of epsilon. If the matrix does not directly exhibit this structure, then we show that putting the matrix into a "generalized kernel form" can be very informative. The resulting formulas are much simpler than classical expressions obtained from complex integrals involving the resolvent. We apply our results to the problem of computing the eigenvalues and eigenvectors of kernel matrices in the "flat limit," a problem that appears in many applications in statistics and approximation theory. In particular, we prove a conjecture from Barthelme'\ and Usevich [SIAM J. Matrix Anal. Appl., 42 (2021), pp. 17--57] which connects the eigenvectors of kernel matrices to multivariate orthogonal polynomials.

Variational inference seeks the best approximation of a target distribution within a chosen family, where "best" means minimizing Kullback-Leibler divergence. When the approximation family is exponential, the optimal approximation satisfies a fixed-point equation. We introduce LSVI (Least Squares Variational Inference), a gradient-free, Monte Carlo-based scheme for the fixed-point recursion, where each iteration boils down to performing ordinary least squares regression on tempered log-target evaluations under the variational approximation. We show that LSVI is equivalent to biased stochastic natural gradient descent and use this to derive convergence rates with respect to the numbers of samples and iterations. When the approximation family is Gaussian, LSVI involves inverting the Fisher information matrix, whose size grows quadratically with dimension $d$. We exploit the regression formulation to eliminate the need for this inversion, yielding $O(d^3)$ complexity in the full-covariance case and $O(d)$ in the mean-field case. Finally, we numerically demonstrate LSVI’s performance on various tasks, including logistic regression, discrete variable selection, and Bayesian synthetic likelihood, showing competitive results with state-of-the-art methods, even when gradients are unavailable.

Exact eigendecomposition of large matrices is very expensive, and it is practically impossible to compute exact eigenvalues. Instead, one may set a more modest goal of approaching the empirical distribution of the eigenvalues, recovering the overall shape of the eigenspectrum. Current approaches to spectral estimation typically work with moments of the spectral distribution. These moments are first estimated using Monte Carlo trace estimators, then the estimates are combined to approximate the spectral density. In this article we show how Kirchhoff forests, which are random forests on graphs, can be used to estimate certain non-linear moments of very large graph Laplacians. We show how to combine these moments into an estimate of the spectral density. If the estimate's desired precision isn't too high, our approach paves the way to the estimation of a graph's spectrum in time sublinear in the number of links.

Given a non-oriented edge-weighted graph, we show how to make some estimation of the associated Laplacian eigenvalues through Monte Carlo evaluation of spectral quantities computed along Kirchhoff random rooted spanning forest trajectories. The sampling cost of this estimation is only linear in the node number, up to a logarithmic factor. By associating a double cover of such a graph with any symmetric real matrix, we can then perform spectral estimation in the same way for the latter.

Variational inference consists in finding the best approximation of a target distribution within a certain family, where `best' means (typically) smallest Kullback-Leiber divergence. We show that, when the approximation family is exponential, the best approximation is the solution of a fixed-point equation. We introduce LSVI (Least-Squares Variational Inference), a Monte Carlo variant of the corresponding fixed-point recursion, where each iteration boils down to ordinary least squares regression and does not require computing gradients. We show that LSVI is equivalent to stochastic mirror descent; we use this insight to derive convergence guarantees. We introduce various ideas to improve LSVI further when the approximation family is Gaussian, leading to a $O(d^3)$ complexity in the dimension $d$ of the target in the full-covariance case, and a $O(d)$ complexity in the mean-field case. We show that LSVI outperforms state-of-the-art methods in a range of examples, while remaining gradient-free, that is, it does not require computing gradients.

Random diffusions are a popular tool in Monte Carlo estimations, with well-established algorithms such as walk-on-spheres (WoS) going back several decades. In this work, we introduce diffusion estimators for the problems of angular synchronization and smoothing on graphs, in the presence of a rotation associated to each edge. Unlike classical WoS algorithms that are pointwise estimators, our diffusion estimators allow for global estimations by propagating along the branches of random spanning subgraphs called multitype spanning forests. Building upon efficient samplers based on variants of Wilson's algorithm, we show that our estimators outperform standard numerical-linear-algebra solvers in challenging instances, depending on the topology and density of the graph.

Node regression consists in predicting the value of a graph label at a node, given observations at the other nodes. To gain some insight into the performance of various estimators for this task, we perform a theoretical study in a context where the graph is random. Specifically, we assume that the graph is generated by a Latent Position Model, where each node of the graph has a latent position, and the probability that two nodes are connected depend on the distance between the latent positions of the two nodes. In this context, we begin by studying the simplest possible estimator for graph regression, which consists in averaging the value of the label at all neighboring nodes. We show that in Latent Position Models this estimator tends to a Nadaraya Watson estimator in the latent space, and that its rate of convergence is in fact the same. One issue with this standard estimator is that it averages over a region consisting of all neighbors of a node, and that depending on the graph model this may be too much or too little. An alternative consists in first estimating the true distances between the latent positions, then injecting these estimated distances into a classical Nadaraya Watson estimator. This enables averaging in regions either smaller or larger than the typical graph neighborhood. We show that this method can achieve standard nonparametric rates in certain instances even when the graph neighborhood is too large or too small.

Determinantal point processes (DPPs) are a class of repulsive point processes, popular for their relative simplicity. They are traditionally defined via their marginal distributions, but a subset of DPPs called "L-ensembles" have tractable likelihoods and are thus particularly easy to work with. Indeed, in many applications, DPPs are more naturally defined based on the L-ensemble formulation rather than through the marginal kernel. The fact that not all DPPs are L-ensembles is unfortunate, but there is a unifying description. We introduce here extended L-ensembles, and show that all DPPs are extended L-ensembles (and vice-versa). Extended L-ensembles have very simple likelihood functions, contain L-ensembles and projection DPPs as special cases. From a theoretical standpoint, they fix some pathologies in the usual formalism of DPPs, for instance the fact that projection DPPs are not L-ensembles. From a practical standpoint, they extend the set of kernel functions that may be used to define DPPs: we show that conditional positive definite kernels are good candidates for defining DPPs, including DPPs that need no spatial scale parameter. Finally, extended L-ensembles are based on so-called ``saddle-point matrices'', and we prove an extension of the Cauchy-Binet theorem for such matrices that may be of independent interest.

We introduce new smoothing estimators for complex signals on graphs, based on a recently studied Determinantal Point Process (DPP). These estimators are built from subsets of edges and nodes drawn according to this DPP, making up trees and unicycles, i.e., connected components containing exactly one cycle. We provide a Julia implementation of these estimators and study their performance when applied to a ranking problem.

Determinantal point processes (DPPs) are repulsive point processes where the interaction between points depends on the determinant of a positive-semi definite matrix. In this paper, we study the limiting process of L-ensembles based on kernel matrices, when the kernel function becomes flat (so that every point interacts with every other point, in a sense). We show that these limiting processes are best described in the formalism of extended L-ensembles and partial projection DPPs, and the exact limit depends mostly on the smoothness of the kernel function. In some cases, the limiting process is even universal, meaning that it does not depend on specifics of the kernel function, but only on its degree of smoothness. Since flat-limit DPPs are still repulsive processes, this implies that practically useful families of DPPs exist that do not require a spatial length-scale parameter.

Discrete Determinantal Point Processes (DPPs) have a wide array of potential applications for subsampling datasets. They are however held back in some cases by the high cost of sampling. In the worst-case scenario, the sampling cost scales as O(n^3) where n is the number of elements of the ground set. A popular workaround to this prohibitive cost is to sample DPPs defined by low-rank kernels. In such cases, the cost of standard sampling algorithms scales as O(np^2 + nm^2) where m is the (average) number of samples of the DPP (usually m << n) and p the rank of the kernel used to define the DPP (m \leq p \leq n). The first term, O(np^2), comes from a SVD-like step. We focus here on the second term of this cost, O(nm^2), and show that it can be brought down to O(nm + m^3 log m) without loss on the sampling's exactness. In practice, we observe very substantial speedups compared to the classical algorithm as soon as n > 1000. The algorithm described here is a close variant of the standard algorithm for sampling continuous DPPs, and uses rejection sampling. In the specific case of projection DPPs, we also show that any additional sample can be drawn in time O(m^3 log m). Finally, an interesting by-product of the analysis is that a realisation from a DPP is typically contained in a subset of size O(m log m) formed using leverage score i.i.d. sampling.

Gaussian process (GP) regression is a fundamental tool in Bayesian statistics. It is also known as kriging and is the Bayesian counterpart to the frequentist kernel ridge regression. Most of the theoretical work on GP regression has focused on a large-$n$ asymptotics, i.e. as the amount of data increases. Fixed-sample analysis is much more difficult outside of simple cases, such as locations on a regular grid. In this work we perform a fixed-sample analysis that was first studied in the context of approximation theory by Driscoll & Fornberg (2002), called the ``flat limit''. In flat-limit asymptotics, the goal is to characterise kernel methods as the length-scale of the kernel function tends to infinity, so that kernels appear flat over the range of the data. Surprisingly, this limit is well-defined, and displays interesting behaviour: Driscoll & Fornberg showed that radial basis interpolation converges in the flat limit to polynomial interpolation, if the kernel is Gaussian. Subsequent work showed that this holds true in the multivariate setting as well, but that kernels other than the Gaussian may have (polyharmonic) splines as the limit interpolant. Leveraging recent results on the spectral behaviour of kernel matrices in the flat limit, we study the flat limit of Gaussian process regression. Results show that Gaussian process regression tends in the flat limit to (multivariate) polynomial regression, or (polyharmonic) spline regression, depending on the kernel. Importantly, this holds for both the predictive mean and the predictive variance, so that the posterior predictive distributions become equivalent. Our results have practical consequences: for instance, they show that optimal GP predictions in the sense of leave-one-out loss may occur at very large length-scales, which would be invisible to current implementations because of numerical difficulties.

Large dimensional least-squares and regularized least-squares problems are expensive to solve. There exist many approximate techniques, some deterministic (like conjugate gradient), some stochastic (like stochastic gradient descent). Among the latter, a new class of techniques uses Determinantal Point Processes (DPPs) to produce unbiased estimators of the solution. In particular, they can be used to perform Tikhonov regularization on graphs using random spanning forests, a specific DPP. While the unbiasedness of these algorithms is attractive, their variance can be high. We show here that variance can be reduced by combining the stochastic estimator with a deterministic gradient-descent step, while keeping the property of unbiasedness. We apply this technique to Tikhonov regularization on graphs, where the reduction in variance is found to be substantial at very small extra cost.

The trace $\tr(q(\ma{L} + q\ma{I})^{-1})$, where $\ma{L}$ is a symmetric diagonally dominant matrix, is the quantity of interest in some machine learning problems. However, its direct computation is impractical if the matrix size is large. State-of-the-art methods include Hutchinson's estimator combined with iterative solvers, as well as the estimator based on random spanning forests (a random process on graphs). In this work, we show two ways of improving the forest-based estimator via well-known variance reduction techniques, namely control variates and stratified sampling. Implementing these techniques is easy, and provides substantial variance reduction, yielding comparable or better performance relative to state-of-the-art algorithms.

Gas sensors lack repeatability over time. They are affected by drift, the result of changes at the sensor level and in the environment. A solution is to design software methods that compensate for the drift. Existing methods are often based on calibration samples acquired at the start of each new measurement session. However, finding a good reference compound is a difficult task and generating calibration samples is time-consuming. We propose a model-based correction method which does not require any calibration sample over time, operating 'blindly'. In this study, we focus on the drift affecting electronic noses. To this end, we built a real data set acquired over 9 months in real-life conditions. By using the proposed method, we show that the drift is partly compensated, thus increasing the reliability of the electronic nose. Besides, we also show that the algorithm can easily adapt if the target compounds are not all sampled during every session.

Kernel matrices are of central importance to many applied fields. In this manuscript, we focus on spectral properties of kernel matrices in the so-called “flat limit,” which occurs when points are close together relative to the scale of the kernel. We establish asymptotic expressions for the determinants of the kernel matrices, which we then leverage to obtain asymptotic expressions for the main terms of the eigenvalues. Analyticity of the eigenprojectors yields expressions for limiting eigenvectors, which are strongly tied to discrete orthogonal polynomials. Both smooth and finitely smooth kernels are covered, with stronger results available in the finite smoothness case.