python - escribir - numpy álgebra lineal de precisión arbitraria

python cientifico pdf

Question

Tengo una matriz 2d numpy [tamaño mediano / grande - digamos 500x500]. Quiero encontrar los valores propios del exponente a nivel de elemento de él. El problema es que algunos de los valores son bastante negativos (-800, -1000, etc.) y sus exponentes tienen un flujo inferior (lo que significa que están muy cerca de cero, por lo que numpy los trata como cero). ¿Hay alguna forma de usar precisión arbitraria en numpy?

La forma en que lo sueño:

import numpy as np

np.set_precision(''arbitrary'') # <--- Missing part
a = np.array([[-800.21,-600.00],[-600.00,-1000.48]])
ex = np.exp(a)  ## Currently warns about underflow
eigvals, eigvecs = np.linalg.eig(ex)

He buscado una solución con gmpy y mpmath en vano. Cualquier idea será bienvenida.

Answer 1

Entonces necesitas 350 dígitos de precisión. No lo conseguirás con los números flotantes IEEE (qué numpy está usando). Puede obtenerlo con el programa 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

Answer 2

No tengo experiencia con Numpy en particular, pero agregar un punto decimal con una cantidad establecida de ceros puede ayudar. Por ejemplo, use 1.0000 en lugar de 1. En las secuencias de comandos python normales donde he tenido este problema, esto ha ayudado, así que a menos que su problema sea causado por una rareza en numpy y no tenga nada que ver con python, esto debería ayudar.

¡Buena suerte!

Answer 3

Por lo que yo sé, Numpy no admite una precisión superior a doble (float64), que es el valor predeterminado si no se especifica.

Intenta usar esto: http://code.google.com/p/mpmath/

Lista de características (entre otras)

Aritmética:

• Números reales y complejos con precisión arbitraria
• Ilimitado exponente de tamaños / magnitudes

Answer 4

SymPy puede calcular la precisión arbitraria:

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

salida:

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

puedes cambiar 100 en N (x, 100) a otra precisión, pero, como probé 1000, el cálculo de eigen vect falló.

Answer 5

En los sistemas de 64 bits, hay un numpy.float128 numpy.float128. (También creo que hay un float96 float96 en los sistemas de 32 bits). Mientras que numpy.linalg.eig no admite flotantes de 128 bits, scipy.linalg.eig (más o menos) sí.

Sin embargo, nada de esto va a importar , a la larga. Cualquier solucionador general para un problema de valor propio va a ser iterativo, en lugar de exacto, ¡por lo que no obtendrás nada manteniendo la precisión extra! np.linalg.eig funciona para cualquier forma, pero nunca devuelve una solución exacta .

Si siempre estás resolviendo matrices de 2x2, es trivial escribir tu propio solucionador que debería ser más exacto. Voy a mostrar un ejemplo de esto al final ...

Independientemente, avanzando en contenedores de memoria sin sentido precisos:

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''

Sin embargo, notarás que lo que obtienes es idéntico a simplemente hacer np.linalg.eig(ex.astype(np.float64) . De hecho, estoy bastante seguro de que eso es lo que hace scipy , mientras que numpy plantea un error que hacerlo en silencio. Sin embargo, podría estar bastante equivocado ...

Si no quiere usar scipy, una solución alternativa es reescalar elementos después de la exponenciación pero antes de resolver los valores propios, convertirlos en flotantes "normales", resolver los valores propios y luego refundir cosas como float128 después y reescalar.

P.ej

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''

Finalmente, si solo resuelve matrices de 2x2 o 3x3, puede escribir su propio solucionador, que devolverá un valor exacto para esas formas de matrices.

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''

Esta devuelve una solución verdaderamente exacta, pero solo funcionará para matrices de 2x2. ¡Sin embargo, es la única solución que realmente se beneficia con la precisión adicional!