Abstract. We propose a new notion of uniqueness for the adversarial Bayes classifier in the setting of binary classification. Analyzing this concept produces a simple procedure for computing all adversarial Bayes classifiers for a well-motivated family of one-dimensional data distributions. This characterization is then leveraged to show that as the perturbation radius increases, certain notions of regularity of adversarial Bayes classifiers improve. These results provide tools for understanding relationships between the Bayes and adversarial Bayes classifiers in one dimension.
We consider scattered data approximation on product regions of equal and different dimensionality. On each of these regions, we assume quasi-uniform but unstructured data sites and construct optimal sparse grids for scattered data interpolation on the product region. For this, we derive new improved error estimates for the respective kernel interpolation error by invoking duality arguments. An efficient algorithm to solve the underlying linear system of equations is proposed. The algorithm is based on the sparse grid combination technique, where a sparse direct solver is used for the elementary anisotropic tensor product kernel interpolation problems. The application of the sparse direct solver is facilitated by applying a samplet matrix compression to each univariate kernel matrix, resulting in an essentially sparse representation of the latter. In this way, we obtain a method that is able to deal with large problems up to billions of interpolation points, especially in case of reproducing kernels of nonlocal nature. Numerical results are presented to qualify and quantify the approach.
This paper considers the problem of estimating a matrix that encodes pairwise distances in a finite metric space (or, more generally, the edge weight matrix of a network) under the barycentric coding model (BCM) with respect to the Gromov-Wasserstein (GW) distance function. We frame this task as estimating the unknown barycentric coordinates with respect to the GW distance, assuming that the target matrix (or kernel) belongs to the set of GW barycenters of a finite collection of known templates. In the language of harmonic analysis, if computing GW barycenters can be viewed as a synthesis problem, this paper aims to solve the corresponding analysis problem. We propose two methods: one utilizing fixed-point iteration for computing GW barycenters, and another employing a differentiation-based approach to the GW structure using a blow-up technique. Finally, we demonstrate the application of the proposed GW analysis approach in a series of numerical experiments and applications to machine learning.
Compressed sensing with subsampled unitary matrices benefits from optimized sampling schemes, which feature improved theoretical guarantees and empirical performance relative to uniform subsampling. We provide, in a first of its kind in compressed sensing, theoretical guarantees showing that the error caused by the measurement noise vanishes with an increasing number of measurements for optimized sampling schemes, assuming that the noise is Gaussian. We moreover provide similar guarantees for measurements sampled with-replacement with arbitrary probability weights. All our results hold on prior sets contained in a union of low-dimensional subspaces. Finally, we demonstrate that this denoising behavior appears in empirical experiments with a rate that closely matches our theoretical guarantees when the prior set is the range of a generative ReLU neural network and when it is the set of sparse vectors.
ReLU matrix decomposition (RMD) is the following problem: given a sparse, nonnegative matrix X and a factorization rank r, identify a rank-r matrix Θ such that X≈max(0,Θ). RMD is a particular instance of nonlinear matrix decomposition (NMD) that finds application in data compression, matrix completion with entries missing not at random, and manifold learning. The standard RMD model minimizes the least squares error, that is, X - max(0,Θ)_F^2. The corresponding optimization problem, Least-Squares RMD (LS-RMD), is nondifferentiable and highly nonconvex. This motivated Saul to propose an alternative model, dubbed Latent-RMD, where a latent variable Z is introduced and satisfies max(0,Z)=X while minimizing Z - Θ_F^2 (“A nonlinear matrix decomposition for mining the zeros of sparse data”, SIAM J. Math. Data Sci., 2022). Our first contribution is to show that the two formulations may yield different low-rank solutions Θ. We then consider a reparametrization of the Latent-RMD, called 3B-RMD, in which Θ is substituted by a low-rank product WH, where W has r columns and H has r rows. Our second contribution is to prove the convergence of a block coordinate descent (BCD) approach applied to 3B-RMD. Our third contribution is a novel extrapolated variant of BCD, dubbed eBCD, which we prove is also convergent under mild assumptions. We illustrate the significant acceleration effect of eBCD compared to eBCD, and also show that eBCD performs well against the state of the art on synthetic and real-world data sets.
Equivariant neural networks play a pivotal role in analyzing datasets with symmetry properties, particularly in complex data structures. However, integrating equivariance with Markov properties presents notable challenges due to the inherent dependencies within such data. Previous research has primarily concentrated on establishing generalization bounds under the assumption of independently and identically distributed data, frequently neglecting the influence of Markov dependencies. In this study, we investigate the impact of Markov properties on generalization performance alongside the role of equivariance within this context. We begin by applying a new McDiarmid's inequality to derive a generalization bound for neural networks trained on Markov datasets, using Rademacher complexity as a central measure of model capacity. Subsequently, we utilize group theory to compute the covering number under equivariant constraints, enabling us to obtain an upper bound on the Rademacher complexity based on this covering number. This bound provides practical insights into selecting low-dimensional irreducible representations, enhancing generalization performance for fixed-width equivariant neural networks.
We introduce a new framework for data denoising, partially inspired by martingale optimal transport. For a given noisy distribution (the data), our approach involves finding the closest distribution to it among all distributions which 1) have a particular prescribed structure (expressed by requiring they lie in a particular domain), and 2) are self-consistent with the data. We show that this amounts to maximizing the variance among measures in the domain which are dominated in convex order by the data. For particular choices of the domain, this problem and a relaxed version of it, in which the self-consistency condition is removed, are intimately related to various classical approaches to denoising. We prove that our general problem has certain desirable features: solutions exist under mild assumptions, have certain robustness properties, and, for very simple domains, coincide with solutions to the relaxed problem. We also introduce a novel relationship between distributions, termed Kantorovich dominance, which retains certain aspects of the convex order while being a weaker, more robust, and easier-to-verify condition. Building on this, we propose and analyze a new denoising problem by substituting the convex order in the previously described framework with Kantorovich dominance. We demonstrate that this revised problem shares some characteristics with the full convex order problem but offers enhanced stability, greater computational efficiency, and, in specific domains, more meaningful solutions. Finally, we present simple numerical examples illustrating solutions for both the full convex order problem and the Kantorovich dominance problem.
Layer factorization has emerged as a widely used technique for training memory-efficient neural networks. However, layer factorization methods face several challenges, particularly a lack of robustness during the training process. To overcome this limitation, dynamical low-rank training methods have been developed, utilizing robust time integration techniques for low-rank matrix differential equations. Although these approaches facilitate efficient training, they still depend on computationally intensive QR and singular value decompositions of matrices with small rank. In this work, we introduce a novel low-rank training method that reduces the number of required QR decompositions. Our approach integrates an augmentation step into a projector-splitting scheme, ensuring convergence to a locally optimal solution. We provide a rigorous theoretical analysis of the proposed method and demonstrate its effectiveness across multiple benchmarks.
We establish the theoretical framework for implementing the maximumn entropy on the mean (MEM) method for linear inverse problems in the setting of approximate (data-driven) priors. We prove a.s. convergence for empirical means and further develop general estimates for the difference between the MEM solutions with different priors μ and ν based upon the epigraphical distance between their respective log-moment generating functions. These estimates allow us to establish a rate of convergence in expectation for empirical means. We illustrate our results with denoising on MNIST and Fashion-MNIST data sets.
Sampling from high-dimensional distributions has wide applications in data science and machine learning but poses significant computational challenges. We introduce Subspace Langevin Monte Carlo (SLMC), a novel and efficient sampling method that generalizes random-coordinate Langevin Monte Carlo and preconditioned Langevin Monte Carlo by projecting the Langevin update onto subsampled eigenblocks of a time-varying preconditioner at each iteration. The advantage of SLMC is its superior adaptability and computational efficiency compared to traditional Langevin Monte Carlo and preconditioned Langevin Monte Carlo. Using coupling arguments, we establish error guarantees for SLMC and demonstrate its practical effectiveness through a few experiments on sampling from ill-conditioned distributions.
Root mean square propagation (abbreviated as RMSProp) is a first-order stochastic algorithm used in machine learning widely. In this paper, a stable gradient-adjusted RMSProp (abbreviated as SGA-RMSProp) with mini-batch stochastic gradient is proposed, and its properties are studied on the linear least squares problem. R-linear convergence of the algorithm is established on the consistent linear least squares problem. The algorithm is also proved to converge R-linearly to a neighborhood of the minimizer for the inconsistent case, with the region of the neighborhood being controlled by the batch size. Furthermore, numerical experiments are conducted to compare the performances of SGA-RMSProp, stochastic gradient descent (abbreviated as SGD), and the original RMSProp with different batch sizes. The faster initial convergence rate of SGA-RMSProp is observed through numerical experiments and an adaptive strategy for switching from SGA-RMSProp to SGD is proposed, which combines the benefits of these two algorithms.
We introduce a machine-learning framework to learn the hyperparameter sequence of first-order methods (e.g., the step sizes in gradient descent) to quickly solve parametric convex optimization problems. Our computational architecture amounts to running fixed-point iterations where the hyperparameters are the same across all parametric instances and consists of two phases. In the first step-varying phase the hyperparameters vary across iterations, while in the second steady-state phase the hyperparameters are constant across iterations. Our learned optimizer is flexible in that it can be evaluated on any number of iterations and is guaranteed to converge to an optimal solution. To train, we minimize the mean square error to a ground truth solution. In the case of gradient descent, the one-step optimal step size is the solution to a least squares problem, and in the case of unconstrained quadratic minimization, we can compute the two and three-step optimal solutions in closed-form. In other cases, we backpropagate through the algorithm steps to minimize the training objective after a given number of steps. We show how to learn hyperparameters for several popular algorithms: gradient descent, proximal gradient descent, and two ADMM-based solvers: OSQP and SCS. We use a sample convergence bound to obtain generalization guarantees for the performance of our learned algorithm for unseen data, providing both lower and upper bounds. We showcase the effectiveness of our method with many examples, including ones from control, signal processing, and machine learning. Remarkably, our approach is highly data-efficient in that we only use 10 problem instances to train the hyperparameters in all of our examples.
We introduce novel estimators for computing the curvature, tangent spaces, and dimension of data from manifolds, using tools from diffusion geometry. Although classical Riemannian geometry is a rich source of inspiration for geometric data analysis and machine learning, it has historically been hard to implement these methods in a way that performs well statistically. Diffusion geometry lets us develop Riemannian geometry methods that are accurate and, crucially, also extremely robust to noise and low-density data. The methods we introduce here are comparable to the existing state-of-the-art on ideal dense, noise-free data, but significantly outperform them in the presence of noise or sparsity. In particular, our dimension estimate improves on the existing methods on a challenging benchmark test when even a small amount of noise is added. Our tangent space and scalar curvature estimates do not require parameter selection and substantially improve on existing techniques.
In a reinforcement learning (RL) setting, the agent's optimal strategy heavily depends on her risk preferences and the underlying model dynamics of the training environment. These two aspects influence the agent's ability to make well-informed and time-consistent decisions when facing testing environments. In this work, we devise a framework to solve robust risk-aware RL problems where we simultaneously account for environmental uncertainty and risk with a class of dynamic robust distortion risk measures. Robustness is introduced by considering all models within a Wasserstein ball around a reference model. We estimate such dynamic robust risk measures using neural networks by making use of strictly consistent scoring functions, derive policy gradient formulae using the quantile representation of distortion risk measures, and construct an actor-critic algorithm to solve this class of robust risk-aware RL problems. We demonstrate the performance of our algorithm on a portfolio allocation example.
This paper studies Kernel density estimation for a high-dimensional distribution $\rho(x)$. Traditional approaches have focused on the limit of large number of data points $n$ and fixed dimension $d$. We analyze instead the regime where both the number $n$ of data points $y_i$ and their dimensionality $d$ grow with a fixed ratio $\alpha=(\log n)/d$. Our study reveals three distinct statistical regimes for the kernel-based estimate of the density $\hat \rho_h^{\mathcal {D}}(x)=\frac{1}{n h^d}\sum_{i=1}^n K\left(\frac{x-y_i}{h}\right)$, depending on the bandwidth $h$: a classical regime for large bandwidth where the Central Limit Theorem (CLT) holds, which is akin to the one found in traditional approaches. Below a certain value of the bandwidth, $h_{CLT}(\alpha)$, we find that the CLT breaks down. The statistics of $\hat \rho_h^{\mathcal {D}}(x)$ for a fixed $x$ drawn from $\rho(x)$ is given by a heavy-tailed distribution (an alpha-stable distribution). In particular below a value $h_G(\alpha)$, we find that $\hat \rho_h^{\mathcal {D}}(x)$ is governed by extreme value statistics: only a few points in the database matter and give the dominant contribution to the density estimator. We provide a detailed analysis for high-dimensional multivariate Gaussian data. We show that the optimal bandwidth threshold based on Kullback-Leibler divergence lies in the new statistical regime identified in this paper. Our findings reveal limitations of classical approaches, show the relevance of these new statistical regimes, and offer new insights for Kernel density estimation in high-dimensional settings.
Estimating nonlinear functionals of probability distributions from samples is a fundamental statistical problem. The "plug-in" estimator obtained by applying the target functional to the empirical distribution of samples is biased. Resampling methods such as the bootstrap derive artificial datasets from the original one by resampling. Comparing the outcome of the plug-in estimator in the original and resampled datasets allows estimating and thus correcting the bias. In the asymptotic setting, iterations of this procedure attain an arbitrarily high order of bias correction, but finite sample results are scarce. This work develops a new theoretical understanding of bootstrap bias correction by viewing it as an iterative linear solver for the combinatorial operation of Möbius inversion. It sharply characterizes the regime of linear convergence of the bootstrap bias reduction for moment polynomials. It uses these results to show its superalgebraic convergence rate for band-limited functionals. Finally, it derives a modified bootstrap iteration enabling the unbiased estimation of unknown order-m moment polynomials in m bootstrap iterations.
A Gaussian Cox process is a popular model for point process data, in which the intensity function is a transformation of a Gaussian process. Posterior inference of this intensity function involves an intractable integral (i.e., the cumulative intensity function) in the likelihood resulting in doubly intractable posterior distribution. Here, we propose a nonparametric Bayesian approach for estimating the intensity function of an inhomogeneous Poisson process without reliance on large data augmentation or approximations of the likelihood function. We propose to jointly model the intensity and the cumulative intensity function as a transformed Gaussian process, allowing us to directly bypass the need of approximating the cumulative intensity function in the likelihood. We propose an exact MCMC sampler for posterior inference and evaluate its performance on simulated data. We demonstrate the utility of our method in three real-world scenarios including temporal and spatial event data, as well as aggregated time count data collected at multiple resolutions. Finally, we discuss extensions of our proposed method to other point processes.
Dimensionality reduction on quadratic manifolds augments linear approximations with quadratic correction terms. Previous works rely on linear approximations given by projections onto the first few leading principal components of the training data; however, linear approximations in subspaces spanned by the leading principal components alone can miss information that are necessary for the quadratic correction terms to be efficient. In this work, we propose a greedy method that constructs subspaces from leading as well as later principal components so that the corresponding linear approximations can be corrected most efficiently with quadratic terms. Properties of the greedily constructed manifolds allow applying linear algebra reformulations so that the greedy method scales to data points with millions of dimensions. Numerical experiments demonstrate that an orders of magnitude higher accuracy is achieved with the greedily constructed quadratic manifolds compared to manifolds that are based on the leading principal components alone.
We introduce a novel approach for learning memory kernels in Generalized Langevin Equations. This approach initially utilizes a regularized Prony method to estimate correlation functions from trajectory data, followed by regression over a Sobolev norm-based loss function with RKHS regularization. Our approach guarantees improved performance within an exponentially weighted L^2 space, with the kernel estimation error controlled by the error in estimated correlation functions. We demonstrate the superiority of our estimator compared to other regression estimators that rely on L^2 loss functions and also an estimator derived from the inverse Laplace transform, using numerical examples that highlight its consistent advantage across various weight parameter selections. Additionally, we provide examples that include the application of force and drift terms in the equation.
