[PYTHON] Performance comparison between 2D matrix calculation and for with numpy

In this article, how do you write two-dimensional matrix calculations in Python numpy like this? That is.

for i in range(0, ni-1):
    for j in range(0, nj-1):
        y[i, j] = x[i, j] + x[i+1, j] - x[i+1, j+1] * x[i, j+1]

Later, we will compare the performance. The Python version is 3.7.4. I am using anaconda3-2019.10.

Preface

COVID-19 + As I changed jobs, I decided to spare time and learn a new language every day when I kept getting stuck at my home in Tokyo. Python seems to be the main thing in the next workplace, so Python. I have to cut my money to get serious about it, so I took a course at Udemy and spent about a week learning the basics. I've been programming for a long time, but Python is a beginner with about two weeks of experience.

After trying it all, I feel like I can do web development with Python, but since I'm doing Python, I'd like to do data analysis and machine learning. So, I rewrote the solution of the one-dimensional advection equation that I did long ago in Python.

▶ 1 Dimensional Flow Calculation by Python

I haven't checked the details, and I've cut corners, but I was satisfied with the results I wanted for the first time. But what a sense of accomplishment? I couldn't feel something like that, so I tried two-dimensional calculation, but the matrix calculation was extremely dull. My Ghost whispers that I shouldn't continue writing anymore. And In python

--If you write for, you lose --Nesting for slows down execution. Cruel.

It's just a rumor. As a solution, I heard that numpy etc. will be used, so I tried using it, but it is just confusing and it seems that I will not be able to learn it immediately, so I will leave it as a memo.

All you can get by reading this article is that you don't have to read "Use numpy to calculate matrices (don't use for)". Other than that, I don't think it will come to your mind even if you read it, so if you understand it, you should actually move your hand to check it and repeat it until it fits in your hand.

Comparison of how to write list and numpy array

Python has `list, tuple, dict```, etc. to represent arrays, and for two-dimensional arrays, list``` comes from other languages. It is easy to understand if it is `` x [i] [j] etc. However, from the mathematical notation, `` `x [i, j] is easier to understand. In the array of numpy, it is written like the latter. As you can see by actually writing it, writing `x [i] [j] ``` many times is annoying and hard to see, so like `x [i, j] ``` Thank you for being able to write.

Then immediately. For example, suppose you have the following formula that calculates x, y in the form ndarray.

y[i, j] = x[i, j] + x[i+1, j] - x[i+1, j+1] * x[i, j+1]

It is not a particularly meaningful formula, but you should often see this kind of calculation not only in numerical systems but also in numerical calculation processing. Nesting this in a loop of i and `` `j```,

for i in range(0, x.shape[0]-1):
    for j in range(0, x.shape[1]-1):
        y[i, j] = x[i, j] + x[i+1, j] - x[i+1, j+1] * x[i, j+1]

Will be. Very easy to understand. If the maximum number of divisions is `ni```, `nj```, etc.

for i in range(0, ni-1):
    for j in range(0, nj-1):
        y[i, j] = x[i, j] + x[i+1, j] - x[i+1, j+1] * x[i, j+1]

It means.

By the way, in numpy, it seems smart to write this calculation as follows.

y[0:-1, 0:-1] = x[0:-1, 0:-1] + x[1:, 0:-1] - x[1:, 1:] * x[0:-1, 1:]

e? What?

If you take a closer look, what is it? Does the expression correspond to the sliced range? I don't feel like that. I don't feel like writing freely.

Try

Let's check if it really matches.

import numpy as np


x = np.random.randint(10, size=(100, 50))
y = np.zeros((100, 50))

for i in range(0, x.shape[0]-1):
    for j in range(0, x.shape[1]-1):
        y[i, j] = x[i, j] + x[i+1, j] - x[i+1, j+1] * x[i, j+1]

z = np.zeros_like(y)
z[0:-1, 0:-1] = x[0:-1, 0:-1] + x[1:, 0:-1] - x[1:, 1:] * x[0:-1, 1:]

