Repository: scikit-learn/scikit-learn
Pull request: [MRG+2] LogisticRegression convert to float64 (newton-cg)
<!--
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/master/CONTRIBUTING.md#Contributing-Pull-Requests
-->
#### Reference Issue
Fixes #8769

#### What does this implement/fix? Explain your changes.
Avoids logistic regression to aggressively cast the data to `np.float64` when `np.float32` is supplied. 

#### Any other comments?
(only for the `newton-cg` case)

<!--
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:
LogisticRegression convert to float64
Looking at the code of LogisticRegression, I have just noticed that it converts automatically the data to float64. I would expect at least the SAG, SAGA, newton-cg and lbfgs solvers to be able to work with float32.

Also, in a similar line of thoughts, the code converts to C-ordered data. How useful / necessary is this for solvers others than liblinear.

I am asking these questions because it seems that we could reduce the memory footprint of LogisticRegression.

Ping @arthurmensch for the SAG part.

Hints from issue discussion:
I can work on this
Great! Handling float32 is not that hard: We need fused type : you can have
a look at cd_fast.pyx. Handling F ordered arrays is possible but requires
to change the way we access some arrays in place of the code where we use
explicit pointers.

On Thu, Apr 20, 2017, 5:19 PM Joan Massich <notifications@github.com> wrote:

> I can work on this
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/scikit-learn/scikit-learn/issues/8769#issuecomment-295775406>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/AD48Y4m29WJ2a8SIOcVzrOsd3E0-AHAsks5rx3dkgaJpZM4NC9c0>
> .
>

Maybe not for L-BFGS:
https://github.com/scipy/scipy/issues/4873

Commit message: Add a test where X.dtype is different of y.dtype

Patch to review:
diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py
index fa7bc7f32e4be..53cc391156d6b 100644
--- a/sklearn/linear_model/logistic.py
+++ b/sklearn/linear_model/logistic.py
@@ -337,10 +337,12 @@ def _multinomial_loss_grad(w, X, Y, alpha, sample_weight):
     n_classes = Y.shape[1]
     n_features = X.shape[1]
     fit_intercept = (w.size == n_classes * (n_features + 1))
-    grad = np.zeros((n_classes, n_features + bool(fit_intercept)))
+    grad = np.zeros((n_classes, n_features + bool(fit_intercept)),
+                    dtype=X.dtype)
     loss, p, w = _multinomial_loss(w, X, Y, alpha, sample_weight)
     sample_weight = sample_weight[:, np.newaxis]
     diff = sample_weight * (p - Y)
+    diff = diff.astype(X.dtype)
     grad[:, :n_features] = safe_sparse_dot(diff.T, X)
     grad[:, :n_features] += alpha * w
     if fit_intercept:
@@ -608,10 +610,10 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True,
     # and check length
     # Otherwise set them to 1 for all examples
     if sample_weight is not None:
-        sample_weight = np.array(sample_weight, dtype=np.float64, order='C')
+        sample_weight = np.array(sample_weight, dtype=X.dtype, order='C')
         check_consistent_length(y, sample_weight)
     else:
-        sample_weight = np.ones(X.shape[0])
+        sample_weight = np.ones(X.shape[0], dtype=X.dtype)
 
     # If class_weights is a dict (provided by the user), the weights
     # are assigned to the original labels. If it is "balanced", then
@@ -624,10 +626,10 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True,
     # For doing a ovr, we need to mask the labels first. for the
     # multinomial case this is not necessary.
     if multi_class == 'ovr':
-        w0 = np.zeros(n_features + int(fit_intercept))
+        w0 = np.zeros(n_features + int(fit_intercept), dtype=X.dtype)
         mask_classes = np.array([-1, 1])
         mask = (y == pos_class)
-        y_bin = np.ones(y.shape, dtype=np.float64)
+        y_bin = np.ones(y.shape, dtype=X.dtype)
         y_bin[~mask] = -1.
         # for compute_class_weight
 
@@ -648,7 +650,7 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True,
             Y_multi = le.fit_transform(y)
 
         w0 = np.zeros((classes.size, n_features + int(fit_intercept)),
-                      order='F')
+                      order='F', dtype=X.dtype)
 
     if coef is not None:
         # it must work both giving the bias term and not
