正文：2 核心篇

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

5.3.3 矩阵的创建

直接生成

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)

间接创建

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]]

改变矩阵的大小

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]]

矩阵元素的存取

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

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

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

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

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

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

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

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

5.3.4 向量的创建

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特殊的矩阵

零矩阵

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

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

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

单位矩阵

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

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

对角矩阵

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

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

上三角矩阵

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]])

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

判断矩阵是否相等

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 矩阵乘法

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

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

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

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次方

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]])

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

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

5.6.1 转置矩阵

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

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 逆矩阵

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.]]

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

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

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

5.7.4 行列式的计算

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

5.8.4 矩阵的秩

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

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

5.9.1 向量的内积

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 向量的长度

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

5.9.4 标准正交基

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

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

解方程组：

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=B$, $X=A^{-1}B$

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

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

2 利用np.linalg.solve()

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

3 利用Sympy库的solve函数

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}

解具有无穷解的方程组

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

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}

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文件

import pandas as pd
import numpy as np

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

5.12 习题

(1)

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

1
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]])

A.I * B
matrix([[ 2. ],
        [-0.5],
        [ 0.5]])
2
np.linalg.solve(A, B)
matrix([[ 2. ],
        [-0.5],
        [ 0.5]])
3
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的随机整数， 计算着两个矩阵的乘积， 求这两个矩阵的秩

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]]))

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

np.linalg.matrix_rank(A)
3

np.linalg.matrix_rank(B)
4

(3)

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

import numpy as np

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

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

(4)

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

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

np.identity(4)
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 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的随机浮点数， 根据其上三角和下三角矩阵创建其对应的对称矩阵

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]])

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]])

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]])

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的特征值和特征向量

$|A-\lambda E|=
\begin{vmatrix}
4- \lambda & 2\\
1 & 5 - \lambda
\end{vmatrix}=0$

$(4 - \lambda)(5 - \lambda) = 0$

$\lambda ^ 2 - 9\lambda + 18 = 0$

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

$A - 3E =
\begin{bmatrix}
1 & 2\\
1 & 2
\end{bmatrix} \to$ 特征向量 $(-2, 1)^T$, 单位化得$(-\frac{2}{5}\sqrt 5,\frac {1}{5}\sqrt 5)^T$

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

import numpy as np

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

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

得特征向量$(-\frac{2}{5}\sqrt 5,\frac {1}{5}\sqrt 5)^T$, $(-\frac{1}{2}\sqrt 2, -\frac{1}{2}\sqrt 2)^T$

6.3 特征值分解

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

6.4 奇异值分解(SVD)

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

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

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

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

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

  • 奇异值是非负实数

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

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

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.      ]]

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

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

使用numpy.linalg.svd进行SVD

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]]

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值, 分别对应取$U$列, $\Sigma$变成k阶方阵,$V^T$取前k行, 对比利用矩阵乘积得到的新矩阵与原始数据的情况

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]])

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对图像进行压缩

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推荐商品

准备工作，调入相关的库

import numpy as np
from numpy import linalg as la
import pandas as pd
import time

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

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\theta = \frac{\mathbf{X} \cdot \mathbf{Y}}{\left | \mathbf{X} \right | * \left | \mathbf{Y} \right |}$

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

"""
以下三种计算方式的参数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()函数
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()函数
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()函数

  • 其中, $P _ {ui}$表示用户$u$对物品$j$的预测值, $item1$代表用户已打分的物品集合, $W _ {ji}$表示物品$j$和物品$i$的相似度, $r _ {ui}$表示用户$u$对物品$i$的打分
  • 该公式的含义: 与用户历史上感兴趣(评分高)的物品(向量)越相似(根据相似度判定算法得出)的物品(向量), 越有可能在用户的推荐列表中获得较高排名.
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()函数

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函数。

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()函数，获得结果。
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 用户——电影评分矩阵的获取

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 
# 输出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

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