[Python] Article qui permet le calcul rapide de matrices éparses

introduction

Une matrice creuse est une matrice dans laquelle la plupart des éléments de la matrice sont 0. Il est couvert dans un large éventail de domaines tels que le calcul numérique, l'apprentissage automatique et la théorie des graphes (je pense). Étant donné que la majeure partie de la matrice est 0, il est possible d'effectuer un calcul matriciel à grande vitesse en ignorant 0. Par exemple, si vous calculez en utilisant uniquement la partie non nulle dans une matrice creuse à grande échelle (1 million de lignes x 1 million de colonnes) où près de 90% est égal à 0, le temps de calcul sera d'environ 1/10 et le calcul peut être effectué rapidement. Il est facile d'imaginer devenir. En d'autres termes, avec une certaine ingéniosité, il est possible de réduire la quantité de calcul par itération de $ O (N ^ 2) $ à $ O (N) $. Dans cet article, j'expliquerai ce calcul de matrice clairsemée à l'aide de Python, Numpy et Scipy. Si vous trouvez quelque chose que vous souhaitez ajouter, nous le mettrons à jour en conséquence. Si vous avez des erreurs, veuillez nous en informer dans les commentaires.

Public cible

• Les personnes qui peuvent gérer Python
• Ceux qui s'intéressent au calcul matriciel

Aperçu

Après avoir résumé les concepts nécessaires pour gérer les calculs matriciels épars en Python, nous nous référerons aux bibliothèques directes et itératives. Enfin, le problème créé dans Article précédent est résolu par une méthode appelée méthode d'itération non stationnaire.

Bibliothèque à utiliser: Scipy, Numpy, matplotlib

Résultat du calcul

Contenu de l'article

