Python数学建模算法与应用（司守奎）--第3章随笔

原创已于 2023-10-02 16:44:07 修改 · 1.1k 阅读

10 ·

CC 4.0 BY-SA版权

文章标签：

#数学建模 #算法 #python

于 2023-09-27 21:33:09 首次发布

本文介绍了线性代数中的特征值与特征向量概念，通过Fibonacci数列和Leslie种群模型示例展示如何运用特征值方法求解，同时探讨了PageRank算法的基本原理及其在实际问题中的应用，包括矩阵奇异值分解在推荐系统和图像压缩中的作用。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

第 3 章线性代数模型

3.1 特征值于特征向量

3.1.1 差分方程

例3.1 Fibonacci 数列的通项

解法一运用特征值和特征向量求通式

# 程序文件 ex3_1_1.py
import sympy as sp
sp.var('k', positive=True, integer=True)
a = sp.Matrix([[0, 1], [1, 1]])
val = a.eigenvals()                 # 求特征值
vec = a.eigenvects()                # 求特征向量
P, D = a.diagonalize()              # 把 a 相似对角化
ak = P @ (D**k) @ (P.inv())
F = ak @ sp.Matrix([1, 1])
s = sp.simplify(F[0])
print(s)
sm = []
for i in range(20):
    sm.append(s.subs(k, i).n())
print(sm)

解法二差分方程的特征根解法

# 程序文件 ex3_1_2.py
import sympy as sp
sp.var('t, c1, c2')
t0 = sp.solve(t**2-t-1)         # 求解特征方程
eq1 = c1+c2-1
eq2 = c1*t0[0]+c2*t0[1]-1
s = sp.solve([eq1, eq2]) 
print('c1=', s[c1])
print('c2=', s[c2])

解法三直接求解

# 程序文件 ex3_1_3.py
import sympy as sp
sp.var('k')
y = sp.Function('y')
f = y(k+2)-y(k+1)-y(k)
s = sp.rsolve(f, y(k), {y(0):1, y(1):1})
print(s)

3.1.2 Leslie 种群模型

$x^{\left( k \right)}=\textit{\textbf{L}}x^{\left( k-1 \right)}=\cdots=\textit{\textbf{L}}^{\left( k \right )}x^{\left ( 0 \right)}, k=1,2,\cdots .$

其中 $\textit{\textbf {L}}=\begin{bmatrix} a_{1} & a_{2} & \cdots & a_{n-1} & a_{n}\\ b_{1} & 0 & \cdots & 0 & 0\\ 0 & b_{2} & \cdots & 0 & 0\\ \vdots & \vdots & & \vdots & \vdots\\ 0 & 0 & \cdots & b_{n-1} & 0 \end{bmatrix}$

例 3.2 某种动物雌性的最大生存年龄为 15 年，以 5 年为间隔，将动物种群分为 3 个年龄组 [0, 5), [5, 10), [10, 15)，利用统计资料，已知 $a_{1}=0,a_{2}=4,a_{3}=3;b_{1}=0.5,b_{2}=0.25$ ，在初始时刻 $t=0$ 时，3 个年龄组的雌性动物个数分别为500，1000，500，对整个种群总量进行估计。

# 程序文件 ex3_2.py
import numpy as np
import sympy as sp
X0 = np.array([500, 1000, 500])
L = np.array([[0, 4, 3], [0.5, 0, 0], [0, 0.25, 0]])
X1 = L @ X0
X2 = L @ X1     # @ 表示矩阵乘法
X3 = L @ X2
Ls = sp.Matrix([[0, 4 , 3], 
    [sp.Rational(1, 2), 0, 0], 
    [0, sp.Rational(1, 4), 0]])  # 定义矩阵
sp.var('lamda')                  # 定义符号变量
p = Ls.charpoly(lamda)           # 计算特征多项式
w1 = sp.roots(p)                 # 计算特征值
w2 = Ls.eigenvals()              # 直接计算特征值
v = Ls.eigenvects()              # 直接计算特征向量
print('特征值为:', w2)
print('特征向量为:\n', v)
P, D = Ls.diagonalize()          # 相似对角化
Pinv = P.inv()                   # 求逆阵
Pinv = sp.simplify(Pinv)
cc = Pinv @ X0
print('P=\n', P)
print('c=', cc[0])

