Graph automorphisms, canonical labeling, nauty

Раскрываем тему канонической нумерации графов. О деревьях была отдельная заметка, теперь переключим внимание на общий случай.

Допустим, мы храним и поддерживаем базу данных молекул. Или путей метаболизма в биологических клетках, или вовсе абстрактных графов, для какой-либо задачи по комбинаторике. В любом случае, речь идёт о графах с дополнительными свойствами.

При помещении в базу очередной записи необходимо проверить, что аналогичной записи в базе ещё нет. Выражаясь математически, нужно убедиться в отсутствии графов, изоморфных добавляемому. Наиболее грубый (и наименее быстрый) способ выглядит так: мы перебираем всю базу и осуществляем проверку изоморфизма каждого графа из базы с новым образцом.

Проверка изоморфизма двух графов — непростая задача. Она очевидно принадлежит к классу NP, но принадлежность её к классам P или NP-полных до сих пор под вопросом. (В комментариях дана ссылка на алгоритм, возможно решающий эту задачу за полиномиальное время.) Пока что для этой задачи придуман собственный класс сложности в NP, под названием GI. Но пусть даже у нас будет эффективная процедура для проверки изоморфизма графов; в любом случае, если записей в базе — миллион, то на проверку присутствия нового образца уйдут многие минуты.

Возможно, перед лицом этой опасности и были изобретены канонические коды графов. Канонический код — это строковое представление графа, не зависящее от порядка нумерации вершин; своего рода хеш-код без промахов. Имея такой код, никогда не потребуется проверять непосредственно изоморфизм графов, достаточно сравнивать две строки. Коды одинаковы тогда и только тогда, когда графы изоморфны.

Таким образом, перед добавлением графа в базу мы вычислим его канонический код, и затем осуществим поиск этого кода (строки) среди имеющихся. Поиск строк, в отличие от изоморфных графов, осуществляется быстро; проблема решена.

Следует отметить, что канонические коды — не единственное средство. Существует немало кодов (см. комментарии), которые вычисляются за полиномиальное время и не являются «каноническими» в том смысле, что могут совпадать у неизоморфных графов. Однако, если коды разные, то графы неизоморфны. Пользуясь этим, можно осуществить отсеивание (screening) большей части неизоморфных графов, и проверять изоморфизм только для оставшихся.

Очевидно, что задача вычисления канонического кода не проще, чем задача проверки изоморфизма; но она интереснее, и канонические коды удобнее в использовании: на один граф вызывается ровно одна специальная процедура.

Любой канонический код состоит из (i) определённой схемы кодирования графа и (ii) процедуры переупорядочения вершин. Нас не интересует схема кодирования; это может быть обход в глубину, или в ширину, или развёрнутая в строку матрица смежности, или просто списки вершин и рёбер... как бы то ни было, кодирование графа подразумевает какой-то порядок его вершин. Наша задача — упорядочить вершины так, чтобы результат не зависел от изначального порядка. Это называется каноническая нумерация (canonical numbering, canonical labeling).

Похоже, что пионерами этой задачи являются Арлазаров, Зуев, Усков и Фараджев, написавшие в 1973 году статью «Алгоритм приведения неориентированных графов к каноническому виду». Насколько большую роль сыграли их идеи — это отдельная тема, выходящая за рамки данной заметки (более того, в оффлайн); мы же останемся в реалиях нынешней сетевой эпохи.

Реалии, однако, мало изменились за последние 15 лет. Австралийский профессор Брендан МакКей (Brendan McKay) в 1980-х годах написал статью «Practical Graph Isomorphism» и соответствующий программный пакет под названием nauty. Долгое время никому не удалось превзойти nauty по быстродействию; хотя для некоторых частных случаев более эффективные алгоритмы всё же появляются, например вот, вот. Особого упоминания заслуживают библиотеки bliss (2007) и saucy2 (2008), которые на больших разреженных графах уверенно обгоняют nauty. Но эти алгоритмы, фактически, не более чем nauty с модификациями.

