[Python] Artikel, der eine schnelle Berechnung der Sparse-Matrix ermöglicht

Einführung

Eine dünn besetzte Matrix ist eine Matrix, in der die meisten Elemente der Matrix 0 sind. Es wird in einer Vielzahl von Bereichen wie numerische Berechnung, maschinelles Lernen und Graphentheorie behandelt (glaube ich). Da der größte Teil der Matrix 0 ist, ist es möglich, eine Matrixberechnung mit hoher Geschwindigkeit durchzuführen, indem 0 ignoriert wird. Wenn Sie beispielsweise nur den Nicht-Null-Teil in einer großen Matrix mit geringer Dichte (1 Million Zeilen x 1 Million Spalten) berechnen, in der fast 90% 0 sind, beträgt die Berechnungszeit etwa 1/10 und die Berechnung kann schnell abgeschlossen werden. Es ist leicht vorstellbar zu werden. Mit anderen Worten, mit einigem Einfallsreichtum ist es möglich, den Rechenaufwand pro Iteration von $ O (N ^ 2) $ auf $ O (N) $ zu reduzieren. In diesem Artikel werde ich eine solche spärliche Matrixberechnung mit Python, Numpy und Scipy erläutern. Wenn Sie etwas finden, das Sie hinzufügen möchten, werden wir es entsprechend aktualisieren. Wenn Sie Fehler haben, teilen Sie uns dies bitte in den Kommentaren mit.

Zielgruppe

Überblick

Nachdem wir die Konzepte zusammengefasst haben, die für die Berechnung spärlicher Matrixberechnungen in Python erforderlich sind, werden wir auf die direkten und iterativen Bibliotheken verweisen. Schließlich wird das in Vorheriger Artikel erstellte Problem durch eine Methode gelöst, die als instationäre Iterationsmethode bezeichnet wird.

Zu verwendende Bibliothek: Scipy, Numpy, matplotlib

Berechnungsergebnis

anim.gif

