Papers

We study the computational properties of two canonical planted average-case problems – noisy planted $k$-XOR and Tensor PCA – by formally unifying them into a family of planted problems parametrized by tensor order $k$, number of entries $m$, and noise level $\delta$. We build a wide range of poly-time average-case reductions within this family, across all regimes $m \in [1, n^k]$. In the denser $m \geq n^{k/2}$ regime, our reductions preserve proximity to the computational threshold, and, as a central application, reduce conjectured-hard $k$-XOR instances with $m \approx n^{k/2}$ to conjectured-hard instances of Tensor PCA. Additionally, we give new order-reducing maps at fixed densities (e.g., $5\to 4$ for $k$-XOR with $m \approx n^{k/2}$ entries and $7\to 4$ for Tensor PCA). In the sparser $m \leq n^{k/2}$ regime, we relate instances of different orders, reducing, for example, $7$-XOR with $m = n^{3.4}$ to the classical setting of $3$-XOR with $m = \widetilde\Theta(n^{1.4})$. Taken together, these results establish a hardness partial order in the space of planted tensor models.

With the advancement of machine learning and artificial intelligence technologies, recommender systems have been increasingly used across a vast variety of platforms to efficiently and effectively match users with items. As application contexts become more diverse and complex, there is a growing need for more sophisticated recommendation techniques. One example is the composite item (for example, fashion outfit) recommendation where multiple levels of user preference information might be available and relevant. In this study, we propose JIMA, a joint interaction modeling approach that uses a single model to take advantage of all data from different levels of granularity and incorporate interactions to learn the complex relationships among lower-order (atomic item) and higher-order (composite item) user preferences as well as domain expertise (e.g., on the stylistic fit). We comprehensively evaluate the proposed method and compare it with advanced baselines through multiple simulation studies as well as with real data in both offline and online settings. The results consistently demonstrate the superior performance of the proposed approach.

The Robust Effect Size Index (RESI) is a recently proposed standardized effect size to quantify association strength across models. However, its confidence interval construction has relied on computationally intensive bootstrap procedures. We establish a general theorem for the asymptotic distribution of the RESI using a Taylor expansion that accommodates a broad class of models. Simulations under various linear and logistic regression settings show that RESI and its CI have smaller bias and more reliable coverage than commonly used effect sizes such as Cohen's d and f. Combining with robust covariance estimation yields valid inference under model misspecification. We use the methods to investigate associations of depression and behavioral problems with sex and diagnosis in Autism spectrum disorders, and demonstrate that the asymptotic approach achieves up to a 50-fold speedup over the bootstrap. Our work provides a scalable and reliable alternative to bootstrap inference, greatly enhancing the applicability of RESI to high-dimensional studies.

The Bootstrap Particle Filter (BPF) and the Ensemble Kalman Filter (EnKF) are two widely used methods for sequential Bayesian filtering: the BPF is asymptotically exact but can suffer from weight degeneracy, while the EnKF scales well in high dimension (typically with localization) yet is exact only in the linear-Gaussian case. We combine these approaches by retaining the EnKF transport step and adding a principled importance-sampling correction. Our first contribution is a general importance-sampling theory for mixture targets and proposals, including variance comparisons between individual- and mixture-based estimators. We then interpret the stochastic EnKF analysis as sampling from explicit Gaussian-mixture proposals obtained by conditioning on the current or previous ensemble, which leads to six self-normalized IS-EnKF schemes. We embed these updates into a broader class of ensemble-based filters and prove consistency and error bounds, including weight-variance comparisons and sufficient conditions ensuring finite-variance importance weights. As a second contribution, we construct transported quasi-Monte Carlo (TQMC) point sets for the Gaussian-mixture laws arising in prediction and analysis, yielding TQMC-enhanced variants that can substantially reduce sampling error without changing the filtering pipeline. Numerical experiments on benchmark models compare the proposed mixture-weighted and TQMC-enhanced filters, showing improved filtering accuracy relative to BPF, EnKF, and the standard weighted EnKF, and that the weighted schemes eliminate the EnKF error plateau often caused by analysis-target mismatch.

