Алгоритм Тарьяна для поиска минимального набора уравнений
Дана система, состоящая из большого количества уравнений (необязательно линейных), где вам необходимо найти всего лишь несколько переменных. Как это сделать эффективно? Какой минимальный набор уравнений вам потребуется? В этой статье мы обсудим графовое представление систем уравнений, применим алгоритм Тарьяна и формализуем процесс на Python.
Постановка задачи
Введём двудольный граф .
— множество вершин, представляющие уравнения,
— вершины, представляющие переменные. Ребро
соединяет пар
, если уравнение
зависит от переменной
.
И будем пользоваться матрицей смежности . По строкам — уравнения, по стобцам — переменные. Если
, то
-ое уравнение зависит от
-переменной.
Найти минимальную подсистему уравнений, решив которую, можно найти переменные
(произвольное подмножество).
Формулировка алгоритма
Количество уравнений и переменных должно совпадать для применения данного алгоритма.
Находится такая перестановка рядов матрицы , что на главной диагонали будут стоять единицы. Другим словами, это такая перенумеровка, что каждая переменная входит в уравнение с тем же номером.
На самом деле, задача состоит в приведении матрицы к блочной нижнетреугольной (БНТ) форме. Это мы сделаем с помощью алгоритма Тарьяна, который найдёт компоненты сильной связности (КСС) в ориентированном графе. КСС — это такой подграф, что для любой пары вершин А и Б, существует путь из А в Б и обратно — из Б в А.
Ключевой момент в том, что в качестве стартовой вершины будет выбрана переменная, для которой хотим найти минимальный набор уравнений (для простоты пока положим, что такая переменная одна). Тогда алгоритм вернёт КСС в том порядке, как они связаны с искомой переменной, неявно приводя таким образом матрицу смежности к БНТ форме.
Hidden text
def tarjan_algorithm_with_start(graph: dict[int, set[int]], start_node: int):
def dfs(v):
nonlocal index
indices[v] = index
lowlink[v] = index
index += 1
stack.append(v)
on_stack[v] = True
for neighbor in graph[v]:
if indices[neighbor] == -1:
dfs(neighbor)
lowlink[v] = min(lowlink[v], lowlink[neighbor])
elif on_stack[neighbor]:
lowlink[v] = min(lowlink[v], indices[neighbor])
if indices[v] == lowlink[v]:
scc = []
while True:
node = stack.pop()
on_stack[node] = False
scc.append(node)
if node == v:
break
strongly_connected_components.append(scc)
index = 0
stack = []
indices = [-1] * len(graph)
lowlink = [-1] * len(graph)
on_stack = [False] * len(graph)
strongly_connected_components = []
dfs(start_node)
return strongly_connected_components
Если искомых переменных несколько, то добавим новую вершину (назовём её корневой), которая будет соединена выходящими ребрами с исследуемыми переменными. Тогда в качестве стартовой вершины следует взять корневую
.
C точки зрения системы уравнений, добавляется вспомогательная переменная , выраженная некой функцией от интересующих нас переменных.
Пример
Три переменных состояния () и две алгебраические переменные
. Введём вспомогательных уравнения, показывающие связь между переменной
и её производной
:
Матрица смежности:
Переставляем строки так, чтобы на главной диагонали были 1. Программа (см. ниже раздел Код, где представлена наивная реализация такого алгоритма) даст перестановку :
В системе восемь переменных, перечислим их здесь для удобства:
Предположим, что нужно найти минимальный набор уравнений, которые нужны для вывода переменной . Вывод программы:
[[3, 0], [4, 1], [6]]
.
Уравнения дадут нужный результат. Здесь мы подействовали на номера
[[3, 0], [4, 1], [6]]
перестановкой {0, 1, 7, 5, 6, 2, 3, 4}: 3 переходит в 5, 0 переходит в 0, 4 переходит в 6 и т.д.
Записав уравнения в заданном порядке получим матрицу сниженнной размерности в блочной нижнетреугольной форме: (переменные / столбцы в таком порядке: ).
Получили блочную нижнетреугольную форму. Вспомним, что уравнения фиктивные, отражающие связь между функцией и её производной.
Пример 2
Найдём уравнения, необходимые для вычисления и
. Введём корневую вершину
(на рисунке обозначена 8):
Рис. 2. Граф с добавленной корневой вершиной
Вывод программы: [[3, 0], [4, 1], [6], [5], [2], [7], [8]].
Потребуются все уравнения в системе, чтобы найти . Однако алгоритм привёл матрицу в удобную БНТ форму.
При подготовке данного материала использовался пример и концепция из статьи Minimal Equation Sets for Output Computation in Object-Oriented Models. Vincenzo Manzoni Francesco Casella Dipartimento di Elettronica e Informazione, Politecnico di Milano.
Код
Hidden text
import itertools
import numpy as np
from numpy.typing import NDArray
from tarjan.utils import tarjan_algorithm_with_start, matr2dct, colorGraph
class MinimalSubsetDetector:
def __init__(self, matrix: NDArray, omitOnes: bool = True):
assert matrix.shape[0] == matrix.shape[1], 'Матрица должна быть квадратной'
permuted = MinimalSubsetDetector.permuteMatrix(matrix)
if permuted is None:
raise ValueError('У переданной матрицы не существует такой перестановки строк, что на главной диагонали стоят единицы.')
self.matrix = permuted
self.size = self.matrix.shape[0]
self.hasRoot = False
if omitOnes:
for i in range(self.size):
self.matrix[i][i] = 0
@staticmethod
def permuteMatrix(matrix) -> NDArray | None:
"""
Находит такую перестановку строк, что на главной диагонали должны стоять только 1.
O(k^N), где k - примерное кол-во 1 в строке.
:param matrix:
:return:
"""
result = []
for row in matrix:
indices = np.where(row == 1)[0]
result.append(indices.tolist())
result_dict = {}
for i, sublist in enumerate(result):
for element in sublist:
if element not in result_dict:
result_dict[element] = [i]
else:
result_dict[element].append(i)
result_list = [result_dict.get(i, []) for i in range(len(result))]
combinations = itertools.product(*result_list)
valid_perm = next(comb for comb in combinations if len(set(comb)) == len(comb))
if valid_perm is None:
print('Такой перестановки не существует.')
return
print(f'Необходимая перестановка: {valid_perm}.')
return matrix[list(valid_perm)]
@classmethod
def initfromTeX(cls, dimension: int, texcode: str, omitOnes: bool = True):
"""
:param omitOnes: опустить единицы на главной диагонали матрицы (убрать петли).
:param dimension: размерность матрицы
:param texcode: код для матрицы
:return:
"""
tex_matrix = texcode.replace("\\", '&').replace("\n", "")
elements = tex_matrix.split("&")
int_elements = [int(element.strip()) for element in elements if len(element)]
return cls(np.array(int_elements).reshape(dimension, dimension), omitOnes=omitOnes)
def addRootNode(self, interested_nodes: list[int]):
"""
Добавление корневой вершины R, которая связана ребрами с interested_nodes. Применяется, когда
нужно узнать минимальный набор уравнений для нахождения группы неизвестных.
:param interested_nodes:
:return:
"""
self.hasRoot = True
row = [0] * self.size
for n in interested_nodes:
row[n] = 1
new_row = np.array(row)
extended_matrix = np.vstack([self.matrix, new_row])
new_column = np.array([0] * (self.size + 1))
new_column = new_column[:, np.newaxis]
self.matrix = np.hstack((extended_matrix, new_column))
def find(self, node: int = None):
if node is None and not self.hasRoot:
raise ValueError('Укажите индекс переменной, которую вы хотите найти.')
return tarjan_algorithm_with_start(matr2dct(self.matrix), start_node=self.size if self.hasRoot else node)
def color_SCC(self, scc):
return colorGraph(self.matrix, scc)
def test1():
detector = MinimalSubsetDetector.initfromTeX(
8,
r'1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\'
r'1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\'
r'0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\'
r'1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\'
r'0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\'
r'1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\'
r'1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\'
r'0 & 0 & 1 & 0 & 0 & 0 & 0 & 1'
)
# Для поиска минимального набора уравнений для переменных у1 и у2.
# detector.addRootNode([6, 7])
answers = detector.find(6) # найти для у1
print(answers)
detector.color_SCC(answers)
if __name__ == '__main__':
test1()
Утилитные функции:
Hidden text
import networkx as nx
import numpy as np
from matplotlib import pyplot as plt
def matr2dct(m):
gr = {}
num_nodes = m.shape[0]
for i in range(num_nodes):
neighbors = set(np.nonzero(m[i])[0])
gr[i] = neighbors
return gr
def colorGraph(matrix, scc):
G = nx.DiGraph(matrix)
pos = nx.spring_layout(G)
colors = ['red', 'green', 'blue', 'yellow', 'orange', 'violet', 'pink']
for component, color in zip(scc, colors):
nx.draw_networkx_nodes(G, pos, nodelist=component, node_color=color, node_size=500)
nx.draw_networkx_edges(G, pos, arrows=True)
nx.draw_networkx_labels(G, pos)
plt.axis('off')
# plt.savefig('graph.png')
plt.show()