Papers

Debiasing group graphical lasso estimates enables statistical inference when multiple Gaussian graphical models share a common sparsity pattern. We analyze the estimation properties of group graphical lasso, establishing convergence rates and model selection consistency under irrepresentability conditions. Based on these results, we construct debiased estimators that are asymptotically Gaussian, allowing hypothesis testing for linear combinations of precision matrix entries across populations. We also investigate regimes where irrepresentibility conditions does not hold, showing that consistency can still be attained in moderately high-dimensional settings. Simulation studies confirm the theoretical results, and applications to real datasets demonstrate the practical utility of the method.

We study the decomposability and the subdifferential of the tensor nuclear norm. Both concepts are well understood and widely applied in matrices but remain unclear for higher-order tensors. We show that the tensor nuclear norm admits a full decomposability over specific subspaces and determine the largest possible subspaces that allow the full decomposability. We derive novel inclusions of the subdifferential of the tensor nuclear norm and study its subgradients in a variety of subspaces of interest. All the results hold for tensors of an arbitrary order. As an immediate application, we establish the statistical performance of the tensor robust principal component analysis, the first such result for tensors of an arbitrary order.

Neural networks are typically optimized with variants of stochastic gradient descent. Under a squared loss, however, the optimal solution to the linear last layer weights is known in closed-form. We propose to leverage this during optimization, treating the last layer as a function of the backbone parameters, and optimizing solely for these parameters. We show this is equivalent to alternating between gradient descent steps on the backbone and closed-form updates on the last layer. We adapt the method for the setting of stochastic gradient descent, by trading off the loss on the current batch against the accumulated information from previous batches. We provide theoretical analyses showing convergence of the method to an optimal solution in the neural tangent kernel regime, as well as quantifying the gains compared to standard SGD in a one-step analysis. Finally, we demonstrate the effectiveness of our approach compared with SGD and Adam on a squared loss in several regression tasks, including neural operators and causal inference.

Wasserstein barycenters provide a principled approach for aggregating probability measures, while preserving the geometry of their ambient space. Existing discrete methods are not because as they assume access to the complete set of samples from the input measures. Meanwhile, neural network approaches do scale well, but rely on complex optimization problems and cannot easily incorporate label information. We address these limitations through gradient flows in the space of probability measures. Through time discretization, we achieve a scalable algorithm that i) relies on mini-batch optimal transport, ii) accepts modular regularization through task-aware functions, and iii) seamlessly integrates supervised information into the ground-cost. We empirically validate our approach on domain adaptation benchmarks that span computer vision, neuroscience, and chemical engineering. Our method establishes a new state-of-the-art Wasserstein barycenter solver, with labeled barycenters consistently outperforming unlabeled ones. Our code at https://github.com/SigmaNova/barycentric-gradient-flows

We propose a reflection-free Langevin framework for sampling and optimization on compact polyhedra. The method is based on the inverse Hessian of the logarithmic barrier, which defines a Dikin–Langevin diffusion whose drift and noise adapt to the local interior-point geometry. We show that trajectories started in the interior remain feasible for all finite times almost surely, so the constrained domain is preserved without reflections or projections. For computation, we discretize the diffusion using the Euler–Maruyama scheme and apply a Metropolis–Hastings correction, yielding a sampler that targets the exact constrained distribution. We also propose an annealed interacting variant for nonconvex optimization. Numerically, the Metropolis-adjusted method outperforms both the Dikin random walk and standard MALA on anisotropic box-constrained Gaussians, and the interacting optimizer escapes suboptimal basins more reliably than the non-interacting method.

The Busemann function has recently found much interest in a variety of geometric machine learning problems, as it naturally defines projections onto geodesic rays of Riemannian manifolds and generalizes the notion of hyperplanes. As several sources of data can be conveniently modeled as probability distributions, it is natural to study this function in the Wasserstein space, which carries a rich formal Riemannian structure induced by Optimal Transport metrics. In this work, we investigate the existence and computation of Busemann functions in Wasserstein space, which admits geodesic rays. We establish closed-form expressions in two important cases: one-dimensional distributions and Gaussian measures. These results enable explicit projection schemes for probability distributions on $\mathbb{R}$, which in turn allow us to define novel Sliced-Wasserstein distances over Gaussian mixtures and labeled datasets. We demonstrate the efficiency of those original schemes on synthetic datasets as well as transfer learning problems.

