%% Chapter 3
% This is a MATLAB Live Editor version of the Fundamentals of Numerical Computation Jupyter Notebooks found at 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,'.','MarkerSize',20)
xlabel('year'), 
ylabel('anomaly ({\circ}C)'), 
axis tight   
set(gca,'fontsize',14)

%
% 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])
set(gca,'fontsize',14)

% 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   
set(gca,'fontsize',14)
hold off

%
% 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])

hold on 
plot(year,y,'.','MarkerSize',20)
xlabel('year'), 
ylabel('anomaly ({\circ}C)'), 
axis tight   
set(gca,'fontsize',14)
hold off

% 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
hold off

%
% 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