Inhalt im Artikel

  1. [Sparse Matrix](# 1-Sparse Matrix)
  2. [Dichte und spärliche Matrix](# 1-1-Dichte und spärliche Matrix)
  3. [Visualisierung der Sparse-Matrix](# 1-2-Visualisierung der Sparse-Matrix)
  4. [Sparse-Matrix-Speicherformat](# 2-Sparse-Matrix-Speicherformat)
  5. Liste der Listen (LIL) (# 2-1-Liste der Listen Liste der Listen - lil)
  6. [Komprimierte Zeilenspeichermethode (CSR)](# 2-2-Komprimierte Zeilenspeichermethode Compressed-Sparse-Row-CSR)
  7. [Komprimierte Spaltenspeichermethode (CSC)](# 2-3-Komprimierte Spaltenspeichermethode Compressed-Sparse-Column-CSC)
  8. [Andere](# 2-4-Andere)
  9. [Vier Regeln usw.](Nr. 3 - Vier Regeln usw.)
  10. [Matrix-Marktformat](# 4-Matrix-Marktformat)
  11. [Direkte Methodenbibliothek](# 5-Direkte Methodenbibliothek)
  12. Iterative Bibliothek (# 6-Iterative Bibliothek)

1. Sparse Matrix

1-1. Dichte und spärliche Matrizen

Dichte Prozession. Eine Matrix mit wenigen Nullen in der Matrix. Vielleicht gibt es viele Leute, die damit besser vertraut sind, wie zum Beispiel mathematische Probleme.

  \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)

Eine spärliche Prozession. Es ist nur so, dass die Matrix viele Nullen enthält. Wenn Sie den Teil der Berechnung in Bezug auf 0 erstellen, können Sie ihn kaum berechnen. Berechnen Sie ihn daher mit einer speziellen Methode. Es wird in diesem Artikel behandelt.

  \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. Visualisierung der spärlichen Matrix

Ich denke, der beste Weg, um spärliche Matrizen zu visualisieren, ist die Verwendung der Spionagefunktion von matplotlib. Dies ist eine Matrix, in der der 0-Teil in Weiß angezeigt wird. Wenn eine Zahl vorhanden ist, wird er in Farbe angezeigt. Angenommen, die obige Sparse-Matrix ist eine $ 100 \ times 100 $ -Matrix, sieht sie folgendermaßen aus:

100100example.png

2. Sparse-Matrix-Speicherformat

Um mit einer Sparse-Matrix umzugehen, die 0 ignoriert, drücken wir eine Sparse-Matrix mit einer speziellen Methode aus. Auf diese Weise ist es möglich, den verwendeten Speicher und den Rechenaufwand erheblich zu reduzieren. Sie sind sich möglicherweise nicht sicher, daher erkläre ich dies anhand der unten gezeigten Matrix. 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)

example.png

2-1. Liste der Listen (LIL)

Speichert (Spalte, Element) Taples für jede Zeile. Wenn Sie Daten zum ersten Mal in einer Matrix speichern, ist die Verwendung der LIL-Methode einfacher. Danach sollte es meiner Meinung nach leicht verständlich sein, auf die unten gezeigte CRS- oder CCS-Methode zu konvertieren. Sobald Sie sich an die Arbeit mit Python-Listen und Numpy-Arrays gewöhnt haben, ist dies nicht mehr so schwierig. Im obigen Beispiel ist

Linie (Säule,Element)
0 (0, a)
1 (0, b)
2 (3, c), (5, d)
4 (4, e)
5 (5, f)

Es wird in Form von gespeichert, und Sie können sehen, dass eine spärliche Matrix mit wenig Speicher beschrieben werden kann. Je größer die dünne Matrix ist, desto besser ist die Speichermethode, bei der diese Nullen ignoriert werden.

scipy.sparse.lil_matrix

Funktion lil_matrix zur Behandlung des von scipy bereitgestellten LIL-Formats. Es wird häufig beim Generieren von Matrizen mit geringer Dichte verwendet, da es einfach ist, Nicht-Null-Elemente der Matrix einzugeben.

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())  #Kehren Sie zur normalen Matrixform zurück
print("LIL type: ")
print(a.rows)  #In welcher Spalte enthält jede Zeile Elemente ungleich Null
print(a.data)  #Nicht-Null-Element in jeder Zeile

Das Ausführungsergebnis ist wie folgt.

(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)

Einige nennen es CRS (Compressed Row Storage). Suchen Sie in Zeilenrichtung und speichern Sie die Elemente der Nicht-Null-Matrix. Wenn es schwer zu verstehen ist, lesen Sie bitte das Material von Professor Nakajima von der Universität Tokio. Ein Beispiel für die Ausführung der CSR-Methode durch Ausschließen diagonaler Komponenten wird gezeigt.

Wenn Sie die CSR-Methode verwenden,

  \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)

Ist

Es kann durch drei Listen von ausgedrückt werden. Wieder können Sie sehen, dass Sie eine dünne Matrix mit weniger Speicher schreiben können. Ich denke, dass dieses CSR-Format häufig bei der Durchführung von Sparse-Matrix-Berechnungen verwendet wird. Wie schnell die Berechnung abgeschlossen ist, wird in Kapitel 3, Vier Regeln beschrieben.

Um zu zeigen, was diese Listen darstellen, ist Indizes der Index für die Elemente der ursprünglichen Matrix, und indptr befindet sich rechts.

  \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)

Um die Idee anhand dieses Beispiels zu erklären,

scipy.sparse.csr_matrix

Funktion csr_matrix zur Behandlung des von scipy bereitgestellten CSR-Formats. Wenn Sie eine beliebige Matrix mtx wie "scipy.sparse.csr_matrix (mtx)" einfügen, wird diese in das CSR-Format konvertiert. Das Format von mtx ist list, numpy.ndarray, numpy.matrix usw. Wenn Sie vom LIL-Format zum CSR-Format wechseln, können Sie tocsr () verwenden.

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:Die Spaltennummer des Elements wird in der Reihenfolge ab der oberen Zeile angezeigt.
[1. 2. 3. 4. 5. 6.]  # data:Die Nicht-Null-Elemente werden nach Indizes angezeigt.
[0 1 2 4 4 5 6]  # indptr:Zählen Sie, wie viele Elemente sich in jeder Zeile befinden, und addieren Sie sie. Anzahl der Zeilen, da sie bei 0 beginnen+Es wird die Größe von 1.

2-3. Compressed Sparse Column (CSC)

Manche nennen es CCS. Suchen Sie in Spaltenrichtung und speichern Sie die Elemente der Nicht-Null-Matrix. Spaltenversion der CSR-Methode.

scipy.sparse.csc_matrix

Funktion csc_matrix zur Verarbeitung des von scipy bereitgestellten CSC-Formats. Wenn Sie eine beliebige Matrix mtx wie "scipy.sparse.csc_matrix (mtx)" einfügen, wird diese in das CSC-Format konvertiert. Das Format von mtx ist list, numpy.ndarray, numpy.matrix usw. Wenn Sie vom LIL-Format zum CSC-Format wechseln, können Sie tocsc () verwenden.

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:Die Zeile des Elements wird in der Reihenfolge der linken Spalte angezeigt.
[1. 2. 3. 5. 4. 6.]  # data:Die Nicht-Null-Elemente werden nach Indizes angezeigt.
[0 2 2 2 3 4 6]  # indptr:Zählen Sie, wie viele Elemente sich in jeder Spalte befinden, und addieren Sie sie. Da es bei 0 beginnt, die Anzahl der Spalten+Es wird die Größe von 1.

2-4 Andere

Es gibt verschiedene andere Matrixspeichermethoden. Wenn Sie mehr wissen möchten, lesen Sie bitte [Wikipedia](https://ja.wikipedia.org/wiki/sparse matrix).

3. Vier Regeln usw.

Bei der Berechnung einer dünn besetzten Matrix wird diese am häufigsten in der CSR-Methode oder der CSC-Methode gespeichert, bevor die Berechnung durchgeführt wird. Es kann auf die gleiche Weise ausgeführt werden wie die bekannte Operation mit vier Regeln.

# 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))

Auf diese Weise kann auch die Multiplikation zwischen dünn besetzten Matrizen berechnet werden, die im CSR-Format gespeichert sind.

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())

Lassen Sie uns auf dieser Grundlage die Wirksamkeit der CSR-Methode (Compressed Row Storage) überprüfen. Versuchen Sie, den folgenden Code auf Ihrem Jupyter-Notebook auszuführen. Sie sehen, dass die ** CSR-Methode die Berechnung erheblich (etwa 1000-mal) schneller abschließt als das Lösen mit einer normalen Matrix. ** ** **

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)  #CSR-Methode
%timeit a_large_dense.dot(b)  #Gewöhnliche Matrix

wichtiger Punkt

Wenn Sie eine csr-Formatmatrix in das Argument der Punktfunktion von numpy wie np.dot (a_csr, b) einfügen,

[[<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>]]

Ich kann nicht gut rechnen.

4. Matrix Market Format

Obwohl in diesem Artikel nicht verwendet, werde ich auch das Matrix Market Format erwähnen, das eine der Darstellungen von spärlichen Matrizen ist. Sparse-Matrix-Daten werden grundsätzlich in diesem Format oder im Rutherford-Boeing-Format verteilt, und es ist erforderlich, dass auf Kenntnisse und Implementierungen im Zusammenhang mit der Sparse-Matrix-Berechnung verwiesen wird. Die SuiteSparse Matrix Collection, früher die Sparse Matrix Collection der Universität von Florida (https://sparse.tamu.edu), ist wahrscheinlich der bekannteste Datensatz für das Problem der spärlichen Matrix.

  \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)

Wird im MM-Format (Matrix Market) ausgedrückt

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

Es wird sein. Sie können es in den Kommentaren oben sehen. Dieser Artikel ist leicht zu verstehen. Vorerst werde ich mit einer spärlichen Prozession namens steam3 spielen, die ich entsprechend ausgewählt habe. Verwenden Sie "scipy.io", wenn Sie das Matrix-Marktformat lesen.

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]])

steam3.png

5. Direkte Methodenbibliothek

Es ist in scipy.sparse.linalg.dsolve. Ich habe das Problem von Ax = b (Einzelheiten siehe Formel unten) implementiert, das in Vorheriger Artikel als Beispiel gelöst wurde (die Implementierung hat sich nicht wesentlich geändert). ).

  \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)