Deep generative models have made significant advances in generating complex content, yet conditional generation remains a fundamental challenge. Existing conditional generative adversarial networks often struggle to balance the dual objectives of assessing authenticity and conditional alignment of input samples within their conditional discriminators. To address this, we propose a novel discriminator design that integrates three key capabilities: unconditional discrimination, matching-aware supervision to enhance alignment sensitivity, and adaptive weighting to dynamically balance all objectives. Specifically, we introduce Sum of Naturalness and Alignment (SONA), which employs separate projections for naturalness (authenticity) and alignment in the final layer with an inductive bias, supported by dedicated objective functions and an adaptive weighting mechanism. Extensive experiments on class-conditional generation tasks show that \ours achieves superior sample quality and conditional alignment compared to state-of-the-art methods. Furthermore, we demonstrate its effectiveness in text-to-image generation, confirming the versatility and robustness of our approach.

Maintaining the predictive performance of pricing models is challenging when insurance portfolios and data-generating mechanisms evolve over time. Focusing on non-life insurance, we adopt the concept-drift terminology from machine learning and distinguish virtual drift from real concept drift in an actuarial setting. Methodologically, we (i) formalize deviance loss and Murphy's score decomposition to assess global and local auto-calibration; (ii) study the Gini score as a rank-based performance measure, derive its asymptotic distribution, and develop a consistent bootstrap estimator of its asymptotic variance; and (iii) combine these results into a statistically grounded, model-agnostic monitoring framework that integrates a Gini-based ranking drift test with global and local auto-calibration tests. An application to a modified motor insurance portfolio with controlled concept-drift scenarios illustrates how the framework guides decisions on refitting or recalibrating pricing models.

In-context learning (ICL) is a key building block of modern large language models, yet its theoretical mechanisms remain poorly understood. It is particularly mysterious how ICL operates in real-world applications where tasks have a common structure. In this work, we address this problem by analyzing a linear attention model trained on low-rank regression tasks. Within this setting, we precisely characterize the distribution of predictions and the generalization error in the high-dimensional limit. Moreover, we find that statistical fluctuations in finite pre-training data induce an implicit regularization. Finally, we identify a sharp phase transition of the generalization error governed by task structure. These results provide a framework for understanding how transformers learn to learn the task structure.

Accurate predictions on tabular data rely on capturing complex, dataset-specific feature interactions. Attention-based methods and graph neural networks, referred to as graph-based tabular deep learning (GTDL), aim to improve predictions by modeling these interactions as a graph. In this work, we analyze how these methods model the feature interactions. Current GTDL approaches primarily focus on optimizing predictive accuracy, often neglecting the accurate modeling of the underlying graph structure. Using synthetic datasets with known ground-truth graph structures, we find that current GTDL methods fail to recover meaningful feature interactions, as their edge recovery is close to random. This suggests that the attention mechanism and message-passing schemes used in GTDL do not effectively capture feature interactions. Furthermore, when we impose the true interaction structure, we find that the predictive accuracy improves. This highlights the need for GTDL methods to prioritize accurate modeling of the graph structure, as it leads to better predictions.

Masked diffusion models have shown promising performance in generating high-quality samples in a wide range of domains, but accelerating their sampling process remains relatively underexplored. To investigate efficient samplers for masked diffusion, this paper theoretically analyzes the MaskGIT sampler for image modeling, revealing its implicit temperature sampling mechanism. Through this analysis, we show that MaskGIT is asymptotically equivalent to a choose-then-sample (CTS) formulation, instantiated as the "moment sampler," which explicitly separates index selection from token sampling. This CTS reformulation is essential: it yields unbiased token sampling and exposes an algorithmic design space for index selection, both of which are inaccessible in MaskGIT's original formulation. Regarding token sampling, we reveal that MaskGIT implicitly adopts a low-temperature sampler, which explains why MaskGIT often degrades with more sampling steps. The CTS reformulation of MaskGIT allows us to fix the temperature sampling to ensure unbiasedness. We also improve the index selection in CTS through two key innovations: a partial caching technique for transformers that approximates longer sampling trajectories without proportional computational cost, and a hybrid approach formalizing the exploration-exploitation trade-off in adaptive unmasking. Experiments in image and text domains demonstrate our theory as well as the efficiency of our proposed methods, advancing both theoretical understanding and practical implementation of masked diffusion samplers.

