正文：2 核心篇

5 将研究对象形式化——线性代数基础

5.3.3 矩阵的创建

直接生成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np

A = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
arr1 = np.array(A)
print("A=", A)
print("arr1=\n", arr1)
B = ((1, 2, 3, 4), (5, 6, 7, 8), (9, 10, 11, 12))
arr2 = np.array(B)
print("B=", B)
print("arr2=\n", arr2)

print("type(A)=", type(A))
print("type(B)=", type(B))
print("type(arr1)=", type(arr1))
print("arr1.shape=", arr1.shape)
A= [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
arr1=
 [[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]
B= ((1, 2, 3, 4), (5, 6, 7, 8), (9, 10, 11, 12))
arr2=
 [[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]
type(A)= <class 'list'>
type(B)= <class 'tuple'>
type(arr1)= <class 'numpy.ndarray'>
arr1.shape= (3, 4)

间接创建

1
2
3
4
5
6
import numpy as np

arr1 = np.random.random((2, 3))
print("arr1=\n", arr1)
arr2 = np.random.randint(3, 30, size=[2, 3])  # 3-30(不包括 30)
print("arr2=\n", arr2)
arr1=
 [[0.7888749  0.29995777 0.79464025]
 [0.48565204 0.77377983 0.84873221]]
arr2=
 [[19 17  7]
 [ 4 25 14]]

改变矩阵的大小

1
2
3
4
5
6
7
8
9
import numpy as np

A = [1, 2, 3, 4, 5, 6]
B = np.array(A)
C1 = B.reshape(2, 3)
C2 = B.reshape(3, 2)
print("B=", B)
print("C1=\n", C1)
print("C2=\n", C2)
B= [1 2 3 4 5 6]
C1=
 [[1 2 3]
 [4 5 6]]
C2=
 [[1 2]
 [3 4]
 [5 6]]

矩阵元素的存取

1
C1[0]  # 第 0 行
array([1, 2, 3])

1
C1[0:2]  #  0-2(不含 2)行
array([[1, 2, 3],
       [4, 5, 6]])

1
C2[[0, 2]]  # 第 0 行和第 2 行
array([[1, 2],
       [5, 6]])

1
C2[:, 1]  # 第 1 列
array([2, 4, 6])

1
C2[:, 0:2]  # 0-2(不含 2)列
array([[1, 2],
       [3, 4],
       [5, 6]])

1
C1[:, [0, 2]]  # 第 0 列和第 2 列
array([[1, 3],
       [4, 6]])

1
C2[2, 1]  # 第 2 行第 1 列
6

1
C2[2][1]  # 第 2 行第 1 列
6

5.3.4 向量的创建

1
2
3
4
5
6
7
8
9
10
import numpy as np

A = [[1, 2, 3, 4, 5]]
B = [[1], [2], [3], [4], [5]]
C = np.array(A)
D = np.array(B)
print("C=\n", C)
print("D=\n", D)
print("C.shape", np.shape(C))
print("D.shape", np.shape(D))
C=
 [[1 2 3 4 5]]
D=
 [[1]
 [2]
 [3]
 [4]
 [5]]
C.shape (1, 5)
D.shape (5, 1)

5.4 特殊的矩阵

零矩阵

1
np.zeros(10)
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

1
np.zeros((2, 4))
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.]])

1
np.array([np.zeros(10)])
array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

单位矩阵

1
np.eye(3)
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

1
np.identity(3)
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

对角矩阵

1
2
arr1 = np.diag([1, 2, 3])
arr1
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])

1
np.diag(arr1)
array([1, 2, 3])

上三角矩阵

1
2
A = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
np.triu(A, 0)
array([[ 1,  2,  3,  4],
       [ 0,  6,  7,  8],
       [ 0,  0, 11, 12],
       [ 0,  0,  0, 16]])

1
np.tril(A, 0)
array([[ 1,  0,  0,  0],
       [ 5,  6,  0,  0],
       [ 9, 10, 11,  0],
       [13, 14, 15, 16]])

判断矩阵是否相等

1
2
3
A = np.array([[1, 2, 3], [4, 5, 6]])
B = np.array([1, 2, 3, 4, 5, 6]).reshape(2, 3)
np.allclose(A, B)
True

5.5 矩阵基本操作

5.5.3 矩阵乘法

1
2
3
A = np.array([[1, 2], [1, 0]])
B = np.diag([2, 2])
A.dot(B)  # 矩阵乘法 AB
array([[2, 4],
       [2, 0]])

1
A * B  # np.array 类型, * 表示对应元素相乘
array([[2, 0],
       [0, 0]])

1
np.multiply(A, B)
array([[2, 0],
       [0, 0]])

1
2
3
A = np.mat([[1, 2], [1, 0]])
B = np.mat(np.diag([2, 2]))
A * B  # np.mat 类型, * 表示矩阵乘法
matrix([[2, 4],
        [2, 0]])

5.5.5 矩阵的乘方

只有方阵才可进行乘方运算.对 array 类型, 矩阵的乘方要经多次 dot 运算得到, 对 matrix 类型, 可以通过 ** 得到, array 类型的 ** 是每个元素的 n 次方

1
2
3
4
A = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
A_array = np.array(A)
A_matrix = np.mat(A)
A_array.dot(A_array).dot(A_array)
array([[ 468,  576,  684],
       [1062, 1305, 1548],
       [1656, 2034, 2412]])

