Кластерный анализ облака точек в Python - определение кластеров из двоичной матрицы

В качестве примера у меня есть следующие входные данные (облако точек, с которым я работаю, более сложное)

Data = [1,1,1,1,1],[1,1,2,1,1],[2,2,2,1,1],[3,3,3,1,1],[4,4,4,1,1],[5,5,5,1,1],[50,50,50,1,1],[95,95,95,1,1],[96,96,96,1,1],[97,97,97,1,1],[98,98,98,1,1],[99,99,99,1,1],[2,2,3,1,1],[2,2,1,1,1],[2,2,4,1,1]

Алгоритм кластеризации дает двоичную верхнетреугольную матрицу (назовем ее матрицей соединений). 1 означает, что две точки соединены. Например. Идентификатор точки 0 (строка 0) связан с самим собой (столбец 0) и 1,2,3,12,13,14. Но в точки 4 и 5 также можно попасть через 3, 12, 13 и 14.

[ 1.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  0.  1.  1.  1.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  0.  0.  1.  1.  0.  0.  0.  0.  0.  0.  1.  0.  1.]
[ 0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.]

Я могу идентифицировать кластеры в каждой строке с помощью кластеризации строк, где s - это двоичная матрица сверху.

def rowclustering(s):
    r = 0
    idx = []
    while r < size(s,0):
        row = []
        for i in range(size(s,1)):
            if s[r][i] == 1:
                row = row + [i]
        r = r + 1
        idx = idx + [row]
    return idx

И возвращаемый idx:

idx = [[0, 1, 2, 3, 12, 13, 14], [1, 2, 3, 12, 13, 14], [2, 3, 4, 12, 13, 14], [3, 4, 5, 12, 13, 14], [4, 5, 12, 14], [5], [6], [7, 8, 9], [8, 9, 10], [9, 10, 11], [10, 11], [11], [12, 13, 14], [13, 14], [14]]

Теперь очевидно, что кластеров меньше, чем 15, потому что некоторые строки связаны через общий идентификатор (например, посмотрите на идентификаторы 4 и 5). Я хочу иметь:

result = [[0, 1, 2, 3, 4, 5, 12, 13, 14], [6], [7, 8, 9, 10, 11]]

Я попытался создать функцию (кластеризацию (idx, f)), которая делает это, где idx является результатом кластеризации строк, а f будет первой строкой в ​​idx, например [0, 1, 2, 3, 12, 13, 14]. Однако эта функция не будет завершена должным образом. Какой будет правильный код для прерывания функции после того, как все соединения (идентификаторов idx) будут выполнены?

def clustering(idx,f):
    for i in f:
        f = f + idx[i]
    f = list(set(f))
    clustering(idx,f)

    return

Проблема, которую я пытаюсь решить, - это своего рода процедура саморазвития. Функция кластеризации должна вызывать сама себя до тех пор, пока не будут выполнены все возможные точечные соединения. Это можно сделать на idx или (может быть, лучше) на матрице соединений (уменьшение матрицы?).

Любая помощь высоко ценится! Дайте мне знать, если мне нужно уточнить свой вопрос. Спасибо.


python cluster-analysis matrix point-clouds
person Michael    schedule 13.12.2016    source источник
comment
Ваша функция кластеризации (idx, f) никогда не может вернуться, она будет рекурсивно повторяться до тех пор, пока не произойдет переполнение стека.   -  person benathon    schedule 13.12.2016
comment
Взгляните на DBSCAN и одноканальный. Но вам нужно connected = 0, а не 1.   -  person Has QUIT--Anony-Mousse    schedule 15.12.2016

Ответы (1)


arrow_upward
2
arrow_downward

Ваша проблема может рассматриваться как поиск подключенных компонентов. Вы можете использовать networkx для получения решения или самостоятельно реализовать BFS (поиск в ширину).

import networkx as nx
import numpy as np
x = """
[ 1.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  0.  1.  1.  1.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  0.  0.  1.  1.  0.  0.  0.  0.  0.  0.  1.  0.  1.]
[ 0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  0.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.]
"""
mat = eval('[' + x.replace('.', '.,').replace(']', '],') + ']')
mat = np.array(mat)
G = nx.from_numpy_matrix(np.array(mat))
print(list(nx.connected_components(G)))
[{0, 1, 2, 3, 4, 5, 12, 13, 14}, {6}, {7, 8, 9, 10, 11}]

