The Hardness of Conditional Independence Testing and the Generalised Covariance Measure

Introduction

Conditional independences lie at the heart of several fundamental concepts such as sufficiency [Fisher20] and ancillarity [Fisher34, Fisher35]; see also Sorensen2017. Dawid79 states that “many results and theorems concerning these concepts are just applications of some simple general properties of conditional independence”. During the last few decades, conditional independence relations have played an increasingly important role in computational statistics, too, since they are the building blocks of graphical models [Lauritzen1996, Koller2009, Pearl2009].

Conditional independence tests also play a key role in causal inference. Constraint-based or independence-based methods [Pearl2009, Spirtes2000, Peters2017book] apply a series of conditional independence tests in order to learn causal structure from observational data. The recently introduced invariant prediction methodology [Peters2016jrssb, Heinze2017] aims to estimate for a specified target variable $Y$, the set of causal variables among potential covariates $X_{1},\ldots,X_{d}$. Given data from different environments labelled by a variable $E$, the method involves testing for each subset $S\subseteq\{1,\ldots,d\}$, the null hypothesis that the environment $E$ is conditionally independent of $Y$, given $X_{S}$.

Given the importance of conditional independence tests in modern statistics, there has been great deal of work devoted to developing conditional independence tests; we review some important examples of tests in Section 1.3. However one issue that conditional independence tests appear to suffer from is that they can fail to control the type I error in finite samples, which can have important consequences in downstream analyses.

As an example application of our results on the GCM, we consider the case where the regressions are performed using kernel ridge regression, and show that provided the conditional expectations are contained in a reproducing kernel Hilbert space, our test statistic has a tractable limit distribution.

The rest of the paper is organised as follows. In Sections 1.1 and 1.2, we first formalise the notion of conditional independence and relevant concepts related to statistical hypothesis testing. In Section 1.3 we review some popular conditional independence tests, after which we set out some notation used throughout the paper. In Section 2 we present our main result on the hardness of conditional independence testing. We introduce the generalised covariance measure in Section 3 first treating the univariate case with $d_{X}=d_{Y}=1$ before extending ideas to the potentially high-dimensional case. In Section 4 we apply the theory and methodology of the previous section to study that particular example of generalised covariance measures based on kernel ridge regression. We present numerical experiments in Section 5 and conclude with a discussion in Section 6. All proofs are deferred to the appendix and supplementary material.

This may be viewed as an extension of the fact that for one-dimensional $X$ and $Y$, the partial correlation coefficient $\rho_{X,Y|Z}$ (the correlation between residuals of linear regressions of $X$ on $Z$ and $Y$ on $Z$) is 0 if and only if X\mbox{{}\perp\mkern-11.0mu\perp{}}Y\mid Z in the case where $(X,Y,Z)$ are jointly Gaussian.

2 Statistical hypothesis testing and notation

is a measurable function whose last argument is reserved for a random variable $U\sim U$ independent of the data which is responsible for the randomness of the test.

Given a sequence of tests $(\psi_{n})_{n=1}^{\infty}$, the following validity properties will be of interest; note the particular names given to these properties differ in literature. Given a level $\alpha\in(0,1)$ and null hypothesis $\mathcal{P}$, we say that the test $\psi_{n}$ has

where the left-hand side is the size of the test; the sequence $(\psi_{n})_{n=1}^{\infty}$ has

In practice, we would like a test to have at least uniformly asymptotic level. Otherwise, even for an arbitrarily large sample size $n$, there can exist null distributions for which the size exceeds the nominal level by some fixed amount.

[lecam1973convergence, barron1989uniformly, Prop. 2, Thm. 3]. To achieve power tending to 1, we need to restrict $\mathcal{Q}$ by imposing certain smoothness conditions for example [balakrishnan2017hypothesis].

The hypothesis testing problem defined by the pair $(\mathcal{P},\mathcal{Q})$ is then said to be untestable [dufour2003identification]. In order to have power at even a single alternative, we need to restrict the null $\mathcal{P}$ in some way. One of the main results of this paper is that conditional independence is untestable.

3 Related work

Our hardness result for conditional independence contributes to an important literature on impossibility results in statistics and econometrics starting with the work of bahadur1956 which shows that there is no non-trivial test for whether a distribution has mean zero. canay2013 shows that certain problems arising in the context of identification of some nonparametric models are not testable. In these examples, the null hypothesis is dense with respect to the TV metric in the alternative hypothesis, a property which implies untestability [romano2004non]. Interestingly, our Proposition 5 shows that conditional independence testing is qualitatively different in that some distributions in the alternative are in fact well-separated from the null. It has been suggested for some time that conditional independence testing is a hard problem (see, e.g., [Bergsma2004], and several talks given by Bernhard Schölkopf). To the best of our knowledge the conjecture that conditional independence is not testable (cf. Corollary 3 with $M=\infty$) is due to Arthur Gretton. We also note that when the conditional distribution of $X$ given $Z$ is known, conditional independence is testable [e.g. berrett2018conditional].

We now briefly review several tests for conditional independence that bear some relation to our proposal here.

Ramsey2014 suggests regressing $X$ on $Z$ and $Y$ on $Z$ and then testing for independence between the residuals. Fan2015 consider this approach in the setting where $Z$ is potentially high-dimensional and under the null hypothesis of X\mbox{{}\perp\mkern-11.0mu\perp{}}Y\mid Z, $X=Z^{T}\beta_{X}+\varepsilon_{X}$, $Y=Z^{T}\beta_{Y}+\varepsilon_{Y}$ with \varepsilon_{X}\mbox{{}\perp\mkern-11.0mu\perp{}}Z and \varepsilon_{Y}\mbox{{}\perp\mkern-11.0mu\perp{}}Z. The following simple example however indicates where such methods can fail.

The Hilbert-Schmidt independence criterion (HSIC) equals the square of the Hilbert–Schmidt norm of the cross-covariance operator, and is used in unconditional independence testing [Gretton2008]. Fukumizu2008 extend this idea to conditional independence testing. To construct a test for continuous variables $Z$, their work requires clustering of the values of $Z$ and permuting $X$ and $Y$ values within the same cluster component. Another extension is proposed by Zhang2011uai. Their kernel conditional independence (KCI) test is stated to yield pointwise asymptotic level control.

4 Notation

No-free-lunch in Conditional Independence Testing

In this section we show that, under certain conditions, no non-trivial test for conditional independence with valid level exists. To state our result, we introduce the following subsets of $\mathcal{E}_{0}$ defined to be the set of all distributions for $(X,Y,Z)$ absolutely continuous with respect to Lebesgue measure.

A proof is given in the appendix. Note that taking $M$ to be finite ensures all the random vectors $(x_{i},y_{i},z_{i})$ are bounded. Thus, for example, averages will converge in distribution to Gaussian limits uniformly over $\mathcal{P}_{0,M}$; however, as the result shows, this does not help in the construction of a non-trivial test for conditional independence. An immediate corollary to Theorem 2 is that there is no non-trivial test for conditional independence with uniformly asymptotic level.

