Repository: scikit-learn/scikit-learn
Pull request: ENH PolynomialFeatures better handle degree=0
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/main/CONTRIBUTING.md
-->

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->
Fixes #19551 .

#### What does this implement/fix? Explain your changes.
This PR makes the following two changes:
1. When `degree=0`, a ValueError will be raised if `include_bias=False`, a single constant column will be returned if `include_bias=True` regardless of the input number of features
2. Fixes an edge case. When passing in a sparse matrix, the output of `degree=0, include_bias=True` used to have n+1 columns with n being number of input features, now it outputs 1 column correctly.

<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
PolynomialFeatures doesn't work correctly when degree=0
#### Describe the bug
PolynomialFeatures.transform fails when degree=0 and the input has more than one column.

#### Steps/Code to Reproduce

```python
from sklearn.preprocessing import PolynomialFeatures
PolynomialFeatures(degree=0, include_bias=True).fit_transform(np.ones((10,1))) # works
PolynomialFeatures(degree=0, include_bias=False).fit_transform(np.ones((10,1))) # works
PolynomialFeatures(degree=0, include_bias=True).fit_transform(np.ones((10,2))) # raises a ValueError due to a broadcasting problem
PolynomialFeatures(degree=0, include_bias=False).fit_transform(np.ones((10,2))) # raises a ValueError due to a broadcasting problem
```
#### Expected Results
No error is raised by any of the examples.

#### Actual Results
Transforming an array with more than one column fails.

#### Versions

System:
    python: 3.6.6 (v3.6.6:4cf1f54eb7, Jun 27 2018, 03:37:03) [MSC v.1900 64 bit (AMD64)]
executable: C:\Program Files (x86)\Microsoft Visual Studio\Shared\Python36_64\python.exe
   machine: Windows-10-10.0.19041-SP0

Python dependencies:
          pip: 21.0.1
   setuptools: 51.0.0
      sklearn: 0.24.1
        numpy: 1.18.5
        scipy: 1.5.4
       Cython: 0.29.14
       pandas: 1.0.5
   matplotlib: None
       joblib: 1.0.0
threadpoolctl: 2.1.0

Built with OpenMP: True

Hints from issue discussion:
Why should we allow degree=0? I suspect we should raise an explicit error for degree <1
Take
Seems pretty easy so I'd like to work on it
> Why should we allow degree=0? I suspect we should raise an explicit error for degree <1

Maybe one could argue that, for educational purpose, it could be interesting to generate a single constant feature so that fitting a linear model (with an intercept) on this would be equivalent to a dummy classifier / regressor predicting the mean class probability / regression value.

But maybe this is too much of a edge case to consider of value in scikit-learn.
+1 for raising a `ValueError` stating that degree should be greater than or equal to 1.
To me, it seems like there's exactly one natural definition for the behavior when degree=0.  Given an array with n columns, when interaction_only is False, for degree d>1 the output will have exactly choose(n+d, n) columns when include_bias is True and one less column when it's False.  Therefore, when degree is 0, this should extend to a single column of ones if include_bias is True or an empty array if it's False.

