ScikitLearn regression: Design matrix X too big for regression. What do I do?

It is a well known issue that can be address using out-of-core learning. By googling the term you will find several ways to address the problem.

For your specific problem, you have first to create a generator that will yield a row (or several of them) of your matrix and than using the partial_fit method of your algorithm.

Standard algorithms of scikit-learn use actually an exact computation of the solution, like sklearn.linear_model.LinearRegression or sklearn.linear_model.LinearRegression.RidgeCV. Other methods are based on batch learning and have a partial_fit methods like sklearn.linear_model.SGDRegressor, allowing to fit only a mini-batch. It is what you are looking for.

The process is: use the generator to yield a mini-batch, apply the partial_fit method, delete the mini-batch from the memory and get a new one.

However, as this method is stochastic and depends of the order of your data and your initialization of the weights, at the opposite of the solution given by the standard regression methods that can fit all the data in the memory. I won't enter into the details but look at gradient descent optimization to understand how it works (http://ruder.io/optimizing-gradient-descent/)

That's not really very large-scale in my opinion and you probably do not need to use out-of-core learning here (although i don't know how much memory you got).

Using out-of-core approaches when not needed, you will pay a price (not robust without heavy tuning)!

It would be also interesting to know if your features are sparse or dense which can make a huge difference (as most solvers can exploit sparse-data!).

Some things to say here:

• RidgeCV uses some problem-tuned cross-validation approach and probably because of this reason, there is no parameter for controlling the underlying solver
  • It's not uncommon to use a problem-tuned CV-method when efficient warm-starting can help the process (performance)
• When doing CV manually, all the tools are available in sklearn, you can chose different solvers
  • Those are very different in terms of the approach and characteristics

solver : {‘auto’, ‘svd’, ‘cholesky’, ‘lsqr’, ‘sparse_cg’, ‘sag’, ‘saga’}

Solver to use in the computational routines:

‘auto’ chooses the solver automatically based on the type of data.

‘svd’ uses a Singular Value Decomposition of X to compute the Ridge coefficients. More stable for singular matrices than ‘cholesky’.

‘cholesky’ uses the standard scipy.linalg.solve function to obtain a closed-form solution.

‘sparse_cg’ uses the conjugate gradient solver as found in scipy.sparse.linalg.cg. As an iterative algorithm, this solver is more appropriate than ‘cholesky’ for large-scale data (possibility to set tol and max_iter).

‘lsqr’ uses the dedicated regularized least-squares routine scipy.sparse.linalg.lsqr. It is the fastest but may not be available in old scipy versions. It also uses an iterative procedure.

‘sag’ uses a Stochastic Average Gradient descent, and ‘saga’ uses its improved, unbiased version named SAGA. Both methods also use an iterative procedure, and are often faster than other solvers when both n_samples and n_features are large. Note that ‘sag’ and ‘saga’ fast convergence is only guaranteed on features with approximately the same scale. You can preprocess the data with a scaler from sklearn.preprocessing.

All last five solvers support both dense and sparse data. However, only ‘sag’ and ‘saga’ supports sparse input when fit_intercept is True.

So i highly recommend to try: sparse_cg or maybe lsqr together with manual CV. If that's working (does for me perfectly), which i would expect, this approach is a lot more stable / robust (compared to out-of-core approaches using SGD) and you don't need to tune their parameters much, which is a huge advantage.

Of course, one can always use sag and sgd, but convergence-theory is based on some strong assumptions in regards to parameter tuning. In very large-scale settings, these two are viable candidates (as others won't work), but here i don't see much merrit (again: i'm not sure how much memory you got). If the above approach is not working, try saga before sgd and respect the rules (standardization + param-tuning). (edit: sag is pretty bad for my test-data!)

Example:

from sklearn.datasets import make_regression
from sklearn.linear_model import Ridge
from time import perf_counter

X, y, = make_regression(n_samples=38000, n_features=7000, n_informative=500,
                        bias=-2.0, noise=0.1, random_state=0)

print(type(X))  # dense!!! if your data is sparse; use that fact!

clf = Ridge(alpha=1.0, solver="lsqr")
start = perf_counter()
clf.fit(X, y)
end = perf_counter()
print('LSQR: used secs: ', end-start)

Output:

LSQR: used secs:  8.489622474064486

So even in the dense-case, that's not that hard to optimize (and uses ~6-8 GB of memory).

And while i would be careful to implicate equivalence between the Ridge-model and the following SGD-based Ridge-model (careful which variable is part of regularization; too lazy to check), here just a demo on how hard it is to tune SGD. Take this with a grain of salt (and maybe don't evaluate the absolute score; but the variance depending on parameters):

Remark: that's a smaller example! Using your original example, none of the SGD-approaches will obtain convergence without manually setting down the learning-rate eta_0 (as internal heuristics are not able to do that for you!)

Partial-Code:

X, y, = make_regression(n_samples=3800, n_features=700, n_informative=500,
                        noise=0.1, random_state=0)

print(type(X))  # dense!!! if your data is sparse; use that fact!

clf = Ridge(alpha=1.0, solver="lsqr", fit_intercept=False)
start = perf_counter()
clf.fit(X, y)
end = perf_counter()
print('LSQR: used secs: ', end-start)
print('train-score: ', clf.score(X, y))

clf = Ridge(alpha=1.0, solver="sparse_cg", fit_intercept=False)
start = perf_counter()
clf.fit(X, y)
end = perf_counter()
print('sparse_cg: used secs: ', end-start)
print('train-score: ', clf.score(X, y))

clf = SGDRegressor(loss='squared_loss', penalty='l2', alpha=1., fit_intercept=False,
                   random_state=0)
start = perf_counter()
clf.fit(X, y)
end = perf_counter()
print('SGD: used secs: ', end-start)
print('train-score: ', clf.score(X, y))

clf = SGDRegressor(loss='squared_loss', penalty='l2', alpha=1., fit_intercept=False,
                   random_state=0, average=True)
start = perf_counter()
clf.fit(X, y)
end = perf_counter()
print('SGD: used secs: ', end-start)
print('train-score: ', clf.score(X, y))

clf = SGDRegressor(loss='squared_loss', penalty='l2', alpha=1., fit_intercept=False,
                   random_state=0, learning_rate="constant", eta0=0.001)
start = perf_counter()
clf.fit(X, y)
end = perf_counter()
print('SGD: used secs: ', end-start)
print('train-score: ', clf.score(X, y))

clf = SGDRegressor(loss='squared_loss', penalty='l2', alpha=1., fit_intercept=False,
                   random_state=0, n_iter=50, average=True)
start = perf_counter()
clf.fit(X, y)
end = perf_counter()
print('SGD: used secs: ', end-start)
print('train-score: ', clf.score(X, y))

Output:

LSQR: used secs:  0.08252486090450709
train-score:  0.999999907282
sparse_cg: used secs:  0.13668818702548152
train-score:  0.999999181151
SGD: used secs:  0.04154542095705427
train-score:  0.743448766459
SGD: used secs:  0.05300238587407993
train-score:  0.774611911034
SGD: used secs:  0.038653031605587
train-score:  0.733585661919
SGD: used secs:  0.46313909066321507
train-score:  0.776444474871

Also: when going for the partial_fit / out-of-memory approach, tuning the size of mini_batches also needs to be done (which is ignored in the above demo -> pure SGD -> one sample at a time)! And again: this is not easy!