Conceptio › Archive › arXiv CS

arXiv CSopen access

Causal Discovery via Transformed Low-Rank Quantile Surfaces

· arxiv_cs

arXiv CS · Papers · License: Open Access

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Causal Discovery via Transformed Low-Rank Quantile Surfaces Ryo Kamimura1,3

arXiv:2609.16931v1 [stat.ME] 15 Sep 2026

Thong Pham2,1,3 *

Shiga University 2 The University of Osaka 3 RIKEN AIP [email protected] [email protected]

Abstract We propose Low-Rank Quantile Surfaces (LRQS), a bivariate causal model in which, in the causal direction, an unknown monotone transformation of the conditional quantile surface admits a low-rank functional decomposition. LRQS subsumes location-scale noise models and post-nonlinear heteroscedastic noise models, while allowing multiple quantile bases to represent changes beyond locationscale effects. We prove generic identifiability of LRQS: the transformed quantile surface is low rank in the causal direction, whereas reverse representability under the corresponding constraints occurs only for exceptional, fine-tuned cause marginals. We provide a simple-yet-powerful causal score using a nonparametric fitting procedure that alternates between rank-constrained approximation of discretized quantile surfaces and isotonic estimation of the unknown monotone transformation. Experiments on synthetic mechanisms with higher-rank distributional shape variation and strong nonlinear distortions, together with standard bivariate benchmarks, show that LRQS is especially effective when conditional distributional shape or observation distortion goes beyond existing location-scale assumptions.

Introduction

Inferring causal direction from observational data requires asymmetry: the conditional distribution in the causal direction should admit a simpler description than the one in the anticausal direction. Classical approaches instantiate this principle through structural restrictions such as additive noise models (ANMs) [Hoyer et al., 2009], post-nonlinear models (PNLs) [Zhang and Hyvärinen, 2009], and location-scale noise models (LSNMs) [Immer et al., 2023]. These models impose location or location-scale structure, either directly or after an invertible transformation. They are identifiable because their assumptions are generally not preserved under reversal, but their expressiveness is limited: in the latent scale, real conditional distributions may vary not only in location and scale, but also in skewness, tail behavior, and other shape features. Quantile-based causal discovery provides a natural way to model distributional asymmetry beyond conditional means. A recent line of work based on quantile partial effects (QPE) assumes that the derivative of the conditional quantile surface with respect to the conditioning variable lies in a finite span of known basis functions [Chen et al., 2026]. This is an expressive observational restriction, but it is imposed on the quantile slope field rather than on the quantile surface itself. Moreover, in their identifiability theory, the finite basis is fixed in advance and is not designed to absorb an unknown monotone observation transformation such as those in PNLs. We propose Low-Rank Quantile Surfaces (LRQS), a bivariate causal model that addresses this limitation. Let QY |X=x (u) denote the conditional quantile function of Y | X = x. LRQS assumes that, * Corresponding author.

Figure 1: Observed and latent quantile surfaces. First row: an LSNM Y = 0.4X + (0.4 + 0.2X 2 )ε. Second row: a PNL-HNM Y = exp{0.4X + (0.4 + 0.2X 2 )ε}. (a) Joint samples of X and Y ; (b) forward observed quantile surface QY |X=x (u); (c) forward latent quantile surface ĥ(QY |X=x (u)); (d) reverse observed quantile surface QX|Y =y (u); (e) reverse latent quantile surface. In the causal direction, the observed surface becomes low rank after a monotone unwarping. in the causal direction, there exists an unknown increasing transformation h = g −1 such that K X bk (x)qk (u). h QY |X=x (u) = a(x) + k=1

Thus the observed quantile surface need not be low rank; instead, it becomes low rank after an unknown monotone unwarping. The case K = 1 recovers a post-nonlinear heteroscedastic noise model (PNL-HNM), while larger K captures richer changes in conditional shape beyond locationscale variation in the latent scale. We prove that this transformed low-rank structure is generically identifiable. In the causal direction, the transformed conditional quantile surface is low rank by construction. For every fixed finite number L of reverse non-intercept components and each prescribed reverse quantile basis, the compatible cause log-densities have local restrictions belonging to a finite-dimensional family, and thus are exceptional. A corresponding result with a free reverse basis holds for L = 1. Neither result requires L = K. Our theory leads to a simple-yet-powerful nonparametric causal score. We estimate conditional quantile matrices in both directions and fit the LRQS structure by alternating between rankconstrained approximation and isotonic estimation of the unknown monotone transformation. The direction with the smaller reconstruction error is selected as causal. Figure 1 gives a visual summary of the LRQS method. Our contributions are: • We introduce LRQS, a transformed low-rank conditional quantile model that extends location-scale and post-nonlinear heteroscedastic noise models. • We establish generic identifiability by characterizing reverse-compatible cause marginals through finite-dimensional local restrictions on their log-densities. This gives, to our knowledge, the first generic identifiability result for the nondegenerate bivariate PNL-HNM class with unknown transformation and reverse noise distribution. • We provide a practical causal discovery algorithm based on low-rank quantile-surface fitting and isotonic estimation of the unknown transformation. • We show through experiments that LRQS is particularly effective on mechanisms with higher-rank conditional shape variation and strong nonlinear observation distortions, while remaining competitive on standard bivariate benchmarks. 2

Related work

Quantile-based causal discovery exploits asymmetries beyond the conditional mean. Bivariate quantile causal discovery uses multiple conditional quantile levels and an independence-of-mechanisms description-length principle [Tagasovska et al., 2020]. The closely related quantile partial effect (QPE) framework assumes that the derivative of the conditional quantile surface with respect to the conditioning variable, i.e. ψ(x, y) = Qx (x, FY |X=x (y)), lies in a fixed finite span in x and y coordinates [Chen et al., 2026]. LRQS differs by imposing low rank on the transformed conditional quantile surface in x and u coordinates: both the quantile bases and the monotone unwarping are learned rather than fixed in advance. The QPE and LRQS model classes overlap nontrivially, and neither theory subsumes the other: (a) ANMs and LSNMs belong to both frameworks; (b) LRQS with an unknown distortion, such as PNL or PNL-HNM, is not contained in QPE with any basis fixed in advance; and (c) conversely, there are QPE models that cannot be expressed as a transformed finite-rank representation for any fixed K in the LRQS framework. Recent alternatives include optimal-transport and velocity-field criteria, such as DIVOT and causal velocity models, which characterize causal asymmetry through transport dynamics or score/velocity equations [Tu et al., 2022, Xi et al., 2025], as well as neural generative or likelihood-based bivariate methods [Goudet et al., 2018, Immer et al., 2023]. These methods are complementary, but LRQS is targeted at an interpretable nonparametric quantile-surface score with explicit transformed low-rank identifiability guarantees. See Appendix A for more detailed discussions.

The model

In the introduction we motivated LRQS as a transformed low-rank structure in the conditional quantile surface. We now give the formal model definition: ! K X Y = g a(X) + bk (X)qk (U ) , (1) k=1

where U ∼ Unif(0, 1), U ⊥ ⊥ X, and g is a continuous, strictly increasing function with inverse h := g −1 . For a fixed K, the transformation g, the basis {qk }K k=1 , the shift function a(x), and the coefficient functions bk (x) are unknown. LRQS generalizes post-nonlinear heteroscedastic noise models. When each qk is a continuous, function, qk is the inverse CDF of some random variable εk , and thus Y = increasing PK g a(X) + k=1 bk (X)εk for some generally dependent noises εk = qk (U ). When K = 1, we get Y = g(a(X) + b(X)ε). This class is sometimes called a post-nonlinear heteroscedastic noise model (PNL-HNM) [Chen et al., 2026]. It contains the LSNM [Immer et al., 2023] when g(z) = z. Beyond location and scale. Multiple quantile bases {qk }K k=1 allow the latent conditional distributions to vary in shape, including skewness, kurtosis, and tail behavior, rather than only in location and scale. See Appendix B for a discussion. The following assumption makes u the conditional quantile level: increasing the latent noise rank increases the response, so the structural function can be read directly as a conditional quantile function. PK Assumption 1. For every fixed x, the function z(x, u) := a(x) + k=1 bk (x)qk (u) is continuous and strictly increasing in u ∈ (0, 1). This assumption holds, for example, if each bk (x) is positive and each basis function qk is continuous and strictly increasing. The following proposition, whose proof is in Appendix C, shows the low-rank structure of the transformed quantile surface. 3

Proposition 1. Denote by QY |X=x (u) the conditional quantile function of Y | X = x. Under Assumption 1, the conditional quantile surface is ! K X Q(x, u) := QY |X=x (u) = g a(x) + bk (x)qk (u) . (2) k=1

PK Equivalently, h Q(x, u) = a(x) + k=1 bk (x)qk (u). Pd A surface H(x, u) has rank at most d if it can be written as H(x, u) = j=1 Aj (x)Bj (u). Proposition 1 therefore implies that h(Q(x, u)) has rank at most K + 1, with the additional component corresponding to the intercept a(x). Evaluating this surface on a grid gives a matrix with the same rank upper bound, which motivates our fitting procedure. The observed surface Q(x, u), by contrast, need not have finite rank. Remark. Y as defined in Eq. (1) is still a valid random variable without Assumption 1, for example, when there exists some x such that z(x, u) is not increasing in u. However, in that case the induced conditional quantile surface need not retain the transformed low-rank structure above.

Identifiability

Fix the forward conditional density r(y | x) = pY |X (y | x) and write ξ(x) = log pX (x). Changing ξ preserves the forward LRQS mechanism but changes the reverse conditional distribution through Bayes’ rule. We study which cause marginals also make the reverse quantile surface Q← (y, u) := QX|Y =y (u) belong to a specified LRQS class. We use K and L for the forward and reverse numbers of non-intercept components, respectively, with no relation imposed between them. The following definition formalizes the set of such tuned marginals: Definition 4.1 (Reverse-compatible marginal set). Fix the forward conditional density r(y | x). For a backward model class M, define n o B(r; M) := ξ : pX,Y (x, y) = eξ(x) r(y | x) and the induced Q← (·, ·) ∈ M . (3) Notion of genericity. For an interval I, write ξ|I for the restriction of ξ to I. Our results place such restrictions of reverse-compatible log-densities in common finite-parameter families. The interval may depend on ξ, but the family depends only on the fixed forward mechanism and the specified reverse class. We use generic identifiability in this local finite-dimensional sense. Here, “local” describes the restriction on the exceptional marginals, not the identification of the causal direction. For each reverse-compatible distribution under consideration, the following two assumptions are required at the same interior point (x0 , y0 ). This point may vary between distributions. Write u0 = FX|Y =y0 (x0 ) and αy0 (x) = ∂y log r(y | x)|y=y0 , with primes on αy0 denoting differentiation in x. Assumption 2 (Local score variation). The conditional y-score has a nonzero derivative in x at x0 : αy′ 0 (x0 ) ̸= 0. Thus αy0 is locally invertible near x0 , allowing us to recover a reverse quantile curve from the Bayes identity used in the proofs. Assumption 3 (Local regularity). In neighborhoods of the corresponding arguments, the densities, conditional quantile maps, and functions defining the reverse representation are C 3 , and the relevant densities are strictly positive. These conditions ensure that the logarithmic derivatives and inverse conditional quantiles used below ← are well-defined locally, and that Q← u (y, u) = 1/pX|Y =y (Q (y, u)) > 0 locally. The general idea in both Theorems 1 and 2 below is that, at a fixed conditioning value y0 , a reverse quantile curve determines the corresponding conditional density of X. Bayes’ rule then recovers the shape of the cause density from this conditional density and the fixed forward mechanism. Thus, restricting one reverse quantile curve at y0 already restricts the cause density on an interval. 4

