LU-разложение, определитель и точная арифметика

Если ты объяснишь, что хранят матрицы L и U, почему вторая правая часть обходится в n², как из разложения получить определитель и чем Гаусс по модулю отличается от Гаусса в дробных числах, - листай до «Теперь сам».

Как думаешь, сколько стоит вторая правая часть

Мост рассчитывают под сотню вариантов нагрузки. Конструкция одна - значит матрица системы одна и та же, меняется только правая часть. Метод Гаусса из первого урока главы решает одну систему из сорока уравнений за 71 430 шагов.

Прежде чем читать дальше, ответь себе: обязательно ли платить столько же за каждую из ста правых частей? Что в работе метода Гаусса зависит от правой части, а что нет?

На что это похоже

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

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

Механизм: множители на месте нулей

Прямой ход для матрицы из первого урока использовал три множителя: 2 и 4 для первого столбца, 3 для второго. После вычитания на их местах в матрице стоят нули, которые больше не нужны. Туда множители и записывают:

исходная матрица A        после прямого хода, множители на месте нулей
2  1  1                   2  1  1
4  3  3                   2  1  1            ← 2 - множитель, дальше строка U
8  7  9                   4  3  2            ← 4 и 3 - множители, 2 - из U

L - множители и единицы на диагонали      U - то, что получилось над диагональю и на ней
1  0  0                                   2  1  1
2  1  0                                   0  1  1
4  3  1                                   0  0  2

L - нижняя треугольная матрица, U - верхняя; от английских lower и upper. Их произведение равно исходной матрице: A = L · U. Проверить можно прямым умножением: третья строка L на столбцы U даёт 4 · 2 = 8, 4 · 1 + 3 · 1 = 7, 4 · 1 + 3 · 1 + 1 · 2 = 9.

Система A x = b превращается в две треугольные: сначала L y = b решается сверху вниз, потом U x = y - снизу вверх. Первая - это те самые вычитания, которые прямой ход делал с правой частью; вторая - обычный обратный ход. Обе стоят n².

Определитель матрицы - число, которое равно нулю ровно тогда, когда у системы нет единственного решения. У треугольной матрицы он равен произведению диагонали. Вычитание строк определитель не меняет, поэтому определитель A - произведение ведущих элементов: 2 · 1 · 2 = 4. Каждая перестановка строк меняет его знак.

Правило

LU-разложение записывает матрицу как произведение нижней треугольной L (множители метода Гаусса, единицы на диагонали) и верхней треугольной U (результат прямого хода). Разложение стоит Θ(n³) и делается один раз; решение системы для каждой правой части - две треугольные системы, Θ(n²). С выбором ведущего элемента дополнительно хранится перестановка строк. Определитель равен произведению диагонали U, умноженному на −1 в степени числа перестановок. Обратную матрицу для решения систем не вычисляют: разложение даёт то же дешевле и точнее.

Разбор: программа и цена повторных решений

def lu(a):
    n = len(a)
    i = 0
    while i < n:
        j = i + 1
        while j < n:
            k = a[j][i] / a[i][i]
            a[j][i] = k
            c = i + 1
            while c < n:
                a[j][c] = a[j][c] - k * a[i][c]
                c = c + 1
            j = j + 1
        i = i + 1
    return a

def reshi_lu(a, b):
    n = len(a)
    y = [0.0] * n
    i = 0
    while i < n:
        s = b[i]
        c = 0
        while c < i:
            s = s - a[i][c] * y[c]
            c = c + 1
        y[i] = s
        i = i + 1
    return obratnyj_hod(a, y)

От gauss функция lu отличается тремя вещами: строка a[j][i] = k сохраняет множитель, вычитание начинается со столбца i + 1 (столбец i занят множителем), и правой части в ней нет вовсе. reshi_lu сначала проводит правую часть через сохранённые множители, потом вызывает обратный ход из первого урока. Перестановки строк здесь опущены, чтобы было видно главное.

Для матрицы из механизма lu возвращает [[2, 1, 1], [2.0, 1.0, 1.0], [4.0, 3.0, 2.0]]. Правая часть [7, 19, 49] даёт [1.0, 2.0, 3.0], а [4, 10, 24] - [1.0, 1.0, 1.0], уже без нового разложения.

n        разложение     одно решение по готовому      метод Гаусса с нуля
10          1 169                431                       1 515
20          8 634              1 451                       9 920
40         66 464              5 291                      71 430

