Instead of binarising discrete variables for NN algorithms, could I attempt to assign them a value by splitting into buckets and sorting? Creating a boolean for each value of each discrete variable explodes the number of features. I have an idea for solving it but I wasn't able to find anything on this on the internet, so perhaps someone could either point me to a resource or let me know whether they think this is possible. And if not, maybe point to another solution (all I could find was binarisation).

So if we have training vectors $X_i$ and similarly labels $y_i$, I propose we partition the training set by some variable $j$, so that for each $X_i$ for which $j_i$ equals some value $q$, $X_i$ will be in $P(q)$.

Then for each $q$ we calculate the mean $y_i$ for $P(q)$ and sort those averages, mapping each $q$ onto its index in the sorted list. This value is used as the feature for $j$. I suspect that this ought to be better than just giving each value a random index, as any remaining non-linearity will be hopefully fixed by the hidden layers. But it seems that these non-linear patterns will be easier to pick up if the function is somewhat quasi-monotonous.

I think this could be improved, potentially instead of just finding the mean, another measure could work better and also for multiple discrete features, I am not sure if this process should be repeated for each variable separately or through a multi-dimensional partitioning.

The main question is of course, would this process be helpful at all, and if not, is there anything similar that could solve this problem?
Just taking the mean seems better, that way the distances between the groups is better represented than just indexing them by order. With the target you should be a bit careful not to include the current row for the mean, since for the test set and the actual predictions you cannot use that target either (for obvious reasons). To get even more information about the category you try to get rid of is to also take the variance or standard deviation of the target within this group and take the mean and variance of important numerical features for these categories. If you have multiple categories you can also take the expected value for features and targets conditioned on multiple categories at once. This is a nice compromise between keeping enough information on high cardinality while keeping the dimensionality low.

Here's an example where the c is the categorical variable, n are numeric and y is your target. We want to get rid of c (because c could have 1000s of values).

xc_1 | xn_1 | xn_2 | y
A    | 2    | 4    | 1
A    | 1    | 2    | 1
A    | 4    | 3    | 0
B    | 3    | 5    | 0
B    | 5    | 4    | 1
B    | 6    | 6    | 0


To keep things a little more simple we will disregard the values of the current row, which is better anyway. For the first row we will look at all other A rows, take the averages (and maybe other statistics) of the numeric features and of y, we will do the same for B, then we get:

xc_1 | xn_1 | xn_2 | xnc_1 | xnc_2 | yc_1 | y
A    | 2    | 4    | 2.5   | 2.5   | 0.5  | 1
A    | 1    | 2    | 3     | 3.5   | 0.5  | 1
A    | 4    | 3    | 1.5   | 3     | 1    | 0 
B    | 3    | 5    | 5.5   | 5     | 0.5  | 0 
B    | 5    | 4    | 4.5   | 5.5   | 0    | 1  
B    | 6    | 6    | 4     | 4.5   | 0.5  | 0


We have captured a lot of information about category 1, which we can now discard.