A fundamental concept in multivariate statistics, the sample correlation matrix, is often used to infer the correlation/dependence structure among random variables, when the population mean and covariance are unknown. A natural block extension of it, the sample block correlation matrix, is proposed to take on the same role, when random variables are generalized to random subvectors. In this paper, we establish a spectral theory of the sample block correlation matrices and apply it to group independent tests and related problems, under the high-dimensional setting. More specifically, we consider a random vector of dimension p, consisting of k subvectors of dimension ${\mathit{p}}_{\mathit{t}}$’s, where ${\mathit{p}}_{\mathit{t}}$’s can vary from 1 to order p. Our primary goal is to investigate the dependence of the k subvectors. We construct a random matrix model called sample block correlation matrix based on N samples for this purpose. The spectral statistics of the sample block correlation matrix include the classical Wilks’ statistic and Schott’s statistic as special cases. It turns out that the spectral statistics do not depend on the unknown population mean and covariance, under the null hypothesis that the subvectors are independent. Further, the limiting behavior of the spectral statistics can be described with the aid of the free probability theory. Specifically, under three different settings of possibly N-dependent k and ${\mathit{p}}_{\mathit{t}}$’s, we show that the empirical spectral distribution of the sample block correlation matrix converges to the free Poisson binomial distribution, free Poisson distribution (Marchenko–Pastur law) and free Gaussian distribution (semicircle law), respectively. We then further derive the CLTs for the linear spectral statistics of the block correlation matrix under a general setting. Our results are established under the general distribution assumption on the random vector. It turns out that the CLTs are universal and do not depend on the 4th cumulants of the vector components, due to a self-normalizing effect of the correlation-type matrices. We further derive the CLT under the alternative hypothesis and discuss the power of our statistics. Based on our theory, real data analysis on stock return data and gene data is also conducted.

The problem of joint sequential detection and isolation is considered in the context of multiple, not necessarily independent, data streams. A multiple testing framework is proposed, where each hypothesis corresponds to a different subset of data streams, the sample size is a stopping time of the observations, and the probabilities of four kinds of error are controlled below distinct, user-specified levels. Two of these errors reflect the detection component of the formulation, whereas the other two the isolation component. The optimal expected sample size is characterized to a first-order asymptotic approximation as the error probabilities go to 0. Different asymptotic regimes, expressing different prioritizations of the detection and isolation tasks, are considered. A novel, versatile family of testing procedures is proposed, in which two distinct, in general, statistics are computed for each hypothesis, one addressing the detection task and the other the isolation task. Tests in this family, of various computational complexities, are shown to be asymptotically optimal under different setups. The general theory is applied to the detection and isolation of anomalous, not necessarily independent, data streams, as well as to the detection and isolation of an unknown dependence structure.

We study nonparametric methods for estimating the value function of an infinite-horizon discounted Markov reward process (MRP). We analyze the kernel-based least-squares temporal difference (LSTD) estimate, which can be understood either as a nonparametric instrumental variables method, or as a projected approximation to the Bellman fixed point equation. Our analysis imposes no assumptions on the transition operator of the Markov chain, but rather only conditions on the reward function and population-level kernel LSTD solutions. Using empirical process theory and concentration inequalities, we establish a nonasymptotic upper bound on the error with explicit dependence on the effective horizon $\mathit{H}={(1-\mathit{\gamma })}^{-1}$ of the Markov reward process, the eigenvalues of the associated kernel operator, as well as the instance-dependent variance of the Bellman residual error. In addition, we prove minimax lower bounds over subclasses of MRPs, which shows that our guarantees are optimal in terms of the sample size n and the effective horizon H. Whereas existing worst-case theory predicts cubic scaling (${\mathit{H}}^3$) in the effective horizon, our theory reveals a much wider range of scalings, depending on the kernel, the stationary distribution, and the variance of the Bellman residual error. Notably, it is only parametric and near-parametric problems that can ever achieve the worst-case cubic scaling.

We present an estimator of the covariance matrix Σ of random d-dimensional vector from an i.i.d. sample of size n. Our sole assumption is that this vector satisfies a bounded ${\mathit{L}}^{\mathit{p}}-{\mathit{L}}^2$ moment assumption over its one-dimensional marginals, for some $\mathit{p}\ge 4$. Given this, we show that Σ can be estimated from the sample with the same high-probability error rates that the sample covariance matrix achieves in the case of Gaussian data. This holds even though we allow for very general distributions that may not have moments of order > p. Moreover, our estimator can be made to be optimally robust to adversarial contamination. This result improves the recent contributions by Mendelson and Zhivotovskiy and Catoni and Giulini, and matches parallel work by Abdalla and Zhivotovskiy (the exact relationship with this last work is described in the paper).

