From 64a7df455279393d9b7fd14533c84c4c0b495399 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Thu, 12 May 2022 18:54:31 +0200
Subject: [PATCH 01/17] fit_intercept for cd solver (using sklearn
 _preprocess_data)

---
 solvers/cd.py | 19 ++++++++++++++++---
 1 file changed, 16 insertions(+), 3 deletions(-)

diff --git a/solvers/cd.py b/solvers/cd.py
index 289a2170..ef5a5274 100644
--- a/solvers/cd.py
+++ b/solvers/cd.py
@@ -5,6 +5,7 @@
     import numpy as np
     from scipy import sparse
     from numba import njit
+    from sklearn.linear_model._base import _preprocess_data
 
 
 if import_ctx.failed_import:
@@ -37,14 +38,22 @@ class Solver(BaseSolver):
     ]
 
     def skip(self, X, y, lmbd, fit_intercept):
-        # XXX - not implemented but this should be quite easy
-        if fit_intercept:
+        # XXX - intercept not implemented for sparse X but it shouldn't be hard
+        if fit_intercept and sparse.issparse(X):
             return True, f"{self.name} does not handle fit_intercept"
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
-        self.y, self.lmbd = y, lmbd
+        # sklearn way of handling intercept: center y and X for dense data
+        if fit_intercept:
+            X, y, X_offset, y_offset, _ = _preprocess_data(
+                X, y, fit_intercept, return_mean=True, copy=True,
+            )
+            self.X_offset = X_offset
+            self.y_offset = y_offset
+
+        self.y, self.lmbd, self.fit_intercept = y, lmbd, fit_intercept
 
         if sparse.issparse(X):
             self.X = X
@@ -66,6 +75,10 @@ def run(self, n_iter):
             L = (self.X ** 2).sum(axis=0)
             self.w = self.cd(self.X, self.y, self.lmbd, L, n_iter)
 
+        if self.fit_intercept:
+            intercept = self.y_offset - self.X_offset @ self.w
+            self.w = np.r_[self.w, intercept]
+
     @staticmethod
     @njit
     def cd(X, y, lmbd, L, n_iter):

From de13e7dda484c5f51ab4642fbe19dc3c77a6a234 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Thu, 12 May 2022 21:03:36 +0200
Subject: [PATCH 02/17] fit_intercept support for python_pgd solver

---
 solvers/python_pgd.py | 27 ++++++++++++++++++++++-----
 1 file changed, 22 insertions(+), 5 deletions(-)

diff --git a/solvers/python_pgd.py b/solvers/python_pgd.py
index 81f600b2..68668bb8 100644
--- a/solvers/python_pgd.py
+++ b/solvers/python_pgd.py
@@ -3,6 +3,7 @@
 with safe_import_context() as import_ctx:
     import numpy as np
     from scipy import sparse
+    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
@@ -22,13 +23,24 @@ class Solver(BaseSolver):
     ]
 
     def skip(self, X, y, lmbd, fit_intercept):
-        # XXX - not implemented but not too complicated to implement
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
+        # XXX - intercept not implemented for sparse X but it shouldn't be hard
+        if fit_intercept and sparse.issparse(X):
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data",
+            )
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
+        # sklearn way of handling intercept: center y and X for dense data
+        if fit_intercept:
+            X, y, X_offset, y_offset, _ = _preprocess_data(
+                X, y, fit_intercept, return_mean=True, copy=True,
+            )
+            self.X_offset = X_offset
+            self.y_offset = y_offset
+
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
 
@@ -40,8 +52,10 @@ def run(self, callback):
         if self.use_acceleration:
             z = np.zeros(n_features)
 
+        intercept = self.y_offset if self.fit_intercept else []
+
         t_new = 1
-        while callback(w):
+        while callback(np.r_[w, intercept]):
             if self.use_acceleration:
                 t_old = t_new
                 t_new = (1 + np.sqrt(1 + 4 * t_old ** 2)) / 2
@@ -53,7 +67,10 @@ def run(self, callback):
                 w -= self.X.T @ (self.X @ w - self.y) / L
                 w = self.st(w, self.lmbd / L)
 
-        self.w = w
+            if self.fit_intercept:
+                intercept = self.y_offset - self.X_offset @ w
+
+        self.w = np.r_[w, intercept]
 
     def st(self, w, mu):
         w -= np.clip(w, -mu, mu)

From a003c8abf140dfeea5e3c3103db80b41e46a4781 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Thu, 12 May 2022 21:06:00 +0200
Subject: [PATCH 03/17] cosmit

---
 solvers/cd.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/solvers/cd.py b/solvers/cd.py
index ef5a5274..11d449a8 100644
--- a/solvers/cd.py
+++ b/solvers/cd.py
@@ -40,7 +40,10 @@ class Solver(BaseSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX - intercept not implemented for sparse X but it shouldn't be hard
         if fit_intercept and sparse.issparse(X):
-            return True, f"{self.name} does not handle fit_intercept"
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data",
+            )
 
         return False, None
 

From 169faa390b86d9150205007105136a227cf4f9c8 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Fri, 13 May 2022 13:03:04 +0200
Subject: [PATCH 04/17] Explicitly excluding sparse case on fit intercept

---
 solvers/cd.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/solvers/cd.py b/solvers/cd.py
index 11d449a8..78765e5f 100644
--- a/solvers/cd.py
+++ b/solvers/cd.py
@@ -78,7 +78,7 @@ def run(self, n_iter):
             L = (self.X ** 2).sum(axis=0)
             self.w = self.cd(self.X, self.y, self.lmbd, L, n_iter)
 
