We propose a new procedure for testing whether two networks are edge-correlated through some latent vertex correspondence. The test statistic is based on counting the cooccurrences of signed trees for a family of nonisomorphic trees. When the two networks are Erdős–Rényi random graphs $\mathcal{G}(\mathit{n},\mathit{q})$ that are either independent or correlated with correlation coefficient ρ, our test runs in ${\mathit{n}}^{2\mathbf{+}\mathit{o}(1)}$ time and succeeds with high probability as $\mathit{n}\to \infty$, provided that $\mathit{n}min\{\mathit{q},1-\mathit{q}\}\ge {\mathit{n}}^{-\mathit{o}(1)}$ and ${\mathit{\rho }}^2>\mathit{\alpha }\approx 0.338$, where α is Otter’s constant so that the number of unlabeled trees with K edges grows as ${(1/\mathit{\alpha })}^{\mathit{K}}$. This significantly improves the prior work in terms of statistical accuracy, running time and graph sparsity.

A seminal result in the ICA literature states that for $\mathit{A}\mathit{Y}=\mathit{\epsilon }$, if the components of ε are independent and at most one is Gaussian, then A is identified up to sign and permutation of its rows (Signal Process. 36 (1994)). In this paper we study to which extent the independence assumption can be relaxed by replacing it with restrictions on higher order moment or cumulant tensors of ε. We document new conditions that establish identification for several nonindependent component models, for example, common variance models, and propose efficient estimation methods based on the identification results. We show that in situations where independence cannot be assumed the efficiency gains can be significant relative to methods that rely on independence.

The extremal dependence structure of a regularly varying random vector $\mathit{X}$ is fully described by its limiting spectral measure. In this paper, we investigate how to recover characteristics of the measure, such as extremal coefficients, from the extremal behaviour of convex combinations of components of $\mathit{X}$. Our considerations result in a class of new estimators of moments of the corresponding combinations for the spectral vector. We show asymptotic normality by means of a functional limit theorem and, focusing on the estimation of extremal coefficients, we verify that the minimal asymptotic variance can be achieved by a plug-in estimator using subsampling bootstrap. We illustrate the benefits of our approach on simulated and real data.

We consider the sequential decision-making problem where the mean outcome is a nonlinear function of the chosen action. Compared with the linear model, two curious phenomena arise in nonlinear models: first, in addition to the “learning phase” with a standard parametric rate for estimation or regret, there is an “burn-in period” with a fixed cost determined by the nonlinear function; second, achieving the smallest burn-in cost requires new exploration algorithms. For a special family of nonlinear functions named ridge functions in the literature, we derive upper and lower bounds on the optimal burn-in cost, and in addition, on the entire learning trajectory during the burn-in period via differential equations. In particular, a two-stage algorithm that first finds a good initial action and then treats the problem as locally linear is statistically optimal. In contrast, several classical algorithms, such as UCB and algorithms relying on regression oracles, are provably suboptimal.

We study the tensor-on-tensor regression, where the goal is to connect tensor responses to tensor covariates with a low Tucker rank parameter tensor/matrix without prior knowledge of its intrinsic rank. We propose the Riemannian gradient descent (RGD) and Riemannian Gauss–Newton (RGN) methods and cope with the challenge of unknown rank by studying the effect of rank over-parameterization. We provide the first convergence guarantee for the general tensor-on-tensor regression by showing that RGD and RGN respectively converge linearly and quadratically to a statistically optimal estimate in both rank correctly-parameterized and over-parameterized settings. Our theory reveals an intriguing phenomenon: Riemannian optimization methods naturally adapt to over-parameterization without modifications to their implementation. We also prove the statistical-computational gap in scalar-on-tensor regression by a direct low-degree polynomial argument. Our theory demonstrates a “blessing of statistical-computational gap” phenomenon: in a wide range of scenarios in tensor-on-tensor regression for tensors of order three or higher, the computationally required sample size matches what is needed by moderate rank over-parameterization when considering computationally feasible estimators, while there are no such benefits in the matrix settings. This shows moderate rank over-parameterization is essentially “cost-free” in terms of sample size in tensor-on-tensor regression of order three or higher. Finally, we conduct simulation studies to show the advantages of our proposed methods and to corroborate our theoretical findings.