We consider estimation of average treatment effects given observational data with high-dimensional pretreatment variables. Existing methods for this problem typically assume some form of sparsity for the regression functions. In this work, we introduce a debiased inverse propensity score weighting (DIPW) scheme for average treatment effect estimation that delivers $\sqrt{\mathit{n}}$-consistent estimates when the propensity score follows a sparse logistic regression model; the outcome regression functions are permitted to be arbitrarily complex. We further demonstrate how confidence intervals centred on our estimates may be constructed. Our theoretical results quantify the price to pay for permitting the regression functions to be unestimable, which shows up as an inflation of the variance of the estimator compared to the semiparametric efficient variance by a constant factor, under mild conditions. We also show that when outcome regressions can be estimated consistently, our estimator achieves semiparametric efficiency. As our results accommodate arbitrary outcome regression functions, averages of transformed responses under each treatment may also be estimated at the $\sqrt{\mathit{n}}$ rate. Thus, for example, the variances of the potential outcomes may be estimated. We discuss extensions to estimating linear projections of the heterogeneous treatment effect function and explain how propensity score models with more general link functions may be handled within our framework. An R package $\mathtt{dipw}$ implementing our methodology is available on CRAN.

The singular subspaces perturbation theory is of fundamental importance in probability and statistics. It has various applications across different fields. We consider two arbitrary matrices where one is a leave-one-column-out submatrix of the other one and establish a novel perturbation upper bound for the distance between the two corresponding singular subspaces. It is well suited for mixture models and results in a sharper and finer statistical analysis than classical perturbation bounds such as Wedin’s theorem. Empowered by this leave-one-out perturbation theory, we provide a deterministic entrywise analysis for the performance of spectral clustering under mixture models. Our analysis leads to an explicit exponential error rate for spectral clustering of sub-Gaussian mixture models. For the mixture of isotropic Gaussians, the rate is optimal under a weaker signal-to-noise condition than that of Löffler et al. (2021).

This paper is concerned with statistical inference for regression coefficients in high-dimensional linear regression models. We propose a new method for testing the coefficient vector of the high-dimensional linear models, and establish the asymptotic normality of our proposed test statistic with the aid of the martingale central limit theorem. We derive the asymptotical relative efficiency (ARE) of the proposed test with respect to the test proposed in Zhong and Chen (J. Amer. Statist. Assoc. 106 (2011) 260–274), and show that the ARE is always greater or equal to one under the local alternative studied in this paper. Our numerical studies imply that the proposed test with critical values derived from its asymptotical normal distribution may retain Type I error rate very well. Our numerical comparison demonstrates the proposed test performs better than existing ones in terms of powers. We further illustrate our proposed method with a real data example.

Permutation tests are widely recognized as robust alternatives to tests based on normal theory. Random permutation tests have been frequently employed to assess the significance of variables in linear models. Despite their widespread use, existing random permutation tests lack finite-sample and assumption-free guarantees for controlling type I error in partial correlation tests. To address this ongoing challenge, we have developed a conformal test through permutation-augmented regressions, which we refer to as PALMRT. PALMRT not only achieves power competitive with conventional methods but also provides reliable control of type I errors at no more than $2\mathit{\alpha }$, given any targeted level α, for arbitrary fixed designs and error distributions. We have confirmed this through extensive simulations.

Compared to the cyclic permutation test (CPT) and residual permutation test (RPT), which also offer theoretical guarantees, PALMRT does not compromise as much on power or set stringent requirements on the sample size, making it suitable for diverse biomedical applications. We further illustrate the differences in a long-Covid study where PALMRT validated key findings previously identified using the t-test after multiple corrections, while both CPT and RPT suffered from a drastic loss of power and failed to identify any discoveries. We endorse PALMRT as a robust and practical hypothesis test in scientific research for its superior error control, power preservation, and simplicity.

