numpy précision arbitraire d'algèbre linéaire

Sur les systèmes 64 bits, il y a un numpy.float128 dtype. (Je crois qu'il y a un float96 dtype sur les systèmes 32 bits) Alors que numpy.linalg.eig ne prend pas en charge 128 bits flotteurs, scipy.linalg.eig (en quelque sorte).

Cependant, rien de tout cela est d'aller à la matière, dans le long terme. Général solveur pour une valeur propre problème va être itératif, plutôt que d'exact, de sorte que vous ne gagnez pas tout en gardant une précision supplémentaire! np.linalg.eig fonctionne pour n'importe quelle forme, mais ne retourne jamais exacte solution.

Si vous êtes toujours à la résolution des matrices 2x2, il est trivial d'écrire votre propre solveur qui devrait être plus exact. Je vais vous montrer un exemple de ce à la fin...

Peu importe, à aller de l'avant dans inutilement précis de la mémoire des conteneurs:

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'

Cependant, vous remarquerez que ce que vous obtenez est identique à juste faire np.linalg.eig(ex.astype(np.float64). En fait, je suis assez sûr que ce scipy est en train de faire, tandis que numpy génère une erreur plutôt que de le faire en silence. Je pourrais être tout à fait tort, mais...

Si vous ne voulez pas utiliser de scipy, une solution de contournement consiste à redimensionner les choses après l'exponentiation mais avant de résolution pour les valeurs propres, les jeter comme "normale", les chars, les résoudre pour les valeurs propres, et puis refonte des choses comme float128 de suite et de le redimensionner.

E. g.

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'

Enfin, si vous êtes seulement à résoudre 2x2 ou 3x3 matrices, vous pouvez écrire vos propres solveur, qui sera retourner une valeur exacte pour ces formes 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'

Celle-ci rentre vraiment une solution exacte, mais ne fonctionne que pour les matrices 2x2. C'est la seule solution qui en fait des avantages de l'une précision supplémentaire, cependant!

OriginalL'auteur Joe Kington

Je n'ai aucune expérience avec numpy en particulier, mais l'ajout d'un point décimal avec un nombre de zéros peut aider. Par exemple, l'utilisation de 1,0000 au lieu de 1. Dans la normale des scripts python où j'ai eu ce problème, ce qui a aidé, donc, sauf si votre problème est causé par une bizarrerie dans numpy et n'a rien à voir avec python, ce qui devrait aider.

Bonne chance!

OriginalL'auteur Sinthet