Python, NumPy и точность вычислений

Пример, когда обычной точности недостаточно...

Любой, кто имел дело с вычислениями на компьютере знает, что они производятся в нем с определенной точностью (имеются ввиду несимвольные операции). Наиболее распространены одинарная (32-битная) и удвоенная (64-битная) формы представления вещественных чисел. Как правило, этой точности вполне достаточно для большинства возникающих задач; однако, проблемы с округлением могут поджидать там, где их совсем не ждут...

В чем собственно проявляется конечная точность вычислений?!

Рассмотрим простой пример вычислений с одинарной и удвоенной точностями, выполненный в консоли Python:

>>> import numpy as np
>>> a = np.float32(1.0 / 10 ** 8)
>>> b = np.float32(1.0)
>>> 
>>> print (b + a) * 10 ** 10
10000000000.0
>>> 
>>> 
>>> a = np.float64(1.0 / 10 ** 8)
>>> b = np.float64(1.0)
>>> 
>>> print (b + a) * 10 ** 10
10000000100.0

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

Переход к удвоенной точности позволяет избежать рокового округления при сложении. Дело в том, что операции при одинарной точности выполняются с точностью до 7 (число значимых цифр можно определить в двоичной системе, в десятичной — оно приближенно) знака; когда мы прибавляем к единице число меньшее \(10^{-7}\) (надежней — \(10^{-8}\)) фактически ничего не происходит: единица остается единицей, и, как результат, мы получаем вычисления с ошибкой. В случае использования удвоенной точности, все проходит нормально: здесь вычисления производятся до 15 значимых десятичных цифр. Библиотека NumPy здесь используется для задания точности вычислений.

Об одном моем опыте с проблемой недостаточной точности  можно прочесть здесь.

Рассмотрим еще один очень поучительный пример, когда недостаточной является не только удвоенная точность, но и так называемая "учетверенная" (128-битная) точность.

Пусть нам необходимо вычислять значения простой функции следующего вида $$f(a,b) = (333.75 - a^2)  b ^6 + a ^2 (11 a^2 b^2 - 121  b ^4  - 2) + 5.5 b ^ 8 + \frac{a}{2.0 b}$$

Реализуем процесс вычислений этой функции при различных точностях на Python:

# coding: utf-8

import numpy as np

