数值分析作业-拟合和最小二乘法

Homework

• 用周期模型 $F_3(t)=c_1+c_2\cos 2\pi t+c_3\sin 2\pi t$ 拟合以下数据:$(0,0),\left(\frac{1}{6},2\right),\left(\frac{1}{3},0\right),\left(\frac{1}{2},-1\right),\left(\frac{2}{3};1\right),\left(\frac{5}{6},1\right)$
• 分别使用法线方程和 QR 分解求解方程的最小二乘解, 并求出二范数误差 $\left[\begin{array}{cc}-4& -4\\ -2& 7\\ 4& -5\end{array}\right]\left[\begin{array}{l}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{l}3\\ 9\\ 0\end{array}\right]$
• 编程: 令 $\left\{x_0,\dots ,x_{20}\right\}$ 为区间 $[2,4]$ 上间隔均匀的点, 并设 $y_i={\sum }_{n=0}^{10}{x}_i^n(0\le i\le 20)$ . 使用法线方程和 QR 分解找出拟合 ${\left\{\left(x_i,y_i\right)\right\}}_{i=0}^{20}$ 的最小二乘多项式 $P(x)={\sum }_{n=0}^{10}{c}_n{x}^n$

$$\left(\begin{array}{ccc}1& 1& 0\\ 1& 1/2& \frac{\sqrt{3}}{2}\\ 1& -1/2& \frac{\sqrt{3}}{2}\\ 1& -1& 0\\ 1& -1/2& \frac{-\sqrt{3}}{2}\\ 1& 1/2& \frac{-\sqrt{3}}{2}\end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right)=\left(\begin{array}{c}0\\ 2\\ 0\\ -1\\ 1\\ 1\end{array}\right)$$ $$\begin{array}{rl}{A}^{\top }AC& =A^{\top }b\\ \left(\begin{array}{ccc}6& 0& 0\\ 0& 3& 0\\ 0& 0& 3\end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right)& =\left(\begin{array}{c}3\\ 2\\ 0\end{array}\right)\\ \left(\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right)& =\left(\begin{array}{c}1/2\\ 2/3\\ 0\end{array}\right)\end{array}$$

$F_3(t)=\frac{1}{2}+\frac{2}{3}cos2\pi t$

---

$$\left(\begin{array}{cc}36& -18\\ -18& 90\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{c}-30\\ 51\end{array}\right)$$

解得: $x_1=\frac{-11}{18},x_2=\frac{4}{9}$ , 误差： $r=b-Ax=\left(\begin{array}{c}\frac{7}{3}\\ \frac{14}{3}\\ \frac{14}{3}\end{array}\right),\Vert r{\Vert }_2=7$

$$A=qr=\left(\begin{array}{ccc}-2/3& -2/3& -1/3\\ -1/3& 2/3& -2/3\\ 2/3& -1/3& -2/3\end{array}\right)\left(\begin{array}{cc}6& -3\\ 0& 9\\ 0& 0\end{array}\right)$$

$q.q^{\top }r=\left(\begin{array}{cc}6& -3\\ 0& 9\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{c}-5\\ 4\end{array}\right)$ 解得 $x_1=\frac{-11}{18},x_2=\frac{4}{9}$ 二范数误差： $r=\Vert q^{\top }.b-r.x{\Vert }_2=7$

---

运行结果

APPENDIX

import math
import numpy as np
 
 
def poly_matr(x):
    y = [[1, x_i, math.pow(x_i, 2), math.pow(x_i, 3), math.pow(x_i, 4), math.pow(x_i, 5), math.pow(x_i, 6), math.pow(x_i, 7), math.pow(x_i, 8), math.pow(x_i, 9), math.pow(x_i, 10)] for x_i in x]  # noqa: E501
    return y
 
 
def QRDecomposition(A):
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))
    for j in range(n):
        v = A[:, j]
        for i in range(j):
            R[i, j] = np.dot(Q[:, i].T, A[:, j])
            v = v - R[i, j] * Q[:, i]
        R[j, j] = np.linalg.norm(v)
        Q[:, j] = v / R[j, j]
    return Q, R
 
 
def least_squares_polynomial():
    x = np.linspace(2,4,num=20,endpoint=True)
    y = [sum(math.pow(x_i, i) for i in range(11)) for x_i in x]
    poly_matrix = poly_matr(x)
    poly_matrix_transpose = np.transpose(poly_matrix)
    poly_matrix_transpose_poly_matrix = np.dot(poly_matrix_transpose, poly_matrix)
    poly_matrix_transpose_poly_matrix_inv = np.linalg.inv(poly_matrix_transpose_poly_matrix)
    poly_matrix_transpose_y = np.dot(poly_matrix_transpose, y)
    a = np.dot(poly_matrix_transpose_poly_matrix_inv, poly_matrix_transpose_y)
    return a
 
 
def QR_least_squares_polynomial():
    x = np.linspace(2,4,num=20,endpoint=True)
    y = [sum(math.pow(x_i, i) for i in range(11)) for x_i in x]
    poly_matrix = poly_matr(x)
    poly_matrix_transpose = np.transpose(poly_matrix)
    Q, R = QRDecomposition(poly_matrix_transpose)
    R_inv = np.linalg.inv(R)
    Q_transpose = np.transpose(Q)
    a = np.dot(np.transpose(np.dot(R_inv, Q_transpose)), y)
    return a
 
 
if __name__ == "__main__":
    print(least_squares_polynomial())
    print("______________________")
    print(QR_least_squares_polynomial())