A Number of computationally efficient statistical meth- ods have been proposed for analyzing massive data sets. Examples include divide-and-conquer approaches [1]- [4]; low-rank approximations: random projection methods [5]- [8], subsampling methods [9]- [11], Nyström approximations [12], [13]; and online learning methods [14]- [16].
An interesting question arising from these new methods is the minimum computational cost required for obtaining statistically satisfactory solutions. This might be viewed as a type of "computational limit" from a statistical perspective. Such an issue has been addressed in certain situations. For divide-and-conquer approaches, [4] derived a sharp upper bound for the number of distributed computing units in the smoothing spline setup, while [17] estimated the quantile regression process under an additional sharp lower bound on the number of quantile levels. For random projection methods, the literature nonetheless only focused on parametric cases such as compressed sensing. For example, [18] showed that the minimum number of random projections is s log n for signal recovery, where n is the number of measurements and s is the number of nonzero components in the true signal. Recently, [7] proposed the randomly sketched kernel ridge regression (KRR) estimator and studied the minimax optimal nonparametric estimation under random projection. To our knowledge, the computational limit for random projection methods remains unknown in nonparametric models.
There are two purposes in this paper: (i) develop an optimal nonparametric testing procedure based on random projection; (ii) explore its computational limit in the kernel ridge regression setup. We remark that classical nonparametric testing methods, e.g., the locally most powerful test, the generalized/penalized likelihood ratio test and the distance-based test [19]- [23], may not be directly applied to big data due to their high computational costs.
Specifically, we consider the following nonparametric model
where x i ∈ X ⊆ R a for a fixed a ≥ 1 are i.i.d. random design points, and i are random noise following Normal distribution with mean zero, variance σfoot_0 . The regression function f belongs to a reproducing kernel Hilbert space (RKHS) H. The hypothesis of interest is
where f 0 is a hypothesized function. Testing optimality has been well studied in literature. [24] [25] established the minimax testing rate for Gaussian sequence model on the Sobolev class or the Besov class. Recently, [26] established the minimax nonparametric testing rate under a general eigen-decaying framework including the polynomial decay and exponential decay kernels. In practice, testing in (2) has wide applications. One motivating example is in the signal detection in cognitive radio and other wireless applications, it is assumed under the null hypothesis that the signal is completely specified, e.g., that no signal is present. Another example is testing the adequacy of a parametric linear model in nonparametric regression, where the null hypothesis assumes f 0 has a linear structure; see [27], [28].
Focusing on the testing in (2), we construct a distancebased test statistic T n,λ = f R -f 0 compared with the classic KRR estimator f n ( [29]) defined as
where f 2 H = f, f H with •, • H the inner product of H, λ > 0 is a smoothing parameter. Solving ( 3) is an ndimensional quadratic program, which involves the computational cost and storage occupation of orders O(n 3 ) and O(n 2 ), respectively. Instead, f R is achieved by randomly projecting the row and column subspaces of the ndimensional kernel matrix to an s-dimensional subspace, reducing (3) to an s-dimensional quadratic program, with time and storage costs as O(s 3 ) and O(s 2 ) under s( n) random projections; see Section 2 for detailed algorithm. The pre-processing step in computing the kernel approximation normally takes O(sn 2 ), and can be easily reduced to O(n 2 (log s)) for suitably chosen random matrices (see [30]), which can be further reduced to O(n 2 (log s)/t) by using t clusters in a parallel fashion. After f R is obtained, T n,λ can be computed in a parallel fashion. Hence, s can be viewed as a simple proxy for computing and storage costs.
In this paper, we reveal a phase transition phenomenon in terms of s to guarantee the testing optimality. Specifically, a sharp lower bound for s is established: when s is above this bound, T n,λ is minimax optimal; otherwise, minimax optimality becomes impossible even when the best possible λ is chosen. We next illustrate more subtle details using the following Figure 1, where the strength of the weakest detectable signals (SWDS) is characterized given any s and λ. In general, we require s > s λ for any λ, where s λ is determined by kernel eigenvalues and λ. An important observation is that the smallest SWDS can be achieved at λ = λ * and s > s λ * := s * (note that when s s * , our testing procedure under a proper λ is still powerful as long as SWDS becomes sufficiently large). Both λ * and s * have precise orders in specific situations. For example, in an morder polynomial decay kernel, the smallest SWDS achieves the minimax optimal rate n -2m 4m+1 established in [24], [25] when λ * = n -4m 4m+1 and s * = n 2 4m+1 . As a by-product, we also derive a sharp lower bound for s for obtaining the minimax optimal estimation. Our results hold for a general class of random projection matrix, such as the sub-Gaussian matrix or certain data-dependent matrix.
It is worth mentioning that the construction of T n,λ crucially relies on the regularity of H, which is often unavailable in practice. Hence, we propose an adaptive test statistic based on the maximum of a sequence of (standardized) nonadaptive test statistics corresponding to various regularities. Based on a recent Gaussian approximation result in [31], we prove that the null limit distribution is an extreme value distribution.
The proofs of main results rely on the behavior of the tail sum of empirical kernel eigenvalues. One technical contribution of this work is to derive upper bounds for a range of tail sums such that nonparametric estimation and testing can now be analyzed in a unified framework. This is obtained by flexibly adjusting the size of the function class associated with the Rademacher average in the local Rademacher complexity theory ( [32]). In simulation studies, we find that the size and power of the proposed non-adaptive and adaptive test statistics are both satisfactory. In particular, the power cannot be further improved as the number of random projections grows beyond some threshold, as predicted by our theory. For an illustration purpose, we also demonstrate that when n = 2 12 , conducting testing based on f R only takes 3.2 seconds in comparison with 42 seconds based on f n . In practice, the smoothing parameter λ can be directly selected via generalized cross validation. We would like to point out that this is an advantage of the random projection method over the divide-and-conquer method [1] in estimation, where the selection of the smoothing parameter is nontrivial; see [33].
The rest of this paper is organized as follows. Section 2 introduces kernel ridge regression together with its approximation based on random projection. Our main results are presented in Section 3: Section 3.1 introduces one primary assumption on random projection; Sections 3.2 and 3.3 study testing consistency and power behaviors in terms of the projection dimension s and the smoothing parameter λ, with specific situations considered in Section 3.5; Section 3.6 proves the lower bound on s given in Section 3.5 to be sharp. An adaptive testing procedure is developed in Section 4. 2 . For a matrix A ∈ R m×n , its operator norm is defined as
A random variable X is said to be sub-Gaussian if there exists a constant σ 2 > 0 such that for any t ≥ 0, P[|X| ≥ t] ≤ 2 exp(-t 2 /(2σ 2 )). The sub-Gaussian norm of X is defined as X ψ2 = inf{t > 0 : E exp(X 2 /t 2 ) ≤ 2}. We will use c, c 1 , c 2 , C to denote generic absolute constants, whose values may vary from line to line.
In this section, we review kernel ridge regression and its variant based on random projection. Suppose that we have n i.i.d. observations {(x i , y i )} n i=1 from (1), where the covariances x i are samples of X from a distribution P X with domain X . Let X be a closed subset of R a , a is fixed,P X a strictly positive Borel measure on X . We recall that a Borel measure P X on X is said to be strictly positive if the measure of every nonempty open subset in X is positive, an example being the Lebesgue measure in R a .
Throughout assume that f ∈ H, where H ⊂ L 2 (P X ) is a reproducing kernel Hilbert space (RKHS) associated with an inner product •, • H and a reproducing kernel function K(•, •) : X × X → R. K is a symmetric positive definite kernel on X satisfying: for any finite set of points {x i } n i=1 in X and real numbers {a i } n i=1 that n i,j=1 a i a j K(x i , x j ) ≥ 0. Assume further that K is a continuous function on X × X and X X K(x, t)dP X (x)dP X (t) < ∞. Then by Mercer's theorem, K has the following spectral expansion:
where µ 1 ≥ µ 2 ≥ • • • ≥ 0 is a sequence of ordered eigenvalues and the eigenfunctions {φ i } ∞ i=1 form a basis in L 2 (P X ). We refer the reader to the standard sources [34], [35] for more details on RKHSs and their properties.
Moreover, for any i, j ∈ N,
Throughout this paper, assume that φ i 's are uniformly bounded, a common condition in literature, e.g., [36], and µ i 's satisfy certain tail sum property.
Assumption A1 is satisfied in two types of commonly used kernels, categorized by the eigenvalue decay rates. The first is µ i i -2m for a constant m > 0, called as polynomial decay kernel (PDK) of order m. Examples of kernels in this class include the m th order Sobolev spaces for some fixed integer m ≥ 1 with Lebesgue measure on a bounded domain; see [35].
The second is µ i exp(-γi p ) for constants γ, p > 0, called as exponential decay kernel (EDK) of order p. Examples of EDK include the Gaussian kernel, which for the Lebesgue measure satisfies such a bound with p = 1 (compact domain) or p = 2 (real line); see [37]. Verification of Assumption A1 with concrete examples is deferred to Section S.5.3 in Supplementary.
Recall the KRR estimator f n from (3). By representer theorem, it has an expression
where ω = ( ω 1 , . . . , ω n ) is a real vector determined by ω = argmin
,j≤n , and I ∈ R n×n is identity. This standard procedure requires storing (K 2 , K, Ky) and inverting K + λI, which requires O(n 2 ) memory usage and O(n 3 ) floating operations.
The above computational and storage constraints become severe for a large sample size, and thus motivate the random projection approach proposed by [7]. Specifically, ω in ( 5) is substituted with S β, where β ∈ R s and S is an s × n real-valued random matrix; see Section 3.1. Then, β is solved as:
Hence, the resulting estimator of f becomes
which requires computing and storing (SK 2 S , SKS , SKy), along with inverting an s × s matrix. The cost in the pre-processing step to compute the kernel approximation normally takes O(sn 2 ), and can be easily reduced to O(n 2 (log s)) for suitably chosen random matrices (see [30]), which can be further reduced to O(n 2 (log s)/t) by using t clusters in a parallel fashion. Furthermore, the memory usage and floating operations are reduced to O(s 2 ) and O(s 3 ), respectively, when s = o(n). On the other hand, s cannot be too small in order to maintain sufficient data information for achieving statistical optimality. Critical lower bounds for s will be derived in Section 3.6.
Consider the nonparametric testing problem (2). For convenience, assume f 0 = 0, i.e., we will test
In general, testing f = f 0 (for an arbitrary known (8) has no loss of generality. Based on f R , we propose the following distancebased test statistic:
In the subsequent sections, we will derive the null limit distribution of T n,λ (Theorems 3.3), and further provide a sufficient and necessary condition in terms of s such that T n,λ is minimax optimal (Section 3.6). As a byproduct, we derive a critical bound in terms of s such that f R is minimax optimal. Proof of such results rely on an exact analysis on the kernel and projection matrices which requires an accurate estimate of the tail sum of the empirical eigenvalues by Lemma 3.1. Our results hold for a general choice of projection matrix which will be discussed in Section 3.
Consider the singular value decomposition K = U DU , where U = (U 1 , U 2 ) with U 1 consisting of the first s λ columns of U and U 2 consisting of the rest n -s λ columns,
Here { µ i } n i=1 are empirical eigenvalues in decreasing order. For any λ > 0, s λ (or s λ ) is defined to be the number of µ i 's (or µ i 's) greater than λ, i.e.,
We have the following assumption on the population eigenvalues through s λ . Assumption A2. s λ diverges as λ → 0.
Assumption A2 is satisfied in various classes of kernels, including PDK and EDK introduced in Section 3.5. In the Supplementary S.5.3, we verify Assumption A2 with two concrete examples:
An accurate upper bound for the tail sum of empirical eigenvalues n i= s λ+1 µ i is needed for studying nonparametric testing and estimation. However, this bound was often assumed to hold in the kernel learning literature, e.g., [38], [39]. The application of concentration inequalities of individual eigenvalues ( [40], [41]) only provides a very loose bound due to accumulative errors. Recently, the local Rademacher complexity theory ( [32]) was employed by [7] to derive a more accurate upper bound that is useful in studying nonparametric estimation. However, this upper bound no longer works for testing problems, due to the improper size of the function class defining Rademacher average. We establish the upper bounds, i.e., Lemma 3.1, for a range of tail sums of empirical eigenvalues in terms of population quantities s λ and µ s λ , with known orders. This result can be applied to both nonparametric estimation and testing, and may be of independent interest. Lemma 3.1. If 1/n < λ → 0, then with probability at least 1 -4e -s λ , n i= s λ+1 µ i s λ µ s λ . Clearly, Lemma 3.1 is a sample analog to the tail sum assumption for µ i in Assumption A1. The proof of Lemma 3.1 is based on the development of a generalized function class and its associated local Rademacher complexity theory as explained in Appendix 7.1.
The following definition of "K-satisfiability" describes a class of matrices that preserve the principal components of the kernel matrix. Definition 1. (K-satisfiability) A matrix S ∈ R s×n is said to be K-satisfiable if there exists a constant c > 0 such that
By Definition 1, a K-satisfiable S will make (SU 1 ) SU 1 "nearly" identity as well as down-weight the tail eigenvalues. Such a matrix will be able to extract the principle information from the kernel matrix. A special case of the above "K-satisfiability" condition was studied in [7] by fixing λ as the optimal estimation rate. However, by choosing a range of λ as threshold to select the leading eigenvalues, our general form of "K-satisfiability" condition allows us to study estimation and testing in a unified framework.
Besides, we need the following definition which will make the statement of our assumptions more precise and concise. Definition 2. An event E is said to be of (a, b)-type for a, b ∈ (0, ∞],
Definition 2 describes events whose probabilities have exponential type lower bounds. It is easy to see that, if E is of (a, b)-type, then P(E)
In particular, E is of (∞, ∞)-type if and only if E occurs almost surely.
Throughout the rest of this paper, assume the following condition on S. Assumption A3. (a) s ≥ qs λ for a sufficiently large constant q > 0. (b) There exist c 1 , c 2 ∈ (0, ∞] such that the event "S is Ksatisfiable" is of (c 1 s, c 2 s λ )-type.
Assumption A3 (a) requires a sufficient amount of random projections to preserve data information. Assumption A3 (b) requires S to be K-satisfiable with high probability which holds in a broad range of situations such as matrix of sub-Gaussian entries (Example 3) and certain data dependent matrix (Example 4). Example 3. Let S be an s × n random matrix of entries S ij / √ s, i = 1, . . . , s, j = 1, . . . , n, where S ij are independent (not necessarily identically distributed) sub-Gaussian variables. Examples of such sub-Gaussian variables include Gaussian variables, bounded variables such as Bernoulli, multinomial, uniform, variables with strongly log-concave density (see [42]), or mixtures of sub-Gaussian variables. The following lemma shows that Assumption A3 (b) holds in all these situations. Lemma 3.2. Let S ij : 1 ≤ i ≤ s, 1 ≤ j ≤ n be independent sub-Gaussian of mean zero and variance one, and λ ∈ (1/n, 1). If s ≥ qs λ for a sufficiently large constant q, then Assumption A3 (b) holds for The proof of Lemma 3.2 relies on the bound of the tail sums of empirical eigenvalues in Lemma 3.1. We point out that obtaining the eigen-decomposition in Example 4 is as burdensome as computing the matrix inverse, which is not preferred in practice. Rather, the purpose of this example is to directly illustrate one situation that Assumption A3 can be satisfied.
In this section, we derive the null limit distribution of (standardized) T n,λ as standard Gaussian, and then extend our result to the case of composite hypothesis testing. Theorem 3.3. Suppose that λ → 0 and s → ∞ as n → ∞.
Suppose Assumption A2 is satisfied. Then under H 0 , we have
Here,
Theorem 3.3 holds once s diverges (no matter how slowly). Theorem 3.3 implies the following testing rule at significance level α:
where z 1-α/2 is the 100 × (1 -α/2)th percentile of N (0, 1).
As an important consequence of Theorem 3.3, we comment that the optimal estimation rate in [7] can also be obtained as a by-product, that is
where r 2 n,λ = λ + µ n,λ . The proof of ( 12) is sketched as follows. Suppose that f 0 ∈ H is the "true" function in (1).
where E is the expectation w.r.t. = ( 1 , • • • , n ) with i the random noise in the regression function defined in (1). By direct examinations, it can be shown that
This completes the proof of (12). The above discussions are summarized in the following corollary.
Corollary 3.4. Suppose that 1/n < λ < 1 and Assumption A1-A3 holds. Then with probability approaching one, it holds that
where r 2 n,λ = λ + µ n,λ and C is an absolute constant. From Corollary 3.4, the best upper bound can be obtained through balancing λ and µ n,λ . Denote λ † the optimizer. This in turn provides a lower bound s † for s according to (10), i.e., s † = s λ † . In Section 3.5, we will show that the upper bound under λ † is minimax optimal, and further provide explicit orders for s † in concrete settings.
In this section, we investigate the power of T n,λ under a sequence of local alternatives by assuming f ∈ B = {f ∈ H : f H ≤ C} for some existing constant C.
For any generic 0-1 valued testing rule φ where φ = 0 if H 0 is preferred and 1 otherwise, define the total error Err(φ, d n ) of φ under a separation rate d n > 0 as
Notice E(φ | H 0 is true) is the probability of making a type I error and E(1 -φ | H 1 is true) is the probability of making a type II error, and the total error represents the maximum possible type I error and type II error. The separation rate d n is used to measure the distance between the null and the alternative hypotheses. Intuitively, the smaller d n is, the harder it is to distinguish the alternative hypothesis from the null. For any ε ∈ (0, 1), define the minimax separation rate d * n (ε) as
where the infimum in (15) is taken over all 0-1 valued testing rules based on samples ((x 1 , y 1 ), . . . , (x n , y n )). d * n (ε) characterizes the smallest separation between the null and local alternatives such that there exists a testing approach with a total error of at most ε. [24] [25] established the minimax separation rate and revealed its difference with optimal estimation rate. Their work are derived based on Gaussian sequence model with focus on the Sobolev class or the Besov class with polynomial decaying eigenvalues. Recently, [26] established the minimax nonparametric testing rate under a general eigen-decaying framework including the polynomial decay and exponential decay kernels. Next, we show that our proposed Wald-type test can achieve the minimax separation rate under appropriate λ and s.
For any f ∈ H, define the squared separation rate
In the following Theorem 3.5, we show that T n,λ can achieve high power provided that s diverges fast enough and the local alternative is separated from the null by at least an amount of d n,λ . It is sufficient to minimize the separation rate d n,λ to achieve optimal testing. We show that our test can achieve the minimax rate of testing by selecting λ to balance the trade-off between the bias of f R and the standard derivation of T n,λ shown in (16). Theorem 3.5. Suppose that 1/n < λ → 0 as n → ∞, Assumption A1-A2 are satisfied, and Assumption A3 holds for c 1 , c 2 ∈ (0, ∞]. Then for any ε > 0, there exist positive constants C ε and N ε such that, with probability greater than 1 -e -c1s -e -c2s λ ,
where d n,λ := λ + σ n,λ and B = {f ∈ H : f H ≤ C} for some existing constant C and P f (•|x, S) is the conditional probability measure under f given x, S.
In view of Theorem 3.5, to maximize the power of T n,λ , one needs to minimize d n,λ = λ + σ n,λ through balancing λ and σ n,λ . Denote λ * the optimizer. The lower bound s * for s is obtained via (10), i.e., s * = s λ * . The explicit forms of λ * and s * varies for different reproducing kernels, and lead to specific optimal testing rate, depending on their eigendecay rate.
In practice, it is often of interest to test certain structure of f , e.g., linearity,
where L(X ) is the class of linear functions over X ⊆ R a . Testing H linear 0 can be easily converted into simple hypothesis testing problem. Intuitively, if the parametric structure is right, then its residuals should be patternless and independent of input features. So we can apply non-parametric smoothing to the parametric residuals and see if their pattern is approximately zero everywhere; see Section 4.2 in [25].
Denote f * n as the least square estimator of f , f R as the randomly projected KRR estimator, and f 0 is a hypothesized "true" parameter with unknown value. Write
It can be shown that T Suppose Assumption A2 is satisfied. Then under H linear 0 , we have
Here, µ * n,λ := E H0 {T * n,λ |x, S} = tr(∆ 2 )/n, σ * 2 n,λ := Var H0 {T * n,λ |x, S} = 2tr(∆ 4 )/n 2 To characterize the power of the composite hypothesis testing problem in (17), we define the test error under a separation rate d n > 0,
where φ * is any desion rule for hypothesis testing problem in (17) and P L (X )f is the projection of f on L(X ). Notice that f -P L (X )f = 0 under the null hypothesis; the magnitude of ||f -P L (X )f || 2 n characterize how far the f is deviated from a linear function. Since the plugin estimate f * n approaches P L(X ) f with 1/n rate, the separation rate for the decision rule for φ
) is the same as φ n,λ given in Theorem 3.5.
Next, we derive the lower bounds for s to achieve optimal estimation and testing in two featured examples: PDK and EDK, based on the main results obtained in Corollary 3.4 and Theorem 3.5. It is easy to check that Assumption A1 and A2 hold for these two examples; see Section S.5.3 in Supplementary.
Theorem 3.7. For the two kinds of eigenvalue decaying rates, suppose Assumption A3 holds. Suppose that 1/n < λ → 0 as n → ∞, then with probability approaches 1, it holds that µ n,λ s λ /n and σ 2 n,λ
Furthermore, we have the following optimal estimation and testing rates by properly choosing the tuning parameters and the lower bound of projection dimension:
• Polynomially decaying kernel (with µ i i -2m ) -When λ n -2m 2m+1 and s n
4m+1 and s n 2 4m+1 with m > 3/2, T n,λ achieves the minimax optimal rate of testing n -2m 4m+1 .
• Exponentionally decaying kernel (with µ i exp(-γi p ))
1 , T n,λ achieves the minimax optimal rate of testing n -1 2 (log n)
.
The derivation of optimal estimation and testing rates is attributed to the accurate characterization of µ n,λ and σ n,λ by employing Lemma 3.1 to bound the tail sums of empirical eigenvalues.
Plugging in µ n,λ s λ /n to Corollary 3.4, we have f R with the convergence rate r 2 n,λ = λ + s λ /n. As shown in (13), λ is the squared bias of f R , and s λ /n quantifies the variance of f R . Hence, the optimal estimation rate r †2 n,λ is achieved via the bias-variance tradeoff as follows
We denote the choice of λ and lower bound of s to achieve the optimal estimation rate as λ † and s † respectively. To find the lower bound for s in achieving optimal testing, by Theorem 3.5, d 2 n,λ λ + √ s λ /n, then the optimal separation rate d * n can be achieved by another type of tradeoff, i.e., the squared bias of f R v.s. the standard derivation of T n,λ , as follows
We use λ * and s * to represent the optimal λ and s λ to achieve d * n 2 . Take PDK as an example, plugging in µ i i -2m to the definition of s λ in (10), we have s λ λ -1 2m . Then Theorem 3.7 can be directly achieved based on the above two types of tradeoffs in estimation and testing. The results for EDK can be achieved similarly. We conclude our findings of this section in the following Table 1.
It is worth emphasizing that λ † , s † are different from λ * , s * due to different types of trade-off discussed above, indicating a fundamental difference between estimation and testing ( [24], [25]). Figure 2 summarizes the two different types of trade-off to achieve the minimax rate in estimation and testing.
In the following Theorem 3.8, we show that the upper bound between f R and f n can fall below the statistical error f n -f 0 n by further increasing the projection dimension.
1. In fact, for EDK, s * (log n -1 2p log log n) 1/p . For simplicity, we keep the main term s * (log n) This result improves upon existing approximation error bound in [7] that is of the same order as the statistic error. When s > s † , the difference between f R and f n is ignorable compared with the difference between f n and f 0 . Theorem 3.8. Let δ > 0 satisfying δ ≤ λ. Define s = argmin{j : µ j < δ }. Suppose s ≥ cs , then with probability approaching 1,
Furthermore, when δ λ and λ → 0 as n → ∞,
In this section, we will show that s * and s † derived in PDK and EDK are actually sharp. For technical convenience, define
Our first result is about the sharpness of s † . Theorem 3.9 shows that when s s † , there exists a true function f such that f R -f 2 n is substantially slower than the optimal estimation rate. Our proof is constructive in the sense we construct the above true function as n i=1 K(x i , •)w i with w i being selected from the orthogonal complement of a subsapce properly generated by S and K; see Appendix 7.2 for details. Theorem 3.9. Suppose s = o(s † ). Then for any s × n random matrix S satisfying Assumption A3, with probability greater than 1 -e -cnδn , it holds that
where c is a constant independent of n.
Our second result is about the sharpness of s * . Theorem 3.10 shows that when s s * , there exists a local alternative f that is not detectable by T n,λ even when it is separated from zero by d * n . In this case, the asymptotic testing power is actually smaller than α. The proof of Theorem 3.10 is similar as that of Theorem 3.9, except that a different true function is constructed; see Appendix 7.3 for details.
where c is a constant independent of n. Recall 1 -α is the significance level.
In view of Theorems 3.5 and 3.10, we observe a subtle phase transition phenomenon for testing signals as shown in Figure 1.
In this section, we focus on the case of PDK as a leading example, and construct an adaptive testing procedure that does not require any exact prior knowledge on m except for m ≥ 2. The adaptive procedure is proven to achieve the minimax rate of testing established by [44] (up to an iterative-logarithmic term).
Consider an RKHS generated by a PDK of order m * ≥ 2, i.e., H = H m * . To reflect the role of m, we modify all previous notation by adding a subscript m. For example, let K m (•, •) be the reproducing kernel function associated with H m , and K m = 1 n [K m (x i , x j )] 1≤i,j≤n be the corresponding empirical kernel matrix. Let S m be an s m × n projection matrix. We will construct the corresponding f R,m (•) based on (7) under S m and λ m . Here
and the corresponding projection dimension s m is an integer satisfying
where q > 0 is a sufficiently large constant.
Given each m, the sketched KRR estimator f R,m = ∆ m y, where
Denote m n (log n) d0 for a constant d 0 ∈ (0, 1/2). Based on T n,m , our adaptive testing procedure is constructed as follows.
Step 1. For any 2 ≤ m ≤ m n → ∞, standardize T n,m as
.
Step 2. Calculate τ * n = max 1≤m≤mn τ m .
Step 3. Find
By allowing m n → ∞, the unknown m * will be eventually covered over a sequence of test statistics. Under the null hypothesis (8), T n,m = 1 n Δ 2 m , and thus τ m is of a standardized quadratic form. Then, τ * n is the maxima of a sequence of dependent τ m 's. Based on a recent Gaussian approximation result in [31], i.e., Lemma S.4 (stated in Supplementary), we prove in the following Theorem 4.1 that the null limit distribution of τ n,mn is some extreme value distribution.
Theorem 4.1. Suppose that m n (log n) d0 for a constant d 0 ∈ (0, 1/2), and, for 2 ≤ m ≤ m n , S m satisfies Assumption A3 (b) with projection dimension s m . Then, under H 0 in (8), for any α ∈ (0, 1), it holds that
Our next result states that the above adaptive testing procedure is asymptotically minimax optimal. Specifically, Theorem 4.2 shows that τ n,mn achieves high power if the local alternative is separated from zero by an order δ(n, m * ) defined as
And, [44] showed that δ(n, m * ) is minimax optimal rate for adaptive testing.
Theorem 4.2. Suppose that m n (log n) d0 for a constant d 0 ∈ (0, 1/2), and S m satisfies Assumption A3 (b) with projection dimension s m . Then, for any ε > 0, there exist positive constants C ε , N ε for any n ≥ N ε , with probability approaching 1,
In the end, we point out that the lower bound for s m given in (20) is slightly smaller than the sharp lower bound for s derived in the non-adaptive case; see Table 1. This is not surprising since the corresponding minimax rate δ(n, m * ), i.e., (23), is larger than the non-adaptive rate, i.e., n -2m * /(4m * +1) .
In this section, we examine the performance of the proposed testing procedure through simulation studies in Sections 5.1 and 5. ∼ N (0, 1) and c is a constant. To fit the model, we consider the periodic Sobolev kernel with eigenvalues μ 2i = μ 2i-1 = (2πi) -2m for i ≥ 1; see [46] for details. Set n = 2 9 , 2 10 , 2 11 , 2 12 , and H 0 : f = 0. The significance level was chosen as 0.05 and the Gaussian random projection matrix was applied in this setting.
We examined the empirical performance of the distancebased test (DT) T n,λ , and adaptive test (AT) τ n,mn . For DT, the projection dimension s was chosen as 2n γ for γ = 1/(4m+1), 2/(4m+1), 3/(4m+1), with m = 2 corresponding to cubic splines. For AT, the projection dimensions s m was chosen as 2n γ (log log n)
We propose a data-adaptive smoothing parameter slection rule to guarantee testing optimality. Specifically, we choose the optimal smoothing parameter λ * satisfying
where σ n,λ is defined in Theorem 3.3 as σ n,λ = 2trace(Δ 4 )/n 2 with Δ = KS (SK 2 S + λSKS ) -1 SK as a projection version of the classical smoothing matrix.
Empirical size was evaluated at c = 0, and power was evaluated at c = 0.01, 0.02, 0.03. Both size and power were calculated based on 500 independent replications. Figure 3 shows that the size of both DT and AT approach the nominal level 0.05 under various choices of (s, n), demonstrating the validity of the proposed testing procedure. Figure 4 displays the power of DT and AT. Under various choices of c and γ, it is not surprising to see from Figure 4 (a), (c), and (e) that the power of DT approaches one as n or c increases. Rather, a key observation is that the power cannot be further improved as γ grows beyond the critical point 2/(4m + 1) when c ≥ 0.02. This is consistent with our theoretical result; see Theorem 3.10. Similar patterns have been observed for the power of AT in Figure 4 (b), (d), and (f ). Of course, the power of AT is usually lower than that of DT under the same setup, especially when the signal strength is weak. This is the price paid for adaptivity.
In this section, we consider a multivariate case and test H 0 : f = 0. Data were generated from where (x i1 , x i2 , x i3 ) follows from N (μ, I 3 ) with μ = (0, 0, 0), i ∼ N (0, 1), and c ∈ {0, 0.05, 0.1, 0.15}. Specifically, we chose the Gaussian kernel
We considered sample sizes n = 2 9 to n = 2 12 and sketch dimensions s = 1.2 log(n), 1.2(log n) 3/2 , 1.2(log n) 2 . For each pair (n, s), experiments were independently repeated 500 times for calculating the size and power. Interpretations for Figure 5 about the size and power are similar to those for Figures 3 and4. Interestingly, we observe that the power increases dramatically as γ increases from 1 to 1.5, while becomes stable near one as γ ≥ 1.5. This is consistent with Theorem 3.7. Figure 6 demonstrates the significant computational advantage of DT (corresponding to γ < 1) over the testing procedure based on standard KRR (corresponding to γ = 1).
In the supplementary, we conduct additional synthetic experiments under the same simulation setup as Sections 5.1 and 5.2 except for using the Bernoulli random matrix. As shown in Figures S.1-S.3, the interpretations remain the same.
The main contribution of this paper is to apply random projection to nonparametric testing, and propose the first "sharp" result regarding projection dimension to guarantee minimax optimal testing respectively for a general class of random projections. In practice, many other random projections also satisfy K-satisfiability, for example, the randomized orthogonal system (ROS) sketches introduced in [47], [7], by randomly sampling and rescaling the rows of a fixed orthonormal matrix. For the ROS examples, Assumption A3 also holds with extra logarithm factor in projection dimension, see Lemma S.8 in supplementary. Stochastic approximation is another computationally efficient method for nonparametric learning. A representative approach is the stochastic gradient descent (SGD) algorithm. In SGD, the total step size within n steps iteration plays the role of the regularization to avoid overfitting, and the SGD estimator ( [16]) can achieve the optimality provided that the total step size has the same order of the effective dimension which is represented by s λ in our work. [48] established an early stopping rule for gradient descent algorithm to achieve optimal nonparametric testing. The result shows that the total step size in gradient descent plays the same role as 1/λ in classic KRR, and the optimal testing rate can be achieved by the same "bias-standard deviation" tradeoff while the bias and the standard deviation are represented as functions of the total step sizes. In sketched KRR, we showed that the projection dimension is required to be greater than the effective dimension s λ to guarantee high power performance. Therefore, although the sketched KRR and SGD are two different computationally efficient approaches, they are connected through the statistical effective dimension.
In this section, we present the main proofs of Lemma 3.1, sharpness properties including Theorem 3.9, Theorem 3.10.
Proofs for the rest of Lemmas and Theorems can be found in the Supplementary.
Before the proof of Lemma 3.1, we first state some definitions and preliminary lemmas. Define κ λ = s λ nλ , for any λ > 0. In fact, κ λ is the variance-to-bias ratio, where λ and s λ /n correspond to (squared-)bias and variance of f R , respectively; see Corollary 3.4 and Lemma S.3 for details. Consider a bundle of function classes indexed by κ λ :
Define a generalized version of local Rademacher complexity function
where σ 1 , . . . , σ n are independent Rademacher random variables, i.e., P(σ i = 1) = P(
When κ λ 1, Ψ λ (•) and Ψ λ (•) become the original local Rademacher complexity functions introduced in [32]. Note that κ λ 1 actually corresponds to the optimal bias vs. variance trade-off required for estimation. Rather, a different type of trade-off is needed for optimal testing as revealed by [22], [24], which corresponds to a different choice of κ λ in F λ as demonstrated in Section 3.
In the following Lemma 5 we represent the generalized local Rademacher complexity function and its empirical version by a function of eigenvalues and κ λ . In Lemma 6, we further show that both Ψ λ and Ψ λ possess unique (positive) fixed points. This fixed point property is crucial in proving Lemma 3.1. We defer the proof of Lemma 5 and Lemma 6 to Section S.5.2 in the Supplementary.
(a) Suppose µ 1 > 1/n. For any λ > 1/n, it holds that
(b) For any λ > 0, it holds that
Lemma 6. There exist uniquely positive r λ and r λ such that Ψ λ (r λ ) = r λ and Ψ λ ( r λ ) = r λ . Furthermore, if λ > 1/n, then r λ s λ /n, and there exists an absolute constant c > 0 such that, with probability at least 1 -e -cs λ , r λ s λ /n.
We are ready to prove Lemma 3.1.
Proof Plugging these fixed points into (26) and (27) in Lemma 5, we have
By Lemma 6, we have r λ s λ /n, leading to r λ /κ λ λ; for the empirical version, with probability at least 1e -cs λ , r λ s λ /n leading to r λ /κ λ λ. Recall that s λ = argmin{i : µ i ≤ λ} -1. Then by (29), with probability at least 1 -e -cs λ , where α = (SK 2 S + λSKS ) -1 SKy. Since the vector (SU D 1/2 ) α belongs to the column space of D 1/2 U S ∈ R n×s , and dim(D 1/2 U S ) = s. Then
By choosing β * orthogonal to the span of U S , so that D 1/2 β * is orthogonal to D 1/2 U S , and then we have A 3 = 0. Applying the minimax characterization of eigenvalues of the n × n matrix D 1/2 , we have with probability greater than 1 -e -cnδn , when s s † ,
where c is a constant. Proof Without loss of generality, here we consider H 0 : f = f 0 with f 0 = 0. We need to prove that when projection dimension s is too small, for any random matrix S, there exists true function f , such that f H ≤ C, f n ≥ d * , our testing rule still cannot detect it. We show the existence of such true function f as follows. We construct the true f (•) = where S = SU . Denote ∆ = (D S( SD 2 S + λ SD S) -1 SD 3/2 ) , then T 1 = n ∆ D 1/2 β * 2 2 . Notice that the span{∆} ⊂ span{D 3/2 S }, 3/2 S ) = s. There always exists β * ∈ R n×1 such that D 1/2 β * is orthogonal to the column space of D 3/2 S , so that ∆ D 1/2 β * = 0 leading to T 1 = 0. Furthermore, based on the above argument, we have T 3 = 0. Then we have T n,λ -µ n,λ σ n,λ = T 2 /n -µ n,λ σ n,λ d -→ N (0, 1).
On the other hand, the signal of such f satisfies
n , where f R is a randomly sketched KRR estimator, and • n is the empirical norm. The sketched KRR estimator f R ([7]) enjoys both theoretical support and computational efficiency, especially Authorized licensed use limited to: New Jersey Institute of Technology. Downloaded on August 01,2021 at 18:34:42 UTC from IEEE Xplore. Restrictions apply.
Authorized licensed use limited to: New Jersey Institute of Technology. Downloaded on August 01,2021 at 18:34:42 UTC from IEEE Xplore. Restrictions apply.