Разложение стоит почти столько же, сколько метод Гаусса целиком, - куб сидит в нём. Зато каждая следующая правая часть при n = 40 обходится в 5 291 шаг вместо 71 430, в тринадцать с половиной раз дешевле, и разрыв растёт вместе с n. Сто нагрузок на мост: 66 464 + 100 · 5 291, около 600 тысяч шагов вместо семи миллионов.

Разбор: тот же алгоритм без погрешности

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

Первый - обыкновенные дроби: числитель и знаменатель длинными целыми. Ответ точен, но числители и знаменатели растут с каждым шагом исключения, и арифметика из первой главы раздела - список разрядов, умножение за квадрат - превращает n³ во что-то гораздо большее.

Второй - считать по простому модулю p. Все действия метода Гаусса - сложение, вычитание, умножение и деление, - а по простому модулю все четыре есть: деление - это умножение на обратный элемент. Множитель a[j][i] / a[i][i] превращается в a[j][i] · a[i][i]⁻¹ mod p, выбор ведущего элемента - в поиск любого ненулевого. Погрешности нет: остатки точны.

та же система 3 × 3 по модулю 1 000 003:         решение [1, 2, 3]
система  3x + y = 1,  x + 2y = 1:
    в дробных числах                              [0.19999999999999998, 0.4]
    по модулю 1 000 003                           [600002, 200001]

Во второй строке ответ - не целые числа, а 1/5 и 2/5. Дробные числа выдали первое с ошибкой в последнем знаке. По модулю те же дроби выглядят как 600 002 и 200 001 (программа - та же gauss, где деление заменено умножением на pow(a[i][i], -1, p), а после каждого действия берётся остаток): это обратный к пяти, 5 · 600 002 = 3 000 010 = 3 · 1 000 003 + 1, и он же, умноженный на два. Число не похоже на дробь, но ведёт себя как она: умножь на пять - получишь единицу.

Граница приёма: ответ по модулю - это остаток. Если настоящее решение - целые числа меньше p, они получены точно. Если дроби или большие числа - их восстанавливают по нескольким модулям китайской теоремой. Так считают точные определители: по модулям размером в машинное слово, без единого длинного числа до самой сборки.

Что ломается без этого

Первое: обратная матрица ради решения системы. Формула x = A⁻¹ b верна, но вычислять A⁻¹ - это решить n систем с правыми частями из одной единицы и нулей - втрое дороже разложения, и ошибок округления больше. Обратную матрицу считают, только когда она сама нужна.

Второе: определитель по школьной формуле. Разложение по строке порождает n! слагаемых: для матрицы 20 на 20 это 2,4 · 10¹⁸. Через исключение - n³.

Третье: определитель как мера «почти вырожденности». Умножь все уравнения на 0,1 - определитель матрицы 100 на 100 уменьшится в 10¹⁰⁰ раз, а система останется той же. У матрицы Гильберта порядка 10 определитель около 10⁻⁵³, у единичной матрицы того же размера (единицы на диагонали, нули вокруг), умноженной на 0,1, - 10⁻¹⁰; плохо обусловлена только первая. Про близость к вырожденности говорит число обусловленности.

Четвёртое: забытый знак. Определитель после перестановки строк меняет знак. Для матрицы со строками 0 2 и 3 4 он равен −6, а произведение ведущих элементов после обмена строк - 3 · 2 = 6.

Зачем это в работе

scipy.linalg.lu_factor и lu_solve - ровно эта пара: разложить один раз, решать много. Любой расчёт с меняющейся нагрузкой, шаг по времени в моделировании, итерационное уточнение ответа пользуются ею.

Умножение матриц стоит те же n³, и ускорения из раздела «Разделяй и властвуй» - Штрассен и его наследники - переносятся на решение систем: разложение можно свести к умножению. На практике выигрывают не они, а блочные версии обычного алгоритма, дружащие с кешем.

Исключение по модулю 2 - отдельный жанр. Сложение по модулю 2 - это исключающее «или», и система превращается в задачу о битах: какие выключатели нажать, чтобы погасли все лампы; как исправить стёртые биты в коде с избыточностью. Строка матрицы помещается в машинное слово, и один шаг исключения обрабатывает 64 неизвестных разом.

Теперь сам

1. Для матрицы со строками 2 1 и 6 4 найди множитель, L, U и определитель.
2. По этому разложению реши систему с правой частью [3, 10].
3. Матрицу 40 на 40 разложили за 66 464 шага. Сколько шагов уйдёт на 50 правых частей вместе с разложением и сколько - у метода Гаусса с нуля?
4. Решение системы по модулю 7 получилось x = 4. Известно, что точное решение - дробь со знаменателем 2 и числителем от 1 до 6. Какая?