Repository: scikit-learn/scikit-learn
Pull request: [MRG] Adds Plotting API to Partial Dependence
<!--
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#pull-request-checklist
-->

#### 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 https://github.com/scikit-learn/scikit-learn/issues/14596


#### What does this implement/fix? Explain your changes.
1. Adds `PartialDependenceDisplay` for `plot_partial_dependence`.
2. Deprecates the `fig` option. The new API can take a single axes and place all the plots in that ax.
3. If a list of axes are given, it must match the number of partial dependence plots. If it does, the partial dependence plots will be place into those axes.
4. The example, `plot_partial_dependence_visualization_api.py`, is used to demonstrate how `plot_partial_dependence` how to use this axes based API.

#### Any other comments?
This PR removes the custom padding between subplots because the padding depends on the size of the figure. When using matplotlib's constrained_layout:

```py
fig, ax = plt.subplots(constrained_layout=True)
disp = plot_partial_dependence(..., ax=ax)
```

The labels of the plots do not overlap. This can not be used in our examples because we are support matplotlib `1.5.1` (which does not have constrained_layout).

<!--
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:
Add Display object for plot_partial_dependence
Seems that Andy and I agree that we don't need a display object for plot_tree, opinions from other core devs are welcomed.
We'll need to return a Display object in plot_partial_dependence, I don't think this breaks current API (because currently we return nothing).
ping @thomasjpfan do you want to finish this?

Hints from issue discussion:
Yup. I’ll finish this. 😆

Commit message: CLN Fixes bugs

Patch to review:
diff --git a/doc/modules/classes.rst b/doc/modules/classes.rst
index e3b0c7f06c419..229501affb5bf 100644
--- a/doc/modules/classes.rst
+++ b/doc/modules/classes.rst
@@ -664,6 +664,12 @@ Plotting
 
 .. currentmodule:: sklearn
 
+.. autosummary::
+   :toctree: generated/
+   :template: class.rst
+
+   inspection.PartialDependenceDisplay
+
 .. autosummary::
    :toctree: generated/
    :template: function.rst
diff --git a/doc/visualizations.rst b/doc/visualizations.rst
index d21b90d0b4171..bef2532370c21 100644
--- a/doc/visualizations.rst
+++ b/doc/visualizations.rst
@@ -59,6 +59,7 @@ values of the curves.
 .. topic:: Examples:
 
     * :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py`
+    * :ref:`sphx_glr_auto_examples_plot_partial_dependence_visualization_api.py`
 
 Available Plotting Utilities
 ============================
@@ -70,6 +71,7 @@ Functions
 
 .. autosummary::
 
+   inspection.plot_partial_dependence
    metrics.plot_roc_curve
 
 
@@ -80,4 +82,5 @@ Display Objects
 
 .. autosummary::
 
+   inspection.PartialDependenceDisplay
    metrics.RocCurveDisplay
diff --git a/examples/inspection/plot_partial_dependence.py b/examples/inspection/plot_partial_dependence.py
index 5d9a07a081e36..09a194c450c10 100644
--- a/examples/inspection/plot_partial_dependence.py
+++ b/examples/inspection/plot_partial_dependence.py
@@ -94,7 +94,7 @@
 fig = plt.gcf()
 fig.suptitle('Partial dependence of house value on non-location features\n'
              'for the California housing dataset, with MLPRegressor')
-fig.tight_layout(rect=[0, 0.03, 1, 0.95])
+fig.subplots_adjust(wspace=0.8, hspace=0.3)
 
 ##############################################################################
 # Partial Dependence computation for Gradient Boosting
@@ -120,7 +120,9 @@
 fig = plt.gcf()
 fig.suptitle('Partial dependence of house value on non-location features\n'
              'for the California housing dataset, with Gradient Boosting')
-fig.tight_layout(rect=[0, 0.03, 1, 0.95])
+fig.subplots_adjust(wspace=0.8, hspace=0.3)
+plt.show()
+
 
 ##############################################################################
 # Analysis of the plots
diff --git a/examples/plot_partial_dependence_visualization_api.py b/examples/plot_partial_dependence_visualization_api.py
new file mode 100644
index 0000000000000..31f58edcafd3d
--- /dev/null
+++ b/examples/plot_partial_dependence_visualization_api.py
@@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
+tht plot with the Visualization API.
+
+.. note::
+
+    See also :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py`
+
+"""
+print(__doc__)
+
+import matplotlib.pyplot as plt
+from sklearn.datasets import load_boston
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.experimental import enable_hist_gradient_boosting  # noqa
+from sklearn.ensemble import HistGradientBoostingRegressor
+from sklearn.inspection import plot_partial_dependence
+
+
+##############################################################################
+# Train models on the boston housing price dataset
+# ================================================
+#
+# First, we load the boston housing price dataset and split the the dataset
+# into a training and test set. Then, we train a histogram gradient boosting
+# decision tree and a multi-layer perceptron on the training set.
+
+boston = load_boston()
+X, y = boston.data, boston.target
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1,
+                                                    random_state=0)
+
+hgbr = make_pipeline(StandardScaler(), HistGradientBoostingRegressor())
+mlp = make_pipeline(StandardScaler(),
+                    MLPRegressor(hidden_layer_sizes=(100, 100),
+                                 tol=1e-2, max_iter=500, random_state=0))
+hgbr.fit(X_train, y_train)
+mlp.fit(X_train, y_train)
+
+##############################################################################
+# Plotting partial dependence of the two models independently
+# ===========================================================
+#
+# Next, we plot a partial dependence curves for features "LSTAT" and "RM" for
+# for the histogram gradient boosting model.
+fig, ax = plt.subplots(figsize=(12, 6))
+ax.set_title("Histogram Gradient Boosting")
+hgbr_disp = plot_partial_dependence(hgbr, X_test, ["LSTAT", "RM"],
+                                    feature_names=boston.feature_names, ax=ax)
+
+##############################################################################
+# The partial depdendence curves can be plotted for the multi-layer perceptron.
+# In this case `line_kw` was passed to
+# `~sklearn.inspection.plot_partial_dependence` to change the color of the
+# curve and `n_cols` was set to 1 to set the number of columns to 1.
+fig, ax = plt.subplots(figsize=(8, 8))
+ax.set_title("Multi-layer Perceptron")
+mlp_disp = plot_partial_dependence(mlp, X_test, ["LSTAT", "RM"],
+                                   feature_names=boston.feature_names, ax=ax,
+                                   n_cols=1, line_kw={"c": "red"})
+
+##############################################################################
+# Plotting partial dependence of the two models together
+# ======================================================
+#
+# The `hgbr_disp` and `mlp_disp`
+# :class:`~sklearn.inspection.PartialDependenceDisplay` objects contain all the
+# computed information needed to recreate the partial dependence curves. This
+# means we can easily create additional plots without needing to recompute the
+# curves.
+#
+# One way to plot the curves is to place them in the same figure, with the
+# curves of each model on each row. First, we create a figure two axes with two
+# rows and one column. The two axes are passed to the
+# :func:`~sklearn.inspection.PartialDependenceDisplay.plot` functions of
+# `hgbr_disp` and `mlp_disp`. The plot funciton will plot the two curves in the
+# space by the passed in axes. The resulting plot places the histogram gradient
+# boosting partial dependence curves on top of the multi-layer perceptron
+# plots.
+
+fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))
+hgbr_disp.plot(ax=ax1)
+ax1.set_title("Histogram Gradient Boosting")
+mlp_disp.plot(ax=ax2, line_kw={"c": "red"})
+ax2.set_title("Multi-layer Perceptron")
+
+##############################################################################
+# Another way to compare the curves is to plot them on top of each other. Here
+# we create a figure with one row and two columns. The axes are passed into the
+# :func:`~sklearn.inspection.PartialDependenceDisplay.plot` function as list,
+# which will plot the partial dependence curves of each model on the same axes.
+
+# sphinx_gallery_thumbnail_number = 4
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
+hgbr_disp.plot(ax=[ax1, ax2],
+               line_kw={"label": "Histogram Gradient Boosting"})
+mlp_disp.plot(ax=[ax1, ax2],
+              line_kw={"label": "Multi-layer Perceptron", "c": "red"})
+ax1.legend()
+ax2.legend()
+plt.show()
diff --git a/sklearn/inspection/__init__.py b/sklearn/inspection/__init__.py
index 6670e4c576c4d..60e0a499fe1cb 100644
--- a/sklearn/inspection/__init__.py
+++ b/sklearn/inspection/__init__.py
@@ -2,9 +2,11 @@
 from .partial_dependence import partial_dependence
 from .partial_dependence import plot_partial_dependence
 from .permutation_importance import permutation_importance
+from ._plot.partial_dependence import PartialDependenceDisplay
 
 __all__ = [
     'partial_dependence',
     'plot_partial_dependence',
-    'permutation_importance'
+    'permutation_importance',
+    'PartialDependenceDisplay'
 ]
diff --git a/sklearn/inspection/_plot/__init__.py b/sklearn/inspection/_plot/__init__.py
new file mode 100644
index 0000000000000..e69de29bb2d1d
diff --git a/sklearn/inspection/_plot/partial_dependence.py b/sklearn/inspection/_plot/partial_dependence.py
new file mode 100644
index 0000000000000..508115ebdb25b
--- /dev/null
+++ b/sklearn/inspection/_plot/partial_dependence.py
@@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
+        a partial dependence curve of one feature. A tuple of two ints will
+        plot a two-way partial dependence curve as a contour plot.
+
+    feature_names : list of str
+        Feature names corrsponding to the indicies in ``features``.
+
+    target_idx : int
+        - In a multiclass setting, specifies the class for which the PDPs
+          should be computed. Note that for binary classification, the
+          positive class (index 1) is always used.
+        - In a multioutput setting, specifies the task for which the PDPs
+          should be computed
+        Ignored in binary classification or classical regression settings.
+
+    pdp_lim : dict
+        Global min and max average predictions. `pdp_lim[1]` is the global min
+        and max for single partial dependence curves. `pdp_lim[2]` is the
+        global min and max for two-way partial dependence curves.
+
+    deciles : dict
+        Deciles for feature indicies in ``features``.
+
+    Attributes
+    ----------
+    bounding_ax_ : matplotlib Axes or None
+        If `ax` is an axes or None, the `bounding_ax_` is the axes where the
+        grid of partial dependence plots are drawn. If `ax` is  list of axes,
+        `bounding_ax_` is None.
+
+    axes_ : dict of matplotlib Axes
+        If `ax` is an axes or None, `axes_[i, j]` is the axes on the ith row
+        and jth column. If `ax` is a list of axes, `axes_[i]` is the ith item
+        in `ax`.
+
+    lines_ : matplotlib Artists
+        If `ax` is an axes or None, `line_[i, j]` is the partial dependence
+        curve on the ith row and jth column. If `ax` is a list of axes,
+        `lines_[i]` is the partial dependence curve corresponding to the ith
+        item in `ax`. Only plots with a single feature will be populated in
+        `lines_`. Contour plots are stored in `contours_`.
+
+    contours_ : matplotlib Artists
+        If `ax` is an axes or None, `contours_[i, j]` is the partial dependence
+        plot on the ith row and jth column. If `ax` is a list of axes,
+        `contours_[i]` is the partial dependence plot corresponding to the ith
+        item in `ax`. Only plots with a two feature will be populated in
+        `contours_`. Line plots are stored in `lines_`.
+
+    figure_ : matplotlib Figure
+        Figure containing partial dependence plots.
+    """
+    def __init__(self, pd_results, features, feature_names, target_idx,
+                 pdp_lim, deciles):
+        self.pd_results = pd_results
+        self.features = features
+        self.feature_names = feature_names
+        self.target_idx = target_idx
+        self.pdp_lim = pdp_lim
+        self.deciles = deciles
+
+    def plot(self, ax=None, n_cols=3, line_kw=None, contour_kw=None, fig=None):
+        """Plot partial dependence plots.
+
+        Parameters
+        ----------
+        ax : Matplotlib axes, list of Matplotlib axes or None, default=None
+            Axes to plot the partial dependence curves. If a single axes is
+            given, it is treated as a bounding axes and a grid of partial
+            depdendence plots will be drawn on that top of it. If a list of
+            axes are passed, the partial dependence plots will be drawn on
+            those axes. By default, a single bounding axes is created and
+            treated as the single axes case.
+
+        n_cols : int, default=3
+            The maximum number of columns in the grid plot.
+
+        line_kw : dict or None, default=None
+            Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
+            For one-way partial dependence plots.
+
+        contour_kw : dict or None, default=None
+            Dict with keywords passed to the ``matplotlib.pyplot.contourf``
+            call for two-way partial dependence plots.
+
+        fig : Matplotlib figure object or None, default=None
+            A figure object onto which the plots will be drawn, after the
+            figure has been cleared. By default, a new one is created.
+
+            .. deprecated:: 0.22
+
+        Returns
+        -------
+        display: :class:`~sklearn.inspection.PartialDependenceDisplay`
+        """
+
+        check_matplotlib_support("plot_partial_dependence")
+        import matplotlib.pyplot as plt  # noqa
+        from matplotlib import transforms  # noqa
+        from matplotlib.ticker import MaxNLocator  # noqa
+        from matplotlib.ticker import ScalarFormatter  # noqa
+        from matplotlib.gridspec import GridSpecFromSubplotSpec  # noqa
+
+        if line_kw is None:
+            line_kw = {}
+        if contour_kw is None:
+            contour_kw = {}
+
+        contour_kw = {**{"alpha": 0.75}, **contour_kw}
+
+        if fig is not None:
+            warnings.warn("The fig parameter is deprecated in version "
+                          "0.22 and will be removed in version 0.24",
+                          DeprecationWarning)
+            fig.clear()
+
+        if ax is None:
+            if fig is None:
+                _, ax = plt.subplots()
+            else:
+                ax = fig.add_subplot(111)
+
+        if isinstance(ax, list):
+            if len(ax) != len(self.features):
+                raise ValueError("Expected len(ax) == len(features), "
+                                 "got len(ax) = {}".format(len(ax)))
+            self.bounding_ax_ = None
+            self.axes_ = dict(enumerate(ax))
+            self.figure_ = ax[0].figure
+            self.lines_ = {}
+            self.contours_ = {}
+
+            def _get_axes(i):
+                return self.axes_[i]
+
+            def _store_line(line, i):
+                self.lines_[i] = line
+
+            def _store_contour(contour, i):
+                self.contours_[i] = contour
+
+        else:  # single axes
+            self.bounding_ax_ = ax
+            self.axes_ = {}
+            self.figure_ = ax.figure
+            self.lines_ = {}
+            self.contours_ = {}
+
+            ax.set_axis_off()
+            n_cols = min(n_cols, len(self.features))
+            n_rows = int(np.ceil(len(self.features) / float(n_cols)))
+
+            gs = GridSpecFromSubplotSpec(n_rows, n_cols,
+                                         subplot_spec=ax.get_subplotspec())
+
+            def _get_axes(i):
+                row_idx, col_idx = i // n_cols, i % n_cols
+                axi = self.figure_.add_subplot(gs[row_idx, col_idx])
+                self.axes_[row_idx, col_idx] = axi
+                return axi
+
+            def _store_line(line, i):
+                row_idx, col_idx = i // n_cols, i % n_cols
+                self.lines_[row_idx, col_idx] = line
+
+            def _store_contour(contour, i):
+                row_idx, col_idx = i // n_cols, i % n_cols
+                self.contours_[row_idx, col_idx] = contour
+
+        # create contour levels for two-way plots
+        if 2 in self.pdp_lim:
+            Z_level = np.linspace(*self.pdp_lim[2], num=8)
+
+        for i, fx, (avg_preds, values) in zip(count(), self.features,
+                                              self.pd_results):
+            axi = _get_axes(i)
+
+            if len(values) == 1:
+                line = axi.plot(values[0], avg_preds[self.target_idx].ravel(),
+                                **line_kw)[0]
+                _store_line(line, i)
+            else:
+                # contour plot
+                XX, YY = np.meshgrid(values[0], values[1])
+                Z = avg_preds[self.target_idx].T
+                CS = axi.contour(XX, YY, Z, levels=Z_level, linewidths=0.5,
+                                 colors='k')
+                contour = axi.contourf(XX, YY, Z, levels=Z_level,
+                                       vmax=Z_level[-1], vmin=Z_level[0],
+                                       **contour_kw)
+                _store_contour(contour, i)
+                axi.clabel(CS, fmt='%2.2f', colors='k', fontsize=10,
+                           inline=True)
+
+            trans = transforms.blended_transform_factory(axi.transData,
+                                                         axi.transAxes)
+            ylim = axi.get_ylim()
+            axi.vlines(self.deciles[fx[0]], 0, 0.05, transform=trans,
+                       color='k')
+            axi.set_xlabel(self.feature_names[fx[0]])
+            axi.set_ylim(ylim)
+
+            if len(values) == 1:
+                axi.set_ylabel('Partial dependence')
+                axi.set_ylim(self.pdp_lim[1])
+            else:
+                # contour plot
+                trans = transforms.blended_transform_factory(axi.transAxes,
+                                                             axi.transData)
+                xlim = axi.get_xlim()
+                axi.hlines(self.deciles[fx[1]], 0, 0.05, transform=trans,
+                           color='k')
+                # hline erases xlim
+                axi.set_ylabel(self.feature_names[fx[1]])
+                axi.set_xlim(xlim)
+
+        return self
diff --git a/sklearn/inspection/partial_dependence.py b/sklearn/inspection/partial_dependence.py
index 692736cb56d70..a559ee1f62295 100644
--- a/sklearn/inspection/partial_dependence.py
+++ b/sklearn/inspection/partial_dependence.py
@@ -5,7 +5,7 @@
 #          Nicolas Hug
 # License: BSD 3 clause
 
