Papers

Recovering dynamical systems from noisy observations is a recurring challenge across scientific domains, including neuroscience and physics. Latent stochastic differential equations (SDEs) address this by modeling the system as an unobserved state that evolves according to a learnable SDE and generates the observations. Variational inference (VI) provides a tractable objective for fitting latent SDEs. Traditional VI algorithms evaluate this objective by numerical simulation over a time discretization, trading fidelity for computational cost. A recent class of algorithms, simulation-free VI, sidesteps this tradeoff by parameterizing the posterior through its instantaneous marginals rather than its drift. In this work, we show that the efficiency of existing simulation-free VI algorithms comes at a price: their parameterizations restrict the approximate posterior to a subset of the SDEs available to simulation-based methods, degrading posterior inference and parameter learning. We propose Helmholtz-SDE, a simulation-free VI algorithm that closes this gap by optimizing over path laws compatible with a prescribed collection of marginals. Helmholtz-SDE recovers dynamics more faithfully than prior simulation-free methods, with the largest gains under high posterior uncertainty. It further matches the performance of simulation-based VI at a fraction of the runtime. Code is available at https://github.com/lindermanlab/helmholtz-sde .

Functional regressors complicate inference in linear regression problems so that the bootstrap can play a useful role in quantifying uncertainty and calibrating intervals. The best bootstrap in practice, though, can depend on factors in the data as well as computational considerations and existing bootstraps can have limitations: residual bootstrap is computationally fast and simple but may fail when the errors are heterogeneous, while paired bootstrap applies more generally in functional linear regression at a cost of much higher computation. To bridge this gap, we develop a wild bootstrap method for functional linear regression, which is akin to a modified version of residual bootstrap but designed to have a wide scope of application like paired bootstrap, including to heteroscedastic errors. Its theoretical consistency is established and numerical studies suggest that wild bootstrap can provide accurate and computationally fast inference. Importantly, we also suggest a practical and effective approach of selecting truncation levels, specifically designed for mean response inference problems. The proposed bootstrap in functional linear regression is further illustrated through a weather data example, and an accompanying R package BTSinFLRM provides numerical implementations.

Longitudinal clinical studies often collect repeated measurements of biomarkers or health-related quality of life together with a time-to-event outcome. These processes are intrinsically linked: longitudinal trajectories may predict event risk, while event occurrence, or its anticipation, can induce informative censoring of the longitudinal process. Joint models provide a principled framework for handling this dependence, but most existing formulations rely on proportional hazards assumptions that may be restrictive and offer limited interpretability on the time scale. We propose a class of semiparametric accelerated failure time joint models that directly model covariate effects on event timing while flexibly capturing longitudinal-event associations. The survival component is specified through an accelerated failure time model with the baseline component represented by a flexible basis expansion, allowing a broad class of smooth baseline specifications. We illustrate the framework using Bernstein polynomial baseline representations and introduce rescaling strategies to improve numerical stability and parameter identifiability under time-warping. Estimation is conducted within a Bayesian framework, enabling joint inference for longitudinal, survival, and association parameters. Simulation studies reflecting realistic longitudinal trajectories, censoring mechanisms, and dependence structures are used to evaluate finite-sample performance. The proposed models show improved recovery of longitudinal treatment effects compared with a standalone linear mixed model when event risk depends on the underlying longitudinal process. Overall, the framework extends existing joint modelling methodology by offering a flexible and interpretable alternative to proportional hazards-based approaches.

Sampling from complex, unnormalized probability densities is a fundamental challenge in Bayesian inference and probabilistic modeling. While Markov chain Monte Carlo (MCMC) methods provide asymptotic guarantees, they often suffer from slow mixing and high computational costs due to fixed or manually tuned trajectory lengths. In this work, we propose a novel framework that treats trajectory termination as a learnable component of the sampling dynamics. By framing MCMC within the theory of non-acyclic generative flow networks (GFlowNets), we train state-dependent neural classifiers to decide when a trajectory has reached a high-density region and should terminate. We theoretically establish the connection between optimal classifiers and the target density via detailed balance conditions and introduce a multilevel training scheme to facilitate exploration in complex geometries. Experimental results across various benchmark densities demonstrate that our approach significantly reduces average trajectory lengths while improving mode coverage and mixing compared to standard MCMC baselines.

