Papers

We study coherence distillation under time-translation-invariant operations: given many copies of a quantum state containing coherence in the energy eigenbasis, the aim is to produce a purer coherent state while respecting the time-translation symmetry. This symmetry ensures that the output remains synchronized with the input and that the process can be realized by energy-conserving unitaries coupling the system to a reservoir initially in an energy eigenstate, thereby modeling thermal operations supplemented by a work reservoir or battery. For qubit systems, we determine the optimal asymptotic fidelity and show that it is governed by the purity of coherence, a measure of asymmetry derived from the right logarithmic derivative (RLD) Fisher information. In particular, we find that the lowest achievable infidelity (one minus fidelity) scales as $1/N$ times the reciprocal of the purity of coherence of each input qubit, where $N$ is the number of copies, giving this quantity a clear operational meaning. We additionally study many other interesting aspects of the coherence distillation problem for qubits, including computing higher-order corrections to the lowest achievable infidelity up to $O(1/N^3)$, and expressing the optimal channel as a boundary value problem that can be solved numerically.

Estimating how individual input variables affect the output of a black-box model is a central task in explainable machine learning. However, existing methods suffer from two key limitations: sensitivity to out-of-distribution (OOD) evaluations, which arises when query points are placed far from the data manifold, and instability under feature correlation, which can lead to unreliable effect estimates in practice. We introduce a unified view of main effect estimation as a design problem, which reveals that all existing methods differ only in their choice of evaluation locations. Building on this formulation, we propose A2D2E, an Estimator based on Accumulated Aggregated D-Optimal Designs, which replaces evaluations with a D-optimal hypercube design to minimize the variance of main effect estimation. A2D2E is model-agnostic, requires no differentiability of the predictor, and admits a closed-form estimator with complexity comparable to existing approaches. We establish that A2D2E is consistent to the same population target as ALE, and extend this result to the realistic setting where only a surrogate model is available. Through extensive simulations across multiple predictive models and dependence settings, we demonstrate that A2D2E outperforms ALE-based methods, with the largest gains under high feature correlation.

Flow models are a cornerstone of modern machine learning. They are generative models that progressively transform probability distributions according to learned dynamics. Specifically, they learn a continuous-time Markov process that efficiently maps samples from a simple source distribution into samples from a complex target distribution. We show that these models are naturally related to the Schrödinger equation, for an unusual Hamiltonian on continuous variables. Moreover, we prove that the dynamics generated by this Hamiltonian can be efficiently simulated on a quantum computer. Together, these results give a quantum algorithm for preparing coherent encodings (a.k.a., qsamples) for a vast family of probability distributions–namely, those expressible by flow models–by reducing the task to an existing classical learning problem, plus Hamiltonian simulation. For statistical problems defined by flow models, such as mean estimation and property testing, this enables the use of quantum algorithms tailored to qsamples, which may offer advantages over classical algorithms based only on samples from a flow model. More broadly, these results reveal a close connection between state-of-the-art machine learning models, such as flow matching and diffusion models, and one of the main expected capabilities of quantum computers: simulating quantum dynamics.

Graph Neural Networks (GNNs) leverage the graph structure to transmit information between nodes, typically through the message-passing mechanism. While these models have found a wide variety of applications, they are known to suffer from over-squashing, where information from a large receptive field of node representations is collapsed into a single fixed sized vector, resulting in an information bottleneck. In this paper, we re-examine the over-squashing phenomenon through the lens of model storage and retrieval capacity, which we define as the amount of information that can be stored in a node's representation for later use. We study some of the limitations of existing tasks used to measure over-squashing and introduce a new synthetic task to demonstrate that an information bottleneck can saturate this capacity. Furthermore, we adapt ideas from the sequence modeling literature on associative memories, fast weight programmers, and the xLSTM model to develop a novel GNN architecture with improved capacity. We demonstrate strong performance of this architecture both on our capacity synthetic task, as well as a range of real-world graph benchmarks.

Cluster-randomized trials (CRTs) are experimental designs where groups or clusters of participants, rather than the individual participants themselves, are randomized to intervention groups. Analyzing CRT requires distinguishing between treatment effects at the cluster level and the individual level, which requires a clear definition of the estimands under the potential outcomes framework. For analyzing survival outcomes, it is common to assess the treatment effect by comparing survival functions or restricted mean survival times between treatment groups. In this article, we formally characterize cluster-level and individual-level treatment effect estimands with right-censored survival outcomes in CRTs and propose doubly robust estimators for targeting such estimands. Under covariate-dependent censoring, our estimators ensure consistency when either the censoring model or the outcome model is correctly specified, but not necessarily both. We explore different modeling options for the censoring and outcome models to estimate the censoring and survival distributions, and investigate a deletion-based jackknife method for variance and interval estimation. Extensive simulations demonstrate that the proposed methods perform adequately in finite samples. Finally, we illustrate our method by analyzing a completed CRT with survival endpoints.

