Numpy之线性代数

作者: 陨星落云 | 来源:发表于2019-10-07 16:12 被阅读0次

矩阵

1. 矩阵初始化

Numpy函数库中存在两种不同的数据类型（矩阵matrix和数组array），都可以用于处理行列表示的数字元素，虽然他们看起来很相似，但是在这两个数据类型上执行相同的数学运算可能得到不同的结果，其中Numpy函数库中的matrix与MATLAB中matrices等价。

getA()是numpy的一个函数，作用是将矩阵转成一个ndarray，getA()函数和mat()函数的功能相反，是将一个矩阵转化为数组。

数组array 矩阵的初始化

import numpy as np

# 全零矩阵

myZero = np.zeros([3,3])
print(myZero)
print(type(myZero))

# 全一矩阵

myOnes = np.ones([3,3])
print(myOnes)
print(type(myOnes))

# 单位矩阵

myEye = np.eye(3)
print(myEye)
print(type(myEye))

# 对称矩阵

a1 = [1.,2.,3.]
myDiag = np.diag(a1)
print(myDiag)
print(type(myDiag))

# 随机矩阵

myRand = np.random.rand(3,3)
print(myRand)
print(type(myRand))

[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
<class 'numpy.ndarray'>
[[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
<class 'numpy.ndarray'>
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
<class 'numpy.ndarray'>
[[1. 0. 0.]
 [0. 2. 0.]
 [0. 0. 3.]]
<class 'numpy.ndarray'>
[[0.8047288  0.08855086 0.01365475]
 [0.27463797 0.61934569 0.59620009]
 [0.01570598 0.54358429 0.5640885 ]]
<class 'numpy.ndarray'>

矩阵matrix矩阵的初始化

通过mat（）函数将目标数据的类型转化成矩阵（matrix）

a1 = [1.,2.,3.]
myDiag = np.mat(np.diag(a1))
print(myDiag)
print(type(myDiag))

[[1. 0. 0.]
 [0. 2. 0.]
 [0. 0. 3.]]
<class 'numpy.matrix'>

直接创建

a = np.mat('1 2 3;4 5 6; 7 8 9')
print(a)
b = np.mat([[1,2,3],[4,5,6],[7,8,9]])
print(b)
c = np.mat(np.arange(9).reshape(3,3))
print(c)
print(type(c))

[[1 2 3]
 [4 5 6]
 [7 8 9]]
[[1 2 3]
 [4 5 6]
 [7 8 9]]
[[0 1 2]
 [3 4 5]
 [6 7 8]]
<class 'numpy.matrix'>

2. 矩阵元素运算

矩阵的元素运算是指矩阵在元素级别的加、减、乘、除运算。
（1）元素相加和相减。
条件：矩阵的行数和列数必须相同。
数学公式： $(A±B)_i,_j=A_i,_j±B_i,_j$

myOnes = np.ones([4,4])
myEye = np.eye(4)
print(myOnes+myEye)
print(myOnes-myEye)

[[2. 1. 1. 1.]
 [1. 2. 1. 1.]
 [1. 1. 2. 1.]
 [1. 1. 1. 2.]]
[[0. 1. 1. 1.]
 [1. 0. 1. 1.]
 [1. 1. 0. 1.]
 [1. 1. 1. 0.]]

（2）矩阵数乘：一个数乘以一个矩阵。

数学公式： $(cA)_i,_j=c·A_i,_j$

mymat = np.arange(9).reshape(3,3)
print(mymat)
print(10*mymat)

[[0 1 2]
 [3 4 5]
 [6 7 8]]
[[ 0 10 20]
 [30 40 50]
 [60 70 80]]

（3）矩阵所有元素求和。
数学公式： $其中\operatorname{sum}(A)=\sum_{i=1}^{m} \sum_{j=1}^{n} A_{i, j}=1 其中, 1<i<m, 1<j<n$

mymat = np.arange(9).reshape(3,3)
print(mymat)
print(np.sum(mymat))

[[0 1 2]
 [3 4 5]
 [6 7 8]]
36

（4）矩阵所有元素之积

数学公式： $其中\operatorname{product}(A)=\prod_{i=1}^{m} \prod_{j=1}^{n} A_{i, j}=1 其中, 1<i<m, 1<j<n$

mymat = np.arange(4).reshape(2,2)+1
print(mymat)
print(np.product(mymat))

[[1 2]
 [3 4]]
24

（5）矩阵各元素的n次幂：n=3。
数学公式： $A_{ij}^3=A_{ij}×A_{ij}×A_{ij}$

mymat = np.arange(9).reshape(3,3)
print(mymat)
print(np.power(mymat,3))

[[0 1 2]
 [3 4 5]
 [6 7 8]]
[[  0   1   8]
 [ 27  64 125]
 [216 343 512]]

3. 矩阵的乘法

（1）矩阵各元素的积：矩阵的点乘同维对应元素的相乘。当矩阵的维度不同时，会根据一定的广播规则将维数扩充到一致的形式。
数学公式： $(A×B)_i,_j =A_i,_j×B_i,_j$

mymat = np.arange(9).reshape(3,3)
print(mymat)
print(np.multiply(mymat,mymat))

[[0 1 2]
 [3 4 5]
 [6 7 8]]
[[ 0  1  4]
 [ 9 16 25]
 [36 49 64]]

（2）矩阵内积

数学公式： $[\boldsymbol{A}, \boldsymbol{B}]_{i, j}=\boldsymbol{A}_{i, 1} \boldsymbol{B}_{1, j}+\boldsymbol{A}_{i, 2} \boldsymbol{B}_{2, j}+\cdots+\boldsymbol{A}_{i n} \boldsymbol{B}_{n j}=\sum_{r=1}^{n} \boldsymbol{A}_{i, r} \boldsymbol{B}_{r, j^{\circ}}$

# 数组array
a = np.arange(9).reshape(3,3)
b = np.ones([3,3])
print(a)
print(b)
print(np.dot(a,b))
# 矩阵matrix
c = np.mat(a)
d = np.mat(b)
print(c*d)

[[0 1 2]
 [3 4 5]
 [6 7 8]]
[[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
[[ 3.  3.  3.]
 [12. 12. 12.]
 [21. 21. 21.]]
[[ 3.  3.  3.]
 [12. 12. 12.]
 [21. 21. 21.]]

（3）向量内积、外积

$x=\left(x_{1}, x_{1}, ...,x_{m}\right)$ $y=\left(y_{1}, y_{2},...,y_{n}\right)$

向量内积： $x^{T} y=\sum_{i=1}^{n} x_{i} y_{i}$

向量外积： $x y^{T}=\left(\begin{array}{ccc}{x_{1} y_{1}} & {\cdots} & {x_{1} y_{n}} \\ {\vdots} & {} & {\vdots} \\ {x_{m} y_{1}} & {\cdots} & {x_{m} y_{n}}\end{array}\right)$

a = [2,1,0]
b = [-1,2,1]
print(a)
print(b)
print(np.dot(a,b))
print(np.outer(a,b))

[2, 1, 0]
[-1, 2, 1]
内积：0
外积：
[[-2  4  2]
 [-1  2  1]
 [ 0  0  0]]

（4）向量叉乘(叉积)：运算结果是一个向量而不是一个标量。并且两个向量的叉积与这两个向量组成的坐标平面垂直。

数学公式： $a=\left(x_{1}, y_{1}, z_{1}\right)$ $b=\left(x_{2}, y_{2}, z_{2}\right)$
$a \times b=\left|\begin{array}{ccc}{\mathrm{i}} & {\mathrm{j}} & {\mathrm{k}} \\ {x_{1}} & {y_{1}} & {z_{1}} \\ {x_{2}} & {y_{2}} & {z_{2}}\end{array}\right|=\left(y_{1} z_{2}-y_{2} z_{1}\right) i-\left(x_{1} z_{2}-x_{2} z_{1}\right) j+\left(x_{1} y_{2}-x_{2} y_{1}\right) k$
其中
$i=(1,0,0) \quad j=(0,1,0) \quad \mathrm{k}=(0,0,1)$
根据i、j、k间关系，有:
$a \times b=\left(y_{1} z_{2}-y_{2} z_{1},-\left(x_{1} z_{2}-x_{2} z_{1}\right), x_{1} y_{2}-x_{2} y_{1}\right)$

例1、已知，a =（2，1，0），b =（-1，2，1），试求（1） $a\times b$ （2） $b\times a$
解:（1） $a\times b$ =（1，-2，5）：（2） $b\times a$ =（-1，2，5）

a = [2,1,0]
b = [-1,2,1]
print()
print(np.cross(a,b))
print(np.cross(b,a))

[ 1 -2  5]
[-1  2 -5]

4. 矩阵的转置

数学公式： $(A^T)_{ij} = A_{ij}$

a = np.arange(9).reshape(3,3)
print(a)
print(a.T)
print(a.transpose())

[[0 1 2]
 [3 4 5]
 [6 7 8]]
[[0 3 6]
 [1 4 7]
 [2 5 8]]
[[0 3 6]
 [1 4 7]
 [2 5 8]]

5. 矩阵对应列行的最大值，最小值，和

# 随机生成一个0-9的3×3矩阵
a = np.random.randint(0,10,(3,3))
print(a)
# 矩阵中的最大值
print(a.max())
# 计算矩阵中最大值的对应索引
print(np.argmax(a))
# 矩阵列中的最小值
print(a.min(axis=0))
# 计算矩阵列中最小值对应的列索引
print(np.argmin(a,0))
# 矩阵列求和
print(a.sum(axis=0))
# 矩阵行求和
print(a.sum(axis=1))

[[3 3 7]
 [9 8 5]
 [4 7 0]]
9
3
[3 3 0]
[0 0 2]
[16 18 12]
[13 22 11]

6. 矩阵的其他操作：行列数、切片、复制、比较

mymat = (np.arange(9)+1).reshape(3,3)
mymat1 = np.random.randint(0,10,(3,3))
print('mymat:',mymat)
print('mymat1:',mymat1)

[m,n] = np.shape(mymat)
print('矩阵的行列数:',m,n)
myscl1 = mymat[1]
print('按行切片：',myscl1)
myscl2 = mymat.T[1]
print('按列切片：',myscl2)
mymatcopy =mymat.copy()
print('复制矩阵:',mymatcopy)
# 比较
print('矩阵元素比较：\n',mymat<mymat1)

mymat: [[1 2 3]
 [4 5 6]
 [7 8 9]]
mymat1: [[1 3 7]
 [7 8 1]
 [6 2 3]]
矩阵的行列数: 3 3
按行切片： [4 5 6]
按列切片： [2 5 8]
复制矩阵: [[1 2 3]
 [4 5 6]
 [7 8 9]]
矩阵元素比较：
 [[False  True  True]
 [ True  True False]
 [False False False]]

7. 矩阵的行列式

import numpy as np
from numpy import linalg

A = np.full((4,4),3)
B = np.random.randint(0,10,(4,4))
print('A:\n',A)
print('B:\n',B)
print('det(A):',linalg.det(A))
print('det(B):',linalg.det(B))

A:
 [[3 3 3 3]
 [3 3 3 3]
 [3 3 3 3]
 [3 3 3 3]]
B:
 [[0 5 2 4]
 [2 2 4 2]
 [8 1 7 5]
 [4 0 4 5]]
det(A): 0.0
det(B): 197.9999999999998

8. 矩阵的逆和伪逆

矩阵的逆

import numpy as np
from numpy import linalg

B = np.random.randint(0,7,(3,3))
print('B:\n',B)
print('inv(B):',linalg.inv(B))

B:
 [[6 2 1]
 [2 3 6]
 [6 6 0]]
inv(B): [[ 0.24       -0.04       -0.06      ]
 [-0.24        0.04        0.22666667]
 [ 0.04        0.16       -0.09333333]]

注意：矩阵不满秩，则会报错

矩阵的伪逆

A = np.full((3,4),4)
print('A:\n',A)
print('A的伪逆：\n',linalg.pinv(A))

A:
 [[4 4 4 4]
 [4 4 4 4]
 [4 4 4 4]]
A的伪逆：
 [[0.02083333 0.02083333 0.02083333]
 [0.02083333 0.02083333 0.02083333]
 [0.02083333 0.02083333 0.02083333]
 [0.02083333 0.02083333 0.02083333]]

9. 矩阵的对称

mymat = np.arange(9).reshape(3,3)
print('mymat:\n',mymat)
print('矩阵的对称：\n',mymat*mymat.T)

mymat:
 [[0 1 2]
 [3 4 5]
 [6 7 8]]
矩阵的对称：
 [[ 0  3 12]
 [ 3 16 35]
 [12 35 64]]

10. 矩阵的秩

import numpy as np
from numpy import linalg

A = np.full((4,4),3)
B = np.random.randint(0,10,(4,4))
print('A:\n',A)
print('B:\n',B)
print('矩阵A的秩:',linalg.matrix_rank(A))
print('矩阵B的秩:',linalg.matrix_rank(B))

A:
 [[3 3 3 3]
 [3 3 3 3]
 [3 3 3 3]
 [3 3 3 3]]
B:
 [[7 1 2 6]
 [0 1 9 5]
 [7 7 4 0]
 [6 6 7 9]]
矩阵A的秩: 1
矩阵B的秩: 4

11. 可逆矩阵的求解

import numpy as np
from numpy import linalg


A = np.random.randint(0,10,(4,4))
b = np.array([1,10,0,1])
print('A:\n',A)
print('b:\n',b)
print('Ax=b,\nx=',linalg.solve(A,b.T))

A:
 [[3 9 2 7]
 [2 5 4 1]
 [8 5 3 8]
 [4 2 3 4]]
b:
 [ 1 10  0  1]
Ax=b,
x= [ 0.95638629  1.18691589  1.0623053  -2.09657321]

12. 矩阵的特征值与特征向量（EVD）

$Av=\lambda v$

这里 $A$ 是实对称矩阵， $v$ 是特征向量， $\lambda$ 是特征值。下面我们使用Numpy求取矩阵的特征值和特征向量。

import numpy as np
from numpy import linalg

A = np.mat(np.array([[8,1,6],[3,5,7],[4,9,2]]))
evals,evecs = linalg.eig(A)
print('特征值：\n',evals,'\n特征向量：\n',evecs)

特征值：
 [15.          4.89897949 -4.89897949] 
特征向量：
 [[-0.57735027 -0.81305253 -0.34164801]
 [-0.57735027  0.47140452 -0.47140452]
 [-0.57735027  0.34164801  0.81305253]]

有了特征值和特征向量，可以还原出原矩阵(公式1-1)：
$A=Q \Sigma Q^{-1}$
A是实对称矩阵，Q是特征向量， $\Sigma$ 是特征值构成的对角矩阵。下面是代码实现：

sigma =evals*np.eye(3)
print(evecs*sigma*linalg.inv(evecs))

[[8. 1. 6.]
 [3. 5. 7.]
 [4. 9. 2.]]

13. 奇异值分解（SVD）

有一个m×n的实数矩阵A，我们想要把它分解成如下的形式(公式2-1)：
$A=U \Sigma V^{T}$
其中 $U$ 和 $V$ 均为单位正交阵，即有 $UU^T=I$ 和 $VV^T=I$ ， $U$ 称为左奇异矩阵， $V$ 称为右奇异矩阵， $Σ$ 仅在主对角线上有值，我们称它为奇异值，其它元素均为0。上面矩阵的维度分别为 $U∈R^{m×m}, Σ∈R^{m×n}, V∈R^{n×n}$ 。

一般地 $Σ$ 有如下形式:
$\Sigma=\left[\begin{array}{ccccc}{\sigma_{1}} & {0} & {0} & {0} & {0} \\ {0} & {\sigma_{2}} & {0} & {0} & {0} \\ {0} & {0} & {\ddots} & {0} & {0} \\ {0} & {0} & {0} & {\ddots} & {0}\end{array}\right]_{m \times n}$

奇异值分解.jpg

对于奇异值分解，我们可以利用上面的图形象表示，图中方块的颜色表示值的大小，颜色越浅，值越大。对于奇异值矩阵 $Σ$ ，只有其主对角线有奇异值，其余均为0。

正常求上面的 $U,V,Σ$ 不便于求，我们可以利用如下性质(公式2-2、2-3):
$\begin{array}{l}{A A^{T}=U \Sigma V^{T} V \Sigma^{T} U^{T}=U \Sigma \Sigma^{T} U^{T}} \\ {A^{T} A=V \Sigma^{T} U^{T} U \Sigma V^{T}=V \Sigma^{T} \Sigma V^{T}}\end{array}$
注：需要指出的是，这里 $ΣΣ^T$ 与 $Σ^TΣ$ 在矩阵的角度上来讲，它们是不相等的，因为它们的维数不同 $ΣΣ^T∈R^{m×m}$ ，而 $Σ^TΣ∈R^{n×n}$ ，但是它们在主对角线的奇异值是相等的，即有
$\Sigma \Sigma^{T}=\left[\begin{array}{cccc}{\sigma_{1}^{2}} & {0} & {0} & {0} \\ {0} & {\sigma_{2}^{2}} & {0} & {0} \\ {0} & {0} & {\ddots} & {0} \\ {0} & {0} & {0} & {\ddots}\end{array}\right] \quad \Sigma^{T} \Sigma=\left[\begin{array}{cccc}{\sigma_{1}^{2}} & {0} & {0} & {0} \\ {0} & {\sigma_{2}^{2}} & {0} & {0} \\ {0} & {0} & {\ddots} & {0} \\ {0} & {0} & {0} & {\ddots}\end{array}\right]_{n \times n}$
可以看到式（2-2）与式（1-1）的形式非常相同，进一步分析，我们可以发现 $AA^T$ 和 $A^TA$ 也是对称矩阵，那么可以利用式（1-1），做特征值分解。利用式（2-2）特征值分解，得到的特征矩阵即为U；利用式（2-3）特征值分解，得到的特征矩阵即为V；对 $ΣΣ^T$ 或 $Σ^TΣ$ 中的特征值开方，可以得到所有的奇异值。

例：对A进行奇异值分解

$A=\left[\begin{array}{ccccc}{1} & {5} & {7} & {6} & {1} \\ {2} & {1} & {10} & {4} & {4} \\ {3} & {6} & {7} & {5} & {2}\end{array}\right]$

A = np.mat([[1,5,7,6,1],
            [2,1,10,4,4],
            [3,6,7,5,2]])

print('A:\n',A)

U,evals,VT = linalg.svd(A)
print('U=\n',U,'\nevals=\n',evals,'\nVT=\n',VT)
sigma = np.mat([[evals[0],0,0,0,0],
                [0,evals[1],0,0,0],
                [0,0,evals[2],0,0]])
#print(sigma)
print(U*sigma*VT)

A:
 [[ 1  5  7  6  1]
 [ 2  1 10  4  4]
 [ 3  6  7  5  2]]
U=
 [[-0.55572489  0.40548161 -0.72577856]
 [-0.59283199 -0.80531618  0.00401031]
 [-0.58285511  0.43249337  0.68791671]] 
evals=
 [18.53581747  5.0056557   1.83490648] 
VT=
 [[-0.18828164 -0.37055755 -0.74981208 -0.46504304 -0.22080294]
 [ 0.01844501  0.76254787 -0.4369731   0.27450785 -0.38971845]
 [ 0.73354812  0.27392013 -0.12258381 -0.48996859  0.36301365]
 [ 0.36052404 -0.34595041 -0.43411102  0.6833004   0.30820273]
 [-0.5441869   0.2940985  -0.20822387 -0.0375734   0.7567019 ]]
[[18.53581747  0.          0.          0.          0.        ]
 [ 0.          5.0056557   0.          0.          0.        ]
 [ 0.          0.          1.83490648  0.          0.        ]]
[[ 1.  5.  7.  6.  1.]
 [ 2.  1. 10.  4.  4.]
 [ 3.  6.  7.  5.  2.]]

参考资料：

参考网站：

https://www.jianshu.com/p/45417c05574c

https://www.cnblogs.com/wj-1314/p/10244807.html

https://www.cnblogs.com/endlesscoding/p/10033527.html

参考书籍：《python科学计算》《机器学习算法原理与编程实实践》

Numpy之线性代数

矩阵

1. 矩阵初始化

2. 矩阵元素运算

3. 矩阵的乘法

4. 矩阵的转置

5. 矩阵对应列行的最大值，最小值，和

6. 矩阵的其他操作：行列数、切片、复制、比较

7. 矩阵的行列式

8. 矩阵的逆和伪逆

9. 矩阵的对称

10. 矩阵的秩

11. 可逆矩阵的求解

12. 矩阵的特征值与特征向量（EVD）

13. 奇异值分解（SVD）

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

数字图像处理与计算机视觉（python）

numpy模块

python数学应用

机器学习爱好者

线性代数

Numpy

Python数据分析三剑客：Pandas、Numpy、Matplotlib

机器学习