Confidence intervals based on the central limit theorem (CLT) are a cornerstone of classical statistics. Despite being only asymptotically valid, they are ubiquitous because they permit statistical inference under weak assumptions and can often be applied to problems even when nonasymptotic inference is impossible. This paper introduces time-uniform analogues of such asymptotic confidence intervals, adding to the literature on confidence sequences (CS)—sequences of confidence intervals that are uniformly valid over time—which provide valid inference at arbitrary stopping times and incur no penalties for “peeking” at the data, unlike classical confidence intervals which require the sample size to be fixed in advance. Existing CSs in the literature are nonasymptotic, enjoying finite-sample guarantees but not the aforementioned broad applicability of asymptotic confidence intervals. This work provides a definition for “asymptotic CSs” and a general recipe for deriving them. Asymptotic CSs forgo nonasymptotic validity for CLT-like versatility and (asymptotic) time-uniform guarantees. While the CLT approximates the distribution of a sample average by that of a Gaussian for a fixed sample size, we use strong invariance principles (stemming from the seminal 1960s work of Strassen) to uniformly approximate the entire sample average process by an implicit Gaussian process. As an illustration, we derive asymptotic CSs for the average treatment effect in observational studies (for which nonasymptotic bounds are essentially impossible to derive even in the fixed-time regime) as well as randomized experiments, enabling causal inference in sequential environments.

Tensor time series, which is a time series consisting of tensorial observations, has become ubiquitous. It typically exhibits high dimensionality. One approach for dimension reduction is to use a factor model structure, in a form similar to Tucker tensor decomposition, except that the time dimension is treated as a dynamic process with a time dependent structure. In this paper, we introduce two approaches to estimate such a tensor factor model by using iterative orthogonal projections of the original tensor time series. These approaches extend the existing estimation procedures and improve the estimation accuracy and convergence rate significantly as proven in our theoretical investigation. Our algorithms are similar to the higher-order orthogonal projection method for tensor decomposition, but with significant differences due to the need to unfold tensors in the iterations and the use of autocorrelation. Consequently, our analysis is significantly different from the existing ones. Computational and statistical lower bounds are derived to prove the optimality of the sample size requirement and convergence rate for the proposed methods. Simulation study is conducted to further illustrate the statistical properties of these estimators.

Segmented regression models offer model flexibility and interpretability as compared to the global parametric and the nonparametric models, and yet are challenging in both estimation and inference. We consider a four-regime segmented model for temporally dependent data with segmenting boundaries depending on multivariate covariates with nondiminishing boundary effects. A mixed integer quadratic programming algorithm is formulated to facilitate the least square estimation of the regression and the boundary parameters. The rates of convergence and the asymptotic distributions of the least square estimators are obtained for the regression and the boundary coefficients, respectively. We propose a smoothed regression bootstrap to facilitate inference on the parameters and a model selection procedure to select the most suitable model within the model class with at most four segments. Numerical simulations and a case study on air pollution in Beijing are conducted to demonstrate the proposed approach, which shows that the segmented models with three or four regimes are suitable for the modeling of the meteorological effects on the PM${}_{2.5}$ concentration.

High-dimensional distributions, especially those with heavy tails, are notoriously difficult for off-the-shelf MCMC samplers: the combination of unbounded state spaces, diminishing gradient information, and local moves results in empirically observed “stickiness” and poor theoretical mixing properties—lack of geometric ergodicity. In this paper, we introduce a new class of MCMC samplers that map the original high-dimensional problem in Euclidean space onto a sphere and remedy these notorious mixing problems. In particular, we develop random-walk Metropolis type algorithms as well as versions of the Bouncy Particle Sampler that are uniformly ergodic for a large class of light and heavy-tailed distributions and also empirically exhibit rapid convergence in high dimensions. In the best scenario, the proposed samplers can enjoy the “blessings of dimensionality” that the convergence is faster in higher dimensions.

Gaussian deterministic approximations are routinely employed in Bayesian statistics to ease inference when the target posterior is intractable. While these approximations are justified, in asymptotic regimes, by Bernstein–von Mises type results, in practice the expected Gaussian behavior might poorly represent the actual shape of the target posterior, thus affecting approximation accuracy. Motivated by these considerations, we derive an improved class of closed-form and valid deterministic approximations of posterior distributions that arise from a novel treatment of a third-order version of the Laplace method yielding approximations within a tractable family of skew-symmetric distributions. Under general assumptions accounting for misspecified models and non-i.i.d. settings, such a family of approximations is shown to have a total variation distance from the target posterior whose convergence rate improves by at least one order of magnitude the one achieved by the Gaussian from the classical Bernstein–von Mises theorem. Specializing this result to the case of regular parametric models shows that the same accuracy improvement can be also established for the posterior expectation of polynomially bounded functions. Unlike available higher-order approximations based on, for example, Edgeworth expansions, our results prove that it is possible to derive closed-form and valid densities which provide a more accurate, yet similarly tractable, alternative to Gaussian approximations of the target posterior, while inheriting its limiting frequentist properties. We strengthen these arguments by developing a practical skew-modal approximation for both joint and marginal posteriors which preserves the guarantees of its theoretical counterpart by replacing the unknown model parameters with the corresponding maximum a posteriori estimate. Simulation studies and real-data applications confirm that our theoretical results closely match the empirical gains observed in practice.