def f(a,b):
    '''
    Взято отсюда: 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)

# -------------- Пробуем различные точности вычислений ---------------
# одинарная точность
a = np.float32(77617)
b = np.float32(33096)

print 'При вычислении с одинарной (float32) точностью:', f(a, b)


# удвоенная точность
a = np.float64(77617)
b = np.float64(33096)

print 'При вычислении с удвоенной (float64) точностью:', f(a, b)


# "учетверенная" точность (IEEE 754 quadruple-precision floating-point format)
a = np.float128(77617)
b = np.float128(33096)

print 'При вычислении с расширенной (float128) точностью:', f(a, b)

Результаты выполнения фрагмента кода будут следующие:

При вычислении с одинарной (float32) точностью: -2.08953762157e+29
При вычислении с удвоенной (float64) точностью: 1.17260394005
При вычислении с расширенной (float128) точностью: 5.76460752303e+17

В документации к C-XSC (библиотеки C++ для интервальных с расширенной точностью) утверждалось, что такие расчеты должны давать одинаковый (1.1726....) результат (почему-то это не удалось воспроизвести на Python с NumPy). Важно, однако, то, что правильный результат совершенно другой, и не превосходит по модулю единицы: -0,827...

Таким образом, казалось бы простая функция... и даже расширенная удвоенная точность является недостаточной! Как же поведет себя результат, если последовательно увеличивать точность вычислений?

Эксперимент с последовательным увеличением точности удобно провести с использованием библиотеки mpmath. Она написана на чистом Python (но может использовать gmpy, если последний установлен).

from mpmath import *

for j in range(6, 50):
    mp.dps = j
    a = mpf(77617)
    b = mpf(33096)
    print 'Decimal precision={0}: f(a, b)={1}'.format(mp.dps, f(a,b))
    

В результате получим следующее:
 

Decimal precision=6: f(a, b)=1.1726
Decimal precision=7: f(a, b)=-7.922816e+28
Decimal precision=8: f(a, b)=1.1726039
Decimal precision=9: f(a, b)=1.17260394
Decimal precision=10: f(a, b)=-7.737125246e+25
Decimal precision=11: f(a, b)=9.6714065569e+24
Decimal precision=12: f(a, b)=1.17260394005
Decimal precision=13: f(a, b)=1.172603940053
Decimal precision=14: f(a, b)=1.1726039400532
Decimal precision=15: f(a, b)=1.17260394005318
Decimal precision=16: f(a, b)=1.172603940053179
Decimal precision=17: f(a, b)=-9.2233720368547758e+18
Decimal precision=18: f(a, b)=-1.15292150460684697e+18
Decimal precision=19: f(a, b)=144115188075855873.2
Decimal precision=20: f(a, b)=1.1726039400531786319
Decimal precision=21: f(a, b)=1125899906842625.1726
Decimal precision=22: f(a, b)=-140737488355326.8273961
Decimal precision=23: f(a, b)=8796093022209.1726039401
Decimal precision=24: f(a, b)=1.17260394005317863185883
Decimal precision=25: f(a, b)=-137438953470.8273960599468
Decimal precision=26: f(a, b)=1.1726039400531786318588349
Decimal precision=27: f(a, b)=1073741825.17260394005317863
Decimal precision=28: f(a, b)=1.172603940053178631858834905
Decimal precision=29: f(a, b)=1.1726039400531786318588349045
Decimal precision=30: f(a, b)=1.17260394005317863185883490452
Decimal precision=31: f(a, b)=1.17260394005317863185883490452
Decimal precision=32: f(a, b)=1.1726039400531786318588349045202
Decimal precision=33: f(a, b)=1.17260394005317863185883490452018
Decimal precision=34: f(a, b)=1.172603940053178631858834904520184
Decimal precision=35: f(a, b)=1.1726039400531786318588349045201837
Decimal precision=36: f(a, b)=-0.827396059946821368141165095479816292
Decimal precision=37: f(a, b)=-0.827396059946821368141165095479816292
Decimal precision=38: f(a, b)=-0.827396059946821368141165095479816292
Decimal precision=39: f(a, b)=-0.827396059946821368141165095479816291999
Decimal precision=40: f(a, b)=-0.827396059946821368141165095479816291999
Decimal precision=41: f(a, b)=-0.82739605994682136814116509547981629199903
Decimal precision=42: f(a, b)=-0.827396059946821368141165095479816291999033
Decimal precision=43: f(a, b)=-0.8273960599468213681411650954798162919990331
Decimal precision=44: f(a, b)=-0.82739605994682136814116509547981629199903312
Decimal precision=45: f(a, b)=-0.827396059946821368141165095479816291999033116
Decimal precision=46: f(a, b)=-0.8273960599468213681411650954798162919990331158
Decimal precision=47: f(a, b)=-0.82739605994682136814116509547981629199903311579
Decimal precision=48: f(a, b)=-0.827396059946821368141165095479816291999033115784
Decimal precision=49: f(a, b)=-0.8273960599468213681411650954798162919990331157844

Удивительно то, что неправильный результат 1.17260394 проявляет достаточную устойчивость на различных точностях. Это вполне может стать причиной того, что он будет ошибочно принят за верный.  Если, например, выполнив вычисления с 8 десятичными знаками, а потом с 30-ю, мы получили одинаковый в первых цифрах результат — вполне вероятно, что он и есть верный! Пример показывает, что и такого рода "основания" могут быть ошибочными.

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

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

Интервальные вычисления в Python могут быть проведены также с использованием mpmath.

for j in range(6, 50):
    iv.dps = j
    a = iv.mpf(77617)
    b = iv.mpf(33096)
    print 'Decimal precision={0}: f(a, b)={1}'.format(mp.dps, f(a,b))
Decimal precision=6: f(a, b)=[-5.0706024009e+30, 6.3382542101e+30]
Decimal precision=7: f(a, b)=[-3.16912650057e+29, 3.16912654779e+29]
Decimal precision=8: f(a, b)=[-5.94211218857e+28, 2.971056097974e+28]
Decimal precision=9: f(a, b)=[-3.7138201178561e+27, 4.9517601582944e+27]
Decimal precision=10: f(a, b)=[-1.54742504910673e+26, 3.09485009825849e+26]
Decimal precision=11: f(a, b)=[-1.934281311383407e+25, 3.86856262277385e+25]
Decimal precision=12: f(a, b)=[-4.8357032784585167e+24, 3.6267774588444373e+24]
Decimal precision=13: f(a, b)=[-3.02231454903657294e+23, 2.26673591177745118e+23]
Decimal precision=14: f(a, b)=[-2.833419889721787128e+22, 2.833419889721790484e+22]
Decimal precision=15: f(a, b)=[-3.5417748621522339103e+21, 3.5417748621522344346e+21]
Decimal precision=16: f(a, b)=[-442721857769029238784.0, 442721857769029246976.0]
Decimal precision=17: f(a, b)=[-27670116110564327424.0, 27670116110564327456.0]
Decimal precision=18: f(a, b)=[-3458764513820540927.0, 2305843009213693953.5]
Decimal precision=19: f(a, b)=[-432345564227567614.828125, 288230376151711745.179688]
Decimal precision=20: f(a, b)=[-36028797018963966.8274231, 18014398509481985.17260742]
Decimal precision=21: f(a, b)=[-3377699720527870.8273963928, 1125899906842625.172604084]
Decimal precision=22: f(a, b)=[-140737488355326.827396061271, 422212465065985.172603942454]
Decimal precision=23: f(a, b)=[-17592186044414.8273960599472, 17592186044417.17260394006735]
Decimal precision=24: f(a, b)=[-1099511627774.8273960599468637, 3298534883329.17260394005325]
Decimal precision=25: f(a, b)=[-137438953470.827396059946822859, 412316860417.172603940053178917]
Decimal precision=26: f(a, b)=[-17179869182.82739605994682137446, 17179869185.17260394005317863941]
Decimal precision=27: f(a, b)=[-2147483646.8273960599468213681761, 2147483649.1726039400531786320407]
Decimal precision=28: f(a, b)=[-268435454.827396059946821368142245, 268435457.172603940053178631864532]
Decimal precision=29: f(a, b)=[-8388606.827396059946821368141165871, 8388609.172603940053178631858840746]
Decimal precision=30: f(a, b)=[-1048574.8273960599468213681411651475, 1048577.1726039400531786318588349559]
Decimal precision=31: f(a, b)=[-131070.827396059946821368141165095792, 131073.172603940053178631858834907439]
Decimal precision=32: f(a, b)=[-8190.827396059946821368141165095483, 1.172603940053178631858834904520184761]
Decimal precision=33: f(a, b)=[-1022.8273960599468213681411650954798445, 1.1726039400531786318588349045201837979]
Decimal precision=34: f(a, b)=[-126.82739605994682136814116509547981678, 1.17260394005317863185883490452018372567]
Decimal precision=35: f(a, b)=[-6.827396059946821368141165095479816292382, 1.172603940053178631858834904520183709123]
Decimal precision=36: f(a, b)=[-0.82739605994682136814116509547981629200549, -0.82739605994682136814116509547981629181741]
Decimal precision=37: f(a, b)=[-0.827396059946821368141165095479816292005489, -0.827396059946821368141165095479816291981979]
Decimal precision=38: f(a, b)=[-0.8273960599468213681411650954798162919996113, -0.827396059946821368141165095479816291998142]
Decimal precision=39: f(a, b)=[-0.82739605994682136814116509547981629199906033, -0.82739605994682136814116509547981629199887666]
Decimal precision=40: f(a, b)=[-0.827396059946821368141165095479816291999037372, -0.827396059946821368141165095479816291999014413]
Decimal precision=41: f(a, b)=[-0.8273960599468213681411650954798162919990345025, -0.8273960599468213681411650954798162919990330676]
Decimal precision=42: f(a, b)=[-0.82739605994682136814116509547981629199903324694, -0.82739605994682136814116509547981629199903306757]
Decimal precision=43: f(a, b)=[-0.827396059946821368141165095479816291999033134834, -0.827396059946821368141165095479816291999033112413]
Decimal precision=44: f(a, b)=[-0.8273960599468213681411650954798162919990331180186, -0.827396059946821368141165095479816291999033115216]
Decimal precision=45: f(a, b)=[-0.82739605994682136814116509547981629199903311591666, -0.82739605994682136814116509547981629199903311574149]
Decimal precision=46: f(a, b)=[-0.827396059946821368141165095479816291999033115785285, -0.827396059946821368141165095479816291999033115763389]
Decimal precision=47: f(a, b)=[-0.8273960599468213681411650954798162919990331157852846, -0.8273960599468213681411650954798162919990331157825476]
Decimal precision=48: f(a, b)=[-0.82739605994682136814116509547981629199903311578442927, -0.82739605994682136814116509547981629199903311578425821]
Decimal precision=49: f(a, b)=[-0.827396059946821368141165095479816291999033115784386505, -0.827396059946821368141165095479816291999033115784365123]

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

Нужно ли теперь всегда проводить интервальный анализ, даже в простых задачах? Во-первых, интервальный анализ выполнить не всегда так просто, как в приведенном примере. Во-вторых, ожидать нежелательных эффектов от влияния ошибок округления по большей части можно только в том случае, когда мы оперируем сильно разнесенными по разрядной сетки числами. Например,  \(1, 10^{-10}, 1000\) — являются "сильно разнесенными" для одинарной точности (поскольку сложение \(1 + 10^{-10} = 1\) для такой машины), следовательно использовать одинарную точность в этом случае не следует. С другой стороны, \(10^{-8}, 10^{-9}, 10^{-10}\) не являются "сильно разнесенными" (для одинарной точности), их порядки отличны не более, чем на количество десятичных знаков (7 десятичных знаков для одинарной точности).

blog comments powered by Disqus