Department of Engineering

Matlab vectorisation tricks

Some basic tips on speeding up matlab code and on exploiting vectorisation are mentioned in the Optimisation section of our matlab page. Some of the tricks below come from newsreader.mathworks.com and the Mathworks site. They are low-level and capable of delivering order-of-magnitude improvements.

I've added the author names (where known). If you have contributions, mail them to tpl@eng.cam.ac.uk.

Indexing using vectors

• If X and V are vectors, then X(V) is [X(V(1)), X(V(2)), ..., X(V(n))]. So for example if X=[2 5 8 11 14 17 20 23 26 29] and V=[4 2 6] then X(V) is [11 5 17]. This idea is often used in association with others. For example, Matlab has a randperm(n) routine that returns a random permutation of the integers from 1 to n. You can use this to produce a random permutation of X by doing X(randperm(length(X)))
• If X and V are the same size and V only consists of true or false elements then MATLAB interprets V as a mask, and returns only the elements of X whose position corresponds to the location of a true in V. For example, if X is an array, then X>6 is an array the same size as X with trues where the corresponding element in X is >6 and falses elsewhere. This array can be used as a "mask" to select only the elements in X which are >6. Try
  X=1:10 V=X>6 X(V)
  or, more succinctly,
  X=1:10 X(X>6)

Creating and manipulating Matrixes

• Array Manipulation - look at flipud, fliplr, rot90, triu (extracts upper triangle), tril (extracts lower triangle), filter.
• Array Creation - look at hankel (try hankel(1:5), for instance), kron, toeplitz, diag, sparse, repmat (replicates and tiles matrices) and find.

Reshaping matrices can help help when vectorising. Fortunately reshaping is a cheap operation - the underlying data isn't moved. It often helps to be able to locate the original location of an element even if the matrix has changed shape. sub2ind and ind2sub are useful in this regard. The following example creates a 3D matrix where one element is 7, makes that matrix into a vector and finds where the 7 is (it's the 3rd element), then works out where that element is in the original matrix.

Extra features of common routines

Some fairly commonly used routines have extra features that are especially useful when vectorising

• ismember - Suppose you have a 2D matrix A and another a 2D matrix B with the same number of columns. To find which rows of A are in B, you can do ismember(A,B,'rows')
• all - to remove from a 2D matrix M all the rows that contain only zeroes, you can do m(all(m==0,2),:)=[] (the ",2" argument to the all command does a per-row comparison)
• unique - this routine returns an array with duplicates removed, but the returned array is sorted too. If you don't want this, try something like
  a=[3 2 2 1 7]; [B,I,J] = unique(a,'first') a(sort(I))
  Note that the final line uses a vector as an 'index'

Vectorising Routines

• bsxfun (Binary Singleton eXpansion FUNction) is quite a new routine. The example that Mathworks offer subtracts the column means from the matrix A
  A = magic(5); A = bsxfun(@minus, A, mean(A))
  Here's another example
  A = magic(5); B= A+3; % The next line is the equivalent of sqrt(A.^2 + B.^2) A = bsxfun(@hypos, A, B)

Examples

These examples are short, so by reading about the functions used and testing with small matrices, you should be able to discover why they work. More than one way is shown to solve some of these questions so that you can compare methods. When big matrices are used, you might find that some methods are hundreds of times faster than others. Remember that the most elegant-looking way may not be the fastest, and faster methods may use a lot more memory. Use flops or tic and toc to assess performance. At the end there's a worked example

• Create a 1000 by 1000 vector full of 7s
  1. X=7*ones(1000,1000);
  2. X=repmat(7,1000,1000);
  3. X(1:1000,1:1000)=7;
  Which method is fastest depends on the version of matlab. For me, currently, method 1 is fastest, with methods 2 and 3 being 50-100% slower. The moral is: test and measure, don't just hypothesise or depend on habit.
• How can you reverse a row vector?
  1. X(end:-1:1)
  2. fliplr(X)
• How can you reverse one column of a matrix
  1. To reverse column c of a matrix M, try M(:,c)=flipud(M(:,c)).
• Remove from a vector all the elements that are equal to the biggest element

     x(x==max(x)) = [] 
  

  (Jos)
• Remove from a 2-D matrix all the rows that contain at least one element less than 3

      A=magic(5)
      A(ismember(A<3,zeros(size(A,2)),'rows'),:)   
  

• How can you swap rows?
  1. M([1,5],:) = M([5,1],:) interchanges rows 1 and 5 of M.
  2. More generally, if V has m components and W has n components, then M(V,W) is the m-by-n matrix formed from the elements of M whose subscripts are the elements of V and W.
