这里写自定义目录标题
练习
Ex1：利用列表推导式写矩阵乘法
一般的矩阵乘法根据公式，可以由三重循环写出：
In [138]: M1 = np.random.rand(2,3)
In [139]: M2 = np.random.rand(3,4)
In [140]: res = np.empty((M1.shape[0],M2.shape[1]))
In [141]: for i in range(M1.shape[0]):
.....: for j in range(M2.shape[1]):
.....: item = 0
.....: for k in range(M1.shape[1]):
.....: item += M1[i][k] * M2[k][j]
.....: res[i][j] = item
.....:
In [142]: ((M1@M2 - res) < 1e-15).all() # 排除数值误差
Out[142]: True
请将其改写为列表推导式的形式。
- 答案：
> res = [M1[i,:] @ M2[:, j] for i in range(M1.shape[0]) for j in range(M2.shape[1])]
> res = np.array(res).reshape(M1.shape[0], M2.shape[1])
- 参考答案：
# 从逐层的角度分析比较清晰
> res = [[sum([M1[i][k] * M2[k][j] for k in range(M1.shape[1])]) for j in range(M2.shape[1])] for i in range(M1.shape[0])]
Ex2：更新矩阵
设矩阵
A
m
×
n
A_{m×n}
Am×n ，现在对
A
A
A 中的每一个元素进行更新生成矩阵
B
B
B ，更新方法是
B
i
j
=
A
i
j
∑
k
=
1
n
1
A
i
k
B_{ij}=A_{ij}\sum_{k=1}^n\frac{1}{A_{ik}}
Bij=Aij∑k=1nAik1 ，例如下面的矩阵为
A
A
A ，则
B
2
,
2
=
5
×
(
1
4
+
1
5
+
1
6
)
=
37
12
B_{2,2}=5\times(\frac{1}{4}+\frac{1}{5}+\frac{1}{6})=\frac{37}{12}
B2,2=5×(41+51+61)=1237 ，请利用 Numpy 高效实现。
A
=
[
1
2
3
4
5
6
7
8
9
]
A=\left[\begin{array}{lll} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array}\right]
A=⎣⎡147258369⎦⎤
> B = A * np.sum(1/A,axis=1).reshape(-1,1)
> B
array([[ 0.54545455, 3.24324324, 7.91623037],
[ 2.18181818, 8.10810811, 15.83246073],
[ 3.81818182, 12.97297297, 23.7486911 ]])
Ex3：卡方统计量
设矩阵
A
m
×
n
A_{m\times n}
Am×n，记
B
i
j
=
(
∑
i
=
1
m
A
i
j
)
×
(
∑
j
=
1
n
A
i
j
)
∑
i
=
1
m
∑
j
=
1
n
A
i
j
B_{ij} = \frac{(\sum_{i=1}^mA_{ij})\times (\sum_{j=1}^nA_{ij})}{\sum_{i=1}^m\sum_{j=1}^nA_{ij}}
Bij=∑i=1m∑j=1nAij(∑i=1mAij)×(∑j=1nAij)，定义卡方值如下：
χ
2
=
∑
i
=
1
m
∑
j
=
1
n
(
A
i
j
−
B
i
j
)
2
B
i
j
\chi^2 = \sum_{i=1}^m\sum_{j=1}^n\frac{(A_{ij}-B_{ij})^2}{B_{ij}}
χ2=i=1∑mj=1∑nBij(Aij−Bij)2
请利用Numpy对给定的矩阵
A
A
A计算
χ
2
\chi^2
χ2
> np.random.seed(0)
> A = np.random.randint(10, 20, (8, 5))
> t1 = A.sum(1).reshape(-1, 1)
> t2 = A.sum(0).reshape(1, -1)
> B = t1 @ t2 / A.sum()
> np.sum(np.power(A-B, 2)/B)
11.842696601945802
Ex4：改进矩阵计算的性能
设 Z Z Z为 m × n m×n m×n的矩阵， B B B和 U U U分别是 m × p m×p m×p和 p × n p×n p×n的矩阵， B i B_i Bi为 B B B的第 i i i行， U j U_j Uj为 U U U的第 j j j列，下面定义 R = ∑ i = 1 m ∑ j = 1 n ∥ B i − U j ∥ 2 2 Z i j \displaystyle R=\sum_{i=1}^m\sum_{j=1}^n\|B_i-U_j\|_2^2Z_{ij} R=i=1∑mj=1∑n∥Bi−Uj∥22Zij，其中 ∥ a ∥ 2 2 \|\mathbf{a}\|_2^2 ∥a∥22表示向量 a a a的分量平方和 ∑ i a i 2 \sum_i a_i^2 ∑iai2。
现有某人根据如下给定的样例数据计算
R
R
R的值，请充分利用Numpy中的函数，基于此问题改进这段代码的性能。
In [145]: np.random.seed(0)
In [146]: m, n, p = 100, 80, 50
In [147]: B = np.random.randint(0, 2, (m, p))
In [148]: U = np.random.randint(0, 2, (p, n))
In [149]: Z = np.random.randint(0, 2, (m, n))
In [150]: def solution(B=B, U=U, Z=Z):
.....: L_res = []
.....: for i in range(m):
.....: for j in range(n):
.....: norm_value = ((B[i]-U[:,j])**2).sum()
.....: L_res.append(norm_value*Z[i][j])
.....: return sum(L_res)
.....:
In [151]: solution(B, U, Z)
Out[151]: 100566
- 答案：
> t1 = np.sum(np.power(B, 2), 1).reshape(-1,1)
> t2 = np.sum(np.power(U, 2), 0)
> t3 = 2 * B@U
> t = np.ones((m,n), dtype=np.int32)
> Y = t1*t + t2*t + t3
> R = np.sum(Y*Z)
> R
308970
- 参考答案
改进方法：
令 Y i j = ∥ B i − U j ∥ 2 2 Y_{ij} = \|B_i-U_j\|_2^2 Yij=∥Bi−Uj∥22，则 R = ∑ i = 1 m ∑ j = 1 n Y i j Z i j \displaystyle R=\sum_{i=1}^m\sum_{j=1}^n Y_{ij}Z_{ij} R=i=1∑mj=1∑nYijZij，这在Numpy中可以用逐元素的乘法后求和实现，因此问题转化为了如何构造Y矩阵。
Y i j = ∥ B i − U j ∥ 2 2 = ∑ k = 1 p ( B i k − U k j ) 2 = ∑ k = 1 p B i k 2 + ∑ k = 1 p U k j 2 − 2 ∑ k = 1 p B i k U k j \begin{aligned} Y_{i j} &=\left\|B_{i}-U_{j}\right\|_{2}^{2} \\ &=\sum_{k=1}^{p}\left(B_{i k}-U_{k j}\right)^{2} \\ &=\sum_{k=1}^{p} B_{i k}^{2}+\sum_{k=1}^{p} U_{k j}^{2}-2 \sum_{k=1}^{p} B_{i k} U_{k j} \end{aligned} Yij=∥Bi−Uj∥22=k=1∑p(Bik−Ukj)2=k=1∑pBik2+k=1∑pUkj2−2k=1∑pBikUkj
从上式可以看出，第一第二项分别为 B B B的行平方和与 U U U的列平方和，第三项是两倍的内积。因此， Y Y Y矩阵可以写为三个部分，第一个部分是 m × n m×n m×n的全 1 1 1矩阵每行乘以 B B B对应行的行平方和，第二个部分是相同大小的全 1 1 1矩阵每列乘以 U U U对应列的列平方和，第三个部分恰为 B B B矩阵与 U U U矩阵乘积的两倍。从而结果如下：
# t1, t2, t3可以用广播的方法直接做加法准算
# 以(B**2).sum(1)代替np.sum(np.power(B, 2),1)更为简略
> (((B**2).sum(1).reshape(-1,1) + (U**2).sum(0) - 2*B@U)*Z).sum()
308970
Ex5：连续整数的最大长度
输入一个整数的Numpy数组，返回其中递增连续整数子数组的最大长度，正向是指递增方向。例如，输入[1,2,5,6,7]，[5,6,7]为具有最大长度的连续整数子数组，因此输出3；输入[3,2,1,2,3,4,6]，[1,2,3,4]为具有最大长度的连续整数子数组，因此输出4。请充分利用Numpy的内置函数完成。（提示：考虑使用nonzero, diff函数）
- 参考答案
# 思路：求出每个连续整数子数组的上下边界，做diff，求长度最大值
> f = lambda x:np.diff(np.nonzero(np.r_[1,np.diff(x)!=1,1])).max()
心得体会
- 这个学习项目的题目质量很好，确实比较费脑。很多时候都自己无法相处最优解，对于numpy掌握还不够。做练习过程中，也发现自己线代水平有待提高！
- 第一次打卡！接下来也要认真完成每一次作业！