Main message. We establish these local restrictions for any prescribed finite reverse quantile basis and, when L = 1, for an unrestricted reverse basis. The reverse monotone transformation is unrestricted in both cases. 4.1

General L, fixed basis

We consider the following backward model class: ( ! ) L X ← ← L ML,fixed = Q : Q (y, u) = g̃ ã(y) + b̃k (y)q̃k (u) for some g̃, ã, {b̃k }k=1 .

(4)

k=1

The reverse basis q̃1 , . . . , q̃L : (0, 1) → R is prescribed, while g̃, ã, b̃1 , . . . , b̃L remain unrestricted, subject to the model assumptions. The key object in the analysis is the backward drift ratio R(y, u) = ∂y Q← (y, u)/∂u Q← (y, u), which measures how the u-th backward quantile of X | Y = y moves as the conditioning value y ← changes, normalized by Q← ∈ ML,fixed , the outer monotone distortion cancels from u . When Q this ratio, so R exposes constraints imposed by the inner low-rank structure: PL ã′ (y) + k=1 b̃′k (y)q̃k (u) g̃ ′ (z ← (y, u))∂y z ← (y, u) = , (5) R(y, u) = ′ ← PL ′ g̃ (z (y, u))∂u z ← (y, u) k=1 b̃k (y)q̃k (u) PL where z ← (y, u) = ã(y) + k=1 b̃k (y) q̃k (u). Theorem 1 (Prescribed reverse quantile basis). Fix the forward conditional density r, a finite integer L ≥ 1, and a reverse quantile basis {q̃k }L k=1 . Every reverse-compatible cause log-density ξ ∈ B(r; ML,fixed ) satisfying both Assumptions 2 and 3 at some interior point agrees, on some nonempty open interval, with a member of a common family parameterized by at most 2L + 3 real numbers. This family depends only on r and the prescribed reverse basis. Proof sketch. Fix an interior point as in the assumptions and write x(u) = Q← (y0 , u). At the fixed conditioning value y0 , Eq. (5) expresses R(y0 , ·) using 2L+1 coefficient values. A common nonzero scaling leaves the ratio unchanged, so at most 2L parameters are needed. Bayes’ rule gives Ru (y0 , u) = c0 − αy0 (x(u)), where c0 = (log pY )′ (y0 ). By Assumption 2, αy0 is locally invertible. Thus those parameters, together with c0 , determine the reverse quantile curve x(u) and its derivative. A second application of Bayes’ rule gives ξ(x(u)) = ℓ0 − log r(y0 | x(u)) − log xu (u), where ℓ0 = log pY (y0 ). Since xu > 0, this determines ξ on the corresponding x-interval. Counting the 2L coefficient parameters, c0 , ℓ0 , and y0 gives the upper bound 2L + 3. The full proof is given in Appendix C.2. Remark 1. Theorem 1 allows an arbitrary unknown reverse distortion g̃, which is not covered by the prescribed-basis QPE theory [Chen et al., 2026]. 4.2

L = 1, free basis

We next consider reverse representations with one non-intercept component, allowing their quantile basis to vary: n o M1,free = Q← : Q← (y, u) = g̃ ã(y) + b̃(y)q̃(u) for some g̃, ã, b̃, q̃ . (6) The proof strategy of Theorem 1 does not apply directly because R(y0 , ·) need not belong to a finite-parameter family when q̃ is unrestricted. The key identity is given in the following proposition: Proposition 2. When L = 1, ∂u R + T (u)R = ρ(y), where T (u) = q̃ ′′ (u)/q̃ ′ (u) and ρ(y) = b̃′ (y)/b̃(y). 5

(7)

Proof. Since R = ã′ (y) + b̃′ (y)q̃(u) / b̃(y)q̃ ′ (u) , ∂u R = −q̃ ′′ (u)/q̃ ′ (u)R + b̃′ (y)/b̃(y). Theorem 2 (Free reverse quantile basis for L = 1). Fix the forward conditional density r. Every reverse-compatible cause log-density ξ ∈ B(r; M1,free ) satisfying both Assumptions 2 and 3 at some interior point agrees, on some nonempty open interval, with a member of a common family parameterized by at most nine real numbers. This family depends only on r; no reverse quantile basis is prescribed. Proof sketch. The key observation is that the unknown reverse basis q̃ enters Eq. (7) only through T (u), which does not depend on y. Fix y0 as in the assumptions and set x(u) = Q← (y0 , u), R(u) = R(y0 , u), and D(u) = Ry (y0 , u). Equation (7) and its y-derivative give Ru + T R = ρ0 and Du + T D = ρ1 , where ρ0 = ρ(y0 ) and ρ1 = ρ′ (y0 ). Multiplying the second identity by R and subtracting D times the first eliminates T : RDu − DRu = ρ1 R − ρ0 D. Bayes’ rule expresses Ru in terms of x. Differentiating the Bayes identity in y before restricting to y0 also expresses Du in terms of x, R, and xu . Substitution into the identity above determines xu on a suitable subinterval, giving a first-order system for (x, R, D) with no unknown quantile-basis function remaining. Under Assumptions 2 and 3, local uniqueness determines the solution from three initial values and the four constants c0 = (log pY )′ (y0 ), c1 = (log pY )′′ (y0 ), ρ0 , and ρ1 . Bayes’ rule then recovers ξ on the corresponding x-interval after specifying ℓ0 = log pY (y0 ). Neither the system nor the reconstruction depends explicitly on u, so the origin of the local u-parameter need not be counted. Including ℓ0 and y0 gives at most 3 + 4 + 2 = 9 real parameters. The full proof is given in Appendix C.3. Remark 2. The elimination in the proof sketch relies on the separation of u and y into T (u) and ρ(y) in Eq. (7). Our argument does not provide a corresponding elimination for multiple reverse components, so the free-basis result is limited to L = 1. Remark 3. To our knowledge, when both candidate directions are modeled as PNL-HNM (K = L = 1), Theorem 2 gives the first generic identifiability result for the nondegenerate PNL-HNM class Y = g(a(X) + b(X)ε). Earlier identifiability results cover important special cases, including ANMs (g(z) = z, b ≡ 1) [Hoyer et al., 2009], PNL models (b ≡ 1) [Zhang and Hyvärinen, 2009], and LSNMs (g(z) = z) [Immer et al., 2023]. Scope of the assumptions. Assumptions 2 and 3 require a common interior point where ∂x ∂y log pX,Y (x0 , y0 ) = αy′ 0 (x0 ) ̸= 0. This condition excludes uniform-noise location-scale representations in either direction, even with varying scale and an unknown monotone transformation, as well as constant-scale exponential and Laplace noise models. In the ANM special case, the same nonvanishing condition appears in classical differential-equation arguments for generic identifiability [Hoyer et al., 2009, Peters et al., 2014]. Our assumptions also exclude pure-scale power-law noise models. These restrictions concern the scope of the theorems and do not imply that every excluded mechanism is nonidentifiable. See Appendix C.4 for details.

Causal scoring by transformed quantile-surface fitting

For a chosen number Kfit of non-intercept components, we consider the following causal score: !!2 Z 1 Kfit X SX→Y = inf EX Q(X, u) − g a(X) + bk (X)qk (u) du, (8) a,bk ,qk , g: increasing, z(x,·): increasing

where z(x, u) := a(x) +

k=1

PKfit

k=1 bk (x)qk (u).

Suppose the forward conditional quantile surface belongs to the LRQS class with K ≤ Kfit . Then the causal population score satisfies SX→Y = 0. 6

In the reverse direction, the same fitted component count corresponds to L = Kfit . Theorem 1 constrains exact reverse representations with a prescribed basis for every fixed finite L, whereas Theorem 2 allows a free reverse basis when L = 1, under the stated assumptions. For a population score whose reverse fitting class is covered by the relevant theorem, assume additionally that zero score implies exact membership in that class. Then the reverse score is positive outside the corresponding reverse-compatible marginal set. The reverse-compatible marginals satisfy the finitedimensional local restrictions established above, so the reverse score is positive for generic cause marginals in the stated sense. The population score above, as well as the theoretical results in the previous section, is defined for continuous conditional quantile surfaces. The discretization described below is introduced only for finite-sample estimation. For a candidate direction X → Y , we first sort the observations by the conditioning variable X and partition the sorted observations into G bins with sizes as equal as possible. Thus, our implementation uses approximately equal-count binning. For the j-th bin Ij , the observed conditional quantile b Y |X∈I (ul ), for j = 1, . . . , G, and l = 1, . . . , B. Given the matrix is defined by [Qobs ]j,l = Q j G × B observed quantile matrix, we estimate the population causal score by alternating between a low-rank approximation of the latent surface and estimation of the monotone transformation. The PKfit shift function a(x) becomes a length-G vector a, while a(x) + k=1 bk (x)qk (u) becomes a G × B matrix Z. Further implementation details are provided in Appendix E.1. The algorithm LowRank computes an empirical approximation to SX→Y , and the algorithm Bivariate-LRQS chooses the direction with the smaller causal score. Algorithm 1 LowRank(Q, Kfit ) Data: G × B matrix Q, number of non-intercept components Kfit , inner and outer loop iterations Tin and Tout Result: Approximation score s Initialize Z = Q + Σ, where Σ is a random Gaussian perturbation matrix. Repeat for t = 1, . . . , Tout : • Find an increasing function g such that g(Z) best approximates Q under the Frobenius norm: flatten the matrices Q and Z, sort the pairs (z, q) by the z-values, and perform isotonic regression. • Update Z by an approximate inverse of g: Z[j, l] ← ĥ(Q[j, l]), where ĥ denotes the inverse mapping approximated by linear interpolation from the isotonic fit. • Repeat this inner loop for τ = 1, . . . , Tin to project Z onto the low-rank structure: + a ← RowMeans(Z) and subtract the intercept function: Zc ← Z − a1⊺ . + Use truncated SVD to update Zc : Zc ← SVDKfit (Zc ), and rebuild Z ← Zc + a1⊺ . + Make each row of Z an increasing sequence by isotonic projection. + Rescale Z so that its overall mean is 0 and variance is 1. Estimate the final increasing function g from the flattened pairs (z, q) of Z and Q by isotonic regression, and calculate the score s = ∥Q − g(Z)∥F . Return the score s.

