Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Authors: Yong Wu
Yanwei Fu
Shouyan Wang
Xinwei Sun
Venue: ICML 2025 poster (2025)
Keywords: causal discovery
unmeasured confounding
proxy variables
integral equation
causal hypothesis testing
negative control
Primary area: general_machine_learning->causality
TL;DR: We propose a proxy-based method for integral solving to identify causal relations when unobserved variables are present.

Abstract: Bivariate causal discovery is challenging when unmeasured confounders exist. To adjust for the bias, previous methods employed the proxy variable (*i.e.*, negative control outcome (NCO)) to test the treatment-outcome relationship through integral equations -- and assumed that violation of this equation indicates the causal relationship. Upon this, they could establish asymptotic properties for causal hypothesis testing. However, these methods either relied on parametric assumptions or required discretizing continuous variables, which may lead to information loss. Moreover, it is unclear when this underlying integral-related assumption holds, making it difficult to justify the utility in practice. To address these problems, we first consider the scenario where only NCO is available. We propose a novel non-parametric procedure, which enjoys asymptotic properties and preserves more information. Moreover, we find that when NCO affects the outcome, the above integral-related assumption may not hold, rendering the causal relation unidentifiable. Informed by this, we further consider the scenario when the negative control exposure (NCE) is also available. In this scenario, we construct another integral restriction aided by this proxy, which can discover causation when NCO affects the outcome. We demonstrate these findings and the effectiveness of our proposals through comprehensive numerical studies.

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Yong Wu 1234
Yanwei Fu Shouyan Wang 123467
Xinwei Sun E
Abstract
1. Introduction
When the causal sufficiency is violated (Spirtes et al., 2001), constraint-based causal discovery is challenging due to the potential bias introduced by unmeasured confounders in conditional independence testing. As foundational step, this paper focuses the bivariate casc, aiming identify causal relationships in the presence of unobserved variables. To adjust for the bias, (Miao et al., 2018) leveraged observed proxy variables (Kuroki & Pearl, 2014; Tehetgen et al., 2024) to examine the existence of causal relations over discrete variables. Later, (Miao et al., 2023; Liu et al., 2023) extended them continuous variables, by investigating the existence of solutions to integral equations.
Bivariate causal discovery is challenging when unmeasured confounders exist. To adjust for the bias, previous methods employed the proxy vari- able (i.e., negative control outcome (NCO)) to test the treatment-outcomme relationship through inte- gral equations - and assumed that violation of this equation indicates the causal relationship. Upon this, they could establish asymptotic properties for causal hypothesis testing. However, these meth- ods either relied on parametric assumptions or required discretizing continuous variables, which may lead to information loss. Moreover, it is unclear when this underlying integral-related as- sumption holds, making difficult justify the utility in practice. To address these problems, we first consider the scenario where only NCO is available. We propose a novel non-pann panaamamttrrcc pro- cedure, which enjoys asymptotic properties and preserves more information. Moreover, we find that when NCO affects the outcome, the above integral-related assumption may not hold, render- ing the causal relation unidentifiable. Informed by this, we further consider the scenario when the negative control exposure (NCE) is also available. In this scenario, we construct another integral re- striction aided by this proxy, which can discover causation when NCO affects the outcome, We demonstrate these findings and the effectiveness of our proposals through comprehensive numeri- cal studies.
Specifically, the integral equation examines whether the re- lationship between the treatment and the outcome can be fully explained by the proxy, a.k.a, negative control outcome (NCO) (Lipsitch et al., 2010). If there no causal relation, the relationship remains consistent across changes in treat- ment. Upon this, (Miao et al., 2023) tested the causal null hypothesis by goodness-of-fit. which relied on parametric assumptions. To allow for nonparametric testing, (Liu et al., 2023) proposed a discretization approach to approximate the equation, and established asymptotic properties. However, may require a large number of bins and substantial data to control the approximation error.
To address these issues when only the NCO is available, we propose a non-parametric testing procedure based on kernel estimator called Proxy Maximum Characteristic Restriction (PMCR), to examine the integral equation. With- out discretization, our approach preserves more information from the original data. After computing the equare residues via PMCR, construct the statisties and estab- lish its asymptotic properties. Compared other first-order moment restriction methods (Mastouri et al., 2021), our procedure leverages the characteristic function to capture all order moments, offering greater power in identifying the causal relation.
'Institute of Science and Technology for Brain- Brain-Inspired In- telligence, Fudan University 'Key Laboratory of Computational Neuroscience and Insinn-nspired Intelligence (Fudan University), Ministry of Education SSttee Key Laboratory of Medical Neurobi- ology and MOE Frontiers Center for Brain Science, Fudan Univer- sity "Zhangiiang Fudan International Innovation Center Schooo of Data Science, Fudan University "Shanghai Engineering Research Center of AI Robotics, Fudan University 'Engineering Research Center of Al Robotics, Ministry of Education, Fudan University, Correspondence to: Xinwei Sun <sunxinwei@ fww.a.cdu..com edu.en">.
Nevertheless, we proceed to note that, regarding the power analysis, both our method and others, are built upon basic assumption that the integral equation does not hold when the causal relation exists. However, by studying the solvability of the integral equation, we surprisingly find that this as- sumption may not hold, if the NCO's effect on the outcome
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
is strong enough to account for the effect of the treatment, making it fail to discover the causal relation with only the NCO. This inspires to leverage the additional negative control exposure (NCE). that was commonly used in causal inference (Miao et al., 2018; Tchetgen et al., 2024). In this scenario, we introduce another integral equation aided by such NCE, which can effectively discriminate the alternative hypothesis from the null, thereby enabling causal identifica- tion when NCO has strong effect the outcome.
3. Preliminary
Problem setup. We consider the problem of testing whether the causal relation X + Y exists, under the unmeasured confounder U. Under Markovian and faithfulness condi- tions, this is equivalent to testing the causal null hypothesis Ho X 1 YJU. To adjust for the confounding bias, we assume the availability of a proxy variable W such that X T W (Kuroki & Pearl, 2014), which also serves as the negative control outcome in causal inference. Fig. I (a) shows the causal diagram over X; Y. U, W. Besides, in some scenarios, we may have access to an additional proxy variable Z (i.e., negative control exposure), which satisfies Z T (W, Y),,,, X) illustrated in Fig. 1 (b). Through- out, we assume X, Y, U, W, Z are continuous variables.
To demonstrate our findings and the effectiveness of our approach, we test it in two scenarios-one in which the NCO influences the outcome and one in which it does not. We find that the procedure using only NCO performs well when the NCO has no effect on the outcome, but fails when it does influence the outcome. In this case, our procedure that additionally leverages the NCE can successfully discover the causal relationship.
Notations. Suppose X, Y, U, W, Z are continuous random variables over the probability space (D, IP). with domains X, 14 W, Z, respectively. For any variable U, we denote [F)) as the space of square-integrable functions with respect to the cumulative distribution function F(u). For any space W, let kw be its positive semi-definite kernel. We denote éw as its associated canonical feature map, i.e., ów(w) == kw (w. .) for any w E W. Besides, we de- note Hw as the corresponding reproducing kernel Hilbert space (RKHS). For any operator A : Hw - Hx. we denote A' as its adjoint operator. For any discrete vari- ables X, Y with respectively i categories, we denote P(y = {p(ylz1), .... p(y) x2)). the probability matrix P(Y[X) = {p(M)IX) ..... [リح x) T
We summarize our contributions as follows:
1. We propose a non-parametric procedure to solve the integral equation, which is more efficient and effective in learning causal relations when only NCO is available.
2. We study the solvability of the integral equation and find that it may not be able to identify the causal relation when the NCO's effect on the outcome is strong enough.
3. When NCE is available, we leverage it to construct an- other integral restriction, which can recover causal rela- tionships when NCO strongly affects the outcome.
4. We thoroughly verify our findings, utility, and effective- ness of our proposed methods on synthetic data.
U
W
Z
W
2. Related works
?
To adjust for the bias caused by unmeasured confounding, (Miao et al., 2018) leveraged proxy variables (Kuroki & Pearl, 2014; Tchetgen et al., 2024) to test the causal null hypothesis over discrete variables. Later, (Miao et al., 2023; Liu et al., 2023) extended their procedure to continuous variables by examining the integral equation. While their methods established asymptotic properties with only sin- gle proxy, their methods either relied parametric assump- tions or the discretization process that may lead to infor- mation loss. To address these issues, we propose a novel con-parametric procedure to examine the integral equation without discretization, hence is more sample efficient in learning causal relationships when only the NCO is avail- able. Additionally, by comprehensively studying the soly- ability of the integral equation, we surprisingly find that previous methods may not be able to identify the causal relation when NCO influences the outcome. When the NCE is available, we further construct another integral equation that can discriminate the alternative hypothesis from the null, enabling causal identification when NCO has strong influence on the outcome.
X
X
(a)
Figure 1. Causal diagrams over X, Y, U, W, Z. W (resp. Z) de- note the negative control outcome (resp. exposure). The dotted line indicates its potential presence or absence.
Previous studies and limitations. For discrete variables, (Miao et al., 2018) proposed to test Ho by examining whether PPWI x) can fully explain the variability of P(yl Z, x) as æ varies. Specifically, given X, Z, W with respectively j k categories, they performed a linear re- gression of 93 == (P(y) Z, x1 P(y/ x;)} on Q := (P(W) Z, 01 ",,,,),,,,,P, x;))), and tested the linearity based on the least-square residues. If Q is the full-column rank when ij k. they derived the null-limiting distribution of the statisties based on these residues.
Later, (Liu et al., 2023) noticed that the rank constraint satisfies when > k. allowing to test the linearity between P(y and P(WIX) without Z. Inspired by this, they extended to continuous variables by discretizing X, Y, W
2

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
and performed the linearity testing. While it could establish asymptotic properties for hypothesis testing, the method was constrained by the need for both large number of bins and substantial sample size to ensure accurate probability estimation. Similarly, (Miao et al., 2023) used goodness-of- fit to examine the following integral equation:
square-integrable function g. we assume ((g(w) 0 almost surely if and only if g(u) 0 almost surely.
Completeness 4.1 is a standard assumption in causal hy- pothesis testing (Miao et al., 2018; 2023; Liu et al., 2023). This condition is widely applicable, as shown by examples provided in (Newey & Powell, 2003; D' H'Hallffouille, 2011; Hu & Shiu, 2018; Andrews, 2017). Here, it means W car- ries all the variability of U. When W, U are discrete with i.j categories (i > j), it reduces to the full-column rank of P(WIU), as used in (Liu et al., 2023) for identification.
p(ulz)= = f hww, p))(uc (ar/://))dd dre for some hiw,, y);
which serves as the theoretical foundation of the discretiza- tion approach in (Liu et al., 2023). However, its goodness- of-fit procedures relied on the parametric assumption.
Proposition 4.2. Under assumption 4.1 and some regularity conditions in I.I, there exists g(w, u) L² (F(())) for all u, such that it solves the following integral equation for all z):
Does (I) holds under H1, ? To discover the causal relation- ship, these methods assumed that (1) does not hold under H1- However, it is unclear when this condition holds, making it hard to justify the validity during practical implementations.
p(ul xz) = ff g(w,, u)p(wlp).
Overview of our contributions. To address problems of existing works when only Win is available, Sec. 4 propose a novel non- non- con-paramttric approach based on Proxy Maxi- mum Characteristic Restriction (PMCR), which can test (I) without discretization, thereby avoiding information loss in causal discovery.
Remark 4.3. When W. U are discrete, this result reduces to PPU X) = ([W]U)--P(W' (X) to establish linear rela- tion under Hg (Liu et al., 2023).
By p(ul=) = J (u))u) under Ho. h(w,y) y) == I g(w, u)p(ylu)du solves (I). To further ensure h(w, y) is square integrable, we require the bounded likelihood ratio condition (Kato et al., 2021) for p(uly)/p(u), which holds except for the extreme dependency between U and Y.
Further, through comprehensively studying the solvability of (1) in Sec. 5, we find that when W + Y (the dotted edge in Fig. (a)), the above integral may have a solution even under HH, making all procedures fail to discover the causal relation using only W. This inspires to consider the scenario where the proxy Z is available, as commonly adopted in causal inference (Tchetgen et al., 2024). In this scenario, Sec. 6 leveraged 2 to construct another integral equation to examine the solution's 's property that involves U. Thanks to this auxiliary information, we can discover the causal relation that may not be identifiable with only W.
Assumption 4.4 (Bounded likelihood ratio). There exists C >> 0 such that 0 < p(uly)
< C almost surely for and y. p(u)
Under these conditions. Theorem 4.5 establishes the exis- tence of a solution to the integral equation in (1).
Theorem 4.5. Suppose assumptions 4.1, 4.4 and the reg- ularity condition B.I hold. If g(w, u) in (2) satisfies Jf almma)) (pp(w) ((l))p())ddd < co, then under Ho,, there exists h(w, y) L² (F(w)) for all y. such that it solves the integral equation in (1) for all (æ.
4. Hypothesis testing with a single proxy
Remark 4.6. Condition g(W, UVM,UU)pp))) E (()
We first consider the scenario when only the proxy W (i.e. NCO) is available. To test the causal null hypothesis, we propose to examine whether the integral equation (1) exists. To this end, Sec. 4.1 first shows that under Hl, the solution exists under some completeness conditions. To solve the equation, Sec. 4.2 introduces a novel estimation method and constructs the testing statistics. Sec. 4.3 establishes its asymptotic level and power for such statistics. Finally, Sec. 4.4 summarizes the algorithm for implementations.
imposes a regularity requirement on the solution gW, U) in (2). In Appx. B, we show that this condition easily holds under linear models.
This result suggests that we can reject Hig when the discrep- ancy between p(ylu, x) and p(ylu) is sufficiently large to make p(w|x)) fails to account for all the variability encoded in p(y(ar). It has been similarly employed for testing Hy (Miao et al., 2023; Liu et al., 2023). However, (Miao et al., 2023) additionally required the equivalence condition for identification, which may not hold beyond factor models.
4.1. Solution existence under the null hypothesis
We show that H, can be tested by examining the integration equation (1), which can hold under Ho as the absence of di- rect effect from X to Y allows p(w|2r) to fully explain away the variability of (uylz) To formally claim this statement, we require the completeness condition.
In particular, (Liu et al., 2023) proposed a discrete approxi- mation of (1) and tested the linearity between P(y) X) and P(WIX), However, this method may suffer from informa- tion loss due to the discretization. In the subsequent section, we propose a novel estimation method to solve the integral equation for testing.
Assumption 4.1 (Completeness of P(UIW).. For any
3

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
4.2. Testing statistics via integral solving
Assumption 4.8 (Smoothness). For H(W, 1) in (4), we assume +) E Hw for all t. This implies that there exists a solution within the RKHS that satisfies (4).
In this section, we construct the test statisties to examine whether the solution exists in (1). To this end, we propose to solve the equation and measure its residue. When the equation (1) holds, previous studies proposed Maxinum Moment Restriction (MMR) (Mastouri et al., 2021: Kallus et al., 2021) for integral solving. In our scenario, it involves solving h(W) from the following moment restriction:
To solve H (w, 1) from (4), we employ recently developed nonparametric methods designed to estimate such condi- tional restrictions (Zhang et al., 2020; Mastouri et al., 2021; Kallus et al., 2021; Ghassami et al., 2022). Unlike these methods, we do not require the completeness of W to ensure the uniqueness of the solution. This distinction arises because most of these methods focus causal inference, where the goal is to accurately identify the bridge function to compute the causal effect. In contrast, our objective is to determine whether a solution exists. Therefore, it suffices for our estimate to approximate any valid solution to achieve this goal. In this paper, we estimate the least norm solution:
By,w {Y- [( = O.
(3)
However, it only leverages the first-order moment in- formation, it will lose the testing power, as illustrated by example B.3 in Appx. B, where (3) holds under H4. To address this, we propose novel estimation method, which leverages the information of all order moments.
Proxy Maximum Characteristic Restriction. To test whether p(y|zr) equals to h(w,, y)p(wlz)dw == -(yl>),
H HW = arg min (H,,W())
H(W.))EHW.o
we consider the following equation:
where Hw.o C Hw contains all solutions in (4). We leave more details in Appx. C.4.
(4) where we choose (YY,t) t) to contain more information than the first-order moment (Stinchcombe & White, 1998), and set H (W, t) as I Tyy, t)h(w, ()p(ylz)dy to make (4) holds. A common choice for t) is exp(ity). where T can be an arbitrarily chosen neighborhood around O. In this case, By {yp(Y,t)) is the characteristic function, and we hence call (4) the Proxy Maximum Characteristic Restric- tion. Since the characteristic function can uniquely deter- mine the probability density and therefore all order mo- ments, solving (4) offers greater utility to identify causal relations. In practice, we can set tY, = sin(ty) and cos(ly), and test whether (4) holds for these choices.
Ey.w tY, (W, x)(x) = OVEET,
Equivalently speaking, given RKHS Hw and Hx C L (Fx() with kernels kw and kx, respectively, the goal is to solve HW, from (- = t(,,t), where A Hw + Hx is a compact operator such that AHW, ())) == E (W, x))xx()); and b(., = EY t)ox (X)}. To solve H, we first note that (4) means, for any g E Hx, we have Wy.w,x [[p(Y, H (W,t))g(X) = for al- most all t. Therefore, similar to (Mastouri et al., 2021), we take g over unit- of Hx, and minimizes:
R(H)= = sup
(E ()) - H(W,
S1
Corollary 4.7 shows that H (w, t) belongs to [2(F(w) for all t, which ensures that (4) is solvable.
Let A(W, == ((Y,)) - HWW,)). (Mastouri et al.,
2021) provided an equivalent form of the risk:
Corollary 4.7. Suppose conditions in Theorem 4.5 hold. For any I, H (w, 1) in (4) exists and belongs to L (F))..)).
R(H)
Y (Xxx
where X',Y' W are independent copies of X, Y, W. Fur- ther, (Zhang et al., 2020) demonstrated that under some conditions for 4x. minimizing R(H ensures us to find the solution. To implement, we propose to minimize the regularized empirical risk:
As demonstrated by (Horowitz, 2012), achieving uniform consistency in testing the existence of a solution to the condi- tional equation is impossible. Therefore, certain smoothness conditions must be imposed to ensure the feasibility of solv- ing the equation. Following existing studies (Mastouri et al., 2021; Ghassami et al., 2022), we solve the solution in the reproducing kernel Hilbert space (RKHS) denoted by Hw.
RA(H) = Σ
A. n²
Kx + H| Hw: H
(5)
Formally, let kw be the reproducing kernel for the RKHS Hw. By the spectral theorem, we can rewrite (w, w') in terms of the eigenvalues and continuous eigenfunctions as kw (w, w') = ΣS ()pp(())))) (w'), where {yyly is the orthonormal basis of (F(F())). We can therefore characterize Hw as Hw
where A ((r, t) H (w,, t) and Kx.y := kx (x, zz). Using the representer theorem (Schölkopf et al., 2001). the estimated function A 1) for a fixed t can be written as A t) = a (w),, where kw(w) = (kw(w;, (w;, E R" and the coefficient a is given by a = (Kw KKWK + n²>Kx) "x KKKKwply.t t). Here, Kx = (1x(2,, E Ruxn and Kw =
H & (F (H.40) < "4
4

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
(Kw (w,, wy) E Ruxn are Gram matrices, and ply.t) 1) = (p(1),)),.. sen via cross-validation.
such that JJ [G((s, () < oo, where G(s, 1) is Gaussian process with zero-mean and covariance:
T)) The optimal À is cho-
t)(s,t),( (s', = EEne.(0mmm (O')),
Constructing the testing statisties. Ideally, if H can well approximate the solution, the equation (4) approximately holds for A/* Therefore, we can assess the validity of Hlo by evaluating the residue of the equation. To this end, we employ the conditional moment test procedure (Bierens, 1982; Bierens& & Ploberger, 1997).
where O' == Y is an independent copy of O; and (ii). p.mm converges weakly to max J [G(,, du(s)
EET
Remark 4.13. For brevity, we only present asymptotic re- sults for Tn (s. t) being real-valued function, or as the real and imaginary parts of a complex- -aluedd function, although they can be trivially extended to complex-valued functions.
To enhance the power, we choose a weight function m(. 8) that transforms the conditional restriction to unconditional one. Commonly chosen weight function includes char- acteristic function, exponential function, sine and cosine functions. According to (Stinchcombe & White, 1998), these functions ensure that for any UWW, Y, :- ((, (W,, t) with (UU Y, x * 0, the set of s E T such that (WW Y, ())} = 0 has Lebesgue mea- sure zero. Let Ü(W, X, = (X, - (W, (). we define
Power analysis. We consider the power performance under two alternatives, under which (4) has no solution. First, we consider the global alternative that has been similarly con- sidered in proximal causal discovery (Liu et al., 2023). That is, for any H (ww, t) E Hw for all t, the global alternative Hii satisfies the following:
HIfx EY, 11 (W))) + 0 for some t E T.
1 a
(s, = Vn Σ U. Vr. f))m(s, s). S.EET. (6)
Besides, consider a sequence of local alternatives Him, There exists H () E Hw for all I, such that:
Since testing (4) for all E T is equivalent to evaluating the maximum of the residuals over T, we define the final statistics for testing Hp as:
HIG EE = ((H + r(X,t) t)
1A
12
where 0 < S and r r(X,1) E Hx. To be a valid alternative, r(X, t)n" can not be written as B{H- for any H E Hw. Theorem 4.14 suggests that our statistics has asymptotic power of one under HIfix and Hi, when a < 2. and has nontrivial power when a = 1/2
/ [T (T, TET
A
p.
max
(7)
where y denotes the measure of 'T (e.g., Lebesgue measure).
4.3. Asymptotic behavior
Theorem 4.14, Suppose assumptions in Theorem 4.12 hold. Besides, we assume E(r(X, < co. Then, we have:
We study asymptotic properties of Apm. We first introduce some regularity conditions.
(i) Global alternative. lim,, maxIET (Ta(s, () = 00
for almost all under Hfx
Assumption 4.9. We assume Ex(m(X,s)|M and Ex (xx,m s) I) are uniformly bounded for all s.
(ii) Local alternative (a < 1/2). lim, MAXIET
(x ()] = 0o for almost all under HIg.
Assumption 4.10. nd -> D0, nd² - O.
(iii) Local alternative (a = 1/2). (s, converges
Assumption 4.11. For any s,t T, E(U(W,Y, <
weakly to G(s, t) + u(X, s, f) such that IJ (G(s, t) + s.X.stt) (ddpusssddiii
oo and ((Im(X, s)= {A(A*
< oc, where
< Do under Hi,, where G(s,t) is defined in Theorem 4.12 and s(X, = E [r(r (A* A ((xx)))m((X)))
(.(+) = (x. s) s) (W)](.)
Assumptions 4.9-4.10 are standard in kernel estimation methods (Darolles et al., 2011; Babii & Florens, 2020; Bey- hum et al., 2024). Asm. 4.9 imposes regularity conditions on the weight function m, while Asm. 4.10 ensures that the regularization bias vanishes asymptotically. Additionally, Asm. 4.11 is required to control the asymptotic variance of the test statistic, which has been similarly assumed in (Huang et al., 2022).
4.4. Implementations
To implement the testing, we need to compute the statistics A p.m in (7) and determine the critical value.
Computing Ap.m. Since Apm involves the integration, generally we should employ Monte-Carle methods for ap- proximation. For computational convenience, we can set m(.,s) s) to the characteristic function and to be symmetric around the origin (e.g., Lebesgue measure). Such choices en- able the integration to be computed in closed form. Besides, according to (Stincheombe & White, 1998), this choice of
Theorem 4.12. Let (O) = U(W, Y, tt)m(X s) -U(W, Y, (A(AAAAAA)) A)- 9o (X), with O == (W, Y.). Suppose assumptions 4.11..11. C.2-C.4, and D. D.2 hold. Under Hi, we have (i), (a 1) converges weakly to G(s,t)
5

---

Bivariate Causal Discovery with Proxy Variables:
Integral Solving and Beyond
m can preserve power when transforming the conditional restriction to the unconditional one. To compute the max- imal value of [ [Tm[ss z)) du((s) over T. we evaluate the process at a grid of equi-distan indices f E [K]] and [ S (s, dd)(ss Corol- lary 4. 4.16 shows that when K is sufficiently large, Apm converges to max ST [G(s, ()(s, dd(ss
equation (I) will also have a solution, rendering the causal relationship non-identifiable.
Proposition 5.1. Suppose U, X, Y, W follow from the linear
Gaussian model, i.e. U = EU,, X a Ogg + Ex: W BUU + 90 +ew,Y =UU + + X + 7x + 7o + sy, where EU, EY, ~ NN 1). When "/ = 0. as long as hxx > 9x Ex(o,,BB (,, the integration equation (1) hus no solution. Further, if [//l > gw Bu, Tv)', (1) has a solution.
estimate Å
Remark 4.15. If Ü (W,Y, .Y)) is a complex function, we can respectively compute Re( (e))). and Im Im(A (A.p).). Our test statistic is given by S = max ( Re( p...) Im( (mp))))))
Remark 5.2. In the proof of Prop. 5.1, we also discuss when the solution exists under IHl.. We find that the confound- ing strength between W and U, specifically Sr. must be sufficiently large to make the solution exist.
Critical value. Since it is difficult to obtain the explicit form of G(s. t), we employ the residue-based wild boot- strap procedure for approximation under the null-limiting distribution. We repeat the procedure for B times. For the b-th time, we first employ the empirical process Th (8, 1) = /n En=1 (wr, ys. ())((..,)) s) to approximate T, (s. t)
To illustrate, we consider example 5.3. We find that when " is sufficiently large, the solution exists even when X + Y, leading to significant drop of power shown in Fig. 2.
1.
for each (s, t), where [w]]]= is a sequence of zero-mean, unit variance variables. Here, we follow (Mammen, 1993)
Example 5.3, Suppose that X, U, W satisfy the linear Gaussian model, L.e. U = EU, X 20 + W 2U + ew, Y = X + U + W + &Y, where Ev, &w, ex ~ (0,, and W, X are standardized. The integral equation (1) has a solution if and only if yw
to set P(u; = 1 - k) = </ s/v5 and PP(w; x) = 1 k/V5 with k= = 15+11 The bootstrapped statistic is given by:
2
-54336/5 0.61. 72+16v5
[T, (x, 42) du))ss
p.
max kE[K]
2.0
Teat Power Rate
Ova samelt
Given the level of significance a, the critical value is com- puted as the
0.8
(1 a)--qantile of { ( AB p. }. de- We then reject the null hypothesis if
2.8
noted by A
1-a p.m.
0.4
Corollary 4.16 shows that the bootstrap p.m. converges maxIeT Sr [G(s, t) dd(ss) Corollary 4.16. Suppose assumptions in Theorem 4.12 hold. If (f((,)) 1) is continuous with respect to L f.for each y, then Ap.m is weakly convergent maxteT [G(s,, di((s) under Hig n, K- + co. Besides, conditional on the origi- nal sample (vi, Wi, Wi, Yi the bootstrapped statistics (8) is also weakly convergent to the maxIeT J [G(ss ()), dd(ss)
B p.m
3pm
2 A 1-o
9.m-
statistics A
54 02
0.0
0.5
0.9
1.4
Figure 2. The change of power across TW in example 5.3.
The result suggests that, when the effect from W + Y is strong enough, using only the NCO (i.e., W) may fail to identify the causal relationship. In the next section, we show that when the additional proxy Z (i.e., NCE) is available, the identification becomes possible again.
Remark 4.17. Many choices of sov, t) satisfy the continuity condition, including characteristic function, sine and cosine functions, etc.
6. Hypothesis testing with two proxies
5. Nonidentiiiaailiit with only NCO
To discover the causal relation when (1) has a solution un- der the alternative hypothesis, we assume Z (i.e., NCE) is available, which has been commonly employed in proximal causal inference (Miao et al., 2018: Cui et al., 2024). When Z is available, we analyze the properties of the solution under the null hypothesis to identify the causal relation.
To identify the causal relation, our method and previous studies in proximal causal discovery (Miao et al., 2023; Liu et al., 2023) rely on the assumption that the integral equation (I) has no solution under Hl,. However, it remains unelear when this condition holds. To clarify this, we explore the condition under the linear model.
Specifically, when h(w, y) satisfies (1), by W 1 X we have p(ylır) = [ [ h( y)p(wlu)dd } p(u(r)du for all (z. v). We will examine the property of the solution J h(w, (w)(w) to test Hlg. To this end, we require the
Proposition 5.1 suggests that the direct effect needs to be sufficiently large for the condition to hold. Additionally, we surprisingly find that when the effect from W to Y is strong enough, even if the effect from X to Y is also strong,
We leave the detailed form of gx, gw in Appx. E.3.
6

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
following completeness condition.
effect. In particular, Appx. E.2 provides a counterexample to illustrate that the causal relation may still exist even when there is no causal effect. A more comprehensive discussion is also provided in Appx. E.2.
Assumption 6.1 (Completeness of ()U For any square-integrable function g. we assume [g)[(]] = 0 almost surely if and only if (()) = 0 almost surely.
Similar to Sec. C, if A can well approximate the solution of (4), the equation
This assumption was also made in (Liu et al., 2023; Miao et al., 2023). It ensures that for any 1( y) that satisfies (1), we must have
Ey,w (ew((((Y)))))))) - (W, [1,,)) OVIET, (1I)
p(ulu)= = / Now, v)p(wo/ udw under THl.-
also approximately holds for all t E T. Therefore, we can assess the validity of Hg based on the residue Ü. Based on U, we can construct the statistics Tx:' (s. f), ppm, with the weight function m(Z, x 8) over Z.
(9)
Testing the null hypothesis through Z. To examine (9) that involves the unobserved confounder U, we employ Z to introduce another restriction. We require the following completeness condition that standard in proximal causal inference (Miao et al., 2018; Tchetgen et al., 2024).
Asymptotic behavior. Theorem 6.6 establishes the weak convergence of T () and A p.m under lil..
(2)
Theorem 6.6. Denote (0,z) = U(WYY, x.].
Assumption 6.2 (Completeness of P(U Z, x)). For any 9 (F))) and any 2,, we assume ((g(u)) 2; = 0 almost surely if and only if g(u) = 0 almost surely.
[m(Z,, 2. ( Emm(Z, x, 22,),,$6w... ())((W))() and O == (W, ).). Sup- pose assumptions in Theorem 4.12 hold. If Asm. 6.1-6.2, and E.6-E.8 hold, under Ho we have (i). T T122(
9. (x)x))) where gs (,()))) ==
2 (s, con- verges weakly to G(s, t) s.t. If [G(s, du(sss(ss))))))ss)s))))s ())p(()) < co, where G(s,I) is a mean-zero Gaussian process with co- variance [) t), (s',t') = (O, (O',z))
Taking expectation over p( p(u|z 2) on both sides of (9), and obtain the following for any x::
p(ulz. x) = f hww 2))(w
x)
(10)
where O' == (W", "Y,Y", X') an independent copy of O; (2) (ii). A p.m converges weakly to max J (ss)))))((())) LET
for all 2). Under assumption 6.2, we can check Ho by first solving h from (1) and examine whether it satisfies p()],,z) f h(w, v)p(w)/,,r)dd We summarize it into the following theorem.
Similarly, we can establish the asymptotic power for A where the global and alternatives are defined accordingly in terms of tY, HW, 0)) x}.. Due to space limit, we leave the result and its proof in Appx. E.4.
Theorem 6.3. Suppose assumptions 6.1, 6.2 hold. For any h(w, y) that satisfies (1), Hy holds if and only if h(w, y) also satisfies the integral equation (10) for any fixed x.
7. Experiments
Remark 6.4. Solving hww,, v) from (10) is different from that of p(y]](())))) x, = J- h(w, y, x x))((w) (w jz, x adw in (Miao et al., 2018). The goal of (Miao et al., 2018) is to solve hww,, y, x) for identifying ply do(x))
In this section, we evaluate our methods on synthetic data. We consider two settings: (i) (Sec. 7.1) W A Y where only
the NCO W is available; (ii) (Sec. 7.2) W + Y where the additional proxy Z (i.e., NCE) is provided?.
J h(w, y, x))((w) () which thus allows/ h to depend on x. In contrast, our goal is testing whether X directly affects Y, so the solution h should be independent of X while ensuring that (10) hold as æ varies.
Compared baselines. For the single-proxy setting, we compare our methods with: (i) Liu (Liu et al., 2023) that de- signed a discretization method for bivariate causal discovery over continuous variables; (ii) KCI the Kernel-based Con- ditional Independence test (KCI) (Zhang et al., 2012) that tested the null hypothesis of X 1 Y|W using kernel matri- ces. For the two-proxy setting, we also conduct (iii) Miao (Miao et al., 2018) that was designed for causal hypothesis testing over discrete variables using W and Z.
When X, Z, and W are discrete, (Miao et al., 2018) pro- posed testing the linear relationship between (P(y) hi and (P(W) Z2 }; across the values of Z, X. Our proce- dure can be seen as the continuous counterpart to this ap- proach, with the difference being that after solving hww, v) from (I) for all x, we only need to examine (10) for all z with a single a. rather than solving h(w, y) from (10) for all pairs of (I, 2).
Implementation details. We set the significance level a to 0.05. We choose o and m to be complex exponential functions. For PMCR estimation, we set K = 100 and follow (Mastouri et al., 2021) to select the optimal A from a
Remark 6.5. One might argue that when W, Z are available, the average causal effect is identifiable. In this regard, our analysis seems redundant as the effect of causal relation can be quantified. However, we would like to mention that the causal discovery conceptually differs from the causal
"ode is available at https: / //tpthb.ttttt////////////t.bbt com/yezichu/ proxi ma causa 1 _dd sdiscovery_cv
7

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
sequence ranging from 4.9 10-6 to 0.25, with step size chosen ensure the sequence contains 50 values. Besides, we use Gaussian kernels with the bandwidth parameters being initialized using the median distance heuristic. For the procedure of Liu (Liu et al., 2023), we follow the paper to set the bin numbers of W and X to Ix = 14, lw = 12, respectively. For the procedure described in Miao, we implement the R code released in the paper and set 1x 3, lw 2,1z 2 by default. For KCI, we adopt the implementations provided in the causallearn packages https: http:///caaaaaaaaa learn. teadthedoccs.on
tiveness of our estimation method PMCR over the MMR, we apply PMCR to the data generated in Example B.3, where we have shown that the solution of the first-moment equa- tion exists under the alternative hypothesis. As shown in Fig. 4, although both methods can asymptotically control the type-I error as n + oo, the power of our procedure approaches while the first-moment method (i.e., MMR) still approximates a 0.05 under Hi..
7.1. Single proxy with W AY
We first evaluate our method in Sec. 4 to the setting where only Wi is provided.
Data generation. We follow (Liu et al., 2023) to generate V E {X,Y,U W via V = fPAA + EV. where PAy and EV respectively denotes the par- ent set and the noise of V. For each V. fv is randomly selected from {linear, tanh, sin, sqrt]. Be- sides, the distribution of εv is randomly chosen from (Gaussian, uniform, exponential, gamma). To mitigate the effect of randomness, we repeat the process 20 times. At each time, we generate 100 replications under each Hlo and Hl, and record the type-I error rate and power rate.
[2]
[1)
Figure 4. Type-I error rate (left) and power rate (right) of our pro- cedure with PMCR and the first-moment method in example B.3.
7.2. Two proxies with W Y
In this section, we apply our method in Sec. to the setting W + Y with both W and Z are available.
Data generation. Following example 5.3*, we set "w 1. which implies there exists h that satisfies the integral equation (1). Similar to the single-proxs setting, we repeat the process 20 times, where at each time we generate 100 replications under H, and [H],
Dut STn
Type: Uirrr Baze
Tett Poseer hate
(1)
Figure 3. T error rate (left) and power rate (right) of our test- ing procedure and baseline methods in the single proxy setting. The solid line reports the average value over 20 times, and the shaded area denotes the region (mean- std, mean std).
4.4
1000
25
400
8/4
180g
Type-I error and power. In Fig. 3, we report the average Iype-I error error rate and power rate for testing procedure and others. As shown, the type-I error rate of our method closely approximates a 0.05 as increases, while other methods fail to control the type-I error. Specifically, conditioning on the proxy W. KCI cannot eliminate the confounding bias, leading to uncontrollable type-I errors; while the additional error in Liu (Liu et al., 2023) may arise from discretization errors with finite bin number or probability estimation error due to limited sample size. Besides, our power approximates to one increases. Compared to previous baselines Liu, these results demonstrate the utility and its ability to make better use of available data in causal discovery.
Figure 5. Type-I error rate (left) and power rate (right) of our pro- cedure and baselines on synthetic data with two proxies.
Type-I error and power We report the average results in Fig. 5. As shown, although our single-proxy procedure can control the type-I it suffers from low power in learning the causal relation, due to the existence of solution under Hi in this example. With additional information provided by Z, the power significantly improves and approaches one as n increases. This verifies our findings in Sec. 5, and demonstrates the utility of employing Z (i.e., NCE) in discovering the causal relation when the effect of W on Y is strong enough to invalidate the procedure with only W.
Comparisons with MMR, To further demonstrate the effec-
We also consider a nonlinear setting, as detailed Appx. F.
8

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
8. Conclusions and discussions
Bierens, H. J. and Ploberger, W. Asymptotic theory of inte- grated conditional moment tests. Econometrica: Journal of the Econometric Society, pp. 1129-1151, 1997.
We introduce non-pann con-aaametric procedure for causal hy- pothesis testing and establish its asymptotic properties. Ad- ditionally, we show that causal relationships may not be identifiable through examining the integral equation with only NCO. We then leverage the additional NCE that can effectively recover the causation when it exists. We believe our findings, supported by theoretical justifications, provide new insights into proximal causal discovery.
Calonico, S., Cattaneo, M. D., and Farrell, M. H. On the effect of bias estimation on coverage accuracy in non- parametric inference. Journal of the American Statistical Association, 1133522):767-779, 2018.
Carrasco, M., Florens, J.-P., and Renault, Linear inverse problems in structural econometrics estimation based on spectral decomposition and regularization. Handbook of econometrics, 6:5633-57555 2007.
Currently, we only analyze the integral-related assumption when NCO has a direct effect on the outcome. In the future, we will extend our analysis to the case when there is a bidirectional edge between them.
Colangelo, K. and Lee, Y.-Y. Double debiased machine learning nonparametric inference with continuous treat- ments, arXiv preprint arXiv:2004.03036 2020.
Impact Statement
This paper presents work whose goal is advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.
Cui, Y., Pu, H., Shi, X., Miao, W., and Tchetgen Tchetgen, Semiparametric proximal causal inference. Journal of the American Statistical Association, 1195546):34881333 2024.
Acknowledgements
Darolles, S., Fan, Y., Florens, J.-P., and Renault, E. Non- parametric instrumental regression. Econometrica, 79(5): 1541-1565, 2011.
This work was supported by Young Scientists Fund of the National Natural Science Foundation of China (Grant No. KRH2305058); the State Key Program of National Natural Science Foundation of China under Grant No. 12331009 STI 2030-Major Major Projects(No. 2021ZD0200407); National Key Research and Development Program of China(No. 2022YFC2405100, The computations in this research were performed using the CFFF platform of Fudan University.
Davis, P. general expression for hermite expansions with applications. The Mathematics Enthusiast, 21(1): 71-87, 2024.
D' H'Haltfocuille, X. On the completeness condition in non- parametric instrumental problems. Econometric Theory, 27(3):660-71,, 2011.
References
Florens, J.-P., Johannes, J., and Van Bellegem, S. Instru- mental regression in partially linear models. The Econo- metrics Journal, 15(2):304-322 2012.
Andrews, D. W. Examples of 12-complete and boundedly- complete distributions. Journal of econometrics, 199(2): 213-220, 2017.
Fukumizu, K., Bach, F. R., and Jordan, M. Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces. Journal of Machine Learning Research, 5 (Jan):73 99.779999919 2004.
Babii, A. and Florens, J.-P. Is completeness necessary? estimation in nonidentified linear models. arXiv preprint arXiv:1709.03477 2017.
Babii, A. and Florens, J.-P. Are unobservables separable? Econometric Theory, pp. 1-33, 2020.
Ghassami, A., Ying, A., Shpitser, I., and Tchetgen, E. T. Minimax kernel machine learning for a class of doubly robust functionals with application to proximal causal inference. In International Conference on Artificial Intel- ligence and Statistics, pp. 7210-7239. PMLR, 2022.
Babusci, D., Dattoli, G., and Quattromini, M. On integrals involving hermite polynomials. Applied Mathematics Letters. 25(8):1157-11666 2012.
Beyhum, J., Florens, J.-P., Lapenta, E., and Keilegom, I. V. Testing for homogeneous treatment effects in linear and nonparametric instrumental variable models. Economet- ric Reviews, pp. 1, 2024.
Guo, R., Cheng, L., Li, J., Hahn, P. R., and Liu, H. A survey of learning causality with data: Problems and methods. ACM Computing Surveys (CSUR), 53(4):1-37, 2020.
Horowitz, J. L. Specification testing in nonparametric in- strumental variable estimation. Journal of Econometrics, 167(2):333-36 2012.
Bierens, H. Consistent model specification tests. Journal of Econometrics, 10(1):105-133 1982.

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Hu, Y. and Shiu, J.-L. Nonparametric identification using instrumental variables: sufficient conditions for complete- ness. Econometric Theory, 34(3):659-69) 2018.
Miao, w., Hu, W., Ogburn, E. L.,, and Zhou, X.-H. Iden- tifying effects of multiple treatments in the presence of unmeasured confounding. Journal of the American Sta- ristical Association, 1 186543):1933-1977 2023.
Huang, W., Linton, O., and Zhang, Z. A unified framework for specification tests of continuous treatment effect mod- els. Journal of Business & Economic Statistics, 40(4): 1817-1830, 2022.
Newey, W. K. and Powell, J. L. Instrumental variable esti- mation of nonparametric models. Econometrica, 71(5): 1565-1578, 2003.
Olver, W. NIST handbook of mathematical functions hardback and CD-ROM. Cambridge university press. 2010.
Kallus, N., Mao, X., and Uehara, M. Causal inference under unmeasured confounding with negative controls: A mini- max learning approach. arXiv preprint arXiv:2103.1 14029, 2021.
Robins, J. M. Confidence intervals for causal parameters. Statistics medicine, 7(7):773-785, 1988.
Kato, M., Imaizumi, M., McAlinn, K., Kakehi, H., and Yasui, S. Learning causal models from conditional mo- ment restrictions by importance weighting. arXiv preprint arXiv:2108.01 0130000..11.22222 2021.
Robins, J. M., Scheines, R., Spirtes, P., and Wasserman, L. Uniform consistency in causal inference. Biometrika, 90 (3):491-515, 2003.
Kress, R. Linear integral equations, 1989.
Schölkopf. B., Herbrich, R., and Smola, A. J. A general- ized representer theorem. In International conference on computational learning theory, pp. 416-426. Springer, 2001.
Kuroki, M. and Pearl, J. Measurement bias and effect restoration in causal inference. Biometrika, 101(2):423- 437, 2014.
Seeger, M. Gaussian processes for machine learning. Inter- national journal of neural systems, 14((2))69-1066 2004.
Lapenta, E. and Lavergne, P. Encompassing tests for non- parametric regressions. Econometric Theory, pp. 1-30, 2022.
Spirtes, P., Glymour, C.,, and Scheines, R. Causation, pre- diction, and search. MIT press, 2001.
Li, Q., Hsiao, C., and Zinn, J. Consistent specification tests for senipraametic/nonpaanntii models based on series estimation methods. Journal of Econometrics, 112(2): 295-325, 2003.
Stinchcombe, M. B. and White, H. Consistent specification testing with nuisance parameters present only under the alternative. Econometric theory, 14(3):295-325, 1998.
Tchetgen, E. J. T., Ying, A., Cui, Y., Shi, X., and Miao, W. An Introduction to Proximal Causal Inference. Statistical Science, 39(3):375 - 390, 2024. doi: 10. 10.1214/23-SS9111 URL https: htt///do..orrr org /10 0. 121 221-/////111111...11111111111111
Lipsitch, M., Tchetgen, E. T., and Cohen, Negative controls: tool for detecting confounding and bias in ob- servational studies. Epidemiology, 21(3):383-388 2010.
Liu, M., Sun, X., Qiao, Y.. and Wang. Causal discovery via conditional independence testing with proxy variables. arXiv preprint arXiv:2305.0528/ 2023.
vd Vaart, A. Asymptotic statistics. cambridge series in statistical and probabilistic mathematics, 1998.
Wellner, J. al. Weak convergence and empirical pro- cesses: with applications to statistics. Springer Science & Business Media, 2013.
Mammen, E. Bootstrap and Wild Bootstrap for High Dimen- sional Linear Models. The Ammals of Staristics, 21(1):255
- 285, 1993. doi: 10.449995 URL https: ///o//. otloo//g////// 2214/2/9/1176349922
Zhang, K., Peters, J., Janzing, D., and Schölkopf, B. Kernel- based conditional independence and application in causal discovery. arXiv preprint arXiv: 2012.
Mastouri, A., Zhu, Y., Gultchin, L.. Korba, A.. Silva, R., Kusner, M., Gretton, A., and Muandet, K. Proximal causal learning with kemell:: Two-stage estimation and moment restriction. International Conference on Ma- chine Learning, pp. 7512-7523. PMLR, 2021.
Zhang, R., Imaizumi, M., Schölkopf, B., and Muandet, K. Maximum moment restriction for instrumental variable regression, arXiv preprint arXiv:2010.07684, 2020.
Miao, W., Geng, Z., and Tchetgen Tchetgen, E. J. Identify- ing causal effects with proxy variables of an unmeasured confounder. Biometrika, 105(4):987-993 2018.
10

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Appendix
Notations
12
Existence of solutions with single proxy
14
B.1 Proof of Theorem 4.5
14
B.2 Proof of Corollary 4.7
16
B.3 Counter-example to the solvability of the first-order moment equation under H1
17
B.4 Verification of [J
g(w,u) ((W)p(()ddddr < co in Theorem 4.5 p(a)
17
Proxy Maximum Characteristic Restriction
20
C.I Bounded linear operator
20
C.2 Hilbert space
20
C,3 Reproducing Kernel Hilbert Space
21
C.4 Validity of optimizing (5)
C.5 Tikhonov regularization
21
C.6 Ill-posed inverse problem and solutions
22
D Proofs of Asymptotic Properties
24
D.1 Proof roadmap and key assumptions
24
D.2 Empirical process
25
D.3 Expected risk difference
26
D.4 Proofs in section 4.3
32
D.5 Auxiliary lemmas
42
E Existence of solutions with two proxies
47
E.1 Proof of Theorem 6.3
47
E.2 Discussions of causal inference and causal discovery
47
E.3 Proof of Proposition 5.1 and example 5.3
49
E,4 Proof of asymptotic properties with proxies
55
F Additional experiments
58

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
A. Notations
introduce notations used throughout the appendix.
Table I: Notations.
Notation
Definition
Z, W
Negative control exposure, negative control outcome and unobserved confounder
P(X)
for any discrete variables with it categories
p(v1)=1)
p(v1 (mk)
for any discrete variables Y, X with 1, categories
p(unkar)
p(u))a))
((y(=1,2))
plul 01, zm)
P(Y y|x, 2)
for any discrete variables X, Z with k m categories
plyl MK, Ik, 21)
plylæx zm)
Hw, Hxx
The reproducing kernel Hilbert spaces (RKHS) defined on the domains of W and X
ow (w), Øx (2)
The canonical feature map defined on the domains of W and
(wo, w'), kx x2)
The reproducing kernel for the RKHS Hw and Hx
R(H)
The population loss function defined in Eq. (22)
(A)
The regularized empirical risk defined in Eq. (5)
(x),(x) = t(x,t)
The operator and right term defined in Eq. (24)
(,b,(z) = t(z,t)
The plugging operator and right term defined in Eq. (26)
A"
The adjoint operator of A and Å defined in Eq. (27) and (28)
(Ag.
The singular value decomposition of the operator A
Hw.o
The set of all solutions of Eq. (22) defined in Eq. (29)
HA = (w, 1)
The population Tikhonov regularization solution defined in Eq. (31)
(A(w) (w) = (W) t)
The empirical Tikhonov regularization solution defined Eq. (32)
(w) (w,)
The least norm solution in Eq. (4)
Ker(A)
Ker(A) (H Ah 0} is the null space of the operator A
Ran(A)
Ran(A) {f Ah = f} is the ranged space of the operator A
The space of square-integrable functions with respect to the cumulative distribution function F(w) and F(x)
(FF(w)),
We say G(s,t) E (2 T, x n) fJ [G ())d(())) < x
m(,)),
The weight function
9. [mm(X, (W))
9.
(W,YY,),&(WW (s,
The residual ((Y,t) H (W, t) and estimated version (Y, (W,
The statistics defined in Eq. (6)
A pm
The statisties defined in Eq. (7)
12

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Notation
Definition
B())
The expectation with respect to both random variable and data
P()))
The expectation with respect to a random variable alone The empirical expectation with respect to a random variable given data
P.(())
|| ||=
The norm with respect to space F
13

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
B. Existence of solutions with single proxy
Let [Fx))) denote the space of all square-integrable functionn of x with respect to a cumulative distribution function F(x), which a Hilbert space with inner product (g), 92) = I 91 (z)ga (e)gg(()pp) (()))))) Let 'T' denote the operator: [2)F(w)) (3(F())), Tg- [ - and let (Ans PP Pn === denote singular value decomposition of T'.
Assumption B.1. Assume the following conditions for all uc:
(1) ( < do and p(pp < oc;
(2) 1 (p(u)x) on,) < oc.
B.1. Proof Theorem 4.5
We first prove that under the conditions Theorem 4.5, there exists a solution g(e, u) E [F([(W)) for all u, such that p(][z) = g(w, ()p(ulz)dd Our proof is based on Picard's Theorem, which is presented below.
Lemma B.2 (Theorem 15.18
18 in (Kress, 1989)), Given Hilbert spaces Hi and H2, a compact operator T H1 - H2 and its adjoint operator T Hg + H1, there exists a singular system (Ays Pnx Qn) t of K with nonzero singular values {)a) and orthogonal sequences (pn, H1), (Φn E H2). Then the equation of the first kind Th with f H2, has solution if and only if
1. E -er(T"*++ where Ker(T")) = {h: T = 0} is the null space of the adjoint operator T*; 2. Σt (3,1(, < +o0.
Proposition 4.2. Under assumption 4.1 and some regularity conditions in B.1, there exists g(w, u) E (F all u, such that it solves the following integral equation for all (u, z):
p(ulz) - =
g(w, u)p(w
(2)
Proof. Our goal is to show the solution exists for the following estimator:
T: 22 [F 23 (F(x)) Tf= =[][(((()) f (F)(
Besides, we consider the following operator:
S [F((x) + (2 Sg = -[o(()] = f (2
By Lemma B.2 and assumption B.1 (2) for p(u|z), the conclusion holds if we show that is the adjoint operator of T, T is compact, and that p(u(z) Ker(S)
(i). S is the adjoint operator of For the operator T, T. E L and Vg E (F(((=)), we have
(TF.0) = Ex -EX(L((W)))))(XX)
Ex LUIX LEX ((W)) U)x())
② Ex XL
Ee,x [LLx:(B(((()))))))xx) (V(BBL((W))0C)](((0))))
Similarly,
(f
= Ew Ew (FW
= Ew B(
= -B//BL(()111)(((()1))1)
14

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where (1) and (2) follow from W L XIU Therefore, have
(T(:0) -(,
(ii). T is compact.
We define the integral kernel
((w, x) p( x)
(12)
(
Then for the operator defined previously, we have the following form:
Tj = / ((w,)))(()(P((()) = (E((()))) S (F())
(13)
Sg /
(14)
By the definition of K in (12), we have:
/f [K(, (ww
/f ((
Under assumption B.1 (I), we can apply Theorem 2.34 in (Carrasco et al., 2007) to obtain that the operator T is Hillert-SShmidt operator, which implies that a compact operator, according to Theorem 2.32 in (Carrasco et al., 2007).
(iii). p(u|2r) E Ker(S)' for any u. For each g E (S),, by iterated expectations, have
E(g(X)WW) = Euw (X) W,U! UH
Eow E [g( 0,
where (1) follows from W 1. X By the completeness assumption, we have
B(g(X)lu} = a.s.
(15)
Then we have
(o.p(ulx)))x(()))) Ex(o(X)p(u())) (p(ulz) (F(((x) is used here.)
/ g(z)p(u,z)dr= / (()) =0 0,
which implies 2(u) E Ker(S)+. Combining the above three steps together, obtain the conclusion.
Theorem 4.5. Suppose assumptions 4.1. 4.4 and the regularity condition B.1 hold. g(w, u) in (2) satisfies < oo, then under lolo, there exists h(w,y) E C²(F(w)) for all y, such that it solves the integral equation in fo) all (x.y)-
Jf g(oo,a)) plu)
Proof. By proposition 4.2, we have g(w, u) satisfies the integral equation p(u|z) g(w, u)p(w])z)dw. Then under Ho. have:
p(y|=) / p(ylu)p(ulz)ddu
U g(w, u)plylu)du } pl p(w|z)dw.
Therefore, we write an expression h(w,y) v) = f (w, u)p(yl du that satisfies solution of the integral Eq. (1).
15

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
We first prove that h(w,,) y) is well-defined for any fixed y. To be specific,
V o(w, w))(ul)) <
lg(w, w)l w(ullddd
g(w,, u)
p(u)- p(uly) ply)du
p(u)
p(u)
(
/
g(w, u) p(u)
Cp(y)
plu)du,
where (1) follows from assumption 4.4. Since we have assumed IJ (g (plppp(u)lpPllppppapaapaaa < co, there holds I (g(w,, u)) do for a.e by Fubini's theorem. Applying the Cauchy-Schwand inequality, we have:
/
g(w, w) p(u) p(u)du S
{
slw; u) 2 p(u)du p(u)
}
}
1/2
(vienduy 1/2
<0.
(16)
Thus, h(w, y) is well-defined for any y-
Next, we show that hww, y) E c (F(w)) for any y. By Cauchy-Schward inequality, we have
I [b(w,)/ v) p( = / { g(w, wpp(olu)dy } p(w)dw
| g(w, w)
p(olu)p(u)du } p(w)dd p(u)
/
g(w, u)
plu)du } { (n(ulw)
} plw)dw
p(u)
g(w, u) p(u)
p(w)p(() (ulpp(()ddddd }
{
p(uly) p(u)
p(v) plu))p(u)dd))
g(w, u) p(u)
p(w)p(u)ddull } { C (Cele)p((u)d
Since [ʃ]
p(w)p(u)dwddu we complete the proof. p())
B.2. Proof of Corollary 4.7
Corollary 4.7, Suppose conditions in Theorem 4.5 hold. For any t, H (ww, f) in (4) exists and belongs to (F(F(w)).
Proof. (i). We first prove that Hww, f) is well-defined. To be specifie, we take op be the complex exponential function eity Then, by H H(w, = f ply, () ddy we have
1f ( dd
I g(w, u)plulu)du dy S lestyy. f lotw, u) plrll)dudy
lem0
/ lettop {
g((w u)
p(u)- p(ul)) () p(u)
plu)duy }
p(u)
(1)
letty- p(y) { glwu, u)
2 C.
plu)du dy.
plu)
where (1) follows from assumption 4.4. Since |esty| = 1, p(y) dy = 1, and by Eq. (16), it follows that H (w, t) is well-defined.
16

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
(ii). We proof (w, 1) (F()). By Cauchy-Schwarz inequality, we have:
/ (w, (P(wwt))30l t)(
J g(w,
} p(w)ddw
= - U g(u, u)
) 2 p(w)dw
e plyl ((ppp))((((l)) (u)ppp)((dddd
p(u)
- / {
g(w, u) p(u)
plu)du }
HI e uu (")p(yluddy p(u)du j p(w)du 12 p(ull)ddy plu)du
(17)
g(w, u) plu)
p(w ((w)p())dodd
The proof is completed
|. (((p((y))dd) < and Jf
g(m..)
p( 2c. p(u)
B.3. Counter-example to the solvability of the first-order moment equation under H₁
Example B.3. Suppose that X, YYU, W satisfy the linear Gaussian model, i.e. X = ex, ax X og + Ev. W
3rU = 3o + sw, Y TUU + X + + Ey, where EU, Ex, are Gaussian noises. Then the solution of X = [( is given by h(w)= = +xxxx w + To + "1000 - (Bo + (l) + puax
poax
Proof. We show that [[Y = Eb. W + bolX), where
(bo, bw)) {70 + wua (go perao) 7x + Tux
7x + xwax
poax
fact, this can achieved by solving (bo, bw.) the following integral equation:
X(Y = E(b W + bol
For the left-hand side, have
X( = 70 + 7x + No B(U
= 70 + 7x + 7u B(as + xxx
(70 a00 + (7x + vx
For the right-hand side, we have
B[o(W)|X) = E(ho + b WIX)
= by - by B(Bu + moU
by bw Bu + My E(ao + ax xxx(X))
(bn + 30 buHwaa) + (b )X,
bw 7x + +x+xxxx
ba - (7x + waxx = (bo + Bo + bwpwva)) (%0 700))
puax
7x + v
by 700 + 7000 (Bo + m000)
puax
B.4. Verification of [f g(w.u) p(u)
(²p(())(())dur < co in Theorem 4.5
Lemma B.4 (Miao et al., 2018)). F N + 8, o3) and Y ~ /V166 + + oz)). then one can verify integral equation p(y[z) = hww, y)p(wu x ddw has a unique solution h(w, y):
hw y)
1
y
Tuz //B1u0
(18)
Tow
17

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where ① is the probability density function (pdf) of the standard normal distribution, Twx 7 -///// 1 and o d3 (n/))/3/ (8))²
Example B.5. We consider the linear Gaussian generation mechanism:
U o7
X = oo + ta0 + 1)
W Bu+ U ((),
Y 7o + +U - 1).
Then h(we, y) given by the following has a unique solution to the integral equation (I);
h(w, v)
1
y
Bu
50
70
1
Besides, at this time, the assumption in Theorem 4.5 has
{ (lwwww) p(u) p(uc)p(u)dddd
02-02k k/2 /o
< 00,
(19)
1 - o2k ko2
where & = (rar))/(o"aa) 1).
Proof. Based on the data generation structure, we can obtain joint distribution
U
o² aa
oaan o²
0² azay o² +1
X
+1
aaun
N
W
So
o² Bu oay
aaavoe
Bu
70
d'au'
o²Bun
oand +1
We now get the conditional distributions p( ((wlz) and p(y|zz) according to the joint distribution
pw + Cov(W, X)
(x - M2), Var(W) (1
Cov2(W,
Var(X)
Var(X) Var(W)
~N o'ayBy o'amBB
(0'ay
o2aat +1 o2a +1 ap Bo, 22 og +1 o²ait +1
Y x== Myy Cov(Y, Var(X)
(x ux), Var(Y) (1
Cov²(V,
Var(X) Var(Y)
N P'anN o'au'u
(o'auno)
o² +1 ao + 7o, 0222 + 0²a0z +
o² +1
By Lemma B.4, the solution to (1) with h(w, y) given above is
kww v)
1
1/
3y
Bo
"Yo
Next, we prove (19). We will first obtain g(w, u). Specifically, we first get p(ular):
U = x ~ N + Cov(U,
Ux), Var(X) (x (x)) Var(U) (1 (1-
Var(X) Var(U)) }
CU,2(U,
oFay a2au
(oax:)
N
22a +1" x- 0²a + 00,0² 0² o² o²a + +15
18

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Similarly, by applying Lemma B.4 with the form of Y|X replaced by that of U U|X above, we have:
glw, u)
Dv
30
(aRaw)2 (1 σ²)
(o 1
(1 0²)
o
We can easily verify that the following integration is finite, by the following:
((w)p(u)dudd o² o2 + k²
/
((w,,)). u) p(u) }
ag
1- ko
where :== (oop) ) o'a +1
19

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
C. Proxy Maximum Characteristic Restriction
For the sake of completeness, we introduce some preliminary concepts that are necessary to understand the theoretical analysis of our PMCR method. First, in section C.I-C.3, we introduce some background knowledge of the linear operators and Reproducing Kernel Hilbert Spaces required in this article. Upon this, we provide details on the derivation of our empirical loss (5) in section C.4. Section C.5 rewrites the loss into the Tikhonov regularized form, which serves the foundation of our theoretical analysis for Theorem 4.12.
C.1. Bounded linear operator
For two normed linear spaces F and & over R. function A F R) is called linear operator it satisfies the following properties:
- 9 (where FF and 9 are both normed linear spaces
1. Homogeneity: A(af) a(Af), for any E f F;;
2. Additivity: A(f + g) Af + Ag. for any f.g F.
Operator Norm and Boundedness. The operator norm of a linear operator A F G is defined
|| A
|| sup fer If
A linear operator A is called bounded if there exists finite constant C' such that for all f E F, we have
| < G|I/||=-
In terms of the operator norm, this condition is equivalent to saying that ||
C.2. Hilbert space
We begin by introducing definitions and basic properties of an inner product space. Based on this, we introduce the Hilbert space.
A function (x R is said to be inner product on F if it satisfies the following three properties
1. (a1fi + a2f2.9)) (1(f1.9) gz(
2. (f-9)==(g.1 ().
3. (f,)) 2 0 and x = 0 if and only if =.
One can always define a norm induced by the inner product: || = 2 For this norm, the following Cauchy- Schwarz inequality holds, i.e., Kf. S -
A Hilbert space is complete inner product space. This means, Hilbert space is inner product space in which every Cauchy sequence (a sequence where the elements get arbitrarily close to each other) converges to element within the space. An orthonormal basis of Hilbert space H is a set of vectors {ei], such that ||e/|| = 1 for each i and (ei, ej)x = 0 for each # j. Besides, each f E H can be expanded in Fourier series:
if. eires
Hilbert adjoint operator. In the context of Hilbert spaces, we can define the adjoint operator. Let Hi and He be Hilbert spaces, and let A H,, + H2 be a linear operator. The adjoint operator A H2 - H, defined by the property that for all
(Af.gina (f. A A* u1
The operator enjoys a number of important properties:
1. If A is bounded, so is A", and ||Allop = || Ilop;
2. (A*)* =A; A;
3. If A is invertible, so A", and (1+)-1 (1-1)*
20

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
C.3. Reproducing Kernel Hilbert Space
For any space W, let kw : W W + be a positive semi-definite kernel. A kernel called characteristic if P Ew-rllww (W; -]] is injective (Fukumizu et al., 2004). We denote by ow its associated canonical feature map éw(w) kw (w,) :) for any W, and Hw its corresponding RKHS of real-valued functions on W. The space Hw is Hilbert space with inner product (v) and norm || || It satisfies two important properties:
1. (w E Hw for all w E W:
2. reproducing property: for all Hw and W, f(w) - (f, (w (w,+) Ho
Since the Reproducing Kernl Hilbert Space (RKHS) is a Hilbert space, it satisfies all properties in secion C.2. Besides, we can define the kernel mean embedding, which helps to take the expectation of a function. Suppose we wish to calculate [{J(W}} for any E Hw. By the reproducing property and linearity the inner product, we have
B{f(W)) = f(w)dP(w) /
/ (f. ow (w)) dP(w) = ( ow (w())d
= (f.
Hw
The object Ww := 6w (w))()(( called the mean embedding of the distribution P(w). This property of RKHS implies that, to calculate the expectation of a function, it suffices to take the inner product between the function and the mean embedding of the corresponding distribution. Following this property, we can also calculate the expectation ((f(W))(X)) for any E Hw:
(Lf(W))(X)) = f f(w))( / (w
(F. m(z))w (w) (w)dd(W,z)
Hw
(20)
Finally, we introduce properties for the norm ||. llww. A function f e Hw if and only if [lf = f) oc. Further, if kw (w, )) is bounded, we have (= & Il To this, note that by Cauchy- Mluchy-Schwarz inequality, for any f E Hw, we get:
(5(w) (kw (we, D S [1ww(w, 11 1
Therefore, have
( & llfllw.w
(21)
C.4. Validity of optimizing (5)
Since (4) implies BY{+(Y,t) - H(W.t))s(X)) = 0 holds for any measurable functions g : X - R, we follow (Zhang al., 2020; Mastouri et al., 2021) to take g over a unit-ball of RKHS Hx with fixed kernel k", and minimizes
R(H)= =
sup VEHx.1ollEs
EE (E((s(Y,t) ((x)))
(22)
(Mastouri et al., 2021) provides an equivalent form of this risk, which is the population version of our empirical loss (5). Lemma C.1 (Lemma 2 in (Mastouri et al., 2021)). Assume that E E({{(Y H(W, () (X, X')] < oc and denote by X' an independent copy of the random variable X. Then R(H) = Bl( H HWW
() X')].
(Zhang et al., 2020; Mastouri et al., 2021) demonstrated that if the kernel function kx derived from the conditional variable X the conditional moment equation (4) integrally strictly positive definite (ISPD defined in Asm. C.4), continuous, and bounded, then the conditional moment equation (4) shares the same solution with R(H). That means, optimizing R R(H) ensures to find the right solution.
C.5. Tikhonov regularization
In this section, we rewrite our loss (5) into the following Tikhonov regularized form, which serves as the foundation to prove Theorem 4.12.
RA [6(2, t) ÂH(w, ()) 1ee + (w, 1)1l
(23)
21

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
This be achieved by reformulating the PMCR into linear ill-posed inverse problem in the RKHS. Specifically, let xx (x: = xxx (2r,-) and ow (w:) kww (w,-) ) be the canonical feature maps. Then, by (0x(()) = 1x(x, (x,)
R(H) of Lemma C.1 can be rewritten in terms of mean square error:
R(H) - [I 6) (X)))))))))
(Xx [1(X,)) (1554--AHUH) AH((,,)(XX))
where
bl,t)) = / ply,
fly,))xxx f)xx (x)p(a, y)drdy, AH(W,)) = I H H(w, f)ox (x)p(a, w w)dxdw.
(24)
According 0x(x) = kx(x,-), we can treat (x as (x(x)(-). Therefore, we have b(x',t) == J p(v, t)xx (x, x' (x, d)dd and (W, ()(z' == Hw, f)kx (z, x' )p(r, w)drdw for all x'.
Thus, we can treat PMCR as linear ill- ill-posed inverse problem the RKHS by the operator To ensure that is a bounded linear operator, we require some standard assumptions (Zhang al., 2020; Mastouri et al., 2021):
Assumption C.2. 3cy oo, Y < cy a.s. and E(Y) cy.
Assumption C.3. (i). (x, and kw (w.-) are continuous and bounded, i.e. there exists x >> 0 such that:
sup (low(w))lxs (w)l <K, ( (x)]lxx S K. w
(ii). Feature maps (W)) and ox (X) are measurable. (iii). ow (W) and øx (X) are characteristic kernels. Assumption C.4. The kernel kxx (z, x²) is integrally strictly positive definite (ISPD), i.e., for any function f that satisfies II (1(2)((())))) < 0o, we have Jʃ (x)) (ar, f'' ((x)dddrr' > 0.
By assumptions C.2 and C.3, b(x, t) E Hx and A is a bounded linear operator from Hw to Hxx Based on the above formulation, we can rewrite R(H) of Lemma C.1 with regularized term as follow:
(H) = [|Mx, t) AH(w, f)(z) lix + (l) f)
(25)
Plugging the estimates of B(a, 1) and A into the loss, we have (23). Based on the ii.d. samples (xi, wi, mila and dx (a:), the estimates A(-,t) and A are given by:
b(x, t) = 1 nE Σ p(u. (x,,,, AH(W,t)(z) = Σ
(26)
In the following, derive the minimizer of the risk (25).
Let A Hx - Hw be an adjoint operator of A such that (Au, u) Hx = (u, A u) for all E Hw and u E Hx. And denote A an adjoint operator of Á. By (Mastouri al., 2021), for any m(w, t) Hw, have:
A" (X, = f
m m(a, k)kw (w, w wx) ws wddddw.
(27)
Since (w((w) = kw (w, :). we have A"1(x, ()(w)) = given by its empirical form:
:of===mmm fm(r, k)ww (w, w p(r, w) for all w'_ The estimate Ā is A ('(() = m(z, k) (w;, w').
(28)
C.6. Ill-posed inverse problem and solutions
Solving R(H) is generally an ill-posed inverse problem, as it may not have a unique solution (Carrasco et al., 2007). We allow the Conditional Characteristic Restrictions (4) to be ill ill-posed and have non-unique solutions. Thus, the set of all solutions is given by Hw.o = (H(,t) E Hw : (:, b(.,)) HD (. + Ker(A),
(29)
22

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where (er(A) ( (H(())) ( 0} is the null space of the adjoint operator A. This solution consists of parts, one is special solution H. Ran(A), and the other is the elements of the null space.
If the solution exists, we can express the solution in the form of the singular value decomposition of A. Let (Aj, (j. be the singular value decomposition of the operator A. Then, if we define the orthogonal projection operator Q H. + Ker(A), (A), we have:
(U( - I (0)
H(::)) - -(,t)-
Thus, we target at the special solution (W,, (), which achieves the least norm, i.e.
H (W,t)
arg min [W,W, H(W.))eW.oo
(30)
By solving for (H) of Eq. (25), we attempt estimate the minimum regularization solutions in respectively the population and in the finite sample regime:
solution (W, in (30) via the Tikhonov
(W,
arg min (x( {(A + (AA'A-A))
(31)
(WW
arg min RA(H) = HWWA)EHW
(32)
23

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
D. Proofs of Asymptotic Properties
In this section, we study the asymptotic properties of the testing statistics A p.m. Since A p.m depends on Tw (8, t) through (7), we first study the asymptotic properties of T, (s,
Notations. For generic random vector W E W, we use (F(())) to denote the space of square integrable functions of W with respect to the cumulative distribution of W. For any s(W), E [Fw)), we denote the -o-morm B(S(W)((W)). We use Hw to denote the
by [[ = / [B[[W)3 and inner product by (J.9)C2(F(w)) reproducing kernel Hilbert spaces of W. For any f(W), E Hw, we denote the Hw-norm by || and inner product by (fr g.w We let IP denote the probability. We let ([f(W)} J (())dP(w) be the expectation with respect to W alone. We differentiate this from [[](W)), which we use to denote full expectation with respect to both W and data W1, ..y Wn- Thus if H depends on the data w1,..., Wa. then P(F(W; } remains a function of H (and thus the data) but [Wf(W; H) is nonrandom scalar. We use both P,, to denote the empirical expectation with respect to W given data W1, ... W P.. (W)) = ; Ei=1 f(W.).
For A, be(w) = b(w, t) defined (24), and A" in (28), the estimates Ā, (.)) == B(w, t) given by (26), and the estimates A is given by (28). Besides, for the operator A, its singular value decomposition is given by (An, Pn, on) We denote H = (w, as the least norm solution is defined in (30). The population Tikhonov regularization solution (w) == (w,, t) and the empirical Tikhonov regularization solution (A(w) == (w, t) are respectively defindd in (31) and (32). Further, recall that
9s (im(X,s)on (W)) (assumption 4.11),
(33)
U(W Y,t) = ((Y - HD (W, (section 4.2), Y,t) (W,YYtt
(34)
(Y, / (W, () 4.2).
(35)
D.1. Proof roadmap and key assumptions
In this section, we introduce the overview of our proof and the required assumptions derive the asymptotic distribution of Ta (s, t). Define UWW, = =)(Y H (W, t), we decompose T, (s, t) as follows:
1
(s,
Ü (wi, yi t)m((a;)s)))
Vn
vnP,
((W.YYt)m((Xs))
VnP, ((x(Y.0) - AM(W,1) m m(X.s))
= -Vil,, [(P(Y,t) (W, + (W,,)) -W,
(36)
= VnP, ()) + VnP (H ())((X))
Expected risk difference
+ Vn(P, - -P) (H HWW (W)()
Empiriaal process
To derive the asymptotic distribution of (s, t), we first investigate the last two terms (36):
Empirical process (Proposition D.3): (Pn, - IP) [1))(W))))) (A(W, = (1122)))
Expected risk difference (Proposition D.4):
VnP (F)P:..))-(AW..))mmx..... vn ==1
Σ
U(w,, ()) + (1),
where Is is defined in (46).
24

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Lastly, show that-
1
v Ei=1 U(w,, ()* O)AAAs 9. plus the remaining term (U,
converges to the zero-man Gaussian process G. i.c.,
1 Σ U(s,, Mi, (AA(AA))) -ta Ga.--
lim VTPP R-000
Since ppm maxieT Is [Tn(s, (Tn(ss dp(s) in (7), we show that 4gm converges to max.T f |Gs,{Pdu((s) in Theo- rem 4.12.
Before proving these properties, we first introduce some regularity conditions. Let Hw denote the function space such that (W,t) E Hww for each f.
Assumption D.I. Let N1 (e, Hw,ll II )w be the bracketing number of size e of Hw- We assume Io' Nile, Hw. II (Ilcca < oo and P(AAAA Hw) 1.
Assumption D.2. Let (Λj, (Py: ob]ر be the singular value decomposition of the operator A described in section C. Then we
assume: (a). For some 17 > 2. Σ; ,3 -- (g., +Py) < oo; (b) For some 0 2 2. E, > 28 (H)) HP))HW < cc.
Assumption D.I restricts the complexity of Hw and ensures Hw is P-Donsker class (vd Vaart, 1998), which was standard assumption to analyze the empirical process (Beyhum et al., 2024; Lapenta & Lavergne, 2022).
Assumption D.2 is the source condition that is commonly assumed in nonparametric regression (Carrasco et al., 2007; Florens et al., 2012). These have also been employed in (Florens et al., 2012: Beyhum et al., 2024) to obtain a faster convergence rate for nonparametric instrumental regression. Here, we require 9s and H; to satisfy the source condition, for establishing the asymptotic properties of the statistic in examining the integral equation. Since gs == Emm(X,, () (W)). the source condition for Is puts requirement on the smoothness for the space Hw when m( s) is chosen properly.
Compared to (Beyhum et al., 2024; Mastouri et al., 2021) that allowed
2. we require 17 > 2 in the source condition for 9.- This is because, after investigating the proof of (Beyhum et al., 2024; Mastouri et al., 2021), we find that their conclusion may not hold for = 2. Specifically, similar to our analysis, a key step in their proofs that related to the source condition is II( + A A)- A A* Ag- - 9l Ho = (1) (in our scemario, it refers to (51)) for ğ E Hw. Here, 9 = Σ3 >o (gg PP){P Ø refers to ge in our scenario), where the source condition was assumed on g (g. our scenario). To prove the property. they
leveraged the property in Lemma 2.5 (d) of (Beyhum et al., 2024) that [[(AI + A) A' A Ag g Kw =0p min(
for ÿ such that |0 W. := Σ, 2 1 < co. Since g satisfies the source condition and . := Ej 2 (g. 9;) to apply the above property, we should take S 17 2 to ensure that ||]| W., < co. That means, we should take r > 2 to min(s.2) = (1). Besides, while this condition is stronger than the case of 1 = 2, we emphasize that requiring 1) > 2 only slightly increases the smoothness requirement of the Hw, making it a modestly stronger condition.
make that 7 1) 2 > 0, in order to ensure that O, x
D.2. Empirical process
Proposition D.3. Under assumptions 4.9, C.2-C.4, and D.1-D.2, the empirical process VT(PP (W (W, s)] = o,(().
Proof. We first proof - ( 3) = (1). ). In fact, we have (W.)) - (x(W.1)3|m(xxx))) (IPP(W, = | (W) -
(W, - FM(W..))m(X.))
_ -
110 (-1))
where (1) follows from assumption 4.9 and (2) follows from (21) by assumption C.3. Since assumptions C.2-C.4 and D.2 are satisfied, we have (w, t) (w, E) = (1) by Lemma 1.17. Therefore, all conditions in Lemma 1.14 are satisfied and we obtain
VT(P - (W, (W, 0) (X,))) = (1).
25

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
The proof is completed.
D.3. Expected risk difference
Proposition D.4 is our main result in this section, whose proof is decomposed into Lemmas D.5-D.9.
Proposition D.4. Under assumptions 4.10, C.2-C.3, and D.2, the expected risk difference term has:
n
VnP ( (FW,,))) ((77,))m(xx8)) m(x.)) =- Vn
Σ
U(w, Us.
(x.) + o1))
Proof. For simplicity, we denote Hf)w) (w, f). b,(x) = b(x, t). Based on the interpretation of PMCR as linear ill- ill-posed problem and the form of Tikhonov regularization solutions in (31) (32), we have the following decomposition (Babii & Florens, 2017; 2020):
6) ( + Gs:
(37)
where
-===AAT A- AHH));
(38)
G2 ->== - A-) - AHH);
(39)
G3 : {(T - -AT AA Ar ĀH!):
(40)
G4 :4:=(A7+ + -'A)-- - (N A) bb:
(41)
G3 OA:=(A+ A'A) H
+
(42)
Therefore, we have
(WW.))))) = Spi (s, t),
1=1
where Sni(s,) t) is define (TP(( s)}. By applying Lemmas D.9, D.5, D.6, D.7 and D.8 to S,,(s, t), (s Suz(s, 1), (n)(s,)) () and Sns(s, 1), respectively, we have:
-WW.:)))))WWWW.((((...). =
Vn Σ
Uww, () + (1)
The proof is completed.
Next, we provide proofs for Lemmas D.5-D.8.
Lemma D.5. Under assumptions C.2, C.3 and D.2, (n2(s,t) = (1()) +
Proof. By the reproducing property that f(w) kw (w. Hw for each f E Hw, we have (AI + A-A)-1(A (b, (w)) = ((AI + A' (Ā* Ab, - AHς), kw(w, :) Therefore, we have the following for (, = VnP(G2m(Xxs))
IP{G2m(X,, s)} = E (AT A) (ā' A ÁH!) m(X, s)
RR. (E(( A- - - -AH!), ow (W)) Hw m( (X,
((AT+ A+A) -1 -1-)6. -AHT), E(m s) (W) Hw
((T) -4- ÂH?), (NT A (X, s/w (W) H
26

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where (1) follows from (20) and (2) follows from {(A + -)= (A7 + A A) in Sec. C.2. Thus, we have:
(1) - II(Ā* (b; AHę))luw |1 + A A (mm))
[P(G2m(X, 8)
(E s)o (W))
< || lopp 11 - 1 [I(AI
s)o (W))
|| - Allop 110g - AHH I|(AT + A A) E () s)w (W) (W)
where (1) follows from the Cauchy-Schwarte inequality and (2) follows from || = ||Allop in See. C.4. By Lem-
mas D.12, we have 18g by 1ex Op(1/vn) and 1|.1 Allop Op(1/vn). With these properties, we have:
1:. AHH Hx llb b, + b, AH, AH, AHę
|b, b, + +AH, - ĀHς 1ax
(43)
I 1b, ÂHд 1lxxx + ||AH, AHl
1h0 bellaex + 1|A Āllop || Op(1/vn).
By assumption D.2 (a) with 9s [mm(X, s)@w (W)), apply Lemma D.10 (d) obtain that II + A (m) (W))l = =p(1 m -.. Combining all the inequalities, we get
Vn (P(G2m(X,))) & vn. Op Vn (21) O,
Vn) (1h)
OptA min(m
(44)
Since we require r) 2 2 assumption D.2 (a), the last term (1))
Lemma D.6. Under assumptions 4.10, C.2, C.3 and D.2, (n,(s, = (1) as n 00.
Proof. By Lemma D.16, we have:
G3 = -
+ Ā-
A AA AAH
= = +
A A)(AT + Ā) - (br AH
- By the reproducing property that f(w) = (f. kw (w; Hy for any f E Hw, we have (AI + A) (A'A- Ā'A)(AI+ + A -) (b. - (w) = ((A 00-+44 A' (A"A A Ā)( + A A) - (b. ΛHT) kow (w, -)) Therefore, we have the following for Sna(s, 1) = VnP(GGm(XXx)
A
IP(Gsm(X,))H = (AI + A'A) 1' + T.. (H))
E ((A( AA' -------- + -AA - -AHH), dow (W)) Hw m((X)) 8)
((XT + A-A) T AA ( - - -ΛH!), A k)dw (W)
@ - -H, Á(XI
Ā AAA) +A' -Bm)) s)S (W)
where (1) follows from (20), and (2) follows from ((V A) (N + A- and {(V + Ā+ A)-
(AI + A' Ā) Thus, we have:
(P(G)m(X,s)) 2 b. AHPIIN |IĀ(A) + - A A) + A" -Bmmm (1)mm(X,)))n)) s)ow (W))llxx
b. AHllIN |IĀ(A/ + Ā -' ||A' A - Ā Alop [[ + A" A)E s)o (W)
where (1) follows from the Cauchy-Schwartz inequality and (2) follows from operator norm inequality.
According to the paragraph above equation (97) in (Mastouri et al., 2021), Ā a compact operator. By Lemma D.10 (c), we have ||T Ā - Op(1/V),, |I Ā Ā) - A" 1lop (p(1/VX). According to assumption D.2
27

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
2 -]. Finally, By Lemmas D.12, (a) and Lemma D.10 (d), we have ( + AA* (X))) = Op(1 mincn 2)
we have 18g- be 1xx Op(1/ Vn) and ||A Allop = Op(//nn). With these properties, we can similarly obtain that
[b, AAH |xx I I|b b + || Allop H |k = Op(1 /T). Combining all the inequalities and 1) 2 2 in
assumption D.2 (a), we get
Vn/PpGsmm(X,))) = vn.
(
Op (12) Og
1 (vT) (V²)
OppA
minmn 11
vn)
(45)
0g
Oy(1).
Vnd
where (+) follows from assumption 4.10.
Lemma D.7. Under assumptions 4.10, C.2, C.3 and D.2, Sna(s, t) = 0,( as +0.
Proof. By the reproducing property, we have {(M + A+ (-A)-11A A AHε (WI + A) b,, kw (w. Therefore, we have the following for Sn1 ts;
Anl (N A'A)- Abbl)) = (((A) +
IX)Gam(X,))] = E [(N + A'A) A' Anς (N + A'A)) m(X, s)
[(AA+ A A) A AH, - (AI + A) 1 b ow (W) mm s)
((AT + -
(AT A A" 1) E(m(X (X, s)w (W)))
where (1) follows from (20). To establish the upper bound for the right-hand side, we derive the form of the operator A A. First, for the operator A Hw - Hx defined in (24), its singular value decomposition given by (Anr (n; tn) For Is (w') == E (r,(X,a 8) k)w (W) (w'), we first proof ||. < oc. In fact, since (m,(X w)]W} and kernel k are assumed to bounded in assumption C.2, we have:
119.11ze [IE(m(X,)) s)ow =
How (W)) How
S CC) = (L Hów (W)) Hw
(46)
-CVB (@w (W), dw (W'))nw C CVE/fw (W, W')} 00.
Thus, I1g=1lww < <0 imply 9. E Hw. which can write the expansion of the basis functions as gs = =E, (9) 9j by sec. C.2. Besides, by the properties of singular value decomposition, we have Apj \jos and A A*Ø; Ajpjj. Therefore, for each 9, Hw. we have:
A" Ag. =
Σ Aa (s.
(47)
Let (u(")) = =E, \J² (go,
2g (5)) Py. By assumption D.2, we have:
1 1/
Σ xj Σn
Pg.. Ps. Σ ΣN₁² 2 Pp)mwep
Σ " Py) Fa R< co.
(48)
Hw
Besides, we have:
ArA = =ΣΣ 232 2 (0. (p:n)))) =
Σ (g.. sssh = Is.
3
Thus, if we define P. = (NI + A² - ĀHH - (AI + be, then by the Cauchy-Schwartz inequality, we have:
IP IP s)} = (Pt, ()m((X)))) (W))) Hw KPP,, A" Ag.)Hw
<
(Pt. (A" - >A*)Aga)) H.
+ (Pr, A Ag.)
(Pe, (A* -Ā*) -A')Aīa) Rew
+ (APt, Aga) He
|| Pelleew -1 - Allopp
|| + |1 -11A0lllxx
28

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where (1) follows from | A'loop | Allop- By AH, by, can decompose Pa as follows by Lemma 1.
P. = (AT A A) An, (AI + A A)- AHÇ
=((I++ A' Ā'A XI) (AI A) (N--AAAA An) H.
- (M ) H²
=X(M + -AA) A 4AA + AH
Therefore, we have:
Poll = ||X(M + A - (Ā A'A)- H| Hew
S ||A(AT + A1 ||A*
Ā Allop II(AT + A) 1l
Since A is a compact operator as shown in the proof of Lemma D.5, by Lemma D.10 (b), we have [1AA(1 + -) min(
1 op S 2.
By assumption D.2 (b), we can apply Lemma D.10 (d) to obtain that || (M + A A' l Aw =pp1 2 -1.. Finally, by Lemma D.13, we have ||A* A A Allop = Op(1/vn). Combining all the inequalities, we get
[Plll 2. (1²)
min(O.2
Next, we prove || < co. In fact, since A is bounded linear operators, we have |Allop < co. Besides, by Eq. (48), we have l[ < oo. Thus, we have |A $ ||Allop ||āallww < oc.
Next we provide the bound for ||APollxx
|ÂPollux |ΛĀÇ(N + A A*A) + - Mx
SA ||Ā(AT X+ A) | Ā - A^Allop II(AT + A A) - Fl
Since À a compact operator, by Lemma D.10 (c), we have |IĀ(A/ + A -1 lop = O,(1/VA), [I A -- -1Al (1/VV)). By assumption D.2 (2), we can apply Lemma D.10 (d) to have [(A + A) - [410l1z || OptA
min
Finally, by Lemma D.13, we have || A Al (p(1///)) Combining all the inequalities, we get:
1 )1)00,,
| A Op
Op
min(6.2)
Therefore, substituting all inequalities to ts) and A 2 2 in assumption D.2 (b), have
vn |P[Gam(X, s)H & Vn. |PPlll 1| - Allop ||A + /n ||APPll || I
min(6.2)
{ x-2-1/2.000 Op
1
min(6.20
< Tn. 12.0 Op
Op Vn + Vn.
= Op + 1/2 (1) 1/2
(49)
By assumption 4.10, the last term is (1))
Lemma D.S. Under assumptions 4.10, C.2, C.3 and D.2, Sn,(s, (1) as n 40.
Proof. By the reproducing property, have {(T A A) A)-1 A HP) (w) ((AT A-A)-
Hi kw (w, H) Thus, we have the following for Sn5 (3,(s,)) f):
[P(G5m((X s)}1
E [[(AT A A) -1 A H.' ml [X,
=|E [E(( + A A) -1 AAb, H. ow
ow (W) Hm. -m(X s)
[(AT + ((A)+AAA)* A - H E(m(X, s)ow (W)))
29

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where (1) follows from (20).
By assumption D.2 (b) and All, b,, we apply Lemma D.11 to obtain that |I(AI + A'A) 11
+ A'A H lww = OptA min(d.2)
Combining this rate with the Cauchy-Schwartz inequality and and
2 2 in assumption D.2 (a), we have
Vn [P([Gmmm(X ((A + A'A) A H1. B/m(XX S [( + A) "A'b; 11 Hw 1ll.sll H
(50)
2 O nx3) 2 o,(1),
where (1) follows from | 00 by Eq. (46) in Lemma D.7 and (2) follows from assumption 4.10.
Lemma D.9. Under assumptions 4.10 C.2 and C.3, (s, = + co, where g-(-) (m,m(X, s)kw (W.)).
1 v -.
+ op(1) as
Proof. By the reproducing property, we have (M + A AHt)(W) ((XI + A'A)-1A- AA. AH!), ow (W)) Therefore, we have the following for (s, t):
P(Gim(X, s)} = (N + A-A)) - (m)(ss)
=B -E(((A++ A-A) Ar - -AH,), ow (W))mw m(X, s)
(IN A' AA"(6 - AHς)) Bmm(X, (W)) Hw
@ (A* (b. AHς) (AI -'A
(A AHP), {(A + A* (A* -) + (A~ - ÄH?), AA 9g))
where (1) follows from (20) and (2) follows from {(V - A- (AT A
We first analyze the second term in RHS. By (26), we obtain:
p(y, txxx(x;) - n
(b - AxT)(x)
Σ
Hun (x)
Ulw,, Vi, ()xxx ()xx
where U is defined in (34). Since A'm, == f ()X W) in (27), we have:
{ (r-(6. (HH3) ()
1
Uwwi, ti / kxx (W)dP(X, (),W)
n
1
Uww, t)A (2x
30

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Therefore, obtain:
Vnl (r AHH). (A*A) 90) Hw
=Vā =V (A'AA)-) E(m(, Hw
( Vnn
(A*A)- - -
Uww,, Vc. t) (xx ()x4 m
=VnE
(A*A)-
Σ
Uww Vi. 1) / (x)
(Z,,
W)
vn
12)
Uw, Vc, t)
(2x, (xx, ())W(W))) Hw
Uww yr 1) (kx( -)
'9)) Hx
3
Uw, yr,
where (1), (2), (3) follows from reproducing property f(x) = (f,kx(z, H2x and g(w) = kw, H)w for each
e Hx and g Nw. Besides, for A (b. AH ]l we have:
A Al = | - A + A. - A AHH
1 |A Abbelluw + 11AAA* AHE A AHH New
I | Ilop -10, bellxex + lopp ||A Allop |1
Since H. E Hw, we must have ||l|] Hw < xc. Besides, according to Sec. C.2, we have lop ||Allop < co since A is
bounded linear operator. Therefore, the last term is Op (+)). By Lemma D.12, we have /T|lAA (b.-AA)))) Hw = Op)). By the Cauchy-Schward inequality, we have:
(VnAr AHP), ((AT A² ()*'))-1))) Hw
S VTl (b. - A ||(((AT + -- - (A*A 9- gell
s Op {II((AT+ + A'A) - A) -)'-))''1 1. }glll
-AAA)-
By A Ag. 9. Lemma D.7, we have:
{(I + A -1
-AI (A'A)-')g, -1AA))1--AAA----
A A" A) -..-
(51)
-1 "A* is Tikhonov regularization scheme, the RHS of the above display converges as A- In fact, according Lemma D.11, we have [[(X1 A) Ag gllir.. S Amiu((,2) since g satisfies assumption D.2 (b). By (48), J. = E, X, 2 (ge, Ppshnw Pp. Thus, similar to the proof of Lemma D.11, we have
Since (AI A)
II{(AT + A-A) -AA "gall
Σ
(g.. Py)
X3
+X
X3
A1
Kg (ps) 2
Σ
S suup j
13
+
If 17 2 4, since the maximum singular value of the operator equals ||Allop < oo, we have
AA
}
X²
}
sup i
x5 + A5
<
s
X3(22
+ A)
S sd'supagn O(2?).
2n
31

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
2 < 17 < 4, we define a := A2. Plugging it into
x2(2/+A) +x) ʃ we define f(x) = x² Noted that f(x) is maximized
AA;
(z+
(by using the first order condition) at x X(2 (1) 4)- Thus, the maximum value of f(x) is
x²² S
(g) - 2) -)) X -
(x + X-2
4
Thus, {(M + +A-A)- (A* (0-A))113 19. - by assumption 4.10. Therefore, we have
1
() o,))
v
D.4. Proof's in section 4.3
Theorem 4.12. Let (s²(0)= == U(W,Y, (1)m(X,)) -U(W ()AA with O (W,Y,X) X). Suppose assumptions 4.9-4.11, C.2-C.4, and D.1- hold. Under THl,, we have (i). (n, t) converges weakly to G(s () such that (Gs < oc, where G(s f) is Gaussian process with zero-mean and covariance:
() (8, (B(ms..O))m (O')),
where O' = (W", Y", X') is an independent copy of O; and (ii). Dgm converges weakly to max (()d(())))
IET
Proof. By (36), we have
(s, VnP, (UU, (W) + (Expected risk difference) + (Empirical process).
By Propositions D.3 and D.4. we have:
(x, = Vn Σ
1
Uww, (6, [m((r,,)) (A(A*A))'s)) (#:) + 01).
Next, apply Lemma D.15 to {U yz. 1) {m((r;,)) A(AAA))= (x:)),, obtain the result. To this end, need to verify UWW Y, (()AAA(()"AAAA)) 9. + m(X, s)] is zero mean and
E [10 Vi, () 8 - (A(A A) 9 ( < <c.
(52)
Notice that the zero-mean is met by ((U( sides, by assumption 4.11, we have Var(U(w,,, Vs, t) s)) - (A = E(U vs, (m,,,, < do for any (s. Therefore,
B[p(Y,t) (W,
0 under Ho. Be-
{A(AA")) A)- 9. (xx)))² < E(U(w,, Vi, [X, + [[m s) exchanging the order of integration, taking the maximum value of f, and setting the measure to be probability measure, we get (52). Thus, we have (s, t) converges weakly to G(s,t) t) in C2 x TT x μ), where G(s, t) a Gaussian process with mean.
(A(AA"A)''g,)am
For any fixed and (s, (T 11), the continuous mapping theorem (Theorem 1.3.6 of (Wellner al., 2013)) to obtain
/ [Tn(s,t)Pduus 4 /
by the continuity of the integral functional. Next, for [T, (s. dd(ss E (T mJ), noted that t in (s) is determined by Uww,y, = p(y. t) (w, t). To ensure that taking the max operation is meaningful, we need to prove that if U (w,, y, 1) E (F y)} for any 1, max[U(w, y, f)) E L (2 By (17), we have: iET
/ (p(y,t) - H(w,t) p(w, yddwyy 2 / Ip(y, p(yddy 22 / H[w,, ()
g(w, u) 1 [/]
p(w)p(u)dddd } KV
et p(ulu)dd) p(u)du } < 00,
2
< 2
p(u)
32

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where the second inequality follows from (a - b)2 S 2n² + 212 and leavp= 1. Thus, taking max operation on both sides, have max J (W, y, p(we, y)dwdy co. EE
Next, prove the continuity of the max functional in metric d. Next, prove the continuity of the max functional. If f2) 8 given any > 0, we have max fi [t) - f2ax)f((t) I max)ff((t) fa(t) d(fi, fe) < 8. Applying the 1E2"
tET" continuous mapping theorem to such continuous metric max, have:
/ [0(s,0)du((s) EET
4 A (1)
max IET
max
The proof is complete.
Theorem 4.14. Suppose assumptions in Theorem 4.12 hold. Besides, we assume BAr(X, ())) < oo. Then, we have:
(i) Global alternative. limm-+00 maXIET (s, t)] for almost all under Hfx
(ii) Local alternative (a < 1/2). lima-+o0 maxTeT [s, oc for almost all s under Hin.
(iii) Local alternative (a = 1/2). (s(st) converges weakly to G(s,t) + p(X,s,t) such that (G(,tt + < oo under Him where G(s,t) is defined in Theorem 4.12 and p(X,s,t) = E (r(( - AA*A) (XA ())(XX)
p(X, (dtt){pp((ssd((((
Proof. (i). The case of HHf
We first decompose Tn (8, 1) as follows
1 "
(s,
vn Σ
(ws, Mes t))((;; s)
vnl,, n
[s((Y.)) - FA(W.0)m(XX.)
= vnln
(v((K,t) - (W.,)) - (FWW.))--WW..)))MC
((p(V.t) .) - -F. ()) - VnP (FF (W,,)))) s)]
VnPn
-VP(P, - [( - -FV(V.))3
where (W,,)) t) is the least squares solution.
We first analyze IP/ [ - (W, s))
For (w, 77 (w, 1), we can decompose b, bran + bker, where Weee Ker(A*) and b, -er(AA)) (An(A). By using a by E RRa(R A) to approximate brcm, we get be bran bf + by + bker. By this decomposition, we have the following for Aw, t) (w,,)) t):
(w,t) (w,+
(w, t)
(AI + A) At, (w,t)
= =(T + Z Ā) -1 (.. - b.) + (AT + A Ā) - (A"A)-1 AAb,
= =AI Å A)- -Aat -b,)+ +(X+ + - (d''A) (A"A)-'A"(8b" ाााा + b, kor)
(b"
= -AT+ Ā- -1 -b,) + ((AT + - - (A'))--"))0 Dj + bj 1
Thus, we have
[[AA(W, (V(W,)))((X.)) = ([[( + - 14)06=
- -13)
m(X,(XX
+P ([(M A' -('AA -A'A))))-A'' m((xs)
+P ([(A++ -'A) (((A++A'A))+ - -A '4776 m(X,s) 8)
+P
[(AT + -) - b))m((X)))
33

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
For the first term, since we can always make lf arbitrarily close iran which causes ba bf < E. Thus, we only need consider the last three terms.
For the second term, since bt E Ker(AA), we have
{(1+ A' A) Ā- (A*A) = (AI + Ā) A.
Then, we have
P { A- - Att m(X. {s}
(1) ((AI + A
[m((X, s)ow (W))) Now
= ((AT+ + A' - A A' -gs + ((AT + A A) 11 A ,9=) Hw
2 ((π) A, (AT + Ā- Ā-)) g)) Hw
((A' - A 1r ,((N + A A) (AI A }ge)ne
(1)
+ (AA Ar) .NI + A) A)
9s) g,)Wk, Hwy
(II)
where (1) follows from (20) and (2) follow from bler Ker(A").
According to Lemma D.16, for the first term, we have
(I) = ((Ā" Ao {(N + A) -1 Ā- A (M + "
By Lemma D.10 (d) and assumption D.2 (b), we have ||(AT + A" I1 -Op(1 min(n.2) -1). Since rj > 2 in
1).
assumption D.2 (a), we have ||(A + A A) " HWw ())) By Lemma D.12 and D.13, have A A
(p(L/Vz) and || Ā Allop Op(1/ Vn). By Lemma D.10 (b), we have II(AI + A A)loop Op(1/A). Combining all these results, we get
(T) $ Ι - A 1lop ||b,or|| |1(AT + A' -) 1 Ā Alopp |(AT + A'A)
-p(1/\m) (,(AA) * (p(1vv)) * (1( = Op nx
By assumption 4.10, we have (I) (1). For (II). by Lemma D.10 (d) and assumption D.2 (b), we have
[[Y + A-A))" Op min(m,21 -1 j. By Lemma D. and, we have A| A lop Op(1/vn).
(II) S 1 Allop !|b wm ||AI + A' 9: Hw
= =pp(1/Vz) Opt1 min(n.2) -1) op(1).
Next, we consider {(A + A-A) 1 (A* - In fact, we have the following decomposition
{(AI + A' A) 1 - S₁ + S2 + S3,
where
S :- = ((AI + A'A) - - (NI + + A) 3b:
S2 (AT A'A)- - 1''b;;;
S3 : {(V+ - A'D
For S2. we have
P[S2m(X, {)} = ((AT + A) -1 ' - A B(m(X X, s)ow W)))R Hw
-(LĀ- A''), (AI + A-AA He
< 7 - A'llop | 1 II(AT + A) -'
34

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
By Lemma D.10 (d) and assumption D.2 (b), we have |[(A1 + A-A) ge Hw
= Op(1 m
2
-1.. By Lemma D.12, we
have || - ollpop = Op(1/vn). Combining all these results and according to assumption D.2 (a) with 7) > 2, we get P(S2m((X s)} = (p(1/vn)== op(1).
For S3, since A" "Ag. Is. we have
[[ {,,,) 8)}
-((AA++ A* AAA)" A'G (m,(X s)ow (W))
-A {(A A-A))=--A'" } g.)
=A
A Ag.)Hw
= (A (AAA A''Aaa - Oz) Hw.
By Lemma D.9, we have I1(AT + 1)=1A) - TallHw = (1) Besides, since ||4" = ||Allpp < oc and ||b ||x < 0o, we have
P(Sam(X.s) s}} 1 ||Allop ||bEllmx II(AI + A) Telluw = op()).
For S1, we have
P[Sim(X,8)] = ([(AT --Ā) - (AT - 4 A- Bfm(X,((XX s)dw (W)))
=([AM + - (M (f
(A')
be. ge)
= ([AAA- - Ā)T + -) 12 (AT + A A) g.)
$ ||A'A Ā' II(M + A) 1 A' 1lop || 11(A1 AA |
By Lemma D.10 (d), we have [[AT + A- ys WN Opp1 minfm.2) -1). By Lemma D.10 (c), have |Â(AT +
Ä A) A)- - llop = (I + A A) A llop 2 1/VA. By Lemma D.13, we have ||A A A Ållop = Op(1//n). Combining
these results and according to assumption D.2 (a) with 17 > 2, we get (PS[mm(X,)) s)} Op(1/VnA) = op()..
Thus, we have
P (((( 4' (1-A) (A-A)-AA4666 m(X,s)) op(1).
For the last term, we have
P P|((Ne+ A) -b))mm(s) = -(AI + -A- A (b. b,). (m(, s)ow (W)HMHw
$ <|(Ar -² -'A")- || -16 bell ll
By Lemma D.10 (c), we have |I( + - S 1/vă. By Lemma D.12, we have 1 bellop Op(l//n). Combining these results, we get P [(AT + Ā² - -1 (6, m)) = =O,(1/VA>)= = (1).
Thus, we have
P((FA(W,)) FF = o))
By the above proof, we similarly have P - - = op(1), and therefore (Pn P) (f,( - ()),,))m((XXx } = (1).
Besides, for the first term of (s, t), there exists f, we have
Pa [[o(Y,)) (W,
P. [( (W, + HW, - (W(,.))m(Xx))
According to the definition of Hifi, there exists r(X,t) such that E(p(Y,t) - (W; = r(X,t), where r r(X,t) cannot be written as E (HH(WW,t) t) - (W,)))) for any
35

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
We first prove (((rr(X,)))))) + ()) F),,, s)]l < co. In fact, we only proof 8) < and (W) (W, (m/)) s)ll < oc.
For E))r(XX s]]} we have
3()((X00)/(X & |r(X, 1x(X,01l2 (|m(XxX 8) s |r(X, |m 1/2
(P(t))))
where the last inequality follows from (21) by assumption C.3. Since ||xxX,))lxxx <o and the selected weight function m satisfies [|mmX,, c)(())))()) < oo, have [[[[(X, t)m(X, s)]3 <0..
For (W,
s)I we have
B[((H((W,)) - (W
1 (W, (JW,))))((((X))))) - C- (W, - (W,
(C. - -F, (W), f)ll
where the second inequality follows from assumption 4.9 and the last inequality follows from (21) by assumption C.3, Since (W, - (W,,) E Hw, we have E (W,) - f (W, w, ()(( s) < <0.
According to the law of large numbers, we know that it converges to (r(X,t) + HW, t) - (W,,,)))((XX We assert that if HIfx holds, then exists 1 satisfies E [r(())) + (W,, (W,
* 0. Proof by contradiction.
r(X, t) = E (w) (W, x) holds for any t, we have:
E tY) (W, (W,) (W))))))
which implies E (,(Y, = (W,
(W, { for any t. This is not consistent with HiI
Combining these results, we have:
lim max [Tn(s, lim +1(0p(1)+00())) = co. n00 IET n-
for almost all under Hfx
(ii). The case of Hy, with 0 < < 1/2.
We first decompose the term into
(fFA(W,E) - ( into Di= Gi. where
G : A AHT)
G2 (N AA) -A') - AH!!);
G3 = ((T -A)+ (N + A² A- (b. A'(b= AHH):
G4 -A
-NT + A² A
G5 ( - -
Gg Go:=H² (w, t) 0 (w, t) = HDdh (w, t) - AA A'b
We define ((x,t) [(r(X,))(x(x));
36

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
For P(Ggm(X, 8)), following Lemma D.5, have
P(G2m(X,s) (X, s)} = ((A + A A) ( (Ā" - (r) - AHl)) ( gs)
(be Anl. (Ā - A)(MI A A-
(b. be + b, - -AH (A (A- (A-A)(AI + AA'A)-
(be be + AH, +((X,))/n - -AHØ, (a-A)(AT + A (1=A)'1'))
t)/X,,)/",, (Ā (A-A)(M) + AA ga) + -6.--b bb, (A A](AI A) '
+ (AHH - -ĀHŞ (A A)(AT + A)A) ga)
By Lemma D.12, we have (|6. bellxx Op(1/Vn) and I Allop Op(1/VTo.. By Lemma D.10 (d), we have
[[A + A A) A OptA 2 -1.. Therefore, we have
P[Ggm(X, s)} - (((X,t)/n", .A A) + A) go)nx + 16; bellmex -1Ā - AAlopp 1 + A A) 00'gellww
+ A Allop || -| Allop 11 A -9ll-w
< (e(X, ()n", "A* - AA + A) ga)Hx +0p(1/n).
For P(Gym(X, s)}, following Lemma D.6, we have
P(Ggm(X,s)) (((2 A (AT - (N+ - (H.)).9)) .9.) Hw
(AT+ A A) " -A-A
-- 7AT バー -A - -AH!), Hw
(br - b bb, AMP - A A) + A A)
+
(Be- - by + AH, + (XX,t)/n AH. A(M
Ā) + A
(e(X,1)/n", A(AI + (A'A- (A)) 9x)Hx
+ (b. -6-,,, bb, A(N + -A'A--AA 4'AA)) ga'ux
+ (AHH ÄH. A(X AAA)) A-Ā)(AT A) -g)
By Lemma D.12, we have 1|6 bellzex (p(1/v7) and || - Allop = =Op(1/VTt.. By Lemma D.10 (d), we have [I(AI + A A) we have
A) goll = Opt1 min(n.)) -1.. By Lemma D.10 (c), we have |IĀ(AL + Ā- "llop Op(L/VX). Therefore,
P(Gam(XX,s) S (((Xtt)\n² A(M + AA A + AA g,)
+ 106 - by bellxx
|IĀ(AI + A ~ ||A'A - A'A'lop |I(AI - A -
+ ΙΙ A Allop | Hw
||Ā(A7 + 7 A) 'llop ||A'A π- Alloop ||(AT + AA*A) gs Hw
- (e(X, t)/n", A(N + -A'A + +A"A)-' g,) + +0, nVx
For G4 and Gs, noted that
G.+ = (AI + A A ĀHς (N A'A) A'b + + A'A)
A'b - -H¹
= (XT + - ' ĀHς (XT + 4'4
+r/Vn) m + (NT A'A) +r//2/r///m/m -HO
(A A A) -1 "" Anς (N + A AT A AH + +AI A -) A" AHl HI,
where G2 and G3 are those defined in Proposition D.4. Therefore, we have P(G s)] = Op (1/n) + 1122/v1 and (C0sm(X00))
P(G5m(X (X, x)} = OA()).
37

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
For P(GIm(X, s), following Lemma D.9, have
P(Gim(X,,)) 8)} =((A+ + A) - AHH),, H.)
((A AHl), (M A" 9s)
= (AA (b. AH,), {(N + A) (A*A)=) g.) + (A- - AHτ), (A*A)-'
'9.) How
(1)
For (II), following Lemma D.9, we have
(Ar AHH), (A* A) gs) Hw
A"
1
clu (xx) 1 Σ H (1ne, t)xx((1))
(A*
n
Hw
A"
U (1, + + ( (2/))/// / ox (xx(x))
(A*A))
Hw
A"
(U (m, AA*A)
Hw
AA* r (z,, () /n /x(()) -
Hw
U(wi, ys. () 9x.)))
A"
(zi, t) // ox (x((x)) (AAA))'90 gg Σ 1=1
Hw
For (1), we have
(A - - (MM -A'A) g.)
(A. - b. + be - AHH), ((N - AA* A)- gz)
(Ar - b, AH. -(x())/="' ((AT + A A* - (A*A) Ho)
(B. be,
(A'AA) ( Hw + (A" - A (A*A)- 9 Hw
+ ()) (A 9m)
By Lemma D.12, we have 16, bellxx Op(1/Vn) and A - Allop = Op(1/vn). By Lemma D.9, we have ([ + A) (A^ g011 (1). Thus, we have
(H..-AH)). ((A + A-A) (A 9o) (o,1///)++((0011111,,,,) ((A -1 (A* (A* A) ga) 9a)Hw Hw
For P(Gem(X, {)}, we have
((A" - AAb, Ho = A(A"A)-'A' (AHς +((X,t)/ax) -Ho g-) Hw
=(-" (I(),t) (X)t))n)),,,. 9.)
38

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Combining these results, we get
P[Gam(X,))}
1
Σ
Uw,, vr. ()AA" + Vn (e(X,0), - A7) + -A
A) "9)px
Vn 1=1
(1)
Vn t) + ' A A) A
+ n
AA'A
(II)
7
+
(A
((, 1) {1V + A- - (A* g.)
a
(II)
+
12
(A*
1
Σ
(2n)) ) (A*A))'9s)) vn na
(E(, (A* (A""))-'g)))
na
By Lemma D.12, we have I0+ bell = Op(1/vn) and || - - Op(1/vn). By Lemma D.13, we have
||A* π A (1(1//n). By Lemma D.10 (d), we have [I A' g) =0p{A min(n.3) -1}. By Lemma D.10
(c), we have |IĀ(AL + A A) 'llop = Op(1/VA). By Lemma D.9, we have { + A) (A* ' gall = (1).
Thus, for the term (I), we have
Vn (1) v
|((x, t)]l -1Ā" - 1op |[ + A) A''A)"'gallluu goll
na
Vn
Vn
n Opp1///) 0z
n2a
For term (II, we have
(II) S Vn
|I((X, -11---I---IAAL | Ā) ||A* - Ā' AlAlo II(AT + A"A)=""alll 2
n
Vn
Op(1/Vn)) op
na
n
For term (III), we have
Vn
III)
|| [|((X, t)ll -1 (1((M + A'A) (ArA) 1)
n
Vn
(
/nn
(1(1) oz
n
n
For the last term, we have
Vn
Vn na
A
r (xx)))((x))
(A*
(A" (I( (A* ()-"A)'''00 -9)) -)
na
71
Hw
vn
Σ r (z..))
E[r (X)) A(A* A)- ge
Hw
By Lemma 11 of (Mastouri et al.. 2021), we have r (X, xx (xx]] Op(1//p) Com- bining all these results, we have:
Vn2
n
s))
Σ
U(w, vi, -) + +0p n
+ o1)
(53)
Vn im1
39

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Besides, for the first term of (5, (), we have
Pn - (W. s)] Pn [
( (W, -FW..))FF(W.)))mmmc (WV
Therefore, combining (53), we have
Tn (x, vnP, (U(,YY,) -
{A(A* [r (X, 8) /n (W, (V - (W,
+ vn
s)] + 0p
+ op(1).
n
For the first term, according to Theorem 4. 42, we have vnP, (U(W, - 9,. converges weakly to C(s, t). For the second term, we have
VnPn fo (X, /n" (W, (W,,)))m(Xxss) VnPn [Tr (X, 1) / (1 A'
| no (Ar=A (1)
vi
According to the law of large numbers, we know that it converges to E (r(x,t)/n*+ -F (W, s)]
We assert that if in holds, then for any E t)(X + (V, - (W) # 0. Proof by contradiction.
r(X, t)/n = (W) - (W, x) for any t. we have
E{ - H (W, x) = (W. H (W,+) for any t.
This implies t(p(Y,t) (W,)))) = 0 for any t, which is not consistent with Hy. Thus, we have (r(Xxt)/n" + (W,t) [) * O. which implies E [(r(X.t) - (A A-rX, (X),) # o. Combining these results, we have
lim Vn --00 na ( 0,(1) + +
//
x2(11) 00.
lim max [Ta(s,))== 0p(1)) 1-000 BET
no
for almost all under Hy, (0 a < 1/2), where (+) follows from the Gaussian process.
(iii). The case of H with 1/2.
Following the proof of Hÿ, with 0 < a < 1/2, we have
(s, = vnP,, (U(W,Y, [m(X. A)-AAA)-'9
(vn) + (1)
+ VnP,, [r (X, t) / + - top na
Taking a 1122, we obtain
(s, = vnl,, (U(WYY,1) (m(x))) (A(A) A)-)))aa))
+ VnP, n [r (X,t) /n (W.t) (W, s)] + op()
For the first term, according to Theorem 4.12, we have VnP, (U(W,Y,1) -(AA(A ()-(A-))- converges weakly to G(s, t). For the second term, we have
A VnPn, (WW.))
H
(W,t)
(W,t))m(XX s)]
= VPP [ (1((X.0)/22²++*PPP (W,,)) A" A'A)- A (W)) +r(1)))((((0)1)))))) (
- (A4-)))-A4))-14>(xx1))
40

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Thus, have
1 (n(s, Vn
(U(ws Mes (I () + m(x),))
() (AtA))))(((AA)""=xxxxxxx)))) op)())
Besides, we have:
E ({U(w,, yo. t)vnt) [m(2;,
(x:)))
= (w, yi. (m((, (A(A' A) gx 9.. +n (mmm "rtta:/""J|mmc....aa///"a""""""" A)aaaa(())))
2n 1/2E (E(U (ws, v4, ttl..)/l21,.......) t))(z ()(A-))- 9x,")31)
(w Ue (m 8) 99 n =E (r(z,, (m(()))) ()
(I)
where (1) follows from EEU wr. 2) = r(,,, f)/ n under Hf. Following 4.12, we have (I) < co. Besides, for
the second term, we have (II) I 2E (r(2,,) +)² 2E (Im(z;,) 5) (A(AAA)) 9. (x:))) < co by inequality a²b² <
(a + /22 and assumption 4.11. Thus. we have EE(U (w,, Yr, t) rrr,,t/////i [m s) (A((AAA) (xr,)) + (Ws Us. H) 8) - {A(A* (((AAA)-1 (ax)))))) as n - 20. Therefore, the first term of (s, t) converges weakly to G(s, 1) in (T T x μu} by Lemma D.15, where G(s, () is defined Theorem 4.12.
Besides, we need to prove [fr(X,t) f) A(A- ")''r(X,,)))m(( s]] co. In fact, we only proof [r(X, t)m(X, s)] Do and E [[A* A)-"A"((XXXX m(X, s)] < oc.
For Bl]r(X, s)]}. we have
(),(XX < 1r(X,011 (M)X
S |r(X, |m 8) s) 1/2 (F())))
where the last inequality follows from (21) by assumption C.3. Since ||((X, 00 and the selected weight function m satisfies [| s) () < co, we have E
<0..
For B
m m(X, +]] we have
B ()-A"))))A(((()) m(X
< C
SC ([(" A)""
where the second inequality follows from assumption 4.9 and the last inequality follows from (21) by assumption C.3. Since (A* ( t) E Hw, we have E [((AAAA)=A''((X t) ())))m((555) s)] < 30.
Therefore, the second term of T, (s, converges weakly to μ(X, s,t) = E [(r(x,t) - A'rXX" By Slutsky's theorem, we have [s, t) converges weakly to G(s, t) + s,t in (T x T, T. μ} under Hig with a a=1/2.
Corollary 4.16. Suppose assumptions in Theorem 4.12 hold. If (yy t) is continuous with respect t for each y, then Apm weakly convergent to maxteT Jr [G(s, (){dd((s) under Iflo, as n, K - cc. Besides, conditional on the original sample {y, wi, x. the bootstrapped statistics (8) is also weakly convergent the maxEET J [G(s, ())))ss)))
Proof. (i). Apm is weakly convergent to maxtet dd)()))))((s)
We denote (t) f (s, t)} du(s) and (t) weakly converges X(t) IT du((s) integral of the Gaussian process G(s,t) still Gaussian process with respect we can obtain the variance f /Var (n(X, W. Y, t)) dy(s). Besides, for the Gaussian process, X(t) is continuous in probability if and only if
dp(s). Since

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
its mean and variance continuous following (Seeger, 2004). Since sly, i) is continuous with respect 4, the vari- ance is continuous. Therefore, (t) is continuous in probability. Assume that we obtains the maximum value at to, Le. max.eT X (t) (to). Since the process (t) is continuous in probability, we have that, for any e > 0, there exists 6 such that as long as jt - tol < 6, (t) - ((fo) > s/3) <.
Since (f1,-.,4k) are evaluated at a grid of equidistant indices, for any to E T, we have limx- ming 1/0 - 4 O. That means, for any > 0, there exists Ko. such that as long as K > Ko. there exists la with 1 < S K, Itx - tol < 8. Further, for any finite ti..... tk. denote TK == {k: ((xx) maxjşK Xt,)) and set óg : X (tka) X (tt:), where tot E Tk and X ((t)) == arg & X (tz). For such K, there exists Nk. such that when n NK. P (fxxn((t) x(tx)] minfe/3, 60/2) < 2K for any le I K. Therefore, for any 0, there exists K > Ko such that minaga 14x - tol < 8, and Nx such that for any > NR. we have:
P( max (txtt) x(to)] E) S P(1
max X (tx) X, (tmo)) > E/3)
KSK
KSK
X(1o0) e/3) + +((X X(1o)) 8/3) s c + P(1 max (tx( - (n(te e/3) + +P((X(tko) X(to) > E/3).
P( ((
KSK
For P([ (( X.((k) - (z > e/33, we have:
P(| max (t&t) - (1oo)))) > E/3) $ P(mmx N,, (1x) & (tto)) KSK
KSK
S 133, & TK max Xn (tx) = (t)(t)))
KSK
Pmaax Xn (tx) (t,)
KSK
I
E
(() X(t) + x(t)) X(tro) + ((kmn) - (n(on)) >0}
JSk
P(Xn((t)) - X(t) + - (tx(tk)) > 2803
JSK
- x((t)) > 50/2) + (t() X(tra)1 > 80/211 Σ (
iSK
Σ
E
+
2k) )c.
=c.
2K
jgk
Denote k == arg minggK Itx - tol- Then for (t(o)) X(to)l > E/3), we have:
(())) - X (to)) > e/3) = P(X((o) - -X(ta.) > E/3]
= P(X((n) - -X(x) + x(1x - -((kk) > e/3}
S (((((1))) - -((4) > 8/3] Se.
Combining these results together, we have lim limx-1o maxeçk X, (1x) =A maxtET X(I).
(ii). Bootstrapped statistics (8) is also weakly convergent to the max:eT f [G(s, dd)
By Theorem 2.9.2 of (Wellner et al., 2013), (s, = t Ei w/Ū Vr t)m( s) is weakly convergent to G(s,t) t) conditional the original sample. Applying the continuous mapping theorem. f (s, du)sss is weakly convergent to [G(s, du))). Using the proof in (i) again, we can obtain that A max fr (7,,,, (s, 4)) is weakly convergent to maxter f (G, ddpss), conditional the original sample.
D.5. Auxiliary lemmas
We provide some auxiliary lemmas needed for the theoretical analysis above. Lemma 1010 (Lemma 2.5 of (Beyhum et al., 2024)). Let (W, ll Ilw) and (X, I1.l be two Hilbert spaces and A W + X be a linear compact operator with singular value decomposition given by (Ans Pn, 22) |l llap be operator norm. Let T: W W be the identity operator. For each X > O, we have the following results:
42

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
(a)
llop || A-A) Al S 1.
(b)
[IA(A + A-A))=1 2.
(c)
1
Ia + AA lopp = IA(N + A'A) 'lopp S 2v
(d) For any y 0 and W such that I19112 = Σ3 M (9 12)12 < co, there holds:
IMAA + A) slll = 0 min(o.2)
Lemma D.11. If (w, 1) the least norm solution to the linear inverse problem and satisfies assumption D.2 (b), then the solution to the Tikhonov regularization Ha [w, t) satisfies that:
|H t) (w, D. S 5Op( (.
Proof. For the operator A Hw + Hx defined in (24), its singular value decomposition given by (An, Pr: ) Thus, we have Hi = Σ3 (HC (2)) Besides, according to Asp. Andn and A on AnsPn we have H;
x²
(A*A+N x (7): Thus, we have
|H>w; (w. 1) Hº (w,, 1
Σ
-1
(Hi
Hw Pj
1
+
-Σ
X
13+1 -1
(H) .....ssuw
22130
1(H Hir 12
Σ (23 + X)2
226
A1
K
spy) 12
< sup ;
Σ
x2 +
26
Applying assumption D.2 (b) for 0 2 and the maximum singular value of the operator equals |Allop < oo, we have
2
2
) S -suupp?-1 0(22).
susup
sup
X3 + 2
X3
+
For 0 < A < 2, we define = X3 and f(x)
Noted that f(z) is maximized (by using the first order condition) at (z+
x²
x = A0(2 - -) Thus, the maximum value of f(x) is
x
& 1 00 - -2 0
S O02)).
(a Aj²
4
The proof is complete.
Lemma D.12 (Lemma 12 of (Mastouri et al., 2021)). Suppose assumptions C.2 and C.3 hold for constants and K, respectively, Define of and o² as follows:
of = E(Il((),
0 = [[]lxx())))))l (W)
43

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
For A, J defined in Eq.(2), and A" in (28), the estimates J given by (26) satisfy the following properties with probability at least 6:
2 log(2/0)
20 log(2/0)
16. bel _
= 0pp
2x 1og(2/8)
20 3 (2/(2/))
||. - Allop
Op
Vni
2x (2/(2/6)
20 2 log(2/0)
| - Aloppp <
Op
Lemma D.13 (Lemma 13 of (Mastouri et al., 2021)) Suppose assumptions C.2 and C.3 hold. For A, A defined respectively (24) and (28), the estimates given by (26) satisfiex:
A A Allop = Op (Tm
Lemma D.14 (Lemma 2.4 of (Beyhum et al., 2024)). For random variables X, W. let m(.) be the function such that [[ bounded. Besides, we denote F as a class of functions of W such that FN -llll < cc. If [( fo)mmllc3F(a.nnn c24p((,00) (1( and Pf F) + 1, then
Vn((P. PP - -0)mm 0y(1).
Lemma D.15 (Lemma 2.1 (Li et al., 2003)). Let ....).... ,,,() be independent and identically distributed zero mean random elements on (S, v) such that B(IlZ(()](<(ss.,,)) : (J, (s)vv((s)) Here, (² v) square integrable
function space with respect to the measure Then n n- 11/2 Σ Z,(-) converges weakly to a mean Gaussian process
with the covariance function given by (2(,,'' = (z) 2z (s')).
Lemma D.16. For operators and A and their adjoint A and A", we have the following transformation:
(Al A Ā A-A) (N + A' (A Ā)() -1
Proof.
(AI+ + A'A)
I- Ā (W 1
=(M (M + - (M (N A'A+ + +
=(XI+ A'A)) - ++AA))((+++
=(AI+ A (A"A A Ā) A'))
Lemma D.17. Suppose that assumptions C.2, C.3, and D.2 hold. The PMCR estimator Aw, 1) satisfies
(w,) - (w, f)ll = Op { { + nd +
minfo.2)
Vnd
In particular, if assumption 4.10 holds, we have [W - H (w, t) t) = op(1)
Proof. We first decompose the estimation bias into two parts:
(w, (w, t) - Hw,, t)ll I (w, (w, t) + (H> -H (w, b)llo
We first consider || | t) - (w, m) In fact, following the decomposition (37), we have Aww t) Hw, 1) = G1 + G2 + G3 + Ga,
44

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where C1, G2 G3, C2 are defined in (38)-(41). For G1, can apply Lemma D.10 (c) have II A A)
Op(1/VA). Besides, according to (43), we have |66 AHH = =Op11//n). Combining these together, we get
Mw & |I(AI + A" ( "lolop 1|b ĀHς lmw = =Op
Vnd
= Op(1/A). Besides, according to (43) and Lemma D.12, we have l|b, AHl || Op(1/vn) and "| Aollop (1(1 (l/vn). Combining these inequalities together, we have:
For G2, we apply Lemma D.10 (b) to obtain that II +
(11 ||22||Hw $ $1(A + A A) 111 A'llop 10 AHęllHw =Op 1nd
For G3. we have:
Gslllw I + -M llo 18: - -AH, AH
=I(A + -
(M + A" (1))=((2) 'A*- - -1b ĀH,
1IAT++ + /' - A a' 11 1
+ +I
+ A" - Ilop -| - Aolp 110 AHf ll Hw
Since À and are compact operators, we can apply Lemma D.10 (b), (c) to obtain that (I A) (N + A A)A)-AA" llop (1/A) and [I(XT + A llop (1/A))). Besides, according to (43) and Lemma D.12, have |bg- AHA New = O,(1/ and | - Aoppop = (1/1/V1)) Combining all the incqualities, we get
1 + +0p
|Gallæw = Opp Vnd
For G4, we have:
|Gal|
|I + X' A AH, lmw
|A(AI + -A:A-
-1
=IIA(AI +
+ AA H2 ||w
$ ||AATT + -'A)
+IA(M+++')) -AAA)) A" A-A)-1 Hw
< |[ + Ā - -'''lopp 1|1-1AA -1A-~llopp |I(A + AA'A) A''))_llop || || He
+ ||A(AT + -/ A) 1lop 1 1opp ||A(AT + A A)- 'lap HH low
where (1) follows from:
A' A' = =[(AT+ (- (AT Āー - (AT A" A(N+ A²A) An]] [ H
(AI +
+ (21 A }) H1
=X(++ - A (AAAAA''A)) + A'A)-
Since À and A are compact operators, we can apply Lemma D.10 (b). (c) obtain that || (10
|[(AT + A* loop = (I/AA)), ||{ + A A' S 2, ||(AT + Á) S 2. Besides, according to Lemma D.12. we have || loop || - Allop Op(1/v/) Combining all the inequalities, we get:
|[GallHw = Op
Vnd
45

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Combining these results for 11 to G4, have
( (110 + nd )
t) (w, (w Op || (w, (w,t)
Vnd
Next, we consider [WH((w,)) 1) (w, 1) Hzw By assumption D.2, we can employ Lemma D.11 to obtain that:
|| |1 (uc, 1) - (w,, f)lleewwwwwww--mmmm))))))) = Op (> min:0.2)
Thus, we have
[w, f) (w. f) = =0p { Vnd
min(o.2)
+ nd
+ x
By assumption 4.10, have nd and A which gives [,/ t) 11 1) oI()..
46

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
E. Existence of solutions with two proxies
E.1. Proof Theorem 6.3
Theorem 6.3. Suppose assumptions 6.1, 6.2 hold. For any h(w, y) that satisfies (1), Ho holds if and only if h(w, y) also satisfies the integral equation (10) for any fixed
Proof. Suppose h(w,)) y) satisfies p(yl=)= = h(w, p)p(w/)))d.. Under Ho, we have X T (W, Y)J,, which leads
/ plll = p(ul=)) = = h(w, v)p(wola)du
= - V h(w, s)p(w)))wwee) (dw plul (l)x)du..
By the completeness in assumption 6.1, h(w, y) solves the following integral equation for all (u, y)-
p(ylu)= = h(w, y)p(wle (op(wpp wlw)ddw.
Since Hy holds, we have Y (Z, X)).. Therefore, for any fixed x, taking expectation over p(u[z, on both sides, we have:
p(ulz,æ) f p(olu)p(ull,,z)du
H h(w, v)p(wl))dw p(ul2, @ f now, 2)
where "(1)" is due to W I (Z, U. That means, h(w, y) solves the integral equation (10). To prove the contrary, note that if h(w, y) (I) satisfies (10), by W 1 (Z,X)JU and 1 ZI(U, X), we have
/ plule, r)p(ulz, x)du = p())2,x)
= = New, v)p(w z, x)dw
H h(w, (ww.v)p(wlwdddd 2()),,,))d. x)du.
By the completeness condition in assumption 6.2, we obtain
x) = / heww,
Since the right side of the equation is independent of æ, we get p(ylu, x) = p(ul u), and thus HH holds.
E.2. Discussions of causal inference and causal discovery
In this section, explore the distinction between causal discovery and causal inference, focusing on why the causal relation cannot be identified solely through the causal effect. We begin by presenting coumter-example that demonstrates that even when the intervention distribution for each x is identical, the independence X JL Y!U may still fail to hold. Following this, we provide an in-depth discussion the differences between causal inference and causal discovery.
We first introduce the notations. For any discrete variables X,Y, with k categories, we denote P(X) == (p(z1),..,,p(en)) P(Y = {p(wlays and IPYY = X| Z) = (p(yl xi, 2:.2
Example E. Suppose X,X,Y Y are binary, and the causal diagram over (U,XYY) is + X,U +YX + Y. The LXX,U U) are given by:
conditional probability matrices P(U), P(XUU), P(Y|X,
(0.4) 0.6
00.2 0.4
(0.5 0.1)
P(U)
P(X
P(Y = 0|X,U)= 0.2 0.3
0.8 0.6
By the definition, we know X 1 Y However, the intervention distribution is the same, iee. Poy|do(X = 0)) = P(y]do(X = 1)) for any y.
47

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Proof. According the backdoor formula, we have
P(Y y yjoo[X = 2)}
Σ
P(Y = u|U u, u, = x)P((U = u).
NE{O,1)
Plugging P(Y = 0|x, U) into the formula, we have:
P[Y Oldo(X 0)} = 0.5 x 0.4 0.1 x 0.6 0.26
P(Y = O|do(X = 1)} = 0.2 x 0.3 x 0.6 0.26
P(Y = 1|do(X = 0)} = 0.5 x 0.4 0.9 x 0.6 = 0.74
P(Y = 1|do(X = 1)} = 0.8 x 0.4 + 0.7 x 0.6 0.74,
which implies intervention distributions are equal. However, through data generation, we know X A Y!U.
Next, we will verify that in this example, P(Y = y| = 2) + E,, P(Y = uIU P) = u|X = z). which implies the example contradicts our assumption that there is no solution in (1) under H2. To this end, we need to obtain probability matrix P(Y]X), and (). First, by P(U) and P(X we can get the probability matrix PPX) and F(U
0.2 0.4)
(0.6) =
(0.4)
(0.32)
(0.75
(0.25 8/17
P(X) = ((X)))P(U)) = 0.8 0.6)
0.68 .(UIX)) =
0.6
0.75 9/17
Besides, we calculate the probability of F()|z) for any y, x. According to the Bayesian formula, we have
P(Y = y = = Σ P(Y = y|x = x,- = = u|X = x)
u
P(X - 2/U
uP
P(Y = y|X = x, = u) Σ
PPX x)
Therefore, we have
P(YIX) = 0.2 43/170
0.8 127/170
According to the Bayesian formula, we have
F(Y u|U =u) = Σ P(Y s/X u,U u) x/U== == =u).
Therefore, have
0.26
0.22
P(YU))
0.74
0.78
Thus, we can verify
P(Y 0|X = 0) 0.2 A 0.23 = 0.26 x 0.25 + 0.22 x 0.75 =
PPY ΣP(Y = OJU uPP u/ = 0)
P(Y = 0|X 1)
43
203
0.26
+.222
X
9
Σ P(Y UFU = uP = u|X =)
170
850
17
17
P(Y 1|X 0) 0.8
+ 0.64
0.22
0.25
0.78 x 0.75
Σ TP(Y IIU u)T(U u/ X 0)
P(Y = 1|X =
127
439
0.22 X
17 +
0.78
x
9
Σ
P(Y = I|U = uPU(U = u|X = 1).
170
850
17
More discussions about casual discovery and causal inference. Causal inference and causal discovery address funda- mentally different problems (Guo et al., 2020). Causal inference focuses on quantifying the effects of interventions, often
48

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
requiring strong assumptions and additional information to ensure accurate estimation. In contrast, causal discovery aims uncover the underlying causal structure, emphasizing the identification of causal relationships rather than their magnitudes.
It may not be feasible to infer the causal relationship from the causal effect. One key reason is that the inference is often complicated by noise in the estimates, making it hard to determine whether nonzero effect arises from an actual causal relationship or random noise perturbing the estimation. Even if we can estimate a confidence interval for the effect at each treatment value (Robins, 1988; Robins et al., 2003; Calonico et al., 2018; Colangelo Lee, 2020), there are no valid statisties determine whether the relation exists. Moreover, as shown in the previous example, causal effect of zero does necessarily imply the existence of the causal relation. Additionally, estimating causal effects often requires satisfying other conditions. For example, proximal causal inference depends additional completeness assumptions (Miao et al., 2018; Mastouri et al., 2021). In our scenario, such conditions are assumed on Z X, W (i.e., for any square-integrable function (x) pr w almost surely if and only if g(2) = 0 almost surely) and (X, WHXX, 2} (Mastouri et al., 2021).
E.3. Proof of Proposition 5.1 and example 5.3
We first prove Proposition 5.1.
Proposition E.2. We consider the linear Gaussian generation mechanism:
U ~ 1))
X ao + ayU + 10),
W = 3o + ByU +00,)
Y = 7o + UUU + 7 + W + +0(00,) 1).
When ?/ 0, long as hxx|
B B
where A 14 a
B
2ne ax
2
20y and
2A
Pg
au
86:
4 (1+a]++08) (1+ +-6)
the integration equation (1) has no solution. Further, if lnw]
|C|+|8 |yxl+ A where
2
C' = 1 - - /..BB and D aoPo
Bu.
7x
96'
(1) has solution.
Proof. Based on the data generation structure, we can obtain joint distribution
0
au + at av
Bo Bu 1 +B Cov(W,Y)
Cov(U, Y) Cov(X, Y) Cov(W, Var(Y)
X
N
ao
au Be Cov(U,
W
Bo Tw 90
Y
70
+ 7x
Cov(X,
where covariance Cov(U, Cov(),,ooXXn Co, Co, and Var(Y) are respectively
Cov(U, Y) = 7U + 7xau + BwBr
Cov(X, Y) au (To + Su 7xar) + 7x
Cov(W,,) = -Bur (rv + auTx + 3w + 7w
Var(Y) = (Tw + 7xar 7w Bu) + + 7/ +
We can therefore derive the explicit form of the conditional distributions ((wl(r) and p(ylzr))
kW + Cov(W,
Cov2 (x - -x)), Var(W) (1 Var(X) Var(W)
W = - ~ N
Var(X)
N {
WIX ux
WIX 2
+ Mo
T MY + Cov(Y, X)
(x px). Var(Y) (1- Co,2(Y,X)
Var(X)
Var(X). Var(Y)
N NIXII VIX
VIX
1x x 10 ,
49

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where Px
WIx W
:P Twwxx
x
XIX
x|x :Po
and xix are defined as follows
WIX
on
Fx
140%
Mo x 30
noau
8u
1+0
O wix = =1+ g8
Bu
1+a0
XIX au au +o Bu +7x 7x
Htae
10 XIX 70 0000:000
Bu +2x ( +4007x
1+
O = (YU + TxaU + Bu)² 2 + 73 + 7iv +
(ove Tw By +2x +xx)
14a0
By applying Lemma B.4, the solution of (1) is given by:
140 - -IXIX 140 /1X wIx)
VIX
vIX WIX
) -IXXX
- /NX /
WIX
y
w
h(w, v)
VorIx - (uxIx) / (1x W ix
2
V (1x) 2
/ (1
where o is the standard normal distribution's 's probability density function (pdf).
(1XIX)
/ (ux
(IX)
For h(w, y) be meaningful, need x
0, which implies
1 20
274 ay
2700
20070
1
at
g
+
ay ayst
+
87
9x
1+
at
2
+
B
apist
92
PU
1x
(54)
au
BU
V
T > 0.
Bu
OU
3u
We discuss the following two cases: (i) + (xx # 0) and W A Y (nw 0); (ii) X + Y (rx #0) (Tw = 0).
and and W A Y
1. 7x # 7w = O.
first rewrite (54) as:
2Yu
Y0
2
+
By
7x
1
0
22
0
2 > 0.
Noting that this is a quadratic function, can get its discriminant
A := -4-CAC-- 4 4(1 + at + (8) (1 + at + +28)
>
az
0. Therefore, this is quadratic function whose discriminant is always positive and opens downward. When 7x satisfies
Besides, we can find 1 +
at
ag
35
-BVA 24
7x
-B-VA
24
(1) will have solution. When
7x -B-VA or -B+VA (1) will have no solution.
Without loss of generality, we consider the case where ay and tu have the same sign. First, we can find B - 277 +
ay
2v1r
20arree R
< 0 since B3 > 0. Thus, we have B v
B + vAJ. Thus, when |7x >
B+vAA
(1) will
24
have no solution. If ou and Yu have the different sign, we have B = -( 2TT +
2 ag
200 A
) > 0 since 88 > 0. Thus,
we have 1 B- V
B A Thus, when |nx|
-B-VA (1) will have no solution. 24
Combining the two cases, we can obtain that as long as |xx
BB+VA the integration equation (1) has no solution.
2A
50

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
2. 7x + 0,7w # O.
We consider the case hx/> >
= BB+A
under YU > 0, since (1) have no solution. We can rewrite (54) as
2A
2
7x_
+
au
30
7x
+
oU
9u
xx +
YW <
1- 2ú
/8
2
og
27u B
+
mox
8
Xx
1
+
2
1
2
00
a
2..
Thus,
if
_x
+
au
Bu
7x
+
ov
Bu
1x
+
26 Av
O, we can obtain that (1) may still have a solution, as long as
1
200
2710 02
20 89
7X
1
+
Dg 7
av
T
>
2
7x
7x
+
Bu ox
7x
+
30 By
We find that if long law] is sufficiently large, the solution to (1) will still exist when hxx
7x
ou 7x + 9u
mg 7x +
TV
< O, the right-hand side of the above inequality is positive. That means, as
Ae
Ag
-B+A
2X B0
7x
+ au
7x +
7 Dy
>
O, we can obtain (1) may still have a solution, as long as
74. 1
2>2 au
270
2ag 90
1
oxpy P7
7x
ai
apPE
Ph.
TW
2
0030
+
Do
1x
+
By ox
7x
+
D6
We find that if 7x + ou 90 Bo
xx + Be ag
Xx +
O, the right-hand side is negative. That also means, as long as |wwl is
sufficiently large, the solution to (I) will still exist when |7x|
BBVA 24
_7x If YU < O, the proof is similar. Besides, in the above cases, long as hw/ |C + with D ou Bo 7x + AC au 7x +
+
2
oo Bu
20 36 (I) has solution.
Remark E.3. If 7x TW- 0, (54) will become 1 - 12 > O. This means that if the strength between W- - U is greater than the confounder strength between W U, (1) will have a solution under HD. Otherwise, similar to the case when 7x A 0, if the effect of W on Y is strong enough (i.e., (ww)) D), the solution exists again. Specifically, if rx = 0, nw o O, 2/ 2 (54) will become 3 Py 7 > 0. If 2 3u Be ?w is large enough, (1) still have solution. If we /Bu 0, we need
TW to be as negative possible; /t// /0u < O, we need yw to be as positive possible.
Next, we prove the claims in example 5.3. We show that as long as the coefficient of W + Y is strong enough in example 5.3, the solution of the integral equation p(y/zr) = f h(w. v)p(wlx)da exists. As an explanation, we will show that a key condition in Picard's theorem B.2 holds, namely, the series Σ > (f(ulz), o,)) converges.
To compute the series, we need the singular value decomposition of the operator T L²(F(w)) + [F(((x),, where Th = Efh(W, y}|=} = p(y(z) for all (z, y). Based on the data duta-generating process in example 5.3, both (3(F(w)) and (F((()) are square-integrable spaces with respect to the standard Gaussian measure. For such spaces, (Carrasco et al., 2007) derived the form of the eigenvectors on. as stated in Lemma E.4. As this result builds upon the concept of generalized Hermite polynomials, we include a brief introduction to facilitate understanding.
We first introduce the concept of Hermite polynomial, which is defined in the square-integrable function space with respect to the standard Gaussian measure. Specifically, we say that function f : RR- + R is square integrable w.r.t. the standard Gaussian measure 7 e -²² 222 if E.~~V(o.) {f²())} < co. We denote by L (2R) the space of all such functions, whose basis functions are characterized by probabilist's Hermite polynomials
He, (x) = (-1) x22 dx dak e
=²/2
The first three Hermite polynomials are Hoo(x) = 1, He (x) = x, Hez(x) = -²- - Let hoy(a)) == Hen[(z) normalized Hermite polynomials, which form a complete orthonormal basis in C²(R). Thus, the Hermite expansion of a
vil denote the
51

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
function / (IR) is given by
f(z) = Σ (x-1(f)hex-1(t)) pk-1(f) = Ex-xva.n)f(((xbe=-1111
Indeed, Hermite polynomials can be equivalently defined by identifying eat- generalized Hermite polynomials Hn (z, y) as those that satisfy
=Ek=
(e,(w) pb 13
We can define the
Hn (z,
k
In the case when y = 2. we have (x -1 -1/2) (x). Generally, the generalized Hermite polynomials Hn (z. (xr.)) related to the Hermite canonical form He, (I) as the following way:
Ha (2g) Hen
(55)
iv/2
We are now ready introduce the eigenvalue system of the operator T' [W) (X)) derived by (Carrasco et al., 2007).
(W - (X), [[J(W)XX - where C2 is square integrable space with respect to the standard Gaussian measure, i.e., (W,X X) is jointly Gaussian with zero mean, unit variance, and correlation pw We have is a self-adjoint operator, and the eigenvalue system for T is given by (P(w)) = he, (w), (x) = he,(s), >, = pw where pwx the correlation coefficient between W and and he, is the normalized Hermite polynomials.
Lemma E.4 (CCarrasco et al., 2007)). Let T
Now prove the result in example 5.3.
Example E.S. Suppose that X, Y. W satisfy the linear Gaussian model, i.e. U- EU, X 2 + EX, W = 20 +
Ew, = X + U + Tw W + where EU, Sy, ~ (N, 1). The integral equation (1) has a solution if and only -15+36vG the series >x= \, (p(ylz)) 2 converges.
-15+36v3 xe 0.61. Besides, when "/ 72+16v5
TW
72+16v5
Proof. We first show that even under Iil,, the integral equation p(y]a) = f hww, coefficient "w is large enough. Specifically, since X and W are normalized, based on the data generation structure, we have
has a solution when the
nw + 1+ y
U
0 0000
v5
X
2
1
7 + + v5 ww - 22 + v5).
N
V
1
v5 TW
+
+1+
v5
Tw
+
+1+
v5
YW
v3
7W 2
(2 + v5 +3+ vb
We can therefore derive the explicit form of the conditional distributions p(w(ar and p(ul(r):
Cov(W,X)
(xc - -xx), Var(W) (1- Var(X) Var(W)
Cov2
WIX - N
Var(X)
4 T;
9
N
25
ur + Cov(Y, X)
(x - ux), Var(Y) (1-
Cov2(V,
Y|X = x x N
Var(X)
Var(X). Var(Y)
((-5
4 ww
2
+ +
I, 2
?W +
(56)
V5
25
5v5"
By applying Lemma B.4, the solution of (1) is given by:
Yw + 2y5-5
1
h( (ww y)
(57)
94
10
ww +
9
20
9
10
YW
+
16
9
20
52

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
For h(w, 3) to be meaningful, we need 9+2v/ 16
10 T +
3
9// > O, which implies "/ > 16v5
20
315+36 R 0.61.
Next, we need to verify the conditions for the series in Picard's theorem B.2, which requires proving that En=o KS on < too for the singular system (An: Pn, 2n) associated with the compact operator Th- = S- In our data generation process, operator I' (W, 7) + (X, 7) satisfies Th = (B((WW, y)[z] = p(y|z) for all (z, y) and is characterized by the integral kernel (12). Thus, by Lemma E.4, we have : C2 + (X,, is self-adjoint he, (z):
operator and the eigenvalue system of operator T' is given by (w) ( is the correlation coefficient between W and X and he, is the normalized Hermite polynomials. Thus, we show that as long 15-+36/5 the following series converges, which can explain why the solution may exist under Hl:
pwx: where pw
as 7/ >
72+16/3
(p(y)zr), hea(x)))
2n
PWX
Let
In (p(ylz), (x))
1 - / 22
p(ylr)ho,
x² 12
1
/
Hc,l =2 dxdx
2m
We will use the generating function method to derive the analytic form of integrals. By equation (6) in (Babusci et al., 2012), we have:
/
Hn (ar + m,e
cx. -x2 dx
v Ve
exp
aa
Hn
b+
aa
m +
a²
4c
2c
4c
Let == nw + 1
2 and 0² == v5
23 Tmu
9
5V5 w
4
+ 8.
By (56), we can obtain
-
2 2
In 227n!
S 20
(n(()e-==
12 dx (1) e
Hu (I, - 11) dx,
V2na²
2n
where (1) follows from Hn( 11/2) = He, (x) by (55). By taking a = 1,b = 0, have:
o²ty²
20
oz
and m -1/2, we
In
20 4n Hn
nyy
v 22mn! (0 + 122)
0 + 112
2( 2 + 11
We consider the case when μ # 0 and μ = 0.
1. When / + 0.
By (55) we have
2
wv +
In
C 200 112)
1
Hen
V 20 + 1
02 + 112
i
12
1
(-1)²
122
+ Hen
!ulv02 02 + 12 }
2z + 122)
2 +
Therefore, the series can be written as
2
1)"
22
2
1 (Tullo oz + we )² (58)
Hen Vn!
Pix
Dixv20 2T
(0² 2 +
12)
0² + n
R=
R=0
According to equation (18.15.27) in (Olver, 2010), it is known that for fixed M- = 0, 1, 2....
M Um (1) cos Wn.m (I)
Hn (x)
0
1
(2n + 1)²m m
(2n + 1)
53

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
where (x) (-1)" de e- is physicist's Hermite polynomials, the coefficients Um (2) are polynomials in x, uo(x) - (,,1(() = 1/6x", Wnm (x) = (2n + 1)2x= { + n) and
T(n + 1)/T +1 if is even,
An
T(n + 2)/ (2m + I) (2 })}
if odd.
By taking M 1 and He, (x:) = = Hn
we have
V2)
(z) = = e cos mi,
+ V2
}1
2 v Ancd +2 +2 An xo
(2n 1) }
(2n + 1)
2 2 Aneda?
Since each term in the series Σ = 12 is positive, its odd-term series and even-term series are both positive term series. If both subseries converge, then 2n 1 also converges.
(i). Even term series. By Stirling's approximation n!
22n (+)
we have:
I(n + 1)
n
v 2m
(e)
2T
2
[( +1]
2
2e
Thus, we can obtain Ho,, S 2 +1 (7)
T for large values of n even. Then for even we have
(e (x)
23
c /2
n!
2nm
Thus, we have
In
122
23
pixx
2x(l²
w²) pivx(o' 22)
22nn
e
== Jn.
We only need to prove the convergence of the series En iveven Jn. By the ratio test, we can obtain:
e 2 n
Jnt2
2 TH2
Pi
(2
lim
lim n
x2+00
Jn
2
Pirx (0 =2)
12
2
72
lim n-000
Piv
2
V n (02 + 12))
This series absolutely converges if and only if
112 Pixx (02 + 112)
which holds and only if yw
-15+33655 72+16V5
2s 0.61 by taking μ = -w + 1 +
v5.
2
0²
26
9
5/5
2w
+
5
6
and
PWx -
(ii). Odd term series. By Stirling's approximation n! ~ / 2TT (2)
we have
[ + 2)
/ 2m 11
ni n
n+1
2m n n+1
2 )
-
(
1
2 2
n2 +
2n + r( +
2n
+1
2e
2n
+
e
Thus, obtain He (I) $
2 23
n+1
m+1
4 for large values of odd. Then for odd n, we have
2n+1
He, (x)
16(n + +1)
2
(n 1)
(n + 1)!
e(2n + 1) 2T
54

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Thus, have
(
In
0C
122 x +2)
2 16 (n (1++1)
20 2o4 ==
Jn.
PW
2 + 12)
2n +1
2 (n +1)
n=0
Again, we only need to show the convergence for the series
En is odd
Sn. By the ratio test, we can obtain:
n+2
11
16(x+3)
2(1
Jnt2
2m(
2n-5
(2a(++3)
lim n 2-000
lim
Jn
n 220
10(n+1)
1
1
25 Oirxlo+ux'
2n+1 2x(
4 Pirx (o² - 12)
(2n 1)² (n + 3)
lim 22000
(2n - 5)² 1I)
We can verify converges if and only if
(2n 10++111 (a+3)
4 16n +13m++113
for positive integers. Thus, this sub-series of odd terms absolutely
(m++))
In 21n +55++25
12
<1,
Pw x (o2 + +22)
2
which holds if and only if pwx G we have yw 25 Yw
-15+36v5
0.61, by taking μ - + 1 +
2 v5"
a²
72+16V3
9
5/5 w
6.
Combining two results, we know that when "w
-15+36v5 72+1645
0.61, the original series converges.
2. When / O.
Since 0, we have "w = 5+2v5 1
15+3665
Thus, the distribution of p(ylx) becomes
72+16v5
1 (29 + 455) 0,
-IX=
16
Thus, we define 22. = 1 (29 + 4V5), we have
16
2e Hn (x)e 22 (2)e 12
He (z)e -2² 222 dz.
I
2Tn!
2TTOy con
2T oo n!
Following Lemma 2.6 of (Davis, 2024), the integral of the stretched Hermite polynomial S, only non-zere for and has the value Sn (n 1)!!(?? 1)m/², We use the above results and take 1, we have I, = 0 for 2 1. Thus, the series is:
1
He (yx)e > /2 is
<2n
PWX
I 2
=
plrx
To
= (Io)
(
20
Eon (Je (fe-21adx
(2)
4.
2T
2TO con
you < 00.
n==
con
where (1) follows from Heo(x: 1 and (2) follows from f 1e= 2T.
E.4. Proof of asymptotic properties with two proxies
Assumption E.6. We assume Ex (m(X, Z, s)W and Ex (Im, Z, 3) are uniformly bounded for all s. co, where
Assumption E.7. For any s, E T, (WU Y.W)), co and [[[m (X, Z, s) (A(A* A) (x)) ""
(.) = Efm((X Z, s))w(W)) (W)](.).
Assumption E.8. Let (Λ₃: 19j, Φyyy be the singular value decomposition of the operator A described in section C. Then
assume: (a). For some 7) 2 2, Σ, 2, 22 (gs, (9)) w |2 < oo; (b) For some 0 > 2, Σ, X; 28 (HO (9)) Hw co.
55

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Theorem 6.6. Denole Te,(
U(W [m((Z (A(A"))))( (X,, a)), where Is(;2)
B(m( s) (W) and O Y,,,YYX), Suppose assumptions in Theorem 4.12 hold. If Asm. 6.1-6.2, and E.6-E.8 hold, under Ho we have, (i). (s.t) converges weakly to G(s,t) s.t. fI (G(s, d){dd(()) < cc. where G(s 1) mean-zero Gaussian process with covariance E((s, (') O' (W', Y',Y,, (') is an independent copy of O: (it). A converges weakly max J [G(s, du))(s
(0, (0', where
IET
Proof. We need to replace the weight function m(x, s) with (ZZ s) over Z. By (36), we have (s, = VnP, (WW, Y, (2, + (Expected risk difference) + (Empirical process).
By Propositions D.3, for fixed x, we can obtain the empirical process
Vn(P, - (w,t) -
= op(1).
By Propositions D.4, for fixed x, we have
Vup { (W,)) - FA(W.0))(22.;)))
1 vn
U(w, ua A) (x,)) + oz)
Therefore, combining all the inequalities, we have
(s, Σ Uww, Ve, t) [m[x,, 21, - (A(A"A) (2(,))) + o,()).
Next, we apply Lemma D.15 to {U(w, ti, t) [m(x, 21, - {A (z,, .)]} to obtain T (s, converges weakly to ((x,4) in (T x T, 1 x fJ,, where G(s,() is a Gaussian process with zero-mean and covariance:
= [((s, t), (1, Z',''',''"'", xx)')
To show ts is zero-mean, noted that
E [U(W,,, (m( (4(""))))))
E [U
E [U1 -x,
0,
where the last equation follows from (11) and (4).
Besides, by assumption E.7, we have (U,; yis [m)(, 8) - 9a (w,, yr, [m( 2,, s) - {A(A*A) A) xz), < 00 for any (z,s, Therefore, by continuous mapping theorem, have A/2), converges weakly max J (G(s,t)
For power analysis, we define the global alternative Hifi and Hi, (0 < a S 1/2) of (11), in terms of ({p(Y,t) H (WW4))Z, x for fixed v.
HIfx E(p(Y,t) H(W, t) 2, x) * for some t E T,
for any H. Hw. For the local alternative IHf there exists 11x E Hw, such that
r(Z,,,t) n
H =(p(Y.)) =
f) +
Vt,
where < a 21/1 1/2, and for any H, r(Z.at) cannot be written as E(H(.. H. (1, 1)|Z, xr} for some L. n
Theorem E.9. Suppose assumptions Theorem 6.6 hold. Besides, we assume E(r(Z, t)*} < for fixed x and any t. Then, we have:
56

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
Global alternative. limp- maxIET (Tx (8, 1) = 00 for almost all s under Hf*.
(ii) Local alternative (ax < 1/2). lim m [T/1] (s, t)l
for almost all under Hin.
(iii) Local alternative (a 1/2). (s, converges weakly to G () ((,, z. 8, t) (T x xT, (ux under Him where G(s,t) defined in Theorem 6.6 and p(Z, X, x. 8, E
8)]
Proof. The proof is similar that of theorem 4.14, with the weight function m(X,s) replaced with m (n1 x. 8).
57

---

Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond
F. Additional experiments
In this section, we evaluate our two-proxy procedure to a nonlinear setting, where W - and both W and Z begin available. Similar to Fig. 5, the two-proxy method outperforms the single-proxy approach by leveraging information provided from the additional proxy, i.e., NCE.
Data generation. We generate U following normal distribution with mean 0 and variance I, denoted by U NV 1). Similarly, we simulate W = (in + &w and % = sui((U) + ez. The treatment assignment mechanism follows the generation process: X = sin) + + Under the alternative hypothesis X T YJU, the outcome is generated from Y = X + sin(U) + 2W + Ey; while under the null hypothesis X T Y the outcome is generated from Y (U) + 2W2 &y. both hypotheses, the noise terms &w,,,,,,,, ey are independently drawn from a standard normal distribution. We repeat the process times, where at each time we generate 100 replications under Hy and [H]
Type-I error and power. The average results presented in Fig. 6. As observed, while single-proxy procedure effectively controls the type-I error, it exhibits low power in identifying causal relationships. In contrast, by incorporating additional information from the NCE, the power improves significantly, demonstrating its efficacy in learning causal connections.
Type Error Rate
Test Power Rate
Ours (single)
Ours (single]
Ours (two)
Ours (two)
0.8-
0.8
0.6
Power
0.6
Type
0.2
0.0
0.0
200
500
900
1200
1500
200
500
900
1200
1500
Sample Size
Sample Size
(a)
(b)
Figure 6. Type-I error rate (left) and power rate (right) of our procedure and baselines in the nonlinear setting with two proxies.
58


Reviews:
[{"review_id": "EoBK8DdCVe", "reviewer": "Reviewer_EyGH", "summary": "The paper proposes a nonparametric procedure for bivariate causal discovery for determining $X \\perp Y \\mid U$, where $U$ is an unmeasured confounder of $X$ and $Y$. It introduces the Proxy Maximum Characteristic Restriction (PMCR) method to solve an integral equation where a proxy variable is available to determine $X \\perp Y \\mid U$. The paper also introduces a second proxy (NCE) and derive an additional integral equation for identifiability. Theoretical results establish asymptotic validity and power, while experiments demonstrate improved type-I error control and power compared to baselines in synthetic settings.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the positive assessment and valuable suggestions about our paper. We address your questions below. We will correct typos and make the number of lemmas consistent in the updated version.\n\n**Q1.** How sensitive is PMCR to kernel and bandwidth choices?\n\n**A.** We would like to clarify that PMCR requires kernels to be bounded, continuous, and integrally strictly positive definite (Assumptions C.2 and C.3). Typical types of kernels satisfying these conditions include Gaussian RBF and Laplacian kernels. We adopt the Gaussian kernel because it is commonly used in the literature. Its bandwidth is typically initialized via the median distance heuristic, ensuring robust performance.\n\n**Q2.** What is the computational efficiency of the proposed approach? How well does it scale to larger datasets and/or high-dimensional variables?\n\n**A.** PMCR has a time complexity of $ O(pn^3) $, where $ n $ is the sample size and $ p $ is the dimension of $ W $. The complexity is linearly proportional to $ p $ because, for multivariate $ W $, the kernel can be constructed as the product of scalar kernels applied to each input dimension $ w $, i.e., $ k(w, w') = \\prod_{j=1}^{p} k(w[:, j], w[:, j]') $, where $ w[:, j] $ denotes the $ j $-th feature. The complexity with respect to $ n $ can be reduced using standard Cholesky techniques when computing the inverse. We record the running time of a single trial as $ n $ varies (when $ W $ is univariate). As shown, the actual running time increases much more slowly than $ O(n^3) $.\n\n| Sample                | 200  | 400  | 600  | 800  | 1000 | 1200 |\n|----------------------|------|------|------|------|------|------|\n| Computational time (s) | 0.39 | 0.97 | 2.64 | 5.62 | 9.65 | 14.84 |\n\n**Q3.** Are the non-identifiability results (Prop. 5.1) generalizable beyond linear models?\n\n**A.** We believe that nonlinear models also suffer from the non-identifiability issue when $W$ has a strong effect on $Y$, but the theoretical analysis is challenging as there is no closed-form solution for non-linear models. \n\nTo empirically verify this point, we consider a non-linear setting: $U \\sim \\mathcal{N}(0,1)$, $W=U +\\varepsilon_W$, $X=U +\\varepsilon_X$, and $Y=X^2+\\gamma_W W^2+U+\\varepsilon_Y$, where $\\varepsilon_W, \\varepsilon_X, \\varepsilon_Y$ are standard normal. We generate $n=800$ samples and record the average power over 20 times. We can see that the power decreases as $\\gamma_W$ increases, which suggests the difficulty of identifying the causal relation when $W$ strongly affects $Y$. \n\n| $\\gamma_W$ | 0    | 0.5  | 1    | 3    | 5    | 7    | 9    | 11   |\n|---------------|------|------|------|------|------|------|------|------|\n| Power        | 0.85 | 0.88 | 0.85 | 0.68 | 0.58 | 0.32 | 0.14 | 0.07 |", "claims_and_evidence": "Yes, claims are backed up by theoretical statements (e.g., Thm. 4.5) and experimental results.", "methods_and_evaluation_criteria": "Synthetic datasets generated under varying conditions (with/without direct proxy effects) are used to evaluate type-I error and power. The evaluation criteria are appropriate for assessing causal discovery methods in the bivariate setting.", "theoretical_claims": "Yes, I checked the proof of Theorem 4.5 in Appendix B.", "experimental_designs_or_analyses": "Yes, the ones in the main text.", "relation_to_broader_scientific_literature": "It extends prior work by addressing the sample efficiency limitations of discretization and by clarifying non-identifiability conditions when using a single proxy.", "essential_references_not_discussed": "None to the best of my knowledge.", "other_strengths_and_weaknesses": "**Strengths**\n- Comprehensive theoretical analysis with detailed asymptotic results.\n- Clear identification and resolution of non-identifiability issues via a second proxy.\n- Thorough experimental results (albeit all synthetic).\n\n**Weaknesses**\n- The paper is very dense and details are hard to grasp without reading parts of the 45+ page appendix!", "other_comments_or_suggestions": "1. In line 289 it should read \"bootstrapped statistic.\"\n2.  The numbering of lemmas in the main text and the appendix is inconsistent (e.g., Theorem 4.5 in the main text is Theorem B.4 in the appendix). Consider using a LaTeX package like `restatable` to maintain consistent numbering.", "questions_for_authors": "1. How sensitive is PMCR to kernel and bandwidth choices?\n2. What is the computational efficiency of the proposed approach? How well does it scale to larger datasets and/or high-dimensional variables?\n3. Are the non-identifiability results (Prop. 5.1) generalizable beyond linear models?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wHHwOWCgO1", "reviewer": "Reviewer_nuHH", "summary": "The paper concerns bivariate causal discovery in the presence of unobserved confounders with the assumption that certain \"proxy variables\" are observed. Existing literature translates the absence of a direct effect of treatment on outcome to the existence of a solution to a certain integral equation and proposes parametric approaches that involve discretization for continuous variables given that a proxy variable that affects the outcome is observed. This paper proposes a nonparametric approach to test the existence. While the integral equation is only a necessary condition (under certain assumptions) for the absence of a direct effect, the paper shows that the integral equation is not sufficient in the linear models case. By assuming that another proxy variable that affects the treatment, they show that a condition involving the existence of a solution for a modified integral equation is necessary and sufficient for the absence of direct effect.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your efforts and suggestions in reviewing our paper. We address your concerns below.\n\n**Q1.** About the type I error level in Liu's paper. \n\n**A.** Liu's type-I error control requires a Lipschitz smooth function, but sqrt does not satisfy this. To make it hold, we conduct experiments with the sigmoid function. As shown, Liu's method approximates the type-I error as $n$ increases. It is also seen that our method has better sample efficiency than theirs. \n\n|Sample|200|400|600|800|1000|1200|\n|-|-|-|-|-|-|-|\n|Our|0.06|0.05|0.05|0.06|0.05|0.05|\n|Liu|0.15|0.13|0.09|0.09|0.08|0.08|\n\n**Q2.** No evidence to support sample efficiency. \n\n**A.** The table above shows our method's superior sample efficiency in type-I error control. As claimed in lines 107-109, this is because Liu's method requires the bin number to go to infinity to make the discretization error vanish. \n\n**Q3.** The power does not approach to one in Fig. 3b. \n\n**A.** We further conduct experiments as $n$ increases. Combined with Fig. 3b, they show a clear trend of power approaching one. When $n=4000$, the power is 0.97. \n\n|Sample|2000|3000|4000|\n|-|-|-|-|\n|Power|0.89|0.94|0.97|\n\n**Q4.** About definitions and claims. \n\n**A.** **About $\\hat{H}^\\lambda(w,t)$, its relation to the paragraph before \"equivalent speaking\".** As mentioned in lines 201-211, $\\hat{H}^\\lambda(w,t)$ is the minimizer of (5), which is designed to solve the least-norm solution as introduced in lines 176-181. From lines 182-208, we explain how the risk (5) is derived. We first introduce the risk in lines 192-193, followed by its equivalent form in lines 197-198. Indeed, the risk (5) is its empirical version, augmented with a regularization term. Additionally, we provide a more detailed explanation in Appx. C.4-C.5 due to space limitations. We will offer further clarifications in the updated version.\n\n**About the motivation of Tikhonov regularization.** The regularization provides a way to solve ill-posedness problems and is commonly adopted in the literature of kernel regression. We have explained in more detail in Appx. C.5. \n\n**About the bridge function.**  The bridge function is the solution to the integral equation used in proximal causal inference, serving a similar role to $H$ in (4). It is a fundamental concept in proximal causal inference, and we have provided references in lines 167-168 due to space limits. \n\n**Q5.** About power analysis after successive approximation. \n\n**A.** Our power analysis focuses on alternatives in (4), which are defined based on characteristic restrictions. By utilizing the characteristic function, our approach demonstrates significantly better power compared to first-order moment methods in the literature (Fig. 4).\n\n**Q6.** Does (4) hold if (1) does not hold under the assumptions in Thm 4.5?\n\n**A.** We are sorry that we're unsure if we fully understand the question. Under assumptions in Thm. 4.5, both (1) and (4) should hold, as (1) is the conclusion of Thm. 4.5, and (4) naturally follows if (1) holds. We would appreciate your clarification if we've misunderstood.\n\n**Q7.** Does assump. 6.1 say that all variability in U is captured by X? That seems to make a \"confounder\" immaterial.\n\n**A.** Yes, it means variability in $U$ is captured by $X$. However, we want to clarify that this assumption is easy to hold and will not make $U$ immaterial. \n\nFirst, the assumption applies to a wide range of models, as long as the dimension of $X$ is greater than that of $U$. Under this condition, [Andrews et al. 2017] demonstrated that completeness generically holds. Even if this condition holds, the confounding problem still matters as it also causes bias since completeness does not imply that $p(x)$ would determine (up to transformation) $p(x,u)$ (this is known as the ``equivalence condition\" in Miao et. al. 2023) and hence the distribution of $p(u)$. In this regard, we still cannot determine whether $p(y|do(x)) \\neq p(y|x)$ when $U$ is unobserved. \n\n**Q8.** About \"there exist\" statement in $H^0(w,t)$\" (line 234). \n\n**A.** Local alternative refers to alternatives that are very close to the null hypothesis. That means, it adds a small deviation $r(X)/n^\\alpha (\\alpha > 0)$ to some $H^0$. To ensure that it does not degenerate to $\\mathbb{H}_0$, we require $r(X)/n^\\alpha$ cannot be written as $E(H-H^0|X)$ for any $H$, as claimed in lines 239-241.\n\n**Q9.** Why need h(w,y) to be square-integrable? \n\n**A.** If $h(w, y)$ is not square-integrable, it cannot be solved using regression methods, as regression approaches minimize the squared loss that involves the second-order moment of the solution. Besides, our method and analysis are built upon the completeness assumption, which is fundamental in proximal causal inference and is imposed on the square-integrable function class.\n\n**Q10.** What does \"feasibility of solutions\" mean in Line 204?\n\n**A.** It refers to the capability of regression methods to achieve consistency for integral solving (Lines 200-203).", "claims_and_evidence": "Most of the theoretical claims are clear in the statement but I have not checked the proofs. There are some unsubstantiated claims.\n1) For e.g. in Section 3, lines 125-127 state that the proposed test is sample-efficient compared to discretization. I don't see any evidence in the form of any comparison for this claim to hold. \n2) \"Besides our power approximates to one as n increases\" (line 422) - Figure 3b does not show the same.", "methods_and_evaluation_criteria": "The paper draws heavily from previous experimental setups in Liu et. al. 2023 which makes comparison easier. However, it is to be noted that these are synthetic datasets and perhaps evaluations on real-world datasets could make the paper stronger.", "theoretical_claims": "No", "experimental_designs_or_analyses": "1) I am not sure why the orange line in Figure 3 does not achieve type I error level in the paper's experiments whereas it does in Liu et. al. 2023. Apart from the sqrt function, the function classes and the noise types are the same in both papers.", "relation_to_broader_scientific_literature": "This work continues the recent thread of research in proximal causal inference where the exchangeability assumption is relaxed by assuming that certain proxy variables are observed. The absence of a direct effect is then translated in terms of conditions on joint distribution of the treatment, outcome and proxy variables. Recently, Miao et. al. 2023 introduced one such condition that concerned the existence of a solution of an integral equation but used a parametric approach to testing. This work introduces a nonparametric approach.", "essential_references_not_discussed": "None that I am familiar with.", "other_strengths_and_weaknesses": "The paper finds an appropriate set of well-motivated assumptions under which a non-parametric test is constructed for bivariate causal discovery assuming access to a proxy variable that affects the outcome. It seems believable that under the nonparametric setup, nonidentifiability is an issue which the paper shows is true even under a linear model. Despite this a condition is proposed under which identifiability is restored. The paper combines multiple ideas that come together for the non-parametric test such as using the characteristic function instead of first order moments like previous work, using a weight function to convert the conditional restriction to an unconditional one. \nRegarding weaknesses, 1) I think the writing of the paper is a major weakness and needs to be improved greatly. Currently, the flow of ideas is abrupt and a few claims/choices don't seem to be well-explained in the text. For e.g., Pg 4, the part before \"Equivalently speaking\" does not seem to be relevant to the final estimate $\\hat{H}^{\\lambda}(w,t)$, the regularization in (5) is not motivated, the term 'bridge function' is never defined in line 172  (see questions for more such instances). \n2) Through each successive approximation, the null hypothesis is being enlarged. While, the paper provides a power analysis, it was not clear to me what whether there are any theoretical guarantees for a type I error control with finite-sample guarantees which makes an earlier claim about sample-efficiency untestable. \n3) Some claims don't seem to be supported by evidence (see claims and evidence section)", "other_comments_or_suggestions": "1. Line 101 - \"allowing test\" should be \"allowing to test\"?\n2. Please specify $\\mathbb{H}_1$ in line 118.", "questions_for_authors": "1) Does (4) hold if (1) does not hold under the assumptions in Thm 4.5? \n2) Like assumption 4.1, does assumption 6.1 say that all variability in U is captured by X? That seems like a very strong assumption and makes a \"confounder\" immaterial. \n3) \"There exists $H^0(w,t)$\" (line 234). Could you explain why this is a \"there exist\" statement. The motivation behind this form of the local alternative seems unclear. \n4) It is not clear to me why in Line 129, you need h(w,y) to be square-integrable? The text prior does not seem to motivate this requriement. \n5) What does \"feasibility of solutions\" mean in Line 204? Feasibility in what sense?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZcfhsejJZk", "reviewer": "Reviewer_g83S", "summary": "This paper aims to test the conditional independence relation 𝑋⊥𝑌∣𝑈 in the presence of an unobserved latent confounder U. Due to the unobservability of U, this conditional independent relation cannot be directly tested. The authors show that conditional independence can still be assessed using a proxy variable Z or W, thereby transforming the hypothesis test into solving an integral equation and evaluating whether the residuals are zero. The paper provides theoretical results on the solution of such an integral equation and an asymptotic analysis of the corresponding test statistic. These theoretical contributions extend previous discretization-based methods to a more robust approach.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your efforts and valuable suggestions in reviewing our paper. We address your concerns below. \n\n**Q1.** The method relies on a partially known prior structure, specifically that $U$ causes both $X$ and $Y$, and that $Z$ and $W$ serve as proxy variables for $U$. This assumption may be restrictive, as in practice, we may not have direct knowledge of the existence of $U$ or whether $Z$ and $W$ are valid proxies.\n\n**A.** This is a standard setting in proximal causal inference. Indeed, since the observational causal inference/discovery typically involves unobserved measurements, it is essential to introduce substitution variables, such as instrumental variables and proxy variables, as considered here. In many scenarios, these assumptions naturally hold and such proxy variables are easy to obtain, for example, they can be noisy measurements of the confounders [1] or as descendants in time-series data [2].\n\n**Q2.** The asymptotic properties of the proposed method depend on a set of assumptions, which may limit its applicability in real-world scenarios.\n\n**A.** The most key assumption is completeness, which is standard in proximal causal inference and can easily hold in our scenario. Others are regularity conditions that are also standard for kernel regression. \n\n**Q3.** The uniqueness of the solution (not only the existence) to the integral equation seems crucial for the hypothesis test. Could the authors provide an intuitive explanation or illustration of its role?\n\n**A.** As claimed in lines 168-176 (right column), uniqueness is not required in our procedure, since our goal is to determine whether a solution exists. Among all solutions to the integral equation, our estimate converges to the least-norm solution (Lemma D.19 in Appendix D.5).  \n\n**Q4.** If $U$ is a vector rather than a single variable, does the proposed theoretical result still hold?\n\n**A.** Our procedure remains valid even when $U$ is a high-dimensional variable, but requires the proxy variable $W$ to satisfy the completeness condition (Assumption 4.1). It implies that the dimension of $W$ is greater than that of $U$. When this holds, the completeness is easy to hold. \n\n**Q5.** The theoretical contribution is intriguing! I am curious whether the proposed method can be viewed as a generalization of Tetrad constraints when considering $Z$ and $W$ as proxy variables.\n\n**A.** Thank you for your insightful question. Our methods can be viewed as a generalization of Tetrad constraints. Specifically, the tetrad constraint was introduced to test whether $X, Y, Z, W$ are conditionally independent given $U$, under linear models. Specifically, the constraint is:\n$$\n\\mathrm{cov}( X,Y ) \\mathrm{cov}( Z,W )= \\mathrm{cov}( X,Z ) \\mathrm{cov}( Y,W )= \\mathrm{cov}( X,W ) \\mathrm{cov}( Y,Z ).\n$$\nWhen proxies $W,Z$ are assumed to satisfy $W \\perp Y|U$ and $Z \\perp X|U$, the constraint degenerates to:\n$$\n\\mathrm{cov}( X,Y ) \\mathrm{cov}( Z,W )= \\mathrm{cov}( X,W ) \\mathrm{cov}( Y,Z ),\n$$\nand it can be used to test $\\mathbb{H}_0: X \\perp Y|U$. In contrast, our constraint for testing $\\mathbb{H}_0$ is formulated by integral equations (1) and (10). Under the linear Gaussian setting with standard normal exogenous noise, our constraint gives rise to the one in \\#. Specifically, suppose $U \\sim \\mathcal{N}(0,1)$, $W=\\mu_UU +\\varepsilon_W$, $Z=\\beta_UU +\\varepsilon_Z$, $X=\\alpha_UU +\\varepsilon_X$, and $Y=\\gamma_UU+\\varepsilon_Y$, where  $\\varepsilon_W,\\varepsilon_Z, \\varepsilon_X, \\varepsilon_Y$ are standard normal. We have $h(W) = \\frac{\\gamma_U}{\\mu_U}W$. Since $\\mathrm{cov}( Z,W ) = \\beta_U \\mu_U$ and $\\mathrm{cov}( Z,Y ) = \\beta_U \\gamma_U$, we have $\\mathrm{cov}( X,Y ) \\mathrm{cov}( Z,W )=\\gamma_U \\alpha_U \\mu_U \\beta_U = \\mathrm{cov}( X,W ) \\mathrm{cov}( Y,Z )$. \n\n[1] Kuroki, M. and Pearl, J. Measurement bias and effect restoration in causal inference. Biometrika, 101(2):423-437, 2014.\n\n[2] Liu, M., Sun, X., Hu, L., and Wang, Y. Causal discovery from subsampled time series with proxy variables. Advances in neural information processing systems, 2023.", "claims_and_evidence": "The claims are generally well-supported by theoretical analysis.", "methods_and_evaluation_criteria": "The proposed method is reasonable and well-motivated.", "theoretical_claims": "The theoretical claims are clearly presented, and the proofs are well-structured.", "experimental_designs_or_analyses": "The experimental design and analyses are generally sound.", "relation_to_broader_scientific_literature": "The proposed method contributes to the problem of CI testing in  the presence of latent confounders.", "essential_references_not_discussed": "The author gives a well-reviewed on the related literature.", "other_strengths_and_weaknesses": "Strengths:\n\n1. This paper addresses an important yet challenging problem in causal discovery involving latent confounders.\n\n2. The proposed method is supported by rigorous theoretical analysis.\n\n3. Experimental results demonstrate that the proposed approach outperforms state-of-the-art methods.\n\n\nWeaknesses:\n\n1. The method relies on a partially known prior structure, specifically that U causes both X and Y, and that Z and W serve as proxy variables for U. This assumption may be restrictive, as in practice, we may not have direct knowledge of the existence of U or whether Z and W are valid proxies.\n\n2. The asymptotic properties of the proposed method depend on a set of assumptions, which may limit its applicability in real-world scenarios.", "other_comments_or_suggestions": "NAN", "questions_for_authors": "1. The uniqueness of the solution (not only the existence) to the integral equation seems crucial for the hypothesis test. Could the authors provide an intuitive explanation or illustration of its role?\n\n2. If U is a vector rather than a single variable, does the proposed theoretical result still hold?\n\n3. The theoretical contribution is intriguing! I am curious whether the proposed method can be viewed as a generalization of Tetrad constraints when considering Z and  W as proxy variables.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper presents a nonparametric method, Proxy Maximum Characteristic Restriction (PMCR), for bivariate causal discovery in the presence of unobserved confounders. By leveraging proxy variables, it reformulates the causal inference problem as an integral equation and proposes a second proxy-based condition for identifiability. 

Pros:

+ Strong theoretical foundation with rigorous asymptotic analysis and well-structured proofs.

+ Innovative use of multiple proxy variables to resolve non-identifiability issues.

+ Experimental results show performance improvements over baseline methods in synthetic settings.

Cons:

+ The writing is dense and very hard to follow, requiring reference to a lengthy appendix.

+ Assumptions such as access to valid proxy variables and certain model structures may limit real-world applicability.
Decision: Accept (poster)