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.

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.

 -- C =: bsxfun (F, A, B)
     Apply a binary function F element-by-element to two array arguments A and B, expanding singleton dimensions in either input argument as necessary.

     F is a function handle, inline function, or string containing the name of the function to evaluate.  The function F must be capable of accepting two column-vector arguments of equal length, or one column vector argument and a scalar.

     The dimensions of A and B must be equal or singleton.  The singleton dimensions of the arrays will be expanded to the same dimensionality as the other array. arrayfun, cellfun

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

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

Here is a real example of the power of broadcasting.  The Floyd-Warshall algorithm is used to calculate the shortest path lengths between every pair of vertices in a graph.  A naive implementation for a graph adjacency matrix of order N might look like this:

     for k = 1:n
       for i = 1:n
         for j = 1:n
           dist(i,j) = min (dist(i,j), dist(i,k) + dist(k,j));
         endfor
       endfor
     endfor

Upon vectorizing the innermost loop, it might look like this:

     for k = 1:n
       for i = 1:n
         dist(i,:) = min (dist(i,:), dist(i,k) + dist(k,:));
       endfor
     endfor

Using broadcasting in both directions, it looks like this:

     for k = 1:n
       dist = min (dist, dist(:,k) + dist(k,:));
     endfor

The relative time performance of the three techniques for a given graph with 100 vertices is 7.3 seconds for the naive code, 87 milliseconds for the singly vectorized code, and 1.3 milliseconds for the fully broadcast code. For a graph with 1000 vertices, vectorization takes 11.7 seconds while broadcasting takes only 1.15 seconds.  Therefore in general it is worth writing code with broadcasting semantics for performance.

However, 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. Because such operation is now valid Octave syntax, this will no longer produce an error.  Instead, the following code should be used:

     if (isequal (size (a), size (b)))
       c = a .* b;
     else
       c = a .* a;
     endif

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.

 -- B =: arrayfun (FCN, A)
 -- B =: arrayfun (FCN, A1, A2, ...)
 -- [B1, B2, ...] =: arrayfun (FCN, A, ...)
 -- B =: arrayfun (..., "UniformOutput", VAL)
 -- B =: arrayfun (..., "ErrorHandler", ERRFCN)

     Execute a function on each element of an array.

     This is useful for functions that do not accept array arguments.  If the function does accept array arguments it is better to call the function directly.

     The first input argument FCN can be a string, a function handle, an inline function, or an anonymous function.  The input argument A can be a logical array, a numeric array, a string array, a structure array, or a cell array.  arrayfun passes all elements of A individually to the function FCN and collects the results.  The equivalent pseudo-code is

          cls = class (FCN (A(1));
          B = zeros (size (A), cls);
          for i = 1:numel (A)
            B(i) = FCN (A(i))
          endfor

     The named function can also take more than two input arguments, with the input arguments given as third input argument A2, fourth input argument A2, ...  If given more than one array input argument then all input arguments must have the same sizes.  For example:

          arrayfun (@atan2, [1, 0], [0, 1])
               => [ 1.57080   0.00000 ]

     If the parameter VAL after a further string input argument "UniformOutput" is set true (the default), then the named function FCN must return a single element which then will be concatenated into the return value and is of type matrix.  Otherwise, if that parameter is set to false, then the outputs are concatenated in a cell array.  For example:

          arrayfun (@(x,y) x:y, "abc", "def", "UniformOutput", false)
          =>
             {
               [1,1] = abcd
               [1,2] = bcde
               [1,3] = cdef
             }

     If more than one output arguments are given then the named function must return the number of return values that also are expected, for example:

          [A, B, C] = arrayfun (@find, [10; 0], "UniformOutput", false)
          =>
          A =
          {
             [1,1] =  1
             [2,1] = [](0x0)
          }
          B =
          {
             [1,1] =  1
             [2,1] = [](0x0)
          }
          C =
          {
             [1,1] =  10
             [2,1] = [](0x0)
          }

     If the parameter ERRFCN after a further string input argument "ErrorHandler" is another string, a function handle, an inline function, or an anonymous function, then ERRFCN defines a function to call in the case that FCN generates an error.  The definition of the function must be of the form

          function [...] = errfcn (S, ...)

     where there is an additional input argument to ERRFCN relative to FCN, given by S.  This is a structure with the elements "identifier", "message", and "index" giving, respectively, the error identifier, the error message, and the index of the array elements that caused the error.  The size of the output argument of ERRFCN must have the same size as the output argument of FCN, otherwise a real error is thrown.  For example:

          function y = ferr (s, x), y = "MyString"; endfunction
          arrayfun (@str2num, [1234],
                    "UniformOutput", false, "ErrorHandler", @ferr)
          =>
             {
               [1,1] = MyString
             }

     spfun, cellfun, structfun

 -- Y =: spfun (F, S)
     Compute f (S) for the nonzero elements of S.

     The input function F is applied only to the nonzero elements of the input matrix S which is typically sparse.  The function F can be passed as a string, function handle, or inline function.

     The output Y is a sparse matrix with the same sparsity structure as the input S.  spfun preserves sparsity structure which is different than simply applying the function F to the sparse matrix S when F (0) != 0.

     Example

     Sparsity preserving spfun versus normal function application

          S = pi * speye (2,2)
          S =

          Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

            (1, 1) -> 3.1416
            (2, 2) -> 3.1416

          y = spfun (@cos, S)
          y =

          Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

            (1, 1) -> -1
            (2, 2) -> -1

          y = cos (S)
          y =

          Compressed Column Sparse (rows = 2, cols = 2, nnz = 4 [100%])

            (1, 1) -> -1
            (2, 1) -> 1
            (1, 2) -> 1
            (2, 2) -> -1

     arrayfun, cellfun, structfun

 -- A =: cellfun ("FCN", C)
 -- A =: cellfun ("size", C, K)
 -- A =: cellfun ("isclass", C, CLASS)
 -- A =: cellfun (@FCN, C)
 -- A =: cellfun (FCN, C)
 -- A =: cellfun (FCN, C1, C2, ...)
 -- [A1, A2, ...] =: cellfun (...)
 -- A =: cellfun (..., "ErrorHandler", ERRFCN)
 -- A =: cellfun (..., "UniformOutput", VAL)

     Evaluate the function named "FCN" on the elements of the cell array C.

     Elements in C are passed on to the named function individually.  The function FCN can be one of the functions

     isempty
          Return 1 for empty elements.

     islogical
          Return 1 for logical elements.

     isnumeric
          Return 1 for numeric elements.

     isreal
          Return 1 for real elements.

     length
          Return a vector of the lengths of cell elements.

     ndims
          Return the number of dimensions of each element.

     numel
     prodofsize
          Return the number of elements contained within each cell element.  The number is the product of the dimensions of the object at each cell element.

     size
          Return the size along the K-th dimension.

     isclass
          Return 1 for elements of CLASS.

     Additionally, cellfun accepts an arbitrary function FCN in the form of an inline function, function handle, or the name of a function (in a character string).  The function can take one or more arguments, with the inputs arguments given by C1, C2, etc.  For example:

          cellfun ("atan2", {1, 0}, {0, 1})
               => [ 1.57080   0.00000 ]

     The number of output arguments of cellfun matches the number of output arguments of the function and can be greater than one.  When there are multiple outputs of the function they will be collected into the output arguments of cellfun like this:

          function [a, b] = twoouts (x)
            a = x;
            b = x*x;
          endfunction
          [aa, bb] = cellfun (@twoouts, {1, 2, 3})
               =>
                  aa =
                     1 2 3
                  bb =
                     1 4 9

     Note that, per default, the output argument(s) are arrays of the same size as the input arguments.  Input arguments that are singleton (1x1) cells will be automatically expanded to the size of the other arguments.

     If the parameter "UniformOutput" is set to true (the default), then the function must return scalars which will be concatenated into the return array(s).  If "UniformOutput" is false, the outputs are concatenated into a cell array (or cell arrays).  For example:

          cellfun ("lower", {"Foo", "Bar", "FooBar"},
                   "UniformOutput", false)
          => {"foo", "bar", "foobar"}

     Given the parameter "ErrorHandler", then ERRFCN defines a function to call in case FCN generates an error.  The form of the function is

          function [...] = errfcn (S, ...)

     where there is an additional input argument to ERRFCN relative to FCN, given by S.  This is a structure with the elements "identifier", "message", and "index" giving respectively the error identifier, the error message, and the index into the input arguments of the element that caused the error.  For example:

          function y = foo (s, x), y = NaN; endfunction
          cellfun ("factorial", {-1,2}, "ErrorHandler", @foo)
          => [NaN 2]

     Use cellfun intelligently.  The cellfun function is a useful tool for avoiding loops.  It is often used with anonymous function handles; however, calling an anonymous function involves an overhead quite comparable to the overhead of an m-file function.  Passing a handle to a built-in function is faster, because the interpreter is not involved in the internal loop.  For example:

          C = {...}
          v = cellfun (@(x) det (x), C); # compute determinants
          v = cellfun (@det, C);         # 40% faster

     arrayfun, structfun, spfun

 -- A =: structfun (FCN, S)
 -- A =: structfun (..., "ErrorHandler", ERRFCN)
 -- A =: structfun (..., "UniformOutput", VAL)
 -- [A, B, ...] =: structfun (...)

     Evaluate the function named NAME on the fields of the structure S.  The fields of S are passed to the function FCN individually.

     structfun accepts an arbitrary function FCN in the form of an inline function, function handle, or the name of a function (in a character string).  In the case of a character string argument, the function must accept a single argument named X, and it must return a string value. If the function returns more than one argument, they are returned as separate output variables.

     If the parameter "UniformOutput" is set to true (the default), then the function must return a single element which will be concatenated into the return value.  If "UniformOutput" is false, the outputs are placed into a structure with the same fieldnames as the input structure.

          s.name1 = "John Smith";
          s.name2 = "Jill Jones";
          structfun (@(x) regexp (x, '(\w+)$', "matches"){1}, s,
                     "UniformOutput", false)
            => scalar structure containing the fields:
                 name1 = Smith
                 name2 = Jones

     Given the parameter "ErrorHandler", ERRFCN defines a function to call in case FCN generates an error.  The form of the function is

          function [...] = errfcn (SE, ...)

     where there is an additional input argument to ERRFCN relative to FCN, given by SE.  This is a structure with the elements "identifier", "message" and "index", giving respectively the error identifier, the error message, and the index into the input arguments of the element that caused the error.  For an example on how to use an error handler, see cellfun.

     cellfun, arrayfun, spfun

Consistent with earlier advice, seek to use Octave built-in functions whenever possible for the best performance.  This advice applies especially to the four functions above.  For example, when adding two arrays together element-by-element one could use a handle to the built-in addition function @plus or define an anonymous function @(x,y) x + y.  But, the anonymous function is 60% slower than the first method. See Operator Overloading, for a list of basic functions which might be used in place of anonymous ones.

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

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

 -- A =: accumarray (SUBS, VALS)
 -- A =: accumarray (SUBS, VALS, SZ)
 -- A =: accumarray (SUBS, VALS, SZ, FCN)
 -- A =: accumarray (SUBS, VALS, SZ, FCN, FILLVAL)
 -- A =: accumarray (SUBS, VALS, SZ, FCN, FILLVAL, ISSPARSE)

     Create an array by accumulating the elements of a vector into the positions defined by their subscripts.

     The subscripts are defined by the rows of the matrix SUBS and the values by VALS.  Each row of SUBS corresponds to one of the values in VALS.  If VALS is a scalar, it will be used for each of the row of SUBS.  If SUBS is a cell array of vectors, all vectors must be of the same length, and the subscripts in the Kth vector must correspond to the Kth dimension of the result.

     The size of the matrix will be determined by the subscripts themselves.  However, if SZ is defined it determines the matrix size.  The length of SZ must correspond to the number of columns in SUBS.  An exception is if SUBS has only one column, in which case SZ may be the dimensions of a vector and the subscripts of SUBS are taken as the indices into it.

     The default action of accumarray is to sum the elements with the same subscripts.  This behavior can be modified by defining the FCN function.  This should be a function or function handle that accepts a column vector and returns a scalar.  The result of the function should not depend on the order of the subscripts.

     The elements of the returned array that have no subscripts associated with them are set to zero.  Defining FILLVAL to some other value allows these values to be defined.  This behavior changes, however, for certain values of FCN.  If FCN is @min (respectively, @max) then the result will be filled with the minimum (respectively, maximum) integer if VALS is of integral type, logical false (respectively, logical true) if VALS is of logical type, zero if FILLVAL is zero and all values are non-positive (respectively, non-negative), and NaN otherwise.

     By default accumarray returns a full matrix.  If ISSPARSE is logically true, then a sparse matrix is returned instead.

     The following accumarray example constructs a frequency table that in the first column counts how many occurrences each number in the second column has, taken from the vector X.  Note the usage of unique  for assigning to all repeated elements of X the same index (see unique).

          X = [91, 92, 90, 92, 90, 89, 91, 89, 90, 100, 100, 100];
          [U, ~, J] = unique (X);
          [accumarray(J', 1), U']
            =>  2    89
                3    90
                2    91
                2    92
                3   100

     Another example, where the result is a multi-dimensional 3-D array and the default value (zero) appears in the output:

          accumarray ([1, 1, 1;
                       2, 1, 2;
                       2, 3, 2;
                       2, 1, 2;
                       2, 3, 2], 101:105)
          => ans(:,:,1) = [101, 0, 0; 0, 0, 0]
          => ans(:,:,2) = [0, 0, 0; 206, 0, 208]

     The sparse option can be used as an alternative to the sparse constructor (see sparse).  Thus

          sparse (I, J, SV)

     can be written with accumarray as

          accumarray ([I, J], SV', [], [], 0, true)

     For repeated indices, sparse adds the corresponding value.  To take the minimum instead, use min as an accumulator function:

          accumarray ([I, J], SV', [], @min, 0, true)

     The complexity of accumarray in general for the non-sparse case is generally O(M+N), where N is the number of subscripts and M is the maximum subscript (linearized in multi-dimensional case).  If FCN is one of @sum (default), @max, @min or @(x) {x}, an optimized code path is used. Note that for general reduction function the interpreter overhead can play a major part and it may be more efficient to do multiple accumarray calls and compute the results in a vectorized manner.

     accumdim, unique, sparse

 -- A =: accumdim (SUBS, VALS)
 -- A =: accumdim (SUBS, VALS, DIM)
 -- A =: accumdim (SUBS, VALS, DIM, N)
 -- A =: accumdim (SUBS, VALS, DIM, N, FCN)
 -- A =: accumdim (SUBS, VALS, DIM, N, FCN, FILLVAL)
     Create an array by accumulating the slices of an array into the positions defined by their subscripts along a specified dimension.

     The subscripts are defined by the index vector SUBS. The dimension is specified by DIM.  If not given, it defaults to the first non-singleton dimension.  The length of SUBS must be equal to size (VALS, DIM).

     The extent of the result matrix in the working dimension will be determined by the subscripts themselves.  However, if N is defined it determines this extent.

     The default action of accumdim is to sum the subarrays with the same subscripts.  This behavior can be modified by defining the FCN function.  This should be a function or function handle that accepts an array and a dimension, and reduces the array along this dimension.  As a special exception, the built-in min and max functions can be used directly, and accumdim accounts for the middle empty argument that is used in their calling.

     The slices of the returned array that have no subscripts associated with them are set to zero.  Defining FILLVAL to some other value allows these values to be defined.

     An example of the use of accumdim is:

          accumdim ([1, 2, 1, 2, 1], [ 7, -10,   4;
                                      -5, -12,   8;
                                     -12,   2,   8;
                                     -10,   9,  -3;
                                      -5,  -3, -13])
          => [-10,-11,-1;-15,-3,5]

     accumarray

Memoization
===========

Memoization is a technique to cache the results of slow function calls and return the cached value when the function is called with the same inputs again, instead of reevaluating it.  It is very common to replace function calls with lookup tables if the same inputs are happening over and over again in a known, predictable way.  Memoization is, at its core, an extension of this practice where the lookup table is extended even during runtime for new arguments not seen previously.  A basic theoretical background can be found on Wikipedia or any undergraduate-level computer science textbook.

Octave's memoize function provides drop-in memoization functionality for any user function or Octave function, including compiled functions.

 -- MEM_FCN_HANDLE =: memoize (FCN_HANDLE)

     Create a memoized version MEM_FCN_HANDLE of function FCN_HANDLE.

     Each call to the memoized version MEM_FCN_HANDLE checks the inputs against an internally maintained table, and if the inputs have occurred previously, then the result of the function call is returned from the table itself instead of evaluating the full function again.  This speeds up the execution of functions that are called with the same inputs multiple times.

     For example, here we take a slow user-written function named slow_fcn and memoize it to a new handle cyc.  The first executions of both versions take the same time, but the subsequent executions of the memoized version returns the previously computed value, thus reducing 2.4 seconds of runtime to only 2.4 milliseconds.  The final check verifies that the same result was returned from both versions.

          >> tic; P = slow_fcn (5040); toc
          Elapsed time is 2.41244 seconds.
          >> tic; P = slow_fcn (5040); toc
          Elapsed time is 2.41542 seconds.

          >> cyc = memoize (@slow_fcn);
          >> tic; R = cyc (5040); toc
          Elapsed time is 2.42609 seconds.
          >> tic; R = cyc (5040); toc
          Elapsed time is 0.00236511 seconds.

          >> all (P == R)
          ans = 1

     clearAllMemoizedCaches

To memoize a function z = foo(x, y), use this general pattern:

     foo2 = memoize (@(X, Y) FOO(X, Y));
     z = foo2 (x, y);

In the above example, the first line creates a memoized version foo2 of the function foo.  For simple functions with only trivial wrapping, this line can also be shortened to:

     foo2 = memoize (@foo);

The second line z = foo2 (x, y); calls that memoized version foo2 instead of the original function, allowing memoize to intercept the call and replace it with a looked-up value from a table if the inputs have occurred before, instead of evaluating the original function again.

Note that this will not accelerate the first call to the function but only subsequent calls.

Note that due to the overhead incurred by memoize to create and manage the lookup tables for each function, this technique is useful only for functions that take at least a couple of seconds to execute.  Such functions can be replaced by table lookups taking only a millisecond or so, but if the original function itself was taking only milliseconds, memoizing it will not speed it up.

Recursive functions can be memoized as well, using a pattern like:

     function z = foo (x, y)
       persistent foo2 = memoize (@foo);
       foo2.CacheSize = 1e6;

       ## Call the memoized version when recursing
       z = foo2 (x, y);
     endfunction

The CacheSize can be optionally increased in anticipation of a large number of function calls, such as from inside a recursive function.  If CacheSize is exceeded, the memoization tables are resized, causing a slowdown.  Increasing the CacheSize thus works like preallocation to speed up execution.

The function clearAllMemoizedCaches clears the memoization tables when they are no longer needed.

 -- clearAllMemoizedCaches: ()
     Clear all memoized caches.

     Memoization maintains internal tables of which functions have been called with which inputs.  This function clears those tables to free memory, or for a fresh start.

     memoize

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) - 1;
          B = zeros (n, 2);
          for i = 1:n
            ## 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.