Algorithm 2 Bivariate-LRQS Data: Data matrix of (X, Y ), number of non-intercept components Kfit , number of initializations m Standardize the data. Build the conditional quantile matrix QX→Y of Y | X = x and compute the minimum score over m independent runs to reduce sensitivity to local optima: sX→Y = min1≤i≤m LowRank(QX→Y , Kfit ). Build the conditional quantile matrix QY →X of X | Y = y and compute the minimum score over m independent runs: sY →X = min1≤i≤m LowRank(QY →X , Kfit ). If sX→Y < sY →X , output X as the parent; otherwise, output Y as the parent.

Experiments

6.1

Experimental setup

Custom benchmarks: To clarify the theoretical limitations of existing methods and evaluate the robustness of our proposed approaches, we generated custom benchmarks consisting of two categories, as summarized in Table 1. For each benchmark variant, we generated 100 independent cause-effect pairs with a sample size of n = 1000. The cause variable x is sampled from a Gaussian distribution N (0, 2). For each model, we employed three types of noise components ε, all normalized to have √ √ mean 0 and variance 1: Gaussian N (0, 1), uniform U (− 3, 3), and beta, where V ∼ Beta(2, 2) √ is scaled as (V − 0.5) × 20. 1. Structural complexity (Kgen > 1): To evaluate cases with complex causal mechanisms without PKgen nonlinear observational distortion, the DGP is y = a(x) + k=1 bk (x)qk (u), where a ∼ GP and bk (x) = |fk (x)| with fk ∼ GP. The noise basis functions qk (u) are standardized odd polynomials. We set Kgen ∈ {2, 3, 4, 5}, resulting in 12 variants (4 component counts × 3 noise types). This benchmark isolates conditional shape variation without an observation-level transformation, allowing comparison of the methods in the absence of nonlinear measurement distortion. 2. Strong nonlinear distortions (PNL): The DGP is the PNL-HNM model y = g(z), where z follows either an ANM, z = f (x) + ε, or an LSNM, z = f (x) + σ(x)ε. We consider five transformations: (a) identity, g(z) = z; (b) cube, g(z) = z 3 ; (c) sigmoid, g(z) = 1/(1 + exp(−z)); (d) exponential, g(z) = exp(z); and (e) hyperbolic tangent, g(z) = tanh(z). This category consists of 30 variants. Table 1 displays representative results for the Gaussian-noise settings, including the exponential and hyperbolic-tangent distortions. The reported Avg. acc. and Time (s), however, are macro-averaged over all 42 custom benchmark variants (12 structural-complexity variants and 30 PNL variants). Complete results for the 42 variants are provided in Appendix D.1. Existing benchmarks: We evaluate all 24 bivariate benchmark datasets used in Chen et al. [2026]. Table 2 reports 12 representative datasets: (i) AN and LS from Tagasovska et al. [2020]; (ii) SIM and SIM-c from Mooij et al. [2016]; (iii) Cha and Net from Guyon et al. [2019]; (iv) Per and Sig from Xi et al. [2025]; (v) Qd-V and NN-V from Chen et al. [2026]; (vi) Tue [Mooij et al., 2016]; and (vii) D4-s1 [Marbach et al., 2009]. Although Table 2 displays these 12 datasets, its Avg. acc. and Time (s) are macro-averaged over all 24 benchmark datasets. Baselines and evaluation protocol: We compare LRQS against ANM [Zheng et al., 2024], DIVOT [Tu et al., 2022], CVEL [Xi et al., 2025], QPE-k, and three QPE-f variants, which are fixed, poly, and lowrank, from Chen et al. [2026]. For Table 1, we report results obtained by running these baselines on the same custom benchmark datasets. For the existing-benchmark comparison in Table 2, we rerun all methods rather than citing benchmark values from prior work. No dataset-specific hyperparameter tuning is performed for this main comparison; instead, a single default configuration for each method is used across all 24 datasets. Implementation and computational details are provided in Appendix E. Proposed method variants: We evaluate two configurations of LRQS: fix-g, which fixes the transformation to the identity mapping (g(z) = z), and est-g, which jointly estimates the unknown increasing transformation g together with the low-rank quantile surface. Unless otherwise stated, both variants use the default settings described in Appendix E.1. 6.2

Main results

Tables 1 and 2 summarize the causal direction identification accuracy and average execution time for the custom and existing benchmarks, respectively. Lightweight. Both fix-g and est-g variants are very fast compared to neural-based methods such as QPE-f. Robustness to structural complexity. As shown in Table 1, QPE-k and the QPE-f variants struggle to achieve consistently high accuracy under multirank structural complexity (Kgen > 1), whereas CVEL performs strongly for larger Kgen . In contrast, both of our proposed methods (fix-g and est8

Table 1: Accuracy comparison on custom benchmarks highlighting theoretical limitations of existing methods. The displayed columns show representative Gaussian-noise settings, while Avg. acc. and Time (s) are averaged over all 42 custom benchmark variants. Structural complexity (g(z) = z, Kgen > 1)

Strong nonlinear distortions (PNL)

Overall

Kgen = 2

Kgen = 3

Kgen = 4

Kgen = 5

AN (Exp)

AN (Tanh)

LS (Exp)

LS (Tanh)

Avg. acc. (all 42)

Time (s)

ANM DIVOT CVEL QPE-k QPE-f (fixed) QPE-f (poly) QPE-f (lowrank)

0.44 0.04 0.59 0.36 0.34 0.33 0.36

0.48 0.00 0.91 0.13 0.28 0.18 0.37

0.43 0.00 0.98 0.09 0.40 0.17 0.63

0.40 0.00 1.00 0.11 0.37 0.15 0.64

0.21 0.00 1.00 0.16 0.68 0.27 0.61

0.27 0.31 0.70 0.10 0.13 0.07 0.31

0.26 0.00 0.94 0.32 0.90 0.53 0.85

0.33 0.08 0.66 0.45 0.47 0.44 0.54

0.421 0.217 0.730 0.561 0.443 0.390 0.535

7.335 1.261 2.055 0.062 34.536 32.588 37.966

LRQS (est-g) LRQS (fix-g)

0.94 0.95

0.97 0.97

0.96 0.99

1.00 1.00

0.91 0.59

1.00 1.00

0.96 0.78

0.952 0.893

2.222 0.050

Method

Table 2: Accuracy and average computational time (seconds per pair) on 12 representative bivariate benchmark datasets using default parameter settings. Avg. acc. and Time (s) are averaged over all 24 benchmark datasets. Method

SIM

SIM-c

Cha

Net

Per

Sig

Qd-V

NN-V

Tue

D4-s1

Avg. acc. (all 24)

Time (s)

ANM DIVOT CVEL QPE-k QPE-f (fixed) QPE-f (poly) QPE-f (lowrank)

1.00 1.00 0.25 0.99 1.00 0.97 0.51

0.42 0.72 0.18 1.00 0.98 0.99 0.52

0.74 0.73 0.64 0.83 0.74 0.75 0.74

0.79 0.70 0.62 0.79 0.70 0.66 0.66

0.73 0.52 0.67 0.60 0.54 0.50 0.55

0.73 0.80 0.51 0.89 0.89 0.91 0.63

0.63 0.90 1.00 0.77 0.95 0.96 0.83

0.28 0.61 0.69 0.89 0.47 0.54 0.73

0.82 0.37 0.83 0.42 0.76 0.64 0.73

0.64 0.40 0.88 0.53 0.75 0.72 0.77

0.55 0.46 0.34 0.54 0.61 0.61 0.48

0.58 0.75 0.42 0.58 0.54 0.54 0.33

0.609 0.612 0.602 0.741 0.739 0.744 0.631

22.354 4.057 2.561 0.067 35.502 31.672 3.862

LRQS (est-g) LRQS (fix-g)

0.98 1.00

0.96 1.00

0.61 0.69

0.52 0.77

0.61 0.70

0.71 0.85

0.64 0.74

0.54 0.52

0.57 0.73

0.61 0.79

0.69 0.78

0.50 0.42

0.654 0.740

2.707 0.053

g) maintain near-perfect accuracy across the evaluated values of Kgen . Notably, fix-g achieves this robustness while operating at a speed comparable to the fastest baselines. Fitted rank and effective complexity. The fitted number of components need not match the generating number. In the Gaussian benchmark with Kgen = 5, both variants achieve 99–100% accuracy across Kfit ∈ {2, 3, 4}. Spectral analysis shows that the first two components capture over 99% of the forward row-centered quantile matrix’s energy in the median case, helping to explain the effectiveness of the default Kfit = 2. Increasing the fitted rank further can reduce accuracy, suggesting that Kfit should be viewed as a regularization parameter (Appendix D.2). Resistance to strong nonlinear distortions. Under strong nonlinear distortions such as exp and tanh, several baseline methods experience substantial degradation in accuracy. While our fix-g also struggles with extreme saturation (e.g., tanh), est-g successfully estimates and unwarps these distortions, sustaining high accuracy. Across all 42 custom benchmark variants, est-g achieves the highest average accuracy, further supporting the benefit of learning the monotone unwarping under strong observation-level distortions. Performance on existing benchmarks. On the standard datasets (Table 2), fix-g achieves perfect accuracy on AN and LS and the highest accuracy among all evaluated methods on the real-world Tue dataset. Across all 24 benchmarks, fix-g remains competitive with the strongest baselines under the common default-parameter protocol. The additional flexibility of est-g is particularly beneficial under strong nonlinear observation distortions. 6.3

Additional experimental analyses

Sensitivity to optimization and estimation settings. Additional initializations yield only modest accuracy gains, and reconstruction scores and directional decisions stabilize after a few outer iterations in a representative example (Appendix D.3). Analyses of sample size, grid resolution, and quantile range show stable performance near the default settings, although overly coarse discretization reduces accuracy (Appendix D.4). 9

Directional score separation. Across the 42 custom benchmark settings, normalized gaps between forward and reverse LRQS scores are typically separated from zero. Incorrect decisions tend to have smaller absolute gaps than correct decisions, indicating weaker directional separation (Appendix D.5). Role of the low-rank constraint. Removing the effective rank constraint yields zero reconstruction scores in both directions for every pair in the two evaluated settings. The low-rank restriction is therefore essential for directional discrimination in these experiments (Appendix D.6). Non-monotone transformations. Both variants achieve perfect accuracy for the tested values Kfit ∈ {1, 2, 3, 4} across the evaluated family of non-monotone transformations. These results demonstrate empirical robustness to this particular violation of the increasing-transformation assumption (Appendix D.7).

