MemoryError when fitting on sparse X as apparently a Hessian matrix is being instantied ?

requirement: pip install libsvmdata -- python utility to download data from LIBSVM website; the first time downloading the data may take 2 mins.

The following script causes a MemoryError on my machine:

from libsvmdata import fetch_libsvm
from glum import GeneralizedLinearRegressor

X, y = fetch_libsvm("finance")
clf = GeneralizedLinearRegressor(family='gaussian', l1_ratio=1, alpha=1).fit(
    X, y)

output:

---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
~/workspace/mem_error.py in <module>
      3 
      4 X, y = fetch_libsvm("finance")
----> 5 clf = GeneralizedLinearRegressor(family='gaussian', l1_ratio=1, alpha=1).fit(
      6     X, y)

~/miniconda3/lib/python3.8/site-packages/glum/_glm.py in fit(self, X, y, sample_weight, offset, weights_sum)
   2448                     )
   2449                 )
-> 2450             coef = self._solve(
   2451                 X=X,
   2452                 y=y,

~/miniconda3/lib/python3.8/site-packages/glum/_glm.py in _solve(self, X, y, sample_weight, P2, P1, coef, offset, lower_bounds, upper_bounds, A_ineq, b_ineq)
    976             # 4.2 coordinate descent ##############################################
    977             elif self._solver == "irls-cd":
--> 978                 coef, self.n_iter_, self._n_cycles, self.diagnostics_ = _irls_solver(
    979                     _cd_solver, coef, irls_data
    980                 )

~/miniconda3/lib/python3.8/site-packages/glum/_solvers.py in _irls_solver(inner_solver, coef, data)
    287     https://www.csie.ntu.edu.tw/~cjlin/papers/l1_glmnet/long-glmnet.pdf
    288     """
--> 289     state = IRLSState(coef, data)
    290 
    291     state.eta, state.mu, state.obj_val, coef_P2 = _update_predictions(

~/miniconda3/lib/python3.8/site-packages/glum/_solvers.py in __init__(self, coef, data)
    529         self.gradient_rows = None
    530         self.hessian_rows = None
--> 531         self.hessian = np.zeros(
    532             (self.coef.shape[0], self.coef.shape[0]), dtype=data.X.dtype
    533         )

MemoryError: Unable to allocate 133. TiB for an array with shape (4272228, 4272228) and data type float64

ping @qb3

Reproduce the failure with the provided libsvmdata script, then inspect glum/_solvers.py at IRLSState.init and the fit/_solve path shown in glum/_glm.py. Trace how the sparse X shape reaches the Hessian allocation and verify that fitting the finance data no longer attempts the impossible allocation or raises MemoryError.