Die Standardeinstellung ist "scipy.sparse.linalg.dsolve (CSR-Format A Matrix, B-Vektor)" mit Super-LU-Zerlegung (http://www.turbare.net/transl/scipy-lecture-notes/advanced/). scipy_sparse / solvers.html) sollte durchgeführt werden.

Das Bild unten zeigt das Vergleichsergebnis mit der Gauß-Seidel-Methode der stationären Wiederholungsmethode.

sp.png

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. Iterative Bibliothek

Eine Bibliothek instationärer Iterationsmethoden ist in scipy implementiert. Es ist in scipy.sparse.linalg.isolve. Wenn Sie mit der instationären Iterationsmethode nicht vertraut sind, lesen Sie Vorheriger Artikel. Einige der Bibliotheken, die verwendet werden können, sind unten aufgeführt. Wenn Sie eine andere Bibliothek als diese verwenden möchten, suchen Sie im Handbuch von Scipy (https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#module-scipy.sparse.linalg) danach. Bitte gib mir.

  1. Symmetrische Matrix
  2. cg Conjugate Gradient-Methode (CG-Methode) - Nur positive symmetrische Matrix mit konstantem Wert
  3. Minimale Restmethode "minres" (MINimum RESidual: MINRES-Methode)
  4. Asymmetrische Matrix
  5. bicg bi-konjugierter Gradient (BiCG-Methode)
  6. "cgs" -Quadrat-Konjugat-Gradienten-Methode (Conjugate Residual Squared: CGS-Methode)
  7. bicgstab BiCG-Stabilisierung (BiCG STAB-Methode)
  8. gmres Generalized Minimal Residual (GMRES-Methode)
  9. "qmr" quasi-minimal residente Methode (QMR-Methode)

Das Folgende ist ein Implementierungsbeispiel der BiCG STAB-Methode und ihres Ausführungsergebnisses für dasselbe Problem, das in der direkten Methode verwendet wird.

anim.gif

sp_isolve.png

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)

Zusammenfassung

Ich denke es ist so.

Verweise

  1. Sparse Matrix in 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

Recommended Posts

[Python] Artikel, der eine schnelle Berechnung der Sparse-Matrix ermöglicht
Artikel, die die Systementwicklung mit Django (Python) ermöglichen _Einführung
Führen Sie mit Python eine Konvertierung in halber und voller Breite mit hoher Geschwindigkeit durch
[Python] Finden Sie die Anzahl der Fibonacci mit hoher Geschwindigkeit (Memoisierung, dynamische Planungsmethode)
[Python] Wie man den Bruchteil einer natürlichen Zahl mit hoher Geschwindigkeit erhält
Skript zum Herunterladen von AWS RDS-Protokolldateien mit hoher Geschwindigkeit
So beschleunigen Sie Python-Berechnungen
Eine kleine Geschichte, die Tabellendaten im CSV-Format mit hoher Geschwindigkeit ausgibt