机器学习之线性代数
Posted ZingpLiu
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了机器学习之线性代数相关的知识,希望对你有一定的参考价值。
说明
题目是优达学城机器学习入门线性代数作业。下面是我的实现。
工具为jupyter notebook,不用该工具请自行导入相关依赖。
完整内容已上传到github:https://github.com/zingp/data-analysis/blob/master/linear_algebra/linear_regression_project.ipynb。
本篇代码中引用的helper.py可到上面github上下载。
1 矩阵运算
1.1 创建一个4*4的单位矩阵
在创建矩阵之前注意选择seed:
# 任意选一个你喜欢的整数,这能帮你得到稳定的结果 seed = 9999
创建矩阵:
# 这个项目设计来帮你熟悉 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 返回矩阵的行数和列数
def shape(M): """返回矩阵的行列""" return len(M), len(M[0])
1.3 每个元素四舍五入到特定的小数位
# 每个元素四舍五入到特定小数数位 # 直接修改参数矩阵,无返回值 def matxRound(M, decPts=4): num_row,num_clo = shape(M) for r in range(num_row): for c in range(num_clo): M[r][c] = round(M[r][c], decPts)
1.4 计算矩阵的转置
def transpose(M): # *M 分解出列表中的子元素(子列表) # zip()将子列表中对应的元素打包成元组,返回包含一个个元组的列表 # 然后用列表推导式...真优雅啊 return [list(col) for col in zip(*M)]
1.5 计算矩阵乘法
# 计算矩阵乘法 AB,如果无法相乘则raise ValueError def matxMultiply(A, B): """矩阵乘法""" row_a, clo_a = shape(A) row_b, clo_b = shape(B) if clo_a == row_b: res = [] for i in range(row_a): res.append([]) for j in range(clo_b): ele_sum = 0 for s in range(clo_a): matx_ele = A[i][s] * B[s][j] if matx_ele is list: print(matx_ele) ele_sum += matx_ele res[i].append(ele_sum) return res else: raise ValueError
以上是我的实现,再看下充分利用列表递推式的实现方式:
def matxMultiply(A,B): _, c = shape(A) r, _ = shape(B) if c != r : raise ValueError Bt = transpose(B) result = [[sum((a*b) for a,b in zip(row,col)) for col in Bt] for row in A] return result
2 高斯消元法
2.1 构建增广矩阵
代码:
# 构造增广矩阵,假设A,b行数相同 def augmentMatrix(A, b): if len(A) != len(b): raise ValueError else: augment_mat = [] for r in range(shape(A)[0]): augment_mat.append([]) for c in range(shape(A)[1]): augment_mat[r].append(A[r][c]) augment_mat[r].append(b[r][0]) return augment_mat
再来看看利用列表递推式和zip函数的实现方式:
def augmentMatrix(A, b): return [ra + rb for ra,rb in zip(A,b)]
2.2 初等行变换
(1)交换两行
# r1 <---> r2 # 直接修改参数矩阵,无返回值 def swapRows(M, r1, r2): if (0 <= r1 < len(M)) and (0 <= r2 < len(M)): M[r1], M[r2] = M[r2], M[r1] else: raise IndexError(\'list index out of range\')
(2)某行乘以标量
# r1 <--- r1 * scale # scale为0是非法输入,要求 raise ValueError # 直接修改参数矩阵,无返回值 def scaleRow(M, r, scale): if not scale: raise ValueError(\'the parameter scale can not be zero\') else: M[r] = [scale*i for i in M[r]]
(3)某行乘以标量加到另一行
# r1 <--- r1 + r2*scale # 直接修改参数矩阵,无返回值 def addScaledRow(M, r1, r2, scale): if not scale: raise ValueError if (0 <= r1 < len(M)) and (0 <= r2 < len(M)): M[r1] = [M[r1][i] + scale * M[r2][i] for i in range(len(M[r2]))] else: raise IndexError(\'list index out of range\')
2.3 高斯消元法求解:Ax = b
(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)推演可逆矩阵
通过这段代码生成矩阵:
from helper import * A = generateMatrix(4,seed,singular=False) b = np.ones(shape=(4,1)) # it doesn\'t matter Ab = augmentMatrix(A.tolist(),b.tolist()) # please make sure you already correct implement augmentMatrix printInMatrixFormat(Ab,padding=4,truncating=0)
得到矩阵:
然后进行初等行变换:
(3)推演奇异矩阵
通过代码生成矩阵:
A = generateMatrix(4,seed,singular=True) b = np.ones(shape=(4,1)) # it doesn\'t matter Ab = augmentMatrix(A.tolist(),b.tolist()) # please make sure you already correct implement augmentMatrix printInMatrixFormat(Ab,padding=4,truncating=0)
得到矩阵:
然后进行初等行变换:
(4)高斯消去法的代码实现
我的low代码:
def gj_Solve(A, b, decPts=4, epsilon=1.0e-16): if len(A) != len(b): raise ValueError elif len(A) != len(A[0]): raise ValueError else: Ab = augmentMatrix(A, b) matxRound(Ab, decPts) num_row, num_clo = shape(Ab) for c in range(num_clo-1): current_max = 0.0 current_row = c for r in range(c, num_row): if abs(Ab[r][c]) > current_max: current_max = abs(Ab[r][c]) current_row = r if abs(current_max) < epsilon: return None else: swapRows(Ab, c, current_row) while abs((Ab[c][c]-1.0)) >= epsilon: scaleRow(Ab, c, 1.0 / Ab[c][c]) for j in range(c): while abs(Ab[j][c]) >= epsilon: addScaledRow(Ab, j, c, -Ab[j][c]) for j in range(c + 1, num_row): while abs(Ab[j][c]) >= epsilon: addScaledRow(Ab, j, c, -Ab[j][c]) res = [] for row in range(num_row): res.append([Ab[row][-1]]) return res
再看看参考答案的实现:
# 实现 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): raise ValueError Ab = augmentMatrix(A,b) for c in range(len(A[0])): AbT = transpose(Ab) col = AbT[c] maxValue = max(col[c:],key=abs) if abs(maxValue) < epsilon: return None maxIndex = col[c:].index(maxValue)+c swapRows(Ab,c,maxIndex) scaleRow(Ab,c,1.0/Ab[c][c]) for i in range(len(A)): if Ab[i][c] != 0 and i != c: addScaledRow(Ab,i,c,-Ab[i][c]) matxRound(Ab) return [[value] for value in transpose(Ab)[-1]
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()
得到样本点如图:
不断修改下面的m和b的值,拟合直线。这里我选去m=3.0, b=7.0
# 请选择最适合的直线 y = mx + b m = 3.0 b = 7.0 # 不要修改这里! plt.xlim((-5,5)) x_vals = plt.axes().get_xlim() y_vals = [m*x+b 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 计算平均平方误差 (MSE)
我们要编程计算所选直线的平均平方误差(MSE), 即数据集中每个点到直线的Y方向距离的平方的平均数,表达式如下:
代码实现:
# 实现以下函数并输出所选直线的MSE def calculateMSE(X,Y,m,b): if len(X) == len(Y) and len(X) != 0: n = len(X) square_li = [(Y[i]-m*X[i]-b)**2 for i in range(n)] return sum(square_li) / float(n) else: raise ValueError print(calculateMSE(X,Y,m,b))
得到的MSE是:1.7601561403444317。
3.3 的到最优参数
可以证明(此处不予证明)求解方程可以找到最优参数。其中向量Y,矩阵X和向量h分别为:
下面看下代码实现:
#实现线性回归 \'\'\' 参数:X, Y 返回:m,b \'\'\' def linearRegression(X, Y): X = [[x, 1] for x in X] Y = [[y] for y in Y] XT = transpose(X) A = matxMultiply(XT, X) b = matxMultiply(XT, Y) ret = gj_Solve(A, b) return ret[0][0], ret[1][0] m,b = linearRegression(X,Y) print(m,b)
# 3.2379 7.1899
最后我们看看得到的回归结果是什么,并用代码画出来:
x1,x2 = -5,5 y1,y2 = x1*m+b, x2*m+b 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.text(1,2,\'y = {m}x + {b}\'.format(m=m,b=b)) plt.show()
最后得到的直线是:
求得的回归结果对当前数据集的MSE是:
print(calculateMSE(X,Y,m,b))
# 1.3549197783872027
本篇就到这里,觉得还行记得点赞哦~~~
以上是关于机器学习之线性代数的主要内容,如果未能解决你的问题,请参考以下文章