For all $M\in(0,\infty]$ and for any sequence $(\psi_{n})_{n=1}^{\infty}$ of tests we have

This result is in stark contrast to unconditional independence testing, where a permutation test can always be used to control the size of any testing procedure. As a consequence, there exist tests with valid level at sample size $n$ and non-trivial power. For example, Hoeffding1948 introduces a rank-based test in the case of univariate random variables and proves that it maintains uniformly asymptotic level and has asymptotic power against each fixed alternative. For the multivariate case, Berrett2017 consider a test based on mutual information and prove level guarantees, as well as uniform power results against a wide class of alternatives. Thus while independence testing remains a hard problem in that it is only possible to have uniform power against certain subsets of alternatives, this is different to conditional independence testing where we can only hope to control the size uniformly over certain subsets of the null hypothesis.

Inspection of the proof shows that Theorem 2 also holds in the case where the variables $X$ and $Y$ have marginal distributions that are absolutely continuous with respect to counting measure, for example. Theorem 2 therefore contains an impossibility result for testing the equality of two conditional distributions (by taking $Y$ to be an indicator specifying the distribution). The continuity of $Z$, however, is necessary. If $Z$ only takes values in $\{1,2\}$, for example, one can reduce the problem of conditional independence testing to unconditional independence testing by combining the tests for X\mbox{{}\perp\mkern-11.0mu\perp{}}Y\mid Z=1 and X\mbox{{}\perp\mkern-11.0mu\perp{}}Y\mid Z=2.

The null hypothesis being dense with respect to TV distance among the alternative hypothesis is a sufficient condition for the problem to be untestable [romano2004non]. Proposition 5, proved in the supplementary material, illustrates that this is not the case here: at least for $M\in(0,\infty)$, there exists an alternative, for which there is no distribution from the null that is arbitrarily close.

For $P,Q\in\mathcal{E}_{0}$, the total variation distance is given by

In Proposition 16 in the appendix, we also show that the null and alternative hypotheses are well-separated in the sense of KL divergence. On the other hand, it is known that if a problem is untestable, the convex closure of the null must contain the alternative [Kraft55, Bertanha2018, Theorem 5 and Corollary 1, respectively]. The problem of conditional independence testing therefore has the interesting property of the null being separated from the alternative, but its convex hull is TV-dense in the alternative.

A practical implication of the negative result of Theorem 2 is that domain knowledge is needed to select a conditional independence test appropriate for the data at hand. However guessing the form of the entire joint distribution in order to apply a test with the appropriate type I error control seems challenging. In Section 3 we introduce a form of test that instead relies on selecting regression methods that have sufficiently low prediction error when regressing $\mathbf{Y}^{(n)}$ and $\mathbf{X}^{(n)}$ on $\mathbf{Z}^{(n)}$, thereby converting the problem of finding an appropriate test to the more familiar task of prediction. Before discussing this methodology, we first sketch some of the main ideas of the proof of Theorem 2 below.

Figure 1 sketches the main components in our construction of $P$, which is laid out formally in Lemmas 13 and 14 in the appendix.

The key idea is as follows. Given $(X,Y,Z)\sim P$, we consider a binary expansion of $(X,Y,Z)$, which we truncate at some point to obtain $(\mathring{X},\mathring{Y},\mathring{Z})$. We then concatenate the digits of $\mathring{X}$ and $\mathring{Z}$ placing the former at the end of the binary expansion, thereby embedding $\mathring{X}$ within $\mathring{Z}$. This way, $\mathring{X}$ can be reconstructed from $\mathring{Z}$, and adding noise gives a distribution that is absolutely continuous with respect to Lebesgue measure. By making the truncation point sufficiently far down the expansions, we can ensure the $\epsilon$ proximity required.

The Generalised Covariance Measure

We have seen how conditional independence testing is not possible without restricting the null hypothesis. In this section we give a general construction for a conditional independence test based on regression procedures for regressing $\mathbf{Y}^{(n)}$ and $\mathbf{X}^{(n)}$ on $\mathbf{Z}^{(n)}$. In the case where $d_{X}=d_{Y}=1$, which we treat in the next section, the basic form of our test statistic is a normalised covariance between the residuals from these regressions. Because of this, we call our test statistic the generalised covariance measure (GCM). In Section 3.2 we show how to extend the approach to handle cases where more generally $d_{X},d_{Y}\geq 1$.

Given a distribution $P$ for $(X,Y,Z)$, we can always decompose

Let $\hat{f}^{(n)}$ and $\hat{g}^{(n)}$ be estimates of the conditional expectations $f_{P}$ and $g_{P}$ formed, for example, by regressing $\mathbf{X}^{(n)}$ and $\mathbf{Y}^{(n)}$ on $\mathbf{Z}^{(n)}$. For $i=1,\ldots,n$, we compute the product between residuals from the regressions:

Here, and in what follows, we have sometimes suppressed dependence on $n$ and $P$ for simplicity of presentation. We then define $T^{(n)}$ to be a normalised sum of the $R_{i}$’s:

Our final test can be based on $|T^{(n)}|$ with large values suggesting rejection. Note that the introduction of notation for the numerator and denominator in the definition of $T^{(n)}$ are for later use in Theorem 8.

In the case where $\hat{f}$ and $\hat{g}$ are formed through linear regressions, the test is similar to one based on partial correlation, and would be identical were the denominator in (3) to be replaced by the the product of the empirical standard deviations of the vectors $(x_{i}-\hat{f}(z_{i}))_{i=1}^{n}$ and $(y_{i}-\hat{g}(z_{i}))_{i=1}^{n}$. This approach however would fail for Example 1 despite $f$ and $g$ being linear (in fact both equal to the zero function) as the product of the variances of the residuals would not in general equal the variance of their product. Indeed, the reader may convince herself using pcor.test from the R package ppcor [ppcor], for example, that common tests for vanishing partial correlation do not yield the correct size in this case.

The following result gives conditions under which when the null hypothesis of conditional independence holds, we can expect the asymptotic distribution of $T^{(n)}$ to be a standard normal.

Applying the Cauchy–Schwarz inequality and Markov’s inequality, we see the requirement that $A_{f}A_{g}=o_{P}(n^{-1})$ is fulfilled if

Note that the rate of convergence requirement on $A_{f}$ and $A_{g}$ is a slower rate of convergence than the rate obtained when estimating Lipschitz regression functions when $d_{Z}=1$, for example. Furthermore, we show in Section 4 that $f$ and $g$ being in a reproducing kernel Hilbert space (RKHS) is enough for them to be estimable at the required rate.

