Loren on the Art of MATLAB

Contents

Note: I am not discussing efficiency in this article. It is highly dependent on the number of elements in the original array and how many will be retained or thrown out. This article focuses on specifying what to keep or replace.

General Setup

Here's the setup for this investigation. I will use a fixed matrix for all the methods and always end up with the same final output. The plan is to show you multiple ways to get the result, since different methods may be appropriate under different circumstances.

A = magic(17);
Result = A;
Result( A < mean(A(:)) ) = 0;

Let's look at the nonzero pattern of Result using spy.

spy(Result)

Method #1 - Using Subscripts of Keepers

Here's a list of the subscripts for the elements to keep unchanged.

[rA,cA] = find(A > (17^2)/2);

Next we convert the subscripts to indices.

Result1 = zeros(size(A));
indices = sub2ind(size(A),rA,cA);
Result1(indices) = A(indices);
isequal(Result, Result1)

ans =
     1

Why did I convert subscripts to indices? Let me illustrate with a very small example.

matrix = [ -1 1 0; 2 0 -2; 0 3 -3]
[rows,cols] = find(matrix==0)

matrix =
    -1     1     0
     2     0    -2
     0     3    -3
rows =
     3
     2
     1
cols =
     1
     2
     3

Now let's see what I get if I use the subscripts to address the selected elements:

matrix(rows,cols)

ans =
     0     3    -3
     2     0    -2
    -1     1     0

I get the full matrix back, even though I selected only 3 elements. This definitely surprised me when I first encountered this. What's happening?

MATLAB matches each row element with each column element. matrix([1 2 3],2) returns the elements from rows 1 through 3 in column 1.

matrix(1:3,2)

ans =
     1
     0
     3

To learn more about indexing in general, you might want to read these posts or search the MATLAB documentation.

Method #2 - Using Indices of Keepers

Here we used the single output form of find which returns indices instead of subscripts.

indA = find(A > (17^2)/2);
Result2 = zeros(size(A));
Result2(indA) = A(indA);
isequal(Result, Result2)

ans =
     1

Method #3 - Using Logical Keepers

We'll try keeping about half of the elements unchanged.

keepA = (A > (17^2)/2);
Result3 = zeros(size(A));
Result3(keepA) = A(keepA);
isequal(Result, Result3)

ans =
     1

keepA is a logical matrix the same size as A. I use logical indexing to populate Result3 with the chosen values from A.

Method #4 - Subscripts for Elements to Set to Zero

If instead we have a list of candidates to set to 0, we have an easier time since we don't need to start off with a matrix of zeros. Instead we start with a copy of A.

Result4 = A;
[rnotA,cnotA] = find(A <= (17^2)/2);

Convert indices to subscripts, as in method #1.

indices = sub2ind(size(A),rnotA,cnotA);

Now zero out the selected matrix elements.

Result4(indices) = 0;
isequal(Result, Result4)

ans =
     1

Method #5 - Indices for Elements to Set to Zero

If we're instead given indices, we simply skip the step of converting subscripts and follow similar logic to that in method #4.

Result5 = A;
indnotA = find(A <= (17^2)/2);
Result5(indnotA) = 0;
isequal(Result, Result5)

ans =
     1

Method #6 - Using Logical Arrays to Specify Zero Elements

Finally, if we have a mask for the values to set to 0, we simply use it to select and set elements.

Result6 = A;
keepnotA = (A <= (17^2)/2);
Result6(keepnotA) = 0;
isequal(Result, Result6)

ans =
     1

Which Method(s) Do You Prefer?

Which method or methods do you naturally find yourself using? Do you ever invert the logic of your algorithm to fit your way of thinking about addressing the data (the ins or the outs)? Please post your thoughts here. I look forward to seeing them.

Get the MATLAB code

Published with MATLAB® 7.6

Contents

Two Conventions

When calling functions in MATLAB that create or reference arrays, there are two different conventions for specifying the first two dimensions of an array:

