Repository: scikit-learn/scikit-learn
Pull request: FIX LinearRegression sparse + intercept + sample_weight
Partial fix for #15438
Fixes https://github.com/scikit-learn/scikit-learn/issues/19578

This PR fixes the issue for LinearRegression only. I intend to check other linear models in separate PR(s).

To account for sample weight, we rescale X and y. When X is sparse we also need to rescale X_offset because the centering is implicit.

ping @agramfort 

Problem statement:
LinearRegression fits wrongly on csr sparse matrix with sample weights
<!--
Before submitting a bug, please make sure the issue hasn't been already
addressed by searching through the past issues.
-->

#### Describe the bug

Fitting a model from a csr sparse matrix while passing sample_weight gives a very wrong fit.

This report is for `LinearRegression`, but I've been experiencing wrong behavior in a similar setup with `Ridge`.

#### Steps/Code to Reproduce

Taken from a small sample of a real life dataset, just copy&paste:

```python
import numpy as np
import scipy as sp

from sklearn.linear_model import LinearRegression

X = np.matrix(
       [[5.33558618e+00, 1.42240293e-03, 3.79195466e-07, 5.33558618e+00,
         1.42240293e-03, 3.79195466e-07, 5.33558618e+00, 1.42240293e-03,
         3.79195466e-07, 1.98842115e+00, 5.30089098e-04, 1.41315361e-07,
         1.98842115e+00, 5.30089098e-04, 1.41315361e-07, 1.98842115e+00,
         5.30089098e-04, 1.41315361e-07, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00],
        [4.02363338e+01, 2.14360316e-03, 1.14201123e-07, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 6.57402127e+00, 3.50233021e-04, 1.86587728e-08,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 4.02363338e+01, 2.14360316e-03,
         1.14201123e-07, 4.02363338e+01, 2.14360316e-03, 1.14201123e-07,
         6.57402127e+00, 3.50233021e-04, 1.86587728e-08, 6.57402127e+00,
         3.50233021e-04, 1.86587728e-08, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00],
        [4.73842584e+01, 6.78046657e-03, 9.70253169e-07, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 5.75883825e+00, 8.24062917e-04, 1.17919563e-07,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 4.73842584e+01, 6.78046657e-03,
         9.70253169e-07, 4.73842584e+01, 6.78046657e-03, 9.70253169e-07,
         5.75883825e+00, 8.24062917e-04, 1.17919563e-07, 5.75883825e+00,
         8.24062917e-04, 1.17919563e-07, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00],
        [6.73948024e+01, 1.51450045e-02, 3.40339540e-06, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 5.09430891e+00, 1.14479646e-03, 2.57259416e-07,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 6.73948024e+01, 1.51450045e-02,
         3.40339540e-06, 6.73948024e+01, 1.51450045e-02, 3.40339540e-06,
         5.09430891e+00, 1.14479646e-03, 2.57259416e-07, 5.09430891e+00,
         1.14479646e-03, 2.57259416e-07, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00],
        [2.20048734e+01, 4.36703062e-03, 8.66669673e-07, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 3.18056958e+00, 6.31207665e-04, 1.25267851e-07,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
         0.00000000e+00, 0.00000000e+00, 2.20048734e+01, 4.36703062e-03,
         8.66669673e-07, 2.20048734e+01, 4.36703062e-03, 8.66669673e-07,
         3.18056958e+00, 6.31207665e-04, 1.25267851e-07, 3.18056958e+00,
         6.31207665e-04, 1.25267851e-07]])

y = np.array([ 0.        , -0.0003287 ,  0.        , -0.0011448 , -0.00063121])
w = np.array([50000.,  2500., 50000.,  2500.,  2500.])

lr = LinearRegression().fit(X, y, sample_weight=w)
print(
    "dense:",
    lr.score(X, y, sample_weight=w),
    lr.coef_[:5],
)

lr = LinearRegression().fit(X, y)
print(
    "dense (no weights):",
    lr.score(X, y),
    lr.coef_[:5],
)

X = sp.sparse.csr_matrix(X)
lr = LinearRegression().fit(X, y, sample_weight=w)
print(
    "sparse:",
    lr.score(X, y, sample_weight=w),
    lr.coef_[:5],
)

lr = LinearRegression().fit(X, y)
print(
    "sparse (no weights):",
    lr.score(X, y),
    lr.coef_[:5],
)
```