In the setting where $Z$ is high-dimensional and $f$ and $g$ are sparse and linear, standard theory for the Lasso [tibshirani96regression, buhlmann2011statistics] shows that it may be used to obtain estimates $\hat{f}$ and $\hat{g}$ satisfying the required properties under appropriate sparsity conditions. In fact, in this case our test statistic is closely related to that involved in the ANT procedure of Ren2015 and the so-called RP test introduced in Shah2018, which amount to a regularised partial correlation. A difference is that the denominator in (3) means the GCM test would not require \varepsilon_{P}\mbox{{}\perp\mkern-11.0mu\perp{}}\xi_{P} unlike the ANT test and the RP test.

We now briefly sketch the reason for the relatively weak requirement on the MSPEs. In the following we suppress dependence on $P$ for simplicity of presentation. We have

The summands in the final term in (5) are i.i.d. with zero mean provided $P\in\mathcal{P}_{0}$, so the central limit theorem dictates that these converge to a standard normal. We also see that the simple form of the GCM gives rise to the term $b$ involving a product of bias-type terms from estimating $f$ and $g$, so each term is only required to converge to 0 at a slow rate such that their product is of smaller order than the variance of the final term. The summands in $\nu_{g}$ are, under the null, mean zero conditional on $(\mathbf{Y}^{(n)},\mathbf{Z}^{(n)})$. This term and similarly $\nu_{f}$ are therefore both relatively well-behaved, and give rise to the weak conditions on $B_{f}$ and $B_{g}$.

Control of the term $b$ in (5) under the alternative can proceed in exactly the same way as under the null. However control of the terms $\nu_{f}$ and $\nu_{g}$ typically requires additional conditions (for example Donsker-type conditions) on the estimators $\hat{f}$ and $\hat{g}$ as under the alternative both the errors $\varepsilon_{i}$ and $\hat{g}$ can depend on $\mathbf{Y}^{(n)}$. A notable exception is when $f$ and $g$ are sparse linear functions; in this setting alternative arguments can be used to show the GCM with Lasso regressions has optimal power when $Z$ has a sparse inverse covariance [Ren2015, Shah2018].

To state a general result avoiding additional conditions, here we will suppose that $\hat{f}$ and $\hat{g}$ have been constructed from an auxiliary training sample, independent of the data $(\mathbf{X}^{(n)},\mathbf{Y}^{(n)},\mathbf{Z}^{(n)})$ (e.g. through sample-splitting); see, for example, [robins2008higher, Zheng2011]. A drawback however, compared to the original GCM, is that the corresponding prediction error terms $A_{f}$ and $A_{g}$ are here out-of-sample prediction errors. These are typically more sensitive to the distribution of $Z$ and larger than the in-sample prediction errors featuring in Theorem 6. For this reason we consider the sample splitting approach to be more of a tool to facilitate theoretical analysis and would usually recommend using the original GCM in practice due to its typically better type I error control.

Then under the conditions of (i) in Theorem 6 we have

with $\tau_{N}^{(n)}$ and $\tau_{D}^{(n)}$ defined as in (3).

Under the conditions of (ii) in Theorem 6 we have

A proof is given in the supplementary material. We see that we achieve optimal $\sqrt{n}$ rates for estimating $\rho_{P}$.

1.2 Relationship to semiparametric models

In viewing $\tau_{N}^{(n)}/\sqrt{n}$ as an estimator of the functional $\rho_{P}$, our GCM test connects to a vast literature in semiparametric statistics. In particular, the requirement of estimating nonparametric quantities (in our case $\sqrt{A_{f}}$ and $\sqrt{A_{g}}$) at a $o(n^{-1/4})$ rate is common for estimators of functionals based on estimating equations involving influence functions [bickel1993efficient]. Our requirement on prediction error necessitates that at least one of $f_{P}$ and $g_{P}$ is Hölder $\beta$-smooth with $\beta/(2\beta+d_{Z})\geq 1/4$. Estimators of the expected conditional covariance functional requiring minimal possible smoothness conditions may be derived using the theory of higher order influence functions [robins2008higher, robins2009semiparametric, li2011higher, robins2017minimax]; these estimators are however significantly more complicated. newey2018cross study another approach to estimation of the functional based on a particular spline-based regression method. The work of chernozhukov2017double uses related ideas to ours here to obtain $1/\sqrt{n}$ convergent estimates and confidence intervals for parameters such as average treatment effects in causal inference settings. A distinguishing feature of our work here is that we only require in-sample prediction error bounds under the null of conditional independence, which is advantageous in our setting for the reasons mentioned in the previous section.

2 Multivariate X𝑋X and Y𝑌Y

We define our aggregated test statistic to be

In order to understand what values of $S_{n}$ indicate rejection, we will compare $S_{n}$ to

Here $\bar{\mathbf{R}}_{jk}$ is the sample mean of the components of $\mathbf{R}_{jk}$.

Let $\hat{G}$ be the quantile function of $\hat{S}_{n}$. This is a random function that depends on the data $(\mathbf{X}^{(n)},\mathbf{Y}^{(n)},\mathbf{Z}^{(n)})$ through $\hat{\boldsymbol{\Sigma}}$. Note that given the $\mathbf{R}_{jk}$, we can approximate $\hat{G}$ to any degree of accuracy via Monte Carlo.

The ground-breaking work of chernozhukov2013gaussian gives conditions under which $\hat{G}$ can well-approximate the quantile function of a version of $S_{n}$ where all bias terms, that is terms corresponding to $b$, $\nu_{g}$ and $\nu_{f}$ are all equal to 0. We will require that those conditions are met by $\varepsilon_{P,j}\xi_{P,k}$ for all $j=1,\ldots,d_{X}$, $k=1,\ldots,d_{Y}$. Below, we lay out these conditions, which take two possible forms. Let $d=\max(d_{X},d_{Y})$; note that $d$ and $\mathcal{P}$ are permitted to change with $n$, though we suppress this in the notation.

The result below shows that under the moment conditions above, provided the prediction error following the regressions goes to zero sufficiently fast, $\hat{G}$ closely approximates the quantile function of $S_{n}$ and therefore may be used to correctly calibrate our test.

Suppose that for $\mathcal{P}\subset\mathcal{P}_{0}$, the following is true: there are constants $C,c,c_{1}>0$ such that for each $n$ and $P\in\mathcal{P}$, there exists $C_{n}\geq 1$ such that one of (A1a) and (A1b) hold, and that (A2) holds. Suppose that

A proof is given in the supplementary material.

If the errors $\{\varepsilon_{P,j}\}_{j=1}^{d_{X}}$ and $\{\xi_{P,k}\}_{k=1}^{d_{Y}}$ are all sub-Gaussian with parameters bounded above by some constant $M$ uniformly across $P\in\mathcal{P}$, we may easily see that both (A1a) and (A1b) are satisfied with $C_{n}$ a constant; see chernozhukov2013gaussian for further discussion.

