Вычислить разреженное транзитивное закрытие scipy разреженной матрицы

Это было поднято в системе отслеживания проблем SciPy. Проблема не столько в формате вывода; реализация Флойда-Уоршалла состоит в том, чтобы начать с матрицы, полной бесконечностей, а затем вставить конечные значения, когда путь найден. Разреженность сразу теряется.

Библиотека networkx предлагает альтернативу с ее all_pairs_shortest_path_length. Его вывод — итератор, возвращающий кортежи вида

(source, dictionary of reachable targets) 

для преобразования в разреженную матрицу SciPy требуется небольшая работа (формат csr здесь естественен). Полный пример:

import numpy as np
import networkx as nx
import scipy.stats as stats
import scipy.sparse as sparse

A = sparse.random(6, 6, density=0.2, format='csr', data_rvs=stats.randint(1, 2).rvs).astype(np.uint8)
G = nx.DiGraph(A)       # directed because A need not be symmetric
paths = nx.all_pairs_shortest_path_length(G)
indices = []
indptr = [0]
for row in paths:
  reachable = [v for v in row[1] if row[1][v] > 0]
  indices.extend(reachable)
  indptr.append(len(indices))
data = np.ones((len(indices),), dtype=np.uint8)
A_trans = A + sparse.csr_matrix((data, indices, indptr), shape=A.shape)
print(A, "\n\n", A_trans)

Причина добавления А обратно следующая. Выход Networkx включает в себя пути длины 0, которые сразу заполнили бы диагональ. Мы не хотим, чтобы это произошло (вы хотели транзитивного закрытия, а не рефлексивно-транзитивного закрытия). Отсюда строка reachable = [v for v in row[1] if row[1][v] > 0]. Но тогда мы вообще не получаем никаких диагональных записей, даже там, где они были у A (пустой путь нулевой длины лучше, чем путь длиной 1, образованный петлей). Поэтому я добавляю A обратно к результату. Теперь в нем есть записи 1 или 2, но имеет значение только тот факт, что они отличны от нуля.

Пример запуска вышеизложенного (я выбираю размер 6 на 6 для удобочитаемости вывода). Исходная матрица:

  (0, 3)    1
  (3, 2)    1
  (4, 3)    1
  (5, 1)    1
  (5, 3)    1
  (5, 4)    1
  (5, 5)    1 

Транзитивное закрытие:

  (0, 2)    1
  (0, 3)    2
  (3, 2)    2
  (4, 2)    1
  (4, 3)    2
  (5, 1)    2
  (5, 2)    1
  (5, 3)    2
  (5, 4)    2
  (5, 5)    1

Вы можете видеть, что это сработало правильно: добавлены записи (0, 2), (4, 2) и (5, 2), все они получены по пути (3, 2).

Кстати, у networkx тоже есть floyd_warshall, но его документация говорит

Этот алгоритм наиболее подходит для плотных графов. Время выполнения — O(n^3), а рабочее пространство — O(n^2), где n — количество узлов в G.

Выход снова плотный. У меня складывается впечатление, что этот алгоритм просто считается плотным по своей природе. Похоже, что all_pairs_shortest_path_length является своего рода алгоритмом Дейкстры.

Транзитивный и рефлексивный

Если вместо транзитивного замыкания (которое является наименьшим транзитивным отношением, содержащим заданное) вы хотите транзитивное и рефлексивное замыкание (наименьшее транзитивное и рефлексивное отношение, содержащее данное) , код упрощается, поскольку мы больше не беспокоимся о путях нулевой длины.

for row in paths:
  indices.extend(row[1])
  indptr.append(len(indices))
data = np.ones((len(indices),), dtype=np.uint8)
A_trans = sparse.csr_matrix((data, indices, indptr), shape=A.shape)

Транзитивный, рефлексивный и симметричный

Это означает нахождение наименьшего отношения эквивалентности, содержащего заданное. Эквивалентно делению вершин на компоненты связности. Для этого вам не нужно заходить на networkx, есть connected_components метода SciPy. Установите directed=False там. Пример:

import numpy as np
import scipy.stats as stats
import scipy.sparse as sparse
import itertools

A = sparse.random(20, 20, density=0.02, format='csr', data_rvs=stats.randint(1, 2).rvs).astype(np.uint8)
components = sparse.csgraph.connected_components(A, directed=False)
nonzeros = []
for k in range(components[0]):
  idx = np.where(components[1] == k)[0]
  nonzeros.extend(itertools.product(idx, idx))
  row = tuple(r for r, c in nonzeros)
  col = tuple(c for r, c in nonzeros)
  data = np.ones_like(row)
B = sparse.coo_matrix((data, (row, col)), shape=A.shape)

Вот как выглядит результат print(B.toarray()) для случайного примера, 20 на 20:

[[1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 1 0 0 0 0 0 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]
 [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 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 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0]
 [0 0 0 0 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 0 0 0 0 1 1 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]]