<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 ★ Least Squares ★
Step1:  4.1 Least Squares and the normal equations
$$
A^TA\bar{X} = A^Tb
$$
for the least squares solution $\bar{x}$ that minimizes the Euclidean length if the residual $r = b -Ax$
Example
Solve the least squares problem
$$
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
\begin{bmatrix}
x_1 \ 
x_2
\end{bmatrix}
=
\begin{bmatrix}
-3 \ 
15 \
9
\end{bmatrix}
$$
Step2:  Example
Find the line that best fits the three data points (t,y)=(1,2),(-1,1) and (1,3) <br/>
The model is $y = c_1 + c_2t$
Step3:  The best line is $y = \frac{7}{4} + \frac{3}{4}t$ 
4.3 QR Factorization
Step4:  Example
Find the reduced QR factorization by applying Gram-Schmidt orthogonalization to the columns of 
$
A = 
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
$
Step5:  Example
Find the full QR factorization of 
$
A = 
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
$
Step6:  Example
Use the full QR factorization to solve the least squares problem
$$
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
\begin{bmatrix}
x_1 \ 
x_2
\end{bmatrix}
=
\begin{bmatrix}
-3 \ 
15 \
9
\end{bmatrix}
$$
Step7:  Modified Gram-Schmidt orthogonalization
Step8:  Householder reflector
Example
Let $x = [3,4]$ and $w =[5,0]$. Find a householder reflector H that satisfies $Hx = w$.
Step9:  Example
Use Householder reflectors to find the QR factorization of 
$A = 
\begin{bmatrix}
3 & 1 \
4 & 3 \
\end{bmatrix}
$
Step10:  Example
Use Householder reflectors to find the QR factorization of
$A =
\begin{bmatrix}
1 & -4 \
2 & 3  \
2 & 2  \
\end{bmatrix}
$
Step11:  4.4 Generalized Minimum Residual (GMRES) Method
Example
Solve $Ax = b$ for the following $A$ and $b = [1,0,0]^T$, using GMRES with $x_0 = [0,0,0]^T$
Step12:  4.5 Nonlinear Least Squares
Gauss-Newton Method
Example
Consider the three circles in the plane with centers $(x_1,y_1) = (-1,0)$,$(x_2,y_2) = (1,1/2)$,$(x_3,y_3) = (1,-1/2)$ <br/>
and radii $R_1 = 1,R_2 = 1/2,R3 = 1/2$,respectively. <br/>
Use the Gauss-Newton Method to find the point for which the sum of the squared distances to the three circles is miniminzed.
Step13:  Levenberg-Marquardt Method
Example
Use Levenberg-Marquardt to fit the model $y = c_1e^{-c_2(t - c_3)^2}$ to the data points $(t_i,y_i) = {(1,3),(2,5),(2,7),(3,5),(4,1)}$
<ASSISTANT_TASK:>
Python Code:
# Import modules
import sys
import math
import numpy as np
from matplotlib import pyplot as plt
from scipy import linalg
from scipy import sparse
Explanation: ★ Least Squares ★
End of explanation
A = np.array([1, -4, 2, 3, 2, 2]).reshape(3, 2)
b = np.array([-3, 15, 9])
x = linalg.lstsq(A, b)
print(x[0])
Explanation: 4.1 Least Squares and the normal equations
$$
A^TA\bar{X} = A^Tb
$$
for the least squares solution $\bar{x}$ that minimizes the Euclidean length if the residual $r = b -Ax$
Example
Solve the least squares problem
$$
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
\begin{bmatrix}
x_1 \ 
x_2
\end{bmatrix}
=
\begin{bmatrix}
-3 \ 
15 \
9
\end{bmatrix}
$$
End of explanation
A = np.array([1, 1, 1, -1, 1, 1]).reshape(3, 2)
b = np.array([2, 1, 3])
x = linalg.lstsq(A, b)
print(x[0])
Explanation: Example
Find the line that best fits the three data points (t,y)=(1,2),(-1,1) and (1,3) <br/>
The model is $y = c_1 + c_2t$
End of explanation
def classical_gram_schmidt_orthogonalization(A):
    Q = np.zeros(A.size).reshape(A.shape)
    R = np.zeros(A.shape[1] ** 2).reshape(A.shape[1], A.shape[1])
    for j in range(A.shape[1]):
        y = A[:,j]
        for i in range(j):
            R[i][j] = np.matmul(Q[:,i], A[:,j])
            y = y - R[i][j] * Q[:,i]
        R[j][j] = linalg.norm(y, 2)
        Q[:,j] = y / R[j][j]
    return Q, R