-        if self.fit_intercept:
+        if self.fit_intercept and not sparse.issparse(self.X):
             intercept = self.y_offset - self.X_offset @ self.w
             self.w = np.r_[self.w, intercept]
 

From 4a6fedabae3d37924df05dc639155d344fdb4861 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Fri, 13 May 2022 14:08:31 +0200
Subject: [PATCH 05/17] fit_intercept for blitz solver (supports sparse data!)

---
 solvers/blitz.py | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/solvers/blitz.py b/solvers/blitz.py
index 09dfdecc..683bab5f 100644
--- a/solvers/blitz.py
+++ b/solvers/blitz.py
@@ -4,6 +4,7 @@
 
 with safe_import_context() as import_ctx:
     import blitzl1
+    import numpy as np
 
 
 class Solver(BaseSolver):
@@ -20,16 +21,11 @@ class Solver(BaseSolver):
         'vol. 37, pp. 1171-1179 (2015)'
     ]
 
-    def skip(self, X, y, lmbd, fit_intercept):
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
-
-        return False, None
-
     def set_objective(self, X, y, lmbd, fit_intercept):
         self.X, self.y, self.lmbd = X, y, lmbd
+        self.fit_intercept = fit_intercept
 
-        blitzl1.set_use_intercept(False)
+        blitzl1.set_use_intercept(self.fit_intercept)
         blitzl1.set_tolerance(0)
         self.problem = blitzl1.LassoProblem(self.X, self.y)
 
@@ -39,7 +35,10 @@ def get_next(previous):
         return previous + 1
 
     def run(self, n_iter):
-        self.coef_ = self.problem.solve(self.lmbd, max_iter=n_iter).x
+        self.sol_ = self.problem.solve(self.lmbd, max_iter=n_iter)
 
     def get_result(self):
-        return self.coef_.flatten()
+        if self.fit_intercept:
+            return np.r_[self.sol_.x, self.sol_.intercept]
+        else:
+            return self.sol_.x.flatten()

From eb98950f34592ec69da4ca6628a2889100d53880 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Fri, 13 May 2022 17:05:41 +0200
Subject: [PATCH 06/17] fit_intercept for skglm solver (only dense data)

---
 solvers/skglm.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/solvers/skglm.py b/solvers/skglm.py
index 05091d4d..62787999 100644
--- a/solvers/skglm.py
+++ b/solvers/skglm.py
@@ -5,6 +5,7 @@
 
 with safe_import_context() as import_ctx:
     import numpy as np
+    from scipy import sparse
     from skglm import Lasso
     from sklearn.exceptions import ConvergenceWarning
 
@@ -25,8 +26,11 @@ class Solver(BaseSolver):
     ]
 
     def skip(self, X, y, lmbd, fit_intercept):
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
+        if fit_intercept and sparse.issparse(X):
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data",
+            )
 
         return False, None
 
@@ -38,7 +42,8 @@ def set_objective(self, X, y, lmbd, fit_intercept):
         n_samples = self.X.shape[0]
         self.lasso = Lasso(
             alpha=self.lmbd / n_samples, max_iter=1, max_epochs=50_000,
-            tol=1e-12, fit_intercept=False, warm_start=False, verbose=False)
+            tol=1e-12, fit_intercept=self.fit_intercept, warm_start=False,
+            verbose=False)
 
         # Cache Numba compilation
         self.run(1)

From eeb0143878dc665593a1b31952f2f764e307c2b5 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Tue, 17 May 2022 13:55:58 +0200
Subject: [PATCH 07/17] fit_intercept for lightning solver (only for dense
 data)

---
 solvers/lightning.py | 25 ++++++++++++++++++++-----
 1 file changed, 20 insertions(+), 5 deletions(-)

diff --git a/solvers/lightning.py b/solvers/lightning.py
index 7c9becdd..3c1ef76b 100644
--- a/solvers/lightning.py
+++ b/solvers/lightning.py
@@ -4,11 +4,11 @@
 
 with safe_import_context() as import_ctx:
     import numpy as np
+    from scipy import sparse
     from lightning.regression import CDRegressor
+    from sklearn.linear_model._base import _preprocess_data    
 
 
-# TODO: lightning always fit an intercept
-#       it is thus not optimizing the same cost function
 class Solver(BaseSolver):
     name = 'Lightning'
 
@@ -25,12 +25,25 @@ class Solver(BaseSolver):
     ]
 
     def skip(self, X, y, lmbd, fit_intercept):
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
+        if fit_intercept and sparse.issparse(X):
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data",
+            )
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
+        # lightning has an attribut intercept_ but it is not handled properly
+        # (as it is simply set to zero). For this reason, we use the sklearn
+        # way of handling intercept: center y and X beforehand for dense data
+        if fit_intercept:
+            X, y, X_offset, y_offset, _ = _preprocess_data(
+                X, y, fit_intercept, return_mean=True, copy=True,
+            )
+            self.X_offset = X_offset
+            self.y_offset = y_offset
+
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
 
@@ -45,6 +58,8 @@ def run(self, n_iter):
 
     def get_result(self):
         beta = self.clf.coef_.flatten()
+
         if self.fit_intercept:
-            beta = np.r_[beta, self.clf.intercept_]
+            intercept = self.y_offset - self.X_offset @ beta
+            beta = np.r_[beta, intercept]
         return beta

