PoissonRegressor lbfgs solver giving coefficients of 0 and Runtime Warning

### Describe the bug

See the following [stack exchange](https://stats.stackexchange.com/questions/622085/sklearn-poissonregressor-giving-all-coefficients-zero) post (the solution to my original issue was to use newton-cholesky solver)

When fitting a Poisson Regression (without regularization) to some dummy data I encounter:

- with lbfgs solver, a Runtime Warning, a non-zero intercept and all coefficients as zero
- with newton-cholesky solver, coefficients as expected

Some people on StackExchange have mentioned it is worth submitting an issue (there was a similar [one](https://github.com/scikit-learn/scikit-learn/issues/24752) faced with Logistic Regression).

### Steps/Code to Reproduce

```python
import statsmodels.api as sm
import statsmodels.formula.api as smf
import pandas as pd
from sklearn.linear_model import PoissonRegressor
from sklearn.preprocessing import (
    OneHotEncoder,
)
data = sm.datasets.get_rdataset('Insurance', package='MASS').data
# Fit Poisson regression using formula interface
formula = "Claims ~ C(District, Treatment(1)) + C(Group, Treatment('<1l')) + C(Age, Treatment('<25')) + Holders"
model_smf = smf.poisson(formula=formula, data=data).fit()
print(type(model_smf))
print(model_smf.summary())
# with sklearn OneHotEncoder

X_train_ohe = OneHotEncoder(sparse_output=False, drop=[1, "<1l", "<25"]).fit(data[["District", "Group", "Age"]])
X_train_ohe = pd.DataFrame(X_train_ohe.transform(data[["District", "Group", "Age"]]), columns=X_train_ohe.get_feature_names_out())

X_train = pd.concat([X_train_ohe, data[["Holders"]]], axis=1)
y_train = data["Claims"]

# one-hot encode the categorical columns, and drop the baseline column
# with lbfgs solver

model_sklearn_lbfgs = PoissonRegressor(alpha=0).fit(X_train, y_train)

print(model_sklearn_lbfgs.intercept_)
print(model_sklearn_lbfgs.coef_)
# with newton-cholesky solver

model_sklearn_nc = PoissonRegressor(alpha=0, solver='newton-cholesky').fit(X_train, y_train)

print(model_sklearn_nc.intercept_)
print(model_sklearn_nc.coef_)
```

### Expected Results

Expected intercept/coefficients (from statsmodels and sklearn with newton-cholesky solver):
```
2.8456859090224444
[-4.47436243e-01 -9.30629212e-01 -1.46547942e+00  1.00603354e+00
  4.71328611e-01 -5.98213625e-01  5.69685530e-01  6.85301286e-01
  2.22504699e+00 -1.58161172e-05]
```

### Actual Results

Result with lbfgs:
```
3.896592058397603
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
/home/akaash/Downloads/statsmodels-play/venv/lib/python3.10/site-packages/sklearn/linear_model/_linear_loss.py:290: RuntimeWarning: invalid value encountered in matmul
  grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights
```

### Versions

```shell
System:
    python: 3.10.9 (main, Mar  1 2023, 18:23:06) [GCC 11.2.0]
executable: /home/akaash/Downloads/statsmodels-play/venv/bin/python
   machine: Linux-5.19.0-46-generic-x86_64-with-glibc2.35

Python dependencies:
      sklearn: 1.3.0
          pip: 23.2.1
   setuptools: 68.0.0
        numpy: 1.25.1
        scipy: 1.11.1
       Cython: None
       pandas: 2.0.3
   matplotlib: 3.7.2
       joblib: 1.3.1
threadpoolctl: 3.2.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 6
         prefix: libopenblas
       filepath: /home/akaash/Downloads/statsmodels-play/venv/lib/python3.10/site-packages/numpy.libs/libopenblas64_p-r0-7a851222.3.23.so
        version: 0.3.23
threading_layer: pthreads
   architecture: Haswell

       user_api: blas
   internal_api: openblas
    num_threads: 6
         prefix: libopenblas
       filepath: /home/akaash/Downloads/statsmodels-play/venv/lib/python3.10/site-packages/scipy.libs/libopenblasp-r0-23e5df77.3.21.dev.so
        version: 0.3.21.dev
threading_layer: pthreads
   architecture: Haswell

       user_api: openmp
   internal_api: openmp
    num_threads: 6
         prefix: libgomp
       filepath: /home/akaash/Downloads/statsmodels-play/venv/lib/python3.10/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
        version: None
```

COMMENT:
@ogrisel 

I am trying the approach you have suggested above - from what I understand using ```np.errstate``` means that instead of the RuntimeWarning, a (```FloatingPointError```) exception is raised, so I have tried the following change: where the matmul was before, I now have:

```python
with np.errstate(all="raise"):
    try:
        grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights
    except FloatingPointError:
        raise ValueError(
            "Overflow detected. Try scaling the target variable or"
            " features, or using a different solver"
        ) from None
```

(I believe the ```None``` surpresses the original FloatingPointErrror, so only the ValueError is returned... below is the output)

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], line 4
      1 # one-hot encode the categorical columns, and drop the baseline column
      2 # with lbfgs solver
----> 4 model_sklearn_lbfgs = PoissonRegressor(alpha=0).fit(X_train, y_train)
      6 print(model_sklearn_lbfgs.intercept_)
      7 print(model_sklearn_lbfgs.coef_)

File ~/Downloads/scikit-learn/sklearn/base.py:1215, in _fit_context..decorator..wrapper(estimator, *args, **kwargs)
   1208     estimator._validate_params()
   1210 with config_context(
   1211     skip_parameter_validation=(
   1212         prefer_skip_nested_validation or global_skip_validation
   1213     )
   1214 ):
-> 1215     return fit_method(estimator, *args, **kwargs)

File ~/Downloads/scikit-learn/sklearn/linear_model/_glm/glm.py:265, in _GeneralizedLinearRegressor.fit(self, X, y, sample_weight)
    262 if self.solver == "lbfgs":
    263     func = linear_loss.loss_gradient
--> 265     opt_res = scipy.optimize.minimize(
    266         func,
    267         coef,
    268         method="L-BFGS-B",
...
    297         ) from None
    298 if self.fit_intercept:
    299     grad[-1] = grad_pointwise.sum()

ValueError: Overflow detected. Try scaling the target variable or features, or using a different solver
```

Is this what you meant? Could something similar could be done for multiclass case too (i.e. when ```self.base_loss.is_multiclass is True```)?

Labels: Bug, help wanted