Explanation: The best line is $y = \frac{7}{4} + \frac{3}{4}t$ 
4.3 QR Factorization
End of explanation
A = np.array([1, -4, 2, 3, 2, 2]).reshape(3, 2)
Q, R = classical_gram_schmidt_orthogonalization(A)
print('Q =')
print(Q)
print('R =')
print(R)
Explanation: Example
Find the reduced QR factorization by applying Gram-Schmidt orthogonalization to the columns of 
$
A = 
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
$
End of explanation
A = np.array([1, -4, 2, 3, 2, 2]).reshape(3, 2)
Q, R = linalg.qr(A)
print('Q =')
print(Q)
print('R =')
print(R)
Explanation: Example
Find the full QR factorization of 
$
A = 
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
$
End of explanation
A = np.array([1, -4, 2, 3, 2, 2]).reshape(3, 2)
b = np.array([-3, 15, 9]).T
Q, R = linalg.qr(A)
lu, piv = linalg.lu_factor(R[:2,:])
x = linalg.lu_solve([lu, piv], np.matmul(Q.T, b).reshape(3, 1)[:2])
print('x = %s' %x.T)
Explanation: Example
Use the full QR factorization to solve the least squares problem
$$
\begin{bmatrix}
1 & -4 \ 
2 & 3  \ 
2 & 2
\end{bmatrix}
\begin{bmatrix}
x_1 \ 
x_2
\end{bmatrix}
=
\begin{bmatrix}
-3 \ 
15 \
9
\end{bmatrix}
$$
End of explanation
def modified_gram_schmidt_orthogonalization(A):
    Q = np.zeros(A.size).reshape(A.shape)
    R = np.zeros(A.shape[1] ** 2).reshape(A.shape[1], A.shape[1])
    for j in range(A.shape[1]):
        y = A[:,j]
        for i in range(j):
            R[i][j] = np.matmul(Q[:,i], y)
            y = y - R[i][j] * Q[:,i]
        R[j][j] = linalg.norm(y, 2)
        Q[:,j] = y / R[j][j]
    return Q, R
Explanation: Modified Gram-Schmidt orthogonalization
End of explanation
x = np.array([3, 4]).reshape(2, 1)
w = np.array([5, 0]).reshape(2, 1)
v = w - x
# Projection matrix
P = np.matmul(v, v.T) / np.matmul(v.T, v)
# Householder reflector
H = np.identity(P.shape[0]) - 2 * P
print('H=\n', H)
def householder_reflector(x):
    w = np.zeros(x.size)
    w[0] = linalg.norm(x, 2)
    v = (w - x).reshape(x.size, 1)
    
    # Projection matrix
    P = np.matmul(v, v.T) / np.matmul(v.T, v)
    
    # Householder reflector
    H = np.identity(P.shape[0]) - 2 * P
    
    return H
Explanation: Householder reflector
Example
Let $x = [3,4]$ and $w =[5,0]$. Find a householder reflector H that satisfies $Hx = w$.
End of explanation
A = np.array([3, 1, 4, 3]).reshape(2, 2)
H1 = householder_reflector(A[:,0])
R = np.matmul(H1, A)
Q = H1
print('Q=\n', Q)
print('R=\n', R)
Explanation: Example
Use Householder reflectors to find the QR factorization of 
$A = 
\begin{bmatrix}
3 & 1 \
4 & 3 \
\end{bmatrix}
$
End of explanation
A = np.array([1, -4, 2, 3, 2, 2]).reshape(3, 2)
H1 = householder_reflector(A[:, 0])
TEMP = np.matmul(H1, A)
H2 = householder_reflector(TEMP[1:, 1])
H2_Ext = np.identity(H1.shape[0])
H2_Ext[-H2_TMP.shape[0]:, -H2_TMP.shape[1]:] = H2
R = np.matmul(np.matmul(H2_Ext, H1), A)
Q = np.matmul(H1, H2_Ext)
print('Q=\n', Q)
print('R=\n', R)
Explanation: Example
Use Householder reflectors to find the QR factorization of
$A =
\begin{bmatrix}
1 & -4 \
2 & 3  \
2 & 2  \
\end{bmatrix}
$
End of explanation
A = np.array([1, 1, 0, 0, 1, 0, 1, 1, 1]).reshape(3, 3)
b = np.array([1, 0, 0]).reshape(3, 1)
x0 = np.zeros(3).reshape(3, 1)
x, info = sparse.linalg.gmres(A, b, x0)
print('x = %s' %x)
A = np.arange(1, 10).reshape(3, 3)
D = np.diag(A.diagonal())
print(D)
print(linalg.inv(D))
Explanation: 4.4 Generalized Minimum Residual (GMRES) Method
Example
Solve $Ax = b$ for the following $A$ and $b = [1,0,0]^T$, using GMRES with $x_0 = [0,0,0]^T$
End of explanation
def R_xy(x, y):
    A = np.zeros(3)
    f = lambda xf, yf, R : math.sqrt(pow(x - xf, 2) + pow(y - yf, 2)) - R
    A[0] = f(-1,   0, 1)
    A[1] = f( 1, 0.5, 0.5)
    A[2] = f( 1,-0.5, 0.5)
    return A
