The kNN estimator returns the complete graph

[kslote1]

With information="knn" on the package's own linear-Gaussian benchmark (n=5, T=500, max_lag=1), every run returns all 25 directed pairs including self-loops, regardless of the ground truth:

seed	true edges	predicted	TP	FP	runtime
0	0	25	0	25	137 s
1	5	25	5	20	131 s
2	3	25	3	22	129 s
3	3	25	3	22	129 s

The unit test confirms this — test_standard_knn in causationentropy/tests/test_data_integration.py fails:

Standard KNN Estimate: TPR: 1.0, FPR: 1.3571428571428572
FAILED (assert fpr == 0)

Two likely contributing defects in knn_conditional_mutual_information (conditional_mutual_information.py:143):

• Line 201: cdist(JS, JS, metric="minkowski", p=k+1) conflates the neighbour count k with the Minkowski order p. With discover_network's default k_means=5 this silently computes an order-6 Minkowski distance. These are independent parameters.
• The docstring states two mutually contradictory (and both incorrect) identities: I(X;Y|Z) = I(X,Y) - I(X,Y;Z) and I(X;Y|Z) = I(X;Y) - I(X;Y|Z).

[chaewon-research]

I dug into this on the current main (ae9868f) and was able to reproduce the complete-graph behavior. I then tried separating a few things that looked like they might be contributing to it: the distance metric, the CMI non-negativity clamp, and whether Z_init is kept during the backward pass.

I used the same general setup as test_standard_knn (T=200, n=5, seed 42, max_lag=2, k=5, alpha 0.01), with 200 shuffles for these diagnostic runs. I kept the existing percentile-based permutation test unchanged.

Configuration	Predicted pairs	TP	FP	FN	Self-loops
Current behavior (Euclidean)	25	6	19	0	5
Chebyshev only	5	5	0	1	0
Chebyshev + unclamped CMI	5	5	0	1	0
Chebyshev + Z_init in backward	6	6	0	0	0
Chebyshev + unclamped CMI + Z_init	6	6	0	0	0
Euclidean + unclamped CMI + Z_init	10	6	4	0	1

The biggest difference came from the metric. Switching from Euclidean to Chebyshev by itself removed all of the false positives in this benchmark, although one true edge (0, 2) was then dropped. Removing the CMI clamp didn’t change the Chebyshev result here.

I also noticed that the standard forward pass uses Z_init, but the backward pass doesn’t keep it in the conditioning set. When I kept Z_init during backward elimination, the missing (0, 2) edge was retained and the benchmark recovered all 6 true edges with no false positives.

So at least for this benchmark, it looks like the complete-graph behavior is mainly driven by the distance metric choice, while the missing Z_init in backward seems to be a separate issue affecting true-edge retention.

I haven’t changed the library source yet. Would you prefer that I keep #14 focused on the kNN metric behavior and open a separate issue for the Z_init/backward-pass problem?

[chaewon-research]

I also reran this using the original benchmark from the issue (T=500, max_lag=1, seeds 0-3) to see if the metric result held up there.

The current Euclidean behavior reproduced the complete-graph failure for all four seeds: all 25 directed pairs were returned each time, including all 5 self-loops.

Using Chebyshev instead of Euclidean eliminated the complete-graph behavior in these runs, but it did not solve the accuracy problem overall:

Seed	True edges	Predicted	TP	FP	FN
0	0	1	0	1	0
1	5	3	3	0	2
2	3	1	1	0	2
3	3	0	0	0	3

So I’d qualify my earlier conclusion a bit. Euclidean seems to be a major contributor to the complete-graph / false-positive behavior, but Chebyshev by itself is not a complete fix because it misses a substantial number of true edges on the original benchmark.

I’m going to hold off on changing the implementation until we decide how to scope this.

[nil-vn]

I summarized the additional diagnostics and experiment results in the attached Markdown file for easier review.

Main findings:

• Chebyshev removes the complete-graph pathology.
• Remaining forward recall is strongly sensitive to k.
• Keeping Z_init in backward improves true-edge retention.
• The observed=0, threshold=0, p=1.0, Pass=True case looks like a separate permutation-test issue.
• Compute_TPR_FPR also appears to have a separate diagonal/self-loop reporting issue.

The full experiment configurations and results are attached below.

issue_14_investigation_summary.md

[Parva-26]

Thanks @chaewon-research and @nil-vn. I reproduced the seed-42 table in section 3 of the summary on ae9868f with the same settings (1000 shuffles): current behaviour 6/19/0 with 5 self-loops, Euclidean + empirical p 1/0/5, Chebyshev + empirical p 5/0/1 (TP/FP/FN). Three small additions that might help with scoping; I'm not planning code changes here.

1. Size of the effect at the test level. I called shuffle_test directly on fully independent X, Y, Z (N=200, 2-D Z, k=5, 100 shuffles, alpha=0.05, 200 repetitions):

estimator / metric	rejected	observed CMI == 0	threshold == 0
knn / euclidean	200/200 (100%)	200/200	200/200
knn / chebyshev	13/200 (6.5%)	111/200	0/200
gaussian	8/200 (4.0%)	0/200	0/200

