Aug 2026· SIAM Journal on Scientific Computing· Vol 48, pp. C1037-C1059· 0 citations· 7 references
Abstract
Abstract.
We propose an adaptive randomized truncation estimator for Krylov subspace methods that optimizes the trade-off between the solution variance and the computational cost while remaining unbiased. The estimator solves a constrained optimization problem to compute the truncation probabilities on the fly with minimal computational overhead. The problem has a closed-form solution when the improvement of the deterministic algorithm satisfies a diminishing returns property. We prove that obtaining the optimal adaptive truncation distribution is impossible in the general case. Without the diminishing return condition, our estimator provides a suboptimal but still unbiased solution. We present experimental results in Gaussian process (GP) hyperparameter training and competitive physics-informed neural networks problems to demonstrate the effectiveness of our approach.
Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/RockyL7/AdaptivelySubsampledKrylov.jl and in the supplementary materials ( AdaptivelySubsampledKrylov_jl-master.zip [14.5KB]). [Formula: see text]
Resolvent Monte Carlo estimates eigenvalues of large matrices by sampling Markov chains and reading the target value off a truncated resolvent quotient, trading exact arithmetic for a stochastic error that the almost-optimal sampling scheme is designed to suppress. This paper studies when that error vanishes outright. An exact closed-form identity is derived for the variance of the moment estimators of a general, possibly signed matrix, and is used to isolate a hierarchy of zero-variance notions ranging from the most local, which constrains only the first draws, through the finite-truncation regime that a practical run can certify, to the global regime in which every moment estimator is deterministic. Determinism of the estimator is separated from correctness of the eigenvalue it reports, and the exact conditions under which each notion holds are exhibited, together with the examples that separate them. A single edgewise condition, termed the eigen-triple condition, forces the truncated quotient to equal the target eigenvalue in finite samples; the associated moment and quotient variances are second order in the maximal edge defect and vanish at the eigen-triple. A linear-time procedure certifies the condition.
Tsvetelina Kostadinov, I. Dimov· Mathematics· 0 citations
In this paper, we propose a novel tamed stochastic gradient Hamiltonian Monte Carlo (tSGHMC) algorithm for sampling and stochastic optimization problems with superlinearly growing stochastic gradients. Under a certain continuity in average condition and a strong convexity condition, we establish a non-asymptotic error bound in Wasserstein-2 distance for tSGHMC with the rate of convergence equal to $1/4$. Then, we derive an upper estimate for the associated expected excess risk, which provides a theoretical guarantee for the performance of tSGHMC. To illustrate the effectiveness of the proposed algorithm, we apply tSGHMC to practical examples, including a newsvendor problem and a Conditional Value-at-Risk minimization problem, using synthetic and real-world datasets. Numerical results support our theoretical findings. Furthermore, we compare tSGHMC with its first-order counterpart, namely, the tamed unadjusted stochastic Langevin algorithm. Simulation results demonstrate that tSGHMC achieves lower root mean square error and expected excess risk across a range of tasks.
Data scarcity poses a fundamental challenge in training generative models to produce initial guesses for parametric optimization problems that are otherwise numerically expensive to solve. We therefore study a $k$-neighborhood data collection strategy that augments datasets of converged solutions with intermediate solver iterates, increasing the amount of training data without additional solver runs. To understand the benefits of this approach, we derive a generalization bound based on Rademacher complexity that reveals the role of the $k$-neighborhoods and related parameters. To achieve this result, we focus on one-sided box-constrained quadratic programs solved by projected gradient descent. We illustrate the behavior of this solver on two examples. The approach proposed in this paper enables a more capable DDDAS paradigm by improving the efficiency of the data-model-optimization loop. We finish by discussing two views of learning solver-iterate data and connect our analysis with GLENS, a new data-efficient global search method.
We consider constrained optimization problems in which input data are affected by estimation errors. In such settings, Wasserstein distributionally robust optimization provides a principled framework to mitigate model risk by optimizing against worst-case distributions within Wasserstein ambiguity sets. However, the numerical resolution of the resulting problems remains challenging, especially in constrained and combinatorial settings. In this paper, we propose a tractable stochastic approach based on two key ingredients: (i) an entropic regularization of the distributionally robust value function, which makes it possible to compute stochastic gradient estimators, and (ii) the combination of these estimators with a stochastic Frank-Wolfe algorithm, allowing us to optimize the regularized robust objective while naturally handling constraints. We illustrate the method, and its interests against empirical risk minimization, on two classical optimization problems, the minimum quadratic spanning tree and the traffic assignment problems. Our approach provides a general, practical way to address Wasserstein distributionally robust formulations in the presence of constraints.
We present an unbiased method for Bayesian posterior means based on kinetic Langevin dynamics that combines advanced splitting methods with enhanced gradient approximations. Our approach avoids Metropolis correction by coupling Markov chains at different discretization levels in a multilevel Monte Carlo approach. Theoretical analysis demonstrates that our proposed estimator is unbiased, attains finite variance, and satisfies a central limit theorem. It can achieve accuracy ϵ>0 for estimating expectations of Lipschitz functions in d dimensions with O(d1/4ϵ−2) expected gradient evaluations, without assuming warm start. We exhibit similar bounds using both approximate and stochastic gradients, and our method’s computational cost is shown to scale independently of the size of the dataset. The proposed method is tested using a multinomial regression problem on the MNIST dataset and a Poisson regression model for soccer scores. Experiments indicate that the number of gradient evaluations per effective sample is independent of dimension, even when using inexact gradients. For product distributions, we give dimension-independent variance bounds. Our results demonstrate that in large-scale applications, the unbiased algorithm we present can be 2–3 orders of magnitude more efficient than the “gold-standard” randomized Hamiltonian Monte Carlo.
Neil K. Chada, B. Leimkuhler, Daniel Paulin et al.· Annals of Statistics· 0 citations
We propose a nonparametric method to learn the Lévy density from data consisting of the process’s probability densities. We recast the problem as identifying the kernel of a nonlocal integral operator from discrete or noisy data, which leads to an ill-posed inverse problem. To regularize it, we construct an adaptive reproducing kernel Hilbert space (RKHS) whose kernel is built directly from the data. Under source and spectral decay conditions, we show that the reconstruction error decays with the mesh size at a near-optimal rate. Importantly, we develop a generalized singular value decomposition-based bilevel optimization algorithm to select the regularization parameter, resulting in efficient and robust computation of the regularized estimator. Numerical experiments for several Lévy densities, drift fields and data types (PDE-based densities and sample ensemble-based kernel density estimation reconstructions) demonstrate that our bilevel RKHS method provides a more stable and competitive alternative to classical L-curve and generalized cross-validation strategies and that the adaptive RKHS norm is more accurate and robust than Lρ2- and ℓ2-norms for regularization.