% Script 3 for Week 3

 
 Vector Spaces (such as R^n)
 
Every vector space has a basis: a set of vectors so that ANY vector can be written in exactly one way as a linear combination of these basis vectors. And any basis of a given vector space has the same "cardinality"--that is, the same number of vectors in it.


Solving Matrix Equations

We want to solve H*X=C where 
H is an mxn matrix 
X is a "variable stack" of needed values in R^n
C is a "constant stack" in R^m

You can think of left multiplication by H as a linear function from R^n to R^m. So for a "target" vector C in R^m we are trying to find X in R^n for which H*X "hits the target." That is, H*X=C.

If there are exact solutions then C must be in the columnspace of H, since H*X is explicitly a linear combination of the columns of H.

And if there is more than one solution then any two solutions differ by a member of the nullspace of H, the subspace of vectors "killed" by H.

And we know that the dimension of the range, i.e. the potential "exactly solvable targets," plus the dimension of the nullspace is the dimension of the domain space: in this case R^n.

H=[4 5 6 9 8 6 9
8 0 4 0 8 1 3]

C=[7;3]

% The columnspace is obviously two dimensional so here we can hit any target we want.

y=H\C

H*y

R=null(H)

% The columns of R form an orthonormal basis of the nullspace of H.

R(:,4)'*R(:,4)

R(:,2)'*R(:,4)

H*R(:,4)

R

% You will note that there are 5 of these columns: the nullspace has dimension 5.

% Any solution X (or Least-Squares Best Approximate Solution if there is no ACTUAL solution) is of the form

%      X = y + generic member of nullspace
      
y+R*[0;33;0;72;0]

H*(y+72*R(:,4)+33*R(:,2))
 
 
% The columns of 
 
 S=rref(R')'
 
% have the same span as the columns of R, but may be easier to work with in some ways.

H*(y+S*[-7;-30;1;22;5])
 
% Back to solutions of systems
 
 A=[2 -3;
 4 -6]
 
 B=[1;3]
 
rref([A B])

% No exact solution!!
  
% Any solution to AX=B must be in the columnspace of A.
 
% A basis for the columnspace is those columns of A corresponding to the pivot columns of rref(A).

% Note: Here the situation is obvious, but in a more general situation, a (possibly) nicer basis of the columnspace is the set of nonzero columns of 

(rref(A'))'.
 
% Anyway, B is not in the columnspace of A so you cannot hit B. There is no solution.

% But you can always find a solution to  A'*A*X=A'*B 

% A'B is ALWAYS in the column space of A'*A.

% This is called the "normal equation" produced from the original equation.

rref([A'*A  A'*B])
 
 A\B
 
% gives an error because A is square and singular. 
 
 g=(A'*A)\(A'*B)
 
 % g gives one particular least-squares solution with a singularity warning.
 
% Better in some ways is: 
 
 h=pinv(A)*B

% which gives another least-squares solution, the solution with smallest norm.
 
% This is one point in the domain which is taken to the exact place [1.4;2.8]

A*g
A*h

% which is the point in the range space closest to the target B=[1;3].

x= lsqminnorm(A,B)

% gives another least squares solution. Use this if you are going to be solving A*x=B for just one B, because it is cheaper to calculate than pinv(A)*B.

% However if you are going to be solving the same equation for many different B targets you should calculate D=pinv(A) once and find your solution as D*B for each B.

% In any case the general least-squares solution is given by any particular solution plus an arbitrary member of the nullspace.

%    g + general member of the nullspace of A
 
null(A)

% This nullspace has dimension 1 so the general solution is

%    g + t*[-0.8321;-0.5547] for arbitrary t

% More generally the solution is 

%    g + null(A)*t 

% where t is a vector of coefficients for the null vectors.

% The times below for a 1000x1000 matrix suggest which methods might be preferred.

A = rand(1000,1000);
B = rand(1000,1);
D = pinv(A);

tic
 g1=A\B;
toc

tic
 g2=(A'*A)\(A'*B);
toc

tic
 g3=inv(A'*A)*(A'*B);
toc

tic
x= lsqminnorm(A,B);
toc

tic
h1=pinv(A)*B;
toc

tic
h2=D*B;
toc

[norm(g1-g2) norm(g1-g3) norm(g1-x) norm(g1-h1) norm(g1-h2)]


% The error in g2 and g3 is caused by roundoff issues in the MANY calculations done to calculate the matrix products and inverse.
 
% If you are going to be doing this for many B targets, predefining D as pinv(A) (the last option) and multiplying D*B is the winner but A\B is close.