Recent multimodal LLMs have shown promise in chart-based visual question answering, but their performance declines sharply on unannotated charts-those requiring precise visual interpretation rather than relying on textual shortcuts. To address this, we introduce ChartAgent, a novel agentic framework that explicitly performs visual reasoning directly within the chart's spatial domain. Unlike textual chain-of-thought reasoning, ChartAgent iteratively decomposes queries into visual subtasks and actively manipulates and interacts with chart images through specialized actions such as drawing annotations, cropping regions (e.g., segmenting pie slices, isolating bars), and localizing axes, using a library of chart-specific vision tools to fulfill each subtask. This iterative reasoning process closely mirrors human cognitive strategies for chart comprehension. ChartAgent achieves state-of-the-art accuracy on the ChartBench and ChartX benchmarks, surpassing prior methods by up to 16.07% absolute gain overall and 17.31% on unannotated, numerically intensive queries. Furthermore, our analyses show that ChartAgent is (a) effective across diverse chart types, (b) achieves the highest scores across varying visual and reasoning complexity levels, and (c) serves as a plug-and-play framework that boosts performance across diverse underlying LLMs. Our work is among the first to demonstrate visually grounded reasoning for chart understanding using tool-augmented multimodal agents.

We study the design of experiments with multiple treatment levels, a setting common in clinical trials and online A/B/n testing. Unlike single-treatment studies, practical analyses of multi-treatment experiments typically first select a winning treatment, and then only estimate the effect therein. Motivated by this analysis paradigm, we propose a design for MUlti-treatment experiments that jointly maximizes the accuracy of winner Selection and effect Estimation (MUSE). Explicitly, we introduce a single objective that balances selection and estimation, and determine the unit allocation to treatments and control by optimizing this objective. Theoretically, we establish finite-sample guarantees and asymptotic equivalence between our proposal and the Neyman allocation for the true optimal treatment and control. Across simulations and a real data application, our method performs favorably in both selection and estimation compared to various standard alternatives.

We survey different perspectives on the stochastic localization process of Eldan, a powerful construction that has had many exciting recent applications in high-dimensional probability and algorithm design. Unlike prior surveys on this topic, our focus is on giving a self-contained presentation of all known alternative constructions of Eldan's stochastic localization, with an emphasis on connections between different constructions. Our hope is that by collecting these perspectives, some of which had primarily arisen within a particular community (e.g., probability theory, theoretical computer science, information theory, or machine learning), we can broaden the accessibility of stochastic localization, and ease its future use.

In classical canonical correlation analysis (CCA), the goal is to determine the linear transformations of two random vectors into two new random variables that are most strongly correlated. Canonical variables are pairs of these new random variables, while canonical correlations are correlations between these pairs. In this paper, we propose and study two generalizations of this classical method: (1) Instead of two random vectors we study more complex data structures that appear in important applications. In these structures, there are $L$ features, each described by $p_l$ scalars, $1 \le l \le L$. We observe $n$ such objects over $T$ time points. We derive a suitable analog of the CCA for such data. Our approach relies on embeddings into Reproducing Kernel Hilbert Spaces, and covers several related data structures as well. (2) We develop an analogous approach for multidimensional random processes. In this case, the experimental units are multivariate continuous, square-integrable functions over a given interval. These functions are modeled as elements of a Hilbert space, so in this case, we define the multiple functional canonical correlation analysis, MFCCA. We justify our approaches by their application to two data sets and suitable large sample theory. We derive consistency rates for the related transformation and correlation estimators, and show that it is possible to relax two common assumptions on the compactness of the underlying cross-covariance operators and the independence of the data.

Data-driven inverse optimization for mixed-integer linear programs (MILPs), which seeks to learn an objective function and constraints consistent with observed decisions, is important for building accurate mathematical models in a variety of domains, including power systems and scheduling. However, to the best of our knowledge, existing data-driven inverse optimization methods primarily focus on learning objective functions under known constraints, and learning both objective functions and constraints from data for MILPs remains largely unexplored. In this paper, we propose a two-stage approach for a class of inverse optimization problems in which the objective is a linear combination of given feature functions and the constraints are parameterized by unknown functions and thresholds. Our method first learns the constraints and then, conditioned on the learned constraints, estimates the objective-function weights. On the theoretical side, we provide finite-sample guarantees for solving the proposed inverse optimization problem. To this end, we develop statistical learning tools for pseudo-metric spaces under sub-Gaussian assumptions and use them to derive a learning-theoretic framework for inverse optimization with both unknown objectives and constraints. On the experimental side, we demonstrate that our method successfully solves inverse optimization problems on scheduling instances formulated as ILPs with up to 100 decision variables.