We propose Kernel Embedding Distributional Reinforcement Learning (KE-DRL) an offline method to estimate kernel mean embeddings of conditional multivariate return distributions from observed offline trajectories with continuous state-action inputs and vector-valued rewards. KE-DRL approximates the conditional embedding and estimates its a finite-dictionary representation coefficients through a maximum mean discrepancy Bellman criterion. For a regular class of return distributions, we show that a Matérn embedding and a Bellman-invariant Sobolev-moment class identifies the unique distributional Bellman fixed point and and derive a Hölder-type bound that controls Wasserstein error by the population embedding Bellman residual. For directly observed responses, we provide finite-sample and uniform error bounds for a regularized conditional mean embedding estimator. Simulation experiments evaluate pointwise embedding recovery against Monte Carlo benchmarks across multiple behavior-target policy pairs. An application to Expedia hotel-search data illustrates conditional multivariate return evaluation under two policies.

Distributed high dimensional mean estimation is a common aggregation routine used often in distributed optimization methods. Most of these applications call for a communication-constrained setting where vectors, whose mean is to be estimated, have to be compressed before sharing. One could independently encode and decode these to achieve compression, but that overlooks the fact that these vectors are often close to each other. To exploit these similarities, recently Suresh et al., 2022, Jhunjhunwala et al., 2021, Jiang et al, 2023, proposed multiple correlation-aware compression schemes. However, in most cases, the correlations have to be known for these schemes to work. Moreover, a theoretical analysis of graceful degradation of these correlation-aware compression schemes with increasing dissimilarity is limited to only the $\ell_2$-error in the literature. In this paper, we propose four different collaborative compression schemes that agnostically exploit the similarities among vectors in a distributed setting. Our schemes are all simple to implement and computationally efficient, while resulting in big savings in communication. The analysis of our proposed schemes show how the $\ell_2$, $\ell_\infty$ and cosine estimation error varies with the degree of similarity among vectors.

Imputing missing node features in graphs is challenging, particularly under high missing rates. Existing methods based on latent representations or global diffusion often fail to produce reliable estimates, and may propagate errors across the graph. We propose FSD-CAP, a two-stage framework designed to improve imputation quality under extreme sparsity. In the first stage, a graph-distance-guided subgraph expansion localizes the diffusion process. A fractional diffusion operator adjusts propagation sharpness based on local structure. In the second stage, imputed features are refined using class-aware propagation, which incorporates pseudo-labels and neighborhood entropy to promote consistency. We evaluated FSD-CAP on multiple datasets. With $99.5\%$ of features missing across five benchmark datasets, FSD-CAP achieves average accuracies of $80.06\%$ (structural) and $81.01\%$ (uniform) in node classification, close to the $81.31\%$ achieved by a standard GCN with full features. For link prediction under the same setting, it reaches AUC scores of $91.65\%$ (structural) and $92.41\%$ (uniform), compared to $95.06\%$ for the fully observed case. Furthermore, FSD-CAP demonstrates superior performance on both large-scale and heterophily datasets when compared to other models.

Popularized by their strong image generation performance, diffusion and related methods for generative modeling have found widespread success in visual media applications. In particular, diffusion methods have enabled new approaches to data compression, where realistic reconstructions can be generated at extremely low bit-rates. This article provides a unifying review of recent diffusion-based methods for generative lossy compression, with a focus on image compression. These methods generally encode the source into an embedding and use a diffusion model to iteratively refine it during decoding, so that the reconstruction approximately follows the true data distribution. The embedding can take various forms and is typically transmitted via an auxiliary entropy model, and recent methods also explore the use of diffusion models themselves for information transmission via channel simulation. We review representative approaches through the lens of rate-distortion-perception theory, highlighting the role of common randomness and connections to inverse problems, and identify open challenges.

