Generating N random unit vectors with their sum equal to 0 (Python)

Viewed 234

I'd like to generate N random 3-dimensional vectors (uniformly) on the unit sphere but with the condition, that their sum is equal to 0. My attempt was to generate N/2 random unit vectors, while the other are just the same vectors with a minus sign. The problem is, as I'm trying to achieve as little correlation as possible, and my idea is obviously not ideal, since half of my vectors are perfectly anti-correlated with their corresponding pair.

2 Answers

Your problem does not really have a solution, but you can generate a set of vectors that are going to be slightly less visibly correlated than your original solution of negating them. To be precise, if you generate N / 2 vectors and negate them, then rotate the negated vectors about their sum by any angle, you can guarantee that the sum will be zero and the correlation will be a more complicated rotation than a negative identity matrix.

import numpy as np
from scipy.spatial.transform import Rotation

N = 10

v1 = np.random.normal(size=(N / 2, 3))
v1 /= np.linalg.norm(v1, axis=1, keepdims=True)

axis = v1.sum(0)
rot = Rotation.from_rotvec(np.random.uniform(2.0 * np.pi) * axis / np.linalg.norm(axis))

v2 = rot.apply(-v1)

result = np.concatenate((v1, v2), axis=0)

This assumes that N is even in all cases. The normal distribution is a fairly standard method to generate points uniformly on a sphere: https://mathworld.wolfram.com/SpherePointPicking.html.

If you had some leeway from the sum being exactly zero, you could align two random sets of N / 2 vectors so that their sums point opposite each other.

In this code, I tried to generate vectors selected from a sphere by converting a theta, phi to x, y, z.

Cartesian System Coordinates

import numpy as np

def vectorize(theta, phi):
    x = np.cos(phi) * np.cos(theta)
    y = np.cos(phi) * np.sin(theta)
    z = np.sin(phi)
    return np.array([x, y, z])

theta_range = np.arange(0, 2 * np.pi, 0.01)
phi_range = np.arange(-np.pi / 2, np.pi / 2, 0.01)
TH, PI = np.meshgrid(theta_range, phi_range)
whole_map = np.vstack((TH.flatten(), PI.flatten())).T

# Number of vectors: 
N = 100

# Selecting N/2 Vectors first at random
v_selected = np.random.choice(range(whole_map.shape[0]), N // 2)

vectors = np.array([vectorize(whole_map[ind][0], whole_map[ind][1]) for ind in v_selected])

# Doubling up the number of vectors by adding the negate of each vector to the vector set
vectors = np.vstack((vectors, -vectors))


print(vectors.sum(axis=0))
# array([1.94289029e-16, 1.17961196e-16, 1.11022302e-16])
# Almost 0, but isn't zero because of floating number precision when converted to binary

Here is the scatter plot of the points generated on the sphere with radius=1: enter image description here

Related