We study the Bayesian density estimation of data living in the offset of an unknown submanifold of the Euclidean space. In this perspective, we introduce a new notion of anisotropic Hölder for the underlying density and obtain posterior rates that are minimax optimal and adaptive to the regularity of the density, to the intrinsic dimension of the manifold, and to the size of the offset, provided that the latter is not too small—while still allowed to go to zero. Our Bayesian procedure, based on location-scale mixtures of Gaussians, appears to be convenient to implement and yields good practical results, even for quite singular data.

We study the performance of the spectral method for the phase synchronization problem with additive Gaussian noises and incomplete data. The spectral method utilizes the leading eigenvector of the data matrix followed by a normalization step. We prove that it achieves the minimax lower bound of the problem with a matching leading constant under a squared ${\ell }_2$ loss. This shows that the spectral method has the same performance as more sophisticated procedures including maximum likelihood estimation, generalized power method, and semidefinite programming, as long as consistent parameter estimation is possible. To establish our result, we first have a novel choice of the population eigenvector, which enables us to establish the exact recovery of the spectral method when there is no additive noise. We then develop a new perturbation analysis toolkit for the leading eigenvector and show it can be well-approximated by its first-order approximation with a small ${\ell }_2$ error. We further extend our analysis to establish the exact minimax optimality of the spectral method for the orthogonal group synchronization.

We develop a theory of asymptotic efficiency in regular parametric models when data confidentiality is ensured by local differential privacy (LDP). Even though efficient parameter estimation is a classical and well-studied problem in mathematical statistics, it leads to several nontrivial obstacles that need to be tackled when dealing with the LDP case. Starting from a regular parametric model $\mathcal{P}={({\mathit{P}}_{\mathit{\theta }})}_{\mathit{\theta }\in \mathrm{\Theta }}$, $\mathrm{\Theta }\subseteq {\mathbb{R}}^{\mathit{p}}$, for the i.i.d. unobserved sensitive data ${\mathit{X}}_1,\dots ,{\mathit{X}}_{\mathit{n}}$, we establish local asymptotic mixed normality (along subsequences) of the model

$${\mathit{Q}}^{(\mathit{n})}\mathcal{P}={\left({\mathit{Q}}^{(\mathit{n})}{\mathit{P}}_{\mathit{\theta }}^{\mathit{n}}\right)}_{\mathit{\theta }\in \mathrm{\Theta }}$$

generating the sanitized observations ${\mathit{Z}}_1,\dots ,{\mathit{Z}}_{\mathit{n}}$, where ${\mathit{Q}}^{(\mathit{n})}$ is an arbitrary sequence of sequentially interactive privacy mechanisms. This result readily implies convolution and local asymptotic minimax theorems. In case $\mathit{p}=1$, the optimal asymptotic variance is found to be the inverse of the supremal Fisher information ${sup}_{\mathit{Q}\in {\mathcal{Q}}_{\mathit{\alpha }}}{\mathit{I}}_{\mathit{\theta }}(\mathit{Q}\mathcal{P})\in \mathbb{R}$, where the supremum runs over all α-differentially private (marginal) Markov kernels. We present an algorithm for finding a (nearly) optimal privacy mechanism $\hat{\mathit{Q}}$ and an estimator ${\hat{\mathit{\theta }}}_{\mathit{n}}({\mathit{Z}}_1,\dots ,{\mathit{Z}}_{\mathit{n}})$ based on the corresponding sanitized data that achieves this asymptotically optimal variance.

We provide nonasymptotic rates of convergence of the Wasserstein Generative Adversarial networks (WGAN) estimator. We build neural networks classes representing the generators and discriminators which yield a GAN that achieves the minimax optimal rate for estimating a certain probability measure μ with support in ${\mathbb{R}}^{\mathit{p}}$. The probability μ is considered to be the push forward of the Lebesgue measure on the d-dimensional torus ${\mathbb{T}}^{\mathit{d}}$ by a map ${\mathit{g}}^{\star }:{\mathbb{T}}^{\mathit{d}}\to {\mathbb{R}}^{\mathit{p}}$ of smoothness $\mathit{\beta }\mathbf{+}1$. Measuring the error with the γ-Hölder Integral Probability Metric (IPM), we obtain up to logarithmic factors, the minimax optimal rate $\mathit{O}({\mathit{n}}^{-\frac{\mathit{\beta }\mathbf{+}\mathit{\gamma }}{2\mathit{\beta }\mathbf{+}\mathit{d}}}\vee {\mathit{n}}^{-\frac{1}{2}})$ where n is the sample size, β determines the smoothness of the target measure μ, γ is the smoothness of the IPM ($\mathit{\gamma }=1$ is the Wasserstein case) and $\mathit{d}\le \mathit{p}$ is the intrinsic dimension of μ. In the process, we derive a sharp interpolation inequality between Hölder IPMs. This novel result of theory of functions spaces generalizes classical interpolation inequalities to the case where the measures involved have densities on different manifolds.

