Погрешность: где метод Гаусса теряет цифры

Если ты объяснишь, почему 0.1 + 0.2 не равно 0.3, как маленький ведущий элемент уничтожает ответ и почему малая невязка не гарантирует верного решения, - листай до «Теперь сам».

Как думаешь, чему равен x

0,00000000000000001 · x + 1 · y = 1
                  1 · x + 1 · y = 2

Коэффициент при x в первом уравнении - 10⁻¹⁷, почти ноль. Если бы там был ноль, получилось бы y = 1 и x = 1. С таким крошечным числом ответ должен измениться на столько же: x и y остаются практически единицами.

Прежде чем читать дальше, ответь себе: что напечатает функция gauss из прошлого урока, если дать ей эту систему как есть? Деления на ноль тут нет.

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

Представь весы, которые показывают вес с точностью до грамма, и на них стоит слон. Положи ему на спину монету - показания не изменятся: монета меньше одного деления шкалы. Сними слона - на весах ноль, а не монета. Сведения о монете не просто округлились, они исчезли в момент, когда её взвесили вместе со слоном.

Дробное число - такие весы: шкала у него относительная, около шестнадцати значащих цифр. Маленькое число, прибавленное к огромному, пропадает целиком.

Механизм: шестнадцать цифр и что с ними происходит

Дробное число хранится в 64 битах: знак, порядок и 53 двоичных разряда значащей части. Это около 16 десятичных цифр при любом порядке величины. Чисел, которые так представляются точно, конечное множество, и результат каждого действия округляется до ближайшего из них.

0.1 + 0.2                    0.30000000000000004       0,1 не представима точно
десять раз по 0.1            0.9999999999999999        не 1.0
10¹⁶ + 1                     10¹⁶                      единица меньше шага шкалы
(10¹⁶ + 1) − 10¹⁶            0.0                       монета исчезла
(1 + 10⁻¹⁵) − 1              1.11 · 10⁻¹⁵              вместо 10⁻¹⁵: ошибка 11 %

Точно в двоичной записи представимы только суммы долей 1/2, 1/4, 1/8 и так далее. Одна десятая такой конечной суммой не является - как одна треть не записывается конечной десятичной дробью, - и хранится округлённой.

Машинный эпсилон - наименьшая относительная прибавка, которую число ещё замечает. Его можно измерить: делить x пополам, пока 1.0 + x отличается от 1.0. Цикл делает 53 витка и останавливается на x ≈ 1,1 · 10⁻¹⁶ - первой прибавке, которую единица уже не замечает; последняя замеченная вдвое больше, 2,2 · 10⁻¹⁶. Эпсилоном называют то одно, то другое число - порядок у них один.

Последняя строка таблицы - главная. Сами по себе слагаемые хранятся с шестнадцатью верными цифрами. Но при вычитании почти равных чисел старшие цифры сокращаются, и в ответе остаются только младшие - те, что уже были округлены. Это называют потерей значащих цифр: относительная ошибка 10⁻¹⁶ превратилась в 10⁻¹.

Теперь видно, что случилось в системе из начала. Множитель для второй строки - 1 / 10⁻¹⁷ = 10¹⁷:

вторая строка:   1 − 10¹⁷ · 1  =  −10¹⁷          единица исчезла рядом с 10¹⁷
правая часть:    2 − 10¹⁷ · 1  =  −10¹⁷          двойка тоже
y = −10¹⁷ / −10¹⁷ = 1
x = (1 − 1 · y) / 10⁻¹⁷ = 0 / 10⁻¹⁷ = 0          верный ответ: x = 1

Вторая строка после вычитания не помнит ни единицу, ни двойку - всё, чем уравнения отличались. Маленький ведущий элемент дал огромный множитель, а огромный множитель утопил вторую строку в первой.

Правило

Перед обработкой столбца на диагональ переставляют строку с наибольшим по модулю элементом в этом столбце - выбор ведущего элемента. Тогда все множители по модулю не больше единицы, и строки не тонут друг в друге. Дробные числа сравнивают не на равенство, а с допуском. Малая невязка - разность Ax − b - не доказывает, что x близко к точному решению: у плохо обусловленной системы ошибка решения может превышать ошибку данных в миллиарды раз.

Разбор: выбор ведущего элемента

В функцию из прошлого урока добавляются поиск наибольшего элемента в столбце и обмен строк - сразу после строки while i < n:, перед j = i + 1:

        luchshaya = i
        r = i + 1
        while r < n:
            if abs(a[r][i]) > abs(a[luchshaya][i]):
                luchshaya = r
            r = r + 1
        if luchshaya != i:
            t = a[i]
            a[i] = a[luchshaya]
            a[luchshaya] = t
            t = b[i]
            b[i] = b[luchshaya]
            b[luchshaya] = t