With Euclidean, the raw CMI for a truly dependent pair (Y = 0.6X + noise, true value 0.154) averaged -0.097 and was negative in 98% of repetitions. So the clamp also turns almost every real dependence into 0, not only the independent cases. With Chebyshev it averaged +0.136.

2. Duplicate points. When a point has at least k identical copies in the joint space, eps is 0, the strict < count gives -1, and digamma(0) is -inf (SciPy 1.17). knn_mutual_information then returns inf and knn_conditional_mutual_information returns nan. This happens, for example, on Poisson(3) count data with k=3, and the wrapper passes both through unchanged (it only clamps finite values). This might be worth a test once the estimator scope is decided.

3. PR #12 (kd_tree). Its KD-tree path uses the Euclidean norm regardless of metric, so whichever option is chosen here (Chebyshev enforced, or Chebyshev as the default) also needs to be applied there. I left the details in a review on #12.

<details> <summary>Script for 1 and 2 (run from the root of main)</summary>

"""
Extra checks for issue #14. Run from the root of a checkout of main (ae9868f):
python issue14_extra_checks.py            # uses the current directory
python issue14_extra_checks.py /path/to/checkout

"""
import importlib

import sys

import warnings
import numpy as np

import scipy

from scipy.special import digamma
sys.path.insert(0, sys.argv[1] if len(sys.argv) > 1 else ".")

warnings.filterwarnings("ignore")

pkg = importlib.import_module("causationentropy")

D = importlib.import_module("causationentropy.core.discovery")

C = importlib.import_module(

"causationentropy.core.information.conditional_mutual_information"

)

M = importlib.import_module("causationentropy.core.information.mutual_information")

print("causationentropy imported from:", pkg.file, "| scipy", scipy.version)
R, N = 200, 200

rng = np.random.default_rng(123)
print("\n[1] Raw kNN CMI (k=5, Z is 2-D, 200 repetitions), mean +- s.e.")

truth = 0.5 * np.log(1 + 0.6**2)

for metric in ("euclidean", "chebyshev"):

indep, dep = [], []

for _ in range(R):

X = rng.standard_normal((N, 1))

Z = rng.standard_normal((N, 2))

Y0 = rng.standard_normal((N, 1))

Y1 = 0.6 * X + rng.standard_normal((N, 1))

indep.append(C.knn_conditional_mutual_information(X, Y0, Z, metric=metric, k=5))

dep.append(C.knn_conditional_mutual_information(X, Y1, Z, metric=metric, k=5))

indep, dep = np.array(indep), np.array(dep)

se = lambda v: v.std(ddof=1) / np.sqrt(R)  # noqa: E731

print(

f"    {metric:9s} independent: {indep.mean():+.3f} +- {se(indep):.3f} (true 0)"

f" | Y=0.6X+noise: {dep.mean():+.3f} +- {se(dep):.3f} (true {truth:.3f}),"

f" negative in {np.mean(dep < 0):.0%} of repetitions"

)
print("\n[2] shuffle_test on fully independent X, Y, Z (100 shuffles, alpha=0.05)")

for info, metric in (

("knn", "euclidean"),

("knn", "chebyshev"),

("gaussian", "euclidean"),

):

passed = zero_obs = zero_thr = 0

for _ in range(R):

X = rng.standard_normal((N, 1))

Y = rng.standard_normal((N, 1))

Z = rng.standard_normal((N, 2))

obs = C.conditional_mutual_information(X, Y, Z, method=info, metric=metric, k=5)

out = D.shuffle_test(

X,

Y,

Z,

obs,

alpha=0.05,

n_shuffles=100,

rng=rng,

information=info,

metric=metric,

k_means=5,

)

passed += bool(out["Pass"])

zero_obs += obs == 0

zero_thr += out["Threshold"] == 0

print(

f"    {info}/{metric:9s}: rejected {passed}/{R} ({passed / R:.1%}); "

f"observed CMI == 0 in {zero_obs}/{R}; threshold == 0 in {zero_thr}/{R}"

)
print("\n[3] Points with >= k exact duplicates in the joint space (eps = 0)")

print("    digamma(0) =", digamma(0))

rng = np.random.default_rng(0)

Xp, Yp, Zp = [rng.poisson(3, (300, d)).astype(float) for d in (1, 1, 2)]

print(

"    Poisson(3) counts, k=3: knn MI =",

M.knn_mutual_information(Xp, Yp, k=3),

"| knn CMI =",

C.knn_conditional_mutual_information(Xp, Yp, Zp, metric="euclidean", k=3),

)

Xd = np.repeat(np.arange(10.0), 20)[:, None]

print(

"    each value repeated 20 times, k=3: knn MI =",

M.knn_mutual_information(Xd, Xd[::-1], k=3),

"| knn CMI =",

C.knn_conditional_mutual_information(Xd, Xd[::-1], Xd, metric="euclidean", k=3),

"| via conditional_mutual_information: Z=None ->",

C.conditional_mutual_information(Xd, Xd[::-1], None, method="knn", k=3),

", with Z ->",

C.conditional_mutual_information(Xd, Xd[::-1], Xd, method="knn", k=3),

)

</details>