These different conventions get used depending on whether the array is strictly an array or when it has a probable physical interpretation in terms of Cartesian grid.

Matrix versus Cartesian Frameworks

In MATLAB, we index into matrices using row and column indices. That is, for matrix

A = magic(3)

A =
     8     1     6
     3     5     7
     4     9     2

I can reference the last element in the second row like this:

A(2,3)

ans =
     7

or

A(2,end)

ans =
     7

and to retrieve the first element of columns 1 and 3, I type this:

A(1,[1 3])

ans =
     8     6

Check out this post to get a fairly complete picture of indexing in MATLAB.

Plotting Matrices

The options xy and ij for axis allow us to control whether to view the data just the way the matrix is displayed (the (1,1) element is in the upper left corner):

colormap(jet(9))
image(A)
colorbar

or we can view the image/matrix as a top-down view of a physical "map" with the y-values now increasing along the vertical axis.

axis xy

In images, the matrices column values plot vertically. In line plots, the columns are plotted vs. their indices, making columns plot horizontally.

h = plot(A);
legend(h,{'Column 1','Column 2','Column 3'})

Creating an N-Dimensional Grid

ndgrid uses matrix indexing since it is designed to suit working in any number of array dimensions.

[ii, jj] = ndgrid(1:3,1:2)

ii =
     1     1
     2     2
     3     3
jj =
     1     2
     1     2
     1     2

Creating an 2- or 3-Dimensional Grid

On the other hand, mesh, interprets its first inputs in a Cartesian framework since it is used primarily for working with 2- and 3-dimensional plots, often depicting something over an area or volume.

[xx, yy] = meshgrid(1:3,1:2)

xx =
     1     2     3
     1     2     3
yy =
     1     1     1
     2     2     2

In the Remarks section of the documentation for meshgrid, you will find a discussion of this, especially with respect to inverting the first two input arguments compared to ndgrid. There is related discussion in the Description section for the function mesh.

Functions designed for two or three dimension, such as contour and quiver, usually work in Cartesian coordinates. Functions that are designed to work in N-dimensions generally work with a matrix indexing framework. Be sure to check the help for functions you use to see which convention applies.

Thoughts?

Any ideas on what we can do, perhaps in examples or documentation (or eslewhere?) to help people out with this potential confusion? If so, please post here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Uses of :

There are several ways that you can use the : operator in MATLAB, and they serve different purposes. Most, but not all of them, have to do with indexing into an array.

Creating a List of Numbers

You can use the : operator to create a vector of evenly-spaced numbers. Here are the integers from -3 to 3.

list1 = -3:3

list1 =
    -3    -2    -1     0     1     2     3

Here are the first few odd positive integers.

list2 = 1:2:10

list2 =
     1     3     5     7     9

Here's how to divide the interval between 0 and pi into equally spaced samples.

nsamp = 5;
sliceOfPi = (0:1/(nsamp-1):1)*pi

sliceOfPi =
         0    0.7854    1.5708    2.3562    3.1416

Note that even though only a few digits of each value in sliceOfPi are printed out, they are double precision numbers. To see more digits, check out the function format.

Collapsing Trailing Dimensions

I have a 4-dimensional array and would like to find the sum of all the elements that are in the final row and column.

b = rand(3,2,6,4);
b32sum = sum(b(3,2,:))

b32sum =
   11.4622

If instead, I sum using all the dimensions, in this case, I get a 1x1x1x4 array instead of scalar output. If I sum these output values, I get my overall sum.

b32sumAll = sum(b(3,2,:,:))
overall = sum(b32sumAll)

b32sumAll(:,:,1,1) =
    2.8279
b32sumAll(:,:,1,2) =
    3.4177
b32sumAll(:,:,1,3) =
    1.9352
b32sumAll(:,:,1,4) =
    3.2815
overall =
   11.4622

Creating a Column Vector

Here's the size of the array b32sumAll.

size(b32sumAll)

ans =
     1     1     1     4