@@ -1203,7 +1205,12 @@ def fit(self, X, y, sample_weight=None):
             raise ValueError("Tolerance for stopping criteria must be "
                              "positive; got (tol=%r)" % self.tol)
 
-        X, y = check_X_y(X, y, accept_sparse='csr', dtype=np.float64,
+        if self.solver in ['newton-cg']:
+            _dtype = [np.float64, np.float32]
+        else:
+            _dtype = np.float64
+
+        X, y = check_X_y(X, y, accept_sparse='csr', dtype=_dtype,
                          order="C")
         check_classification_targets(y)
         self.classes_ = np.unique(y)
diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py
index 01f470a0e44cc..efb1c55783610 100644
--- a/sklearn/linear_model/tests/test_logistic.py
+++ b/sklearn/linear_model/tests/test_logistic.py
@@ -1136,3 +1136,38 @@ def test_saga_vs_liblinear():
                 liblinear.fit(X, y)
                 # Convergence for alpha=1e-3 is very slow
                 assert_array_almost_equal(saga.coef_, liblinear.coef_, 3)
+
+
+def test_dtype_match():
+    # Test that np.float32 input data is not cast to np.float64 when possible
+
+    X_ = np.array(X).astype(np.float32)
+    y_ = np.array(Y1).astype(np.float32)
+
+    for solver in ['newton-cg']:
+        for multi_class in ['ovr', 'multinomial']:
+
+            # Check type consistency
+            lr_32 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_32.fit(X_, y_)
+            assert_equal(lr_32.coef_.dtype, X_.dtype)
+
+            # Check accuracy consistency
+            lr_64 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_64.fit(X, Y1)
+            assert_almost_equal(lr_32.coef_, lr_64.coef_.astype(np.float32))
+
+
+def test_dtype_missmatch_to_profile():
+    # Test that np.float32 input data is not cast to np.float64 when possible
+
+    X_ = np.array(X).astype(np.float32)
+    y_ = np.array(Y1).astype(np.float32)
+
+    for solver in ['newton-cg']:
+        for multi_class in ['ovr', 'multinomial']:
+
+            # Check type consistency
+            lr_32 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_32.fit(X_, y_)
+            assert_equal(lr_32.coef_.dtype, X_.dtype)
diff --git a/sklearn/utils/class_weight.py b/sklearn/utils/class_weight.py
index 353f0023fc686..dea3459f3e2e9 100644
--- a/sklearn/utils/class_weight.py
+++ b/sklearn/utils/class_weight.py
@@ -41,12 +41,17 @@ def compute_class_weight(class_weight, classes, y):
     # Import error caused by circular imports.
     from ..preprocessing import LabelEncoder
 
+    if y.dtype == np.float32:
+        _dtype = np.float32
+    else:
+        _dtype = np.float64
+
     if set(y) - set(classes):
         raise ValueError("classes should include all valid labels that can "
                          "be in y")
     if class_weight is None or len(class_weight) == 0:
         # uniform class weights
-        weight = np.ones(classes.shape[0], dtype=np.float64, order='C')
+        weight = np.ones(classes.shape[0], dtype=_dtype, order='C')
     elif class_weight == 'balanced':
         # Find the weight of each class as present in y.
         le = LabelEncoder()
@@ -55,11 +60,11 @@ def compute_class_weight(class_weight, classes, y):
             raise ValueError("classes should have valid labels that are in y")
 
         recip_freq = len(y) / (len(le.classes_) *
-                               bincount(y_ind).astype(np.float64))
+                               bincount(y_ind).astype(_dtype))
         weight = recip_freq[le.transform(classes)]
     else:
         # user-defined dictionary
-        weight = np.ones(classes.shape[0], dtype=np.float64, order='C')
+        weight = np.ones(classes.shape[0], dtype=_dtype, order='C')
         if not isinstance(class_weight, dict):
             raise ValueError("class_weight must be dict, 'balanced', or None,"
                              " got: %r" % class_weight)


Code review comments:
diff_hunk: @@ -1136,3 +1136,38 @@ def test_saga_vs_liblinear():
                 liblinear.fit(X, y)
                 # Convergence for alpha=1e-3 is very slow
                 assert_array_almost_equal(saga.coef_, liblinear.coef_, 3)