High-dimensional linear contextual bandit problems remain a significant challenge due to the curse of dimensionality. Existing methods typically consider either the model parameters to be sparse or the eigenvalues of context covariance matrices to be (approximately) sparse, lacking general applicability due to the rigidity of conventional reward estimators. To overcome this limitation, a powerful pointwise estimator is introduced in this work that adaptively navigates both kinds of sparsity. Based on this pointwise estimator, a novel algorithm, termed HOPE, is proposed. Theoretical analyses demonstrate that HOPE not only achieves improved regret bounds in previously discussed homogeneous settings (i.e., considering only one type of sparsity) but also, for the first time, efficiently handles two new challenging heterogeneous settings (i.e., considering a mixture of two types of sparsity), highlighting its flexibility and generality. Experiments corroborate the superiority of HOPE over existing methods across various scenarios.

Reinsurance treaty pricing must satisfy stringent regulatory standards, yet current quoting practices remain opaque and difficult to audit. We introduce ClauseLens, a clause-grounded reinforcement learning framework that produces transparent, regulation-compliant, and risk-aware treaty quotes. ClauseLens models the quoting task as a Risk-Aware Constrained Markov Decision Process (RA-CMDP). Statutory and policy clauses are retrieved from legal and underwriting corpora, embedded into the agent's observations, and used both to constrain feasible actions and to generate clause-grounded natural language justifications. Evaluated in a multi-agent treaty simulator calibrated to industry data, ClauseLens reduces solvency violations by 51%, improves tail-risk performance by 27.9% (CVaR_0.10), and achieves 88.2% accuracy in clause-grounded explanations with retrieval precision of 87.4% and recall of 91.1%. These findings demonstrate that embedding legal context into both decision and explanation pathways yields interpretable, auditable, and regulation-aligned quoting behavior consistent with Solvency II, NAIC RBC, and the EU AI Act.

We identify and analyze a surprising phenomenon of Latent Diffusion Models (LDMs) where the final steps of the diffusion can degrade sample quality. In contrast to conventional arguments that justify early stopping for numerical stability, this phenomenon is intrinsic to the dimensionality reduction in LDMs. We provide a principled explanation by analyzing the interaction between latent dimension and stopping time. Under a Gaussian framework with linear autoencoders, we characterize the conditions under which early stopping is needed to minimize the distance between generated and target distributions. More precisely, we show that lower-dimensional representations benefit from earlier termination, whereas higher-dimensional latent spaces require later stopping time. We further establish that the latent dimension interplays with other hyperparameters of the problem such as constraints in the parameters of score matching. Experiments on synthetic and real datasets illustrate these properties, underlining that early stopping can improve generative quality. Together, our results offer a theoretical foundation for understanding how the latent dimension influences the sample quality, and highlight stopping time as a key hyperparameter in LDMs.

In this paper we will give a characterization of the learnability of forgiving 0-1 loss functions in the multiclass setting with effectively finite cardinality of the output and label space. To do this, we create a new combinatorial dimension that is based off of the Natarajan Dimension and we show that a hypothesis class is learnable in our setting if and only if this Generalized Natarajan Dimension is finite. We also show how this dimension characterizes other known learning settings such as a vast amount of instantiations of learning with set-valued feedback and a modified version of list learning.

Micro-randomized trials (MRTs) are increasingly used to evaluate mobile health interventions with binary proximal outcomes. Standard inverse probability weighting (IPW) estimators are unbiased but unstable in small samples or under extreme randomization. Estimated mean excursion effect (EMEE) improves efficiency but lacks double robustness. We propose a doubly robust EMEE (DR-EMEE) with stabilized and truncated weights, combining per-decision IPW and outcome regression. We prove double robustness, asymptotic efficiency, and provide finite-sample variance corrections, with extensions to machine learning nuisance estimators. In simulations, DR-EMEE reduces root mean squared error, improves coverage, and achieves up to twofold efficiency gains over IPW and five to ten percent over EMEE. Applications to HeartSteps, PAMAP2, and mHealth datasets confirm stable and efficient inference across both randomized and observational settings.