Domain generalization (DG) is the problem of generalizing from several distributions (or domains), for which labeled training data are available, to a new test domain for which no labeled data is available. For the prevailing benchmark datasets in DG, there exists a single classifier that performs well across all domains. In this work, we study a fundamentally different regime where the domains satisfy a posterior drift assumption, in which the optimal classifier might vary substantially with domain. We establish a decision-theoretic framework for DG under posterior drift, and investigate the practical implications of this framework through experiments on language and vision tasks.

We present spd-metrics-id, a Python package for computing distances and divergences between symmetric positive-definite (SPD) matrices. Unlike traditional toolkits that focus on specific applications, spd-metrics-id provides a unified, extensible, and reproducible framework for SPD distance computation. The package supports a wide variety of geometry-aware metrics, including Alpha-z Bures-Wasserstein, Alpha-Procrustes, affine-invariant Riemannian, log-Euclidean, and others, and is accessible both via a command-line interface and a Python API. Reproducibility is ensured through Docker images and Zenodo archiving. We illustrate usage through a connectome fingerprinting example, but the package is broadly applicable to covariance analysis, diffusion tensor imaging, and other domains requiring SPD matrix comparison. The package is openly available at https://pypi.org/project/spd-metrics-id/.

We introduce the Divergence Phase Index (DPI), a novel framework for quantifying phase differences in one and multidimensional signals, grounded in harmonic analysis via the Riesz transform. Based on classical Hilbert Transform phase measures, the DPI extends these principles to higher dimensions, offering a geometry-aware metric that is invariant to intensity scaling and sensitive to structural changes. We applied this method on both synthetic and real-world datasets, including intracranial EEG (iEEG) recordings during epileptic seizures, high-resolution microscopy images, and paintings. In the 1D case, the DPI robustly detects hypersynchronization associated with generalized epilepsy, while in 2D, it reveals subtle, imperceptible changes in images and artworks. Additionally, it can detect rotational variations in highly isotropic microscopy images. The DPI's robustness to amplitude variations and its adaptability across domains enable its use in diverse applications from nonlinear dynamics, complex systems analysis, to multidimensional signal processing.

Survival analysis provides statistical methods to model the time until an event occurs. Reporting delays arise when event times are not observed at their occurrence but are only revealed upon reporting. This issue is particularly critical for timely risk evaluation when the observation window is short due to administrative censoring. In this study, we incorporate right-censored reporting delays by jointly modeling parametric hazards for the event and reporting processes. We then construct a consistent estimator for the model parameters and develop a Monte Carlo expectation-maximization algorithm to compute it. To address the challenges posed by administrative censoring, we leverage these findings and propose a transfer-learning procedure. Experimental results demonstrate that our method improves the accuracy of timely risk evaluation under administrative censoring.

Conformal prediction offers finite-sample coverage guarantees under minimal assumptions. However, existing methods treat the entire modeling process as a black box, overlooking opportunities to exploit and understand modular structure. We introduce a conformal prediction framework for two-stage sequential models, where an upstream predictor generates intermediate representations for a downstream model. By decomposing the overall prediction residual into stage-specific components, our method enables practitioners to attribute uncertainty to specific pipeline stages. We develop a risk-controlled parameter selection procedure using family-wise error rate (FWER) control to calibrate stage-wise scaling parameters, and introduce an adaptive extension for non-stationary settings. Experiments on synthetic distribution shifts, as well as real-world supply chain and stock market data, demonstrate that our approach improves coverage under structural, stage-wise shifts compared to standard conformal methods, while identifying stage-wise error contribution. This framework offers diagnostic advantages and robust coverage that standard conformal methods lack.

This paper is a follow-up to our earlier study, Natural Disasters in Canada (2017). We analyze the Canadian Disaster Database (CDD) to examine the frequency and severity of various natural disasters over the past 120 years and to identify emerging trends. We generate annual loss distributions for individual disaster types, as well as an aggregate annual loss distribution across all event types. Our analysis provides evidence that Canada is experiencing warmer and wetter conditions and indicates a substantial likelihood of extreme national-level losses.