3.1.3 PageRank 算法

1. 基本原理

Google 搜索引擎对检索结果的一种排序算法，基本思想主要来自于文献引文分析。即文献的质量和重要性由其被引数量和引文质量衡量。

Google 给出页面排序遵循两条标准：（1）链接数目（2）链出页面的重要程度

PageRank 利用网络的拓扑结构判断网页的重要性。

（1）邻接矩阵 $\textit{\textbf{W}}=\left( \omega_{ij} \right )_{N \times N}$ ，若 $i$ 到 $j$ 有超链接，则 $\omega_{ij}=1$ ，否则为0。

（2）记 $\textit{\textbf{W}}$ 的行和为 $r_{i}=\sum\limits_{j=1}^{N}\omega_{ij}$ ，即 $i$ 的链出链接数目。

（3）定义矩阵 $\textit{\textbf{P}}=\left(p_{ij} \right )_{N \times N}, p_{ij}=\frac{\omega_{ij}}{r_{i}} \ \ \ i, j=1, 2, \cdots, N$ ， $\textit{\textbf{P}}$ 为马尔可夫链的状态转移概率矩阵， $p_{ij}$ 为 $i$ 到 $j$ 的概率。

（4）马尔科夫链的平稳分布 $\textit{\textbf{x}}=\left [x_{1},x_{2},\cdots,x_{N} \right ]^{T}$ 满足 $\textit{\textbf{P}}^{T}\textit{\textbf{x}}=\textit{\textbf{x}},\ \ \sum\limits_{i=1}^{N}x_{i}=1$ ， $x_{k}=\sum\limits_{i=1}^{N}p_{ik}x_{i}=\sum\limits_{i}^{N}\frac{\omega_{ik}}{r_{i}}x_{i}$

平稳分布（PageRank 值）是状态转移概率矩阵 $\textit{\textbf{P}}$ 的转置矩阵 $\textit{\textbf{P}}^{T}$ 的最大特征值（=1）所对应的归一化特征向量。

