% reviewtut.m
% A review tutorial for the first lab session of Inverse Theory ESS523,
% to get students brushed up on both Matlab usage and linear algebra.
% Includes material from Appendix A of Aster, Borchers, & Thurber,
% chapter 1 of Gelb, and the very dredges of Andy's mind...
% Andy Ganse, 1 Oct 2005
% ---------------------------------------------------------------------
% JUST START READING AND FOLLOW THE INSTRUCTIONS AS YOU GO...
% ---------------------------------------------------------------------
% Okay all you inverse theory fans out there!
% Let's get started with the ultra ultra basics of Matlab:
% a.) Notice that these Matlab "script" files always end in ".m"
% b.) You can run them in Matlab by entering the filename without the ".m"
% at the Matlab command prompt. Alternatively you can open them in
% the editor of the GUI version of Matlab and click the "Run" button.
% c.) If you want to "comment out" a line in the script, i.e. use it for
% writing rather than for code that gets run, put a "%" as the first
% character. As you'll see below, Matlab actually considers anything
% to the right of a "%" as a comment (note that includes the beginning
% of the line).
%%% %% By the way, when "commenting out" a block of code that already has
%%%% % comments in it you may end up with several %s in there - no problem.
% d.) There are several help mechanisms:
% - to get info about a single function like cos() for example, enter
% "help cos" at the Matlab prompt.
% - to get a list of commands which relate to a subject like
% convolution for example), enter "lookfor convolution" at the
% Matlab prompt.
% - online there is general documentation at:
% http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.shtml
% and a users' guide at:
% http://www.mathworks.com/access/helpdesk/help/techdoc/math/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% As you go thru the tutorial in the rest of this file, you will "comment
% out" and "uncomment out" various lines, and change some of the variable
% values, and then rerun the script to note the changes in computations and
% output. You can do this anywhere you like, but the lines labeled "HEY!"
% in the file point out particular places to do so, to emphasize key
% points. Deciphering which lines are code to uncomment, and which are my
% notes and explanations, is perhaps part of the learning process - but
% please don't hesitate to ask for clarification and help!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SETTING VARIABLES
% -----------------
a=2.718 % stored as a double precision floating point number
b=4e-6 % ditto
c='yeehaw' % this is a "string", ie a word rather than number.
% A string could actually contain numerical digits but it can't
% be used in a numerical calculation (there are a few weird
% exceptions to that but please don't worry about those for now)
% Notice that without a semicolon at the end of the statement, Matlab echos
% the values at the command prompt.
% When we do vectors and matrices below, placing that semicolon at the end
% of the line is crucial, so that (for example) a vector of length 10,000
% doesn't take 5 minutes to scroll down your screen.
%% HEY!: Put semicolons after the statements above now and rerun the
%% script to compare what happens.
% A FEW STANDARD MATHEMATICAL OPERATORS AND FUNCTIONS
% --------------------------------------------------
%% HEY!: keep the above lines specifying a, b, & c uncommented, and now
%% uncomment the lines below:
% d=a+b
% e=3*a+4/b
% disp('howdy'); % a handy function for outputting one-line notes
% f=strcat(c,' man') % note there are bunch of these "str___" commands
% % for manipulating strings of characters
% g=4*log(b)+sin(pi*a) % note: pi is built-in, and sin() expects radians,
% % and log() is ln... for base10 use log10().
% disp('Note the following produces an error: why?');
% stringnumber='1.5708'
% h=sin(stringnumber) % why does this line produce an error?
%% HEY!: what can you do to the previous two lines to fix the error?
%% Lastly, remember you can always look up the use of other math functions
%% with "help" & "lookfor"
% COMPLEX NUMBERS
% ---------------
% Matlab accepts both "i" and "j" as the imaginary number.
%% HEY!: now comment out everything above, and uncomment the lines below
% a=1.3+4.5i;
% b=6.7+2.1j;
% c=a+b
% phase=angle(c)
% amplitude=abs(c)
% foo=real(c)
% boo=imag(c)
% VECTORS AND MATRICES
% (when starting this section, comment out all above, and gradually
% uncomment lines a few at a time in this section, running each time)
% -------------------------------------------------------------------
% Matlab will both calculate and display row and column vectors in the
% correct directions.
% rowvector=[1,2,3,4,5] % or also simply [1 2 3 4 5]
% colvector=[1;2;3;4;5]
% matrix=[1,2,3; 2,4,9; 7,11,13]
% Recall what a matrix transpose is. Matlab uses the ' mark as a
% transpose operator.
%% HEY!: note the output of the transposes below to remind yourself what
%% a transpose is (uncomment the next 3 lines):
% rowvector'
% colvector'
% matrix'
%
% We can also create vectors by specifying a range:
% rowvector=1:5
% colvector=(1:5)'
% rowvector2=2:.1:2.6 % you can specify a step increment too.
% colvectorlong=(1:1000)'; % what happens if you remove the semicolon here?
% We can access just one element of a vector or matrix if needed:
% p=rowvector(1)
% q=colvector(3)
% r=matrix(2,3) % (first arg here is row, 2nd arg is column)
% We can specify "all" of a row or column like this:
% s=matrix(2,:) % 2nd row, all columns
% Recall we can multiply vectors and matrices, but unlike scalars, here
% order is important
% innerproduct=rowvector*colvector % also called dot product
% outerproduct=colvector*rowvector
%% HEY!: compare the innerproduct and outerproduct of vectors - which is
%% a scalar and which is a matrix? What are the dimensions of the matrix?
%% What is the geometric interpretation of a dot product of two vectors?
%
% Most of the other mathematical operators and functions work on vectors &
% matrices too, just with some twists:
% rowvector+rowvector
% rowvector+2
% rowvector*2
% rowvector.*rowvector % the "." on * and / and ^ means element-by-element
% disp('the following will cause an error');
% rowvector*rowvector % HEY!: why does this cause an error?
% rowvector.^2
% matrix*2
% matrix.^2
% matrix^2
%% HEY!: closely investigate the output of the previous two lines to
%% notice the different between using ".^" and "^" for the matrix.
%% What is matrix^0 ? And matrix.^0 ?
%
% We can also multiply a vector and a matrix. This is the foundation
% for a linear problem such as Ax=b:
% x=[1,2,3]'
% A=matrix % let's reuse the earlier "matrix". (is it uncommented?)
% b=A*x
% The above represents 3 equations and 3 variables in a linear problem.
%% HEY!: what are the 3 equations? is the solution unique or could there
%% be more than one to this problem?
%
% To take the inverse (if it exists) of a matrix, use the inv() function.
% If the inverse doesn't exist or if the matrix is even "almost" singular
% Matlab will tell you in a warning message.
% Ainv=inv(A)
% So if you have a linear problem like the above, and as long as the matrix
% A is "nice", ie invertible (square, full rank, and so on), you could
% solve for x like this if you didn't have it:
% xsolve=inv(A)*b
%% HEY!: Now let's mess up the matrix A and see what happens. Change its
%% 2nd row to equal the first and try running it again:
% A(2,:)=A(1,:)
% b=A*x
% xsolve=inv(A)*b
%% HEY!: So what happened? Why'd the computation mess up this time?
%
% Now say we had four equations and three variables. This is an over-
% determined problem; we have more pieces of information than variables
% to solve for. In this case our matrix A has more rows than columns.
% There may not be a solution that exactly satisfies all the equations so
% we find a least-squares (LS) solution.
% x=[1,2,3]';
% A=[1,2,3; 2,4,9; 7,11,13; 8,14,2]; % 4 rows, 3 cols
% b=A*x+[.1,-.1,.1,-.1]'; % (sortof artificially making up the problem here!)
% We can use our class material to solve this:
% xbest=inv(A'*A)*A'*b
% HEY!: how close is xbest to your original x?
%
% HEY!: Now what if we had three equations and four variables? This is
% an underdetermined problem, where you don't have enough information to
% uniquely solve the problem, so you have an infinite number of solutions.
% An example might be "fit the best plane to two given points".
%
% Anyways, many math functions work as usual on vectors and matrices, just
% repeated on each element:
% sin(rowvector)
% exp(matrix)
% log(5-matrix)
%% HEY!: where do the complex numbers come from in that matrix output?
%
% Matlab has a function to create an identity matrix of whatever size:
% I=eye(4)
% And two functions to create matrices of ones or zeros:
% M=ones(5,3)
% N=zeros(2,4)
% A function to create a vector or matrix of uniformly-dist random numbers:
% P=rand(1000,1); % this is a col vector, row vector would have been (1,50)
% And a vector or matrix of normally-dist (Gaussian) random numbers:
% Q=randn(1000,1); % this is a col vector, row vector would have been (1,50)
% We'll compare the values in P and Q in a few moments when we learn how
% plot a histogram. Note the semicolons there as these vectors are so
% long, in order to make sense of their histograms. Ok, back to this soon.
%
%% HEY!: now to remind yourself of some linear algebra identities:
% let's load up the "good" matrix again (the one we didn't mess up):
% A=matrix;
%
% samematrix=eye(3)*A
% the following two results should be the same:
% Ainv'
% inv(A')
% the following two results should be the same:
% Ainv*A % is the result exactly an identity matrix? or at least close?
% A*Ainv % why might it be a teeny tiny bit off?
% Let's bring in another matrix:
% B=[2,8,11; 1,7,12; 6,10,4]
% the following two results should be the same:
% (A*B)'
% B'*A'
% the following two results should be the same:
% inv(A*B)
% inv(B)*inv(A)
% HEY!: (does the order matter on those multiplications?)
% PLOTTING - FINALLY!
% See, we had to learn how to create and manipulate vectors and matrices
% before we could make sense of this section...
% ------------------------------------------------------------
% (You can comment out everything in earlier sections now, except for the
% P & Q random vectors - we'll finally plot those histograms.)
%
% myXvector=1:50;
% myYvector=sin(.2*myXvector);
% plot(myXvector,myYvector); % this will open a figure window
% grid on; % put a grid on the plot
% let's open another figure window and show how to fancy up a plot:
% figure;
% plot(myXvector,myYvector,'r--',...
% myXvector,myYvector.*myXvector/25,'ks',...
% 'LineWidth',2,...
% 'MarkerEdgeColor','k',...
% 'MarkerFaceColor','g',...
% 'MarkerSize',10);
% title('This is my Plot','FontWeight','bold','FontSize',12);
% xlabel('FooFoos');
% ylabel('BooBoos');
% Check "help plot" for some info on plot options, and if you want to
% really get into the rest, learn about "handle graphics". You can check
% and set every option for a plot via its "handle" which can be optionally
% returned like this:
% h=plot(myXvector,myYvector);
% set(h) % tells you all the options available for that plot
% get(h) % tells you the current values of all properties of that plot
% set(h,'LineWidth',1,'Marker','o'), % changes these two properties
%
% And here's a handy way to plot matrix contents:
% figure;
% C=randn(30); % 30x30 matrix of Gaussian values
% imagesc(C); % plot with colors scaled between max & min elements
% colorbar; % add a colorbar legend to the plot
% a handy option for this type of plot, to move the origin to the bottom
% left (default is top left).
% axis xy
% A quick digression, simply to make a cool-looking plot(!). One way to
% smooth an image or surface like that random matrix is to filter it by
% convolving it with a smoothing kernel. Matlab has a built-in convolution
% function:
% x=linspace(-1,1,15)'; % 15 points from -1 to 1 for the bell-curve below:
% g=exp(-x.^2/2/(.4)^2)/sqrt(2*pi)/.4; % points on a Gaussian bell-curve
% kernel=g*g'; % HEY!: what is the dimension of kernel here: matrix? scalar?
% D=conv2(C,kernel);
% imagesc(D);
% colorbar;
% HEY!: try changing the number of points in the bell-curve above - ie the
% size of the smoothing kernel, and look at the effect in the rerun plot.
% How's it look with 5 points, 10 points, 30 points?
%
% Ok, let's finally look at those histograms for P & Q, the uniformly-
% distributed and normally-distributed (respecitively) random vectors we
% created earlier. (Are they still uncommented?)
% Matlab has a handy histogram plotting function that does it all for you:
% hist(P,20); % specifying 20 bins to histogram over. default is 10.
% figure; % start a 2nd plot
% hist(Q,20); % specifying 20 bins to histogram over. default is 10.
% HEY!: okay, see the difference between being uniformly distributed
% between 0 & 1, and being normally distributed about 0 with stdev 1?
% Try this: shrink and expand the 20 bins to say 4 and to 200. How's that
% affect the histogram plots, and your interpretation of the distributions?
%
%
%
%
%
%
%
%
% YOUR HOMEWORK FOR THIS LAB:
% -------------------------------
% Use the example bits in this file to piece together a script of your
% own which demonstrates the gravity problem discussed in class and the
% book, where the ball's elevation over time is a parabola. (This is
% example 1.1 in the Aster book.)
% Do as above in this file, where we made up some x ("m" in the book, ie
% the initial position, initial velocity, and acceleration of gravity).
% Create a vector of times t_i, say maybe 1:.1:10 (or whatever range).
% Create an A matrix ("G" in the book) with the 3 columns based on the t_i.
% Calculate the trajectory y_i of elevations over times t_i, a parabola.
% Now add some Gaussian noise to y using the randn() function (play with
% how big its standard deviation is).
% Use the normal equations to find the LS soln "mhat" for the initial pos,
% initial vel, and gravitational accel.
% Calculate the trajectory yhat of the best fit to the data, from mhat.
%
% Print out and turn in a plot showing the noisy data y and the best fit
% yhat to the data, on top of each other. And show your original m and
% your mhat as well for comparison.
%
% Don't hesitate to holler if any of this is not clear or if you're having
% Matlab trouble...
%
% Cheers!
% -Andy