Papers

Recent works in econometrics have shown that there can be issues with applying a constant treatment effect model in longitudinal settings with staggered treatment and heterogeneous treatment effects. We focus on the issue that the estimated constant treatment effect may be a weighted average, with some negative weights, of treatment effects that are heterogeneous across treatment durations. When this issue arises, the estimated constant treatment effect and estimated heterogeneous treatment effects may result in conflicting results. Through the example of estimating the effect of drug-induced homicide (DIH) prosecutions reported by media on unintentional drug-overdose deaths in the United States, we illustrate how the negative weighting issue can lead to conflicting results in practice. Moreover, although research has shown that the negative weight issue may arise in linear regression models, we show this issue may also arise in logistic regression models. Using a linear link, we estimated a constant treatment effect risk ratio of 0.977 (95% CI:(0.866, 1.101)) and an average risk ratio of 0.728 (range: 0.507-0.979) over different treatment durations. Using a logistic link, we estimated a constant treatment risk ratio effect of 1.064 (95% CI: (0.972, 1.165)) and an average risk ratio of 0.739 (range: 0.538-1.008) over different treatment durations. Under both models, the estimated constant treatment effect is either smaller in magnitude or has a different sign than almost all estimated heterogeneous treatment effects, suggesting a negative weighting issue is present. Our results suggest additional care is needed when applying constant treatment effect models in longitudinal settings.

We propose a methodology based on the standard ReLU Deep Neural Networks (DNN) to make predictions and quantify their uncertainty. Classically, people rely on linear, non-linear, or non-parametric kernel methods to fit and then predict the time series. As the universal approximation ability was revealed for DNN, its application has become more and more popular for prediction tasks in various scientific areas. However, the corresponding uncertainty quantification has not been studied thoroughly. Particularly, the uncertainty in prediction will consist of two parts: (1) the future variability; (2) the estimation variability within training data. To capture both variabilities, we build the so-called pertinent prediction interval (PPI) with the DNN model estimator. We first explore the consistency property of the DNN estimator with beta-mixing dependent data. Subsequently, we show that the implied forward bootstrap series is still beta-mixing and possesses the same stationary distribution as the original time series in probability, which is a key condition to enable the PPI. Lastly, the desired PPI is built after imposing minimal conditions on the limiting distribution of predictive roots. Simulations and real-data analysis are deployed to challenge our approach with standard non-parametric methods.

Nominal LoRA rank is a design parameter; calibrated spectral evidence is a separate inferential quantity. This article develops a finite-sample framework for inferring effective rank structure in public foundation-model adapters. The theoretical core is an exact chi-square divergence for the fixed-dimensional Gaussian rank-one reference experiment, with an unknown signal direction integrated under a rotation-invariant reference prior. The resulting series yields a computable finite-sample Le Cam bound at concrete layer sizes, an explicit remainder bound for numerical truncation, and the rectangular Baik-Ben Arous-Peche (BBP) limit. A compact-manifold Laplace expansion shows that finite-sample likelihood evidence also depends on leading spectral gaps through the factor $s_1^{|m-n|}\prod_{i\ge2}(s_1^2-s_i^2)$, motivating joint calibration of clustered singular values. Building on these results, we introduce an empirical-null workflow for PEFT LoRA adapters: factor reconstruction, Monte Carlo $p$-values, stagewise and block testing, and module-wise and corpus-level BH reporting. In an audit of 26 public adapters, 684 modules, six architecture families, and 31,770 public-checkpoint spectra rows, calibrated effective rank is typically much smaller than nominal rank and differs systematically from 95% energy retention. A measured RoBERTa-RTE slice on $n=24$ examples illustrates the measurement path from calibrated ranks to task evaluation, without treating the slice as a utility study. The main empirical finding is that calibrated effective rank is usually far below nominal rank, and that energy retention and statistical surprise answer different questions.

