Math 584, Mathematics of Medical Imaging

Matlab Lessons

This is a simple tutorial for using Matlab for Math 584, Mathematics of Medical Imaging.

The black text in bold is what you should type at the Matlab command prompt. The non-bold black text which follows is the expected Matlab output. The red text will explain what you should be learning. Please experiment by changing the commands so you can get a feel for the language of Matlab.

0. Some basic facts:

Type help <command> at the Matlab command prompt to see useful information about certain Matlab commands.
>> help zeros

ZEROS Zeros array.
    ZEROS(N) is an N-by-N matrix of zeros.

    ZEROS(M,N) or ZEROS([M,N]) is an M-by-N matrix of zeros.

    ZEROS(M,N,P,...) or ZEROS([M N P ...]) is an M-by-N-by-P-by-... array of
    zeros.

    ZEROS(SIZE(A)) is the same size as A and all zeros.

    ZEROS with no arguments is the scalar 0.

    ZEROS(M,N,...,CLASSNAME) or ZEROS([M,N,...],CLASSNAME) is an
    M-by-N-by-... array of zeros of class CLASSNAME.

    Example:
       x = zeros(2,3,'int8');

    See also eye, ones.

To supress output terminate a command with a semicolon ;

>> X= -1:0.1:1

X =

Columns 1 through 4

-1.00000000000000 -0.90000000000000 -0.80000000000000 -0.70000000000000

Columns 5 through 8

-0.60000000000000 -0.50000000000000 -0.40000000000000 -0.30000000000000

Columns 9 through 12

-0.20000000000000 -0.10000000000000                  0   0.10000000000000

Columns 13 through 16

   0.20000000000000   0.30000000000000   0.40000000000000   0.50000000000000

Columns 17 through 20

   0.60000000000000   0.70000000000000   0.80000000000000   0.90000000000000

Column 21

   1.00000000000000

versus

>> X=-1:0.1:1;
>>

This syntax produces no output in the command window.

Matlab Lesson 1: Linear Algebra

Defining matrices, and vectors

Define a 3x3 zero matrix and store it in the variable M

>> M=zeros(3,3)

M =

0 0 0

Similarly, the command ones(m,n) defines an mxn matrix of all ones.

Now, let's set the diagonal entries of M equal to one.

Type "M" at the end of the statement to display the matrix

>> M(1,1)=1; M(2,2)=1; M(3,3)=1; M

M =

1 0 0

0 1 0

0 0 1

This is one way to define the 3x3 identity matrix. Here is another way:

>> B=diag([1 1 1])

B =

1 0 0

0 1 0

0 0 1

Actually there's an even quicker way to define the 3x3 identity matrix:

>> B2=eye(3)

B2 =

1 0 0

0 1 0

0 0 1

The expression [1 1 1] represents a 1x3 row vector. To verify this, type

>> size([1 1 1])

ans =

1 3

To define a column vector, use semicolons between the entries:

>> [1; 2; 3; 4]

ans =

Or you can turn a row vector into a column vector by using the transpose operator '

>> row_vector=[1 2 3 4 5]; column_vector=row_vector.'

column_vector =

Of course, the transpose operator works on matrices of arbitrary size.

You can also define a row vector using an arithmetic sequence. For example:

>> X =0:0.2:pi;

produces a vector X starting at 0 and increasing each time by 0.2 until it reaches pi. Many functions defined in Matlab can be applied to vectors with the output a vector of the values of the function. For example

>> VCos=cos(0:0.2:pi)

VCos =

Columns 1 through 7

1.0000 0.9801 0.9211 0.8253 0.6967 0.5403 0.3624

Columns 8 through 14

0.1700 -0.0292 -0.2272 -0.4161 -0.5885 -0.7374 -0.8569

Columns 15 through 16

-0.9422 -0.9900

The expression 0:0.2:pi represents the row vector consisting of sixteen evenly spaced numbers between 0 and pi (starting at 0 with step size 0.2).

You can plot vectors in Matlab. Try typing

>> figure; plot(0:0.2:pi,cos(0:0.2:pi))

In general, if X=(x₁,....,x_n) and Y=(y₁,....,y_n) are two vectors containing real numbers of the same length and type (either both columns or both rows), then the command plot(X,Y) produces a plot of the pairs { (x_j.y_j) j=1,...n}. Matlab also has some built in functions for defining special classes of matrices. For example, the 4x4 Hilbert matrix:

>> hilb(4)

ans =

1.0000 0.5000 0.3333 0.2500

0.5000 0.3333 0.2500 0.2000

0.3333 0.2500 0.2000 0.1667

0.2500 0.2000 0.1667 0.1429

For more information, try the following Matlab help topics: zeros, ones, size, length, transpose, eye, hilb, toeplitz, hankel, plot

Adding, subtracting and multiplying matrices

You can add two matrices of the same dimension (define M and B as above):

