Высшая математика

Глобальная оптимизация: генетические алгоритмы и эволюционные вычисления

Ловушка локальных минимумов

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

Для преодоления этой проблемы были созданы методы глобальной стохастической оптимизации. Одним из самых мощных, универсальных и красивых классов таких методов являются генетические алгоритмы (ГА), которые черпают свое вдохновение напрямую из биологической теории эволюции Чарльза Дарвина.

Биологическая метафора: хромосомы, популяции и отбор

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

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

Кроссовер и мутация: двигатели математического прогресса

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

Однако селекция и кроссовер могут привести к тому, что популяция быстро выродится и станет состоять из одинаковых особей, застряв в локальном экстремуме. Чтобы поддерживать генетическое разнообразие и исследовать новые, неизведанные участки пространства, применяется оператор мутации. С небольшой вероятностью (например, 1%) алгоритм случайно изменяет один или несколько генов (битов) у потомков. Мутация позволяет алгоритму случайно «перепрыгивать» через холмы целевой функции, вырываясь из ловушек локальных минимумов. Повторяя эти циклы (оценка, селекция, скрещивание, мутация) на протяжении сотен поколений, генетический алгоритм стягивает популяцию к глобальному оптимуму даже для недифференцируемых, разрывных и зашумленных функций, с которыми не справляется ни один градиентный метод.

Подробнее

Численное решение краевых задач: метод стрельбы и метод прогонки

Отличие краевой задачи от задачи Коши

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

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

Метод стрельбы: искусство прицеливания в математике

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

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

Метод конечных разностей и алгоритм прогонки

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

В результате дифференциальное уравнение превращается в систему линейных алгебраических уравнений. Особенность этой системы в том, что ее матрица является трехдиагональной (каждое уравнение связывает только три соседних узла: левый, центральный и правый). Для решения таких СЛАУ применяется специальный, невероятно быстрый и элегантный алгоритм Томаса (в русскоязычной литературе известный как метод прогонки). Прогонка состоит из прямого хода (последовательного исключения поддиагональных элементов) и обратного хода (нахождения неизвестных). Этот метод требует всего O(N) арифметических операций и является абсолютно устойчивым для задач с диагональным преобладанием, что делает его индустриальным стандартом в задачах теплопроводности и сопротивления материалов.

Подробнее

Дискретное вейвлет-преобразование: альтернатива Фурье в анализе сигналов

Анализ сигналов: выход за пределы синусоид

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

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

Вейвлеты: маленькие волны с конечной энергией

Этим новым инструментом стало вейвлет-преобразование, активно развивавшееся в 1980-х годах (труды Морле, Гроссмана, Ив Мейер и Ингрид Добеши). Слово «вейвлет» (wavelet) переводится как «маленькая волна» или «всплеск». В отличие от синусоиды, вейвлет отличен от нуля только на коротком отрезке времени, а за его пределами быстро затухает до нуля.

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

Дискретное преобразование и индустриальные стандарты

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

Благодаря своей способности концентрировать энергию сигнала в небольшом количестве коэффициентов (большинство коэффициентов ДВП для реальных сигналов оказываются близки к нулю), вейвлеты совершили революцию в сжатии данных. Именно вейвлет-преобразование (конкретно, вейвлеты Коэна-Добеши-Фово) лежит в основе графического стандарта JPEG 2000, обеспечивая сжатие изображений высокого разрешения без характерных для обычного JPEG раздражающих квадратных «артефактов блочности». Также вейвлеты стали золотым стандартом для очистки медицинских и акустических сигналов от шумов (алгоритмы wavelet denoising/shrinkage).

Подробнее

Итерационные методы подпространств Крылова для гигантских матриц

Как решать системы с миллиардами уравнений?

При дискретизации трехмерных дифференциальных уравнений в частных производных (например, в аэродинамике или гидрогеологии) методом конечных элементов возникают колоссальные системы линейных алгебраических уравнений (СЛАУ). Размерность таких матриц может достигать сотен миллионов и даже миллиардов строк. Использование прямого метода Гаусса для таких гигантов невозможно: даже суперкомпьютеру не хватит оперативной памяти для хранения матрицы, не говоря уже о годах, необходимых для расчетов.

