Алгоритм Тарьяна для поиска минимального набора уравнений

Дана система, состоящая из большого количества уравнений (необязательно линейных), где вам необходимо найти всего лишь несколько переменных. Как это сделать эффективно? Какой минимальный набор уравнений вам потребуется? В этой статье мы обсудим графовое представление систем уравнений, применим алгоритм Тарьяна и формализуем процесс на Python.

Постановка задачи

Введём двудольный граф $G = (V_1 \cup V_2, \ E)$ . V_1 — множество вершин, представляющие уравнения, V_2 — вершины, представляющие переменные. Ребро $e \in E$ соединяет пар, если уравнение зависит от переменной .

И будем пользоваться матрицей смежности . По строкам — уравнения, по стобцам — переменные. Если $M_{ij} = 1$ , то -ое уравнение зависит от -переменной.

$S:\left\{\begin{array}{cccc} f_1\left(x_1, x_2, \ldots, x_n\right) & = & 0 & e_1 \\ f_2\left(x_1, x_2, \ldots, x_n\right) & = & 0 & e_2 \\ \vdots & & & \\ f_n\left(x_1, x_2, \ldots, x_n\right) & = & 0 & e_n \end{array},\right.$

Найти минимальную подсистему $v_1 \subset V_1$ уравнений, решив которую, можно найти переменные $v_2 = \{x_i, x_j, \dots\} \subset V_2$ (произвольное подмножество).

Формулировка алгоритма

Количество уравнений и переменных должно совпадать для применения данного алгоритма.

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

$M^{'} = \begin{pmatrix} 1 & a_{12} & \dots & a_{1n}\\ a_{21} & 1 & \dots & a_{2n} \\ \vdots & \vdots & \dots & a_{n - 1, n} \\ \dots & \dots & \dots & 1 \end{pmatrix}$

На самом деле, задача состоит в приведении матрицы к блочной нижнетреугольной (БНТ) форме. Это мы сделаем с помощью алгоритма Тарьяна, который найдёт компоненты сильной связности (КСС — это такой подграф, что для любой пары вершин А и Б, существует путь из А в Б и обратно — из Б в А.

Ключевой момент в том, что в качестве стартовой вершины будет выбрана переменная, для которой хотим найти минимальный набор уравнений (для простоты пока положим, что такая переменная одна). Тогда алгоритм вернёт КСС в том порядке, как они связаны с искомой переменной, неявно приводя таким образом матрицу смежности к БНТ форме.

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 точки зрения системы уравнений, добавляется вспомогательная переменная $r = f_r(x_1, \dots x_h)$ , выраженная некой функцией от интересующих нас переменных.

Пример

$\left\{\begin{array}{l}\dot{x_1}(t)=-x_1(t) \\\dot{x_2}(t)=x_1(t)-x_2(t) \\\dot{x_3}(t)=x_1(t) \\y_1(t)=3 x_2(t)+x_1(t) \\y_2(t)=2 x_3(t)\end{array} \right.$

Три переменных состояния () и две алгебраические переменные . Введём вспомогательных уравнения, показывающие связь между переменной и её производной $\dot{x_i}$ :

$\left\{\begin{array}{lll}\dot{x_1}+x_1 & =0 & e_0 \\\dot{x_2}-x_1+x_2 & =0 & e_1 \\\dot{x_3}-x_1 & =0 & e_2 \\y_1-3 x_2-x_1 & =0 & e_3 \\y_2-2 x_3 & =0 & e_4 \\x_1-f\left(\dot{x_1}\right) & =0 & e_5 \\x_2-g\left(\dot{x_2}\right) & =0 & e_6 \\x_3-h\left(\dot{x_3}\right) & =0 & e_7\end{array}\right.$

Матрица смежности:

$\begin{aligned}\begin{array}{l} e_0 \\ e_1 \\ e_2 \\ e_3 \\ e_4 \\ e_5 \\ e_6 \\ e_7 \end{array}\left(\begin{array}{llllllll} 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \end{array}\right) \\ & \end{aligned}$

Переставляем строки так, чтобы на главной диагонали были 1. Программа (см. ниже раздел Код, где представлена наивная реализация такого алгоритма) даст перестановку $\{0, 1, 7, 5, 6, 2, 3, 4\}$ :

$\begin{aligned}& \begin{array}{l}e_0 \\e_1 \\e_7 \\e_5 \\e_6 \\e_2 \\e_3 \\e_4\end{array}\left(\begin{array}{llllllll}1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & 0 & 1\end{array}\right) \\&\end{aligned}$

В системе восемь переменных, перечислим их здесь для удобства:

$\begin{array}{llllllllll} x_1 & x_2 & x_3 & \dot{x}_1 & \dot{x_2} & \dot{x}_3 & y_1 & y_2 \end{array}$

Предположим, что нужно найти минимальный набор уравнений, которые нужны для вывода переменной . Вывод программы: [[3, 0], [4, 1], [6]].

Рис. 1. Граф с выделенными КСС, необходимых для нахождения

Уравнения дадут нужный результат. Здесь мы подействовали на номера [[3, 0], [4, 1], [6]]перестановкой {0, 1, 7, 5, 6, 2, 3, 4}: 3 переходит в 5, 0 переходит в 0, 4 переходит в 6 и т.д.

Записав уравнения в заданном порядке получим матрицу сниженнной размерности в блочной нижнетреугольной форме: (переменные / столбцы в таком порядке: $x_1, \dot{x_1}, x_2, \dot{x_2}, y_1$ ).

$\begin{aligned}& \begin{array}{l}e_5 \\e_0 \\e_6 \\e_1 \\e_3\end{array}\end{aligned}\left(\begin{array}{lllll}1 & 1 & 0 & 0 & 0 \\1 & 1 & 0 & 0 & 0 \\0 & 0 & 1 & 1 & 0 \\0 & 1 & 1 & 1 & 0 \\0 & 1 & 0 & 1 & 1\end{array}\right)$

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

$\left\{\begin{array}{l}\dot{x_1}(t)=-x_1(t) \\\dot{x_2}(t)=x_1(t)-x_2(t) \\y_1(t)=3 x_2(t)+x_1(t)\end{array} \right.$

Пример 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()