[Medium] Fast solvers for situations with many columns.

Currently, the solvers all build a Fisher matrix/Hessian. If the data matrix is n x m, this Fisher matrix will have shape m x m. If m is less than a few thousand, this is fine. However, if m gets very large, possibly due to lots of categorical variables, the current solvers will be inadequate because they are constructing a very large dense matrix.

Currently, we can avoid that matrix construction using the diag_fisher=True option. However, this is not particularly well supported and performance might not be very good. We should consider adding a benchmark with more than 10,000 columns to better track performance in these cases.

As part of this effort, it would be nice to separate the diag_fisher behavior and simply use a different solver function for high dimensional problems.

Start by reviewing the current solver functions and the diag_fisher=True path described in the issue. Measure performance and memory use on data with more than 10,000 columns, then determine how a separate high-dimensional solver should be exposed and benchmarked. Done means high-column problems avoid constructing a dense m x m Fisher matrix and have performance coverage.