本次作业是CME 193 Introduction to Python一书中关于Numpy库的练习部分。
第一题:
为了完成本题,我们需要使用的Numpy或者Scipy函数有:
1.使用np.randorm.normal来生成具有高斯分布的A矩阵的元素。
2.使用np.zeros函数来生成指定大小的零矩阵。
3.使用scipy.linalg.toeplitz函数来生成拓普里兹矩阵B。
4.矩阵的乘法使用np.matmul来计算
5.矩阵的转置通过np.array类的内置函数transpose来获得。
知道了需要用哪些函数后,可以立马写出本题的代码,如下所示:
from scipy.linalg import toeplitz
import numpy as np
def foo(A, B, lam):
dim_B = B.shape
return np.matmul(A, B - lam * np.identity(dim_B[0]))
n = 200
m = 500
#generate metrix A with random Gasussian entries(loc = 1, scale = 2)
A = np.zeros([n*m,1])
for idx in range(len(A)):
A[idx] = np.random.normal(1, 2)
A = A.reshape(n, m)
#generate Toeplitz matrix B, its first row equal A's first row
B = toeplitz(A[0, :])
#A+A
print(A+A)
#three multiplications
print(np.matmul(A, A.transpose()))
print(np.matmul(A.transpose(), A))
print(np.matmul(A, B))
#test the function foo
print( foo(A, B, 0.007) )
运行该代码后我们可以得到四次矩阵乘法运算的结果。由于结果数据量较大,并且每次的结果都是随机生成的,就不放在这篇文章里了。
第二题:
为了完成本题,我们需要使用的Numpy或者Scipy函数有:
1.使用np.linalg.solve函数来解线性方程组。
2.使用np.allclose来比较两个矩阵之间的各元素之间的相等关系。我们使用这个函数来验证所求得的线性方程组的解是否正确。
本题的代码如下:
from scipy.linalg import toeplitz
import numpy as np
n = 200
m = 500
#generate metrix A with random Gasussian entries(loc = 1, scale = 2)
A = np.zeros([n*m,1])
for idx in range(len(A)):
A[idx] = np.random.normal(1, 2)
A = A.reshape(n, m)
#generate Toeplitz matrix B, its first row equal A's first row
B = toeplitz(A[0, :])
#generate vector x with the second row of A, and get b with Bx
x = A[1,:]
b = np.matmul(B, x)
x = np.linalg.solve(B, b)
print(np.allclose(np.matmul(B, x), b))
运行该代码,输出如下:
True
证明我们使用以上的方法求解出的方程解是正确的。
第三题:
为了完成本题,我们需要使用的Numpy或者Scipy函数有:
1.使用np.linalg.norm函数来获得矩阵或向量的范数。通过该函数的第二个参数,我们可以设定求各种范数。例如,使用np.linalg.norm(A, 'fro')我们就可以求得A矩阵的F范数。
2.使用scipy.linalg.svdvals来获得某个矩阵的奇异值。
3.使用np.array类的内置函数max和min来获得矩阵或者向量中的最大或最小值。在这里,我们使用这个函数来获得B矩阵的最大和最小奇异值。
本题代码如下:
from scipy.linalg import toeplitz
from scipy.linalg import svdvals
import numpy as np
n = 200
m = 500
#generate metrix A with random Gasussian entries(loc = 1, scale = 2)
A = np.zeros([n*m,1])
for idx in range(len(A)):
A[idx] = np.random.normal(1, 2)
A = A.reshape(n, m)
#generate Toeplitz matrix B, its first row equal A's first row
B = toeplitz(A[0, :])
#get different kinds of norms of A and B
A_norm_F = np.linalg.norm(A, 'fro')
B_norm_inf = np.linalg.norm(B, np.inf)
print(A_norm_F)
print(B_norm_inf)
#
print(svdvals(B).max())
print(svdvals(B).min())
运行该代码,我们可以得到输出:
707.240508048
968.495347708
551.073144544
0.151395334649
从上到下分别是:A的F范数、B的∞范数,以及B的最大和最小奇异值。
第四题:
本题使用power iteration来求解矩阵的最大特征值以及特征矩阵。
该算法的解法与有效性的证明请见:power_method
按照题意,我把n在[5, 20, 100, 200, 500, 1000]中变换,并记录求解的时间、迭代次数,在结果中进行输出。
本题的代码如下:
from scipy.linalg import toeplitz
import numpy as np
import time
for n in [5, 20, 100, 200, 500, 1000]:
s_time = time.clock()
#generate metrix A with random Gasussian entries(loc = 1, scale = 2)
Z = np.zeros([n*n,1])
for idx in range(len(Z)):
Z[idx] = np.random.normal(1, 2)
Z = Z.reshape(n, n)
v = u = np.ones([n, 1])
lam = 0
iter = 0
while(True):
v = np.matmul(Z, u)
pre_lam = lam
lam = v.max()
u = v / lam
if ( abs(pre_lam-lam) < 1e-100 ):
break;
iter += 1
#print(u)
print('lambda:' + str(lam))
print('iter:' + str(iter))
e_time = time.clock()
lapse = e_time - s_time
print('n: ' + str(n) + ' with time: ' + str(lapse))
输出如下:
lambda:7.5188091439
iter:51
n: 5 with time: 0.0025315535553157094
lambda:14.0884754581
iter:161
n: 20 with time: 0.003170367865388353
lambda:96.9368425841
iter:25
n: 100 with time: 0.02972639626556344
lambda:202.634475773
iter:20
n: 200 with time: 0.1407568270563342
lambda:496.67729975
iter:20
n: 500 with time: 0.7963087041511473
lambda:997.340717821
iter:20
n: 1000 with time: 2.3361704010505475
可以看出,随着n的增加:
1.迭代次数不一定增加,即算法的收敛速度稳定不变。
2.算法运行的时间呈平方上升的趋势。
第五题:
为了完成本题,我们需要使用的Numpy或者Scipy函数有:
1.使用np.randorm.sample函数来获取[0,1)之间的一个随机值。本题中我使用这个方法来使得我们可以在p的概率下给A的元素赋值为1。
2.使用scipy.linalg.svdvals来获得某个矩阵的奇异值。
3.使用np.array类的内置函数max和min来获得矩阵或者向量中的最大或最小值。在这里,我们使用这个函数来获得B矩阵的最大和最小奇异值。
根据题意,我把n在[3, 10, 50, 100, 500]之间变化,同时把p在[0, 0.2, 0.4, 0.6, 0.8, 1]之间变化,观察C的最大特征值情况。
本题代码如下:
from scipy.linalg import svdvals
import numpy as np
for n in [3, 10, 50, 100, 500]:
for p in [0, 0.2, 0.4, 0.6, 0.8, 1]:
#generate metrix C
C = np.zeros([n*n, 1])
for idx in range(len(C)):
C[idx] = 1 if np.random.sample(1) < p else 0
lsv = svdvals(C).max()
print('n:' + str(n) + " p:" + str(p) + '\n\t' + 'largest value:' + str(lsv))
输出如下:
n:3 p:0
largest value:0.0
n:3 p:0.2
largest value:1.41421356237
n:3 p:0.4
largest value:2.2360679775
n:3 p:0.6
largest value:2.44948974278
n:3 p:0.8
largest value:2.44948974278
n:3 p:1
largest value:3.0
n:10 p:0
largest value:0.0
n:10 p:0.2
largest value:4.472135955
n:10 p:0.4
largest value:5.83095189485
n:10 p:0.6
largest value:7.81024967591
n:10 p:0.8
largest value:8.60232526704
n:10 p:1
largest value:10.0
n:50 p:0
largest value:0.0
n:50 p:0.2
largest value:23.2808934536
n:50 p:0.4
largest value:31.9530906173
n:50 p:0.6
largest value:38.4447655735
n:50 p:0.8
largest value:45.0444225182
n:50 p:1
largest value:50.0
n:100 p:0
largest value:0.0
n:100 p:0.2
largest value:44.9110231458
n:100 p:0.4
largest value:62.6976873577
n:100 p:0.6
largest value:77.2075125878
n:100 p:0.8
largest value:89.5823643358
n:100 p:1
largest value:100.0
n:500 p:0
largest value:0.0
n:500 p:0.2
largest value:223.45916853
n:500 p:0.4
largest value:316.572266631
n:500 p:0.6
largest value:387.565478339
n:500 p:0.8
largest value:447.269493706
n:500 p:1
largest value:500.0
可以总结出其中的规律:
1.当C为全0矩阵的时候,其最大特征值为0。
2.当C为全1矩阵的时候,其最大特征值为n。
3.当C中1占到p,0占到(1-p)的时候,其最大特征值随着p的增大而增大。
第六题:
为了完成本题,我们需要使用的Numpy或者Scipy函数有:
1.使用np.argmin来获取array中最小值的坐标。
2.使用abs来获取绝对值。
为了验证所写函数的正确性,我设置A = [1,3,7,9],z = 4.5。这里,距离z最近的A值为3。
本题代码如下:
import numpy as np
def nearest_neighbor(A, z):
dist = abs(A-z)
min_arg = np.argmin(dist)
return A[min_arg]
A = np.array([1,3,7,9])
z = 4.5
print(nearest_neighbor(A, z))
输出如下:
3
证明了我们函数的正确性。