【CME 193】 Numpy

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(