Generalized estimating equations (GEE) are widely used for population-averaged inference on correlated binary responses, but ordinary GEE can fail under separation, a situation that is more likely in small-sample, sparse, or rare-event settings, leading to nonconvergence, infinite or extreme estimates, and unreliable inference. Existing penalized GEE (PGEE) approaches mitigate some of these problems but do not generally guarantee finite estimates under nonindependence working structures and often rely on correlation-coefficient parameterizations whose admissible range shrinks as fitted probabilities approach zero or one, forcing the working association toward independence under separation. We propose a PGEE framework that combines a Jeffreys-prior penalty with marginalized odds-ratio working parameterizations. The odds-ratio parameterization avoids this failure, while the penalty, with tunable strength $\delta$ and default $\delta = 1/2$, stabilizes estimation under separation. Under working independence, PGEE reduces to the Jeffreys-prior penalized maximum-likelihood estimator, yielding finite estimates for logit, probit, complementary log-log, and cauchit links. Under nonindependence odds-ratio structures, where a formal finiteness guarantee is unavailable, PGEE achieves near-complete empirical convergence even in separated settings. We also propose one-step and hybrid variants, OPGEE and HPGEE, that reduce computational cost. Simulations show that all three variants substantially outperform ordinary GEE under separation while retaining the performance of ordinary GEE in regular settings. We illustrate the method using a respiratory-illness trial in which ordinary GEE fails, and provide an implementation in the R package geer.

Generalized estimating equations (GEE) are widely used for correlated data, but with small to moderate numbers of independent clusters the ordinary GEE regression estimators can be substantially biased. We develop a first-order bias-reduction principle for GEE by viewing the estimator as a clustered-data $M$-estimator and deriving an adjustment to the estimating equations that targets the leading bias term while accounting for the dependence of the working covariance on the mean parameters. The resulting class includes three bias-reduced estimators and three one-step bias-corrected analogs, nesting the bias-corrected estimator of Lunardon and Scharfstein (2017) and the bias-reduced and bias-corrected estimators of Paul and Zhang (2014) as special cases. The framework applies to general response types through correlation-coefficient parameterizations for the association structure and extends to correlated binary data through pairwise odds-ratio parameterizations, yielding the first bias-reduced and bias-corrected GEE estimators under this parameterization, for which the marginal-mean compatibility constraints are far less restrictive than those of correlation-coefficient parameterizations, making them better suited for small-sample settings. Under standard regularity conditions, all six estimators share the same asymptotic distribution as the ordinary GEE. Simulation studies show that the proposed estimators reduce bias while maintaining efficiency and coverage close to those of ordinary GEE across a range of settings, and a clinical trial analysis illustrates the proposed estimators in practice. Software is available in the R package geer.

Interpreting machine-learning models has attracted increasing attention, particularly in the physical sciences, where one often seeks to understand the underlying mechanisms rather than merely make predictions. Multiple linear regression is often regarded as an interpretable alternative to more complex models, such as deep neural networks, because its predictions are expressed as explicit weighted sums of input features. However, when input features are strongly correlated, namely in the presence of multicollinearity, the learned weights can exhibit large dataset-to-dataset fluctuations and oscillatory behavior across physically similar features, making their interpretation difficult or even impossible. Although the instability of the weights under multicollinearity is well known in statistics, its consequences for physical interpretation, in particular its connection to oscillatory weights across physically similar features, have not been systematically clarified. Here, we theoretically discuss the mechanism behind this loss of interpretability by analyzing the eigenmodes of the feature correlation matrix. We show that small-eigenvalue modes associated with multicollinearity amplify fluctuations in the weights and generate oscillatory patterns that do not necessarily reflect meaningful contributions. We test this theoretical picture numerically on physics datasets and show that Ridge regularization suppresses these unstable modes, although the resulting weights must still be interpreted with caution. We further confirm the generality of our findings beyond physics by analyzing a diverse collection of publicly available datasets. Our results clarify why, in the presence of multicollinearity, physical interpretation can remain difficult even for linear regression models.

