Grosse matrice symétrique - comment stocker et utiliser intelligemment - Python

Question

J'ai une matrice symétrique. Maintenant, le problème est que j'ai besoin de remplir une telle matrice de dimensions (32**3) x (32**3). La raison pour laquelle j'ai besoin de remplir la matrice est parce que dans mon programme je l'utilise ensuite pour divers calculs: je l'inverse, je le multiplie avec d'autres matrices ... et il me semble que pour effectuer ces diverses calculs que vous avez besoin de stocker réellement la matrice complète, vous ne pouvez pas utiliser par exemple seulement la moitié de cela (mais je peux me tromper, auquel cas s'il vous plaît dites-moi comment je devrais faire).

Le problème est qu'une telle matrice est tout simplement trop grand pour mon ordinateur et je reçois l'erreur suivante:

Traceback (most recent call last): 
    File "program.py", line 191, in <module> 
     A = zeros((n_x*n_y*n_z. n_x*n_y*n_z), float) 
MemoryError 

Ici, n_x = 32. Alors, comment puis-je résoudre ce problème? Y a-t-il un moyen de stocker une matrice aussi grande, ou une façon intelligente d'éviter de la stocker? Les deux façons seraient bien pour moi, à condition que je puisse les utiliser sans faire d'erreurs dans les calculs.

Pour être complet, je déclare dans ce qui suit comment la matrice A est construite:

n_x = n_y = n_z = 32 
    L_x = L_y = L_z = n_x 
    A = zeros((n_x*n_y*n_z , n_x*n_y*n_z), float) 
    P_0 = 50.0 
    sigma_x = sigma_y = sigma_z = 0.9 
    sigma_root = np.sqrt(sigma_x**2 + sigma_y**2 + sigma_z**2) 
    twosigmasquared = 2.*sigma_root**2 
    for l in range(n_x*n_y*n_z): 
     for m in range(n_x*n_y*n_z): 
      A[l][m] = P_0*(L_x/(np.sqrt(2.*np.pi)*sigma_root*n_x**2)) * (L_y/(np.sqrt(2.*np.pi)*sigma_root*n_y**2)) * (L_z/(np.sqrt(2.*np.pi)*sigma_root*n_z**2))*np.exp((-((x[l]-x[m])**2)-((y[l]-y[m])**2)-((z[l]-z[m])**2))/twosigmasquared) 
      A[m][l] = A[l][m] 

Answer 1

Pour économiser 50% d'espace, vous pouvez utiliser une matrice lil_sparse de scipy.

from scipy import sparse as S 
A = S.lil_matrix((n_x*n_y*n_z , n_x*n_y*n_z), float) 

n_x = n_y = n_z = 32 
L_x = L_y = L_z = n_x 
P_0 = 50.0 
sigma_x = sigma_y = sigma_z = 0.9 
sigma_root = np.sqrt(sigma_x**2 + sigma_y**2 + sigma_z**2) 
twosigmasquared = 2.*sigma_root**2 
for l in range(n_x*n_y*n_z): 
    for m in range(l, n_x*n_y*n_z): # Filling only the top half 
     A[l][m] = P_0*(L_x/(np.sqrt(2.*np.pi)*sigma_root*n_x**2)) * (L_y/(np.sqrt(2.*np.pi)*sigma_root*n_y**2)) * (L_z/(np.sqrt(2.*np.pi)*sigma_root*n_z**2))*np.exp((-((x[l]-x[m])**2)-((y[l]-y[m])**2)-((z[l]-z[m])**2))/twosigmasquared) 

Alors, au lieu d'accéder à la matrice elle-même, vous pouvez écrire une fonction d'assistance:

def getA(i, j): 
    if i < j: 
     return A[j, i] 
    else: 
     return A[i, j] 

---

Cependant, cela ne va pas vous aider beaucoup si vous voulez calculer l'inverse de la matrice en utilisant des approches standard, veulent multiplier la matrice efficacement ou faire des opérations du tout. Sauvegarder la matrice entière en mémoire est probablement un meilleur choix.