Печально другое: мало кто из программистов сумел хотя бы подойти близко к пониманию того, как работает nauty. Если мы посмотрим на современные программные пакеты, то окажется, что они либо напрямую используют эту библиотеку: Magma, MuPAD, GRAPE, Planarity, Cyberaide, igraph (использует bliss), либо используют более примитивные алгоритмы: OpenBabel, Marvin. В любом случае, об осмыслении nauty речи не идёт. Приятным исключением является InChI, в которой присутствует алгоритм, явно сделанный с оглядкой на nauty. Настоящая популярность к алгоритму МакКея придёт тогда, когда в каждом пакете он будет реализован по-своему.


Перестановки и упорядоченные разбиения

Говоря здесь о разбиениях, мы подразумеваем упорядоченные разбиения множества вершин графа, т.е. представления этого множества вершин виде упорядоченного набора его попарно непересекающихся непустых подмножеств. [a b c | d e] и [a c b | e d] — одно и то же разбиение, но [d e | a b c] от него отличается. Подмножества [a b c] и [d e] называются ячейками (cells).

Если каждая ячейка разбиения π₁ является подмножеством какой-либо ячейки разбиения π₂, то π₁ называется подразбиением π₂, а π₂ — надразбиением π₁. Также говорят, что π₁ мельче π₂, а π₂ крупнее π₁. Мы подразумеваем, что при подразбиении сохраняется порядок ячеек. Например, [a b | c | d e] является подразбиением [a b c | d e] и [a b | c d e], но не [c | a b d e].

Разбиение, в котором все подмножества состоят из одного элемента, называется дискретным (discrete partition). Дискретное разбиение очевидным образом образует перестановку (permutation) вершин. Например, дискретное разбиение [b | c | a | d | e] образует перестановку (b, c, a, d, e).

Разбиение, состоящее из одной ячейки, называется элементарным (unit partition). Любое разбиение является подразбиением элементарного.

Покажем процедуру перебора всех перестановок вершин, основанную на переборе разбиений. При своей крайней примитивности, она является базой всего алгоритма канонической нумерации. Для получения всех перестановок в процедуру следует передать элементарное разбиение.

function AllPermutations(π = [V₀|...|V_k])
   if π — дискретное разбиение:
      yield π
   else:
      выбрать V_i ∈ π : |V_i| ≠ 1
      for v ∈ V_i:
         AllPermutations([V₀|...|{v}|V_i\v|...|V_k])

На этом месте у читателя может возникнуть вопрос: зачем нужны разбиения, если наша цель — перебор перестановок? Да, для перебора всех перестановок не нужно вводить понятие разбиений. Но на самом деле, все перестановки нам не нужны, а нужна одна перестановка — каноническая, которую мы выбираем как «лучшую» среди какого-то класса перестановок (не обязательно всех). Очевидно, что чем меньше этот класс, тем лучше, но есть условие: состав этого класса перестановок не должен зависеть от исходной нумерации вершин. Именно здесь нам поможет идея разбиений: мы будем не отбрасывать перестановки, а уточнять разбиения, участвующие в переборе.


Годные разбиения

Мы называем разбиение π годным (equitable), если для него выполняется следующее свойство: для любых двух ячеек V₁, V₂ ∈ π (не обязательно различных) и для любых v₁, v₂ ∈ V₁ выполняется равенство d(v₁, V₂)=d(v₂, V₂), где d(v,V) — количество вершин в V, смежных с v. Очевидно, что все дискретные разбиения являются годнымы.

В работе МакКея доказано, что для каждого разбиения π существует единственное его наикрупнейшее подразбиение, являющееся годным. Процедура уточнения (refine) разбиения находит это годное подразбиение. В упрощённом виде она выглядит так:

function Refine(π = [V₀|...|V_k])
   for Vi ∈ π :
      for Vj ∈ π:
         подразбить V_j, сгруппировав элементы v ∈ V_j по d(v, V_i) и отсортировав по возрастанию d(v, V_i)
         заменить V_j в π на результат этого подразбиения, удлиняя цикл по V_i, но не удлиняя цикл по V_j
   return π

Следовательно, чтобы перебрать все годные подразбиения и все происходящие от них перестановки, нам нужно всего лишь вставить вызов процедуры Refine в перебор:

function AllEquitablePermutations(π = [V₀|...|V_k])
   Refine(π)
   if π — дискретное разбиение:
      yield π
   else:
      выбрать V_i ∈ π : |V_i| ≠ 1
      for v ∈ V_i:
         AllPermutations([V₀|...|{v}|V_i\v|...|V_k])

Процесс перебора перестановок можно изобразить как дерево разбиений, которое обходится слева направо. Среди листьев первым при обходе встречается самый левый. Обозначим его ξ.


(Правая половина дерева показана не целиком. Вызовы Refine, которые не уточняют подразбиение, также опущены.)


Сравнения перестановок и автоморфизмы

Наша задача — найти каноническую перестановку множества вершин исходного графа. Для этого нужно обзавестить функцией сравнения перестановок cmp, которая бы определяла, которая из двух «лучше» (или они одинаковые). Функция должна задавать транзитивное отношение, а её результат не должен зависеть от первоначальной нумерации вершин графа. Пример такой функции: по исходному графу G и перестановкам π₁ и π₂ построить графы G₁ и G₂, вычислить их матрицы смежности, и в качестве cmp(π₁, π₂) вернуть результат лексикографического сравнения матриц, развёрнутых в строки. Таким образом, канонической перестановкой будет считаться та, для которой матрица смежности наибольшая; если таких перестановок две и более, то они считаются неразличимыми и в качестве канонической можно взять любую из них.

Какой смысл для нас имеют «неразличимые» перестановки? Пусть cmp(π₁, π₂)=0. Несложно доказать, что в этом случае перестановка γ=π₂^-1∘π₁ (композиция π₁ и перестановки, обратной π₂) является автоморфизмом вершин исходного графа.

Использование автоморфизмов позволяет существенно (иногда — на порядки) сократить перебор при нахождении канонической перестановки. Поэтому в nauty и производных алгоритмах обязательно осуществляется нахождение всей группы автоморфизмов графа, тогда как нахождение собственно канонической нумерации сделано опцией.

Есть следующий способ найти генераторы группы автоморфизмов графа: при обработке первого встретившегося дискретного разбиения ξ оно сохраняется. Каждое последующее дискретное разибение π сравнивается с ξ: если перестановка ξ^-1∘π является автоморфизмом графа, то эта перестановка записывается в группу генераторов.


Сокращение перебора 1: автоморфизм ξ^-1∘π

После обнаружения автоморфизма ξ^-1∘π, обход можно «поднять» по дереву до общего предка π и ξ, поскольку поддерево этого предка, в котором находится π, «изоморфно» поддереву, в котором находится ξ, и не содержит новой информации. Строгое доказательство приведено в работе МакКея.


Сокращение перебора 2: автоморфизм ρ^-1∘π

Как было сказано, поиск канонической нумерации — только опция в алгоритме нахождения автоморфизмов. Когда эта опция включена, мы храним ещё одну перестановку ρ, которая является «лучшей» на данный момент (best so far, BSF). При обнаружении первой перестановки ξ мы копируем её в ρ. Каждую последующую перестановку π мы сравниваем с ρ. Если cmp(π, ρ) > 0, то мы копируем π в ρ. Если cmp(π, ρ) < 0, то ничего не делаем. Если же cmp(π, ρ) = 0, то ρ^-1∘π является новым автоморфизмом, и обход можно «поднять» по дереву до общего предка π и ρ.


Сокращение перебора 3: орбиты

Орбита вершины v — это подмножество вершин, в которые может перейти v в результате действия автоморфизмов. Например, на следующем рисунке показаны орбиты вершин графа:


