Дана система, состоящая из большого количества уравнений (необязательно линейных), где вам необходимо найти всего лишь несколько переменных. Как это сделать эффективно? Какой минимальный набор уравнений вам потребуется? В этой статье мы обсудим графовое представление систем уравнений, применим алгоритм Тарьяна и формализуем процесс на 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):
Вывод программы: [[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()
Комментарии (6)
wataru
16.07.2024 10:46+1И вообще очень странно описано. Сначала у вас граф двудольный, с вершинам уравнениями и вершинами переменными. В этом графе вообще компонент сильной связности нет (кроме тривиальных). А потом вы отождествляете вершину-уравнение и вершину-переменную и уже в этом графе ищите КСС. Почему это работает вообще?
Кстати, ваш наивный алгоритм перестановки строк будет очень медленным. Там фактически O(n^n). Можно заметить, что эта задача эквивалентна задаче поиска максимального паросочетания в двудольном графе и решается за O(n^3) алгоритмом Куна.
y_mur
16.07.2024 10:46Не вникал в код и пока вообще не понял что здесь и о чем.
Просто возникла мысль, что если матрицу сделать в стиле С++, то перестановку строк можно будет делать очень быстро через перестановку ссылок на строки.
wataru
16.07.2024 10:46Мысль интересная, но на практике не важная.
Во-первых, перестановка строк - это доли процента от основной работы (перебор перестановок у автора. И даже если реализовать паросочетания, то все равно, перестановка O(n^2), когда как основная часть будет O(n^3)).
Во-вторых, если массив не расплющивать в один линейный, то любое обращение к нему будет через 2 адресации, вместо одной, что заметно медленнее. Так что не факт что ваша оптимизация вообще в плюс, хоть и минимальный.
y_mur
16.07.2024 10:46Во-вторых, если массив не расплющивать в один линейный, то любое обращение к нему будет через 2 адресации, вместо одной, что заметно медленнее. Так что не факт что ваша оптимизация вообще в плюс, хоть и минимальный.
В моих тестах на алгоритме Гаусса (решение системы линейных уравнений) сравнение двойной адресации и расплющивания в линейный массив дают:
на clang: производительность одинаково высокая (как в варианте с двойной адресацией на msvc).
на msvc: "расплющивание" мешает компилятору применять SSE / AVX векторизацию из-за чего вариант с двойной адресацией работает в 2-3 раза быстрее.
И это тесты без использования выигрыша от возможности быстрой перестановки строк.
Пруфы выложу позже, когда наиграюсь с этими эффектами и оформлю их в статью.
y_mur
16.07.2024 10:46Дана система, состоящая из большого количества уравнений (необязательно линейных), где вам необходимо найти всего лишь несколько переменных.
А можно привести примеры, где такое можно встретить на практике? Зачем оно может понадобиться?
wataru
Вот что непонятно, а что делать с тождественными уравнениями? Это же нигде не проверяется.
В простейшем случае линейных уравнений, может быть так, что одно уравнение есть линейная комбинация других. Если случайно взять их все, то одно из них будет бесполезным. И может так получится, что выбранный набор недостаточен для нахождения нужной переменной.