1
A_array ** 3
array([[  1,   8,  27],
       [ 64, 125, 216],
       [343, 512, 729]], dtype=int32)

1
A_matrix ** 3
matrix([[ 468,  576,  684],
        [1062, 1305, 1548],
        [1656, 2034, 2412]])

5.6.1 转置矩阵

利用转置矩阵创建对称矩阵

1
2
3
4
5
6
7
8
9
import numpy as np

arr1 = np.random.randint(1, 16, size=[3, 3])
arr2 = np.triu(arr1)
arr2 += arr2.T - np.diag(np.diag(arr2))

print("arr1=\n", arr1)
print("arr2=\n", arr2)
print(np.allclose(arr2, arr2.T))
arr1=
 [[ 8 14 15]
 [13 13 12]
 [ 3  8 12]]
arr2=
 [[ 8 14 15]
 [14 13 12]
 [15 12 12]]
True

5.6.2 逆矩阵

1
2
3
4
5
6
7
8
9
10
import numpy as np

A = [[1, 2], [2, 5]]
C1 = np.array(A)
C2 = np.mat(A)
C1_inverse = np.linalg.inv(C1)  # array 类型不能使用 C1.I
C2_inverse = C2.I
print("C1_inverse=\n", C1_inverse)
print("C2_inverse=\n", C2_inverse)
print("C1.dot(C1_inverse)=\n", C1.dot(C1_inverse))
C1_inverse=
 [[ 5. -2.]
 [-2.  1.]]
C2_inverse=
 [[ 5. -2.]
 [-2.  1.]]
C1.dot(C1_inverse)=
 [[1. 0.]
 [0. 1.]]

1
C2 ** (-1)
matrix([[ 5., -2.],
        [-2.,  1.]])

“导致降维”, 行列式为 0 的矩阵不可逆

1
2
3
4
try:
    np.mat([[0, 0], [0, 1]]) ** (-1)
excerpt excerption as e:
    print(e)
Singular matrix

5.7.4 行列式的计算

1
np.linalg.det(np.diag([1, 2, 3, 4]))  # np.linalg 线性代数
23.999999999999993

5.8.4 矩阵的秩

1
np.linalg.matrix_rank(np.diag([1, 2, 3]))
3

1
np.linalg.matrix_rank(np.diag([1, 2, 0]))
2

5.9.1 向量的内积

1
2
3
4
5
A = [[1, 2, 3]]
B = [[4, 5, 6]]
C1 = np.array(A).reshape(3, 1)  # 转换成列向量
C2 = np.array(B).reshape(3, 1)
np.dot(C1.T, C2)
array([[32]])

5.9.2 向量的长度

1
2
A = np.array([[0, 3, 4]])
np.linalg.norm(A)
5.0

5.9.4 标准正交基

定义 5.17 在复数范围内满足ATA=AAT=EA^TA=AA^T=E, 则称AA为酉矩阵
酉矩阵常用于奇异值分解(SVD)中, 正交矩阵是实数特殊化的酉矩阵

综合实例——线性代数在实际问题中的应用

解方程组：

