Repository: scikit-learn/scikit-learn
Issue category: Feature Request

Issue:
GradientBoostingClassifier.fit accepts sparse X, but .predict does not
I have a sparse dataset that is too large for main memory if I call `X.todense()`.  If I understand correctly, `GradientBoostingClassifier.fit` will accept my sparse `X`, but it is not currently possible to use `GradientBoostingClassifier.predict` on the results.  It would be great if that were not the case.

Here is a minimal example of the issue:

```
from scipy import sparse
from sklearn.datasets.samples_generator import make_classification
from sklearn.ensemble import GradientBoostingClassifier

X, y = make_classification(n_samples=20, n_features=5, random_state=0)
X_sp = sparse.coo_matrix(X)

clf = GradientBoostingClassifier()
clf.fit(X,y)
clf.predict(X)  # works

clf.fit(X_sp, y)  # works
clf.predict(X_sp)  # fails with TypeError: A sparse matrix was passed, but dense data is required.
```



Hints:
Confirmed. I think small rework of https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/ensemble/_gradient_boosting.pyx#L39 is needed.
Looks like function _predict_regression_tree_inplace_fast was written to somehow optimize prediction speed. But i see some problems in it, its main loop https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/ensemble/_gradient_boosting.pyx#L92 differs from the main loop here https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/tree/_tree.pyx#L784, and additionally this function doesn't  work with sparse matrices.

Maybe it's better to rely on predict method of Tree class, but some performance tests is needed.


Patch:
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
index b5f10da91d28f..a4b775ec66d0a 100644
--- a/doc/whats_new.rst
+++ b/doc/whats_new.rst
@@ -36,6 +36,11 @@ Enhancements
      (`#6288 <https://github.com/scikit-learn/scikit-learn/pull/6288>`_) by
      `Jake VanderPlas`_.
 
+   - :class:`ensemble.GradientBoostingClassifier` and :class:`ensemble.GradientBoostingRegressor`
+     now support sparse input for prediction.
+     (`#6101 <https://github.com/scikit-learn/scikit-learn/pull/6101>`_)
+     By `Ibraim Ganiev`_.
+
 Bug fixes
 .........
 
diff --git a/sklearn/ensemble/_gradient_boosting.pyx b/sklearn/ensemble/_gradient_boosting.pyx
index 9e6e9f6d29c0e..71371f5c24a48 100644
--- a/sklearn/ensemble/_gradient_boosting.pyx
+++ b/sklearn/ensemble/_gradient_boosting.pyx
@@ -8,13 +8,22 @@
 
 cimport cython
 
+from libc.stdlib cimport free
+from libc.string cimport memset
+
 import numpy as np
 cimport numpy as np
 np.import_array()
 
+from scipy.sparse import issparse
+from scipy.sparse import csr_matrix
+
 from sklearn.tree._tree cimport Node
 from sklearn.tree._tree cimport Tree
-
+from sklearn.tree._tree cimport DTYPE_t
+from sklearn.tree._tree cimport SIZE_t
+from sklearn.tree._tree cimport INT32_t
+from sklearn.tree._utils cimport safe_realloc
 
 ctypedef np.int32_t int32
 ctypedef np.float64_t float64
@@ -27,24 +36,19 @@ from numpy import bool as np_bool
 from numpy import float32 as np_float32
 from numpy import float64 as np_float64
 
-# Define a datatype for the data array
-DTYPE = np.float32
-ctypedef np.float32_t DTYPE_t
-ctypedef np.npy_intp SIZE_t
-
 
 # constant to mark tree leafs
-cdef int LEAF = -1
-
-cdef void _predict_regression_tree_inplace_fast(DTYPE_t *X,
-                                                Node* root_node,
-                                                double *value,
-                                                double scale,
-                                                Py_ssize_t k,
-                                                Py_ssize_t K,
-                                                Py_ssize_t n_samples,
-                                                Py_ssize_t n_features,
-                                                float64 *out):
+cdef SIZE_t TREE_LEAF = -1
+
+cdef void _predict_regression_tree_inplace_fast_dense(DTYPE_t *X,
+                                                      Node* root_node,
+                                                      double *value,
+                                                      double scale,
+                                                      Py_ssize_t k,
+                                                      Py_ssize_t K,
+                                                      Py_ssize_t n_samples,
+                                                      Py_ssize_t n_features,
+                                                      float64 *out):
     """Predicts output for regression tree and stores it in ``out[i, k]``.
 
     This function operates directly on the data arrays of the tree
@@ -84,22 +88,108 @@ cdef void _predict_regression_tree_inplace_fast(DTYPE_t *X,
         shape ``(n_samples, K)``.
     """
     cdef Py_ssize_t i
