Доклады Российской академии наук. Математика, информатика, процессы управления, 2023, T. 514, № 2, стр. 109-117

1. ВВЕДЕНИЕ

Начиная с SINDy [1] и PDE-FIND [2], поиск уравнений стал привлекательным инструментом машинного обучения для получения интерпретируемых и лаконичных моделей из физических данных. Было предложено множество расширений SINDy [3, 4]. Результатом работы каждого из описанных алгоритмов (кроме случая [4], который мы обсудим позже) является набор коэффициентов при слагаемых из относительно большого набора-библиотеки. Отсутствие произвола в итоговых коэффициентах справедливо для всех методов, основанных на регрессии. Это может быть разреженная регрессия [1], регрессия с помощью нейронной сети [5], символьная регрессия на графах с помощью генетического программирования [6], символьная регрессия на лесе (наборе графов слагаемых) с помощью более продвинутого генетического алгоритма [7].

Алгоритмы, основанные на SINDy [1], используют символьную (разреженную) регрессию с предопределенной библиотекой слагаемых. Каждое слагаемое является дифференциальным оператором, произведением дифференциальных операторов, примененных к входному полю данных или другой отдельно заданной нелинейностью. Ограничением и по сути недостатком этого подхода является предположение о том, что уравнение может быть выражено исключительно с использованием слагаемых, присутствующих в предопределенной библиотеке. Количество слагаемых в библиотеке может быть интерактивно сокращено для получения краткого выражения. Таким образом, требуется знать структуру нелинейности для каждого набора данных, что приводит нас к классической постановке обратной задачи, осложненной большим числом слагаемых в уравнении.

В случае, когда процесс, управляющий данными, неизвестен полностью или частично, постановка задачи, содержащая заранее определенный набор слагаемых или форму уравнения, становится непригодной. Например, нельзя восстановить волновое уравнение, ожидая уравнение параболического типа. В таких ситуациях используется вторая группа методов обнаружения уравнений, основанная на эволюционной оптимизации. Такие алгоритмы [6–9] менее требовательны к количеству априорных предположений и обеспечивают более широкое пространство поиска, позволяя получать различные локальные оптимумы в ходе многократных запусков.

В контексте поиска уравнений часто акцент делается на проблеме восстановления известного уравнения, которое интегрируется численно для получения данных. Близость восстановленного уравнения к известному ответу является основной метрикой качества в литературе [1, 4, 6, 7, 10]. Для ответа на этот вопрос в общем случае нам необходимо оценить, насколько хорошо полученное уравнение описывает данные статистически, что включает учет неопределенности.

В статье [4] используются разнообразные стратегии бутстрэппинга, включая применение бутстрэпа к данным и предварительно определенной библиотеке слагаемых. В результате, даже с использованием методов градиентного спуска, можно получить разные уравнения от запуска к запуску. Основная идея заключается в том, что чем чаще мы получаем одно и то же уравнение, используя случайную часть доступной информации, тем выше вероятность того, что оно точно описывает процесс. Однако фиксированная форма ${{u}_{t}} = F(x,{{u}_{x}},{{u}_{{xx}}},...)$ снижает вероятность нахождения кардинально новых уравнений, вместо этого приводя к изменениям в основном в коэффициентах.

В качестве альтернативы мы можем собрать все уравнения, полученные в результате эволюционной оптимизации, и оценить изменчивость коэффициентов каждого слагаемого и формы уравнения, т.е. провести анализ совместного распределения коэффициентов. Для построения сложных совместных распределений можно использовать байесовские сети, которые позволяют исследовать эти зависимости. Структурное обучение байесовских сетей выходит за рамки данной работы. Мы рекомендуем читателю обратиться к обзорной статье [11]. В регрессионном подходе байесовские сети также могут быть использованы для определения взаимовстречаемых слагаемых в уравнении, что поможет выявить устойчивые структуры.

Для понимания того, как неопределенность в коэффициентах или выбор устойчивой подструктуры на основе совместного распределения отражается на описании данных и решении, нам необходимо иметь возможность решать каждое уравнение, с которым мы сталкиваемся в процессе поиска. В классической постановке, как показано в [4], даже при бутстрэппинге данных и библиотеки мы получаем близкое уравнение, которое может отличаться только коэффициентами. Это упрощает решение уравнения и последующее сравнение с данными, поскольку мы можем настроить конкретный солвер для определенного типа уравнения и изменять только коэффициенты. Любые члены, не соответствующие заданному уравнению, при этом исключаются из рассмотрения вне зависимости от их значимости. К сожалению, такой подход не применим в общих случаях поиска, когда уравнение неизвестно.

