Repository: scikit-learn/scikit-learn
Pull request: FIX compute `y_std` properly with multi-target in GPR
#### Reference Issues/PRs

Fixes [#17394](https://github.com/scikit-learn/scikit-learn/issues/17394). Fixes [#18065](https://github.com/scikit-learn/scikit-learn/issues/18065).



#### What does this implement/fix? Explain your changes.

GaussianProcessRegressor could not retrieve `y_std` 
for `predict(X, return_std=True)` when n_targets were bigger than 1. 
This happened because line 415 in file
"sklearn/gaussian_process/_gpr.py" tried to multiply
`y_var * self._y_train_std ** 2` using simple multiplication ( a1 * a2).

However, it fails when `self._y_train_std` has more than one feature
(when `n_targets` is more than 1), so we need to implement this
multiplication using `np.outer` product, because it will handle
the conventional scalar-array multiplication for EACH output feature
(`self._y_train_std` contains one normalization rate for each output
feature).

After this fix, each output feature will receive its respective 
normalized variance when `return_std == True`.

A simple example of how to train gpr and reproduce the error:
```
#=====breaking-code-before-fix============================================

import numpy as np

X_train = np.random.rand((11, 10))
y_train = np.random.rand((11, 6)) # 6 target features == ERROR

X_test = np.random.rand((4, 10))

import sklearn.gaussian_process as gp

kernel = gp.kernels.ConstantKernel(1.0, (1e-1, 1e3)) * gp.kernels.RBF(10.0, (1e-3, 1e3))

model = gp.GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, alpha=0.1, normalize_y=True)

model.fit(X_train, y_train)
params = model.kernel_.get_params()

y_pred, std = model.predict(X_test, return_std=True)

#=====end-of-breaking-code-before-fix======================================


#=====well-running-code-before-fix============================================

import numpy as np

X_train = np.random.rand((11, 1))
y_train = np.random.rand((11, 1)) # 1 target feature == FINE

X_test = np.random.rand((4, 1))

import sklearn.gaussian_process as gp

kernel = gp.kernels.ConstantKernel(1.0, (1e-1, 1e3)) * gp.kernels.RBF(10.0, (1e-3, 1e3))

model = gp.GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, alpha=0.1, normalize_y=True)

model.fit(X_train, y_train)
params = model.kernel_.get_params()

y_pred, std = model.predict(X_test, return_std=True)

#=====end-of-well-running-code-before-fix======================================
```

Line 415 in file "sklearn/gaussian_process/_gpr.py":
```
# before
y_var = y_var * self._y_train_std ** 2

# after
y_var = np.outer(y_var, self._y_train_std ** 2)
```




<!--

#### Any other comments?
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:
GaussianProcessRegressor cannot correctly predict std in multi-target scene
I used GaussianProcessRegressor like this:
```python
from sklearn.gaussian_process.kernels import Matern
from sklearn.gaussian_process import GaussianProcessRegressor
import numpy as np

ls = [
    ([2, 3, 4, 5], [1, 1, 3]),
    ([5, 3, 1, 3], [1, 4, 3]),
    ([2, 7, 2, 5], [8, 2, 3]),
    ([-2, 3, 4, 5], [1, 0, 4]),
    ([2, 8, 7, 5], [-1, 2, 3]),
]

X = np.array([ele[0] for ele in ls])
y = np.array([ele[1] for ele in ls])

gp = GaussianProcessRegressor(
    kernel=Matern(nu=2.5),
    alpha=1e-6,
    normalize_y=True,
    n_restarts_optimizer=5,
    random_state=np.random.RandomState(None),
)

gp.fit(X, y)
# x = np.array([200, 300, 600, 500.]).reshape(1, -1)
x = np.array([
    [2, 3, 6, 5.],
    [2, 3, 3, 5.],
    [5, 3, 1, 3.00001],
    [5, 3, 1, 3]
])

v, d = gp.predict(x, return_std=True)

print(v, d)

_, c = gp.predict(x, return_cov=True)

print(c)
```
But I got this Error:

```python
ValueError: operands could not be broadcast together with shapes (4,) (3,) 
```

I scaned the code of `GaussianProcessRegressor`,  and I found out what happened in function `predict(self, X, return_std=False, return_cov=False)`:
- code below: y_var.shape==(4,) but self._y_train_std.shape==(3,)
```
# undo normalisation
y_var = y_var * self._y_train_std ** 2
```
I think `y_var = y_var * self._y_train_std ** 2` will work well if there is only one target.
In multi-target scene, we should change it like this:
```python
# undo normalisation
# y_var = y_var * self._y_train_std ** 2
y_var = y_var.reshape((-1, 1))
y_var = np.einsum("ij,j->ij", y_var, self._y_train_std ** 2)
```
PS:
I think we need to add an `EPS` to the `self._y_train_std` to avoid zero division error.
```python
# Normalize target value
if self.normalize_y:
    self._y_train_mean = np.mean(y, axis=0)
    # self._y_train_std = np.std(y, axis=0)
    self._y_train_std = np.std(y, axis=0) + self.EPS # avoid zero division error

    # Remove mean and make unit variance
    y = (y - self._y_train_mean) / self._y_train_std
```

GaussianProcessRegressor doesn't work with multidemensional output when normalize_y=True
The GaussianProcessRegressor doesn't work with multidemensional output when normalize_y=True. In this example, the code runs fine when `normalize_y` is False, but breaks when it is true:

```py
    import numpy as np
    from sklearn.gaussian_process import GaussianProcessRegressor
    from sklearn.gaussian_process.kernels import RBF
    
    def f(x): return(np.array([np.sin(7 * x), x ** 4]))
    
    kernel = RBF()
    gp = GaussianProcessRegressor(kernel=RBF(length_scale=15.7), n_restarts_optimizer=50,
                             normalize_y=True) # (works when normalize_y is False)
    
    X = np.linspace(0, 5, 5)
    gp.fit(np.atleast_2d(X).T, f(X).T)
    
    newx = np.atleast_2d([1, 2, 3, 4]).T
    gp.predict(newx, return_std=False)
    gp.predict(newx, return_std=True)
```

Hints from issue discussion:
Hello, can I work on this ?
what version are you using?

I was able to reproduce this in the following environment.

```
>>> import sklearn
>>> sklearn.show_versions() 

System:
    python: 3.6.9 (default, Apr 18 2020, 01:56:04)  [GCC 8.4.0]
executable: /usr/bin/python3
   machine: Linux-4.19.104+-x86_64-with-Ubuntu-18.04-bionic

Python dependencies:
          pip: 19.3.1
   setuptools: 47.3.1
      sklearn: 0.23.1
        numpy: 1.18.5
        scipy: 1.4.1
       Cython: 0.29.20
       pandas: 1.0.5
   matplotlib: 3.2.2
       joblib: 0.15.1
threadpoolctl: 2.1.0

Built with OpenMP: True
```
I tried to run this code and its working fine for me.
@JK186 what environment are you using?
take
Hello. Can anyone consult me, whether it is correct that multitarget gpr should provide std of form (n_samples, n_output_dims) instead of (n_samples,), whereas cov should be of form (n_samples, n_samples, n_output_dims) instead of (n_samples, n_samples)?
I discovered it in Rasmussen and Williams' [book](http://www.gaussianprocess.org/gpml/chapters/RW.pdf): when you try to work with multiple outputs in GP and the classes are independent from each other, then you just need to work with each class separately: 

> "For example in the robot inverse dynamics problem described in section 2.5 there are really seven torques to be predicted. A simple approach is to model each output variable as independent from the others and treat them separately." (Rasmussen and Williams, 2004. p. 190) 

If so, you can easily divide your `y` into three arrays (per class) and work with each of them separately. I am currently working on the pull request where I just work with each target class separately, assuming that each class is independent from each other. 

However, as was stated in (Rasmussen and Williams, 2004. p. 190) there might be correlations caused by the noise process correlation or if the prior already have such structure so that cross-correlation between classes exists. Do I need to handle this also? 
take
Is it correct that multitarget regression predicts y_std of form (n_samples, n_output_dims) and y_cov of form (n_samples, n_samples, n_output_dims)?

Commit message: fix changelog

Patch to review:
diff --git a/doc/whats_new/v1.0.rst b/doc/whats_new/v1.0.rst
index 40076c0b275dd..d9686ac26af7a 100644
--- a/doc/whats_new/v1.0.rst
+++ b/doc/whats_new/v1.0.rst
@@ -43,6 +43,14 @@ Fixed models
   between sparse and dense input. :pr:`21195`
   by :user:`Jérémie du Boisberranger <jeremiedbb>`.
 
+:mod:`sklearn.gaussian_process`
+...............................
+
+- |Fix| Compute `y_std` properly with multi-target in
+  :class:`sklearn.gaussian_process.GaussianProcessRegressor` allowing
+  proper normalization in multi-target scene.
+  :pr:`20761` by :user:`Patrick de C. T. R. Ferreira <patrickctrf>`.
+
 :mod:`sklearn.linear_model`
 ...........................
 
diff --git a/sklearn/gaussian_process/_gpr.py b/sklearn/gaussian_process/_gpr.py
index cf9a349f1b074..90b313d7b1fcf 100644
--- a/sklearn/gaussian_process/_gpr.py
+++ b/sklearn/gaussian_process/_gpr.py
@@ -349,7 +349,7 @@ def predict(self, X, return_std=False, return_cov=False):
         y_mean : ndarray of shape (n_samples,) or (n_samples, n_targets)
             Mean of predictive distribution a query points.
 
-        y_std : ndarray of shape (n_samples,), optional
+        y_std : ndarray of shape (n_samples,) or (n_samples, n_targets), optional
             Standard deviation of predictive distribution at query points.
             Only returned when `return_std` is True.
 
@@ -424,7 +424,11 @@ def predict(self, X, return_std=False, return_cov=False):
                     y_var[y_var_negative] = 0.0
 
                 # undo normalisation
-                y_var = y_var * self._y_train_std ** 2
+                y_var = np.outer(y_var, self._y_train_std ** 2)
+
+                # if y_var has shape (n_samples, 1), reshape to (n_samples,)
+                if y_var.shape[0] == y_var.size:
+                    y_var = np.ravel(y_var)
 
                 return y_mean, np.sqrt(y_var)
             else:
diff --git a/sklearn/gaussian_process/tests/test_gpr.py b/sklearn/gaussian_process/tests/test_gpr.py
index b641be30a824a..8b3bc0dccd534 100644
--- a/sklearn/gaussian_process/tests/test_gpr.py
+++ b/sklearn/gaussian_process/tests/test_gpr.py
@@ -652,3 +652,31 @@ def test_gpr_predict_error():
     err_msg = "At most one of return_std or return_cov can be requested."
     with pytest.raises(RuntimeError, match=err_msg):
         gpr.predict(X, return_cov=True, return_std=True)
+
+
+def test_y_std_with_multitarget_normalized():
+    """Check the proper normalization of `y_std` in multi-target scene.
+
+    Non-regression test for:
+    https://github.com/scikit-learn/scikit-learn/issues/17394
+    https://github.com/scikit-learn/scikit-learn/issues/18065
+    """
+    rng = np.random.RandomState(1234)
+
+    n_samples, n_features, n_targets = 12, 10, 6
+
+    X_train = rng.randn(n_samples, n_features)
+    y_train = rng.randn(n_samples, n_targets)
+    X_test = rng.randn(n_samples, n_features)
+
+    # Generic kernel
+    kernel = WhiteKernel(1.0, (1e-1, 1e3)) * C(10.0, (1e-3, 1e3))
+
+    model = GaussianProcessRegressor(
+        kernel=kernel, n_restarts_optimizer=10, alpha=0.1, normalize_y=True
+    )
+    model.fit(X_train, y_train)
+    y_pred, y_std = model.predict(X_test, return_std=True)
+
+    assert y_pred.shape == (n_samples, n_targets)
+    assert y_std.shape == (n_samples, n_targets)


Code review comments:
diff_hunk: @@ -424,7 +424,11 @@ def predict(self, X, return_std=False, return_cov=False):
                     y_var[y_var_negative] = 0.0
 
                 # undo normalisation
-                y_var = y_var * self._y_train_std ** 2
+                y_var = np.outer(y_var, self._y_train_std ** 2)
+
+                # if y_var has shape (n_samples, 1), reshape to (n_samples,)
+                if y_var.shape[0] == y_var.size:
+                    y_var = np.ravel(y_var)
path: sklearn/gaussian_process/_gpr.py
text: @user1:
np.squeeze(y_var, axis=1) will do that automatically

@user1:
well not that easily since it will raise an error if the given axis has a size > 1 :/

@author:
I just saw the message now, but thanks for adding this feature.