%% Chapter 3
% _This is a MATLAB Live Editor version of the *Fundamentals of Numerical Computation* 
% Jupyter Notebooks found at_ <https://tobydriscoll.net/project/fnc/ _https://tobydriscoll.net/project/fnc/_>_._ 
% Example 3.1.1 Worldwide Temperature
% Here are 5-year averages of the worldwide temperature anomaly as compared 
% to the 1951-1980 average (source: NASA).

clear
t = (1955:5:2000)';
y = [ -0.0480; -0.0180; -0.0360; -0.0120; -0.0040;
    0.1180; 0.2100; 0.3320; 0.3340; 0.4560 ];
plot(t,y,'.')
xlabel('year'), ylabel('anomaly ({\circ}C)'), axis tight   % ignore this line
%% 
% A polynomial interpolant can be used to fit the data. Here we build one using 
% a Vandermonde matrix. First, though, we express time as decades since 1950, 
% as it improves the condition number of the matrix. 

t = (t-1950)/10;  
V = t.^0;      % vector of ones
n = length(t);
for j = 1:n-1
    V(:,j+1) = t.*V(:,j);
end
c = V\y;
p = @(x) polyval(flipud(c),(x-1950)/10);
hold on, fplot(p,[1955 2000])
%% 
% As you can see, the interpolant does represent the data, in a sense. However 
% it's a crazy-looking curve for the application. Trying too hard to reproduce 
% all the data exactly is known as _overfitting_.
% 
% *Example 3.1.2 Temperature Fit*
% 
% Here are the 5-year temperature averages again.

year = (1955:5:2000)';
y = [ -0.0480; -0.0180; -0.0360; -0.0120; -0.0040;
    0.1180; 0.2100; 0.3320; 0.3340; 0.4560 ];
%% 
% The standard best-fit line results from using a linear polynomial that meets 
% the least squares criterion. 

t = year - 1955;             % better matrix conditioning later
V = [ t.^0 t ];    % Vandermonde-ish matrix
c = V\y;
f = @(x) polyval(c(end:-1:1),x-1955);
fplot(f,[1955 2000])
hold on, plot(year,y,'.')
xlabel('year'), ylabel('anomaly ({\circ}C)'), axis tight   % ignore this line
%% 
% If we use a global cubic polynomial, the points are fit more closely.

V = [ t.^0 t t.^2 t.^3];    % Vandermonde-ish matrix
c = V\y;   f = @(x) polyval(c(end:-1:1),x-1955);
fplot(f,[1955 2000])
%% 
% If we were to continue increasing the degree of the polynomial, the residual 
% at the data points would get smaller, but overfitting would increase.
% 
% *Example 3.1.3 Pi Rate*
% 
% Finding numerical approximations to $\pi$ has fascinated people for millenia. 
% One famous formula is
% 
% $$ \frac{\pi^2}{6} = 1 + \frac{1}{2^2} + \frac{1}{3^2} + \cdots.$$
% 
% Say $s_k$ is the sum of the first $k$ terms of the series above, and $p_k 
% = \sqrt{6s_k}$. Here is a fancy way to compute these sequences in a compact 
% code.

k = (1:100)';
s = cumsum( 1./k.^2 );   % cumulative summation
p = sqrt(6*s);
plot(k,p,'.-')
xlabel('k'), ylabel('p_k')
title('Sequence converging to \pi')
%% 
% This graph suggests that $p_k\to\pi$ but doesn't give much information about 
% the rate of convergence. Let $\epsilon_k = |\pi - p_k|$ be the sequence of errors. 
% By plotting the error sequence on a log-log scale, we can see a nearly linear 
% relationship.

ep = abs(pi-p);    % error sequence
loglog(k,ep,'.'), title('log-log convergence')
xlabel('k'), ylabel('error'),  axis tight
%% 
% This suggests a power-law relationship where $\epsilon_k\approx a k^b$, or 
% $\log \epsilon_k \approx b (\log k) + \log a$.