-    cdef int32 node_id
     cdef Node *node
     for i in range(n_samples):
         node = root_node
         # While node not a leaf
-        while node.left_child != -1 and node.right_child != -1:
+        while node.left_child != TREE_LEAF:
             if X[i * n_features + node.feature] <= node.threshold:
                 node = root_node + node.left_child
             else:
                 node = root_node + node.right_child
         out[i * K + k] += scale * value[node - root_node]
 
+def _predict_regression_tree_stages_sparse(np.ndarray[object, ndim=2] estimators,
+                                           object X, double scale,
+                                           np.ndarray[float64, ndim=2] out):
+    """Predicts output for regression tree inplace and adds scaled value to ``out[i, k]``.
+
+    The function assumes that the ndarray that wraps ``X`` is csr_matrix.
+    """
+    cdef DTYPE_t* X_data = <DTYPE_t*>(<np.ndarray> X.data).data
+    cdef INT32_t* X_indices = <INT32_t*>(<np.ndarray> X.indices).data
+    cdef INT32_t* X_indptr = <INT32_t*>(<np.ndarray> X.indptr).data
+
+    cdef SIZE_t n_samples = X.shape[0]
+    cdef SIZE_t n_features = X.shape[1]
+    cdef SIZE_t n_stages = estimators.shape[0]
+    cdef SIZE_t n_outputs = estimators.shape[1]
+
+    # Initialize output
+    cdef float64* out_ptr = <float64*> out.data
+
+    # Indices and temporary variables
+    cdef SIZE_t sample_i
+    cdef SIZE_t feature_i
+    cdef SIZE_t stage_i
+    cdef SIZE_t output_i
+    cdef Node *root_node = NULL
+    cdef Node *node = NULL
+    cdef double *value = NULL
+
+    cdef Tree tree
+    cdef Node** nodes = NULL
+    cdef double** values = NULL
+    safe_realloc(&nodes, n_stages * n_outputs)
+    safe_realloc(&values, n_stages * n_outputs)
+    for stage_i in range(n_stages):
+        for output_i in range(n_outputs):
+            tree = estimators[stage_i, output_i].tree_
+            nodes[stage_i * n_outputs + output_i] = tree.nodes
+            values[stage_i * n_outputs + output_i] = tree.value
+
+    # Initialize auxiliary data-structure
+    cdef DTYPE_t feature_value = 0.
+    cdef DTYPE_t* X_sample = NULL
+
+    # feature_to_sample as a data structure records the last seen sample
+    # for each feature; functionally, it is an efficient way to identify
+    # which features are nonzero in the present sample.
+    cdef SIZE_t* feature_to_sample = NULL
+
+    safe_realloc(&X_sample, n_features)
+    safe_realloc(&feature_to_sample, n_features)
+
+    memset(feature_to_sample, -1, n_features * sizeof(SIZE_t))
+
+    # Cycle through all samples
+    for sample_i in range(n_samples):
+        for feature_i in range(X_indptr[sample_i], X_indptr[sample_i + 1]):
+            feature_to_sample[X_indices[feature_i]] = sample_i
+            X_sample[X_indices[feature_i]] = X_data[feature_i]
+
+        # Cycle through all stages
+        for stage_i in range(n_stages):
+            # Cycle through all trees
+            for output_i in range(n_outputs):
+                root_node = nodes[stage_i * n_outputs + output_i]
+                value = values[stage_i * n_outputs + output_i]
+                node = root_node
+
+                # While node not a leaf
+                while node.left_child != TREE_LEAF:
+                    # ... and node.right_child != TREE_LEAF:
+                    if feature_to_sample[node.feature] == sample_i:
+                        feature_value = X_sample[node.feature]
+                    else:
+                        feature_value = 0.
+
+                    if feature_value <= node.threshold:
+                        node = root_node + node.left_child
+                    else:
+                        node = root_node + node.right_child
+                out_ptr[sample_i * n_outputs + output_i] += (scale
+                    * value[node - root_node])
+
+    # Free auxiliary arrays
+    free(X_sample)
+    free(feature_to_sample)
+    free(nodes)
+    free(values)
+
 
