Производительность: Matlab против Python
Я недавно перешел с Matlab в Python, Преобразуя один из моих длинных кодов, я с удивлением обнаружил, Python быть очень медленным Я профилировал и проследил проблему с одной функцией, тянущей время. Эта функция вызывается из разных мест в моем коде (будучи частью других функций, которые вызываются рекурсивно). Профилировщик предполагает, что в обе функции было сделано 300 вызовов Matlab а также Python,
Вкратце, следующие коды кратко описывают проблему:
MATLAB
Класс, содержащий функцию:
classdef ExampleKernel1 < handle
methods (Static)
function [kernel] = kernel_2D(M,x,N,y)
kernel = zeros(M,N);
for i= 1 : M
for j= 1 : N
% Define the custom kernel function here
kernel(i , j) = sqrt((x(i , 1) - y(j , 1)) .^ 2 + ...
(x(i , 2) - y(j , 2)) .^2 );
end
end
end
end
end
и скрипт для вызова test.m:
xVec=[
49.7030 78.9590
42.6730 11.1390
23.2790 89.6720
75.6050 25.5890
81.5820 53.2920
44.9680 2.7770
38.7890 78.9050
39.1570 33.6790
33.2640 54.7200
4.8060 44.3660
49.7030 78.9590
42.6730 11.1390
23.2790 89.6720
75.6050 25.5890
81.5820 53.2920
44.9680 2.7770
38.7890 78.9050
39.1570 33.6790
33.2640 54.7200
4.8060 44.3660
];
N=size(xVec,1);
kex1=ExampleKernel1;
tic
for i=1:300
K=kex1.kernel_2D(N,xVec,N,xVec);
end
toc
Дает вывод
clear all
>> test
Elapsed time is 0.022426 seconds.
>> test
Elapsed time is 0.009852 seconds.
PYTHON 3.4
Класс, содержащий функцию CustomKernels.py:
from numpy import zeros
from math import sqrt
class CustomKernels:
"""Class for defining the custom kernel functions"""
@staticmethod
def exampleKernelA(M, x, N, y):
"""Example kernel function A"""
kernel = zeros([M, N])
for i in range(0, M):
for j in range(0, N):
# Define the custom kernel function here
kernel[i, j] = sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2)
return kernel
и скрипт для вызова test.py:
import numpy as np
from CustomKernels import CustomKernels
from time import perf_counter
xVec = np.array([
[49.7030, 78.9590],
[42.6730, 11.1390],
[23.2790, 89.6720],
[75.6050, 25.5890],
[81.5820, 53.2920],
[44.9680, 2.7770],
[38.7890, 78.9050],
[39.1570, 33.6790],
[33.2640, 54.7200],
[4.8060 , 44.3660],
[49.7030, 78.9590],
[42.6730, 11.1390],
[23.2790, 89.6720],
[75.6050, 25.5890],
[81.5820, 53.2920],
[44.9680, 2.7770],
[38.7890, 78.9050],
[39.1570, 33.6790],
[33.2640, 54.7200],
[4.8060 , 44.3660]
])
N = xVec.shape[0]
kex1 = CustomKernels.exampleKernelA
start=perf_counter()
for i in range(0,300):
K = kex1(N, xVec, N, xVec)
print(' %f secs' %(perf_counter()-start))
Дает вывод
%run test.py
0.940515 secs
%run test.py
0.884418 secs
%run test.py
0.940239 secs
РЕЗУЛЬТАТЫ
Сравнивая результаты, кажется, Matlab примерно в 42 раза быстрее после clear all "вызывается, а затем в 100 раз быстрее, если скрипт запускается несколько раз без вызова" clear all ". По крайней мере, на порядок, если не на два порядка быстрее. Это очень удивительный результат для меня. Я ожидал, что результат будет наоборот.
Может кто-нибудь, пожалуйста, пролить свет на это?
Может кто-нибудь предложить более быстрый способ сделать это?
ПРИМЕЧАНИЕ
Я также пытался использовать numpy.sqrt что ухудшает производительность, поэтому я использую math.sqrt в Python,
РЕДАКТИРОВАТЬ
for Циклы для вызова функций являются чисто фиктивными. Они там просто для " симуляции " 300 обращений к функции. Как я описал ранее, функции ядра (kernel_2D в Matlab а также kex1 в Python) вызываются из разных мест в программе. Чтобы сделать проблему короче, я " симулирую " 300 вызовов, используя for петля. for Циклы внутри функций ядра являются существенными и неизбежными из-за структуры матрицы ядра.
РЕДАКТИРОВАТЬ 2
Вот большая проблема: https://github.com/drfahdsiddiqui/bbfmm2d-python
4 ответа
Вы хотите избавиться от тех, for петли. Попробуй это:
def exampleKernelA(M, x, N, y):
"""Example kernel function A"""
i, j = np.indices((N, M))
# Define the custom kernel function here
kernel[i, j] = np.sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2)
return kernel
Вы также можете сделать это с вещанием, которое может быть еще быстрее, но немного менее интуитивно понятным из MATLAB,
После дальнейшего расследования я обнаружил, что с помощью indices как указано в ответе еще медленнее.
Решение: использовать meshgrid
def exampleKernelA(M, x, N, y):
"""Example kernel function A"""
# Euclidean norm function implemented using meshgrid idea.
# Fastest
x0, y0 = meshgrid(y[:, 0], x[:, 0])
x1, y1 = meshgrid(y[:, 1], x[:, 1])
# Define custom kernel here
kernel = sqrt((x0 - y0) ** 2 + (x1 - y1) ** 2)
return kernel
Результат: очень, очень быстро, в 10 раз быстрее, чем indices подход. Я получаю времена, которые ближе к C.
Однако: Использование meshgrid с Matlab биения C а также Numpy будучи в 10 раз быстрее, чем оба.
Все еще задаюсь вопросом, почему!
Matlab использует коммерческую библиотеку MKL. Если вы используете бесплатный дистрибутив Python, проверьте, используется ли в Python MKL или другая высокопроизводительная библиотека blas или это библиотека по умолчанию, которая может быть намного медленнее.
Сравнение Jit-компиляторов
Было упомянуто, что Matlab использует внутренний Jit-компилятор, чтобы получить хорошую производительность при выполнении таких задач. Давайте сравним jit-компилятор Matlabs с jit-компилятором Python (Numba).
Код
import numba as nb
import numpy as np
import math
import time
#If the arrays are somewhat larger it makes also sense to parallelize this problem
#cache ==True may also make sense
@nb.njit(fastmath=True)
def exampleKernelA(M, x, N, y):
"""Example kernel function A"""
#explicitly declaring the size of the second dim also improves performance a bit
assert x.shape[1]==2
assert y.shape[1]==2
#Works with all dtypes, zeroing isn't necessary
kernel = np.empty((M,N),dtype=x.dtype)
for i in range(M):
for j in range(N):
# Define the custom kernel function here
kernel[i, j] = np.sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2)
return kernel
def exampleKernelB(M, x, N, y):
"""Example kernel function A"""
# Euclidean norm function implemented using meshgrid idea.
# Fastest
x0, y0 = np.meshgrid(y[:, 0], x[:, 0])
x1, y1 = np.meshgrid(y[:, 1], x[:, 1])
# Define custom kernel here
kernel = np.sqrt((x0 - y0) ** 2 + (x1 - y1) ** 2)
return kernel
@nb.njit()
def exampleKernelC(M, x, N, y):
"""Example kernel function A"""
#explicitly declaring the size of the second dim also improves performance a bit
assert x.shape[1]==2
assert y.shape[1]==2
#Works with all dtypes, zeroing isn't necessary
kernel = np.empty((M,N),dtype=x.dtype)
for i in range(M):
for j in range(N):
# Define the custom kernel function here
kernel[i, j] = np.sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2)
return kernel
#Your test data
xVec = np.array([
[49.7030, 78.9590],
[42.6730, 11.1390],
[23.2790, 89.6720],
[75.6050, 25.5890],
[81.5820, 53.2920],
[44.9680, 2.7770],
[38.7890, 78.9050],
[39.1570, 33.6790],
[33.2640, 54.7200],
[4.8060 , 44.3660],
[49.7030, 78.9590],
[42.6730, 11.1390],
[23.2790, 89.6720],
[75.6050, 25.5890],
[81.5820, 53.2920],
[44.9680, 2.7770],
[38.7890, 78.9050],
[39.1570, 33.6790],
[33.2640, 54.7200],
[4.8060 , 44.3660]
])
#compilation on first callable
#can be avoided with cache=True
res=exampleKernelA(xVec.shape[0], xVec, xVec.shape[0], xVec)
res=exampleKernelC(xVec.shape[0], xVec, xVec.shape[0], xVec)
t1=time.time()
for i in range(10_000):
res=exampleKernelA(xVec.shape[0], xVec, xVec.shape[0], xVec)
print(time.time()-t1)
t1=time.time()
for i in range(10_000):
res=exampleKernelC(xVec.shape[0], xVec, xVec.shape[0], xVec)
print(time.time()-t1)
t1=time.time()
for i in range(10_000):
res=exampleKernelB(xVec.shape[0], xVec, xVec.shape[0], xVec)
print(time.time()-t1)
Спектакль
exampleKernelA: 0.03s
exampleKernelC: 0.03s
exampleKernelB: 1.02s
Matlab_2016b (your code, but 10000 rep., after few runs): 0.165s
Я получил ~5-кратное улучшение скорости по сравнению с сеткой, используя только трансляцию:
def exampleKernelD(M, x, N, y):
return np.sqrt((x[:,1:] - y[:,1:].T) ** 2 + (x[:,:1] - y[:,:1].T) ** 2)