Когда речь идет о повседневных арифметических операциях, проблемы с конечной точностью вычислений не выглядят столь пугающими. И наилучшей проверкой того, что результат получен правильно, является сравнение значений полученных на различных точностях.
Если, например, вычисления, полученные на одинарной и удвоенной точностях совпадают, то создается чувство уверенности в результате, по крайней мере с точностью сопоставимой с одинарной. Здесь, я бы хотел привести один интересный пример, демонстрирующий, что даже в сравнительно несложной арифметической задаче подобная устойчивость при переменной точности представления чисел не может служить основанием для такой уверенности.
Необходимо вычислить значение функции двух переменных при определенных значениях (даны ниже) ее аргументов:
Этот пример был мною замечен, когда я разбирался с библиотекой C-XSC (система классов языка C для (преимущественно) интервальных научных расчетов с произвольной точностью). Данная библиотека отлично подходит для практического исследования вычислительной устойчивости различных численных алгоритмов.
Для эмуляции вычислений с плавающей точкой на Python установим пакет mpmath. В принципе, если ограничиться точностями IEEE 754, можно было бы использовать и NumPy, но нам нужно показать, что результаты, получаемые в рамках IEEE 754 неверны.
Вычислим значение функции f(a, b) при a = 77617 и b = 33096.
# coding: utf-8
from mpmath import *
def f(a,b):
'''
From: Ramp S.M. Algorithms for verified inclusions -- theory and practice, USA, NY, 1988.
'''
return (333.75 - a * a) * b ** 6 + a * a * (11 * a * a * b * b - 121 * b ** 4 - 2) + 5.5 * b ** 8 + a/(2.0*b)
# Одинарная точность
mp.dps = 8
a = mpf(77617)
b = mpf(33096)
print 'Значение f(a, b) с одинарной точностью:', f(a, b)
# Удвоенная точность
mp.dps = 16
a = mpf(77617)
b = mpf(33096)
print 'Значение f(a, b) с удвоенной точностью:', f(a, b)
Читатель может заметить, что задание точности через dps в mpmath это не совсем правильный подход при эмуляции реальной удвоенной точности. Если речь идет об удвоенной и\или одинарной точности в рамках IEEE то, пожалуй, более корректно описывать их характеристики через двоичную систему. Здесь, однако, это не так важно; зато использование mp.dps более проще интерпретируется — т.е. как количество десятичных значащих цифр в представлении числа.
Выполняя код, получим:
Значение f(a, b) с одинарной точностью: 1.1726039
Значение f(a, b) с удвоенной точностью: 1.172603940053179
Вполне убедительно, не так ли? Только значение-то неправильное! Правильное значение при данных a и b вообще меньше единицы!
Дополним пример следующими вычислениями:
for i in range(8, 40):
mp.dps = i
a = mpf(77617)
b = mpf(33096)
print 'Точность dps={0}, f(a,b)={1}'.format(mp.dps, f(a,b))
Получим (некоторые строки опущены):
Точность dps=8, f(a,b)=1.1726039
Точность dps=9, f(a,b)=1.17260394
Точность dps=10, f(a,b)=-7.737125246e+25
...
Точность dps=13, f(a,b)=1.172603940053
Точность dps=14, f(a,b)=1.1726039400532
Точность dps=15, f(a,b)=1.17260394005318
Точность dps=16, f(a,b)=1.172603940053179
Точность dps=17, f(a,b)=-9.2233720368547758e+18
...
Точность dps=28, f(a,b)=1.172603940053178631858834905
Точность dps=29, f(a,b)=1.1726039400531786318588349045
Точность dps=30, f(a,b)=1.17260394005317863185883490452
...
Точность dps=36, f(a,b)=-0.827396059946821368141165095479816292
Точность dps=37, f(a,b)=-0.827396059946821368141165095479816292
Точность dps=38, f(a,b)=-0.827396059946821368141165095479816292
Точность dps=39, f(a,b)=-0.827396059946821368141165095479816291999
Видно, что несмотря на устойчивость, некоторые значения все-таки существенно отличаются от обычных, что уже заставляет задуматься.
Скажу сразу, значения, получаемые при точности dps=36 и выше, правильные. Но откуда знать, что при дальнейшем увеличении точности не произойдет еще какого-либо скачка, ведь, как было видно со значением 1.17260..., даже устойчивость результата при различных точностях не может гарантировать его правильность.
К счастью, и на этот вопрос можно ответить, используя пакет mpmath. Данный пакет позволяет производить интервальные вычисления, что дает возможность оценить область возможных значений функции при представлении ее аргументов с различной точностью.
Выполним расчеты, используя аппарат интервальной арифметики mpmath:
for j in range(6, 40):
iv.dps = j
a = iv.mpf(77617)
b = iv.mpf(33096)
print 'Точность={0}: f(a, b)={1}'.format(mp.dps, f(a,b))
Получим следующее:
Точность dps=6: f(a, b)=[-5.0706024009e+30, 6.3382542101e+30]
Точность dps=7: f(a, b)=[-3.16912650057e+29, 3.16912654779e+29]
Точность dps=8: f(a, b)=[-5.94211218857e+28, 2.971056097974e+28]
Точность dps=9: f(a, b)=[-3.7138201178561e+27, 4.9517601582944e+27]
Точность dps=10: f(a, b)=[-1.54742504910673e+26, 3.09485009825849e+26]
Точность dps=11: f(a, b)=[-1.934281311383407e+25, 3.86856262277385e+25]
Точность dps=12: f(a, b)=[-4.8357032784585167e+24, 3.6267774588444373e+24]
Точность dps=13: f(a, b)=[-3.02231454903657294e+23, 2.26673591177745118e+23]
Точность dps=14: f(a, b)=[-2.833419889721787128e+22, 2.833419889721790484e+22]
Точность dps=15: f(a, b)=[-3.5417748621522339103e+21, 3.5417748621522344346e+21]
Точность dps=16: f(a, b)=[-442721857769029238784.0, 442721857769029246976.0]
Точность dps=17: f(a, b)=[-27670116110564327424.0, 27670116110564327456.0]
Точность dps=18: f(a, b)=[-3458764513820540927.0, 2305843009213693953.5]
Точность dps=19: f(a, b)=[-432345564227567614.828125, 288230376151711745.179688]
Точность dps=20: f(a, b)=[-36028797018963966.8274231, 18014398509481985.17260742]
Точность dps=21: f(a, b)=[-3377699720527870.8273963928, 1125899906842625.172604084]
Точность dps=22: f(a, b)=[-140737488355326.827396061271, 422212465065985.172603942454]
Точность dps=23: f(a, b)=[-17592186044414.8273960599472, 17592186044417.17260394006735]
Точность dps=24: f(a, b)=[-1099511627774.8273960599468637, 3298534883329.17260394005325]
Точность dps=25: f(a, b)=[-137438953470.827396059946822859, 412316860417.172603940053178917]
Точность dps=26: f(a, b)=[-17179869182.82739605994682137446, 17179869185.17260394005317863941]
Точность dps=27: f(a, b)=[-2147483646.8273960599468213681761, 2147483649.1726039400531786320407]
Точность dps=28: f(a, b)=[-268435454.827396059946821368142245, 268435457.172603940053178631864532]
Точность dps=29: f(a, b)=[-8388606.827396059946821368141165871, 8388609.172603940053178631858840746]
Точность dps=30: f(a, b)=[-1048574.8273960599468213681411651475, 1048577.1726039400531786318588349559]
Точность dps=31: f(a, b)=[-131070.827396059946821368141165095792, 131073.172603940053178631858834907439]
Точность dps=32: f(a, b)=[-8190.827396059946821368141165095483, 1.172603940053178631858834904520184761]
Точность dps=33: f(a, b)=[-1022.8273960599468213681411650954798445, 1.1726039400531786318588349045201837979]
Точность dps=34: f(a, b)=[-126.82739605994682136814116509547981678, 1.17260394005317863185883490452018372567]
Точность dps=35: f(a, b)=[-6.827396059946821368141165095479816292382, 1.172603940053178631858834904520183709123]
Точность dps=36: f(a, b)=[-0.82739605994682136814116509547981629200549, -0.82739605994682136814116509547981629181741]
Точность dps=37: f(a, b)=[-0.827396059946821368141165095479816292005489, -0.827396059946821368141165095479816291981979]
Точность dps=38: f(a, b)=[-0.8273960599468213681411650954798162919996113, -0.827396059946821368141165095479816291998142]
Точность dps=39: f(a, b)=[-0.82739605994682136814116509547981629199906033, -0.82739605994682136814116509547981629199887666]
Интересно проследить эволюцию интервала возможных значений функции: он стабилизируется только при использовании точностей от 36 и более значимых десятичных цифр, хотя и постепенно сужается.
Из интервальных расчетов становится вполне ясным, что доверять следует только результатам, получаемым при точностях от dps=36 и более десятичных цифр.
Этот пример — наглядная демонстрация того, насколько нужно быть осторожным осуществляя вычисления с плавающей точкой, и что даже 128-битная (например, numpy.float128 если говорить о Python и NumPy) точность может быть недостаточной. Он также показывает, что нельзя доверять и устойчивости результата, полученного на различных точностях. Применение аппарата интервальных вычислений может быть одним из вариантов решения в этом случае, которое позволяет оценить необходимую точность для получения адекватного результата.
Комментарии (36)
gbg
08.07.2016 12:30+46Делай раз — смотрим на формулу: у вас 3*10^4 возводится в степень 8.
Делай два — результат этого возведения — 3*10^32 — упс, вышли за приделы точности double
Делай три — открываем тетрадочку с лекциями по численным методам и читаем «Для обеспечения устойчивого счета, требуется заложить минимум 4 порядка точности арифметики». Суммируем 32+4=36. Ой, что-то теория работает.
Это я к чему. Вычислительная математика работает строго по принципу GIGO — мусор на входе — мусор на выходе. Вы заранее засунули компьютеру ерунду, которую он обсчитать не может, и удивляетесь, что он выдал вам ерунду. Так математический софт не пишут.
burrdarr
08.07.2016 16:20+2Плюсанусь, однако замечу, что не столь важен максимальный порядок, сколь разница между максимальным и минимальным. Собственно, вангую, что не будь «а/2бе», который таки порядка 10^0, то и устойчивость достигалась бы на меньшей точности, чем 36.
gbg
08.07.2016 16:29+6«Устойчивость» вообще термин здесь малоуместный.
«Устойчивость» означает, что малые изменения параметров системы (обычно — коэффициентов диффура, или граничных условий), ведут к малым изменениям решения.
И этот термин относится к алгоритму счета, там и критерий специальный есть, и теория, которая этот критерий строго доказывает, и так далее.
Здесь же у нас чисто арифметическая задача «сколько брать знаков на расчет» и ответ к ней — давным давно эмпирически предложен — запас должен быть 4 (это как самый минимум!) знака.
Вы совершенно правы, указывая на то, что нужно брать разницу между самой большой и самой меньшей степенями десятки.
Однако, при численном решении реальных задач (в матфизике, например), задачу сначала нормируют (чтобы самое большое число по модулю было равно 1), а затем уже с нормированной и безразмерной задачей работают. Тогда критерий становится простым — «хотим N знаков точности — берем 4+N знака арифметики (а лучше — больше)».scidam
10.07.2016 08:40Прежде всего, f(а, b) — отображение, заданное в виде арифметического выражения, относительно простого, да. Но оно очень чувствительно к изменениям входящих параметров; это — фактически — явление неустойчивости.
gbg
10.07.2016 09:18+2Ваша формула — полином. Это гладкая функция, непрерывная вместе со всеми своими производными. Изменение параметров на бесконечно малую величину в данном случае приведет к изменению результата также на бесконечно малую величину.
В очередной раз повторяю, что термин «устойчивость» относится к конкретному алгоритму вычислений, а не к формуле. Вы же не приводите алгоритм (как конкретно интерпретатор пайтона это все считал? Прямо с точностью до действий?).
Например, для задачи об охлаждении неоднородного стержня, алгоритм прямого счета является условно устойчивым (нужно придерживаться определенного соотношения шагов разностной схемы, иначе на выходе — мусор), а схема Кранка-Николсон (не НиколсонА, Филлис Николсон — женщина!) — абсолютно устойчива и работает для любого соотношения шагов.
dbagaev
08.07.2016 16:34+2Вычислительная математика работает и по другим принципам, иногде получая вполне корректные данные на входе, но выдавая мусор на выходе, и дело часто отнюдь не только в том, что где-то мы вылазим за пределы допустимых значений или теряем разряды. Проблема в том, что некоторые вычисления являются вычислительно нестабильными, например элементарное вычитание двух чисел либо суммирование ряда.
Я не буду сейчас приводить и объяснять конкретные примеры, но в одном из них выполнение серии безобидных операций без выходов за пределы чисел на выходе дает число, в котором все разряды кроме старшего неправильные. Есть и контр-примеры, в которых изменение порядка вычислений сохраняет вычислительную стабильность даже при переполнениях.mbait
08.07.2016 21:54вылазим за пределы допустимых значений или теряем разряды
Есть ещё ошибки округления. Но это вроде бы полный набор. Какие ещё нестабильности вы знаете?gbg
08.07.2016 22:01-ошибки входных данных (неточности измерений физических величин)
-ошибки, вносимые математической моделью (несоответствие между моделью и реальностью)
-погрешность численного метода (неустойчивость из-за «плохой» геометрии, погрешность аппроксимации численной схемы и так далее)
-погрешность арифметики
scidam
10.07.2016 08:28-2Не совсем правильно считать данный пример ерундой; Для иллюстрации влияния погрешностей округления, действительно, можно привести пример гораздо проще, но смысл данного примера в другом: допустим мы решили посчитать данное выражение с одинарной точностью, посчитали; мы понимаем, что здесь степени и возможны влияния погрешностей округления. Для проверки результата мы можем пересчитать все это в удвоенной точности, и получим такой же результат! Вычисляя выражение с учетверенной точностью (все в рамках IEEE 754), снова получаем тот же результат. Так легко и поверить, что он правильный, но он-то оказывается не верным!
Вывод из примера, что получение мало отличающихся результатов при увеличении точности не гарантирует, что результат правильный.gbg
10.07.2016 08:46+1Прежде чем что-либо считать, всегда можно априорно оценить результат.
Более того, всегда можно заставить компьютер работать с символьными выражениями, то есть сохранять все знаки (симулируя старательного школьника, который решает пример).
В очередной раз повторяю — никакие реальные вычисления не делаются без априорных теоретических оценок. Доверять можно только сочетанию из формального доказательства и результата работы программы. Иначе — получите GIGO. В вычислительной практике никто не опирается только на то, что насчитал компьютер.
ARad
08.07.2016 14:07Wolfram Alpha считает что ответ примерно равен -4.27229?10^20
playermet
08.07.2016 20:05+1А мне Wolfram Alpha сказал, что -54767/66192. Это как раз примерно -0.82739.
ARad
09.07.2016 04:31Я вставлял ссылку WA, но ее почемуто удалили ;( не знаю почему у нас разные ответы
playermet
09.07.2016 14:23Я вставлял в форму WA вот такую строку: «a = 77617, b = 33096, c = (333.75 — a * a) * b ** 6 + a * a * (11 * a * a * b * b — 121 * b ** 4 — 2) + 5.5 * b ** 8 + a/(2.0*b)». Возможно это не совсем правильный способ, но по другому я не знаю как. Ответ получил как писал выше.
Еще я пытался вводить строку с уже подставленными значениями a и b: "(333.75 — 77617 * 77617) * 33096 ** 6 + 77617 * 77617 * (11 * 77617 * 77617 * 33096 * 33096 — 121 * 33096 ** 4 — 2) + 5.5 * 33096 ** 8 + 77617/(2.0*33096)". В этом случае WA ответил без дробей, сразу числом из поста.
Dum_spiro_spero
10.07.2016 08:41Мне такое сказал Maple — но не сразу. Пришлось представить 333.75 как 333+3/4, и т.п…
Вообще пример хороший — например в прикладном софте такое запросто можно проморгать. Формула есть? Есть. Выглядит вроде адекватно — никаких мелких знаменателей нет. Восьмая степень — это все ж не двадцать восьмая. Уравнение состояния Бенедикта — Вебба — Рубина содержит шестую степень — и ничего. Так что вопросы откуда взялась такая странная формула возникнут не у всех, а значит «тестирования» формулы может и не быть.gbg
10.07.2016 12:09«Восьмая степень» — при использовании double, это ограничение на аргумент — не стоит допускать, чтобы аргумент под этой 8 степенью был больше 10 по модулю. При использовании float, вообще никакие расчеты не пойдут — результат не вместится.
michael_vostrikov
08.07.2016 15:27+8Ну прям секрет открыли. А еще unsigned int нельзя использовать для вычислений с результатом, который выходит за пределы 2^32.
Если посчитать первый множитель в обычном калькуляторе, то получается результат -7,9171109033773850490791882372801e+36, откуда видно, что для точных вычислений надо как минимум 37 разрядов.
Siper
08.07.2016 22:51-1Плюс-минус 1 — это офигительная точность для того диапазона значений, которые принимает данная функция))
Мы ж не деньги считаем (их в 8ю степень не возводят).
scidam
10.07.2016 08:51-1Благодарю всех за комментарии, но, к сожалению, практически все здесь увидели пример, лишь как обычную иллюстрацию влияния погрешностей округления. Смысл примера — в другом. Мы получаем практически одинаковые неправильные результаты вычислений последовательно повышая точность представления чисел, и это интересно и показательно, т.к. получение мало отличающихся результатов на различных разрядных сетках является «доказательством» его правильности.
michael_vostrikov
10.07.2016 11:55+2Я бы сказал, это смысл комментов в другом. То, что вы повысили точность, не означает, что она стала достаточной. Не знаю, почему вы считаете, что разные разрядные сетки что-то доказывают. Нельзя получить 1 в выражении 1000000000877 — 876 — 1000000000000, если разрядная сетка 1, 2, 4, или 8 знаков, и младшие разряды в первом числе превращаются в 0. Смысл не в количестве разрядов, а в том, помещается ли туда результат.
Поэтому фраза
Но откуда знать, что при дальнейшем увеличении точности не произойдет еще какого-либо скачка
технически абсолютно неграмотна. Оттуда и знать, что точность теперь достаточная. Почему достаточная? Потому что порядок множителей меньше или равен 10^36.scidam
11.07.2016 02:18-1Не знаю, почему вы считаете, что разные разрядные сетки что-то доказывают...
Конечно не доказывают, но могут придавать уверенности, особенно в сложных вычислениях. (реально сталкивался с тем, что люди выполняют «такую проверку», осуществляя расчеты в академических учреждениях...)
Если бы наша задача состояла в вычислении собственных чисел какой-либо несимметричной матрицы, процедура вычислений не была бы столь прозрачна, и было бы очень замачиво «проверить» результат, вычислив его, например, на удвоенной и расширенной удвоенной точностях. И вот, мы проделали вычисления, и результаты «почти» совпали; вдобавок еще провели вычисления на учетвернной точности, получили «такие же» значения… так недолго и поверить, что результат «правильный» (понимаю, что вычисления на различных разрядных сетках — не основание для проверки результатов вычислений собственных значений несим. матриц, для этого есть другие подходы — например, построение их спектральных портретов — псевдоспектров).
Пример показывает, что «проверка» на различных разрядных сетках ничего не значит, т.к. почти одинаковые неправильные значения могут получаться и при увеличении точности.
Про наличие еще «одного скачка» — применительно к примеру, согласен, что фраза не учитывает, что в этом простом выражении можно легко оценить необходимое число разрядов; а если бы процедура вычислений была бы достаточно сложна и получение оценки требуемого числа разрядов было бы проблематично…
gbg
11.07.2016 02:26Если бы наша задача состояла в вычислении собственных чисел какой-либо несимметричной матрицы, процедура вычислений не была бы столь прозрачна, и было бы очень замачиво «проверить» результат, вычислив его, например, на удвоенной и расширенной удвоенной точностях.
Совершенно бессмысленная затея. В ситуации с матрицами и собственными числами нужно оценивать число обусловленности матрицы, а не гонять компьютер почем зря. Более того, для сколько-нибудь «интересной» задачи расчет без аппаратной поддержки арифметики займет чрезвычайно длительное время, а на выходе даст пшик.scidam
12.07.2016 03:23Число обусловленности здесь не особо подходит, например, для вырожденой матрицы число обусловленности (спектральное) есть бесконечность, и что делать дальше, если собственные значения нужно найти (одно из этих значений будет 0, а остальные)… Но здесь не об этом…
В MatLab есть тестовая матрица gallery(5), следующего вида:
А = [[-9, 11, -21, 63, -252], [70, -69, 141, -421, 1684], [-575, 575, -1149, 3451, -13801], [3891, -3891, 7782, -23345, 93365], [1024, -1024, 2048, -6144, 24572]]. Нужно найти ее собственные значения (несимвольно);
Есть даже ее совсем неправильный спетральный портрет, который подталкивает поверить в неправильно вычисленные собственные значения.
Вот и неправильный результат: eigs (A) = ( -0.03697710+0.02749668j, -0.03697710-0.02749668j, 0.01470225+0.04265918j, 0.01470225-0.04265918j, 0.04454971+0.j ) — (вычислено при помощи numpy.linalg.eig); MatLab раньше (R2008a, например) давал также неправильный результат, сейчас — не знаю, т.к. не пользуюсь MatLab давно.
Какой же результат правильный?! — 0.0 (т.е. 5-ти кратное нулевое соб. значение)
Правда этот пример не такой красивый (для игр с точностями), т.к. на различных точностях дает различные наборы собственных значений.
gbg
12.07.2016 08:20+1Будем последовательны. Число обусловленности прекрасно «работает» в случае вырожденной матрицы — оно сообщает нам, что матрица поганая, а следовательно, стандартный метод для нее — не подходит.
В вычислительной математик1е нет одного универсального метода, который бы оптимально подходил для всех задач, к каждому методу жестко прибит гвоздями-сотками критерий его применимости (как и вообще в науке, кстати). Вы в данной публикации критерии применимости систематически игнорируете, в результате закономерно получаете лажу.
gnomeby
11.07.2016 13:50-2Всё правильно вы сделали. В современном мире люди без спец математического образования могут, что-то считать по формулам и не знать ничего про коварные точности.
А все ваши критики в этой статье просто матемаческие снобы, считающие людей с не таким образованием просто более низшей кастой.
buriy
11.07.2016 21:08+2Мне кажется, вы переусложняете без необходимости.
У каждого программиста есть более простой пример перед глазами:
(int8) (2 ** 36)совпадает с
(int16) (2 ** 36)и с
(int32) (2 ** 36).
И что это доказывает? Что 2 ** 36 равно нулю, поскольку все ваши вычисления дали нуль?
Так и здесь: поскольку вы отнимаете друг от друга числа порядка 10**36, то для точного результата мантисса должна вмещать в себя по крайней мере 36 десятичных разрядов. Сколько разрядов вмещает мантисса у double, помните? Всего лишь 15-17.
Это даже в википедии написано: English Wikipedia: Double-precision_floating-point_format
napa3um
Сайт, посвящённый проблеме точности вычислений в разных языках: http://0.30000000000000004.com/