-from itertools import count
+from itertools import chain
 import numbers
 from collections.abc import Iterable
 
@@ -23,6 +23,7 @@
 from ..ensemble.gradient_boosting import BaseGradientBoosting
 from sklearn.ensemble._hist_gradient_boosting.gradient_boosting import (
     BaseHistGradientBoosting)
+from ._plot.partial_dependence import PartialDependenceDisplay
 
 
 __all__ = ['partial_dependence', 'plot_partial_dependence']
@@ -374,7 +375,7 @@ def plot_partial_dependence(estimator, X, features, feature_names=None,
                             target=None, response_method='auto', n_cols=3,
                             grid_resolution=100, percentiles=(0.05, 0.95),
                             method='auto', n_jobs=None, verbose=0, fig=None,
-                            line_kw=None, contour_kw=None):
+                            line_kw=None, contour_kw=None, ax=None):
     """Partial dependence plots.
 
     The ``len(features)`` plots are arranged in a grid with ``n_cols``
@@ -388,19 +389,23 @@ def plot_partial_dependence(estimator, X, features, feature_names=None,
         A fitted estimator object implementing `predict`, `predict_proba`,
         or `decision_function`. Multioutput-multiclass classifiers are not
         supported.
+
     X : array-like, shape (n_samples, n_features)
         The data to use to build the grid of values on which the dependence
         will be evaluated. This is usually the training data.
+
     features : list of {int, str, pair of int, pair of str}
         The target features for which to create the PDPs.
         If features[i] is an int or a string, a one-way PDP is created; if
         features[i] is a tuple, a two-way PDP is created. Each tuple must be
         of size 2.
         if any entry is a string, then it must be in ``feature_names``.
+
     feature_names : seq of str, shape (n_features,), optional
         Name of each feature; feature_names[i] holds the name of the feature
         with index i. By default, the name of the feature corresponds to
         their numerical index.
+
     target : int, optional (default=None)
         - In a multiclass setting, specifies the class for which the PDPs
           should be computed. Note that for binary classification, the
@@ -408,6 +413,7 @@ def plot_partial_dependence(estimator, X, features, feature_names=None,
         - In a multioutput setting, specifies the task for which the PDPs
           should be computed
         Ignored in binary classification or classical regression settings.
+
     response_method : 'auto', 'predict_proba' or 'decision_function', \
             optional (default='auto') :
         Specifies whether to use :term:`predict_proba` or
@@ -417,14 +423,18 @@ def plot_partial_dependence(estimator, X, features, feature_names=None,
         and we revert to :term:`decision_function` if it doesn't exist. If
         ``method`` is 'recursion', the response is always the output of
         :term:`decision_function`.
+
     n_cols : int, optional (default=3)
         The maximum number of columns in the grid plot.
+
     grid_resolution : int, optional (default=100)
         The number of equally spaced points on the axes of the plots, for each
         target feature.
+
     percentiles : tuple of float, optional (default=(0.05, 0.95))
         The lower and upper percentile used to create the extreme values
         for the PDP axes. Must be in [0, 1].
+
     method : str, optional (default='auto')
         The method to use to calculate the partial dependence predictions:
 
@@ -448,23 +458,37 @@ def plot_partial_dependence(estimator, X, features, feature_names=None,
         - 'auto':
           - 'recursion' is used for estimators that supports it.
           - 'brute' is used for all other estimators.
+
     n_jobs : int, optional (default=None)
         The number of CPUs to use to compute the partial dependences.
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
         for more details.
+
     verbose : int, optional (default=0)
         Verbose output during PD computations.
+
     fig : Matplotlib figure object, optional (default=None)
         A figure object onto which the plots will be drawn, after the figure
         has been cleared. By default, a new one is created.
+
+        .. deprecated:: 0.22
+
     line_kw : dict, optional
         Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
         For one-way partial dependence plots.
+
     contour_kw : dict, optional
-        Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
+        Dict with keywords passed to the ``matplotlib.pyplot.contourf`` call.
         For two-way partial dependence plots.
 
+    ax : Matplotlib axes, list of Matplotlib axes or None, (default=None)
+        By default,
+
+    Returns
+    -------
+    sklearn.inspection.PartialDependenceDisplay
+
     Examples
     --------
     >>> from sklearn.datasets import make_friedman1
@@ -538,7 +562,7 @@ def convert_feature(fx):
         if isinstance(fxs, (numbers.Integral, str)):
             fxs = (fxs,)
         try:
-            fxs = [convert_feature(fx) for fx in fxs]
+            fxs = tuple(convert_feature(fx) for fx in fxs)
         except TypeError:
             raise ValueError('Each entry in features must be either an int, '
                              'a string, or an iterable of size at most 2.')
@@ -550,21 +574,19 @@ def convert_feature(fx):
 
     features = tmp_features
 
-    names = []
-    try:
-        for fxs in features:
-            names_ = []
-            # explicit loop so "i" is bound for exception below
-            for i in fxs:
-                names_.append(feature_names[i])
-            names.append(names_)
-    except IndexError:
-        raise ValueError('All entries of features must be less than '
-                         'len(feature_names) = {0}, got {1}.'
-                         .format(len(feature_names), i))
+    if isinstance(ax, list):
+        if len(ax) != len(features):
+            raise ValueError("Expected len(ax) == len(features), "
+                             "got len(ax) = {}".format(len(ax)))
+
+    for i in chain.from_iterable(features):
+        if i >= len(feature_names):
+            raise ValueError('All entries of features must be less than '
+                             'len(feature_names) = {0}, got {1}.'
+                             .format(len(feature_names), i))
 
     # compute averaged predictions
-    pd_result = Parallel(n_jobs=n_jobs, verbose=verbose)(
+    pd_results = Parallel(n_jobs=n_jobs, verbose=verbose)(
         delayed(partial_dependence)(estimator, X, fxs,
                                     response_method=response_method,
                                     method=method,
@@ -577,7 +599,7 @@ def convert_feature(fx):
     # Also note: as multiclass-multioutput classifiers are not supported,
     # multiclass and multioutput scenario are mutually exclusive. So there is
     # no risk of overwriting target_idx here.
-    avg_preds, _ = pd_result[0]  # checking the first result is enough
+    avg_preds, _ = pd_results[0]  # checking the first result is enough
     if is_regressor(estimator) and avg_preds.shape[0] > 1:
         if target is None:
             raise ValueError(
@@ -587,9 +609,9 @@ def convert_feature(fx):
                 'target must be in [0, n_tasks], got {}.'.format(target))
         target_idx = target
 
-    # get global min and max values of PD grouped by plot type
+    # get global min and max average predictions of PD grouped by plot type
     pdp_lim = {}
-    for avg_preds, values in pd_result:
+    for avg_preds, values in pd_results:
         min_pd = avg_preds[target_idx].min()
         max_pd = avg_preds[target_idx].max()
         n_fx = len(values)
@@ -598,71 +620,12 @@ def convert_feature(fx):
         max_pd = max(max_pd, old_max_pd)
         pdp_lim[n_fx] = (min_pd, max_pd)
 
-    # create contour levels for two-way plots
-    if 2 in pdp_lim:
-        Z_level = np.linspace(*pdp_lim[2], num=8)
-
-    if fig is None:
-        fig = plt.figure()
-    else:
-        fig.clear()
-
-    if line_kw is None:
-        line_kw = {'color': 'green'}
-    if contour_kw is None:
-        contour_kw = {}
-
-    n_cols = min(n_cols, len(features))
-    n_rows = int(np.ceil(len(features) / float(n_cols)))
-    axs = []
-    for i, fx, name, (avg_preds, values) in zip(
-            count(), features, names, pd_result):
-        ax = fig.add_subplot(n_rows, n_cols, i + 1)
-
-        if len(values) == 1:
-            ax.plot(values[0], avg_preds[target_idx].ravel(), **line_kw)
-        else:
-            # make contour plot
-            assert len(values) == 2
-            XX, YY = np.meshgrid(values[0], values[1])
-            Z = avg_preds[target_idx].T
-            CS = ax.contour(XX, YY, Z, levels=Z_level, linewidths=0.5,
-                            colors='k')
-            ax.contourf(XX, YY, Z, levels=Z_level, vmax=Z_level[-1],
-                        vmin=Z_level[0], alpha=0.75, **contour_kw)
-            ax.clabel(CS, fmt='%2.2f', colors='k', fontsize=10, inline=True)
-
-        # plot data deciles + axes labels
-        deciles = mquantiles(X[:, fx[0]], prob=np.arange(0.1, 1.0, 0.1))
-        trans = transforms.blended_transform_factory(ax.transData,
-                                                     ax.transAxes)
-        ylim = ax.get_ylim()
-        ax.vlines(deciles, [0], 0.05, transform=trans, color='k')
-        ax.set_xlabel(name[0])
-        ax.set_ylim(ylim)
-
-        # prevent x-axis ticks from overlapping
-        ax.xaxis.set_major_locator(MaxNLocator(nbins=6, prune='lower'))
-        tick_formatter = ScalarFormatter()
-        tick_formatter.set_powerlimits((-3, 4))
-        ax.xaxis.set_major_formatter(tick_formatter)
-
-        if len(values) > 1:
-            # two-way PDP - y-axis deciles + labels
-            deciles = mquantiles(X[:, fx[1]], prob=np.arange(0.1, 1.0, 0.1))
-            trans = transforms.blended_transform_factory(ax.transAxes,
-                                                         ax.transData)
-            xlim = ax.get_xlim()
-            ax.hlines(deciles, [0], 0.05, transform=trans, color='k')
-            ax.set_ylabel(name[1])
-            # hline erases xlim
-            ax.set_xlim(xlim)
-        else:
-            ax.set_ylabel('Partial dependence')
-
-        if len(values) == 1:
-            ax.set_ylim(pdp_lim[1])
-        axs.append(ax)
+    deciles = {}
+    for fx in chain.from_iterable(features):
+        if fx not in deciles:
+            deciles[fx] = mquantiles(X[:, fx], prob=np.arange(0.1, 1.0, 0.1))
 
-    fig.subplots_adjust(bottom=0.15, top=0.7, left=0.1, right=0.95, wspace=0.4,
-                        hspace=0.3)
+    display = PartialDependenceDisplay(pd_results, features, feature_names,
+                                       target_idx, pdp_lim, deciles)
+    return display.plot(ax=ax, n_cols=n_cols, line_kw=line_kw,
+                        contour_kw=contour_kw, fig=fig)
diff --git a/sklearn/inspection/tests/test_partial_dependence.py b/sklearn/inspection/tests/test_partial_dependence.py
index 19399224e07ba..83eb540536a37 100644
--- a/sklearn/inspection/tests/test_partial_dependence.py
+++ b/sklearn/inspection/tests/test_partial_dependence.py
@@ -7,7 +7,6 @@
 
 import sklearn
 from sklearn.inspection import partial_dependence
-from sklearn.inspection import plot_partial_dependence
 from sklearn.inspection.partial_dependence import (
     _grid_from_X,
     _partial_dependence_brute,
@@ -22,7 +21,7 @@
 from sklearn.linear_model import LogisticRegression
 from sklearn.linear_model import MultiTaskLasso
 from sklearn.tree import DecisionTreeRegressor
-from sklearn.datasets import load_boston, load_iris
+from sklearn.datasets import load_iris
 from sklearn.datasets import make_classification, make_regression
 from sklearn.cluster import KMeans
 from sklearn.metrics import r2_score
@@ -438,152 +437,3 @@ def test_partial_dependence_pipeline():
         values_pipe[0],
         values_clf[0] * scaler.scale_[features] + scaler.mean_[features]
     )
-
-
-def test_plot_partial_dependence(pyplot):
-    # Test partial dependence plot function.
-    boston = load_boston()
-    clf = GradientBoostingRegressor(n_estimators=10, random_state=1)
-    clf.fit(boston.data, boston.target)
-
-    grid_resolution = 25
-    plot_partial_dependence(clf, boston.data, [0, 1, (0, 1)],
-                            grid_resolution=grid_resolution,
-                            feature_names=boston.feature_names)
-    fig = pyplot.gcf()
-    axs = fig.get_axes()
-    assert len(axs) == 3
-    assert all(ax.has_data for ax in axs)
-
-    # check with str features and array feature names
-    plot_partial_dependence(clf, boston.data, ['CRIM', 'ZN',
-                                               ('CRIM', 'ZN')],
-                            grid_resolution=grid_resolution,
-                            feature_names=boston.feature_names)
-
-    fig = pyplot.gcf()
-    axs = fig.get_axes()
-    assert len(axs) == 3
-    assert all(ax.has_data for ax in axs)
-
-    # check with list feature_names
-    feature_names = boston.feature_names.tolist()
-    plot_partial_dependence(clf, boston.data, ['CRIM', 'ZN',
-                                               ('CRIM', 'ZN')],
-                            grid_resolution=grid_resolution,
-                            feature_names=feature_names)
-    fig = pyplot.gcf()
-    axs = fig.get_axes()
-    assert len(axs) == 3
-    assert all(ax.has_data for ax in axs)
-
-
-def test_plot_partial_dependence_multiclass(pyplot):
-    grid_resolution = 25
-    clf = GradientBoostingClassifier(n_estimators=10, random_state=1)
-    iris = load_iris()
-
-    # Test partial dependence plot function on multi-class input.
-    clf.fit(iris.data, iris.target)
-    plot_partial_dependence(clf, iris.data, [0, 1],
-                            target=0,
-                            grid_resolution=grid_resolution)
-    fig = pyplot.gcf()
-    axs = fig.get_axes()
-    assert len(axs) == 2
-    assert all(ax.has_data for ax in axs)
-
-    # now with symbol labels
-    target = iris.target_names[iris.target]
-    clf.fit(iris.data, target)
-    plot_partial_dependence(clf, iris.data, [0, 1],
-                            target='setosa',
-                            grid_resolution=grid_resolution)
-    fig2 = pyplot.gcf()
-    axs2 = fig2.get_axes()
-    assert len(axs2) == 2
-    assert all(ax.has_data for ax in axs2)
-
-    # check that the pd plots are the same for 0 and "setosa"
-    assert all(axs[0].lines[0]._y == axs2[0].lines[0]._y)
-    # check that the pd plots are different for another target
-    clf.fit(iris.data, iris.target)
-    plot_partial_dependence(clf, iris.data, [0, 1],
-                            target=1,
-                            grid_resolution=grid_resolution)
-    fig3 = pyplot.gcf()
-    axs3 = fig3.get_axes()
-    assert any(axs[0].lines[0]._y != axs3[0].lines[0]._y)
-
-
-def test_plot_partial_dependence_multioutput(pyplot):
-    # Test partial dependence plot function on multi-output input.
-    (X, y), _ = multioutput_regression_data
-    clf = LinearRegression()
-    clf.fit(X, y)
-
-    grid_resolution = 25
-    plot_partial_dependence(clf, X, [0, 1],
-                            target=0,
-                            grid_resolution=grid_resolution)
-    fig = pyplot.gcf()
-    axs = fig.get_axes()
-    assert len(axs) == 2
-    assert all(ax.has_data for ax in axs)
-
-    plot_partial_dependence(clf, X, [0, 1],
-                            target=1,
-                            grid_resolution=grid_resolution)
-    fig = pyplot.gcf()
-    axs = fig.get_axes()
-    assert len(axs) == 2
-    assert all(ax.has_data for ax in axs)
-
-
-@pytest.mark.parametrize(
-    "data, params, err_msg",
-    [(multioutput_regression_data[0], {"target": None, 'features': [0]},
-      "target must be specified for multi-output"),
-     (multioutput_regression_data[0], {"target": -1, 'features': [0]},
-      r'target must be in \[0, n_tasks\]'),
-     (multioutput_regression_data[0], {"target": 100, 'features': [0]},
-      r'target must be in \[0, n_tasks\]'),
-     (make_classification(random_state=0),
-     {'features': ['foobar'], 'feature_names': None},
-     'Feature foobar not in feature_names'),
-     (make_classification(random_state=0),
-     {'features': ['foobar'], 'feature_names': ['abcd', 'def']},
-      'Feature foobar not in feature_names'),
-     (make_classification(random_state=0), {'features': [(1, 2, 3)]},
-      'Each entry in features must be either an int, '),
-     (make_classification(random_state=0), {'features': [1, {}]},
-      'Each entry in features must be either an int, '),
-     (make_classification(random_state=0), {'features': [tuple()]},
-      'Each entry in features must be either an int, '),
-     (make_classification(random_state=0),
-      {'features': [123], 'feature_names': ['blahblah']},
-      'All entries of features must be less than '),
-     (make_classification(random_state=0),
-      {'features': [0, 1, 2], 'feature_names': ['a', 'b', 'a']},
-      'feature_names should not contain duplicates')]
-)
-def test_plot_partial_dependence_error(pyplot, data, params, err_msg):
-    X, y = data
-    estimator = LinearRegression().fit(X, y)
-
-    with pytest.raises(ValueError, match=err_msg):
-        plot_partial_dependence(estimator, X, **params)
-
-
-def test_plot_partial_dependence_fig(pyplot):
-    # Make sure fig object is correctly used if not None
-    (X, y), _ = regression_data
-    clf = LinearRegression()
-    clf.fit(X, y)
-
-    fig = pyplot.figure()
-    grid_resolution = 25
-    plot_partial_dependence(
-        clf, X, [0, 1], target=0, grid_resolution=grid_resolution, fig=fig)
-
-    assert pyplot.gcf() is fig
diff --git a/sklearn/inspection/tests/test_plot_partial_dependence.py b/sklearn/inspection/tests/test_plot_partial_dependence.py
new file mode 100644
index 0000000000000..0fabc201963eb
--- /dev/null
+++ b/sklearn/inspection/tests/test_plot_partial_dependence.py
@@ -0,0 +1,325 @@
+import numpy as np
+from scipy.stats.mstats import mquantiles
+
+import pytest
+from numpy.testing import assert_allclose
+
+from sklearn.datasets import load_boston
+from sklearn.datasets import load_iris
+from sklearn.datasets import make_classification, make_regression
+from sklearn.ensemble import GradientBoostingRegressor
+from sklearn.ensemble import GradientBoostingClassifier
+from sklearn.linear_model import LinearRegression
+from sklearn.inspection import plot_partial_dependence
+
+
+boston = load_boston()
+
+
+@pytest.fixture(scope="module")
+def clf_boston():
+    clf = GradientBoostingRegressor(n_estimators=10, random_state=1)
+    clf.fit(boston.data, boston.target)
+    return clf
+
+
+@pytest.mark.parametrize("grid_resolution", [10, 20])
+def test_plot_partial_dependence(pyplot, clf_boston, grid_resolution):
+    # Test partial dependence plot function.
+    feature_names = boston.feature_names
+    disp = plot_partial_dependence(clf_boston, boston.data,
+                                   [0, 1, (0, 1)],
+                                   grid_resolution=grid_resolution,
+                                   feature_names=feature_names,
+                                   contour_kw={"cmap": "jet"})
+    fig = pyplot.gcf()
+    axs = fig.get_axes()
+    assert disp.figure_ is fig
+    assert len(axs) == 4
+
+    assert disp.bounding_ax_ is not None
+    assert len(disp.axes_) == 3
+    assert len(disp.lines_) == 2
+    assert len(disp.contours_) == 1
+    assert disp.features == [(0, ), (1, ), (0, 1)]
+    assert np.all(disp.feature_names == feature_names)
+    assert len(disp.deciles) == 2
+    for i in [0, 1]:
+        assert_allclose(disp.deciles[i],
+                        mquantiles(boston.data[:, i],
+                                   prob=np.arange(0.1, 1.0, 0.1)))
+
+    single_feature_positions = [(0, 0), (0, 1)]
+
+    for i, pos in enumerate(single_feature_positions):
+        ax = disp.axes_[pos]
+        assert ax.get_xlabel() == boston.feature_names[i]
+        assert ax.get_ylabel() == "Partial dependence"
+        assert_allclose(ax.get_ylim(), disp.pdp_lim[1])
+
+        line = disp.lines_[pos]
+
+        avg_preds, values = disp.pd_results[i]
+        assert avg_preds.shape == (1, grid_resolution)
+        target_idx = disp.target_idx
+
+        line_data = line.get_data()
+        assert_allclose(line_data[0], values[0])
+        assert_allclose(line_data[1], avg_preds[target_idx].ravel())
+
+    # two feature position
+    ax = disp.axes_[0, 2]
+    coutour = disp.contours_[0, 2]
+    expected_levels = np.linspace(*disp.pdp_lim[2], num=8)
+    assert_allclose(coutour.levels, expected_levels)
+    assert coutour.get_cmap().name == "jet"
+    assert ax.get_xlabel() == boston.feature_names[0]
+    assert ax.get_ylabel() == boston.feature_names[1]
+
+
+def test_plot_partial_dependence_str_features(pyplot, clf_boston):
+    grid_resolution = 25
+    # check with str features and array feature names and single column
+    disp = plot_partial_dependence(clf_boston, boston.data,
+                                   [('CRIM', 'ZN'), 'ZN'],
+                                   grid_resolution=grid_resolution,
+                                   feature_names=boston.feature_names,
+                                   n_cols=1, line_kw={"alpha": 0.8})
+    fig = pyplot.gcf()
+    axs = fig.get_axes()
+    assert len(axs) == 3
+
+    assert disp.figure_ is fig
+    assert len(disp.axes_) == 2
+    assert len(disp.lines_) == 1
+    assert len(disp.contours_) == 1
+
+    # line
+    ax = disp.axes_[1, 0]
+    assert ax.get_xlabel() == "ZN"
+    assert ax.get_ylabel() == "Partial dependence"
+
+    line = disp.lines_[1, 0]
+    avg_preds, values = disp.pd_results[1]
+    target_idx = disp.target_idx
+    assert line.get_alpha() == 0.8
+
+    line_data = line.get_data()
+    assert_allclose(line_data[0], values[0])
+    assert_allclose(line_data[1], avg_preds[target_idx].ravel())
+
+    # contour
+    ax = disp.axes_[0, 0]
+    coutour = disp.contours_[0, 0]
+    expect_levels = np.linspace(*disp.pdp_lim[2], num=8)
+    assert_allclose(coutour.levels, expect_levels)
+    assert ax.get_xlabel() == "CRIM"
+    assert ax.get_ylabel() == "ZN"
+
+
+def test_plot_partial_dependence_custom_axes(pyplot, clf_boston):
+    grid_resolution = 25
+    fig, (ax1, ax2) = pyplot.subplots(1, 2)
+    feature_names = boston.feature_names.tolist()
+    disp = plot_partial_dependence(clf_boston, boston.data,
+                                   ['CRIM', ('CRIM', 'ZN')],
+                                   grid_resolution=grid_resolution,
+                                   feature_names=feature_names, ax=[ax1, ax2])
+    assert fig is disp.figure_
+    assert disp.bounding_ax_ is None
+    assert len(disp.axes_) == 2
+    assert disp.axes_[0] is ax1
+    assert disp.axes_[1] is ax2
+
+    ax = disp.axes_[0]
+    assert ax.get_xlabel() == "CRIM"
+    assert ax.get_ylabel() == "Partial dependence"
+
+    line = disp.lines_[0]
+    avg_preds, values = disp.pd_results[0]
+    target_idx = disp.target_idx
+
+    line_data = line.get_data()
+    assert_allclose(line_data[0], values[0])
+    assert_allclose(line_data[1], avg_preds[target_idx].ravel())
+
+    # contour
+    ax = disp.axes_[1]
+    coutour = disp.contours_[1]
+    expect_levels = np.linspace(*disp.pdp_lim[2], num=8)
+    assert_allclose(coutour.levels, expect_levels)
+    assert ax.get_xlabel() == "CRIM"
+    assert ax.get_ylabel() == "ZN"
+
+
+def test_plot_partial_dependence_incorrent_num_axes(pyplot, clf_boston):
+    grid_resolution = 25
+    fig, (ax1, ax2, ax3) = pyplot.subplots(1, 3)
+
+    msg = r"Expected len\(ax\) == len\(features\), got len\(ax\) = 3"
+    with pytest.raises(ValueError, match=msg):
+        plot_partial_dependence(clf_boston, boston.data,
+                                ['CRIM', ('CRIM', 'ZN')],
+                                grid_resolution=grid_resolution,
+                                feature_names=boston.feature_names,
+                                ax=[ax1, ax2, ax3])
+
+    disp = plot_partial_dependence(clf_boston, boston.data,
+                                   ['CRIM', ('CRIM', 'ZN')],
+                                   grid_resolution=grid_resolution,
+                                   feature_names=boston.feature_names)
+
+    with pytest.raises(ValueError, match=msg):
+        disp.plot(ax=[ax1, ax2, ax3])
+
+
+def test_plot_partial_dependence_multiclass(pyplot):
+    grid_resolution = 25
+    clf_int = GradientBoostingClassifier(n_estimators=10, random_state=1)
+    iris = load_iris()
+
+    # Test partial dependence plot function on multi-class input.
+    clf_int.fit(iris.data, iris.target)
+    disp_target_0 = plot_partial_dependence(clf_int, iris.data, [0, 1],
+                                            target=0,
+                                            grid_resolution=grid_resolution)
+    assert disp_target_0.figure_ is pyplot.gcf()
+    assert len(disp_target_0.axes_) == 2
+    assert len(disp_target_0.lines_) == 2
+    assert len(disp_target_0.contours_) == 0
+    assert disp_target_0.target_idx == 0
+
+    # now with symbol labels
+    target = iris.target_names[iris.target]
+    clf_symbol = GradientBoostingClassifier(n_estimators=10, random_state=1)
+    clf_symbol.fit(iris.data, target)
+    disp_symbol = plot_partial_dependence(clf_symbol, iris.data, [0, 1],
+                                          target='setosa',
+                                          grid_resolution=grid_resolution)
+    assert disp_symbol.figure_ is pyplot.gcf()
+    assert len(disp_symbol.axes_) == 2
+    assert len(disp_symbol.lines_) == 2
+    assert len(disp_symbol.contours_) == 0
+    assert disp_symbol.target_idx == 0
+
+    for int_result, symbol_result in zip(disp_target_0.pd_results,
+                                         disp_symbol.pd_results):
+        avg_preds_int, values_int = int_result
+        avg_preds_symbol, values_symbol = symbol_result
+        assert_allclose(avg_preds_int, avg_preds_symbol)
+        assert_allclose(values_int, values_symbol)
+
+    # check that the pd plots are different for another target
+    disp_target_1 = plot_partial_dependence(clf_int, iris.data, [0, 1],
+                                            target=1,
+                                            grid_resolution=grid_resolution)
+    target_0_data_y = disp_target_0.lines_[0, 0].get_data()[1]
+    target_1_data_y = disp_target_1.lines_[0, 0].get_data()[1]
+    assert any(target_0_data_y != target_1_data_y)
+
+
+multioutput_regression_data = make_regression(n_samples=50, n_targets=2,
+                                              random_state=0)
+
+
+@pytest.fixture(scope="module")
+def clf_multioutput():
+    clf = LinearRegression()
+    return clf.fit(*multioutput_regression_data)
+
+
+@pytest.mark.parametrize("target", [0, 1])
+def test_plot_partial_dependence_multioutput(pyplot, target, clf_multioutput):
+    # Test partial dependence plot function on multi-output input.
+    X, y = multioutput_regression_data
+    clf = LinearRegression()
+    clf.fit(X, y)
+
+    grid_resolution = 25
+    disp = plot_partial_dependence(clf_multioutput, X, [0, 1], target=target,
+                                   grid_resolution=grid_resolution)
+    fig = pyplot.gcf()
+    axs = fig.get_axes()
+    assert len(axs) == 3
+    assert disp.target_idx == target
+    assert disp.bounding_ax_ is not None
+
+    positions = [(0, 0), (0, 1)]
+
+    for i, pos in enumerate(positions):
+        ax = disp.axes_[pos]
+        assert ax.get_ylabel() == "Partial dependence"
+        assert ax.get_xlabel() == "{}".format(i)
+
+
+dummy_classification_data = make_classification(random_state=0)
+
+
+@pytest.mark.parametrize(
+    "data, params, err_msg",
+    [(multioutput_regression_data, {"target": None, 'features': [0]},
+      "target must be specified for multi-output"),
+     (multioutput_regression_data, {"target": -1, 'features': [0]},
+      r'target must be in \[0, n_tasks\]'),
+     (multioutput_regression_data, {"target": 100, 'features': [0]},
+      r'target must be in \[0, n_tasks\]'),
+     (dummy_classification_data,
+     {'features': ['foobar'], 'feature_names': None},
+     'Feature foobar not in feature_names'),
+     (dummy_classification_data,
+     {'features': ['foobar'], 'feature_names': ['abcd', 'def']},
+      'Feature foobar not in feature_names'),
+     (dummy_classification_data, {'features': [(1, 2, 3)]},
+      'Each entry in features must be either an int, '),
+     (dummy_classification_data, {'features': [1, {}]},
+      'Each entry in features must be either an int, '),
+     (dummy_classification_data, {'features': [tuple()]},
+      'Each entry in features must be either an int, '),
+     (dummy_classification_data,
+      {'features': [123], 'feature_names': ['blahblah']},
+      'All entries of features must be less than '),
+     (dummy_classification_data,
+      {'features': [0, 1, 2], 'feature_names': ['a', 'b', 'a']},
+      'feature_names should not contain duplicates')]
+)
+def test_plot_partial_dependence_error(pyplot, data, params, err_msg):
+    X, y = data
+    estimator = LinearRegression().fit(X, y)
+
+    with pytest.raises(ValueError, match=err_msg):
+        plot_partial_dependence(estimator, X, **params)
+
+
+@pytest.mark.parametrize("params, err_msg", [
+    ({'target': 4, 'features': [0]},
+     'target not in est.classes_, got 4'),
+    ({'target': None, 'features': [0]},
+     'target must be specified for multi-class'),
+    ({'target': 1, 'features': [4.5]},
+     'Each entry in features must be either an int,'),
+])
+def test_plot_partial_dependence_multiclass_error(pyplot, params, err_msg):
+    iris = load_iris()
+    clf = GradientBoostingClassifier(n_estimators=10, random_state=1)
+    clf.fit(iris.data, iris.target)
+
+    with pytest.raises(ValueError, match=err_msg):
+        plot_partial_dependence(clf, iris.data, **params)
+
+
+def test_plot_partial_dependence_fig_deprecated(pyplot):
+    # Make sure fig object is correctly used if not None
+    X, y = make_regression(n_samples=50, random_state=0)
+    clf = LinearRegression()
+    clf.fit(X, y)
+
+    fig = pyplot.figure()
+    grid_resolution = 25
+
+    msg = ("The fig parameter is deprecated in version 0.22 and will be "
+           "removed in version 0.24")
+    with pytest.warns(DeprecationWarning, match=msg):
+        plot_partial_dependence(
+            clf, X, [0, 1], target=0, grid_resolution=grid_resolution, fig=fig)
+
+    assert pyplot.gcf() is fig
diff --git a/sklearn/setup.py b/sklearn/setup.py
index 53f6d3f6eb30c..7f956ed5f5758 100644
--- a/sklearn/setup.py
+++ b/sklearn/setup.py
@@ -32,6 +32,7 @@ def configuration(parent_package='', top_path=None):
     config.add_subpackage('impute')
     config.add_subpackage('impute/tests')
     config.add_subpackage('inspection')
+    config.add_subpackage('inspection/_plot')
     config.add_subpackage('inspection/tests')
     config.add_subpackage('mixture')
     config.add_subpackage('mixture/tests')


Code review comments:
diff_hunk: @@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
path: examples/plot_partial_dependence_visualization_api.py
text: @user1:
```suggestion
example, we show how to plot partial dependence plots and how to quickly customize
```

diff_hunk: @@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
+tht plot with the Visualization API.
+
+.. note::
+
+    See also :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py`
+
+"""
+print(__doc__)
+
+import matplotlib.pyplot as plt
+from sklearn.datasets import load_boston
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.experimental import enable_hist_gradient_boosting  # noqa
+from sklearn.ensemble import HistGradientBoostingRegressor
+from sklearn.inspection import plot_partial_dependence
+
+
+##############################################################################
+# Train models on the boston housing price dataset
+# ================================================
+#
+# First, we load the boston housing price dataset and split the the dataset
+# into a training and test set. Then, we train a histogram gradient boosting
+# decision tree and a multi-layer perceptron on the training set.
+
+boston = load_boston()
+X, y = boston.data, boston.target
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1,
+                                                    random_state=0)
+
+hgbr = make_pipeline(StandardScaler(), HistGradientBoostingRegressor())
+mlp = make_pipeline(StandardScaler(),
+                    MLPRegressor(hidden_layer_sizes=(100, 100),
+                                 tol=1e-2, max_iter=500, random_state=0))
+hgbr.fit(X_train, y_train)
+mlp.fit(X_train, y_train)
+
+##############################################################################
+# Plotting partial dependence of the two models independently
+# ===========================================================
+#
+# Next, we plot a partial dependence curves for features "LSTAT" and "RM" for
+# for the histogram gradient boosting model.
+fig, ax = plt.subplots(figsize=(12, 6))
+ax.set_title("Histogram Gradient Boosting")
+hgbr_disp = plot_partial_dependence(hgbr, X_test, ["LSTAT", "RM"],
+                                    feature_names=boston.feature_names, ax=ax)
+
+##############################################################################
+# The partial depdendence curves can be plotted for the multi-layer perceptron.
+# In this case `line_kw` was passed to
+# `~sklearn.inspection.plot_partial_dependence` to change the color of the
+# curve and `n_cols` was set to 1 to set the number of columns to 1.
+fig, ax = plt.subplots(figsize=(8, 8))
+ax.set_title("Multi-layer Perceptron")
+mlp_disp = plot_partial_dependence(mlp, X_test, ["LSTAT", "RM"],
+                                   feature_names=boston.feature_names, ax=ax,
+                                   n_cols=1, line_kw={"c": "red"})
+
+##############################################################################
+# Plotting partial dependence of the two models together
+# ======================================================
+#
+# The `hgbr_disp` and `mlp_disp`
+# :class:`~sklearn.inspection.PartialDependenceDisplay` objects contain all the
+# computed information needed to recreate the partial dependence curves. This
+# means we can easily create additional plots without needing to recompute the
+# curves.
+#
+# One way to plot the curves is to place them in the same figure, with the
+# curves of each model on each row. First, we create a figure two axes with two
+# rows and one column. The two axes are passed to the
+# :func:`~sklearn.inspection.PartialDependenceDisplay.plot` functions of
+# `hgbr_disp` and `mlp_disp`. The plot funciton will plot the two curves in the
path: examples/plot_partial_dependence_visualization_api.py
text: @user1:
funciton -> function

