Когда порядок производной становится дробным числом

В классическом математическом анализе, заложенном Ньютоном и Лейбницем, мы привыкли оперировать производными только целого порядка: первая производная (скорость), вторая производная (ускорение), третья производная (рывок) и так далее. Однако еще в 1695 году маркиз де Лопиталь задал Лейбницу в письме интригующий вопрос: «А что будет, если порядок производной n сделать равным 1/2?». Лейбниц ответил, что это приведет к удивительным парадоксам, из которых однажды будут извлечены великие следствия. Так зародилась совершенно невероятная ветвь математики — дробное исчисление (Fractional Calculus).

На протяжении веков этот раздел оставался чистой математической абстракцией, лишенной физического смысла. Но в конце 20-го века произошел настоящий взрыв: выяснилось, что дифференциальные уравнения с дробными производными феноменально точно описывают сложнейшие физические процессы с эффектами «долговременной памяти» и пространственной нелокальности. Моделирование вязкоупругих материалов (полимеров), аномальной диффузии в пористых средах, электрических фрактальных цепей и даже динамики финансовых рынков оказалось невозможным без использования дробных операторов.

Операторы Римана-Лиувилля и Капуто

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

Однако производная Римана-Лиувилля от константы не равна нулю, что вызывает огромные трудности при постановке начальных физических условий (инженеры не понимают, как измерить «дробное начальное значение» температуры или скорости). Поэтому для решения прикладных дифференциальных уравнений чаще используется определение Капуто. Производная Капуто сначала дифференцирует функцию до ближайшего целого числа, а затем берет от нее дробный интеграл. Это позволяет использовать стандартные, понятные классические начальные условия Коши, делая уравнения физически осмысленными.

Методы Грюнвальда-Летникова и проблема вычислительной памяти

Поскольку аналитически решить дробные дифференциальные уравнения удается лишь в тривиальных случаях, возникает острая необходимость в эффективных численных методах. Самым интуитивным подходом к дискретизации является схема Грюнвальда-Летникова. Она напрямую обобщает классические конечные разности. Если первая производная вычисляется по двум точкам, вторая — по трем, то для вычисления дробной производной по Грюнвальду-Летникову требуется бесконечный биномиальный ряд.

В этом кроется главная вычислительная сложность («проклятие») дробного исчисления. Классическая производная локальна: чтобы найти скорость автомобиля прямо сейчас, нам нужно знать его положение только за последнюю миллисекунду. Дробная же производная глобальна. Это математический оператор с эффектом «затухающей памяти». Чтобы вычислить состояние дробной системы на n-ном временном шаге, компьютеру необходимо учесть всю историю состояния системы, начиная с самого первого шага t=0! Это приводит к тому, что время вычислений растет пропорционально квадрату количества шагов (O(N^2)), а требования к оперативной памяти растут линейно. Для длительного интегрирования таких систем математики вынуждены разрабатывать хитрые методы «короткой памяти» (усечения хвоста) и применять алгоритмы быстрых дискретных сверток на базе БПФ.

Почему не подходят символьные и численные методы?

Для оптимизации сложных функций (например, при обучении глубоких нейронных сетей с миллиардами параметров) нам необходимо быстро и сверхточно вычислять градиенты. Классическая вычислительная математика исторически предлагала два подхода, но оба оказались непригодны для задач современных масштабов. Первый подход — символьное дифференцирование (которым занимаются системы компьютерной алгебры типа Mathematica). Программа аналитически применяет жесткие правила дифференцирования к исходной формуле. К сожалению, для сложных алгоритмов с циклами ветвления этот метод приводит к феномену «экспоненциального разбухания выражений», когда итоговая формула производной занимает гигабайты памяти. Второй подход — численное дифференцирование конечными разностями. Как мы обсуждали ранее, этот метод страдает от катастрофической потери точности из-за неизбежных ошибок округления и требует N+1 вызовов функции для функции от N переменных, что при миллионах переменных замедлит расчеты до бесконечности.

Третий, поистине революционный путь — это Автоматическое Дифференцирование (АД). Это не символьная аналитика и не численное приближение. АД основывается на том простом факте, что любая, даже самая сложная компьютерная программа в конечном итоге состоит из последовательности примитивных арифметических операций (сложение, умножение) и базовых функций (синус, экспонента). АД применяет цепное правило дифференцирования (chain rule) к этим элементарным шагам прямо в процессе выполнения программы, сохраняя идеальную машинную точность и не создавая раздутых математических формул.

Прямой проход и алгебра дуальных чисел

Автоматическое дифференцирование имеет два основных режима работы. Прямой режим (Forward mode) концептуально основан на использовании так называемых дуальных чисел — специальной алгебраической структуры вида a + b*epsilon, где число epsilon обладает уникальным свойством: его квадрат строго равен нулю (по аналогии с комплексной мнимой единицей, квадрат которой равен -1). Применяя арифметику дуальных чисел к коду программы, мы одновременно за один единственный вычислительный проход получаем и итоговое значение функции, и ее точную направленную производную. Прямой режим невероятно эффективен и быстр, если мы анализируем функцию с малым числом входов и огромным числом выходов.

Обратный проход (Reverse Mode) и алгоритм Backpropagation

Однако в машинном обучении ситуация строго обратная: функция потерь (Loss function) — это всего лишь одно единственное скалярное число (один выход), которое сложнейшим образом зависит от миллионов весовых коэффициентов (входов). Для решения таких задач прямой режим потребовал бы миллионов повторных проходов. Здесь на сцену выходит обратный режим АД (Reverse mode), который в индустрии искусственного интеллекта больше известен под термином «алгоритм обратного распространения ошибки» (Backpropagation).

Обратный режим АД работает в два математических этапа. Сначала выполняется обычный прямой проход, во время которого алгоритм строит в памяти гигантский направленный вычислительный граф и бережно запоминает промежуточные значения абсолютно всех переменных (что требует колоссального объема оперативной видеопамяти VRAM). Затем алгоритм начинает двигаться по этому графу в обратном направлении, от выхода к входам, передавая сигналы чувствительности (градиенты). Благодаря гениальному применению цепного правила с конца, обратный режим позволяет вычислить точный градиент функции сразу по всем миллионам параметров всего за один обратный проход графа! Без этого потрясающего вычислительного алгоритма современный бум генеративных нейронных сетей был бы абсолютно невозможен.

Когда одиночный выстрел уходит в бесконечность

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

Суть проблемы заключается в жесточайшей чувствительности дифференциальной системы к начальным данным. Если в задаче присутствуют быстрорастущие экспоненциальные компоненты, то даже микроскопическая ошибка (в миллиардных долях) при выборе начального угла «орудия» приведет к тому, что к середине расчетного интервала численное решение улетит в математическую бесконечность из-за переполнения регистров памяти компьютера. Алгоритм просто не сможет дойти до правого края, чтобы оценить ошибку промаха. Невозможно прицелиться, если любое малейшее дрожание руки приводит к отклонению снаряда на миллионы километров.

Метод многократной (параллельной) стрельбы

Чтобы обуздать экспоненциальный взрыв и восстановить вычислительную стабильность, математиками Келлером, Осборном и Деуфлхардом был разработан метод многократной стрельбы (Multiple Shooting Method). Идея этого метода заключается в том, чтобы не пытаться «прострелить» весь длинный и сложный интервал одним выстрелом.

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

Сшивка решений и матрица Якоби

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

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