% A first example of least squares approximation and
% projections in a function setting. One key objective
% is to casually identify functions with long column
% vectors of sample points.
%
% The number of data points is chosen to just fill the
% screen from top to bottom, but still be visible.
% The inner product does not perfectly reflect a Riemann
% sum of the integral, the standrad inner product as it
% weighs both endpoints the same as every interior point.
% Still, the orthogonality of the chosen trig functions
% unsuspectingly pops up in the unexpected orthogonality
% of A when calculating A'*A.
% Immediate improvements for higher order approximations
% come up.
% Very nicely, c(2)=0 is readily explained geometrically
% from the symmetry properties.
% Caution: for periodic extensions either consider the
% odd extension, or include a constant c0*ones(size(xx)).
%
% July 2008, M. Kawski, http://math.asu.edu/~kawski
%
xx=[0:.1:2]';
yy=1-abs(xx-1);
f1=sin(pi/2*xx);
f2=sin(2*pi/2*xx);
f3=sin(3*pi/2*xx);
[xx,yy,f1,f2,f3]
A=[f1,f2,f3];
% Just rename to connect with notation of the 3D picture.
b=yy;
% Verify the sizes of the matrices. Invertibility.
A'*A
A'*b
% Recall that A\b is three times faster than inv(A)*b.
c=A'*A\A'*b
% Here parenntheses are not needed, (but recommended)
% as MATLAB reads left 2 rite
[A'*A\A'*b,(A'*A)\A'*b]
% Calculate the projection from c in two separate steps.
p=A*c
% But be careful when doing it all in one step in MATLAB.
% In the next command parentheses are essential!
p=A*((A'*A)\A'*b)
% Enjoy the built-in MATLAB notation for pseudo inverse
[A\b,(A'*A)\A'*b]
[xx,yy,p,f1,f2,f3]
plot(xx,yy,'-b',xx,yy,'ob',xx,p,'-r',...
     xx,c(1)*f1,'-g',xx,c(2)*f2,'-g',xx,c(3)*f3,'-g',...
     xx,f1,'-c',xx,f2,'-c',xx,f3,'-c')
axis([-.2,2.2,-1.2,1.2])