We propose a family of control variate estimators for variance reduction in design-based survey sampling and causal inference, with and without interference. In these settings, inverse probability weighting (IPW) estimators are widely used, but may have large variance when sampling, treatment, or exposure probabilities are small. Building on the observation that several common estimators, including the Hajek estimator, the normalized estimator, the augmented inverse probability weighting (AIPW) estimator, and the targeted maximum likelihood estimator (TMLE), all correct the Horvitz-Thompson estimator by canceling part of its randomness, we provide a unified interpretation of these estimators as special cases of a general control variate estimator. We then construct optimal control variates that can reduce the finite sample variance compared to these common estimators. We parameterize the proposed control variates by their bases and characterize the optimal bases through a stochastic optimization formulation. In survey sampling and causal inference without interference, the optimal bases are characterized by leading eigenvectors of matrices that depend on both the design-based sampling structure and the model-based outcome uncertainty. In causal inference under network interference, the optimal bases solve a nonconvex quadratic optimization problem; we provide a $\frac{1}{2}$-approximate solution and an alternating local search heuristic. We apply the control variate estimators to the Swiss Environmental Panel survey data and the Chinese social network data, and conduct extensive simulations to show that the proposed control variate estimators can achieve substantial variance reduction.

Constraint-based causal discovery like PC and FCI depends on its conditional independence test. Partial correlation and the Generalised Covariance Measure (GCM) detect only the conditional covariance of residuals, so they miss dependence in the mean's nonlinear part, the scale, and the tails. Tests that detect more are biased inside PC, not scalable, only continuous, or not aimed at the tails. Our Generalised Feature Covariance Measure (GFCM) is valid, sensitive beyond covariance, robust inside PC, and applicable to mixed-type data. It runs the GCM template on a configurable set of residual features with conditional mean zero (centered moments and conditional quantile indicators), pooled in blocks and combined by the Cauchy rule, with a growing-knot spline nuisance at regression cost. We contribute (i) a centering result making the scale feature Neyman orthogonal, where the uncentered version is biased; (ii) the orientation asymmetry the mean-quantile construction creates inside PC, and its fix; (iii) a Phi-faithfulness theory under which PC with GFCM recovers the CPDAG of the set's detection class; and (iv) a benchmark of CI tests sensitive beyond covariance on synthetic data, semi-synthetic tail injections, and PC discovery on random DAGs. Under size-corrected power, GFCM recovers the scale and tail edges the covariance family misses and alone keeps power at the deep conditioning sets PC issues. It stays calibrated as n grows, whereas FFCI, the boosted GCM, and the partial copula test do not, and it handles mixed-type data directly. Inside PC at scale it attains the lowest skeleton SHD among tests that stay calibrated, while the others inflate false edges. Validity rests on an additive nuisance, and the tail advantage is shown on simulated and semi-synthetic data, as no fully real benchmark with both heavy tails and known structure exists.

Decorrelating a multivariate system need not destroy source-variable identity. We introduce Constrained Oblique Rotation with Anchored Loadings (CORAL), which minimizes residual cross-correlation while guaranteeing a declared minimum correlation between each transformed variable and its designated source. For a p-variable correlation matrix $R$, we show that every exact decorrelator can be written as $R^{-1/2}Q$ for some orthogonal matrix $Q$, and define $\rho_\star(R)$ as the maximum common source fidelity compatible with exact decorrelation. Constructive lower bounds and rigorous analytical upper bounds tightly bracket $\rho_\star$ at [0.972,0.976], [0.959,0.962], and [0.949,0.950] in simulations with p={6,18,50}, respectively, compared with PCA's largest achievable minimum correlation between distinct principal components and matched source variables of 0.358, 0.329, and 0.172. Corresponding intervals are [0.751,0.758] for World Development Indicators and [0.826,0.835] for wine chemistry data sets. Thus, loss of source-variable identity is not inherent to exact decorrelation but depends on the decorrelator selected. CORAL uses constrained Riemannian optimization and extends to exact support restrictions.

In this paper, we propose SHOPCA (Shape Operator-based Principal Component Analysis), a novel method for unsupervised metric learning and dimensionality reduction that incorporates differential geometric information into the covariance structure of classical PCA. SHOPCA regularizes the global covariance matrix using the mean shape operator, defined as the average of the absolute local shape operators estimated from the data manifold, steering principal components toward directions of both maximum variance and informative curvature. A single trace-normalized mixing coefficient $\alpha$ controls the regularization, recovering standard PCA at $\alpha = 0$ and a curvature-driven embedding as $\alpha \to \infty$. We further introduce a fully unsupervised criterion for selecting $\alpha$ based on the spectral eigengap of the regularized covariance matrix, maximizing the relative separation between the top-$d$ and remaining eigenvalues without using class labels. We evaluate SHOPCA on more than 50 real-world benchmark datasets, comparing it with PCA, ISOMAP, and UMAP using Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), Fowlkes-Mallows index (FM), and V-measure. Results show that SHOPCA consistently improves clustering quality over PCA across a broad range of datasets and surpasses UMAP on small-sample settings, where iterative neighborhood-based manifold estimation can degrade. SHOPCA is computationally tractable, parameter-efficient, and applicable to domains requiring fully unsupervised, geometry-aware dimensionality reduction.

