Numpy 练习

先设置随机数种子：

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() 返回的是将数组平铺后，最小值的下标。