If additionally we have $A_{f,j},A_{g,k}=o_{\mathcal{P}}\left(\log(d)^{-1}\min\{n^{-1/2},\log(d)^{-4}\}\right)$, (6), (7) and (8) will all be satisfied.

GCM Based on Kernel Ridge Regression

We now apply the results of the previous section to a GCM based on estimating the conditional expectations via kernel ridge regression. For simplicity, we consider only the univariate case where $d_{X}=d_{Y}=1$. In the following, we make use of the notation introduced in Section 3.1.

Consider forming estimates $\hat{f}=\hat{f}^{(n)}$ and $\hat{g}=\hat{g}^{(n)}$ through kernel ridge regressions of $\mathbf{X}^{(n)}$ and $\mathbf{Y}^{(n)}$ on $\mathbf{Z}^{(n)}$ in the following way. For $\lambda>0$, let

We will consider selecting a final tuning parameter $\hat{\lambda}$ in the following data-dependent way:

The term minimised on the RHS is an upper bound on the mean-squared prediction error omitting constant factors depending on $\sigma^{2}$ (defined below in Theorem 11) and $\|f_{P}\|_{\mathcal{H}}^{2}$ or $\|g_{P}\|_{\mathcal{H}}^{2}$. Because of the hidden dependence on these quantities, this is not necessarily a practically effective way of selecting $\lambda$: our use of it here is simply to facilitate theoretical analysis. Finally define $\hat{f}=\hat{f}_{\hat{\lambda}}$, and define $\hat{g}$ analogously. We will write $T^{(n)}$ for the test statistic formed as in (3) with these choices of $\hat{f}$ and $\hat{g}$.

Let $\mathcal{P}$ be such that $u_{P}(z),v_{P}(z)\leq\sigma^{2}$ for all $z$ and $P\in\mathcal{P}$.

A proof is given in the supplementary material.

An application of the dominated convergence theorem shows that a sufficient condition for (10) to hold is that $\sum_{j=1}^{\infty}\sup_{P\in\mathcal{P}}\mu_{P,j}<\infty$.

Experiments

Section 3 proposes the generalised covariance measure (GCM). Although we provide detailed computations for kernel ridge regression in Section 4, the technique can be combined with any regression method. In practice, the choice may depend on external knowledge of the specific application the user has in mind. In this section, we study the empirical performance of the GCM with boosted regression trees as the regression method. In particular, we use the R package xgboost [chen2018xgboost, chen2016xgboost] with a ten-fold cross-validation scheme over the parameter maxdepth.

Theorem 2 states that if a conditional independence test has power against an alternative at a given sample size, then there is a distribution from the null that is rejected with probability larger than the significance level. We now illustrate the no-free-lunch theorem empirically.

Let us fix an RKHS $\mathcal{H}$ that corresponds to a Gaussian kernel with bandwidth $\sigma=1$. We now compute for different sample sizes the rejection rates for data sets generated from the following model: $Z=N_{Z}$, $Y=f_{a}(Z)+N_{Y}$, and $X=f_{a}(Z)+N_{X}$, with $N_{X},N_{Y},N_{Z}\sim\mathcal{N}(0,1)$, i.i.d., and $f_{a}(z):=\exp(-z^{2}/2)\sin(az)$ defining a function $f_{a}\in\mathcal{H}$. Figure 2 shows a plot of $f_{a}$ for $a=6$ and $a=18$.

Clearly, for any $a$, we have X\mbox{{}\perp\mkern-11.0mu\perp{}}Y\mid Z, but for large values of $a$ the independence will be harder to detect from data. We now fix three different sample sizes $n=100$, $n=1000$, and $n=10000$. For any of such sample size $n$, we can find an $a$, i.e., a distribution from the null, such that the probability of (falsely) rejecting X\mbox{{}\perp\mkern-11.0mu\perp{}}Y\mid Z is larger than the prespecified level $\alpha$. Figure 3 shows the results for the GCM test with boosted regression trees and the significance level $\alpha=0.05$: for any sample size, there exists a distribution from the null, for which the test rejects the null hypothesis of conditional independence.

For $n=100$, we can choose $a=6$, for $n=1000$, we choose $a=12$, and for $n=10000$, $a=18$.

is the Fourier transform of $f_{a}$. Equation (11) shows that a null hypothesis $\mathcal{P}$ containing all of the above models for $a>0$, violates one of the assumptions in Theorem 11: for this choice of RKHS and null hypothesis there is no $M$ such that $\sup_{P\in\mathcal{P}}\max(\|f_{P}\|_{\mathcal{H}},\|g_{P}\|_{\mathcal{H}})<M$. (Note that not all sequences of functions with growing RKHS norm also yield a violation of level guarantees: some functions with large RKHS norm, e.g., modifications of constant functions, can be easily learned from data.) Other conditional independence tests fail on the examples in Figure 3, too, for a similar reason. However most of these other methods are less transparent in the underlying assumptions, since they do not come with uniform level guarantees.

2 On Level and Power

It is of course impossible to provide an exhaustive simulation-based level and power analysis. We therefore concentrate on a small choice of distributions from the null and the alternative. In the following, we compare the GCM with three other conditional independence tests: KCI [Zhang2011uai] with its implementation from CondIndTests [Heinze2017], and the residual prediction test [Heinze2017, Shah2018]. We also compare to a test that performs the same regression as GCM, but then tests for independence between the residuals, rather than vanishing correlation, using HSIC [Gretton2008]. (This procedure is similar to the one that Fan2015 propose to use in the case of additive noise models.) As we discuss in Example 1, we do not expect this test to hold level in general. We then consider the following distributions from the null:

$Z\sim\mathcal{N}(0,1)$, $X=f_{a}(Z)+0.3\mathcal{N}(0,1)$, $Y=f_{a}(Z)+0.3\mathcal{N}(0,1)$, $a=2$;

$Z_{1},Z_{2}\sim\mathcal{N}(0,1)$ independent, $X=f_{1}(Z_{1})-f_{1}(Z_{2})+0.3\mathcal{N}(0,1)$, $Y=f_{1}(Z_{1})+f_{1}(Z_{2})+0.3\mathcal{N}(0,1)$;

$Z\sim\mathcal{N}(0,1)$, $X_{1}=f_{1}(Z)+0.3\mathcal{N}(0,1)$, $X_{2}=f_{1}(Z)+X_{1}+0.3\mathcal{N}(0,1)$, $Y_{1}=f_{1}(Z)+0.3\mathcal{N}(0,1)$, $Y_{2}=f_{1}(Z)+Y_{1}+0.3\mathcal{N}(0,1)$; and

$Z\sim\mathcal{N}(0,1)$, $Y=f_{2}(Z)\cdot\mathcal{N}(0,1)$, $X=f_{2}(Z)\cdot\mathcal{N}(0,1)$.