• Vectorise the following, where data is 200x400

  for i = 1:200
     for j = 1:400
       if data(i,j) > 0
             data(i,j) = 0;
       end
     end
  end
  

  data(data > 0) = 0;
• Vectorise the following threshold function

  function A = threshold(Matrix, max_value, min_value)
      [a, b] = size(Matrix);
      for i = 1:a
         for j = 1:b
                c =  Matrix(i,j);
                if ( c > max_value )
                   Matrix(i,j) = max_value;
                elseif ( c < min_value )
                   Matrix(i,j) = min_value;
                end
         end
      end
  

  By Peter J. Acklam
  A = min( max( Matrix, min_value ), max_value );
• z is a vector of 100 numbers. Produce a vector p where p(1) sums the 1st 5 elements of z, p(2) sums the next 5 and so on.

  This can be done using a loop inside a loop, but by reshaping the data to match the way matlab's sum command works, you can avoid explicit loops.

  p = sum(reshape(z,[5 20]))
  (Anh Huy Phan)
• Sum all the elements of a 2D matrix m that are greater than 23
  m=magic(5); % create a sample matrix sum(m(m(:)>23))
• Obtain the mean of each column of x. Ignore elements <=0.
  keepers = (x>0); colSums = sum(x .* keepers); counts = sum(keepers); means = colSums ./ counts;
• Vectorise

    for i = 1:100
      for j = 1:100
        r(i,j) = sqrt(i^2+j^2);
      end
    end
  

  [i,j]=meshgrid(1:100,1:100); r = sqrt(i.^2+j.^2); % The following alternative works with new versions bsxfun(@hypos, i, j)
• Vectorise the following, where elements depend on previous ones

  n=1000;
  x(1)=1;
  for j=1:n-1,
     x(j+1) = x(j) + n - j;
  end
  

  n=1000; x(1)=1; j=1:n-1; x(j+1) = n - j; cumsum(x);
• Subtract 3 from each element of x which is greater than 3
  x(x>3)=x(x>3)-3;
  (Yuri Strukov)
• From A = [3,4,3] and B = [1,2,3,4,5,6,7,8,9,10] produce a vector C where C(1) is the sum of the first A(1) elements of B, C(2) is the sum of the next A(2) elements of B, etc.
  By Ulrich Elsner
  foo=cumsum(B); C=diff([0 foo(cumsum(A))])
