Inconsistent Results For Logistic Regressions across multiple computers

### Describe the bug

Hey all - 

I'm working on teaching some students some logistic regressions and noticed different computers can produce slightly different intercepts/coefs. At first I thought it was maybe environment differences, but I have been able to reproduce the variances when accounting for the various packages/interpreters. 

I can create a new environment on computer 1 (windows 10, intel 8th gen CPU), computer 2 (windows 10, intel 11th gen CPU), and a coworker's M1 MacBook - all three produce different results. I thought maybe minute differences in Numpy's Openblas or MKL could be the culprit but that yielded the same varying results. I've tried not splitting the data, different random states, different C values, different solvers... 

For the purposes of reproducibility, I ran this is exact code & uploaded the data to Google sheets:

```python
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

bank_df = pd.read_csv('https://docs.google.com/spreadsheets/d/e/2PACX-1vSGfjG4mq1_4HS4iwRN7EZK6YHzPDi8HpB_giY7kiqbDZRsRNjbfhuQ2J6xkHGk1YVYN9H0TxOf2tgw/pub?gid=1909291157&single=true&output=csv')
bank_df.drop(columns=['ID', 'ZIP_Code'], inplace=True)

bank_df['Education'] = bank_df['Education'].astype('category')
new_categories = {1: 'Undergrad', 2: 'Graduate', 3: 'Advanced/Professional'}
bank_df.Education.cat.rename_categories(new_categories, inplace=True)
bank_df = pd.get_dummies(bank_df, prefix_sep='_', drop_first=True)

y = bank_df['Personal_Loan']
X = bank_df.drop(columns=['Personal_Loan'])

# partition data
train_X, valid_X, train_y, valid_y = train_test_split(X, y, test_size=0.4, random_state=1)

# fit a logistic regression (set penalty=l2 and C=1e42 to avoid regularization)
logit_reg = LogisticRegression(penalty="l2", C=1e42, solver='liblinear', random_state=0)
logit_reg.fit(train_X, train_y)
```

All versions of all other packages are held constant. I received this results across the three computers: 


| values | Computer 1 with Sklearn 1.1.2 | Computer 2 with Sklearn 1.1.2 | M1 Macbook with Sklearn 1.1.1 |
|---|---|---|---|
|Intercept|	-12.4919578542609|	-12.6123809362746	| -12.4934360611768|
|Age coef|	-0.037738|	-0.032577	|-0.037685|
|Experience coef|	0.039255|	0.034171|0.039202|
|Income coef|	0.058843|	0.058837	|0.058844|
|Family coef|	0.612243|	0.613192	|0.612251|

I have tried many other things that all produce slightly different results that are not captured in the above table. Multicollinearity is present within the data, but I would still expect consistent results.

I understand the results are effectively the same and this produces little real world impacts. However, when a student is comparing results and they see 12.4 vs 12.6 which is a decently large difference. I think some of these differences are too large to just be contributed to floating-point computations -- especially for size of the data. I assume I'm missing something, but I'm reaching the end of my rope on troubleshooting. Thoughts / expected behavior? 

### Steps/Code to Reproduce

```python
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression

bank_df = pd.read_csv('https://docs.google.com/spreadsheets/d/e/2PACX-1vSGfjG4mq1_4HS4iwRN7EZK6YHzPDi8HpB_giY7kiqbDZRsRNjbfhuQ2J6xkHGk1YVYN9H0TxOf2tgw/pub?gid=1909291157&single=true&output=csv')
bank_df.drop(columns=['ID', 'ZIP_Code'], inplace=True)

bank_df['Education'] = bank_df['Education'].astype('category')
new_categories = {1: 'Undergrad', 2: 'Graduate', 3: 'Advanced/Professional'}
bank_df.Education.cat.rename_categories(new_categories, inplace=True)
bank_df = pd.get_dummies(bank_df, prefix_sep='_', drop_first=True)

y = bank_df['Personal_Loan']
X = bank_df.drop(columns=['Personal_Loan'])

# partition data
train_X, valid_X, train_y, valid_y = train_test_split(X, y, test_size=0.4, random_state=1)

# fit a logistic regression (set penalty=l2 and C=1e42 to avoid regularization)
logit_reg = LogisticRegression(penalty="l2", C=1e42, solver='liblinear')
logit_reg.fit(train_X, train_y)
```

### Expected Results

| values | Computer 1 with Sklearn 1.1.2 | Computer 2 with Sklearn 1.1.2 | M1 Macbook with Sklearn 1.1.1 |
|---|---|---|---|
|Intercept|	-12.4919578542609|	-12.6123809362746	| -12.4934360611768|
|Age coef|	-0.037738|	-0.032577	|-0.037685|
|Experience coef|	0.039255|	0.034171|0.039202|
|Income coef|	0.058843|	0.058837	|0.058844|
|Family coef|	0.612243|	0.613192	|0.612251|

### Actual Results

n/a

### Versions

```shell
System:
    python: 3.9.13 (main, Aug 25 2022, 23:51:50) [MSC v.1916 64 bit (AMD64)]
executable: C:\Users\Nick\.conda\envs\sktest\python.exe
   machine: Windows-10-10.0.19044-SP0

Python dependencies:
      sklearn: 1.1.2
          pip: 22.2.2
   setuptools: 63.4.1
        numpy: 1.23.1
        scipy: 1.9.1
       Cython: None
       pandas: 1.4.4
   matplotlib: None
       joblib: 1.1.0
threadpoolctl: 2.2.0

Built with OpenMP: True

threadpoolctl info:
       filepath: C:\Users\Nick\.conda\envs\sktest\Lib\site-packages\sklearn\.libs\vcomp140.dll
         prefix: vcomp
       user_api: openmp
   internal_api: openmp
        version: None
    num_threads: 12

       filepath: C:\Users\Nick\.conda\envs\sktest\Library\bin\mkl_rt.1.dll
         prefix: mkl_rt
       user_api: blas
   internal_api: mkl
        version: 2021.4-Product
    num_threads: 6
threading_layer: intel
```

COMMENT:
> What happens with C=np.inf instead of C=1e42?

It doesn't converge :)

With LBFGS, since you have the verbose, you can see small numerical differences (in term of objective and gradient) at each iteration that leads the coefficients to diverge after a while.

Labels: Bug, module:linear_model, Needs Investigation