The 2024 July Revolution in Bangladesh represents a landmark event in the study of civil resistance. This study investigates the central paradox of the success of this student-led civilian uprising: how state violence, intended to quell dissent, ultimately fueled the movement's victory. We employ a mixed-methods approach. First, we develop a qualitative narrative of the conflict's timeline to generate specific, testable hypotheses. Then, using a disaggregated, event-level dataset, we employ a multi-method quantitative analysis to dissect the complex relationship between repression and mobilisation. We provide a framework to analyse explosive modern uprisings like the July Revolution. Initial pooled regression models highlight the crucial role of protest momentum in sustaining the movement. To isolate causal effects, we specify a Two-Way Fixed Effects panel model, which provides robust evidence for a direct and statistically significant local suppression backfire effect. Our Vector Autoregression (VAR) analysis provides clear visual evidence of an immediate, nationwide mobilisation in response to increased lethal violence. We further demonstrate that this effect was non-linear. A structural break analysis reveals that the backfire dynamic was statistically insignificant in the conflict's early phase but was triggered by the catalytic moral shock of the first wave of lethal violence, and its visuals circulated around July 16th. A complementary machine learning analysis (XGBoost, out-of-sample R$^{2}$=0.65) corroborates this from a predictive standpoint, identifying "excessive force against protesters" as the single most dominant predictor of nationwide escalation. We conclude that the July Revolution was driven by a contingent, non-linear backfire, triggered by specific catalytic moral shocks and accelerated by the viral reaction to the visual spectacle of state brutality.

Arterial aneurysm (Fig.1) is a bulb-shape local expansion of human arteries, the rupture of which is a leading cause of morbidity and mortality in US. Therefore, the prediction of arterial aneurysm rupture is of great significance for aneurysm management and treatment selection. The prediction of aneurysm rupture depends on the analysis of the time series of aneurysm growth history. However, due to the long time scale of aneurysm growth, the time series of aneurysm growth is not always accessible. We here proposed a method to reconstruct the aneurysm growth time series directly from patient parameters. The prediction is based on data pairs of [patient parameters, patient aneurysm growth time history]. To obtain the mapping from patient parameters to patient aneurysm growth time history, we first apply autoencoder to obtain a compact representation of the time series for each patient. Then a mapping is learned from patient parameters to the corresponding compact representation of time series via a five-layer neural network. Moving average and convolutional output layer are implemented to explicitly taking account the time dependency of the time series. Apart from that, we also propose to use prior knowledge about the mechanism of aneurysm growth to improve the time series reconstruction results. The prior physics-based knowledge is incorporated as constraints for the optimization problem associated with autoencoder. The model can handle both algebraic and differential constraints. Our results show that including physical model information about the data will not significantly improve the time series reconstruction results if the training data is error-free. However, in the case of training data with noise and bias error, incorporating physical model constraints can significantly improve the predicted time series.

Meaningful comparison between sets of observations often necessitates alignment or registration between them, and the resulting optimization problems range in complexity from those admitting simple closed-form solutions to those requiring advanced and novel techniques. We compare different Procrustes problems in which we align two sets of points after various perturbations by minimizing the norm of the difference between one matrix and an orthogonal transformation of the other. The minimization problem depends significantly on the choice of matrix norm; we highlight recent developments in nonsmooth Riemannian optimization and characterize which choices of norm work best for each perturbation. We show that in several applications, from low-dimensional alignments to hypothesis testing for random networks, when Procrustes alignment with the spectral or robust norm is the appropriate choice, it is often feasible to replace the computationally more expensive spectral and robust minimizers with their closed-form Frobenius-norm counterpart. Our work reinforces the synergy between optimization, geometry, and statistics.

Variable-ratio matching is a flexible alternative to conventional $1$-to-$k$ matching for designing observational studies that emulate a target randomized controlled trial (RCT). To achieve fine balance – that is, matching treated and control groups to have the same marginal distribution on selected covariates – conventional approaches typically partition the data into strata based on estimated entire numbers and then perform a series of $1$-to-$k$ matches within each stratum, with $k$ determined by the stratum-specific entire number. This “divide-and-conquer" strategy has notable limitations: (1) fine balance typically does not hold in the final pooled sample, and (2) more controls may be discarded than necessary. To address these limitations, we propose a one-shot variable-ratio matching algorithm. Our method produces designs with exact fine balance on selected covariates in the matched sample, mimicking a hypothetical RCT where units are first grouped into sets of different sizes and one unit within each set is assigned to treatment while others to control. Moreover, our method achieves comparable or superior balance across many covariates and retains more controls in the final matched design, compared to the “divide-and-conquer" approach. We demonstrate the advantages of the proposed design over the conventional approach via simulations and using a dataset studying the effect of right heart catheterization on mortality among critically ill patients. The algorithm is implemented in the R package match2C.