High-throughput single-cell and spatial transcriptomic technologies provide high-resolution snapshots of heterogeneous cellular states, but their destructive nature prevents repeated measurements of the same cells over time. Consequently, temporal and spatial dynamics must be inferred from independently sampled, unaligned cell populations, making it challenging to reconstruct developmental trajectories. Optimal transport (OT) offers a geometric framework for aligning cell populations and inferring developmental trajectories, but many existing approaches focus on modeling the evolution of distributions of cells in gene expression space rather than the relational structure encoded by gene expression networks. To address this limitation, we introduce a geometric framework for analyzing the spatiotemporal evolution of gene expression networks through embeddings in Gromov–Wasserstein (GW) space. By representing each developmental stage as a graph combining gene expression and spatial proximity, our approach enables comparisons of network structure across time, continuous interpolation between developmental stages via GW geodesics, and quantification of network-level changes using Ollivier-Ricci curvature. We evaluate our framework on a spatiotemporal transcriptomic Drosophila dataset and show that GW geodesic interpolations reproduce main trends in curvature dynamics observed in empirical gene expression networks. Agreement with higher-order Co-Optimal Transport (COOT) distances, which jointly represent spatial and temporal information, further validates the framework and suggests that hypernetwork representations successfully record salient biological changes across time. In general, our approach provides a unified geometric approach to study dynamically evolving biological networks.

The increasing availability of large and complex datasets across many scientific disciplines has led to widespread adoption of machine learning (ML) for prediction. However, most ML algorithms focus on point estimation and provide limited information about predictive uncertainty or the conditional distribution of the response, restricting their ability to characterize rare or extreme outcomes. We develop QXGB, a quantile-based gradient boosting framework, and introduce a convolution smoothed loss within it that estimates conditional quantiles for constructing dense cumulative distribution functions (CDFs), exceedance probabilities, and tail behaviour relevant to extreme outcomes. This approach preserves the computational efficiency of extreme gradient boosting while restoring the Hessian information XGBoost relies on for tree splitting, in turn providing interpretable measures of extreme value and exceedance probability predictions. We derive the gradients and Hessians needed to integrate convolution smoothed quantile loss with different kernel specifications into XGBoost, and with simulated data, benchmark this approach against alternative smoothed quantile regression losses, the native quantile objective in the XGBoost Python package, and independent versus multi-output tree estimation. The practical relevance is illustrated in an application predicting fine particulate matter (PM$_{2.5}$) in northern California, including periods where levels were elevated due to wildfire smoke. Our results show that convolution smoothed QXGB, particularly when paired with multi-output trees, delivers accurate predictions with near-zero quantile crossing, well-calibrated CDF and exceedance probability estimates, and useful tail characterization for extreme values. Interval estimation is also evaluated as a measure of data spread.

Power-law anisotropy in internal representations has been observed across a wide range of biological and artificial neural systems, from state-of-the-art language models to the mouse cerebral cortex. This anisotropy is a key geometric property of high-dimensional information processing and underlies a variety of theoretical analyses. However, the mechanism by which it emerges from input structure and task-driven learning has remained unclear. Here, we characterize this formation process by exactly solving the learning dynamics of a wide two-layer linear neural network in a teacher–student setting with power-law input and teacher structures. We show that, in the feature-learning regime, the local power-law exponent of the internal-representation spectrum evolves nonmonotonically over the course of training and exhibits up to four distinct asymptotic regimes across modes and training times. By contrast, in the lazy regime, the exponent remains essentially unchanged. We further demonstrate numerically that similar exponent dynamics arise in more realistic nonlinear networks. Together, these results suggest a general mechanism by which the dynamic interaction between input statistics and task structure gives rise to power-law internal representations.

