用numpy快速计算平方欧氏距离的数值问题及优化方案咨询

Great question—this is a classic tradeoff between raw speed and numerical stability in distance calculations, and you’ve nailed exactly why that formula has quirks! Let’s break this down.

Why the Numerical Errors Happen #

Mathematically, the formula ||x - y||² = ||x||² + ||y||² - 2*(x @ y.T) is perfect. But when working with floating-point numbers (like float64), we run into catastrophic cancellation: when two very similar large numbers are subtracted, the small differences get amplified, leading to tiny errors (or even negative values) where we expect zero. This is especially bad for diagonal elements (distance from a point to itself) or nearly identical points.

Fixes That Keep Speed (Almost) Unchanged #

You don’t have to give up the 10x speed boost—here are simple, fast fixes:

1. Force Diagonal Elements to Zero
   The diagonal of your result matrix must be zero (distance from a point to itself). This is a guaranteed correction, and it’s a single fast vectorized operation:

   np.fill_diagonal(res, 0.0)
   

2. Truncate Tiny Errors to Zero
   Any value with an absolute magnitude smaller than a reasonable threshold (like 1e-10) is just floating-point noise. We can zero these out with boolean indexing (another fast numpy operation):

   res[np.abs(res) < 1e-10] = 0.0
   

3. Safe Square Root Calculation
   If you need the actual Euclidean distance (not squared), take the absolute value first to avoid nan from negative noise:

   sqrt_res = np.sqrt(np.abs(res))
   

   For valid distances, this won’t change anything—only fixes the spurious negatives.

4. Upgrade to Higher Precision
   If your data is in float32, switching to float64 (numpy’s default for rand) reduces error magnitude significantly. For extreme cases, you could even use np.longdouble, though this has minor memory/performance costs.

Improved Example Code #

Here’s how this looks in practice:

import numpy as np
M = np.random.rand(1000, 10)
M_norm = np.sum(M**2, axis=1)
res = M_norm[:, np.newaxis] + M_norm[np.newaxis, :] - 2. * M @ M.T

# Apply fixes
np.fill_diagonal(res, 0.0)
res[np.abs(res) < 1e-10] = 0.0
sqrt_res = np.sqrt(np.abs(res))

unique_diag = np.unique(np.diag(res))
sqrt_unique_diag = np.sqrt(unique_diag)
print(unique_diag)
print(sqrt_unique_diag)

This will output exactly [0.] for the diagonal, with no nan values.

Why Scipy’s pdist Doesn’t Use This Method #

You’re exactly right—numerical stability is the reason! scipy.pdist calculates squared Euclidean distance directly as sum((x-y)**2) for each pair. This is slower, but it avoids catastrophic cancellation because it computes the difference first (small numbers) before squaring, which preserves precision for similar points. Scipy prioritizes correctness over raw speed for its core functions, which makes sense for general-purpose use.

These fixes let you keep the speed of the numpy formula while mitigating the worst numerical issues—perfect when you need performance but can tolerate a tiny bit of post-processing!

内容的提问来源于stack exchange，提问作者no_use123