-@cython.nonecheck(False)
 def predict_stages(np.ndarray[object, ndim=2] estimators,
-                   np.ndarray[DTYPE_t, ndim=2, mode='c'] X, double scale,
+                   object X, double scale,
                    np.ndarray[float64, ndim=2] out):
     """Add predictions of ``estimators`` to ``out``.
 
@@ -112,25 +202,31 @@ def predict_stages(np.ndarray[object, ndim=2] estimators,
     cdef Py_ssize_t K = estimators.shape[1]
     cdef Tree tree
 
-    for i in range(n_estimators):
-        for k in range(K):
-            tree = estimators[i, k].tree_
+    if issparse(X):
+        _predict_regression_tree_stages_sparse(estimators, X, scale, out)
+    else:
+        if not isinstance(X, np.ndarray):
+            raise ValueError("X should be in np.ndarray or csr_matrix format,"
+                             "got %s" % type(X))
+
+        for i in range(n_estimators):
+            for k in range(K):
+                tree = estimators[i, k].tree_
 
-            # avoid buffer validation by casting to ndarray
-            # and get data pointer
-            # need brackets because of casting operator priority
-            _predict_regression_tree_inplace_fast(
-                <DTYPE_t*> X.data,
-                tree.nodes, tree.value,
-                scale, k, K, X.shape[0], X.shape[1],
-                <float64 *> (<np.ndarray> out).data)
-            ## out += scale * tree.predict(X).reshape((X.shape[0], 1))
+                # avoid buffer validation by casting to ndarray
+                # and get data pointer
+                # need brackets because of casting operator priority
+                _predict_regression_tree_inplace_fast_dense(
+                    <DTYPE_t*> (<np.ndarray> X).data,
+                    tree.nodes, tree.value,
+                    scale, k, K, X.shape[0], X.shape[1],
+                    <float64 *> (<np.ndarray> out).data)
+                ## out += scale * tree.predict(X).reshape((X.shape[0], 1))
 
 
-@cython.nonecheck(False)
 def predict_stage(np.ndarray[object, ndim=2] estimators,
                   int stage,
-                  np.ndarray[DTYPE_t, ndim=2] X, double scale,
+                  object X, double scale,
                   np.ndarray[float64, ndim=2] out):
     """Add predictions of ``estimators[stage]`` to ``out``.
 