К счастью, такие матрицы являются крайне разреженными. Более 99% их элементов — это нули. Для их хранения используются специальные форматы (например, CSR — Compressed Sparse Row), в которых хранятся только ненулевые элементы и их координаты. Однако прямые методы в процессе исключения неизбежно превращают нули в ненулевые элементы (проблема заполнения). Поэтому для решения таких огромных разреженных СЛАУ применяются исключительно современные итерационные методы, вершиной эволюции которых являются алгоритмы на основе подпространств Крылова.

Подпространства Крылова: выжимание максимума из умножения матрицы на вектор

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

Если мы возьмем начальный вектор невязки (r) и начнем последовательно умножать его на матрицу A, мы получим последовательность векторов: r, Ar, A^2r, A^3r и так далее. Линейная оболочка (пространство, натянутое на эти векторы) называется подпространством Крылова. Гениальность подхода заключается в том, что вместо решения системы в исходном пространстве с миллионами измерений, мы ищем приближенное решение (проекцию) в подпространстве Крылова, размерность которого мала (например, 50 или 100). На каждом шаге (итерации Арнольди или Ланцоша) к подпространству добавляется новый ортогональный вектор, и приближение становится все точнее.

GMRES и Метод сопряженных градиентов (CG)

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

Однако если матрица несимметричная (что часто бывает в задачах конвекции-диффузии и гидродинамики), метод CG не работает. Для таких матриц в 1986 году Саадом и Шульцем был разработан метод GMRES (Generalized Minimal Residual method). Метод GMRES минимизирует норму невязки на всем подпространстве Крылова. В отличие от CG, он вынужден хранить в памяти все базисные векторы, поэтому на практике применяется с «рестартами» (GMRES(m)): после достижения m итераций накопленные векторы сбрасываются, и процесс начинается заново с текущего наилучшего приближения. Сегодня алгоритмы подпространств Крылова, оснащенные мощными предобуславливателями (preconditioners) типа ILU или Algebraic Multigrid (AMG), являются абсолютным индустриальным стандартом при суперкомпьютерном моделировании.

Подробнее

Численное решение интегральных уравнений: методы Фредгольма и Вольтерры

Когда неизвестная функция скрыта под знаком интеграла

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

Интегральные уравнения принято классифицировать по двум признакам. Если пределы интегрирования фиксированы (постоянны), то это уравнение Фредгольма. Если верхний предел является переменной интегрирования (то есть интеграл зависит от текущего положения), то это уравнение Вольтерры, описывающее процессы эволюции с памятью. Если неизвестная функция находится только под знаком интеграла, это уравнение первого рода; если она присутствует еще и вне интеграла — это уравнение второго рода. Как и в случае с ОДУ, аналитически решить интегральные уравнения удается лишь в редких частных случаях, что делает применение численных методов абсолютно необходимым.

Метод квадратур: сведение к СЛАУ

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

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

Метод последовательных приближений (итерации Пикара)

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

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

Подробнее

Адаптивные методы интегрирования ОДУ: алгоритмы Рунге-Кутты-Фельберга

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

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

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

Вложенные методы и оценка локальной погрешности

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

Если бы мы делали это «в лоб», объем вычислений возрос бы в два раза. Гениальность идеи Фельберга (и последующих авторов, таких как Дорман и Принс) состоит в том, что коэффициенты для методов 4-го и 5-го порядков подбираются специальным образом. Они используют одни и те же промежуточные вычисления (те самые коэффициенты наклона касательных внутри шага). То есть, вычислив набор наклонов для метода 5-го порядка, мы можем путем простого перемножения с другими весами получить и результат для 4-го порядка практически бесплатно! Разность между результатами 5-го и 4-го порядка служит превосходной оценкой локальной ошибки (погрешности усечения) на данном шаге.

Алгоритм управления шагом (Step Size Control)

Имея оценку локальной ошибки, программа-решатель включает встроенный ПИД-регулятор шага. Пользователь заранее задает желаемую допустимую погрешность (tolerance). Алгоритм сравнивает вычисленную оценку ошибки с заданным допуском.

Если ошибка оказалась больше допуска, это означает, что шаг был слишком велик. Алгоритм отбрасывает вычисленное значение, уменьшает шаг по специальной формуле (обычно умножая на коэффициент около 0.8 от расчетного оптимального, чтобы иметь запас надежности) и пересчитывает этот участок заново. Если же ошибка оказалась значительно меньше допуска, алгоритм принимает шаг как успешный, но для следующего шага автоматически увеличивает размер интервала, ускоряя симуляцию. Метод Рунге-Кутты-Фельберга (RKF45) и его более современная модификация, метод Дормана-Принса (DOPRI5), стали мировым стандартом адаптивного интегрирования и используются по умолчанию во встроенных функциях MATLAB (ode45) и библиотеках SciPy в Python.