• Find the maximum of sin(x) * cos(y) where x is 1,2..7 and y is 1,2 .. 5.
  x=1:7 [rx,cx]=size(x); y=1:5 [ry,cy]=size(y); X=repmat(x,cy,1) Y=repmat(y',1,cx) % Doing max(max( sin(X).*cos(Y))) would find the max in one % line. To know the location of the maximum too, use XY=sin(X).*cos(Y) [r,c]=max(XY); [r2,c2]=max(max(XY)); sprintf('The biggest element in XY is %d at XY(%d,%d)',r2,c(c2), ... c2)
• Count how many times the components of a matrix are repeated.
  1. bins = min(min(X)):max(max(X)); [numTimesInMatrix, Number] = hist(X(:),bins);
     From Mathworks - "Hist can only partially vectorize this problem, so it uses a loop as well. Therefore, the algorithm runs with order O(length(bins)). Also, this algorithm reports the number of occurrences of all numbers in bin, including those that appear zero times. To get rid of these, add this":
     ind=find(~numTimesInMatrix); Number(ind)=[];numTimesInMatrix(ind)=[];
  2. x = sort(X(:)); difference = diff([x;max(x)+1]); count = diff(find([1;difference])); y = x(find(difference));
     Note that we use the (:) operation to ensure that x is a vector.
  3. count = sparse(1,X,1);
     (read about sparse if, like me, you were surprised by this)
• Retrieve those elements that are shared in matrices A and B
  1. intersect(A(:),B(:))
     (Sijmen-de.Jong)
  2. For positive integers:
     a=A(:); % make A a row b=B(:); m=max([a;b]); x1=zeros(m,1); x2=zeros(m,1); x1(a)=ones(length(a),1); x2(b)=ones(length(b),1); find(x1&x2)
  3. For equally sized A and B:
     a = A(:)'+1; b = B(:)'+1; simple = [find(sparse(1,a,1)>0);find(sparse(1,b,1)>0)]; values = find(full(sparse(1,simple',1))>1)-1
• Scale all the rows of each column by the data in col 300.
  [nr,nc] = size(A); B = A(:,300); A = A ./ B( :, ones(nc,1) );
• Write an m-file which will replace each element of a matrix with a 4x4 matrix of that element.
  1. From pete@electrosystems.com
     x=[1 2 3; 4 5 6]; y=kron(x,ones(4))
  2. A=randn(100); A(ceil((1:(size(A,1)*4))/4),ceil((1:(size(A,2)*4))/4));
• A vector x contains some 0s. Create y s.t. y[i] is 0 if x[i] is 0, otherwise y[i] is log(x[i])
  y=zeros(size(x)); y(find(x))=log(x(find(x)))
• Duplicate each element of a vector.
  1. By Ulrich Elsner
     rep=2; foo=repmat(b,rep,1); b=foo(:)' % reshape into row
• Multiply every column of matrix X of size (1000,500) by another vector V which has dimension (1000).
  1. A non-vectorised way to do this is:
     X2 = zeros(1000,500) for k = 1:500 X2(:,k) = V.*X(:,k); end
  2. X2 = diag(V)*X;
  3. rows=1000; columns = 500; Gcolvec = [1:rows]'; % Be sure that G is a column vector Gmatrix = repmat(Gcolvec,1,columns); R = ones(rows, columns); tic, GxR = Gmatrix.*R; toc
  4. By marco@chinook.physics.utoronto.ca
     Ra = R .* (G * ones(1,500)); Rb = R .* repmat(G,1,500); Rc = R .* G(:,ones(1,500));
  5. The solution proposed in the Mathworks "How Do I Vectorize My Code?" is the fastest
     X2 = diag(sparse(V))*X;
• If you have two matrices, which together give x and y coordinates into another matrix, eg.

   xc = [1 2 1; 3 2 4 ];
   yc = [1 1 1; 2 2 1; ];
  
   M = [ 6  7  8  9;
        10 11 12 13];
  

  produce the 2x3 matrix R which contains the values of M at the locations given in xc and yc. Thus:

   R = [ 6  7 6;
        12 11 9];
  

  1. R=M(sub2ind(size(M),yc,xc));
  2. a = [1 4 3 2] b = [2 3 1 4] ind = sub2ind([4 4],a,b) A = rand(4,4) A(ind)
• Create a matrix X, where each column is a shifted copy of the vector v
  v = (1:5)'; X = toeplitz(v, v([1,length(v):-1:2]))
  To shift up rather than down, use
  % By Robert Richter hankel(v,v([end,1:end-1]))
• Unique sort (i.e. duplicates removed)
  • x = sort(x(:)); unique(x)
  • x = sort(x(:)); difference = diff([x;NaN]); y = x(difference~=0);
    From Mathworks - "You may be wondering why we didn't use the find function, rather than ~=0 - especially since this function is specifically mentioned above. This is because the find function does not return indices for NaN elements, and the last element of difference as we defined it is a NaN."
• Given a 2D matrix remove all duplicate rows (where a row is a considered a duplicate of another if it has the same elements)
  a=[1 2 5 2 6 7 1 2 9 1 3 6 5 1 2 3 6 7 7 2 6] [result_sorted,INX] = unique(sort(a,2),'rows'); a(INX,:)
  (by Jos)
• Create a vector vi which is made by interleaving elements from an array v1 of length N and an array v2 of length N-1
  % By Peter J. Acklam N=10000 v1 = 1:N; v2 = 1:N-1; vi = zeros(1, 2*N-1); vi(1:2:end) = v1; vi(2:2:end) = v2; % or vi = [ v1 ; v2 0 ]; vi = vi(1:end-1);
• Write a function d = nearney(x, y) where x, y and d are column vectors of the same length and where d(i) is the smallest distance on a plane from the point [x(i) y(i)] to [x(j) y(j)] for all i ~= j; that is the distance from [x(i) y(i)] to its nearest neighbor.
  1. From moler@mathworks.com - "This gets good marks for terseness, but has lousy time and space complexity"
     function d = nearney(x,y) [X,Y] = meshgrid(x+i*y); d = min(abs(X-Y) + realmax*eye(length(x)))';
  2. From Sijmen de Jong
     z = [x, y]; d = sum(z'.^2); [m,n] = size(z); d = sqrt(min(d(ones(m,1),:)+d(ones(m,1),:)'-2*z*z'+realmax*eye(m)));
• Another 'find the nearest' question Write the code to do this

  function out = nearest(y,x1,x2);
  % function out = nearest(y,x1,x2);
  % y,x1,x2, and out are row vectors of the same length
  % out(i) = x1(i) or x2(i) depending on which one is closer to y(i)
  %
  % for example if:
  %        y = [3    5    7  9    11]
  %       x1 = [3.1  6    7  8.9  12]
  %       x2 = [3.2  5.9  0  0    10.9]
  % then out = [3.1  5.9  7  8.9  10.9]
  

  1. a = [x1; x2]; [c d] = min(abs([y; y] - a)); out = diag(a(d,:))';
  2. out = x1; i2 = find(abs(x1-y) > abs(x2-y)); out(i2) = x2(i2);
  3. More generally
     function out = nearest(y,x) % function out =nearest(y,x) % y is a row-vector and x is composed of row-vectors of the same length as y. % % out(i)=the element of x(:,i) closest to y(i) % % For example if: % y= [3 5 7 9 11] % x=[[3.1 6 7 8.9 12]; % [3.2 5.9 0 0 10.9]] % then out=[3.1 5.9 7 8.9 10.9] % [c d]=min(abs(y(ones(size(x,1),1),:)-x)); out=x(d+(0:size(x,2)-1)*size(x,1));
• Given a matrix of points, find the distance between each pair of points
  points = [ 3 4 5; 6 5 7; 6 7 1; 7 6 0; 3 2 4; 1 0 3; 5 6 2; 9 0 1] m=8; p=3; % Here's a one-liner by Peter Acklam. Results are in D D= sqrt(sum(abs(repmat(permute(points, [1 3 2]), [1 m 1]) ... - repmat(permute(points, [3 1 2]), [m 1 1]) ).^2,3)); % Here's another solution by Peter Acklam which exploits symmetry [ i j ] = find(triu(ones(m),1)); E= zeros(m,m); E(i+m*(j-1) ) = sqrt(sum(abs(points(i,:) -points(j,:)).^2,2)); E(j + m*(i-1)) = E(i+m*(j-1)); % And here's a solution by Jyri Kivinen based on the notes in % Borg and Groenen's "Modern Multidimensional Scaling: Theory and % Applications" K=points*points'; d=diag(K); one=ones(length(d),1); D=sqrt(d*one'+one*d'-2*K);
• Suppose have a 4-column matrix and a row vector [0 1 0 1]. How can I create new matrix containing the same columns data as the original data when the corresponding row vector element is 1, and a column of 0s otherwise?
  Given
  A= [1 2 3 4 6 7 1 2 1 2 4 3 1 2 3 1]; v=[0 1 0 1];
  1. [M,N] =size(A); B=A.*repmat(v,M,1)
  2. B=A*diag(v)
     (by /PB)
• Given a vector how can each sequence of more than 1 zero be reduced to a single zero? E.g. given [0 0 1 5 6 0 7 0 7 0 0 9 8 0 0] produce [0 1 5 6 0 7 0 7 0 9 8 0]
  
  a = [0 0 1 5 6 0 7 0 7 0 0 9 8 0 0] a( (a(1:end-1)==0) & (a(2:end)==0) )=[]
  (by Alan B)
• Suppose you have a vector x=[2 6 99 3 10] and a vector y=['a' 'b' 'c' 'd' 'e']. How could you re-order y so that if the elements of x were re-ordered in the same way, the elements would be in order? (in this example the re-ordered y would be "adbec")
  
  [dummy I]=sort(x); y=y(I)
• For each row of an array how do you find the longest sequence of numbers that are less than the mean for that row? (the actual question involved analysing an image of a vertical crack, measuring the width of the crack).
  
  % create a test array m=rand(10,10) % find mean of each row means=mean(m,2); for row=1:10 [Y,I] =max(diff([0 (find(~(m(row,:)<means(row)))) numel(m(1,:))+1]) -1); disp(sprintf('Row %d biggest crack (width %d) at %d ',row, Y, I)) end
  This uses code that "kinor" on a discussion board described as follows. Suppose you wanted to find out about sequences of numbers less than 1.

      diff(find(~(A<1)) )
  

  will give you the number of steps from one element >=1 to the next element >=1, the number of steps being one plus the number of elements being <1 on the way, so the number you are interested in is

      diff(find(~(A<1)) )-1
  

  You then have to take account of the first and final values, where you'll need to decide whether sequences can begin at an edge, etc.

A worked example

A bubblesort using 1 loop

Each time the program goes round the loop it will look at the (1st 2nd), (3rd 4th) , .... pairs, swapping the numbers if necessary, then it will look at the (2nd 3rd), (4th 5th) ... pairs, swapping the numbers if necessary. The resulting code might not be very fast ...

See Also

• How Do I Vectorize My Code? (from Mathworks)
• Peter Acklam's MATLAB array manipulation tips and tricks