sklearn gives different trees for scaled and unscaled data

I am not very convinced.  I used the same model with exact same hyperparameters and setting the same random state for a synthetically generated dataset for three scenarios. Firstly, made the tree with unscaled data, then using maxmin scaling and then scaling each input arbitrarily to close to double digits. The tree output was the same for the maxmin scaled and arbitrarily scaled data but different when it comes to the unscaled data. I am enclosing the code used for reference. This issue is not observed when I tried the same approach with other machine learning packages like MATLAB. I am also enclosing the zip file with the jupyter notebook and the output trees for reference.
[test_tree.zip](https://github.com/scikit-learn/scikit-learn/files/10844429/test_tree.zip)

```python
import pandas as pd
import numpy as np
from sklearn import tree
import dtreeviz

dfe  = np.random.exponential(scale=1e-7, size=(100000,1))
mld  = np.random.exponential(scale=1e2, size=(100000,1))
nh4  = np.random.exponential(scale=1e-5, size=(100000,1))
no3  = np.random.exponential(scale=1e-5, size=(100000,1))
po4  = np.random.exponential(scale=1e-5, size=(100000,1))
rsn  = np.random.exponential(scale=1e2, size=(100000,1))
si   = np.random.exponential(scale=1e-4, size=(100000,1))
sos  = np.random.exponential(scale=1e2, size=(100000,1))
tos  = np.random.exponential(scale=1e2, size=(100000,1))
wo5  = np.random.exponential(scale=1e-7, size=(100000,1))
phyc = np.random.exponential(scale=1e-4, size=(100000,1))

columns=['dfe', 'mld', 'nh4', 'no3', 'po4', 'rsn','si',
            'sos', 'tos', 'wo5','phyc']
arr = np.concatenate((dfe,mld,nh4,no3,po4,rsn,si,sos,tos,wo5,phyc), axis=1)
da = pd.DataFrame(arr, columns=columns)
df = da.copy(deep=True)

from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

X = df.drop('phyc', axis=1)
y = df['phyc']
X_train, X_test, y_train, y_test = train_test_split(X, y,test_size=0.1,random_state=101,shuffle=False)

import graphviz
from sklearn.tree import DecisionTreeRegressor
model = DecisionTreeRegressor(
    criterion='squared_error',
    splitter='best',
    min_samples_leaf=5,
    max_depth=5,
    # min_samples_split=2,
    random_state=101
    )

#Unscaled Tree

unscaled_tree=model.fit(X_train,y_train)
y_pred = unscaled_tree.predict(X_test)
rmse_unscaled = mean_squared_error(y_test,y_pred, squared=False)
print('rmse :',rmse_unscaled)
print(rmse_unscaled /df['phyc'].std())

from sklearn.tree import export_graphviz
estimator=unscaled_tree
export_graphviz(estimator,
                out_file='test_unscaled.dot',
                feature_names=X.columns,
                class_names='phyc',
                rounded=True,proportion=False,
                precision=2, filled=True,)

from subprocess import call
call(['dot','-Tpng', 'test_unscaled.dot', '-o', 'test_unscaled.png','-Gdpi=600'])


from IPython.display import Image
Image(filename = 'test_unscaled.png')

#Tree maxmin scaled

from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler(feature_range=(0,1))
dfs = scaler.fit_transform(da)
dfs = pd.DataFrame(dfs, columns=['dfe', 'mld', 'nh4', 'no3', 'po4', 'rsn','si','sos', 'tos', 'wo5','phyc'])
Xs = dfs.drop(['phyc'], axis=1)
ys = dfs['phyc']

Xs_train, Xs_test, ys_train, ys_test = train_test_split(Xs, ys, shuffle=False,test_size=0.1, random_state=101)

scaled_tree= model.fit(Xs_train,ys_train)
ys_pred = scaled_tree.predict(Xs_test)

from sklearn.tree import export_graphviz
estimator=scaled_tree
export_graphviz(estimator,
                out_file='test_scaled.dot',
                feature_names=X.columns,
                class_names='phyc',
                rounded=True,proportion=False,
                precision=2, filled=True,)
from subprocess import call
call(['dot','-Tpng', 'test_scaled.dot', '-o', 'test_scaled.png','-Gdpi=600'])
from IPython.display import Image
Image(filename = 'test_scaled.png')

#Tree arbitrary scaled

dfr = da.copy(deep=True)
dfr['dfe']=dfr['dfe']*1e9
dfr['nh4']=dfr['nh4']*1e7
dfr['no3']=dfr['no3']*1e7
dfr['phyc']=dfr['phyc']*1e5
dfr['po4']=dfr['po4']*1e6
dfr['si']=dfr['si']*1e5
dfr['wo5']=dfr['wo5']*1e9

Xr = dfr.drop('phyc', axis=1)
yr = dfr['phyc']
Xr_train, Xr_test, yr_train, yr_test = train_test_split(Xr, yr,test_size=0.1,random_state=101,shuffle=False)

arscaled_tree=model.fit(Xr_train,yr_train)
yr_pred = arscaled_tree.predict(Xr_test)

from sklearn.tree import export_graphviz
estimator=arscaled_tree

export_graphviz(estimator,
                out_file='test_arscaled.dot',
                feature_names=X.columns,
                class_names='phyc',
                rounded=True,proportion=False,
                precision=2, filled=True,)


from subprocess import call
call(['dot','-Tpng', 'test_arscaled.dot', '-o', 'test_arscaled.png','-Gdpi=600'])

from IPython.display import Image
Image(filename = 'test_arscaled.png')
```

_Originally posted by @sandupal in https://github.com/scikit-learn/scikit-learn/issues/25495#issuecomment-1447103075_

COMMENT:
We agree that it *can* be produced by rounding error. However, we have proved to our own satisfaction this is not due to rounding errors.

What we find is that a variable whose range is 1e-8 to 1e-5 is showing up as unimportant when it is used without scaling. However, it becomes the most important predictor when it is used with scaling- so that the variable does in fact explain such errors. Roundoff error would not produce such a change as the values, when scaled, would have essentially the same. It might produce changes in reproducibility, but should not fundamentally start finding different relationships. Again, we do not find this behavior using the same values when we use Matlab.  This tells us there is some problem in the implmentation of scikitlearn.

Labels: Needs Triage