Token-level knowledge distillation (KD) matches two conditional distributions per position, yet the standard objectives compare them pointwise: a Kullback-Leibler gradient is blind to which wrong token receives probability mass. We develop a distributional view in which the teacher is represented not by a single softened output but by a family of multi-temperature views - marginals of the annealing path of its logits - and the student is trained against a geometry-aware aggregate of these views under an embedding-based ground cost. We formalize the resulting design space (mixtures, log-linear pooling, entropic Wasserstein barycenters, and a debiased Sinkhorn-divergence flagship in hub and path forms), prove an exact collapse result showing log-linear pooling of tempered views is equivalent to a single temperature, and give a multi-marginal Schrodinger-bridge reading that yields falsifiable predictions. On instruction-tuned Pythia pairs, experiments yield three empirical laws: (i) dispersion law - the benefit of multi-temperature aggregation grows monotonically with the effective temperature dispersion of the views, not with their number; (ii) dispersed views unlock the aggregation operator - the barycenter separates from the arithmetic mixture exactly when transport-based aggregation starts to beat averaging; and (iii) two-regime picture governed by the ceiling gap $\Gamma=\mathrm{PPL}_{\mathrm{SFT}}-\mathrm{PPL}_{T}$: when the fine-tuned teacher barely beats a supervised student the gentle transport objective is the best KD loss but no KD beats supervised fine-tuning, whereas at a real ceiling the ranking inverts - and the sign of the fidelity-generalization correlation flips. We argue that "which distillation loss is the best" is not a fixed property of the loss but a function of $\Gamma$.

Parameter identification of stochastic dynamical systems driven by mixed noises is challenging due to intractable likelihood functions. We propose PENN-GMD, a parameter estimation neural network that maps partially observed trajectories to a Gaussian mixture distribution (GMD) over the system parameters. Unlike conventional uncertainty estimates, the GMD employs full covariance matrices to explicitly reveal parameter couplings and multi-modal likelihood structures. The network is trained by minimizing the negative log-likelihood via a surjective parameterization that hard-encodes all GMD constraints, thereby approximating the true likelihood. We validate the method on five numerical examples with increasing complexity, including systems driven by fractional Gaussian and Lévy noises, oscillators with colored noise, coupled neurons under different observability, and an aeroelastic airfoil with unidentifiable stochastic disturbances. Results demonstrate that PENN-GMD accurately recovers likelihood distributions, captures parameter couplings, and naturally diagnoses non-identifiability through variance broadening or mode splitting. These capabilities establish PENN-GMD as a practical tool for uncertainty-aware parameter identification in complex stochastic systems where conventional likelihood-based methods are infeasible.

In this paper we show how to extend quantile-stratified sampling to produce simulations from various multivariate distributions. These simulations have desirable space-filling and coverage properties relative to simulation using IID sampling. We examine the coverage performance of these simulations against IID sampling by looking at plots of ordered log-density values from the simulations.

Population-based adaptive importance sampling (AIS) methods use a set of proposal densities to approximate complex target distributions. Their performance is commonly assessed through effective sample size (ESS) and related weight-based diagnostics, which measure the concentration of normalized importance weights. However, a large ESS only indicates that the normalized sample weights are not strongly concentrated; it does not describe how the proposal components are arranged in the sampling space. In population-based AIS, several proposal components may generate samples in the same region of the target, so the sample weights can appear well balanced even though the effective number of distinct proposal components is small. This letter introduces the effective number of proposals (ENP), a similarity-aware proposal-level diagnostic for population-based AIS. ENP combines the total normalized weight assigned to each proposal with a redundancy measure computed from similarities among target-weighted samples, estimating the number of non-redundant empirical proposal contributions to the approximation. We establish basic effective-number properties and show that ENP detects proposal collapse and duplication missed by standard ESS. We also illustrate its use as a targeted feedback signal for proposal rejuvenation.

Pretrained diffusion models represent image distributions through a continuum of progressively smoothed distributions. This multiscale structure organizes generation from global structure to fine detail and supports high-quality, diverse samples. We exploit the same multiscale diffusion prior for linear imaging inverse problems. Rather than using the pretrained model only as a denoiser in an outer iteration, we define a surrogate likelihood whose center is aligned with the clean-image coordinate and whose covariance accounts for residual diffusion uncertainty. This construction defines an explicit surrogate posterior path, from which we derive continuous posterior dynamics. A tunable Langevin component supports target tracking and allows the amount of posterior exploration to be adapted to the application. We prove endpoint consistency and a finite-horizon tracking bound and, in the exact-score setting, first-order weak accuracy. For computation, we derive the Posterior-Dynamics Implicit–Explicit sampler (PD-IMEX), a stable method using one score evaluation per diffusion scale and an implicit data-consistency update. Experiments on deblurring, super-resolution, and inpainting show strong reconstruction quality at 100 score evaluations, coarse-grid stability, and controllable fidelity–diversity behavior.

