Есть такой алгоритм. В целях масштабирования, я хочу его немного перепилить, но что-то под ночь совсем голова закипает вот на этом куске:
I = np.where(np.diag(A + R) > 0)[0]
K = I.size # Identify exemplars
if K > 0:
c = np.argmax(S[:, I], axis=1)
c[I] = np.arange(K) # Identify clusters
# Refine the final set of exemplars and clusters and return results
for k in range(K):
ii = np.where(c == k)[0]
j = np.argmax(np.sum(S[ii[:, np.newaxis], ii], axis=0))
I[k] = ii[j]
c = np.argmax(S[:, I], axis=1)
c[I] = np.arange(K)
labels = I[c]
# Reduce labels to a sorted, gapless, list
cluster_centers_indices = np.unique(labels)
labels = np.searchsorted(cluster_centers_indices, labels)
S, A, R - квадратные несимметричные матрицы NxN
I = np.where(np.diag(A + R) > 0)[0]
K = I.size # Identify exemplars
c = np.argmax(S[:, I], axis=1)
S[:, I] - матрица, составленная из столбцов, индексы которых мы получили на предыдущем шаге (Причем, так как K <=N, то это уже не NxN, а NxK)
c - одномерный массив c индексами максимумов в строках матрицы NxK
А вот дальше совсем беда.
Если поможет, могу показать ссылку на исходный код на матлабе, с которого переписывался код по ссылке, и даже статью.
Заранее благодарю за помощь.