In the remainder of this section, we refer to these settings as (a) “$a=2$”, (b) “$a=4$”, (c) “biv. $Z$”, (d) “biv. $X,Y$”, and (e) “multipl. noise”, respectively. For each of the sample sizes $50$, $100$, $200$, $300$, and $400$, we first generate $100$ data sets, and then compute rejection rates of the considered conditional independence tests. The results are shown in Figure 4. For rejection rates below $0.11$ the hypothesis “the size of the test is less than $0.05$” is not rejected at level 0.01 (pointwise). The GCM indeed has promising behaviour in terms of type I error control. As expected, however, it requires the sample size to be big enough to obtain a reliable estimate for the conditional mean.

We then investigate the tests’ power by altering the data generating processes (a)–(e), described above. Each equation for $Y$ receives an additional term $+0.2X$, which yields X\mbox{{}\not\!\perp\mkern-11.0mu\perp{}}Y\mid Z (for (d), we add the term $+0.2X_{2}$ to the equation of $Y_{2}$). Figure 5 shows empirical rejection rates.

All methods, except for RPT, are able to correctly reject the hypothesis that the distribution is from the null, particularly with increasing sample size. In our experimental setup, it is the level analysis, that poses a greater challenge for the methods other than GCM.

Discussion

A key result of this paper is that conditional independence testing is hard: non-trivial tests that maintain valid level over the entire class of distributions satisfying conditional independence and that are absolutely continuous with respect to Lebesgue measure cannot exist. In unconditional independence testing, control of type I error is straightforward and research efforts have focussed on power properties of tests. Our result indicates that in conditional independence testing, the basic requirement of type I error control deserves further attention. We argue that as domain knowledge is necessary in order to select a conditional independence test appropriate for a particular setting, there is a need to develop conditional independence tests whose suitability is reasonably straightforward to judge.

In this work we have introduced the GCM framework to address this need. The ability for the GCM to maintain the correct level relies almost exclusively on the predictive properties of the regression procedures upon which it is based. Selecting a good regression procedure, whilst mathematically an equally impossible problem, can at least be usefully informed by domain knowledge. We hope to see further applications of GCM-based tests in the future. On the theoretical side, it would be interesting to understand more precisely the tradeoff between the type I and type II errors in conditional independence testing. Often, work on testing fixes a null and then considers what sorts of classes of alternative distributions it is possible, or impossible to maintain power against. In the context of conditional independence testing, the problem set is even richer, in that one must also consider subclasses of null distributions, and can then study power properties associated with that null.

Acknowledgements

We thank Kacper Chwialkowski, Kenji Fukumizu, Arthur Gretton, Bernhard Schölkopf, and Ilya Tolstikhin for helpful discussions, initiated by BS, on the hardness of conditional independence testing, and KC and AG for helpful comments on the manuscript. BS has raised the point that for finitely many data conditional independence testing may be arbitrarily hard in several of his talks, e.g., at the Machine Learning Summer School in Tübingen in 2013. We are very grateful to Matey Neykov, Anton Lundborg and Cyrill Scheidegger for kindly pointing out some errors in earlier versions of this manuscript, and suggesting potential fixes. We also thank Peter Bühlmann for helpful discussions regarding the aggregation of tests via taking the maximum test statistic. Finally, we thank four anonymous referees and an associate editor for helpful comments that have improved the manuscript.

Appendix A Proof of Theorem 2

Let $L=L(\eta)$ be as defined in Lemma 13 (taking $\delta=\eta$). From Lemma 15 applied to $\check{R}$, we know there exists a finite union $R^{\sharp}$ of hypercubes each of the form

such that $\mu(R^{\sharp}\triangle\check{R})<\eta/\max(L,M_{1}^{n})$. Now on the region $B_{M_{1}}^{c}$ defining $\Omega_{1}$ we know that the density of $(\mathbf{V},U)$ is bounded above by $M_{1}^{n}$. Thus we have that

Then since $R^{r}\uparrow R^{\sharp}$ as $r\downarrow 0$, there exists $r_{0}>0$ such that $\mu(R^{\sharp}\setminus R^{r_{0}})<\eta/M_{1}^{n}$.

A.2 Auxilliary Lemmas

Step 2: We can now apply Lemma 14 with $W=(\check{Y},\check{Z})$ and $N=\check{X}$. This gives us $K!$ random vectors $\mathring{V}^{(1)},\ldots,\mathring{V}^{(K!)}$ where $K>2^{r}>n$; for each $m=1,\ldots,K!$, $\mathring{V}^{(m)}=(\mathring{X}^{(m)},\mathring{Y}^{(m)},\mathring{Z}^{(m)})$ satisfies

$\mathring{X}^{(m)}$ may be recovered from $\mathring{Z}^{(m)}$ via $\mathring{X}^{(m)}=\mathring{g}_{m}(\mathring{Z}^{(m)})$ for some function $\mathring{g}_{m}$;

From (a’), by the triangle inequality we have that

For $j\in\{1,\ldots,K\}$ let $G_{j}:=\cup_{k=1}^{K}G_{jk}$. Let $\mathcal{J}$ be the set of $n$-tuples $(j_{1},\ldots,j_{n})$ of distinct elements of $\{1,\ldots,K\}$. Now define

and let $D=C\times$. Fix $(j_{l})_{l=1}^{n}\in\mathcal{J}$ and $(k_{l})_{l=1}^{n}\in\{1,\ldots,K\}^{n}$, and set $G:=\prod_{l=1}^{n}G_{j_{l}k_{l}}$. Then $G$ has non-empty intersection with a given $S_{m}^{n}$ if and only if $k_{l}=\pi_{m}(j_{l})$ for all $l$. Thus if $(k_{l})_{l=1}^{n}\notin\mathcal{J}$, $G$ is disjoint from all $S_{m}^{n}$. On the other hand if $(k_{l})_{l=1}^{n}\in\mathcal{J}$, the number of $S_{m}^{n}$ that intersect $G$ is $(K-n)!$, the number of permutations of $\{1,\ldots,K\}$ whose outputs are fixed at $n$ points. We therefore have that all but at most $(K-n)!$ of the support sets $T_{m}$ are disjoint from $G\times$, whence

Now the set $C$ is the disjoint union of all sets $\prod_{l=1}^{n}G_{j_{l}k_{l}}$ with $(j_{l})_{l=1}^{n}\in\mathcal{J}$ and $(k_{l})_{l=1}^{n}\in\{1,\ldots,K\}^{n}$. Thus summing over all such sets we obtain

This gives that there exists at least one $m=m^{*}$ with

for every $m$. Putting things together, we have that there must exist an $m^{*}$ with

for $K$ sufficiently large, which can be arranged by taking $r$ sufficiently large. ∎

Moreover, defining $\mathring{W}^{(m)}=(W_{1},\ldots,W_{l-1},\mathring{W}_{l}^{(m)})$,