{x+2y+z=72x−y+3z=73x+y+2z=18\left\{\begin{matrix} x + 2y + z = 7 \\ 2x - y + 3z = 7 \\ 3x + y + 2z =18 \end{matrix}\right. 

1
2
3
4
5
6
import numpy as np

A = np.mat([[1, 2, 1], [2, -1, 3], [3, 1, 2]])
B = np.mat([7, 7, 18]).reshape(3, 1)
AB = np.hstack((A, B))
AB
matrix([[ 1,  2,  1,  7],
        [ 2, -1,  3,  7],
        [ 3,  1,  2, 18]])

1 利用逆矩阵求AX=BAX=B, X=A−1BX=A^{-1}B

1
2
X = A.I * B
X
matrix([[ 7.],
        [ 1.],
        [-2.]])

1
A * X  # 验证
matrix([[ 7.],
        [ 7.],
        [18.]])

2 利用 np.linalg.solve()

1
np.linalg.solve(A, B)
matrix([[ 7.],
        [ 1.],
        [-2.]])

3 利用 Sympy 库的 solve 函数

1
2
3
4
5
6
7
import sympy
from sympy.abc import x, y, z

eq = [x + 2 * y + z - 7, 
      2 * x - y + 3 * z - 7, 
      3 * x + y + 2 * z - 18]
sympy.solve(eq, [x, y, z])
{x: 7, y: 1, z: -2}

解具有无穷解的方程组

{x+2y+z−2w=02x+3y−w=0x−y−5z+7w=0\left\{\begin{matrix} x + 2y + z - 2w = 0 \\ 2x + 3y - w = 0 \\ x - y - 5z + 7w = 0 \end{matrix}\right. 

1
2
3
4
5
6
import numpy as np

A = np.array([[1, 2, 1, -2],
             [2, 3, 0, -1],
             [1, -1, -5, 7]])
np.linalg.matrix_rank(A)
2

1
2
3
4
5
6
import sympy
from sympy.abc import x, y, z, w

eq = [x + 2 * y + z - 2 * w, 2 * x + 3 * y - w, x - y - 5 * z + 7 * w]
result = sympy.solve(eq, [x, y, z, w])
result
{x: -4*w + 3*z, y: 3*w - 2*z}

1
2
3
4
A = {z:1, w:2}
x = float(result[x].evalf(subs=A))
y = float(result[y].evalf(subs=A))
x, y
(-5.0, 4.0)

例 5.18 获取数据集的数据 读取 csv 文件

1
2
3
4
5
import pandas as pd
import numpy as np

dataset = pd.read_csv("iris.csv")
data = np.array(dataset)

5.12 习题

(1)

分别使用三种方法求线性方程组的解

{x+y+z=2x+2y+4z=3x+3y+9z=5\left \{ \begin{matrix} x+y+z=2 \\x+2y+4z=3 \\x+3y+9z=5 \end{matrix}\right .

1
1
2
3
4
5
6
7
8
import numpy as np

A = np.mat([[1, 1, 1],
            [1, 2, 4],
            [1, 3, 9]])
B = np.mat([2, 3, 5]).reshape(3, 1)
AB = np.hstack((A, B))
AB
matrix([[1, 1, 1, 2],
        [1, 2, 4, 3],
        [1, 3, 9, 5]])

1
A.I * B
matrix([[ 2. ],
        [-0.5],
        [ 0.5]])
2
1
np.linalg.solve(A, B)
matrix([[ 2. ],
        [-0.5],
        [ 0.5]])
3
1
2
3
4
5
import sympy
from sympy.abc import x, y, z

eq = [x + y + z - 2, x + 2 * y + 4 * z - 3, x + 3 * y + 9 * z - 5]
sympy.solve(eq, [x, y, z])
{x: 2, y: -1/2, z: 1/2}

(2)

分别创建 3 * 4 阶和 4 * 5 阶的矩阵，元素值为 1-20 的随机整数，计算着两个矩阵的乘积，求这两个矩阵的秩

1
2
3
4
5
import numpy as np

A = np.random.randint(1, 20, size=[3, 4])
B = np.random.randint(1, 20, size=[4, 5])
A, B
(array([[ 2, 18, 16,  6],
        [19,  1, 12, 19],
        [11,  9, 18, 17]]),
 array([[ 7, 14, 14, 14,  6],
        [11, 10, 11,  8, 10],
        [ 4,  3,  8,  4, 11],
        [11, 10, 14,  1,  3]]))

1
A.dot(B)
array([[342, 316, 438, 242, 386],
       [401, 502, 639, 341, 313],
       [435, 468, 635, 315, 405]])

1
np.linalg.matrix_rank(A)
3

1
np.linalg.matrix_rank(B)
4

(3)

求下面矩阵的逆矩阵，求逆矩阵时要先求行列式，行列式不为 0 时逆矩阵存在，之后再进行求逆操作

A=[123221343]A=\begin{bmatrix} 1 & 2 & 3 \\ 2 & 2 & 1 \\ 3 & 4 & 3 \end{bmatrix} 

1
2
3
4
5
6
import numpy as np

A = np.mat([[1, 2, 3],
            [2, 2, 1],
            [3, 4, 3]])
np.linalg.det(A)
1.9999999999999993

1
A.I
matrix([[ 1. ,  3. , -2. ],
        [-1.5, -3. ,  2.5],
        [ 1. ,  1. , -1. ]])

(4)

分别创建四阶零矩阵和四阶单位矩阵，以及对角线元素分别为 1，2，3，4 的对角矩阵

1
np.zeros((4, 4))
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])

1
np.identity(4)
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

1
np.diag([1, 2, 3, 4])
array([[1, 0, 0, 0],
       [0, 2, 0, 0],
       [0, 0, 3, 0],
       [0, 0, 0, 4]])

(5)

创建一个四阶方阵，元素值为 1~20 的随机浮点数，根据其上三角和下三角矩阵创建其对应的对称矩阵

1
2
A = np.random.random(16).reshape(4, 4) * 20
A
array([[14.08629999,  3.99280371, 13.3828697 , 13.32337192],
       [11.80703601, 11.72735218,  7.43395284, 14.95671889],
       [18.40134451, 17.69715348, 14.43624563,  4.31954707],
       [11.07252469,  7.31171911, 16.52018505, 11.79353701]])

1
np.triu(A)
array([[14.08629999,  3.99280371, 13.3828697 , 13.32337192],
       [ 0.        , 11.72735218,  7.43395284, 14.95671889],
       [ 0.        ,  0.        , 14.43624563,  4.31954707],
       [ 0.        ,  0.        ,  0.        , 11.79353701]])

1
np.triu(A).T + np.triu(A) - np.diag(np.diag(A))
array([[14.08629999,  3.99280371, 13.3828697 , 13.32337192],
       [ 3.99280371, 11.72735218,  7.43395284, 14.95671889],
       [13.3828697 ,  7.43395284, 14.43624563,  4.31954707],
       [13.32337192, 14.95671889,  4.31954707, 11.79353701]])

1
np.tril(A).T + np.tril(A) - np.diag(np.diag(A))
array([[14.08629999, 11.80703601, 18.40134451, 11.07252469],
       [11.80703601, 11.72735218, 17.69715348,  7.31171911],
       [18.40134451, 17.69715348, 14.43624563, 16.52018505],
       [11.07252469,  7.31171911, 16.52018505, 11.79353701]])

6 从数据中提取重要信息——特征值与矩阵分解

6.1.4 特征值的实现

已知

A=[4215]A= \begin{bmatrix} 4 & 2 \\ 1 & 5 \end{bmatrix}

求 A 的特征值和特征向量

∣A−λE∣=∣4−λ215−λ∣=0|A-\lambda E|= \begin{vmatrix} 4- \lambda & 2\\ 1 & 5 - \lambda \end{vmatrix}=0

(4−λ)(5−λ)=0(4 - \lambda)(5 - \lambda) = 0

λ2−9λ+18=0\lambda ^ 2 - 9\lambda + 18 = 0

特征值: λ1=3\lambda _ {1} = 3, λ2=6\lambda _ {2} = 6

A−3E=[1212]→A - 3E = \begin{bmatrix} 1 & 2\\ 1 & 2 \end{bmatrix} \to 特征向量 (−2,1)T(-2, 1)^T, 单位化得(−255,155)T(-\frac{2}{5}\sqrt 5,\frac {1}{5}\sqrt 5)^T

A−6E=[−221−1]→A - 6E = \begin{bmatrix} -2 & 2\\ 1 & -1 \end{bmatrix} \to 特征向量 (1,1)T(1, 1)^T, 单位正交化得(−122,−122)T(-\frac{1}{2}\sqrt 2, -\frac{1}{2}\sqrt 2)^T

1
2
3
4
import numpy as np

A = np.array([[4, 2], [1, 5]])
eig_val, eig_vex = np.linalg.eig(A)  # eig()函数分别求解特征值和特征向量矩阵
1
eig_val  # 特征值
array([3., 6.])

1
eig_vex  # 特征向量矩阵(已标准化)
array([[-0.89442719, -0.70710678],
       [ 0.4472136 , -0.70710678]])

得特征向量(−255,155)T(-\frac{2}{5}\sqrt 5,\frac {1}{5}\sqrt 5)^T, (−122,−122)T(-\frac{1}{2}\sqrt 2, -\frac{1}{2}\sqrt 2)^T

6.3 特征值分解

特征值分解是将矩阵AA分解成A=QΣQ−1A=Q\Sigma Q^{-1}的形式, 要求矩阵必须是 n 维方阵, 将特征值分解算法推广到所有矩阵之上, 就是更加通用的奇异值分解(SVD)

6.4 奇异值分解(SVD)

  • 这个知识点本科和考研居然没有学过……大概是因为计算量太大了吧

  • 参考教程【学长小课堂】什么是奇异值分解 SVD–SVD 如何分解时空矩阵

  • SVD 可以将任意的矩阵分解为三个矩阵相乘, 简记为A=UΣVTA=U\Sigma V^T

  • 其中UU和VV是正交矩阵, 分别对应AATAA^T和ATAA^TA的特征向量,Σ\Sigma是一个对角矩阵, 对角元素为AA的特征值(按降序排列), 奇异值2^2=特征值

  • 可以把AA看作是一个线性变换,将 A 分解为: UU旋转, Σ\Sigma拉伸, VTV^T旋转

  • 奇异值是非负实数

  • 定理：设 A 和 B 是两个矩阵（注意不限定大小，但是必须满足做矩阵积 AB 和 BA 都有意义），则 AB 与 BA 拥有相同的非零特征值，并且保持相同重数，并且属于非零特征值的特征向量间存在制约关系。为什么 A^TA 与 AA^T 有相同特征值？

Am∗n=Um∗mΣm∗nVm∗nTA _ {m * n} = U _ {m * m}\Sigma _ {m * n}V _ {m * n}^T

[011110]=[1/61/21/32/60−1/31/6−1/21/3][300100][1/21/2−1/21/2]\begin{bmatrix} 0 & 1\\ 1 & 1\\ 1 & 0 \end{bmatrix}= \begin{bmatrix} 1/\sqrt6 & 1/\sqrt2 & 1/\sqrt3\\ 2/\sqrt6 & 0 & -1/\sqrt3\\ 1/\sqrt6 & -1/\sqrt2 & 1/\sqrt3 \end{bmatrix} \begin{bmatrix} \sqrt3 & 0\\ 0 & 1\\ 0 & 0 \end{bmatrix} \begin{bmatrix} 1/\sqrt2 & 1/\sqrt2\\ -1/\sqrt2 & 1/\sqrt2 \end{bmatrix} 

1
2
3
4
5
import numpy as np

# A = np.array([[1, 5, 7, 6, 1], [2, 1, 10, 4, 4], [3, 6, 7, 5, 2]])
A = np.array([[0, 1], [1, 1], [1, 0]])  # https://zhuanlan.zhihu.com/p/29846048
A
array([[0, 1],
       [1, 1],
       [1, 0]])

使用定义进行 SVD

1
2
3
4
5
6
7
8
9
10
11
12
sigmal_val, U = np.linalg.eigh(A.dot(A.T))  # 计算特征值, 左奇异矩阵
sigmal_sort_id = np.argsort(sigmal_val)[::-1]  # 让特征值降序排列
sigmal_val = np.sort(sigmal_val)[::-1]  # 特征值对应的特征向量也对应排序
U = U[:, sigmal_sort_id]  # 左奇异矩阵按这个 id 排序
sigmal = np.diag(np.sqrt(sigmal_val))  # 计算奇异值矩阵(方阵)
sigmal_inv = np.linalg.inv(sigmal)  # 奇异值矩阵的逆矩阵
V_part_T = sigmal_inv.dot(U.T).dot(A)  # 右奇异矩阵

print("AAT=:\n", AAT)
print("左奇异矩阵: \n", U)
print("奇异值矩阵: \n", np.round(sigmal, 2))
print("右奇异矩阵的转置(部分):\n", np.round(V_part_T, 6))
AAT=:
 [[1 1 0]
 [1 2 1]
 [0 1 1]]
左奇异矩阵: 
 [[ 4.08248290e-01 -7.07106781e-01  5.77350269e-01]
 [ 8.16496581e-01 -7.58447699e-17 -5.77350269e-01]
 [ 4.08248290e-01  7.07106781e-01  5.77350269e-01]]
奇异值矩阵: 
 [[1.73 0.   0.  ]
 [0.   1.   0.  ]
 [0.   0.   0.  ]]
右奇异矩阵的转置(部分):
 [[ 0.707107  0.707107]
 [ 0.707107 -0.707107]
 [ 0.        0.      ]]

1
np.round(U.dot(sigmal).dot(V_part_T), 2)
array([[0., 1.],
       [1., 1.],
       [1., 0.]])

不研究了, 费劲, 怪不得考研不考, 电脑都鼓捣老半天, 别说手算了

使用 numpy.linalg.svd 进行 SVD

1
2
3
4
5
6
7
U, s, VT = np.linalg.svd(A)
Sigma = np.zeros(np.shape(A))  # 先创建一个零矩阵
Sigma[:len(s), :len(s)] = np.diag(s)  # 将该矩阵的对角元素替换为奇异值
print("左奇异矩阵: \n", U)
print("奇异值:\n", s)
print("奇异值矩阵: \n", Sigma)
print("右奇异矩阵的转置:\n", VT)
左奇异矩阵: 
 [[-4.08248290e-01  7.07106781e-01  5.77350269e-01]
 [-8.16496581e-01  5.55111512e-17 -5.77350269e-01]
 [-4.08248290e-01 -7.07106781e-01  5.77350269e-01]]
奇异值:
 [1.73205081 1.        ]
奇异值矩阵: 
 [[1.73205081 0.        ]
 [0.         1.        ]
 [0.         0.        ]]
右奇异矩阵的转置:
 [[-0.70710678 -0.70710678]
 [-0.70710678  0.70710678]]

1
np.round(U.dot(Sigma).dot(VT), 2)
array([[ 0.,  1.],
       [ 1.,  1.],
       [ 1., -0.]])

6.5.4 利用 SVD 进行矩阵近似

在下述 A 中, SVD 得到的特征值前两项较大(343, 25),最后一项较小(3.37). 通过取不同的 k 值, 分别对应取UU列, Σ\Sigma变成 k 阶方阵,VTV^T取前 k 行, 对比利用矩阵乘积得到的新矩阵与原始数据的情况

1
2
3
4
import numpy as np

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

1
2
3
4
5
6
7
U, s, VT = np.linalg.svd(A)
Sigma = np.zeros(np.shape(A))  # 先创建一个零矩阵
Sigma[:len(s), :len(s)] = np.diag(s)  # 将该矩阵的对角元素替换为奇异值

for k in range(3, 0, -1):
    D = U[:, :k].dot(Sigma[:k, :k].dot(VT[:k, :]))
    print("k=", k, "压缩后的矩阵: \n", np.round(D, 1))
k= 3 压缩后的矩阵: 
 [[ 1.  5.  7.  6.  1.]
 [ 2.  1. 10.  4.  4.]
 [ 3.  6.  7.  5.  2.]]
k= 2 压缩后的矩阵: 
 [[ 2.   5.4  6.8  5.3  1.5]
 [ 2.   1.  10.   4.   4. ]
 [ 2.1  5.7  7.2  5.6  1.5]]
k= 1 压缩后的矩阵: 
 [[1.9 3.8 7.7 4.8 2.3]
 [2.1 4.1 8.2 5.1 2.4]
 [2.  4.  8.1 5.  2.4]]

可以选取合适的 k 值, 保留较大的奇异值和特征向量, 实现用较少的数据量达到较好的矩阵近似效果.

6.6 综合实例 1 ——利用 SVD 对图像进行压缩

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt


def get_approx_SVD1(data, percent):
    """
    SVD 并还原压缩后的数据, 取整体奇异值的百分比
    :param data: 原始矩阵
    :param percent: 奇异值总和百分比
    """
    U, s, VT = np.linalg.svd(data)
    Sigma = np.zeros(np.shape(data))
    Sigma[:len(s), :len(s)] = np.diag(s)
    count = (int)(sum(s)) * percent  # 取整体奇异值的百分比
    k = -1  # k 是奇异值总和的百分比的个数
    curSum = 0  # 初始为第一个奇异值
    while curSum <= count:  # 当取的奇异值之和 <= 应取得奇异值之和 时:
        k += 1
        curSum += s[k]
    D = U[:, :k].dot(Sigma[:k, :k].dot(VT[:k, :]))  # 取 U 的前 k 列, Sigma 的前 k 阶方阵, VT 的前 k 行
    D[D < 0] = 0  # 钳制颜色范围
    D[D > 255] = 255
    # numpy.rint()将数组的元素四舍五入为最接近的整数
    return np.rint(D).astype("uint8")


def get_approx_SVD2(data, percent):
    """
    SVD 并还原压缩后的数据, 取奇异值之和的百分比
    :param data: 原始矩阵
    :param percent: 奇异值个数百分比
    """
    U, s, VT = np.linalg.svd(data)
    Sigma = np.zeros(np.shape(data))
    Sigma[:len(s), :len(s)] = np.diag(s)
    k = (int)(percent * len(s))  # 取最大的前 k 个奇异值
    D = U[:, :k].dot(Sigma[:k, :k].dot(VT[:k, :]))
    D[D < 0] = 0
    D[D > 255] = 255
    return np.rint(D).astype("uint8")


def rebuild_img(filename, p, get_approx_SVD, index):
    """
    导入图像, 进行 SVD 压缩, 并重构图像
    :param filename: 文件名
    :param p: 百分比
    :param get_approx_SVD: 调用的 SVD 筛选方法
    :param index: 显示图像的序号
    """
    img = Image.open(filename, 'r')
    a = np.array(img)
    R0 = a[:, :, 0]  # 红色
    G0 = a[:, :, 1]  # 绿色
    B0 = a[:, :, 2]  # 蓝色
    R = get_approx_SVD(R0, p)
    G = get_approx_SVD(G0, p)
    B = get_approx_SVD(B0, p)
    I = np.stack((R, G, B), 2)

    # 加载图像
    plt.subplot(2, 5, index)
    if get_approx_SVD == get_approx_SVD1:
        plt.title("SVD1 " + str((int)(p * 100)) + '%')
    else:
        plt.title("SVD2 " + str((int)(p * 100)) + '%')
    plt.imshow(I)
    plt.axis('off')


if __name__ == "__main__":
    filename = "lenna.bmp"
    plt.figure(figsize=(10, 5))
    index = 0
    for p in np.arange(0.2, 1.2, 0.2):
        index += 1
        rebuild_img(filename, p, get_approx_SVD1, index)
    for p in np.arange(0.2, 1.2, 0.2):
        index += 1
        rebuild_img(filename, p, get_approx_SVD2, index)
    plt.show()

png

6.7 综合实例 2——利用 SVD 推荐商品

准备工作，调入相关的库

1
2
3
4
import numpy as np
from numpy import linalg as la
import pandas as pd
import time

1 加载测试数据集，形成 loadExData()函数

1
2
3
4
5
6
7
8
9
10
11
12
13
def loadExData():
    """
    行: 代表用户
    列: 代表电影
    值: 代表用户对电影的评分, 0 表示未评分
    :return: 将 csv 转换成的 np.array 数组
    """
    """
    pd.read_csv():
    sep：读取 csv 文件时指定的分隔符，默认为逗号。
    header：设置导入 DataFrame 的列名称(第一行)，默认为 "infer"
    """
    return np.array(pd.read_csv("评分.csv", sep=",", header=None))

2 计算相似度, 形成 ecludSim()欧式距离, pearsSim()皮尔逊相关系数, cosSim() 余弦相似度 三个判断相似度的函数

  • 欧氏距离就是两点间距离公式, 值域为 R

  • 皮尔逊相关系数应该就是概率论与数理统计的那个 r?我记得初中的时候利用计算器的这个 bug 爆机, 值域为[-1, 1]

  • 余弦相似度就是判断两个向量的夹角余弦值, 为 0 就说明两个向量垂直, 不相关, 值域为[-1, 1]

  • 欧氏距离: $dist(\mathbf{X} \cdot \mathbf{Y}) = \left | \mathbf{X} - \mathbf{Y} \right | $

  • 皮尔逊相关系数: $r = \frac{( \mathbf{X-\bar{X}} ) \cdot ( \mathbf{Y-\bar{Y}} )}
    {\left | \mathbf{X-\bar{X}} \right | * \left | \mathbf{Y-\bar{Y}} \right |} $

  • 余弦相似度: cosθ=X⋅Y∥X∥∗∥Y∥cos\theta = \frac{\mathbf{X} \cdot \mathbf{Y}}{\left \| \mathbf{X} \right \| * \left \| \mathbf{Y} \right \|}

  • 最大最小归一化: x′=x−min(x)max(x)−min(x)x' = \frac{x - min(x)}{max(x) - min(x)}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
"""
以下三种计算方式的参数 X 和 Y 都采用一维数组(向量), 返回的值越接近 1, 相关度越强
"""


def ecludSim(X, Y):
    """
    利用欧氏距离计算相似度, 并归一化
    """
    return 1.0 / (1.0 + la.norm(X - Y))


def pearsSim(X, Y):
    """
    利用皮尔逊相关系数计算相似度, 并归一化
    """
    if len(X < 3):
        return 1.0
    else:
        return 0.5 + 0.5 * np.corrcoef(X, Y, rowvar=1)[0][1]
    
    
def cosSim(X, Y):
    """
    利用余弦相似度计算相似度, 并归一化
    """
    return 0.5 + 0.5 * (float(X.dot(Y)) / (la.norm(X) * la.norm(Y)))

3 对矩阵降维处理, 形成 svd_item()函数

1. 计算选取的奇异值数目 k 值, 形成 SigmaPct()函数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
def SigmaPct(sigma, percentage):
    """
    按照前 k 个奇异值的平方和占总奇异值的平方和的百分比 percentage 确定 k 的值
    后续计算 SVD 时需要将 item 原始矩阵降维
    """
    sum_sigma = sum(sigma ** 2)
    sum_sigma1 = sum_sigma * percentage  # 求所有奇异值 sigma 的平方和的百分比
    sum_sigma2 = 0
    k = 0
    for i in sigma:
        sum_sigma2 += i * i
        k += 1
        if sum_sigma2 >= sum_sigma1:
            return k
2. 降维处理, 形成 svd_item()函数
1
2
3
4
5
6
7
8
9
10
11
def svd_item(data, percentage):
    """
    降维处理
    :return: 降维的物品数据
    """
    n = np.shape(data)[1]  # 物品种类数据
    U, s, VT = la.svd(data)  # 数据集进行奇异值分解，返回的 s 为对角线上的值
    k = SigmaPct(s, percentage)  # 确定了 k 的值,前 k 个已经包含了 percentage 的能力
    Sigma = np.eye(k) * s[:k]  # 构建对角矩阵
    # 将数据转换到 k 维空间(低维)，构建转换后的物品
    return data.T.dot(U[:, :k].dot(la.inv(Sigma)))  # 返回降维的物品数据

4 在已经降维的数据中, 对用户未打分的一个物品进行评分预测, 形成 svd_predict()函数

Puj=Σi∈item1WjiruiΣi∈item1WjiP _ {uj} = \frac {\Sigma _ {i \in item1} W _ {ji} r _ {ui}}{\Sigma _ {i \in item1}W _ {ji}}

  • 其中, PuiP _ {ui}表示用户uu对物品jj的预测值, item1item1代表用户已打分的物品集合, WjiW _ {ji}表示物品jj和物品ii的相似度, ruir _ {ui}表示用户uu对物品ii的打分
  • 该公式的含义: 与用户历史上感兴趣(评分高)的物品(向量)越相似(根据相似度判定算法得出)的物品(向量), 越有可能在用户的推荐列表中获得较高排名.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def svd_predict(data, user, simMeas, FormedItems, item):
    """
    基于 item 的相似性对用户未评过分的物品进行预测评分
    :param data: 数组矩阵
    :param user: 用户编号
    :param simMeas: 相似度算法
    :param FormedItems: 数组矩阵的行对应用户
    :param item: 列对应物品
    :return: 对物品的预测评分，返回后用于分数的排序
    """
    n = np.shape(data)[1]  # 得到数据集中的物品种类数据(各个列)
    Totalsim = 0.0  # 初始化两个评分值
    TotalratSim = 0.0  # 相似性总和变量初始化
    # 遍历给定的用户行中的每个物品
    # 即（对用户评过分的物品进行遍历，并将它与其他物品进行比较），计算相似度
    for j in range(n):
        # 得到给定的用户 user 对商品的评分
        Rating_user = data[user, j]
        # 只对 评价过的商品 和 不是自己的商品 求相似度
        if Rating_user != 0 and j != item:
            # 计算 svd 转换过后矩阵的相似度,物品 item 与物品 j 之间的相似度
            # 相似度的计算方法也会作为一个参数传递给该函数
            Similarity = simMeas(FormedItems[item, :], FormedItems[j, :])
            Totalsim += Similarity  # 对相似度不断累加求和
            TotalratSim += Similarity * Rating_user  # 对相似度及对应评分值乘积求和
    if Totalsim == 0:
        return 0
    else:
        return TotalratSim / Totalsim  # 得到对物品的预测评分，返回后用于分数的排序

5 产生前 N 个评分值高的物品，返回物品编号以及预测评分值，形成 recommend()函数

1
2
3
4
5
6
7
8
9
10
11
12
13
def recommend(data, user, FormedItems, N, simMeas):
    # 为未评价的物品建立一个用户未评分 item 的列表
    unratedItems = np.array(np.nonzero(data[user, :] == 0))[0]
    if len(unratedItems) == 0:
        return "你已评价完所有物品"  # 若都已经评过分，则退出
    Scoresitem = []
    for item in unratedItems:  # 对未评分的物品 item，都计算其预测评分
        # 计算评价值
        estimatedScore = svd_predict(data, user, simMeas, FormedItems,
                                     item)
        Scoresitem.append((item, estimatedScore))  # 记录商品及评价值
    Scoresitem = sorted(Scoresitem, key=lambda x: x[1], reverse=True)  # 按照得分逆序排序
    return Scoresitem[:N]  # 返回前 N 个评分的物品名

6 对指定用户进行商品推荐，形成 recommend_predict 函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def recommend_predict():
    user_item = loadExData()  # 读取数据
    percentage = 0.9  # 奇异值平方和的百分比
    user = 1  # 预测的用户
    n = 4  # 推荐个数
    FormedItems = svd_item(user_item, percentage)  # 获得 SVD 降维后的物品
    
    s_t = time.time()
    simMeas = cosSim  # 相似度
    r1 = recommend(user_item, user, FormedItems, n, simMeas)
    print('利用余弦相似度计算距离，进行的奇异值分解推荐：')
    print("按相似度推荐的物品编号为：", r1, "\n 用时:", time.time() - s_t)
    
    s_t = time.time()
    simMeas = ecludSim
    r2 = recommend(user_item, user, FormedItems, n, simMeas)
    print('利用欧氏距离计算距离，进行的奇异值分解推荐：')
    print("按相似度推荐的物品编号为：", r2, "\n 用时:", time.time() - s_t)
    
    s_t = time.time()
    simMeas = pearsSim
    r3 = recommend(user_item, user, FormedItems, n, simMeas)
    print('利用皮尔逊相关系数计算距离，进行的奇异值分解推荐：')
    print("按相似度推荐的物品编号为：", r3, "\n 用时:", time.time() - s_t)
7 调用 recommend_predict()函数，获得结果。
1
2
if __name__ == "__main__":
    recommend_predict()
利用余弦相似度计算距离，进行的奇异值分解推荐：
按相似度推荐的物品编号为：[(1605, 3.7280993014041743), (935, 3.72778748418618), (366, 3.7274892395964807), (812, 3.7260645630753846)] 
用时: 1.6344106197357178
利用欧氏距离计算距离，进行的奇异值分解推荐：
按相似度推荐的物品编号为：[(8, 3.690142726938304), (514, 3.689943318086311), (14, 3.6894721117751406), (172, 3.689287056157151)] 
用时: 1.240149736404419
利用皮尔逊相关系数计算距离，进行的奇异值分解推荐：
按相似度推荐的物品编号为：[(1, 3.7096774193548385), (2, 3.7096774193548385), (3, 3.7096774193548385), (4, 3.7096774193548385)] 
用时: 0.8249673843383789

6.8.2 用户——电影评分矩阵的获取

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np
import pandas as pd

# 1、读入数据集 header = [用户 id, 电影 id, 评分(取整), 用户提交时间所用秒数]
header = ["user_id","item_id","rating","timestamp"]
data = pd.read_csv("u.data", sep="\t", names=header)  # 以 csv 文件格式读取 u.data, 分隔符是\t, 列名是 header
# 2、生成用户—物品评分矩阵
# 检查是否有重复的用户物品打分记录(数据清洗)
# data.duplicated()   # 返回布尔型数据，告诉重复值的位置
data.duplicated(subset = ["user_id","item_id"]).sum()  # 返回重复值的个数
item_id_user = data.groupby("item_id").count()["user_id"]  # 返回表格(物品 id, 评分用户次数)
print(item_id_user)
# 构建用户物品矩阵
users_num = data.user_id.max()  # 选择最大的 id, 可以剔除重复数据
items_num = data.item_id.max()
user_item_rating = np.zeros((users_num,items_num))
for line in data.itertuples():  # itertuples()：将 DataFrame 迭代成元组
    # 以元组的方式赋值
    user_item_rating[line[1] - 1, line[2] - 1] = line[3]
# 生成用户——电影评分矩阵存储成文本文件, 数据为 user_item_rating, 分隔符为，
np.savetxt("评分 2.csv", user_item_rating, delimiter = ",")
# 输出 u.data 中的前 5 行内容
data.head()
item_id
1       452
2       131
3        90
4       209
5        86
       ... 
1678      1
1679      1
1680      1
1681      1
1682      1
Name: user_id, Length: 1682, dtype: int64
user_id item_id rating timestamp
0 196 242 3 881250949
1 186 302 3 891717742
2 22 377 1 878887116
3 244 51 2 880606923
4 166 346 1 886397596 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 输出 u.data 中的大小
print("数据集的大小", data.shape)
# 输出客户数和电影数
print("客户数 =", users_num)
print("电影数 =", items_num)
print("客户数 * 电影数:", users_num * items_num)
# 查看生成的 user_item_rating 非零元素
print("user_item_rating 中的非零元素个数", len(user_item_rating.nonzero()[1]))
print("user_item_rating 中的非零元素", user_item_rating.nonzero()[1])
# 查看生成的 user_item_rating 矩阵的稀疏性
sparsity = round(len(user_item_rating.nonzero()[1]) / float(users_num * items_num), 3)
print("user_item_rating 矩阵的稀疏性：", sparsity)
print("user_item_rating 矩阵的大小", user_item_rating.shape)
# 以表格形式显示用户—电影评分表矩阵, 列名为电影, 行名为用户
pd.DataFrame(user_item_rating)
数据集的大小 (100000, 4)
客户数 = 943
电影数 = 1682
客户数 * 电影数: 1586126
user_item_rating 中的非零元素个数 100000
user_item_rating 中的非零元素 [   0    1    2 ... 1187 1227 1329]
user_item_rating 矩阵的稀疏性：0.063
user_item_rating 矩阵的大小 (943, 1682)
0 1 2 3 4 5 6 7 8 9 ... 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681
0 5.0 3.0 4.0 3.0 3.0 5.0 4.0 1.0 5.0 3.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1 4.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4 4.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
938 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 5.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
939 0.0 0.0 0.0 2.0 0.0 0.0 4.0 5.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
940 5.0 0.0 0.0 0.0 0.0 0.0 4.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
941 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
942 0.0 5.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

943 rows × 1682 columns

Icon
Icon
Icon
Icon
上一篇：Python-人工智能数学基础(7-9)
下一篇：Python-人工智能数学基础(2-4)