diff_hunk: @@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
+tht plot with the Visualization API.
path: examples/plot_partial_dependence_visualization_api.py
text: @user1:
```suggestion
the plot with the Visualization API.
```

diff_hunk: @@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
+tht plot with the Visualization API.
+
+.. note::
+
+    See also :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py`
+
+"""
+print(__doc__)
+
+import matplotlib.pyplot as plt
+from sklearn.datasets import load_boston
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.experimental import enable_hist_gradient_boosting  # noqa
+from sklearn.ensemble import HistGradientBoostingRegressor
+from sklearn.inspection import plot_partial_dependence
+
+
+##############################################################################
+# Train models on the boston housing price dataset
+# ================================================
+#
+# First, we load the boston housing price dataset and split the the dataset
+# into a training and test set. Then, we train a histogram gradient boosting
+# decision tree and a multi-layer perceptron on the training set.
+
+boston = load_boston()
+X, y = boston.data, boston.target
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1,
+                                                    random_state=0)
+
+hgbr = make_pipeline(StandardScaler(), HistGradientBoostingRegressor())
path: examples/plot_partial_dependence_visualization_api.py
text: @user1:
There is no need to use the StandardScaler. It will also speed-up the partial-dependence computation since we will use the recursive method included in the HistGradientBoostingRegressor which is not available when using a Pipeline.

diff_hunk: @@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
+tht plot with the Visualization API.
+
+.. note::
+
+    See also :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py`
+
+"""
+print(__doc__)
+
+import matplotlib.pyplot as plt
+from sklearn.datasets import load_boston
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.experimental import enable_hist_gradient_boosting  # noqa
+from sklearn.ensemble import HistGradientBoostingRegressor
+from sklearn.inspection import plot_partial_dependence
+
+
+##############################################################################
+# Train models on the boston housing price dataset
+# ================================================
+#
+# First, we load the boston housing price dataset and split the the dataset
+# into a training and test set. Then, we train a histogram gradient boosting
+# decision tree and a multi-layer perceptron on the training set.
+
+boston = load_boston()
+X, y = boston.data, boston.target
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1,
+                                                    random_state=0)
+
+hgbr = make_pipeline(StandardScaler(), HistGradientBoostingRegressor())
+mlp = make_pipeline(StandardScaler(),
+                    MLPRegressor(hidden_layer_sizes=(100, 100),
+                                 tol=1e-2, max_iter=500, random_state=0))
+hgbr.fit(X_train, y_train)
+mlp.fit(X_train, y_train)
+
+##############################################################################
+# Plotting partial dependence of the two models independently
+# ===========================================================
+#
+# Next, we plot a partial dependence curves for features "LSTAT" and "RM" for
+# for the histogram gradient boosting model.
+fig, ax = plt.subplots(figsize=(12, 6))
+ax.set_title("Histogram Gradient Boosting")
+hgbr_disp = plot_partial_dependence(hgbr, X_test, ["LSTAT", "RM"],
+                                    feature_names=boston.feature_names, ax=ax)
+
+##############################################################################
+# The partial depdendence curves can be plotted for the multi-layer perceptron.
+# In this case `line_kw` was passed to
+# `~sklearn.inspection.plot_partial_dependence` to change the color of the
+# curve and `n_cols` was set to 1 to set the number of columns to 1.
+fig, ax = plt.subplots(figsize=(8, 8))
+ax.set_title("Multi-layer Perceptron")
+mlp_disp = plot_partial_dependence(mlp, X_test, ["LSTAT", "RM"],
+                                   feature_names=boston.feature_names, ax=ax,
+                                   n_cols=1, line_kw={"c": "red"})
+
+##############################################################################
+# Plotting partial dependence of the two models together
+# ======================================================
+#
+# The `hgbr_disp` and `mlp_disp`
+# :class:`~sklearn.inspection.PartialDependenceDisplay` objects contain all the
+# computed information needed to recreate the partial dependence curves. This
+# means we can easily create additional plots without needing to recompute the
+# curves.
+#
+# One way to plot the curves is to place them in the same figure, with the
+# curves of each model on each row. First, we create a figure two axes with two
+# rows and one column. The two axes are passed to the
+# :func:`~sklearn.inspection.PartialDependenceDisplay.plot` functions of
+# `hgbr_disp` and `mlp_disp`. The plot funciton will plot the two curves in the
+# space by the passed in axes. The resulting plot places the histogram gradient
path: examples/plot_partial_dependence_visualization_api.py
text: @user1:
"space by the passed in axes" -> there is something wrong in the sentence

@user1:
"The given axes will be used by the plotting function to draw the PD"?

diff_hunk: @@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
+tht plot with the Visualization API.
+
+.. note::
+
+    See also :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py`
+
+"""
+print(__doc__)
+
+import matplotlib.pyplot as plt
+from sklearn.datasets import load_boston
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.experimental import enable_hist_gradient_boosting  # noqa
+from sklearn.ensemble import HistGradientBoostingRegressor
+from sklearn.inspection import plot_partial_dependence
+
+
+##############################################################################
+# Train models on the boston housing price dataset
+# ================================================
+#
+# First, we load the boston housing price dataset and split the the dataset
+# into a training and test set. Then, we train a histogram gradient boosting
+# decision tree and a multi-layer perceptron on the training set.
+
+boston = load_boston()
+X, y = boston.data, boston.target
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1,
+                                                    random_state=0)
+
+hgbr = make_pipeline(StandardScaler(), HistGradientBoostingRegressor())
+mlp = make_pipeline(StandardScaler(),
+                    MLPRegressor(hidden_layer_sizes=(100, 100),
+                                 tol=1e-2, max_iter=500, random_state=0))
+hgbr.fit(X_train, y_train)
+mlp.fit(X_train, y_train)
+
+##############################################################################
+# Plotting partial dependence of the two models independently
+# ===========================================================
+#
+# Next, we plot a partial dependence curves for features "LSTAT" and "RM" for
+# for the histogram gradient boosting model.
+fig, ax = plt.subplots(figsize=(12, 6))
+ax.set_title("Histogram Gradient Boosting")
+hgbr_disp = plot_partial_dependence(hgbr, X_test, ["LSTAT", "RM"],
+                                    feature_names=boston.feature_names, ax=ax)
+
+##############################################################################
+# The partial depdendence curves can be plotted for the multi-layer perceptron.
+# In this case `line_kw` was passed to
+# `~sklearn.inspection.plot_partial_dependence` to change the color of the
+# curve and `n_cols` was set to 1 to set the number of columns to 1.
+fig, ax = plt.subplots(figsize=(8, 8))
+ax.set_title("Multi-layer Perceptron")
+mlp_disp = plot_partial_dependence(mlp, X_test, ["LSTAT", "RM"],
+                                   feature_names=boston.feature_names, ax=ax,
+                                   n_cols=1, line_kw={"c": "red"})
+
+##############################################################################
+# Plotting partial dependence of the two models together
+# ======================================================
+#
+# The `hgbr_disp` and `mlp_disp`
+# :class:`~sklearn.inspection.PartialDependenceDisplay` objects contain all the
+# computed information needed to recreate the partial dependence curves. This
+# means we can easily create additional plots without needing to recompute the
+# curves.
+#
+# One way to plot the curves is to place them in the same figure, with the
+# curves of each model on each row. First, we create a figure two axes with two
+# rows and one column. The two axes are passed to the
+# :func:`~sklearn.inspection.PartialDependenceDisplay.plot` functions of
+# `hgbr_disp` and `mlp_disp`. The plot funciton will plot the two curves in the
+# space by the passed in axes. The resulting plot places the histogram gradient
+# boosting partial dependence curves on top of the multi-layer perceptron
path: examples/plot_partial_dependence_visualization_api.py
text: @user1:
what about:

The partial dependence plot of the GBDT and MLP models are overlaid on each other.

@user1:
Oh seeing the plot, you might want to use "above" instead of "top" :)

@user2:
"places the gbrt plots in the first row and the mlp plots in the second row"?

diff_hunk: @@ -120,7 +120,9 @@
 fig = plt.gcf()
 fig.suptitle('Partial dependence of house value on non-location features\n'
              'for the California housing dataset, with Gradient Boosting')
-fig.tight_layout(rect=[0, 0.03, 1, 0.95])
+fig.subplots_adjust(wspace=0.8, hspace=0.3)
+plt.show()
path: examples/inspection/plot_partial_dependence.py
text: @user1:
This will block the execution of the script until the figure is closed, isn't it?

@author:
Yup. I'll remove it. (looks like we usually just call it once)

diff_hunk: @@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
Indicies -> Indices

@user1:
int -> integer

diff_hunk: @@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
+        a partial dependence curve of one feature. A tuple of two ints will
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
ints -> integers

diff_hunk: @@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
+        a partial dependence curve of one feature. A tuple of two ints will
+        plot a two-way partial dependence curve as a contour plot.
+
+    feature_names : list of str
+        Feature names corrsponding to the indicies in ``features``.
+
+    target_idx : int
+        - In a multiclass setting, specifies the class for which the PDPs
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
PDP should be defined at the first occurence

