I’m not quite sure what you mean, but I think part of what your getting at is that such “huge” problems do indeed tend to have “huge” reductions during pre-solving (if not, sometimes you’re just screwed, even a MILP is NP-hard after all). However, it’s really hard to do really good pre-solving before getting it into standard form. Getting it into standard form can be surprisingly hard, which certainly seems counter-intuitive for a purely linear problem. Consider that, if it was convenient for you to express your variable as a matrix rather than a vector, the most general form of a constraint is AXB \le C. However, each element of AXB is an entire array (since it is in fact a symbolic expression in terms of the elements of X). For small enough matrices, this is fine, you can think of it as a rank-3 array and in principle you can even stack-allocate this (although in practice you might not find that you have to do some work to even get here, since there isn’t exactly an abundance of libraries that does this sort of thing for you). However, as the problems get very large (thousands or millions of rows and columns in the matrices), it becomes obligatory to take advantage of the sparsity of the matrices (if they are mostly dense, probably the problem is just intractable in the first place). This is the sort of thing that’s beautifully easy in Julia and a nightmare in Python. Even if numpy provides you with a type of sparse array that’s appropriate here it’s not going to do you much good because it has to interoperate with the (very slow, python-implemented) symbolic variables that are the elements of at least one of these matrices. numpy pretty much falls back to Python code when it hits something it doesn’t already have library functions for written in C. In Julia all of this is quite nice, JuMP variables and expressions work just fine with whatever type of array you can dream up. You might have to define your own array operations to get it really efficient, but this only winds up being a few lines of code ususally, and you can easily parallelize those operations.
Working on such a large problem with gigantic sparse matrices might seem like a very esoteric use case, but those of us who use commercial solvers in industry often find ourselves dealing with such large problems. Yes, I might be able to find ways to re-formulate a lot of those so that the variable is no longer an unwieldy matrix, but that’s a big pain in the ass, and it’s very hard to do that in a generic way for every problem, so this is really the sort of thing I want my solver interface to do for me.
.