V = [ k.^0, log(k) ];    % fitting matrix
c = V \ log(ep)          % coefficients of linear fit
%% 
% In terms of the parameters $a$ and $b$ used above, we have 

a = exp(c(1)),  b = c(2),
%% 
% It's tempting to conjecture that $b=-1$ asymptotically. Here is how the numerical 
% fit compares to the original convergence curve. 

hold on, loglog(k,a*k.^b,'r'), title('power-law fit')
axis tight    % ignore this line
%% 
% *Example 3.2.1 Normal Instability*
% 
% Because the functions $\sin^2(t)$, $\cos^2(t)$, and $1$ are linearly dependent, 
% we should find that the following matrix is somewhat ill-conditioned.

t = linspace(0,3,400)';
A = [ sin(t).^2, cos((1+1e-7)*t).^2, t.^0 ];
kappa = cond(A)
%% 
% Now we set up an artificial linear least squares problem with a known exact 
% solution that actually makes the residual zero.

x = [1;2;1];
b = A*x;
%% 
% Using backslash to find the solution, we get a relative error that is about 
% $\kappa$ times machine epsilon.

x_BS = A\b;
observed_err = norm(x_BS-x)/norm(x)
max_err = kappa*eps
%% 
% If we formulate and solve via the normal equations, we get a much larger relative 
% error. With $\kappa^2\approx 10^{14}$, we may not be left with more than about 
% 2 accurate digits.

N = A'*A;
x_NE = N\(A'*b);
observed_err = norm(x_NE-x)/norm(x)
digits = -log10(observed_err)
%% 
% *Example 3.3.1 QR Factorization*
% 
% MATLAB gives direct access to both the full and thin forms of the QR factorization. 

A = magic(5);
A = A(:,1:4);
[m,n] = size(A)
%% 
% Here is the full form:

[Q,R] = qr(A);
szQ = size(Q)
szR = size(R)
%% 
% We can test that $\mathbf{Q}$ is orthogonal.

QTQ = Q'*Q
norm(QTQ - eye(m))
%% 
% With a second input argument given, the thin form is returned.

[Q,R] = qr(A,0);
szQ = size(Q)
szR = size(R)
%% 
% Now $\mathbf{Q}$ cannot be an orthogonal matrix, because it is not even square, 
% but it is still ONC.

Q'*Q - eye(4)
%% 
% *Example 3.4.1 Householder QR Factorization*
% 
% We will use Householder reflections to produce a QR factorization of the matrix

A = magic(6);
A = A(:,1:4);
[m,n] = size(A);
%% 
% Our first step is to introduce zeros below the diagonal in column 1. Define 
% the vector 

z = A(:,1);
%% 
% Applying the Householder definitions gives us

v = z - norm(z)*eye(m,1);
P = eye(m) - 2/(v'*v)*(v*v');   % reflector
%% 
% By design we can use the reflector to get the zero structure we seek:

P*z
%% 
% Now we let 

A = P*A
%% 
% We are set to put zeros into column 2. We must not use row 1 in any way, lest 
% it destroy the zeros we just introduced. To put it another way, we can repeat 
% the process we just did on the smaller submatrix

A(2:m,2:n)
%%
z = A(2:m,2);
v = z - norm(z)*eye(m-1,1);
P = eye(m-1) - 2/(v'*v)*(v*v');   % reflector
%% 
% We now apply the reflector to the submatrix.

A(2:m,2:n) = P*A(2:m,2:n)
%% 
% We need two more iterations of this process.

for j = 3:n
    z = A(j:m,j);
    k = m-j+1;
    v = z - norm(z)*eye(k,1);
    P = eye(k) - 2/(v'*v)*(v*v');
    A(j:m,j:n) = P*A(j:m,j:n);
end
%% 
% We have now reduced the original $\mathbf{A}$ to an upper triangular matrix 
% using four orthogonal Householder reflections:

R = A