Sufficient dimension reduction (SDR) seeks the minimal subspace of the predictors that captures the full conditional distribution of the response, which is known as the central subspace (CS). When the response is multivariate, the problem becomes considerably more challenging, particularly when the sample size is limited. Existing methods face different limitations:inverse regression approaches rely on strong distributional assumptions and matrix inversion, and their multi-response extensions suffer from severe slice sparsity; forward regression methods depend on computationally intensive iterative smoothing whose cost grows with the response dimension; and deep learning-based approaches demand large amounts of labeled data. To circumvent these shortcomings, we propose an SDR framework based on the generalized Stein's lemma. Our method constructs a cross-moment matrix between the multivariate response and the marginal score function of the predictors, and recovers the CS via its singular value decomposition. The proposed method does not rely on the linearity condition, avoids matrix inversion as well as iterative smoothing, and can leverage unlabeled data when available. We establish convergence guarantees for the proposed estimator under standard regularity conditions. Moreover, we propose a practical rank-selection algorithm to estimate the dimension of the CS. Extensive simulation studies and a real data application demonstrate that the proposed methods consistently outperform existing approaches across a variety of settings, particularly in moderate-dimensional, label-scarce scenarios with high noise levels.

Diffusion models learn to reverse a predefined corruption process, but sampling still requires a costly time discretization and depends on the chosen noise schedule. We study these two issues for variance-preserving diffusions with matrix-valued schedules. Our analysis transfers reverse-time discretization errors to the forward corruption law and treats two numerical schemes within a common framework. The first freezes the score and yields, through a matrix-sensitive local comparison and forward information dissipation, an ambient-dimensional step complexity with leading factor $d/\varepsilon^2$ for KL accuracy $\varepsilon^2$. The second keeps the known Gaussian drift exact and freezes the posterior mean. For data of metric-entropy dimension $k$, a forward Markov identity, an anisotropic covering estimate, and Stieltjes integration by parts give the corresponding factor $k\log k/\varepsilon^2$. In both cases, the proof identifies a local error, accumulates it through the forward evolution, and inserts the result into a common KL decomposition. The local errors further provide directional criteria for matrix schedules and an asymptotically optimal square-root adaptive grid. A high-dimensional Gaussian-mixture experiment illustrates the resulting schedule and grid improvements.

In this paper, we study the effects of employing multiplier bootstrap to analyze the asymptotic distributions of the largest eigenvalues of high-dimensional sample covariance matrices in both spiked and non-spiked models. Our findings demonstrate that the multiplier bootstrap establishes several phase transitions in the limiting edge distributions of both unconditional and conditional bootstrapped covariance matrices, provided the different classes of multipliers. In the nonspiked setting, unbounded multipliers lead to Frechet or Gumbel limits for the largest eigenvalue of the bootstrapped covariance matrix, both conditionally on the observed data and unconditionally. For bounded multipliers, the unconditional model exhibits transitions among Tracy-Widom, Gaussian, or Weibull limits, determined jointly by the aspect ratio p/n, the upper-endpoint behavior of the multipliers, and the population covariance matrix. The conditional model displays analogous Gaussian and Weibull regimes; in contrast, the conditional counterpart of the unconditional Tracy-Widom regime collapses to a point mass. In the spiked setting, under suitable signal-strength conditions, the leading eigenvalues of both the unconditional and conditional bootstrapped sample covariance matrices are asymptotically Gaussian for bounded as well as unbounded multipliers, under some mild assumptions. Our theoretical results also clarify the feasibility and adaptability of the multiplier bootstrap for spectral inference in high-dimensional sample covariance models. Numerical simulations confirm the accuracy of our results and the effectiveness of the proposed spectral inference procedures, which may be of independent interest.

