```python
# 任意选一个你喜欢的整数，这能帮你得到稳定的结果
seed = 2333 # todo
```

# 欢迎来到线性回归项目

若项目中的题目有困难没完成也没关系，我们鼓励你带着问题提交项目，评审人会给予你诸多帮助。

所有选做题都可以不做，不影响项目通过。如果你做了，那么项目评审会帮你批改，也会因为选做部分做错而判定为不通过。

其中非代码题可以提交手写后扫描的 pdf 文件，或使用 Latex 在文档中直接回答。

# 1 矩阵运算

## 1.1 创建一个 4*4 的单位矩阵


```python
# 这个项目设计来帮你熟悉 python list 和线性代数
# 你不能调用任何NumPy以及相关的科学计算库来完成作业


# 本项目要求矩阵统一使用二维列表表示，如下：
A = [[1,2,3], 
     [2,3,3], 
     [1,2,5]]

B = [[1,2,3,5], 
     [2,3,3,5], 
     [1,2,5,1]]

# 向量也用二维列表表示
C = [[1],
     [2],
     [3]]

#TODO 创建一个 4*4 单位矩阵
I = [[1,0,0,0],
     [0,1,0,0],
     [0,0,1,0],
     [0,0,0,1]]
```

## 1.2 返回矩阵的行数和列数


```python
# TODO 返回矩阵的行数和列数
def shape(M):
    return len(M), len(M[0])
```


```python
# 运行以下代码测试你的 shape 函数
%run -i -e test.py LinearRegressionTestCase.test_shape
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.001s
    
    OK


## 1.3 每个元素四舍五入到特定小数数位


```python
# TODO 每个元素四舍五入到特定小数数位
# 直接修改参数矩阵，无返回值
def matxRound(M, decPts=4):
    for row in range(len(M)):
        for col in range(len(M[row])):
            M[row][col] = round(M[row][col], decPts)
```


```python
# 运行以下代码测试你的 matxRound 函数
%run -i -e test.py LinearRegressionTestCase.test_matxRound
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.010s
    
    OK


## 1.4 计算矩阵的转置


```python
# TODO 计算矩阵的转置
def transpose(M):
    rows, cols = shape(M)
    new = []
    for col in range(cols):
        new.append([])
        for row in range(rows):
            new[col].append(0)
    for row in range(rows):
        for col in range(cols):
            new[col][row] = M[row][col]
    return new

#     Pythonic approach
#     return [list(col) for col in zip(*M)]
```


```python
# 运行以下代码测试你的 transpose 函数
%run -i -e test.py LinearRegressionTestCase.test_transpose
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.030s
    
    OK


## 1.5 计算矩阵乘法 AB


```python
# TODO 计算矩阵乘法 AB，如果无法相乘则raise ValueError
def dot(a, b):
    return sum(i * j for i, j in zip(a, b))

def matxMultiply(A, B):
    n, m = shape(A)
    p, q = shape(B)
    if m != p:
        raise ValueError('dimension of matrices don\'t match')
    result = []
    for row in range(n):
        result.append([])
        for col in range(q):
            # ith row jth column = dot product of ith row of A and jth column of B
            # im not using [row][col] to index but rather just appending...
            # just because i can
            result[row].append(dot(A[row], [i[col] for i in B]))
    return result
```


```python
# 运行以下代码测试你的 matxMultiply 函数
%run -i -e test.py LinearRegressionTestCase.test_matxMultiply
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.124s
    
    OK


---

# 2 Gaussian Jordan 消元法

## 2.1 构造增广矩阵

$ A = \begin{bmatrix}
    a_{11}    & a_{12} & ... & a_{1n}\\
    a_{21}    & a_{22} & ... & a_{2n}\\
    a_{31}    & a_{22} & ... & a_{3n}\\
    ...    & ... & ... & ...\\
    a_{n1}    & a_{n2} & ... & a_{nn}\\
\end{bmatrix} , b = \begin{bmatrix}
    b_{1}  \\
    b_{2}  \\
    b_{3}  \\
    ...    \\
    b_{n}  \\
\end{bmatrix}$

