Как написать векторизованные функции в MATLAB

Question

Я только изучаю MATLAB, и мне трудно понять факторы производительности циклов по сравнению с векторизованными функциями.

В моем предыдущем вопросе: Вложено для циклов очень медленноMATLAB (предварительно выделенный) Я понял, что использование векторизованной функции и 4 вложенных циклов дает 7-кратное различие во времени выполнения .

В этом примере вместо циклического прохождения всех измерений4-мерный массив и вычисление медианы для каждого вектора, было намного чище и быстрее просто вызвать медиану (стек, n), где n означало рабочее измерение медианной функции.

Но медиана - это очень простопример и Мне просто повезло что у него реализован этот параметр измерения .

Мой вопрос заключается в том, как вы сами пишете функцию, котораяработает так же эффективно, как и тот, в котором реализован этот диапазон измерений ?

Например, у вас есть функция my_median_1D, которая работает толькона одномерном векторе и возвращает число.

Как написать функцию my_median_nD, которая действует как медиана MATLAB, принимая n-мерный массив и «рабочее измерение» параметр?

Обновление

Я нашел код для вычисления медианы в более высоких измерениях

% In all other cases, use linear indexing to determine exact location
% of medians.  Use linear indices to extract medians, then reshape at
% end to appropriate size.
cumSize = cumprod(s);
total = cumSize(end);            % Equivalent to NUMEL(x)
numMedians = total / nCompare;

numConseq = cumSize(dim - 1);    % Number of consecutive indices
increment = cumSize(dim);        % Gap between runs of indices
ixMedians = 1;

y = repmat(x(1),numMedians,1);   % Preallocate appropriate type

% Nested FOR loop tracks down medians by their indices.
for seqIndex = 1:increment:total
  for consIndex = half*numConseq:(half+1)*numConseq-1
    absIndex = seqIndex + consIndex;
    y(ixMedians) = x(absIndex);
    ixMedians = ixMedians + 1;
  end
end

% Average in second value if n is even
if 2*half == nCompare
  ixMedians = 1;
  for seqIndex = 1:increment:total
    for consIndex = (half-1)*numConseq:half*numConseq-1
      absIndex = seqIndex + consIndex;
      y(ixMedians) = meanof(x(absIndex),y(ixMedians));
      ixMedians = ixMedians + 1;
    end
  end
end

% Check last indices for NaN
ixMedians = 1;
for seqIndex = 1:increment:total
  for consIndex = (nCompare-1)*numConseq:nCompare*numConseq-1
    absIndex = seqIndex + consIndex;
    if isnan(x(absIndex))
      y(ixMedians) = NaN;
    end
    ixMedians = ixMedians + 1;
  end
end

Не могли бы вы объяснить мне, , почему этот код так эффективен по сравнению с простым вложеннымпетли ?Он имеет вложенные циклы, как и другие функции.

Я не понимаю как это может быть в 7 раз быстрее , а также, что почему это так сложно .

Обновление2

Я понял, что использование медианы не является хорошим примером, поскольку сама по себе это сложная функция, требующая сортировки массива или других простых приемов.Вместо этого я повторил тесты со средним значением, и результаты стали еще более сумасшедшими: 19 секунд против 0,12 секунды. Это означает, что , построенный для суммирования, в 160 раз быстрее, чем вложенные циклы .

Мне действительно трудно понять, как ли ведущий в отрасли язык может иметь такую ​​крайнюю разницу в производительности в зависимости от стиля программирования, но я вижу моменты, упомянутые в ответах ниже.

Answer 1

Обновление 2 (для ответа на ваш обновленный вопрос)