Подробнее

Методы многомерной оптимизации: градиентный спуск и метод Ньютона

Поиск наилучшего решения в мире множества вариантов

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

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

Градиентный спуск: движение по наискорейшему склону

Самым известным и концептуально простым алгоритмом является метод градиентного спуска (Gradient Descent), предложенный еще Огюстеном Луи Коши. Градиент — это вектор, указывающий направление наискорейшего возрастания функции. Следовательно, антиградиент (вектор, взятый со знаком минус) указывает направление наискорейшего убывания.

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

Метод Ньютона и учет кривизны пространства

Чтобы ускорить сходимость и избежать зигзагов в оврагах, необходимо использовать информацию не только о наклоне функции (первой производной), но и о ее кривизне. Эту информацию предоставляет матрица вторых производных, называемая матрицей Гессе (или гессианом). Использование гессиана лежит в основе многомерного метода Ньютона.

Метод Ньютона строит на каждом шаге квадратичную аппроксимацию целевой функции (в виде параболоида) и мгновенно прыгает в вершину этого параболоида. В окрестности минимума метод Ньютона обладает фантастической квадратичной скоростью сходимости. Однако вычисление, хранение и, что самое главное, обращение матрицы Гессе на каждом шаге (требующее решения огромной СЛАУ) делает классический метод Ньютона абсолютно неприменимым для задач с большим количеством параметров (например, при обучении нейросетей). Поэтому на практике используют квазиньютоновские методы (такие как популярный BFGS) или метод сопряженных градиентов, которые хитроумными способами оценивают кривизну пространства, не вычисляя сам гессиан в явном виде.

Подробнее

Численная линейная алгебра: вычисление собственных значений и векторов

Скрытые резонансы: проблема спектрального анализа

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

Физический и инженерный смысл этой задачи колоссален. Собственные значения описывают частоты свободных (резонансных) колебаний механических систем: от гитарной струны и подвески автомобиля до мостов и небоскребов. Совпадение внешней нагрузки с этими частотами вызывает разрушительный резонанс. В квантовой механике собственные значения оператора Гамильтона определяют дискретные уровни энергии атомов (отсюда происходит термин «спектр»). В теории графов и анализе социальных сетей анализ спектра матриц смежности позволяет находить скрытые кластеры и самых влиятельных участников.

Степенной метод для поиска доминирующего значения

Аналитически собственные значения находятся путем решения так называемого характеристического уравнения (поиска корней полинома, степень которого равна размерности матрицы). Однако, как мы знаем из теоремы Абеля-Руффини, для матриц размером 5x5 и больше точных формул корней не существует. Поэтому в реальности собственные значения всегда ищут только итерационными численными методами.

Самым простым подходом является степенной метод (итерации фон Неймана). Он позволяет найти наибольшее по модулю (доминирующее) собственное значение. Идея поразительно проста: мы берем случайный начальный вектор и начинаем многократно умножать его на нашу матрицу. На каждом шаге вектор будет поворачиваться в сторону собственного вектора, соответствующего максимальному собственному значению. После достаточного количества итераций (и периодической нормализации вектора, чтобы числа не ушли в бесконечность) процесс сходится к искомому главному собственному вектору. Именно этот простой принцип лежит в основе первоначального алгоритма PageRank, созданного Ларри Пейджем и Сергеем Брином для ранжирования веб-страниц в поисковой системе Google, где матрицей выступает граф ссылок всего интернета.

Алгоритм QR-разложения для нахождения всего спектра

Степенной метод хорош, если нам нужно только одно максимальное значение. Но что делать, если инженеру нужен весь спектр частот (все собственные значения)? Для этого в 1961 году Джоном Фрэнсисом и Верой Кублановской был независимо разработан QR-алгоритм, который считается одним из 10 самых важных численных алгоритмов 20-го века.

