LocalOutlierFactor with Mahalanobis distance returns different results based on `n_jobs` parameter

### Describe the bug

I encountered the following bug while doing an outlier analysis on a large dataset:

To detect outliers in a dataset with known outliers, I followed these steps:

1) compute the covariance for a set of preselected numerical features with either `MinCovDet` or `EmpiricalCovariance` and extracted the inverse matrix VI using the `precision_` parameter.
2) I initialized the LocalOutlierFactor class with the metric `mahalanobis`, passed `VI`, and set n_jobs to 1
3) I used `fit_predict()`to get the outlier predictions

Everything looked as expected.

Later on, I tried speeding up the computation by adding a few cores to `n_jobs` and had to realise that no outliers were detected anymore. After some debugging, I found that:

1) outlier detection seems not to work when using multiple cores (all samples are marked as inliers)
2) this effect seems to be tied to `mahalanobis` distance, because the bug does not appear with minkowski (p=2) distance regardless of the number of cores.

One workaround seems to be to precompute the mahalanobis distance and passing it to the LOF, but this is again using only one core for the computation and is not feasible for large datasets due to the size of the full pairwise distance matrix.

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.covariance import MinCovDet
from sklearn.neighbors import LocalOutlierFactor
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt

RNG = np.random.RandomState(42)


def make_data(n_inliers=200, n_outliers=10, n_features=10):
    # Inliers: standard normal in 10D
    X_in = RNG.normal(loc=0.0, scale=1.0, size=(n_inliers, n_features))
    # Outliers: shifted by +8 in all dimensions
    X_out = RNG.normal(loc=8.0, scale=1.0, size=(n_outliers, n_features))
    X = np.vstack([X_in, X_out])
    y_true = np.hstack([np.ones(n_inliers, dtype=int), -np.ones(n_outliers, dtype=int)])
    return X, y_true


def robust_inverse_covariance(X, random_state=0):
    # Fit robust covariance estimator
    mcd = MinCovDet(random_state=random_state).fit(X)
    cov = mcd.covariance_  # shape (n_features, n_features)
    # Compute inverse covariance (VI) expected by 'mahalanobis' metric
    VI = mcd.precision_ 
    # ensure dtype=float64 and C-contiguous so it serializes safely across processes
    VI = np.ascontiguousarray(VI, dtype=np.float64)
    return VI, cov


def run_lof(X, VI, n_neighbors=20, n_jobs=1):
    lof = LocalOutlierFactor(
        n_neighbors=n_neighbors,
        metric='mahalanobis',
        metric_params={'VI': VI},
        n_jobs=n_jobs
    )
    # fit_predict returns -1 for outliers and 1 for inliers
    labels = lof.fit_predict(X)
    return labels


def run_lof_precomputed(X, VI, n_neighbors=20, n_jobs=1):
    # compute full pairwise Mahalanobis distance matrix in the main process
    D = cdist(X, X, metric='mahalanobis', VI=VI)
    lof = LocalOutlierFactor(
        n_neighbors=n_neighbors,
        metric='precomputed',
        n_jobs=n_jobs
    )
    labels = lof.fit_predict(D)
    return labels


def stats(labels):
    unique, counts = np.unique(labels, return_counts=True)
    return dict(zip(unique.tolist(), counts.tolist()))


# Executing example
## make data
X, y_true = make_data(n_inliers=200, n_outliers=10, n_features=10)

## plot data
plt.figure(figsize=(8, 6))
plt.scatter(X[:, 0], X[:, 1], c=y_true, cmap='coolwarm', edgecolor='k')
plt.title("Toy dataset with inliers (blue) and outliers (red)")
plt.xlabel("Feature 0")
plt.ylabel("Feature 1")
plt.show()

## Fit MinCovDet on the whole dataset 
VI, cov = robust_inverse_covariance(X, random_state=42)

## Run LOF with a single core (n_jobs=1)
labels_n1 = run_lof(X, VI, n_neighbors=20, n_jobs=1)
print("n_jobs=1 unique counts:", stats(labels_n1))
## expected result: n_jobs=1 unique counts: {-1: 11, 1: 199}