Foundation models are now used in settings where the prompts they receive can change quickly. Users change, topics change, policies change, and the model may suddenly face a kind of request that was rare in the calibration data. This makes fixed calibration risky. Conformal prediction and conformal risk control give model-agnostic ways to control error, but they work best when the calibration data still look like the future data. This paper develops PromptShift CRC, a drift-aware conformal risk control method for foundation-model outputs under prompt and domain shift. The method embeds prompts and responses, measures how far the current prompt stream has moved from the calibration pool, gives more weight to relevant or recent calibration examples, and updates the risk level online after observed violations. It reports three practical diagnostics: realized risk error, prompt drift, and effective calibration size. We give conditions under which the method controls risk up to terms for distribution mismatch and weighted quantile uncertainty. In a synthetic prompt-shift benchmark, static conformal risk control fails sharply after drift, while PromptShift-CRC gives the best coverage among the adaptive baselines considered. We then evaluate the same calibration layer on public benchmark derived streams for question answering, toxicity, summarization factuality, and long-context hallucination risk

Finding D-optimal designs for generalized linear models (GLMs) is challenging due to the dependence of the Fisher information matrix on unknown parameters and the lack of closed-form solutions, particularly when input factors include both discrete and continuous variables. Although classical algorithms and recent metaheuristic approaches have offered partial solutions, there remains a need for robust and computationally efficient methods. In this paper, we propose a penalized Particle Swarm Optimization (PSO) approach, named $p$-PSO. Here we introduce a new, general-purpose penalty formulation for constrained optimization and demonstrate its effectiveness in optimal design problems. The formulation is algorithm-agnostic and applicable to a broad class of black-box optimization methods. Results show that the method is highly efficient, with its primary contribution being a penalty formulation that enables the direct use of an off-the-shelf PSO algorithm and extends naturally to more general constrained optimization tasks.

Conformal prediction provides distribution-free prediction intervals under exchangeability, but many modern data streams are neither independent nor stable. They exhibit recurring regimes, changing seasonal frequencies, abrupt shifts, and gradual drift. We propose drift-aware spectral conformal prediction (DASC), a streaming uncertainty quantification framework for structured non-exchangeable data subject to distributional drift. DASC forms conformal prediction intervals using calibration residuals weighted by local spectral similarity, while a transport-based drift score monitors whether the current test distribution has moved away from past calibration regimes. When drift is mild, DASC borrows calibration residuals from structurally similar historical windows; when drift is severe, it contracts or reweights the calibration pool and updates the target miscoverage level online. The method also reports an effective sample size diagnostic that warns when a weighted conformal quantile is statistically fragile. We establish an approximate coverage bound that decomposes coverage loss into drift, residual mismatch, and weighted effective sample size. In synthetic experiments and five stress-test regimes, DASC maintains near-nominal coverage after drift where rolling, recency-weighted, and spectral-only conformal methods can under-cover. In real electricity and weather streams, DASC reduces average interval width by approximately 28% and 42%, respectively, relative to the best calibrated non-DASC baseline, while preserving calibrated or conservative coverage. A financial volatility example shows a more nuanced regime in which spectral-only calibration is competitive, but DASC retains near-nominal coverage and adds drift diagnostics.

Conformal prediction gives prediction intervals with finite-sample coverage when the data are exchangeable. Many time-indexed datasets are not exchangeable: they have seasons, recurring regimes, changing frequencies, or other forms of structured dependence. This paper studies a simple way to use that structure. We propose spectral adaptive conformal prediction, a method that forms weighted conformal quantiles using local spectral similarity and then updates the target miscoverage level online. The spectral weights choose calibration residuals that look relevant to the current test point. The adaptive update corrects the long-run miss rate when uncertainty changes over time. The theory makes both parts controllable. We give an approximate coverage bound that splits the error into a spectral mismatch term and an effective-sample-size term, prove that kernel spectral weighting never increases the mismatch term relative to uniform weighting, show that a bandwidth of order N^(-1/(d+2)) balances the two terms, and establish an unconditional long-run calibration bound for the adaptive update that holds for every sample path without independence or stationarity. Simulations with recurring regimes and slowly changing frequencies, together with four real-data examples spanning monthly, weekly, and daily U.S. and European series, show when the hybrid method improves on strong adaptive baselines and when it does not, and an effective-sample-size safeguard, computable at prediction time without outcomes, detects and repairs the one observed failure.