Limitations and conclusions

Limitations. LRQS is a bivariate method and does not address multivariate graphs, latent confounding, or selection bias. Its generic-identifiability analysis characterizes exceptional cause marginals through local restrictions. The theory covers any prescribed finite reverse quantile basis, while the free-basis result is proved for one reverse component. The practical score also depends on conditional quantile estimation, the chosen number of non-intercept components Kfit , discretization, and a nonconvex alternating fit. Conclusion. We introduced Low-Rank Quantile Surfaces, a model in which the causal-direction conditional quantile surface becomes low rank after an unknown monotone transformation. LRQS extends location-scale and post-nonlinear heteroscedastic noise models, possesses generic identifiability, and leads to a simple nonparametric causal score. Experiments show strong performance under higher-rank shape variation and nonlinear observation distortions, suggesting latent low-rank quantile structure as an effective inductive bias for causal discovery.

Acknowledgment This work was partially supported by the Japan Science and Technology Agency (JST) under CREST Grant Number JPMJCR22D2 and by the Japan Society for the Promotion of Science (JSPS) under KAKENHI Grant Number JP24K20741. References Yikang Chen, Xingzhe Sun, and Dehui Du. Causal discovery via quantile partial effect. In The Fourteenth International Conference on Learning Representations, 2026. Anish Dhir, Samuel Power, and Mark Van Der Wilk. Bivariate causal discovery using Bayesian model selection. In Proceedings of the 41st International Conference on Machine Learning, ICML’24. JMLR.org, 2024. Olivier Goudet, Diviyan Kalainathan, Philippe Caillou, Isabelle Guyon, David Lopez-Paz, and Michele Sebag. Learning functional causal models with generative neural networks. In Explainable and interpretable models in computer vision and machine learning, pages 39–80. Springer, 2018. Isabelle Guyon, Alexander Statnikov, and Berna Bakir Batu. Cause effect pairs in machine learning. Springer, 2019. Patrik O. Hoyer, Dominik Janzing, Joris Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems 21, pages 689–696. Curran Associates Inc., 2009. Alexander Immer, Christoph Schultheiss, Julia E Vogt, Bernhard Schölkopf, Peter Bühlmann, and Alexander Marx. On the identifiability and estimation of causal location-scale noise models. In International Conference on Machine Learning, pages 14316–14332, 2023. 10

Mariyam Khan, Shohei Shimizu, and Thong Pham. Beyond additivity: Causal discovery in locationscale noise models with hidden variables. arXiv preprint, 2026. doi: 10.48550/arXiv.2606.08196. David Lopez-Paz, Krikamol Muandet, and Benjamin Recht. The randomized causation coefficient. The Journal of Machine Learning Research, 16(1):2901–2907, 2015. Daniel Marbach, Thomas Schaffter, Claudio Mattiussi, and Dario Floreano. Generating realistic in silico gene networks for performance assessment of reverse engineering methods. Journal of computational biology, 16(2):229–239, 2009. Joris M Mooij, Jonas Peters, Dominik Janzing, Jakob Zscheischler, and Bernhard Schölkopf. Distinguishing cause from effect using observational data: methods and benchmarks. Journal of Machine Learning Research, 17(32):1–102, 2016. Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 15:2009–2053, 2014. Thong Pham, Shohei Shimizu, Hideitsu Hino, and Tam Le. Scalable counterfactual distribution estimation in multivariate causal models. In Proceedings of CLeaR 2024, 2024. Thong Pham, Takashi Nicholas Maeda, and Shohei Shimizu. Causal additive models with unobserved causal paths and backdoor paths. In Proceedings of AISTATS 2026, 2026. Paul Rolland, Volkan Cevher, Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear additive noise models. In Proceedings of the 39th International Conference on Machine Learning, pages 18741–18753, 2022. Hirofumi Suzuki, Kentaro Kanamori, Takuya Takagi, Thong Pham, Takashi Nicholas Maeda, and Shohei Shimizu. I-CAM-UV: Integrating causal graphs over non-identical variable sets using causal additive models with unobserved variables. In Proceedings of AAAI 2026, 2026. Natasa Tagasovska, Valérie Chavez-Demoulin, and Thibault Vatter. Distinguishing cause from effect using quantiles: bivariate quantile causal discovery. In Proceedings of the 37th International Conference on Machine Learning, ICML’20. JMLR.org, 2020. Ruibo Tu, Hedvig Kjellstrom, Kun Zhang, and Cheng Zhang. Optimal transport for causal discovery. In International Conference on Learning Representations, 2022. Johnny Xi, Hugh Dance, Peter Orbanz, and Benjamin Bloem-Reddy. Distinguishing cause from effect with causal velocity models. In Proceedings of the 42nd International Conference on Machine Learning, ICML’25. JMLR.org, 2025. Hiroshi Yokoyama, Ryusei Shingaki, Kaneharu Nishino, Shohei Shimizu, and Thong Pham. Causaldiscovery-based root-cause analysis and its application in time-series prediction error diagnosis. In 2025 International Joint Conference on Neural Networks (IJCNN), pages 1–10, 2025. K. Zhang and A. Hyvärinen. On the identifiability of the post-nonlinear causal model. In Proc. 25th Conference on Uncertainty in Artificial Intelligence (UAI2009), pages 647–655, 2009. Yujia Zheng, Biwei Huang, Wei Chen, Joseph Ramsey, Mingming Gong, Ruichu Cai, Shohei Shimizu, Peter Spirtes, and Kun Zhang. Causal-learn: Causal discovery in Python. Journal of Machine Learning Research, 25(60):1–8, 2024.

Detailed related work

Bivariate cause-effect identification. Because X → Y and Y → X impose no distinct conditional-independence constraints, bivariate observational discovery requires additional asymmetry assumptions [Mooij et al., 2016, Guyon et al., 2019]. Transport, velocity, and score-based approaches. DIVOT connects functional causal models to optimal transport and derives directional criteria from transport dynamics [Tu et al., 2022]. Causal velocity models (CVEL) treat the cause as a time-like parameter, relating counterfactual velocity fields to the joint score function [Xi et al., 2025]. For multivariate nonlinear additive Gaussian-noise models, SCORE recovers a causal ordering from the Hessian of the joint log-density [Rolland et al., 2022]. Generative, supervised, and likelihood-based methods. Causal generative neural networks (CGNNs) compare generative fits using maximum mean discrepancy [Goudet et al., 2018], whereas the randomized causation coefficient learns causal decisions from labeled cause-effect pairs [LopezPaz et al., 2015]. LOCI provides feature-map and neural-network estimators for LSNMs [Immer et al., 2023], and QPE-f estimates quantile partial effects using normalizing flows [Chen et al., 2026]. Bayesian model selection offers a complementary approach: Dhir et al. [2024] compare flexible bivariate models by marginal likelihood under priors expressing independent causal mechanisms. Broader connections. Additive models address unobserved causal and backdoor paths [Pham et al., 2026] and graph integration across non-identical variable sets [Suzuki et al., 2026]; locationscale models additionally accommodate heteroscedasticity with hidden variables [Khan et al., 2026]. These approaches build on regression and residual-independence criteria. Related downstream tasks include root-cause analysis of time-series prediction errors [Yokoyama et al., 2025] and optimaltransport-based counterfactual distribution estimation [Pham et al., 2024], which targets outcome distributions rather than causal direction.

Beyond latent location-scale mechanisms

For K = 1, h(Q(x, u)) = a(x) + b(x)q(u): after removing location and scale, every latent conditional quantile curve u 7→ h(Q(x, u)) has the same standardized shape q(u). Hence conditional skewness, tail ratios, and other location-scale invariants (when the corresponding moments exist) cannot vary with the cause x. (A nonlinear g can still make the observed quantile curve u 7→ Q(x, u) change shape with x, but only in the restricted way expressible by a PNL-HNM, a location-scale family pushed through one fixed nonlinearity.) For K = 2, h(Q(x, u)) = a(x) + b1 (x)q1 (u) + b2 (x)q2 (u): the relative coefficient λ(x) = b2 (x)/b1 (x) can vary with x, so the standardized latent conditional quantile curve u 7→ h(Q(x, u)) can change by shape variations beyond location and scale on this latent scale. For example: • if q1 is a symmetric Gaussian quantile and q2 is an asymmetric Gumbel quantile, varying λ(x) changes conditional skewness; • if q1 is Gaussian and q2 a heavy-tailed Student-t quantile, varying λ(x) changes relative tail thickness.

Theoretical details

C.1

Proof of Proposition 1

PK For fixed x, define Zx (u) = a(x) + k=1 bk (x)qk (u). By Assumption 1, Zx is strictly increasing. Since U ∼ Unif(0, 1), the conditional quantile of Zx (U ) is Zx (u). Because g is strictly increasing, monotone transformations preserve quantile order, so QY |X=x (u) = g(Zx (u)). 12

C.2

Proof of Theorem 1

We first establish two consequences of Bayes’ rule used in both theorem proofs. The first relates the reverse drift ratio to the forward conditional score; the second reconstructs the cause log-density from a reverse quantile curve. Lemma 3. Write c(y) = ∂y log pY (y), and α(x, y) = ∂y log r(y | x). Then Ru (y, u) = c(y) − α(Q← (y, u), y), ←

←

ξ(Q (y, u)) = log pY (y) − log r(y | Q

(y, u)) − log Q← u (y, u).

(9) (10)