Алгоритм основан на факторизации (разложении) матрицы. На каждом шаге исходная матрица A представляется в виде произведения ортогональной матрицы Q и верхнетреугольной матрицы R (A = QR). Затем эти матрицы умножаются в обратном порядке: строится новая матрица A1 = RQ. Удивительный математический факт состоит в том, что новая матрица A1 имеет точно такие же собственные значения, что и исходная матрица A (они подобны). Если повторять этот процесс итеративно (находя QR-разложение и перемножая в обратном порядке), то матрица постепенно стремится к диагональному (или блочно-диагональному) виду, где все ее собственные значения просто стоят на главной диагонали. С применением предварительного приведения к форме Хессенберга и стратегий сдвигов, QR-алгоритм работает феноменально быстро и стабильно, являясь сегодня ядром библиотек линейной алгебры, таких как LAPACK, встроенных в Python, MATLAB и R.

Подробнее

Метод Монте-Карло: стохастический подход к вычислениям и интегрирование

Случайность на службе строгой математики

Численные методы, которые мы обсуждали до сих пор, являются детерминированными: они следуют строгим формулам, и при одинаковых входных данных всегда дают абсолютно идентичный результат с точностью до бита. Однако в середине 20-го века, во время работы над Манхэттенским проектом, физики Станислав Улам, Джон фон Нейман и Николас Метрополис предложили радикально иной подход, названный кодовым словом «Монте-Карло» в честь знаменитого казино. Этот метод использует генерацию случайных чисел для решения сугубо детерминированных математических задач.

Суть метода Монте-Карло заключается в проведении огромного числа виртуальных стохастических экспериментов. Классический пример — вычисление площади фигуры сложной формы или числа Пи. Если вписать фигуру в квадрат с известной площадью, а затем «бросать» случайные точки в этот квадрат, то отношение числа точек, попавших внутрь фигуры, к общему числу брошенных точек будет стремиться к отношению их площадей. Закон больших чисел теории вероятностей гарантирует, что при стремлении числа испытаний к бесконечности мы получим точное математическое решение.

Проклятие размерности в задачах интегрирования

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

Если мы попытаемся использовать классические детерминированные методы (например, метод Симпсона) для вычисления D-мерного интеграла, мы столкнемся с «проклятием размерности». Если на каждую ось мы поместим по 10 узлов сетки, то для 3-мерного объема потребуется 10^3 = 1000 узлов. Но для 100-мерного интеграла потребуется 10^100 вычислений функции — число, превышающее количество атомов во Вселенной. Ни один суперкомпьютер с этим не справится. Ошибка классических методов зависит от размерности пространства. Метод Монте-Карло же обладает уникальным свойством: скорость его сходимости (убывания ошибки) равна O(1/sqrt(N)), где N — количество испытаний, и эта скорость абсолютно не зависит от размерности пространства! Это делает его единственным рабочим инструментом для задач высокой размерности.

Алгоритмы уменьшения дисперсии и цепи Маркова (MCMC)

Главный недостаток базового метода Монте-Карло — его медленная сходимость. Из-за корня в формуле сходимости, чтобы уменьшить ошибку в 10 раз, нужно увеличить количество испытаний в 100 раз. Чтобы бороться с этим, математики разработали изощренные техники снижения дисперсии (variance reduction techniques).

Одним из мощнейших инструментов является метод выборки по значимости (importance sampling). Вместо того чтобы разбрасывать случайные точки равномерно по всему объему, алгоритм концентрирует выборку в тех областях, где подынтегральная функция вносит наибольший вклад в итоговый результат. Другим прорывом стали алгоритмы Монте-Карло по схеме марковских цепей (Markov Chain Monte Carlo, MCMC), в частности алгоритм Метрополиса-Гастингса. В этом случае случайные блуждания не являются независимыми: следующая точка выбирается на основе положения предыдущей, формируя интеллектуальный поиск в многомерном пространстве параметров. Алгоритмы MCMC произвели революцию в байесовском машинном обучении, искусственном интеллекте и расшифровке генома, доказав, что контролируемая случайность — один из самых мощных инструментов познания.

Подробнее

Решение жестких систем обыкновенных дифференциальных уравнений (ОДУ)

Феномен жесткости: когда классические методы дают сбой

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

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

Неявные методы и А-устойчивость

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

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

Формулы дифференцирования назад (BDF) и методы Гира

Простейшие неявные методы обладают низким порядком точности (первым или вторым). Для точного моделирования жестких систем требуются методы более высоких порядков. Американский математик Чарльз Гир в начале 1970-х годов разработал семейство жестко-устойчивых алгоритмов, известных как формулы дифференцирования назад (Backward Differentiation Formulas, BDF).

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

Подробнее