Repository: dask/dask-ml

Issue:
Behavior of regression metrics when multioutput=None is inconsistent with scikit-learn
**What happened**:

Using `multioutput=None` with the regression metrics in `dask_ml.metrics.regression` results in an error.

**What you expected to happen**:

I expected the behavior to be the same as the equivalent `scikit-learn` metrics functions, where `multioutput=None` means "all elements have the same weight".

**Minimal Complete Verifiable Example**:

In `scikit-learn`, the value of `multioutput` is passed through to `np.average()`.

https://github.com/scikit-learn/scikit-learn/blob/877c6e6db42006445ccf0695c0dde3294ff4dd4a/sklearn/metrics/_regression.py#L195

`np.average()` treats the value `None` as "equally-weighted", which is the same as passing `multioutput = "uniform_average"`. From https://numpy.org/doc/stable/reference/generated/numpy.average.html

> If weights=None, then all data in a are assumed to have a weight equal to one. 

```python
import numpy as np
from sklearn.metrics import mean_squared_error

a = np.random.uniform(size=(100, 3))
b = np.random.uniform(size=(100, 3))

raw_output = mean_squared_error(a, b, multioutput="raw_values")
print(raw_output)

# [0.17118814 0.15964742 0.13095381]

mean_squared_error(a, b, multioutput=None)

# 0.15392978995168105
```

In `dask-ml`, passing `multioutput=None` results in an error.

```python
import dask.array as da
from dask_ml.metrics.regression import mean_squared_error

a = da.random.uniform(size=(100, 3))
b = da.random.uniform(size=(100, 3))

raw_output = mean_squared_error(a, b, multioutput="raw_values")
print(raw_output)
# dask.array<mean_agg-aggregate, shape=(3,), dtype=float64, chunksize=(3,), chunktype=numpy.ndarray>

mean_squared_error(a, b, multioutput=None)
# ValueError: Weighted 'multioutput' not supported.
```

**Anything else we need to know?**:

**Environment**:

- Dask version:

<details><summary>output of 'conda list | grep -E "dask|distributed"'</summary>

```text
dask                      2021.3.0           pyhd3eb1b0_0  
dask-core                 2021.3.0           pyhd3eb1b0_0  
dask-glm                  0.2.0                    pypi_0    pypi
dask-saturn               0.2.2                    pypi_0    pypi
dask-sphinx-theme         1.3.5                    pypi_0    pypi
distributed               2021.3.0         py37h06a4308_0  
```

</details>

- Python version, Operating System

<details><summary>output of 'conda info'</summary>

```text

     active environment : None
       user config file : /home/jlamb/.condarc
 populated config files : 
          conda version : 4.9.2
    conda-build version : not installed
         python version : 3.7.6.final.0
       virtual packages : __glibc=2.27=0
                          __unix=0=0
                          __archspec=1=x86_64
       base environment : /home/jlamb/miniconda3  (writable)
           channel URLs : https://repo.anaconda.com/pkgs/main/linux-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/r/linux-64
                          https://repo.anaconda.com/pkgs/r/noarch
          package cache : /home/jlamb/miniconda3/pkgs
                          /home/jlamb/.conda/pkgs
       envs directories : /home/jlamb/miniconda3/envs
                          /home/jlamb/.conda/envs
               platform : linux-64
             user-agent : conda/4.9.2 requests/2.22.0 CPython/3.7.6 Linux/5.4.0-70-generic ubuntu/18.04.4 glibc/2.27
                UID:GID : 1000:1000
             netrc file : None
           offline mode : False
```

</details>

- Install method (conda, pip, source): installed from source as of https://github.com/dask/dask-ml/commit/f5e5bb4d4d21782b0e93f0f6541e6f2501b0c06c, using `python setup.py install`.


Discussion:

Patch:
diff --git a/ci/posix.yaml b/ci/posix.yaml
index 8379f907..2904cd12 100644
--- a/ci/posix.yaml
+++ b/ci/posix.yaml
@@ -10,13 +10,13 @@ jobs:
     matrix:
       linux37:
         envFile: 'ci/environment-3.7.yaml'
-        SKLARN_DEV: "no"
+        SKLEARN_DEV: "no"
       linux38:
         envFile: 'ci/environment-3.8.yaml'
-        SKLARN_DEV: "no"
+        SKLEARN_DEV: "no"
       earliest:
         envFile: 'ci/environment-3.6.yaml'
-        SKLARN_DEV: "no"
+        SKLEARN_DEV: "no"
       sklearnDev:
         envFile: 'ci/environment-3.7.yaml'
         SKLEARN_DEV: "yes"
