Lasso loss is (mathematically) invariant wrt simultaneous rotation of X and y, but the solver outcome is not

### Describe the bug

For any given set of coefficients and intercept, the Lasso loss is invariant with respect to replacing X and y with U@X and U@y, where U is any orthonormal matrix (proof see below). However in practice, when I run Lasso with fit_intercept=True, that is not the case, not just the coefficient values, but the set of nonzero coefficients is different after applying such a rotation if the rotation is sufficiently far from the unit matrix. 

Setting tol to very low values doesn't seem to help.

Here the invariance proof:

y: n x 1
X: n x m
a: m x 1

err(a) = y - X @ a
loss(a) = err(a).T @ err(a) + c * sum(abs(a))

If U is any orthonormal matrix, that is, U.T @ U == np.eye(len(U))
and we define

Xnew := U @ X
ynew := U @ y

Then 
err_new(a) = ynew - Xnew @ a = U @ y - U @ X @ a = U @ err(a)
err_new(a).T @ err_new(a) = err(a).T @ U.T @ U @ err(a) = err(a).T @ err(a)
thus 
loss_new(a) = loss(a) for any a

So the Lasso solver looking for an optimal a should return the same coefficients for the rotated and unrotated versions.

### Steps/Code to Reproduce

```
import pickle
import numpy as np
from sklearn.linear_model import Lasso
from scipy.linalg import expm

N= 1000
M = 300
np.random.seed(42)
X = np.random.randn(N, M)
y = np.random.randn(N)

def make_U(size: int, eps: float):
    """
    Make an orthonormal matrix
    @param eps : controls closeness to unit matrix
    """
    randmat = np.random.randn(size, size)
    U = expm(eps * (randmat.T - randmat))
    maxerr = np.max(np.abs(U.T.dot(U) - np.eye(len(U))))
    assert maxerr < 1e-10
    print("Biggest U deviation from eye:", np.max(np.abs(U - np.eye(len(U)))))
    print("Biggest U.T@U deviation from eye:", maxerr)
    return U

def loss(lasso, X, y, alpha):
    err = lasso.predict(X) - y
    loss = (1/(2*len(y))) * (err*err).sum() + np.abs(lasso.coef_).sum()*alpha
    return loss

alpha = 0.0447

lasso_args = {
    "max_iter": int(1e5),
    "fit_intercept": True,
    "tol": 1e-12,
}

for eps in [0.0, 0.1]:
    print("eps:", eps)
    U = make_U(len(y), eps)
    Xnew = U.dot(X)
    ynew = U.dot(y)
    lasso = Lasso(alpha=alpha, **lasso_args)
    lasso.fit(Xnew, ynew)

    nonzeros = np.nonzero(lasso.coef_)[0]
    print(nonzeros)
    print("loss:", loss(lasso, X, y, alpha))
```

### Expected Results

The loss value and the set of nonzero coefficients is the same for eps=0.0 and eps=0.1


### Actual Results

Output running against the random data above:

```
eps: 0.0
Biggest U deviation from eye: 0.0
Biggest U.T@U deviation from eye: 0.0
[ 15  17  21  24  28  53  56  69  79  81  86  93 111 125 142 156 179 184
 185 194 195 198 204 211 221 224 233 236 254 261 263 270 273 283 285 286
 289 299]
loss: 0.44069201663480084
eps: 0.1
Biggest U deviation from eye: 1.0497406241802238
Biggest U.T@U deviation from eye: 6.661338147750939e-15
[ 15  17  19  21  24  28  53  56  69  79  81  86  93 111 125 142 156 179
 184 185 188 194 195 198 204 211 221 224 233 236 254 261 263 270 271 273
 283 285 286 289 299]
loss: 0.4433901503481833
```

Output running against my data (happy to provide CSVs), the difference in set of nonzero coefficients is much bigger:
```
eps: 0.0
Biggest U deviation from eye: 0.0
Biggest U.T@U deviation from eye: 0.0
[ 12  24  26 119 126 236]
loss: 20.46512689044932
eps: 0.1
Biggest U deviation from eye: 1.0662997612721612
Biggest U.T@U deviation from eye: 7.993605777301127e-15
[  0   1   2   3   4   5   6   7   9  10  11  12  13  14  15  16  17  18
  19  20  21  22  23  24  26 119 126 245]
loss: 20.844058180320488
```

### Versions

```shell
System:
    python: 3.8.16 (default, Mar  2 2023, 03:18:16) [MSC v.1916 64 bit (AMD64)]
executable: C:\...\.conda\envs\pizza3.8\python.exe
   machine: Windows-10-10.0.19044-SP0

Python dependencies:
      sklearn: 1.2.1
          pip: 23.0.1
   setuptools: 65.6.3
        numpy: 1.21.6
        scipy: 1.9.0
       Cython: None
       pandas: 1.3.5
   matplotlib: 3.5.3
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Users\...\.conda\envs\pizza3.8\Lib\site-packages\numpy\.libs\libopenblas.XWYDX2IKJW2NMTWSFYNGFUWKQU3LYTCZ.gfortran-win_amd64.dll
        version: 0.3.17
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12

       user_api: openmp
   internal_api: openmp
         prefix: vcomp
       filepath: C:\Users\...\.conda\envs\pizza3.8\Lib\site-packages\sklearn\.libs\vcomp140.dll
        version: None
    num_threads: 12

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Users\...\.conda\envs\pizza3.8\Lib\site-packages\scipy\.libs\libopenblas.PZA5WNOTOH6FZLB2KBVKAURAKVTFSNNU.gfortran-win_amd64.dll
        version: 0.3.18
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12

System:
    python: 3.8.16 (default, Mar  2 2023, 03:18:16) [MSC v.1916 64 bit (AMD64)]
executable: C:\Users\...\.conda\envs\pizza3.8\python.exe
   machine: Windows-10-10.0.19044-SP0

Python dependencies:
      sklearn: 1.2.1
          pip: 23.0.1
   setuptools: 65.6.3
        numpy: 1.21.6
        scipy: 1.9.0
       Cython: None
       pandas: 1.3.5
   matplotlib: 3.5.3
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Users\...\.conda\envs\pizza3.8\Lib\site-packages\numpy\.libs\libopenblas.XWYDX2IKJW2NMTWSFYNGFUWKQU3LYTCZ.gfortran-win_amd64.dll
        version: 0.3.17
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12

       user_api: openmp
   internal_api: openmp
         prefix: vcomp
       filepath: C:\Users\...\.conda\envs\pizza3.8\Lib\site-packages\sklearn\.libs\vcomp140.dll
        version: None
    num_threads: 12

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Users\...\.conda\envs\pizza3.8\Lib\site-packages\scipy\.libs\libopenblas.PZA5WNOTOH6FZLB2KBVKAURAKVTFSNNU.gfortran-win_amd64.dll
        version: 0.3.18
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12
```

COMMENT:
Why is it not? As per the math above, I'm multiplying X and y with U _from the left_ while keeping the coefficient vector constant. As the coefficient vector doesn't change, neither does the L1-penalty and thus neither does the loss. 

What you must be thinking of is multiplying X with an orthonormal matrix from the right and leaving y untouched, in which case indeed the loss and the optimal solution would change.

Also, when I run the solver with `fit_intercept=False` the expected invariance does hold in every case I tried.

Labels: Bug, Needs Triage