1. [Matrice clairsemée](matrice # 1-clairsemée)
2. [Matrice dense et clairsemée](# 1-1-Matrice dense et clairsemée)
3. [Visualisation de la matrice clairsemée](# 1-2-Visualisation de la matrice clairsemée)
4. [Format de stockage de matrice éparse](format de stockage de matrice # 2-Sparse)
5. Liste des listes (LIL) (# 2-1-List of Lists list-of-lists - lil)
6. [Méthode de stockage de lignes compressées (CSR)](# 2-2-Méthode de stockage de lignes compressées compressées-sparse-row-csr)
7. [Méthode de stockage sur colonne compressée (CSC)](# 2-3-Méthode de stockage sur colonne compressée compressée-sparse-column-csc)
8. [Autre](# 2-4-Autre)
9. [Quatre règles, etc.](# 3-Quatre règles, etc.)
10. [Format Matrix Market](# 4-matrix-market-format)
11. [Bibliothèque de méthodes directes](# 5-Bibliothèque de méthodes directes)
12. Bibliothèque itérative (# 6-Bibliothèque itérative)

1. Matrice clairsemée

1-1. Matrices denses et éparses

• Matrice dense

Procession dense. Une matrice avec quelques 0 dans la matrice. Il y a peut-être beaucoup de gens qui sont plus familiers avec cela, comme les problèmes de mathématiques.

  \left(
    \begin{array}{cccc}
      2  & 1 &  1 & \ldots  & \ldots & 1 \\
      1 & 2  & -1 & 1 & \ldots & 1 \\
      1  & 1 & 2  & 1 & 1 & \vdots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \vdots\\
      1  & 1 & \ddots & 1 & 2  & 1 \\
      1  & 1 & \ldots & 1 & 1 & 2 
    \end{array}
  \right)

• Matrice clairsemée

Une procession clairsemée. C'est juste qu'il y a beaucoup de 0 dans la matrice. Si vous concevez la partie du calcul liée à 0, vous pouvez difficilement la calculer, alors calculez-la par une méthode spéciale. Il sera traité dans cet article.

  \left(
    \begin{array}{cccc}
      2  & -1 &  0 & \ldots  & \ldots & 0 \\
      -1 & 2  & -1 & 0 & \ldots & 0 \\
      0  &-1 & 2  & -1 & \ddots & \vdots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & 0\\
      0  & 0 & \ddots & -1 & 2  & -1 \\
      0  & 0 & \ldots & 0 & -1 & 2 
    \end{array}
  \right)

1-2. Visualisation d'une matrice clairsemée

Je pense que la meilleure façon de visualiser des matrices éparses est d'utiliser la fonction d'espionnage de matplotlib. Il s'agit d'une matrice, où la partie 0 est affichée en blanc, et s'il y a un nombre, elle est affichée en couleur. En supposant que la matrice creuse ci-dessus est une matrice de 100 $ \ fois 100 $, cela ressemble à ceci:

2. Format de stockage de matrice éparse

Afin de gérer la matrice creuse ignorant 0, nous exprimons la matrice creuse en utilisant une méthode spéciale. En faisant cela, il est possible de réduire considérablement la mémoire utilisée et la quantité de calcul. Vous n'êtes peut-être pas sûr, alors je vais vous expliquer en utilisant la matrice ci-dessous. Cependant, a = 1, b = 2, c = 3, d = 4, e = 5, f = 6.

  \left(
    \begin{array}{cccc}
      a & 0 & 0 & 0 & 0 & 0 \\
      b & 0 & 0 & 0 & 0 & 0 \\
      0 & 0 & 0 & c & 0 & d \\
      0 & 0 & 0 & 0 & 0 & 0 \\
      0 & 0 & 0 & 0 & e & 0 \\
      0 & 0 & 0 & 0 & 0 & f 
    \end{array}
  \right)

2-1. Liste des listes (LIL)

Stocke un tuple de (colonnes, éléments) pour chaque ligne. Lors du premier stockage de données dans une matrice, il est plus facile d'utiliser la méthode LIL. Après cela, je pense qu'il devrait être facile de comprendre la conversion à la méthode CRS ou CCS ci-dessous. Une fois que vous vous êtes habitué à travailler avec des listes Python et des tableaux Numpy, ce n'est pas si difficile. Dans l'exemple ci-dessus,

Il est stocké sous la forme de, et vous pouvez voir qu'une matrice creuse peut être décrite avec une petite quantité de mémoire. Plus la matrice clairsemée est grande, meilleure est la méthode de stockage qui ignore ces zéros.

scipy.sparse.lil_matrix

Fonction lil_matrix pour gérer le format LIL fourni par scipy. Il est souvent utilisé lors de la génération de matrices creuses car il est facile de saisir des éléments non nuls de la matrice.

from scipy.sparse import lil_matrix
 
# Create sparse matrix
a=lil_matrix((6,6))
 
# Set the non-zero values
a[0,0]=1.; a[1,0]=2.; a[2, 3]=3.; a[2, 5]=4.; a[4, 4]=5.; a[5, 5]=6.

print("(row, column) value")
print(a, "\n")
print("Normal type: ")
print(a.todense())  #Retour à la forme matricielle normale
print("LIL type: ")
print(a.rows)  #Dans quelle colonne chaque ligne contient des éléments non nuls
print(a.data)  #Élément différent de zéro dans chaque ligne

Le résultat de l'exécution est le suivant.

(row, column) value
  (0, 0)	1.0
  (1, 0)	2.0
  (2, 3)	3.0
  (2, 5)	4.0
  (4, 4)	5.0
  (5, 5)	6.0
  
Normal type:
[[1. 0. 0. 0. 0. 0.]
 [2. 0. 0. 0. 0. 0.]
 [0. 0. 0. 3. 0. 4.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 5. 0.]
 [0. 0. 0. 0. 0. 6.]]

LIL type: 
[list([0]) list([0]) list([3, 5]) list([]) list([4]) list([5])]
[list([1.0]) list([2.0]) list([3.0, 4.0]) list([]) list([5.0]) list([6.0])]

2-2. Compressed Sparse Row (CSR)

Certains l'appellent CRS (Compressed Row Storage). Recherchez dans le sens de la ligne et stockez les éléments de la matrice non nulle. S'il est difficile à comprendre, veuillez lire le Matériel du professeur Nakajima de l'Université de Tokyo. Un exemple d'exécution de la méthode CSR en excluant des composants diagonaux est présenté.

Si vous utilisez la méthode CSR,

  \left(
    \begin{array}{cccc}
      a & 0 & 0 & 0 & 0 & 0 \\
      b & 0 & 0 & 0 & 0 & 0 \\
      0 & 0 & 0 & c & 0 & d \\
      0 & 0 & 0 & 0 & 0 & 0 \\
      0 & 0 & 0 & 0 & e & 0 \\
      0 & 0 & 0 & 0 & 0 & f 
    \end{array}
  \right)

Est

• indices: IA = [0, 0, 3, 5, 4, 5]
• Affiche le numéro de colonne de l'élément dans l'ordre à partir de la ligne supérieure. La taille des indices est le nombre d'éléments non nuls.
• data: A = [a, b, c, d, e, f]
• Indique les éléments non nuls selon les indices. les données auront la même taille que les indies.
• indptr: IN = [0, 1, 2, 4, 4, 5, 6]
• Comptez le nombre d'éléments dans chaque ligne et affichez la valeur totale. Puisqu'il commence à 0, la taille est le nombre de lignes + 1.

Il peut être exprimé par trois listes de. Encore une fois, vous pouvez voir que vous pouvez écrire une matrice creuse avec moins de mémoire. Je pense que ce format CSR est souvent utilisé lors de la réalisation de calculs matriciels clairsemés. La rapidité avec laquelle le calcul est effectué est décrite au chapitre 3, Quatre règles.

Pour montrer ce que ces listes représentent, indices est l'indice des éléments de la matrice d'origine, et indptr est à droite.

  \left(
    \begin{array}{cccccc|c}
      col0    & col1   & col2  &  col3   & col4  & col5 & indptr  \\
          &   &   &     &   & & 0 \\
      a_0 & 0 & 0 & 0   & 0 & 0 & 1 \\
      b_0 & 0 & 0 & 0   & 0 & 0 & 2 \\
      0   & 0 & 0 & c_3 & 0 & d_5& 4 \\
      0   & 0 & 0 & 0   & 0 & 0 & 4 \\
      0   & 0 & 0 & 0   & e_4 & 0 & 5 \\
      0   & 0 & 0 & 0   & 0 & f_5 & 6 
    \end{array}
  \right)

Pour expliquer l'idée à l'aide de cet exemple,

• Etat initial
  • indices: IA = []
  • data: A = []
• ʻIndptr: IN = [0, 0, 0, 0, 0, 0, 0] `(taille du nombre de lignes + 1).
• Ligne 0
• Recherche les éléments non nuls dans l'ordre à partir de la 0ème colonne.
• Découvrez a! Ajoutez l'élément a à A.
• Ajoutez 0 à IA car a est dans la colonne 0.
• Sur la ligne 0, nous avons trouvé un élément, alors ajoutez 1 à ʻIN [1] `.
  • IA = [0], A = [a], IN = [0, 1, 0, 0, 0, 0, 0]
• La première ligne
• Recherche les éléments non nuls dans l'ordre à partir de la 0ème colonne.
• Découvrez b! Ajoutez l'élément b à A.
• b est dans la 0ème colonne, donc ajoutez 0 à IA.
• Dans la première ligne, nous avons trouvé un élément, alors ajoutez 1 + IN [1] à ʻIN [2] `.
  • IA = [0, 0], A = [a, b], IN = [0, 1, 2, 0, 0, 0, 0]
• Répétez ceci à partir de la deuxième ligne.
  • IA = [0, 0, 3, 5, 4, 5], A = [a, b, c, d, e, f], IN = [0, 1, 2, 4, 4, 5, 6]

scipy.sparse.csr_matrix

Fonction csr_matrix pour gérer le format CSR fourni par scipy. Si vous mettez une matrice arbitraire mtx comme scipy.sparse.csr_matrix (mtx), elle sera convertie au format CSR. Le format de mtx est list, numpy.ndarray, numpy.matrix, etc. Lors du passage du format LIL au format CSR, vous pouvez utiliser tocsr ().

a_csr = a.tocsr()
print("Normal type: ")
print(a_csr.todense())
print("csr type:")
print(a_csr.indices)
print(a_csr.data)
print(a_csr.indptr)

Normal type:
[[1. 0. 0. 0. 0. 0.]
 [2. 0. 0. 0. 0. 0.]
 [0. 0. 0. 3. 0. 4.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 5. 0.]
 [0. 0. 0. 0. 0. 6.]]

csr type:
[0 0 3 5 4 5]  # indices:Le numéro de colonne de l'élément est affiché dans l'ordre à partir de la ligne supérieure.
[1. 2. 3. 4. 5. 6.]  # data:Les éléments non nuls sont affichés selon des indices.
[0 1 2 4 4 5 6]  # indptr:Comptez le nombre d'éléments dans chaque ligne et additionnez-les. Nombre de lignes car il commence à 0+Il devient la taille de 1.

2-3. Colonne clairsemée compressée (CSC)

Certains l'appellent CCS. Recherchez dans la direction des colonnes et stockez les éléments de la matrice non nulle. Version en colonne de la méthode CSR.

scipy.sparse.csc_matrix

Fonction csc_matrix pour gérer le format CSC fourni par scipy. Si vous mettez une matrice arbitraire mtx comme scipy.sparse.csc_matrix (mtx), elle sera convertie au format CSC. Le format de mtx est list, numpy.ndarray, numpy.matrix, etc. Lors du passage du format LIL au format CSC, vous pouvez utiliser tocsc ().

a_csc = a.tocsc()
print("Normal type: ")
print(a_csc.todense())
print("csc type:")
print(a_csc.indices)
print(a_csc.data)
print(a_csc.indptr)

Normal type:
[[1. 0. 0. 0. 0. 0.]
 [2. 0. 0. 0. 0. 0.]
 [0. 0. 0. 3. 0. 4.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 5. 0.]
 [0. 0. 0. 0. 0. 6.]]

csc type:
[0 1 2 4 2 5]  # indices:La ligne de l'élément est affichée dans l'ordre à partir de la colonne de gauche.
[1. 2. 3. 5. 4. 6.]  # data:Les éléments non nuls sont affichés selon des indices.
[0 2 2 2 3 4 6]  # indptr:Comptez le nombre d'éléments dans chaque colonne et additionnez-les. Puisqu'il commence à 0, le nombre de colonnes+Il devient la taille de 1.

2-4. Autre

Il existe diverses autres méthodes de stockage matriciel, donc si vous voulez en savoir plus, veuillez consulter [Wikipedia](https://ja.wikipedia.org/wiki/sparse matrix).

3. Quatre règles, etc.

Lors du calcul d'une matrice creuse, il est plus courant de la stocker dans la méthode CSR ou la méthode CSC avant d'effectuer le calcul. Elle peut être effectuée de la même manière que l'opération classique à quatre règles numpy.

# csr example
# Create sparse matrix
a = scipy.sparse.lil_matrix((6,6))
# Set the non-zero values
a[0,0]=1.; a[1,0]=2.; a[2, 3]=3.; a[2, 5]=4.; a[4, 4]=5.; a[5, 5]=6.
a_csr = a.tocsr()
a_matrix = a.todense()

b = np.arange(6).reshape((6,1))
print(b)
print(a_matrix * b)
print(a_csr * b)
print(a_matrix.dot(b))
print(a_csr.dot(b))

La multiplication entre des matrices éparses stockées au format CSR peut également être calculée comme ceci.

b_sparse = scipy.sparse.lil_matrix((6,6))
b_sparse[0,1]=1.; b_sparse[3,0]=2.; b_sparse[3, 3]=3.; b_sparse[5, 4]=4.
b_csr = b_sparse.tocsr()
print(b_csr.todense())
print(a_matrix * b_csr)
print((a_csr * b_csr).todense())

Sur cette base, vérifions l'efficacité de la méthode CSR (stockage de lignes compressées). Essayez d'exécuter le code suivant sur votre bloc-notes Jupyter. Vous pouvez voir que la méthode ** CSR complète le calcul considérablement (environ 1000 fois) plus rapidement que la résolution à l'aide d'une matrice ordinaire. ** **

import scipy.sparse
import numpy as np
import matplotlib.pyplot as plt
num = 10**4
# Creat LIL type sparse matrix
a_large = scipy.sparse.lil_matrix((num,num))
# Set the non-zero values
a_large.setdiag(np.ones(num)*2)
a_large.setdiag(np.ones(num-1)*-1, k=1)
a_large.setdiag(np.ones(num-1)*-1, k=-1)
a_large_csr = a_large.tocsr()
a_large_dense = a_large.todense()
plt.spy(a_large_csr)
b = np.ones(num)
%timeit a_large_csr.dot(b)  #Méthode CSR
%timeit a_large_dense.dot(b)  #Matrice ordinaire

point important

Si vous mettez une matrice de format csr dans l'argument de la fonction point de numpy comme np.dot (a_csr, b),

[[<6x6 sparse matrix of type '<class 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Row format>]
 [<6x6 sparse matrix of type '<class 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Row format>]
 [<6x6 sparse matrix of type '<class 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Row format>]
 [<6x6 sparse matrix of type '<class 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Row format>]
 [<6x6 sparse matrix of type '<class 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Row format>]
 [<6x6 sparse matrix of type '<class 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Row format>]]

Je ne peux pas bien calculer.

4. Format Matrix Market

Bien que non utilisé dans cet article, je mentionnerai également Matrix Market Format, qui est l'une des représentations des matrices clairsemées. Les données matricielles clairsemées sont essentiellement distribuées dans ce format ou format Rutherford Boeing, et il est nécessaire de se référer aux documents et implémentations liés au calcul matriciel clairsemé. La SuiteSparse Matrix Collection anciennement la Sparse Matrix Collection de l'Université de Floride est probablement l'ensemble de données le plus connu pour le problème de la matrice clairsemée.

  \left(
    \begin{array}{cccc}
      a & 0 & 0 & 0 & 0 & 0 \\
      b & 0 & 0 & 0 & 0 & 0 \\
      0 & 0 & 0 & c & 0 & d \\
      0 & 0 & 0 & 0 & 0 & 0 \\
      0 & 0 & 0 & 0 & e & 0 \\
      0 & 0 & 0 & 0 & 0 & f 
    \end{array}
  \right)

Est exprimé au format Matrix Market (MM)

%%MatrixMarket matrix coordinate real general
%=================================================================================
%
% This ASCII file represents a sparse MxN matrix with L 
% nonzeros in the following Matrix Market format:
%
% +----------------------------------------------+
% |%%MatrixMarket matrix coordinate real general | <--- header line
% |%                                             | <--+
% |% comments                                    |    |-- 0 or more comment lines
% |%                                             | <--+         
% |    M  N  L                                   | <--- rows, columns, entries
% |    I1  J1  A(I1, J1)                         | <--+
% |    I2  J2  A(I2, J2)                         |    |
% |    I3  J3  A(I3, J3)                         |    |-- L lines
% |        . . .                                 |    |
% |    IL JL  A(IL, JL)                          | <--+
% +----------------------------------------------+   
%
% Indices are 1-based, i.e. A(1,1) is the first element.
%
%=================================================================================
  6  6  6
    1     1   a
    2     1   b
    3     4   c
    5     5   e
    3     6   d
    6     6   f

Ce sera. Vous pouvez le voir dans les commentaires ci-dessus. Cet article est facile à comprendre. Pour le moment, je jouerai avec une procession clairsemée appelée steam3 que j'ai choisi de manière appropriée. Utilisez scipy.io lors de la lecture du format de marché Matrix.

• Lancer scipy.io.mminfo vous indiquera le nombre de lignes et de colonnes de la matrice et le nombre d'éléments non nuls.
• Vous pouvez lire une matrice au format Matrix Market avec scipy.io.mmread. La valeur renvoyée est numpy.array ou format coo_matrix (similaire au format de stockage csr).

import scipy.io as io
print(io.mminfo("steam3/steam3.mtx"))
steam3 = io.mmread("steam3/steam3.mtx").tocsr()
print(steam3.todense())
plt.spy(steam3)

(80, 80, 928, 'coordinate', 'real', 'general')
matrix([[-3.8253876e+05,  0.0000000e+00,  0.0000000e+00, ...,
          0.0000000e+00,  0.0000000e+00,  0.0000000e+00],
        [ 9.1147205e-01, -9.5328382e+03,  0.0000000e+00, ...,
          0.0000000e+00,  0.0000000e+00,  0.0000000e+00],
        [ 0.0000000e+00,  0.0000000e+00, -2.7754141e+02, ...,
          0.0000000e+00,  0.0000000e+00,  0.0000000e+00],
        ...,
        [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
         -4.6917325e+08,  0.0000000e+00, -1.2118734e+05],
        [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
          0.0000000e+00, -1.3659626e+07,  0.0000000e+00],
        [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
          2.1014810e+07, -3.6673643e+07, -1.0110894e+04]])

5. Bibliothèque de méthodes directes

Il se trouve dans scipy.sparse.linalg.dsolve. J'ai implémenté le problème d'Ax = b (voir la formule ci-dessous pour plus de détails) résolu dans Article précédent à titre d'exemple (l'implémentation n'a pas beaucoup changé). ).

  \left(
    \begin{array}{cccc}
      2 \left(\frac{1}{d} + 1 \right) & -1 &  0 & \ldots  & \ldots & 0 \\
      -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots & 0 \\
      0  &-1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \ddots\\
      0  & 0 & \ddots & -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 \\
      0  & 0 & \ldots & 0 & -1 & 2 \left(\frac{1}{d} + 1 \right)
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      T_1^{n+1}  \\
      T_2^{n+1}  \\
      T_3^{n+1}    \\
      \vdots \\
      T_{M-1}^{n+1} \\
      T_M^{n+1} 
    \end{array}
  \right)
   = 
   \left(
    \begin{array}{c}
      T_2^{n} + 2 \left(\frac{1}{d} - 1 \right) T_1^{n} + \left(T_0^n + T_0^{n+1}\right) \\
      T_3^{n} + 2 \left(\frac{1}{d} - 1 \right) T_2^{n} + T_1^n  \\
      T_4^{n} + 2 \left(\frac{1}{d} - 1 \right) T_3^{n} + T_2^n  \\
      \vdots \\
      T_M^{n} + 2 \left(\frac{1}{d} - 1 \right) T_{M-1}^{n} + T_{M-2}^n  \\
      \left(T_{M+1}^n + T_{M+1}^{n+1}\right) + 2 \left(\frac{1}{d} - 1 \right) T_{M}^{n} + T_{M-1}^n
    \end{array}
  \right)

Le paramètre par défaut est scipy.sparse.linalg.dsolve (format csr A matrice, vecteur b) avec décomposition Super LU (http://www.turbare.net/transl/scipy-lecture-notes/advanced/ scipy_sparse / solvers.html) doit être fait.

L'image ci-dessous montre le résultat de la comparaison avec la méthode Gauss-Seidel de la méthode de répétition stationnaire.

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse
# Make stencils
# Creat square wave
Num_stencil_x = 101
x_array = np.float64(np.arange(Num_stencil_x))
temperature_array = x_array + 100
temperature_lower_boundary = 150
temperature_upper_boundary = 150
Time_step = 100
Delta_x = max(x_array) / (Num_stencil_x-1)
C = 1
Delta_t = 0.2
kappa = 0.5
d = kappa * Delta_t / Delta_x**2
total_movement = C * Delta_t * (Time_step+1)
exact_temperature_array = (temperature_upper_boundary - temperature_lower_boundary) / (x_array[-1] - x_array[0]) * x_array + temperature_lower_boundary
plt.plot(x_array, temperature_array, label="Initial condition")
print("Δx:", Delta_x, "Δt:", Delta_t, "d:", d)

temperature_sp = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_sp)) * 2 *(1/d+1) \
                - np.eye(len(temperature_sp), k=1) \
                - np.eye(len(temperature_sp), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_sp), temperature_upper_boundary)
    b_array = 2 * (1/d - 1) * temperature_sp + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    a_csr = scipy.sparse.csr_matrix(a_matrix)
    temperature_sp = spla.dsolve.spsolve(a_csr, b_array)
plt.plot(x_array, temperature_sp, label="SuperLU")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

6. Bibliothèque itérative

Une bibliothèque de méthodes d'itération non stationnaires est implémentée dans scipy. Il se trouve dans scipy.sparse.linalg.isolve. Si vous n'êtes pas familier avec la méthode d'itération non stationnaire, consultez Article précédent. Certaines des bibliothèques qui peuvent être utilisées sont répertoriées ci-dessous. Si vous souhaitez utiliser une bibliothèque autre que celles-ci, recherchez-la dans le manuel de Scipy (https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#module-scipy.sparse.linalg). S'il te plait donne moi.

1. Matrice symétrique
2. Méthode du gradient conjugué cg (méthode CG) - Matrice symétrique à valeur constante positive uniquement
3. méthode résiduelle minimale «minres» (MINimum RESidual: méthode MINRES)
4. Matrice asymétrique
5. Gradient bi-conjugué bicg (méthode BiCG)
6. Méthode du gradient conjugué au carré cgs (Conjugate Residual Squared: méthode CGS)
7. Stabilisation BiCG bicgstab (méthode BiCG STAB)
8. Résiduel minimal généralisé gmres (méthode GMRES)
9. Méthode résidente quasi-minimale qmr (méthode QMR)

Voici un exemple d'implémentation de la méthode BiCGSTAB et de son résultat d'exécution pour le même problème utilisé dans la méthode directe.

temperature_sp = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_sp)) * 2 *(1/d+1) \
                - np.eye(len(temperature_sp), k=1) \
                - np.eye(len(temperature_sp), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_sp), temperature_upper_boundary)
    b_array = 2.0 * (1/d - 1.0) * temperature_sp + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    a_csr = scipy.sparse.csr_matrix(a_matrix)
    temperature_sp = spla.isolve.bicgstab(a_csr, b_array)[0]
plt.plot(x_array, temperature_sp, label="Scipy.isolve (1000 steps)")

plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

Résumé

• La matrice clairsemée a une méthode de stockage qui ignore 0
• La méthode LIL est souvent utilisée lors de la création d'une matrice creuse.
• La méthode CSR et la méthode CSC sont souvent utilisées lors de l'exécution de calculs matriciels épars.
• Il a été confirmé que le calcul par la méthode CSR est extrêmement rapide que le calcul sous la forme d'une matrice ordinaire.
• Les matrices éparses sont souvent distribuées en utilisant la méthode Matrix Market.
• Les fonctions itératives directes et non stationnaires sont implémentées dans Scipy et sont faciles à utiliser.

Je pense que c'est comme ça.

Les références

1. Matrice clairsemée dans Scipy
2. https://www.atmarkit.co.jp/ait/articles/1706/21/news018_2.html
3. http://xxxxxeeeee.hatenablog.com/entry/20100825/1282743025
4. https://qiita.com/Hase8388/items/369fbd57318d8a976999
5. Matrix Market File Format