【CME 193】 Numpy

最新推荐文章于 2021-04-18 11:28:07 发布

KiwiHe

最新推荐文章于 2021-04-18 11:28:07 发布

阅读量355

点赞数

分类专栏：【Python学习作业】

本文链接：https://blog.youkuaiyun.com/KiwiHe/article/details/80394421

版权

【Python学习作业】专栏收录该内容

17 篇文章

订阅专栏

这篇博客介绍了Numpy库的使用，包括利用np.random.normal生成随机矩阵，使用scipy.linalg.toeplitz创建 toeplitz 矩阵，通过numpy的基本矩阵操作如dot、A.T和eye，以及np.linalg.solve解决线性方程组。还讨论了矩阵的Frobenius范数、特征值计算（通过SVD和幂法）、奇异值计算和二项分布。最后提到了np.argmin函数的应用。

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

9 Numpy

Generate matrices A, with random Gaussian entries, B, a Toeplitz matrix, where A 2 Rnm and B 2 Rmm,
for n = 200, m = 500.

本题中使用np.random.normal 生成A，其原理为：

使用scipy.linalg.toeplitz 生成B，其原理为:

import numpy as np
from scipy.linalg import toeplitz
import time 

#N(0,1) 200*500
n = 200
m = 500
A = np.random.normal(loc = 0.0, scale = 1.0, size = (n,m))

#toeplitz 500*500
B = toeplitz(list(range(1,501)))

以下9-1~9-3程序均以此为基础实现：

本题只需使用numpy中基本矩阵操作 --dot求积，A.T求逆，eye生成单位矩阵等：

#9-1
print("9-1")
print("A + A = ")
print(A+A)

print("A * A^T = ")
print(np.dot(A, A.T))

print("A^T * A = ")
print(np.dot(A.T, A))

print("A * B = ")
print(np.dot(A, B))

print("A* (B - lamda * I) = ")
#eye -> I
lam = float(input("Please input a la:"))
tmp = B - lam * np.eye(500)
print(np.dot(A, tmp))

本题中使用np.linalg.solve求解线性方程组：

#9-2
print()
print("9-2")
b = np.array(list(range(1,501)))  
x = np.linalg.solve(B,b)  
print("x = ")
print(x)

本题中涉及求矩阵的Frobenius范数及无范数，经由不同参数调用np.linalg.norm即可，原理如下:

求最大/最小特征值时对矩阵进行SVD分解即可得到奇异值的对角矩阵，再对该矩阵求最大/最小元即可，原理如下：

#9-3
print()
print("9-3")
print("the Frobenius norm of A: ")  
A_fro = np.linalg.norm(A, 'fro')
print(A_fro)
print("the infinity norm of B: ")  
B_inf = np.linalg.norm(B, np.inf) 
print(B_inf)

#SVD for singular value
u, s, vt = np.linalg.svd(B,full_matrices=True)  
print("The largest singular value is:")
print(max(s))  
print("The smallest singular value is:")
print(min(s))

本题用幂法求特征值和特征向量，其原理见：

（https://baike.baidu.com/item/幂法求矩阵特征值/3082486?fr=aladdin）

精度设置为1e-10，返回值为：lam：主特征值；u：主特征向量；iter: 迭代次数；tim：使用时间

#9-4
print()
print("9-4")
for nx in[100,200,500]:
	print("When n = " + str(nx))
	start = time.clock()
	Z = np.zeros([nx*nx,1])  
	for idx in range(len(Z)):  
		Z[idx] = np.random.normal(1, 2)  
	Z = Z.reshape(nx, nx)  
	#Z = np.random.standard_normal((nx, nx)) 
	v = u = np.ones([nx, 1])  
	lam = 0
	iter = 0
	while(True):
		v = np.matmul(Z, u)  
		pre = lam  
		lam = v.max()  
		u = v / lam  
		if ( abs(pre - lam) < 1e-10 ):  
			break;  
		iter += 1  
	print('lambda:' + str(lam))  
	print('iter:' + str(iter))  
	end = time.clock()  
	tim = end - start
	print('time used: ' + str(tim))  
	print()

当n逐渐增大时，程序的迭代次数并不随之增大，可以看出该算法的收敛是稳定的：

本题中也涉及奇异值的计算，其原理与9.3类似，但要先生成C，其中涉及函数np.random.binomial，原理如下：

二项分布（binomial distribution）：

P (N) = (n N) p N (1 - p) n - N

numpy给出的api是：

numpy.random.RandomState.binomial(n, p, size=None)

表示对一个二项分布进行采样（size表示采样的次数，draw samples from a binomial distribution.），参数中的n, p分别对应于公式中的n,p，函数的返回值表示n

中成功（success）的次数（也即N）。

#9-5
print()
print("9-5")
for n in [50,100,200]:  
    for p in (0.2,0.5,0.8):  
        print("When n = %d, p = %.3f" % (n,p))  
        C = np.random.binomial(1,p,(n,n))  
        u, s, v = np.linalg.svd(C,full_matrices=True)  
        maxvalue = max(s)  
        print("The largest singular value is: " + str(maxvalue)) 
    print()

从程序输出结果可以看出，随着n/p的增大，矩阵的最大奇异值也逐渐增大：

按照本题的提示，我们使用np.argmin，其原理如下：

#9-6
print()
print("9-6")
print("Consume A is an array range in (1,100)")
A = np.array(range(1,100))
z = (float)(input("Please input a z:"))
D = abs(A - z)  
index = np.argmin(D) 
print("The element closed to z is: " + str(A[index]))

2018.05.21