diff_hunk: @@ -0,0 +1,110 @@
+"""
+=========================================
+Partial Dependence with Visualization API
+=========================================
+The :func:`~sklearn.inspection.plot_partial_dependence` function returns a
+:class:`~sklearn.inspection.PartialDependenceDisplay` object that can be used
+for plotting without needing to recalculate the partial dependence. In this
+example we should how to plot partial dependence plots and quickly customize
+tht plot with the Visualization API.
+
+.. note::
+
+    See also :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py`
+
+"""
+print(__doc__)
+
+import matplotlib.pyplot as plt
+from sklearn.datasets import load_boston
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.experimental import enable_hist_gradient_boosting  # noqa
+from sklearn.ensemble import HistGradientBoostingRegressor
+from sklearn.inspection import plot_partial_dependence
+
+
+##############################################################################
+# Train models on the boston housing price dataset
+# ================================================
+#
+# First, we load the boston housing price dataset and split the the dataset
+# into a training and test set. Then, we train a histogram gradient boosting
+# decision tree and a multi-layer perceptron on the training set.
+
+boston = load_boston()
+X, y = boston.data, boston.target
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1,
+                                                    random_state=0)
+
+hgbr = make_pipeline(StandardScaler(), HistGradientBoostingRegressor())
+mlp = make_pipeline(StandardScaler(),
+                    MLPRegressor(hidden_layer_sizes=(100, 100),
+                                 tol=1e-2, max_iter=500, random_state=0))
+hgbr.fit(X_train, y_train)
+mlp.fit(X_train, y_train)
+
+##############################################################################
+# Plotting partial dependence of the two models independently
+# ===========================================================
+#
+# Next, we plot a partial dependence curves for features "LSTAT" and "RM" for
+# for the histogram gradient boosting model.
+fig, ax = plt.subplots(figsize=(12, 6))
+ax.set_title("Histogram Gradient Boosting")
+hgbr_disp = plot_partial_dependence(hgbr, X_test, ["LSTAT", "RM"],
+                                    feature_names=boston.feature_names, ax=ax)
+
+##############################################################################
+# The partial depdendence curves can be plotted for the multi-layer perceptron.
+# In this case `line_kw` was passed to
+# `~sklearn.inspection.plot_partial_dependence` to change the color of the
+# curve and `n_cols` was set to 1 to set the number of columns to 1.
+fig, ax = plt.subplots(figsize=(8, 8))
+ax.set_title("Multi-layer Perceptron")
+mlp_disp = plot_partial_dependence(mlp, X_test, ["LSTAT", "RM"],
+                                   feature_names=boston.feature_names, ax=ax,
+                                   n_cols=1, line_kw={"c": "red"})
+
+##############################################################################
+# Plotting partial dependence of the two models together
+# ======================================================
+#
+# The `hgbr_disp` and `mlp_disp`
+# :class:`~sklearn.inspection.PartialDependenceDisplay` objects contain all the
+# computed information needed to recreate the partial dependence curves. This
+# means we can easily create additional plots without needing to recompute the
+# curves.
+#
+# One way to plot the curves is to place them in the same figure, with the
+# curves of each model on each row. First, we create a figure two axes with two
path: examples/plot_partial_dependence_visualization_api.py
text: @user1:
"a figure two axes" -> is it missing a "with"?