Q-learning and SARSA are foundational reinforcement learning algorithms whose practical success depends critically on step-size calibration. Step-sizes that are too large can cause numerical instability, while step-sizes that are too small can lead to slow progress. We propose implicit variants of Q-learning and SARSA that reformulate their iterative updates as fixed-point equations. This yields an adaptive step-size adjustment that scales inversely with feature norms, providing automatic regularization without manual tuning. Our non-asymptotic analyses demonstrate that implicit methods maintain stability over significantly broader step-size ranges. Under favorable conditions, it permits arbitrarily large step-sizes while achieving comparable convergence rates. Empirical validation across benchmark environments spanning discrete and continuous state spaces shows that implicit Q-learning and SARSA exhibit substantially reduced sensitivity to step-size selection, achieving stable performance with step-sizes that would cause standard methods to fail.

As generative models continue to evolve, detecting AI-generated images remains a critical challenge. While effective detection methods exist, they often lack formal interpretability and may rely on implicit assumptions about fake content, potentially limiting robustness to distributional shifts. In this work, we introduce a rigorous, statistically grounded framework for fake image detection that focuses on producing a probability score interpretable with respect to the real-image population. Our method leverages the strengths of multiple existing detectors by combining training-free statistics. We compute p-values over a range of test statistics and aggregate them using classical statistical ensembling to assess alignment with the unified real-image distribution. This framework is generic, flexible, and training-free, making it well-suited for robust fake image detection across diverse and evolving settings.

Differential item functioning (DIF) or measurement invariance (MI) testing for single-item assessments has previously been impossible. Part of the issue is that there are no conditioning variables to serve as a proxy for the latent variable–regression-based DIF methods. Another reason is that factor-analytic approaches require multiple items to estimate parameters. In this technical working paper, I propose an approach for evaluating DIF/MI in a single-item assessment of a construct. The current methods should NOT replace using multiple-indicator MG-CFA/IRT analyses of DIF/MI or regression mased methods when possible. More items generally provide significantly better construct coverage and provide more rigorous DIF/MI evaluation.

Explainable boosting machines (EBMs) are popular "glass-box" models that learn a set of univariate functions using boosting trees. These achieve explainability through visualizations of each feature's effect. However, unlike linear model coefficients, uncertainty quantification for the learned univariate functions requires computationally intensive bootstrapping, making it hard to know which features truly matter. We provide an alternative using recent advances in statistical inference for gradient boosting, deriving methods for statistical inference as well as end-to-end theoretical guarantees. Using a moving average instead of a sum of trees (Boulevard regularization) allows the boosting process to converge to a feature-wise kernel ridge regression. This produces asymptotically normal predictions that achieve the minimax-optimal MSE for fitting Lipschitz GAMs with $p$ features of $O(p n^{-2/3})$, successfully avoiding the curse of dimensionality. We then construct prediction intervals for the response and confidence intervals for each learned univariate function with a runtime independent of the number of datapoints, enabling further explainability within EBMs. Code is available at https://github.com/hetankevin/ebm-inference.

Unsupervised text segmentation is crucial because boundary labels are expensive, subjective, and often fail to transfer across domains and granularity choices. We propose Embed-KCPD, a training-free method that represents sentences as embedding vectors and estimates boundaries by minimizing a penalized KCPD objective. Beyond the algorithmic instantiation, we develop, to our knowledge, the first dependence-aware theory for KCPD under $m$-dependent sequences, a finite-memory abstraction of short-range dependence common in language. We prove an oracle inequality for the population penalized risk and a localization guarantee showing that each true change point is recovered within a window that is small relative to segment length. To connect theory to practice, we introduce an LLM-based simulation framework that generates synthetic documents with controlled finite-memory dependence and known boundaries, validating the predicted scaling behavior. Across standard segmentation benchmarks, Embed-KCPD often outperforms strong unsupervised baselines. A case study on Taylor Swift's tweets illustrates that Embed-KCPD combines strong theoretical guarantees, simulated reliability, and practical effectiveness for text segmentation.

