Вычисление определителей и обращение матриц: численные аспекты
Разрыв между теоретической алгеброй и вычислительной реальностью
В курсе высшей алгебры студенты учатся решать системы линейных алгебраических уравнений (СЛАУ) вида Ax = b с помощью правила Крамера (через определители) или путем явного вычисления обратной матрицы (x = A^-1 * b). Эти методы красивы в теории и отлично подходят для ручного решения простейших систем 3x3 на листке бумаги. Однако в реальной вычислительной математике при написании программного обеспечения инженеры избегают этих методов как огня. Почему же аналитические фавориты терпят крах при встрече с кремниевыми процессорами?
Все дело в вычислительной сложности. Вычисление определителя по классической формуле разложения по строке (формула Лейбница) требует факториального количества операций (O(N!)). Для матрицы размером всего 20x20 вычисление определителя таким способом потребует 20! операций. Если суперкомпьютер будет выполнять миллиард операций в секунду, ему потребуется несколько сотен тысяч лет, чтобы решить эту крошечную задачу. Использование правила Крамера для СЛАУ требует вычисления N+1 таких определителей, что делает этот метод абсолютно бесполезным для практических расчетов в физике и инженерии.
Почему обращать матрицы явно — плохая идея?
Вторая распространенная ошибка новичков в программировании — явное вычисление обратной матрицы A^-1 для решения системы уравнений. Во-первых, вычисление обратной матрицы методом Гаусса-Жордана требует в 3 раза больше арифметических операций, чем решение самой системы методом исключения Гаусса (или с помощью LU-разложения). Мы тратим драгоценное время процессора впустую.
Во-вторых, обращение разреженных матриц приводит к катастрофическому заполнению памяти. Как мы помним из предыдущих статей, в методе конечных элементов матрица A на 99% состоит из нулей. Однако ее обратная матрица A^-1 почти всегда является полностью плотной (все элементы не равны нулю). Если исходная матрица из миллиона уравнений в формате CSR занимала несколько десятков мегабайт, то обратная матрица займет терабайты оперативной памяти. Именно поэтому золотое правило вычислительной линейной алгебры гласит: «Никогда не вычисляйте обратную матрицу явно, если вам нужно только решить систему уравнений». Вычисляйте LU-разложение или используйте итерационные методы.
Обусловленность и устойчивость обратной матрицы
Еще одной колоссальной проблемой является накопление ошибок округления. Каждая математическая матрица характеризуется числом обусловленности (Condition Number) — грубо говоря, это мера того, насколько матрица близка к вырожденной (к матрице с нулевым определителем).
Если число обусловленности велико (матрица плохо обусловлена), то микроскопические ошибки в исходных данных или ошибки округления с плавающей запятой внутри процессора при вычислении обратной матрицы будут многократно усилены (иногда в миллионы раз). Полученная компьютерная «обратная матрица» при умножении на исходную даст результат, крайне далекий от единичной матрицы. Анализ собственных и сингулярных чисел (SVD), а также использование алгоритмов предобуславливания (Preconditioning) стали обязательными инструментами для стабилизации процесса решения СЛАУ и борьбы с физической некорректностью инженерных моделей.
