CHAPTER 2: MATRIX COMPUTATIONS AND LINEAR ALGEBRA

Lecture 2.1: Two-dimensional arrays (matrices) and matrix operations

In linear algebra, a matrix is a rectangular array of elements (numbers). If a matrix A = [ ai,j ] has m rows and n columns, it has m-by-n elements ai,j  where the first index i  runs for rows of A and the second index j runs for columns of A. MATLAB matrices are direct and convenient numerical realization of matrices in linear algebra.

A = [ 0.1 0.2 0.3 0.4 ; -0.1,-0.2,-0.3,-0.4 ; 0 1 0 1 ]

% Notice that rows in MATLAB matrices are separated by the sign ";"

% Notice that elements at each row can be separated by spaces and commas

A =

0.1000    0.2000    0.3000    0.4000

-0.1000   -0.2000   -0.3000   -0.4000

0    1.0000         0    1.0000

% Number of elements in each row must be identical

??? Error using ==> vertcat

All rows in the bracketed expression must have the same

number of columns.

The colon sign ":" represents all elements of a row or a column of a MATLAB matrix.

ans =

0.1000

-0.1000

0

ans =

0.4000

-0.4000

1.0000

ans =

0.1000    0.2000    0.3000    0.4000

ans =

0     1     0     1

ans =

0.2000

ans =

0

Matrix operations:

·         all vector operations, such as pointwise additions, subtractions, multiplications, divisions, and powers of matrices of the same size

A = [ 1 2 3 4 ; -1,-2,-3,-4 ];

B = [ 4 3 2 1 ; -4,-3,-2,-1 ];

C1 = A + B, C2 = A.*B, C3 = A./B, C4 = A.^B

% Pointwise matrix operations are preferable compared to double loops in i,j

% because of (i) clear mathematical notations, (ii) computational efficiency

C1 =

5     5     5     5

-5    -5    -5    -5

C2 =

4     6     6     4

4     6     6     4

C3 =

0.2500    0.6667    1.5000    4.0000

0.2500    0.6667    1.5000    4.0000

C4 =

1.0000    8.0000    9.0000    4.0000

1.0000   -0.1250    0.1111   -0.2500

·         size: outputs numbers of rows and columns

nRow =     2

nCol =     4

·         transpose: A', outputs a transposed matrix AT

B =  1    -1

2    -2

3    -3

4    -4

ans =

4     2

·         elementary matrices:

• zeros: outputs a null matrix with all elements being zeros

zeros(2,6) % a null 2-by-6 matrix

zeros(3) % a null 3-by-3 matrix

ans =    0     0     0     0     0     0

0     0     0     0     0     0

ans =    0     0     0

0     0     0

0     0     0

• ones: outputs an matrix with all elements being ones

ones(2,6) % a 2-by-6 matrix of ones

ones(3) % a 3-by-3 matrix of ones

ans =

1     1     1     1     1     1

1     1     1     1     1     1

ans =

1     1     1

1     1     1

1     1     1

• eye: outputs an identity matrix (a diagonal matrix)

ans =

1     0     0     0     0     0

0     1     0     0     0     0

ans =

1     0     0

0     1     0

0     0     1

• rand, randn: outputs a random matrix with all elements being random numbers

rand(2,6) % random numbers are uniformly distributed in the interval [0,1]

randn(3) % random numbers are normally distributed with mean zero and variance one

ans =

0.9501    0.6068    0.8913    0.4565    0.8214    0.6154

0.2311    0.4860    0.7621    0.0185    0.4447    0.7919

ans =

-0.4326    0.2877    1.1892

-1.6656   -1.1465   -0.0376

0.1253    1.1909    0.3273

·         special matrices:

• hilb: outputs the Hilbert matrix with elements ai,j = 1/(i+j-1)

C1 =   1.0000    0.5000    0.3333    0.2500    0.2000

0.5000    0.3333    0.2500    0.2000    0.1667

0.3333    0.2500    0.2000    0.1667    0.1429

0.2500    0.2000    0.1667    0.1429    0.1250

0.2000    0.1667    0.1429    0.1250    0.1111

