TensorFlow. Решение систем нелинейных уравнений

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

Идея решения

Для получения решения необходимо выполнить следующие действия:

  1. Определить область поиска решений/сетку начальных условий
  2. Построение графа, реализующего систему
  3. Выбрать начальные условия
  4. Решить оптимизационную задачу
  5. Перейти на п.3, если не найдены все решения/не перебраны все выбранные начальные условия
  6. Объединить эквивалентные решения
  7.  Profit

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

Пример решения задачи

Для большей ясности изложения решения, рассмотрим его на примере следующей системы:

\[
\begin{cases}
x^2 — 2y^2 — xy + 2x — y + 1 = 0 \\
2x^2 — y^2 +xy + 3y — 5 = 0
\end{cases}
\]

Инициализация переменных (начальных условий)

Импортируем пакет tensorflow

import tensorflow as tf

Создаем интерактивную сессию

sess = tf.InteractiveSession()

Далее объявляем инициализаторы. Одним из простейших вариантов является использование случайного равномерного распределения для инициализации переменных.

x_init = tf.initializers.random_uniform(minval=-3.0, maxval=3.0)
y_init = tf.initializers.random_uniform(minval=-3.0, maxval=3.0)

Теперь создаем переменные и передаем им объект инициализатора

x = tf.get_variable('x', shape=1 ,initializer=x_init)
y = tf.get_variable('y', shape=1 ,initializer=y_init)

Запустим для демонстрации 5 раз и выведем начальные значения переменных. Для запуска выполняем инициализацию всех переменных sess.run(tf.global_variables_initializer()) и вычисление начальных значений sess.run([x,y])

for i in range(5):   
    sess.run(tf.global_variables_initializer())
    sess.run([x,y])
    print(x.eval(), y.eval())

# Output:
# [-0.65766764] [1.8217659]
# [2.3475676] [-2.5066323]
# [-0.06060576] [-0.32884097]
# [2.2769918] [1.5647278]
# [1.4144545] [-0.12284446]

Заметим, что вывод, скорее всего, будет другой, так как значения инициализируются случайным образом.

Построение графа системы нелинейных уравнений

Следующим шагом является создание графа, реализующего левую часть системы (в правой части должны быть 0). Граф формируется на базе созданных ранее переменных \(x\) и  \(y\). При необходимости использования математических функций, их можно найти в пакете tf.math, например, квадратный корень.

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

eq1 = x*x - 2*y*y - x*y + 2*x - y + 1
eq2 = 2*x*x - y*y  + x*y + 3*y - 5

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

print(eq1.eval(), eq2.eval())

# Output
# [6.096011] [-1.5560191]

Поиск одиночного решения

Теперь можно перейти к процессу поиска решения системы уравнений. По определению, необходимо, чтобы значений выражений левых частей (в нашем случае eq1 и eq2) были равны 0, а на практике имели минимальное отклонение от 0.

Первым делом необходимо задать функцию потерь \(E\) и выбрать тип оптимизатора. Будем использовать среднеквадратичное отклонение в качестве функции потери и градиентный спуск в качестве оптимизатора.

loss = tf.losses.mean_squared_error([eq1, eq2],[(0.0,), (0.0,)]) # функция потерь
optimizer = tf.train.GradientDescentOptimizer(0.01).minimize(loss) # метод оптимизации

Зададим \(\epsilon\) меньше которого должна быть ошибка решения \(E < \epsilon\)

epsilon = 0.0001

Предварительно инициализируем переменными новыми значениями и выполним оптимизацию.

sess.run(tf.global_variables_initializer())
sess.run([x,y])
loss_v = sess.run([loss])
step=0
print("Step: " + str(step) + " Loss: " + str(loss_v))
print("x = %f, x = %f" % (x.eval(), y.eval()))
while loss_v[0] >= epsilon:
    step = step+1
    sess.run(optimizer)
    sess.run([x, y])
    loss_v = sess.run([loss])
    print("Step: " + str(step) + " Loss: " + str(loss_v))
    print("x = %f, x = %f" % (x.eval(), y.eval()))

# Output:
# ...
# Step: 121 Loss: [9.734182e-05]
# x = -1.499883, x = 0.508757

Поиск всех решений

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

Для удобства сделаем функцию поиска корня, возвращающую историю поиска:

def find_root():
    sess.run([x,y])
    loss_v = sess.run([loss])
    step=0
    x_root = []
    y_root = []
    while loss_v[0] >= epsilon:
        step = step+1
        sess.run(optimizer)
        sess.run([x, y])
        loss_v = sess.run([loss])        
        x_root.append(x.eval())
        y_root.append(y.eval())
    return (x_root, y_root)

Тогда для каждого корня можно получить траекторию:

sess.run(tf.global_variables_initializer())
(x_root, y_root) = find_root()

import matplotlib.pyplot as plt
%matplotlib inline  

plt.plot(x_root, y_root, 'r')
plt.grid()
plt.show()

Теперь запустим поиск множество (100) раз, сохраняя найденные корни для дальнейшего анализа и построим траектории для получаемых корней.

import numpy as np
x_roots = []
y_roots = []
for _ in range(100):
    sess.run(tf.global_variables_initializer())
    (x_root, y_root) = find_root()
    x_roots.append(x_root[-1])
    y_roots.append(y_root[-1])
    plt.plot(x_root, y_root, 'r')
plt.grid()
plt.show()

На данный момент мы имеем 100 пар корней, среди которых, естественно, есть повторяющиеся. Теперь необходимо найти повторяющиеся корни.

Предварительно выполним округление до десятых, т.к. шаг алгоритма оптимизации и целевая ошибка была 1 сотая. Кроме того, сразу сделаем NumPy массив из корней.

x_roots = np.around(x_roots, decimals=1)
y_roots = np.around(y_roots, decimals=1)
roots = np.hstack((x_roots, y_roots))

Следующей командой найдем уникальные корни, получив тем самым ответ:

np.unique(roots, axis=0)

# Output:
#array([[-1.7, -0.3],
#       [-1.5,  0.5],
#       [ 1. ,  1. ]], dtype=float32)

Таким образом у системы имеется три решения в рассматриваемой области:

\([-1.7, -0.3], [-1.6, 0.5], [1, 1]\)

Заключение

Приведенные способ хоть и решает поставленную задачу, он не является единственно верным и оптимальным. Так, например, для решения систем с одним корнем не обязательно делать перебор. Кроме того, случайная инициализация вектора начальных условий не всегда является оптимальной.

В целом подобным способом можно решать нелинейные системы алгебраических уравнений практически любой сложности с соответствующими оговорками.

Код и статья в формате Jupiter Notebook.