Vectorization and Faster Code Execution
***************************************

Vectorization is a programming technique that uses vector operations instead of element-by-element loop-based operations.  Besides frequently producing more succinct Octave code, vectorization also allows for better optimization in the subsequent implementation.  The optimizations may occur either in Octave's own Fortran, C, or C++ internal implementation, or even at a lower level depending on the compiler and external numerical libraries used to build Octave.  The ultimate goal is to make use of your hardware's vector instructions if possible or to perform other optimizations in software.

Vectorization is not a concept unique to Octave, but it is particularly important because Octave is a matrix-oriented language.  Vectorized Octave code will see a dramatic speed up (10X--100X) in most cases.

This chapter discusses vectorization and other techniques for writing faster code.

Basic Vectorization
===================

To a very good first approximation, the goal in vectorization is to write code that avoids loops and uses whole-array operations.  As a trivial example, consider

     for i = 1:n
       for j = 1:m
         c(i,j) = a(i,j) + b(i,j);
       endfor
     endfor

compared to the much simpler

     c = a + b;

This isn't merely easier to write; it is also internally much easier to optimize.  Octave delegates this operation to an underlying implementation which, among other optimizations, may use special vector hardware instructions or could conceivably even perform the additions in parallel.  In general, if the code is vectorized, the underlying implementation has more freedom about the assumptions it can make in order to achieve faster execution.

This is especially important for loops with "cheap" bodies.  Often it suffices to vectorize just the innermost loop to get acceptable performance.  A general rule of thumb is that the "order" of the vectorized body should be greater or equal to the "order" of the enclosing loop.

As a less trivial example, instead of

     for i = 1:n-1
       a(i) = b(i+1) - b(i);
     endfor

write

     a = b(2:n) - b(1:n-1);

This shows an important general concept about using arrays for indexing instead of looping over an index variable.  See Index Expressions. Also use boolean indexing generously.  If a condition needs to be tested, this condition can also be written as a boolean index.  For instance, instead of

     for i = 1:n
       if (a(i) > 5)
         a(i) -= 20
       endif
     endfor

write

     a(a>5) -= 20;

which exploits the fact that a > 5 produces a boolean index.

Use elementwise vector operators whenever possible to avoid looping (operators like .* and .^).  See Arithmetic Ops.  For simple inline functions, the vectorize function can do this automatically.

(vectorize)

Also exploit broadcasting in these elementwise operators both to avoid looping and unnecessary intermediate memory allocations. See Broadcasting.

Use built-in and library functions if possible.  Built-in and compiled functions are very fast.  Even with an m-file library function, chances are good that it is already optimized, or will be optimized more in a future release.

For instance, even better than

     a = b(2:n) - b(1:n-1);

is

     a = diff (b);

Most Octave functions are written with vector and array arguments in mind.  If you find yourself writing a loop with a very simple operation, chances are that such a function already exists.  The following functions occur frequently in vectorized code:

  * Index manipulation

       * find

       * sub2ind

       * ind2sub

       * sort

       * unique

       * lookup

       * ifelse / merge

  * Repetition

       * repmat

       * repelems

  * Vectorized arithmetic

       * sum

       * prod

       * cumsum

       * cumprod

       * sumsq

       * diff

       * dot

       * cummax

       * cummin

  * Shape of higher dimensional arrays

       * reshape

       * resize

       * permute

       * squeeze

       * deal

Broadcasting
============

Broadcasting refers to how Octave binary operators and functions behave when their matrix or array operands or arguments differ in size.  Since version 3.6.0, Octave now automatically broadcasts vectors, matrices, and arrays when using elementwise binary operators and functions. Broadly speaking, smaller arrays are ``broadcast'' across the larger one, until they have a compatible shape.  The rule is that corresponding array dimensions must either

 1. be equal, or

 2. one of them must be 1.

In case all dimensions are equal, no broadcasting occurs and ordinary element-by-element arithmetic takes place.  For arrays of higher dimensions, if the number of dimensions isn't the same, then missing trailing dimensions are treated as 1.  When one of the dimensions is 1, the array with that singleton dimension gets copied along that dimension until it matches the dimension of the other array.  For example, consider

     x = [1 2 3;
          4 5 6;
          7 8 9];

     y = [10 20 30];

     x + y

Without broadcasting, x + y would be an error because the dimensions do not agree.  However, with broadcasting it is as if the following operation were performed:

     x = [1 2 3
          4 5 6
          7 8 9];

     y = [10 20 30
          10 20 30
          10 20 30];

     x + y
     =>    11   22   33
           14   25   36
           17   28   39

