Repository: scikit-learn/scikit-learn
Pull request: ENH  Optimize runtime for IsolationForest
<!--
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 #19275
This is a follow-up PR of #23149.

#### What does this implement/fix? Explain your changes.
As discussed in the [comment](https://github.com/scikit-learn/scikit-learn/pull/23149#issuecomment-1107512467), the indexing operation in `_bagging.py: _parallel_build_estimators` is quite expensive, so when `bootstrap_features` is set to **False** and `max_features_` is equal to number of features, we can skip indexing to get better performance.

#### Benchmark result
Code used for profiling:
```python
from sklearn.datasets import make_classification
from scipy.sparse import csc_matrix, csr_matrix
from sklearn.ensemble import IsolationForest

X, y = make_classification(n_samples=50000, n_features=1000)
X = csc_matrix(X)
X.sort_indices()
IsolationForest(n_estimators=10, max_samples=256, n_jobs=1).fit(X)
```
**Sparse Input:**
Before (4.97s)
![image](https://user-images.githubusercontent.com/16646940/164909922-70a89054-a81e-4931-acca-46e4779d8eba.png)


After (0.241s)
![image](https://user-images.githubusercontent.com/16646940/164909946-4c070410-0176-43ed-a94b-cd9e08afb1af.png)


**Dense Input:**
Before (7.99s)
![image](https://user-images.githubusercontent.com/16646940/164909872-3cec6573-efd9-4baf-83a2-d47caaec7e14.png)

After (5.19s)
![image](https://user-images.githubusercontent.com/16646940/164909792-52b711b6-0280-44cc-8566-036ba0043162.png)


#### Any other comments?
`boostrap_features` is set to **False** already for IsolationForest, but users can still set `max_features` freely, which will cause a worse runtime when it's not 1.0. 

Since in IsolationForest, we're randomly picking a single feature to split, shall we deprecate the `max_features` argument and always set it to 1.0 for better runtime ?

<!--
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:
IsolationForest extremely slow with large number of columns having discrete values
The following example takes an unreasonable amount of time to run:
```python
import numpy as np
from sklearn.ensemble import IsolationForest
from sklearn.datasets import fetch_rcv1
from scipy.sparse import csc_matrix, csr_matrix

X, y = fetch_rcv1(return_X_y=True)
X = csc_matrix(X)
X.sort_indices()
iso = IsolationForest(n_estimators=100, max_samples=256).fit(X)
```

In theory it should be very fast, since each sub-sample it takes is a small sparse matrix in which most columns will have only zeros.

Using `n_jobs>1` also makes it use a very unreasonable amount of memory for some reason.

If the input is passed as dense, the running time still looks worse than it should. I guess the issue from a quick glance at the code is that it doesn't remember which columns are already not possible to split in a given node:
```python
X_sample = csr_matrix(X)[:1000,:]
X_sample = X_sample.toarray()
iso = IsolationForest(n_estimators=100, max_samples=256).fit(X_sample)
```

Hints from issue discussion:
Here is an edited version of your reproducer to get more precise measurements:

```python
import sys
import os
import pickle
import numpy as np
from time import perf_counter
from sklearn.ensemble import IsolationForest
from sklearn.datasets import fetch_rcv1
from scipy.sparse import csc_matrix, csr_matrix

X, y = fetch_rcv1(return_X_y=True)
X = csc_matrix(X)
X.sort_indices()

data_size = X.data.nbytes + X.indices.nbytes + X.indptr.nbytes
print(f"data size: {data_size / 1e6:.1f} MB")

n_jobs = int(sys.argv[1])
print(f"Running IsolationForest with n_jobs={n_jobs}...")
tic = perf_counter()
iso = IsolationForest(n_estimators=100, max_samples=256, n_jobs=n_jobs).fit(X)
print(f"duration: {perf_counter() - tic:.1f} s")


fname = "/tmp/tmp_model.pkl"
with open(fname, "wb") as f:
    pickle.dump(iso, f)

model_size = os.stat(fname).st_size
print(f"final model size: {model_size / 1e6:.1f} MB")
```

I can then use [memory_profiler](https://pypi.org/project/memory-profiler/) to monitory memory usage over time. Here is what I observe:

### n_jobs = 1

```
(dev) ogrisel@arm64-apple-darwin20 /tmp % mprof run iso.py 1
mprof: Sampling memory every 0.1s
running new process
data size: 731.2 MB
Running IsolationForest with n_jobs=1...
duration: 24.1 s
final model size: 38.0 MB
(dev) ogrisel@arm64-apple-darwin20 /tmp % mprof plot
```

![isoforest_rcv1_csc_n_jobs_1](https://user-images.githubusercontent.com/89061/105971909-5522d980-608b-11eb-99a3-a487d7a8d398.png)


### n_jobs = 4

```

(dev) ogrisel@arm64-apple-darwin20 /tmp % mprof run iso.py 4
mprof: Sampling memory every 0.1s
running new process
data size: 731.2 MB
Running IsolationForest with n_jobs=4...
duration: 18.3 s
final model size: 38.0 MB
(dev) ogrisel@arm64-apple-darwin20 /tmp % mprof plot 
```
![isoforest_rcv1_csc_n_jobs_4](https://user-images.githubusercontent.com/89061/105971928-594ef700-608b-11eb-8807-962f0161b940.png)

so a peak memory usage of 5 GB which indeed seems quite large but maybe expected: the resampling from the RF bagging performed by each model in parallel is causing the dataset to be replicated once per parallel tree. As the original data is already  731.2 MB, 5 GB in total memory usage does not sound too catastrophic to me: 1 x 731 (original data) + 4 x 731 (concurrent resampled copies) ~ 3GB. Assuming that the Python GC lags a bit, getting a few extra GB.

Also I am not sure what is going one when you row-wise resample a scipy CSC matrix but I suspect there might be intermediate data structures being materialized.

As for execution time, I don't really know if ~25s is abnormally slow or not to fit 100 isolation trees on this sparse data.

> I guess the issue from a quick glance at the code is that it doesn't remember which columns are already not possible to split in a given node:

That sounds like an interesting optimization avenue, indeed. A PR would be appreciated, I think :)
Actually my analysis above is probably wrong: `bootstrap` is `False` by default for Isolation Forests (contrary to RF). So there is indeed something fish going on. This deserves further investigation.
Hum I just noticed the `max_samples=256`: indeed each tree should be trained on a very small subsample of the dataset so the memory usage seems far to high for this value, especially in parallel mode. Using memory_profiler with the line by line report could help identify what's going on.

BTW, since there is row-wise resampling, converting the original data to CSC is probably useless and even counterproductive. Keeping the CSR representation is much more efficient for row-wise resampling. But this is not the main cause of the performance problem here, I think.
> BTW, since there is row-wise resampling, converting the original data to CSC is probably useless and even counterproductive. Keeping the CSR representation is much more efficient for row-wise resampling.

I am not so sure about that either. Will need line by line profiling to understand what's happening.
To provide some context about the running time, fitting the same model takes <200ms in a different implementation I have ([here](https://www.github.com/david-cortes/isotree)) which uses an array of sub-sample row indices instead of copying the data. It's a very big difference.
Hi @CheriseCodes the labels 'module: something' identify the scikit-learn module the issue is referring to.
This is specified in the [Triaging](https://scikit-learn.org/stable/developers/bug_triaging.html?highlight=triaging#a-typical-workflow-for-triaging-issues) page. This kind of labeling is useful to track which module is the more used or the less troubleshooted by maintainers.
Thanks @cmarmo. I appreciate the resource link. (for anyone confused, I deleted my comment from earlier)
After a little profiling, it seems like the problem has something to do with `scipy.sparse._sparsetools.csr_row_index` (although I haven't figured out exactly where it's used)

#### IsoForest with 10 estimators:
![image](https://user-images.githubusercontent.com/16646940/113248383-70aa9a80-92ef-11eb-8287-87e05fc7ec45.png)
`csr_row_index` get called 10 times, and takes up 40% of the run time

#### RandomForest with 10 estimators:
![image](https://user-images.githubusercontent.com/16646940/113236388-54026880-92d7-11eb-84f5-db71019d4f92.png)
only get called once
I think one of the reasons is that DecisionTree has a `check_input` which is set to True by default. `RandomForest` is derived from `BaseForest` which uses `_parallel_build_trees` to construct underlying estimators, and it sets `check_input` to False during the process:
```python
# sklearn.ensemble._forest.py
def _parallel_build_trees(tree, forest, X, y, sample_weight, tree_idx, n_trees,
                          verbose=0, class_weight=None,
                          n_samples_bootstrap=None):
    """
    Private function used to fit a single tree in parallel."""
   ......
        tree.fit(X, y, sample_weight=curr_sample_weight, check_input=False)
    else:
        tree.fit(X, y, sample_weight=sample_weight, check_input=False)
    return tree
```

meanwhile `IsoForest` is derived from `BaseBagging` which uses `_parallel_build_estimators`, and it does not set `check_input` to False, so we will be checking X repetitively:
```python
# sklearn.ensemble._bagging.py
def _parallel_build_estimators(n_estimators, ensemble, X, y, sample_weight,
                               seeds, total_n_estimators, verbose):
    """Private function used to build a batch of estimators within a job."""
   ......
    for i in range(n_estimators):
        random_state = seeds[i]
        estimator = ensemble._make_estimator(append=False,
                                             random_state=random_state)

        # Draw random feature, sample indices
        features, indices = _generate_bagging_indices(random_state,
                                                      bootstrap_features,
                                                      bootstrap, n_features,
                                                      n_samples, max_features,
                                                      max_samples)
         .......
            estimator.fit(X[:, features], y, sample_weight=curr_sample_weight)

        else:
            estimator.fit((X[indices])[:, features], y[indices])

        estimators.append(estimator)
        estimators_features.append(features)

    return estimators, estimators_features
```

The slicing operation `X[:, features]` is also pretty slow, fitting 10 estimators with n_jobs=1 takes 10s, adding `check_input=False` takes around 6s, and removing the slicing operation (incorrectly) takes under 1s. 

Maybe we should change `IsolationForest` to derive from `BaseForest`? 
I had a chance to look at the code some this week. I think @ogrisel is on to something regarding the materializing of alternative sparse matrices. Specifically, if you call fit with a row major (csr) matrix, when IsolationForest fit method there is a call to _validate_data (defined in base.py) that appears to only accept column major sparse data (csc). see below:

```
#_iforest.py:236
def fit(self, X, y=None, sample_weight=None):       
  X = self._validate_data(X, accept_sparse=["csc"])
  if issparse(X):
      # Pre-sort indices to avoid that each individual tree of the
      # ensemble sorts the indices.
      X.sort_indices()

```

This will cause the code to convert the data from csr to csc. I tried adding "csr" to the list to understand the performance impact and that made things much worse performance wise. This is because BaseDecisionTree only accepts csc type sparse arrays. The code in iForest makes the conversion happen once. Accepting a csr in iforest causes repeated conversion on the creation of each estimator within BaseDecisionTree. (see below). This makes me think that the documentation should advise against using csr matrices with this class.

```   
#_classes.py:150
def fit(
        self, X, y, sample_weight=None, check_input=True, X_idx_sorted="deprecated"
    ):

        random_state = check_random_state(self.random_state)

        if self.ccp_alpha < 0.0:
            raise ValueError("ccp_alpha must be greater than or equal to 0")

        if check_input:
            # Need to validate separately here.
            # We can't pass multi_ouput=True because that would allow y to be
            # csr.
            check_X_params = dict(dtype=DTYPE, accept_sparse="csc")
```

Last thought. If the data being passed was really only 256 rows (per the opening message of this thread), that shouldn't matter much. In my testing though, it seems that sample weights are being used at times and the entire large dataset being converted. Something else to think on.

Commit message: fix merge conflict

Patch to review:
diff --git a/doc/whats_new/v1.2.rst b/doc/whats_new/v1.2.rst
index 5f45cdf6a7168..e409461bccc77 100644
--- a/doc/whats_new/v1.2.rst
+++ b/doc/whats_new/v1.2.rst
@@ -44,6 +44,12 @@ Changelog
 - |Enhancement| :class:`cluster.Birch` now preserves dtype for `numpy.float32`
   inputs. :pr:`22968` by `Meekail Zain <micky774>`.
 
+:mod:`sklearn.ensemble`
+.......................
+
+- |Efficiency| Improve runtime performance of :class:`ensemble.IsolationForest`
+  by avoiding data copies. :pr:`23252` by :user:`Zhehao Liu <MaxwellLZH>`.
+
 Code and Documentation Contributors
 -----------------------------------
 
diff --git a/sklearn/ensemble/_bagging.py b/sklearn/ensemble/_bagging.py
index adac7e063191f..9c8faa783c788 100644
--- a/sklearn/ensemble/_bagging.py
+++ b/sklearn/ensemble/_bagging.py
@@ -88,6 +88,8 @@ def _parallel_build_estimators(
     bootstrap_features = ensemble.bootstrap_features
     support_sample_weight = has_fit_parameter(ensemble.base_estimator_, "sample_weight")
     has_check_input = has_fit_parameter(ensemble.base_estimator_, "check_input")
+    requires_feature_indexing = bootstrap_features or max_features != n_features
+
     if not support_sample_weight and sample_weight is not None:
         raise ValueError("The base estimator doesn't support sample weight")
 
@@ -135,10 +137,11 @@ def _parallel_build_estimators(
                 not_indices_mask = ~indices_to_mask(indices, n_samples)
                 curr_sample_weight[not_indices_mask] = 0
 
-            estimator_fit(X[:, features], y, sample_weight=curr_sample_weight)
-
+            X_ = X[:, features] if requires_feature_indexing else X
+            estimator_fit(X_, y, sample_weight=curr_sample_weight)
         else:
-            estimator_fit(X[indices][:, features], y[indices])
+            X_ = X[indices][:, features] if requires_feature_indexing else X[indices]
+            estimator_fit(X_, y[indices])
 
         estimators.append(estimator)
         estimators_features.append(features)
diff --git a/sklearn/ensemble/_iforest.py b/sklearn/ensemble/_iforest.py
index 4be74d2873b9e..e125e177bdcae 100644
--- a/sklearn/ensemble/_iforest.py
+++ b/sklearn/ensemble/_iforest.py
@@ -8,6 +8,7 @@
 from warnings import warn
 
 from ..tree import ExtraTreeRegressor
+from ..tree._tree import DTYPE as tree_dtype
 from ..utils import (
     check_random_state,
     check_array,
@@ -254,7 +255,7 @@ def fit(self, X, y=None, sample_weight=None):
         self : object
             Fitted estimator.
         """
-        X = self._validate_data(X, accept_sparse=["csc"])
+        X = self._validate_data(X, accept_sparse=["csc"], dtype=tree_dtype)
         if issparse(X):
             # Pre-sort indices to avoid that each individual tree of the
             # ensemble sorts the indices.
diff --git a/sklearn/ensemble/tests/test_iforest.py b/sklearn/ensemble/tests/test_iforest.py
index 55046527b13c3..976c552e78890 100644
--- a/sklearn/ensemble/tests/test_iforest.py
+++ b/sklearn/ensemble/tests/test_iforest.py
@@ -20,7 +20,7 @@
 from sklearn.ensemble import IsolationForest
 from sklearn.ensemble._iforest import _average_path_length
 from sklearn.model_selection import train_test_split
-from sklearn.datasets import load_diabetes, load_iris
+from sklearn.datasets import load_diabetes, load_iris, make_classification
 from sklearn.utils import check_random_state
 from sklearn.metrics import roc_auc_score
 
@@ -347,3 +347,12 @@ def test_n_features_deprecation():
 
     with pytest.warns(FutureWarning, match="`n_features_` was deprecated"):
         est.n_features_
+
+
+def test_iforest_with_n_jobs():
+    """This is a non-regression test for PR #23252"""
+    X, _ = make_classification(n_samples=50000, n_features=100)
+    X = csc_matrix(X)
+    X.sort_indices()
+    for n_jobs in [2, 3, 4]:
+        IsolationForest(n_estimators=10, max_samples=256, n_jobs=n_jobs).fit(X)


Code review comments:
diff_hunk: @@ -347,3 +347,12 @@ def test_n_features_deprecation():
 
     with pytest.warns(FutureWarning, match="`n_features_` was deprecated"):
         est.n_features_
+
+
+def test_iforest_with_n_jobs():
+    """This is a non-regression test for PR #23252"""
+    X, _ = make_classification(n_samples=50000, n_features=100)
+    X = csc_matrix(X)
+    X.sort_indices()
+    for n_jobs in [2, 3, 4]:
+        IsolationForest(n_estimators=10, max_samples=256, n_jobs=n_jobs).fit(X)
path: sklearn/ensemble/tests/test_iforest.py
text: @user1:
I think testing for `n_jobs=2` is enough for a non-regression test.

```suggestion
    IsolationForest(n_estimators=10, max_samples=256, n_jobs=2).fit(X)
```

diff_hunk: @@ -347,3 +347,12 @@ def test_n_features_deprecation():
 
     with pytest.warns(FutureWarning, match="`n_features_` was deprecated"):
         est.n_features_
+
+
+def test_iforest_with_n_jobs():
+    """This is a non-regression test for PR #23252"""
+    X, _ = make_classification(n_samples=50000, n_features=100)
+    X = csc_matrix(X)
+    X.sort_indices()
path: sklearn/ensemble/tests/test_iforest.py
text: @user1:
Do not need to sort here.

```suggestion
```

diff_hunk: @@ -347,3 +347,12 @@ def test_n_features_deprecation():
 
     with pytest.warns(FutureWarning, match="`n_features_` was deprecated"):
         est.n_features_
+
+
+def test_iforest_with_n_jobs():
+    """This is a non-regression test for PR #23252"""
+    X, _ = make_classification(n_samples=50000, n_features=100)
path: sklearn/ensemble/tests/test_iforest.py
text: @user1:
For me, the bug segfault was triggered around here:

```suggestion
    X, _ = make_classification(n_samples=85_000, n_features=100, random_state=0)
```

I think the bug triggers when joblib's [Parallel](https://joblib.readthedocs.io/en/latest/generated/joblib.Parallel.html) uses memmapped, which is triggered when the data is big enough.