## Run LOF allowing multiple cores (n_jobs=4)
labels_n4 = run_lof(X, VI, n_neighbors=20, n_jobs=4)
print("n_jobs=4 unique counts:", stats(labels_n4))
## exprected result: n_jobs=4 unique counts: {1: 210}


## fix attempt: precompute Mahalanobis with single core and LOF with multiple
labels_precomputed = run_lof_precomputed(X, VI, n_neighbors=20, n_jobs=4)
print("precomputed (n_jobs=4) unique counts:", stats(labels_precomputed))
## expected result: precomputed (n_jobs=4) unique counts: {-1: 11, 1: 199}

## Print outlier labels
pred_outliers = np.where(labels_n1 == -1)[0]
print("Indices detected as outliers by n_jobs=1", pred_outliers)
## expected output: Indices detected as outliers by n_jobs=1: [116 200 201 202 203 204 205 206 207 208 209]
y_true_outliers = np.where(y_true == -1)[0]
print("Indices of true outliers:", y_true_outliers)
## expected output: Indices of true outliers: [200 201 202 203 204 205 206 207 208 209]


```

### Expected Results

It shouldn't matter how many cores are passed to the 'n_jobs' parameters. The results should be comparable, or at least detect some outliers.

### Actual Results

Results for LOF should be comparable, but are wildly different based on changes to the 'n_jobs' parameter, with all other parameters staying the same:

Results for n_jobs=1:
n_jobs=1 unique counts: {-1: 11, 1: 199}

Results for n_jobs=4:
n_jobs=4 unique counts: {1: 210}

### Versions

```shell
Python dependencies:
      sklearn: 1.7.2
          pip: 25.2
   setuptools: 80.9.0
        numpy: 2.3.4
        scipy: 1.15.2
       Cython: None
       pandas: 2.2.3
   matplotlib: 3.10.6
       joblib: 1.5.2
threadpoolctl: 3.6.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 5
         prefix: libscipy_openblas
       filepath: /home/USER/environments/conda_environments/MY_ENV/lib/python3.11/site-packages/numpy.libs/libscipy_openblas64_-8fb3d286.so
        version: 0.3.30
threading_layer: pthreads
   architecture: Haswell

       user_api: blas
   internal_api: openblas
    num_threads: 5
         prefix: libopenblas
       filepath: /home/USER/environments/conda_environments/MY_ENV/lib/libopenblasp-r0.3.30.so
        version: 0.3.30
threading_layer: pthreads
   architecture: Zen

       user_api: openmp
   internal_api: openmp
    num_threads: 5
         prefix: libgomp
       filepath: /home/USER/environments/conda_environments/MY_ENV/lib/libgomp.so.1.0.0
        version: None
```

COMMENT:
The non-thread-safe object seems to be the `self.buffer` for the Mahalanobis Distance implementation in `sklearn/metrics/_dist_metrics.pyx.tp` excerpted below.
I imagine the class would have to be rewritten without that array and with multithreading in mind?

``` python
### sklearn\metrics\_dist_metrics.pyx.tp | line 1600
#------------------------------------------------------------
# Mahalanobis Distance
#  d = sqrt( (x - y)^T V^-1 (x - y) )
cdef class MahalanobisDistance{{name_suffix}}(DistanceMetric{{name_suffix}}):
	""" ... """
    cdef float64_t[::1] buffer

    def __init__(self, V=None, VI=None):
        if VI is None:
            if V is None:
                raise ValueError("Must provide either V or VI "
                                 "for Mahalanobis distance")
            VI = np.linalg.inv(V)
        if VI.ndim != 2 or VI.shape[0] != VI.shape[1]:
            raise ValueError("V/VI must be square")

        self.mat = np.asarray(VI, dtype=np.float64, order='C')

        self.size = self.mat.shape[0]

        # We need to create a buffer to store the vectors' coordinates' differences
        self.buffer = np.zeros(self.size, dtype=np.float64)
```

Labels: Bug, Needs Investigation