Ответ 1
Фраза "стабильность" относится к алгоритму. Если ваш алгоритм нестабилен, чтобы начать с того, что увеличение точности или уменьшение ошибки округления шагов компонента не будет многого.
Более сложные процедуры numpy, такие как "решать", являются оболочками для подпрограмм ATLAS/BLAS/LAPACK. Вы можете ссылаться на документацию там, например, "dgesv" решает систему вещественных линейных уравнений с использованием декомпозиции LU с частичной разворачиванием и перестановкой строк: базовые документы документа Fortran для LAPACK можно увидеть здесь http://www.netlib.org/lapack/explore-html/, но http://docs.scipy.org/doc/numpy/user/install.html указывает, что доступно множество различных версий стандартных реализаций - оптимизация скорости и точность будут варьироваться между ними.
В ваших примерах нет большого округления, "+" не имеет ненужного округления, точность зависит исключительно от округления неявного в типе данных с плавающей точкой, когда меньшее число имеет младшие биты, которые не могут быть представлены в ответе. Сумма и точка зависят только от порядка оценки. Cumsum не может быть легко перенастроен, поскольку он выводит массив.
Для кумулятивного округления во время функции "cumsum" или "dot" у вас есть выбор:
В Linux 64-разрядный numpy обеспечивает доступ к высокоточному "длинному двойному" типу float128, который можно использовать для снижения потери точности промежуточных вычислений за счет производительности и памяти.
Однако на моем Win64 установите "numpy.longdouble" карты на "numpy.float64" для обычного C двойного типа, поэтому ваш код не является межплатформенным, проверьте "finfo". (На платформе Canopy Express Win64 нет ни float96, ни float128 с действительно более высокой точностью)
log2(finfo(float64).resolution)
> -49.828921423310433
actually 53-bits of mantissa internally ~ 16 significant decimal figures
log2(finfo(float32).resolution)
> -19.931568 # ~ only 7 meaningful digits
Так как sum() и dot() уменьшают массив до одного значения, максимизация точности выполняется со встроенными функциями:
x = arange(1, 1000000, dtype = float32)
y = map(lambda z : float32(1.0/z), arange(1, 1000000))
sum(x) # 4.9994036e+11
sum(x, dtype = float64) # 499999500000.0
sum(y) # 14.357357
sum(y, dtype = float64) # 14.392725788474309
dot(x,y) # 999999.0
einsum('i,i', x, y) # * dot product = 999999.0
einsum('i,i', x, y, dtype = float64) # 999999.00003965141
- обратите внимание на однократное округление округления в пределах "точки" в этом случае, поскольку каждое почти целое число округляется до точного целого числа
Оптимизация округления зависит от того, что вы суммируете - добавление многих небольших чисел сначала может помочь отложить округление, но не избежать проблем, когда существуют большие числа, но отменяют друг друга, поскольку промежуточные вычисления по-прежнему вызывают потерю точности
пример, показывающий зависимость порядка заказа...
x = array([ 1., 2e-15, 8e-15 , -0.7, -0.3], dtype=float32)
# evaluates to
# array([ 1.00000000e+00, 2.00000001e-15, 8.00000003e-15,
# -6.99999988e-01, -3.00000012e-01], dtype=float32)
sum(x) # 0
sum(x,dtype=float64) # 9.9920072216264089e-15
sum(random.permutation(x)) # gives 9.9999998e-15 / 2e-15 / 0.0