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:
Choose a random point from all the \$10^6\$ points, add it to the chosen set of points
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
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.
data
example, they are not the same length, so points in data aren't part of a sphere. Is the example data not good or is there something that doesn't work? \$\endgroup\$