Proof. Let F (x, y) = FX|Y =y (x) and p← (x | y) = pX|Y =y (x). Differentiating F (Q← (y, u), y) = ← ← u gives Q← | y) and R = −Fy (Q← , y), where subscripts on F and p← denote partial u = 1/p (Q ← derivatives before evaluation at x = Q← (y, u). Consequently, Ru = −p← | y)/p← (Q← | y). y (Q ← ξ(x) Substituting Bayes’ formula p (x | y) = e r(y | x)/pY (y) yields Eq. (9); taking logarithms in the identity for Q← u gives Eq. (10). Fix an interior point (x0 , y0 ) satisfying Assumptions 2 and 3, and set x(u) = Q← (y0 , u). Choose an interval Jx around x0 on which αy′ 0 remains nonzero, and a quantile interval Iu around u0 such that x(Iu ) ⊂ Jx . Set Ix = x(Iu ). Below, αy−1 denotes the inverse of αy0 restricted to Jx . At y0 , 0 the reverse LRQS formula for R in Eq. (5) involves 2L + 1 coefficients. Its value is unchanged by a common nonzero scaling of these coefficients. Since at least one denominator coefficient is nonzero, normalizing that coefficient leaves at most 2L real parameters. Denote them collectively by θ and the resulting function by Rθ (u). Set c0 = c(y0 ). Equation (9) and the local invertibility of αy0 give x(u) = αy−1 c0 − Rθ′ (u) , u ∈ Iu . 0

(11)

Thus both x(u) and its derivative xu (u) are determined by (θ, c0 ) for fixed y0 . With ℓ0 = log pY (y0 ), Eq. (10) becomes ξ(x(u)) = ℓ0 − log r(y0 | x(u)) − log xu (u). Since xu > 0, the map u 7→ x(u) is invertible locally, so this identity determines ξ|Ix without introducing another unknown function. The local inverse branches are determined by the fixed mechanism r. Counting the at most 2L coefficients in θ, the two scalars c0 , ℓ0 , and the anchor y0 gives the claimed upper bound of 2L + 3 real parameters. C.3

Proof of Theorem 2

Proof. Fix ξ ∈ B(r; M1,free ). Choose an interior point satisfying Assumptions 2 and 3, and let Iu be a sufficiently small quantile interval around the corresponding u0 . By Eq. (9), Ruu (y0 , u) = −αy′ 0 (Q← (y0 , u))Q← u (y0 , u) is nonzero on Iu . Hence R(y0 , ·) is not identically zero there. After choosing a smaller interval, which need not contain u0 , we may therefore assume R(y0 , u) ̸= 0 throughout. Set Ix = Q← (y0 , Iu ). Use the functions x(u) = Q← (y0 , u), R(u) = R(y0 , u), and D(u) = Ry (y0 , u). Set c0 = (log pY )′ (y0 ), c1 = (log pY )′′ (y0 ), ρ0 = ρ(y0 ), and ρ1 = ρ′ (y0 ). Also define βy0 (x) = ∂y2 log r(y | x) y=y . Both αy0 and βy0 are fixed by r and y0 ; primes on these functions denote 0 differentiation in x. To obtain an equation for Du , we differentiate Eq. (9) with respect to y while holding u fixed. Since y enters both arguments of α(Q← (y, u), y), the chain rule gives ← ∂y α(Q← (y, u), y) = αx (Q← (y, u), y) Q← y (y, u) + αy (Q (y, u), y). ′ At y = y0 , the quantile derivative is Q← y (y0 , u) = R(u)xu (u), while αx (x, y0 ) = αy0 (x) and αy (x, y0 ) = βy0 (x). Thus, evaluating the Bayes identity and its y-derivative at y = y0 , and using Du (u) = Ryu (y0 , u), gives

Ru = c0 − αy0 (x),

(12)

Du = c1 − βy0 (x) − αy′ 0 (x)Rxu .

(13)

On the other hand, Eq. (7) and its y-derivative give Ru + T R = ρ0 ,

Du + T D = ρ1 ,

because T depends only on u. Together with Eq. (12), the first identity yields T = (ρ0 − c0 + αy0 (x))/R. Substituting this into the second identity and comparing with Eq. (13) gives R c1 − βy0 (x) − ρ1 + ρ0 − c0 + αy0 (x) D , xu = αy′ 0 (x)R2 (14) Ru = c0 − αy0 (x), ρ0 − c0 + αy0 (x) D Du = ρ1 − . R Thus no unknown quantile-basis function remains in the system. Its right-hand side is C 1 wherever R ̸= 0 and αy′ 0 (x) ̸= 0, as ensured by Assumptions 2 and 3. For fixed y0 and a chosen u∗ ∈ Iu , local uniqueness therefore determines (x, R, D) from the three initial values (x∗ , R∗ , D∗ ) = (x, R, D)(u∗ ) and the four scalars (c0 , c1 , ρ0 , ρ1 ). On the actual solution, xu = Q← u (y0 , u) > 0. With ℓ0 = log pY (y0 ), Eq. (10) becomes ξ(x(u)) = ℓ0 − log r(y0 | x(u)) − log xu (u). Since u 7→ x(u) is locally invertible, this determines ξ|Ix , after shrinking Iu and Ix = x(Iu ) if necessary. Neither the system nor this reconstruction depends explicitly on u. Hence changing the origin of the local u-parameter only reparametrizes the same function of x, and u∗ contributes no additional parameter. Counting (x∗ , R∗ , D∗ , c0 , c1 , ρ0 , ρ1 , ℓ0 , y0 ) gives at most nine real parameters. Allowing their admissible values defines a family depending only on r and containing all the required restrictions. Additional conditions for a valid reverse-compatible joint density can only restrict this family. C.4

Scope and exclusions of the local assumptions

We explain the roles of Assumptions 2 and 3 and derive several cases outside their joint scope. All density derivatives below are evaluated in open neighborhoods where the relevant densities are positive and the required derivatives exist. Implications of Assumption 3. Purely discrete, deterministic, or singular conditional laws are outside this setting. The assumption is local: bounded support or nonsmoothness away from the selected neighborhoods does not by itself violate it. Nor is boundedness of quantile derivatives as u → 0 or u → 1 required. Assumption 2 is symmetric. Write α(x, y) = ∂y log r(y | x), so that αx (x0 , y0 ) = αy′ 0 (x0 ). The two factorizations of the joint density give αx (x, y) = ∂x ∂y log pX,Y (x, y) = ∂y ∂x log pX|Y =y (x).

(15)

Thus Assumption 2 is symmetric in X and Y , despite being stated through the forward conditional density. Moreover, differentiating Eq. (9) in u gives Ruu (y, u) = −αx (Q← (y, u), y) Q← u (y, u).

(16)

Because Q← u > 0 at regular points, the nondegeneracy condition αx (x0 , y0 ) ̸= 0 is equivalent to Ruu (y0 , u0 ) ̸= 0, i.e., nonzero curvature of the reverse drift ratio in the quantile coordinate at the corresponding point. The condition detects dependence within smooth regions of the density. Indeed, on an open rectangle, αx ≡ 0 implies log pX,Y (x, y) = A(x) + B(y) by integration. If the density is positive and C 2 throughout the interior of a rectangular support, this factorization implies independence. Consequently, every dependent distribution in that setting has a point where Assumption 2 holds. For PNL-HNMs. Consider Y = g(a(X) + b(X)ε), with ε ⊥⊥ X, b > 0, and h = g −1 . Work locally where the functions are sufficiently differentiable and h′ > 0. Set e = (h(y) − a(x))/b(x) 14

and ν = log pε . The conditional density satisfies r(y | x) = pε (e)h′ (y)/b(x). Differentiating its logarithm with respect to y, using ey = h′ (y)/b(x), gives α(x, y) =

h′ (y) ′ h′′ (y) ν (e) + ′ . b(x) h (y)

(17)

Differentiating with respect to x, using ex = −(a′ (x) + b′ (x)e)/b(x), then gives αx (x, y) = −

′′ h′ (y) ′ ′ ′ ′ a (x) + b (x)e ν (e) + b (x)ν (e) . b(x)2

(18)

The outer transformation contributes only the positive factor h′ (y), so it cannot remove the degeneracies identified below. By Eq. (15), the same conclusions apply to reverse representations after exchanging the variables. For Gaussian noise. Standard Gaussian noise gives αx (x, y) = h′ (y)(a′ (x) + 2b′ (x)e)/b(x)2 . At a conditioning value where a′ and b′ are not both zero, this expression is nonzero for some noise value. Uniform noise, including varying scale. For uniform noise on any nondegenerate interval, ν ′ = ν ′′ = 0 throughout the support interior. Equation (18) therefore gives αx = 0, regardless of the location function, scale function, or monotone transformation. Hence no point satisfies both assumptions. In the reverse direction, for uniform reverse noise, q̃ is affine, so T = q̃ ′′ /q̃ ′ = 0. Equation (7) therefore gives Ru = ρ(y), hence Ruu = 0. Equation (16) again gives αx = 0. Thus a forward mechanism satisfying Assumption 2 cannot admit a regular uniform-noise reverse PNL-HNM. Constant-scale exponential and Laplace noise. For standard exponential noise, ν ′ (e) = −1 and ν ′′ (e) = 0 on e > 0. Thus αx (x, y) = h′ (y)b′ (x)/b(x)2 . When b is constant, Assumption 2 fails throughout the support interior, even if a is nonconstant. For standard Laplace noise, ν(e) = −|e| − log 2, so ν ′′ (e) = 0 away from e = 0. With constant b, Eq. (18) again gives αx = 0 at all smooth points. At e = 0, the density has a kink and the required local regularity fails. Hence constant-scale Laplace models have no point satisfying both assumptions. Unlike uniform noise, these exclusions need not persist under varying scale. For Laplace noise, away from its kink, αx (x, y) = h′ (y)b′ (x) sign(e)/b(x)2 . Thus exponential or Laplace noise can satisfy Assumption 2 at a regular point where b′ (x) ̸= 0. Pure-scale power-law noise models. For a pure-scale mechanism with a ≡ const, suppose the noise log-density has the form ν(t) = const+γ log t on its positive support. Then tν ′′ (t)+ν ′ (t) = 0, and Eq. (18) gives αx = 0 for every scale function b. Examples include ε ∼ Beta(η, 1), with density pε (t) = ηtη−1 on 0 < t < 1, and Pareto noise with density pε (t) = ηt−η−1 on t > 1, where η > 0. For nonconstant b, these mechanisms can describe dependent variables, but their interior log-densities remain separable in x and y. They therefore fall outside Assumption 2, including after a monotone transformation. Piecewise-affine reverse quantile bases. The exclusion also extends beyond one-component models. Suppose all prescribed reverse bases are piecewise affine with finitely many knots. On any interval avoiding their combined knots, the reverse latent surface has the form z ← (y, u) = A(y) + B(y)u. At regular points B(y) > 0, and R(y, u) = (A′ (y) + B ′ (y)u)/B(y) is affine in u. Hence Eq. (16) gives αx = 0. Genuine knots violate the required smoothness; if a knot disappears in the resulting surface, continuity still gives Ruu = 0 there. Thus exact reverse representations built from such bases cannot satisfy the two assumptions jointly. Connection to classical ANM results. For an ANM Y = a(X) + ε, with ε ⊥⊥ X and ν = log pε , direct differentiation gives αx (x, y) = −a′ (x)ν ′′ (y − a(x)). The differential-equation arguments of Hoyer et al. [2009, Theorem 1] and Peters et al. [2014, Condition 19 and Proposition 21] use this nonvanishing product; their finite-dimensional genericity statements additionally require it to 15

be nonzero for all but countably many cause values along some fixed conditioning slice. These particular genericity results do not cover uniform, exponential, or Laplace noise: ν ′′ vanishes on every smooth piece of the positive-density region, while the required differentiability fails at the Laplace kink. Our condition requires only one regular point satisfying both assumptions, but shares this nondegeneracy ingredient.

Additional experimental results

D.1

Additional results

We provide the complete experimental results across all 42 custom mechanism variants and the full suite of 24 existing bivariate benchmark datasets. For the existing benchmarks, we report both a controlled comparison using a single default configuration for each method and a complementary comparison under dataset-specific hyperparameter tuning. Comprehensive evaluation on structural complexity. Table 3 details the performance under increasing structural complexity (Kgen ∈ {2, 3, 4, 5}) across Gaussian, uniform, and beta noise distributions. LRQS remains consistently strong across the different complexity levels and noise distributions. CVEL is also highly competitive in several higher-complexity settings, whereas QPE-k and the QPE-f variants show larger performance degradation in a number of multirank settings. Overall, the results support the robustness of the low-rank quantile surface representation to substantial conditional shape variation. Consistent robustness against nonlinear distortions. Tables 4 and 5 present the complete results under nonlinear observation distortions for both post-nonlinear additive noise models (PNL-AN) and post-nonlinear heteroscedastic noise models (PNL-HNM). The est-g variant remains highly accurate across a broad range of noise distributions and nonlinear transformations, with particularly strong performance under severe nonlinear distortions. While several competing methods are effective in specific settings, the overall results demonstrate the benefit of explicitly estimating the unknown monotone transformation when observation-level distortions are substantial. Extended results on existing benchmarks. Tables 6 and 7 report the complete results on all 24 bivariate benchmark datasets using the same default-parameter protocol as in Table 2. A single default configuration for each method is used across all datasets, without dataset-specific tuning. Under this controlled setting, LRQS (fix-g) remains competitive with the strongest baselines while retaining its low computational cost. For completeness, Tables 8 and 9 additionally report a dataset-specific tuning comparison following the evaluation style of Chen et al. [2026]. Baseline values are cited from Chen et al. [2026], while LRQS is evaluated using dataset-specific hyperparameter tuning. This complementary comparison shows the performance achievable when hyperparameters are adapted to individual datasets, whereas the default-parameter comparison above provides a more controlled assessment under a common evaluation protocol. Table 3: Detailed accuracy comparison on multirank structural-complexity benchmarks. Gaussian, uniform, and beta denote the noise distributions. Method

Gaussian

ANM DIVOT CVEL QPE-k QPE-f (fixed) QPE-f (poly) QPE-f (lowrank) LRQS (est-g) LRQS (fix-g)

D.2

Uniform

Beta

Time (s)

Kgen = 2

Kgen = 3

Kgen = 4

Kgen = 5

Kgen = 2

Kgen = 3

Kgen = 4

Kgen = 5

Kgen = 2

Kgen = 3

Kgen = 4

Kgen = 5

0.44 0.04 0.59 0.36 0.34 0.33 0.36 0.94 0.95

0.48 0.00 0.91 0.13 0.28 0.18 0.37 0.97 0.97

0.43 0.00 0.98 0.09 0.40 0.17 0.63 0.96 0.99

0.40 0.00 1.00 0.11 0.37 0.15 0.64 1.00 1.00

0.39 0.18 0.43 0.98 0.53 0.54 0.45 0.96 0.98

0.49 0.03 0.80 0.77 0.30 0.35 0.30 0.93 0.84

0.59 0.00 0.98 0.47 0.30 0.31 0.53 0.82 0.69

0.60 0.00 1.00 0.23 0.29 0.37 0.64 0.71 0.73

0.46 0.04 0.62 0.69 0.46 0.49 0.38 0.95 0.94

0.58 0.00 0.94 0.22 0.25 0.29 0.39 0.88 0.88

0.50 0.00 0.98 0.07 0.27 0.30 0.57 0.90 0.89

0.54 0.00 1.00 0.04 0.20 0.21 0.58 0.91 0.98

8.420 1.348 1.994 0.057 37.048 27.550 38.597 2.038 0.046

Sensitivity to fitted rank and effective spectral rank

Sensitivity to the fitted rank. We first examine the sensitivity of LRQS to the fitted number of non-intercept basis functions, denoted by Kfit . We use the same 100 cause-effect pairs from the structural complexity benchmark with Gaussian noise and Kgen = 5, and vary Kfit ∈ {2, 3, 4, 5} 16

Table 4: Detailed accuracy comparison on PNL-AN benchmarks. Gaussian, uniform, and beta denote the noise distributions. Id, Sig, Exp, and Tanh denote identity, sigmoid, exponential, and hyperbolic tangent transformations. Method ANM DIVOT CVEL QPE-k QPE-f (fixed) QPE-f (poly) QPE-f (lowrank) LRQS (est-g) LRQS (fix-g)

Gaussian

Uniform

Beta

Time (s)

Cube

Sig

Exp

Tanh

Cube

Sig

Exp

Tanh

Cube

Sig

Exp

Tanh

1.00 0.96 0.02 0.90 0.83 0.87 0.77 0.90 0.94

0.23 0.00 1.00 0.34 0.35 0.14 0.43 0.99 0.82

0.63 0.89 0.04 0.75 0.24 0.15 0.48 0.95 0.96

0.21 0.00 1.00 0.16 0.68 0.27 0.61 1.00 1.00

0.27 0.31 0.70 0.10 0.13 0.07 0.31 0.91 0.59

1.00 1.00 0.76 1.00 0.27 0.37 0.60 0.93 0.98

0.20 0.00 0.99 0.74 0.13 0.12 0.22 0.95 0.56

0.58 0.91 0.87 1.00 0.30 0.34 0.64 0.94 0.91

0.20 0.00 1.00 0.49 0.55 0.49 0.72 0.98 0.97

0.17 0.24 1.00 0.43 0.23 0.30 0.51 0.93 0.69

1.00 1.00 0.19 1.00 0.33 0.41 0.51 0.96 0.93

0.20 0.00 1.00 0.60 0.26 0.14 0.29 0.96 0.58

0.59 0.93 0.52 0.97 0.17 0.17 0.50 0.95 0.95

0.17 0.00 1.00 0.29 0.70 0.49 0.69 0.98 1.00

0.24 0.33 0.99 0.25 0.24 0.28 0.52 0.89 0.74

5.795 1.127 2.178 0.067 32.326 34.431 37.496 2.356 0.051

Table 5: Detailed accuracy comparison on PNL-HNM benchmarks. Gaussian, uniform, and beta denote the noise distributions. Id, Sig, Exp, and Tanh denote identity, sigmoid, exponential, and hyperbolic tangent transformations. Method ANM DIVOT CVEL QPE-k QPE-f (fixed) QPE-f (poly) QPE-f (lowrank) LRQS (est-g) LRQS (fix-g)

Gaussian

Uniform

Beta

Time (s)

Cube

Sig

Exp

Tanh

Cube

Sig

Exp

Tanh

Cube

Sig

Exp

Tanh

0.36 0.25 0.18 0.94 0.92 0.90 0.59 1.00 1.00

0.34 0.00 0.99 0.48 0.68 0.55 0.73 1.00 0.97

0.35 0.36 0.20 0.79 0.72 0.65 0.56 0.99 0.95

0.26 0.00 0.94 0.32 0.90 0.53 0.85 1.00 1.00

0.33 0.08 0.66 0.45 0.47 0.44 0.54 0.96 0.78

0.35 0.30 0.24 1.00 0.45 0.50 0.36 1.00 1.00

0.36 0.00 0.99 0.76 0.73 0.58 0.75 1.00 0.95

0.34 0.32 0.39 0.94 0.31 0.30 0.38 1.00 0.97

0.32 0.02 0.95 0.65 0.79 0.66 0.76 0.98 1.00

0.36 0.10 0.74 0.60 0.34 0.40 0.46 0.99 0.83

0.41 0.36 0.17 0.98 0.60 0.62 0.41 0.99 0.99

0.31 0.00 0.98 0.63 0.69 0.58 0.82 0.98 0.88

0.36 0.34 0.22 0.91 0.46 0.41 0.37 0.99 0.97

0.32 0.04 0.92 0.39 0.79 0.62 0.75 1.00 1.00

0.34 0.09 0.76 0.53 0.34 0.36 0.52 0.97 0.77

8.006 1.324 1.982 0.060 34.736 34.774 37.930 2.234 0.051

Table 6: Detailed accuracy comparison on the first group of bivariate benchmark datasets. All methods are evaluated using a single default configuration across datasets. Method

AN-s

LS-s

MNU

SIM

SIM-c

SIM-g

SIM-ln

Tue

Cha

Net

ANM DIVOT CVEL QPE-k QPE-f (fixed) QPE-f (poly) QPE-f (lowrank)

1.00 1.00 0.25 0.99 1.00 0.97 0.51

1.00 1.00 0.09 0.88 0.94 1.00 0.17

0.42 0.72 0.18 1.00 0.98 0.99 0.52

0.25 0.34 0.25 0.78 1.00 1.00 0.56

0.25 0.91 0.70 1.00 0.93 0.99 0.89

0.74 0.73 0.64 0.83 0.74 0.75 0.74

0.79 0.70 0.62 0.79 0.70 0.66 0.66

0.68 0.68 0.80 0.83 0.54 0.61 0.55

0.72 0.60 0.65 0.68 0.79 0.84 0.72

0.55 0.46 0.34 0.54 0.61 0.61 0.48

0.73 0.52 0.67 0.60 0.54 0.50 0.55

0.73 0.80 0.51 0.89 0.89 0.91 0.63

LRQS (est-g) LRQS (fix-g)

0.98 1.00

0.81 0.93

0.96 1.00

0.73 0.92

0.74 0.84

0.61 0.69

0.52 0.77

0.71 0.72

0.56 0.64

0.69 0.78

0.61 0.70

0.71 0.85

Table 7: Detailed accuracy comparison on the second group of bivariate benchmark datasets. All methods are evaluated using a single default configuration across datasets. Method

Multi

D4-s1

D4-s2a

D4-s2b

D4-s2c

Per

Sig

Vex

Qd-V

Sig-V

RbF-V

NN-V

ANM DIVOT CVEL QPE-k QPE-f (fixed) QPE-f (poly) QPE-f (lowrank)

0.60 0.36 0.87 0.88 0.76 0.81 0.70

0.58 0.75 0.42 0.58 0.54 0.54 0.33

0.62 0.60 0.47 0.67 0.54 0.54 0.66

0.61 0.59 0.41 0.61 0.53 0.57 0.53

0.58 0.54 0.47 0.64 0.47 0.50 0.55

0.63 0.90 1.00 0.77 0.95 0.96 0.83

0.28 0.61 0.69 0.89 0.47 0.54 0.73

0.17 0.05 0.91 0.63 0.93 0.96 0.97

0.82 0.37 0.83 0.42 0.76 0.64 0.73

0.76 0.46 0.90 0.67 0.69 0.64 0.67

0.47 0.59 0.90 0.68 0.68 0.61 0.69

0.64 0.40 0.88 0.53 0.75 0.72 0.77

LRQS (est-g) LRQS (fix-g)

0.84 0.86

0.50 0.42

0.53 0.64

0.51 0.54

0.51 0.53

0.64 0.74

0.54 0.52

0.78 0.70

0.57 0.73

0.60 0.79

0.44 0.65

0.61 0.79

while keeping the other settings fixed (n = 1000, G = B = 10, Tin = 5, and Tout = 10; for est-g, m = 5). In addition to causal-direction accuracy, we report the normalized directional score gap. For sreverse + sforward > 0, define sreverse − sforward δ= . (19) sreverse + sforward 17

Table 8: Detailed accuracy comparison on the first group of bivariate benchmark datasets under dataset-specific tuning. Baseline values are cited from Chen et al. [2026], while LRQS results are obtained using dataset-specific hyperparameter tuning. Method

AN-s

LS-s

MNU

SIM

SIM-c

SIM-g

SIM-ln

Tue

Cha

Net

ANM DIVOT CVEL QPE-k QPE-f

0.43 0.62 1.00 0.99 1.00

0.47 0.69 0.98 0.88 1.00

0.46 0.45 0.98 1.00 1.00

0.45 0.69 0.93 0.78 0.99

0.40 1.00 0.94 1.00 1.00

0.45 0.68 0.63 0.83 0.88

0.49 0.47 0.72 0.79 0.88

0.41 0.60 0.90 0.83 0.86

0.46 0.63 0.76 0.68 0.92

0.65 0.38 0.64 0.54 0.70

0.41 0.44 0.68 0.60 0.85

0.47 0.49 0.62 0.89 0.86

LRQS (est-g) LRQS (fix-g)

1.00 1.00

0.78 0.99

0.98 1.00

0.75 0.80

0.74 0.80

0.82 0.81

0.83 0.82

0.72 0.82

0.66 0.70

0.84 0.90

Table 9: Detailed accuracy comparison on the second group of bivariate benchmark datasets under dataset-specific tuning. Baseline values are cited from Chen et al. [2026], while LRQS results are obtained using dataset-specific hyperparameter tuning. Method

Multi

D4-s1

D4-s2a

D4-s2b

D4-s2c

Per

Sig

Vex

Qd-V

Sig-V

RbF-V

NN-V

ANM DIVOT CVEL QPE-k QPE-f

0.48 0.34 0.97 0.88 0.96

0.50 0.50 0.67 0.58 0.79

0.48 0.57 0.51 0.67 0.71

0.46 0.55 0.58 0.61 0.62

0.48 0.55 0.58 0.64 0.60

0.49 0.97 1.00 0.77 1.00

0.44 0.82 0.84 0.89 0.90

0.39 0.05 0.96 0.63 0.91

0.49 0.32 0.91 0.42 0.91

0.50 0.44 0.94 0.67 0.91

0.43 0.63 0.92 0.68 0.94

0.48 0.47 0.87 0.53 0.90

LRQS (est-g) LRQS (fix-g)

0.86 0.93

0.50 0.67

0.59 0.73

0.55 0.61

0.54 0.57

0.99 0.99

0.78 0.81

0.87 0.70

0.61 0.76

0.75 0.84

0.62 0.69

0.62 0.79

Positive values favor the true forward direction, negative values favor the reverse direction, and values close to zero indicate weak directional separation. Table 10 summarizes the results. Table 10: Sensitivity of LRQS to the fitted rank Kfit on the same 100 structural complexity pairs with Gaussian noise and Kgen = 5. The score gap is reported as median [Q1 , Q3 ]. “Decision flips” denotes the fraction of pairs whose predicted direction differs from that obtained with Kfit = 2. Method

Kfit

Accuracy

Median δ [Q1 , Q3 ]

Exact ties

Decision flips

fix-g

2 3 4 5

1.00 1.00 1.00 1.00

0.499 [0.402, 0.604] 0.563 [0.462, 0.665] 0.589 [0.502, 0.690] 0.597 [0.492, 0.684]

0 0 0 0

– 0% 0% 0%

est-g

2 3 4 5

1.00 1.00 0.99 0.79

0.543 [0.423, 0.665] 0.616 [0.469, 0.737] 0.652 [0.463, 0.822] 0.871 [0.466, 1.000]†

0 0 0 9

– 0% 1% 21%

† For est-g with K fit = 5, the score-gap summary is computed over the 91 non-tied pairs; δ is undefined for the nine exact

ties with sforward = sreverse = 0.

The fix-g variant is insensitive to the fitted rank over the evaluated range: all 100 pairs are correctly oriented for every Kfit ∈ {2, 3, 4, 5}, with no decision flips relative to Kfit = 2. The est-g variant is also stable for moderate changes in rank, retaining 99–100% accuracy for Kfit ∈ {2, 3, 4}. However, increasing the fitted rank to Kfit = 5 reduces the accuracy to 79%, produces nine exact ties, and changes the predicted direction for 21% of the pairs relative to Kfit = 2. These results show that increasing Kfit does not necessarily improve causal identification. In finitesample estimation, Kfit is therefore better viewed as a regularization parameter controlling the flexibility of the fitted quantile surface than as an estimate of the number of generating components. A conservatively small fitted rank preserves the directional asymmetry in this experiment, whereas excessive flexibility can weaken causal discrimination. Effective spectral rank. We next examine the effective rank of the empirical quantile surfaces in the same structural complexity setting with Gaussian noise and Kgen = 5. For each of the 100 18

1.00

Forward Reverse

Cumulative spectral energy E(J)

0.99 0.98 0.97 0.96 0.95 0.94 0.93

3 Approximation rank J

Figure 2: Cumulative spectral energy of the row-centered empirical quantile matrices for the structural complexity benchmark with Gaussian noise and Kgen = 5. Curves show the median across 100 cause-effect pairs, and error bars indicate the interquartile range. cause-effect pairs, we construct the 10 × 10 empirical quantile matrix in each direction and subtract its row means. Let σ1 ≥ σ2 ≥ · · · denote the singular values of the resulting row-centered matrix. We define the cumulative spectral energy explained by the first J components as PJ 2 j=1 σj E(J) = P 2 . j σj Because this benchmark uses g = id, this analysis is performed directly on the empirical quantile matrices before any rank projection and does not use an estimated unwarping transformation. As shown in Figure 2, the spectral energy is strongly concentrated in the first few components. At J = 2, the median cumulative energy is 0.9927 [0.9902, 0.9944] in the forward direction and 0.9845 [0.9786, 0.9892] in the reverse direction. The median paired difference ∆(2) = Eforward (2) − Ereverse (2) is 0.0070 [0.0027, 0.0139], with ∆(2) > 0 for 85% of the pairs. These results clarify that Kgen = 5 denotes the number of components in the generating mechanism, rather than the effective rank of the empirical quantile matrix. In this setting, the row-centered forward quantile matrix is empirically close to rank two, helping to explain why a small Kfit remains effective under a five-component generating mechanism. Both directions exhibit strong spectral concentration, with a modest but consistent advantage in the forward direction. D.3

Sensitivity to iterations and initialization

Sensitivity to outer iterations. Algorithm 1 uses a finite number of alternating isotonicregression, inverse-update, and rank-projection steps. To examine its empirical sensitivity to the number of outer iterations, we use one representative pair from the LSNM-tanh Gaussian setting and run the est-g variant with 50 different initializations for Tout ∈ {0, 1, 5, 10}. As shown in Figure 3, both forward and reverse reconstruction scores decrease substantially over the first few outer iterations. The median forward score decreases from 3.165 at Tout = 0 to 0.271 19

4.0

Forward Reverse

Reconstruction score

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0

5 Number of outer iterations Tout

Figure 3: Sensitivity of LRQS (est-g) to the number of outer iterations Tout on one representative LSNM-tanh Gaussian pair. Curves show the median reconstruction score across 50 initializations, and error bars indicate the interquartile range.

Table 11: Sensitivity of LRQS (est-g) to the number of initializations m on 100 LSNM-tanh Gaussian pairs. m

Accuracy

Exact ties

1 5 10 50

0.94 0.95 0.96 0.97

0 0 0 0

at Tout = 5, while the corresponding reverse score decreases from 3.745 to 0.590. The changes from Tout = 5 to 10 are comparatively small. The numbers of initializations selecting the correct direction are 45/50, 47/50, 50/50, and 50/50 for Tout = 0, 1, 5, 10, respectively. These results indicate that, in this representative example, the reconstruction scores and directional decisions empirically stabilize after a modest number of outer iterations, although Algorithm 1 does not come with a general convergence guarantee.

Sensitivity to the number of initializations. We next evaluate sensitivity to the number of initializations m using 100 pairs from the same LSNM-tanh Gaussian setting. For each direction, we retain the minimum reconstruction score over the first m initializations, as in Algorithm 2. As shown in Table 11, increasing m yields only modest improvements in accuracy. In particular, increasing the number of initializations tenfold from the default m = 5 to m = 50 improves accuracy by two percentage points, from 95% to 97%. Thus, in this experiment, the default initialization count provides a reasonable balance between accuracy and additional computation. 20

1.00

Causal-direction accuracy

0.95 0.90 0.85 0.80 n = 250 n = 1000 n = 4000

0.75 0.70

20 Number of bins G

Figure 4: Sensitivity of LRQS (est-g) to the sample size n and the number of conditioning-variable bins G on 100 LSNM-tanh Gaussian pairs. The number of quantile levels is fixed to B = 10 with umin = 0.05. Each point reports causal-direction accuracy over the 100 pairs.

D.4

Sensitivity to sample size and quantile-surface discretization

We examine the sensitivity of LRQS to the sample size and the discretization used to construct the empirical conditional quantile surface. All experiments in this subsection use the est-g variant on 100 LSNM-tanh Gaussian cause-effect pairs. Sample size and number of bins. Let umin denote the smallest quantile level used in the grid. We first vary the sample size n ∈ {250, 1000, 4000} and the number of conditioning-variable bins G ∈ {5, 10, 20, 30}, while fixing B = 10 and umin = 0.05. As shown in Figure 4, performance degrades when the discretization along the conditioning variable is too coarse. With G = 5, accuracy ranges from 76% to 83% across the evaluated sample sizes, with 3, 5, and 6 exact ties for n = 250, 1000, 4000, respectively. In contrast, for G ∈ {10, 20, 30} and n ≥ 1000, accuracy remains between 97% and 100%. Thus, the default setting n = 1000 and G = 10 lies within a broader region of high empirical accuracy, although overly coarse binning can degrade directional discrimination. Quantile-grid resolution and quantile range. We next examine the discretization along the quantile axis. Under the default construction ul = (2l−1)/(2B), changing B also changes the outermost quantile level umin = 1/(2B). We therefore interpret this experiment as a comparison of quantilegrid resolutions rather than as an isolated effect of B. With n = 1000 and G = 10, accuracy increases from 87% at B = 5 to 97%, 98%, and 100% at B = 10, 20, and 30, respectively. To separately examine sensitivity to the quantile range, we fix G = B = 10 and vary umin ∈ {0.10, 0.05, 0.025, 0.01}. Table 12 summarizes these analyses. Across the explicitly varied quantile ranges, accuracy remains between 89% and 93% for n = 250 and between 91% and 97% for n = 1000, with no exact ties. A joint higher-resolution setting with n = 1000 and G = B = 30 also achieves 100% accuracy. Overall, these results indicate that 21

Table 12: Sensitivity of LRQS (est-g) to the quantile-grid resolution and quantile range on LSNMtanh Gaussian pairs. Each accuracy is computed over 100 cause-effect pairs. Quantile-grid resolution (n = 1000, G = 10) B

umin

Accuracy

Exact ties

5 10 20 30

0.1000 0.0500 0.0250 0.0167

0.87 0.97 0.98 1.00

2 0 0 0

Quantile-range sensitivity (G = B = 10) n

umin = 0.10

umin = 0.05

umin = 0.025

umin = 0.01

250 1000

0.89 0.91

0.91 0.97

0.90 0.96

0.93 0.96

fix-g est-g

1.75 1.50

Density

1.25 1.00 0.75 0.50 0.25 0.00

0.4

0.2

0.0 0.2 0.4 Normalized directional score gap

0.6

0.8

1.0

Figure 5: Distribution of the normalized directional score gap δ for LRQS (fix-g) and LRQS (est-g) across the 42 custom benchmark settings in a separate diagnostic experiment. Each setting contains 100 cause-effect pairs. The solid vertical line denotes δ = 0, and the dashed lines indicate the descriptive near-tie region |δ| < 0.05. performance is stable over a reasonable neighborhood of the default discretization, while excessively coarse binning or quantile grids can reduce accuracy in this benchmark. D.5

Directional score separation

In a separate diagnostic experiment, we examine the normalized directional score gap δ defined in Eq. (19) across the same 42 custom benchmark settings used in Table 1, with 100 cause-effect pairs per setting. Figure 5 shows that the score-gap distributions are predominantly shifted toward positive values for both variants. The median normalized gap is 0.270 for fix-g and 0.422 for est-g. For descriptive purposes, we define a near tie as |δ| < 0.05. Under this criterion, near ties account for 9.2% of the pairs for fix-g and 4.6% for est-g, with no exact score ties observed for either variant. 22

Table 13: Directional score separation for LRQS in the diagnostic experiment. The near-tie threshold |δ| < 0.05 is used only as a descriptive diagnostic. Method LRQS (fix-g) LRQS (est-g)

Median δ

Near ties

Median |δ| (correct)

Median |δ| (incorrect)

0.270 0.422

9.2% 4.6%

0.302 0.441

0.081 0.075

Table 14: Ablation of the low-rank constraint for LRQS (est-g). With Kfit = 9, the rank truncation is effectively removed for the 10 × 10 empirical quantile matrices. Setting Structural complexity (Kgen = 5, Gaussian) AN-tanh (Gaussian)

Kfit = 2

Kfit = 9

100% correct 88% correct

100/100 exact ties 100/100 exact ties

As shown in Table 13, correctly identified pairs tend to exhibit substantially larger absolute score gaps than incorrectly identified pairs. For fix-g, the median |δ| is 0.302 among correct decisions and 0.081 among incorrect decisions; for est-g, the corresponding values are 0.441 and 0.075. Thus, incorrect decisions in this experiment tend to occur when the forward and reverse reconstruction scores are relatively close, whereas correct decisions generally exhibit clearer directional separation. D.6

Ablation of the low-rank constraint

To assess whether the directional discrimination of LRQS arises from the low-rank constraint rather than from the flexibility of the monotone unwarping alone, we perform an ablation using the est-g variant. With G = B = 10, the row-centered quantile matrix has rank at most 9. We therefore compare the standard setting Kfit = 2 with Kfit = 9, which effectively removes the rank truncation. We evaluate two representative settings from the custom benchmarks: structural complexity with Gaussian noise and Kgen = 5, and AN-tanh with Gaussian noise. Each setting contains the same 100 cause-effect pairs used in the corresponding experiments. As shown in Table 14, removing the effective rank constraint eliminates directional discrimination in both settings. For every pair with Kfit = 9, both directions are perfectly reconstructed, sforward = sreverse = 0, so the causal direction cannot be determined. In contrast, the standard low-rank setting retains clear directional information. This ablation shows that the performance of est-g in these experiments cannot be attributed to the flexibility of the monotone unwarping alone. Rather, the low-rank restriction provides an essential inductive bias for preserving asymmetry between the two candidate directions. D.7

Stress test under monotonicity violations

The identifiability analysis assumes that the observation-level transformation g is strictly increasing. To examine empirical behavior outside this assumed regime, we conduct a stress test using the nonmonotone transformation family 4 z , z ≥ 0, gκ (z) = |z|4+κ , z < 0, where κ ∈ {0, 0.5, 1, 1.5, 2}. The transformation folds the latent variable around zero and is therefore non-monotone for every value of κ. At κ = 0, the transformation is symmetric, whereas increasing κ introduces greater asymmetry between the negative and positive sides. We evaluate this stress test on the structural complexity benchmark with Gaussian noise, Kgen = 5, and n = 1000, using 100 cause-effect pairs. We set G = B = 10 and vary the fitted rank over Kfit ∈ {1, 2, 3, 4, 5}. As shown in Table 15, fix-g achieves perfect directional accuracy across all evaluated values of κ and Kfit . The est-g variant also achieves perfect accuracy for Kfit ≤ 4 across the entire transformation 23

Table 15: Causal-direction accuracy under violations of the increasing-transformation assumption. Results are reported over 100 cause-effect pairs for each setting. Method

Kfit

κ=0

κ = 0.5

κ=1

κ = 1.5

κ=2

LRQS (fix-g)

1 2 3 4 5

1.00 1.00 1.00 1.00 1.00

LRQS (est-g)

1 2 3 4 5

1.00 1.00 1.00 1.00 0.92

1.00 1.00 1.00 1.00 0.93

1.00 1.00 1.00 1.00 0.95

1.00 1.00 1.00 1.00 0.90

1.00 1.00 1.00 1.00 0.86

family. Performance decreases only for the more flexible Kfit = 5 setting, where accuracy ranges from 86% to 95%. Thus, LRQS remains empirically effective under this particular family of non-monotone observation transformations, especially with conservative fitted ranks.

Computational resources and implementation details

We provide the complete details of our computational environments, software versions, and specific implementation patches applied to the baseline codes. E.1

Implementation details of LRQS

Construction of the empirical quantile matrix. For a candidate direction X → Y , we sort the observations by the conditioning variable X and partition the sorted indices into G groups using numpy.array split. This produces approximately equal-count bins whose sample sizes are as equal as possible; when the sample size n is not divisible by G, the bin sizes differ by at most one. Unless explicitly specified otherwise, the B quantile levels are ul =

2l − 1 , 2B

l = 1, . . . , B.

For example, B = 10 gives ul ∈ {0.05, 0.15, . . . , 0.95}. Within each bin, the corresponding empirical response quantiles are computed using numpy.quantile without overriding its default quantile rule. No special tie-aware binning rule is applied. Observations with identical values of the conditioning variable are partitioned according to their positions after sorting and are not explicitly constrained to remain in the same bin. Ties in the response variable are passed directly to numpy.quantile. For nonempty bins, no minimum-size rule, bin merging, or additional smoothing is applied. If an empty bin occurs, its quantile row is copied from the preceding bin; if the first bin is empty, its row is set to zero. The same construction is applied in the reverse direction after exchanging the conditioning and response variables. Default settings. Unless otherwise stated, we use G = 10 bins and B = 10 quantile levels to construct the empirical quantile surfaces, and set the number of non-intercept basis functions to Kfit = 2. In the iterative optimization process (Algorithm 1), the inner loop for rank-constrained projection and the outer loop for alternating isotonic-regression and inverse-update steps are repeated Tin = 5 and Tout = 10 times, respectively. For each additional initialization, we add Gaussian perturbations Σij ∼ N (0, 0.52 ) to the initial surface to reduce sensitivity to initialization during the alternating estimation procedure. 24

E.2

Computational environments

The experiments in this study were conducted across two distinct computing environments depending on the hardware requirements of the evaluated methods. 1. Local CPU environment: All experiments for our proposed LRQS method (both fix-g and est-g) and the CPU-based baselines (ANM, DIVOT, QPE-k) were executed on a local workstation. • OS: Windows 11 • CPU: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80 GHz (4 cores, 8 threads) • Memory (RAM): 8.0 GB (3200 MT/s) Software stack: Python 3.13.6. The core dependencies for executing our methods and local baselines include NumPy (v2.4.2) and scikit-learn (v1.8.0). 2. GPU computing cluster: The experiments involving the GPU-based baselines, CVEL and the three QPE-f variants (fixed, poly, and lowrank), were conducted on a Docker computing cluster to leverage hardware acceleration. • GPU: Single NVIDIA RTX 2080 Ti (11 GB VRAM) Software stack: We utilized the official NVIDIA PyTorch Docker container image (nvcr.io/nvidia/pytorch:23.12-py3), which provides a highly optimized environment containing Python 3.10 and PyTorch 2.2.0. For the baseline methods, we used the official implementations released by the original authors. The evaluation protocol for Table 2 is described in Section 6. E.3

Implementation patches for QPE-f baseline

To ensure a fair and reproducible evaluation of the three QPE-f variants (fixed, poly, and lowrank), we utilized the official source code provided by Chen et al. [2026]. However, we encountered severe numerical instability and compatibility issues when running the original implementation in our GPU environment. Specifically, the model frequently produced −∞ for intermediate QPE scores, leading to identical and trivial predictions across multiple benchmark datasets. To conduct a valid comparative study, we introduced two minimal, mathematically equivalent patches to the official codebase: 1. Type-hint compatibility fix: We updated outdated type hints (from torch.types import Tensor, Device) that caused import errors in newer PyTorch environments. 2. Numerical stabilization for Jacobian computation: We traced the −∞ issue to the denominator computation of the QPE term (∂u/∂y). In the original code, this was computed using the Jacobian-vector product (jvp), which frequently collapsed to zero in our en[ = −(∂u/∂x)/(∂u/∂y) to diverge. We patched this vironment, causing the term QPE by replacing the jvp-based calculation with log abs det jacobian(...).exp(), which retrieves the exact same analytical Jacobian through the normalizing flow’s native method. This alternative implementation is mathematically equivalent, resolving the −∞ collapse and yielding dataset-specific performance scores consistent with expectations.

Record · ID 919419 · SHA-256 979f28c5baddccf5
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.