例 3.3 已知一个 $N=6$ 的网络， $\textit{\textbf{W}}=\left [ \begin{matrix} 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1\\ 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1\\ 1 & 0 & 0 & 0 & 0 & 0 \end{matrix} \right ]$ ， $\textit{\textbf{P}}=\left[ \begin{matrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{3} & \frac{1}{3} & \frac{1}{3}\\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 & 0 \end{matrix} \right]$ ，求它的 PageRank 取值

# 程序文件 ex3_3.py
import numpy as np
from scipy.sparse.linalg import eigs
import pylab as plt
L = [(1, 2), (2, 3), (2, 4), (3, 4), (3, 5),
    (3, 6), (4, 1), (5, 6), (6, 1)]
w = np.zeros((6, 6))                        # 邻接矩阵初始化
for i in range(len(L)):
    w[L[i][0]-1, L[i][1]-1] = 1
r = np.sum(w, axis=1, keepdims=True)
P = w/r                                     # 这里利用矩阵广播
val, vec = eigs(P.T, 1)
V = vec.real
V = V.flatten()                             # 展开成 (n, ) 形式的数组
V = V/V.sum()
print('V=', np.round(V, 4)) 
plt.bar(range(1, len(w)+1), V, width=0.6, color='b')
plt.show()

2. 随机冲浪模型的 PageRank 值

用公式 $\tilde{\textit{\textbf{P}}}=\frac{\left(1-d \right )}{N}\textit{\textbf{e}}\textit{\textbf{e}}^{T}+d\textit{\textbf{P}}$ 表示转移概率矩阵，其中，

$\textit{\textbf{e}}$ 为分量全为 1 的 $N$ 维列向量

$\textit{\textbf{e}}\textit{\textbf{e}}^{T}$ 为全 1 矩阵

$d \in \left(0, 1 \right )$ 为阻尼因子（damping factor），Google 取 $d=0.85$

例 3.4 用随机冲浪模型计算例 3.3 所示有向图中各顶点的 PageRank 值

# 程序文件 ex3_4.py
import numpy as np
from scipy.sparse.linalg import eigs
import pylab as plt
L = [(1, 2), (2, 3), (2, 4), (3, 4), (3, 5),
    (3, 6), (4, 1), (5, 6), (6, 1)]
w = np.zeros((6, 6))
for i in range(len(L)):
    w[L[i][0]-1, L[i][1]-1] = 1
r = np.sum(w, axis=1, keepdims=True)
P = (1-0.85)/w.shape[0]+0.85*w/r            # 这里利用矩阵广播
val, vec = eigs(P.T, 1)
V = vec.real
V = V.flatten()                             # 展开成 (n, ) 形式的数组
V = V/V.sum()
print('V=', np.round(V, 4)) 
plt.bar(range(1, len(w)+1), V, width=0.6, color='b')
plt.show()

3.2 矩阵的奇异值分解及应用

3.2.1 矩阵的奇异值分解

定理 3.1 设 $\textit{\textbf{A}}$ 是一个秩为 $r$ 的 $m \times n$ 矩阵，则存在正交矩阵 $\textbf{\textit{U}}$ 和 $\textbf{\textit{V}}$ ，使得

$\textit{\textbf{U}}^{T}\textit{\textbf{A}}\textit{\textbf{V}}=\left[\begin{matrix} \sum & 0 \\ 0 & 0 \end{matrix} \right ]$ （3.11）

式中： $\sum=diag\left \{ \sigma_{1},\sigma_{2},\cdots,\sigma_{r} \right \}$ ， $\sigma_{1} \geqslant \sigma_{2} \geqslant \cdots \geqslant \sigma_{r} \geqslant 0, \ \sigma_{1},\sigma_{2},\cdots,\sigma_{r}$ 为矩阵 $\textit{\textbf{A}}^{T}\textit{\textbf{A}}$ 对应的正特征值。称

$\textit{\textbf{A}}=\textit{\textbf{U}}\left[ \begin{matrix} \sum & 0& \\ 0 & 0 & \end{matrix} \right ]\textit{\textbf{V}}^{T}$ （3.12）

是 $\textit{\textbf{A}}$ 的奇异值分解， $\sigma_{i} \left( i=1, 2,\cdots,r \right )$ 称为 $\textit{\textbf{A}}$ 的奇异值。

矩阵的 $F\left(\text{Frobenius} \right )$ 范数定义为

$\left \| \textit{\textbf{A}} \right \|_{F}^{2}=\text{tr}\left(\textit{\textbf{A}}^{T}\textit{\textbf{A}} \right )=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}a_{ij}^{2}$ (3.13)

由于

$\text{tr}\left(\textit{\textbf{A}}^{T}\textit{\textbf{A}} \right ) = \text{tr} \left(\textit{\textbf{V}} \left[\begin{matrix} \sum & 0 \\ 0 & 0 \end{matrix} \right] ^{T} \textit{\textbf{U}}^{T} \textit{\textbf{U}} \left[\begin{matrix} \sum & 0 \\ 0 & 0 \end{matrix} \right] \textit{\textbf{V}}^{T} \right )$

$= \text{tr}\left(\textit{\textbf{V}} \left[\begin{matrix} \sum^{2} & 0 \\ 0 & 0 \end{matrix} \right] \textit{\textbf{V}}^{T} \right ) = \sum\limits_{i=1}^{r}\sigma_{i}^{2}$

所以

$\left \| \textit{\textbf{A}} \right \|_{F}^{2}=\sum\limits_{i=1}^{r}\sigma_{i}^{2}$ (3.14)

式 3.14 表明矩阵的 $F$ 范数的平方等于矩阵的所有奇异值的平方和。

