Speeding-up builtin Matlab functions – part 2

Last week I showed how we can speed-up built-in Matlab functions, by creating local copies of the relevant m-files and then optimizing them for improved speed using a variety of techniques.Accelerating MATLAB Performance Today I will show another example of such speed-up, this time of the Financial Toolbox’s maxdrawdown function, which is widely used to estimate the relative risk of a trading strategy or asset. One might think that such a basic indicator would be optimized for speed, but experience shows otherwise. In fact, this function turned out to be the main run-time performance hotspot for one of my clients. The vast majority of his code was optimized for speed, and he naturally assumed that the built-in Matlab functions were optimized as well, but this was not the case. Fortunately, I was able to create an equivalent version that was 30-40 times faster (!), and my client remains a loyal Matlab fan.

In today’s post I will show how I achieved this speed-up, using different methods than the ones I showed last week. A combination of these techniques can be used in a wide range of other Matlab functions. Additional speed-up techniques can be found in other performance-related posts on this website, as well in my book Accelerating MATLAB Performance.

Profiling

As I explained last week, the first step in speeding up any function is to profile its current run-time behavior using Matlab’s builtin Profiler tool, which can either be started from the Matlab Editor toolstrip (“Run and Time”) or via the profile function.

The profile report for the client’s function showed that it had two separate performance hotspots:

  1. Code that checks the drawdown format (optional 2nd input argument) against a set of allowed formats
  2. Main code section that iteratively loops over the data-series values to compute the maximal drawdown

In order top optimize the code, I copied %matlabroot%/toolbox/finance/finance/maxdrawdown.m to a local folder on the Matlab path, renaming the file (and the function) maxdrawdn, in order to avoid conflict with the built-in version.

Optimizing input args pre-processing

The main problem with the pre-processing of the optional format input argument is the string comparisons, which are being done even when the default format is used (which is by far the most common case). String comparison are often more expensive than numerical computations. Each comparison by itself is very short, but when maxdrawdown is run in a loop (as it often is), the run-time adds up.

Here’s a snippet of the original code:

if nargin < 2 || isempty(Format)
    Format = 'return';
end
if ~ischar(Format) || size(Format,1) ~= 1
    error(message('finance:maxdrawdown:InvalidFormatType'));
end
choice = find(strncmpi(Format,{'return','arithmetic','geometric'},length(Format)));
if isempty(choice)
    error(message('finance:maxdrawdown:InvalidFormatValue'));
end

An equivalent code, which avoids any string processing in the common default case, is faster, simpler and more readable:

if nargin < 2 || isempty(Format)
    choice = 1;
elseif ~ischar(Format) || size(Format,1) ~= 1
    error(message('finance:maxdrawdown:InvalidFormatType'));
else
    choice = find(strncmpi(Format,{'return','arithmetic','geometric'},length(Format)));
    if isempty(choice)
        error(message('finance:maxdrawdown:InvalidFormatValue'));
    end
end

The general rule is that whenever you have a common case, you should check it first, avoiding unnecessary processing downstream. Moreover, for improved run-time performance (although not necessarily maintainability), it is generally preferable to work with numbers rather than strings (choice rather than Format, in our case).

Vectorizing the main loop

The main processing loop uses a very simple yet inefficient iterative loop. I assume that the code was originally developed this way in order to assist debugging and to ensure correctness, and that once it was ready nobody took the time to also optimize it for speed. It looks something like this:

MaxDD = zeros(1,N);
MaxDDIndex = ones(2,N);
...
if choice == 1   % 'return' format
    MaxData = Data(1,:);
    MaxIndex = ones(1,N);
    for i = 1:T
        MaxData = max(MaxData, Data(i,:));
        q = MaxData == Data(i,:);
        MaxIndex(1,q) = i;
        DD = (MaxData - Data(i,:)) ./ MaxData;
        if any(DD > MaxDD)
            p = DD > MaxDD;
            MaxDD(p) = DD(p);
            MaxDDIndex(1,p) = MaxIndex(p);
            MaxDDIndex(2,p) = i;
        end
    end
else             % 'arithmetic' or 'geometric' formats
    ...

This loop can relatively easily be vectorized, making the code much faster, and arguably also simpler, more readable, and more maintainable:

if choice == 3
    Data = log(Data);
end
MaxDDIndex = ones(2,N);
MaxData = cummax(Data);
MaxIndexes = find(MaxData==Data);
DD = MaxData - Data;
if choice == 1	% 'return' format
    DD = DD ./ MaxData;
end
MaxDD = cummax(DD);
MaxIndex2 = find(MaxDD==DD,1,'last');
MaxIndex1 = MaxIndexes(find(MaxIndexes<=MaxIndex2,1,'last'));
MaxDDIndex(1,:) = MaxIndex1;
MaxDDIndex(2,:) = MaxIndex2;
MaxDD = MaxDD(end,:);

Let’s make a short run-time and accuracy check – we can see that we achieved a 31-fold speedup (YMMV), and received the exact same results:

>> data = rand(1,1e7);
 
>> tic, [MaxDD1, MaxDDIndex1] = maxdrawdown(data); toc  % builtin Matlab function
Elapsed time is 7.847140 seconds.
 
>> tic, [MaxDD2, MaxDDIndex2] = maxdrawdn(data); toc  % our optimized version
Elapsed time is 0.253130 seconds.
 
>> speedup = round(7.847140 / 0.253130)
speedup =
    31
 
>> isequal(MaxDD1,MaxDD2) && isequal(MaxDDIndex1,MaxDDIndex2)  % check accuracy
ans =
  logical
   1

