Main menu

Метод дифференциальной эволюции: векторная мутация в глобальной оптимизации

Оптимизация в непрерывных недифференцируемых пространствах

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

Для таких задач были созданы стохастические генетические алгоритмы (ГА). Классические ГА работают путем кодирования параметров в виде бинарных строк (хромосом) из нулей и единиц. Однако при оптимизации непрерывных вещественных параметров (например, длина волны лазера, концентрация реагента в молях, напряжение в вольтах) бинарное кодирование порождает массу проблем, связанных с потерей точности и эффектом Хэмминга. В 1995 году Кеннет Прайс и Райнер Сторн разработали революционный алгоритм — Дифференциальную эволюцию (Differential Evolution, DE). Этот алгоритм изначально создавался для работы напрямую с векторами чисел с плавающей запятой, устраняя необходимость в громоздких бинарных преобразованиях.

Векторная разность как направляющая сила мутации

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

Для того чтобы создать нового «мутанта» из текущей особи (целевого вектора), алгоритм случайным образом выбирает из популяции еще три других, абсолютно независимых вектора: A, B и C. Затем он вычисляет векторную разность между B и C. Эта математическая разность масштабируется (умножается на постоянный коэффициент мутации F) и прибавляется к вектору A. Полученный вектор становится основой для нового потомка. В чем гениальность этого подхода? Векторная разность автоматически адаптирует алгоритм к топологии ландшафта! В начале работы, когда популяция разбросана далеко друг от друга, разности огромны, и алгоритм делает гигантские разведывательные шаги, перепрыгивая через любые локальные минимумы. Когда алгоритм нащупывает глобальную долину, популяция стягивается в плотный клубок, разности между векторами становятся микроскопическими, и алгоритм переходит в режим ювелирной, сверхточной настройки параметров, не требуя ручного изменения радиуса мутации (в отличие от методов имитации отжига).

Скрещивание, селекция и индустриальное применение

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

Всего три простых математических операции — сложение векторов, масштабирование разности и условная замена. Алгоритм дифференциальной эволюции умещается в 20 строк программного кода, требует настройки всего двух интуитивных параметров (F и CR) и феноменально легко распараллеливается на многоядерных процессорах. Благодаря своей беспрецедентной способности находить глобальный минимум сложных невыпуклых функций, метод DE стал основным алгоритмом для настройки гиперпараметров нейронных сетей, автоматического синтеза фрактальных Wi-Fi антенн сложной формы и проектирования молекул новых лекарственных препаратов в хемоинформатике.

Подробнее

Методы Ланцоша и Арнольди: поиск собственных значений гигантских матриц

Охота за экстремальными резонансами

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

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

Подпространство Крылова и ортогональный базис

Оба метода базируются на элегантной концепции подпространств Крылова, которую мы упоминали в контексте решения СЛАУ. Алгоритм выбирает случайный начальный вектор q1 и начинает многократно умножать его на гигантскую матрицу A. В результате формируется последовательность векторов: q1, A*q1, A^2*q1, A^3*q1. Эти векторы несут в себе колоссальную математическую информацию о доминирующих собственных векторах матрицы (так как при многократном умножении вектор вытягивается вдоль главного собственного направления).

Проблема в том, что эти векторы очень быстро становятся почти коллинеарными (параллельными), что уничтожает их вычислительную полезность (матрица их скалярных произведений становится вырожденной). Гениальность метода Арнольди (Уолтер Арнольди, 1951) заключается в применении процесса ортогонализации Грама-Шмидта прямо на лету! Каждый новый сгенерированный вектор тут же жестко проецируется и очищается от компонентов всех предыдущих векторов. В результате мы получаем идеальный ортонормированный базис подпространства Крылова. В этом новом, маленьком базисе исходная гигантская матрица проецируется (сжимается) в крошечную матрицу Хессенберга размера m x m (где m — число итераций, например, 50 или 100).

Метод Ланцоша и проблема потери ортогональности

Если исходная матрица системы была симметричной (как в задачах сопротивления материалов), процесс Арнольди превращается в Метод Ланцоша (Корнелий Ланцош, 1950). Из-за свойств симметрии матрица Хессенберга магическим образом схлопывается в симметричную трехдиагональную матрицу. Это означает, что для генерации каждого нового ортогонального вектора алгоритму Ланцоша не нужно вспоминать и вычитать все предыдущие 100 векторов — ему достаточно использовать короткую рекуррентную формулу, опирающуюся только на два последних вектора! Это делает метод Ланцоша феноменально быстрым и экономным по памяти.

