Yeo-Johnson Power Transformer gives Numpy warning (and raises scipy.optimize._optimize.BracketError in some cases)

### Describe the bug

When I use a power transformer with yeo-johnson method I get this warning in numpy:

`../lib/python3.10/site-packages/numpy/core/_methods.py:235: RuntimeWarning: overflow encountered in multiply`

Strangely I can't surpress this warning with filter.warnings()

**Tl;dr: The eroor seems to appear when multiple instances of the transformer are called, in this example because n_jobs=2 is set in the preprocessor. But the bug also appears with n_jobs=1 when you call this preprocessor in a grid search for example.**

### Steps/Code to Reproduce

The following code snippet will raise a warning:

```python
import numpy as np
from sklearn.preprocessing import PowerTransformer

rng = np.random.default_rng(0)
X = rng.exponential(size=(50,))
X -= X.max()
X = -X
X += 100
X = X.reshape(-1, 1)

method = "yeo-johnson"
transformer = PowerTransformer(method=method, standardize=False)
transformer.fit_transform(X)
```

### Expected Results

No Warnings

### Actual Results
```pytb
/Users/glemaitre/mambaforge/envs/dev/lib/python3.8/site-packages/numpy/core/_methods.py:233: RuntimeWarning: overflow encountered in multiply
  x = um.multiply(x, x, out=x)

```

### Versions

```shell
System:
    python: 3.8.12 | packaged by conda-forge | (default, Sep 16 2021, 01:38:21)  [Clang 11.1.0 ]
executable: /Users/glemaitre/mambaforge/envs/dev/bin/python
   machine: macOS-12.3.1-arm64-arm-64bit

Python dependencies:
      sklearn: 1.2.dev0
          pip: 21.3
   setuptools: 58.2.0
        numpy: 1.21.6
        scipy: 1.8.0
       Cython: 0.29.24
       pandas: 1.4.2
   matplotlib: 3.4.3
       joblib: 1.0.1
threadpoolctl: 2.2.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /Users/glemaitre/mambaforge/envs/dev/lib/libopenblas_vortexp-r0.3.18.dylib
        version: 0.3.18
threading_layer: openmp
   architecture: VORTEX
    num_threads: 8

       user_api: openmp
   internal_api: openmp
         prefix: libomp
       filepath: /Users/glemaitre/mambaforge/envs/dev/lib/libomp.dylib
        version: None
    num_threads: 8
```

COMMENT:
I think there are two distinct problems to handle here:
1. Make the expressions that use `np.power` more robust against overflow.
2. Avoid overflow in `x_trans.var()` by regularizing the output's exponent. Without regularization, the output's exponent can blow up for certain inputs for very little gain in the loss function.

Here is a patch that addresses both of the above issues:
```python
import numpy as np
from scipy import optimize
from sklearn.preprocessing import PowerTransformer as OriginalPowerTransformer


class PowerTransformer(OriginalPowerTransformer):
    """Apply a power transform featurewise to make data more Gaussian-like."""

    def _yeo_johnson_inverse_transform(self, x, lmbda):
        """Return inverse-transformed input x following Yeo-Johnson inverse transform."""
        x_inv = np.zeros_like(x)
        pos = x >= 0

        # when x >= 0
        if abs(lmbda) < np.spacing(1.0):
            x_inv[pos] = np.exp(x[pos]) - 1
        else:  # lmbda != 0
            # x_inv[pos] = np.power(x[pos] * lmbda + 1, 1 / lmbda) - 1  # Original line
            x_inv[pos] = np.exp(np.log1p(x[pos] * lmbda) / lmbda) - 1

        # when x < 0
        if abs(lmbda - 2) > np.spacing(1.0):
            # x_inv[~pos] = 1 - np.power(-(2 - lmbda) * x[~pos] + 1, 1 / (2 - lmbda))  # Original line
            x_inv[~pos] = 1 - np.exp(np.log1p(-(2 - lmbda) * x[~pos]) / (2 - lmbda))
        else:  # lmbda == 2
            x_inv[~pos] = 1 - np.exp(-x[~pos])

        return x_inv

    def _yeo_johnson_transform(self, x, lmbda):
        """Return transformed input x following Yeo-Johnson transform."""
        out = np.zeros_like(x)
        pos = x >= 0  # binary mask

        # when x >= 0
        if abs(lmbda) < np.spacing(1.0):
            out[pos] = np.log1p(x[pos])
        else:  # lmbda != 0
            # out[pos] = (np.power(x[pos] + 1, lmbda) - 1) / lmbda  # Original line
            out[pos] = (np.exp(np.log1p(x[pos]) * lmbda) - 1) / lmbda

        # when x < 0
        if abs(lmbda - 2) > np.spacing(1.0):
            # out[~pos] = -(np.power(-x[~pos] + 1, 2 - lmbda) - 1) / (2 - lmbda)  # Original line
            out[~pos] = -(np.exp(np.log1p(-x[~pos]) * (2 - lmbda)) - 1) / (2 - lmbda)
        else:  # lmbda == 2
            out[~pos] = -np.log1p(-x[~pos])

        return out

    def _yeo_johnson_optimize(self, x):
        """Find and return optimal lambda parameter of the Yeo-Johnson transform."""
        x_tiny = np.finfo(np.float64).tiny

        def _neg_log_likelihood(lmbda):
            """Compute the negative log likelihood of the observed data x."""
            x_trans = self._yeo_johnson_transform(x, lmbda)
            x_trans_var = x_trans.var()

            # Reject transformed data that would raise a RuntimeWarning in np.log
            if x_trans_var < x_tiny:
                return np.inf

            log_var = np.log(x_trans_var)
            # loglike = -n_samples / 2 * log_var  # Original line
            # loglike += (lmbda - 1) * (np.sign(x) * np.log1p(np.abs(x))).sum()  # Original line
            loglike = -0.5 * log_var
            loglike += (lmbda - 1) * (np.sign(x) * np.log1p(np.abs(x))).mean()
            loglike -= 0.01 * np.mean(np.abs(np.log(np.abs(x_trans_var[x_trans_var != 0]))))  # Regularize the exponent

            return -loglike

        # the computation of lambda is influenced by NaNs so we need to
        # get rid of them
        x = x[~np.isnan(x)]
        # choosing bracket -2, 2 like for boxcox
        lmbda = cast(float, optimize.brent(_neg_log_likelihood, brack=(-2, 2)))
        return lmbda
```

Labels: Bug, module:preprocessing