#### Expected Results

Both weighted results are the same.

Both unweighted results are the same.

#### Actual Results

Weighted results differ. The one for the sparse matrix is clearly wrong.

Unweighted results are the same.

#### Versions
<!--
Please run the following snippet and paste the output below.
For scikit-learn >= 0.20:
import sklearn; sklearn.show_versions()
For scikit-learn < 0.20:
import platform; print(platform.platform())
import sys; print("Python", sys.version)
import numpy; print("NumPy", numpy.__version__)
import scipy; print("SciPy", scipy.__version__)
import sklearn; print("Scikit-Learn", sklearn.__version__)
import imblearn; print("Imbalanced-Learn", imblearn.__version__)
-->

```
System:
    python: 3.8.6 (default, Nov 20 2020, 23:57:10)  [Clang 12.0.0 (clang-1200.0.32.27)]
executable: /Users/carlos/Base/Venvs/Default/bin/python3
   machine: macOS-11.2-x86_64-i386-64bit

Python dependencies:
          pip: 21.0.1
   setuptools: 49.2.1
      sklearn: 0.24.1
        numpy: 1.20.1
        scipy: 1.6.0
       Cython: 0.29.21
       pandas: 1.2.2
   matplotlib: 3.1.2
       joblib: 1.0.0
threadpoolctl: 2.1.0

Built with OpenMP: True
```


<!-- Thanks for contributing! -->

Hints from issue discussion:
OTOH, here is a simple example that does work:

```python
X = sp.sparse.csr_matrix(np.array([[1, 0]] * 100 + [[0, 1]] * 100))
y = np.array([1] * 80 + [0] * 20 + [1] * 50 + [0] * 50)
w = np.array([1] * 80 + [4] * 20 + [1] * 50 + [2] * 50)
LinearRegression(fit_intercept=False).fit(X, y, sample_weight=w).coef_

=> array([0.5       , 0.33333333])
```
Does it look like the scale of the weights might play a part? Can you check
how well it is working with weights of different scale?

No, it doesn't. Dividing them by 100 or 10000 still throws different results for the sparse/weighted case. Using an all ones weight vector works fine but as soon as I change some of the weights to, say, 2, the sparse/weighted case fails again.
Another fact to take into account is that the matrix X'X is very ill conditioned, so if some kind of iterative solver is chosen for sparse input (I don't know if that's the case), then you have problems in both fronts: observations wildly differ in weight while features have very different "geometries", then it might be very difficult to accommodate a step schedule that converges.
Hi @memeplex I believe your findings are similar to #15438 in particular as said in @lorentzenchr [comment](https://github.com/scikit-learn/scikit-learn/issues/15438#issuecomment-550432207).
If so, we can close?

Commit message: address comments

Patch to review:
diff --git a/doc/whats_new/v1.1.rst b/doc/whats_new/v1.1.rst
index f64a6bda6ea95..7786870dc523f 100644
--- a/doc/whats_new/v1.1.rst
+++ b/doc/whats_new/v1.1.rst
@@ -594,6 +594,10 @@ Changelog
   :class:`linear_model.ARDRegression` now preserve float32 dtype. :pr:`9087` by
   :user:`Arthur Imbert <Henley13>` and :pr:`22525` by :user:`Meekail Zain <micky774>`.
 
+- |Fix| The `intercept_` attribute of :class:`LinearRegression` is now correctly
+  computed in the presence of sample weights when the input is sparse.
+  :pr:`22891` by :user:`Jérémie du Boisberranger <jeremiedbb>`. 
+
 :mod:`sklearn.manifold`
 .......................
 
