Operations on matrices
Matlab stands for 'matrix laboratory'. Not surprisingly, matrices, vectors and multidimensional arrays are at the heart of the language. Here we describe how to create, access, modify and otherwise manipulate matrices - the bread and butter of the Matlab programmer.
Contents
- Creating Matrices
- The Size of a Matrix
- Transposing a Matrix
- Sums and Means
- Concatenating Matrices
- Basic Indexing
- Logical Indexing
- Assignment
- Deletion
- Expansion
- Linear Indexing
- Reshaping and Replication
- Element-Wise Matrix Arithmetic
- Matrix Multiplication
- Solving linear systems
- More Linear Algebra
- Multidimensional Arrays
- Sparse Matrices
- Other numeric data types
- Other Useful Functions
Creating Matrices
There are a number of ways to create a matrix in Matlab. We begin by simply entering data directly. Entries on each row are separated by a space or comma and rows are separated by semicolons, (or newlines). We say that this matrix is of size 4-by-3 indicating that it has 4 rows and 3 columns. We, (and Matlab) always refer to rows first and columns second.
A = [1 3 5 ; 2 4 1 ; 3 3 3 ; 2 1 9]
A =
1 3 5
2 4 1
3 3 3
2 1 9
We can often exploit patterns in the entries to create matrices more succinctly.
A = 1:10 % start at 1, increment by 1, stop at 10 B = 1:2:10 % start at 1, increment by 2, stop at 10 C = 10:-1:3 % start at 10, decrement by 1, stop at 3
A =
1 2 3 4 5 6 7 8 9 10
B =
1 3 5 7 9
C =
10 9 8 7 6 5 4 3
We can also create an empty matrix.
D = [];
Alternatively, there are several functions that will generate matrices for us.
A = zeros(4,5) % 4-by-5 matrix of all zeros B = ones (2,3) % 2-by-3 matrix of all ones C = rand(3,3) % 3-by-3 matrix of uniform random numbers in [0,1] D = randn(2,5) % 2-by-5 matrix of standard normally distributed numbers [E,F] = meshgrid(1:5) % 5-by-5 grids of numbers G = eye(4) % 4-by-4 identity matrix H = diag(1:4) % 4-by-4 diagonal matrix I = logspace(0,2,6) % 1-by-6 matrix of log-spaced numbers from 10^0 to 10^2 J = blkdiag(rand(2,2),ones(3,2))% 5-by-4 block diagonal matrix K = tril(ones(3,4)) % 3-by-4 matrix whose lower triangular part is all ones. L = triu(ones(3,4)) % 3-by-4 matrix whose upper triangular part is all ones. M = magic(6) % Create a magic square in which sum(A,1) = sum(A,2) = trace(A) = trace(rot90(A))
A =
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
B =
1 1 1
1 1 1
C =
0.8147 0.9134 0.2785
0.9058 0.6324 0.5469
0.1270 0.0975 0.9575
D =
2.7694 3.0349 -0.0631 -0.2050 1.4897
-1.3499 0.7254 0.7147 -0.1241 1.4090
E =
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
F =
1 1 1 1 1
2 2 2 2 2
3 3 3 3 3
4 4 4 4 4
5 5 5 5 5
G =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
H =
1 0 0 0
0 2 0 0
0 0 3 0
0 0 0 4
I =
1.0000 2.5119 6.3096 15.8489 39.8107 100.0000
J =
0.9595 0.0357 0 0
0.6557 0.8491 0 0
0 0 1.0000 1.0000
0 0 1.0000 1.0000
0 0 1.0000 1.0000
K =
1 0 0 0
1 1 0 0
1 1 1 0
L =
1 1 1 1
0 1 1 1
0 0 1 1
M =
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
The functions true() and false(), act just like ones() and zeros() but create logical arrays whose entries take only 1 byte each rather than 32.
clear all
A = false(3,3)
B = zeros(3,3)
whos
A =
0 0 0
0 0 0
0 0 0
B =
0 0 0
0 0 0
0 0 0
Name Size Bytes Class Attributes
A 3x3 9 logical
B 3x3 72 double
The Size of a Matrix
We can determine the size of a matrix by using the size() command
[nrows,ncols] = size(A)
nrows =
3
ncols =
3
and the number of elements by using the numel() command.
n = numel(A)
n =
9
We refer to dimensions of size 1 as singleton dimensions. The length() command gives the number of elements in the first non-singleton dimension, and is frequently used when the input is a row or column vector; however, it can make code less readable as it fails to make the dimensionality of the input explicit.
We can also determine the size along a specific dimension with size().
n = size(A,2)
n =
3
Transposing a Matrix
A m-by-n matrix can be transposed into a n-by-m matrix by using the transpose operator '.
A = [1 2 3 4 ; 5 6 7 8] B = A'
A =
1 2 3 4
5 6 7 8
B =
1 5
2 6
3 7
4 8
Sums and Means
You can use the sum() and mean() functions to sum up or take the average of entries along a certain dimension.
C = sum(A,1) % Sum out dimension 1, (rows) D = sum(A,2) % Sum out dimension 2, (cols) E = mean(A,1) % Take the average along dimension 1, (rows)
C =
6 8 10 12
D =
10
26
E =
3 4 5 6
The [] argument to the min and max functions indicates that you will specify a dimension.
F1 = max(A, 2) % Larger of A and 2 elementwise F = max(A,[],2) % Find the max of each row (collapse dim 2) G = min(A,[],1) % Find the min of each column (collapse dim 1)
F1 =
2 2 3 4
5 6 7 8
F =
4
8
G =
1 2 3 4
Concatenating Matrices
Matrices can be concatenated by enclosing them inside of square brackets and using either a space or semicolon to specify the dimension. Care must be taken that the matrices are of the right size or Matlab will return an error.
A = [[1 2 3],ones(1,3)] % concatenate [1 2 3], [1 1 1] along columns B = [[1 2 3];ones(1,3)] % concatenate [1 2 3], [1 1 1] along rows C = [99 A 42] % add a number at the beginning or end of an array D = [A ; A] % duplicate the whole row
A =
1 2 3 1 1 1
B =
1 2 3
1 1 1
C =
99 1 2 3 1 1 1 42
D =
1 2 3 1 1 1
1 2 3 1 1 1
Basic Indexing
Individual entries can be extracted from a matrix by simply specifying the indices inside round brackets. We can also extract several entries at once by specifying a matrix, or matrices of indices or use the : operator to extract all entries along a certain dimension. The 'end' statement stands for the last index of a dimension.
A = magic(6); B = A(3,5); % extract the entry 3 rows down, 5 cols over C = A([1,2,3],4); % extract the entries (1,4) ; (2,4) ; (3,4) D = A(4,[1,1,1]); % extract the entry (4,1) three times E = A([2,5],[3,1]); % extract the entries (2,3) ; (2,1) ; (5,3) ; (5,1) F = A(:,4); % extract the fourth column G = A(4,:); % extract the fourth row H = A(:); % extract every entry as a column vector I = A(end,3); % extract the entry in the last row, 3rd column J = A(end-1,end-1); % extract the entry in the second to last row & col K = A(end-4:end,1); % extract the last three entries from the first col L = A(2:end,2:end); % extract everything except the first row and col M = A(end:-1:1,:); % extract everything with the order of the rows reversed. N = diag(A); % extract the main diagonal of A O = diag(rot90(A)); % extract the counter diagonal of A P = diag(A,-2) ; % extract the diagonal entries two diagonals left and below the main
Logical Indexing
We can also extract entries using a bit pattern, i.e. a matrix of logical values. Only the entries corresponding to true are returned. This can be particularly useful for selecting elements that satisfy some logical criteria such as being larger than a certain value. We can create a logical matrix by relating a numeric matrix to either a scalar value or matrix of the same size via one of the logical operators, < > <= >= == ~= or by a binary function such as isprime() or isfinite().
B = A > 30
B =
1 0 0 0 0 0
0 1 0 0 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 0 1 0 0 0
0 1 0 0 0 0
We can then use this logical matrix to extract elements from A. In the following line, we repeat the call to A > 30 but pass the result directly in, without first storing the interim result.
B1 = A(A > 30) % get all elements in A greater than 30 B = A(isprime(A) & (A > 30)) % get all prime elements in A greater than 30
B1 =
35
31
32
36
33
34
B =
31
We could also achieve the same result using the find() function, which returns the indices of all of the non-zero elements in a matrix. While this command is useful when the indices themselves are of interest, using find() can be slightly slower than logical indexing although it is a very common code idiom.
B2 = A(find(A > 30)) % same result as A(A>30) but calculated differently
B2 =
35
31
32
36
33
34
We can check that two matrices are equal, (i.e. the same size with the same elements) with the isequal() function. Using the == relation returns a matrix of logical values, not a single value.
test = isequal(B1,B2) test2 = all(B1==B2)
test =
1
test2 =
1
Assignment
Assignment operations, in which we change a value or values in a matrix, are performed in a very similar way to the indexing operations above. Both parallel and logical indexing can be used. We indicate which entries will be changed by performing an indexing operation on the left hand side and then specify the new values on the right hand side. The right must be either a scalar value, or a matrix with the same dimensions as the resulting indexed matrix on the left. Matlab automatically expands scalar values on the right to the correct size.
A(3,2) = 999; % assign 999 to entry (3,2) A(:,1:3:end) = 999; % assign 999 to every third column A(:,1) = [2;3;5;9;8;7]; % assign new values to the first column. A(A == 999) = 444; % assign all entries equal to 999 the value 444
We can assign every value at once by using the colon operator. The following command temporarily converts A to a column vector, assigns the values on the right hand side and converts back to the original dimensions.
A(:) = 1:36
A =
1 7 13 19 25 31
2 8 14 20 26 32
3 9 15 21 27 33
4 10 16 22 28 34
5 11 17 23 29 35
6 12 18 24 30 36
Recall from the indexing section that indices can be repeated returning the corresponding entry multiple times as in A([1,1,1],3). You can also repeat indices in assignments but the results are not what you might expect.
A = ones(3,5); A([1,2,3,1,1],1) = A([1,2,3,1,1],1) + 1
A =
2 1 1 1 1
2 1 1 1 1
2 1 1 1 1
You may have expected the entry A(1,1) to now have a value of 4 instead of 2 since we indexed entry A(1,1) three times. Matlab calculates the right hand side completely before assigning the values and so the value of 2 is simply assigned to A(1,1) three times.
Deletion
Assigning [] deletes the corresponding entries from the matrix. Only deletions that result in a rectangular matrix are allowed.
A([1,3],:) = [] % delete the first and third rows from A A(:,end) = [] % delete the last column from A
A =
2 1 1 1 1
A =
2 1 1 1
Expansion
When the indices in an assignment operation exceed the size of the matrix, Matlab, rather than giving an error, quietly expands the matrix for you. If necessary, it pads the matrix with zeros. Using this feature is somewhat inefficient, however, as Matlab must reallocate a sufficiently large chunk of contiguous memory and copy the array. It is much faster to preallocate the maximum desired size with the zeros command first, whenever the maximum size is known in advance: see here for details.
[nrows,ncols] = size(A) A(4,10) = 222
nrows =
1
ncols =
4
A =
2 1 1 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 222
A = 3 A(1:5,1:5) = 3
A =
3
A =
3 3 3 3 3
3 3 3 3 3
3 3 3 3 3
3 3 3 3 3
3 3 3 3 3
Linear Indexing
When only one dimension is specified in an indexing or assignment operation, Matlab performs linear indexing by counting from top to bottom and then left to right so that the last entry in the first column comes just before the first entry in the second column.
A = zeros(3,5); % create an empty matrix A(4) = 99 % assign 99 to the fourth entry (row 3, col 1) A(1:15) = 1:15 % assign vals 1:15 to their corresponding entries.
A =
0 99 0 0 0
0 0 0 0 0
0 0 0 0 0
A =
1 4 7 10 13
2 5 8 11 14
3 6 9 12 15
The functions ind2sub() and sub2ind() will convert from a linear index to regular indices and vice versa, respectively. In both cases, you must specify the size of the underlying matrix.
[rowNDX, colNDX] = ind2sub(size(A),12) linearNDX = sub2ind(size(A),rowNDX,colNDX)
rowNDX =
3
colNDX =
4
linearNDX =
12
Sometimes, when dealing with multi-dimensional arrays, it is annoying that these functions return/ require multiple separate arguments. We therefore provide the following alternative functions: ind2subv and subv2ind
ndx = ind2subv(size(A),12) linearNDX = subv2ind(size(A),ndx)
Undefined function or method 'ind2subv' for input arguments of type 'double'. Error in ==> matrixOperations at 217 ndx = ind2subv(size(A),12)
Reshaping and Replication
It is sometimes useful to reshape an array of size m-by-n to size p-by-q where m*n = p*q. The reshape() function lets you do just that. The elements are placed in such a way so as to preserve the order induced by linear indexing. In other words, if a(3) = 3 before reshaping, a(3) will still equal 3 after reshaping.
A = zeros(5,6); A(1:30) = 1:30 B = reshape(A,3,10) check = A(11) == B(11)
Further, the repmat() function can be used to tile an array m-by-n times.
A = [1 2 ; 3 4]
B = repmat(A,3,6) % copy A vertically 3 times and horizontally 6 times
Element-Wise Matrix Arithmetic
We can perform the arithmetical operations, addition, subtraction, multiplication, division, and exponentiation on every element of a matrix. In fact, we can also use functions such as sin(), cos(), tan(), log() , exp() , etc to operate on every entry but we will focus on the former list for now.
If one operand is a scalar value, an element-wise operation is automatically performed. However, if both operands are matrices, a dot must precede the operator as in .* , .^ , ./, and further, both matrices must be the same size.
A = [1 2 3; 4 5 6]; B = A + 1; % Add or subtract a scalar value from every entry C = 3.*A; % Multiply every entry by a scalar value D = A ./ 3; % Divide every entry by a scalar value E = A .^ 3; % Exponentiate every entry by a scalar value F = A - [2 4 8 ; 9 1 2]; % Add or subtract two matrices of the same size, (element-wise) G = A ./ B; % Divide every entry in A by the corresponding entry in B H = A .* B; % Multiply every entry in A by the corresponding entry in B I = A .^ B; % Exponentiate every entry in A by the corresponding entry in B
Matlab also has the .\ operator which is the same as the ./ operator with the order of the operands reversed so that (A ./ B) = (B .\ A). As this is infrequently used, it should be avoided for the sake of clarity.
Matrix Multiplication
We can also perform matrix multiplication of an m-by-n matrix and an n-by-p matrix yielding an m-by-p matrix. Suppose we multiple A*B = C, then C(i,j) = A(i,:)*B(:,j) , that is, the dot product of the ith row from A and the jth column from B. The dot, (or inner) product of a and b is just sum(a.*b). Further, if A is a square matrix, we can multiply A by itself k times by the matrix exponentiation A^k.
A = [1 2 3 4 ; 5 6 7 8] B = 2*ones(4,3) C = A*B % Matrix multiplication of A,B D = (A * A') / 100 % Matrix multiply A and A', divide every entry by 100 E = D ^ 4 % Matrix multiply D by itself 4 four times
Solving linear systems
Suppose we have the matrix equation Y = XW where X is an n-by-d matrix, W is a d-by-k matrix and thus Y is an n-by-k matrix. If X is invertible, we could solve for W= inv(X)*Y. The inv() function returns the inverse of a matrix. If n ~= d, however, we can still solve for the least squares estimate of W by taking the pseudo inverse of X, namely inv(X'*X)*X', or more concisely using the Matlab function, pinv(X). Matlab allows you to solve more directly, (and efficiently) for W, (i.e. the lsq estimate of W), however, by using the matrix division X \ Y. Both X and Y must have the same number of rows. (Below we specify a seed to the random number generators so that they return the same values every time this demo is run).
seed = 1; rand('twister',seed);randn('state',seed); X = randn(100, 8); w = randn(8,1); % loading vector, Wtrue y = X*w; % Target or output variable, (no noise)
Wlsq = X \ y; % recommended way Wlsq1 = inv(X'*X)*X'*y; % normal eqns - Not the recommended way Wlsq2 = pinv(X)*y; % Not the recommended way [w'; Wlsq'; Wlsq1'; Wlsq2'] % In this simple example, all same
Matlab also supports matrix right division such that X \ Y = (Y' / X')' but as this is infrequently used, it should be avoided for the sake of clarity.
More Linear Algebra
Matlab was original designed primarily as a linear algebra package, and this remains its forte. Here we only list a few functions for brevity, do not display the results. Many functions can also take additional arguments: type doc svd for instance to see documentation for the svd() function.
A = magic(5); B = det(A); % The determinant of A C = rank(A); % The rank of A E = trace(A); % sum of diagonal entries G = orth(A); % an orthonormal basis for range of A H = null(A); % an orthonormal basis for the nullspace of A I = chol(A*A'); % Cholesky decomposition s.t. if R = chol(X) then R'*R = X [evecs,evals] = eig(A); % eigen vectors and values of A, (use eigs on sparse matrices) [U,S,V] = svd(A); % singular value decomposition s.t. A = U*S*V' [Q,R] = qr(A); % QR decomposition of A
Multidimensional Arrays
Numeric matrices in Matlab can extend to an arbitrary number of dimensions, not just 2. We can use the zeros(), ones(), rand(), randn() functions to create n-dimensional matrices by simply specifying n parameters. We can also use repmat() to replicate matrices along any number of dimensions, or the cat() function, which is a generalization of the [] concatenation we saw earlier. Indexing, assignment,and extension work just as before, only with n indices, as opposed to just two. Finally, we can use functions like sum(), mean(), max() or min() by specifying the dimension over which we want the function to operate. sum(A,3) for example, sums over, or marginalizes out, the 3rd dimension.
A = ones(3,5,9,2,2); % a 5 dimensional array B = cat(3,[1 2 ; 4 5],[3 2 ; 1 1]) % concatenate two matrices along the third dimension C = repmat([1 2 ; 3 4],[2,2,2,2]); % tile a matrix in 4D, twice per dimension D = C(1,1,1,1) % retrieve the entry (1,1,1,1) A(1,2,2,1) = 99; % assign an entry F = mean(C,4); % Take the mean along the 4th dimension G = max(A,[],3); % Take the maximum along the 3rd dimension
Taking the mean of say a 4-by-4-by-2-by-2 matrix along the 3rd dimension results in a matrix of size 4-by-4-by-1-by-2. If we want to remove the 3rd singleton dimension, (which is only acting now as a place holder) we can use the squeeze() function.
size(C)
E1 = mean(C,3); size(E1) E2 = squeeze(E1); size(E2)
The ndims() functions indicates how many dimensions an array has. Final singleton dimensions are ignored but singleton dimensions occurring before non-singleton dimensions are not.
F = ndims(A) G = ndims(squeeze(mean(A,2)))
The meshgrid() function we saw earlier extends to 3 dimensions. If you need to grid n-dimensional space, use the ndgrid() function but keep in mind that the number of elements grows exponentially with the dimension.
a = 1:10;
n = numel(ndgrid(a,a,a,a,a,a)) % 1 million entries.
Sparse Matrices
When dealing with large matrices containing many zeros, you can save a great deal of space by using Matlab's sparse matrix construct. Sparse matrices can be used just like ordinary matrices but can be slower depending on the operation. The functions full() and sparse() convert back and forth. Currently Matlab supports double and logical sparse matrices.
A = zeros(100,100); A([1,4,8],[7,9,33]) = reshape(1:9,3,3); n = nnz(A) % The number of non-zero entries in A nzeros = nonzeros(A); % All of the non-zero entries in one big column vector A = sparse(A) % Convert to a sparse matrix check = issparse(A) % Is is really a sparse data type? B = A * rand(100,1); % Perform operations as you would with a non sparse matrix C = full(A); % Convert back to a full matrix.
The spy() function can be used to visualize the sparsity pattern of a matrix.
The spalloc() function can be used to preallocate space for a sparse matrix. The following command creates a 100-by-100 matrix with room currently for 10 non-zero elements. More than 10 non-zero elements can be added later but this can be slow as Matlab will need to find a larger chunk of memory and copy the non-zero elements.
A = spalloc(100,100,10)
n = nzmax(A) % how many non-zero elements do I have room for?
Other numeric data types
Matlab has limited support for 11 numeric data types similar to those in the C programming language. Below we create matrices of each type and show the space each matrix requires. Matrices can also be created by using the commands int8() , single() , int64() etc. The cast() command, converts from one data type to another. You can determine the class of a variable with the class() command and the maximum or minimum values each class is able to represent with the intmax() , intmin() , realmax() , and realmin() functions. The uint classes are unsigned and not able to represent negative numbers. Unfortunately many Matlab functions do not support types other than double or logical. Functions such as sum() have an optional parameter 'native', which performs summation without automatically casting to double. To perform variable precision arithmetic, check out the vpa() function available in the symbolic math toolbox.
clear A = zeros(100,100,'double'); % same as zeros(100,100) B = zeros(100,100,'int64'); % signed 64 bit integer C = zeros(100,100,'uint64'); % unsigned 64 bit integer D = zeros(100,100,'single'); % single precision number E = zeros(100,100,'int32'); % signed 32 bit integer F = zeros(100,100,'uint32'); % unsigned 32 bit integer G = zeros(100,100,'int16'); % signed 16 bit integer H = zeros(100,100,'uint16'); % unsigned 16 bit integer I = zeros(100,100,'int8'); % signed 8 bit integer J = zeros(100,100,'uint8'); % unsigned 8 bit integer K = false(100,100); % logical array whos % display size of variables AA = realmax('double'); % max sizes representable by different types BB = intmax('int64'); CC = intmax('uint64'); DD = realmax('single'); EE = intmax('int32'); FF = intmax('uint32'); GG = intmax('int16'); HH = intmax('uint16'); II = intmax('int8'); JJ = intmax('uint8'); K = int32(100); % construct an int32 directly check1 = class(K); % check its class L = cast(K,'double'); % cast to a double() - also see typecast() check2 = class(L); % check its class M = int8([3,1,2,1,4]); % create several int8s N = sum(M,'native'); % sum ints in 'native' mode, i.e. don't cast to double O = sum(M); % don't sum in 'native' mode to see the difference class1 = class(N); % check the class type class2 = class(O); % check the class type
Other Useful Functions
The cumsum() and cumprod() functions can be useful for generating a running sum or product of an array. The diff() function returns the differences between consecutive elements. You can specify the dimension over which you want them to operate. If you leave this blank, they operate over the first non-singleton dimension.
A = cumsum(1:6) C = cumprod(1:6) D = diff(A)
The histc function is useful for, (among other things) counting the number of occurrences of numbers in an array.
A = sort(floor(10*rand(1,10))+1) % random ints from 1 to 10 counts = histc(A,1:10) % count how often each int occurs
The filter() function can be used to calculate values that depend on previous values in an array. While it is quite a complicated function, here is an easy way to calculate the points halfway between each consecutive point in an array. The first result is just half the value of the first element. You can calculate a running average in which only a window of k elements are included with filter(ones(1,k)/k,1,data).
A = 1:10 B = filter([0.5,0.5],1,A);