From d2849629caa1b03166c8ba0c5be38ff38dc1b73c Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Tue, 17 May 2022 17:20:16 +0200
Subject: [PATCH 08/17] fit_intercept for noncvx_pro solver (only for dense
 data for now)

---
 solvers/noncvx_pro.py | 25 ++++++++++++++++++++++---
 1 file changed, 22 insertions(+), 3 deletions(-)

diff --git a/solvers/noncvx_pro.py b/solvers/noncvx_pro.py
index 8aa4ece6..cb593e26 100644
--- a/solvers/noncvx_pro.py
+++ b/solvers/noncvx_pro.py
@@ -6,6 +6,7 @@
     from numpy.linalg import norm
     import scipy.optimize as sciop
     from scipy.sparse import issparse
+    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
@@ -14,14 +15,28 @@ class Solver(BaseSolver):
     stopping_strategy = 'iteration'
 
     def set_objective(self, X, y, lmbd, fit_intercept):
+        # sklearn way of handling intercept: center y and X for dense data
+        # when X is sparse, X_offset is computed but X is not centered
+        if fit_intercept:
+            X, y, X_offset, y_offset, _ = _preprocess_data(
+                X, y, fit_intercept, return_mean=True, copy=True,
+            )
+            self.X_offset = X_offset
+            self.y_offset = y_offset
+
         self.X, self.y, self.lmbd = X, y, lmbd
+        self.fit_intercept = fit_intercept
 
     def skip(self, X, y, lmbd, fit_intercept):
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
         # XXX: make this solver work with sparse matrices.
         if issparse(X):
             return True, f"{self.name} does not support sparse design matrices"
+        # XXX: even if sparse support is added, the test below should be kept
+        # unless fit_intercept is properly handled for sparse matrices
+        # (by manually considering X_offset in calculations)
+        if fit_intercept and issparse(X):
+            return True, \
+                f"{self.name} doesn't handle fit_intercept with sparse matrix",
         return False, None
 
     def run(self, n_iter):
@@ -60,7 +75,7 @@ def nabla_f(v):
                 g = u * (Cx - Xty) / lmbd + v
                 return f, g
 
-        opts = {'gtol': 1e-8, 'maxiter': n_iter, 'maxcor': 100, 'ftol': 0}
+        opts = {'gtol': 1e-30, 'maxiter': n_iter, 'maxcor': 100, 'ftol': 1e-30}
         u0 = np.ones(n_features)
 
         lbfgs_res = sciop.minimize(
@@ -70,5 +85,9 @@ def nabla_f(v):
 
         self.w = v * u_opt(v)
 
+        if self.fit_intercept and not issparse(self.X):
+            intercept = self.y_offset - self.X_offset @ self.w
+            self.w = np.r_[self.w, intercept]
+
     def get_result(self):
         return self.w.flatten()

From 2b5bf12a9f2922906a97375ac382ce8d77c649d1 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Tue, 17 May 2022 17:53:41 +0200
Subject: [PATCH 09/17] Fix linting

---
 solvers/lightning.py  | 2 +-
 solvers/noncvx_pro.py | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/solvers/lightning.py b/solvers/lightning.py
index 3c1ef76b..4eb13b9e 100644
--- a/solvers/lightning.py
+++ b/solvers/lightning.py
@@ -6,7 +6,7 @@
     import numpy as np
     from scipy import sparse
     from lightning.regression import CDRegressor
-    from sklearn.linear_model._base import _preprocess_data    
+    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
diff --git a/solvers/noncvx_pro.py b/solvers/noncvx_pro.py
index 143e37cd..15193d52 100644
--- a/solvers/noncvx_pro.py
+++ b/solvers/noncvx_pro.py
@@ -21,7 +21,7 @@ class Solver(BaseSolver):
 
     def set_objective(self, X, y, lmbd, fit_intercept):
         # sklearn way of handling intercept: center y and X for dense data
-        # when X is sparse, X_offset is computed but X is not centered
+        # when X is sparse, X_offset is computed but X is not centered
         if fit_intercept:
             X, y, X_offset, y_offset, _ = _preprocess_data(
                 X, y, fit_intercept, return_mean=True, copy=True,
@@ -31,11 +31,11 @@ def set_objective(self, X, y, lmbd, fit_intercept):
 
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
-        
+
         if X.shape[0] >= X.shape[1]:
             self.C = X.T @ X
             if issparse(self.C):
-                self.C = self.C.toarray()        
+                self.C = self.C.toarray()
 
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX: make this solver work with sparse matrices.

From 3a365959b6e05d0b3f1aceaf9cab24cfee351e5a Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Wed, 18 May 2022 15:45:07 +0200
Subject: [PATCH 10/17] fit_intercept support on solver L-BFGS-B (only dense
 data for now)

---
 solvers/l_bfgs_b.py | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/solvers/l_bfgs_b.py b/solvers/l_bfgs_b.py
index 2410c352..187a0b76 100644
--- a/solvers/l_bfgs_b.py
+++ b/solvers/l_bfgs_b.py
@@ -5,6 +5,8 @@
     import numpy as np
     from numpy.linalg import norm
     from scipy.optimize import fmin_l_bfgs_b
+    from scipy.sparse import issparse
+    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
@@ -28,13 +30,22 @@ class Solver(BaseSolver):
     ]
 
     def skip(self, X, y, lmbd, fit_intercept):
-        # XXX - not implemented but this should be quite easy
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
+        # XXX - intercept not implemented for sparse X for now
+        if fit_intercept and issparse(X):
+            return True, \
+                f"{self.name} doesn't handle fit_intercept with sparse data",
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
+        # sklearn way of handling intercept: center y and X for dense data
+        if fit_intercept:
+            X, y, X_offset, y_offset, _ = _preprocess_data(
+                X, y, fit_intercept, return_mean=True, copy=True,
+            )
+            self.X_offset = X_offset
+            self.y_offset = y_offset
+
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
 
@@ -65,5 +76,9 @@ def gradf(w):
 
         self.w = w_hat
 
+        if self.fit_intercept and not issparse(self.X):
+            intercept = self.y_offset - self.X_offset @ self.w
+            self.w = np.r_[self.w, intercept]
+
     def get_result(self):
         return self.w.flatten()

From a8cb44ec520bd258a5463fa19643fcea1e71d9b8 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Wed, 18 May 2022 17:12:19 +0200
Subject: [PATCH 11/17] fit_intercept support for solver Julia-PGD (only dense
 for now)

---
 solvers/julia_pgd.py | 23 ++++++++++++++++++++---
 1 file changed, 20 insertions(+), 3 deletions(-)

diff --git a/solvers/julia_pgd.py b/solvers/julia_pgd.py
index 999a419f..322a2abb 100644
--- a/solvers/julia_pgd.py
+++ b/solvers/julia_pgd.py
@@ -6,6 +6,9 @@
 from benchopt.helpers.julia import assert_julia_installed
 
 with safe_import_context() as import_ctx:
+    import numpy as np
+    from scipy.sparse import issparse
+    from sklearn.linear_model._base import _preprocess_data
     assert_julia_installed()
 
 
@@ -27,15 +30,25 @@ class Solver(JuliaSolver):
         'algorithm for linear inverse problems", SIAM J. Imaging Sci., '
         'vol. 2, no. 1, pp. 183-202 (2009)'
     ]