diff --git a/sklearn/linear_model/_base.py b/sklearn/linear_model/_base.py
index f0322c3924426..b0c708d3516c3 100644
--- a/sklearn/linear_model/_base.py
+++ b/sklearn/linear_model/_base.py
@@ -325,11 +325,14 @@ def _preprocess_data(
 # sample_weight makes the refactoring tricky.
 
 
-def _rescale_data(X, y, sample_weight):
+def _rescale_data(X, y, sample_weight, square_root=True):
     """Rescale data sample-wise by square root of sample_weight.
 
     For many linear models, this enables easy support for sample_weight.
 
+    Set square_root=False if the square root of the sample weights has already been
+    done prior to calling this function.
+
     Returns
     -------
     X_rescaled : {array-like, sparse matrix}
@@ -340,7 +343,8 @@ def _rescale_data(X, y, sample_weight):
     sample_weight = np.asarray(sample_weight)
     if sample_weight.ndim == 0:
         sample_weight = np.full(n_samples, sample_weight, dtype=sample_weight.dtype)
-    sample_weight = np.sqrt(sample_weight)
+    if square_root:
+        sample_weight = np.sqrt(sample_weight)
     sw_matrix = sparse.dia_matrix((sample_weight, 0), shape=(n_samples, n_samples))
     X = safe_sparse_dot(sw_matrix, X)
     y = safe_sparse_dot(sw_matrix, y)
@@ -676,10 +680,9 @@ def fit(self, X, y, sample_weight=None):
             X, y, accept_sparse=accept_sparse, y_numeric=True, multi_output=True
         )
 
-        if sample_weight is not None:
-            sample_weight = _check_sample_weight(
-                sample_weight, X, dtype=X.dtype, only_non_negative=True
-            )
+        sample_weight = _check_sample_weight(
+            sample_weight, X, dtype=X.dtype, only_non_negative=True
+        )
 
         X, y, X_offset, y_offset, X_scale = _preprocess_data(
             X,
@@ -691,9 +694,9 @@ def fit(self, X, y, sample_weight=None):
             return_mean=True,
         )
 
-        if sample_weight is not None:
-            # Sample weight can be implemented via a simple rescaling.
-            X, y = _rescale_data(X, y, sample_weight)
+        # Sample weight can be implemented via a simple rescaling.
+        sample_weight = np.sqrt(sample_weight)
+        X, y = _rescale_data(X, y, sample_weight, square_root=False)
 
         if self.positive:
             if y.ndim < 2:
@@ -708,10 +711,10 @@ def fit(self, X, y, sample_weight=None):
             X_offset_scale = X_offset / X_scale
 
             def matvec(b):
-                return X.dot(b) - b.dot(X_offset_scale)
+                return X.dot(b) - sample_weight * b.dot(X_offset_scale)
 
             def rmatvec(b):
-                return X.T.dot(b) - X_offset_scale * np.sum(b)
+                return X.T.dot(b) - X_offset_scale * b.dot(sample_weight)
 
             X_centered = sparse.linalg.LinearOperator(
                 shape=X.shape, matvec=matvec, rmatvec=rmatvec
diff --git a/sklearn/linear_model/tests/test_base.py b/sklearn/linear_model/tests/test_base.py
index 6b89b7e25e899..f59a5f6f08c5c 100644
--- a/sklearn/linear_model/tests/test_base.py
+++ b/sklearn/linear_model/tests/test_base.py
@@ -13,7 +13,6 @@
 
 from sklearn.utils._testing import assert_array_almost_equal
 from sklearn.utils._testing import assert_array_equal
-from sklearn.utils._testing import assert_almost_equal
 from sklearn.utils._testing import assert_allclose
 from sklearn.utils import check_random_state
 
@@ -26,6 +25,7 @@
 from sklearn.datasets import make_regression
 from sklearn.datasets import load_iris
 from sklearn.preprocessing import StandardScaler
+from sklearn.preprocessing import add_dummy_feature
 
 rng = np.random.RandomState(0)
 rtol = 1e-6
@@ -55,45 +55,42 @@ def test_linear_regression():
     assert_array_almost_equal(reg.predict(X), [0])
 
 
-def test_linear_regression_sample_weights():
-    # TODO: loop over sparse data as well
-
+@pytest.mark.parametrize("array_constr", [np.array, sparse.csr_matrix])
+@pytest.mark.parametrize("fit_intercept", [True, False])
+def test_linear_regression_sample_weights(array_constr, fit_intercept):
     rng = np.random.RandomState(0)
 
     # It would not work with under-determined systems
-    for n_samples, n_features in ((6, 5),):
+    n_samples, n_features = 6, 5
 
-        y = rng.randn(n_samples)
-        X = rng.randn(n_samples, n_features)
-        sample_weight = 1.0 + rng.rand(n_samples)
+    X = array_constr(rng.normal(size=(n_samples, n_features)))
+    y = rng.normal(size=n_samples)
 
-        for intercept in (True, False):
+    sample_weight = 1.0 + rng.uniform(size=n_samples)
 
-            # LinearRegression with explicit sample_weight
-            reg = LinearRegression(fit_intercept=intercept)
-            reg.fit(X, y, sample_weight=sample_weight)
-            coefs1 = reg.coef_
-            inter1 = reg.intercept_
+    # LinearRegression with explicit sample_weight
+    reg = LinearRegression(fit_intercept=fit_intercept)
+    reg.fit(X, y, sample_weight=sample_weight)
+    coefs1 = reg.coef_
+    inter1 = reg.intercept_
 
-            assert reg.coef_.shape == (X.shape[1],)  # sanity checks
-            assert reg.score(X, y) > 0.5
+    assert reg.coef_.shape == (X.shape[1],)  # sanity checks
+    assert reg.score(X, y) > 0.5
 
-            # Closed form of the weighted least square
-            # theta = (X^T W X)^(-1) * X^T W y
-            W = np.diag(sample_weight)
-            if intercept is False:
-                X_aug = X
-            else:
-                dummy_column = np.ones(shape=(n_samples, 1))
-                X_aug = np.concatenate((dummy_column, X), axis=1)
+    # Closed form of the weighted least square
+    # theta = (X^T W X)^(-1) @ X^T W y
+    W = np.diag(sample_weight)
+    X_aug = X if not fit_intercept else add_dummy_feature(X)
 
-            coefs2 = linalg.solve(X_aug.T.dot(W).dot(X_aug), X_aug.T.dot(W).dot(y))
+    Xw = X_aug.T @ W @ X_aug
+    yw = X_aug.T @ W @ y
+    coefs2 = linalg.solve(Xw, yw)
 
-            if intercept is False:
-                assert_array_almost_equal(coefs1, coefs2)
-            else:
-                assert_array_almost_equal(coefs1, coefs2[1:])
-                assert_almost_equal(inter1, coefs2[0])
+    if not fit_intercept:
+        assert_allclose(coefs1, coefs2)
+    else:
+        assert_allclose(coefs1, coefs2[1:])
+        assert_allclose(inter1, coefs2[0])
 
 
 def test_raises_value_error_if_positive_and_sparse():


Code review comments:
diff_hunk: @@ -708,10 +711,10 @@ def fit(self, X, y, sample_weight=None):
             X_offset_scale = X_offset / X_scale
 
             def matvec(b):
-                return X.dot(b) - b.dot(X_offset_scale)
+                return X.dot(b) - sample_weight * b.dot(X_offset_scale)
 
             def rmatvec(b):
-                return X.T.dot(b) - X_offset_scale * np.sum(b)
+                return X.T.dot(b) - X_offset_scale * b.dot(sample_weight)
path: sklearn/linear_model/_base.py
text: @user1:
```suggestion
                return X.T.dot(b) - X_offset_scale * b.dot(sample_weight_sqrt)
```

diff_hunk: @@ -708,10 +711,10 @@ def fit(self, X, y, sample_weight=None):
             X_offset_scale = X_offset / X_scale
 
             def matvec(b):
-                return X.dot(b) - b.dot(X_offset_scale)
+                return X.dot(b) - sample_weight * b.dot(X_offset_scale)
path: sklearn/linear_model/_base.py
text: @user1:
Then we use the sqrt variant explicitly only when needed:

```suggestion
                return X.dot(b) - sample_weight_sqrt * b.dot(X_offset_scale)
```