Numpy произвольная точность линейной алгебры

Ответ 1

SymPy может вычислять произвольную точность:

from sympy import exp, N, S
from sympy.matrices import Matrix

data = [[S("-800.21"),S("-600.00")],[S("-600.00"),S("-1000.48")]]
m = Matrix(data)
ex = m.applyfunc(exp).applyfunc(lambda x:N(x, 100))
vecs = ex.eigenvects()
print vecs[0][0] # eigen value
print vecs[1][0] # eigen value
print vecs[0][2] # eigen vect
print vecs[1][2] # eigen vect

выход:

-2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807463648841185335e-261
2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807466621962539464e-261
[[-0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999994391176386872]
[                                                                                                      1]]
[[1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000560882361313]
[                                                                                                    1]]

вы можете изменить 100 в N (x, 100) на другую точность, но, как я пробовал 1000, вычисление собственного vect не удалось.

Ответ 2

В 64-битных системах существует numpy.float128 dtype. (Я полагаю, что там есть и float96 dtype на 32-битных системах). Хотя numpy.linalg.eig не поддерживает 128-битные поплавки, scipy.linalg.eig (сорт) делает.

Однако, в конечном счете, ничто из этого не имеет значения. Любой общий решатель для проблемы с собственным значением будет итеративным, а не точным, поэтому вы ничего не набираете, сохраняя дополнительную точность! np.linalg.eig работает для любой формы, но никогда не возвращает точное решение.

Если вы всегда решаете матрицы 2x2, тривиально написать собственный решатель, который должен быть более точным. Я покажу пример этого в конце...

Независимо, продвигаясь вперед в бессмысленно точные контейнеры памяти:

import numpy as np
import scipy as sp
import scipy.linalg

a = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)
ex = np.exp(a)
print ex

eigvals, eigvecs = sp.linalg.eig(ex)

# And to test...
check1 = ex.dot(eigvecs[:,0])
check2 = eigvals[0] * eigvecs[:,0]
print 'Checking accuracy..'
print check1, check2
print (check1 - check2).dot(check1 - check2), '<-- Should be zero'

Однако вы заметите, что то, что вы получаете, идентично просто выполнению np.linalg.eig(ex.astype(np.float64). На самом деле, я уверен, что делает scipy, а numpy вызывает ошибку, а не делает это молча. Я мог бы ошибаться, хотя...

Если вы не хотите использовать scipy, одним из способов обхода является перемасштабирование вещей после экспоненциальности, но прежде чем решать для собственных значений, отбросить их как "нормальные" поплавки, решить для собственных значений, а затем переделать вещи как float128 после и отмасштабировать.

например.

import numpy as np

a = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)
ex = np.exp(a)
factor = 1e300
ex_rescaled = (ex * factor).astype(np.float64)

eigvals, eigvecs = np.linalg.eig(ex_rescaled)
eigvals = eigvals.astype(np.float128) / factor

# And to test...
check1 = ex.dot(eigvecs[:,0])
check2 = eigvals[0] * eigvecs[:,0]
print 'Checking accuracy..'
print check1, check2
print (check1 - check2).dot(check1 - check2), '<-- Should be zero'

Наконец, если вы решаете только матрицы 2x2 или 3x3, вы можете написать свой собственный решатель, который вернет точное значение для этих форм матриц.

import numpy as np

def quadratic(a,b,c):
    sqrt_part = np.lib.scimath.sqrt(b**2 - 4*a*c)
    root1 = (-b + sqrt_part) / (2 * a)
    root2 = (-b - sqrt_part) / (2 * a)
    return root1, root2

def eigvals(matrix_2x2):
    vals = np.zeros(2, dtype=matrix_2x2.dtype)
    a,b,c,d = matrix_2x2.flatten()
    vals[:] = quadratic(1.0, -(a+d), (a*d-b*c))
    return vals

def eigvecs(matrix_2x2, vals):
    a,b,c,d = matrix_2x2.flatten()
    vecs = np.zeros_like(matrix_2x2)
    if (b == 0.0) and (c == 0.0):
        vecs[0,0], vecs[1,1] = 1.0, 1.0
    elif c != 0.0:
        vecs[0,:] = vals - d
        vecs[1,:] = c
    elif b != 0:
        vecs[0,:] = b
        vecs[1,:] = vals - a
    return vecs

def eig_2x2(matrix_2x2):
    vals = eigvals(matrix_2x2)
    vecs = eigvecs(matrix_2x2, vals)
    return vals, vecs

a = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)
ex = np.exp(a)
eigvals, eigvecs =  eig_2x2(ex) 

# And to test...
check1 = ex.dot(eigvecs[:,0])
check2 = eigvals[0] * eigvecs[:,0]
print 'Checking accuracy..'
print check1, check2
print (check1 - check2).dot(check1 - check2), '<-- Should be zero'

Это возвращает действительно точное решение, но будет работать только для 2x2-матриц. Это единственное решение, которое действительно выигрывает от дополнительной точности, однако!

Ответ 3

Насколько мне известно, numpy не поддерживает более высокую, чем двойную точность (float64), которая по умолчанию не указана.

Попробуйте использовать это: http://code.google.com/p/mpmath/

Список функций (среди прочих)

Арифметика:

• Реальные и комплексные числа с произвольной точностью
• Неограниченные размеры/величины экспоненты

Ответ 4

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

Удачи!

Ответ 5

Итак, вам нужно 350 цифр точности. Вы не получите этого с номерами с плавающей запятой IEEE (что используется numpy). Вы можете получить его с помощью программы bc:

$ bc -l
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
scale=350
e(-800)
.<hundreds of zeros>00366