机器学习之线性代数

说明

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

  工具为jupyter notebook,不用该工具请自行导入相关依赖。

  完整内容已上传到github:https://github.com/ZingP/machine-learning/tree/master/linear_algebra

  本篇代码中引用的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

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

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏漫漫深度学习路

tensorflow学习笔记(三十):tf.gradients 与 tf.stop_gradient() 与 高阶导数

gradient tensorflow中有一个计算梯度的函数tf.gradients(ys, xs),要注意的是,xs中的x必须要与ys相关,不相关的话,会报错...

1.4K9
来自专栏利炳根的专栏

学习笔记CB013: TensorFlow、TensorBoard、seq2seq

tensorflow基于图结构深度学习框架,内部通过session实现图和计算内核交互。

4397
来自专栏决胜机器学习

有趣的算法(一)——n阶层尾部有几个0

有趣的算法(一)——n阶层尾部有几个0 (原创内容,转载请注明来源,谢谢) 最近在网上看到好几次这个题目,觉得挺有意思,则准备用PHP进行实现。 1、题目 ...

3596
来自专栏祥子的故事

tensorflow | 重新学习 | 了解graph 和 Session

3648
来自专栏月色的自留地

《连连看》算法c语言演示(自动连连看)

2619
来自专栏Bay的专栏

如何使用 scikit-learn 为机器学习准备文本数据

文本数据需要特殊处理,然后才能开始将其用于预测建模。

2265
来自专栏技术沉淀

命令行工具:cut列操作

1977
来自专栏鸿的学习笔记

写给开发者的机器学习指南(十三)

在我们实际使用支持向量机(SVM)之前,我先简要介绍一下SVM是什么。 基本SVM是一个二元分类器,它通过选取代表数据点之间最大间隔的超平面将数据集分成2部分。...

611
来自专栏xingoo, 一个梦想做发明家的程序员

动态规划基本要素

动态规划性质: 1  最优子结构性质  2 子问题重叠性质 ----->该问题可用动态规划算法求解的基本要素 1 最优子结构 当问题的最优解包含了其子问题的最优...

20710
来自专栏人工智能LeadAI

使用TensorFlow实现手写识别(Softmax)

准备工作 由于将TensorFlow安装到了Conda的tensorflow环境,虽然可以用Jupyter notebook打开,但是没有提示,写代码不方便,所...

3885

扫码关注云+社区