Evaluating the quality of search, ranking and RAG systems traditionally requires a significant number of human relevance annotations. In recent times, several deployed systems have explored the usage of Large Language Models (LLMs) as automated judges for this task while their inherent biases prevent direct use for metric estimation. We present a statistical framework extending Prediction-Powered Inference (PPI) that combines minimal human annotations with LLM judgments to produce reliable estimates of metrics which require sub-instance annotations. Our method requires as few as 100 human-annotated queries and 10,000 unlabeled examples, reducing annotation requirements significantly compared to traditional approaches. We formulate our proposed framework (PRECISE) for inference of relevance uplift for an LLM-based query reformulation application, extending PPI to sub-instance annotations at the query-document level. By reformulating the metric-integration space, we reduced the computational complexity from O(2^|C|) to O(2^K), where |C| represents corpus size (in order of millions). Detailed experiments across prominent retrieval datasets demonstrate that our method reduces the variance of estimates for the business-critical Precision@K metric, while effectively correcting for LLM bias in low-resource settings.

Live win probabilities have given sports collapses a common numerical language: the highest win probability attained by the team that eventually lost. Yet this familiar number is a selected pathwise extreme. The losing path is identified by the terminal outcome and then searched for its most favorable point. Under ideal sequential calibration, we derive an exact continuous-path benchmark for this statistic, a conservative bound for discretely reported paths, and a probability-integral-transform diagnostic for collections of games. Unlike fixed-time calibration summaries and proper scores, the diagnostic asks whether a published feed produces severe losing paths as often as a calibrated sequential model should. We apply the benchmark separately to public regular-season NFL and NBA feeds from 2018-2024, using a season-stratified dyadic bootstrap to account for recurring teams. We detect no global departure in the NFL and a systematic excess of extreme losing-team peaks in the NBA. In the upper 5% benchmark tail, the NBA excess is 1.3% (95% interval 0.4% to 2.3%); first crossings of a published win probability of 0.95 show the same positive gap between observed and implied loss rates.

This paper investigates when two general ridge estimators coincide under spatial error models. First, in the general linear model, we derive a necessary and sufficient condition based on the commutativity of an extended dispersion matrix and an orthogonal projector, thereby extending the classical result of Zyskind. Next, we establish parameter-free conditions for first-order spatial autoregressive and spatial moving average processes. Also, we obtain a necessary and sufficient condition valid for all values of the spatial correlation coefficient in the specified parameter range and characterize all penalty matrices for which the two general ridge estimators coincide. A numerical experiment illustrates the theoretical results and shows that the proposed conditions can simplify the two-step estimation procedure by eliminating the need to estimate the spatial correlation coefficient.

Modern generative modeling methods have demonstrated strong performance in learning complex data distributions from clean samples. In many scientific and imaging applications, however, clean samples are unavailable, and only noisy or linearly corrupted measurements can be observed. Moreover, latent structures, such as manifold geometries, present in the data are important to extract for further downstream scientific analysis. In this work, we introduce Riemannian AmbientFlow, a framework for simultaneously learning a probabilistic generative model and the underlying, nonlinear data manifold directly from corrupted observations. Building on the variational inference framework of AmbientFlow, our approach incorporates data-driven Riemannian geometry induced by normalizing flows, enabling the extraction of manifold structure through pullback metrics and Riemannian Autoencoders. We establish theoretical guarantees showing that, under appropriate geometric regularization and measurement conditions, the learned model recovers the underlying data distribution up to a controllable error and yields a smooth, bi-Lipschitz manifold parametrization. We further show that the resulting smooth decoder can serve as a principled generative prior for inverse problems with recovery guarantees. We empirically validate our approach on low-dimensional synthetic manifolds and on MNIST.

