Дана система, состоящая из большого количества уравнений (необязательно линейных), где вам необходимо найти всего лишь несколько переменных. Как это сделать эффективно? Какой минимальный набор уравнений вам потребуется? В этой статье мы обсудим графовое представление систем уравнений, применим алгоритм Тарьяна и формализуем процесс на 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]].

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

Уравнения дадут нужный результат. Здесь мы подействовали на номера [[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()