A huge symmetric matrix - how it is convenient to store and use - Python

Question

I have a symmetric matrix. Now the problem is that I need to fill in such a matrix of dimensions (32**3) x (32**3) . The reason I need to fill in the matrix is ​​because in my program I use it for various calculations: I invert it, I multiply it by other matrices ... and it seems to me that you have to actually store these different ones the full matrix, you can’t use, for example, only half of this (but I can be wrong, in this case, please tell me how I should do this).

The problem is that such a matrix is ​​simply too large for my computer, and I get the following error:

 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 

Here n_x = 32 . So how can I solve this problem? Is there a way to store such a large matrix or a smart way to avoid storing it? Both methods would be good for me if I could use them without making mistakes in the calculations.

For completeness, I report the following how the matrix A constructed:

  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

To save 50% of space, you can use the lil_sparse matrix from 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) 

Then, instead of accessing the matrix itself, you can write an auxiliary function:

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

---

However, this will not help you if you want to calculate the inverse matrix using standard approaches, want to efficiently multiply the matrix, or perform any operations at all.

Storing the entire matrix in memory is probably the best choice.