We study function estimation in the empirical Bayes setting for Poisson and normal means. Specifically, given observations $X_i\sim f(\cdot; \theta_i)$ with latent parameters $\theta_i\sim \pi$, the goal is to estimate $\mathbb{E}_{\pi}[\ell(\theta)|X = x]$. This task lies between classical deconvolution (recovering the full prior $\pi$), and standard empirical Bayes mean estimation. While the minimax risk for estimating $\pi$ in the Wasserstein distance is known to decay only logarithmically, we show that estimating the corresponding posterior smooth functionals admits dramatically faster rates. In particular, for polynomial functions of degree $k$ in the Poisson model, we establish a tight total regret bound of $\Theta((\frac{\log n}{\log \log n})^{k+1})$ and $\Theta((\log n)^{2k+1})$ for bounded and subexponential priors, respectively, attainable by estimators mimicking those that achieve optimal regret for the mean estimation problem (Robbins, minimum distance, ERM). In the normal means model, we establish tight total regret bound of $\Theta((\frac{\log n}{\log \log n})^{k+1})$ for bounded priors, and bounds that match up to a polylogarithmic factor for subgaussian priors. Our analysis identifies the approximation-theoretic origin of this improvement: smooth functions can be well-approximated by low-degree polynomials, whereas Lipschitz functions have only $O(\frac{1}{k})$ degree-$k$ polynomial approximation error. The results reveal a sharp hierarchy in the difficulty of empirical Bayes problems: ranging from slow, logarithmic deconvolution to near-parametric convergence for smooth posterior functionals, and establish new connections between nonparametric empirical Bayes theory, polynomial approximation, and statistical inverse problems.

The marginal likelihood, or Bayesian evidence, is a crucial quantity for Bayesian model comparison but its computation can be challenging for complex models, even in parameters space of moderate dimension. The learned harmonic mean estimator has been shown to provide accurate and robust estimates of the marginal likelihood simply using posterior samples. It is agnostic to the sampling strategy, meaning that the samples can be obtained using any method. This enables marginal likelihood calculation and model comparison with whatever sampling is most suitable for the task. However, the internal density estimators considered previously for the learned harmonic mean can struggle with highly multimodal posteriors. In this work we introduce flow matching-based continuous normalizing flows as a powerful architecture for the internal density estimation of the learned harmonic mean. We demonstrate the ability to handle challenging multimodal posteriors, including an example in 20 parameter dimensions, showcasing the method's ability to handle complex posteriors without the need for fine-tuning or heuristic modifications to the base distribution.

We investigate the detection of weak complex-valued signals immersed in non-Gaussian, range-varying interference, with emphasis on maritime radar scenarios. The proposed methodology exploits a Complex-valued Variational AutoEncoder (CVAE) trained exclusively on clutter-plus-noise to perform Out-Of-Distribution detection. By operating directly on in-phase / quadrature samples, the CVAE preserves phase and Doppler structure and is assessed in two configurations: (i) using unprocessed range profiles and (ii) after local whitening, where per-range covariance estimates are obtained from neighboring profiles. Using extensive simulations together with real sea-clutter data from the CSIR maritime dataset, we benchmark performance against classical and adaptive detectors (MF, NMF, AMF-SCM, ANMF-SCM, ANMF-Tyler). In both configurations, the CVAE yields a higher detection probability Pd at matched false-alarm rate Pfa, with the most notable improvements observed under whitening. We further integrate the CVAE with the ANMF through a weighted log-p fusion rule at the decision level, attaining enhanced robustness in strongly non-Gaussian clutter and enabling empirically calibrated Pfa control under H0. Overall, the results demonstrate that statistical normalization combined with complex-valued generative modeling substantively improves detection in realistic sea-clutter conditions, and that the fused CVAE-ANMF scheme constitutes a competitive alternative to established model-based detectors.

