andrey.core.ci_kernel.KCI

class andrey.core.ci_kernel.KCI(data, *, ridge=0.001, rank=None)

Kernel conditional-independence test (KCI) with the gamma approximation of its null.

Each query compares Gaussian-kernel features of x and y after kernel ridge regression on the conditioning set removes what it explains. The statistic is matched to a gamma law through its null mean and variance, which gives the p-value in closed form.

Parameters:

• data (ndarray) – Observed data, one row per sample. Each column is standardized to unit variance.

• ridge (float) – Regularization of the kernel ridge regression on the conditioning set. Smaller values fit the conditioning set more closely.

• rank (int | None) – With a conditioning set, replace each kernel matrix by a pivoted incomplete Cholesky factorization of at most this many columns, so a query’s time grows with the square of the number of rows instead of its cube. The p-value then approximates the exact one. None computes the exact kernels.

Raises:

ValueError – If data is not 2-D, has fewer than 2 rows, contains NaN or infinity, or has values too large or too small to standardize, ridge is not positive and finite, or rank is below 1.

Notes

How it works. With Kx the centered Gaussian kernel matrix of a variable and n rows, an empty conditioning set gives T = tr(Kx Ky) / n, with null mean tr(Kx) tr(Ky) / n^2 and variance 2 tr(Kx^2) tr(Ky^2) / n^4: the gamma HSIC test of hsic_test(). A conditioning set S joins x and S into one variable xs and residualizes with R = ridge (Ks + ridge I)^-1: A = R Kxs R, B = R Ky R, T = tr(A B) / n. Under the null, T is a weighted sum of chi-squared variables whose mean sum_t A_tt B_tt / n and variance 2 sum_st (A_st B_st)^2 / n^2 come from the elementwise product of A and B. The gamma law with that mean and variance gives the p-value [1].

Every kernel is Gaussian, exp(-|a - b|^2 / (2 w^2)), on standardized columns. Without a conditioning set, w is 0.8 below 200 rows, 0.5 below 1,200, and 0.3 beyond. With k conditioning variables, the kernels of xs, y, and S share w = c sqrt(k), c 1.2 below 200 rows, 0.7 below 1,200, and 0.4 beyond, and S enters xs at half scale. x joins the conditioning set and y does not, so of the two variables the one with the lower tie_ranks() rank, a rank computed from each column’s own data, takes the role of x, and the conditioning set’s columns enter its kernels in rank order: the p-value depends on neither the argument order nor the column order.

With rank, each kernel of a conditional query is a pivoted incomplete Cholesky factor L with K ~ L L' [2], centered, and R = I - Ls (Ls' Ls + ridge I)^-1 Ls', so A and B are products of n x rank factors and no n x n system is solved. The factorization stops before rank columns when L L' reproduces the kernel to rounding.

An exact query with a conditioning set runs on a CUDA GPU, in float64, under the cuda backend, or under auto from ANDREY_KCI_GPU_THRESHOLD rows squared. Its p-value then matches the CPU’s to rounding.

Limitations. Measured in PC on simulated data from sparse random graphs (mean degree 2), each variable a weighted sum of z + sin(z), or of cos(z), over its parents plus Gaussian noise, with 20 datasets per size on 4 cores. From 5 to 20 variables with 500 to 2,000 rows, and at 5 variables with 5,000 rows, its median errors (SHD) stayed below the empty graph’s at every size. The median fit took 4 to 18 seconds at 500 rows, 37 seconds to 11 minutes at 2,000 rows, and 13 minutes at 5 variables with 5,000 rows; at 20 variables and 2,000 rows, on shared nodes, 4 of 40 fits passed a half-hour cap. Each query with a conditioning set costs two products of n x n matrices, plus an eigendecomposition for a conditioning set not seen before and two more products for a new y with it; the test keeps up to 256 MB of each kind of matrix, at least one. So time grows as the cube of the number of rows and memory as its square: thousands of rows are practical, tens of thousands are not. The kernel widths follow the row count and the size of the conditioning set, not the data. On simulated draws with up to three conditioning variables and 500 rows, it rejected a true independence at rates between alpha and twice alpha; a test that rejects too often makes PC and FCI keep edges.

References

[1]

Zhang, K., Peters, J., Janzing, D., and Scholkopf, B. (2011). Kernel-based conditional independence test and application in causal discovery. UAI, 804-813.

[2]

Fine, S., and Scheinberg, K. (2001). Efficient SVM training using low-rank kernel representations. Journal of Machine Learning Research 2, 243-264.

property pin_blas: bool

on fewer than 2,000 rows.

Type:

Whether the skeleton scans run the test on one BLAS thread