>> C=M+B

C =

2 0 0

0 2 0

0 0 2

Example of adding vectors:

>> V=[1 2 3 4] + [0 -1 0 -2]

V =

1 1 3 2

You can multiply matrices if they have compatible dimensions:

>> w=C*[1; 2; 3] % Here C is a 3x3 matrix

w =

The "*" operator is used for matrix multiplication. If you need to multiply two matrices (or vectors) element-by-element use the ".*" operator. For example:

>> [1 2 3 4 5].*[1 -1 1 -1 1]

ans =

1 -2 3 -4 5

Similarly, "./" can be used for element-by-element division of matrices:

>> [1; 2; 3; 4]./[5; 1; 5; 1]

ans =

0.2000

2.0000

0.6000

4.0000

You can always multiply a matrix by a complex scalar:

>> (1-i)*eye(4)

ans =

1.0000 - 1.0000i 0 0 0

0 1.0000 - 1.0000i 0 0

0 0 1.0000 - 1.0000i 0

0 0 0 1.0000 - 1.0000i

For more information, try the following Matlab help topics: mtimes, times, divide, plus, minus

Solving matrix equations

Define

>> A = [1 2 3 4; 2 3 4 1; 3 4 1 2; 4 1 2 3]

A =

1 2 3 4

2 3 4 1

3 4 1 2

4 1 2 3

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

B =

To solve the matrix equation Ax=B, simply type

>> x=A\B

x =

-0.2250

0.0250

0.2750

The symbol "\" denotes left matrix division. To check the solution, type

>> A*x

ans =

In this case, the matrix A is invertible, because it has a non-zero determinant:

>> det(A)

ans =

160

This implies that the equation Ax=B has a unique solution. Matlab uses Gaussian elimination to solve the equation. When A is not a square matrix (so the system is over- or under-determined), Matlab finds a solution (or an approximate solution) according to some least-squares criteria.

Even when A is an invertible square matrix, the system of equations may be numerically unstable. That is, small fluctuations in B could cause significant changes in x. One way to test the numeric stability of a system of equations is to compute the condition number:

>> cond(A)

ans =

5.0000

This number denotes the ratio between the largest and smallest singular values. Large condition numbers indicate numerical instability.

For more information, try the following Matlab help topics: mldivide, mrdivide, det, cond

Other matrix functions

View the "help text" on singular value decomposition.

>> help svd

SVD Singular value decomposition.

[U,S,V] = SVD(X) produces a diagonal matrix S, of the same

dimension as X and with nonnegative diagonal elements in

decreasing order, and unitary matrices U and V so that

X = U*S*V'.

S = SVD(X) returns a vector containing the singular values.

[U,S,V] = SVD(X,0) produces the "economy size"

decomposition. If X is m-by-n with m > n, then only the

first n columns of U are computed and S is n-by-n.

Exercises

(a) Define a row vector of length 100 whose n^th entry is (1-n)/(1+n). Hint: 1. Create vectors of length 100 with entries 0,1,2,3,4,... and 2,3,4,5,... and use the ./ syntax.

(b) (extra credit) Plot the entries of this vector using the Matlab plot command.

2. What matrix is produced by the Matlab expression (1:10)*(1:10).' ? Explain.

3. Compute the condition number of the nxn matrix consisting of all ones except zeros on the diagonal. Do this for various values of n. What is the pattern? Can you give an explanation? (note: each computation requires only a short Matlab expression, try >> help ones).

(a) Use Matlab to compute the condition number of the 8x8 Hilbert matrix (see above or type "help hilb").

(b) What does this large number tell you about solving linear equations involving the Hilbert matrix? Explain.

(c) (extra credit) Use Matlab to find a vector B such that the numerical solution to Ax=B is particularly unstable (here A is the 8x8 Hilbert matrix). Hint: compute the singular value decomposition of A.

Matlab tips and tricks

Always end Matlab statements with a semicolon, unless you want Matlab to display the resulting value of the statement (or expression). You can get the value of any variable by simply typing its name (press enter without a semicolon).
You can use the up and down arrows on your keyboard to retrieve previously executed Matlab statements. For more efficient retrieval, type the first few characters of the statement before using the up arrow.
From the Matlab prompt, type help [function name] to get help on any built-in Matlab function.
It is sometimes tempting to use "for loops" to define and manipulate matrices and arrays. However this should be avoided whenever possible because Matlab is very inefficient when implementing such loops. There are certain Matlab operators and functions which can often be used in place of "for loops". See the following help topics: times, divide, colon, find.
You can combine a sequence of Matlab statements into a single Matlab-executable file called an m-file. See help function.
You can save large data objects in files (for future use). See the help topics: save, load.
You can have more than one statement on a single line: just separate the statements by semicolons.
Have you discovered other useful tips and tricks? Send them to cle@math.upenn.edu .