If I want to transform this into a column vector, I can use reshape, permute or : on the right-hand side.

b32vecR = reshape(b32sumAll,[],1);
b32vecP = permute(b32sumAll, [4 1:3]);
b32vec = b32sumAll(:)
allthesame = isequal(b32vec, b32vecP, b32vecR)

b32vec =
    2.8279
    3.4177
    1.9352
    3.2815
allthesame =
     1

Retaining Array Shape During Assignment

Now let me pour some new values into b32sumAll, the array that looks almost like a vector.

b32sumAll(:) = [1 2; 3 5]

b32sumAll(:,:,1,1) =
     1
b32sumAll(:,:,1,2) =
     3
b32sumAll(:,:,1,3) =
     2
b32sumAll(:,:,1,4) =
     5

Notice that I only have to have the same number of elements on both the left- and right-hand sides. The values are poured in from the right-hand side ordered as if that array had been turned into a column vector. Note that a scalar for the right-hand side is also acceptable and then MATLAB performs a scalar expansion to fill the left-hand side.

Working with All the Entries in Specified Dimensions

If I want to manipulate values in some specific dimensions, I can also use the : operator to specify the dimensions I'd like to leave alone. For example, suppose I want to perform a left shift on the values in the second dimension of my 3-D array. Let me first create an array for illustration.

a3 = zeros(2,3,2);
a3(:) = 1:numel(a3)

a3(:,:,1) =
     1     3     5
     2     4     6
a3(:,:,2) =
     7     9    11
     8    10    12

Now let's shift the column values all over to the left, and have the right-most one become the last column. Note that columns are dimension 2. Here's a way to do this.

a3r1 = a3(:,[2:size(a3,2) 1],:)

a3r1(:,:,1) =
     3     5     1
     4     6     2
a3r1(:,:,2) =
     9    11     7
    10    12     8

How do I do this if the array could be any number of dimensions, not just 3? I can't index with :, because I don't know how many : to use, since I don't knwo the dimension. In this case, use the string value ':', create the indices in a cell array, and use the notion of the comma-separated list to perform the actual indexing. For a3, here's what this looks like.

indices = {':', [2:size(a3,2) 1],':'}
a3rnew = a3(indices{:})
sameshifts = isequal(a3r1, a3rnew)

indices = 
    ':'    [1x3 double]    ':'
a3rnew(:,:,1) =
     3     5     1
     4     6     2
a3rnew(:,:,2) =
     9    11     7
    10    12     8
sameshifts =
     1

Now, more generally, when the number of dimensions is not known, I simply create a cell array of all ':' string values, replace the specific dimensions I care about with some other indexing expression, and then I can use the index cell array to do the operation.

a5 = zeros(2,5,2,3,4);
a5(:) = 1:numel(a5);
indices = repmat({':'}, 1, ndims(a5));
myshift = [2:size(a5,2) 1]
indices{2} = myshift
a5new = a5(indices{:});

myshift =
     2     3     4     5     1
indices = 
    ':'    [1x5 double]    ':'    ':'    ':'

Let's compare two "rows" and see that the values have shifted left by 1 column.

a5(1,:,2,2,3)
a5new(1,:,2,2,3)

ans =
   151   153   155   157   159
ans =
   153   155   157   159   151

Related Reference Material

For more in depth understanding of the various ways in which the : operator can be used, look in the MATLAB documention

and in related blog articles

: is Complicated

The : operator behaves differently in MATLAB depending on its usage. This can lead to confusion. Have I missed anything major on this topic? Let me know.

Get the MATLAB code

Published with MATLAB® 7.3

Contents

Dimensions

Arrays in MATLAB are N-dimensional, with an infinite number of trailing singleton dimensions. Trailing singleton dimensions past the second are not displayed or reported on, e.g., with size. No array has fewer than two dimensions. Empty arrays are the logical extension of general arrays but have at least one size 0 dimension. The syntax [] denotes the empty array of size 0x0 and is treated as an exception occasionally (for backward compatibility); it is also used specially during assignment.