Following the wide-spread adoption of machine learning models in real-world applications, the phenomenon of performativity, i.e. model-dependent shifts in the test distribution, becomes increasingly prevalent. Unfortunately, since models are usually trained solely based on samples from the original (unshifted) distribution, this performative shift may lead to decreased test-time performance. In this paper, we study the question of whether and when performative binary classification problems are learnable, via the lens of the classic PAC (Probably Approximately Correct) learning framework. We motivate several performative scenarios, accounting in particular for linear shifts in the label distribution, as well as for more general changes in both the labels and the features. We construct a performative empirical risk function, which depends only on data from the original distribution and on the type performative effect, and is yet an unbiased estimate of the true risk of a classifier on the shifted distribution. Minimizing this notion of performative risk allows us to show that any PAC-learnable hypothesis space in the standard binary classification setting remains PAC-learnable for the considered performative scenarios. We also conduct an extensive experimental evaluation of our performative risk minimization method and showcase benefits on synthetic and real data.

Gene expression levels, hormone secretion, and internal body temperature each oscillate over an approximately 24-hour cycle, or display circadian rhythms. Many circadian biology studies have investigated how these rhythms vary across cohorts, uncovering associations between atypical rhythms and diseases such as cancer, metabolic syndrome, and sleep disorders. A challenge in analyzing circadian biology data is that the oscillation peak and trough times for a phenomenon differ across individuals. If these individual-level differences are not accounted for in trigonometric regression, which is prevalent in circadian biology studies, then estimates of the population-level amplitude parameters can suffer from attenuation bias. This attenuation bias could lead to inaccurate study conclusions. To address attenuation bias, we propose a refined two-stage (RTS) method for trigonometric regression given longitudinal data obtained from each individual participating in a study. In the first stage, the parameters of individual-level models are estimated. In the second stage, transformations of these individual-level estimates are aggregated to produce population-level parameter estimates for inference. Simulation studies show that our RTS method mitigates bias in parameter estimation, obtains greater statistical power, and maintains appropriate type I error control when compared to the standard two-stage (STS) method, which ignores individual-level differences in peak and trough times. The only exception for parameter estimation and statistical power occurs when the oscillation amplitudes are weak relative to random variability in the data and the sample size is small. Illustrations with cortisol level data and heart rate data show that our RTS method obtains larger population-level amplitude parameter estimates and smaller $p$-values for multiple hypothesis tests when compared to the STS method.

Exposure to diverse non-genetic factors, known as the exposome, is a critical determinant of health outcomes. However, analyzing the exposome presents significant methodological challenges, including: high collinearity among exposures, the longitudinal nature of repeated measurements, and potential complex interactions with individual characteristics. In this paper, we address these challenges by proposing a novel statistical framework that extends Bayesian profile regression. Our method integrates profile regression, which handles collinearity by clustering exposures into latent profiles, into a linear mixed model (LMM), a framework for longitudinal data analysis. This profile-LMM approach effectively accounts for within-person variability over time while also incorporating interactions between the latent exposure clusters and individual characteristics. We validate our method using simulated data, demonstrating its ability to accurately identify model parameters and recover the true latent exposure cluster structure. Finally, we apply this approach to a large longitudinal data set from the Lifelines cohort to identify combinations of exposures that are significantly associated with diastolic blood pressure.

We address the open question of counterfactual identification for high-dimensional multivariate outcomes from observational data. Pearl (2000) argues that counterfactuals must be identifiable (i.e., recoverable from the observed data distribution) to justify causal claims. A recent line of work on counterfactual inference shows promising results but lacks identification, undermining the causal validity of its estimates. To address this, we establish a foundation for multivariate counterfactual identification using continuous-time flows, including non-Markovian settings under standard criteria. We characterise the conditions under which flow matching yields a unique, monotone, and rank-preserving counterfactual transport map with tools from dynamic optimal transport, ensuring consistent inference. Building on this, we validate the theory in controlled scenarios with counterfactual ground-truth and demonstrate improvements in axiomatic counterfactual soundness on real images.

