讲解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))