Получив эту крошечную матрицу проекции, мы просто решаем ее стандартным QR-алгоритмом за миллисекунды. Ее собственные значения (называемые числами Ритца) с фантастической скоростью сходятся к истинным экстремальным собственным значениям гигантской исходной матрицы. Однако у метода Ланцоша есть вычислительный недуг: из-за конечной разрядности компьютеров (ошибок округления float) генерируемые векторы постепенно теряют свою взаимную перпендикулярность. Это порождает появление фиктивных (призрачных) собственных значений, которые являются математической галлюцинацией. Для борьбы с этим современная вычислительная алгебра использует методы полной или выборочной переортогонализации, которые встроены в ядро таких сверхмощных индустриальных пакетов, как ARPACK (используемых функциями `eigs` в MATLAB и SciPy).

Подробнее

Векторизация вычислений (SIMD): аппаратное ускорение численных алгоритмов

Когда закон Мура перестает ускорять код

На протяжении десятилетий программисты и математики привыкли к простому правилу: чтобы программа для численного интегрирования или умножения матриц заработала быстрее, нужно просто подождать пару лет, пока Intel или AMD выпустят процессор с более высокой тактовой частотой. Этот феномен, подпитываемый законом Мура, позволял писать неоптимальный код и все равно получать рост производительности. Однако в начале 2000-х годов рост тактовых частот уперся в непреодолимый физический барьер — тепловой предел кремниевых чипов (частота застыла на уровне 3-5 ГГц). Чтобы продолжать наращивать вычислительную мощь, создателям процессоров пришлось менять архитектуру, внедряя многоядерность и, что еще более важно для математики, векторные инструкции.

Концепция SIMD (Single Instruction, Multiple Data — одна инструкция, множество данных) является фундаментом высокопроизводительных вычислений (HPC). В классической, скалярной модели процессора (SISD), чтобы сложить два массива из четырех чисел (A1+B1, A2+B2, A3+B3, A4+B4), процессору нужно выполнить цикл: загрузить два числа, сложить их, записать результат, и повторить это 4 раза. Архитектура SIMD внедряет в процессор сверхширокие регистры. Загрузив сразу 4 числа из массива A в один 256-битный регистр, и 4 числа из массива B во второй регистр, процессор за один единственный такт времени выполняет сложение сразу всех четырех пар чисел одновременно! Это обеспечивает мгновенный 4-кратный (или даже 8-кратный) прирост чистой математической производительности.

Наборы инструкций AVX и выравнивание памяти

Современные процессоры оснащены мощнейшими блоками векторной математики. Инструкции SSE (128 бит) позволяют складывать 4 числа типа float за такт. Инструкции AVX и AVX2 (256 бит) удваивают этот показатель до 8 чисел. Новейшие серверные процессоры с инструкциями AVX-512 оперируют колоссальными 512-битными регистрами, выполняя по 16 операций сложения или умножения 32-битных чисел с плавающей запятой за один микроскопический такт процессора. Более того, существует специальная инструкция FMA (Fused Multiply-Add), которая за один такт выполняет сразу две операции (A * B + C), удваивая пиковую производительность алгоритма матричного умножения.

Но эта фантастическая аппаратная мощь недоступна программам автоматически. Чтобы векторные блоки заработали, данные в оперативной памяти компьютера должны быть расположены (упакованы) строго последовательно, без разрывов, и желательно быть выровненными (aligned) по границам 32-байтовых или 64-байтовых кеш-линий. Если вы написали алгоритм метода конечных разностей, который скачет по памяти (собирая данные из разных классов или объектов), процессор просто не сможет загрузить их широким вектором. Он откатится к медленным скалярным вычислениям (потеряв до 90% скорости).

Компиляторы, раскрутка циклов и библиотеки BLAS

Чтобы заставить численный код летать, разработчики применяют агрессивную раскрутку циклов (Loop Unrolling) и прагмы для компиляторов (например, #pragma omp simd в OpenMP). Важнейшим требованием для векторизации является отсутствие логических ветвлений (операторов if/else) внутри тяжелого математического цикла, так как ветвление разрушает векторный конвейер.

В мире профессиональной вычислительной математики никто не пишет циклы умножения матриц вручную. Вся базовая алгебра (нахождение скалярных произведений, умножение матрицы на вектор, разложение матриц) делегируется стандартизированным библиотекам BLAS (Basic Linear Algebra Subprograms) и LAPACK. Такие реализации, как Intel MKL (Math Kernel Library) или OpenBLAS, написаны на глубоком ассемблерном уровне. Они досконально учитывают размер L1/L2 кэша процессора и микроархитектуру AVX регистров, нарезая матрицы на микро-блоки. Именно благодаря идеальной SIMD-векторизации внутри библиотек BLAS, скрипты на языке Python (использующие NumPy) умудряются работать со скоростью скомпилированного языка C++ при решении гигантских систем дифференциальных уравнений.

Подробнее

Соц. сети