例 3.5 试求矩阵 $\textit{\textbf{A}}=\left[\begin{matrix} 1 & 0 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 0 \end{matrix} \right ]$ 的奇异值分解

# 程序文件 ex3_5.py
import numpy as np 
from numpy.linalg import svd
a = np.array([[1, 0, 1], [0, 1, 1], [0, 0, 0]])
u, s, vt = svd(a)       # a = u @ np.diag(s) @ vt
print(u); print(s); print(vt)

3.2.2 奇异值的分解应用

例 3.6 推荐系统的评分

根据已有评分，给出为评分项的评分估计

表3.1 原始的打分数据（单位：分）
	叉烧粉肠	新疆手抓饭	四川火锅	粤式烧鹅饭	大盘鸡盘面	东北饺子	重庆辣子鸡	广东虾饺	剁椒鱼头	兰州拉面	烤羊排
丁一	5	2	1	4	0	0	2	4	0	0	0
刘二	0	0	0	0	0	0	0	0	0	3	0
张三	1	0	5	2	0	0	3	0	3	0	1
李四	0	5	0	0	4	0	1	0	0	0	0
王五	0	0	0	0	0	4	0	0	0	4	0
马六	0	0	1	0	0	0	1	0	0	5	0
陈七	5	0	2	4	2	1	0	3	0	1	0
胡八	0	4	0	0	5	4	0	0	0	0	5
赵九	0	0	0	0	0	0	4	0	4	5	0
钱十	0	0	0	4	0	0	1	5	0	0	0
孙甲	0	0	0	0	4	5	0	0	0	0	3
周乙	4	2	1	4	0	0	2	4	0	0	0
吴丙	0	1	4	1	2	1	5	0	5	0	0
郑丁	0	0	0	0	0	4	0	0	0	4	0
冯戊	2	5	0	0	4	0	0	0	0	0	0
储己	5	0	0	0	0	0	0	4	2	0	0
魏庚	0	2	4	0	4	3	4	0	0	0	0
高辛	0	3	5	1	0	0	4	1	0	0	0

解法一非压缩数据的模型一

相似系数：指标变量 $x_{1},x_{2},\cdots,x_{11}$ 分别表示叉烧肠粉、新疆手抓饭、...、烤羊排，用 $i=1,2,\cdots,18$ 分别表示用户丁一...高辛。记用户 $i\left( i=1,2,\cdots,18 \right )$ 关于指标变量 $x_{j} \left( j=1,2,\cdots,11 \right )$ 的观测值为 $a_{ij}$ , 构造数据矩阵 $\textit{\textbf{A}}=\left( a_{ij} \right )_{18\times 11}$ 。

计算 $x_{i},x_{j}$ 间的相关系数

$r_{ij}=\frac{\sum\limits_{k=1}^{18}\left(a_{ki}-\mu_{i} \right ) \left(a_{kj}-\mu_{j} \right )}{ \sqrt{ \sum\limits_{k=1}^{18} \left(a_{ki}-\mu_{i} \right )^{2} } \sqrt{ \sum\limits_{k=1}^{18} \left(a_{kj}-\mu_{j} \right )^{2} } }$ (3.15)

式中： $\mu_{j}=\frac{1}{18}\sum_{k=1}^{18}a_{kj}$

相关系数取值范围为 $\left[-1,1 \right ]$ ，可对其归一化，记 $\tilde{r}_{ij}=0.5+0.5r_{ij}$

评分估计：基于菜品之间相似度估计某顾客未购菜品评分。欲估计菜品为 $x_{k}$ ，已评分菜品为 $x_{j} \ \left(j=j_{1},j_{2},\cdots,j_{n_{k}} \right )$ ，评分分别为 $s_{j} \ \left(j=j_{1},j_{2},\cdots,j_{n_{k}} \right )$ ， $x_{k}$ 与 $x_{j}$ 相似度为 $\tilde{r}_{ij} \ \left(j_{1},j_{2},\cdots,j_{n_{k}} \right )$ ，利用相似度加权估计 $x_{k}$ 的评分