diff_hunk: @@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
+        a partial dependence curve of one feature. A tuple of two ints will
+        plot a two-way partial dependence curve as a contour plot.
+
+    feature_names : list of str
+        Feature names corrsponding to the indicies in ``features``.
+
+    target_idx : int
+        - In a multiclass setting, specifies the class for which the PDPs
+          should be computed. Note that for binary classification, the
+          positive class (index 1) is always used.
+        - In a multioutput setting, specifies the task for which the PDPs
+          should be computed
+        Ignored in binary classification or classical regression settings.
+
+    pdp_lim : dict
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
Shall we precise the type of thing we have in pdp_lim.
The description `pdp_lim[2]` looks like we have a list instead.

@user1:
Oh I see that this is a dict of where keys are int. It looks a bit weird.

@author:
Yea it is key by the number of feature, which means it is a dict with at most two elements.

diff_hunk: @@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
+        a partial dependence curve of one feature. A tuple of two ints will
+        plot a two-way partial dependence curve as a contour plot.
+
+    feature_names : list of str
+        Feature names corrsponding to the indicies in ``features``.
+
+    target_idx : int
+        - In a multiclass setting, specifies the class for which the PDPs
+          should be computed. Note that for binary classification, the
+          positive class (index 1) is always used.
+        - In a multioutput setting, specifies the task for which the PDPs
+          should be computed
+        Ignored in binary classification or classical regression settings.
+
+    pdp_lim : dict
+        Global min and max average predictions. `pdp_lim[1]` is the global min
+        and max for single partial dependence curves. `pdp_lim[2]` is the
+        global min and max for two-way partial dependence curves.
+
+    deciles : dict
+        Deciles for feature indicies in ``features``.
+
+    Attributes
+    ----------
+    bounding_ax_ : matplotlib Axes or None
+        If `ax` is an axes or None, the `bounding_ax_` is the axes where the
+        grid of partial dependence plots are drawn. If `ax` is  list of axes,
+        `bounding_ax_` is None.
+
+    axes_ : dict of matplotlib Axes
+        If `ax` is an axes or None, `axes_[i, j]` is the axes on the ith row
+        and jth column. If `ax` is a list of axes, `axes_[i]` is the ith item
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
jth -> j-th