Generalized linear mixed models (GLMMs) are widely used for analyzing correlated data, such as longitudinal and multilevel data. With over 15 $\texttt{R}$ packages available on $\texttt{CRAN}$ for fitting GLMMs, practitioners face a difficult choice regarding which package yields accurate estimates, converges reliably, and offers reasonable computational speed. Existing comparisons are either limited to methods within a single package or focus on narrow criteria such as speed alone. To address this gap, we systematically compared seven representative $\texttt{R}$ packages – $\texttt{lme4}$, $\texttt{GLMMadaptive}$, $\texttt{glmmTMB}$, $\texttt{MASS}$, $\texttt{hglm}$, $\texttt{brms}$, and $\texttt{rstanarm}$ – that implement different estimation frameworks. By using Monte Carlo simulations across 24 scenarios, we evaluated each package in terms of convergence ratios, computational time, estimation accuracy, and hypothesis testing performance. Our results showed that $\texttt{lme4_AGQ}$ and $\texttt{GLMMadaptive}$ yield the highest accuracy and convergence ratios, although $\texttt{GLMMadaptive}$ becomes slower under complex random-effect structures. $\texttt{lme4_LA}$ and $\texttt{glmmTMB}$ are computationally fast but exhibit lower convergence ratios and larger bias, especially for variance components. $\texttt{MASS}$ and $\texttt{hglm}$ are also fast, but $\texttt{MASS}$ yields liberal univariate tests and $\texttt{hglm}$ lacks support for correlated random effects and multivariate testing. Between two Bayesian packages, $\texttt{rstanarm}$ converges reliably and produces valid univariate tests, whereas $\texttt{brms}$ is extremely slow, limiting its practical utility. Based on these findings, we provide practical recommendations for choosing GLMM tool in applied research.

Flow matching is a powerful generative modeling framework, valued for its simplicity and strong empirical performance. However, its standard formulation treats signals on structured spaces, such as fMRI data on brain graphs, as points in Euclidean space, overlooking the rich topological features of their domains. To address this, we introduce topological flow matching, a topology-aware generalization of flow matching. We interpret flow matching as a framework for solving a degenerate Schrödinger bridge problem and inject topological information by augmenting the reference process with a Laplacian-derived drift. This principled modification captures the structure of the underlying domain while preserving the desirable properties of flow matching: a stable, simulation-free objective and deterministic sample paths. As a result, our framework serves as a drop-in replacement for standard flow matching. We demonstrate its effectiveness on diverse structured datasets, including brain fMRIs, ocean currents, seismic events, and traffic flows.

We introduce biarchetype analysis for the first time in the context of univariate functional data. This unsupervised methodology extends archetype analysis by simultaneously identifying archetypal structures across both the cases (countries, in our application) and the temporal argument. Both cases and time points are expressed as mixtures of biarchetypes, yielding a concise and highly interpretable representation of complex functional observations. Although biarchetype analysis is not intended as a clustering technique, it offers superior interpretability compared with biclustering approaches, as it is based on extreme, representative patterns rather than average centroids, thereby enhancing human comprehension. We apply the proposed method to 10-year government bond yields of European countries over the period 2001-2025. The results identify three distinct time regimes (the pre-crisis period, the euro-area sovereign debt crisis, and the post-crisis period), and reveal Germany, Greece, and Hungary as country archetypes.

A recent methodology for extracting archetypal profiles from three-way asymmetric proximity data is applied to a dataset comprising 23 x 23 x 3 two-mode, three-way asymmetric dissimilarity matrices. The asymmetric dissimilarities are based on the oriented wavelet squared coherence of 10-year government bond yields among 23 European countries across three time intervals from 2001 to the present. First, the h-plot is computed for the unconditional two-mode three-way data, which is then represented in a unified Euclidean space, providing an intuitive and interpretable visualization. Subsequently, archetypoid analysis is performed. This unsupervised methodology identifies the archetypal countries and expresses all remaining countries as mixtures of these archetypal instances. Additionally, the degree of asymmetry for each country is calculated, offering further insight into the structure of the dissimilarities. The dataset and code are provided to support reproducible research.

