Numpy overflow with discrete data

[Matyasch]

Dear maintainers, thank you for the great package.

Package version:
0.1.3.8

Problem desciprtion:
I have encountered the following issue when using the chisq CI test. If the number of nodes and/or the number of discrete classes is sufficiently large, then integer representations of numpy (int64) can overflow, resulting in an error. In particular, np.prod(cardSXY) at line 291 in the _Fill3DCountTable function can overflow, returning the result as negative (related stackoverflow question). Then, in the function _Fill3DCountTableByBincount, the product cardS * cardX * cardY at line 248 also overflows, resulting in a negative number being passed to the minlength argument of np.bincount. This throws the following error: ValueError: 'minlength' must not be negative.

Possible solution:
Python uses "unlimited integers", making them better suited for this scenario. Thus, using prod from the math package instead, at line 291, as prod(cardSXY.tolist()), solves this issue. I did not notice any performance issues yet with this solution.

Let me know what you think,
Mátyás

[MarkDana]

Dear @Matyasch, thank you for your insightful observation and even providing a possible solution!

Indeed, when e.g., cardinalities like cardS * cardX * cardY are sufficiently large overflow can occur. However, I'm not entirely sure this is the root cause. In a dataset with N samples, the maximum value for cardS, cardX, cardY would be N (assuming no repeated values at all, e.g., when chisq is mistakenly run on a continuous dataset). In this case cardS * cardX * cardY = N^3, while np.cbrt(np.iinfo(np.int64).max) = 2097152, meaning that this overflow should only be triggered with over 2 million samples. If that's the situation, instead of modifying the code, we might want to warn users that the data isn't suitable for a chisq test.

Could you kindly check if this is the case with your data? If not, could you please provide a minimal code snippet and dataset so we can reproduce the issue, or print the corresponding cardS, cardX, cardY values? We'll investigate further -- it's possible the overflow is occurring earlier, e.g., during the np.bincount step. It's also interesting to hear that using math.prod solves the problem in your case.

Many thanks :)

[jdramsey]

Just a side comment--I've solved this problem elsewhere using AD trees.

[Matyasch]

Dear @MarkDana, I am confused about how the number of samples relates to cardS, cardX or cardY. I can reproduce the issue with a small amount of binary samples as long as the conditioning set size is large enough, causing cardS to explode:

import numpy as np
from causallearn.utils.cit import CIT

p = 63  # number of variables, 62 would still work as 2*2*2^60 is still representable by np.int64
N = 100  # number of samples, can be even smaller
# generate random binary data
np.random.seed(0)
data = np.random.binomial(1, 0.5, (N, p))
# perform test
cit = CIT(data, "chisq")
cit(0, 1, list(range(2, p)))  # test if X0 _||_ X1 | X2, ..., Xp

Printing np.prod(cardSXY) at line 291 writes -9223372036854775808, which is exactly the negative of 2^63, the number we want.

[MarkDana]

Hi @Matyasch thank you for your example, and apologies for the late reply.

Last time I was wrong as cardS here is a list of numbers, not a single large number.

The issue is now fixed by restricting 0 < np.prod < ...' (https://github.com/py-why/causal-learn/pull/200). As you can check from _Fill3DCountTableByBincount' and `_Fill3DCountTableByUnique', the large number of conditioned variables is originally supported.