As to why to allow it, I see a few reasons.  One is that it seems weird to disallow edge cases just because they're edge cases; I'm allowed to call `scipy.special.comb(5,0)` even though if that threw an exception I could always work around it by just special casing it and plugging in 1 by myself when I want to call comb.  Another is composability: as a library writer I may be using PolynomialFeatures as a part of something bigger (maybe it's a configurable piece of a ColumnTransformer or something) where my caller is in control of the degree and where 0 makes sense as part of the bigger picture.  Finally, if it were supported, `PolynomialFeatures(degree=0)` would arguably be the simplest built-in transformer for generating a single constant feature, which is often useful.
> As to why to allow it, I see a few reasons. One is that it seems weird to disallow edge cases just because they're edge cases

Having a transformer that outputs an empty array with 0 column (if `include_bias=False`) is weird and therefore surprising. All scikit-learn estimators will raise an error when fed with such a weird datastructure. Unsuspecting users might not understand what in their pipeline is causing such a crash in the downstream estimator. When having to make API design decisions we often try to follow the [principle of least astonishment](https://en.wikipedia.org/wiki/Principle_of_least_astonishment).

> Another is composability: as a library writer I may be using PolynomialFeatures as a part of something bigger (maybe it's a configurable piece of a ColumnTransformer or something) where my caller is in control of the degree and where 0 makes sense as part of the bigger picture.

I get the composability argument but it would be stronger with a concrete case where this is actually useful.

> PolynomialFeatures(degree=0) would arguably be the simplest built-in transformer for generating a single constant feature, which is often useful.

Interesting that you raise this now, since we were just having the dual discussion about warning the users in the presence of constant features which are mostly useless for machine learning (linear models have a built-in biais, tree-based models will ignore that feature, ...): #19547. I would be interested in a non-synthetic use case where feeding a constant feature to a predictive pipeline is practically useful.
To me we should at least raise a `ValueError` when fitting a `PolynomialFeatures` with `degree=0` and `include_bias=False`.

I could accept to have `PolynomialFeatures(degree=0, include_bias=True`) that outputs a constant column of ones for the sake of composability but I don't really see the practical value, so I would not not be strongly opposed to just reject `degree=0` irrespective of the value of `include_bias` for the sake of simplicity of both the docstring and the code :)
> Interesting that you raise this now, since we were just having the dual discussion about warning the users in the presence of constant features which are mostly useless for machine learning (linear models have a built-in biais, tree-based models will ignore that feature, ...): #19547. I would be interested in a non-synthetic use case where feeding a constant feature to a predictive pipeline is practically useful.

I can describe the use case that resulted in my opening this issue, but it's a bit involved.  We're working on a library for computing heterogeneous treatment effects.  One of our estimators works like this:
1. Users give us one first stage model for estimating an outcome from features and another for estimating a treatment from features
2. We use the first stage models to compute outcome and treatment residuals, which we store for later use
3. The user also gives us a featurizer and a final model
4. We use the final model to regress the outcome residuals against the cross-product of the treatment residuals and the featurized features to estimate a treatment effect

We can run those last two steps multiple times against the residuals stored in step 2.  Then if we use a constant featurizer, we'll compute the average treatment effect, while if we use an identity featurizer that lets us compute heterogenous effects (and featurizers that keep only a subset of the features let us consider heterogeneity across only that subset, etc.).

As I tried to allude to in my compositionality comment, the constant featurizer might not be super useful as the only direct input to a model, but when used as a configurable part of a bigger system it makes sense (e.g. here we're multiplying each featurizer output by each column of another array).  We have also had other cases (although I can't think of a concrete example off hand) where a featurizer returning a zero-width array is useful because we're concatenating the results of different featurizers which are provided by users and a user might want to omit some data from a later step and can use a featurizer returning an empty array to do so.
> > As to why to allow it, I see a few reasons. One is that it seems weird to disallow edge cases just because they're edge cases
> 
> Having a transformer that outputs an empty array with 0 column (if `include_bias=False`) is weird and therefore surprising. All scikit-learn estimators will raise an error when fed with such a weird datastructure. Unsuspecting users might not understand what in their pipeline is causing such a crash in the downstream estimator. When having to make API design decisions we often try to follow the [principle of least astonishment](https://en.wikipedia.org/wiki/Principle_of_least_astonishment).

I get what you're saying, but my colleague and I were both surprised by this _not_ working.  It seems unlikely to me that there are users setting the degree to 0 with or without a bias and then being surprised by the results since there is no other reasonable behavior when degree is 0, but even if that is a concern, why not warn instead of raising an exception so that the (admittedly uncommon) valid uses for such a featurizer are still possible?

Commit message: raise error for PolynomialFeature when include_bias=False and degree=0

Patch to review:
diff --git a/sklearn/preprocessing/_polynomial.py b/sklearn/preprocessing/_polynomial.py
index 09306cc56ab11..f8b67ad2b1cea 100644
--- a/sklearn/preprocessing/_polynomial.py
+++ b/sklearn/preprocessing/_polynomial.py
@@ -291,6 +291,13 @@ def fit(self, X, y=None):
                 raise ValueError(
                     f"degree must be a non-negative integer, got {self.degree}."
                 )
+
+            if self.degree == 0 and not self.include_bias:
+                raise ValueError(
+                    "Setting degree to zero and include_bias to False will result to an"
+                    " empty dataframe."
+                )
+
             self._min_degree = 0
             self._max_degree = self.degree
         elif (
@@ -380,7 +387,7 @@ def transform(self, X):
                 to_stack.append(
                     sparse.csc_matrix(np.ones(shape=(n_samples, 1), dtype=X.dtype))
                 )
-            if self._min_degree <= 1:
+            if self._min_degree <= 1 and self._max_degree > 0:
                 to_stack.append(X)
             for deg in range(max(2, self._min_degree), self._max_degree + 1):
                 Xp_next = _csr_polynomial_expansion(
@@ -442,6 +449,9 @@ def transform(self, X):
             else:
                 current_col = 0
 
+            if self._max_degree == 0:
+                return XP
+
             # degree 1 term
             XP[:, current_col : current_col + n_features] = X
             index = list(range(current_col, current_col + n_features))
diff --git a/sklearn/preprocessing/tests/test_polynomial.py b/sklearn/preprocessing/tests/test_polynomial.py
index 2129247125d6c..ba8b6d14cf25f 100644
--- a/sklearn/preprocessing/tests/test_polynomial.py
+++ b/sklearn/preprocessing/tests/test_polynomial.py
@@ -898,3 +898,21 @@ def test_get_feature_names_deprecated(Transformer):
     msg = "get_feature_names is deprecated in 1.0"
     with pytest.warns(FutureWarning, match=msg):
         poly.get_feature_names()
+
+
+def test_polynomial_features_behaviour_on_zero_degree():
+    """Check that PolynomialFeatures raises error when degree=0 and include_bias=False,
+    and output a single constant column when include_bias=True
+    """
+    X = np.ones((10, 2))
+    poly = PolynomialFeatures(degree=0, include_bias=False)
+    err_msg = (
+        "Setting degree to zero and include_bias to False will result to an empty"
+        " dataframe."
+    )
+    with pytest.raises(ValueError, match=err_msg):
+        poly.fit_transform(X)
+
+    for _X in [X, sparse.csr_matrix(X), sparse.csc_matrix(X)]:
+        poly = PolynomialFeatures(degree=0, include_bias=True)
+        assert poly.fit_transform(_X).shape[1] == 1


Code review comments:
diff_hunk: @@ -291,6 +291,13 @@ def fit(self, X, y=None):
                 raise ValueError(
                     f"degree must be a non-negative integer, got {self.degree}."
                 )
+
+            if self.degree == 0 and not self.include_bias:
path: sklearn/preprocessing/_polynomial.py
text: @user1:
```suggestion
            elif self.degree == 0 and not self.include_bias:
```

diff_hunk: @@ -291,6 +291,13 @@ def fit(self, X, y=None):
                 raise ValueError(
                     f"degree must be a non-negative integer, got {self.degree}."
                 )
+
+            if self.degree == 0 and not self.include_bias:
+                raise ValueError(
+                    "Setting degree to zero and include_bias to False will result to an"
+                    " empty dataframe."
path: sklearn/preprocessing/_polynomial.py
text: @user1:
```suggestion
                    "Setting degree to zero and include_bias to False would result in an"
                    " empty output array."
```

diff_hunk: @@ -898,3 +898,21 @@ def test_get_feature_names_deprecated(Transformer):
     msg = "get_feature_names is deprecated in 1.0"
     with pytest.warns(FutureWarning, match=msg):
         poly.get_feature_names()
+
+
+def test_polynomial_features_behaviour_on_zero_degree():
+    """Check that PolynomialFeatures raises error when degree=0 and include_bias=False,
+    and output a single constant column when include_bias=True
+    """
+    X = np.ones((10, 2))
+    poly = PolynomialFeatures(degree=0, include_bias=False)
+    err_msg = (
+        "Setting degree to zero and include_bias to False will result to an empty"
+        " dataframe."
+    )
+    with pytest.raises(ValueError, match=err_msg):
+        poly.fit_transform(X)
+
+    for _X in [X, sparse.csr_matrix(X), sparse.csc_matrix(X)]:
+        poly = PolynomialFeatures(degree=0, include_bias=True)
+        assert poly.fit_transform(_X).shape[1] == 1
path: sklearn/preprocessing/tests/test_polynomial.py
text: @user1:
This should be a stronger assert statement, something like `assert_array_equal(_X, np.ones(...))`.

@author:
Updated as suggested :)