The weak convergence of the quantile processes, which are constructed based on different estimators of the finite population quantiles, is shown under various well-known sampling designs based on a superpopulation model. The results related to the weak convergence of these quantile processes are applied to find asymptotic distributions of the smooth L-estimators and the estimators of smooth functions of finite population quantiles. Based on these asymptotic distributions, confidence intervals are constructed for several finite population parameters like the median, the α-trimmed means, the interquartile range and the quantile based measure of skewness. Comparisons of various estimators are carried out based on their asymptotic distributions. We show that the use of the auxiliary information in the construction of the estimators sometimes has an adverse effect on the performances of the smooth L-estimators and the estimators of smooth functions of finite population quantiles under several sampling designs. Further, the performance of each of the above-mentioned estimators sometimes becomes worse under sampling designs, which use the auxiliary information, than their performances under simple random sampling without replacement (SRSWOR).

This paper is concerned with high-dimensional two-sample mean problems, which receive considerable attention in recent literature. To utilize the correlation information among variables for enhancing the power of two-sample mean tests, we consider the setting in which the precision matrix of high-dimensional data possesses a linear structure. Thus, we first propose a new precision matrix estimation procedure with considering its linear structure, and further develop regularization methods to select the true basis matrices and remove irrelevant basis matrices. With the aid of an estimated precision matrix, we propose a new test statistic for the two-sample mean problems by replacing the inverse of sample covariance matrix in the Hotelling test by the estimated precision matrix. The proposed test is applicable for both the low-dimensional setting and high-dimensional setting even if the dimension of the data exceeds the sample size. The limiting null distributions of the proposed test statistic under both null and alternative hypotheses are derived. We further derive the asymptotical power function of the proposed test and compare its asymptotic power with some existing tests. We find the estimation error of the precision matrix does not have impact on the asymptotical power function. Moreover, asymptotic relative efficiency of the proposed test to the classical Hotelling test tends to infinity when the ratio of the dimension of data to the sample size tends to 1. We conduct a Monte Carlo simulation study to assess the finite sample performance of the proposed precision matrix estimation procedure and the proposed high-dimensional two-sample mean test. Our numerical results imply that the proposed regularization method is able to effectively remove irrelevant basis matrices. The proposed test performs well compared with the existing methods especially when the elements of the vector have unequal variances. We also illustrate the proposed methodology by an empirical analysis of a real-world data set.

Given a composite null $\mathcal{P}$ and composite alternative $\mathcal{Q}$, when and how can we construct a p-value whose distribution is exactly uniform under the null, and stochastically smaller than uniform under the alternative? Similarly, when and how can we construct an e-value whose expectation exactly equals one under the null, but its expected logarithm under the alternative is positive? We answer these basic questions, and other related ones, when $\mathcal{P}$ and $\mathcal{Q}$ are convex polytopes (in the space of probability measures). We prove that such constructions are possible if and only if $\mathcal{Q}$ does not intersect the span of $\mathcal{P}$. If the p-value is allowed to be stochastically larger than uniform under $\mathit{P}\in \mathcal{P}$, and the e-value can have expectation at most one under $\mathit{P}\in \mathcal{P}$, then it is achievable whenever $\mathcal{P}$ and $\mathcal{Q}$ are disjoint. More generally, even when $\mathcal{P}$ and $\mathcal{Q}$ are not polytopes, we characterize the existence of a bounded nontrivial e-variable whose expectation exactly equals one under any $\mathit{P}\in \mathcal{P}$. The proofs utilize recently developed techniques in simultaneous optimal transport. A key role is played by coarsening the filtration: sometimes, no such p-value or e-value exists in the richest data filtration, but it does exist in some reduced filtration, and our work provides the first general characterization of this phenomenon. We also provide an iterative construction that explicitly constructs such processes, and under certain conditions it finds the one that grows fastest under a specific alternative Q. We discuss implications for the construction of composite nonnegative (super)martingales, and end with some conjectures and open problems.