diff --git a/dask_ml/metrics/__init__.py b/dask_ml/metrics/__init__.py
index 05349d10..32b85f72 100644
--- a/dask_ml/metrics/__init__.py
+++ b/dask_ml/metrics/__init__.py
@@ -6,6 +6,7 @@ from .pairwise import (  # noqa
 )
 from .regression import (  # noqa
     mean_absolute_error,
+    mean_absolute_percentage_error,
     mean_squared_error,
     mean_squared_log_error,
     r2_score,
diff --git a/dask_ml/metrics/regression.py b/dask_ml/metrics/regression.py
index d1849ef8..c48c7629 100644
--- a/dask_ml/metrics/regression.py
+++ b/dask_ml/metrics/regression.py
@@ -16,7 +16,7 @@ def _check_sample_weight(sample_weight: Optional[ArrayLike]):
 def _check_reg_targets(
     y_true: ArrayLike, y_pred: ArrayLike, multioutput: Optional[str]
 ):
-    if multioutput != "uniform_average":
+    if multioutput is not None and multioutput != "uniform_average":
         raise NotImplementedError("'multioutput' must be 'uniform_average'")
 
     if y_true.ndim == 1:
@@ -40,12 +40,12 @@ def mean_squared_error(
     _check_sample_weight(sample_weight)
     output_errors = ((y_pred - y_true) ** 2).mean(axis=0)
 
-    if isinstance(multioutput, str):
+    if isinstance(multioutput, str) or multioutput is None:
         if multioutput == "raw_values":
-            return output_errors
-        elif multioutput == "uniform_average":
-            # pass None as weights to np.average: uniform mean
-            multioutput = None
+            if compute:
+                return output_errors.compute()
+            else:
+                return output_errors
     else:
         raise ValueError("Weighted 'multioutput' not supported.")
     result = output_errors.mean()
@@ -67,12 +67,75 @@ def mean_absolute_error(
     _check_sample_weight(sample_weight)
     output_errors = abs(y_pred - y_true).mean(axis=0)
 
-    if isinstance(multioutput, str):
+    if isinstance(multioutput, str) or multioutput is None:
         if multioutput == "raw_values":
-            return output_errors
-        elif multioutput == "uniform_average":
-            # pass None as weights to np.average: uniform mean
-            multioutput = None
+            if compute:
+                return output_errors.compute()
+            else:
+                return output_errors
+    else:
+        raise ValueError("Weighted 'multioutput' not supported.")
+    result = output_errors.mean()
+    if compute:
+        result = result.compute()
+    return result
+
+
+def mean_absolute_percentage_error(
+    y_true: ArrayLike,
+    y_pred: ArrayLike,
+    sample_weight: Optional[ArrayLike] = None,
+    multioutput: Optional[str] = "uniform_average",
+    compute: bool = True,
+) -> ArrayLike:
+    """Mean absolute percentage error regression loss.
+
+    Note here that we do not represent the output as a percentage in range
+    [0, 100]. Instead, we represent it in range [0, 1/eps]. Read more in
+    https://scikit-learn.org/stable/modules/model_evaluation.html#mean-absolute-percentage-error
+
+    Parameters
+    ----------
+    y_true : array-like of shape (n_samples,) or (n_samples, n_outputs)
+        Ground truth (correct) target values.
+    y_pred : array-like of shape (n_samples,) or (n_samples, n_outputs)
+        Estimated target values.
+    sample_weight : array-like of shape (n_samples,), default=None
+        Sample weights.
+    multioutput : {'raw_values', 'uniform_average'} or array-like
+        Defines aggregating of multiple output values.
+        Array-like value defines weights used to average errors.
+        If input is list then the shape must be (n_outputs,).
+        'raw_values' :
+            Returns a full set of errors in case of multioutput input.
+        'uniform_average' :
+            Errors of all outputs are averaged with uniform weight.
+    compute : bool
+        Whether to compute this result (default ``True``)
+
+    Returns
+    -------
+    loss : float or array-like of floats in the range [0, 1/eps]
+        If multioutput is 'raw_values', then mean absolute percentage error
+        is returned for each output separately.
+        If multioutput is 'uniform_average' or ``None``, then the
+        equally-weighted average of all output errors is returned.
+        MAPE output is non-negative floating point. The best value is 0.0.
+        But note the fact that bad predictions can lead to arbitarily large
+        MAPE values, especially if some y_true values are very close to zero.
+        Note that we return a large value instead of `inf` when y_true is zero.
+    """
+    _check_sample_weight(sample_weight)
+    epsilon = np.finfo(np.float64).eps
+    mape = abs(y_pred - y_true) / da.maximum(y_true, epsilon)
+    output_errors = mape.mean(axis=0)
+
+    if isinstance(multioutput, str) or multioutput is None:
+        if multioutput == "raw_values":
+            if compute:
+                return output_errors.compute()
+            else:
+                return output_errors
     else:
         raise ValueError("Weighted 'multioutput' not supported.")
     result = output_errors.mean()
@@ -90,7 +153,7 @@ def r2_score(
     compute: bool = True,
 ) -> ArrayLike:
     _check_sample_weight(sample_weight)
-    _, y_true, y_pred, multioutput = _check_reg_targets(y_true, y_pred, multioutput)
+    _, y_true, y_pred, _ = _check_reg_targets(y_true, y_pred, multioutput)
     weight = 1.0
 
     numerator = (weight * (y_true - y_pred) ** 2).sum(axis=0, dtype="f8")
diff --git a/docs/source/modules/api.rst b/docs/source/modules/api.rst
index 601357ba..7a2d4d06 100644
--- a/docs/source/modules/api.rst
+++ b/docs/source/modules/api.rst
@@ -245,6 +245,7 @@ Regression Metrics
    :toctree: generated/
 
    metrics.mean_absolute_error
+   metrics.mean_absolute_percentage_error
    metrics.mean_squared_error
    metrics.mean_squared_log_error
    metrics.r2_score


Test patch:
diff --git a/tests/metrics/test_regression.py b/tests/metrics/test_regression.py
index dfdc5480..62f0a9d7 100644
--- a/tests/metrics/test_regression.py
+++ b/tests/metrics/test_regression.py
@@ -1,18 +1,32 @@
 import numbers
 
 import dask.array as da
+import numpy as np
 import pytest
 import sklearn.metrics
+from dask.array.utils import assert_eq
 
 import dask_ml.metrics
+from dask_ml._compat import SK_024
 
+_METRICS_TO_TEST = [
+    "mean_squared_error",
+    "mean_absolute_error",
+    "r2_score",
+]
 
-@pytest.fixture(params=["mean_squared_error", "mean_absolute_error", "r2_score"])
+# mean_absolute_percentage_error() was added in scikit-learn 0.24.0
+if SK_024:
+    _METRICS_TO_TEST.append("mean_absolute_percentage_error")
+
+
+@pytest.fixture(params=_METRICS_TO_TEST)
 def metric_pairs(request):
     """Pairs of (dask-ml, sklearn) regression metrics.
 
     * mean_squared_error
     * mean_absolute_error
+    * mean_absolute_percentage_error (if scikit-learn >= 0.24.0)
     * r2_score
     """
     return (
@@ -60,3 +74,53 @@ def test_mean_squared_log_error():
     result = m1(a, b)
     expected = m2(a, b)
     assert abs(result - expected) < 1e-5
+
+
+@pytest.mark.parametrize("multioutput", ["uniform_average", None])
+def test_regression_metrics_unweighted_average_multioutput(metric_pairs, multioutput):
+    m1, m2 = metric_pairs
+
+    a = da.random.uniform(size=(100,), chunks=(25,))
+    b = da.random.uniform(size=(100,), chunks=(25,))
+
+    result = m1(a, b, multioutput=multioutput)
+    expected = m2(a, b, multioutput=multioutput)
+    assert abs(result - expected) < 1e-5
+
+
+@pytest.mark.parametrize("compute", [True, False])
+def test_regression_metrics_raw_values(metric_pairs, compute):
+    m1, m2 = metric_pairs
+
+    if m1.__name__ == "r2_score":
+        pytest.skip("r2_score does not support multioutput='raw_values'")
+
+    a = da.random.uniform(size=(100, 3), chunks=(25, 3))
+    b = da.random.uniform(size=(100, 3), chunks=(25, 3))
+
+    result = m1(a, b, multioutput="raw_values", compute=compute)
+    expected = m2(a, b, multioutput="raw_values")
+
+    if compute:
+        assert isinstance(result, np.ndarray)
+    else:
+        assert isinstance(result, da.Array)
+
+    assert_eq(result, expected)
+    assert result.shape == (3,)
+
+
+def test_regression_metrics_do_not_support_weighted_multioutput(metric_pairs):
+    m1, _ = metric_pairs
+
+    a = da.random.uniform(size=(100, 3), chunks=(25, 3))
+    b = da.random.uniform(size=(100, 3), chunks=(25, 3))
+    weights = da.random.uniform(size=(3,))
+
+    if m1.__name__ == "r2_score":
+        error_msg = "'multioutput' must be 'uniform_average'"
+    else:
+        error_msg = "Weighted 'multioutput' not supported."
+
+    with pytest.raises((NotImplementedError, ValueError), match=error_msg):
+        _ = m1(a, b, multioutput=weights)