if (z == y).all():
    print('Matches')
else:
    print('Wrong')

The execution result is ... correct. 0 in range is not necessary, but I dare to understand it. Other

--From 0 to ni, nj (` `x.shape [0], x.shape [1]` `) -- 1to ni, nj(`` `x.shape [0], x.shape [1]` `) --From 1 to` `ni-1, nj-1 (`` x.shape [0] -1, x.shape [1] -1```) --From ``` 2``` to ni-3, nj-3``` ( `x.shape [0] -3, x.shape [1] -3```)

And so on.

#0 to ni, nj (x.shape[0], x.shape[1])
for i in range(0, x.shape[0]):
    for j in range(x.shape[1]):
        y[i, j] = x[i, j] + x[i, j]

z[0:, 0:] = x[0:, 0:] + x[0:, 0:]

#1 to ni, nj (x.shape[0], x.shape[1])
for i in range(1, x.shape[0]):
    for j in range(1, x.shape[1]):
        y[i, j] = x[i, j] + x[i-1, j] - x[i-1, j-1] * x[i, j-1]

z[1:, 1:] = x[1:, 1:] + x[0:-1, 1:] - x[0:-1, 0:-1] * x[1:, 0:-1]

#1 to ni-1, nj-1 (x.shape[0]-1, x.shape[1]-1)
for i in range(1, x.shape[0]-1):
    for j in range(1, x.shape[1]-1):
        y[i, j] = x[i, j] + x[i-1, j] - x[i+1, j-1] * x[i, j+1]

z[1:-1, 1:-1] = x[1:-1, 1:-1] + x[0:-2, 1:-1] - x[2:, 0:-2] * x[1:-1, 2:]

#2 to ni-3, nj-3 (x.shape[0]-3, x.shape[1]-3)
for i in range(2, x.shape[0]-3):
    for j in range(2, x.shape[1]-3):
        y[i, j] = x[i, j] + x[i-1, j] - x[i+1, j-1] * x[i, j+1]

z[2:-3, 2:-3] = x[2:-3, 2:-3] + x[1:-4, 2:-3] - x[3:-2, 1:-4] * x[2:-3, 3:-2]

Summary of how to write

Hmmm, is that something like this?

# i,j from n to ni-m, nj-Loop to m
# x[i, j], x[i-ln, j], x[i+ln, j-ln], x[i, j+ln]And so on

for i in range(n, x.shape[0]-m):
    for j in range(n, x.shape[1]-m):
        y[i, j] = x[i, j] + x[i-ln, j] - x[i+ln, j-ln] * x[i, j+ln]

z[n:-m, n:-m] = x[n:-m, n:-m] + x[n-ln:-m-ln, n:-m] - \
    x[n+ln:-m+ln, n-ln:-m-ln] * x[n:-m, n+ln:-m+ln

A lot of summary clutter. It's extremely difficult to read, so I'll make it a table.

※ ni = x.shape[0]

Is it correct? It's hard to understand no matter how you write it. As a Python beginner, my frank impression is that it should not be used from the standpoint of maintainability. At least at this point, I thought so.

Comparison of processing speed

I looped and measured so that the calculation would rotate 100 million times. In the first one-dimensional calculation, the mesh was made finer and calculated about 20 million times, so it seems that there are many 100 million times, and in the actual calculation there are not many at all. It was troublesome to match the fractions, so leave it as it is. The measurement result is the average of 10 executions each.

Overwhelming joiko! !! !! Pepe

It is not at the level of maintainability. Roughly finishing at about 1% means that a calculation that takes an hour will finish in 36 seconds. What about one day? For a week? ?? 1 month? ?? ?? For readability, there is absolutely no option to use for. Please write in the comments.

However, whether or not this numpy calculation is fast is another matter. It seems that numpy is implemented in C and Fortran, but it's faster to nest do loops in Fortran! I think there is something like that. Not verified.

that's all. I think there is something that should be done more like this, but once this is done.