That is, the smaller array of size [1 3] gets copied along the singleton dimension (the number of rows) until it is [3 3].  No actual copying takes place, however.  The internal implementation reuses elements along the necessary dimension in order to achieve the desired effect without copying in memory.

Both arrays can be broadcast across each other, for example, all pairwise differences of the elements of a vector with itself:

     y - y'
     =>    0   10   20
         -10    0   10
         -20  -10    0

Here the vectors of size [1 3] and [3 1] both get broadcast into matrices of size [3 3] before ordinary matrix subtraction takes place.

A special case of broadcasting that may be familiar is when all dimensions of the array being broadcast are 1, i.e., the array is a scalar.  Thus for example, operations like x - 42 and max (x, 2) are basic examples of broadcasting.

For a higher-dimensional example, suppose img is an RGB image of size [m n 3] and we wish to multiply each color by a different scalar.  The following code accomplishes this with broadcasting,

     img .*= permute ([0.8, 0.9, 1.2], [1, 3, 2]);

Note the usage of permute to match the dimensions of the [0.8, 0.9, 1.2] vector with img.

For functions that are not written with broadcasting semantics, bsxfun can be useful for coercing them to broadcast.

(bsxfun)

Broadcasting is only applied if either of the two broadcasting conditions hold.  As usual, however, broadcasting does not apply when two dimensions differ and neither is 1:

     x = [1 2 3
          4 5 6];
     y = [10 20
          30 40];
     x + y

This will produce an error about nonconformant arguments.

Besides common arithmetic operations, several functions of two arguments also broadcast.  The full list of functions and operators that broadcast is

           plus      +  .+
           minus     -  .-
           times     .*
           rdivide   ./
           ldivide   .\
           power     .^  .**
           lt        <
           le        <=
           eq        ==
           gt        >
           ge        >=
           ne        !=  ~=
           and       &
           or        |
           atan2
           hypot
           max
           min
           mod
           rem
           xor

           +=  -=  .+=  .-=  .*=  ./=  .\=  .^=  .**=  &=  |=

Beware of resorting to broadcasting if a simpler operation will suffice. For matrices A and B, consider the following:

     C = sum (permute (A, [1, 3, 2]) .* permute (B, [3, 2, 1]), 3);

This operation broadcasts the two matrices with permuted dimensions across each other during elementwise multiplication in order to obtain a larger 3-D array, and this array is then summed along the third dimension. A moment of thought will prove that this operation is simply the much faster ordinary matrix multiplication, C = A*B;.

A note on terminology: ``broadcasting'' is the term popularized by the Numpy numerical environment in the Python programming language.  In other programming languages and environments, broadcasting may also be known as binary singleton expansion (BSX, in MATLAB, and the origin of the name of the bsxfun function), recycling (R programming language), single-instruction multiple data (SIMD), or replication.

Broadcasting and Legacy Code
----------------------------

The new broadcasting semantics almost never affect code that worked in previous versions of Octave.  Consequently, all code inherited from MATLAB that worked in previous versions of Octave should still work without change in Octave.  The only exception is code such as

     try
       c = a.*b;
     catch
       c = a.*a;
     end_try_catch

that may have relied on matrices of different size producing an error. Due to how broadcasting changes semantics with older versions of Octave, by default Octave warns if a broadcasting operation is performed.  To disable this warning, refer to its ID (see warning_ids):

     warning ("off", "Octave:broadcast");

If you want to recover the old behavior and produce an error, turn this warning into an error:

     warning ("error", "Octave:broadcast");

For broadcasting on scalars that worked in previous versions of Octave, this warning will not be emitted.

Function Application
====================

As a general rule, functions should already be written with matrix arguments in mind and should consider whole matrix operations in a vectorized manner.  Sometimes, writing functions in this way appears difficult or impossible for various reasons.  For those situations, Octave provides facilities for applying a function to each element of an array, cell, or struct.

(arrayfun)

(spfun)

(cellfun)

(structfun)

Accumulation
============

Whenever it's possible to categorize according to indices the elements of an array when performing a computation, accumulation functions can be useful.

(accumarray)

(accumdim)

JIT Compiler
============

Vectorization is the preferred technique for eliminating loops and speeding up code.  Nevertheless, it is not always possible to replace every loop.  In such situations it may be worth trying Octave's experimental Just-In-Time (JIT) compiler.

A JIT compiler works by analyzing the body of a loop, translating the Octave statements into another language, compiling the new code segment into an executable, and then running the executable and collecting any results.  The process is not simple and there is a significant amount of work to perform for each step.  It can still make sense, however, if the number of loop iterations is large.  Because Octave is an interpreted language every time through a loop Octave must parse the statements in the loop body before executing them. With a JIT compiler this is done just once when the body is translated to another language.