+
+
+def test_dtype_match():
+    # Test that np.float32 input data is not cast to np.float64 when possible
+
+    X_ = np.array(X).astype(np.float32)
path: sklearn/linear_model/tests/test_logistic.py
text: @user1:
X_32 = ... astype(32)
X_64 = ... astype(64)

diff_hunk: @@ -41,12 +41,17 @@ def compute_class_weight(class_weight, classes, y):
     # Import error caused by circular imports.
     from ..preprocessing import LabelEncoder
 
+    if y.dtype == np.float32:
+        _dtype = np.float32
path: sklearn/utils/class_weight.py
text: @user1:
why not `_dtype=y.dtype`...

is it so you can have `y.dtype` to be `int` and `weight` to be of `float`?

diff_hunk: @@ -1136,3 +1136,38 @@ def test_saga_vs_liblinear():
                 liblinear.fit(X, y)
                 # Convergence for alpha=1e-3 is very slow
                 assert_array_almost_equal(saga.coef_, liblinear.coef_, 3)
+
+
+def test_dtype_match():
+    # Test that np.float32 input data is not cast to np.float64 when possible
+
+    X_ = np.array(X).astype(np.float32)
+    y_ = np.array(Y1).astype(np.float32)
+
+    for solver in ['newton-cg']:
+        for multi_class in ['ovr', 'multinomial']:
+
+            # Check type consistency
+            lr_32 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_32.fit(X_, y_)
+            assert_equal(lr_32.coef_.dtype, X_.dtype)
+
+            # Check accuracy consistency
+            lr_64 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_64.fit(X, Y1)
+            assert_almost_equal(lr_32.coef_, lr_64.coef_.astype(np.float32))
+
+
+def test_dtype_missmatch_to_profile():
+    # Test that np.float32 input data is not cast to np.float64 when possible
+
path: sklearn/linear_model/tests/test_logistic.py
text: @user1:
and this newline too

diff_hunk: @@ -1136,3 +1136,38 @@ def test_saga_vs_liblinear():
                 liblinear.fit(X, y)
                 # Convergence for alpha=1e-3 is very slow
                 assert_array_almost_equal(saga.coef_, liblinear.coef_, 3)
+
+
+def test_dtype_match():
+    # Test that np.float32 input data is not cast to np.float64 when possible
+
+    X_ = np.array(X).astype(np.float32)
+    y_ = np.array(Y1).astype(np.float32)
+
+    for solver in ['newton-cg']:
+        for multi_class in ['ovr', 'multinomial']:
+
+            # Check type consistency
+            lr_32 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_32.fit(X_, y_)
+            assert_equal(lr_32.coef_.dtype, X_.dtype)
+
+            # Check accuracy consistency
+            lr_64 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_64.fit(X, Y1)
path: sklearn/linear_model/tests/test_logistic.py
text: @user1:
Can you ensure (maybe using `astype`?) `X, Y1` are of `float64` before this test? (If in future it is changed, this test will still pass)

diff_hunk: @@ -1136,3 +1136,38 @@ def test_saga_vs_liblinear():
                 liblinear.fit(X, y)
                 # Convergence for alpha=1e-3 is very slow
                 assert_array_almost_equal(saga.coef_, liblinear.coef_, 3)
+
+
+def test_dtype_match():
+    # Test that np.float32 input data is not cast to np.float64 when possible
+
+    X_ = np.array(X).astype(np.float32)
+    y_ = np.array(Y1).astype(np.float32)
+
+    for solver in ['newton-cg']:
+        for multi_class in ['ovr', 'multinomial']:
+
+            # Check type consistency
+            lr_32 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_32.fit(X_, y_)
+            assert_equal(lr_32.coef_.dtype, X_.dtype)
+
+            # Check accuracy consistency
+            lr_64 = LogisticRegression(solver=solver, multi_class=multi_class)
+            lr_64.fit(X, Y1)
+            assert_almost_equal(lr_32.coef_, lr_64.coef_.astype(np.float32))
+
+
+def test_dtype_missmatch_to_profile():
path: sklearn/linear_model/tests/test_logistic.py
text: @user1:
This test can be removed...