如何在Python中高效生成满足欧氏距离约束的随机向量对？

Your original approach works but is extremely inefficient because it relies on reject sampling—generating random pairs and discarding almost all of them when T is small. For T=0.05 and d=4, the probability that two random [-1,1]^d vectors are within T of each other is tiny, so your while loop runs thousands of iterations per pair. This is completely infeasible for N=1e6.

The Better Approach: Generate y Directly Around x #

Instead of guessing pairs, generate each x first, then create a y that's guaranteed to be within distance T of x (and still within [-1,1] for all components). Here's how to do it efficiently with vectorized NumPy operations (no slow Python loops!):

1. Generate all x vectors at once: Use NumPy's vectorized random functions to create N x d vectors in one go—this is way faster than looping.
2. Sample delta vectors uniformly within a d-dimensional ball of radius T: This ensures ||delta|| ≤ T, so y = x + delta will automatically satisfy ||x - y|| ≤ T.
3. Adjust y to stay within [-1,1]: Since T is small, most x + delta will already be within bounds. For the rare cases where a component goes out of range, we re-sample delta until it fits.

Implementation Code #

import numpy as np

def generate_close_vectors(N, d, T):
    # Step 1: Generate all x vectors (N x d) with components in [-1, 1]
    x = np.random.uniform(low=-1, high=1, size=(N, d))
    
    # Step 2: Generate delta vectors uniformly within radius T
    # Standard method for uniform d-dimensional ball sampling:
    # 1. Generate standard normal vectors
    delta = np.random.normal(size=(N, d))
    # 2. Normalize to unit vectors
    norms = np.linalg.norm(delta, axis=1, keepdims=True)
    delta_unit = delta / norms
    # 3. Scale by T * U(0,1)^(1/d) to get uniform distribution in the ball
    r = T * np.random.uniform(size=(N, 1)) ** (1/d)
    delta = delta_unit * r
    
    # Step 3: Compute y and fix boundary violations
    y = x + delta
    
    # Re-sample delta for any y components outside [-1,1]
    out_of_bounds = (y < -1) | (y > 1)
    while np.any(out_of_bounds):
        # Get indices of pairs needing rework
        idx = np.where(out_of_bounds.any(axis=1))[0]
        # Regenerate delta for these indices
        delta_new = np.random.normal(size=(len(idx), d))
        norms_new = np.linalg.norm(delta_new, axis=1, keepdims=True)
        delta_unit_new = delta_new / norms_new
        r_new = T * np.random.uniform(size=(len(idx), 1)) ** (1/d)
        delta[idx] = delta_unit_new * r_new
        # Update y for these indices
        y[idx] = x[idx] + delta[idx]
        # Check again for bounds
        out_of_bounds = (y < -1) | (y > 1)
    
    # Return as list of pairs (or keep as numpy arrays for maximum efficiency)
    return list(zip(x, y))

# Example usage
d = 4
T = 0.05
N = 10**6  # Works smoothly even for 1 million pairs!
pairs = generate_close_vectors(N, d, T)

# Verify a few pairs
for i in range(5):
    x, y = pairs[i]
    print(f"Pair {i}: Distance = {np.linalg.norm(x-y):.6f}")

Why This Is Fast #

• Vectorized Operations: All random generation happens in bulk with NumPy, which uses optimized C backend code instead of slow Python loops.
• No Wasted Samples: We directly generate valid y vectors instead of discarding invalid pairs. The only looping is for rare boundary cases, which is negligible when T is small.
• Sound Statistical Sampling: The delta generation method ensures uniform distribution within the radius-T ball around each x, so your pairs are statistically valid.

Performance Note #

For N=1e6, this code will run in seconds (not hours like your original approach). If you can work directly with the NumPy arrays x and y instead of converting to a list of pairs, skip the list(zip(...)) step for even better speed.

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