+    support_sparse = False
 
     def skip(self, X, y, lmbd, fit_intercept):
-        # XXX - fit intercept is not yet implemented in julia.jl
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
+        # XXX - fit intercept is not yet implemented in julia.jl for sparse X
+        if fit_intercept and issparse(X):
+            return True, \
+                f"{self.name} doesn't handle fit_intercept with sparse data",
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
+        # sklearn way of handling intercept: center y and X for dense data
+        if fit_intercept:
+            X, y, X_offset, y_offset, _ = _preprocess_data(
+                X, y, fit_intercept, return_mean=True, copy=True,
+            )
+            self.X_offset = X_offset
+            self.y_offset = y_offset
+
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
 
@@ -45,5 +58,9 @@ def set_objective(self, X, y, lmbd, fit_intercept):
     def run(self, n_iter):
         self.beta = self.solve_lasso(self.X, self.y, self.lmbd, n_iter)
 
+        if self.fit_intercept and not issparse(self.X):
+            intercept = self.y_offset - self.X_offset @ self.beta
+            self.beta = np.r_[self.beta.ravel(), intercept]
+
     def get_result(self):
         return self.beta.ravel()

From 2394eb266d9bbf5e00815dd888f2029c50a4c66f Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Thu, 19 May 2022 23:46:04 +0200
Subject: [PATCH 12/17] fit_intercept support for r-pgd solver (only dense data
 for now)

---
 solvers/r_pgd.py | 20 +++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/solvers/r_pgd.py b/solvers/r_pgd.py
index 54c65421..d1c1abf0 100644
--- a/solvers/r_pgd.py
+++ b/solvers/r_pgd.py
@@ -5,6 +5,8 @@
 
 with safe_import_context() as import_ctx:
     import numpy as np
+    from scipy.sparse import issparse
+    from sklearn.linear_model._base import _preprocess_data
 
     from rpy2 import robjects
     from rpy2.robjects import numpy2ri
@@ -34,12 +36,20 @@ class Solver(BaseSolver):
     ]
 
     def skip(self, X, y, lmbd, fit_intercept):
-        if fit_intercept:
-            return True, f"{self.name} does not handle fit_intercept"
+        # rpy2 does not directly support sparse matrices (workaround exists)
+        if fit_intercept and issparse(X):
+            return True, \
+                f"{self.name} doesn't handle fit_intercept with sparse data"
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
+        # sklearn way of handling intercept: center y and X (for dense data)
+        if fit_intercept:
+            X, y, self.X_offset, self.y_offset, _ = _preprocess_data(
+                X, y, fit_intercept, copy=True, return_mean=True
+            )
+
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
         self.r_pgd = robjects.r['proximal_gradient_descent']
@@ -52,4 +62,8 @@ def run(self, n_iter):
         self.w = np.array(as_r(coefs, "vector"))
 
     def get_result(self):
-        return self.w.flatten()
+        if self.fit_intercept:
+            intercept = self.y_offset - self.X_offset @ self.w
+            return np.r_[self.w.flatten(), intercept]
+        else:
+            return self.w.flatten()

From 75177bf9eb403e77701e6de425ea03d37e410303 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Tue, 31 May 2022 15:40:14 +0200
Subject: [PATCH 13/17] Adding sklearn requirements for fit_intercept support

---
 solvers/cd.py         | 2 +-
 solvers/julia_pgd.py  | 1 +
 solvers/lightning.py  | 1 +
 solvers/noncvx_pro.py | 1 +
 solvers/python_pgd.py | 1 +
 solvers/r_pgd.py      | 2 +-
 solvers/skglm.py      | 2 +-
 7 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/solvers/cd.py b/solvers/cd.py
