Вычислительная линейная алгебра

Линейная алгебра — это основа вычислительной науки. Как матанализ — основа математики в каком-то смысле, то вычислительная линейная алгебра — это основа вычислительной науки. Может показаться, что это что-то простое. Что такое вообще линейная алгебра? Это наука о векторных пространствах. Мы можем складывать векторы, есть матрицы или операторы, мы можем умножать матрицу на вектор, решать линейные системы. Выглядит просто.

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

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

Чем занимается вычислительная линейная алгебра? Опять же история из моей жизни. Когда я был на третьем курсе, встретил своего будущего научного руководителя, он мне пытался рассказать, чем занимается. Он говорит: «Я занимаюсь умножением матрицы на вектор». В этот момент мой энтузиазм упал до нуля, потому что это какой-то идиотизм — умножать матрицу на вектор. Это мы умеем. То есть решать линейные системы — просто, умножать матрицу на вектор — тоже вроде тривиально: строчка на столбец, строчка на столбец. Что тут может быть сложного?

Но есть одна небольшая проблемка: эти матрицы очень большие, это матрицы миллион на миллион, миллиард на миллиард. Если прикинуть, сколько нужно операций, чтобы это сделать, это очень много операций. Вы даже эту матрицу целиком не сможете сохранить в памяти компьютера. Если это матрица миллион на миллион, то у вас получается 1012 элементов, каждый элемент занимает в двойной точности 8 байт, 8 на 1012. Сейчас есть параллельные компьютеры, большие диски, вы можете это сохранить. Но вы хотите это делать на обычном компьютере, на ноутбуке, в своей обычной жизни решать системы такими матрицами. Это нормальная задача. Если вы рассмотрите, скажем, двумерную поверхность, разобьете ее здесь на тысячу частей, здесь на тысячу частей — у вас как раз получится миллион элементов, и матрица будет иметь порядок 106.

Здесь собака зарыта: как умножать матрицу на вектор, если матрица большая и, не дай бог, еще и неразреженная.

Разреженная матрица — это когда у вас матрица с большим количеством нулей, то есть все дифференциальные операторы связаны с разрядами, а если у вас оператор является интегральным, то есть матрица плотная, все элементы не нули, все их надо вроде как считать, — эта задача кажется неподъемной. Такого рода задачи научились решать не так давно, где-то в конце 80-х годов.

Есть список лучших достижений XX века. В области вычислительных методов, в области линейной алгебры это быстрый мультипольный метод Грингарда и Рохлина. Грингард сейчас директор Курантовского института в Нью-Йорке, активный ученый, Рохлин тоже активный ученый. Они придумали метод, как можно то, что делается за n2 операций, где n — это число элементов в векторе, тот самый миллион, сделать за О (n*logn). Это был прорыв. Единственное отличие, что теперь мы вычисляем матрично-векторные произведения не точно, а с какой-то заданной точностью, мы задаем некоторый порог точности. Скажем, мы хотим три знака после запятой и получаем ответ с точностью три знака после запятой. Потом мы говорим: хотим пять знаков после запятой. И нам надо посчитать уже подольше. Тем не менее это все равно О (n*logn), умноженный на какой-то множитель, который зависит от точности.

Это смена парадигмы, что нужно те операции, к которым мы привыкли, которым учат в университете, которые выполняются точно, с точностью до машинной погрешности, выполнять приближенно с некоторой точностью, которая нужна в приложении. Часто, например, в приложении не нужно 16 знаков после запятой, достаточно иметь 2–3 знака, и это хорошо.