the probability that $(N,\mathring{W}^{(m)})$ takes any value is bounded above by $K^{-1}2^{-(m+d)r}M$;

Define the random variable $\mathring{N}$ by

This is a concatenation of the binary expansions of $2^{r}N_{j}\in\{0,1,2,\ldots,2^{t}-1\}$ for $j=1,\ldots,d$. Observe that $\mathring{N}\in\{0,1,\ldots,K-1\}$ and that $N_{j}$ may be recovered from $\mathring{N}$ by examining its binary expansion. Indeed, $2^{r}N_{j}$ is the residue modulo $2^{t}$ of $\left\lfloor\mathring{N}/2^{t(j-1)}\right\rfloor$.

One important feature of this construction is that we can recover $E$ from $\mathring{N}_{m}$ (and $m$) via $\pi_{m}(E)=\left\lfloor\mathring{N}_{m}/K\right\rfloor$, and thereby determine $E$, which then reveals $\mathring{N}$ and each of the individual $N_{j}$. In summary, this gives us $K!$ different embeddings of the vector $N$ into a single random variable.

using the independence of $E$ in the second line above. This gives (iv).

The following well-known result appears for example in weaver2013measure.

such that $\mu(B\triangle B^{\sharp})\leq\epsilon$, where $\mu$ denotes Lebesgue measure and $\triangle$ denotes the symmetric difference operator.

Appendix B KL-Separation of null and alternative

References

Appendix C Proof of additional results in Sections 2 and B

In this section we provide proofs of Propositions 5 and 16.

The proof of Proposition 5 makes frequent use of the following lemma.

Let probability measures $P$ and $Q$ be defined on $[-M,M]$. Then

Let $\left\|Q-P\right\|_{\text{TV}}=\delta$. Then

where we have used Lemma 17 in the second and last lines above. We now apply Lemma 17 once more to give

C.1.2 Proof of Proposition 16

where $I_{Q}(X,Y)$ denotes the mutual information between $X$ and $Y$, and $\rho_{Q}$ is the correlation coefficient of the bivariate Gaussian $(X,Y)$, both under $Q$.

Straightforward calculations (Section 10.1.2 of Bishop2006) show that $p^{*}_{1}(x|z)=p^{*}_{1}(x)$ and $p^{*}_{2}(y|z)=p^{*}_{2}(y)$ are Gaussian densities with mean zero and variances $\Sigma_{11}:=(\Sigma^{-1}_{Q})_{11}^{-1}=\sigma^{2}/(\sigma^{2}+1)$ and $\Sigma_{22}:=(\Sigma^{-1}_{Q})_{22}^{-1}=\sigma^{2}$, respectively, where $\Sigma_{Q}$ is the covariance matrix of the bivariate distribution for $(X,Y)$ under $Q$. It then follows from (17) that $p_{3}^{*}(z)=q_{Z}(z)$. Thus we have,

Appendix D Proofs of Results in Section 3

In this section, we will use $a\lesssim b$ as shorthand for $a\leq Cb$ for some constant $C\geq 0$, where what $C$ is constant with respect to will be clear from the context.

We begin by proving (i). We shall suppress the dependence on $P$ and $n$ at times to lighten the notation. Recall the decomposition,

Thus the summands $\varepsilon_{i}\xi_{i}$ in the final term of (18) are i.i.d. mean zero with finite variance, so the central limit theorem dictates that this converges to a standard normal distribution. By the Cauchy–Schwarz inequality, we have

We now turn to $\nu_{f}$ and $\nu_{g}$. Conditional on $\mathbf{Y}$ and $\mathbf{Z}$, $\nu_{g}$ is a sum of mean-zero independent terms and

Using Slutsky’s lemma, we may conclude that $\tau_{N}\stackrel{{\scriptstyle d}}{{\to}}\mathcal{N}(0,1)$. We now argue that the denominator $\tau_{D}$ will converge to 1 in probability, which will give us $T_{n}\stackrel{{\scriptstyle d}}{{\to}}\mathcal{N}(0,1)$ again by Slutsky’s lemma.

First note that from the above we have in particular that $(b+\nu_{f}+\nu_{g})/n\stackrel{{\scriptstyle P}}{{\to}}0$. Thus $\sum_{i=1}^{n}R_{i}/n\stackrel{{\scriptstyle P}}{{\to}}0$ by the weak law of large numbers (WLLN). It suffices therefore to show that $\sum_{i=1}^{n}R_{i}^{2}/n\stackrel{{\scriptstyle P}}{{\to}}1$. Now

Multiplying out and using the inequality $2|ab|\leq a^{2}+b^{2}$ we have

by the same argument as used to show $\nu_{g}\stackrel{{\scriptstyle P}}{{\to}}0$. Similarly, we also have that the corresponding term involving $f$, $\hat{f}$ and $\xi_{i}^{2}$ tends to 0 in probability. For the final term in $\text{I}_{i}$, we have

The uniform result (ii) follows by an analogous argument to the above, the only differences being that all convergence in probability statements must be uniform, and the convergence in distribution via the central limit theorem must also be uniform over $\mathcal{P}$. These stronger properties follow easily from the stronger assumptions given in the statement of the result; that they suffice for uniform versions of the central limit theorem, WLLN and the particular applications of Slutsky’s lemma required here to hold is shown in Lemmas 18, 19 and 20 below.

D.2 Uniform convergence results

For each $n$, let $P_{n}\in\mathcal{P}$ satisfy

By the central limit theorem for triangular arrays [vaart_1998, Proposition 2.27], we have

thus taking limits in (22) immediately yields the result. ∎

Also, by Markov’s inequality and then the triangle inequality, we have for each $n$ that

using Hölder’s inequality. By Markov’s inequality, we have

We prove (a) first. Given $\epsilon>0$, let $N$ be such that for all $n\geq N$ and for all $P\in\mathcal{P}$

Thus for all $n\geq N$ and for all $P\in\mathcal{P}$,

Then for all $n\geq N$ and for all $P\in\mathcal{P}$, for $t\geq 0$

D.3 Proof of Theorem 8

The proof of this result is very similar to that of Theorem 6, and we will adopt the same notation here. We shall suppress the dependence on $P$ and $n$ at times to lighten the notation. We shall denote the auxiliary dataset by $\mathbf{A}$. We begin by proving (i). We have

Thus the summands $\varepsilon_{i}\xi_{i}-\rho$ in the final term are i.i.d. mean zero with finite variance, so the central limit theorem dictates that this converges to a standard normal distribution.

Control of the term $b$ is identical to that in the proof of Theorem 6. Turning to $\nu_{f}$ and $\nu_{g}$, Conditional on $\mathbf{Z}$ and the auxiliary dataset $\mathbf{A}$, $\nu_{g}$ is a sum of mean-zero independent terms and

