Code Vectorization Guide：个人极力推荐！

zhouyang664

文章引自：http://www.mathworks.com/support/tech-notes/1100/1109.html

文章介绍了代码向量化的常用技术，比较经典，希望大家可以多多学习！
Tools and Techniques

• MATLAB Indexing or Referencing (Subscripted, Linear, and Logical)
• Array Operations vs. Matrix Operations
• Boolean Array Operations
• Constructing Matrices from Vectors
• Utility Functions
  

Extended Examples

• Matrix Functions of Two Vectors
• Ordering, Setting, and Counting Operations
• Sparse Matrix Structures
• Additional Examples
  

More Resources

• Matrix Indexing and Manipulation
• Matrix Memory Preallocation
• Maximizing MATLAB Performance
  

---

Tools and Techniques　
Section 1: MATLAB Indexing or ReferencingMATLAB 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 IndexingIn 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 IndexingIn 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 IndexingWith 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; endYou 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 Vectorsy(:,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)*10The 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 itA = 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.

zhouyang664

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(-x2 - y2)

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.