先设置随机数种子:
import numpy as np
import time
seed = np.int64(time.time())
print(seed)
np.random.seed(seed)
Generate matrices A, with random Gaussian entries, B, a Toeplitz matrix, where A ∈Rn×m and B ∈Rm×m, for n = 200, m = 500.
A = np.random.normal(size=[n, m])
print(np.mean(A))
print(np.var(A))
以上代码可以生成一个n行m列的服从标准正态分布的矩阵,通过输出矩阵的均值和方差可以验证
def toeplitz(vectorX, vectorY):
toeplitz_matrix = np.zeros([len(vectorX), len(vectorY)])
for i in range(len(vectorX)):
for j in range(len(vectorY) - i):
toeplitz_matrix[j][j+i] = vectorY[i]
for j in range(i, len(vectorX)):
toeplitz_matrix[j][j-i] = vectorX[i]
return toeplitz_matrix
传进行向量和列向量,生成toeplitz矩阵,当行向量的第一个元素和列向量的第一个元素不一样时,矩阵第一行第一个元素取列向量第一个元素
vectorX = np.random.normal(size = [m])
vectorY = np.random.normal(size = [m])
B = toeplitz(vectorX, vectorY)
生成toeplitz矩阵
验证生成toeplitz矩阵函数的正确性:
生成一个4*4的toeplitz矩阵:
vectorX = np.random.normal(size=[4])
vectorY = np.random.normal(size=[4])
B = toeplitz(vectorX, vectorY)
print(vectorX)
print("___")
print(vectorY)
print("___")
print(B)
输出如图。可见函数正确
sum_A_A = A + A
mul_A_AT = np.matmul(A, A.transpose())
mul_AT_A = np.matmul(A.transpose(), A)
mul_A_B = np.matmul(A, B)
def mul_A_B_I(A, B, constant):
sub = np.eye([len(B), len(B)])
return np.matmul(A, B - constant*sub)
def solve_linsys(B, b):
x = np.linalg.solve(B, b)
print("The result is {}".format(np.allclose(np.dot(B, x), b)))
return x
测试代码:
b = np.random.random(size=[m])
solve_linsys(B, b)
输出如图,函数正确
inf_norm = np.linalg.norm(B, np.inf)
f_norm = np.linalg.norm(A, "fro")
_, s, _ = np.linalg.svd(B)
max_singular = np.max(s)
min_singular = np.min(s)
print(inf_norm)
print(f_norm)
print(max_singular)
print(min_singular)
输出如下:
def power_iteration(Z, stop):
u = np.random.normal(size=[len(Z)]).transpose()
u = u/np.max(u)
v = u.copy()
eigen_change = np.inf
eigen_val = np.max(v)
iteration = 0
s_time = time.clock()
while abs(eigen_change) > stop:
v = np.matmul(Z, u)
max_v = np.max(v)
eigen_change = eigen_val - max_v
eigen_val = max_v
u = v/max_v
iteration += 1
print("The result is {}".format(np.allclose(np.matmul(Z, u), eigen_val*u)))
print("Total iteration: {}".format(iteration))
print("Total time: {}s".format(time.clock() - s_time))
return [eigen_val, u, iteration, time.clock() - s_time]
传进矩阵以及阈值,当特征值变化量小于stop时停止迭代,返回结果
测试代码:
Z = np.random.normal(size=[20, 20])
stop = 1e-8
print(power_iteration(Z, stop))
随机生成的矩阵在迭代中可能出现不收敛现象,再加上计算量的考虑我尝试了20*20的矩阵
最后输出如下:
n = 1000
p_ = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
for p in p_:
C = np.random.random(size=[n, n])
for i in range(n):
for j in range(n):
C[i][j] = 1 if C[i][j] > p else 0
_, s, _ = np.linalg.svd(C)
max_singular = np.max(s)
print("n = {}, p = {}, max_singular = {}".format(n, p, max_singular))
输出如下:
在实验中,n, p, 以及最大奇异值max_singular的关系为 n*(1-p) ≈ max_singular
def nearest_neighbor(matrixA, val_z):
idx = np.argmin(np.abs(matrixA - val_z))
return matrixA[idx//np.shape(matrixA)[0]][idx%np.shape(matrixA)[1]]
思路:将矩阵的每个值减去目标值,然后取绝对值,最小元素的下标即是和目标值最接近的值的下标,然后根据下标直接返回矩阵中该元素即可
测试代码:
test_matrix = np.random.random(size=[50, 50])
res = nearest_neighbor(test_matrix, 0.5)
flag = True
for i in range(50):
for j in range(50):
if np.abs(test_matrix[i][j] - 0.5) < np.abs(res-0.5):
flag = False
if not flag:
print("The function is wrong")
else:
print("The function is right")
函数编写正确。特别要注意的是对一个多维数组使用np.argmin() 返回的是将数组平铺后,最小值的下标。