Product Support

1109 - Code Vectorization Guide

This technical note provides an introduction to vectorization techniques. In order to understand some of the possible techniques, an introduction to MATLAB referencing is provided. Then several vectorization examples are discussed.

This technical note examines how to identify situations where vectorized techniques would yield a quicker or cleaner algorithm. Vectorization is often a smooth process; however, in many application-specific cases, it can be difficult to construct a vectorized routine. Understanding the tools and examples presented in this technical note is a good introduction to the topic, but it is by no means an exhaustive resource of all the vectorization techniques available.

Tools and Techniques

Section 1: MATLAB Indexing or Referencing

MATLAB enables you to select subsets of an array or matrix. There are three basic types of indexing available in MATLAB: subscripted, linear, and logical. These methods are explained in the MATLAB Digest article, Matrix Indexing in MATLAB. While this article goes into detail on MATLAB indexing, the remainder of this section is a good introduction. If you are not familiar with MATLAB indexing, you should read the MATLAB Digest article, and/or the documentation on colon, paren (round(), square[], and curly{} braces) and end.

Subscripted Indexing

In subscripted indexing, the values of the subscripts are the indices of the matrix where the matrix's elements are desired. That is, elements = matrix(subscripts). Thus, if A = 6:10, then A([3,5]) denotes the third and fifth elements of matrix A:

A = 6:10; 
     
A([3,5])

ans=
8 10

For multidimensional arrays, multiple index parameters are used for subscripted indexing:

A = [11 14 17; ... 
     12 15 18; ... 
     13 16 19]; 
     
A(2:3,2)

ans=
15
16

Linear Indexing

In linear, or absolute indexing, every element of a matrix can be referenced by its subscripted indices (as described above), or by its single, columnwise, linear index:

A = [11 14 17; ... 
     12 15 18; ... 
     13 16 19]; 
     
A(6)

ans=
16

A([3,1,8])

ans=
13 11 18

See the functions SUB2IND and IND2SUB for more information.

Also, notice how the returned elements of the indexed matrix preserve the shape that was specified by the index variable. In the previous example, the index variable was a vector of size 1-by-3, and so the result was also of size 1-by-3. If an index variable of size 3-by-1 was used, a preserved shape would have been seen in the results:

A([3;1;8])

ans=
13
11
18

Logical Indexing

With logical, or Boolean, indexing, the index parameter is a logical matrix that is the same size as A and contains only 0's and 1's.

The elements of A that are selected have a '1' in the corresponding position of the logical indexing matrix. For example, if A = 6:10, then A(logical([0 0 1 0 1])) denotes the third and fifth elements of A:

A = 6:10; 
    
A(logical([0 0 1 0 1]))

ans=
8 10

For more information on matrix indexing and matrix manipulation, see the resources listed under Section 10 of this technical note.

Section 2: Array Operations vs. Matrix Operations

y(i) = fcn(x1(i), x2(i), ...)

The simplest type of vector operations, array operations, can be thought of as bulk processing. In this approach, the same operation is performed for each corresponding element in a data set, which may include more than one matrix. The operation is performed on an element-by-element basis.

Matrix operations, or linear algebra operations, operate according to the rules of linear algebra. This type of operation is sometimes useful for vectorization as illustrated in Section 6 and Section 8 of this technical note.

For example, you may have data measurements from an experiment that measures the volume of a cone, (V = 1/12 * pi * D^2 * H), by recording the diameter (D) of the end, and the height (H). If you had run the experiment once, you would just have one value for each of the two variables; however, D and H are scalars. Here is the calculation you will need to run in MATLAB:

V = 1/12*pi*(D^2)*H;

Now, suppose that you actually run the experiment 100 times. Now D and H are vectors of length 100, and you want to calculate the corresponding vector V, which represents the volume for each run.

In most programming languages, you would set up a loop, the equivalent of the following MATLAB code:

for n = 1:100
    V(n) = 1/12*pi*(D(n)^2)*H(n);
end

With MATLAB, you can ignore the fact that you ran 100 experiments. You can perform the calculation element-by-element over each vector, with (almost) the same syntax as the first equation, above.

% Provide data so that the code can be executed 
% Use of negative values will be apparent later 
D = [-0.2 1.0 1.5 3.0 -1.0 4.2 3.1]; 
H = [ 2.1 2.4 1.8 2.6 2.6 2.2 1.8]; 


% Perform the vectorized calculation 
V = 1/12*pi*(D.^2).*H;

The only difference is the use of the .* and ./ operators. These differentiate array operators (element-by-element operators) from the matrix operators (linear algebra operators), * and /.

For more information, refer to the Matrices and Linear Algebra section of the MATLAB Mathematics documentation.

Section 3: Boolean Array Operations

y = bool(x1, x2, ...)

A logical extension of the bulk processing technique is to vectorize comparisons and decision making. MATLAB comparison operators, like the array operators above, accept vector inputs and produce vector outputs.

Suppose that after running the above experiment, you find that negative numbers have been measured for the diameter. As these values are clearly erroneous, you decide that those runs for the experiment are invalid.

You can find out which runs are valid using the >= operator on the vector D:

D = [-0.2 1.0 1.5 3.0 -1.0 4.2 3.14]; 
D >= 0

ans=
0 1 1 1 0 1 1

Now you can exploit the logical indexing power of MATLAB to remove the erroneous values:

Vgood = V(D>=0);

This selects the subset of V for which the corresponding elements of D are nonnegative.

Suppose that if all values for masses are negative, you want to display a warning message and exit the routine. You need a way to condense the vector of Boolean values into a single value. There are two vectorized Boolean operators, ANY and ALL, which perform Boolean AND and OR functions over a vector. Thus you can perform the following test:

if all(D < 0) 
   warning('All values of diameter are negative.'); 
   return; 
end

You can also compare two vectors of the same size using the Boolean operators, resulting in expressions such as:

(V >= 0) & (D > H)

ans=
0 0 0 1 0 1 1

The result is a vector of the same size as the inputs.

Additionally, since MATLAB uses IEEE arithmetic, there are special values to denote overflow, underflow, and undefined operations: INF , -INF , and NaN , respectively.INF and -INF are supported by relational operators (i.e., Inf==Inf returns true, and Inf<1 returns false). However, by the IEEE standard, NaN is never equal to anything, even itself (NaN==NaN returns false). Therefore, there are special Boolean operators, ISINF and ISNAN, to perform logical tests for these values. For example, in some applications, it's useful to exclude NaNs in order to make valid computations:

x = [2 -1 0 3 NaN 2 NaN 11 4 Inf]; 
    
xvalid = x(~isnan(x))

xvalid =
2 -1 0 3 2 11 4 Inf

Also, if you try to compare two matrices of different sizes, an error will occur. You need to check sizes explicitly if the operands might not be the same size. See the documentation on SIZE, ISEQUAL, and ISEQUALWITHEQUALNANS(R13 and later) for more information.

Section 4: Constructing Matrices from Vectors

y(:,i) = fcn(x(:))

Creating simple matrices, such as constant matrices, is easy in MATLAB. The following code creates a size 5-by-5 matrix of all 10s:

A = ones(5,5)*10

The multiplication here is not necessary; you can achieve the same result using the indexing technique described in the next example.

Another common application involves selecting specified elements of a vector to create a new matrix. This is often a simple application of the indexing power of MATLAB. The next example creates a matrix corresponding to elements [2, 3; 5, 6; 8, 9; 11, 12; ...] of a vector x:

% Create a vector from 1 to 51 
x = 1:51; 

% Reshape it such that it wraps through 3 rows 
x = reshape(x, 3, length(x)/3); 

% Select the 2nd and 3rd rows, and transpose it
A = x(2:3,:)';

In this example, indexing is used to exploit the pattern of the desired matrix. There are several tools available to exploit different types of patterns and symmetry. Relevant functions are summarized in Section 5 of this technical note.

There is a well-known technique for duplicating a vector of size M-by-1 n times, to create a matrix of size M-by-N. In this method, known as Tony's Trick, the first column of the vector is indexed (or referenced) n times:

v = (1:5)' n = 3; 
    
M = v(:,ones(n,1))

M =
1 1 1
2 2 2
3 3 3
4 4 4
5 5 5

Now it is apparent how the first matrix, size 5-by-5 constant matrix of all 10's, can be created without the array multiplication operation:

A = 10; 
    
A = A(ones(5,5))

A =
10 10 10 10 10
10 10 10 10 10
10 10 10 10 10
10 10 10 10 10
10 10 10 10 10

The same technique can be applied to row vectors by switching the subscripts. It can also duplicate specific rows or columns of a matrix. For example,

B = magic(3)

B =
8 1 6
3 5 7
4 9 2

N = 2; B(:,N(ones(1,3)))

ans=
1 1 1
5 5 5
9 9 9

While it is good to understand how this technique works, usually you do not need to use it explicitly. Instead, you can use the functions REPMAT and MESHGRID, which implement this technique. Refer to the documentation and Section 6 for more information and examples.

Section 5: Utility Functions

The following table lists functions that are useful for creating patterned matrices and symmetric matrices. Some have already been used in this technical note. Clicking on a function will take you to the documentation for more information.

CUMPROD	Provide a cumulative product of elements
CUMSUM	Provide a cumulative sum of elements
DIAG	Create diagonal matrices and diagonals of a matrix
EYE	Produce an identity matrix
FILTER	Create one-dimensional digital filter (good for vectors where each element depends on the previous element)
GALLERY	Provides a gallery of different types of matrices
HANKEL	Create a Hankel matrix
MESHGRID	Create `X` and `Y` arrays for 3-D plots
NDGRID	Generate of arrays for N-D functions and interpolation. (Note: uses different coordinate system than `MESHGRID`)
ONES	Create an array of ones
PASCAL	Create a Pascal matrix
REPMAT	Replicate and tile an array
RESHAPE	Change size of array
SPDIAGS	Create a sparse matrix formed from diagonals
TOEPLITZ	Create a Toeplitz matrix
ZEROs	Create an array of zeros

Extended Examples

Section 6: Matrix Functions of Two Vectors

y(i,j) = fcn(x1(i), x2(j))

Suppose you want to evaluate a function F of two variables:

F(x,y) = x*exp(-x² - y²)

You need to evaluate the function at every point in vector x, and for each point in x, at every point in vector y. In other words, you need to define a grid of values for F, given vectors x and y.

You can duplicate x and y to create an output vector of the desired size using MESHGRID. This allows you to use the techniques from Section 2 to compute the function.

x = (-2:.2:2); 
y = (-1.5:.2:1.5)';
[X,Y] = meshgrid(x, y); 
F = X .* exp(-X.^2 - Y.^2);

In some cases, you can use the matrix multiplication operator in order to avoid creating the intermediate matrices. For example, if

F(x,y) = x*y

where x and y are vectors, then you can simply take the outer product of x and y:

x = (-2:2); 
y = (-1.5:.5:1.5); 
     
x'*y

      3.0000 2.0000 1.0000      0 -1.0000 -2.0000 -3.0000
      1.5000 1.0000 0.5000      0 -0.5000 -1.0000 -1.5000
           0       0       0      0       0       0       0
     -1.5000 -1.0000 -0.5000      0 0.5000 1.0000 1.5000
     -3.0000 -2.0000 -1.0000      0 1.0000 2.0000 3.0000

When creating matrices from two vectors, there are also cases where sparse matrices make use of more efficient use of storage space, and also make use of very efficient algorithms. An example is discussed further in Section 8.

Section 7: Ordering, Setting, and Counting Operations

In the examples discussed so far, any calculations done on one element of a vector have been independent of other elements in the same vector. However, in many applications, the calculation that you are trying to do depends heavily on these other values. For example, suppose you are working with a vector x which represents a set. You do not know without looking at the rest of the vector whether a particular element is redundant and should be removed. It is not obvious at first how to remove these redundant values without resorting to loops.

This area of vectorization requires a fair amount of ingenuity. There are many functions available that you can use to vectorize code:

DIFF	Acts as difference operator: diff(X), for a vector `X`, is: `[X(2) - X(1), X(3) - X(2), ... X(n) - X(n-1)]`
FIND	Finds indices of the nonzero, non-NaN elements
INTERSECT	Finds the set intersection
MAX	Find largest component
MIN	Find smallest component
SETDIFF	Finds the set difference
SETXOR	Finds the set exclusive `OR`
SORT	Sort in ascending order
UNION	Finds the set union
UNIQUE	Find unique elements of a set

Now, proceed with this example, eliminating redundant elements of a vector. Note that once a vector is sorted, any redundant elements are adjacent. In addition, any equal adjacent elements in a vector create a zero entry in the DIFF of that vector. This suggests the following implementation for the operation. You are attempting to select the elements of the sorted vector that correspond to nonzero differences.

% First try. NOT QUITE RIGHT!! 
x = sort(x(:)); 
difference = diff(x); 
y = x(difference~=0);

This is almost correct, but you have forgotten to take into account the fact that the DIFF function returns a vector that has one fewer element than the input vector. In your first algorithm, the last unique element is not accounted for. You can fix this by adding one element to the vector x before taking the difference. You need to make sure that the added element is always different than the previous element. One way to do this is to add a NaN.

% Final version. 
x = sort(x(:)); 
difference = diff([x;NaN]); 
y = x(difference~=0);

The (:) operation was used to ensure that x is a vector. The ~=0 operation was used rather than the FIND function because the FIND function does not return indices for NaN elements, and the last element of difference, as defined, is a NaN. Alternatively, this example can be accomplished with the code

y=unique(x);

but the point of the above exercise is to exploit vectorization and demonstrate writing your own code for efficiency. Additionally, the UNIQUE function provides more functionality than necessary, and this can slow down the execution of the code.

Suppose that you do not just want to return the set x. You want to know how many instances of each element in the set occurred in the original matrix. Once you have the sorted x, you can use the FIND function to determine the indices where the distribution changes. The difference between subsequent indices will indicate the number of occurrences for that element. This is the "diff of find of diff" trick. Building on the above example:

% Find the redundancy in a vector x 
x = sort(x(:)); 
difference = diff([x;max(x)+1]); 
count = diff(find([1;difference]));
y = x(find(difference)); 
plot(y,count)

This plots the number of occurrences of each element of x. Note that we avoided using NaN here, because FIND does not return indices for NaN elements. However, the number of occurrences of NaNs and Infs are easily computed as special cases:

count_nans = sum(isnan(x(:))); 
count_infs = sum(isinf(x(:)));

Another trick for vectorizing summing and counting operations is to exploit the way sparse matrices are created. This is discussed in greater detail in an example in Section 9.

Section 8: Sparse Matrix Structures

You can use sparse matrices to increase efficiency in some cases. Often vectorization is easier if you construct a large intermediate matrix. In some cases, you can take advantage of a sparse matrix structure to vectorize code without requiring large amounts of storage space for this intermediate matrix.

Suppose that in the last example, you knew beforehand that the domain of the set y is a subset of the integers, {k+1, k+2, ..., k+n}; that is,

y = 1:n + k

These might represent indices in a colormap. You can then count the occurrences of each element of the set. This is an alternative method to the "diff of find of diff" trick from the last Section.

Construct a large m-by-n matrix A, where m is the number of elements in the original vector x, and n is the number of elements in the set y.

if x(i) = y(j) 
   A(i,j) = 1
else A(i,j) = 0

Looking back at Sections 3 and 4, you might suspect that you need to construct matrices from x and y. This would work, but it would require a lot of storage space. You can do better by exploiting the fact that most of the matrix A consists of 0's, with only one 1 value per element of x.

Here is how to construct the matrix (note that y(j) = k+j, from the above formula):

x = sort(x(:)); 
A = sparse(1:length(x), x+k, 1, length(x), n);

Now you can perform a sum over the columns of A to get the number of occurrences.

count = sum(A);

In this case you do not have to form the sorted vector y explicitly, since you know beforehand that y = 1:n + k.

The key here is to use the data (i.e., x) to control the structure of the matrix A. Since x takes on integer values in a known range, you are able to construct the matrix more efficiently.

Suppose you want to multiply each column of a very large matrix by the same vector. There is a way to do this using sparse matrices that saves space and can also be faster than matrix construction using the tools outlined in Section 5. Here's how it works:

F = rand(1024,1024); 
x = rand(1024,1);   % Point-wise multiplication of each row of F .
Y = F * diag(sparse(x)); % Point-wise multiplication of each column of F. 
Y = diag(sparse(x)) * F;

The matrix multiplication operator handles the bulk of the work, while sparse matrices prevent the temporary variable from being too large.

You can also use sparse matrices for counting and binning operations. This involves manipulating the way sparse matrices are created. This is discussed in the next section of examples.

Section 9: Additional Examples

The following examples use several techniques discussed in this technical note, as well as a few other relevant techniques. Try using tic and toc (or t=cputime and cputime-t) around the examples to see how they speed up the code.

Vectorizing a double FOR loop that creates a matrix by computation:

Tools:

array multiplication

Before:

A = magic(100); 
B = pascal(100); 
for j = 1:100 
   for k = 1:100; 
     X(j,k) = sqrt(A(j,k)) * (B(j,k) - 1); 
   end 
end

After:

A = magic(100); 
B = pascal(100); 
X = sqrt(A).*(B-1);

Vectorizing a loop that creates a vector whose elements depend on the previous element:

Tools:

FILTER, CUMSUM, CUMPROD

Before:

A = 1; 
L = 1000; 
for i = 1:L 
  A(i+1) = 2*A(i)+1; 
end

After:

L = 1000; 
A = filter([1],[1 -2],ones(1,L+1));

In addition, if your vector construction only uses addition or multiplication, you can use the CUMSUM or CUMPROD function.

Before:

n=10000; 
V_B=100*ones(1,n); 
V_B2=100*ones(1,n); 
ScaleFactor=rand(1,n-1); 
for i = 2:n 
  V_B(i) = V_B(i-1)*(1+ScaleFactor(i-1));
end 
for i=2:n 
  V_B2(i) = V_B2(i-1)+3; 
end

After:

n=10000; 
V_A=100*ones(1,n); 
V_A2 = 100*ones(1,n); 
ScaleFactor=rand(1,n-1); 
V_A=cumprod([100 1+ScaleFactor]);
V_A2=cumsum([100 3*ones(1,n-1)]);

Note that if you run the two examples one after another, V_B and V_A will differ. This is because you recreate it before computing V_A. If you use the same ScaleFactor matrix for both examples, V_B and V_A will agree.

Vectorizing code that finds the cumulative sum of a vector at every fifth element:

Tools:

CUMSUM, vector indexing

Before:

% Use an arbitrary vector, x 
x = 1:10000; 
y = []; 
for n = 5:5:length(x) 
  y = [y sum(x(1:n))]; 
end

After, using preallocation:

x = 1:10000; 
ylength = (length(x) - mod(length(x),5))/5; 
% Avoid using ZEROS function during preallocation
y(1:ylength) = 0; 
for n = 5:5:length(x) 
  y(n/5) = sum(x(1:n)); 
end

After, using vectorization (preallocation is no longer needed):

x = 1:10000; 
cums = cumsum(x); 
y = cums(5:5:length(x));

Vectorizing code that repeats a vector value when the following value is zero:

Tools:

FIND, CUMSUM, DIFF

Ideally, you want to replace the zeros within a vector with values and zeros by repeating the previous value. For example, convert this:

a=2; 
b=3; 
c=5; 
d=15; 
e=11;
x = [a 0 0 0 b 0 0 c 0 0 0 0 d 0 e 0 0 0 0 0];

into this:

x = [a a a a b b b c c c c c d d e e e e e e];

Solution:

valind = find(x); 
x(valind(2:end)) = diff(x(valind)); 
x = cumsum(x);

Vectorizing code that accumulates a sum at designated indices:

Tools:

SPARSE

Before:

% The values are summed at designated indices 
values = [20 15 45 50 75 10 15 15 35 40 10]; 
% The indices associated with the values are summed cumulatively 
indices = [2 4 4 1 3 4 2 1 3 3 1]; 
totals = zeros(max(indices),1);
for n = 1:length(indices) 
  totals(indices(n)) = totals(indices(n)) + values(n); 
end

After:

% The values are summed at designated indices 
values = [20 15 45 50 75 10 15 15 35 40 10]; 
% The indices associated with the values are summed cumulatively 
indices = [2 4 4 1 3 4 2 1 3 3 1]; 
totals = full(sparse(indices,1,values));

This method exploits the way sparse matrices are created. When using the SPARSE function to create a sparse matrix, any values that are assigned to the same index are summed rather than replacing the existing value. This is called "vector accumulation," and is the way that MATLAB handles sparse matrices.

More Resources

Section 10: Matrix Indexing and Manipulation

The MATLAB Digest article Matrix Indexing in MATLAB provides much more detail on indexing than Section 1 of this technical note.

The Getting Started documentation link will guide you through matrix manipulation in MATLAB. It includes matrix creation, indexing, manipulation, array operation, matrix operations, as well as other topics. (*Click on the Manipulating Matrices link)

Peter Acklam has built a Web page devoted to MATLAB array manipulation tips and tricks.

Section 11: Matrix Memory Preallocation

For more information on preallocation to speed up computations, refer to the following online resolution:
How do I preallocate memory when using MATLAB?
The following technical note is a comprehensive guide for memory management in MATLAB:
The Technical Support Guide to Memory Management

Section 12: Maximizing MATLAB Performance

While this technical note generally discusses how to vectorize code, often times the motivation is to speed up the performance. Therefore, this section provides some additional resources for maximizing the performance of your code:

General tips on speeding up code:
How do I increase the speed or performance of MATLAB?

Conclusion

There are conventional methods for vectorizing code and reducing run-time as introduced throughout this technical note, but there are also techniques (such as the last example in Section 9) that exploit the functionality of the SPARSE function. The above discussion is by no means an exhaustive list of all the useful techniques available in MATLAB. It is meant to be an introduction to ideas and theories that can be applied to vectorizing code.

Products & Services	Industries	Academia	Support	User Community	Company