• Python--线性代数篇


    讲解Python在线性代数中的应用,包括:

    一、矩阵创建

    先导入Numpy模块,在下文中均采用np代替numpy

    1 import numpy as np

    矩阵创建有两种方法,一是使用np.mat函数或者np.matrix函数,二是使用数组代替矩阵,实际上官方文档建议我们使用二维数组代替矩阵来进行矩阵运算;因为二维数组用得较多,而且基本可取代矩阵。

     1 >>> a = np.mat([[1, 2, 3], [4, 5, 6]])   #使用mat函数创建一个2X3矩阵
     2 >>> a
     3 matrix([[1, 2, 3],
     4         [4, 5, 6]])
     5 >>> b = np.matrix([[1, 2, 3], [4, 5, 6]])#np.mat和np.matrix等价
     6 >>> b
     7 matrix([[1, 2, 3],
     8         [4, 5, 6]])
     9 >>> a.shape     #使用shape属性可以获取矩阵的大小
    10 (2, 3)
    1 >>> c = np.array([[1, 2, 3], [4, 5, 6]]) #使用二维数组代替矩阵,常见的操作通用
    2 >>> c#注意c是array类型,而a是matrix类型
    3 array([[1, 2, 3],
    4        [4, 5, 6]])

    单位阵的创建

    1 >>> I = np.eye(3)
    2 >>> I
    3 array([[ 1.,  0.,  0.],
    4        [ 0.,  1.,  0.],
    5        [ 0.,  0.,  1.]])

    矩阵元素的存取操作:

    1 >>> a[0]#获取矩阵的某一行
    2 matrix([[1, 2, 3]])
    3 >>> a[:, 0].reshape(-1, 1)#获取矩阵的某一列
    4 matrix([[1],
    5         [4]])
    6 >>> a[0, 1]#获取矩阵某个元素
    7 2

    二、矩阵乘法和加法

    矩阵类型,在满足乘法规则的条件下可以直接相乘

     1 >>> A = np.mat([[1, 2, 3], [3, 4, 5], [6, 7, 8]])#使用mat函数
     2 >>> B = np.mat([[5, 4, 2], [1, 7, 9], [0, 4, 5]])
     3 >>> A   #注意A, B都是matrix类型,可以使用乘号,如果是array则不可以直接使用乘号
     4 matrix([[1, 2, 3],
     5         [3, 4, 5],
     6         [6, 7, 8]])
     7 >>> B
     8 matrix([[5, 4, 2],
     9         [1, 7, 9],
    10         [0, 4, 5]])
    11 >>> A * B#学过线性代数的都知道:A * B != B * A
    12 matrix([[  7,  30,  35],
    13         [ 19,  60,  67],
    14         [ 37, 105, 115]])
    15 >>> B * A
    16 matrix([[ 29,  40,  51],
    17         [ 76,  93, 110],
    18         [ 42,  51,  60]])

    如果是使用数组代替矩阵进行运算则不可以直接使用乘号,应使用dot()函数。dot函数用于矩阵乘法,对于二维数组,它计算的是矩阵乘积,对于一维数组,它计算的是内积。

     1 >>> C = np.array([[1, 2, 3], [3, 4, 5], [6, 7, 8]])
     2 >>> D = np.array([[5, 4, 2], [1, 7, 9], [0, 4, 5]])
     3 >>> C          #C, D都是array类型,不能直接使用乘号,应该使用dot()函数
     4 array([[1, 2, 3],
     5        [3, 4, 5],
     6        [6, 7, 8]])
     7 >>> D
     8 array([[5, 4, 2],
     9        [1, 7, 9],
    10        [0, 4, 5]])
    11 #>>> C * D, Error, 注意这不是矩阵乘法!!!
    12 >>> np.dot(C, D)#正确的写法,得到的结果和上一段代码的第11行的结果的一样的。
    13 array([[  7,  30,  35],
    14        [ 19,  60,  67],
    15        [ 37, 105, 115]])

    如何理解对于一维数组,它计算的是内积???

    注意:在线性代数里面讲的维数和数组的维数不同,如线代中提到的n维行向量在Python中是一维数组,而线代中的n维列向量在Python中是一个shape为(n, 1)的二维数组!

    第16行,第18行:F是一维数组,G是二维数组,维数不同,个人认为相乘没有意义,但是16行没有错误,18行报错。关于dot()的乘法规则见:NumPy-快速处理数据--矩阵运算

     1 >>> E = np.array([1, 2, 3])
     2 >>> F = np.array([4, 3, 9])
     3 >>> E.shape#E,F都是一维数组
     4 (3,)
     5 >>> np.dot(E, F)
     6 37
     7 >>> np.dot(F, E)
     8 37
     9 >>> G = np.array([4, 3, 9]).reshape(-1, 1)
    10 >>> G
    11 array([[4],
    12        [3],
    13        [9]])
    14 >>> G.shape
    15 (3, 1)
    16 >>> np.dot(F, G)#因此dot(F, G)不再是内积,而是一个只有一个元素的数组
    17 array([106])
    18 >>> np.dot(G, F)#ValueError: shapes (3,1) and (3,) not aligned: 1 (dim 1) != 3 (dim 0)
    19 >>> E.shape = (1, -1)#把E改为二维数组
    20 >>> E
    21 array([[1, 2, 3]])
    22 >>> E.shape
    23 (1, 3)
    24 >>> np.dot(G, E)#3×1的G向量乘以1×3的E向量会得到3×3的矩阵
    25 array([[ 4,  8, 12],
    26        [ 3,  6,  9],
    27        [ 9, 18, 27]])

    矩阵的加法运算

    1 >>> A + B#矩阵的加法对matrix类型和array类型是通用的
    2 matrix([[ 6,  6,  5],
    3         [ 4, 11, 14],
    4         [ 6, 11, 13]])
    5 >>> C + D
    6 array([[ 6,  6,  5],
    7        [ 4, 11, 14],
    8        [ 6, 11, 13]])

    矩阵的数乘运算

    1 >>> 2 * A#矩阵的数乘对matrix类型和array类型是通用的
    2 matrix([[ 2,  4,  6],
    3         [ 6,  8, 10],
    4         [12, 14, 16]])
    5 >>> 2 * C
    6 array([[ 2,  4,  6],
    7        [ 6,  8, 10],
    8        [12, 14, 16]])

    三、矩阵的转置

     1 >>> A = np.array([[1, 2, 3], [3, 4, 5], [6, 7, 8]])
     2 >>> B = np.array([[5, 4, 2], [1, 7, 9], [0, 4, 5]])
     3 >>> A
     4 array([[1, 2, 3],
     5        [3, 4, 5],
     6        [6, 7, 8]])
     7 >>> A.T  #A的转置
     8 array([[1, 3, 6],
     9        [2, 4, 7],
    10        [3, 5, 8]])
    11 >>> A.T.T#A的转置的转置还是A本身
    12 array([[1, 2, 3],
    13        [3, 4, 5],
    14        [6, 7, 8]])

    验证矩阵转置的性质:(A±B)'=A'±B'

    1 >>> (A + B).T
    2 array([[ 6,  4,  6],
    3        [ 6, 11, 11],
    4        [ 5, 14, 13]])
    5 >>> A.T + B.T
    6 array([[ 6,  4,  6],
    7        [ 6, 11, 11],
    8        [ 5, 14, 13]])

    验证矩阵转置的性质:(KA)'=KA'

    1 >>> 10 * (A.T)
    2 array([[10, 30, 60],
    3        [20, 40, 70],
    4        [30, 50, 80]])
    5 >>> (10 * A).T
    6 array([[10, 30, 60],
    7        [20, 40, 70],
    8        [30, 50, 80]])

    验证矩阵转置的性质:(A×B)'= B'×A'

    1 >>> np.dot(A, B).T
    2 array([[  7,  19,  37],
    3        [ 30,  60, 105],
    4        [ 35,  67, 115]])
    5 >>> np.dot(B.T, A.T)
    6 array([[  7,  19,  37],
    7        [ 30,  60, 105],
    8        [ 35,  67, 115]])

    四、方阵的迹

    方阵的迹就是主对角元素之和,使用trace()函数获得方阵的迹:

     1 >>> A
     2 array([[1, 2, 3],
     3        [3, 4, 5],
     4        [6, 7, 8]])
     5 >>> B
     6 array([[5, 4, 2],
     7        [1, 7, 9],
     8        [0, 4, 5]])
     9 >>> np.trace(A)  # A的迹等于A.T的迹
    10 13
    11 >>> np.trace(A.T)
    12 13
    13 >>> np.trace(A+B)# 和的迹 等于 迹的和
    14 30
    15 >>> np.trace(A) + np.trace(B)
    16 30

    五、计算行列式

    1 >>> A
    2 array([[1, 2],
    3        [1, 3]])
    4 >>> np.linalg.det(A)
    5 1.0

    六、逆矩阵/伴随矩阵

    若A存在逆矩阵(满足det(A) != 0,或者A满秩),使用linalg.inv求得方阵A的逆矩阵

     1 import numpy as np
     2 >>> A = np.array([[1, -2, 1], [0, 2, -1], [1, 1, -2]])
     3 >>> A
     4 array([[ 1, -2,  1],
     5        [ 0,  2, -1],
     6        [ 1,  1, -2]])
     7 >>> A_det = np.linalg.det(A)      #求A的行列式,不为零则存在逆矩阵
     8 >>> A_det
     9 -3.0000000000000004
    10 >>> A_inverse = np.linalg.inv(A)  #求A的逆矩阵
    11 >>> A_inverse
    12 array([[ 1.        ,  1.        ,  0.        ],
    13        [ 0.33333333,  1.        , -0.33333333],
    14        [ 0.66666667,  1.        , -0.66666667]])
    15 >>> np.dot(A, A_inverse)          #A与其逆矩阵的乘积为单位阵
    16 array([[ 1.,  0.,  0.],
    17        [ 0.,  1.,  0.],
    18        [ 0.,  0.,  1.]])
    19 >>> A_companion = A_inverse * A_det  #求A的伴随矩阵
    20 >>> A_companion
    21 array([[-3., -3., -0.],
    22        [-1., -3.,  1.],
    23        [-2., -3.,  2.]])

    七、解一元线性方程

    使用np.linalg.solve()解一元线性方程组,待解方程为:

     x + 2y +  z = 7
    2x -  y + 3z = 7
    3x +  y + 2z =18
     1 >>> import numpy as np
     2 >>> A = np.array([[1, 2, 1], [2, -1, 3], [3, 1, 2]])
     3 >>> A    #系数矩阵
     4 array([[ 1,  2,  1],
     5        [ 2, -1,  3],
     6        [ 3,  1,  2]])
     7 >>> B = np.array([7, 7, 18])
     8 >>> B
     9 array([ 7,  7, 18])
    10 >>> x = np.linalg.solve(A, B)
    11 >>> x
    12 array([ 7.,  1., -2.])
    13 >>> np.dot(A, x)#检验正确性,结果为B
    14 array([  7.,   7.,  18.])

    使用np.allclose()检测两个矩阵是否相同:

    1 >>> np.allclose(np.dot(A, x), B)#检验正确性
    2 True

    使用 help(np.allclose) 查看 allclose() 的用法:

    allclose(a, b, rtol=1e-05, atol=1e-08)
        Parameters
        ----------
        a, b : array_like
            Input arrays to compare.
        rtol : float
            The relative tolerance parameter (see Notes).
        atol : float
            The absolute tolerance parameter (see Notes).
        
        Returns
        -------
        allclose : bool
            Returns True if the two arrays are equal within the given
            tolerance; False otherwise.

    八、计算矩阵距离

    矩阵的距离,这里是的是欧几里得距离,其他距离表示方法我们以后再谈,这里说一下如何计算两个形状相同矩阵之间的距离。

     1 >>> A = np.array([[0, 1], [1, 0]])#先创建两个矩阵
     2 >>> B = np.array([[1, 1], [1, 1]])
     3 >>> C = A - B       #计算距离矩阵C
     4 >>> C
     5 array([[-1,  0],
     6        [ 0, -1]])
     7 >>> D = np.dot(C, C)#距离矩阵的平方
     8 >>> E = np.trace(D) #计算矩阵D的迹
     9 >>> E
    10 2
    11 >>> E ** 0.5        #将E开平方得到距离
    12 1.4142135623730951

    关于计算矩阵距离我也不理解。网上看的帖子,先记下来

    九、矩阵的秩

    numpy包中的linalg.matrix_rank方法计算矩阵的秩:

     1 >>> import numpy as np
     2 >>> I = np.eye(3)#先创建一个单位阵
     3 >>> I
     4 array([[ 1.,  0.,  0.],
     5        [ 0.,  1.,  0.],
     6        [ 0.,  0.,  1.]])
     7 >>> np.linalg.matrix_rank(I)#
     8 3
     9 >>> I[1, 1] = 0#将该元素置为0
    10 >>> I
    11 array([[ 1.,  0.,  0.],
    12        [ 0.,  0.,  0.],
    13        [ 0.,  0.,  1.]])
    14 >>> np.linalg.matrix_rank(I)#此时秩变成2
    15 2

    十、求方阵的特征值特征向量

     1 >>> import numpy as np
     2 >>> x = np.diag((1, 2, 3))#创建一个对角矩阵!
     3 >>> x
     4 array([[1, 0, 0],
     5        [0, 2, 0],
     6        [0, 0, 3]])
     7 >>> a,b = np.linalg.eig(x)#特征值保存在a中,特征向量保存在b中
     8 >>> a
     9 array([ 1.,  2.,  3.])
    10 >>> b
    11 array([[ 1.,  0.,  0.],
    12        [ 0.,  1.,  0.],
    13        [ 0.,  0.,  1.]])

    根据公式 Ax = λx 检验特征值与特征向量是否正确:

     1 for i in range(3):#方法一
     2     if np.allclose(np.dot(a[i], b[:, i]), x[:, i]):#np.allclose()方法在第七节提到过
     3         print 'Right'
     4     else:
     5         print 'Error'
     6 
     7 for i in range(3):#方法二
     8     if (np.dot(a[i], b[:, i]) == x[:, i]).all():
     9         print 'Right'
    10     else:
    11         print 'Error'

    注意,如果写成 if np.dot(a[i], b[:, i]) == x[:, i]: 是错误的:(矩阵包含有多个值,应该使用a.any()或者a.all()判断)

     ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() 

    十一、判断正定矩阵

    设M是n阶方阵,如果对任何非零向量z,都有z'Mz> 0,其中z' 表示z的转置,就称M正定矩阵。

    判定定理1:对称阵A为正定的充分必要条件是:A的特征值全为正。
    判定定理2:对称阵A为正定的充分必要条件是:A的各阶顺序主子式都为正。
    判定定理3:任意阵A为正定的充分必要条件是:A合同于单位阵。

    下面用定理1判断对称阵是否为正定阵

     1 >>> import numpy as np
     2 >>> A = np.arange(16).reshape(4, 4)
     3 >>> A
     4 array([[ 0,  1,  2,  3],
     5        [ 4,  5,  6,  7],
     6        [ 8,  9, 10, 11],
     7        [12, 13, 14, 15]])
     8 >>> A = A + A.T             #将方阵转换成对称阵
     9 >>> A
    10 array([[ 0,  5, 10, 15],
    11        [ 5, 10, 15, 20],
    12        [10, 15, 20, 25],
    13        [15, 20, 25, 30]])
    14 >>> B = np.linalg.eigvals(A)#求B的特征值,注意:eig()是求特征值特征向量
    15 >>> B
    16 array([  6.74165739e+01 +0.00000000e+00j,
    17         -7.41657387e+00 +0.00000000e+00j,
    18          2.04219701e-15 +3.94306094e-15j,
    19          2.04219701e-15 -3.94306094e-15j])
    20 
    21 if np.all(B>0):             #判断是不是所有的特征值都大于0,用到了all函数,显然对称阵A不是正定的
    22     print 'Yes'

    创建一个对角元素都为正的对角阵,它一定是正定的:

    1 >>> A = np.diag((1, 2, 3))#创建对角阵,其特征值都为正
    2 >>> B = np.linalg.eigvals(A)#求特征值
    3 >>> B
    4 array([ 1.,  2.,  3.])
    5 >>> if np.all(B>0):#判断特征值是否都大于0
    6     print 'Yes'

    网上查到更简便的方法是对对称阵进行cholesky分解,如果像这样没有提示出错,就说明它是正定的。如果提示出错,就说明它不是正定矩阵,你可以使用try函数捕获错误值:

     1 # -*- coding: utf-8 -*-
     2 import numpy as np
     3 
     4 A = np.arange(16).reshape(4, 4)
     5 A = A + A.T
     6 print A
     7 try:
     8     B = np.linalg.cholesky(A)
     9 except :
    10     print ('不是正定矩阵,不能进行cholesky分解。')

    当不能进行cholesky分解时,出现的异常是: LinAlgError: Matrix is not positive definite ,但是但是LinAlgError不是Python标准异常,因此不能使用这条语句。

    1 except LinAlgError as reason:
    2     print ('不是正定矩阵,不能进行cholesky分解。
    出错原因是:' + str(reason))
    http://www.cnblogs.com/moon1992/
  • 相关阅读:
    [RK3288][Android6.0] U-boot 启动流程小结【转】
    学习笔记二十三——字符函数库cctype【转】
    【Git学习笔记】用git pull取回远程仓库某个分支的更新,再与本地的指定分支自动merge【转】
    Git 少用 Pull 多用 Fetch 和 Merge 【已翻译100%】【转】
    git 拉取和获取 pull 和 fetch 区别【转】
    setprecision、fixed、showpoint的用法总结(经典!!超经典!!)【转】
    Android休眠唤醒机制简介(二)
    获取元素个数的函数
    返回两个时间范围内的一个随机时间
    全角半角转换函数
  • 原文地址:https://www.cnblogs.com/moon1992/p/4960700.html
Copyright © 2020-2023  润新知