Choosing evenly distributed points from a million points on a sphere

Question

I have \$10^5\$ to \$10^6\$ points on a sphere, and want to choose some points from them which are as close as uniformly distributed as possible. For that reason, I do the following: at each step I add the point in the data which is farthest away from any point which I have already chosen. In order to do this, I find the closest distance from each point in the data to any of the points that are already chosen. Then I pick the point from the data which has the largest minimum distance to any of the points that are already chosen. The algorithm is the following:

1. Choose a random point from all the \$10^6\$ points, add it to the chosen set of points

2. For each of the \$10^6\$ points:

   2.1 Find the distances to all the points which are chosen already

   2.2 Find the min distance to a point which is chosen already

   2.3 Find the point which is furthest away from any point which is chosen already

   2.4 Add this point to the chosen points

3. Repeat 2 until I have enough points (say 1000 or 10000 or 100000).

data is \$10^6\$ by 3 array of points on a sphere, e.g.

[[ 0.26750522 -0.92342735  0.27517791]
 [-0.26753053  0.9228284  -0.27715548]
 [ 0.34058837  0.35873472  0.86908513]
 [-0.33563178 -0.35794416 -0.8713365 ]
 [-0.8945523  -0.15682272  0.41854847]
 [ 0.906739    0.15103321 -0.39371734]
 [ 0.49138659 -0.64830133 -0.581588  ]
 [-0.87922161  0.24492971 -0.40863039]
 [ 0.18062012  0.39852148 -0.89919798]
 [-0.49103509  0.65872966  0.57005243]
 [-0.55615839 -0.82248645 -0.11918007]
 [ 0.85802915 -0.25207255  0.44748789]
 [-0.16990087 -0.40390349  0.89888579]]

I have implemented this in Python:

import random
import numpy

N = 1000

starting_point_ind = random.SystemRandom().randint (1 , len(data))

points_from_data = numpy.array(numpy.zeros((N, 3)), dtype=float, order='C')
points_from_data[0] = data[starting_point_ind]

for i in range(1, N):
    distancesToClosestPoints = numpy.array(numpy.zeros(len(data)), dtype=float)

    for j, point in enumerate(data):
        distancesToEstablishedPoints = numpy.linalg.norm (point - points_from_data[0:i] , axis=1)

        distancesToClosestPoints[j] = distancesToEstablishedPoints.min()


    k = distancesToClosestPoints.argmax()

    points_from_data[i] = data[k]

This doesn't run very efficiently for large N, so I also tried to optimise it by doing

import random
from sklearn.metrics import pairwise_distances
import time

t0 = time.time()

N = 1000

starting_point_ind = random.SystemRandom().randint (1 , len(data))

points_from_data = numpy.array([data[starting_point_ind]])

for i in range(1, N):
    pairwise_distances_to_data = pairwise_distances (data, Y=points_from_data, metric='euclidean', n_jobs=-1)
    pairwise_distances_to_data = numpy.array(pairwise_distances_to_data)

    min_distances_to_data = numpy.amin(pairwise_distances_to_data, axis=1)

    k = min_distances_to_data.argmax()

    points_from_data = numpy.append(points_from_data, [data[k]], axis=0)

print(points_from_data)

But this still has the issue that the pairwise distances are computed each time, and has obvious memory issues for large N. I would like to hear suggestions how I can improve the performance of the algorithm.

Answer 1

Since we can assume all your points are on a sphere, there is a better way to compute distance than using the L2 Norm (Euclidean distance, which you're using).

If we transform your data points to Spherical coordinates, we could use a simple Manhattan distance on the theta and rho coordinates, which is less expensive to compute. Since you compute this distance about \$Nm\$ times (6 million points * 1000 sample points), I'm sure it would make a good performance difference.

Small note on the Spherical Coordinates :

Instead of having x,y,z coordinates, we use r, theta, rho, which correspond to :

• r : the length of the vector
• theta : angle on the XY plane
• phi : angle of inclination starting from theta

So, to change data to spherical coordinates, we'd do the following (I assume data is arranged in a x,y,z manner). The conversion is very fast, but it could be further optimized if needed :

def convert_to_spherical(data):
    spherical_data = np.zeros_like(data)
    for i,d in enumerate(data):
        r = math.sqrt(d[0]**2 + d[1]**2 + d[2]**2)
        theta = math.acos(d[2] / r)
        phi = math.atan(d[1] / d[0])

        spherical_data[i] = np.array([r, theta, phi])

    return spherical_data

Then, you could use this instead of data and change the metric used by pairwise_distances to manhattan. You'll need to dig a little bit to understand why this works, but consider that since the parameter r is always the same, you can simply compute the difference between the angles.

Since my laptop is pretty bad, I couldn't generate data to sample 1000 points out of 1000000. These times were generated sampling 300 points from 100000 points, I started timing after generating the points but I included the coordinates transformation in the times :

Time using your optimized solution : 213.1844208240509s
Time using the spherical coordinates : 140.83173203468323s

Considering this pretty big time difference, I'm pretty sure performance would be much better on your data sample.

If afterwards, you need to get your coordinates back to the cartesian plane, the formula is pretty easy and you can find it online easily.

I've generated the points using thanks to this answer :

def sample_spherical(npoints, ndim=3):
    vec = numpy.random.randn(ndim, npoints)
    vec /= numpy.linalg.norm(vec, axis=0)
    return numpy.transpose(vec)