@thomasjpfan Thanks for your help!
Indeed, I tried to remove the casting of inputs in int64 from this PR and run @jianlingzhong 's example (here: #16803 (comment)), and the problem happens when forming the extension because that one has very high dimensionality (so needs int64)

The problem I have is that in this example, no error is thrown by the cython code, but wrong negative indices are created silently instead (a ValueError: column index exceeds matrix dimensions error is raised after the call of _csr_polynomial_expansion , at concatenation)

Also, there is something I didn't get when tracing what happens in the cython code: I added two print lines here:

print(d)    
expanded_dimensionality = int((d**2 + d) / 2 - interaction_only*d)
print(expanded_dimensionality)

d is 120006, however expanded_dimensionality is 758209071, when it should be 7200660015 ... And I thought it was due to the int32 overflow so I tried to cast 7200660015 into int32 and it returns -2147483648 not 758209071 so I have trouble understanding what's going on here...
=> In fact what happens is that 758209071 = int((np.int32(120006**2) + np.int32(120006))/2 - np.int32(120006)), with an overflow coming from the square power.

It might not be a good idea to use a try ... except... approach as the expansion would take a long time when the number of dimensions is high. It may be a good idea to calculate the resultant dimension number and check if int32 is sufficient.

The problem I have is that in this example, no error is thrown by the cython code, but wrong negative indices are created silently instead (a ValueError: column index exceeds matrix dimensions error is raised after the call of _csr_polynomial_expansion , at concatenation)

We would need to update the Cython to properly error when it notices that int64 is not good enough.

It might not be a good idea to use a try ... except... approach as the expansion would take a long time when the number of dimensions is high. It may be a good idea to calculate the resultant dimension number and check if int32 is sufficient.

I agree. If there is a fast way to compute the dimension, we can quickly decide the output format of the indices from there. At a glance, I think it can be done by using _deg2_column or _deg3_column to compute the largest possible indices. (We would need to explicitly use the int64 versions of these functions)

It might not be a good idea to use a try ... except... approach as the expansion would take a long time when the number of dimensions is high. It may be a good idea to calculate the resultant dimension number and check if int32 is sufficient.

I agree. If there is a fast way to compute the dimension, we can quickly decide the output format of the indices from there. At a glance, I think it can be done by using _deg2_column or _deg3_column to compute the largest possible indices. (We would need to explicitly use the int64 versions of these functions)

I agree, I'll try to compute the largest possible index as you suggested @thomasjpfan :