При рассмотрении решателей дифференциальных уравнений, подходящих для случая, описанного выше, требуется довольно сильно ослабить требования на точность итогового решения в угоду возможности решать несколько классов уравнений без дополнительной настройки внутри одного запуска алгоритма. Алгоритмов, удовлетворяющих таким требованиям, немного. Решатель, основанный на деревьях решений и написанный на языке программирования Julia [12], а также решатели DeepXDE и TEDEouS, основанные на нейронных сетях и написанные на языке Python [13, 14].

Таким образом, для наиболее полного решения задачи устойчивого поиска дифференциального уравнения по данным требуется совместное использование трех различных алгоритмов – алгоритма поиска дифференциальных уравнений по данным (в нашем случае, EPDE [9]), алгоритма автоматизированного решения дифференциальных уравнений, основанного на нейронных сетях (TEDEouS [14]), и алгоритма определения структуры совместного распределения в виде байесовской сети (BAMT [15]). Использование данных алгоритмов в отличие от аналогичных решений, например [4], позволило отказаться от заранее заданной библиотеки слагаемых – единственным ограничением, в сущности, является максимальный порядок дифференциального уравнения, а так же не использовать заранее заданную форму уравнения. Реализация и численные эксперименты с использованием алгоритма устойчивого поиска описаны в последующем тексте.

2. ПОСТАНОВКА ЗАДАЧИ УСТОЙЧИВОГО ПОИСКА ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ПО ДАННЫМ

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

С точки зрения машинного обучения, задача поиска уравнения сводится к классическим проблемам символьной регрессии. Наша цель состоит в решении задачи оптимизации (1):

где $M(\theta )$ – параметризованная модель, данные D и $L(M(\theta ),D)$ – функция потерь, определяющая соответствие между данными и моделью. При этом, ввиду постановки устойчивой задачи, оценивается математическое ожидание $\mathbb{E}[ \cdot ]$ функции потерь, так как набор параметров $\theta $ является случайной величиной.

В задаче поиска дифференциальных уравнений по данным сами данные D обычно представляют собой наблюдения на сетке $\bar {x}$ (пространственной, временной) или с использованием сетки других параметров. Исследование различных параметризаций может быть выгодно в некоторых случаях. Например, в случае различных классов уравнений. Поэтому число параметров ${\text{|}}\theta {\text{|}}$ может изменяться в ходе оптимизации. При поиске дифференциальных уравнений так же возможно учесть физические граничные условия в процессе оптимизации. Например, можно зафиксировать граничные условия и включить их в функцию потерь $L(M(\theta ),D)$ в качестве параметра регуляризации.

Алгоритм поиска уравнения. Использование эволюционного алгоритма позволяет избежать ограничений, порождаемых заранее заданной библиотекой. В процессе эволюционной оптимизации слагаемые собираются из параметризованных семейств элементарных операций $F$, где каждое семейство содержит элементы ${{f}_{j}}(p_{1}^{{(j)}},p_{2}^{{(j)}},...,\bar {x})\, \in \,F$. Обычно рассматриваются только дифференциальные операторы ${{F}_{{diff}}} = \left\{ {\frac{{{{\partial }^{{{{p}_{{n + 1}}}}}}u}}{{\partial x_{1}^{{{{p}_{1}}}}...\partial x_{n}^{{{{p}_{n}}}}}}} \right\}$. Кроме того, мы можем включить обратную координатную функцию ${{F}_{{inv}}} = \frac{1}{{{{x}_{i}}}}$ и другие параметризованные дискретные поля ${{F}_{{func}}} = v({{p}_{1}},...,{{p}_{k}})$. Входными данными алгоритма является выбранный набор элементарных операций $F = \bigcup\limits_j^{} {{{F}_{j}}} $, из которых создаются слагаемые.

Мы используем общепринятые процедуры эволюционной оптимизации независимо от выбранных семейств токенов для оптимизации. А именно, мы обучаем модель в форме

В уравнении (2) $C = \{ {{c}_{1}},...{{c}_{L}}\} \in {{R}^{L}}$ представляет собой набор констант (амплитуд слагаемых), a_i(P_i, $\bar {x}) = \prod\limits_{j = 1}^{j = {{N}_{{tokens}}}} {{{f}_{j}}(p_{1}^{{(i)}},p_{2}^{{(i)}},...,\bar {x})} $ обозначает произведения токенов ${{f}_{j}} \in F$, выбранных из семейств токенов $F$, ${{P}_{i}} = \{ p_{1}^{{(i)}},p_{2}^{{(i)}},...\} $ представляет собой набор параметров для слагаемого ${{a}_{i}}$, и P представляет мультииндекс параметров. Количество параметров (и размер хромосомы) ${\text{|}}{\mathbf{P}}{\text{|}}$ может изменяться в процессе эволюционной оптимизации. В большинстве случаев рассматриваются только точки наблюдений, т.е. $\bar {x} \in X$. Однако сетки данных и модели могут быть разделены.

