Numpy - Clustering - Distance - Vectoriation

Question

Я собрал образец данных (400 k образцов, размер = 205, 200 кластеров), используя sklearn Kmeans.

Я хочу знать для каждого кластера максимальное расстояние между центром кластера и самым отдаленным образцом кластера, чтобы узнать что-то о «размере» кластера. Вот мой код:

import numpy as np 
import scipy.spatial.distance as spd 
diam = np.empty([200]) 
for i in range(200): 
    diam[i] = spd.cdist(seed[np.newaxis, i, 1:], data[data[:, 0]==i][:,1:]).max() 

"семена" являются кластерные центры (200x206). В первом столбце «семя» содержится количество выборок внутри кластера (здесь неактуально).

«данные» - это образцы (400kx206). Первый столбец данных содержит номер кластера.

Вопрос: Это делается с использованием петли (не так «numpy»). Можно ли «векторизовать» его?

Answer 1

Вот один Векторизованных подход -

# Sort data w.r.t. col-0 
data_sorted = data[data[:, 0].argsort()] 

# Get counts of unique tags in col-0 of data and repeat seed accordingly. 
# Thus, we would have an extended version of seed that matches data's shape. 
seed_ext = np.repeat(seed,np.bincount(data_sorted[:,0]),axis=0) 

# Get euclidean distances between extended seed version and sorted data 
dists = np.sqrt(((data_sorted[:,1:] - seed_ext[:,1:])**2).sum(1)) 

# Get positions of shifts in col-0 of sorted data 
shift_idx = np.append(0,np.nonzero(np.diff(data_sorted[:,0]))[0]+1) 

# Final piece of puzzle is to get tag based maximum values from dists, 
# where each tag is unique number in col-0 of data 
diam_out = np.maximum.reduceat(dists,shift_idx) 

время выполнения тестов и проверить выход -

Определение функции:

def loopy_cdist(seed,data): # from OP's solution 
    N = seed.shape[0] 
    diam = np.empty(N) 
    for i in range(N): 
     diam[i]=spd.cdist(seed[np.newaxis,i,1:], data[data[:,0]==i][:,1:]).max() 
    return diam 

def vectorized_repeat_reduceat(seed,data): # from this solution 
    data_sorted = data[data[:, 0].argsort()] 
    seed_ext = np.repeat(seed,np.bincount(data_sorted[:,0]),axis=0) 
    dists = np.sqrt(((data_sorted[:,1:] - seed_ext[:,1:])**2).sum(1)) 
    shift_idx = np.append(0,np.nonzero(np.diff(data_sorted[:,0]))[0]+1) 
    return np.maximum.reduceat(dists,shift_idx) 

def vectorized_indexing_maxat(seed,data): # from @moarningsun's solution 
    seed_repeated = seed[data[:,0]] 
    dist_to_center = np.sqrt(np.sum((data[:,1:]-seed_repeated[:,1:])**2, axis=1)) 
    diam = np.zeros(len(seed)) 
    np.maximum.at(diam, data[:,0], dist_to_center) 
    return diam 

Verify выходы:

In [417]: # Inputs 
    ...: seed = np.random.rand(20,20) 
    ...: data = np.random.randint(0,20,(40000,20)) 
    ...: 

In [418]: np.allclose(loopy_cdist(seed,data),vectorized_repeat_reduceat(seed,data)) 
Out[418]: True 

In [419]: np.allclose(loopy_cdist(seed,data),vectorized_indexing_maxat(seed,data)) 
Out[419]: True 

Runtimes:

In [420]: %timeit loopy_cdist(seed,data) 
10 loops, best of 3: 35.9 ms per loop 

In [421]: %timeit vectorized_repeat_reduceat(seed,data) 
10 loops, best of 3: 28.9 ms per loop 

In [422]: %timeit vectorized_indexing_maxat(seed,data) 
10 loops, best of 3: 24.1 ms per loop 

Answer 2

Очень похожая идея, как @Divakar, но без сортировки:

seed_repeated = seed[data[:,0]] 
dist_to_center = np.sqrt(np.sum((data[:,1:]-seed_repeated[:,1:])**2, axis=1)) 

diam = np.zeros(len(seed)) 
np.maximum.at(diam, data[:,0], dist_to_center) 

ufunc.at, как известно, медленно, хотя, так что было бы интересно посмотреть, что быстрее.

Answer 3

Мы можем быть немного умнее об индексации и экономить около ~ 4 раза в цене.

Первая позволяет строить некоторые данные правильной формы:

seed = np.random.randint(0, 100, (200,206)) 
data = np.random.randint(0, 100, (4e5,206)) 
seed[:, 0] = np.arange(200) 
data[:, 0] = np.random.randint(0, 200, 4e5) 
diam = np.empty(200) 

Время оригинального ответа:

%%timeit 
for i in range(200): 
    diam[i] = spd.cdist(seed[np.newaxis, i, 1:], data[data[:, 0]==i][:,1:]).max() 

1 loops, best of 3: 1.35 s per loop 

ответ moarningsun в:

%%timeit 
seed_repeated = seed[data[:,0]] 
dist_to_center = np.sqrt(np.sum((data[:,1:]-seed_repeated[:,1:])**2, axis=1)) 
diam = np.zeros(len(seed)) 
np.maximum.at(diam, data[:,0], dist_to_center) 

1 loops, best of 3: 1.33 s per loop 

ответ Divakar в:

%%timeit 
data_sorted = data[data[:, 0].argsort()] 
seed_ext = np.repeat(seed,np.bincount(data_sorted[:,0]),axis=0) 
dists = np.sqrt(((data_sorted[:,1:] - seed_ext[:,1:])**2).sum(1)) 
shift_idx = np.append(0,np.nonzero(np.diff(data_sorted[:,0]))[0]+1) 
diam_out = np.maximum.reduceat(dists,shift_idx) 

1 loops, best of 3: 1.65 s per loop 

Как мы видим, на самом деле ничего не получилось с помощью векторизованных решений, помимо большего объема памяти.Чтобы избежать этого, мы должны вернуться к первоначальному ответу, который действительно правильный способ сделать это, и вместо того, чтобы попытаться уменьшить количество индексации:

%%timeit 
idx = data[:,0].argsort() 
bins = np.bincount(data[:,0]) 
counter = 0 
for i in range(200): 
    data_slice = idx[counter: counter+bins[i]] 
    diam[i] = spd.cdist(seed[None, i, 1:], data[data_slice, 1:]).max() 
    counter += bins[i] 

1 loops, best of 3: 281 ms per loop 

Дважды проверьте ответ:

np.allclose(diam, dam_out) 
True 

Это проблема с допущением, что петли python плохи. Они часто бывают, но не во всех ситуациях.