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

- yes I am in fact getting the same output with ``` verbose = 100 ``` as you
- ... and the same issue of infinite values for ``` grad_pointwise ```

following your suggestion to print the value at each iteration, I added a line ```print(gradient_pointwise)``` at [Line 290 of _linear_loss.py](https://github.com/scikit-learn/scikit-learn/blob/d90dd114ab85b85e49facdb7b37242b7362696c8/sklearn/linear_model/_linear_loss.py#L290). Without any ```verbose``` value specified, I found that the first value printed had reasonable coordinate values, then the second value had very large/infinite coordinate values, then it reverted to the first value of ```grad_pointwise```, and then printed the intercept (3.896...) and coefficients (all zeroes)

@lorentzenchr 

I have now tried fitting the GLM with LBFGS, now on the same X_train but with the "Holders" column scaled, i.e.
```python
X_train_scaled = X_train.copy()
X_train_scaled["Holders"] = X_train_scaled["Holders"]/X_train_scaled["Holders"].max()
```
and it does indeed converge now! And to the same coefficients as from ```solver='newton-cholesky'``` (aside from the scaling of 'Holders' coefficient due to the scaling of the feature). Scaling the feature works for both statsmodels and scikit-learn now.

Labels: Bug, help wanted