Uncertainty in the solution of an inverse problem and in the tasks performed on it is quantified by posterior expectations, each an average of an integrand over $M$ posterior samples. While designed quadratures improve on the $O(M^{-1/2})$ error of Monte-Carlo estimation, they solve an optimization problem, often against the posterior density, for every new observation, which can be computationally costly. To address this limitation, we introduce the quadrature field, a set-equivariant network that maps an observation and its $M$ posterior samples to an $M$-node signed-weight quadrature in one forward pass. Trained once on a family of posteriors to minimize the worst-case integration error over a class of functions, it serves any observation, any $M$ and any integrand in that class with no further optimization. We show that, with high probability and up to a computable slack, the resulting quadrature is never worse than the Monte-Carlo estimate built from the same samples. We validate the quadrature field on closed-form and on learned posteriors, one constrained by a partial differential equation, where it improves on the Monte-Carlo estimate in median at every node count, often by orders of magnitude.

Deep neural networks (DNNs) frequently fail to generalize to out-of-distribution (OOD) medical images because of variations in scanners and acquisition protocols. Retraining DNN models to address these distribution shifts is often impractical due to the high cost of acquiring and annotating new medical datasets. To address this, we introduce VarDeepPCA, a novel lightweight variational DNN framework designed to restore/refine degraded segmentation maps by leveraging intrinsic geometric priors. Unlike existing approaches that require target-domain data or extensive pre-training, our VarDeepPCA explicitly learns a distribution of valid anatomical geometries using only small in-distribution (ID) datasets. Theoretically, our novel variational learning framework leverages a reinterpretation of the softmax mapping to implicitly perform exact distribution modeling, thereby enabling computationally efficient, sampling-free learning and inference. This also enables VarDeepPCA to provide uncertainty estimates associated with its restored segmentation maps. We empirically validate our framework across 4 distinct clinical applications, using 14 publicly available datasets, involving segmentation of the myocardium, neuroretinal rim, prostate, and fetal head. Comparisons against 15 existing methods demonstrate that VarDeepPCA consistently restores segmentation maps produced by the existing methods on OOD data to (i) significantly improve anatomical plausibility of geometries and clinical utility of the segmentations, and (ii) significantly reduce errors, without needing any more training data than that used by existing methods.

A single observed network reflects several mechanisms at once: communities, hubs, and clustering coexist in one graph, each a different model. We treat the network as a combination of candidate mechanisms and study, from a single graph, how strongly each mechanism contributes and how they combine. We address two questions. The first is how to measure each mechanism's contribution when the mechanisms must themselves be estimated from the graph: fitting the mechanisms and their strengths from the same data biases the strengths toward zero, and a correction removes this bias and yields valid confidence intervals. The second is whether the rule of combination is itself recoverable: when a graph is generated by two mechanisms acting together, the graph alone determines whether they combine additively or interact, exactly when the graph is dense enough, a sharp threshold below which no test can decide. The estimate calibrates the candidate mechanisms against the observed edges. We establish matching minimax rate, against a known-design benchmark and the estimated-design problem itself, confirm the methods in simulation, and apply them to real networks, where the signed coefficients recover known structure and, in one case, a confidence interval excludes any positive contribution from a candidate mechanism.

This paper explores policy gradient algorithms for training stochastic policies to sample from structured discrete probability distributions under the Generative Flow Network (GFlowNet) framework. Building on extensive theoretical connections between GFlowNets and entropy-regularized reinforcement learning, we derive equivalents of standard policy gradient algorithms for training GFlowNets, as well as experimentally explore their various methodological aspects, including baseline training and advantage estimation. Most importantly, our work is the first to derive and successfully apply proximal policy optimization to GFlowNets, showing its improved convergence speed and data efficiency compared to standard GFlowNet training objectives on benchmarks ranging from synthetic energies to molecular graph generation.

