Part 2. Tarjan's algorithm for reducing a nonlinear system of equations to a computable sequence of substitutions

This is the second part of the article, devoted to the issue of applying Tarjan's algorithm to solving systems of equations. The first part considered the problem of finding a minimal set of equations, from which we will need only a few definitions.

Now we will study another problem. A system of nonlinear equations is given. It is necessary to divide the system into two subsystems – “equations” and “substitutions” using Tarjan's algorithm. The substitution subsystem must be formally computable. By substituting the equations of the second category into the first, we obtain a system of lower dimension, which will be easier to solve by numerical methods.

Formalization of the problem and an introductory example

Let's introduce a new definition: computable sequence of formal substitutions – each element of the sequence (equation) contains on the right side only already calculated unknowns, i.e. expressed through known variables.

Let's define a matrix $S_{n \times n}$ for the right part according to the principle from the previous article: $s_{ij} = 1$ If -variable is included in -th equation.

If the sequence is computable, then the matrix – lower triangular with zeros on the diagonal. Let's explain with an example:

$\begin{cases} x_1= A_1 \\ x_2 = \cos{x_1} + A_2 \\ x_3 = x_1^2 + x_1 + x_2 \\ x_4 = A_3 \end{cases}$

The corresponding matrix :

$S = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{pmatrix}$

The declared computability is that the known x_1 can be substituted into the equation that follows, finding x_2 famous x_1, x_2 are substituted into the equation for x_3 etc.

