Improving Matlab function within big simulation

I have a very big Matlab simulation project in my hands, which I wanted to optimize, since I'm running it many times to tune parameters and the like.

Using Matlab's profile I identified one function that is eating up most of my time, specifically the output(i,1) assignment line.

This function is called a LOT, where input is a 10x1 double passed as an argument, and output is also a 10x1 vector.

function output = my_function(input)

a = size(input,1);
output = input*0;
dens = density(input);

% for each i, output(i) is the maximum between output(i+1) and mean(output(i+1:end))
for i = 1:a-1
    output(i,1)= max(mean(dens(i+1:a,1)),dens(i+1,1));
end
output(a,1) = dens(a,1);

end

density() function:

function dens = density(input)

dens = 1000*(1 - (input+288.9414)./(508929.2.*(input+68.12963)).*(input-3.9863).^2);

end

My ideas:

• I think vectorization would maybe help to get rid of the loop, but I'm not familiar at all with the technique.

• Is there a faster/alternative way to calculate the mean (maybe without Matlab's built-in function call?)

Vectorization

If you look at all the calls to mean, you see that what you are calculating in each loop is a "Cumulative moving average", but in the oposite direction of what one usually would do.

A cumulative moving average can be calculated quite simply using the cumulative sum, cumsum, and dividing by the number of elements. Since you're starting with all elements and up with only one, we need to flip it.

The running mean, (or cumulative moving average) can be calculated like this:

• First flip dens upside down

• Take the cumulative sum of it

• Divide each element by the number of elements counted

• Flip it back up

Thus:

running_mean = flipud(cumsum(flipud(dens))./(1:numel(dens)).');

Now, we need to find which element is larger, the cumulative sum, or the original density variable. You are not including the first element of the cumulative sum, and you are appending dens(end) at the end of the vector.

The maximum value of each row of a matrix can be found using max(A,[],2), where the number two represent the second dimension (columns). So, the values we want to use are:

max([running_mean(2:end), dens(2:end)],[],2)

In addition, we need to append the last element, so:

output = [max([running_mean(2:end), dens(2:end)],[],2); dens(end)];

To summarize, the entire operation requires no loops, one call to cumsum, and one call to max. In addition we need to flip it twice, but that takes virtually no time.

In total, your function can be written like this:

function output = my_function(vector)

n = numel(vector); % numel, not length
output = zeros(n,1) % prior to R2015b, use: output(n,1) = 0

% Create an anonymous function called density:
density = @(vector) 1000(1 - (vector+288.9414)./(508929.2.(vector+68.12963)).*(vector-3.9863).^2);

% Calculate the density:
dens = density(vector);

% Calculate the running mean, and create the final output vector:
running_mean = flipud(cumsum(flipud(dens))./(1:numel(dens)).');
output = [max([running_mean(2:end), dens(2:end)],[],2); dens(end)];

end


Time in Octave shows the improvement, as a function of vector size. For vectors with only ten elements, the execution time is 25% of the original code.

Review of your code:

-
Don't use input as a variable name! input is used to ask users for a command line input. Using it as a variable name makes the function useless, but will also potentially create strange errors or warnings. (I called it vector in this review. You should pick a more descriptive name)

The same goes with max, mean, length etc. Don't use function names as variable names.

-
Your input is a vector, so no need to specify that you always want the first column. Use linear indexing instead of subscripts. It's much faster, and you won't mess up the dimensions (oops, it was a column vector, not a row vector).

-
Don't use size(vector,1), and don't use length(vector), use numel(vector). It's faster and more robust.

-
You are pre-allocating output. Very good! That's a common bottleneck.

-
output = vector*0; is clever. Prior to Matlab R2015b zeros(n,1) was actually far from the fastest way to create a vector of zeros! In fact, output(n,1) = 0 was up to 500 times faster! This is because zeros creates a matrix, then fills the elements with zeros, where as the latter just creates a matrix, where the values are zero by default.

-
In your loop, don't use i as the iterator, it is used for the imaginary unit (sqrt(-1)). It's common to use for instance ii, jj etc.

-
If your density function is only what you have posted there, then you can simply create an anonymous function inside the main function:

density = @(vector) 1000*(1 - (vector+288.9414)./(508929.2.*(vector+68.12963)).*(vector-3.9863).^2);

Call it like:

dens = density(vector);

Review of reviews:

There are a few things I don't like about the other reviews that I want to point out:

-
"You don't need to specify that 1 when indexing a vector, MATLAB knows". Many Matlab users use max(A), when they want to find the maximum value of all columns in a matrix. It's much simpler than max(A, [], 1). You might however come into trouble if A has varying number of rows (and matrices often do). Because the moment that matrix has only one row, max(A) will give the maximum value of the entire vector, instead of each (one row) column. This can give you a major headache!

-
"Why not use length?" Well, in this case, length is much better than size, but numel should still be the preferred choice. length is better than size for two reasons. Using size limits the function to only work on column vectors only, and will fail on row vectors. Second, numel is much faster than length for large vectors, and both are much faster than size for all vectors.

-
"Don't use the profiler, use tic/toc". True, the profiler includes the time it takes to run itself, but it's still the easiest wa

StackExchange Code Review Q#71788, answer score: 7