We study anomaly detection in temporally dependent network sequences. Methods based only on adjacency matrices, which are widely used for static networks, can miss changes in the way edges evolve. We instead represent each pair of consecutive networks by separate formation and dissolution event matrices, and embed the resulting matrix sequences using unfolded adjacency spectral embedding. Comparing these embeddings against a stationary baseline yields vertex- and network-level statistics that detect an anomalous transition and identify whether it involves formation, dissolution, or both. For autoregressive random dot product graphs, we establish uniform rowwise consistency of the transition-event embeddings and derive high-probability guarantees for vertex- and network-level detection and exact recovery of the anomalous vertex set. The theory separates a vertex's own displacement from interference caused by other vertices and from autoregressive memory. We quantify the memory left by an earlier anomaly, show that it decays geometrically after the transition mechanism returns to baseline, and give conditions under which it is negligible. A degree-corrected stochastic block model extension gives exact community recovery and provides community-reassignment, centre-shift, split, and merge anomaly statistics with high-probability detection and recovery guarantees. Simulations and applications to an international trade dataset and a primary-school contact network illustrate the performance of the proposed method, revealing a dissolution-driven trade decline followed by formation-driven recovery during COVID-19 and temporary merge-type mixing between school classes.

We prove a quantitative central limit theorem for linear functionals of regularized empirical-risk minimizers in the proportional-dimensional regime \(p=O(n)\). The data columns are independent, not necessarily identically distributed, and satisfy a uniform columnwise Poincaré inequality. Under uniform curvature and smoothness assumptions, and for a quadratic regularizer, we show that every nondegenerate statistic \(\sqrt n\,u^\top\hat\theta\), centered by its expectation and normalized by its standard deviation, converges to a standard normal random variable in Wasserstein distance, with rate \(O((\log n)^7n^{-1/4})\). The proof is based on moment and stability bounds for the minimizer, a second-order leave-one-out expansion, and a perturbative normal-approximation argument for functions of independent variables. We also prove the variance upper bound \(\Var(u^\top\hat\theta)\le C\norm{u}_2^2/n\), identifying the \(\sqrt n\) fluctuation scale.

This paper presents a framework for constructing structure-preserving representations of complex systems through discrete approximation, and demonstrates its use in studying the vertical temperature and salinity structures in the mesopelagic zone across the global ocean using the ARGO dataset. Clustering serves as a means of organizing complexity into a finite set of structures that approximate the overall oceanic conditions, and a color encoding design then integrates these structures into a coherent map, with the three color components derived from interpretable geometric features of a profile: its initial level, its magnitude of variation, and its shape. Instead of focusing on specific depth levels or computing zonal averages within selected regions, our approach preserves the full vertical structure of individual profiles and incorporates each profile in the global ocean, capturing both fine-scale profile detail and large-scale spatial variability. By clustering over one million profiles collected over a decade, we identify and characterize representative profile shapes, which form the basis for a visualization strategy that provides an integrated, comprehensive, and interpretable presentation of the large-scale spatial distributions of these oceanic vertical patterns.

Transformers applied to spatial imperfect-information games must represent map geometry while tracking hidden entities through time. We ask whether geometry-aware positional encodings improve these capabilities, without claiming a new positional encoding. We construct a four-level benchmark on a hexagonal naval pursuit game: controlled geometry and topology probes, an exact-Bayes hidden-target tracking task, offline policy imitation at 1k and 10k games, and 7,200 fixed-seed games against three legacy opponents. Across matched Transformer backbones, HexRoPE reduces exact-belief posterior cross-entropy relative to no positional encoding by 0.278 on D6-transformed test orbits and 0.329 on a larger map; both hierarchical-bootstrap confidence intervals exclude zero, and both Holm-adjusted p-values are below 0.001. At 1k games, HexRoPE improves policy action accuracy by 4.63 percentage points over no encoding and 2.05 points over rectangular relative bias; the gains shrink to 1.55 and 0.41 points at 10k games. However, HexRoPE does not improve aggregate gameplay win rate: its paired effect over no encoding is -1.56 percentage points (95% CI [-4.50, 1.17]). Rectangular relative bias is strongest on D6 belief consistency but fails sharply when extrapolating from radius 3 to radius 4, while graph bias provides only a small blocked-edge gain. The results show that geometric inductive bias improves belief estimation and data-efficient imitation, but those representation gains do not automatically produce stronger closed-loop play.

