标签：learning div 优雅 3.3 表示 .com 项目 pytho false
题目是优达学城机器学习入门线性代数作业。下面是我的实现。
工具为jupyter notebook，不用该工具请自行导入相关依赖。
完整内容已上传到github：https://github.com/ZingP/machine-learning/tree/master/linear_algebra
本篇代码中引用的helper.py可到上面github上下载。
在创建矩阵之前注意选择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]]
def shape(M):
"""返回矩阵的行列"""
return len(M), len(M[0])
# 每个元素四舍五入到特定小数数位
# 直接修改参数矩阵，无返回值
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)
def transpose(M):
# *M 分解出列表中的子元素（子列表）
# zip()将子列表中对应的元素打包成元组，返回包含一个个元组的列表
# 然后用列表推导式...真优雅啊
return [list(col) for col in zip(*M)]
# 计算矩阵乘法 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.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)]
（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‘)
（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]
用代码生成随机样本点：
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()
得到直线如下图：
我们要编程计算所选直线的平均平方误差(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。
可以证明（此处不予证明）求解方程
可以找到最优参数。其中向量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
本篇就到这里，觉得还行记得点赞哦~~~
标签：learning div 优雅 3.3 表示 .com 项目 pytho false
原文地址：http://www.cnblogs.com/zingp/p/8011295.html