When developing clinical prediction models, it can be challenging to balance between global models that are valid for all patients and personalized models tailored to individuals or potentially unknown subgroups. To aid such decisions, we propose a diagnostic tool for contrasting global regression models and patient-specific (local) regression models. The core utility of this tool is to identify where and for whom a global model may be inadequate. We focus on regression models and specifically suggest a localized regression approach that identifies regions in the predictor space where patients are not well represented by the global model. As localization becomes challenging when dealing with many predictors, we propose modeling in a dimension-reduced latent representation obtained from an autoencoder. Using such a neural network architecture for dimension reduction enables learning a latent representation simultaneously optimized for both good data reconstruction and for revealing local outcome-related associations suitable for robust localized regression. We illustrate the proposed approach with a clinical study involving patients with chronic obstructive pulmonary disease. Our findings indicate that the global model is adequate for most patients but that indeed specific subgroups benefit from personalized models. We also demonstrate how to map these subgroup models back to the original predictors, providing insight into why the global model falls short for these groups. Thus, the principal application and diagnostic yield of our tool is the identification and characterization of patients or subgroups whose outcome associations deviate from the global model.

Previous work revealed associations between flood exposure and adverse health outcomes during and in the aftermath of flood events. Floods are highly heterogeneous events, largely owing to vast differences in flood durations, i.e., flash-floods versus slow-moving floods. However, little to no work has incorporated exposure duration into the modeling of flood-related health impacts or has investigated duration-driven effect heterogeneity. To address this gap, we propose an exposure duration varying coefficient modeling (EDVCM) framework for estimating exposure day-specific health effects of consecutive-day environmental exposures that vary in duration. We develop the EDVCM within an area-level self-matched study design to eliminate time-invariant confounding followed by conditional Poisson regression modeling for exposure effect estimation and adjustment of time-varying confounders. Using a Bayesian framework, we introduce duration- and exposure day-specific exposure coefficients within the conditional Poisson model and assign them a two-dimensional Gaussian process prior to allow for sharing of information across both duration and exposure day. This approach enables highly-resolved insights into duration-driven effect heterogeneity while ensuring model stability through information sharing. Through simulations, we demonstrate that the EDVCM out-performs conventional approaches in terms of both effect estimation and uncertainty quantification. We apply the EDVCM to nationwide, multi-decade Medicare claims data linked with high-resolution flood exposure measures to investigate duration-driven heterogeneity in flood effects on musculoskeletal system disease hospitalizations.

Generating quantum data by learning the underlying quantum distribution poses challenges in both theoretical and practical scenarios, yet it is a critical task for understanding quantum systems. A fundamental question in quantum machine learning (QML) is the universality of approximation: whether a parameterized QML model can approximate any quantum distribution. We address this question by proving a universality theorem for the Many-body Projected Ensemble (MPE) framework, a method for quantum state design that uses a single many-body wave function to prepare random states. This demonstrates that MPE can approximate any distribution of pure states within a 1-Wasserstein distance error. This theorem provides a rigorous guarantee of universal expressivity, addressing key theoretical gaps in QML. For practicality, we propose an Incremental MPE variant with layer-wise training to improve the trainability. Numerical experiments on clustered quantum states and quantum chemistry datasets validate MPE's efficacy in learning complex quantum data distributions.

Natural gradients have long been studied in deep reinforcement learning due to their fast convergence properties and covariant weight updates. However, computing natural gradients requires inversion of the Fisher Information Matrix (FIM) at each iteration, which is computationally prohibitive in nature. In this paper, we present an efficient and scalable natural policy optimization technique that leverages a rank-1 approximation to full inverse-FIM. We theoretically show that under certain conditions, a rank-1 approximation to inverse-FIM converges faster than policy gradients and, under some conditions, enjoys the same sample complexity as stochastic policy gradient methods. We benchmark our method on a diverse set of environments and show that it achieves superior performance to standard actor-critic and trust-region baselines.