We consider nonparametric regression when the association between a response and its covariates changes across an unknown partition of a spatial domain. The proposed estimator learns the partition and the cluster-specific regression functions jointly. A neural network depending only on location determines cluster membership, while separate neural networks describe the covariate–response relationship within the clusters. An annealed softmax relaxation permits gradient-based estimation of the otherwise discrete assignments. Graph-Laplacian and occupancy penalties are used to discourage fragmented regions and degenerate solutions. We establish identifiability up to label permutation, bound partition error under a margin condition, and decompose prediction risk into regression and assignment components. The resulting rate agrees with that of an oracle estimator when the partition is estimated sufficiently accurately. Simulations show that joint estimation is useful when regression surfaces change abruptly across spatial boundaries, including settings with nonlinear effects, unequal region sizes, preferential sampling, and spatially correlated errors. Finally, a real data analysis is provided to demonstrate the validity and effectiveness of the proposed method.

Low-rank matrix decompositions can uncover patterns and structure in data and have a number of different applications across many disciplines. Extensions to "joint" low-rank decompositions have been proposed to link datasets from different modalities. While these methods enable the discovery of common patterns across modalities, they require that all the multimodal data share one or more dimensions. We propose a new analysis method, PathFinder, that enables co-analysis of datasets that do not necessarily all share a dimension. The key insight is that as long as pairs or subgroups of matrices do share some dimension, and that there are one or more paths that link across the data matrices, a global joint decomposition can be sought out. This enables the joint estimation of common patterns across different modalities, species, or scales, where a one-to-one mapping across all data along some dimension is not necessarily available. We show that PathFinder is a general umbrella under which many matrix decomposition methods fall as special cases. It can be used to discover common patterns across disparate datasets and to make predictions for missing data or modalities.

Watermarking provides a principled way to authenticate text generated by large language models (LLMs). In practice, however, the final text may be mixed-source, with watermark evidence surviving at only a subset of token positions after rewriting, insertion, deletion, or paraphrasing. Although prior work has studied global detection of watermark signals, when such signals can be localized remains unclear. We formulate watermark localization as a token-level multiple-testing problem based on pivotal statistics, with a latent indicator recording whether watermark dependence survives at each position. Under an asymptotic regime indexed by exponents for signal sparsity, next-token concentration, and effective-vocabulary growth, we derive a sharp boundary for global detection and phase transitions for discovery and classification within the class of coordinatewise pivot-based localization rules. We show that discovery is strictly harder than detection and that consistent classification is impossible across the parameter regime within this class. We then develop an adaptive thresholding method that does not require knowledge of the exponents or time-varying next-token distributions, but uses a data-driven estimate of the surviving watermark fraction. The method attains the optimal discovery boundary and near-optimal discovery power relative to homogeneous pivot-based rules. Simulations support the theoretical phase transitions, while experiments on model-generated texts demonstrate practical localization performance under common edit mechanisms.

Thinned count autoregressions are popular for modelling infectious disease surveillance data due to their flexibility and interpretability. However, such a model is challenging to fit since it involves high-dimensional and serially correlated integer-valued unknowns. One solution is to consider an analogous continuous-valued surrogate model whose values are post-hoc mapped to integers. This procedure produces biased estimates as inference is performed using samples from such a surrogate model and not the thinned count autoregression itself. In this work, we propose a sequential importance sampling procedure to correct this misspecified model. We demonstrate its validity in a simulation study and its applicability for epidemic curve reconstruction using rotavirus data from Germany and meningococcus data from France.

The Widely Applicable Information Criterion (WAIC) is a cornerstone of Bayesian model selection, but its conditional independence assumption renders it inappropriate for sequential and spatially correlated data - a limitation that leads to systematically underestimated model complexity and over-optimistic predictive assessments. We introduce CC-WAIC, a principled generalization of WAIC that explicitly incorporates the full posterior covariance structure of log-likelihood contributions. CC-WAIC reduces exactly to WAIC when independence holds, making it a natural extension rather than an ad-hoc modification. To overcome the prohibitive computational cost of the full covariance matrix, we develop a linear-time implementation using a banded covariance approximation that reduces complexity, with theoretical guarantees on the approximation error. We further introduce an effective sample size correction to mitigate finite-sample bias in MCMC estimation. Through extensive simulations on Hidden Markov Models and real-world applications to Old Faithful geyser data and S&P 500 volatility modelling, we demonstrate that CC-WAIC substantially outperforms WAIC, LOO-CV, iWAIC, and WAICNF, particularly under strong temporal dependence and limited sample sizes. The proposed criterion offers a computationally scalable and theoretically grounded tool for Bayesian model selection in the dependent data settings that are ubiquitous across modern science. Limitations include reliance on exponential mixing and exact conditional likelihoods; extensions to long-memory processes and approximate inference are discussed.

