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:
I gave this another shot and made a suggestion in https://github.com/scikit-learn/scikit-learn/pull/29681, based on @akaashp2000 's PR, but including tests.

However, I found another edge case that is not detectable (see https://github.com/scikit-learn/scikit-learn/pull/29681/files#r1719809543), to repro:

``` python
import sys
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

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())

# NOTE: Use just the last colum here
X_train = data[["Holders"]]
y_train = data["Claims"]

model_sklearn_lbfgs = PoissonRegressor(alpha=0, tol = 1e-8, max_iter = 10000, verbose = True).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', tol = 1e-8, max_iter = 1000).fit(X_train, y_train)
print(model_sklearn_nc.intercept_)
print(model_sklearn_nc.coef_)
```

which gives

```
RUNNING THE L-BFGS-B CODE

           * * *

Machine precision = 2.220D-16
 N =            2     M =           10
 This problem is unconstrained.

           * * *

Tit   = total number of iterations
Tnf   = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip  = number of BFGS updates skipped
Nact  = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F     = final function value

           * * *

   N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
    2      1      3      1     0     0   4.300D+04  -1.426D+02

CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH             
3.896592058397603
[0.]
3.2885810888136406
[0.00087977]
```

How to proceed?

Labels: Bug, help wanted