# 任意选一个你喜欢的整数，这能帮你得到稳定的结果
seed = 81 矩阵运算
1.1 创建一个 4*4 的单位矩阵
# 这个项目设计来帮你熟悉 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):
return len(M),len(M[0])# 运行以下代码测试你的 shape 函数
%run -i -e test.py LinearRegressionTestCase.test_shape1.3 每个元素四舍五入到特定小数数位¶
# TODO 每个元素四舍五入到特定小数数位
# 直接修改参数矩阵，无返回值
def matxRound(M, decPts=4):
for i in range(len(M)):
for j in range(len(M[i])):
M[i][j] = round(M[i][j],decPts)
pass# 运行以下代码测试你的 matxRound 函数
%run -i -e test.py LinearRegressionTestCase.test_matxRound1.4 计算矩阵的转置¶
# TODO 计算矩阵的转置
def transpose(M):
return [list(row) for row in zip(*M)]# 运行以下代码测试你的 transpose 函数
%run -i -e test.py LinearRegressionTestCase.test_transpose1.5 计算矩阵乘法 AB¶
# TODO 计算矩阵乘法 AB，如果无法相乘则raise ValueError
def matxMultiply(A, B):
multiply = []
if len(A[0]) != len(B):
raise ValueError
result = [list(row) for row in zip(*B)]
for Al in range(len(A)):
row =[]
for Bl in range(len(result)):
num = 0
for Br in range(len(result[0])):
num += A[Al][Br] * result[Bl][Br]
row.append(num)
multiply.append(row)
return multiply# 运行以下代码测试你的 matxMultiply 函数
%run -i -e test.py LinearRegressionTestCase.test_matxMultiply2 Gaussign Jordan 消元法
2.1 构造增广矩阵
# TODO 构造增广矩阵，假设A，b行数相同
def augmentMatrix(A, b):
return [AA + bb for AA, bb in zip(A,b)]# 运行以下代码测试你的 augmentMatrix 函数
%run -i -e test.py LinearRegressionTestCase.test_augmentMatrix2.2 初等行变换
- 交换两行
- 把某行乘以一个非零常数
- 把某行加上另一行的若干倍：
# TODO r1 <---> r2 # 直接修改参数矩阵，无返回值 def swapRows(M, r1, r2): M[r1],M[r2] = M[r2],M[r1] pass# 运行以下代码测试你的 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: raise ValueError for i in range(len(M[r])): M[r][i] *= scale pass# 运行以下代码测试你的 scaleRow 函数 %run -i -e test.py LinearRegressionTestCase.test_scaleRow# TODO r1 <--- r1 + r2*scale # 直接修改参数矩阵，无返回值 def addScaledRow(M, r1, r2, scale): for i in range(len(M[r1])): M[r1][i] += M[r2][i] * scale pass# 运行以下代码测试你的 addScaledRow 函数 %run -i -e test.py LinearRegressionTestCase.test_addScaledRow2.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为奇异矩阵两种情况。
-
-
推演有以下要求:
- 展示每一列的消元结果, 比如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}来渲染分数，如下： 
-
# 不要修改这里！ 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}来渲染分数，如下： 
-
2.3.3 实现 Gaussian Jordan 消元法
from decimal import Decimal from fractions import Fraction # 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 len(A)!= len(b) : return None # 构造扩增矩阵 result = augmentMatrix(A, b) # 转换为阶梯型矩阵 for j in range(len(result[0])-1): row,maxNum = 0,0 for i in range(len(result)): if i >= j: if abs(result[i][j]) > maxNum: maxNum = abs(result[i][j]) row = i # 返回奇异矩阵 if(abs(maxNum) < epsilon): return None else: # 找出最大放到最上面 if(row != j): swapRows(result,j,row) row,maxNum = 0,0 scaleRow(result,j,(float(1)/float(result[j][j]))) # 消除下面 for dis in range(len(result)): if dis != j: if abs(-result[dis][j]) > epsilon: addScaledRow(result,dis,j,-result[dis][j]) newResult =[] for i in range(len(result)): row =[] row.append(result[i][len(result)]) newResult.append(row) return newResult# 运行以下代码测试你的 gj_Solve 函数 %run -i -e test.py LinearRegressionTestCase.test_gj_Solve3 线性回归
-
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()
3.2 拟合一条直线
3.2.1 猜测一条直线¶
#TODO 请选择最适合的直线 y = mx + b m1 = 3.7 b1 = 14.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方向距离的平方的平均数，表达式如下：
???=1?∑?=1?(??−???−?)2
# TODO 实现以下函数并输出所选直线的MSE def calculateMSE(X,Y,m,b): return sum([(y-m*x -b)**2 for x,y in zip(X,Y)])/len(X) print(calculateMSE(X,Y,m1,b1))3.2.3 调整参数 ?,?m,b 来获得最小的平方平均误差
你可以调整3.2.1中的参数 ?1,?1m1,b1 让蓝点均匀覆盖在红线周围，然后微调 ?1,?1m1,b1 让MSE最小。
3.3 (选做) 找到参数 ?,?m,b 使得平方平均误差最小
这一部分需要简单的微积分知识( (?2)′=2?(x2)′=2x )。因为这是一个线性代数项目，所以设为选做。
刚刚我们手动调节参数，尝试找到最小的平方平均误差。下面我们要精确得求解 ?,?m,b 使得平方平均误差最小。
定义目标函数 E 为
-
-
因为 ?=?2???E=n2MSE, 所以 ?E 取到最小值时，???MSE 也取到最小值。要找到 ?E 的最小值，即要找到 ?,?m,b 使得 ?E 相对于 ?m, ?E 相对于 ?b 的偏导数等于0.
因此我们要解下面的方程组。
∂?∂?=0∂?∂?=0{∂E∂m=0∂E∂b=0
3.3.1 计算目标函数相对于参数的导数
首先我们计算两个式子左边的值
证明/计算：
∂?∂?=∑?=1?−??(??−???−?)∂E∂m=∑i=1n−xi(yi−mxi−b)
∂?∂?=∑?=1?−(??−???−?)∂E∂b=∑i=1n−(yi−mxi−b)
TODO 证明:
3.3.2 实例推演
现在我们有了一个二元二次方程组
∑?=1?−??(??−???−?)=0∑?=1?−(??−???−?)=0{∑i=1n−xi(yi−mxi−b)=0∑i=1n−(yi−mxi−b)=0
为了加强理解，我们用一个实际例子演练。
我们要用三个点 (1,1),(2,2),(3,2)(1,1),(2,2),(3,2) 来拟合一条直线 y = m*x + b, 请写出
- 目标函数 ?E,
- 二元二次方程组，
- 并求解最优参数 ?,?m,b
-
3.4 求解 ???ℎ=???XTXh=XTY
在3.3 中，我们知道线性回归问题等价于求解 ???ℎ=???XTXh=XTY (如果你选择不做3.3，就勇敢的相信吧，哈哈)
TODO 写出目标函数，方程组和最优参数
3.3.3 将方程组写成矩阵形式
我们的二元二次方程组可以用更简洁的矩阵形式表达，将方程组写成矩阵形式更有利于我们使用 Gaussian Jordan 消元法求解。
请证明
[∂?∂?∂?∂?]=???ℎ−???[∂E∂m∂E∂b]=XTXh−XTY
其中向量 ?Y, 矩阵 ?X 和 向量 ℎh 分别为 :
?=?1?2...??,?=?1?2...??11...1,ℎ=[??]Y=[y1y2...yn],X=[x11x21......xn1],h=[mb]
TODO 证明:
至此我们知道，通过求解方程 ???ℎ=???XTXh=XTY 来找到最优参数。这个方程十分重要，他有一个名字叫做 Normal Equation，也有直观的几何意义。你可以在 子空间投影 和 投影矩阵与最小二乘 看到更多关于这个方程的内容。
-
3.4 求解 ???ℎ=???XTXh=XTY
在3.3 中，我们知道线性回归问题等价于求解 ???ℎ=???XTXh=XTY (如果你选择不做3.3，就勇敢的相信吧，哈哈)
import numpy as np ''' 参数：X, Y 存储着一一对应的横坐标与纵坐标的两个一维数组 返回：m，b 浮点数 ''' def linearRegression(X,Y): x_num = [[1]] * len(X) newX =augmentMatrix (map(list,zip(X)),x_num) newY = map(list,zip(Y)) A = matxMultiply (transpose(newX),newX) b = matxMultiply(transpose(newX) ,newY) result = gj_Solve(A,b,decPts=4,epsilon=1.0e-16) print result[0][0],result[1][0] return result[0][0],result[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"你求得的回归结果是什么？ 请使用运行以下代码将它画出来。
-
# 请不要修改下面的代码 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()print(calculateMSE(X,Y,m2,b2))
5252
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