A significant gap exists between theory and practice in deep learning. Generalization and approximation error bounds are often derived for simplified models or are too loose to be informative. Many rely on the manifold hypothesis and on geometric regularity such as intrinsic dimension, curvature, and reach. Progress requires insight into data-manifold geometry and suitable benchmarks, yet existing options are polarized: analytic manifolds with known geometry but limited applicability, or real-world datasets where geometry is only coarsely estimable. We introduce a benchmarking framework for studying data geometry. We repurpose and extend dSprites and COIL-20 with additional transformation dimensions and dense, axis-aligned sampling, and pair them with finite-difference estimators that recover curvature, reach, and volume at near-ground-truth accuracy in a regime where general-purpose estimators are unreliable or difficult to deploy. The framework is intended as a controlled testbed, useful as a calibration environment for geometric estimators and a sandbox for probing theoretical assumptions. To illustrate its use, we present two application studies, namely assessing the scaling behavior of the bounds of Genovese et al. and Fefferman et al., and tracking the layer-wise geometry of a $\beta$-VAE, highlighting the behavior of current bounds and the value of controlled benchmarks for guiding and validating future theory. A reference implementation is available at https://github.com/koulakis/manifold-microscope.

We introduce the Multiplex Network Hawkes model, which extends the network Hawkes framework of Linderman & Adams (2014) by allowing multiple excitation layers whose weights depend on observed edge and node covariates. We use the model to investigate how contagion in financial networks is affected by different transmission channels. The multiplex structure separates channel-specific contributions within a single inferred transmission network, allowing candidate propagation mechanisms to be compared directly rather than being absorbed into one homogeneous excitation layer. Covariate-dependent excitation allows us to investigate sources of transmission. We make posterior inference about the inferred directed network and its excitation dynamics using an MCMC sampler. The application uses a broad cross-industry credit default swap (CDS) dataset of 99 North American and European firms, including banks, insurers and non-financial firms over 2004-2022. We evaluate three candidate contagion channels associated with asset similarity, solvency and profitability. The results indicate sparse contagion pathways, with systemic-risk transmission concentrated in outward flows from a small number of influential institutions rather than in mutual feedback between institutions. The channel results show that industry similarity is the most consistently supported asset-similarity effect, while aggregate layer contributions indicate that asset-similarity, solvency and profitability channels all contribute to inferred excitation.

Defining and estimating causal effects for ordinal data is challenging. Standard average treatment effects are not appropriate for ordinal scales, and alternative estimands, such as the probabilities that the treatment outcome exceeds or does not worsen the control outcome, are generally not identifiable. Existing work provides sharp bounds for these quantities based only on marginal distributions. Motivated by a previous observation showing that bounds obtained under an independence working assumption can be substantially tighter, we investigate conditions under which such bounds are valid. We show that commonly used notions of positive dependence, including positive quadrant dependence and positive regression dependence, are not sufficient to justify these bounds. We then propose a new dependence condition, diagonal tail dominance (DTD), under which the independence-based bounds are guaranteed to hold. We explain why this condition is quite strong and may not be appropriate in many settings, limiting the justification for using the independence-based bounds. However, local DTD may be plausible in many applications, and we derive improved bounds that exploit an independence working assumption on selected parts of the probability table. Through theoretical results, numerical examples, and an analysis of data from a clinical trial of a new treatment for acute ischemic stroke, we illustrate the properties of the bounds and the role of the proposed conditions.

In this paper we investigate the efficacy of the score-based martingale posteriors (SMP) (Cui & Walker, 2025; Fong et al., 2023) in the context of modern and large-scale machine learning problems and its potential for meaningful uncertainty quantification. SMPs work with a stochastic gradient ascent-type recursion on the parameter space of stochastic models and construct a martingale on the parameter space. Under simple mathematical assumptions, the recursion can be built so that the parameters form a martingale sequence which possesses a limiting, in time, random variable, the latter of which can be simulated very quickly, in contrast to Monte Carlo-based methods such as Markov chain Monte Carlo. In this expository paper we explore the SMP for inferring the parameters of deep neural networks (DNNs) and, where feasible, compare our results to the state-of-the-art Monte Carlo methods aimed at inferring conventional Bayesian posteriors.

