Intercepts of Newton-Cholesky logistic regression get corrupted when warm starting

### Describe the bug

When using multinomial logistic regression with warm starts from a previous iteration, the final coefficients in the model are correct, but the intercepts somehow get filled with incorrect numbers somewhere.

As a result, predictions from a warm-started model differ from those of a cold-start model that has more iterations on the same data.

The issue appears to have been introduced recently as it works fine with version 1.5, but not with 1.6 or 1.7.

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
X, y = load_iris(return_X_y=True)

model1 = LogisticRegression(
    solver="newton-cholesky",
    max_iter=2
).fit(X, y)
model2 = LogisticRegression(
    solver="newton-cholesky",
    max_iter=1,
    warm_start=True
).fit(X, y).fit(X, y)

np.testing.assert_almost_equal(
    model1.coef_,
    model2.coef_
)

np.testing.assert_almost_equal(
    model1.predict_proba(X[:5]),
    model2.predict_proba(X[:5])
)
```

### Expected Results

Intercepts should be the same, up to shifting by a constant if needed.

### Actual Results

Intercepts are different, as are predicted probabilities

### Versions

```shell
System:
    python: 3.12.6 | packaged by conda-forge | (main, Sep 22 2024, 14:16:49) [GCC 13.3.0]
executable: /home/david/miniforge3/bin/python
   machine: Linux-6.12.33+deb12-amd64-x86_64-with-glibc2.36

Python dependencies:
      sklearn: 1.7.1
          pip: 24.2
   setuptools: 74.1.2
        numpy: 2.0.1
        scipy: 1.14.1
       Cython: 3.1.0
       pandas: 2.2.3
   matplotlib: 3.9.2
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 20
         prefix: libscipy_openblas
       filepath: /home/david/.local/lib/python3.12/site-packages/numpy.libs/libscipy_openblas64_-99b71e71.so
        version: 0.3.27
threading_layer: pthreads
   architecture: Haswell

       user_api: blas
   internal_api: mkl
    num_threads: 14
         prefix: libmkl_rt
       filepath: /home/david/miniforge3/lib/libmkl_rt.so.2
        version: 2023.2-Product
threading_layer: gnu

       user_api: openmp
   internal_api: openmp
    num_threads: 20
         prefix: libgomp
       filepath: /home/david/miniforge3/lib/libgomp.so.1.0.0
        version: None
```

COMMENT:
> Yes, this step is done after the optimization procedure (using `n_iter_` steps) to make all unpenalized parameters unique, i.e., sum to zero.
> 
> As I explained, this destroys the equivalence of same number of steps with warm start vs all steps at once. It is even worse: `LogisticRegression(max_iter=6, solver='newton-cholesky')` converges, but calling 6 or more time `model2.fit(X, y)` does not converge.
> 
> So what is more important:
> 
> * Making repeated warm-start calls to `fit` equivalent to one `fit` with many steps.
> * Keeping current symmetrizing behavior.
> 
> The newton-cholesky solver is not a SGD solver nor a minibatch solver and has no `partial_fit` method. You should give it at least some number of iterations and if you prefer an even higher accuracy later, you can use a warm start fit to improve.
> 
> For me, this behavior is fine because I don't have a use case to call repeated fit with small number of iterations.

I don't follow.

If it were the case that the _only_ thing going on is the centering of intercepts, then the predicted probabilities would be the same regardless of whether the intercepts sum to zero or to some other number, so long as all of them are off by the same constant.

For example:
```python
import numpy as np
from scipy.special import softmax
rng = np.random.default_rng(seed=123)
k = 3
n = 10
m = 5
X = rng.standard_normal(size=(m,n))
random_coefs = rng.standard_normal(size=(n, k))
random_intercepts = rng.standard_normal(size=(1,k))
centered_intercepts = random_intercepts - random_intercepts.sum()

pred_uncentered = X@random_coefs + random_intercepts
pred_centered = X@random_coefs + centered_intercepts

np.testing.assert_almost_equal(
    softmax(pred_uncentered),
    softmax(pred_centered)
)
```

The _only_ way in which an exact Newton update would give the same result for coefficients given the same start, is if predictions from which gradients and Hessians were the same, which is the case whenever the intercepts are equivalent up to an arbitrary constant shift - i.e.:
```python
np.testing.assert_almost_equal(
    random_intercepts - random_intercepts.sum(),
    centered_intercepts - centered_intercepts.sum()
)
```

Hence, centering of the intercepts should not change the end result in terms of predicted probabilities.

Yet in this case, the coefficients are matching, which means that the intercepts are being used correctly prior to the Newton update, but are then being incorrectly modified in a non-equivalent way afterwards.

If they received the same update and get centered after the update, the result would still be equivalent between the warm and cold start, at the very least in terms of predicted probabilities.

Labels: Bug, module:linear_model