Deep neural networks have achieved tremendous success due to their representation power and adaptation to low-dimensional structures. Their potential for estimating structured regression functions has been recently established in the literature. However, most of the studies require the input dimension to be fixed, and consequently, they ignore the effect of dimension on the rate of convergence and hamper their applications to modern big data with high dimensionality. In this paper, we bridge this gap by analyzing a k-way nonparametric interaction model in both growing dimension scenarios (d grows with n but at a slower rate) and in high dimension ($\mathit{d}\gtrsim \mathit{n}$). In the latter case, sparsity assumptions and associated regularization are required to obtain optimal convergence rates. A new challenge in diverging dimension setting is in calculation mean-square error; the covariance terms among estimated additive components are an order of magnitude larger than those of the variances and can deteriorate statistical properties without proper care. We introduce a critical debiasing technique to amend the problem. We show that under certain standard assumptions, debiased deep neural networks achieve a minimax optimal rate both in terms of $(\mathit{n},\mathit{d})$. Our proof techniques rely crucially on a novel debiasing technique that makes the covariances of additive components negligible in the mean-square error calculation. In addition, we establish the matching lower bounds.

The “sample amplification” problem formalizes the following question: Given n i.i.d. samples drawn from an unknown distribution P, when is it possible to produce a larger set of $\mathit{n}\mathbf{+}\mathit{m}$ samples which cannot be distinguished from $\mathit{n}\mathbf{+}\mathit{m}$ i.i.d. samples drawn from P? In this work, we provide a firm statistical foundation for this problem by deriving generally applicable amplification procedures, lower bound techniques and connections to existing statistical notions. Our techniques apply to a large class of distributions including the exponential family, and establish a rigorous connection between sample amplification and distribution learning.

Under the usual nonparametric regression model with Gaussian errors, Least Squares Estimators (LSEs) over natural subclasses of convex functions are shown to be suboptimal for estimating a d-dimensional convex function in squared error loss when the dimension d is 5 or larger. The specific function classes considered include: (i) bounded convex functions supported on a polytope (in random design), (ii) Lipschitz convex functions supported on any convex domain (in random design) and (iii) convex functions supported on a polytope (in fixed design). For each of these classes, the risk of the LSE is proved to be of the order ${\mathit{n}}^{-2/\mathit{d}}$ (up to logarithmic factors) while the minimax risk is ${\mathit{n}}^{-4/(\mathit{d}\mathbf{+}4)}$, when $\mathit{d}\ge 5$. In addition, the first rate of convergence results (worst case and adaptive) for the unrestricted convex LSE are established in fixed design for polytopal domains for all $\mathit{d}\ge 1$. Some new metric entropy results for convex functions are also proved, which are of independent interest.

Estimation problems with constrained parameter spaces arise in various settings. In many of these problems, the observations available to the statistician can be modelled as arising from the noisy realization of the image of a random linear operator; an important special case is random design regression. We derive sharp rates of estimation for arbitrary compact elliptical parameter sets and demonstrate how they depend on the distribution of the random linear operator. Our main result is a functional that characterizes the minimax rate of estimation in terms of the noise level, the law of the random operator, and elliptical norms that define the error metric and the parameter space. This nonasymptotic result is sharp up to an explicit universal constant, and it becomes asymptotically exact as the radius of the parameter space is allowed to grow. We demonstrate the generality of the result by applying it to both parametric and nonparametric regression problems.