The JIT compiler is a very new feature in Octave and not all valid Octave statements can currently be accelerated.  However, if no other technique is available it may be worth benchmarking the code with JIT enabled.  The function jit_enable is used to turn compilation on or off.  The function jit_startcnt sets the threshold for acceleration.  Loops with iteration counts above jit_startcnt will be accelerated.  The function debug_jit is not likely to be of use to anyone not working directly on the implementation of the JIT compiler.

(jit_enable)

(jit_startcnt)

(debug_jit)

Miscellaneous Techniques
========================

Here are some other ways of improving the execution speed of Octave programs.

  * Avoid computing costly intermediate results multiple times. Octave currently does not eliminate common subexpressions.  Also, certain internal computation results are cached for variables.  For instance, if a matrix variable is used multiple times as an index, checking the indices (and internal conversion to integers) is only done once.

  * Be aware of lazy copies (copy-on-write). When a copy of an object is created, the data is not immediately copied, but rather shared.  The actual copying is postponed until the copied data needs to be modified.  For example:

          a = zeros (1000); # create a 1000x1000 matrix
          b = a; # no copying done here
          b(1) = 1; # copying done here

     Lazy copying applies to whole Octave objects such as matrices, cells, struct, and also individual cell or struct elements (not array elements).

     Additionally, index expressions also use lazy copying when Octave can determine that the indexed portion is contiguous in memory.  For example:

          a = zeros (1000); # create a 1000x1000 matrix
          b = a(:,10:100);  # no copying done here
          b = a(10:100,:);  # copying done here

     This applies to arrays (matrices), cell arrays, and structs indexed using ().  Index expressions generating comma-separated lists can also benefit from shallow copying in some cases.  In particular, when A is a struct array, expressions like {a.x}, {a(:,2).x} will use lazy copying, so that data can be shared between a struct array and a cell array.

     Most indexing expressions do not live longer than their parent objects.  In rare cases, however, a lazily copied slice outlasts its parent, in which case it becomes orphaned, still occupying unnecessarily more memory than needed.  To provide a remedy working in most real cases, Octave checks for orphaned lazy slices at certain situations, when a value is stored into a "permanent" location, such as a named variable or cell or struct element, and possibly economizes them.  For example:

          a = zeros (1000); # create a 1000x1000 matrix
          b = a(:,10:100);  # lazy slice
          a = []; # the original "a" array is still allocated
          c{1} = b; # b is reallocated at this point

  * Avoid deep recursion. Function calls to m-file functions carry a relatively significant overhead, so rewriting a recursion as a loop often helps.  Also, note that the maximum level of recursion is limited.

  * Avoid resizing matrices unnecessarily. When building a single result matrix from a series of calculations, set the size of the result matrix first, then insert values into it.  Write

          result = zeros (big_n, big_m)
          for i = over:and_over
            ridx = ...
            cidx = ...
            result(ridx, cidx) = new_value ();
          endfor

     instead of

          result = [];
          for i = ever:and_ever
            result = [ result, new_value() ];
          endfor

     Sometimes the number of items can not be computed in advance, and stack-like operations are needed.  When elements are being repeatedly inserted or removed from the end of an array, Octave detects it as stack usage and attempts to use a smarter memory management strategy by pre-allocating the array in bigger chunks.  This strategy is also applied to cell and struct arrays.

          a = [];
          while (condition)
            ...
            a(end+1) = value; # "push" operation
            ...
            a(end) = []; # "pop" operation
            ...
          endwhile

  * Avoid calling eval or feval excessively. Parsing input or looking up the name of a function in the symbol table are relatively expensive operations.

     If you are using eval merely as an exception handling mechanism, and not because you need to execute some arbitrary text, use the try statement instead.  See The try Statement.

  * Use ignore_function_time_stamp when appropriate. If you are calling lots of functions, and none of them will need to change during your run, set the variable ignore_function_time_stamp to "all".  This will stop Octave from checking the time stamp of a function file to see if it has been updated while the program is being run.

Examples
========

The following are examples of vectorization questions asked by actual users of Octave and their solutions.

  * For a vector A, the following loop

          n = length (A);
          B = zeros (n, 2);
          for i = 1:length (A)
            ## this will be two columns, the first is the difference and
            ## the second the mean of the two elements used for the diff.
            B(i,:) = [A(i+1)-A(i), (A(i+1) + A(i))/2)];
          endfor

     can be turned into the following one-liner:

          B = [diff(A)(:), 0.5*(A(1:end-1)+A(2:end))(:)]

     Note the usage of colon indexing to flatten an intermediate result into a column vector.  This is a common vectorization trick.