MATLAB оптимизирован для работы с массивами. Как только вы к этому привыкнете, на самом деле очень приятно просто набрать одну строку и сделать так, чтобы MATLAB самостоятельно выполнял весь цикл 4D, не беспокоясь об этом. MATLAB часто используется для прототипирования / однократных вычислений, поэтому имеет смысл сэкономить время на программирование человека и отказаться от некоторой гибкости C [++ | #].

Вот почему MATLAB внутренне выполняет некоторые циклы действительно хорошо - часто кодируя их как скомпилированную функцию.

Фрагмент кода, который вы даете , на самом деле не содержит соответствующей строки кода, которая выполняет основную работу, а именно

% Sort along given dimension
x = sort(x,dim);

Другими словами, код, который вы показываете, должен получить доступ только к медианным значениям по их правильному индексу в теперь отсортированном многомерном массиве x (который не занимает много времени). Фактическая работа с доступом ко всем элементам массива была выполнена с помощью sort, которая является встроенной (то есть скомпилированной и высоко оптимизированной) функцией.

Оригинальный ответ (о том, как создавать свои собственные быстрые функции, работающие с массивами )

На самом деле существует довольно много встроенных модулей, которые принимают параметр измерения: min(stack, [], n), max(stack, [], n), mean(stack, n), std(stack, [], n), median(stack,n), sum(stack, n) ... вместе с тем, что другие встроенные функции, такие как exp(), sin(), автоматически работают с каждым элементом всего массива (т. е. sin(stack) автоматически выполняет четыре вложенных цикла для вас, если stack равен 4D), вы можете создать много функций, которые могут вам понадобиться, просто полагаться на существующие встроенные модули .

Если этого недостаточно для конкретного случая, вам следует взглянуть на repmat, bsxfun, arrayfun и accumarray, которые являются очень мощными функциями для выполнения действий «путем MATLAB». Просто ищите в SO вопросы (или, скорее, ответы) , используя one из этих , таким образом я узнал много о сильных сторонах MATLAB .

В качестве примера , скажем, вы хотите реализовать p-норму стека вдоль измерения n, вы можете написать

function result=pnorm(stack, p, n)
result=sum(stack.^p,n)^(1/p);

... где вы эффективно повторно используете «способность к размеру» sum.

Обновление

Как отмечает Макс в комментариях, взгляните также на оператор двоеточия (:) , который является очень мощным инструментом для выбора элементов из массива (или даже изменения его формы, что более обычно делается с reshape).

В общем, посмотрите в разделе Операции с массивами в справке - он содержит repmat и др. упомянутый выше, но также cumsum и некоторые более неясные вспомогательные функции, которые вы должны использовать в качестве строительных блоков.

Answer 2

Векторизация

В дополнение к тому, что уже было сказано, вы также должны понимать, что векторизация включает в себя распараллеливание, то есть выполнение параллельных операций с данными в отличие от последовательного выполнения (например, инструкции SIMD), и даже использование преимуществ потоков и мультипроцессоры в некоторых случаях ...

MEX-файлы

Теперь, хотя точка «интерпретированный и скомпилированный» уже обсуждалась, никто не упомянул, что вы можете расширить MATLAB, написав MEX-файлы, которые являются скомпилированными исполняемыми файлами, написанными на C, которые можно вызывать напрямую как обычную функцию изнутри. MATLAB. Это позволяет вам реализовать компоненты, критичные к производительности, используя язык более низкого уровня, такой как C.

Столбец-майор порядка

Наконец, при попытке оптимизировать некоторый код, всегда помните, что MATLAB хранит матрицы в порядке старших столбцов. Доступ к элементам в этом порядке может привести к значительным улучшениям по сравнению с другими произвольными порядками.

Например, в своем предыдущем связанном вопросе вы вычисляли median набора сложенных изображений вдоль некоторого измерения. Теперь порядок, в котором упорядочены эти размеры, сильно влияет на производительность. Иллюстрация:

%# sequence of 10 images
fPath = fullfile(matlabroot,'toolbox','images','imdemos');
files = dir( fullfile(fPath,'AT3_1m4_*.tif') );
files = strcat(fPath,{filesep},{files.name}');      %'

I = imread( files{1} );

%# stacked images along the 1st dimension: [numImages H W RGB]
stack1 = zeros([numel(files) size(I) 3], class(I));
for i=1:numel(files)
    I = imread( files{i} );
    stack1(i,:,:,:) = repmat(I, [1 1 3]);   %# grayscale to RGB
end

%# stacked images along the 4th dimension: [H W RGB numImages]
stack4 = permute(stack1, [2 3 4 1]);

%# compute median image from each of these two stacks
tic, m1 = squeeze( median(stack1,1) ); toc
tic, m4 = median(stack4,4); toc
isequal(m1,m4)

Разница во времени была огромна:

Elapsed time is 0.257551 seconds.     %# stack1
Elapsed time is 17.405075 seconds.    %# stack4

Answer 3

Не могли бы вы объяснить мне, почему этот код так эффективен по сравнению с простыми вложенными циклами?Он имеет вложенные циклы, как и другие функции.

Проблема с вложенными циклами заключается не в самих вложенных циклах.Это операции, которые вы выполняете внутри.

Каждый вызов функции (особенно для не встроенной функции) создает небольшую нагрузку;более того, если функция выполняет, например, проверку ошибок, которая занимает одинаковое количество времени независимо от размера ввода.Таким образом, если функция имеет только 1 мс служебной информации, если вы вызовете ее 1000 раз, вы потеряете секунду.Если вы можете вызвать его один раз для выполнения векторизованного расчета, вы платите только один раз.

Кроме того, JIT-компилятор (pdf) может помочь векторизовать простые for-циклы, где вы, например, выполняете только основные арифметические операции.Таким образом, циклы с простыми вычислениями в вашем посте значительно ускоряются, а циклы, вызывающие median, - нет.

Answer 4

В этом случае

M = median(A,dim) returns the median values for elements along the dimension of A specified by scalar dim

Но с помощью общей функции вы можете попробовать разделить ваш массив с mat2cell (который может работать с массивами nD, а не только с матрицами) и применить вашу функцию my_median_1D черезcellfun.Ниже я буду использовать median в качестве примера, чтобы показать, что вы получаете эквивалентные результаты, но вместо этого вы можете передать его любой функции, определенной в m-файле, или анонимной функции, определенной с нотацией @(args).

>> testarr = [[1 2 3]' [4 5 6]']

testarr =

     1     4
     2     5
     3     6

>> median(testarr,2)

ans =

    2.5000
    3.5000
    4.5000

>> shape = size(testarr)

shape =

     3     2

>> cellfun(@median,mat2cell(testarr,repmat(1,1,shape(1)),[shape(2)]))

ans =

    2.5000
    3.5000
    4.5000