@user1:
ith -> i-th

diff_hunk: @@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
+        a partial dependence curve of one feature. A tuple of two ints will
+        plot a two-way partial dependence curve as a contour plot.
+
+    feature_names : list of str
+        Feature names corrsponding to the indicies in ``features``.
+
+    target_idx : int
+        - In a multiclass setting, specifies the class for which the PDPs
+          should be computed. Note that for binary classification, the
+          positive class (index 1) is always used.
+        - In a multioutput setting, specifies the task for which the PDPs
+          should be computed
+        Ignored in binary classification or classical regression settings.
+
+    pdp_lim : dict
+        Global min and max average predictions. `pdp_lim[1]` is the global min
+        and max for single partial dependence curves. `pdp_lim[2]` is the
+        global min and max for two-way partial dependence curves.
+
+    deciles : dict
+        Deciles for feature indicies in ``features``.
+
+    Attributes
+    ----------
+    bounding_ax_ : matplotlib Axes or None
+        If `ax` is an axes or None, the `bounding_ax_` is the axes where the
+        grid of partial dependence plots are drawn. If `ax` is  list of axes,
+        `bounding_ax_` is None.
+
+    axes_ : dict of matplotlib Axes
+        If `ax` is an axes or None, `axes_[i, j]` is the axes on the ith row
+        and jth column. If `ax` is a list of axes, `axes_[i]` is the ith item
+        in `ax`.
+
+    lines_ : matplotlib Artists
+        If `ax` is an axes or None, `line_[i, j]` is the partial dependence
+        curve on the ith row and jth column. If `ax` is a list of axes,
+        `lines_[i]` is the partial dependence curve corresponding to the ith
+        item in `ax`. Only plots with a single feature will be populated in
+        `lines_`. Contour plots are stored in `contours_`.
+
+    contours_ : matplotlib Artists
+        If `ax` is an axes or None, `contours_[i, j]` is the partial dependence
+        plot on the ith row and jth column. If `ax` is a list of axes,
+        `contours_[i]` is the partial dependence plot corresponding to the ith
+        item in `ax`. Only plots with a two feature will be populated in
+        `contours_`. Line plots are stored in `lines_`.
+
+    figure_ : matplotlib Figure
+        Figure containing partial dependence plots.
+    """
+    def __init__(self, pd_results, features, feature_names, target_idx,
+                 pdp_lim, deciles):
+        self.pd_results = pd_results
+        self.features = features
+        self.feature_names = feature_names
+        self.target_idx = target_idx
+        self.pdp_lim = pdp_lim
+        self.deciles = deciles
+
+    def plot(self, ax=None, n_cols=3, line_kw=None, contour_kw=None, fig=None):
+        """Plot partial dependence plots.
+
+        Parameters
+        ----------
+        ax : Matplotlib axes, list of Matplotlib axes or None, default=None
+            Axes to plot the partial dependence curves. If a single axes is
+            given, it is treated as a bounding axes and a grid of partial
+            depdendence plots will be drawn on that top of it. If a list of
+            axes are passed, the partial dependence plots will be drawn on
+            those axes. By default, a single bounding axes is created and
+            treated as the single axes case.
+
+        n_cols : int, default=3
+            The maximum number of columns in the grid plot.
+
+        line_kw : dict or None, default=None
+            Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
+            For one-way partial dependence plots.
+
+        contour_kw : dict or None, default=None
+            Dict with keywords passed to the ``matplotlib.pyplot.contourf``
+            call for two-way partial dependence plots.
+
+        fig : Matplotlib figure object or None, default=None
+            A figure object onto which the plots will be drawn, after the
+            figure has been cleared. By default, a new one is created.
+
+            .. deprecated:: 0.22
+
+        Returns
+        -------
+        display: :class:`~sklearn.inspection.PartialDependenceDisplay`
+        """
+
+        check_matplotlib_support("plot_partial_dependence")
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
Since this is a pattern that we will have in each object for plotting, could we make a base class to handle that in a specific private method.

For instance a base class could call `check_matplotlib_support` and it could be an abstractmethod. The child class would have to implement it to import all matplotlib stuff. A naming could be somethign like `_validate_mpl_import`?

@author:
I would tend toward having a `sklearn.utils._plot._validate_mpl_import` and call that. (Trying to avoid inheritance for now)

What do you see `_validate_mpl_import` returning?

diff_hunk: @@ -0,0 +1,243 @@
+from itertools import count
+import warnings
+
+import numpy as np
+
+from ...utils import check_matplotlib_support
+
+
+class PartialDependenceDisplay:
+    """Partial Dependence visualization.
+
+    If is recommended to use
+    :func:`~sklearn.inspection.plot_partial_dependence` to create a
+    :class:`~sklearn.inspection.PartialDependenceDisplay`. All parameters are
+    stored as attributes.
+
+    Read more in
+    :ref:`sphx_glr_auto_examples_plot_roc_curve_visualization_api.py` and the
+    :ref:`User Guide <visualizations>`.
+
+    Parameters
+    ----------
+    pd_results : list of (ndarray, ndarray)
+        Results of `sklearn.inspection.partial_dependence` for ``features``.
+
+    features : list of {(int, ), (int, int)}
+        Indicies of features for a given plot. A tuple of one int will plot
+        a partial dependence curve of one feature. A tuple of two ints will
+        plot a two-way partial dependence curve as a contour plot.
+
+    feature_names : list of str
+        Feature names corrsponding to the indicies in ``features``.
+
+    target_idx : int
+        - In a multiclass setting, specifies the class for which the PDPs
+          should be computed. Note that for binary classification, the
+          positive class (index 1) is always used.
+        - In a multioutput setting, specifies the task for which the PDPs
+          should be computed
+        Ignored in binary classification or classical regression settings.
+
+    pdp_lim : dict
+        Global min and max average predictions. `pdp_lim[1]` is the global min
+        and max for single partial dependence curves. `pdp_lim[2]` is the
+        global min and max for two-way partial dependence curves.
+
+    deciles : dict
+        Deciles for feature indicies in ``features``.
+
+    Attributes
+    ----------
+    bounding_ax_ : matplotlib Axes or None
+        If `ax` is an axes or None, the `bounding_ax_` is the axes where the
+        grid of partial dependence plots are drawn. If `ax` is  list of axes,
+        `bounding_ax_` is None.
+
+    axes_ : dict of matplotlib Axes
+        If `ax` is an axes or None, `axes_[i, j]` is the axes on the ith row
+        and jth column. If `ax` is a list of axes, `axes_[i]` is the ith item
+        in `ax`.
+
+    lines_ : matplotlib Artists
+        If `ax` is an axes or None, `line_[i, j]` is the partial dependence
+        curve on the ith row and jth column. If `ax` is a list of axes,
+        `lines_[i]` is the partial dependence curve corresponding to the ith
+        item in `ax`. Only plots with a single feature will be populated in
+        `lines_`. Contour plots are stored in `contours_`.
+
+    contours_ : matplotlib Artists
+        If `ax` is an axes or None, `contours_[i, j]` is the partial dependence
+        plot on the ith row and jth column. If `ax` is a list of axes,
+        `contours_[i]` is the partial dependence plot corresponding to the ith
+        item in `ax`. Only plots with a two feature will be populated in
+        `contours_`. Line plots are stored in `lines_`.
+
+    figure_ : matplotlib Figure
+        Figure containing partial dependence plots.
+    """
+    def __init__(self, pd_results, features, feature_names, target_idx,
+                 pdp_lim, deciles):
+        self.pd_results = pd_results
+        self.features = features
+        self.feature_names = feature_names
+        self.target_idx = target_idx
+        self.pdp_lim = pdp_lim
+        self.deciles = deciles
+
+    def plot(self, ax=None, n_cols=3, line_kw=None, contour_kw=None, fig=None):
+        """Plot partial dependence plots.
+
+        Parameters
+        ----------
+        ax : Matplotlib axes, list of Matplotlib axes or None, default=None
+            Axes to plot the partial dependence curves. If a single axes is
+            given, it is treated as a bounding axes and a grid of partial
+            depdendence plots will be drawn on that top of it. If a list of
+            axes are passed, the partial dependence plots will be drawn on
+            those axes. By default, a single bounding axes is created and
+            treated as the single axes case.
+
+        n_cols : int, default=3
+            The maximum number of columns in the grid plot.
+
+        line_kw : dict or None, default=None
+            Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
+            For one-way partial dependence plots.
+
+        contour_kw : dict or None, default=None
+            Dict with keywords passed to the ``matplotlib.pyplot.contourf``
+            call for two-way partial dependence plots.
+
+        fig : Matplotlib figure object or None, default=None
+            A figure object onto which the plots will be drawn, after the
+            figure has been cleared. By default, a new one is created.
+
+            .. deprecated:: 0.22
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
It is weird to have something deprecated in a new object.