By allowing the effects of $p$ covariates in a linear regression model to vary as functions of $R$ additional effect modifiers, varying-coefficient models (VCMs) strike a compelling balance between interpretable-but-rigid parametric models popular in classical statistics and flexible-but-opaque methods popular in machine learning. But in high-dimensional settings where $p$ and/or $R$ exceed the number of observations, existing approaches to fitting VCMs fail to identify which covariates have a non-zero effect and which effect modifiers drive these effects. We propose sparseVCBART, a fully Bayesian model that approximates each coefficient function in a VCM with a regression tree ensemble and encourages sparsity with a global–local shrinkage prior on the regression tree leaf outputs and a hierarchical prior on the splitting probabilities of each tree. We show that the sparseVCBART posterior contracts at a near-minimax optimal rate, automatically adapting to the unknown sparsity structure and smoothness of the true coefficient functions. Compared to existing state-of-the-art methods, sparseVCBART achieved competitive predictive accuracy and substantially narrower and better-calibrated uncertainty intervals, especially for null covariate effects. We use sparseVCBART to investigate how the effects of interpersonal conversations on prejudice could vary according to the political and demographic characteristics of the respondents.

We consider a problem of covariance estimation from a sample of i.i.d. high-dimensional random vectors. To avoid the curse of dimensionality, we impose an additional assumption on the structure of the covariance matrix $\Sigma$. To be more precise, we study the case when $\Sigma$ can be approximated by a sum of double Kronecker products of smaller matrices in a tensor train (TT) format. Our setup naturally extends widely known Kronecker sum and CANDECOMP/PARAFAC models but admits richer interaction across modes. We suggest an iterative polynomial time algorithm based on TT-SVD and higher-order orthogonal iteration (HOOI) adapted to Tucker-2 hybrid structure. We derive non-asymptotic dimension-free bounds on the accuracy of covariance estimation taking into account hidden Kronecker product and tensor train structures. The efficiency of our approach is illustrated with numerical experiments.

Predicting species distributions using occupancy models accounting for imperfect detection is now commonplace in ecology. Recently, modeling spatial and temporal autocorrelation was proposed to alleviate the lack of replication in occupancy data, which often prevents model identifiability. However, how such models perform in highly heterogeneous datasets where missing or single-visit data dominates remains an open question. Motivated by a heterogeneous fine-scale butterfly occupancy dataset, we evaluate the performance of a multi-season occupancy model with spatial and temporal random effects to a skewed (Poisson) distribution of the number of surveys per site, overlap of covariates between occupancy and detection submodels, and spatiotemporal clustering of observations. Results showed that the model is robust to heterogeneous data and covariate overlap. However, when spatiotemporal gaps were added, site occupancy was biased towards the average occupancy, itself overestimated. Random effects did not correct the influence of gaps, due to identifiability issues of variance and autocorrelation parameters. Occupancy analysis of two butterfly species further confirmed these results. Overall, multi-season occupancy models with autocorrelation are robust to heterogeneous data and covariate overlap, but still present identifiability issues and are challenged by severe data gaps, which contaminate predictions even in data-rich areas.

Despite the progress in the development of generative models, their usefulness in creating synthetic data that improve prediction performance of classifiers has been put into question. Besides heuristic principles such as "synthetic data should be close to the real data distribution", it is actually not clear which specific properties affect the generalization error. Our paper addresses this question through the lens of high-dimensional regression. Theoretically, we show that, for linear models, the covariance shift between the target distribution and the distribution of the synthetic data affects the generalization error but, surprisingly, the mean shift does not. Furthermore we prove that, in some settings, matching the covariance of the target distribution is optimal. Remarkably, the theoretical insights from linear models carry over to deep neural networks and generative models. We empirically demonstrate that the covariance matching procedure (matching the covariance of the synthetic data with that of the data coming from the target distribution) performs well against several recent approaches for synthetic data selection, across training paradigms, architectures, datasets and generative models used for augmentation.

We address the problem of estimating a high-dimensional matrix from linear measurements, with a focus on designing optimal rank-adaptive algorithms. These algorithms infer the matrix by estimating its singular values and the corresponding singular vectors up to an effective rank, adaptively determined based on the data. We establish instance-specific lower bounds for the sample complexity of such algorithms, uncovering fundamental trade-offs in selecting the effective rank: balancing the precision of estimating a subset of singular values against the approximation cost incurred for the remaining ones. Our analysis identifies how the optimal effective rank depends on the matrix being estimated, the sample size, and the noise level. We propose an algorithm that combines a Least-Squares estimator with a universal singular value thresholding procedure. We provide finite-sample error bounds for this algorithm and demonstrate that its performance nearly matches the derived fundamental limits. Our results rely on an enhanced analysis of matrix denoising methods based on singular value thresholding. We validate our findings with applications to multivariate regression and linear dynamical system identification.

