The data type of input data for LinearRegression class will affect the results

### Describe the bug

Our team just used the class sklearn.linear_model.LinearRegression to do multi-linear regression. And, we found out that the same data, which means the values are identical for each element, with different data formats will return different coefficients. The array we use are all integer values but saved as float32 and float64, separately. However, the array in float64 as input get us the expected result, but float32 doesn't.

The code below is a simple reproducible example.

### Steps/Code to Reproduce

```python
In [1]: from sklearn.linear_model import LinearRegression
    ...: import numpy as np
    ...: rng = np.random.default_rng(0)
    ...: x_1_array = rng.integers(low=1000, high=10000, size=10_000_000)
    ...: x_2_array = rng.integers(low=1000, high=10000, size=10_000_000)
    ...: x_3_array = rng.integers(low=1000, high=10000, size=10_000_000)
    ...: y_array = 3 * x_1_array + x_2_array - 5 * x_3_array
    ...: intercept_data = rng.integers(low=-3, high=3)
    ...: print(f'True intercept {intercept_data}')
    ...: y_array += intercept_data
    ...: X_train = np.stack([x_1_array,x_2_array,x_3_array],axis=-1)
    ...: y_train = y_array
    ...: linreg = LinearRegression()
    ...: model = linreg.fit(X_train.astype(np.float64), y_train)
    ...: print(f'Float 64 Intercept: {linreg.intercept_}')
    ...: print(f'Float 64 Coefficient: {linreg.coef_}')
    ...: model = linreg.fit(X_train.astype(np.float32), y_train)
    ...: print(f'Float 32 Intercept: {linreg.intercept_}')
    ...: print(f'Float 32 Coefficient: {linreg.coef_}')
```

### Expected Results

![image](https://user-images.githubusercontent.com/51532833/168492124-d98f9663-4899-47ee-83b2-4952462c37a6.png)

### Actual Results

![image](https://user-images.githubusercontent.com/51532833/168492136-4e0f4100-9e48-4f58-aa2a-97e826f06c61.png)


### Versions

```shell
System:
    python: 3.7.9 (v3.7.9:13c94747c7, Aug 15 2020, 01:31:08)  [Clang 6.0 (clang-600.0.57)]
executable: /Users/lazy/PycharmProjects/projects/bin/python
   machine: Darwin-21.4.0-x86_64-i386-64bit

Python dependencies:
          pip: 21.3
   setuptools: 58.2.0
      sklearn: 1.0.2
        numpy: 1.21.4
        scipy: 1.7.2
       Cython: None
       pandas: 1.3.4
   matplotlib: 3.4.3
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True
```

COMMENT:
Hello everyone ! I worked on this bug on my side and here are some feedbacks in case this issue is still relevant :

1. Something I found strange is that when I did float64 --> float32, the error was bigger than float64 --> float16. I don't understand why since the rounding should be more important, hence more error ?

2. I tested on the Ridge and I had the same behaviour indeed but with Lasso I had very good results :

Float 64 Intercept: 20.000815382406472
Float 64 Coefficient: [ 2.99999985  0.99999985 -4.99999985]

Float 32 Intercept: 19.98779296875
Float 32 Coefficient: [ 3.         1.0000001 -5.       ]

I believe that it is due to the use of the following line for LinearRegression and Ridge : coef_, _, rank_, singular_ = linalg.lstsq(X_train, y_train) whereas Lasso is more sophisticated, hence better results. I guess it might be linked to a more precise calculation of the coefficients thanks to the l1 norm.

3. Regarding the influence of the number of data I found out that the value of the maximum is a huge factor of error. Its influence is bigger than the size of the possible value range in the data :

X_train from 10000 to 10100 : 502377 % (relative error)
X_train from 4900 to 5000 : 247330 % (relative error)

Bigger numbers have a big impact on the intercept calculation and I believe that it is one of the main issue. We have important roundings causing the X_offset value to be very high when using float32. And here comes the huge factor of error propagating by the dot product operation y_offset-np.dot(X_offset, self.coef_.T).

4. Using a Standard Scaler didn't make the calculation of the coefficients that bad on my side compared to what was presented. Here are my results :

Float 64 Intercept: -5477.377062199998
Float 64 Coefficient: [  7793.46267406   2597.7576475  -12992.64703011]

Float 32 Intercept: -5477.37744140625
Float 32 Coefficient: [  7793.746    2597.8232 -12985.981 ]

Except for the last coefficient we almost have a perfect. I used the given data so I don't understand why Mr. glemaitre has such a big difference compare to me. I precise that I used the version of sklearn mentionned in the issue.

To me it seems that it is a pretty good solution, but not perfect since we still have a small variation and we force the normalization.

5. Since the calculation of the intercept depends on the values of the coefficients it may be interesting to consider this calculation as an another source of error for the intercept value. Maybe we should go back to the linalg.lstsq() function in order to make it more robust to type-casting and normalization since we saw a small error ? Mr. glemaitre mentionned this too as a potential error source.

I hope my humble contribution might interest/help you a bit.

Thank you in advance for reading and giving me some feedback.

Labels: Bug, module:linear_model, float32