• For the "polynomial expansion" case, I think we can just find the maximum column index col_max in the indices, and indeed just compute _deg2_column(col_max, col_max) or _deg3_column(col_max, col_max, col_max), with the int64 version of those functions: if it's bigger than 2147483647 (np.iinfo(np.int32)), then we need to convert indices to int64 before doing the expansion
• For the "interaction expansion" case, I think we would need to find the maximum column index col_max as above, then among all rows that contain a column col_max, find the max column number that they contain except col_max, col_max_2 (and then similarly find the subsequent col_max_3), and then decide based on the int64 computation of _deg2_column(col_max, col_max_2) (or `_deg3_column(col_max, col_max_2, col_max_3)``. (I need to double check this)

(A simpler idea (and which wouldn't need to go through all the elements in indices), would be to just look at the expanded_dimensionality and see if it's larger than 2147483647. But it may sometimes be worse than computing the largest possible index as described above, in some weird cases for instance if many last columns of the expanded matrix are empty... I don't know if those cases would arise often however, and if it's worth going through all elements of indices, just to be optimal in those rare cases))

(A more efficient solution instead of going through all indices to find the max, would be to solve on i the equation _deg2_column(i, i)=2147483647 and see if indices contains any index bigger than this, But I don't know if it's worth too)

After that, I think we also need to check whether the indptr should be casted to int64 too: I think in theory the biggest number in indptr is its last element and should be equal to total_nnz (the the number of nonzero elements of the expansion), but I need to double check in case I'm missing some edge cases (for example I'm not sure yet how the code deals with sparse matrix that have duplicates)

But first, I already pushed some casting to int64 needed to avoid overflow in the computation of the expanded dimensionality and the number of nonzero elements (I'm not 100% sure of the syntax for casting in cython, this is how I saw it online and it compiled, but I may be mistaken)

I also see that all current tests pass but the example from @jianlingzhong is not passing, so I'll add it as a test

I'm going for the simplest solution described above (("A simpler idea (and which wouldn't need to go through all the elements in indices), would be to just look at the expanded_dimensionality and see if it's larger than 2147483647"), which I guess is enough to perform well in most cases ?

Regarding the possibility for multiple values that are pointed to the same place in the sparse matrix, there is no such pb because _validate_data takes care of that (if there are such multiple values (duplicates) it sums them, as is standard for sparse matrices)

Also, I'm wondering if there can be problems with the data inside the array (for instance if it's float32 but the product of two float32 is too big to be a float32). Since it seems that overflow fails silently in cython code this could be bad for ML algorithms it would overflow without telling the user... I'll look at this and raise an issue if needed.

I'm also thinking there could be a silent problem if the input array indices are int64, but the expanded dimensionality is even bigger than the max int64.

Also, to test because of some bug in the fit function of PolynomialFeatures, it will be very long to fit_transform in @jianlingzhong's example, but this PR should fix it so I'll add the test after it's merged #19734

I didn't have time to advance yet but I'll get back to it soon

fyi, #19734 was just merged so you should be good to add that test.

fyi, #19734 was just merged so you should be good to add that test.

Thanks for the heads up @frrad ! I'll get started on adding the test

I'm still in progress, but I'm having trouble to explicitly cast a variable with fused type in say an int64 in cython: it seems the syntax var = <np.int64_t>var doesn't work, I'll try to find how to do it, but it would help in case anyone knows

@thomasjpfan thanks a lot for the help ! It's much clearer now how casting works in cython,
Indeed, factoring out the code into its own fused function sounds good, I'll try to do that this week

in 90e12b8, similarly to @thomasjpfan 's suggestion I separated:

• the part which decides whether we need casting or not on one hand, in the python code (in _cast_to_int64_if_needed)
• the part that does the actual expansion in the fused _csr_polynomial_expansion function on the other hand, , in the cython code.

(@thomasjpfan in #19676 (comment) you suggested to do the expansion outside of _csr_polynomial_expansion but I started extracting the casting check outside of it and then it just felt natural to keep _csr_polynomial_expansion as it is, let me know what you think)

I went for a simple rule for casting: if the expanded dimensionality needs to be an int64, then I convert all indices and inptr to int64. This may not be optimal since maybe the actual expanded matrix is so sparse that the expanded indptr (related to the number of nonzeros elements in each row), could be an int32 array, however there are actually no way to know this in advance: to know that we would need to actually compute and sum how many nonzeros elements there are in each expanded row, so we might find halfway that we should have started with int64 from the beginning, which is not ideal... I think this method is a good heuristic that prevents overflow errors while still allowing people to use int32 in most cases where they can.

Let me know what you think !

Tests are green now so it's encouraging, I'll just need to add a what's new entry, and maybe I was thinking add a bit more cases to the test to be sure the code works in every case

My initial thought was to start the computation and when it discovers that it should be int64 to begin with, create a int64 container, copy over the results, and continue from there.

Give that, I see that the heuristic we have here is simple. My concern is that there would be a memory regression for input that did not actual need the int64 conversion.

Thanks a lot for the review @thomasjpfan ! And sorry for the delay in my answer

Regarding the strategy to start the computation, I agree it would be best, but I couldn't find any way to do that efficiently: the problem I have is that I don't know how it could find on the fly that we should have used int64 to begin with: indeed it seems cython code, when doing operations that will overflow, doesn't raise any warning/error (see this notebook for an example https://colab.research.google.com/drive/18KW2Zk1eq5d-WkTBf_l23UF2zNs0USpk?usp=sharing), contrary to numpy for instance, so we can't catch such error to know we should have done an int64 conversion.

So I guess the only way to know if say c=a+b will overflow (with a and b both np.int32), is to convert them to np.int64, and compare c with the max possible value for int32 no ? So I guess this would slow down the computations in some cases, whereas with the heuristic above, in many cases the computations can be done without any conversions ?

But I agree that it could lead to memory regressions for input that did not actual need the int64 conversion.

Maybe we could do the following, by mixing the two strategies, what do you think ? :

• if the heuristic above finds that we can keep the data as int32, then we do so, this way the computation is fast (no conversions and comparisons with the max np.int32 at each operation), and keeps a low memory
• if the heuristic above finds that the expanded dimensionality is too big, then we do checks on the fly as you said i.e. for each operation we do it with int64 converted indices, check whether the result can hold on int32, if yes we return the int32, if not we copy the previous results in int64 and keep going with int64 computations for the next results as you said ?

Let me know how that sounds, if it's a good idea I'll start implementing it

(Btw, I have limited time these days, so if you or anyone feels like taking over the PR, go ahead, it could go faster, otherwise I'll try to continue from time to time when I have some time, sorry for the delay)

@wdevazelhes I have labeled as 'help wanted' if this is ok with you. If you find some time to continue the PR you can always... help yourself... :)

@wdevazelhes I have labeled as 'help wanted' if this is ok with you. If you find some time to continue the PR you can always... help yourself... :)

@cmarmo thanks a lot, that works :)

Hello, scikit-learn community. I have some time and the great interest to make contribution to the project. I investigated history of this issue and read current PR and now I have a working solution for this problem. Could I take this issue to work on it and to create a new PR?

I slightly rewrote code of _csr_polynomial_expansion.pyx to move logic for array allocation outside of the c-code.

Closing as superseded by #23731. Thanks @wdevazelhes for your work!