Параметры модели P определяются с использованием эволюционной оптимизации. Процесс оптимизации включает операцию скрещивания, которая обменивает слагаемые между особями, и мутацию, которая может включать изменения токенов или параметров. Этот процесс аналогичен стандартной символьной регрессии. Для завершения формулировки задачи эволюционной оптимизации необходимо определить функцию приспособленности.

Мы включаем разреженную регрессию в шаг вычисления функции пригодности. Мы используем градиентный спуск для поиска набора коэффициентов ${{C}_{{opt}}}$, которые минимизируют расхождение. А именно, мы решаем задачу оптимизации, показанную в уравнении

В уравнении (3) случайное слагаемое выбирается в качестве объясняемой переменной, а остальные слагаемые в данной особи используются как признак. С помощью такой формулировки задачи мы отходим от фиксированной формы уравнения u_t = = $F(u,{{u}_{x}},{{u}_{{xx}}}...)$ и допускаем произвольные уравнения. Кроме того, мы вводим ${{l}_{1}}$-регуляризацию с использованием ${{l}_{1}}$-нормы ${\text{||}}\, \cdot \,{\text{|}}{{{\text{|}}}_{1}}$. Слагаемые с коэффициентами ${{c}_{i}} \in c$, ${{c}_{i}} < \varepsilon $, где $\varepsilon $ – заданный порог, исключаются из модели как незначимые.

Автоматизированный инструмент для построения байесовской сети [15] используется для автоматического формирования совместного распределения. Впоследствии уравнения сэмплируются из байесовской сети, начиная с произвольной левой части с коэффициентом, равным единице. Каждое выбранное/сэмплированное уравнение затем решается c помощью нейросетевого метода [14] для получения неопределенностей в решении. Байесовская сеть служит в качестве устойчивой структуры мета-уравнения, из которой можно сэмплировать отдельные уравнения. Пример байесовской сети для дифференциального уравнения показан в Приложении A .

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

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

3. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

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

Системы обыкновенных дифференциальных уравнений. Мы вновь обращаемся к набору данных Lynx-Hare [16], который ранее использовался в [4] для восстановления модели Лотки-Вольтерры. Мы используем его, поскольку он прост и, следовательно, нагляден.

Модель Лотки-Вольтерра в [4] для этих данных¹¹ была получена с использованием фиксированного набора слагаемых для каждого уравнения. А именно, была использована модель уравнения (4).

При оптимизации используются упрощенные модели, уравнение (4), полученные путем случайных исключений слагаемых. Такой процесс в статье называется бутстраппинг библиотеки. Неопределенность слагаемых была получена с помощью 1000 моделей. Полученное распределение предельных коэффициентов имеет вид уравнение (5).

В уравнении (5) $\epsilon = O{{(10}^{{ - 5}}})$ и коэффициенты представлены средними значениями вместе с двухсторонними доверительными интервалами, построенными с использованием 95-го перцентиля стандартного распределения.

Отметим, что в отличие от [4] мы не нормализуем данные по дисперсии. Таким образом, коэффициенты слагаемых $uv$ должны быть умножены на 16.565 и 21.414 для первого и второго уравнений соответственно, чтобы достичь соответствия с результатами статьи. Уравнения в оригинальной статье были решены с использованием метода Рунге–Кутты с помощью решателя ode45 в Matlab с фиксированной формой уравнения (4).

В итоге, для обнаружения модели Лотки-Вольтерра из данных с использованием группы методов разреженной регрессии, необходимо знать форму уравнения, настроить решатель и, возможно, настроить метод бутстрэпа для получения различных результатов при каждом запуске. В качестве преимуществ мы получаем коэффициенты, близкие к теоретически полученным [17], которые показаны в уравнении (6) (остальные коэффициенты равны нулю).

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

Для оценки подхода, описываемого в статье, мы рассматриваем два случая обнаружения системы: (а) система, ограниченная уравнениями первого порядка, и (б), ограниченная уравнениями второго порядка. В отличие от заранее определенных библиотечных алгоритмов, каждое уравнение рассматривается независимо, что означает, что система состоит из нескольких взаимосвязанных уравнений, а не записывается в векторной форме.

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

После анализа байесовской сети мы можем определить одну или несколько компонент графа. В нашем случае была обнаружена одна компонента. Затем мы выбираем произвольный начальный узел, например, $\dot {u}$ для первого уравнения и ${\dot {v}}$ для второго уравнения. После выбора начального узла мы сэмплируем 30 уравнений для оценки неопределенности коэффициентов и решаем их для определения неопределенности в области данных.

