Геометрия скрытых структур в больших данных

В эпоху Big Data исследователи часто работают с наборами данных (Point Clouds), которые представляют собой миллионы точек в пространствах очень высокой размерности (сотни и тысячи координат). Это могут быть профили экспрессии генов, данные функциональной МРТ мозга или финансовые транзакции. Традиционные методы машинного обучения (кластеризация, PCA) опираются на метрику — точные расстояния между точками. Но в многомерных пространствах метрика теряет свою силу, шумы искажают расстояния, и методы дают сбои. Нужен был инструмент, устойчивый к растяжениям, сжатиям и непрерывным деформациям данных.

Этим инструментом стала алгебраическая топология, адаптированная для компьютеров под названием Топологический Анализ Данных (TDA). Топология изучает форму объектов (связность, наличие дыр, туннелей и пустот). Например, топологически кофейная кружка и бублик — это один и тот же объект, так как у обоих есть ровно одна дырка. Главный численный алгоритм TDA называется Устойчивыми гомологиями (Persistent Homology). Он позволяет компьютеру находить глобальные топологические структуры в хаосе дискретных многомерных точек.

Симплициальные комплексы и фильтрация Виеториса-Рипса

Как из разрозненных точек собрать сплошной топологический объект? Программа создает математические связи (симплексы). Точка — это симплекс 0-мерный. Отрезок между двумя точками — 1-мерный. Треугольник (с заливкой) — 2-мерный, тетраэдр — 3-мерный. Алгоритм начинает «раздувать» вокруг каждой точки данных гиперсферу радиуса R.

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

Вычисление штрих-кодов и матричная редукция

С вычислительной точки зрения процесс отслеживания гомологий сводится к тяжелой линейной алгебре. Для каждого шага фильтрации компьютер строит матрицу границ (Boundary matrix), которая описывает, какие вершины образуют какие грани. Элементы этой матрицы лежат в поле вычетов по модулю 2 (цифры только 0 и 1).

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

Плотно заполненные матрицы: проклятие интегральных операторов

При численном решении интегральных уравнений (например, уравнений Фредгольма) или при использовании метода граничных элементов (МГЭ) в электродинамике инженеры сталкиваются со страшной проблемой. В отличие от дифференциальных операторов, которые связывают только соседние узлы сетки и порождают удобные разреженные матрицы с нулями, интегральные операторы глобальны. Каждая точка на поверхности антенны электромагнитно взаимодействует со всеми остальными точками антенны. В результате получается СЛАУ, матрица которой заполнена на 100% (dense matrix). Если сетка содержит миллион узлов, матрица будет содержать триллион ненулевых элементов, что потребует 8 терабайт оперативной памяти и сделает матричное умножение катастрофически медленным.

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

От плотной матрицы к разреженной через пороговую фильтрацию

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

После преобразования матрица меняет свое физическое лицо. Теперь ее элементы — это не интенсивность взаимодействия узлов, а коэффициенты взаимодействия различных вейвлет-масштабов. Поскольку оператор (например, функция Грина) гладкий вдали от сингулярности, подавляющее большинство (более 95-99%) этих новых вейвлет-коэффициентов оказываются микроскопически малыми числами, очень близкими к нулю. Алгоритм просто обнуляет все коэффициенты, которые по модулю меньше заданного порога отсечения (Thresholding). Плотно заполненная матрица волшебным образом превращается в сильно разреженную (Sparse matrix) без потери физической точности решения! Теперь для ее хранения требуется ничтожный объем памяти, а итерационные решатели (например, GMRES) выполняют умножение матрицы на вектор за время O(N log N) вместо O(N^2), превращая нерешаемые электромагнитные задачи в повседневные инженерные расчеты.

Аппроксимация непериодических краевых задач

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

Для непериодических задач идеальным базисом являются ортогональные полиномы Чебышева. Метод спектральной коллокации Чебышева (Chebyshev Pseudospectral Method) берет лучшее от метода конечных разностей (простота работы в узлах) и спектральных методов (экспоненциальная сходимость). Ключевой особенностью метода является расположение расчетных узлов. Мы не можем использовать равномерную сетку из-за феномена Рунге (дикой раскачки полиномов на краях). Вместо этого используются точки Гаусса-Лобатто-Чебышева, которые сгущаются к границам отрезка по закону косинуса: x_j = cos(pi * j / N). Такое сгущение идеально компенсирует поведение многочленов, обеспечивая абсолютно стабильную и фантастически точную аппроксимацию.

Матрицы спектрального дифференцирования

Гениальность метода коллокации заключается в том, как он вычисляет производные. Если в методе конечных разностей мы берем значения функции в трех соседних точках и делим их на шаг h, то в спектральном методе Чебышева мы берем значения функции во всех N узлах сетки одновременно! Это глобальный подход.

Алгоритм заранее строит так называемую матрицу спектрального дифференцирования D_N. Эта матрица полностью заполнена (плотная). Чтобы найти производную функции во всех узлах сетки, мы просто умножаем вектор значений функции на эту матрицу. Элементы этой матрицы рассчитываются аналитически (или через Быстрое преобразование Фурье) и зависят от взаимного расположения узлов Чебышева. Если нам нужна вторая производная (для уравнения Пуассона), мы просто возводим матрицу D_N в квадрат! После подстановки граничных условий задача мгновенно сводится к решению СЛАУ. Скорость сходимости потрясает воображение: если решение является гладкой (аналитической) функцией, то ошибка убывает быстрее любой степени числа узлов. Там, где методу конечных разностей потребовались бы сотни тысяч точек, методу Чебышева достаточно 30-40 узлов для достижения машинной точности (10^-15). Этот алгоритм незаменим в теории гидродинамической устойчивости (уравнение Орра-Зоммерфельда) и физике плазмы.