cross_entropy_loss value and gradient change between calls with 4+ threads, so multinomial LogisticRegression stops early #3820

Describe the bug

optimization_solver_cross_entropy_loss returns different values and gradients for the same X, y and beta when it runs on 4 or more threads. With 1 or 2 threads, repeated calls are bit-identical and match a NumPy reference.

This breaks sklearnex's multinomial LogisticRegression(solver="lbfgs"). scipy's L-BFGS-B sees an objective that changes from call to call and stops on CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH. That is status 0, so there is no ConvergenceWarning. The model silently ends up much worse, and the fit looks 25x to 50x faster than it should.

The binary logistic_loss is not affected: it is bit-identical at every thread count I tried.

To Reproduce

import numpy as np
import daal4py

rng = np.random.default_rng(0)
n_samples, n_features, n_classes = 200_000, 100, 7
X = rng.standard_normal((n_samples, n_features))
y = rng.integers(0, n_classes, (n_samples, 1)).astype(np.float64)
beta = 0.1 * rng.standard_normal((n_classes * (n_features + 1), 1))

loss = daal4py.optimization_solver_cross_entropy_loss(
    nClasses=n_classes,
    numberOfTerms=n_samples,
    interceptFlag=True,
    resultsToCompute="value|gradient",
)
loss.setup(X, y, beta)


def evaluate(n_threads):
    daal4py.daalinit(n_threads)
    res = loss.compute(X, y, beta)
    return res.valueIdx[0, 0], res.gradientIdx.ravel().copy()


value_ref, grad_ref = evaluate(1)
for n_threads in [1, 2, 4, 8, 14]:
    # Same X, y and beta every time: results should only differ by rounding.
    results = [evaluate(n_threads) for _ in range(20)]
    values = np.array([v for v, _ in results])
    grads = np.array([g for _, g in results])
    print(
        f"{n_threads:2d} threads: "
        f"max rel. diff to 1 thread: value {np.abs(values - value_ref).max() / abs(value_ref):.1e}, "
        f"gradient {np.abs(grads - grad_ref).max() / np.abs(grad_ref).max():.1e}"
    )

Output on an Intel Core Ultra 5 225H:

 1 threads: max rel. diff to 1 thread: value 0.0e+00, gradient 0.0e+00
 2 threads: max rel. diff to 1 thread: value 0.0e+00, gradient 0.0e+00
 4 threads: max rel. diff to 1 thread: value 1.5e-05, gradient 3.6e-03
 8 threads: max rel. diff to 1 thread: value 1.1e-04, gradient 4.9e-03
14 threads: max rel. diff to 1 thread: value 1.1e-04, gradient 7.2e-03

The single-threaded value matches -log_softmax(X @ W.T + b)[range(n), y].mean() computed with NumPy/SciPy. Differences of 1e-4 on the value and 1e-2 on the gradient are far above what reordering a float64 sum would give (~1e-15), so I think this is a race rather than rounding.

End to end with sklearnex:

import numpy as np
from sklearn.datasets import make_classification
from sklearn.metrics import log_loss
from sklearnex.linear_model import LogisticRegression
import daal4py

X, y = make_classification(n_samples=200_000, n_features=100, n_informative=50,
                           n_classes=7, class_sep=0.5, random_state=0)
X *= np.logspace(-1, 1, X.shape[1])  # ill-conditioned: lbfgs needs many iterations
for n_threads in [1, 4, 14]:
    daal4py.daalinit(n_threads)
    clf = LogisticRegression(C=10.0, max_iter=1000).fit(X, y)
    print(f"{n_threads:2d} threads: n_iter={clf.n_iter_[0]}, train log loss={log_loss(y, clf.predict_proba(X)):.4f}")

 1 threads: n_iter=1000, train log loss=1.7761
 4 threads: n_iter=923, train log loss=1.7761
14 threads: n_iter=54, train log loss=1.7877

The 1-thread fit also warns that it hit max_iter. The 14-thread fit raises no warning.

Expected behavior

Repeated calls with the same inputs give the same value and gradient up to rounding, whatever the thread count. L-BFGS then converges to the same model as with 1 thread.

Impact on a real dataset

I found this in scikit-learn-benchmarks. The case is covtype (7 classes, 200k train rows, 100 Nystroem features), LogisticRegression(C=10, max_iter=1000), with test metrics averaged over 5 splits. On an Intel Core Ultra X7 358H laptop:

Stock scikit-learn on the same machine gets 0.607 at every core count. On this data, evaluating the loss twice at the same point inside the L-BFGS loop gave values up to ~2% apart and gradients up to ~60% apart.

It depends on the thread count, not on the P/E core mix. 8 E-cores alone show it, and daal4py.daalinit(4) with 14 CPUs available avoids it.

The same benchmark on an Intel Xeon 6787P shows no problem: 1000 iterations and identical metrics from 1 to 172 cores. So this may be specific to the AVX2 code path, but I haven't run the reproducer on that machine.

Environment:

• OS: Ubuntu, Linux 6.17
• CPU: Intel Core Ultra 5 225H (reproducer), Intel Core Ultra X7 358H (benchmark), both AVX2 without AVX-512
• scikit-learn-intelex 2026.1.0, daal 2026.1.0 (link version 20260100b'P'_20260601), tbb 2023.1.0, all from PyPI
• scikit-learn 1.8.0, numpy 2.5.3, scipy 1.18.1, Python 3.12.14