$\hat{s}_{k} = \frac{ \sum\limits_{i=1}^{n_{k}} \tilde{r}_{kj_{i}} s_{j_{i}}} {\sum\limits_{i=1}^{n_{k}} \tilde{r}_{kj_{i}} }$ （3.16）

# 程序文件 ex3_6_1.py
import numpy as np
import pandas as pd
a = np.loadtxt('data3_6_1.txt')
b = 0.5*np.corrcoef(a.T)+0.5                # 归一化相似度
c = pd.DataFrame(b)
c.to_excel('data3_6_2.xlsx', index=False)
print('请输入人员编号1-18')
user = int(input())
n = a.shape[1]                              # 变量的个数
no = np.where(a[user-1,:]==0)[0]            # 未评分编号
yb = set(range(n))-set(no)                  # 已评分编号
yb = list(yb)
ys = a[user-1, yb]                          # 已评分分数
sc = np.zeros(len(no))                      # 初始化
for i in range(len(no)):
    sim = b[no[i], yb]
    sc[i] = ys @ sim/sum(sim)
print('未评分项的编号为：', no+1)
print('未评分项的分数为：', np.round(sc, 4))

解法二基于奇异值分解压缩数据的模型二

基于奇异值分解对原始数据进行压缩。对 $\textit{\textbf{A}}$ 进行奇异值分解，得

$\textit{\textbf{A}}=\textit{\textbf{U}}_{18\times18}\left[ \begin{matrix} \sum_{11} \\ 0 \end{matrix} \right] \textit{\textbf{V}}_{11\times11}^{T}$

式中， $\sum_{11\times11}=\text{diag}\left(\sigma_{1},\sigma_{2},\cdots,\sigma_{11} \right )$ ，这里， $\sigma_{1},\sigma_{2},\cdots,\sigma_{11}$ 为矩阵 $\textit{\textbf{A}}$ 的正奇异值。

有 $\frac{\sum\limits_{i=1}^{6}\sigma_{i}^{2}}{\sum\limits_{i=1}^{11}\sigma_{i}^{2}}\geqslant0.9$ ，故只需 6 个奇异值，即可达到主成分贡献率的 90%。可通过压缩的方式，将原始分数矩阵的行由 18 维压缩到 6 维，避免稀疏矩阵的情况。压缩后的数据矩阵

$\textit{\textbf{B}}_{6\times11}=\text{diag}\left(\sigma_{1},\sigma_{2},\cdots,\sigma_{6} \right ) \textit{\textbf{U}} ^{T} \left[ :\ 6,: \ \right ]\textit{\textbf{A}}$

式中， $\textit{\textbf{U}} ^{T} \left[ :\ 6,: \ \right ]$ 为 $\textit{\textbf{U}} ^{T}$ 的前 6 行所组成的 6 行 8 列矩阵。

余弦相似度：两个指定向量的余弦相似度就是用二者夹角 $\theta$ 的余弦值 $\text{cos}\theta=\frac{\nu_{1}\cdot \nu_{2}}{\left |\nu_{1} \right |\cdot \left |\nu_{2} \right |}$ 表示。

然后，如解法一，进行归一化，并估计未购菜品的评分。

# 程序文件 ex3_6_1.py
import numpy as np
import pandas as pd
a = np.loadtxt('data3_6_1.txt')
u, sigma, vt = np.linalg.svd(a)
print(sigma)
cs = np.cumsum(sigma**2)
rate = cs/cs[-1]                                # 计算信息累积贡献率
ind = np.where(rate>=0.9)[0][0]+1               # ind 为奇异值的个数，满足信息提取率达到 90%
b = np.diag(sigma[:ind]) @ u.T[:ind, :] @ a     # 得到降维数据
c = np.linalg.norm(b, axis=0, keepdims=True)    # 逐列求范数
d = 0.5*b.T @ b/(c.T @ c)+0.5                   # 求相似度
#d = 0.5*np.corrcoef(b.T)+0.5
dd = pd.DataFrame(d)
dd.to_excel('data3_6_3.xlsx', index=False)
print('请输入人员编号1-18')
user = int(input())
n = a.shape[1]                                  # 变量的个数
no = np.where(a[user-1,:]==0)[0]                # 未评分编号
yb = set(range(n))-set(no)                      # 已评分编号
yb = list(yb)
ys = a[user-1, yb]                              # 已评分分数
sc = np.zeros(len(no))                          # 初始化
for i in range(len(no)):
    sim = d[no[i], yb]
    sc[i] = ys @ sim/sum(sim)