def DR_xy(x, y):
    A = np.zeros(6).reshape(3, 2)
    fx = lambda xf, yf : (x - xf) / math.sqrt(pow(x - xf, 2) + pow(y - yf, 2))
    fy = lambda xf, yf : (y - yf) / math.sqrt(pow(x - xf, 2) + pow(y - yf, 2))
    A[0][0] = fx(-1, 0)
    A[0][1] = fy(-1, 0)
    A[1][0] = fx(1, 0.5)
    A[1][1] = fy(1, 0.5)
    A[2][0] = fx(1, -0.5)
    A[2][1] = fy(1, -0.5)
    return A
def gauss_newton_method(x0, y0, k):
    
    xk = np.array([x0, y0])
    
    for _ in range(k):
        x = xk[0]
        y = xk[1]
        A = DR_xy(x, y)
        r = R_xy(x, y)
        v = np.matmul(linalg.inv(np.matmul(A.T, A)), -np.matmul(A.T, r))
        xk = xk + v
        
    return xk
x = gauss_newton_method(0, 0, 8)
print('x = %s' %x)
Explanation: 4.5 Nonlinear Least Squares
Gauss-Newton Method
Example
Consider the three circles in the plane with centers $(x_1,y_1) = (-1,0)$,$(x_2,y_2) = (1,1/2)$,$(x_3,y_3) = (1,-1/2)$ <br/>
and radii $R_1 = 1,R_2 = 1/2,R3 = 1/2$,respectively. <br/>
Use the Gauss-Newton Method to find the point for which the sum of the squared distances to the three circles is miniminzed.
End of explanation
def R_xy(c):
    c1 = c[0]
    c2 = c[1]
    c3 = c[2]
    
    r = np.zeros(5)
    f = lambda t, y : c1 * math.exp( -c2 * pow(t - c3, 2) ) - y
    
    r[0] = f(1, 3)
    r[1] = f(2, 5)
    r[2] = f(2, 7)
    r[3] = f(3, 5)
    r[4] = f(4, 1)
    
    return r
def DR_xy(data, c):
    c1 = c[0]
    c2 = c[1]
    c3 = c[2]
    
    DR = np.zeros(15).reshape(5, 3)
    
    f0 = lambda t : math.exp( -c2 * pow(t - c3, 2) )
    f1 = lambda t : -c1 * pow(t - c3, 2) * math.exp( -c2 * pow(t - c3, 2) )
    f2 = lambda t : 2 * c1 * c2 * (t - c3) * math.exp( -c2 * pow(t - c3, 2) )
    
    for i in range(5):
        t = data[i][0]
        DR[i][0] = f0(t)
        DR[i][1] = f1(t)
        DR[i][2] = f2(t)
    
    return DR
def levenberg_marquardt_method(data, c, la, k):
    ck = c
    for _ in range(k):
        A = DR_xy(data, ck)
        r = R_xy(ck)
        mAr = -np.matmul(A.T, r)
        invA = np.linalg.inv(np.matmul(A.T, A) + la * np.diag(np.matmul(A.T, A).diagonal()))
        v = np.matmul(invA, mAr)
        ck = ck + v
    return ck
data = np.array([(1, 3), (2, 5), (2, 7), (3, 5), (4, 1)])
c = np.array([1, 1, 1])
c = levenberg_marquardt_method(data, c, 50, 1200)
f = lambda t, c1, c2, c3 : c1 * np.exp( -c2 * np.power(t - c3, 2) )
X = np.linspace(0, 5, 100)
Y = f(X, *c)
plt.plot(X, Y, color='cyan')
plt.plot(data[:,0], data[:,1], linestyle='', markersize=8, marker='.', color='blue')
plt.show()
Explanation: Levenberg-Marquardt Method
Example
Use Levenberg-Marquardt to fit the model $y = c_1e^{-c_2(t - c_3)^2}$ to the data points $(t_i,y_i) = {(1,3),(2,5),(2,7),(3,5),(4,1)}$
End of explanation
<END_TASK>