# 任意选一个你喜欢的整数，这能帮你得到稳定的结果
seed = 98661
欢迎来到线性回归项目¶
若项目中的题目有困难没完成也没关系，我们鼓励你带着问题提交项目，评审人会给予你诸多帮助。
所有选做题都可以不做，不影响项目通过。如果你做了，那么项目评审会帮你批改，也会因为选做部分做错而判定为不通过。
其中非代码题可以提交手写后扫描的 pdf 文件，或使用 Latex 在文档中直接回答。
# 这个项目设计来帮你熟悉 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 返回矩阵的行数和列数¶
# TODO 返回矩阵的行数和列数
def shape(M):
rows = len(M)
columns = len(M[0])
return rows,columns
# 运行以下代码测试你的 shape 函数
%run -i -e test.py LinearRegressionTestCase.test_shape
1.3 每个元素四舍五入到特定小数数位¶
# TODO 每个元素四舍五入到特定小数数位
# 直接修改参数矩阵，无返回值
def matxRound(M, decPts=4):
for i in range(shape(M)[0]):
for j in range(shape(M)[1]):
# print(i, j)
# print(M[i][j])
M[i][j] = round(M[i][j], decPts)
# 运行以下代码测试你的 matxRound 函数
%run -i -e test.py LinearRegressionTestCase.test_matxRound
1.4 计算矩阵的转置¶
# TODO 计算矩阵的转置
def transpose(M):
M_trans = []
for j in range(shape(M)[1]):
M_trans_rows = []
for i in range(shape(M)[0]):
M_trans_rows.append(M[i][j])
M_trans.append(M_trans_rows)
return M_trans
# 运行以下代码测试你的 transpose 函数
%run -i -e test.py LinearRegressionTestCase.test_transpose
1.5 计算矩阵乘法 AB¶
# TODO 计算矩阵乘法 AB，如果无法相乘则raise ValueError
def matxMultiply(A, B):
if shape(A)[1] == shape(B)[0]:
matx_multi = []
for i_A in range(shape(A)[0]):
matx_multi_rows = []
for j_B in range(shape(B)[1]):
matx_multi_items = []
for j_A in range(shape(A)[1]):
matx_multi_items.append(A[i_A][j_A] * B[j_A][j_B])
matx_multi_rows.append(sum(matx_multi_items))
matx_multi.append(matx_multi_rows)
return matx_multi
else:
raise ValueError("Matrix A's column number doesn't equal to Matrix b's row number")
# 运行以下代码测试你的 matxMultiply 函数
%run -i -e test.py LinearRegressionTestCase.test_matxMultiply
2 Gaussign 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}$
# TODO 构造增广矩阵，假设A，b行数相同
def augmentMatrix(A, b):
aug_matix = []
for i in range(shape(A)[0]):
aug_matix.append(A[i] + b[i])
return aug_matix
# 运行以下代码测试你的 augmentMatrix 函数
%run -i -e test.py LinearRegressionTestCase.test_augmentMatrix
2.2 初等行变换¶
- 交换两行
- 把某行乘以一个非零常数
- 把某行加上另一行的若干倍：
# TODO r1 <---> r2
# 直接修改参数矩阵，无返回值
def swapRows(M, r1, r2):
M[r1], M[r2] = M[r2], M[r1]
# 运行以下代码测试你的 swapRows 函数
%run -i -e test.py LinearRegressionTestCase.test_swapRows
# TODO r1 <--- r1 * scale
# scale为0是非法输入，要求 raise ValueError
# 直接修改参数矩阵，无返回值
def scaleRow(M, r, scale):
if scale != 0:
for j in range(shape(M)[1]):
M[r][j] = M[r][j] * scale
else:
raise ValueError("scale为0是非法输入")
# 运行以下代码测试你的 scaleRow 函数
%run -i -e test.py LinearRegressionTestCase.test_scaleRow
# TODO r1 <--- r1 + r2*scale
# 直接修改参数矩阵，无返回值
def addScaledRow(M, r1, r2, scale):
for j in range(shape(M)[1]):
M[r1][j] = M[r1][j] + M[r2][j] * scale
# 运行以下代码测试你的 addScaledRow 函数
%run -i -e test.py LinearRegressionTestCase.test_addScaledRow
2.3 Gaussian Jordan 消元法求解 Ax = b¶
2.3.1 算法¶
步骤1 检查A，b是否行数相同
步骤2 构造增广矩阵Ab
步骤3 逐列转换Ab为化简行阶梯形矩阵 中文维基链接
对于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}$
推演有以下要求:¶
- 展示每一列的消元结果, 比如3*3的矩阵, 需要写三步
- 用分数来表示
- 分数不能再约分
- 我们已经给出了latex的语法,你只要把零改成你要的数字(或分数)即可
- 检查你的答案, 可以用这个, 或者后面通过单元测试后的
gj_Solve
你可以用python的 fractions 模块辅助你的约分
以下开始你的尝试吧!¶
# 不要修改这里！
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)
请按照算法的步骤3，逐步推演可逆矩阵的变换。
在下面列出每一次循环体执行之后的增广矩阵。
要求：
- 做分数运算
- 使用
\frac{n}{m}来渲染分数，如下：- $\frac{n}{m}$
- $-\frac{a}{b}$
$ Ab = \begin{bmatrix} -7 & 9 & 0 & 1 \\ -10 & 7 & 6 & 1 \\ -4 & 3 & -5 & 1 \end{bmatrix}$
$ --> \begin{bmatrix} -10 & 7 & 6 & 1 \\ -7 & 9 & 0 & 1 \\ -4 & 3 & -5 & 1 \end{bmatrix}$
$ --> \begin{bmatrix} 1 & -\frac{7}{10} & -\frac{3}{5} & -\frac{1}{10} \\ -7 & 9 & 0 & 1 \\ -4 & 3 & -5 & 1 \end{bmatrix}$
$ --> \begin{bmatrix} 1 & -\frac{7}{10} & -\frac{3}{5} & -\frac{1}{10} \\ 0 & \frac{41}{10} & -\frac{21}{5} & \frac{3}{10} \\ 0 & \frac{1}{5} & -\frac{37}{5} & \frac{3}{5} \end{bmatrix}$
$ --> \begin{bmatrix} 1 & -\frac{7}{10} & -\frac{3}{5} & -\frac{1}{10} \\ 0 & 1 & -\frac{42}{41} & \frac{3}{41} \\ 0 & 0 & -\frac{295}{41} & \frac{24}{41} \end{bmatrix}$
$ --> \begin{bmatrix} 1 & 0 & -\frac{54}{41} & -\frac{2}{41} \\ 0 & 1 & -\frac{42}{41} & \frac{3}{41} \\ 0 & 0 & -\frac{295}{41} & \frac{24}{41} \end{bmatrix}$
$ --> \begin{bmatrix} 1 & 0 & -\frac{54}{41} & -\frac{2}{41} \\ 0 & 1 & -\frac{42}{41} & \frac{3}{41} \\ 0 & 0 & 1 & -\frac{24}{295} \end{bmatrix}$
$ --> \begin{bmatrix} 1 & 0 & -\frac{54}{41} & -\frac{2}{41} \\ 0 & 1 & 0 & -\frac{3}{295} \\ 0 & 0 & 1 & -\frac{24}{295} \end{bmatrix}$
$ --> \begin{bmatrix} 1 & 0 & 0 & -\frac{46}{295} \\ 0 & 1 & 0 & -\frac{3}{295} \\ 0 & 0 & 1 & -\frac{24}{295} \end{bmatrix}$
# 不要修改这里！
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)
请按照算法的步骤3，逐步推演奇异矩阵的变换。
在下面列出每一次循环体执行之后的增广矩阵。
要求：
- 做分数运算
- 使用
\frac{n}{m}来渲染分数，如下：- $\frac{n}{m}$
- $-\frac{a}{b}$
$ Ab = \begin{bmatrix} 6 & -7 & 0 & 1 \\ 1 & -10 & 0 & 1 \\ -6 & -2 & 0 & 1 \end{bmatrix}$
$ --> \begin{bmatrix} 1 & -\frac{7}{6} & 0 & \frac{1}{6} \\ 0 & -\frac{53}{6} & 0 & \frac{5}{6} \\ 0 & -9 & 0 & 2 \end{bmatrix}$
$ --> \begin{bmatrix} 1 & -\frac{7}{6} & 0 & \frac{1}{6} \\ 0 & -9 & 0 & 2 \\ 0 & -\frac{53}{6} & 0 & \frac{5}{6} \end{bmatrix}$
$ --> \begin{bmatrix} 1 & -\frac{7}{6} & 0 & \frac{1}{6} \\ 0 & 1 & 0 & -\frac{2}{9} \\ 0 & -\frac{53}{6} & 0 & \frac{5}{6} \end{bmatrix}$
$ --> \begin{bmatrix} 1 & 0 & 0 & -\frac{5}{54} \\ 0 & 1 & 0 & -\frac{2}{9} \\ 0 & 0 & 0 & -\frac{61}{54} \end{bmatrix}$
2.3.3 实现 Gaussian Jordan 消元法¶
# 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 为奇异矩阵
"""
def gj_Solve(A, b, decPts=4, epsilon = 1.0e-16):
if shape(A)[0] != shape(b)[0] or shape(A)[0] != shape(A)[1]:
return None # A、b 高度不同，返回 None
else:
Ab = augmentMatrix(A, b) # 构造增广矩阵Ab
for j in range(shape(A)[1]):
for i in range(j, shape(A)[0]):
if abs(Ab[i][j]) > abs(Ab[j][j]):
swapRows(Ab, j, i) # 使用第一个行变换，将绝对值最大值所在行交换到对角线元素所在行
if abs(Ab[j][j]) < epsilon: # 判断行交换后的对角线元素是否为0
return None # 奇异矩阵，返回 None
else:
scaleRow(Ab, j, 1 / Ab[j][j]) # 使用第二个行变换，将列j的对角线元素缩放为1
for i in range(shape(A)[0]):
if i != j:
addScaledRow(Ab, i, j, -Ab[i][j]) # 使用第三个行变换，将列j除对角线元素外的其他元素消为0
matxRound(Ab, decPts) # 每个元素四舍五入到特定小数数位
x = []
for i in range(shape(A)[0]):
x.append([Ab[i][-1]]) # 取最后一列作为返回列向量x
return x
# 运行以下代码测试你的 gj_Solve 函数
%run -i -e test.py LinearRegressionTestCase.test_gj_Solve
(选做) 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 证明：
3 线性回归¶
3.1 随机生成样本点¶
# 不要修改这里！
# 运行一次就够了！
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()
#TODO 请选择最适合的直线 y = mx + b
m1 = 4.3
b1 = 7.6
# 不要修改这里！
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} $$
# TODO 实现以下函数并输出所选直线的MSE
def calculateMSE(X,Y,m,b):
SE = 0
MSE = 0
for i in range(len(X)):
# print(X[i], Y[i])
SE = SE + pow((Y[i] - m * X[i] - b), 2)
MSE = SE / len(X)
return MSE
print(calculateMSE(X,Y,m1,b1))
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 证明:
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 = m*x + b, 请写出
- 目标函数 $E$,
- 二元二次方程组，
- 并求解最优参数 $m, b$
TODO 写出目标函数，方程组和最优参数
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 证明:
3.4 求解 $X^TXh = X^TY$¶
在3.3 中，我们知道线性回归问题等价于求解 $X^TXh = X^TY$ (如果你选择不做3.3，就勇敢的相信吧，哈哈)
# TODO 实现线性回归
'''
参数：X, Y 存储着一一对应的横坐标与纵坐标的两个一维数组
返回：m，b 浮点数
'''
def linearRegression(X, Y):
matx_X = []
matx_Y = []
for i in range(len(X)):
matx_X.append([X[i], 1]) # 把列表X 表示成矩阵形式(末尾增加各个元素为1 的一列)的matx_X
matx_Y.append([Y[i]]) # 把列表Y 表示成列向量形式的matx_Y
matx_X_T = transpose(matx_X) # 求matx_X 的转置矩阵
h = gj_Solve(matxMultiply(matx_X_T, matx_X), matxMultiply(matx_X_T, matx_Y)) # 求解向量h
return h[0][0], h[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)
你求得的回归结果是什么？ 请使用运行以下代码将它画出来。
# 请不要修改下面的代码
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是多少？
print(calculateMSE(X,Y,m2,b2))