Tokenization – the process of decomposing a given text into a sequence of subwords called tokens – is one of the key components in the development of language models. Particularly, auto-regressive language models generate texts token by token, i.e., by predicting the next-token distribution given the previous ones, and thus tokenization directly affects their efficiency in text generation. Since each language model has their own vocabulary as a set of possible tokens, they struggle to cooperate with each other at the level of next-token distributions such as model ensemble. In this paper, we establish a theoretical framework of lossless vocabulary reduction, which efficiently converts a given auto-regressive language model into the one with an arbitrarily small vocabulary without any loss in accuracy. This framework allows language models with different tokenization to cooperate with each other efficiently by reduction to their maximal common vocabulary. Specifically, we empirically demonstrate its applicability to model ensemble with different tokenization.

Synthetic data can improve generalization when real data is scarce, but excessive reliance may introduce distributional mismatches that degrade performance. In this paper, we present a learning-theoretic framework to quantify the trade-off between synthetic and real data. Our approach leverages algorithmic stability to derive generalization error bounds, characterizing the optimal synthetic-to-real data ratio that minimizes expected test error as a function of the Wasserstein distance between the real and synthetic distributions. We motivate our framework in the setting of kernel ridge regression with mixed data, offering a detailed analysis that may be of independent interest. Our theory predicts the existence of an optimal ratio, leading to a U-shaped behavior of test error with respect to the proportion of synthetic data. Empirically, we validate this prediction on CIFAR-10 and a clinical brain MRI dataset. Our theory extends to the important scenario of domain adaptation, showing that carefully blending synthetic target data with limited source data can mitigate domain shift and enhance generalization. We conclude with practical guidance for applying our results to both in-domain and out-of-domain scenarios.

Suppose (standardized) measurements or statistics are monitored to raise an alarm when a threshold is exceeded. Often, the underlying population is heterogenous with respect to important discrete variables and thus samples may consist of imbalanced classes. We propose to use thresholds which depend on such covariates to boost the sensitivity for rare classes, which otherwise tend to be ignored. Under mild conditions, we identify optimal threshold functions and develop a feasible procedure for their computation. Further, for the proportional rule a nonparametric estimator of the threshold function is proposed and a central limit theorem is shown, including the case that conditional mean and variance used for standardization are estimated. For feasible uncertainty quantification a bootstrap scheme is proposed. The approach is illustrated and evaluated by a real data analysis.

Normalizing flows with a Gaussian base provide a computationally efficient way to approximate posterior distributions in Bayesian inference, but they often struggle to capture complex posteriors with multimodality and heavy tails. We propose a stick-breaking mixture base with component-wise tail adaptation (StiCTAF) for posterior approximation. The method first learns a flexible mixture base to mitigate the mode-seeking bias of reverse KL divergence through a weighted average of component-wise ELBOs. It then estimates local tail indices of unnormalized densities and finally refines each mixture component using a shared backbone combined with component-specific tail transforms calibrated by the estimated indices. This design enables accurate mode coverage and anisotropic tail modeling while retaining exact density evaluation and stable optimization. Experiments on synthetic posteriors demonstrate improved tail recovery and better coverage of multiple modes compared to benchmark models. We also present a real-data analysis illustrating the practical benefits of our approach for posterior inference.

Let $Z_1, \cdots, Z_n$ denote the eigenvalues of the product $\prod_{j=1}^{k_n} \boldsymbol{A}_j$, where $\{\boldsymbol{A}_j\}_{1 \le j \le k_n}$ are independent $n\times n$ complex Ginibre matrices. Define $\alpha = \lim\limits_{n \to \infty} \frac{n}{k_n}$. We prove that $X_n,$ a suitably rescaled version of $\max_{1 \le j \le n} |Z_j|^2,$ converges weakly as follows: to a non-trivial distribution $\Phi_\alpha$ for $\alpha \in (0, +\infty)$, to the Gumbel distribution when $\alpha = +\infty$, and to the standard normal distribution when $\alpha = 0$. This result reveals a phase transition at the boundaries of $\alpha$. Furthermore, we establish the exact rates of convergence in each regime.