Right-Hand Side

A = rand(2,3,2,1,1,1)
ndA = ndims(A)
sizeA = size(A)
B = A(1,2:3,2)
sizeB = size(B)

A(:,:,1) =
    0.2974    0.6932    0.9830
    0.0492    0.6501    0.5527
A(:,:,2) =
    0.4001    0.6252    0.3759
    0.1988    0.7334    0.0099
ndA =
     3
sizeA =
     2     3     2
B =
    0.6252    0.3759
sizeB =
     1     2

A(1,2,1,1)

ans =
    0.6932

try
    A(1,1,1,3)
catch
    s = lasterror();
    s.message
end

ans =
Error using ==> evalin
Index exceeds matrix dimensions.

A = 1:(3*4*7);
A = reshape(A,1,3,4,1,7);
sizeA = size(A)

sizeA =
     1     3     4     1     7

Show the middle chunk and index into it for one value.

A(:,:,10:14)
A(1,2,12)

ans(:,:,1) =
    28    29    30
ans(:,:,2) =
    31    32    33
ans(:,:,3) =
    34    35    36
ans(:,:,4) =
    37    38    39
ans(:,:,5) =
    40    41    42
ans =
    35

A = magic(3)
B = A(:,2)
C = A(1:2,:)

A =
     8     1     6
     3     5     7
     4     9     2
B =
     1
     5
     9
C =
     8     1     6
     3     5     7

V = 1:3;
V(end-1) = 7

V =
     1     7     3

A = magic(5)
A(1:3,1:2:5)

A =
    17    24     1     8    15
    23     5     7    14    16
     4     6    13    20    22
    10    12    19    21     3
    11    18    25     2     9
ans =
    17     1    15
    23     7    16
     4    13    22

A = 1:10
B = [1 3; 2 4]
C = A(B)

A =
     1     2     3     4     5     6     7     8     9    10
B =
     1     3
     2     4
C =
     1     3
     2     4

When both A and B are vectors, the number of elements in C is the number of elements in B and with orientation of A.

row = 1:4
column = (1:4)'
row(column)
column(row)

row =
     1     2     3     4
column =
     1
     2
     3
     4
ans =
     1     2     3     4
ans =
     1
     2
     3
     4

A = magic(6)
ind1 = [1 4]
ind2 = [1 2; 2 6]
C = A(ind1,ind2)

A =
    35     1     6    26    19    24
     3    32     7    21    23    25
    31     9     2    22    27    20
     8    28    33    17    10    15
    30     5    34    12    14    16
     4    36    29    13    18    11
ind1 =
     1     4
ind2 =
     1     2
     2     6
C =
    35     1     1    24
     8    28    28    15

A = magic(3)
B = [true false; true true]
C = A(B)
Arow = A(:)'
C = Arow(B)

A =
     8     1     6
     3     5     7
     4     9     2
B =
     1     0
     1     1
C =
     8
     3
     1
Arow =
     8     3     4     1     5     9     6     7     2
C =
     8     3     1

Left-Hand Side

A = magic(3)
A(:,2) = 17
A(1,:) = 2:2:6

A =
     8     1     6
     3     5     7
     4     9     2
A =
     8    17     6
     3    17     7
     4    17     2
A =
     2     4     6
     3    17     7
     4    17     2

A = 1:3
A(3,[2 4]) = 17

A =
     1     2     3
A =
     1     2     3     0
     0     0     0     0
     0    17     0    17

A = single(magic(2))
A(1) = 17
class(A)

A =
     1     3
     4     2
A =
    17     3
     4     2
ans =
single

A = magic(3)
A([1 3],:) = []

A =
     8     1     6
     3     5     7
     4     9     2
A =
     3     5     7

A = 1:9;
A(:) = magic(3)

A =
     8     3     4     1     5     9     6     7     2

References

Here are some references to documentation and other articles related to indexing.

Let me know if it's helpful having all the indexing information in one place.

Published with MATLAB® 7.2