The substitutions can be expressed as an equation x = Sx + C . Solving the equation (e.g., by Gauss's method), we will establish whether the permutations are a computable sequence.

Algorithm

Reduce each equation of the system to the form $x_i = f(x_1, \dots, x_n, t)$ . should not be included in the right-hand side of the expression. Otherwise, assign this expression to the 'equations' group.
Find strongly connected components using Tarjan's algorithm.
Remove algebraic loops. Rerun Tarjan's algorithm.
Size Components $1 \times 1$ are substitutions, the rest are equations.
Order and obtain a formal computable sequence.

We see that, as in the previous article, the strength of Tarjan's algorithm is that it reduces the matrix to block lower triangular form. Only point 3 about removing algebraic loops is unclear, which we will discuss below.

Example

This example and the idea of the algorithm were taken from the article About detection substitutions in nonlinear algebraic equations with help of Tarjan's algorithm. Isakov A. A., Senichenkov Yu. B.

System:

$\left\{\begin{array}{l} x_{1}=-7x_{7}+x_{8}-2 \\ x_{2}=3x_{4} - 2x_{6}+1 \\ x_ {3}=\frac{x_{5}}{2}+8 \\ x_{4}=\frac{x_{6}}{x_{1}}+12 \\ x_{5}=(x_{ 7}+x_{8}+x_{3})^2-23 \\ x_{6}=x_{3}+x_{2}-3 \\ x_{7}=2 x_{3}+2 \ \ x_{8}=(x_{7}-2)^3+2 \end{array}\right.$

Visualization of the graph of relationships between equations and variables:

Fig. 1. Indices are shifted by 1 to the left.

On Removing Algebraic Loops

An algebraic loop is simple cycle in our graph, which is found by the NetworkX library method nx.simple_cycles. In fact, this is when, as a result of the substitution, it turns out that the variable depends on itself, violating the principle stated in paragraph 1 of the algorithm.

$1. \ x_3 \Longrightarrow x_5 \Longrightarrow x_7 \Longrightarrow x_8 \Longrightarrow x_3$ $2. \x_2 \Longrightarrow x_4 \Longrightarrow x_6 \Longrightarrow x_2$

Definition cyclic order (cycling order) := {number of incoming edges from the vertices of a strongly connected component (SCC); number of outgoing edges to the vertices of the current SCC}.

It is necessary to select a vertex with the maximum cyclic order, assign it to the category of 'equations' and delete from the graph. CO(x_5) = 3 equal to the maximum, we will remove it.
If there are several vertices with the maximum cyclic order, then we will choose the one that was encountered first in Tarjan's algorithm. In this case, we will delete the vertex x_2 with order 2, which was encountered first.

Results of the algorithm

Let's move on to programming the algorithm. Data input:

M = np.array([
    [0, 0, 0, 0, 0, 0, 1, 1],
    [0, 0, 0, 1, 0, 1, 0, 0],
    [0, 0, 0, 0, 1, 0, 0, 0],
    [1, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 1, 0, 0, 0, 1, 1],
    [0, 1, 1, 0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 1, 0],
])

tsk = SubstitutionDetector(dim=8)
tsk.set_matrix(M)
tsk.break_loops()
tsk.get_answer()

Program output:

Для цикла [2, 4, 7, 6] будет удалена вершина: 4. Причина: Максимальный цикловый порядок
Для цикла [1, 3, 5] будет удалена вершина: 1. Причина: Максимальный цикловый порядок, первым встретился в алгоритме Тарьяна
Индексы уравнений: [4, 1]
Индексы(!) подстановок: [2, 6, 7, 0, 5, 3]
Индексы(!) подстановок, которые действительно есть в системе: [2, 6, 7, 0, 5, 3]

Therefore, the “equations” will be expressions on the left side of which is x_5, x_2 and the “substitutions” will be and remember that order matters.

Equations:

$\left\{\begin{array}{l}x_{2}=3x_{4} - 2x_{6}+1 \\x_{5}=(x_{7}+x_{8}+x_{3} )^2-23 \\\end{array}\right.$

Substitutions:

$\left\{\begin{array}{l}x_{3}=\frac{x_{5}}{2}+8 \\x_{7}=2 x_{3}+2 \\x_{8} =(x_{7}-2)^3+2 \\x_{1}=-7x_{7}+x_{8}-2 \\x_{6}=x_{3}+x_{2}-3 \\x_{4}=\frac{x_{6}}{x_{1}}+12 \\\end{array}\right.$

Note the order of the substitutions: now the computability condition is satisfied. After successive substitution of one equation into the next one, we obtain the system:

$\left\{\begin{array}{l} x_{3}=f(x_5) \\ x_{7}=f_1(x_5) \\ x_{8}=f_2(x_5) \\ x_{1}= f_3(x_5) \\ x_{6}=g(x_2, x_5) \\ x_{4}=g_1(x_2, x_5) \\ \end{array}\right.$

Now we can solve the first subsystem, which was classified as “equations”. When x_2, x_5 are found, then through them the remaining unknown ones can be found.

Code

Utility functions:

Hidden text

def tarjan(graph: dict[int, set[int]]):
    """
    Пример представления графа.
    graph = {
    0: {1},
    1: {3},
    2: {1, 4},
    3: {0},
    4: {2, 5},
    5: set()
    }
    :param graph:
    :return:
    """
    n = len(graph)
    sccs = []
    index = [0]
    indexes = [-1] * n
    lows = [float('Inf')] * n
    S = []

    def strongconnect(v):
        indexes[v] = index[0]
        lows[v] = index[0]
        index[0] += 1
        S.append(v)
        for chld in graph[v]:
            if indexes[chld] == -1:
                strongconnect(chld)
                lows[v] = min(lows[v], lows[chld])
            elif chld in S:
                lows[v] = min(lows[v], lows[chld])

        if lows[v] == indexes[v]:
            scc = [v]
            w = S.pop()
            while w != v:
                if w not in scc:
                    scc.append(w)
                w = S.pop()
            sccs.append(scc)

    for v in graph.keys():
        if indexes[v] == -1:
            strongconnect(v)

    return sccs

  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

SubstitutionDetector class:

Hidden text

import networkx as nx
import numpy as np
from tarjan.utils import tarjan, matr2dct


class SubstitutionDetector:
    def __init__(self, dim):
        self.translator = dict()
        self.inverted_translator = dict()
        self.graph = None
        self.scc = None
        self.matrix = np.zeros((dim, dim))
        self.equations_indices = []
        self.substitution_indices = []
        self.in_system = []

    def add_equation(self, varnum: int, dependnums: list[int]):
        for i in dependnums:
            self.matrix[varnum - 1][i - 1] = 1
            self.in_system.append(varnum - 1)

    def set_matrix(self, matr):
        self.matrix = matr
        for i in range(self.matrix.shape[0]):
            if any(elem != 0 for elem in self.matrix[i]):
                self.in_system.append(i)
        self.initialize()

    def cycling_order(self, vertex_idx):
        """
        max(исходящие ребра, выходящие ребра) в вершину той же компоненты сильной связности
        :param vertex_idx:
        :return:
        """
        group = self.get_component_by_vertex(vertex_idx)
        ones_in_row = 0
        ones_in_column = 0

        for i, r in enumerate(self.matrix[vertex_idx, :]):
            if r == 1 and i in group:
                ones_in_row += 1

        for i, r in enumerate(self.matrix[:, vertex_idx]):
            if r == 1 and i in group:
                ones_in_column += 1

        return max(ones_in_row, ones_in_column)

    def _strong_connected_components(self):
        return tarjan(matr2dct(self.matrix))

    def initialize(self):
        self.scc = self._strong_connected_components()
        self.graph = nx.DiGraph(self.matrix)

    def get_cycles(self):
        result = []
        nodes = set()
        lst = list(nx.simple_cycles(self.graph))
        lst.sort(key=len, reverse=True)
        for cycle in lst:
            if all((n not in nodes for n in cycle)):
                result.append(cycle)
                for n in cycle:
                    nodes.add(n)
        return result

    def get_component_by_vertex(self, vertex_idx):
        cycles = self.get_cycles()
        for cycle in cycles:
            if vertex_idx in cycle:
                return cycle

    def find_to_delete(self, cycle):
        dct = dict(zip(cycle, [self.cycling_order(c) for c in cycle]))
        inv_dct = {}
        for k, v in dct.items():
            inv_dct[v] = [key for key in dct if dct[key] == v]
        maxKey = max(inv_dct.keys())
        candidates = inv_dct[maxKey]
        if len(candidates) == 1:
            return candidates[0], False
        group = self.get_component_by_vertex(candidates[0])
        return min(candidates, key=group.index), True

    def break_loops(self):
        N = self.matrix.shape[0]
        self.translator = dict(zip(list(range(N)), list(range(N))))
        for cycle in self.get_cycles():
            num, reason = self.find_to_delete(cycle)
            self.equations_indices.append(num)
            print(f'Для цикла {cycle} будет удалена вершина: {num}. Причина: Максимальный цикловый порядок' + ', первым встретился в алгоритме Тарьяна' * reason)
            self.matrix = np.delete(self.matrix, num, axis=0)
            self.matrix = np.delete(self.matrix, num, axis=1)
            del self.translator[num]
            for k in self.translator:
                if k > num:
                    self.translator[k] -= 1
        self.inverted_translator = {v: k for k, v in self.translator.items()}
        self.initialize()

    def get_answer(self):
        print(f'Индексы уравнений: {self.equations_indices}')
        for sublist in self.scc:
            if len(sublist) > 1:
                print(f'{sublist} не подстановка, а видимо уравнение. Добавить в eq_indices?')
                continue
            self.substitution_indices.append(self.inverted_translator[sublist[0]])
        print(f'Индексы(!) подстановок: {self.substitution_indices}')
        print(f'Индексы(!) подстановок, которые действительно есть в системе: {[i for i in self.substitution_indices if i in self.in_system]}')

Additional example

Finally, I would like to draw attention to another way of initializing the system:

  tsk = SubstitutionDetector(8)
  tsk.add_equation(1, [2, 5])
  tsk.add_equation(3, [1, 6])
  tsk.add_equation(2, [4, 7])
  tsk.add_equation(4, [3, 8])
  tsk.add_equation(6, [7])

  tsk.initialize()
  tsk.break_loops()
  tsk.get_answer()

Для цикла [0, 1, 3, 2] будет удалена вершина: 0. Причина: Максимальный цикловый порядок, первым встретился в алгоритме Тарьяна.

Индексы уравнений: [0] Индексы(!) подстановок: [6, 5, 2, 7, 3, 1, 4] Индексы(!) подстановок, которые действительно есть в системе: [5, 2, 3, 1]

In this case, the system does not contain equations describing the variables. x_7, x_8, x_5 but the algorithm works by ordering and classifying existing equations.