numpy实战(高级编程技术week 11)
做题前
在这里,我使用了这样的代码来产生符合题目条件的矩阵
import numpy as np
from scipy.linalg import toeplitz
A = np.mat(np.random.normal(size=(200,500)))
B = np.mat(toeplitz(np.random.normal(size=500),np.random.normal(size=500)))
A,B
第一题
这里还是一些比较简单正常的计算,代码如下
a1 = A + A
a2 = A.dot(A.T)
a3 = A.T.dot(A)
a4 = A.dot(B)
a1,a2,a3,a4
结果如下:
.png)
同时,根据题目要求,我还写了这样的一个函数,并且测试了一下。
def t1(lamda):
return A.dot(B-lamda*np.eye(500))
t1(4)
结果如图所示
第二题
首先使用np.random.normal
函数生成这样的向量,为了求解该方程,可以考虑方程左右两边同时乘上矩阵
的逆,即
。
b = np.random.normal(size=(500,1))
x = B.I.dot(b)
x
结果如下
问题三:范数
这一道题主要考虑的是np.linalg.norm
函数的使用。
查阅文档得知,该函数的原型为numpy.linalg.norm(x, ord=None, axis=None, keepdims=False)
同时参数的功能如下表:
由此,可以写出这道题的代码
np.linalg.norm(A, ord='fro') # Frobenius norm
np.linalg.norm(B, ord=np.inf) # infinity norm
np.linalg.norm(B, ord=2) # biggest singular values
np.linalg.norm(B, ord=-2) # smallest singular values
结果可见
问题四
关于这个幂迭代法求解矩阵的特征向量,我使用了以下代码来实现。
def power_iteration(A, error):
v_k = np.random.rand(A.shape[1])
v_k1 = np.random.rand(A.shape[1])
count = 0
while np.linalg.norm(v_k - v_k1) >= error:
v_k = v_k1
u_k1 = np.dot(A, v_k)
m_k1 = np.max(u_k1)
v_k1 = u_k1/m_k1
count += 1
return (np.max(u_k1), v_k1, count)
power_iteration(np.array([[-4, 14, 0], [-5, 13, 0], [-1, 0, 2]]), 0.00000001)
#Z = np.random.normal(size=(200,200))
#power_iteration(Z, 0.000000001)
在进一步的测试中,验证了该程序的正确性。
我计算了 该矩阵的特征值和特征向量,与自己算的值是一致的。见下图:
.png)
问题五
这里需要探寻
,
,还有最大奇异值之间的关系。从scipy
的文档中找到关于求奇异值的函数:scipy.linalg.svd(C)
,其中第二个返回值就是奇异值列表。
from scipy import linalg
def get_singular_value(size, p):
t = np.random.rand(size, size) > p
C = np.zeros((size,size))
C[t]=1
s_v = linalg.svd(C)[1]
l_s_v = np.max(s_v)
return l_s_v
get_singular_valut(200,0.6)
t = []
for n in range(5,300):
t.append(get_singular_value(n,0.6))
plt.plot(range(5,300), t)
t = []
for n in np.linspace(0,1,30):
t.append(get_singular_value(100,n))
plt.plot(np.linspace(0,1,30),t)
plt.show()
通过实验,得知,当n渐渐变大时,最大奇异值也在变大
当p变小的时候,最大奇异值变小,如下图
问题六
这里要写一个函数,去寻找数据A中最接近 的数,函数需要返回最接近的值,以下是代码实现。
def find_nearest(A, z):
return A[np.argmin(np.abs(A-z))]
以下是运行截图: