Vectorization

A computation is vectorized by taking advantage of vector operations. A variety of programming situations can be vectorized, and often improving speed to 10 times faster or even better. Vectorization is one of the most general and effective techniques for writing fast M-code.

5.1 Vectorized Computations

Many standard Matlab functions are “vectorized”: they can operate on an array as if the function had been applied individually to every element.

>> sqrt([1,4;9,16])

>> abs([0,1,2,−5,−6,−7])

ans =

1 2

3 4

ans =

0 1 2 5 6 7

Consider the following function:

function d = minDistance(x,y,z)

% Find the min distance between a set of points and the origin

nPoints = length(x);

d = zeros(nPoints,1); % Preallocate

for k = 1:nPoints % Compute distance for every point d(k) = sqrt(x(k)ˆ2 + y(k)ˆ2 + z(k)ˆ2);

end

d = min(d); % Get the minimum distance

For every point, its distance from the origin is computed and stored in d. For speed, array d is preallocated (see Section 3). The minimum distance is then found with min. To vectorize the distance

computation, replace the for loop with vector operations:

function d = minDistance(x,y,z)

% Find the min distance between a set of points and the origin

d = sqrt(x.ˆ2 + y.ˆ2 + z.ˆ2); % Compute distance for every point d = min(d); % Get the minimum distance

The modified code performs the distance computation with vector operations. The x, y and z arrays are first squared using the per-element power operator, .^ (the per-element operators for multiplication and division are .* and ./). The squared components are added with vector addition. Finally, the square root of the vector sum is computed per element, yielding an array of distances. (A further improvement: it is equivalent to compute d = sqrt(min(x.^2 + y.^2 + z.^2))).

9

The first version of the minDistance program takes 0.73 seconds on 50000 points. The vectorized version takes less than 0.04 seconds, more than 18 times faster.

Some useful functions for vectorizing computations:

min, max, repmat, meshgrid, sum, cumsum, diff, prod, cumprod, accumarray

5.2 Vectorized Logic

The previous section shows how to vectorize pure computation. But bottleneck code often involves conditional logic. Like computations, Matlab’s logic operators are vectorized:

>> [1,5,3] < [2,2,4]

ans =

1 0 1

Two arrays are compared per-element. Logic operations return “logical” arrays with binary values. How is this useful? Matlab has a few powerful functions for operating on logical arrays:

>> find([1,5,3] < [2,2,4]) % Find indices of nonzero elements

ans =

1 3

>> any([1,5,3] < [2,2,4]) % True if any element of a vector is nonzero

% (or per−column for a matrix)

ans =

1

>> all([1,5,3] < [2,2,4]) % True if all elements of a vector are nonzero

% (or per−column for a matrix)

ans =

0

Vectorized logic also works on arrays.

>> find(eye(3) == 1) % Find which elements == 1 in the 3x3 identity matrix ans =

1 % (element locations given as indices)

5

9

By default, find returns element locations as indices (see page 16 for indices vs. subscripts).

10

Example 1: Vector Normalization

To normalize a single vector v to unit length, one can use v = v/norm(v). However, to use norm for set of vectors v(:,1), v(:,2), . . . requires computing v(:,k)/norm(v(:,k)) in a loop. Alternatively,

vMag = sqrt(sum(v.ˆ2));

v = v./vMag(ones(1,size(v,1)),:);

Example 2: Removing elements

The situation often arises where array elements must be removed on some per-element condition. For example, this code removes all NaN and infinite elements from an array x:

i = find(isnan(x) | isinf(x)); % Find bad elements in x x(i) = []; % and delete them

Alternatively,

i = find(˜isnan(x) & ˜isinf(x)); % Find elements that are not NaN and not infinite x = x(i); % Keep those elements

Both of these solutions can be further streamlined by using logical indexing:

x(isnan(x) | isinf(x)) = []; % Delete bad elements

or

x = x(˜isnan(x) & ˜isinf(x)); % Keep good elements

Example 3: Piecewise functions

The sinc function has a piecewise definition,

sinc(x) =

( sin(x)/x, x = 0

1, x = 0

This code uses find with vectorized computation to handle the two cases separately:

function y = sinc(x)

% Computes the sinc function per−element for a set of x values.

y = ones(size(x)); % Set y to all ones, sinc(0) = 1 i = find(x ˜= 0); % Find nonzero x values

y(i) = sin(x(i)) ./ x(i); % Compute sinc where x ˜= 0

A concise alternative is y = (sin(x) + (x == 0))./(x + (x == 0)).

11

Example 4: Drawing images with meshgrid

The meshgrid function takes two input vectors and converts them to matrices by replicating the first over the rows and the second over the columns.

>> [x,y] = meshgrid(1:5,1:3)

The matrices above work like a map for a width 5, height 3 image. For each pixel, the x-location can be read from x and the y-location from y. This may seem like a gratuitous use of memory as x and y simply record the column and row positions, but this is useful. For example, to draw an ellipse,

% Create x and y for a width 150, height 100 image

[x,y] = meshgrid(1:150,1:100);

% Ellipse with origin (60,50) of size 15 x 40

Img = sqrt(((x−60).ˆ2 / 15ˆ2) + ((y−50).ˆ2 / 40ˆ2)) > 1;

% Plot the image

imagesc(Img); colormap(copper);

axis image, axis off

Drawing lines is just a change in the formula.

[x,y] = meshgrid(1:150,1:100);

% The line y = x*0.8 + 20

Img = (abs((x*0.8 + 20) − y) > 1);

imagesc(Img); colormap(copper);

axis image, axis off

Polar functions can be drawn by first converting x and y variables with the cart2pol function.

[x,y] = meshgrid(1:150,1:100);

[th,r] = cart2pol(x − 75,y − 50); % Convert to polar

% Spiral centered at (75,50) Img = sin(r/3 + th);

imagesc(Img); colormap(hot);

axis image, axis off

12

Example 5: Polynomial interpolation

Given n points x1 , . . . xn and corresponding function values y1 , . . . yn , the coefficients c0 , . . . , cn−1 of the interpolating polynomial

can be found by solving

P (x) = cn−1 x

n−1

+ · · · + c1 x + c0

 x1 n−1 x1 n−2 · · · x1 2 x1 1   cn 1   1 

 n−1 x n−2 · · · x 2 x

1   c

  y 

 x2 2

2 2  

n−2 

 2 

 . .  

.  =  . 

 . .  

  

  . 

  

xn n−1 xn n−2 · · · xn 2 xn 1 c0 yn

function c = polyint(x,y)

% Given a set of points and function values x and y,

% computes the interpolating polynomial.

x = x(:); % Make sure x and y are both column vectors y = y(:);

n = length(x); % n = Number of points

% Construct the matrix on the left−hand−side

xMatrix = repmat(x, 1, n); % Make an n by n matrix with x on every column powMatrix = repmat(n−1:−1:0, n, 1); % Make another n by n matrix of exponents

A = xMatrix .ˆ powMatrix; % Compute the powers

c = A\y; % Solve matrix equation for coefficients

The strategy to construct the left-hand side matrix is to first make two n×n matrices of bases and exponents and then put them together using the per element power operator, .^ . The repmat function (“replicate matrix”) is used to make the base matrix xMatrix and the exponent matrix powMatrix.

 x(1) x(1) · · · x(1) 

 x(2) x(2) · · · x(2) 

 n − 1 n − 2 · · · 0 

 n − 1 n − 2 · · · 0 

xMatrix = 

 .

 powMatrix =  

.   . . 

x(n) x(n) · · · x(n)

n − 1 n − 2 · · · 0