Stochastic trace estimation is a standard tool for approximating the trace of a large-scale matrix available only through matrix-vector products. However, in tensor-structured settings, unstructured Gaussian or Rademacher test vectors may be prohibitively expensive to store and compute with, while cheaper rank-one tensor-product vectors can require sample complexities that grow exponentially with the tensor order. This work studies Gaussian random tensor train vectors as a structured alternative for stochastic trace estimation. We show that, with a suitable choice of the tensor train rank, random tensor train vectors recover dimension-independent guarantees for the Girard–Hutchinson estimator. In particular, a median-of-means variant with tensor train rank $r \geq d-1$ achieves the same dependence on the accuracy $\varepsilon$ and failure probability $\delta$ as the classical estimator based on unstructured Gaussian vectors. We further prove an oblivious subspace injection result for sketches formed from independent Gaussian random tensor train vectors: tensor train rank $r\geq d-1$ and $\mathcal{O}(\varepsilon^{-2}(k+\log(1/\delta)))$ samples suffice for a $k$-dimensional target subspace. Finally, we investigate the use of such sketches within the Nyström++ framework. We show that the resulting estimator can achieve the desired $\mathcal{O}(\varepsilon^{-1})$ sample complexity under an additional spectral-tail condition. These results provide clarififcation on both the potential and the limitations of random tensor train vectors in stochastic trace estimation.

This paper studies the information gap between mixture detection and label recovery in binomial logistic mixtures. Standard likelihood-based criteria such as the Bayesian information criterion (BIC) can detect the presence of two components, but this does not guarantee that the corresponding labels are recoverable. We show that this gap is intrinsic to binomial logistic mixtures with a fixed number of trials: observed-data evidence for mixture structure and per-observation information for label recovery have different local orders in the component separation, and only the former accumulates with the sample size. As a result, there exists a detectable-but-unrecoverable regime in which BIC selects two components while the posterior labels remain essentially uninformative. To address this issue, we propose two feasibility-aware inference procedures: a recoverability-aware BIC with a posterior-entropy penalty and an entropy-regularized estimator that mitigates the tendency of the maximum likelihood estimator to produce overly separated components and overly concentrated posterior responsibilities. Numerical experiments confirm the predicted gap and demonstrate that the proposed methods avoid misleading component selections and improve the calibration of posterior label probabilities.

This paper proposes a novel nonparametric test to assess the uncorrelatedness between two high-dimensional random vectors. We develop our test by generalizing the random integration proposed by Jiang et al. (2023, 2024), and the resulting test statistic estimates a weighted squared $\mathscr{L}_2$ norm of the covariance matrix. Asymptotic properties of the test statistic are derived by letting both the sample size $n$ and the dimension $p$ diverge to infinity. Under the null hypothesis of uncorrelatedness, our proposed test statistic is asymptotically normal with zero mean and unit variance, without requiring any specification of the relative magnitude regarding $n$ and $p$. Monte Carlo simulations demonstrate the good finite-sample performance of our proposed methods. Compared with many existing tests, our test statistic is more powerful at detecting “weak but pervasive” dependence while maintaining a comparable empirical size. The advantages of the proposed methods are further illustrated by an empirical analysis that assesses the correlation between DNA methylation and gene expression.

In high-dimensional semi-supervised linear regression, prediction-powered inference (PPI) corrects an external predictor with a rectifier estimated from the labeled data. In a linear model, however, this rectifier cancels the predictor: PPI and PPI++ reduce to ordinary least squares and can inflate variance when the predictor is close to the oracle. We propose the Debiased External-model-Assisted Lasso (DEAL), which routes the external estimator and the unlabeled covariates into the variance of a debiased estimator, with a bias-aware, cross-fitted shrinkage step that adapts across target-only, near-oracle, and biased-but-informative regimes. We prove coordinate-wise asymptotic normality with an adaptive variance, extend validity to the projection parameter under misspecification and nonlinear labelers, and show that, at a common unlabeled budget, DEAL intervals are shorter than those of debiased Lasso, PPI, and PPI++; a shift-aware variant preserves coverage under covariate shift. In simulations, DEAL intervals are 0.49-0.87 of the debiased-Lasso length, and across six real-data applications spanning astronomy, chemistry, proteomics, and oncology, the last using a large-language-model oracle, they tighten in every case, with median length ratios of 0.23-0.53.