Policy targets evolve faster than the Coupled Model Intercomparison Project cycles, complicating adaptation and mitigation planning that must often contend with outdated projections. Climate model output emulators address this gap by offering inexpensive surrogates that can rapidly explore alternative futures while staying close to Earth System Model (ESM) behavior. The focus is on emulators designed to provide inputs to impact models. Using monthly ESM fields of near-surface temperature, precipitation, relative humidity, and wind speed, it is shown that deep generative models have the potential to model the joint distribution of variables relevant for impacts. The specific model proposed uses score-based diffusion on a spherical mesh and runs on a single mid-range graphical processing unit. A thorough suite of diagnostics is introduced to compare emulator outputs with their parent ESMs, including their probability densities, cross-variable correlations, time of emergence, or tail behavior. The emulator performance is evaluated across three distinct ESMs in both pre-industrial and forced regimes. The results show that the emulator produces distributions that closely match the ESM outputs and captures key forced responses. They also reveal important failure cases, notably for variables with a strong regime shift in the seasonal cycle. Although not a perfect match to the ESM, the inaccuracies of the emulator are small relative to the magnitude of internal variability in ESM projections. This suggests that the generative emulators can be useful in supporting impact assessment. Priorities for future development toward daily resolution, finer spatial scales, and bias-aware training are discussed. Code is made available at https://github.com/shahineb/climemu.

Choosing an appropriate learning rate remains a key challenge in scaling depth of modern deep networks. The classical maximal update parameterization ($\mu$P) enforces a fixed per-layer update magnitude, which is well suited to homogeneous multilayer perceptrons (MLPs) but becomes ill-posed in heterogeneous architectures where residual accumulation and convolutions introduce imbalance across layers. We introduce Arithmetic-Mean $\mu$P (AM-$\mu$P), which constrains not each individual layer but the network-wide average one-step pre-activation second moment to a constant scale. Combined with a residual-aware He fan-in initialization - scaling residual-branch weights by the number of blocks ($\mathrm{Var}[W]=c/(K\cdot \mathrm{fan\text{-}in})$) - AM-$\mu$P yields width-robust depth laws that transfer consistently across depths. We prove that, for one- and two-dimensional convolutional networks, the maximal-update learning rate satisfies $\eta^\star(L)\propto L^{-3/2}$; with zero padding, boundary effects are constant-level as $N\gg k$. For standard residual networks with general conv+MLP blocks, we establish $\eta^\star(L)=\Theta(L^{-3/2})$, with $L$ the minimal depth. Empirical results across a range of depths confirm the $-3/2$ scaling law and enable zero-shot learning-rate transfer, providing a unified and practical LR principle for convolutional and deep residual networks without additional tuning overhead.

Traditional conformal prediction methods construct prediction sets such that the true label falls within the set with a user-specified coverage level. However, poorly chosen coverage levels can result in uninformative predictions, either producing overly conservative sets when the coverage level is too high, or empty sets when it is too low. Moreover, the fixed coverage level cannot adapt to the specific characteristics of each individual example, limiting the flexibility and efficiency of these methods. In this work, we leverage recent advances in e-values and post-hoc conformal inference, which allow the use of data-dependent coverage levels while maintaining valid statistical guarantees. We propose to optimize an adaptive coverage policy by training a neural network using a leave-one-out procedure on the calibration set, allowing the coverage level and the resulting prediction set size to vary with the difficulty of each individual example. We support our approach with theoretical coverage guarantees and demonstrate its practical benefits through a series of experiments.

Statistical methods for metric spaces provide a general and versatile framework for analyzing complex data types. We introduce a novel approach for constructing confidence regions around new predictions from any bagged regression algorithm with metric-space-valued responses. This includes the recent extensions of random forests for metric responses: Fréchet random forests (Capitaine et al., 2024), random forest weighted local constant Fréchet regression (Qiu et al., 2024), and metric random forests (Bulté and Sørensen, 2024). Our prediction regions leverage out-of-bag observations generated during a single forest training, employing the entire data set for both prediction and uncertainty quantification. We establish asymptotic guarantees of out-of-bag prediction balls for four coverage types under certain regularity conditions. Moreover, we demonstrate the superior stability and smaller radius of out-of-bag balls compared to split-conformal methods through extensive numerical experiments where the response lies on the Euclidean space, sphere, hyperboloid, and space of positive definite matrices. A real data application illustrates the potential of the confidence regions for quantifying the uncertainty in the study of solar dynamics and the use of data-driven non-isotropic distances on the sphere.