High-dimensional multi-way (tensor) data pose significant challenges for covariance estimation due to the curse of dimensionality. We introduce a unified framework for scalable estimation of tensor covariances based on a Kronecker-structured sparse inverse Cholesky (KSIC) projection. Our approach is grounded in the geometry of information projection, defining the estimator as the moment-matching projection of a target distribution onto a manifold characterized by sparse, Kronecker-factored inverse Cholesky factors. By leveraging physical or data-driven nearest-neighbor sparsity, KSIC provides a geometry-aware representation that is both statistically interpretable and computationally efficient. Our framework integrates two estimation regimes: a nonparametric estimator that projects the empirical covariance directly onto the manifold, utilizing the KSIC structure to implicitly regularize rank-deficient data; and a parametric estimator that fits generative covariance models (e.g., Matérn) by maximizing the likelihood of their KSIC projections, formulated as a nested double forward Kullback-Leibler minimization. Theoretically, we establish the conditions for the existence of the KSIC projection and finite-sample concentration rates for the nonparametric regime, proving that the KSIC estimator gainfully exploits cross-mode information and is robust to data scarcity. Numerical experiments demonstrate that the proposed KSIC estimators achieve state-of-the-art accuracy and scalability, particularly in settings with high dimensionality and limited sample sizes. We apply KSIC to spatiotemporal temperature anomalies and functional MRI data, demonstrating its broad applicability across diverse multi-way data domains.

Modeling multivariate distributions with nonlinear dependence, multimodality, and tractable analytical structure for downstream applications is a central challenge in uncertainty quantification. Sliced Normal (SN) distributions were introduced in prior works at the National Aeronautics and Space Administration (NASA) to address this need by representing densities through polynomial feature maps. This construction provides a compact algebraic alternative to more opaque generative models, while retaining the ability to capture nonlinear parameter dependencies and multi-modal behavior. In this paper, we build on the SN framework and develop several improvements that make the approach more reliable and scalable. First, we reformulate SN parameter estimation as a convex optimization problem over a positive semidefinite matrix, replacing the original nonconvex likelihood search with a formulation amenable to standard optimization tools. Second, we clarify the expressive power of the SN class by connecting polynomial log-density modeling to a Stone–Weierstrass-type universal approximation argument on compact domains. Third, we propose a high-dimensional fitting procedure that partitions variables into approximately independent groups, fits SN models within each subgroup, and then assembles the subgroup models through a cross-block completion step to recover residual dependence. We demonstrate the resulting SN modeling pipeline on NASA loss-of-control flight data, where the method captures nonlinear dependence patterns in both low-dimensional slices and a higher-dimensional block-assembled model.

Objective: Small-sample molecular classification requires feature selectors that identify predictive, stable, and nonredundant subsets for binary and multiclass outcomes. We propose ARISE (Adaptive Residual-Informed Stability Ensemble), which integrates complementary relevance signals, class-balanced stability assessment, residual-informed redundancy control, and multiclass pairwise coverage. Methods: ARISE combines seven percentile-normalized relevance components through 15 predefined profiles, adaptively weighted by nested inner cross-validation. It was evaluated on five molecular datasets, eight feature-set sizes, three fixed classifiers (k-nearest neighbours, support vector machine, and random forest), and six filter comparators. Generalization was estimated by five-fold outer cross-validation repeated 50 times using balanced accuracy, macro-F1, and Cohen's kappa. Results: Across 210,000 held-out assessments, ARISE ranked first in all 15 dataset-metric combinations. Equal-dataset means were 0.793 for balanced accuracy, 0.776 for macro-F1, and 0.725 for kappa, exceeding the strongest aggregate comparator by 0.022, 0.023, and 0.028, respectively. Performance remained strong across compact feature sets, although the optimal budget differed by dataset. Conclusion: ARISE provides a transparent, adaptive framework that jointly addresses relevance, stability, redundancy, and multiclass discrimination. Its consistent results across datasets, classifiers, metrics, and feature-set sizes support further evaluation for small-sample molecular classification.