That $\nu_{g}\stackrel{{\scriptstyle P}}{{\to}}0$ follows exactly as in the argument preceding (20), and similarly for $\nu_{f}$. Using Slutsky’s lemma, we may conclude that $\tau_{N}-\sqrt{n}\rho_{P}\stackrel{{\scriptstyle d}}{{\to}}\mathcal{N}(0,1)$.

The argument that $\tau_{D}\stackrel{{\scriptstyle P}}{{\to}}\sigma$ proceeds similarly to that in the proof of Theorem 6, but with conditioning on $\mathbf{X}$ or $\mathbf{Y}$ replaced by conditioning on $\mathbf{A}$. The uniform result (ii) follows by an analogous argument, see the comments at the end of the proof of Theorem 6.

D.4 Proof of Theorem 9

The proof of Theorem 9 relies heavily on results from chernozhukov2013gaussian which we state in the next section for convenience, after which we present the proof Theorem 9.

$C_{n}^{2}(\log(pn))^{7}/n\leq Cn^{-c}$ for some constants $C,c>0$.

The labels of the corresponding results in chernozhukov2013gaussian are given in brackets. A slight difference between our presentation of these results here and the statements in chernozhukov2013gaussian is that we consider the maximum absolute value rather than the maximum.

There exists an absolute constant $C^{\prime}>0$ such that for all $t\geq 0$,

Then there exists a constant $c^{\prime}>0$ such that

The following result includes a slight variant of Lemma 3.2 of chernozhukov2013gaussian whose proof follows in exactly the same way.

Writing $G_{\boldsymbol{\Sigma}}$ and $G_{\boldsymbol{\Theta}}$ for the quantile functions of $V$ and $\max_{j=1,\ldots,p}|U_{j}|$ respectively,

Note the second inequality does not appear in chernozhukov2013gaussian but follows easily in a similar manner to the first inequality.

D.4.2 Proof of Theorem 9

We will assume, without loss of generality, that ${\mathbf{var}}_{P}(\varepsilon_{P,j}\xi_{P,k})=1$ for all $P\in\mathcal{P}$. Furthermore, we will suppress dependence on $P$ and $n$ at times in order to lighten the notation. We will use $C^{\prime}$ to denote a positive constant that may change from line to line.

in terms of $\kappa_{P}$, and later bound $\kappa_{P}$ itself. Fixing $P\in\mathcal{P}$ and suppressing dependence on this, we have

Now let $\Omega$ be the event that $\max_{j,k}|\delta_{jk}|\leq u_{\delta}$ and $\max_{j,k}|\Delta_{jk}|\leq u_{\Delta}$.

From Theorem 22, we have $\text{I}=o(1)$. Lemma 23 and Lemma 21 give

We thus see that writing $a_{n}=\log(d)^{-2}$, if $\max_{j,k}|\delta_{jk}|=o_{\mathcal{P}}(a_{n}^{1/4})$, $\max_{j,k}|\Delta_{jk}|=o_{\mathcal{P}}(a_{n})$ and $\|\boldsymbol{\Sigma}-\hat{\boldsymbol{\Sigma}}\|_{\infty}=o_{\mathcal{P}}(a_{n})$, then we will have $\sup_{P\in\mathcal{P}}\sup_{\alpha\in(0,1)}v_{P}(\alpha)\to 0$. These remaining properties are shown in Lemma 26.

D.5 Auxiliary Lemmas

Consider the setup of Theorem 9 and its proof (Section D.4.2). Let $a_{n}=\log(d)^{-2}$. We have that

$\max_{j,k}|\delta_{jk}|=o_{\mathcal{P}}(a_{n}^{1/4})$;

$\max_{j,k}|\Delta_{jk}|=o_{\mathcal{P}}(a_{n})$;

$\|\boldsymbol{\Sigma}-\hat{\boldsymbol{\Sigma}}\|_{\infty}=o_{\mathcal{P}}(a_{n})$.

The arguments here are similar to those in the proof of Theorem 6, but with the added complication of requiring uniformity over expressions corresponding to different components of $X$ and $Y$. We will at times suppress the dependence of quantities on $P$ to lighten notation.

We begin by showing (i). Let us decompose each $\delta_{jk}$ as $\delta_{jk}=b_{jk}+\nu_{g,jk}+\nu_{f,jk}$, these terms being defined as the analogues of $b$, $\nu_{g}$ and $\nu_{f}$ but corresponding to the regression of $\mathbf{X}_{j}^{(n)}$ and $\mathbf{Y}_{j}^{(n)}$ on to $\mathbf{Z}^{(n)}$.

By the Cauchy–Schwarz inequality, we have $b_{jk}\leq\sqrt{n}A_{f,j}^{1/2}A_{g,k}^{1/2}=o_{\mathcal{P}}(a_{n}^{1/4})$ using (6). Let us write $\omega_{ik}=g_{k}(z_{i})-\hat{g}_{k}(z_{i})$. In order to control $\max_{j,k}|\nu_{g,jk}|$ we will use Lemma 29. Given $\epsilon>0$, we have

whence $\max_{j,k}|\nu_{g,jk}|=o_{\mathcal{P}}(a_{n}^{1/4})$. Similarly $\max_{j,k}|\nu_{f,jk}|=o_{\mathcal{P}}(a_{n}^{1/4})$, which completes the proof of (i).

Lemma 27 shows that the first term on the RHS is $o_{\mathcal{P}}(a_{n})$. For the second term we have

from (i) and Lemma 24, noting that (A2) implies in particular that $\log(d)^{3}=o(n)$. Thus applying Lemma 28, we have that $\max_{j,k}|\Delta_{jk}|=o_{\mathcal{P}}(a_{n})$.

We already know that $\max_{j,k}|\Delta_{jk}|=o_{\mathcal{P}}(a_{n})$ so applying Lemma 28, we see that $\max_{j,k}|(1+\Delta_{jk})^{-1}-1|=o_{\mathcal{P}}(a_{n})$. It is then straightforward to see that (26) holds. This completes the proof of (iii). ∎

Consider the setup of Theorem 9 and its proof (Section D.4.2) as well as that of Lemma 26. We have that

Fix $i$ and consider $R_{jk,i}R_{lm,i}-\varepsilon_{ij}\xi_{ik}\varepsilon_{il}\xi_{im}$. Writing $\eta_{j}=f_{j}(z_{i})-\hat{f}_{j}(z_{i})$ and $\omega_{k}=g_{k}(z_{i})-\hat{g}_{k}(z_{i})$, and suppressing dependence on $i$ (so e.g. $\varepsilon_{ij}=\varepsilon_{j}$) we have