Для случая (a) мы выделяем две устойчивые системы, представленные в уравнении (7). Для наглядности все коэффициенты представлены в виде доверительных интервалов как для маргинальных распределений. Однако мы отмечаем, что каждый коэффициент – многомерная случайная величина, для которой учитывается совместное распределение.

Первая система в уравнении (7) представляет собой систему Лотки-Вольтерра, которая является в основном правильной, с дополнительными шумовыми членами, определяющими дисперсию. Второе уравнение представляет собой процесс более высокого масштаба, а именно $\dot {u},\dot {v} = {{C}_{1}},{{C}_{2}}$, т.е. локально постоянное увеличение или уменьшение скорости роста популяции.

На рис. 2 показаны решения выбранных уравнений, полученные с использованием решателя на основе нейронных сетей.

В случае (б) размерность пространства поиска возрастает с учетом слагаемых второго порядка, но мы по-прежнему в основном получаем уравнения первого порядка. Итоговые уравнения для случая (б) показаны в (8).

Решение сэмплированных уравнений для случая (б) показано на рис. 3.

По сравнению с уравнением (7), коэффициенты (8) менее похожи на теоретические коэффициенты (6), но более похожи на полученные с помощью E-SINDy (уравнение (5)). Однако мы все еще можем извлечь устойчивую подмодель и более точно подобрать их при необходимости.

Мы собрали процентную ошибку для всех моделей по сравнению с результатами, полученными в работе [17], в табл. 1. Однако, поскольку это несинтетический случай, у нас нет “правильного ответа”, и метрики ошибки не дают никакой информации.

Метод, предложенный в данной работе, позволяет обнаруживать уравнения, не полагаясь на фиксированную структуру или заранее заданную форму (ср. [4] и (4)). В случае модели Лотки-Вольтерра входными данными являются значения $u$ и $v$ – популяций рыси и зайца соответственно, а также их численные производные до определенного порядка. Порядок и степень членов являются единственными ограничениями для алгоритма, которые могут быть легко расширены, как показано при переходе от случая (а) к случаю (б).

Представленный эксперимент представляет собой значительное расширение классического подхода с заранее определенными библиотеками. С помощью этого подхода мы приобретаем возможность обнаруживать уравнения, независимо от экспертного мнения. Следует отметить, что время, затраченное экспертом, значительно сокращается в пользу времени выполнения алгоритма. Для классического алгоритма [4] мы говорим о секундах, тогда как для предложенного подхода в сложных случаях речь идет о десятках минут.

Уравнения в частных производных. Рассмотрим процесс поиска уравнения в частных производных на примере невязкого уравнения Бюргерса, описанного также в работе [4]. Это уравнение представляет собой нелинейное однородное уравнение первого порядка:

Мы берем в качестве точных данных аналитическое решение, формирующее ударную волну:

где $\alpha = 0.5$, $A = 1000$. Данные без шума вычисляются в точках $({{x}_{i}},{{t}_{i}}) = ( - 4000 + i\Delta x,\;j\Delta t)$ c $\Delta x$ = = 31.25 и $\Delta t = 0.0157$ для $1 \geqslant i$, $j \geqslant 256$. Характеристики решения, используемого для определения невязкого уравнения Бюргерса, заключаются в том, что ударная волна формируется в момент времени t = 2 и движется вдоль линии x = = $500(t - 2)$.

Применяя предложенный алгоритм, было успешно получено уравнение, описывающее динамику рассматриваемой системы:

Отметим, что предельные коэффициенты перед членами уравнения (10) были оценены аналогичным образом. Кроме того, на рис. 4 представлены результаты решения сэмплированных уравнений на основе байесовской сети. Заметно, что применение алгоритма и в этом случае позволяет достичь хорошего соответствия между “средним” решением и исходными данными.

4. ЗАКЛЮЧЕНИЕ

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

$ \bullet $ Предопределенная библиотека слагаемых может быть заменена на алгоритм эволюционной оптимизации, который создает слагаемые из элементарных операций.

$ \bullet $ Оценка неопределенности должна учитывать совместные распределения, чтобы обеспечить более полное понимание взаимосвязей между слагаемыми, в частности их взаимовстречаемость.

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

Устойчивая постановка превращает задачу поиска уравнений в мощный инструмент, который применим не только к известным процессам и уравнениям, но также подходит для частично известных процессов (например, известное уравнение Навье–Стокса и неизвестное уравнение состояния) или даже полностью “черный ящик”. Это соответствует принципам машинного обучения, позволяя извлекать новые физические законы из данных и давая возможность компьютерам творчески предлагать новые уравнения и открывать действительно новые законы природы.

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

А. ПРИМЕР БАЙЕСОВСКОЙ СЕТИ