• pascal: outputs the Pascal matrix with integer entries, made up from Pascal's triangle: ai,j = ai-1,j + ai,j-1

C2 = pascal(6)

% Pascal matrices produce a table of binomial coefficients

% Coefficients at the k-th anti-diagonal of Pascal matrices produce

% coefficients of the binomial expansion: (x + y)^(k-1)

C2 =      1     1     1     1     1     1

1     2     3     4     5     6

1     3     6    10    15    21

1     4    10    20    35    56

1     5    15    35    70   126

1     6    21    56   126   252

• vander: outputs the Vandermonde matrix whose columns are powers of a vector, i.e. ai,j = vin-j

C3 =       1           1           1           1           1           1

32          16           8           4           2           1

243          81          27           9           3           1

1024         256          64          16           4           1

3125         625         125          25           5           1

7776        1296         216          36           6           1

• toeplitz: outputs the Toeplitz matrices

r = [ 1 2 3 4 5 6 ]; s = [ -1,-2,-3,-4, -5];

C4 = toeplitz(r,s) % r and s are the column and row vectors of C4

C5 = toeplitz(r)   % the symmetric Toeplitz matrix

rInv = r(length(r):-1:2); rInv = [r(length(r)),rInv];

C6 = toeplitz(r,rInv)' % the circulant matrix

Warning: Column wins diagonal conflict.

> In c:\progra~1\matlab6\toolbox\matlab\elmat\toeplitz.m at line 18

C4 =

1    -2    -3    -4    -5

2     1    -2    -3    -4

3     2     1    -2    -3

4     3     2     1    -2

5     4     3     2     1

6     5     4     3     2

C5 =

1     2     3     4     5     6

2     1     2     3     4     5

3     2     1     2     3     4

4     3     2     1     2     3

5     4     3     2     1     2

6     5     4     3     2     1

Warning: Column wins diagonal conflict.

> In c:\progra~1\matlab6\toolbox\matlab\elmat\toeplitz.m at line 18

C6 =

1     2     3     4     5     6

6     1     2     3     4     5

5     6     1     2     3     4

4     5     6     1     2     3

3     4     5     6     1     2

2     3     4     5     6     1

• hadamard, hankel: more specialized matrices

·         matrix-vector multiplications

A =  1     2     3     1

0     1     4     2

3     0     2     3

x =   1

2

1

1

% the number of columns of A must be equal to the number of rows of x

y1 = A*x  % MATLAB matrix multiplication operator

% a row-oriented loop:

[n,m] = size(A);

for k = 1:n,     y2(k) = A(k,:)*x; end

y2

% a column-oriented loop

y3 = zeros(n,1);

for j = 1 : m,   y3 = y3 + A(:,j)*x(j); end

y3

y1 =  9

8

8

y2 =     9     8     8

y3 =   9

8

8

·         matrix-matrix mutiplications:

If A is a n-by-m matrix, B is a m-by-q matrix, then the matrix multiplication makes sense and C = A*B is a n-by-q matrix with elements: ci,j = ai,k bk,j . The other product B*A is defined if n = q. For a square matrices of the same size, both products are defined but they are not equal in a general case: A*B B*A.

A = [ -2 1 0 ; 1,-2,1 ; 0,1,-2];

B = [1 0 1 ; 0 1 0; 1 0 1];

C1 = A*B, C2 = B*A

C1 =

-2     1    -2

2    -2     2

-2     1    -2

C2 =

-2     2    -2

1    -2     1

-2     2    -2

n = 3;

for k = 1 : n

for j = 1 : n

C3(k,j) = A(k,:)*B(:,k); % the dot product algorithm

end

end

C3

C4 = zeros(n,n);

for j = 1 : n

for p = 1 : n

C4(:,j) = C4(:,j) + A(:,p)*B(k,j); % the column-oriented algorithm

end

end

C4

C5 = zeros(n,n);

for p = 1 : n

C5 = C5 + A(:,k)*B(k,:); % the outer-product oriented algorithm

end

C5

C3 =

-2    -2    -2

-2    -2    -2

-2    -2    -2

C4 =

-1     0    -1

0     0     0

-1     0    -1

C5 =

0     0     0

3     0     3

-6     0    -6