說明
題目是優達學城機器學習入門線性代數作業。下面是我的實現。
工具為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
本篇就到這里，覺得還行記得點贊哦~~~