机器学习之线性代数

说明

　　题目是优达学城机器学习入门线性代数作业。下面是我的实现。

　　工具为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

本篇就到这里，觉得还行记得点赞哦~~~