返回 $ Ab = \begin{bmatrix}
    a_{11}    & a_{12} & ... & a_{1n} & b_{1}\\
    a_{21}    & a_{22} & ... & a_{2n} & b_{2}\\
    a_{31}    & a_{22} & ... & a_{3n} & b_{3}\\
    ...    & ... & ... & ...& ...\\
    a_{n1}    & a_{n2} & ... & a_{nn} & b_{n} \end{bmatrix}$


```python
# TODO 构造增广矩阵，假设A，b行数相同
def augmentMatrix(A, b):
#     assert n == shape(b)[0] # this is in the assumption!
    new = [A[row] + b[row] for row in range(len(A))]
    return new
```


```python
# 运行以下代码测试你的 augmentMatrix 函数
%run -i -e test.py LinearRegressionTestCase.test_augmentMatrix
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.004s
    
    OK


## 2.2 初等行变换
- 交换两行
- 把某行乘以一个非零常数
- 把某行加上另一行的若干倍：


```python
# TODO r1 <---> r2
# 直接修改参数矩阵，无返回值
def swapRows(M, r1, r2):
    if r1 != r2:
        M[r1], M[r2] = M[r2], M[r1]
```


```python
# 运行以下代码测试你的 swapRows 函数
%run -i -e test.py LinearRegressionTestCase.test_swapRows
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.004s
    
    OK



```python
## TODO r1 <--- r1 * scale
# scale为0是非法输入，要求 raise ValueError
# 直接修改参数矩阵，无返回值
def scaleRow(M, r, scale):
    if scale == 0:
        raise ValueError('cannot scale a matrix by zero')
    for col in range(len(M[r])):
        M[r][col] *= scale
```


```python
# 运行以下代码测试你的 scaleRow 函数
%run -i -e test.py LinearRegressionTestCase.test_scaleRow
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.005s
    
    OK



```python
# TODO r1 <--- r1 + r2*scale
# 直接修改参数矩阵，无返回值
def addScaledRow(M, r1, r2, scale):
    if scale == 0:
        raise ValueError('cannot scale a matrix by zero')
    for col in range(len(M[r1])):
        M[r1][col] += scale * M[r2][col]
        
```


```python
# 运行以下代码测试你的 addScaledRow 函数
%run -i -e test.py LinearRegressionTestCase.test_addScaledRow
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 0.002s
    
    OK


## 2.3  Gaussian Jordan 消元法求解 Ax = b

### 2.3.1 算法

步骤1 检查A，b是否行数相同

步骤2 构造增广矩阵Ab

