<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Example 3.1.1
Here are 5-year averages of the worldwide temperature anomaly as compared to the 1951-1980 average (source
Step1:  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.
Step2:  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
Here are the 5-year temperature averages again.
Step3:  The standard best-fit line results from using a linear polynomial that meets the least squares criterion.
Step4:  If we use a global cubic polynomial, the points are fit more closely.
Step5:  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
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  terms of the series above, and $p_k = \sqrt{6s_k}$. Here is a fancy way to compute these sequences in a compact code.
Step6:  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.
Step7:  This suggests a power-law relationship where $\epsilon_k\approx a k^b$, or $\log \epsilon_k \approx b (\log k) + \log a$.
Step8:  In terms of the parameters $a$ and $b$ used above, we have
Step9:  It's tempting to conjecture that $b\to -1$ asymptotically. Here is how the numerical fit compares to the original convergence curve.
Step10:  Example 3.2.1
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.
Step11:  Now we set up an artificial linear least squares problem with a known exact solution that actually makes the residual zero.
Step12:  Using backslash to find the solution, we get a relative error that is about $\kappa$ times machine epsilon.
Step13:  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.
Step14:  Example 3.3.1
We can access both the thin and full forms of the QR factorization. Thin is the default (we use the numpy version here, not the scipy version).
Step15:  Here is a standard call
Step16:  We can test that $\mathbf{Q}$ has orthonormal columns.
Step17:  Here's the full or "complete" factorization.
Step18:  Example 3.4.1
We will use Householder reflections to produce a QR factorization of the matrix
Step19:  Our first step is to introduce zeros below the diagonal in column 1. Define the vector
Step20:  Applying the Householder definitions gives us
Step21:  By design we can use the reflector to get the zero structure we seek
Step22:  Now we let
Step23:  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
Step24:  We now apply the reflector to the submatrix.
Step25:  We need two more iterations of this process.
Step26:  We have now reduced the original  to an upper triangular matrix using four orthogonal Householder reflections
<ASSISTANT_TASK:>
Python Code:
year = arange(1955,2005,5)
y = array([ -0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
    0.1180, 0.2100, 0.3320, 0.3340, 0.4560 ])
fig,ax = subplots()
ax.scatter(year,y,color="k",label="data")
xlabel("year")
ylabel("anomaly (degrees C)")
title("World temperature anomaly");
Explanation: Example 3.1.1
Here are 5-year averages of the worldwide temperature anomaly as compared to the 1951-1980 average (source: NASA).
End of explanation
t = (year-1950)/10
V = vander(t)
c = solve(V,y)
p = poly1d(c)
tt = linspace(1955,2000,500)
ax.plot(tt,p((tt-1950)/10),label="interpolant")
ax.legend();
fig
Explanation: 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.
End of explanation
year = arange(1955,2005,5)
y = array([ -0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
    0.1180, 0.2100, 0.3320, 0.3340, 0.4560 ])
t = (year-1950)/10
Explanation: 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
Here are the 5-year temperature averages again.
End of explanation
V = [ [t[i],1] for i in range(t.size) ]    # Vandermonde-ish matrix
c,res,rank,sv = lstsq(V,y)
fig,ax = subplots()
ax.scatter(year,y,color="k",label="data")
p = poly1d(c)
tt = linspace(1955,2000,500)
ax.plot(tt,p((tt-1950)/10),label="linear fit")
xlabel("year")
ylabel("anomaly (degrees C)")
title("World temperature anomaly");
ax.legend();
Explanation: The standard best-fit line results from using a linear polynomial that meets the least squares criterion.
End of explanation
V = [ [t[i]**3,t[i]**2,t[i],1] for i in range(t.size) ]    # Vandermonde-ish matrix
c,res,rank,sv = lstsq(V,y,rcond=None)
p = poly1d(c)
ax.plot(tt,p((tt-1950)/10),label="cubic fit")
fig
Explanation: If we use a global cubic polynomial, the points are fit more closely.
End of explanation
a = array([1/(k+1)**2 for k in range(100)])
s = cumsum(a)        # cumulative summation
p = sqrt(6*s)
plot(range(100),p,"o")
xlabel("$k$") 
ylabel("$p_k$") 
title("Sequence convergence");
Explanation: 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
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  terms of the series above, and $p_k = \sqrt{6s_k}$. Here is a fancy way to compute these sequences in a compact code.
End of explanation
ep = abs(pi-p)    # error sequence
loglog(range(100),ep,"o")
xlabel("$k$") 
ylabel("error") 
title("Sequence convergence");
Explanation: 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.
End of explanation
V = [ [1,log(k+1)] for k in range(100) ]     # fitting matrix
c = lstsq(V,log(ep),rcond=None)[0]           # coefficients of linear fit
print(c)
Explanation: This suggests a power-law relationship where $\epsilon_k\approx a k^b$, or $\log \epsilon_k \approx b (\log k) + \log a$.
End of explanation
a,b = exp(c[0]),c[1]
print("b:",b)
Explanation: In terms of the parameters $a$ and $b$ used above, we have
End of explanation
loglog(range(100),ep,"o",label="sequence")
k = arange(1,100)
plot(k,a*k**b,"--",label="power fit")
xlabel("$k$");  ylabel("error"); 
legend(); title("Sequence convergence");
Explanation: It's tempting to conjecture that $b\to -1$ asymptotically. Here is how the numerical fit compares to the original convergence curve.
End of explanation
t = linspace(0,3,400)
A = array([ [ sin(t)**2,cos((1+1e-7)*t)**2,1 ] for t in t ])
kappa = cond(A)
print(kappa)
Explanation: Example 3.2.1
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.
End of explanation
x = array([1.,2,1])
b = A@x
Explanation: Now we set up an artificial linear least squares problem with a known exact solution that actually makes the residual zero.
End of explanation
x_BS = lstsq(A,b,rcond=None)[0]
print("observed error:",norm(x_BS-x)/norm(x))
print("max error:",kappa/2**52)
Explanation: Using backslash to find the solution, we get a relative error that is about $\kappa$ times machine epsilon.
End of explanation
N = A.T@A
x_NE = solve(N,A.T@b)
print("observed error:",norm(x_NE-x)/norm(x))
print("digits:",-log10(norm(x_NE-x)/norm(x)))
Explanation: 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.
End of explanation
A = 1.0 + floor(9*rand(6,4))
A.shape
Explanation: Example 3.3.1
We can access both the thin and full forms of the QR factorization. Thin is the default (we use the numpy version here, not the scipy version).
End of explanation
Q,R = qr(A)
print("Q:",Q)
print("R:",R)
Explanation: Here is a standard call:
End of explanation
norm(Q.T@Q - eye(4))
Explanation: We can test that $\mathbf{Q}$ has orthonormal columns.
End of explanation
Q,R = qr(A,"complete")
print(Q.shape)
norm(Q.T@Q - eye(6))
Explanation: Here's the full or "complete" factorization.
End of explanation
A = 1.0 + floor(9*rand(6,4))
m,n = A.shape
Explanation: Example 3.4.1
We will use Householder reflections to produce a QR factorization of the matrix
End of explanation
z = A[:,0]
Explanation: Our first step is to introduce zeros below the diagonal in column 1. Define the vector
End of explanation
v = z - norm(z)*hstack([1,zeros(m-1)])
P = eye(m) - (2/dot(v,v))*outer(v,v)   # reflector
Explanation: Applying the Householder definitions gives us
End of explanation
P@z
Explanation: By design we can use the reflector to get the zero structure we seek:
End of explanation
A = P@A
print(A)
Explanation: Now we let
End of explanation
A[1:,1:]
z = A[1:,1]
v = z - norm(z)*hstack([1,zeros(m-2)])
P = eye(m-1) - (2/dot(v,v))*outer(v,v)
Explanation: 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
End of explanation
A[1:,1:] = P@A[1:,1:]
print(A)
Explanation: We now apply the reflector to the submatrix.
End of explanation
for j in [2,3]:
    z = A[j:,j]
    v = z - norm(z)*hstack([1,zeros(m-j-1)])
    P = eye(m-j) - (2/dot(v,v))*outer(v,v)
    A[j:,j:] = P@A[j:,j:]
Explanation: We need two more iterations of this process.
End of explanation
R = A
print(R)
Explanation: We have now reduced the original  to an upper triangular matrix using four orthogonal Householder reflections:
End of explanation
<END_TASK>