@@ -216,7 +312,7 @@ cpdef _partial_dependence_tree(Tree tree, DTYPE_t[:, ::1] X,
             stack_size -= 1
             current_node = node_stack[stack_size]
 
-            if current_node.left_child == LEAF:
+            if current_node.left_child == TREE_LEAF:
                 out[i] += weight_stack[stack_size] * value[current_node - root_node] * \
                           learn_rate
                 total_weight += weight_stack[stack_size]
diff --git a/sklearn/ensemble/gradient_boosting.py b/sklearn/ensemble/gradient_boosting.py
index c9a36aac1bd99..fcfeb45a09157 100644
--- a/sklearn/ensemble/gradient_boosting.py
+++ b/sklearn/ensemble/gradient_boosting.py
@@ -1158,8 +1158,10 @@ def _staged_decision_function(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
@@ -1169,7 +1171,7 @@ def _staged_decision_function(self, X):
             Regression and binary classification are special cases with
             ``k == 1``, otherwise ``k==n_classes``.
         """
-        X = check_array(X, dtype=DTYPE, order="C")
+        X = check_array(X, dtype=DTYPE, order="C",  accept_sparse='csr')
         score = self._init_decision_function(X)
         for i in range(self.estimators_.shape[0]):
             predict_stage(self.estimators_, i, X, self.learning_rate, score)
@@ -1479,8 +1481,10 @@ def decision_function(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
@@ -1490,7 +1494,7 @@ def decision_function(self, X):
             Regression and binary classification produce an array of shape
             [n_samples].
         """
-        X = check_array(X, dtype=DTYPE, order="C")
+        X = check_array(X, dtype=DTYPE, order="C",  accept_sparse='csr')
         score = self._decision_function(X)
         if score.shape[1] == 1:
             return score.ravel()
@@ -1504,8 +1508,10 @@ def staged_decision_function(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
@@ -1524,8 +1530,10 @@ def predict(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
@@ -1544,8 +1552,10 @@ def staged_predict(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
@@ -1561,8 +1571,10 @@ def predict_proba(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Raises
         ------
@@ -1589,8 +1601,10 @@ def predict_log_proba(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Raises
         ------
@@ -1614,8 +1628,10 @@ def staged_predict_proba(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
@@ -1840,15 +1856,17 @@ def predict(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
         y : array of shape = [n_samples]
             The predicted values.
         """
-        X = check_array(X, dtype=DTYPE, order="C")
+        X = check_array(X, dtype=DTYPE, order="C",  accept_sparse='csr')
         return self._decision_function(X).ravel()
 
     def staged_predict(self, X):
@@ -1859,8 +1877,10 @@ def staged_predict(self, X):
 
         Parameters
         ----------
-        X : array-like of shape = [n_samples, n_features]
-            The input samples.
+        X : array-like or sparse matrix, shape = [n_samples, n_features]
+            The input samples. Internally, it will be converted to
+            ``dtype=np.float32`` and if a sparse matrix is provided
+            to a sparse ``csr_matrix``.
 
         Returns
         -------
diff --git a/sklearn/tree/_tree.pyx b/sklearn/tree/_tree.pyx
index 4e8160f7dc780..4f11cff0569e6 100644
--- a/sklearn/tree/_tree.pyx
+++ b/sklearn/tree/_tree.pyx
@@ -839,8 +839,8 @@ cdef class Tree:
         # which features are nonzero in the present sample.
         cdef SIZE_t* feature_to_sample = NULL
 
-        safe_realloc(&X_sample, n_features * sizeof(DTYPE_t))
-        safe_realloc(&feature_to_sample, n_features * sizeof(SIZE_t))
+        safe_realloc(&X_sample, n_features)
+        safe_realloc(&feature_to_sample, n_features)
 
         with nogil:
             memset(feature_to_sample, -1, n_features * sizeof(SIZE_t))
@@ -985,8 +985,8 @@ cdef class Tree:
         # which features are nonzero in the present sample.
         cdef SIZE_t* feature_to_sample = NULL
 
-        safe_realloc(&X_sample, n_features * sizeof(DTYPE_t))
-        safe_realloc(&feature_to_sample, n_features * sizeof(SIZE_t))
+        safe_realloc(&X_sample, n_features)
+        safe_realloc(&feature_to_sample, n_features)
 
         with nogil:
             memset(feature_to_sample, -1, n_features * sizeof(SIZE_t))
diff --git a/sklearn/tree/_utils.pxd b/sklearn/tree/_utils.pxd
index d11880908c318..fce3abcb734db 100644
--- a/sklearn/tree/_utils.pxd
+++ b/sklearn/tree/_utils.pxd
@@ -10,6 +10,7 @@
 
 import numpy as np
 cimport numpy as np
+from _tree cimport Node 
 
 ctypedef np.npy_float32 DTYPE_t          # Type of X
 ctypedef np.npy_float64 DOUBLE_t         # Type of y, sample_weight
@@ -36,6 +37,9 @@ ctypedef fused realloc_ptr:
     (unsigned char*)
     (WeightedPQueueRecord*)
     (DOUBLE_t*)
+    (DOUBLE_t**)
+    (Node*)
+    (Node**)
 
 cdef realloc_ptr safe_realloc(realloc_ptr* p, size_t nelems) except *
 


Test patch:
diff --git a/sklearn/ensemble/tests/test_gradient_boosting.py b/sklearn/ensemble/tests/test_gradient_boosting.py
index 1ebf82fa3ae44..634bc259a1167 100644
--- a/sklearn/ensemble/tests/test_gradient_boosting.py
+++ b/sklearn/ensemble/tests/test_gradient_boosting.py
@@ -1050,6 +1050,9 @@ def check_sparse_input(EstimatorClass, X, X_sparse, y):
     assert_array_almost_equal(sparse.feature_importances_,
                               auto.feature_importances_)
 
+    assert_array_almost_equal(sparse.predict(X_sparse), dense.predict(X))
+    assert_array_almost_equal(dense.predict(X_sparse), sparse.predict(X))
+
     if isinstance(EstimatorClass, GradientBoostingClassifier):
         assert_array_almost_equal(sparse.predict_proba(X),
                                   dense.predict_proba(X))
@@ -1061,6 +1064,14 @@ def check_sparse_input(EstimatorClass, X, X_sparse, y):
         assert_array_almost_equal(sparse.predict_log_proba(X),
                                   auto.predict_log_proba(X))
 
+        assert_array_almost_equal(sparse.decision_function(X_sparse),
+                                  sparse.decision_function(X))
+        assert_array_almost_equal(dense.decision_function(X_sparse),
+                                  sparse.decision_function(X))
+
+        assert_array_almost_equal(
+            np.array(sparse.staged_decision_function(X_sparse)),
+            np.array(sparse.staged_decision_function(X)))
 
 @skip_if_32bit
 def test_sparse_input():