This paper presents four theoretical contributions that improve the usability of risk certificates for neural networks based on PAC-Bayes bounds. First, two bounds on the KL divergence between Bernoulli distributions enable the derivation of the tightest explicit bounds on the true risk of classifiers across different ranges of empirical risk. The paper next focuses on the formalization of an efficient methodology based on implicit differentiation that enables the introduction of the optimization of PAC-Bayesian risk certificates inside the loss/objective function used to fit the network/model. The last contribution is a method to optimize bounds on non-differentiable objectives such as the 0-1 loss. These theoretical contributions are complemented with an empirical evaluation on the MNIST and CIFAR-10 datasets. In fact, this paper presents the first non-vacuous generalization bounds on CIFAR-10 for neural networks. Code to reproduce all experiments is available at github.com/Diegogpcm/pacbayesgradients.

Engineering design involves demanding models encompassing many decision variables and uncontrollable parameters. In addition, unavoidable aleatoric and epistemic uncertainties can be very impactful and add further complexity. The state-of-the-art adopts two steps, uncertainty quantification and design optimization, to optimize systems under uncertainty by means of robust or stochastic metrics. However, conventional scenario-based, surrogate-assisted, and mathematical programming methods are not sufficiently scalable to be affordable and precise in large and complex cases. Here, a multi-level approach is proposed to accurately optimize resource-intensive, high-dimensional, and complex engineering problems under uncertainty with minimal resources. A non-intrusive, fast-scaling, Kriging-based surrogate is developed to map the combined design/parameter domain efficiently. Multiple surrogates are adaptively updated by hierarchical and orthogonal decomposition to leverage the fewer and most uncertainty-informed data. The proposed method is statistically compared to the state-of-the-art via an analytical testbed and is shown to be concurrently faster and more accurate by orders of magnitude.

The Median-of-Means (MoM) is a robust estimator widely used in machine learning that is known to be (minimax) optimal in scenarios where samples are i.i.d. In more grave scenarios, samples are contaminated by an adversary that can inspect and modify the data. Previous work has theoretically shown the suitability of the MoM estimator in certain contaminated settings. However, the (minimax) optimality of MoM and its limitations under adversarial contamination remain unknown beyond the Gaussian case. In this paper, we present upper and lower bounds for the error of MoM under adversarial contamination for multiple classes of distributions. In particular, we show that MoM is (minimax) optimal in the class of distributions with finite variance, as well as in the class of distributions with infinite variance and finite absolute $(1+r)$-th moment. We also provide lower bounds for MoM's error that match the order of the presented upper bounds, and show that MoM is sub-optimal for light-tailed distributions.

We study the problem of estimating a continuous ability parameter from sequential binary responses by actively asking questions with varying difficulties, a setting that arises naturally in adaptive testing and online preference learning. Our goal is to certify that the estimate lies within a desired margin of error, using as few queries as possible. We propose a simple algorithm that adaptively selects questions to maximize Fisher information and updates the estimate using a method-of-moments approach, paired with a novel test statistic to decide when the estimate is accurate enough. We prove that this Fisher-tracking strategy achieves optimal performance in both fixed-confidence and fixed-budget regimes, which are commonly invested in the best-arm identification literature. Our analysis overcomes a key technical challenge in the fixed-budget setting – handling the dependence between the evolving estimate and the query distribution – by exploiting a structural symmetry in the model and combining large deviation tools with Ville's inequality. Our results provide rigorous theoretical support for simple and efficient adaptive testing procedures.

We propose a test for a change in the mean for a sequence of functional observations that are only partially observed on subsets of the domain, with no information available on the complement. The framework accommodates important scenarios, including both abrupt and gradual changes. The significance of the test statistic is assessed via a permutation test. In addition to the classical permutation approach with a fixed number of permutation samples, we also discuss a variant with controlled resampling risk that relies on a random (data-driven) number of permutation samples. The small sample performance of the proposed methodology is illustrated in a Monte Carlo simulation study and an application to real data.

Spatial prediction refers to the estimation of unobserved values from spatially distributed observations. Although recent advances have improved the capacity to model diverse observation types, adoption in practice remains limited in industries that demand interpretability. To mitigate this gap, surrogate models that explain black-box predictors provide a promising path toward interpretable decision making. In this study, we propose a graph partitioning problem to construct spatial segments that minimize the sum of within-segment variances of individual predictions. The assignment of data points to segments can be formulated as a mixed-integer quadratic programming problem. While this formulation potentially enables the identification of exact segments, its computational complexity becomes prohibitive as the number of data points increases. Motivated by this challenge, we develop an approximation scheme that leverages the structural properties of graph partitioning. Experimental results demonstrate the computational efficiency of this approximation in identifying spatial segments.