We study the graph alignment problem for correlated Gaussian Orthogonal Ensemble (GOE) matrices, where the goal is to recover a hidden vertex permutation given two correlated symmetric Gaussian matrices $(A, B)$ with correlation $1/\sqrt{1+\sigma^2}$. While the maximum likelihood estimator is information-theoretically optimal, its computation, which reduces to a quadratic assignment problem, is intractable. Motivated by this, we analyze convex relaxations based on minimizing $\|AX - XB\|_F$ over the set of doubly stochastic matrices and the unit hypercube. We show that when the correlation parameter satisfies $\sigma = o(n^{-1/2}/\log^4 n)$, the solution of either relaxation $(X^\star)$ concentrates around the ground-truth permutation matrix $(\Pi^\star)$, i.e., $\|X^\star-\Pi^\star\|_F^2 = o(n)$, implying recovery of all but a vanishing fraction of vertices after simple post-processing. Combined with existing lower bounds, our results precisely characterize that $\|X^\star-\Pi^\star\|_F^2$ transitions from $o(n)$ for $\sigma = \tilde{o}(n^{-1/2})$ to $\Omega(n)$ for $\sigma = \tilde{\Omega}(n^{-1/2})$. In doing so, our analysis significantly tightens prior results and extends them beyond doubly stochastic relaxations.

Test-time training (TTT) adapts a pretrained model to each prompt via parameter updates, improving accuracy under pretraining-to-test distribution shifts. Yet, its performance often suffers from instability and sensitivity to hyperparameters such as update steps and subspace. We explain this behavior through a decision-theoretic lens, treating TTT as implicit Bayesian inference in the kernel regime. Under a Gaussian process benchmark, we show that TTT reduces prediction error when updates are spectrally matched to the prompt's signal-to-noise ratio and aligned with query-relevant eigen-directions. This perspective underpins the following results: (1) we show when fixed update steps and subspaces fail under distribution shifts, motivating adaptive strategies; (2) we prove that selecting update steps via prompt evidence admits a PAC-Bayes guarantee against overfitting; and (3) we characterize the Bayes-optimal update subspace under a linear-Gaussian correction model, yielding a scoring rule for selecting Transformer blocks and heads. Our theory helps explain the empirical instability of TTT, taking a step toward principled guidance for when, how far, and which directions to adapt.

In large language model (LLM) serving, each request accumulates persistent graphics processing unit (GPU) memory during service as its key-value cache grows with every generated token. Under high concurrency, aggregate memory usage therefore increases endogenously over time: the service process itself creates future capacity pressure. When memory capacity is exceeded, systems evict active requests, discarding cached state and restarting them later, which wastes computation and reduces throughput. We develop a discrete-time dynamical model of memory-constrained LLM inference that captures admission, memory growth, and eviction under continuous batching. In the saturated-input regime, the system admits both eviction-free fixed points and limit cycles with evictions. For homogeneous workloads, we show that the eviction-free equilibrium is unstable and that, except for a Lebesgue-measure-zero exact-capture set, the system converges to a unique worst-case limit cycle that is asymptotically stable outside this exceptional set, with throughput losses as large as 50%. For heterogeneous workloads, we prove a stability criterion in the two-class common-input setting and explain how the survival-polynomial mechanism generalizes to multiple classes and heterogeneous-input lengths. Under an input-dominated scaling regime, coprime decoding lengths stabilize the eviction-free equilibrium, while non-coprime lengths create synchronized modes that drive instability. These results characterize when workload heterogeneity desynchronizes completions and helps stabilize memory-constrained serving. More broadly, we identify service-induced congestion as a structural instability mechanism and derive scheduling design principles for sustaining high throughput.