В начале работы nauty каждая вершина принадлежит своей собственной орбите. По мере открывания новых автоморфизмов орбиты увеличиваются, а их количество — сокращается. Например, если найден автоморфизм (a, b, c, d, e) → (b, a, c, e, d), то список орбит из [a][b][c][d][e] становится [a b] [c] [d e]. Если после этого найден автоморфизм (a, b, c, d, e) → (a, c, b, d, e), то список орбит становится [a b c] [d e]

Рассмотрим «линию предков» первого дискретного разбиения ξ. Предположим, мы в процессе обхода дерева вернулись к какому-либо из предков ξ (назовём его ν), чтобы продолжить обход других потомков ν, перебирая другие элементы ячейки V_i. Если рассматриваемый элемент v_k ∈ V_i лежит на одной орбите с уже рассмотренным элементом v_j ∈ V_i, то существует автоморфизм, открытый в процессе обхода более ранних узлов, который переводит v_j в v_k. В работе МакКея доказано, что в этом случае нет смысла рассматривать всё поддерево, произведённое отделением v_k: все листья этого поддерева будут изоморфны ранее открытым листьям, и все автоморфизмы, которые будут найдены при обходе этого поддерева, будут повторять предыдущие.

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


Для всех графов, однако, верно то, что разбиение по орбитам является подразбиением уточнения тривиального разбиения.


Сокращение перебора 4: fix и mcr

С узлами, которые не являются предками ξ, предыдущее сокращение перебора делать нельзя. Для них есть менее сильный, но тоже действенный критерий отсечения потомков. Для применения этого критерия необходимо сохранять информацию о каждом автоморфизме в отдельности (вместо множества орбит, которое накапливает информацию о всех автоморфизмах в совокупности). Для каждого найденного автоморфизма γ запоминаются два множества: fix(γ) и mcr(γ). fix(γ) — это множество вершин, которое при действиии автоморфизма остаются на своих местах. mcr(γ) — это множество вершин, которые являются минимальными по номеру в своих орбитах. Имеются в виду орбиты, составленные по одному автоморфизму γ, а не те, которые используются в предыдущем способе отсечения.

Рассмотрим множество V_i, элементы которого мы перебираем для продолжения обхода дерева из какого-то узла, глубины l. В процессе пути от корня дерева до этого узла l вершин были «зафиксированы», когда мы отделяли их от их ячеек. Назовём множество этих вершин fixed.

В работе МаКея доказано следующее: если fixed ⊂ fix(γ), то можно присвоить в V_i пересечение V_i ∩ mcr(γ) без ущерба для алгоритма. Фактически, речь идёт об автоморфизмах, которые «не затрагивают» вершины, которые мы зафиксировали. Их орбиты можно эксплуатировать так же, как и в предыдущем разделе, что и делается путём пересечения V_i и mcr(γ). Означенная операция выполняется для всех γ, fix и mcr которых мы успели сохранить на данный момент. Следует отметить, что генераторов группы автоморфизмов графа может быть очень много, и сохранять все fix и mcr не всегда возможно. МакКей рекомендует ограничивать количество сохраняемых fix и mcr числом, пропорциональным количеству вершин в графе, с коэффициентом пропорции больше 1 (у него почему-то 50/32).


Алгоритм целиком

В коде имеются глобальные переменные:

• getcanon — опция получения канонической нумерации в ρ
• gca_first — уровень общего предка текущего узла и ξ
• gca_canon — уровень общего предка текущего узла и ρ
• orbits — массив индексов орбит; содержит для каждой вершины номер минимальной вершины в её орбите
• Ψ — списов посчитанных множеств fix и mcr для найденных автоморфизмов

function AutomorphismSearch(π)
   FirstNode(1, π)