print('未评分项的编号为：', no+1)
print('未评分项的分数为：', np.round(sc, 4))

例 3.7 利用 SVD 进行图像压缩

用矩阵 $\textit{\textbf{A}}$ 表示图像的 $\textit{mn}$ 个像素，对其奇异值分解，得到 $\textit{\textbf{A}}= \textit{\textbf{U}}\textit{\textbf{S}}\textit{\textbf{V}}^{T}$ ，其中，奇异值按照从大到小的顺序。若从中选择 $k$ 个大奇异值以及这些奇异值对应的左和右奇异向量逼近原图像，便可以使用 $m\times k +k+k\times n=k\left(m+n+1 \right )$ 个数值代替原来的 $m \times n$ 个图像数据。这 $k\left(m+n+1 \right )$ 个被选择的新数据是矩阵 $\textit{\textbf{A}}$ 的前 $k$ 个奇异值、 $m \times m$ 左奇异向量矩阵 $\textit{\textbf{U}}$ 的前 $k$ 列和右奇异向量矩阵 $\textit{\textbf{V}}$ 的前 $k$ 列的元素。

压缩比率为 $\rho=\left(1-\frac{k\left(m+n+1 \right )}{mn} \right ) \times 100\%$

在接收端，接收到奇异值 $\sigma_{1},\sigma_{2},\cdots,\sigma_{k}$ 以及左奇异向量 $u_{1},u_{2},\cdots,u_{k}$ 和右奇异向量 $\nu_{1},\nu_{2},\cdots,\nu_{k}$ 后，即可通过截尾的奇异值分解公式重构原图像：

$\hat{\textit{\textbf{A}}}=\sum\limits_{i=1}^{k}\sigma_{i}u_{i}\nu_{i}^{T}$

# 程序文件 ex3_7.py
import numpy as np
from numpy import linalg as LA 
from PIL import Image
import pylab as plt                 # 加载 Matplotlib 的 pylab 接口
plt.rc('font', size=13) 
plt.rc('font',family='SimHei')
a = Image.open('Lena.bmp')          # 返回一个 PIL 图像对象
if a.mode!='L':
    a = a.convert('L')              # 转换为灰度图像
b = np.array(a).astype(float)       # 把图像对象转换为数组
[p, d, q] =  LA.svd(b)
m, n = b.shape
R = LA.matrix_rank(b)               # 图像矩阵的秩
plt.figure(0)
plt.plot(np.arange(1, len(d)+1), d, 'k.')
plt.ylabel('奇异值'); plt.xlabel('序号')
plt.title('图像矩阵的奇异值')
CR = []
for K in range(1, int(R/4), 10):
    plt.figure(K)
    plt.subplot(121)
    plt.title('原图')
    plt.imshow(b, cmap='gray')
    I = p[:, :K+1] @ (np.diag(d[:K+1])) @ (q[:K+1, :])
    plt.subplot(122)
    plt.title('图像矩阵的秩='+str(K))
    plt.imshow(I, cmap='gray')
    src = m*n
    compress = K*(m+n+1)
    ratio = (1-compress/src)*100    # 计算压缩比率
    CR.append(ratio)
    print('Rank=%d:K=%d个:ratio=%5.2f'%(R, K ,ratio))
plt.figure();
plt.plot(range(1, int(R/4), 10), CR, 'bo-')
plt.title('奇异值个数与压缩比率的关系')
plt.xlabel('奇异值个数')
plt.ylabel('压缩比率')
plt.show()