步骤3 逐列转换Ab为化简行阶梯形矩阵 [中文维基链接](https://zh.wikipedia.org/wiki/%E9%98%B6%E6%A2%AF%E5%BD%A2%E7%9F%A9%E9%98%B5#.E5.8C.96.E7.AE.80.E5.90.8E.E7.9A.84-.7Bzh-hans:.E8.A1.8C.3B_zh-hant:.E5.88.97.3B.7D-.E9.98.B6.E6.A2.AF.E5.BD.A2.E7.9F.A9.E9.98.B5)
    
    对于Ab的每一列（最后一列除外）
        当前列为列c
        寻找列c中 对角线以及对角线以下所有元素（行 c~N）的绝对值的最大值
        如果绝对值最大值为0
            那么A为奇异矩阵，返回None (你可以在选做问题2.4中证明为什么这里A一定是奇异矩阵)
        否则
            使用第一个行变换，将绝对值最大值所在行交换到对角线元素所在行（行c） 
            使用第二个行变换，将列c的对角线元素缩放为1
            多次使用第三个行变换，将列c的其他元素消为0
            
步骤4 返回Ab的最后一列

**注：** 我们并没有按照常规方法先把矩阵转化为行阶梯形矩阵，再转换为化简行阶梯形矩阵，而是一步到位。如果你熟悉常规方法的话，可以思考一下两者的等价性。

### 2.3.2 算法推演

为了充分了解Gaussian Jordan消元法的计算流程，请根据Gaussian Jordan消元法，分别手动推演矩阵A为***可逆矩阵***，矩阵A为***奇异矩阵***两种情况。

#### 推演示例 


$Ab = \begin{bmatrix}
    -7 & 5 & -1 & 1\\
    1 & -3 & -8 & 1\\
    -10 & -2 & 9 & 1\end{bmatrix}$

$ --> $
$\begin{bmatrix}
    1 & \frac{1}{5} & -\frac{9}{10} & -\frac{1}{10}\\
    0 & -\frac{16}{5} & -\frac{71}{10} & \frac{11}{10}\\
    0 & \frac{32}{5} & -\frac{73}{10} & \frac{3}{10}\end{bmatrix}$

$ --> $
$\begin{bmatrix}
    1 & 0 & -\frac{43}{64} & -\frac{7}{64}\\
    0 & 1 & -\frac{73}{64} & \frac{3}{64}\\
    0 & 0 & -\frac{43}{4} & \frac{5}{4}\end{bmatrix}$

$ --> $
$\begin{bmatrix}
    1 & 0 & 0 & -\frac{3}{16}\\
    0 & 1 & 0 & -\frac{59}{688}\\
    0 & 0 & 1 & -\frac{5}{43}\end{bmatrix}$
    

#### 推演有以下要求:
1. 展示每一列的消元结果, 比如3*3的矩阵, 需要写三步
2. 用分数来表示
3. 分数不能再约分
4. 我们已经给出了latex的语法,你只要把零改成你要的数字(或分数)即可
5. 检查你的答案, 可以用[这个](http://www.math.odu.edu/~bogacki/cgi-bin/lat.cgi?c=sys), 或者后面通过单元测试后的`gj_Solve`

_你可以用python的 [fractions](https://docs.python.org/2/library/fractions.html) 模块辅助你的约分_

#### 以下开始你的尝试吧!


```python
# 不要修改这里！
from helper import *
A = generateMatrix(3,seed,singular=False)
b = np.ones(shape=(3,1),dtype=int) # it doesn't matter
Ab = augmentMatrix(A.tolist(),b.tolist()) # 请确保你的增广矩阵已经写好了
printInMatrixFormat(Ab,padding=3,truncating=0)
```

    -10,  9,  5 ||  1 
     -4,  3, -4 ||  1 
     -2,  3,  5 ||  1 


请按照算法的步骤3，逐步推演***可逆矩阵***的变换。

在下面列出每一次循环体执行之后的增广矩阵。

要求：
1. 做分数运算
2. 使用`\frac{n}{m}`来渲染分数，如下：
 - $\frac{n}{m}$
 - $-\frac{a}{b}$


$ Ab = \begin{bmatrix}
    -10 & 9 & 5 & 1 \\
    -4 & 3 & -4 & 1 \\
    -2 & 3 & 5 & 1 \end{bmatrix}$

$ --> \begin{bmatrix}
    1 & -\frac{9}{10} & -\frac{5}{10} & -\frac{1}{10} \\
    0 & -\frac{3}{5} & -6 & \frac{3}{5} \\
    0 & \frac{6}{5} & 4 & \frac{4}{5} \end{bmatrix}$
    
$ --> \begin{bmatrix}
    1 & 0 & \frac{5}{2} & \frac{1}{2} \\
    0 & 1 & \frac{10}{3} & \frac{2}{3} \\
    0 & 0 & -4 & 1 \end{bmatrix}$
    
$ --> \begin{bmatrix}
    1 & 0 & 0 & \frac{9}{8} \\
    0 & 1 & 0 & \frac{3}{2} \\
    0 & 0 & 1 & -\frac{1}{4} \end{bmatrix}$


```python
# 不要修改这里！
A = generateMatrix(3,seed,singular=True)
b = np.ones(shape=(3,1),dtype=int)
Ab = augmentMatrix(A.tolist(),b.tolist()) # 请确保你的增广矩阵已经写好了
printInMatrixFormat(Ab,padding=3,truncating=0)
```

      6,  1,  9 ||  1 
      0,  1,  1 ||  1 
     -6,  1, -7 ||  1 


请按照算法的步骤3，逐步推演***奇异矩阵***的变换。

在下面列出每一次循环体执行之后的增广矩阵。

要求：
1. 做分数运算
2. 使用`\frac{n}{m}`来渲染分数，如下：
 - $\frac{n}{m}$
 - $-\frac{a}{b}$


$ Ab = \begin{bmatrix}
    6 & 1 & 9 & 1 \\
    0 & 1 & 1 & 1 \\
    -6 & 1 & -7 & 1 \end{bmatrix}$

$ --> \begin{bmatrix}
    1 & \frac{1}{6} & \frac{3}{2} & \frac{1}{6} \\
    0 & 1 & 1 & 1 \\
    0 & 2 & 2 & 2 \end{bmatrix}$
    
$ --> \begin{bmatrix}
    1 & 0 & \frac{4}{3} & 0 \\
    0 & 1 & 1 & 0 \\
    0 & 0 & 0 & 0 \end{bmatrix}$
    

### 2.3.3 实现 Gaussian Jordan 消元法


```python
# TODO 实现 Gaussain Jordan 方法求解 Ax = b

""" Gaussian Jordan 方法求解 Ax = b.
    参数
        A: 方阵 
        b: 列向量
        decPts: 四舍五入位数，默认为4
        epsilon: 判读是否为0的阈值，默认 1.0e-16
        
    返回列向量 x 使得 Ax = b 
    返回None，如果 A，b 高度不同
    返回None，如果 A 为奇异矩阵
"""

from pprint import pprint

def gj_Solve(A, b, decPts=4, epsilon=1.0e-16):
    size = len(A)
    if size != len(b):
        return None
    C = augmentMatrix(A, b)
    for c in range(size):
        absCol = [abs(row[c]) for row in C[c:]]
        maxNum = max(absCol)
        if maxNum < epsilon: # singular matrix
            return None
        swapRows(C, absCol.index(maxNum) + c, c)
        scaleRow(C, c, 1 / C[c][c])
        for r in range(size):
            if r == c or abs(C[r][c]) < epsilon:
                continue
            addScaledRow(C, r, c, -C[r][c])
    solution = [[row[-1]] for row in C]
    return solution
```


```python
# 运行以下代码测试你的 gj_Solve 函数
%run -i -e test.py LinearRegressionTestCase.test_gj_Solve
```

    .
    ----------------------------------------------------------------------
    Ran 1 test in 4.122s
    
    OK


###### (选做) 2.4 算法正确判断了奇异矩阵：

在算法的步骤3 中，如果发现某一列对角线和对角线以下所有元素都为0，那么则断定这个矩阵为奇异矩阵。

我们用正式的语言描述这个命题，并证明为真。

证明下面的命题：

**如果方阵 A 可以被分为4个部分: ** 

$ A = \begin{bmatrix}
    I    & X \\
    Z    & Y \\
\end{bmatrix} , \text{其中 I 为单位矩阵，Z 为全0矩阵，Y 的第一列全0}$，

**那么A为奇异矩阵。**

提示：从多种角度都可以完成证明
- 考虑矩阵 Y 和 矩阵 A 的秩
- 考虑矩阵 Y 和 矩阵 A 的行列式
- 考虑矩阵 A 的某一列是其他列的线性组合

**TODO** 证明：

由于 $A,\ I$ 为方阵，则 $Y$ 为方阵，易知，$X$ 的行数与 $Z$ 的列数一致。  
设 $n$ 为 $I$ 的大小，并设 $x_1,\ x_2,\dots,\ x_n$ 为矩阵 $X$ 的第一列的元素；设 $A_i$ 为矩阵 $A$ 的第 $i$ 列。

由于 $Y$ 的第一列均为 0，那么 $A_{n+1}$（$Y$ 所在的这一列）可以被表示为
$$ A_{n+1} = x_1 A_1 + x_2 A_2 + \cdots + x_n A_n = \sum_{i=1}^n x_i A_i $$

因此 $A_{n+1}$ 这一列并不是线性独立的（是 $A$ 前 $n$ 列的线性组合），又因 $A$ 为方阵，则得出 $A$ 不满秩，即 $A$ 是奇异矩阵，证毕。

# 3  线性回归

## 3.1 随机生成样本点


```python
# 不要修改这里！
# 运行一次就够了！
from helper import *
from matplotlib import pyplot as plt
%matplotlib inline

X,Y = generatePoints(seed,num=100)

## 可视化
plt.xlim((-5,5))
plt.xlabel('x',fontsize=18)
plt.ylabel('y',fontsize=18)
plt.scatter(X,Y,c='b')
plt.show()
```

## 3.2 拟合一条直线

### 3.2.1 猜测一条直线


```python
#TODO 请选择最适合的直线 y = mx + b
m1 = 2/5
b1 = 13.8

# 不要修改这里！
plt.xlim((-5,5))
x_vals = plt.axes().get_xlim()
y_vals = [m1*x+b1 for x in x_vals]
plt.plot(x_vals, y_vals, '-', color='r')

plt.xlabel('x',fontsize=18)
plt.ylabel('y',fontsize=18)
plt.scatter(X,Y,c='b')

plt.show()
```

### 3.2.2 计算平均平方误差 (MSE)

我们要编程计算所选直线的平均平方误差(MSE), 即数据集中每个点到直线的Y方向距离的平方的平均数，表达式如下：
$$
MSE = \frac{1}{n}\sum_{i=1}^{n}{(y_i - mx_i - b)^2}
$$


```python
# TODO 实现以下函数并输出所选直线的MSE

def calculateMSE(X,Y,m,b):
    error = 0
    for x, y in zip(X, Y):
        diff = (y - m * x - b) ** 2
        error += diff
    error /= len(X)
    return error

#     Pythonic approach
#     return sum([(y - m * x - b) ** 2 for x, y in zip(X, Y)]) / len(Y)

print(calculateMSE(X,Y,m1,b1))
```

    1.489379289068421


### 3.2.3 调整参数 $m, b$ 来获得最小的平方平均误差

你可以调整3.2.1中的参数 $m1,b1$ 让蓝点均匀覆盖在红线周围，然后微调 $m1, b1$ 让MSE最小。

## 3.3 (选做) 找到参数 $m, b$ 使得平方平均误差最小

**这一部分需要简单的微积分知识(  $ (x^2)' = 2x $ )。因为这是一个线性代数项目，所以设为选做。**

刚刚我们手动调节参数，尝试找到最小的平方平均误差。下面我们要精确得求解 $m, b$ 使得平方平均误差最小。

定义目标函数 $E$ 为
$$
E = \frac{1}{2}\sum_{i=1}^{n}{(y_i - mx_i - b)^2}
$$

因为 $E = \frac{n}{2}MSE$, 所以 $E$ 取到最小值时，$MSE$ 也取到最小值。要找到 $E$ 的最小值，即要找到 $m, b$ 使得 $E$ 相对于 $m$, $E$ 相对于 $b$ 的偏导数等于0. 

因此我们要解下面的方程组。

$$
\begin{cases}
\displaystyle
\frac{\partial E}{\partial m} =0 \\
\\
\displaystyle
\frac{\partial E}{\partial b} =0 \\
\end{cases}
$$

### 3.3.1 计算目标函数相对于参数的导数
首先我们计算两个式子左边的值

证明/计算：
$$
\frac{\partial E}{\partial m} = \sum_{i=1}^{n}{-x_i(y_i - mx_i - b)}
$$

$$
\frac{\partial E}{\partial b} = \sum_{i=1}^{n}{-(y_i - mx_i - b)}
$$

**TODO** 证明:

设 $u_i = (y_i - mx_i - b)$，我们可以重写目标函数 $E$ 为
$$ E = \frac{1}{2} \sum_{i=1}^n u_i^2 $$
设 $E_i = \frac{1}{2} u_i^2$，我们可以再重写目标函数 $E$ 为（trivially）
$$ E = \sum_{i=1}^n E_i$$
则，我们有
$$ \frac{\partial E}{\partial m} = \sum_{i=1}^n \frac{\partial E_i}{\partial m} \quad\textrm{和}\quad \frac{\partial E}{\partial b} = \sum_{i=1}^n \frac{\partial E_i}{\partial b} $$

我们先求 $E_i$ 关于 $m$ 的偏导数，即 $\frac{\partial E_i}{\partial m}$，根据链式法则，我们有，对于任意的 $i \in \{1,2,\dots,n\}$
$$ \frac{\partial E_i}{\partial m} = \frac{\mathrm{d} E_i}{\mathrm{d} u_i} \cdot \frac{\partial u_i}{\partial m} $$
在 $u_i$ 中，只有一项 $-mx_i$ 是和 $m$ 相关的，因此 $\frac{\partial u_i}{\partial m} = -x_i$，即
$$ \frac{\partial E_i}{\partial m} = u_i\cdot (-x_i) = -x_i(y_i - mx_i - b) $$
即
$$ \frac{\partial E}{\partial m} = \sum_{i=1}^n -x_i(y_i - mx_i - b) $$

同理，我们可得
$$ \frac{\partial E_i}{\partial b} = \frac{\partial E_i}{\partial u_i}\cdot\frac{\partial u_i}{\partial b} = u_i\cdot (-1) = -(y_i - mx_i - b) $$
即
$$ \frac{\partial E}{\partial b} = \sum_{i=1}^n -(y_i - mx_i - b) $$
证毕。

### 3.3.2 实例推演

现在我们有了一个二元二次方程组

$$
\begin{cases}
\displaystyle
\sum_{i=1}^{n}{-x_i(y_i - mx_i - b)} =0 \\
\displaystyle
\sum_{i=1}^{n}{-(y_i - mx_i - b)} =0 \\
\end{cases}
$$

为了加强理解，我们用一个实际例子演练。

我们要用三个点 $(1,1),\ (2,2),\ (3,2)$ 来拟合一条直线 $y = mx + b$, 请写出

- 目标函数 $E$, 
- 二元二次方程组，
- 并求解最优参数 $m, b$

**TODO** 写出目标函数，方程组和最优参数

1. 目标函数 $E$

\begin{align}
    E &= \frac{1}{2}\sum_{i=1}^n(y_i - mx_i - b)^2 \\
      &= \frac{1}{2}\left((1 - m - b)^2 + (2 - 2m - b)^2 + (3 - 2m - b)^2\right) \\
      &= \frac{1}{2}\left(9m^2 - 22m + 10mb - 12b + 3b^2 + 14\right) \\
      &= \frac{9}{2}m^2 - 11m + 5mb - 6b + \frac{3}{2}b^2 + 7
\end{align}

2. 二元二次方程组

$$\begin{cases}
    \displaystyle \sum_{i=1}^n -x_i(y_i - mx_i - b) = 0 \\
    \displaystyle \sum_{i=1}^n -(y_i - mx_i - b) = 0 \\
\end{cases} \; \Rightarrow \;
\begin{cases}
    \displaystyle -(1 - m - b) - 2(2 - 2m - b) - 3(2 - 3m - b) = 0 \\
    \displaystyle -(1 - m - b) - (2 - 2m - b) - (2 - 3m - b) = 0 \\
\end{cases} \; \Rightarrow \;
\begin{cases}
    \displaystyle -11 + 14m + 6b = 0 \\
    \displaystyle -5 + 8m + 3b = 0 \\
\end{cases}$$

3. 求解最优参数 $m,\ b$

\begin{cases}
    -11 + 14m + 6b = 0 \\
    -5 + 8m + 3b = 0 \\
\end{cases}

$$ (1) + -2 \times (2) $$
$$ \Downarrow $$

\begin{cases}
    \displaystyle -1 - 2m = 0 \quad \\
    \displaystyle -5 + 8m + 3b = 0 \quad \\
\end{cases}

$$ \therefore\ m = -\frac{1}{2} $$

$$ \therefore\ b = \frac{11 - 14m}{6} = \frac{18}{6} = 3 $$

$$\therefore\begin{cases}
    \displaystyle m = -\frac{1}{2}\\
    \displaystyle b = 3\\
\end{cases}$$

### 3.3.3 将方程组写成矩阵形式

我们的二元二次方程组可以用更简洁的矩阵形式表达，将方程组写成矩阵形式更有利于我们使用 Gaussian Jordan 消元法求解。

请证明 
$$
\begin{bmatrix}
    \frac{\partial E}{\partial m} \\
    \frac{\partial E}{\partial b} 
\end{bmatrix} = X^TXh - X^TY
$$

其中向量 $Y$, 矩阵 $X$ 和 向量 $h$ 分别为 :
$$
Y =  \begin{bmatrix}
    y_1 \\
    y_2 \\
    ... \\
    y_n
\end{bmatrix}
,
X =  \begin{bmatrix}
    x_1 & 1 \\
    x_2 & 1\\
    ... & ...\\
    x_n & 1 \\
\end{bmatrix},
h =  \begin{bmatrix}
    m \\
    b \\
\end{bmatrix}
$$

**TODO** 证明:

首先，
$$ X^\top X = \begin{bmatrix}
    \displaystyle \sum_{i=1}^n x_i^2 & \displaystyle \sum_{i=1}^n x_i \\
    \displaystyle \sum_{i=1}^n x_i & \displaystyle \sum_{i=1}^n 1
\end{bmatrix}$$
因此，
$$ X^\top Xh = \begin{bmatrix}
    \displaystyle \sum_{i=1}^n mx_i^2 + \sum_{i=1}^n bx_i \\
    \displaystyle \sum_{i=1}^n mx_i + \sum_{i=1}^n b
\end{bmatrix}$$
另外，
$$ X^\top Y = \begin{bmatrix}
    \displaystyle \sum_{i=1}^n x_i y_i \\
    \displaystyle \sum_{i=1}^n y_i
\end{bmatrix}$$
因此，
$$ X^\top Xh - X^\top Y = \begin{bmatrix}
    \displaystyle \sum_{i=1}^n mx_i^2 + \sum_{i=1}^n bx_i - \sum_{i=1}^n x_i y_i \\
    \displaystyle \sum_{i=1}^n mx_i + \sum_{i=1}^n b - \sum_{i=1}^n y_i
\end{bmatrix} = \begin{bmatrix}
    \displaystyle \sum_{i=1}^n (mx_i^2 + bx_i - x_i y_i) \\
    \displaystyle \sum_{i=1}^n (mx_i + b - y_i)
\end{bmatrix} = \begin{bmatrix}
    \displaystyle \sum_{i=1}^n -x_i(y_i - mx_i - b) \\
    \displaystyle \sum_{i=1}^n -(y_i - mx_i - b)
\end{bmatrix} = \begin{bmatrix}
    \displaystyle \frac{\partial E}{\partial m} \\
    \displaystyle \frac{\partial E}{\partial b} 
\end{bmatrix}$$
证毕。

至此我们知道，通过求解方程 $X^TXh = X^TY$ 来找到最优参数。这个方程十分重要，他有一个名字叫做 **Normal Equation**，也有直观的几何意义。你可以在 [子空间投影](http://open.163.com/movie/2010/11/J/U/M6V0BQC4M_M6V2AJLJU.html) 和 [投影矩阵与最小二乘](http://open.163.com/movie/2010/11/P/U/M6V0BQC4M_M6V2AOJPU.html) 看到更多关于这个方程的内容。

### 3.4 求解 $X^TXh = X^TY$ 

在3.3 中，我们知道线性回归问题等价于求解 $X^TXh = X^TY$ (如果你选择不做3.3，就勇敢的相信吧，哈哈)


```python
# TODO 实现线性回归
'''
参数：X, Y 存储着一一对应的横坐标与纵坐标的两个一维数组
返回：m，b 浮点数
'''
def linearRegression(X,Y):
    X = [[i, 1] for i in X]
    Y = [[i] for i in Y]
    A = matxMultiply(transpose(X), X)
    b = matxMultiply(transpose(X), Y)
    x = gj_Solve(A, b)
    return x[0][0], x[1][0]

m2,b2 = linearRegression(X,Y)
assert isinstance(m2,float),"m is not a float"
assert isinstance(b2,float),"b is not a float"
print(m2,b2)
```

    0.30111301039028665 13.867374963487162


你求得的回归结果是什么？
请使用运行以下代码将它画出来。


```python
# 请不要修改下面的代码
x1,x2 = -5,5
y1,y2 = x1*m2+b2, x2*m2+b2

plt.xlim((-5,5))
plt.xlabel('x',fontsize=18)
plt.ylabel('y',fontsize=18)
plt.scatter(X,Y,c='b')
plt.plot((x1,x2),(y1,y2),'r')
plt.title('y = {m:.4f}x + {b:.4f}'.format(m=m2,b=b2))
plt.show()
```

你求得的回归结果对当前数据集的MSE是多少？


```python
print(calculateMSE(X,Y,m2,b2))
```

    1.402604878826932


cbfbe7869a2df9bb93fa554b1ce81d5482748e0c
ipynb
Jupyter Notebook
udacity-dl/3-linear-algebra/Linear Regression Project.ipynb
yu-george/ml
8eedc62df9e7f37312bba39fa45bb9e9c6361028
Unlicense
udacity-dl/3-linear-algebra/Linear Regression Project.ipynb
yu-george/ml
8eedc62df9e7f37312bba39fa45bb9e9c6361028
Unlicense
udacity-dl/3-linear-algebra/Linear Regression Project.ipynb
yu-george/ml
8eedc62df9e7f37312bba39fa45bb9e9c6361028
Unlicense
Qwen/Qwen-72B
1. YES
2. YES