Какая математика зашита за этими быстрыми алгоритмами? Если вернуться к университетскому курсу — и это беда не только российских университетов, это беда многих мировых университетов, — линейную алгебру учат неправильно. Например, там встречается такое замечательное понятие, как «определитель». Определитель учат, но определитель в линейной алгебре не используется нигде и никогда, это чисто теоретическое понятие. Жуткая неустойчивость с вычислительной точки зрения, когда мы хотим что-то сделать устойчивым. Что значит устойчивость? У вас есть какое-то вычисление, вы выполняете его на компьютере в конечной точности и хотите получить результат с гарантированной погрешностью. Если алгоритм неустойчив, вы этот результат получите неправильно. У вас алгоритм работает правильно в точной арифметике, а на компьютере в условиях машинной погрешности он работает неточно. Решать ту же линейную систему по правилу Крамера, если 2*2, значит, нужно составить определитель, составить так называемые миноры, — это делать не надо.

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

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

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

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

Есть замечательное матричное разложение — сингулярное разложение. Это самое важное матричное разложение, на нем основано практически все.

Оно позволяет заданную матрицу, которая определяется n2 параметров, во многих случаях приближать так называемыми матрицами малого ранга. То есть матрицы малого ранга определяются небольшим количеством параметров. Есть некий параметр, который называется рангом, и мы хотим, чтобы этот ранг был небольшим. В исходной матрице, если она квадратная, n2 чисел, это понятно. Малоранговая матрица определяется 2*n*r параметров. Если ранг у вас маленький, вы можете существенно уменьшить количество параметров.

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

Если возвращаться к теореме, что матриц малого ранга не бывает, если вы возьмете случайную матрицу, она будет не матрицей малого ранга, тут можно ответить: если вы берете случайную матрицу, она с вероятностью 1 будет плотной, то есть все элементы не 0, и с вероятностью 1 она будет почти ортогональной, то есть она хорошо вращается. На практике, если вы берете практически любую задачу, матрица будет разреженная и плохо обусловленная. Практические задачи и случайные матрицы не имеют ничего между собой общего в этом смысле. Можно на философском уровне объяснить, почему матрицы малого ранга возникают во многих приложениях, то есть они возникают и в факторном анализе, они возникают и в различных приложениях, связанных с решением дифференциальных интегральных уравнений, и во многих-многих других. В некотором смысле сингулярное разложение идеально с обеих точек зрения. У вас бывают алгоритмы быстрые, которые задачи решают быстро, но неустойчивые, неточные, а бывают алгоритмы точные, но при этом медленные. Сингулярное разложение — оно и быстрое, и точное.

Линейная алгебра в последние годы, где-то, наверное, в 80-х, переживает второе рождение, связанное как раз с малоранговыми матрицами, с матрицами, которые имеют специальную структуру. Например, такая вещь, как тёплицева матрица. Это означает, что элементы матриц находятся вдоль диагонали — диагональная матрица. Структура матриц обусловлена приложениями. Если это тёплицева матрица, значит, у вас какая-то величина зависит только от расстояния между двумя объектами — например, взаимодействие, которое зависит только от радиуса.

Матрицы малого ранга связаны с некоторыми свойствами геометрического разделения. Это приводит к новым задачам, к новым теоремам и к новой науке.

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

Есть классическая задача при умножении матриц: строчка на столбец. Есть большое количество пакетов, которые нацелены на то, чтобы сделать это наиболее эффективно на современных компьютерах. В большом проекте, который долго развивался, в ATLAS’е, — там один из лидеров был Джек Донгарра, — пытались настроить какие-то параметры алгоритма на архитектуры. Его много лет развивали, и где-то в 2003 году один человек, японец по фамилии Гото, предложил новый алгоритм. Обычно говорят, что не надо писать свой алгоритм для базовых задач умножения матриц, это уже всегда сделано в пакетах, не пытайтесь даже. Но если у вас есть идея, как это улучшить, это иногда стоит попробовать. Один человек написал пакет, и он оказался быстрее на 20–30%, чем пакеты, которые разрабатывались большими коммерческими компаниями, большими консорциумами. У него была одна маленькая идея, связанная с программистскими вещами, с алгеброй, он смог это сделать. Даже те алгоритмы, которые, кажется, уже в бронзе отлиты, — это далеко не так. Если у вас есть идея, вы можете попробовать. А если нет, то берите готовое.