This paper considers a multi-environment linear regression model in which data from multiple experimental settings are collected. The joint distribution of the response variable and covariates may vary across different environments, yet the conditional expectations of the response variable, given the unknown set of important variables, are invariant. Such a statistical model is related to the problem of endogeneity, causal inference, and transfer learning. The motivation behind it is illustrated by how the goals of prediction and attribution are inherent in estimating the true parameter and the important variable set. We construct a novel environment invariant linear least squares (EILLS) objective function, a multi-environment version of linear least squares regression that leverages the above conditional expectation invariance structure and heterogeneity among different environments to determine the true parameter. Our proposed method is applicable without any additional structural knowledge and can identify the true parameter under a near-minimal identification condition related to the heterogeneity of the environments. We establish nonasymptotic ${\ell }_2$ error bounds on the estimation error for the EILLS estimator in the presence of spurious variables. Moreover, we further show that the ${\ell }_0$ penalized EILLS estimator can achieve variable selection consistency in high-dimensional regimes. These nonasymptotic results demonstrate the sample efficiency of the EILLS estimator and its capability to circumvent the curse of endogeneity in an algorithmic manner without any additional prior structural knowledge. To the best of our knowledge, this paper is the first to realize statistically efficient invariance learning in the general linear model.

Statistical inference for time series such as curve estimation for time-varying models or testing for existence of a change point have garnered significant attention. However, these works are generally restricted to the assumption of independence and/or stationarity at its best. The main obstacle is that the existing Gaussian approximation results for nonstationary processes only provide an existential proof, and thus they are difficult to apply. In this paper, we provide two clear paths to construct such a Gaussian approximation for nonstationary series. While the first one is theoretically more natural, the second one is practically implementable. Our Gaussian approximation results are applicable for a very large class of nonstationary time series, obtain optimal rates and yet have good applicability. Building on such approximations, we also show theoretical results for change-point detection and simultaneous inference in presence of nonstationary errors. Finally, we substantiate our theoretical results with simulation studies and real data analysis.

Graphon estimation has been one of the most fundamental problems in network analysis and has received considerable attention in the past decade. From the statistical perspective, the minimax error rate of graphon estimation has been established by (Ann. Statist. 43 (2015) 2624–2652) for both stochastic block model (SBM) and nonparametric graphon estimation. The statistical optimal estimators are based on constrained least squares and have computational complexity exponential in the dimension. From the computational perspective, the best-known, polynomial-time estimator is based on universal singular value thresholding (USVT), but it can only achieve a much slower estimation error rate than the minimax one. It is natural to wonder if such a gap is essential. The computational optimality of the USVT or the existence of a computational barrier in graphon estimation has been a long-standing open problem. In this work, we take the first step toward it and provide rigorous evidence for the computational barrier in graphon estimation via low-degree polynomials. Specifically, in SBM graphon estimation, we show that for low-degree polynomial estimators, their estimation error rates cannot be significantly better than that of the USVT under a wide range of parameter regimes and in nonparametric graphon estimation, we show low-degree polynomial estimators achieve estimation error rates strictly slower than the minimax rate. Our results are proved based on the recent development of low-degree polynomials by (Ann. Statist. 50 (2022) 1833–1858), while we overcome a few key challenges in applying it to the general graphon estimation problem. By leveraging our main results, we also provide a computational lower bound on the clustering error for community detection in SBM with a growing number of communities and this yields a new piece of evidence for the conjectured Kesten–Stigum threshold for efficient community recovery. Finally, we extend our computational lower bounds to sparse graphon estimation and biclustering with additive Gaussian noise, and provide discussion on the optimality of our results.

Elliptical distribution is a basic assumption underlying many multivariate statistical methods. For example, in sufficient dimension reduction and statistical graphical models, this assumption is routinely imposed to simplify the data dependence structure. Before applying such methods, we need to decide whether the data are elliptically distributed. Currently existing tests either focus exclusively on spherical distributions, or rely on bootstrap to determine the null distribution, or require specific forms of the alternative distribution. In this paper, we introduce a general nonparametric test for elliptical distribution based on kernel embedding of the probability measure that embodies the two properties that characterize an elliptical distribution: namely, after centering and rescaling, (1) the direction and length of the random vector are independent, and (2) the directional vector is uniformly distributed on the unit sphere. We derive the asymptotic distributions of the test statistic via von Mises expansion, develop the sample-level procedure to determine the rejection region, and establish the consistency and validity of the proposed test. We also develop the concentration bounds of the test statistic, allowing the dimension to grow with the sample size, and further establish the consistency in this high-dimension setting. We compare our method with several existing methods via simulation studies, and apply our test to a SENIC dataset with and without a transformation aimed to achieve ellipticity.

