How to write vectorized functions in MATLAB

Question

I am just learning MATLAB, and it's hard for me to understand performance factors for loops and vectorized functions.

In my previous question: Nested for loops are very slow in MATLAB (preallocated) I realized that using a vectorized function against 4 nested loops made a 7x time difference .

In this example, instead of iterating over all dimensions of a 4-dimensional array and calculating the median for each vector, it was much cleaner and faster to just call the median (stack, n), where n denoted the working size of the median function.

But the median is a very simple example, and I was just lucky that he had this parameter implemented . .

My question is, how do you write a function yourself that works as efficiently as one that has this size range ?

For example, you have a function my_median_1D that only works with a one-dimensional vector and returns a number.

How do you write the function my_median_nD , which acts as the median of MATLAB, taking an n-dimensional array and the "working size" parameter?

Refresh

I found code to calculate the median in higher dimensions

 % 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 

Could you explain to me that why is this code so efficient compared to simple nested loops ? It has nested loops just like other functions.

I don’t understand how it can be 7 times faster , and also why it is so difficult .

Update 2

I realized that using the median is not a good example, as it is a complex function that requires sorting an array or other neat tricks. I repeated the tests with an average and the results are even crazier: 19 seconds versus 0.12 seconds. This means that the built-in method for the sum is 160 times faster than nested loops .

I find it very difficult to understand how such an extreme performance difference in an industry-leading language depends on the programming style, but I see the points mentioned in the answers below.

Answer 1

Update 2 (to answer your updated question)

MATLAB is optimized for working with arrays. Once you get used to it, it’s really nice to just type one line and MATLAB to do the whole 4D loop without worrying about it. MATLAB is often used for prototype / one-time calculations, so it makes sense to save time for human coding and abandon some C flexibility [++ | #].

This is why MATLAB internally does some loops really well - often by encoding them as a compiled function.

The code snippet you give doesn’t actually contain the corresponding line of code that does the bulk of the work, namely

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

In other words, the code you are showing requires only access to the median values ​​at their correct index in the now sorted multidimensional array x (which does not take much time). The actual work of accessing all elements of the array was done using sort , which is a built-in (i.e., compiled and optimized) function.

Original answer (on how to create your own fast functions that work with arrays)

In fact, there are quite a few built-in modules that take a measurement parameter: min(stack, [], n) , max(stack, [], n) , mean(stack, n) , std(stack, [], n) , median(stack,n) , sum(stack, n) ... while other built-in functions, such as exp() , sin() , automatically work on every element of your entire array (i.e. sin(stack) automatically executes four nested loops for you if stack is 4D), you can create many functions that you might need, just rely on existing built-in functions .

If this is not enough for a specific case, you should look at repmat , bsxfun , arrayfun and accumarray , which are very powerful functions for performing MATLAB path actions. Just do a SO search for questions (or rather answers) using one , I learned a lot about MATLABs strengths.

As an example , let's say you wanted to implement a p-norm stack of size n , you can write

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

... where you effectively reuse the "dimension" sum .

Update

As Max points out in the comments, also see the colon operator ( : , which is a very powerful tool for selecting elements from an array (or even changing its shape, which is more commonly done with reshape ).

In general, look in the Array Operations section in help - it contains repmat et al. mentioned above, but also cumsum and some more obscure auxiliary functions that you should use as building blocks.

Answer 2

Could you explain to me why this code is so efficient compared to simple nested loops? It has nested loops just like other functions.

The problem with nested loops is not nested loops. These are the operations you perform internally.

Each function call (especially for a non-built-in function) generates a bit of overhead; moreover, if the function performs, for example, error checking, which takes the same time regardless of the size of the input. Thus, if the function has only 1 ms overhead, if you call it 1000 times, you will lose the second. If you can call it once to perform a vectorized calculation, you pay overhead only once.

In addition, the JIT compiler (pdf) can help in vectorizing simple for-loops, where you, for example, perform only basic arithmetic operations. Thus, loops with simple calculations in your message are accelerated a lot, while loops that call median are not.

Answer 3

Vectorization

In addition to what has already been said, you should also understand that vectorization involves parallelization, that is, performing parallel operations with data as opposed to sequential execution (think about SIMD instructions), and even in some cases use threads and multiprocessors ...

Mex files

Now, although the “interpreted or compiled” argument has already been argued, no one has mentioned that you can extend MATLAB by writing MEX files that are compiled by executable files written in C that can be called directly as a normal function from within MATLAB. This allows you to implement performance-critical components using a lower-level language such as C.

Column order

Finally, when trying to optimize some code, always remember that MATLAB stores matrices in column order. Access to items in this order can provide significant improvements over other custom orders.

For example, in your previous related question, you calculated a median set of stacked images along a certain size. Now the order in which these sizes are ordered greatly affects performance. Illustration:

 %# 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 HW 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: [HW 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) 

The time difference was huge:

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

Answer 4

In this case

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

But with a generic function, you can try splitting your array with mat2cell (which can work with nD arrays, not just matrices) and apply your my_median_1D function via cellfun . Below I use median as an example to show that you are getting equivalent results, but instead you can pass it any function defined in the m file, or an anonymous function defined using the @(args) notation.

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