function FirstNode(level, π = [V₀|...|V_k])
   π ← Refine(π)
   if π — дискретное:
      ξ ← π
      gca_first ← level
      if getcanon:
         ρ ← π
         gca_canon ← level
      return level-1
   Vi ← первая ячейка π, в которой больше 1 вершины
   for v ∈ Vi
      if orbits[v] ≠ v:
         continue (сокращение 3)
      π'←[V₀|...|{v}|V_i\v|...|V_k]
      fixed ← fixed ∪ v
      if v — первая вершина в V_i:
         rtnlevel ← FirstNode(level+1, π')
         gca_first ← level
      else:
         rtnlevel ← OtherNode(level+1, π')
      fixed ← fixed \ v
      if rtnlevel < level:
         return rtnlevel
      if level < gca_canon:
         gca_canon ← level
   return level-1

function OtherNode(level, π = [V₀|...|V_k])
   π ← Refine(π)
   if π — дискретное:
      if cmp(π,ξ)=0:
         AddFixMcr(ξ^-1∘π)
         JoinOrbits(ξ^-1∘π)
         return gca_first (сокращение 1)
      if getcanon:
         comp ← cmp(π, ρ)
         if comp=0:
            AddFixMcr(ρ^-1∘π)
            JoinOrbits(ρ^-1∘π)
            return gca_canon (сокращение 2)
         if comp>0:
            ρ ← π
      return level-1
   Vi ← первая ячейка π, в которой больше 1 вершины
   for (fix, mcr) ∈ Ψ:
      if fixed ⊂ fix:
         V_i ← V_i ∩ mcr (сокращение 4)
   for v ∈ Vi:
      π'←[V₀|...|{v}|V_i\v|...|V_k]
      fixed ← fixed ∪ v
      rtnlevel ← OtherNode(level+1, π')
      fixed ← fixed \ v
      if rtnlevel < level:
         return rtnlevel
      if level < gca_canon:
         gca_canon ← level
   return level-1


Детали реализации

Авторская имплементация nauty содержит трюки для увеличения быстродействия и уменьшения расхода памяти:

• Для хранения графа и множеств используются битовые массивы
• Вся цепочка подразбиений от корня дерева к текущему узлу кодируется всего двумя массивами (lab и ptn). При возвращении к родительскому узлу разбиение восстанавливается функцией recover

Подробности можно узнать в nauty User’s Guide и в коде самой nauty.


Применение

Прелесть изложенного алгоритма в том, что в нём не много завязок на то, что исходное множество вершин — это вершины простого графа. Фактически, можно находить автоморфизмы и каноническую нумерацию для чего угодно, например:

• Для ориентированных графов. Никаких изменений в алгоритме не требуется. МакКей в статье упоминает ориентированные графы, и в nauty есть соответствующая опция digraph
• Для графов с раскрашенными вершинами. Всё, что требуется — начинать не с элементарного разбиения, а сгруппировать вершины по цветам. Эта опция также описана у МакКея и имеется в nauty
• Для графов с раскрашенными рёбрами. Всё, что требуется — изменить функцию cmp так, чтобы она учитывала цвета рёбер. Также, можно более эффективно написать процедуру Refine (уточнение разбиения)
• Для молекул со стереохимическими свойствами. Опять же, модификации требуются только в функции cmp

Модификации алгоритма

Не было упомянуто об огромном количестве дополнительных оптимизаций, имеющихся в коде nauty. МакКей признаёт, что у него самого не все они задокументированы. Вот только некоторые:

• «Даровые» автоморфизмы (cheap automorphisms, implicit automorphisms), для обнаружения которых не требуется доходить до дискретного разбиения. Работают в nauty только для неориентированных графов (цвета вершин допускаются)
• Проверка автоморфизма cmp(π,ξ) может делаться более эффективно без вызова cmp, с учётом того, что ξ никогда не меняется и нас интересует только то, есть ли изоморфизм
• Псевдо-канонические коды разбиений, благодаря которым можно выявить отсутствие автоморфизма, не проверяя его явно. Они же используются для предварительного определения знака функции cmp
• При выборе ячейки V_i для подразбиения можно применять различные эвристики, с целью сокращения перебора