Disclaimer: The code above seems to work (quite well in fact) for a 1D data vector. You’d need to modify it a bit to handle 2D data – the returned maximal drawdown are still computed correctly but not the returned indices, due to their computation using the find function. This modification is left as an exercise for the reader…

Very similar code could be used for the corresponding maxdrawup function. Although this function is used much less often than maxdrawdown, it is in fact widely used and very similar to maxdrawdown, so it is surprising that it is missing in the Financial Toolbox. Here is the corresponding code snippet:

% Code snippet for maxdrawup (similar to maxdrawdn)
MaxDUIndex = ones(2,N);
MinData = cummin(Data);
MinIndexes = find(MinData==Data);
DU = Data - MinData;
if choice == 1	% 'return' format
    DU = DU ./ MinData;
end
MaxDU = cummax(DU);
MaxIndex = find(MaxDD==DD,1,'last');
MinIndex = MinIndexes(find(MinIndexes<=MaxIndex,1,'last'));
MaxDUIndex(1,:) = MinIndex;
MaxDUIndex(2,:) = MaxIndex;
MaxDU = MaxDU(end,:);

Similar vectorization could be applied to the emaxdrawdown function. This too is left as an exercise for the reader…

Conclusions

Matlab is composed of thousands of internal functions. Each and every one of these functions was meticulously developed and tested by engineers, who are after all only human. Whereas supreme emphasis is always placed with Matlab functions on their accuracy, run-time performance sometimes takes a back-seat. Make no mistake about this: code accuracy is almost always more important than speed (an exception are cases where some accuracy may be sacrificed for improved run-time performance). So I’m not complaining about the current state of affairs.

But when we run into a specific run-time problem in our Matlab program, we should not despair if we see that built-in functions cause slowdown. We can try to avoid calling those functions (for example, by reducing the number of invocations, or limiting the target accuracy, etc.), or optimize these functions in our own local copy, as I’ve shown last week and today. There are multiple techniques that we could employ to improve the run time. Just use the profiler and keep an open mind about alternative speed-up mechanisms, and you’d be half-way there.

Let me know if you’d like me to assist with your Matlab project, either developing it from scratch or improving your existing code, or just training you in how to improve your Matlab code’s run-time/robustness/usability/appearance. I will be visiting Boston and New York in early June and would be happy to meet you in person to discuss your specific needs.

Categories: Low risk of breaking in future versions, Stock Matlab function, Toolbox

Tags: , ,

Bookmark and SharePrint Print

6 Responses to Speeding-up builtin Matlab functions – part 2

  1. John says:

    Hi Yair

    just wondering if the find function is really required
    for example is: MaxIndexes = MaxData==Data; equivalent to MaxIndexes = find(MaxData==Data);

    best regards

  2. Marshall says:

    Hi Yair,

    One of the lines in your code

    MaxIndex2 = find(MaxDD==DD,1,'last');

    is one that is incredibly common but is undoubtedly inefficient because the == compares every element, despite the fact that our call to find requests that we stop. As an example, the following two are roughly comparable in time:

    % compare using vectorization
    X = zeros(1e9,1); X(5e8) = 9;
    tic; find(X==9,1); t_vec = toc;
     
    % compare using a for-loop
    tic;
    for i = 1:length(X)
        if(X(i)==9), break; end
    end
    t_iter = toc;
     
    fprintf('Vectorized: %2.2f\n',t_vec);
    fprintf('Iterative: %2.2f\n',t_iter);
     
    Vectorized: 3.94
    Iterative: 3.58

    However, in the following scenario, the issue becomes obvious:

    X = zeros(1e9,1); X(1) = 9;
    tic; find(X==9,1); t_vec = toc;
     
    % compare using a for-loop
    tic;
    for i = 1:length(X)
        if(X(i)==9), break; end
    end
    t_iter = toc;
     
    fprintf('Vectorized: %2.2f\n',t_vec);
    fprintf('Iterative: %2.2f\n',t_iter);
     
    Vectorized: 2.19
    Iterative: 0.00

    From this we can deduce that the call to find() in the first example takes about 2 seconds, with the other 2 seconds due to the comparison X==9. However, in the second example, the entirety of X is still compared to 9 before the find is called

    A smart compiler would optimize away the common code pattern of find(X==y), but it appears Matlab’s JIT doesn’t. Is there a more obvious idiomatic Matlab method of avoiding performing a full comparison when it’s not necessary?

    • @Marshal – I accept that additional speedup is possible for the line you noted, in cases involving 1e9 data elements (as in your example). As I noted in my post, the new function is ~40 times faster than Matlab’s built-in function, to the point where additional speedups might not have been very cost-effective, especially since my use-case involved up to a few million (not billions of) data elements. In such cases, shaving off a few dozen millisecs did not seem to be worth the extra coding time and risk.

      If you are interested in extra speedup, an additional tip would be to enclose the entire section that computes MaxIndex1, MaxIndex2 and MaxDDIndex with a check that the MaxDDIndex output arg is in fact requested – in many use cases it is not, and in such cases there is no need to compute these indexes:

      if nargout > 1   MaxIndex2 = find(MaxDD==DD,1,'last');
         MaxIndex1 = MaxIndexes(find(MaxIndexes<=MaxIndex2,1,'last'));
         MaxDDIndex(1,:) = MaxIndex1;
         MaxDDIndex(2,:) = MaxIndex2;
      end
    • Marshall says:

      Thanks Yair–I wasn’t trying to say “look you could have sped this up!”, I was instead asking if there was a better way of performing the loop I made in the second example, using vectorization instead of a for loop, which is the un-Matlab way to do it. I was mainly disappointed that Matlab wasn’t optimizing away vectorization of X==9 in conjunction with find(X==9,1).

Leave a Reply

Your email address will not be published. Required fields are marked *