Strict stationarity is an assumption commonly used in time-series analysis in order to derive asymptotic distributional results for second-order statistics, like sample autocovariances and sample autocorrelations. Focusing on weak stationarity, this paper derives the asymptotic distribution of the maximum of sample autocovariances and sample autocorrelations under weak conditions by using Gaussian approximation techniques. The asymptotic theory for parameter estimators obtained by fitting a (linear) autoregressive model to a general weakly stationary time series is revisited and a Gaussian approximation theorem for the maximum of the estimators of the autoregressive coefficients is derived. To perform statistical inference for the aforementioned second-order parameters of interest, a bootstrap algorithm, the so-called second-order wild bootstrap is applied. Consistency of the bootstrap procedure is proven without imposing strict stationary conditions or structural process assumptions, like linearity. The good finite sample performance of the second-order wild bootstrap is demonstrated by means of simulations.

We study a multivariate version of trend filtering, called Kronecker trend filtering or KTF, for the case in which the design points form a lattice in d dimensions. KTF is a natural extension of univariate trend filtering (Int. J. Comput. Vis. 70 (2006) 214–255; SIAM Rev. 51 (2009) 339–360; Ann. Statist. 42 (2014) 285–323), and is defined by minimizing a penalized least squares problem whose penalty term sums the absolute (higher-order) differences of the parameter to be estimated along each of the coordinate directions. The corresponding penalty operator can be written in terms of Kronecker products of univariate trend filtering penalty operators, hence the name Kronecker trend filtering. Equivalently, one can view KTF in terms of an ${\ell }_1$-penalized basis regression problem where the basis functions are tensor products of falling factorial functions, which is a piecewise polynomial (discrete spline) basis that underlies univariate trend filtering.

This paper is a unification and extension of the results in (In Advances in Neural Information Processing Systems (2016); in Advances in Neural Information Processing Systems (2017)). We develop a complete set of theoretical results that describe the behavior of kth-order Kronecker trend filtering in d dimensions, for every $\mathit{k}\ge 0$ and $\mathit{d}\ge 1$. This reveals a number of interesting phenomena, including the dominance of KTF over linear smoothers in estimating heterogeneously smooth functions, and a phase transition at $\mathit{d}=2(\mathit{k}\mathbf{+}1)$, a boundary past which (on the high dimension-to-smoothness side) linear smoothers fail to be consistent entirely. We also leverage recent results on discrete splines from (Tibshirani (2020)), in particular, discrete spline interpolation results that enable us to extend the KTF estimate to any off-lattice location in constant-time (independent of the size of the lattice n).

We study the problem of community recovery and detection in multi-layer stochastic block models, focusing on the critical network density threshold for consistent community structure inference. Using a prototypical two-block model, we reveal a computational barrier for such multilayer stochastic block models that does not exist for its single-layer counterpart: When there are no computational constraints, the density threshold depends linearly on the number of layers. However, when restricted to polynomial-time algorithms, the density threshold scales with the square root of the number of layers, assuming correctness of a low-degree polynomial hardness conjecture. Our results provide a nearly complete picture of the optimal inference in multiple-layer stochastic block models and partially settle the open question in (J. Amer. Statist. Assoc. 118 (2023) 2433–2445) regarding the optimality of the bias-adjusted spectral method.

This paper proposes a Gaussian process (GP) approach for testing conditional moment restrictions. Tests are based on squared Neyman orthogonal function-parametric processes integrated with respect to a GP distribution. This methodology leads to a general unified framework of kernel-based tests having the following properties: (i) bootstrap tests are easy to implement in the presence of nuisance parameters (they are simple quadratic forms, and there is no need to reestimate the nuisance parameters in each bootstrap replication); and (ii) the new tests are valid under general conditions, including higher-order conditional moments of unknown form, regularized estimators (e.g., Lasso) or parameters at the boundary of the parameter space. Novel applications include distance kernel tests for zero conditional treatment effects. The paper introduces Neyman orthogonal kernels, a new asymptotic theory and a detailed local power analysis. Monte Carlo experiments and a real data application illustrate the sensitivity of tests to the dimension of covariates and to the mean and covariance kernel of the GP.