РЕДАКТИРОВАТЬ:

Собственно, что-то в этой проблеме заставило меня вспомнить кое-что, что я читал некоторое время назад. Фактически это можно вычислить, используя только матричные операции. Это очень здорово. Ваша исходная матрица - это матрица смежности (A), и нам также нужно указать матрицу степеней (D), которая содержит степень каждого узла на диагонали. Мы можем использовать их для определения матрицы Лапласа (L), а затем использовать немного теории спектральных графов. (ура!)

# Make the undirected version of the graph (no self loops)
A = (mat + mat.T) * (1 - np.eye(mat.shape[0]))
# Make the degree matrix
D = np.diag(A.sum(axis=1) + A.sum(axis=0)) / 2
# thats all we need to define the laplacian
L = D - A

# The number of zeros eigenvalues of the Laplacian is exactly the number of CCs
np.isclose(np.linalg.eigvals(L), 0).sum()

3 

# The connected compoments themselves are identified by rows that have the same nullspace vector
u, s, vh = np.linalg.svd(L)
ns = vh[(s >= 1e-13).sum():].conj().T

array([[-0.32684842, -0.06239247, -0.0197079 ],
   [-0.32684842, -0.06239247, -0.0197079 ],
   [-0.32684842, -0.06239247, -0.0197079 ],
   [-0.32684842, -0.06239247, -0.0197079 ],
   [-0.32684842, -0.06239247, -0.0197079 ],
   [-0.32684842, -0.06239247, -0.0197079 ],
   [-0.19222441,  0.97663659,  0.09607676],
   [-0.01778075, -0.04721352,  0.44435878],
   [-0.01778075, -0.04721352,  0.44435878],
   [-0.01778075, -0.04721352,  0.44435878],
   [-0.01778075, -0.04721352,  0.44435878],
   [-0.01778075, -0.04721352,  0.44435878],
   [-0.32684842, -0.06239247, -0.0197079 ],
   [-0.32684842, -0.06239247, -0.0197079 ],
   [-0.32684842, -0.06239247, -0.0197079 ]])

Теперь мы вычислили ответ! Это просто немного странно интерпретировать. Небольшая обработка может превратить это в желаемое представление.

# the following is a little numpy trick to find unique rows 
# chopping off the last few decimal places to account for numerical errors
ns_ = np.ascontiguousarray(np.round(ns, 8)).view(np.dtype((np.void, ns.dtype.itemsize * ns.shape[1])))
ns_basis, row_to_cc_id = np.unique(ns_, return_inverse=True)
# Finally we can just use this to convert to the desired output format
groups = [[] for _ in range(len(ns_basis))]
for row, id in enumerate(row_to_cc_id):
    groups[id].append(row)
print(groups)
[[0, 1, 2, 3, 4, 5, 12, 13, 14], [6], [7, 8, 9, 10, 11]]
person Erotemic    schedule 13.12.2016
comment
Ваш ответ отличный! Спасибо, это именно то, что я искал. Я взял вашу отредактированную версию, и она работает как шарм. Не могли бы вы мне сказать, как масштабируется этот алгоритм, что, если у меня будет 10 000 или 1 000 000 баллов? - person Michael; 14.12.2016
comment
Первый ответ - просто использовать алгоритм поиска в ширину для поиска связанных компонентов en.wikipedia.org/ wiki / Connected_component_ (graph_theory). Этот алгоритм имеет линейное время O (n) и очень быстрый. Возможно, вы захотите рассмотреть возможность его повторной реализации, чтобы вам не приходилось выполнять преобразование в граф networkx. Второй алгоритм больше для развлечения. Вызов SVD (разложение по сингулярным значениям) занимает примерно ~ O (n ^ 3) времени, поэтому лучше не использовать эту версию на практике. - person Erotemic; 14.12.2016
comment
Кроме того, вам следует подумать о хранении ваших данных в разреженном формате вместо плотной матрицы. Список смежности (en.wikipedia.org/wiki/Adjacency_list) идеально подходит для BFS. алгоритм. (Или вы также можете использовать поиск в глубину (DFS), они делают более или менее то же самое) - person Erotemic; 14.12.2016