Numpy之线性代数

矩阵

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}
其中UV均为单位正交阵,即有UU^T=IVV^T=IU称为左奇异矩阵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^TA^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.]]

参考资料:

参考网站:

//www.greatytc.com/p/45417c05574c

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

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

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,743评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,296评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,285评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,485评论 1 283
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,581评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,821评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,960评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,719评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,186评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,516评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,650评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,329评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,936评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,757评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,991评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,370评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,527评论 2 349

推荐阅读更多精彩内容