04_numpy学习笔记（下）：线性代数
文章目录
数据科学中通常会有大量计算，如果仅有标量积算，是远远不能满足计算需求的，而矩阵计算刚好能解决这个问题。Numpy 定义了
matrix 类型，使用该 matrix类型创建的是矩阵对象，它们的加减乘除运算缺省采用矩阵方式计算，因此用法和Matlab十分类似。但是由于 NumPy 中同时存在 ndarray 和 matrix对象，因此用户很容易将两者弄混。这有违 Python 的“显式优于隐式”的原则，因此官方并不推荐在程序中使用 matrix 。
一、矩阵和向量积
numpy.dot(a, b[, out])计算两个矩阵的乘积，如果是一维数组则是它们的内积。
import numpy as np
x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 3, 4, 5, 6])
z = np.dot(x, y)
print(z) # 70
x = np.array([[1, 2, 3], [3, 4, 5], [6, 7, 8]])
print(x)
# [[1 2 3]
# [3 4 5]
# [6 7 8]]
y = np.array([[5, 4, 2], [1, 7, 9], [0, 4, 5]])
print(y)
# [[5 4 2]
# [1 7 9]
# [0 4 5]]
z = np.dot(x, y)
print(z)
# [[ 7 30 35]
# [ 19 60 67]
# [ 37 105 115]]
z = np.dot(y, x)
print(z)
# [[ 29 40 51]
# [ 76 93 110]
# [ 42 51 60]]
注意：在线性代数里面讲的维数和数组的维数不同，如线代中提到的n维行向量在 Numpy 中是一维数组，而线性代数中的n维列向量在 Numpy 中是一个shape为(n, 1)的二维数组。
二、矩阵特征值与特征向量
numpy.linalg.eig(a). 计算方阵的特征值和特征向量。numpy.linalg.eigvals(a). 计算方阵的特征值。
（1）求方阵的特征值和特征向量
import numpy as np
# 创建一个对角矩阵！
x = np.diag((1, 2, 3))
print(x)
# [[1 0 0]
# [0 2 0]
# [0 0 3]]
print(np.linalg.eigvals(x))
# [1. 2. 3.]
a, b = np.linalg.eig(x)
# 特征值保存在a中，特征向量保存在b中
print(a)
# [1. 2. 3.]
print(b)
# [[1. 0. 0.]
# [0. 1. 0.]
# [0. 0. 1.]]
# 检验特征值与特征向量是否正确
for i in range(3):
if np.allclose(a[i] * b[:, i], np.dot(x, b[:, i])):
print('Right')
else:
print('Error')
# Right
# Right
# Right
（2）判断对称阵是否为正定阵（特征值是否全部为正）。
import numpy as np
A = np.arange(16).reshape(4, 4)
print(A)
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]
# [12 13 14 15]]
A = A + A.T # 将方阵转换成对称阵
print(A)
# [[ 0 5 10 15]
# [ 5 10 15 20]
# [10 15 20 25]
# [15 20 25 30]]
B = np.linalg.eigvals(A) # 求A的特征值
print(B)
# [ 6.74165739e+01 ‐7.41657387e+00 1.82694656e‐15 ‐1.72637110e‐15]
# 判断是不是所有的特征值都大于0，用到了all函数，显然对称阵A不是正定的
if np.all(B > 0):
print('Yes')
else:
print('No')
# No
三、矩阵分解
1. 奇异值分解
SVD（Singular Value Decomposition）可以理解为：将一个比较复杂的矩阵用更小更简单的3个子矩阵的相乘来表示，这3个小矩阵描述了大矩阵重要的特性。
定义：矩阵的奇异值分解是指将一个秩为的实矩阵分解为三个实矩阵乘积的形式：
A
m
×
n
=
U
m
×
m
Σ
m
×
n
V
n
×
n
T
≈
U
m
×
k
Σ
k
×
k
V
k
×
n
T
A_{m\times n} = U_{m \times m} \Sigma_{m \times n} V ^ { T }_{n\times n} \approx U_{m \times k} \Sigma_{k \times k} V ^ { T }_{k\times n}
Am×n=Um×mΣm×nVn×nT≈Um×kΣk×kVk×nT
其中
U
U
U 是
m
m
m阶正交矩阵（
U
U
U的列向量称为左奇异向量），
V
V
V是
n
n
n阶正交矩阵**（
V
V
V的列向量称为右奇异向量），
∑
\sum
∑是矩形对角矩阵，称为奇异值矩阵**，对角线上的元素称为奇异值。
Σ
=
[
D
r
×
r
0
0
0
]
m
×
n
\Sigma = \begin{bmatrix} D_{r\times r}&0\\ 0&0\\ \end{bmatrix}_{m\times n}
Σ=[Dr×r000]m×n
D
D
D是一个
r
×
r
r \times r
r×r的对角阵，
D
D
D的对角线元素是
A
A
A的前
r
r
r个奇异值
σ
1
≥
σ
2
≥
⋯
≥
σ
r
>
0
\sigma _ { 1 } \geq \sigma _ { 2 } \geq \cdots \geq \sigma _ { r } > 0
σ1≥σ2≥⋯≥σr>0（非负，降序）。
知识点：任意一个实矩阵可以由其外积展开式表示
A
=
σ
1
u
1
v
1
T
+
σ
2
u
2
v
2
T
+
⋯
+
σ
r
u
r
v
r
T
A = \sigma _ { 1 } u _ { 1 } v _ { 1 } ^ { T } + \sigma _ { 2 } u _ { 2 } v _ { 2 } ^ { T } + \cdots + \sigma _ { r } u _ { r } v _ { r } ^ { T }
A=σ1u1v1T+σ2u2v2T+⋯+σrurvrT
其中 u k v k T u _ { k } v _ { k } ^ { T } ukvkT为矩阵，是列向量 u k u _ { k } uk和行向量 v k T v _ { k } ^ { T } vkT的外积， σ k \sigma_k σk为奇异值， u k , v k T , σ k u _ { k } , v _ { k } ^ { T } , \sigma _ { k } uk,vkT,σk通过矩阵 A A A的奇异值分解得到。
知识点：奇异值在矩阵中按照从大到小排列，在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上的比例。我们可以用最大的个奇异值的矩阵和相乘来近似描述矩阵，从而实现了降维、减少数据存储、提升计算性能等效果。
-
u, s, v = numpy.linalg.svd(a, full_matrices=True, compute_uv=True, hermitian=False) 奇异
值分解。-
a 是一个形如(M,N)矩阵
-
full_matrices 的取值是为False或者True，默认值为True，这时u 的大小为(M,M)， v 的大小为(N,N)。否则u 的大小为(M,K)， v 的大小为(K,N) ，K=min(M,N)。
-
compute_uv 的取值是为False或者True，默认值为True，表示计算u,s,v 。为False的时候只计算s 。
-
总共有三个返回值u,s,v ， u 大小为(M,M)， s 大小为(M,N)， v 大小为(N,N)， a =usv 。
-
其中s 是对矩阵a 的奇异值分解。s 除了对角元素不为0 ，其他元素都为0 ，并且对角元素从大到小排列。s 中有n 个奇异值，一般排在后面的比较接近0，所以仅保留比较大的r个奇异值。
-
注：Numpy中返回的v 是通常所谓奇异值分解a=usv’ 中v 的转置。
import numpy as np
A = np.array([[4, 11, 14], [8, 7, ‐2]])
print(A)
# [[ 4 11 14]
# [ 8 7 ‐2]]
u, s, vh = np.linalg.svd(A, full_matrices=False)
print(u.shape) # (2, 2)
print(u)
# [[‐0.9486833 ‐0.31622777]
# [‐0.31622777 0.9486833 ]]
print(s.shape) # (2,)
print(np.diag(s))
# [[18.97366596 0. ]
# [ 0. 9.48683298]]
print(vh.shape) # (2, 3)
print(vh)
# [[‐0.33333333 ‐0.66666667 ‐0.66666667]
# [ 0.66666667 0.33333333 ‐0.66666667]]
a = np.dot(u, np.diag(s))
a = np.dot(a, vh)
print(a)
2. QR分解
QR（正交三角）分解法是求一般矩阵全部特征值的最有效并广泛应用的方法，一般矩阵先经过正交相似变化成为Hessenberg矩阵，然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R，所以称为QR分解法，与此正规正交矩阵的通用符号Q有关。
如果实（复）非奇异矩阵A能够化成正交（酉）矩阵Q与实（复）非奇异上三角矩阵R的乘积，即A=QR，则称其为A的QR分解。
使用qr分解有助于加快解方程或求解速度即收敛速度。
-
q,r = numpy.linalg.qr(a, mode=‘reduced’) 计算矩阵a 的QR分解。
-
a 是一个(M, N)的待分解矩阵。
-
mode = reduced ：返回(M, N)的列向量两两正交的矩阵q ，和(N, N)的三角阵r （Reduced QR分解）。
-
mode = complete ：返回(M, M)的正交矩阵q ，和(M, N)的三角阵r （Full QR分解）。
-
import numpy as np
A = np.array([[1, 1], [1, ‐2], [2, 1]])
print(A)
# [[ 1 1]
# [ 1 ‐2]
# [ 2 1]]
q, r = np.linalg.qr(A, mode='complete')
print(q.shape) # (3, 3)
print(q)
# [[‐0.40824829 0.34503278 ‐0.84515425]
# [‐0.40824829 ‐0.89708523 ‐0.16903085]
# [‐0.81649658 0.27602622 0.50709255]]
print(r.shape) # (3, 2)
print(r)
# [[‐2.44948974 ‐0.40824829]
# [ 0. 2.41522946]
# [ 0. 0. ]]
print(np.dot(q, r))
# [[ 1. 1.]
# [ 1. ‐2.]
# [ 2. 1.]]
a = np.allclose(np.dot(q, q.T), np.eye(3))
print(a) # True
3. Cholesky分解
Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的。Cholesky分解法又称平方根法，是当A为实对称正定矩阵时，LU三角分解法的变形。
重要性质：
-
若A对称正定，则 A − 1 A^{-1} A−1亦对称正定，且 a i i > 0 a_{ii} \gt 0 aii>0；
-
A A A的顺序主子阵 A k A_k Ak亦对称正定；
-
A的特征值 λ i > 0 \lambda_i>0 λi>0；
-
A的全部顺序主子式 d e t ( A k ) > 0 det(A_k) \gt 0 det(Ak)>0。（A能够作Cholesky分解的充要条件）
import numpy as np
A = np.array([[1, 1, 1, 1], [1, 3, 3, 3],
[1, 3, 5, 5], [1, 3, 5, 7]])
print(A)
# [[1 1 1 1]
# [1 3 3 3]
# [1 3 5 5]
# [1 3 5 7]]
print(np.linalg.eigvals(A))
# [13.13707118 1.6199144 0.51978306 0.72323135]
L = np.linalg.cholesky(A)
print(L)
# [[1. 0. 0. 0. ]
# [1. 1.41421356 0. 0. ]
# [1. 1.41421356 1.41421356 0. ]
# [1. 1.41421356 1.41421356 1.41421356]]
print(np.dot(L, L.T))
# [[1. 1. 1. 1.]
# [1. 3. 3. 3.]
# [1. 3. 5. 5.]
# [1. 3. 5. 7.]]
四、范数和其它数字
1. 矩阵的范数
矩阵的F范数：设矩阵
A
=
(
a
i
,
j
)
m
×
n
\mathbf A=(a_{i,j})_{m\times n}
A=(ai,j)m×n ，则其F 范数为：
∣
∣
A
∣
∣
F
=
∑
i
,
j
a
i
,
j
2
||\mathbf A||_F=\sqrt{\sum_{i,j}a_{i,j}^{2}}
∣∣A∣∣F=∑i,jai,j2 。它是向量的
L
2
L_2
L2范数的推广。
-
numpy.linalg.norm(x, ord=None, axis=None, keepdims=False) 计算向量或者矩阵的范数。
根据ord 参数的不同，计算不同的范数：
| ord | norm for matrices | norm for vectors |
|---|---|---|
| None | Frobenius norm | 2-norm |
| ‘fro’ | Frobenius norm | - |
| ‘nuc’ | nuclear norm | - |
求向量的范数：
import numpy as np
A = np.array([[1, 2, 3, 4], [2, 3, 5, 8],
[1, 3, 5, 7], [3, 4, 7, 11]])
print(A)
# [[ 1 2 3 4]
# [ 2 3 5 8]
# [ 1 3 5 7]
# [ 3 4 7 11]]
# 求向量的范数
print(np.linalg.norm(A, ord=1)) # 30.0
print(np.max(np.sum(A, axis=0))) # 30
print(np.linalg.norm(A, ord=2))
# 20.24345358700576
print(np.max(np.linalg.svd(A, compute_uv=False)))
# 20.24345358700576
print(np.linalg.norm(A, ord=np.inf)) # 25.0
print(np.max(np.sum(A, axis=1))) # 25
print(np.linalg.norm(A, ord='fro'))
# 20.273134932713294
print(np.sqrt(np.trace(np.dot(A.T, A))))
# 20.273134932713294
2. 方阵的行列式
- numpy.linalg.det(a) 计算行列式。
import numpy as np
x = np.array([[1, 2], [3, 4]])
print(x)
# [[1 2]
# [3 4]]
print(np.linalg.det(x))
# ‐2.0000000000000004
3. 矩阵的秩
- numpy.linalg.matrix_rank(M, tol=None, hermitian=False) 返回矩阵的秩。
import numpy as np
I = np.eye(3) # 先创建一个单位阵
print(I)
# [[1. 0. 0.]
# [0. 1. 0.]
# [0. 0. 1.]]
r = np.linalg.matrix_rank(I)
print(r) # 3
I[1, 1] = 0 # 将该元素置为0
print(I)
# [[1. 0. 0.]
# [0. 0. 0.]
# [0. 0. 1.]]
r = np.linalg.matrix_rank(I) # 此时秩变成2
print(r) # 2
4. 矩阵的迹
- numpy.trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None) 方阵的迹就是主对角元素之和。
import numpy as np
x = np.array([[1, 2, 3], [3, 4, 5], [6, 7, 8]])
print(x)
# [[1 2 3]
# [3 4 5]
# [6 7 8]]
y = np.array([[5, 4, 2], [1, 7, 9], [0, 4, 5]])
print(y)
# [[5 4 2]
# [1 7 9]
# [0 4 5]]
print(np.trace(x)) # A的迹等于A.T的迹
# 13
print(np.trace(np.transpose(x)))
# 13
print(np.trace(x + y)) # 和的迹 等于 迹的和
# 30
print(np.trace(x) + np.trace(y))
# 30
五、解方程和逆矩阵
1. 逆矩阵（inverse matrix）
设
A
A
A 是数域上的一个
n
n
n 阶矩阵，若在相同数域上存在另一个
n
n
n 阶矩阵
B
B
B，使得：
A
B
=
B
A
=
E
AB=BA=E
AB=BA=E （
E
E
E 为单
位矩阵），则我们称
B
B
B 是
A
A
A 的逆矩阵，而
A
A
A 则被称为可逆矩阵。
numpy.linalg.inv(a)计算矩阵a 的逆矩阵（矩阵可逆的充要条件： det(a) != 0 ，或者a 满
秩）。
计算矩阵的逆矩阵。
import numpy as np
A = np.array([[1, ‐2, 1], [0, 2, ‐1], [1, 1, ‐2]])
print(A)
# [[ 1 ‐2 1]
# [ 0 2 ‐1]
# [ 1 1 ‐2]]
# 求A的行列式，不为零则存在逆矩阵
A_det = np.linalg.det(A)
print(A_det)
# ‐2.9999999999999996
A_inverse = np.linalg.inv(A) # 求A的逆矩阵
print(A_inverse)
# [[ 1.00000000e+00 1.00000000e+00 ‐1.11022302e‐16]
# [ 3.33333333e‐01 1.00000000e+00 ‐3.33333333e‐01]
# [ 6.66666667e‐01 1.00000000e+00 ‐6.66666667e‐01]]
x = np.allclose(np.dot(A, A_inverse), np.eye(3))
print(x) # True
x = np.allclose(np.dot(A_inverse, A), np.eye(3))
print(x) # True
A_companion = A_inverse * A_det # 求A的伴随矩阵
print(A_companion)
# [[‐3.00000000e+00 ‐3.00000000e+00 3.33066907e‐16]
# [‐1.00000000e+00 ‐3.00000000e+00 1.00000000e+00]
# [‐2.00000000e+00 ‐3.00000000e+00 2.00000000e+00]]
2. 求解线性方程组
numpy.linalg.solve(a, b)求解线性方程组或矩阵方程。
# x + 2y + z = 7
# 2x ‐ y + 3z = 7
# 3x + y + 2z =18
import numpy as np
A = np.array([[1, 2, 1], [2, ‐1, 3], [3, 1, 2]])
b = np.array([7, 7, 18])
x = np.linalg.solve(A, b)
print(x) # [ 7. 1. ‐2.]
x = np.linalg.inv(A).dot(b)
print(x) # [ 7. 1. ‐2.]
y = np.allclose(np.dot(A, x), b)
print(y) # True
扫一扫
773
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
