Я получил этот код, где расстояние - это треугольная треугольная матрица, определяемая следующим образом:

distance = np.tril(scipy.spatial.distance.cdist(points, points))  
def make_them_touch(distance):
    """
    Return the every distance where two points touched each other. See example below.
    """
    thresholds = np.unique(distance)[1:] # to avoid 0 at the beginning, not taking a lot of time at all
    result = dict()
    for t in thresholds:
            x, y = np.where(distance == t)
            result[t] = [i for i in zip(x,y)]
    return result

Моя проблема в том, что np.where довольно медленный с большой матрицей (например, 2000 * 100).
Как я могу ускорить этот код, улучшив np.where или изменив алгоритм?

EDIT: как указывал MaxU, лучшая оптимизация здесь НЕ должна генерировать квадрат матрицы и использовать итераторы.

Пример:

points = np.array([                                                                        
...: [0,0,0,0],                                                            
...: [1,1,1,1],         
...: [3,3,3,3],              
...: [6,6,6,6]                             
...: ])  

In [106]: distance = np.tril(scipy.spatial.distance.cdist(points, points))

In [107]: distance
Out[107]: 
array([[ 0.,  0.,  0.,  0.],
   [ 2.,  0.,  0.,  0.],
   [ 6.,  4.,  0.,  0.],
   [12., 10.,  6.,  0.]])

In [108]: make_them_touch(distance)
Out[108]: 
{2.0: [(1, 0)],
 4.0: [(2, 1)],
 6.0: [(2, 0), (3, 2)],
 10.0: [(3, 1)],
 12.0: [(3, 0)]}