Обмен строк - обмен двух ссылок в списке, одно действие, а не n. Та же система, оба варианта:

ведущий элемент      без выбора                        с выбором
10⁻¹⁷                [0.0, 1.0]                        [1.0, 1.0]
10⁻¹⁰                [1.000000082740371, 0.9999999999]  [1.0000000001, 0.9999999999]
0                    деление на ноль                   [1.0, 1.0]

Ответ на вопрос из начала: без выбора функция печатает x = 0. Ни ошибки, ни предупреждения - просто неверное число. Вторая строка показывает, что беда не в одном экзотическом значении: при ведущем элементе 10⁻¹⁰ в ответе осталось семь верных цифр из шестнадцати. Чем меньше ведущий элемент, тем больше цифр теряется.

С выбором на диагональ встаёт единица из второго уравнения, множитель равен 10⁻¹⁷, и теперь первая строка почти не меняет вторую - как и должно быть. Цена выбора - один проход по столбцу на каждом шаге, n² на всё: на системе из 80 уравнений 551 050 шагов вместо 541 250.

Разбор: система, которую перестановки не спасают

Матрица Гильберта: в строке i, столбце j стоит 1 / (i + j + 1), счёт с нуля.

1      1/2    1/3
1/2    1/3    1/4
1/3    1/4    1/5

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

n       ошибка решения        невязка
3       1 · 10⁻¹⁴             2 · 10⁻¹⁶
5       6 · 10⁻¹³             4 · 10⁻¹⁶
7       3 · 10⁻⁸              2 · 10⁻¹⁶
9       2 · 10⁻⁵              4 · 10⁻¹⁶
11      0,01                  4 · 10⁻¹⁶
13      4,2                   9 · 10⁻¹⁶

При n = 13 в ответе вместо единиц стоят числа, которые отличаются от них на несколько единиц. А невязка - шестнадцатый знак: подставь этот неверный ответ в уравнения, и они сойдутся с машинной точностью. Проверка подстановкой говорит «всё верно».

Дело в самой системе. Строки матрицы Гильберта почти пропорциональны друг другу, уравнения почти повторяют одно и то же, и очень разные x дают почти одинаковые b. Округление правой части в шестнадцатом знаке - законная ошибка данных - сдвигает решение на единицы. Такую систему называют плохо обусловленной, а меру этого свойства - числом обусловленности: во сколько раз относительная ошибка решения может превысить относительную ошибку данных. У матрицы Гильберта размера 13 оно около 10¹⁸ - больше, чем шестнадцать значащих цифр способны покрыть.

Алгоритм тут ни при чём, и лучший алгоритм не поможет. Помогает другая постановка задачи или другая арифметика.

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

Первое: if x == 0.3. Сравнение дробных на равенство почти всегда ошибка. Сравнивают так: abs(x - 0.3) < 1e-9, с допуском по смыслу задачи.

Второе: деньги в дробных числах. Десять раз по 0,1 - не единица. Суммы хранят в целых копейках или в десятичном типе.

Третье: сложение большого с маленькими. Тысяча единиц, прибавленных к 10¹⁶ по одной, не меняют его; те же единицы, сложенные сначала между собой, дают верные 10¹⁶ + 1000. Порядок суммирования влияет на ответ: начинать лучше с маленьких.

Четвёртое: вера в невязку. «Подставил - сошлось» проверяет, что решена какая-то близкая система, а не что найдено решение этой.

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

Библиотечный решатель всегда выбирает ведущий элемент, и вместе с ответом у него можно спросить оценку числа обусловленности - numpy.linalg.cond. Если она порядка 10¹², в ответе верны четыре цифры из шестнадцати, как бы аккуратно ни была написана программа.

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

Машинное обучение сознательно работает с числами одинарной и половинной точности - 32 и 16 бит - и вынуждено следить за тем же: где суммы теряют малые слагаемые и где вычитание съедает цифры.

А в прошлых главах раздела погрешности не было вовсе: целые числа и остатки точны. Быстрое преобразование Фурье из раздела «Разделяй и властвуй» как раз стоит на границе: умножение длинных чисел через него идёт в дробных числах, и длину блоков выбирают так, чтобы накопленная ошибка не дошла до половины единицы.

Теперь сам

1. Что напечатает print(0.1 + 0.2 == 0.3) и почему?
2. К числу 10¹⁶ прибавили 1.0, потом вычли 10¹⁶. Что получилось? А если прибавляли 2.0?
3. В столбце под диагональю стоят числа 0.001, −5 и 2, на диагонали - 0.01. Какую строку поставит на диагональ выбор ведущего элемента и каким будет наибольший по модулю множитель?
4. Решатель вернул ответ, невязка 10⁻¹⁵. Число обусловленности матрицы около 10¹⁰. Сколько верных цифр можно ждать в ответе?