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 = 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.