index 78765e5f..62f589d1 100644
--- a/solvers/cd.py
+++ b/solvers/cd.py
@@ -27,7 +27,7 @@ class Solver(BaseSolver):
     name = "cd"
 
     install_cmd = 'conda'
-    requirements = ['numba']
+    requirements = ['numba', "scikit-learn"]
     references = [
         'W. J. Fu, "Penalized Regressions: the Bridge versus the Lasso", '
         'J. Comput. Graph. Statist., vol.7, no. 3, pp. 397-416, '
diff --git a/solvers/julia_pgd.py b/solvers/julia_pgd.py
index 322a2abb..84ef695d 100644
--- a/solvers/julia_pgd.py
+++ b/solvers/julia_pgd.py
@@ -21,6 +21,7 @@ class Solver(JuliaSolver):
     # Config of the solver
     name = 'Julia-PGD'
     stopping_strategy = 'iteration'
+    requirements = ["scikit-learn"]
     references = [
         'I. Daubechies, M. Defrise and C. De Mol, '
         '"An iterative thresholding algorithm for linear inverse problems '
diff --git a/solvers/lightning.py b/solvers/lightning.py
index 4eb13b9e..1ef51caf 100644
--- a/solvers/lightning.py
+++ b/solvers/lightning.py
@@ -15,6 +15,7 @@ class Solver(BaseSolver):
     install_cmd = 'conda'
     requirements = [
         'cython',
+        'scikit-learn',
         'pip:git+https://github.com/scikit-learn-contrib/lightning.git'
     ]
     references = [
diff --git a/solvers/noncvx_pro.py b/solvers/noncvx_pro.py
index b39e95ad..17066a2b 100644
--- a/solvers/noncvx_pro.py
+++ b/solvers/noncvx_pro.py
@@ -13,6 +13,7 @@ class Solver(BaseSolver):
     name = "noncvx-pro"
 
     stopping_strategy = 'iteration'
+    requirements = ["scikit-learn"]
     references = [
         "Clarice Poon and Gabriel Peyré, "
         "'Smooth Bilevel Programming for Sparse Regularization', "
diff --git a/solvers/python_pgd.py b/solvers/python_pgd.py
index 68668bb8..fb23160b 100644
--- a/solvers/python_pgd.py
+++ b/solvers/python_pgd.py
@@ -8,6 +8,7 @@
 
 class Solver(BaseSolver):
     name = 'Python-PGD'  # proximal gradient, optionally accelerated
+    requirements = ['scikit-learn']
     stopping_strategy = "callback"
 
     # any parameter defined here is accessible as a class attribute
diff --git a/solvers/r_pgd.py b/solvers/r_pgd.py
index d1c1abf0..8dd1a661 100644
--- a/solvers/r_pgd.py
+++ b/solvers/r_pgd.py
@@ -22,7 +22,7 @@ class Solver(BaseSolver):
     name = "R-PGD"
 
     install_cmd = 'conda'
-    requirements = ['r-base', 'rpy2']
+    requirements = ['r-base', 'rpy2', 'scikit-learn']
     stopping_strategy = 'iteration'
     support_sparse = False
     references = [
diff --git a/solvers/skglm.py b/solvers/skglm.py
index 5ef691e8..073b43bf 100644
--- a/solvers/skglm.py
+++ b/solvers/skglm.py
@@ -15,7 +15,7 @@ class Solver(BaseSolver):
 
     install_cmd = 'conda'
     requirements = [
-        'pip:skglm'
+        'pip:skglm', 'scikit-learn'
     ]
     references = [
         'Q. Bertrand and Q. Klopfenstein and P.-A. Bannier and G. Gidel'

From 2b634fcd7d297676fe9f00bf48360966ad813bef Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Wed, 20 Jul 2022 11:10:37 +0200
Subject: [PATCH 14/17] Removing sklearn dependency wherever possible. Handling
 fit_intercept manually instead.

---
 solvers/cd.py         | 16 +++++++---------
 solvers/julia_pgd.py  | 19 ++++++++-----------
 solvers/l_bfgs_b.py   | 18 ++++++++----------
 solvers/lightning.py  | 17 +++++++----------
 solvers/noncvx_pro.py | 20 ++++++++------------
 solvers/r_pgd.py      | 18 +++++++++---------
 6 files changed, 47 insertions(+), 61 deletions(-)

diff --git a/solvers/cd.py b/solvers/cd.py
index 62f589d1..6abc46a0 100644
--- a/solvers/cd.py
+++ b/solvers/cd.py
@@ -5,7 +5,6 @@
     import numpy as np
     from scipy import sparse
     from numba import njit
-    from sklearn.linear_model._base import _preprocess_data
 
 
 if import_ctx.failed_import:
@@ -27,7 +26,7 @@ class Solver(BaseSolver):
     name = "cd"
 
     install_cmd = 'conda'
-    requirements = ['numba', "scikit-learn"]
+    requirements = ['numba']
     references = [
         'W. J. Fu, "Penalized Regressions: the Bridge versus the Lasso", '
         'J. Comput. Graph. Statist., vol.7, no. 3, pp. 397-416, '
@@ -48,13 +47,12 @@ def skip(self, X, y, lmbd, fit_intercept):
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
-        # sklearn way of handling intercept: center y and X for dense data
-        if fit_intercept:
-            X, y, X_offset, y_offset, _ = _preprocess_data(
-                X, y, fit_intercept, return_mean=True, copy=True,
-            )
-            self.X_offset = X_offset
-            self.y_offset = y_offset
+        # Handling intercept: center y and X (dense data only)
+        if fit_intercept and not sparse.issparse(self.X):
+            self.X_offset = np.average(X, axis=0)
+            X -= self.X_offset
+            self.y_offset = np.average(y, axis=0)
+            y -= self.y_offset
 
         self.y, self.lmbd, self.fit_intercept = y, lmbd, fit_intercept
 
diff --git a/solvers/julia_pgd.py b/solvers/julia_pgd.py
index 84ef695d..58081e3b 100644
--- a/solvers/julia_pgd.py
+++ b/solvers/julia_pgd.py
@@ -8,7 +8,6 @@
 with safe_import_context() as import_ctx:
     import numpy as np
     from scipy.sparse import issparse
-    from sklearn.linear_model._base import _preprocess_data
     assert_julia_installed()
 
 
@@ -21,7 +20,6 @@ class Solver(JuliaSolver):
     # Config of the solver
     name = 'Julia-PGD'
     stopping_strategy = 'iteration'
-    requirements = ["scikit-learn"]
     references = [
         'I. Daubechies, M. Defrise and C. De Mol, '
         '"An iterative thresholding algorithm for linear inverse problems '
@@ -36,19 +34,18 @@ class Solver(JuliaSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX - fit intercept is not yet implemented in julia.jl for sparse X
         if fit_intercept and issparse(X):
-            return True, \
-                f"{self.name} doesn't handle fit_intercept with sparse data",
+            return (True,
+                f"{self.name} doesn't handle fit_intercept with sparse data")
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
-        # sklearn way of handling intercept: center y and X for dense data
-        if fit_intercept:
-            X, y, X_offset, y_offset, _ = _preprocess_data(
-                X, y, fit_intercept, return_mean=True, copy=True,
-            )
-            self.X_offset = X_offset
-            self.y_offset = y_offset
+        # Handling intercept: center y and X (dense data only)
+        if fit_intercept and not issparse(self.X):
+            self.X_offset = np.average(X, axis=0)
+            X -= self.X_offset
+            self.y_offset = np.average(y, axis=0)
+            y -= self.y_offset
 
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
diff --git a/solvers/l_bfgs_b.py b/solvers/l_bfgs_b.py
index 187a0b76..ce9d2767 100644
--- a/solvers/l_bfgs_b.py
+++ b/solvers/l_bfgs_b.py
@@ -6,7 +6,6 @@
     from numpy.linalg import norm
     from scipy.optimize import fmin_l_bfgs_b
     from scipy.sparse import issparse
-    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
@@ -32,19 +31,18 @@ class Solver(BaseSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX - intercept not implemented for sparse X for now
         if fit_intercept and issparse(X):
-            return True, \
-                f"{self.name} doesn't handle fit_intercept with sparse data",
+            return (True,
+                f"{self.name} doesn't handle fit_intercept with sparse data")
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
-        # sklearn way of handling intercept: center y and X for dense data
-        if fit_intercept:
-            X, y, X_offset, y_offset, _ = _preprocess_data(
-                X, y, fit_intercept, return_mean=True, copy=True,
-            )
-            self.X_offset = X_offset
-            self.y_offset = y_offset
+        # Handling intercept: center y and X (dense data only)
+        if fit_intercept and not issparse(self.X):
+            self.X_offset = np.average(X, axis=0)
+            X -= self.X_offset
+            self.y_offset = np.average(y, axis=0)
+            y -= self.y_offset
 
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
diff --git a/solvers/lightning.py b/solvers/lightning.py
index 9c444021..e8f1bf4c 100644
--- a/solvers/lightning.py
+++ b/solvers/lightning.py
@@ -6,7 +6,6 @@
     import numpy as np
     from scipy import sparse
     from lightning.regression import CDRegressor
-    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
@@ -15,7 +14,6 @@ class Solver(BaseSolver):
     install_cmd = 'conda'
     requirements = [
         'cython',
-        'scikit-learn',
         'pip:git+https://github.com/scikit-learn-contrib/lightning.git'
     ]
     references = [
@@ -36,14 +34,13 @@ def skip(self, X, y, lmbd, fit_intercept):
 
     def set_objective(self, X, y, lmbd, fit_intercept):
         # lightning has an attribut intercept_ but it is not handled properly
-        # (as it is simply set to zero). For this reason, we use the sklearn
-        # way of handling intercept: center y and X beforehand for dense data
-        if fit_intercept:
-            X, y, X_offset, y_offset, _ = _preprocess_data(
-                X, y, fit_intercept, return_mean=True, copy=True,
-            )
-            self.X_offset = X_offset
-            self.y_offset = y_offset
+        # (as it is simply set to zero). For this reason, we handle intercept
+        # manually: center y and X beforehand (for dense data only)
+        if fit_intercept and not sparse.issparse(self.X):
+            self.X_offset = np.average(X, axis=0)
+            X -= self.X_offset
+            self.y_offset = np.average(y, axis=0)
+            y -= self.y_offset
 
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
diff --git a/solvers/noncvx_pro.py b/solvers/noncvx_pro.py
index 396f063f..ef442da2 100644
--- a/solvers/noncvx_pro.py
+++ b/solvers/noncvx_pro.py
@@ -7,7 +7,6 @@
     from numpy.linalg import norm
     import scipy.optimize as sciop
     from scipy.sparse import issparse
-    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
@@ -15,7 +14,6 @@ class Solver(BaseSolver):
 
     stopping_strategy = 'iteration'
     stopping_criterion = SufficientDescentCriterion(eps=1e-10, patience=5)
-    requirements = ["scikit-learn"]
     references = [
         "Clarice Poon and Gabriel Peyré, "
         "'Smooth Bilevel Programming for Sparse Regularization', "
@@ -23,14 +21,12 @@ class Solver(BaseSolver):
     ]
 
     def set_objective(self, X, y, lmbd, fit_intercept):
-        # sklearn way of handling intercept: center y and X for dense data
-        # when X is sparse, X_offset is computed but X is not centered
-        if fit_intercept:
-            X, y, X_offset, y_offset, _ = _preprocess_data(
-                X, y, fit_intercept, return_mean=True, copy=True,
-            )
-            self.X_offset = X_offset
-            self.y_offset = y_offset
+        # Handling intercept: center y and X (dense data only)
+        if fit_intercept and not issparse(self.X):
+            self.X_offset = np.average(X, axis=0)
+            X -= self.X_offset
+            self.y_offset = np.average(y, axis=0)
+            y -= self.y_offset
 
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept
@@ -48,8 +44,8 @@ def skip(self, X, y, lmbd, fit_intercept):
         # unless fit_intercept is properly handled for sparse matrices
         # (by manually considering X_offset in calculations)
         if fit_intercept and issparse(X):
-            return True, \
-                f"{self.name} doesn't handle fit_intercept with sparse matrix",
+            return (True,
+                f"{self.name} doesn't handle fit_intercept with sparse matrix")
         return False, None
 
     def run(self, n_iter):
diff --git a/solvers/r_pgd.py b/solvers/r_pgd.py
index 8dd1a661..c168983a 100644
--- a/solvers/r_pgd.py
+++ b/solvers/r_pgd.py
@@ -6,7 +6,6 @@
 with safe_import_context() as import_ctx:
     import numpy as np
     from scipy.sparse import issparse
-    from sklearn.linear_model._base import _preprocess_data
 
     from rpy2 import robjects
     from rpy2.robjects import numpy2ri
@@ -22,7 +21,7 @@ class Solver(BaseSolver):
     name = "R-PGD"
 
     install_cmd = 'conda'
-    requirements = ['r-base', 'rpy2', 'scikit-learn']
+    requirements = ['r-base', 'rpy2']
     stopping_strategy = 'iteration'
     support_sparse = False
     references = [
@@ -38,17 +37,18 @@ class Solver(BaseSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # rpy2 does not directly support sparse matrices (workaround exists)
         if fit_intercept and issparse(X):
-            return True, \
-                f"{self.name} doesn't handle fit_intercept with sparse data"
+            return (True,
+                f"{self.name} doesn't handle fit_intercept with sparse data")
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
-        # sklearn way of handling intercept: center y and X (for dense data)
-        if fit_intercept:
-            X, y, self.X_offset, self.y_offset, _ = _preprocess_data(
-                X, y, fit_intercept, copy=True, return_mean=True
-            )
+        # Handling intercept: center y and X (dense data only)
+        if fit_intercept and not issparse(self.X):
+            self.X_offset = np.average(X, axis=0)
+            X -= self.X_offset
+            self.y_offset = np.average(y, axis=0)
+            y -= self.y_offset
 
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept

From fbf7f9539c2694336ea557fe72b3c80bb415247d Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Wed, 20 Jul 2022 11:21:45 +0200
Subject: [PATCH 15/17] Linter

---
 solvers/julia_pgd.py  | 5 ++---
 solvers/l_bfgs_b.py   | 3 +--
 solvers/noncvx_pro.py | 3 +--
 solvers/r_pgd.py      | 3 +--
 4 files changed, 5 insertions(+), 9 deletions(-)

diff --git a/solvers/julia_pgd.py b/solvers/julia_pgd.py
index 58081e3b..ba14f4e2 100644
--- a/solvers/julia_pgd.py
+++ b/solvers/julia_pgd.py
@@ -34,14 +34,13 @@ class Solver(JuliaSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX - fit intercept is not yet implemented in julia.jl for sparse X
         if fit_intercept and issparse(X):
-            return (True,
-                f"{self.name} doesn't handle fit_intercept with sparse data")
+            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
         # Handling intercept: center y and X (dense data only)
-        if fit_intercept and not issparse(self.X):
+        if fit_intercept and not sparse.issparse(self.X):
             self.X_offset = np.average(X, axis=0)
             X -= self.X_offset
             self.y_offset = np.average(y, axis=0)
diff --git a/solvers/l_bfgs_b.py b/solvers/l_bfgs_b.py
index ce9d2767..b4a33e71 100644
--- a/solvers/l_bfgs_b.py
+++ b/solvers/l_bfgs_b.py
@@ -31,8 +31,7 @@ class Solver(BaseSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX - intercept not implemented for sparse X for now
         if fit_intercept and issparse(X):
-            return (True,
-                f"{self.name} doesn't handle fit_intercept with sparse data")
+            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
 
         return False, None
 
diff --git a/solvers/noncvx_pro.py b/solvers/noncvx_pro.py
index ef442da2..d8e51ec8 100644
--- a/solvers/noncvx_pro.py
+++ b/solvers/noncvx_pro.py
@@ -44,8 +44,7 @@ def skip(self, X, y, lmbd, fit_intercept):
         # unless fit_intercept is properly handled for sparse matrices
         # (by manually considering X_offset in calculations)
         if fit_intercept and issparse(X):
-            return (True,
-                f"{self.name} doesn't handle fit_intercept with sparse matrix")
+            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
         return False, None
 
     def run(self, n_iter):
diff --git a/solvers/r_pgd.py b/solvers/r_pgd.py
index c168983a..da979985 100644
--- a/solvers/r_pgd.py
+++ b/solvers/r_pgd.py
@@ -37,8 +37,7 @@ class Solver(BaseSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # rpy2 does not directly support sparse matrices (workaround exists)
         if fit_intercept and issparse(X):
-            return (True,
-                f"{self.name} doesn't handle fit_intercept with sparse data")
+            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
 
         return False, None
 

From d81d5b79482076a6760976aa642caef4a21203e3 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Wed, 20 Jul 2022 11:26:39 +0200
Subject: [PATCH 16/17] Linter

---
 solvers/julia_pgd.py  | 7 +++++--
 solvers/l_bfgs_b.py   | 5 ++++-
 solvers/noncvx_pro.py | 5 ++++-
 solvers/r_pgd.py      | 5 ++++-
 4 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/solvers/julia_pgd.py b/solvers/julia_pgd.py
index ba14f4e2..8d64de29 100644
--- a/solvers/julia_pgd.py
+++ b/solvers/julia_pgd.py
@@ -34,13 +34,16 @@ class Solver(JuliaSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX - fit intercept is not yet implemented in julia.jl for sparse X
         if fit_intercept and issparse(X):
-            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data"
+            )
 
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
         # Handling intercept: center y and X (dense data only)
-        if fit_intercept and not sparse.issparse(self.X):
+        if fit_intercept and not issparse(self.X):
             self.X_offset = np.average(X, axis=0)
             X -= self.X_offset
             self.y_offset = np.average(y, axis=0)
diff --git a/solvers/l_bfgs_b.py b/solvers/l_bfgs_b.py
index b4a33e71..498bd5ab 100644
--- a/solvers/l_bfgs_b.py
+++ b/solvers/l_bfgs_b.py
@@ -31,7 +31,10 @@ class Solver(BaseSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # XXX - intercept not implemented for sparse X for now
         if fit_intercept and issparse(X):
-            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data"
+            )
 
         return False, None
 
diff --git a/solvers/noncvx_pro.py b/solvers/noncvx_pro.py
index d8e51ec8..54992d45 100644
--- a/solvers/noncvx_pro.py
+++ b/solvers/noncvx_pro.py
@@ -44,7 +44,10 @@ def skip(self, X, y, lmbd, fit_intercept):
         # unless fit_intercept is properly handled for sparse matrices
         # (by manually considering X_offset in calculations)
         if fit_intercept and issparse(X):
-            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data"
+            )
         return False, None
 
     def run(self, n_iter):
diff --git a/solvers/r_pgd.py b/solvers/r_pgd.py
index da979985..29031b8f 100644
--- a/solvers/r_pgd.py
+++ b/solvers/r_pgd.py
@@ -37,7 +37,10 @@ class Solver(BaseSolver):
     def skip(self, X, y, lmbd, fit_intercept):
         # rpy2 does not directly support sparse matrices (workaround exists)
         if fit_intercept and issparse(X):
-            return True, f"{self.name} doesn't handle fit_intercept with sparse data"
+            return (
+                True,
+                f"{self.name} doesn't handle fit_intercept with sparse data"
+            )
 
         return False, None
 

From fb38f68777e240035bf5513a5d3ae02567ebd525 Mon Sep 17 00:00:00 2001
From: "Cassio F. Dantas" <cassiofragadantas@gmail.com>
Date: Wed, 20 Jul 2022 13:46:03 +0200
Subject: [PATCH 17/17] removing scikit-learn dependency in python-pgd solver

---
 solvers/python_pgd.py | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/solvers/python_pgd.py b/solvers/python_pgd.py
index fb23160b..f7bff705 100644
--- a/solvers/python_pgd.py
+++ b/solvers/python_pgd.py
@@ -3,12 +3,10 @@
 with safe_import_context() as import_ctx:
     import numpy as np
     from scipy import sparse
-    from sklearn.linear_model._base import _preprocess_data
 
 
 class Solver(BaseSolver):
     name = 'Python-PGD'  # proximal gradient, optionally accelerated
-    requirements = ['scikit-learn']
     stopping_strategy = "callback"
 
     # any parameter defined here is accessible as a class attribute
@@ -34,13 +32,12 @@ def skip(self, X, y, lmbd, fit_intercept):
         return False, None
 
     def set_objective(self, X, y, lmbd, fit_intercept):
-        # sklearn way of handling intercept: center y and X for dense data
-        if fit_intercept:
-            X, y, X_offset, y_offset, _ = _preprocess_data(
-                X, y, fit_intercept, return_mean=True, copy=True,
-            )
-            self.X_offset = X_offset
-            self.y_offset = y_offset
+        # Handling intercept: center y and X (dense data only)
+        if fit_intercept and not sparse.issparse(self.X):
+            self.X_offset = np.average(X, axis=0)
+            X -= self.X_offset
+            self.y_offset = np.average(y, axis=0)
+            y -= self.y_offset
 
         self.X, self.y, self.lmbd = X, y, lmbd
         self.fit_intercept = fit_intercept