How does $\chi^2$ feature selection work? I can't find the information how $\chi^2$ are used to select numerical features for a model.

Fro instance, If I employ the sklearn library:

from sklearn.datasets import load_iris
from sklearn.feature_selection import chi2

iris = load_iris()
X, y = iris.data, iris.target # X contains 4 features, y does 3 classes


I can build the following table for the dataset

pivot_table = np.zeros((X.shape[1], len(np.unique(y))))
pivot_table[:, 0] = X[y==0].sum(axis=0)
pivot_table[:, 1] = X[y==1].sum(axis=0)
pivot_table[:, 2] = X[y==2].sum(axis=0)

print(pivot_table.T)
array([[250.3, 171.4,  73.1,  12.3],
       [296.8, 138.5, 213. ,  66.3],
       [329.4, 148.7, 277.6, 101.3]])


Then applying  the contingency table to it:

from scipy.stats import chi2_contingency
chi2_contingency(pivot_table, correction=False)


Output:

125.58397740261773
1.0924329599765022e-24
6
[[213.82265358 301.31664021 361.36070621]
 [111.8757204  157.6540915  189.0701881 ]
 [137.51492279 193.78458652 232.40049069]
 [ 43.88670323  61.84468177  74.168615  ]]


But this only provides the information about dependence of variables. And they are dependent because p-value is quite low.

Then I use the other function and I get what I want:

chi2(X,y)
(array([ 10.81782088,   3.7107283 , 116.31261309,  67.0483602 ]),
 array([4.47651499e-03, 1.56395980e-01, 5.53397228e-26, 2.75824965e-15]))


But how does this work? What math lies at the heart of it?
Just go inside the chi2 function

It calculates chi2 values for Y's based on the contingency table then computes p-values for each feature