[Python] Articles that enable sparse matrix calculations at high speed

Introduction

A sparse matrix is a matrix in which most of the elements of the matrix are 0. It is covered in a wide range of fields such as numerical calculation, machine learning, and graph theory (I think). Since most of the matrix is 0, it is possible to perform matrix calculation at high speed by ignoring 0. For example, if you calculate using only the non-zero part in a large-scale sparse matrix (1 million rows x 1 million columns) where nearly 90% is 0, the calculation time will be about 1/10 and the calculation can be completed quickly. You can easily imagine becoming. In other words, with some ingenuity, it is possible to reduce the amount of calculation per iteration from $ O (N ^ 2) $ to $ O (N) $. In this article, I will explain such sparse matrix calculation using Python, Numpy, and Scipy. If you find something you want to add, we will update it accordingly. If you have any mistakes, please let us know in the comments.

Target audience

Overview

After summarizing the concepts needed to handle sparse matrix calculations in Python, we will refer to the direct and iterative library. Finally, the problem created in Previous article is solved by a method called non-stationary iterative method.

Library to use: Scipy, Numpy, matplotlib

Calculation result

anim.gif

Table of contents in the article

  1. [Sparse matrix](# 1-Sparse matrix)
  2. [Dense and sparse matrix](# 1-1-Dense and sparse matrix)
  3. [Visualization of sparse matrix](# 1-2-Visualization of sparse matrix)
  4. [Sparse matrix storage format](# 2-Sparse matrix storage format)
  5. List of Lists (LIL) (# 2-1-List of Lists list-of-lists--lil)
  6. [Compressed Row Storage Method (CSR)](# 2-2-Compressed Row Storage Method compressed-sparse-row-csr)
  7. [Compressed column storage method (CSC)](# 2-3-Compressed column storage method compressed-sparse-column-csc)
  8. [Other](# 2-4-Other)
  9. [Four arithmetic operations, etc.](# 3-Four arithmetic operations, etc.)
  10. [Matrix Market format](# 4-matrix-market-format)
  11. [Direct Method Library](# 5-Direct Method Library)
  12. Iterative Library (# 6-Iterative Library)

1. Sparse matrix

1-1. Dense and sparse matrices

Dense procession. A matrix with few 0s in the matrix. Perhaps many people are more familiar with math problems.

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

A sparse matrix. It's just that there are a lot of 0s in the matrix. If you devise the part of the calculation related to 0, you can hardly calculate it, so calculate it by a special method. It will be dealt with in this 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. Visualization of sparse matrix

I think the best way to visualize sparse matrices is to use matplotlib's spy function. This is a matrix where the 0 part is displayed in white, and if there is any number, it is displayed in color. Assuming the above sparse matrix is a $ 100 \ times 100 $ matrix, it looks like this:

100100example.png

2. Sparse matrix storage format

In order to handle sparse matrices ignoring 0, we use a special method to represent sparse matrices. By doing this, it is possible to significantly reduce the memory used and the amount of calculation. You may not be sure, so I will explain using the matrix shown below. However, 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. List of Lists (LIL)

Stores tuples of (columns, elements) for each row. When storing data in the matrix for the first time, it is easier to use the LIL method. After that, I think it should be easy to understand to convert to the CRS or CCS method shown below. Once you get used to working with Python lists and Numpy arrays, it's not that difficult. In the above example,

line (Column,element)
0 (0, a)
1 (0, b)
2 (3, c), (5, d)
4 (4, e)
5 (5, f)

It is stored in the form of, and you can see that a sparse matrix can be described with a small amount of memory. The larger the sparse matrix, the better the storage method that ignores these zeros.

scipy.sparse.lil_matrix

Function lil_matrix for handling LIL format provided by scipy. It is often used when generating sparse matrices because it makes it easy to input non-zero elements of a matrix.

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())  #Return to normal matrix form
print("LIL type: ")
print(a.rows)  #In what column each row contains a nonzero element
print(a.data)  #Non-zero element in each row

The execution result is as follows.

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

Some call it CRS (Compressed Row Storage). Search in the row direction and store the elements of the nonzero matrix. If it is difficult to understand, please read the Material by Professor Nakajima of the University of Tokyo. An example of executing the CSR method by excluding diagonal components is shown.

If you use the CSR method,

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

Is

It can be expressed by three lists of. Again, you can see that you can write a sparse matrix with less memory. I think that this CSR format is often used when performing sparse matrix calculations. How quickly the calculation is completed is described in Chapter 3 of the four arithmetic operations.

To show what these lists represent, indices is the subscript of the elements of the original matrix, and indptr is on the right.

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

To explain the idea using this example,

scipy.sparse.csr_matrix

Function csr_matrix for handling the CSR format provided by scipy. If you put an arbitrary matrix mtx like scipy.sparse.csr_matrix (mtx), it will be converted to CSR format. The format of mtx is list, numpy.ndarray, numpy.matrix, etc. When changing from LIL format to CSR format, you can use 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:The column number of the element is displayed in order from the top row.
[1. 2. 3. 4. 5. 6.]  # data:The non-zero elements are shown according to indices.
[0 1 2 4 4 5 6]  # indptr:Count how many elements there are in each row and add them up. Number of lines because it starts from 0+It becomes the size of 1.

2-3. Compressed Sparse Column (CSC)

Some call it CCS. Search in the column direction and store the elements of the nonzero matrix. Column version of CSR method.

scipy.sparse.csc_matrix

Function csc_matrix for handling CSC format provided by scipy. If you put an arbitrary matrix mtx like scipy.sparse.csc_matrix (mtx), it will be converted to CSC format. The format of mtx is list, numpy.ndarray, numpy.matrix, etc. When changing from LIL format to CSC format, you can use 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:The row of the element is displayed in order from the left column.
[1. 2. 3. 5. 4. 6.]  # data:The non-zero elements are shown according to indices.
[0 2 2 2 3 4 6]  # indptr:Count how many elements there are in each column and add them up. Since it starts from 0, the number of columns+It becomes the size of 1.

2-4. Other

There are various other matrix storage methods, so if you want to know more, please see [Wikipedia](https://ja.wikipedia.org/wiki/sparse matrix).

3. Four arithmetic operations, etc.

When calculating a sparse matrix, it is most common to store it in the CSR or CSC method before performing the calculation. It can be performed in the same way as the familiar numpy arithmetic operations.

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

Multiplication between sparse matrices stored in CSR format can also be calculated like this.

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

Based on these, let's check the effectiveness of the CSR (compressed row storage) method. Try running the following code on your Jupyter notebook. You can see that the calculation is completed considerably (about 1000 times) faster with the ** CSR method than with the solution using an ordinary 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 method
%timeit a_large_dense.dot(b)  #Ordinary matrix

important point

If you put a csr format matrix in the argument of the dot function of numpy like 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>]]

I can't calculate well.

4. Matrix Market format

Although not used in this article, I will also mention one of the representations of sparse matrices, the Matrix Market format (https://math.nist.gov/MatrixMarket/formats.html). Sparse matrix data is basically distributed in this format or Rutherford Boeing format, and it is necessary knowledge to refer to papers and implementations related to sparse matrix calculation. The SuiteSparse Matrix Collection formerly the University of Florida Sparse Matrix Collection (https://sparse.tamu.edu) is probably the most famous dataset for sparse matrix problems.

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

Is expressed in Matrix Market (MM) format

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

It will be. You can see it in the comments above. This article is easy to understand. For the time being, I will play with a sparse matrix called steam3 that I chose appropriately. Use scipy.io when reading the Matrix market format.

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. Direct method library

It's in scipy.sparse.linalg.dsolve. I implemented the problem of Ax = b (see the formula below for details) solved in Previous article as an example (implementation has not changed much). ).

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

The default setting is scipy.sparse.linalg.dsolve (csr format A matrix, b vector) with SuperLU Decomposition (http://www.turbare.net/transl/scipy-lecture-notes/advanced/ scipy_sparse / solvers.html) should be done.

The image below shows the comparison result with the Gauss-Seidel method of the stationary iterative method.

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 library

A library of non-stationary iterative methods is implemented in scipy. It is in scipy.sparse.linalg.isolve. If you are not familiar with the transient iterative method, see Previous article. Some of the libraries that can be used are listed below. If you want to use a library other than these, look for it in Scipy's Manual (https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#module-scipy.sparse.linalg). Please give me.

  1. Symmetric matrix
  2. cg Conjugate Gradient method (CG method) --Only definite symmetric matrix
  3. minres minimum residual method (MINimum RESidual: MINRES method)
  4. Asymmetric matrix
  5. bicg bi-conjugate Gradient (BiCG method)
  6. cgs square conjugate gradient method (Conjugate Residual Squared: CGS method)
  7. bicgstab biconjugate gradient stabilization method (BiCG Stabilization: BiCG STAB method)
  8. gmres Generalized Minimal Residual (GMRES method)
  9. qmr quasi-minimal resident method (QMR method)

For the same problem used in the direct method, the implementation example in the BiCG STAB method and its execution result are shown below.

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)

Summary

I think it's like that.

References

  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] Articles that enable sparse matrix calculations at high speed
Articles that enable system development with Django (Python) _Introduction
Perform half-width / full-width conversion at high speed with Python
[Python] Find Fibonacci numbers at high speed (memoization, dynamic programming)
[Python] How to get divisors of natural numbers at high speed
A script that downloads AWS RDS log files at high speed
How to speed up Python calculations
A small story that outputs table data in CSV format at high speed