Joint models for longitudinal and time-to-event data are widely used in many disciplines. Nonetheless, existing model comparison criteria do not indicate whether a model adequately fits the data or which components may be misspecified. We introduce a Bayesian posterior predictive checks framework for assessing a joint model's fit to the longitudinal and survival processes and their association. The framework supports multiple settings, including existing subjects, new subjects with only covariates, dynamic prediction at intermediate follow-up times, and cross-validated assessment. For the longitudinal component, goodness-of-fit is assessed through the mean, variance, and correlation structure, while the survival component is evaluated using empirical cumulative distributions and probability integral transforms. The association between processes is examined using time-dependent concordance statistics. We apply these checks to the Bio-SHiFT heart failure study, and a simulation study demonstrates that they can identify model misspecification that standard information criteria fail to detect. The proposed methodology is implemented in the freely available R package JMbayes2.

The generalised balanced power diagram (GBPD) is regarded in the literature as a suitable geometric model for describing polycrystalline microstructures with curved grain boundaries. This article compiles properties of GBPDs with regard to affine transformations and flat sections. Furthermore, it extends an algorithm known for power diagrams for generating digital images, which is more efficient than the usual brute force approach, on GBPDs.

We review vaccine efficacy (VE) estimands for susceptibility in individual randomized trials with natural (unmeasured) exposure, where individual responses are measured as time from vaccination until an event (e.g., disease from the infectious agent). Common VE estimands are written as $1-\theta$, where $\theta$ is some ratio effect measure (e.g., ratio of incidence rates, cumulative incidences, hazards, or odds) comparing outcomes under vaccination versus control. Although the ratio effects are approximately equal with low control event rates, we explore the quality of that approximation using a nonparametric formulation. Traditionally, the primary endpoint VE estimands are full immunization (or biological) estimands that represent a subset of the intent-to-treat population, excluding those that have the event before the vaccine has been able to ramp-up to its full effect, requiring care for proper causal interpretation. Besides these primary VE estimands that summarize an effect of the vaccine over the full course of the study, we also consider local VE estimands that measure the effect at particular time points. We discuss interpretational difficulties of local VE estimands (e.g., depletion of susceptibles bias), and using frailty models as sensitivity analyses for the individual-level causal effects over time.

Predicting fuel assembly bow in pressurized water reactors requires solving tightly coupled fluid-structure interaction problems, whose direct simulations can be computationally prohibitive, making large-scale uncertainty quantification (UQ) very challenging. This work introduces a general mathematical framework for coupling Gaussian process (GP) surrogate models representing distinct physical solvers, aimed at enabling rigorous UQ in coupled multiphysics systems. A theoretical analysis establishes that the predictive variance of the coupled GP system remains bounded under mild regularity and stability assumptions, ensuring that uncertainty does not grow uncontrollably through the iterative coupling process. The methodology is then applied to the coupled hydraulic-structural simulation of fuel assembly bow, enabling global sensitivity analysis and full UQ at a fraction of the computational cost of direct code coupling. The results demonstrate accurate uncertainty propagation and stable predictions, establishing a solid mathematical basis for surrogate-based coupling in large-scale multiphysics simulations.

Ranking or assessing centrality in multivariate and non-Euclidean data is difficult because there is no canonical order and many depth notions become computationally fragile in high-dimensional or structured settings. We introduce a preference-based notion of centrality defined through population proximity comparisons with respect to a random reference draw, yielding a metric-intrinsic statistical functional that is well-defined on general metric spaces. Because the induced pairwise preferences may be non-transitive, we map them to a coherent one-dimensional score via a Bradley–Terry–Luce cross-entropy projection, viewed as a calibrated aggregation device rather than a correctly specified model. We develop two finite-sample estimators a convex M-estimator and a fast spectral estimator based on a comparison operator, and establish identifiability and consistency under mild conditions. Simulations and real-data examples, including high-dimensional and functional observations, illustrate that the proposed scores provide stable, interpretable rankings aligned with the underlying preference centrality.

We show that the percentile-percentile (P-P) process constructed from an independent and identically distributed sample of pairs converges in distribution in $L^1[0,1]$ if and only if the associated P-P curve is absolutely continuous. When this condition holds, the limiting distribution is Gaussian and the process admits a valid bootstrap approximation.