@user2:
I agree. We need it for backward-compatibility though. We could make it private. I'm not sure if that's less weird?

@user3:
Should we just not publicly document it then?

diff_hunk: @@ -598,71 +620,12 @@ def convert_feature(fx):
         max_pd = max(max_pd, old_max_pd)
         pdp_lim[n_fx] = (min_pd, max_pd)
 
-    # create contour levels for two-way plots
-    if 2 in pdp_lim:
-        Z_level = np.linspace(*pdp_lim[2], num=8)
-
-    if fig is None:
-        fig = plt.figure()
-    else:
-        fig.clear()
-
-    if line_kw is None:
-        line_kw = {'color': 'green'}
-    if contour_kw is None:
-        contour_kw = {}
-
-    n_cols = min(n_cols, len(features))
-    n_rows = int(np.ceil(len(features) / float(n_cols)))
-    axs = []
-    for i, fx, name, (avg_preds, values) in zip(
-            count(), features, names, pd_result):
-        ax = fig.add_subplot(n_rows, n_cols, i + 1)
-
-        if len(values) == 1:
-            ax.plot(values[0], avg_preds[target_idx].ravel(), **line_kw)
-        else:
-            # make contour plot
-            assert len(values) == 2
-            XX, YY = np.meshgrid(values[0], values[1])
-            Z = avg_preds[target_idx].T
-            CS = ax.contour(XX, YY, Z, levels=Z_level, linewidths=0.5,
-                            colors='k')
-            ax.contourf(XX, YY, Z, levels=Z_level, vmax=Z_level[-1],
-                        vmin=Z_level[0], alpha=0.75, **contour_kw)
-            ax.clabel(CS, fmt='%2.2f', colors='k', fontsize=10, inline=True)
-
-        # plot data deciles + axes labels
-        deciles = mquantiles(X[:, fx[0]], prob=np.arange(0.1, 1.0, 0.1))
-        trans = transforms.blended_transform_factory(ax.transData,
-                                                     ax.transAxes)
-        ylim = ax.get_ylim()
-        ax.vlines(deciles, [0], 0.05, transform=trans, color='k')
-        ax.set_xlabel(name[0])
-        ax.set_ylim(ylim)
-
-        # prevent x-axis ticks from overlapping
-        ax.xaxis.set_major_locator(MaxNLocator(nbins=6, prune='lower'))
-        tick_formatter = ScalarFormatter()
-        tick_formatter.set_powerlimits((-3, 4))
-        ax.xaxis.set_major_formatter(tick_formatter)
-
-        if len(values) > 1:
-            # two-way PDP - y-axis deciles + labels
-            deciles = mquantiles(X[:, fx[1]], prob=np.arange(0.1, 1.0, 0.1))
-            trans = transforms.blended_transform_factory(ax.transAxes,
-                                                         ax.transData)
-            xlim = ax.get_xlim()
-            ax.hlines(deciles, [0], 0.05, transform=trans, color='k')
-            ax.set_ylabel(name[1])
-            # hline erases xlim
-            ax.set_xlim(xlim)
-        else:
-            ax.set_ylabel('Partial dependence')
-
-        if len(values) == 1:
-            ax.set_ylim(pdp_lim[1])
-        axs.append(ax)
+    deciles = {}
+    for fx in chain.from_iterable(features):
+        if fx not in deciles:
+            deciles[fx] = mquantiles(X[:, fx], prob=np.arange(0.1, 1.0, 0.1))
 
-    fig.subplots_adjust(bottom=0.15, top=0.7, left=0.1, right=0.95, wspace=0.4,
-                        hspace=0.3)
+    display = PartialDependenceDisplay(pd_results, features, feature_names,
+                                       target_idx, pdp_lim, deciles)
+    return display.plot(ax=ax, n_cols=n_cols, line_kw=line_kw,
+                        contour_kw=contour_kw, fig=fig)
path: sklearn/inspection/partial_dependence.py
text: @user1:
Regarding `fig`, could we set `display.figure_` instead avoiding to have the `fig` in the new class. We can raise the deprecation warning in this function directly as well.

@user1:
Regarding the introduction of `fig`

diff_hunk: @@ -448,23 +458,37 @@ def plot_partial_dependence(estimator, X, features, feature_names=None,
         - 'auto':
           - 'recursion' is used for estimators that supports it.
           - 'brute' is used for all other estimators.
+
     n_jobs : int, optional (default=None)
         The number of CPUs to use to compute the partial dependences.
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
         for more details.
+
     verbose : int, optional (default=0)
         Verbose output during PD computations.
+
     fig : Matplotlib figure object, optional (default=None)
         A figure object onto which the plots will be drawn, after the figure
         has been cleared. By default, a new one is created.
+
+        .. deprecated:: 0.22
+
     line_kw : dict, optional
         Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
         For one-way partial dependence plots.
+
     contour_kw : dict, optional
-        Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
+        Dict with keywords passed to the ``matplotlib.pyplot.contourf`` call.
         For two-way partial dependence plots.
 
+    ax : Matplotlib axes, list of Matplotlib axes or None, (default=None)
+        By default,
path: sklearn/inspection/partial_dependence.py
text: @user1:
Missing the rest of the sentence ;)

diff_hunk: @@ -0,0 +1,325 @@
+import numpy as np
+from scipy.stats.mstats import mquantiles
+
+import pytest
+from numpy.testing import assert_allclose
+
+from sklearn.datasets import load_boston
+from sklearn.datasets import load_iris
+from sklearn.datasets import make_classification, make_regression
+from sklearn.ensemble import GradientBoostingRegressor
+from sklearn.ensemble import GradientBoostingClassifier
+from sklearn.linear_model import LinearRegression
+from sklearn.inspection import plot_partial_dependence
+
+
+boston = load_boston()
path: sklearn/inspection/tests/test_plot_partial_dependence.py
text: @user1:
it would be better to name it BOSTON.

Maybe this is time to use a fixture for the data here instead?

diff_hunk: @@ -448,23 +458,37 @@ def plot_partial_dependence(estimator, X, features, feature_names=None,
         - 'auto':
           - 'recursion' is used for estimators that supports it.
           - 'brute' is used for all other estimators.
+
     n_jobs : int, optional (default=None)
         The number of CPUs to use to compute the partial dependences.
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
         for more details.
+
     verbose : int, optional (default=0)
         Verbose output during PD computations.
+
     fig : Matplotlib figure object, optional (default=None)
         A figure object onto which the plots will be drawn, after the figure
         has been cleared. By default, a new one is created.
+
+        .. deprecated:: 0.22
+
     line_kw : dict, optional
         Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
         For one-way partial dependence plots.
+
     contour_kw : dict, optional
-        Dict with keywords passed to the ``matplotlib.pyplot.plot`` call.
+        Dict with keywords passed to the ``matplotlib.pyplot.contourf`` call.
         For two-way partial dependence plots.
 
+    ax : Matplotlib axes, list of Matplotlib axes or None, (default=None)
path: sklearn/inspection/partial_dependence.py
text: @user1:
Let's go for option 3

```python
ax : Matplotlib axes, list of Matplotlib axes, default=None
```