This article presents two new multilevel hierarchical matrices to perform fast kernel matrix-vector products in $d$ dimensions $(d>1)$, namely $\mathcal{H}_{hw}^{\prime}$, $\mathcal{H}_{hw}^{\prime \prime}$. We introduce a novel hybrid weak admissibility condition to construct these hierarchical matrices. More precisely, at coarser levels of the tree (i.e., larger blocks), we use a weak admissibility condition, wherein \textbf{far-field} and \textbf{vertex-sharing} clusters are considered admissible, while at the finer levels (i.e., smaller blocks), the \textbf{adjacent} clusters are treated as admissible. The admissible blocks at coarser levels of $\mathcal{H}_{hw}^{\prime}$ and $\mathcal{H}_{hw}^{\prime \prime}$ are represented using semi-nested and fully nested bases, respectively. In both hierarchical matrices, the finer levels at which adjacent clusters are treated as admissible are constructed in a non-nested fashion. Thus, they can be viewed as a hybrid of $\mathcal{H}$ and $\mathcal{H}^2$ matrices. We use adaptive cross approximation (ACA) and nested cross approximation (NCA) to construct non-nested and nested bases, respectively. The fast matrix-vector products are performed using the proposed hierarchical matrix formats. To evaluate the performance and usefulness of the proposed hierarchical matrices, we conduct some numerical experiments. The numerical results show that they outperform their uniform admissibility counterparts in terms of storage and matrix-vector product cost, demonstrating the effectiveness of the proposed hybrid formats. We also show that the proposed formats are competitive with the standard (or strong) admissibility-based $\mathcal{H}^2$-matrices. Furthermore, we showcase some interesting applications of the proposed hierarchical matrices.
Jan Blechschmidt, Oliver G. Ernst, Moritz Poguntke, Björn Sprungk
Sampling from the posterior is the central task of computational Bayesian inverse problems. The standard workhorse in Bayesian inference---Markov chain Monte Carlo (MCMC)---is sequential, yields correlated samples, and must be rerun for each new observation. Conditional flow matching offers an amortized alternative: a transport map, trained once on joint samples of parameter and data, that yields independent approximate posterior samples for any observation at negligible online cost, without new likelihood evaluations. We give a careful, MCMC-literate assessment of flow matching for PDE-based inverse problems with function-valued parameters. Exploiting the flow's tractable density, we derive computable accuracy estimates for the learned posterior in total variation distance and Kullback--Leibler divergence and, moreover, propose a hybrid sampler that is made asymptotically exact by Metropolization. We validate the accuracy estimates and demonstrate amortization in several numerical examples, including electrical impedance tomography and a likelihood-free Lotka--Volterra model. Flow matching reproduces posterior means and credible bands closely at a fraction of the forward solves required by MCMC, but its accuracy degrades for concentrated posteriors and for observations that are atypical under the prior.
We study Gibbs averaging of finite Markov chains with a prescribed stationary distribution $π$. Given a transition matrix $P$ and the Gibbs kernel $G$ induced by a state-space partition, we prove that $GPG$ is an information projection for every admissible Bregman generator $φ$, with divergences defined using transition probabilities relative to $π$. A Pythagorean identity characterizes this projection and recovers the relative entropy and weighted Frobenius cases. A dual characterization identifies $GPG$ as a kernel closest to stationarity among those with the same projection chain as $P$; for relative entropy, it maximizes the entropy rate. The identities also give monotonicity under partition coarsening. For every reversible ergodic baseline $P$, Gibbs averaging cannot increase the squared weighted Frobenius distance to stationarity at any positive integer time, and the inequality is strict unless $P=GPG$. The corresponding mixing time therefore cannot increase. To choose the partition, we identify the divergence of $G$ from stationarity as the exact worst-case divergence of $GPG$ over all $π$-stationary baselines $P$. This objective depends only on the block masses. We give an exact algorithm for concave block costs, including those arising from relative entropy and power generators $φ(u)=u^p$ with $1<p<2$, and a greedy approximation algorithm with a provable guarantee for $p>2$. For power generators $φ(u)=u^p$ with $p>1$, the objective is proportional to Tsallis entropy of order $2-p$: its minimization favors imbalanced block masses when $1<p<2$ and balanced block masses when $p>2$, while all $k$-block partitions have the same objective value when $p=2$. A numerical experiment compares the resulting projected chains with a local Metropolis sampler on a multi-modal $π$.
Mame Diarra Toure, Quang Hung Pham, Mohamed Ryad Zemouri, David Stephens, Luc Vouligny, Martin Gagnon
Industrial systems continuously generate time series that vary in sampling frequency, scale, dynamics, and operating regime. Consequently a single time-series foundation model may not be able to represent all regimes effectively and retraining a large model for every new dataset or regime shift is computationally expensive. Modular foundation models address this issue by reusing specialised encoders trained on diverse datasets and by allowing the expansion of capabilities through the addition of new experts. This approach requires compact experts to minimize storage cost and a proper selection mechanism of experts for unseen data. We propose a modular foundation model framework that combines end-to-end Bayesian encoder training with uncertainty based routing. Using Singular Bayesian Neural Networks (SBNNs), we build a collection of compact Bayesian encoders whose epistemic uncertainty guides label-free selection. At the selected ranks, our Bayesian collection uses 63.8\% less storage than its deterministic counterpart. We evaluate this framework on time-series data through distribution shift, encoder selection, and probabilistic forecasting. Across seven incoming datasets, uncertainty-based routing is competitive with distance-based routing in terms of forecast MSE, while reducing average ECE from 0.086 to 0.072 and average NLL from 0.907 to 0.836. It also increases PICP90 from 82.28\% to 89.03\%. In a separate model-family comparison, we evaluate the SBNN pipeline against full-rank Bayes by Backprop, deep ensembles, and MC Dropout. The SBNN pipeline achieves an average MSE of 0.390, comparable with 0.396 for deep ensembles, 0.389 for MC dropout, and 0.443 for the full-rank Bayesian model. It uses 63.9\% fewer parameters than MC dropout, 81.9\% fewer than the full-rank model, and 92.8\% fewer than the deep ensemble.
Particle Gibbs samplers with backward sampling are widely used for sampling the smoothing distribution of state-space models, but, with few particles, the new trajectory often coincides with the current one at many time steps. A known but rarely used remedy replaces the random selection of the final particle by a Metropolis--Hastings move which avoids the current particle. We show that, with backward sampling, this remedy only helps at the last time steps, because its effect on the update rates decays geometrically with the distance to the final time. We propose instead Metropolized backward sampling, which makes such a move at every step of the backward pass, forcibly switching trajectories whenever possible. Additionally, we introduce a coupled version which prevents proposing the current trajectory, maximizing the chances of an actual switch versus the current trajectory and performing perfect sampling when the state-space model has no observations, a feature not shared by backward sampling. We prove that, at every time step, both methods dominate backward sampling as they update the state more often and have smaller lag-one autocorrelations than it does. Our methods are drop-in replacements for backward sampling, and can be used to improve any particle Gibbs sampler.
We introduce fastC-GC and fastC-GC$^*$, computationally efficient methods for time-series graph estimation using the two-test rule motivated by Reichenbach's common cause principle (RCCP) and the conditioning sets of causalised Granger causality (c-GC and c-GC$^*$). The methods use standardized data, ridge-residual correlation, and complete circular shift enumeration evaluated by Fourier transforms. Shared precision calculations and Schur downdates avoid repeated regression fits. On a single-threaded benchmark, fastC-GC was 68-277 times faster than the compared c-GC implementation, and fastC-GC$^*$ was 192-886 times faster than c-GC$^*$. Across four suites with reference graphs, aggregate mean area under the precision-recall curve (AUPRC) differed from the corresponding baseline by at most 0.019, but F1 was lower at the selected thresholds in every suite. The computational gains therefore accompany a change in estimator and selected graphs. Among eight methods, Bagged PCMCI+ had the highest mean AUPRC on sparse VAR, an autocorrelated structural process, and DREAM3; CUTS+ led on Lorenz-96. Wind-tunnel and river-discharge applications illustrate differences between inferred graphs without establishing causal accuracy. We prove the computational identities and a restricted population adjacency result for fastC-GC$^*$. Causal interpretation requires the stated graphical assumptions; validity of the implemented tests additionally requires a suitable shift-invariant null law. Empirical null rejection rates informed the conditional threshold $β= 0.005$, without establishing general test validity or graph-wide error control.
Zhenyu Wang, Molei Liu, Jing Lei, Francis Bach, Zijian Guo
When synthesizing multi-source high-dimensional data, a key objective is to extract low-dimensional representations that effectively approximate the original features across different sources. Such representations facilitate the discovery of transferable structures and help mitigate systematic biases such as batch effects. We introduce Stable Principal Component Analysis (StablePCA), a distributionally robust framework for constructing stable latent representations by maximizing the worst-case explained variance over multiple sources. A primary challenge in extending classical PCA to the multi-source setting lies in the nonconvex rank constraint, which renders the StablePCA formulation a nonconvex optimization problem. To overcome this challenge, we conduct a convex relaxation of StablePCA and develop an efficient Mirror-Prox algorithm to solve the relaxed problem, with global convergence guarantees. Since the relaxed problem generally differs from the original formulation, we further introduce a data-dependent certificate to assess how well the algorithm solves the original nonconvex problem and establish the condition under which the relaxation is tight. Finally, we explore alternative distributionally robust formulations of multi-source PCA based on different loss functions.
The ratio of two densities provides a direct characterization of their differences. We consider the two-sample comparison problem by estimating this ratio given i.i.d. observations from two distributions. To this end, we propose additive tree models for density ratio estimation along with efficient algorithms using a new loss function, the balancing loss. The loss allows tree-based models to be trained using several algorithms originally designed for supervised learning, such as forward-stagewise optimization and gradient boosting. Moreover, the balancing loss resembles an exponential family kernel, and it can serve as a pseudo-likelihood with conjugate priors. This property enables generalized Bayesian inference on the density ratio using backfitting samplers designed for Bayesian additive regression trees (BART). Our Bayesian strategy provides uncertainty quantification for the inferred density ratio, which is critical for applications involving high-dimensional and data-limited distributions with potentially substantial uncertainty. We further show connections of the balancing loss to the exponential loss in binary classification and to the variational form of f-divergence, particularly the squared Hellinger distance. Numerical experiments demonstrate that the proposed boosting algorithms achieve high accuracy at low computational cost, and that the Bayesian additive tree is competitive in accuracy while additionally providing uncertainty quantification, which the compared estimators do not offer. Finally, we demonstrate its application to assessing the quality of generative models for microbiome compositional data.
Penalized maximum likelihood is a widely used approach to estimate a sparse precision or a sparse covariance matrix. A zero in the precision matrix encodes a conditional independence and a zero in the covariance matrix a marginal independence, and the two estimators have been developed separately, each with its own algorithm. We show that these two problems can be unified via a single reparametrization of the last column of either matrix. This reparametrization decouples the likelihood function for each problem to a one-dimensional minimization followed by a penalized multiple regression. Our approach applies to any penalty that is separable in the off-diagonal entries and linear on the diagonal. When the penalty is the $\ell_1$ norm, it unifies the graphical Lasso and the covariance graphical Lasso algorithms.
Papers in statistical theory rest on long proofs that few readers check line by line. We present KAPPA (Kit of Agents for Proof-checking Papers with a proof Assistant), open-source software that organizes language model agents, run by Claude Code, to translate the proofs of a LaTeX paper into the proof assistant Lean 4 with its library Mathlib and check them there, result by result. Each numbered result receives one of four verdicts: verified unconditionally; verified conditionally, with every cited or unproved input stated as a named hypothesis; issue recorded, with its place in the paper and the side at fault; or unfinished, with the reason. The kit follows the paper's own proof step by step, rather than finding a new one. Separate agents write and review, each statement is compared with the paper before it is proved and fingerprinted afterwards, and every verdict rests on evidence that can be rerun without a language model. The kit runs unattended for hours, waits out usage limits, re-checks only what a revision changed and reports its cost by role, model and result. We describe the design, walk through a complete run, and apply the kit to four example papers with answer keys and to three published statistical papers, checked blind in their versions from before public corrections. In all three, the kit found the published error at the step that the correction names, together with further problems that we confirmed by hand. Given the same papers directly, the same language model working with Lean classified fewer results correctly on the longest paper and reported some statements narrower than the paper's as proved; no verdict of the kit claimed more than its Lean evidence showed.
Structured sparse learning is a leading paradigm for high-dimensional data modeling, as it enables rigorous integration of domain-validated structural priors into the learning process. Among these priors, pre-defined grouped structure and piecewise smoothness are ubiquitous in real-world applications. However, most existing structured regularization methods are built on squared loss and convex penalties, rarely jointly model these two complementary structural priors, and suffer from two key limitations. Their optimal tuning parameters are sensitive to the unknown noise level $σ$ due to the squared loss. Meanwhile, convex penalties induce inherent estimation bias and severe over-shrinkage of large magnitudes. To address these issues, we propose a knowledge-driven Nonconvex Structured Smooth Adjustment for Square-Root Regularization (NS$^3$R), which unifies the $σ$-free square-root loss with a joint nonconvex $\ell_{2,q}$ penalty ($0 < q < 1$) to simultaneously encode group sparsity and structured smoothness priors. Theoretically, we establish rigorous non-asymptotic error bounds for the NS$^3$R estimator under a mild $(2,q)$-group restricted compatibility condition. Computationally, we develop an efficient double-screened IRLS with thresholding (DS-IRLS-th) algorithm, which adopts a nested double-screening strategy to reduce computational cost and improve numerical stability, with provable global convergence for the resulting nonconvex optimization problem. Extensive experiments on synthetic and real-world data demonstrate that NS$^3$R consistently outperforms state-of-the-art methods in finite-sample settings.
The Half-Trek Criterion (HTC) is the primary graphical tool for determining generic identifiability of causal effect coefficients in linear structural equation models (SEMs) with latent confounders. However, HTC is inherently node-wise: it simultaneously resolves all incoming edges of a node, leaving a gap of "inconclusive" causal effects (15-23% in moderate graphs). We introduce Iterative Identification Closure (IIC), a general framework that decouples causal identification into two phases: (1) a seed function S_0 that identifies an initial set of edges from any external source of information (instrumental variables, interventions, non-Gaussianity, prior knowledge, etc.); and (2) Reduced HTC propagation that iteratively substitutes known coefficients to reduce system dimension, enabling identification of edges that standard HTC cannot resolve. The core novelty is iterative identification propagation: newly identified edges feed back to unlock further identification -- a mechanism absent from all existing graphical criteria, which treat each edge (or node) in isolation. This propagation is non-trivial: coefficient substitution alters the covariance structure, and soundness requires proving that the modified Jacobian retains generic full rank -- a new theoretical result (Reduced HTC Theorem). We prove that IIC is sound, monotone, converges in O(|E|) iterations (empirically <=2), and strictly subsumes both HTC and ancestor decomposition. Exhaustive verification on all graphs with n<=5 (134,144 edges) confirms 100% precision (zero false positives); with combined seeds, IIC reduces the HTC gap by over 80%. The propagation gain is gamma~4x (2 seeds identifying ~3% of edges to 97.5% total identification), far exceeding gamma<=1.2x of prior methods that incorporate side information without iterative feedback.
Vansh Bansal, Pradip Sharma Poudel, James G. Scott
Simulation-based inference often encodes observations $Y$ into a fixed-dimensional summary $S$, with the posterior conditioning only on $S$. While exact preservation often requires infinite information, any such summary imposes a finite information budget on inference which is rarely an explicit choice or measurement. In this paper, we make the budget explicit by modeling the summary as a lossy channel from $Y$ to $S$, and ask two questions: how large must the budget be, and how best is it spent? First, we show for parametric models that this budget, measured by the mutual information $I(Y;S)$, grows unboundedly as the posterior concentrates, even when fixed-dimensional sufficient statistics exist. More generally, the required budget is set by the distribution of the random posterior, and we bound it from above by quantizing posteriors and from below by converses that weight each posterior state by its probability. Second, we show that a finite budget is best spent stochastically. Under a common-support condition on the simulator, covering the continuous observation spaces typical of simulators, a deterministic summary necessarily loses something: a stochastic summary achieves \emph{strictly} lower distortion than every deterministic summary and every time-shared mixture of deterministic summaries at every positive budget too small for exact preservation. Motivated by these results, we insert a binary bottleneck into neural posterior estimators, measure each trained summary's rate, and compare summaries at matched rate. Across four simulator benchmarks and several encoder families, stochastic bottlenecks achieve favorable rate--loss tradeoffs, and architecture rankings change at different information budgets.
Slender concrete columns reinforced with fiber-reinforced polymer (FRP) bars fail at medium to high axial load levels either by concrete crushing or by loss of stability, while the tensile strains in the reinforcement remain far below the rupture strain. Each non-linear finite element (FE) analysis yields only the capacity of the governing failure mode, and the governing mode changes over a broad zone in which both capacities are nearly equal. Both characteristics conflict with active-learning Kriging-assisted Monte Carlo simulation (AK-MCS), which assumes a single, smooth limit-state function. This paper proposes a piecewise, two-branch AK-MCS framework in which each failure mode is represented by its own Kriging model, trained only on the FE analyses in which that mode governs, and the governing prediction is the minimum of both. Safeguards against data starvation keep rarely governing branches usable. After convergence, additional FE analyses are placed based on a leave-one-out calibration check and on the width of the confidence interval on the reliability index. A design-point search on the surrogate and importance sampling at the design points provide the reliability index even where direct Monte Carlo simulation finds no failure. For an FRP-reinforced concrete column example with mixed failure, 79% and 62% of the training points lie near the mode boundary. The reliability index from the design-point search agrees with the Monte Carlo estimate to within 0.05, and an FE analysis at the design point confirms the surrogate to within 0.6%.