Thanks for the review @adrinjalali. I updated this PR with the following changes:

1. I changed x_trans_exp = np.log(np.abs(x_trans)) to x_trans_exp = np.log(np.abs(x_trans[x_trans != 0])) to avoid RuntimeWarning: divide by zero encountered in log.
2. I added a changelog entry as requested.
3. I rebased the PR on the latest upstream main.

Hmm, this is changing some values in a docstring we have (hence fail in CI)

@lsorber please avoid force pushing so that we can see the history of changes. At the end we squash and merge anyway.

@adrinjalali Apologies, I thought I had rebased on the latest upstream main, but found that wasn't the case and so I did another rebase and force push. I'll avoid doing that from now.

It seems that the doctests fail because they expect to see specific lambda values for the given example. However, because of the new regularisation term on the exponent of the transformed variables, the lamdbas are slightly different from the expected lambdas (see below).

How should we proceed from here?

3076 
3077     Examples
3078     --------
3079     >>> import numpy as np
3080     >>> from sklearn.preprocessing import PowerTransformer
3081     >>> pt = PowerTransformer()
3082     >>> data = [[1, 2], [3, 2], [4, 5]]
3083     >>> print(pt.fit(data))
3084     PowerTransformer()
3085     >>> print(pt.lambdas_)
Expected:
    [ 1.386... -3.100...]
Got:
    [ 1.3761854 -3.091368 ]

/home/vsts/work/1/s/sklearn/preprocessing/_data.py:3085: DocTestFailure

maybe @lorentzenchr would have a better idea how to proceed here.

Can we use scipy.stats.yeojohnson (requiring scipy 1.2 should fine) and contribute such fixes upstream in scipy?

I reviewed scipy's implementation and it's very similar to sklearn's implementation in that it is prone to both of the sources of overflow that this PR addresses.

I am open to applying the improvements of this PR to scipy and then updating this PR to use scipy's implementation. Are you asking me to do this @lorentzenchr? Should we check with the scipy maintainers whether they are open to such a change? The new regularisation term is non-standard and so may meet some resistance.

I think the idea is to not have our own implementation, and use scipy's when available, and suggest this fix on the scipy side.

I opened #26308 to rely on scipy.stats.yeojohnson.
I think the improvement of this PR should be contributed upstream in scipy. Therefore I close.