Testing the significance of a variable or group of variables X for predicting a response Y, given additional covariates Z, is a ubiquitous task in statistics. A simple but common approach is to specify a linear model, and then test whether the regression coefficient for X is nonzero. However, when the model is misspecified, the test may have poor power, for example, when X is involved in complex interactions, or lead to many false rejections. In this work, we study the problem of testing the model-free null of conditional mean independence, that is, that the conditional mean of Y given X and Z does not depend on X. We propose a simple and general framework that can leverage flexible nonparametric or machine learning methods, such as additive models or random forests, to yield both robust error control and high power. The procedure involves using these methods to perform regressions, first to estimate a form of projection of Y on X and Z using one-half of the data, and then to estimate the expected conditional covariance between this projection and Y on the remaining half of the data. While the approach is general, we show that a version of our procedure using spline regression achieves what we show is the minimax optimal rate in this nonparametric testing problem. Numerical experiments demonstrate the effectiveness of our approach both in terms of maintaining Type I error control, and power, compared to several existing approaches.

Random matrix theory has become a widely useful tool in high-dimensional statistics and theoretical machine learning. However, random matrix theory is largely focused on the proportional asymptotics in which the number of columns grows proportionally to the number of rows of the data matrix. This is not always the most natural setting in statistics where columns correspond to covariates and rows to samples.

With the objective to move beyond the proportional asymptotics, we revisit ridge regression (${\ell }_2$-penalized least squares) on i.i.d. data $({\mathit{x}}_{\mathit{i}},{\mathit{y}}_{\mathit{i}})$, $\mathit{i}\le \mathit{n}$, where ${\mathit{x}}_{\mathit{i}}$ is a feature vector and ${\mathit{y}}_{\mathit{i}}=⟨\mathit{\beta },{\mathit{x}}_{\mathit{i}}⟩\mathbf{+}{\mathit{\epsilon }}_{\mathit{i}}\in \mathbb{R}$ is a response. We allow the feature vector to be high-dimensional, or even infinite-dimensional, in which case it belongs to a separable Hilbert space, and assume either ${\mathit{z}}_{\mathit{i}}:={\mathrm{\Sigma }}^{-1/2}{\mathit{x}}_{\mathit{i}}$ to have i.i.d. entries, or to satisfy a certain convex concentration property.

Within this setting, we establish nonasymptotic bounds that approximate the bias and variance of ridge regression in terms of the bias and variance of an “equivalent” sequence model (a regression model with diagonal design matrix). The approximation is up to multiplicative factors bounded by $(1\pm \mathrm{\Delta })$ for some explicitly small Δ.

Previously, such an approximation result was known only in the proportional regime and only up to additive errors: in particular, it did not allow to characterize the behavior of the excess risk when this converges to 0. Our general theory recovers earlier results in the proportional regime (with better error rates). As a new application, we obtain a completely explicit and sharp characterization of ridge regression for Hilbert covariates with regularly varying spectrum. Finally, we analyze the overparametrized near-interpolation setting and obtain sharp “benign overfitting” guarantees.

This paper considers the problem of testing and estimation of change point where signals after the change point can be highly irregular, which departs from the existing literature that assumes signals after the change point to be piecewise constant or vary smoothly. A two-step approach is proposed to effectively estimate the location of the change point. The first step consists of a preliminary estimation of the change point that allows us to obtain unknown parameters for the second step. In the second step, we use a new procedure to determine the position of the change point. We show that, under suitable conditions, the desirable ${\mathcal{O}}_{\mathbb{P}}(1)$ rate of convergence of the estimated change point can be obtained. We apply our method to analyze the Baidu search index of COVID-19 related symptoms and find December 8, 2019, to be the starting date of the COVID-19 pandemic.

This paper considers decentralized Federated Learning (FL) under heterogeneous distributions among distributed clients or data blocks for the M-estimation. The mean squared error and consensus error across the estimators from different clients via the decentralized stochastic gradient descent algorithm are derived. The asymptotic normality of the Polyak–Ruppert (PR) averaged estimator in the decentralized distributed setting is attained, which shows that its statistical efficiency comes at a cost as it is more restrictive on the number of clients than that in the distributed M-estimation. To overcome the restriction, a one-step estimator is proposed which permits a much larger number of clients while still achieving the same efficiency as the original PR-averaged estimator in the nondistributed setting. The confidence regions based on both the PR-averaged estimator and the proposed one-step estimator are constructed to facilitate statistical inference for decentralized FL.

This paper presents asymptotic results for the maximum likelihood and restricted maximum likelihood (REML) estimators within a two-way crossed mixed effect model, when the number of rows, columns, and the number of observations per cell tend to infinity. The relative growth rate for the number of rows, columns, and cells is unrestricted, whether considered pairwise or collectively. Under very mild conditions (which include moment conditions instead of requiring normality for either the random effects or errors), the estimators are proven to be asymptotically normal, with a structured covariance matrix. We also discuss the case where the number of observations per cell is fixed at 1.