Я получил ~ 5-кратное улучшение скорости по сравнению с сеткой, используя только трансляцию:

авно перешел с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

 norok228 сент. 2017 г., 21:14
Учитывая, что у вас уже есть доступ к коду C ++ (согласно вашемуРЕДАКТИРОВАТЬ 2), Я хотел бы рассмотреть возможность создания привязок этого кода к Python, а не переводить его, если только вы не выполняете этот перевод по определенным причинам, помимо наличия алгоритма, доступного в Python.
 Martin Beckett28 сент. 2017 г., 19:34
Как правило, не пытайтесь перебрать массив в Python. Вызовите операции для всего массива (ов), используя numpy, чтобы фактическое вычисление для каждого элемента выполнялось внутри библиотеки.
 Daniel F28 сент. 2017 г., 19:35
Силаnumpy это способность избавиться от тех,for петли
 norok228 сент. 2017 г., 19:50
Если проблема заключается в петле, по которой вы вызываетеexampleKernelA функционировать 300 раз, вы должны, вероятно, рассмотретьnumba«s@jit, В целом, циклы в Python медленны по сравнению с скомпилированными языками точно в срок (или, конечно, раньше времени), такими как современныеMATLAB distrubutions.
 Fahd Siddiqui28 сент. 2017 г., 19:36
Я понимаю, что вы говорите, это верно и для Matlab. Но структура матрицы ядра в этом случае неизбежна. В любом случае, почему вызов функции настолько дорогой в Python и менее дорогой в Matlab?

Ответы на вопрос(5)

Сравнение 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
 Cris Luengo16 июн. 2018 г., 21:18
Поменяйте циклы OP, код MATLAB будет значительно быстрее. Кроме того, fastmath не следует показывать в этом сравнении.
 Cris Luengo16 июн. 2018 г., 22:54
да, вы правы, это небольшой массив, он, вероятно, помещается в кэш. Неважно. :)
 max911116 июн. 2018 г., 22:40
@Cris Luengo Я уже пытался переключать циклы без эффекта (возможно, из-за небольшого размера массива), я попробую это без fastmath и добавлю результаты. Для действительно честного сравнения следует использовать новейшую версию 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 раз быстрее, чем оба.

Все еще задаюсь вопросом, почему!

Решение Вопроса

Вы хотите избавиться от тех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.

 Luis Mendo29 сент. 2017 г., 19:46
@ Percusse Я не следую. Можете ли вы привести пример в Octave или Numpy, где вещание не для бинарного оператора (т.е. с двумя входами)?
 Luis Mendo28 сент. 2017 г., 22:14
Почему вещаниенемного менее интуитивно понятный из Matlab? Matlab имел трансляцию (с другим названием)с 2007 годаи это происходит неявнос 2017 года
 Fahd Siddiqui28 сент. 2017 г., 20:04
Благодарю. Предлагаемая форма не работает.line 41, in exampleKernelA kernel[i, j] = sqrt((x[i, 0] - y[j, 0]) ** 2 + (x[i, 1] - y[j, 1]) ** 2) TypeError: only length-1 arrays can be converted to Python scalars
 Andras Deak29 сент. 2017 г., 23:00
@percusse для разумного обсуждения этого вопроса, вы должны сначала определить вещание, потому что я должен согласиться с Луисом, что я не понимаю вашего различия. Кроме того, я не считаю, что вещание интуитивно понятно, если вы не понимаете, как ведет себя bsxfun.
 Daniel F29 сент. 2017 г., 07:30
Прости мой последнийMATLAB опыт есть. , Некоторое время назад Человек я чувствую себя старым сейчас.

используя только трансляцию:

def exampleKernelD(M, x, N, y):
    return np.sqrt((x[:,1:] - y[:,1:].T) ** 2 + (x[:,:1] - y[:,:1].T) ** 2)

платный дистрибутив Python, проверьте, используется ли в Python MKL или другая высокопроизводительная библиотека blas или это библиотека по умолчанию, которая может быть намного медленнее.

 max911116 июн. 2018 г., 17:37
MKL релевантно, если вызываются подпрограммы BLAS, что в данном примере не актуально. Здесь важен только Jit-компилятор.

Ваш ответ на вопрос