We see that the sum on the RHS contains terms of four different types of which $\eta_{j}\omega_{k}\eta_{l}\omega_{m}$, $\eta_{j}\omega_{k}\eta_{l}\xi_{m}$, $\eta_{j}\omega_{k}\varepsilon_{l}\xi_{m}$ and $\eta_{j}\xi_{k}\varepsilon_{l}\xi_{m}$ are representative examples. We will control the sizes of each of these when summed up over $i$. Turning first to $\eta_{j}\omega_{k}\eta_{l}\omega_{m}$, note that $2|\eta_{j}\omega_{k}\eta_{l}\omega_{m}|\leq\eta_{j}^{2}\omega_{k}^{2}+\eta_{l}^{2}\omega_{m}^{2}$.

The argument of (19) combined with (6) shows that

Next we have $2|\eta_{j}\omega_{k}\eta_{l}\xi_{m}|\leq\eta_{j}^{2}\omega_{k}^{2}+\eta_{l}^{2}\xi_{m}^{2}$.

Arguing as in (25), we have for any $\epsilon>0$,

using Markov’s inequality in the final line. Next by (8), $\max_{i,m}\xi_{im}^{2}=O_{\mathcal{P}}(\tau_{f,n}^{2})$, so

Again by (8), this is $o_{\mathcal{P}}(1)$, so by bounded convergence, we have that

Considering the third term, we have $2|\eta_{j}\omega_{k}\varepsilon_{l}\xi_{m}|\leq\eta_{j}^{2}\xi_{m}^{2}+\omega_{k}^{2}\varepsilon_{l}^{2}$, so this term may be controlled in the same way.

using (27). This completes the proof of the result. ∎

Let $\epsilon,\delta>0$. As $f^{\prime}$ is continuous at 0, it is bounded on a sufficiently small interval $(-\delta^{\prime},\delta^{\prime})\subseteq D$. Let $M=\sup_{x\in(-\delta^{\prime},\delta^{\prime})}|f^{\prime}(x)|$ and set $\eta=\min(\delta^{\prime},\delta/M)$. Note by the mean-value theorem we have the inequality $|f(x)-c|\leq M|x|\leq\delta$ for all $x\in(-\eta,\eta)$. Thus

We will apply a symmetrisation argument to the inner conditional expectation. To this end, introduce $W^{\prime}$ such that $W^{\prime}$ and $W$ have the same distribution conditional on $V$ and such that W^{\prime}\mbox{{}\perp\mkern-11.0mu\perp{}}W\mid V. In addition, let $S_{1},\ldots,S_{n}$ be i.i.d. Rademacher random variables independent of all other quantities. The RHS of the last display is equal to

Appendix E Proof of Theorem 11

We will prove (ii) first. From Theorem 6 and Remark 7, it is enough to show that

and an analogous result for $\hat{g}$. We know from Lemma 30 that

for a constant $C>0$. Note that the first inequality in the last display allows us to effectively move from a fixed design with a design-dependent tuning parameter $\lambda$ to a random design but where $\lambda$ is fixed since the minimum is outside the expectation. For $P\in\mathcal{P}$, let $\phi_{P}:[0,\infty)\to[0,\infty)$ be given by

Observe that $\phi_{P}$ is increasing and $\lim_{\lambda\downarrow 0}\sup_{P\in\mathcal{P}}\phi_{P}(\lambda)=0$ by (10). Let $\lambda_{P,n}=n^{-1/2}\sqrt{\phi_{P}(n^{-1/2})}$ so $\sup_{P\in\mathcal{P}}\lambda_{P,n}=o(n^{-1/2})$. Thus for $n$ sufficiently large $\phi_{P}(\lambda_{P,n})\leq\phi_{P}(n^{-1/2})$, whence for such $n$ we have

To show (i), set $\mathcal{P}=\{P\}$ in the preceding argument and note that $\lim_{\lambda\downarrow 0}\phi_{P}(\lambda)=0$ by dominated convergence theorem using the summability of the eigenvalues i.e. (10) always holds.

The following result gives a bound on the prediction error of kernel ridge regression with fixed design. The arguments are similar to those used in the analysis of regular ridge regression, see for example JMLR:v14:dhillon13a.

Let $z_{1},\ldots,z_{n}\in\mathcal{Z}$ (for some input space $\mathcal{Z}$) be deterministic and suppose

Let $\mathbf{X}=(x_{1},\ldots,x_{n})^{T}$. We know from the representer theorem [Kimeldorf1970, Schoelkopf2001] that

Let $V=\text{span}\{k(\cdot,z_{1}),\ldots,k(\cdot,z_{n})\}\subseteq\mathcal{H}$ and write $f=u+v$ where $u\in V$ and $v\in V^{\perp}$. Then

where $\langle\cdot,\cdot\rangle$ denotes the inner product of $\mathcal{H}$. Write $u=\sum_{i=1}^{n}\alpha_{i}k(\cdot,z_{i})$. Then

where $K_{i}$ is the $i$th column (or row) of $K$. Thus K\alpha=\Big{(}f(z_{1}),\ldots,f(z_{n})\Big{)}^{T}. By Pythagoras’ theorem

Now let the eigendecomposition of $K$ be given by $K=UDU^{T}$ with $D_{ii}=\hat{\mu}_{i}$ and define $\theta=U^{T}K\alpha$. We see that $n$ times the left-hand side of (29) is

Now as $\theta=DU^{T}\alpha$ note that $\theta_{i}=0$ when $d_{i}=0$. Let $D^{+}$ be the diagonal matrix with $i$th diagonal entry equal to $D_{ii}^{-1}$ if $D_{ii}>0$ and 0 otherwise. Then

using the inequality $(a+b)^{2}\geq 4ab$ in the final line. Finally note that

Putting things together gives the result. ∎

Consider the setup of Theorem 11. For all $r>0$,

The proof below is not included in the published version of this paper, where instead some results in a book are cited in order to establish this result. However, it was subsequently discovered that the results in the book were incorrect. Below is a complete proof of the result.

By the Cauchy–Schwarz inequality, for all $z$ and $z^{\prime}$,

so $K-\Phi\Phi^{T}/n$ is positive semi-definite.

By Weyl’s inequality, noting that the nonzero eigenvalues of $\Phi^{T}\Phi$ and $\Phi\Phi^{T}$ coincide, we have, for all $i$,

where $(\cdot)_{+}:=\max(\cdot,0)$ denotes the positive part. This will prove concavity of $f$ as

so subtracting (31) yields $f(tA+(1-t)B)\geq tf(A)+(1-t)f(B)$ as desired.

Certainly (31) holds when $r\geq\lambda_{1}(tA+(1-t)B)$. Now by Lidskii’s inequality, for each $j=1,\ldots,d$,

For convenience, let us set $\lambda_{d+1}(tA+(1-t)B)=0$. Then for any $j=1,\ldots,d$, if $\lambda_{j+1}(tA+(1-t)B)\leq r\leq\lambda_{j}(tA+(1-t)B)$, we have

using (32) for the first inequality. We thus have that (31) holds whatever the value of $r$, and so $f$ is concave, which completes the proof. ∎

References