exercise 9 Numpy
事先声明,这里用到函数前提是安装好numpy和scipy包并在开头引用如下:
from numpy import *
from scipy.linalg import toeplitz
import time
Generate matrices A, with random Gaussian entries, B, a Toeplitz matrix, where A ∈Rn × m and B ∈Rm × m, for n = 200, m = 500.
用随机数构造矩阵可以用numpy下的random.normal(size)函数构造符合标准正态分布的随机数的矩阵,size是要构造的随机列表的大小,如下几个例子
A=mat(random.normal(size=(200,500)))
b0=mat(random.normal(size=500))
上面一个是构造一个2维的列表,行数是200,列数是500,下面是构造一个一维的列表,长度为500.
而构造toeplitz矩阵可以用scipy.linalg下的toeplitz(c,r)函数,这个函数接收两个一维的列表,以前者c为第一行,后者r为第二行构造tooeplitz矩阵,如果两个列表c,r的第一个数不同,默认使用前者c的第一个数。
构造A,B的代码如下:
A=mat(random.normal(size=(200,500)))
b0=mat(random.normal(size=500))
b1=mat(random.normal(size=500))
B=toeplitz(b0,b1)
练习9-1
Matrix operations
Calculate A + A, AAT,ATA and AB. Write a function that computes A(B − λI) for any λ.
这里加减法用numpy下的矩阵加减法就可以了,而乘法就需要用numpy的dot(A,B)函数完成,A,B分别代表要相乘的两个矩阵。转置使用A.T即可完成。需要注意做乘法时要满足矩阵乘法的规则。而I单位矩阵可以用eye(size)函数构造,size代表的是单位矩阵的是size*size大的。
代码如下
def fun_e1(A,B,lamda):
return dot(A,B-lamda*eye(500))
def e1(A,B):
print("A+A=",A+A)
print("A*At=",dot(A,A.T))
print("At*A=",dot(A.T,A))
lamda=5
ans = fun_e1(A,B,lamda)
print("A(B-lamda*I)=",ans)
return ans
练习9-2
Solving a linear system
Generate a vector b with m entries and solve Bx = b.
这题直接用numpy里的linalg的solve(B,b)函数就能解决,代码如下
def e2(B):
b=mat(random.normal(size=(500,1)))
print("b:",b)
x=linalg.solve(B,b)
print("Bx=b,the x is:",x)
练习9-3
Norms
Compute the Frobenius norm of A and the infinity norm of B. Also find the largest and smallest singular values of B.
这题用numpy里的linalg的norm(A,choose)函数就可以了,A代表要求范数的矩阵,无choose或choose为'fro'时,是计算frobenius范数,choose为numpy.inf时求无穷范数,为2求最大奇异值,-2求最小奇异值。代码如下:
def e3(A,B):
print("Frobenius norm of A",linalg.norm(A,'fro'))
print("Infiinity norm of B",linalg.norm(B,inf))
print("Largest singular of B",linalg.norm(B,2))
print("Smallest singular of B",linalg.norm(B,-2))
练习9-4
Exercise 9.4: Power iteration
Generate a matrix Z, n × n, with Gaussian entries, and usethe power iteration to find the largest
eigenvalue and corresponding eigenvector of Z. How many iterations are neededtill convergence?
Optional: use the time.clock() method to compare computation time when varying n.
生成矩阵可以用前面用过的random.normal(200,200)的函数,算largest eigenvalue和对应的eigenvector和迭代可以用幂迭代法来计算,所以需要生成全零和全一矩阵,可以分别用numpy里的zeros()和ones()函数来做用相邻两次迭代的向量的范数的差的绝对值作为判停条件,这里认定当它小于0.00005时为收敛,同时这里还要算迭代的次数,用time.clock()算经历了多少时间。代码如下:
def e4():
times=0
Z=random.standard_normal((200,200))
Z0=ones(200)
Z1=zeros(200)
Znorm=0
sta=time.clock()
while(1):
Z1=dot(Z,Z0)
tem=Znorm
Znorm=linalg.norm(Z1)
Z0=Z1/Znorm
times+=1
if(abs(Znorm-tem)<0.00005):
break
sto=time.clock()
print("Largest eigenvalue",Znorm)
print("Corresponding eigenvector",Z0)
print("Iterations",times)
print("Time needed",sto-sta)
练习9-5
Exercise 9.5: Singular values
Generate an n × n matrix, denoted by C, where each entry is 1 withprobability p and 0 otherwise. Use
the linear algebra library of Scipy to compute the singular values of C. What can you say about the relationship between n, p and the largest singular value?
C矩阵的生成可以用numpy里的binomial(1,p,(200,200))函数生成最大最小奇异值上面已经说了,这里直接给出代码。
def e5(): p=0.5 C=random.binomial(1,p,(200,200)) print("Largest singular of C",linalg.norm(C,2)) print("Smallest singular of C",linalg.norm(C,-2)) print("The relationship between largest singular and n,p,n*p=",200*p,"≈largestsingular=",linalg.norm(C,2))
练习9-6
Exercise 9.6: Nearest neighbor
Write a function that takes a value z and an array A and finds the element in A that is closest to z. The
function should return the closest value, not index.
Hint: Use the built-in functionality of Numpy rather than writing code to findthis value manually. In particular, use brackets and argmin.
首先生成数组A,也是用到random.normal的函数。然后可以用A[A>z]和A[A<=z]以z为界限分开A为两部分,然后可以用argmax()和argmin()找数组里最大数最小数对应的索引,通过一系列条件判断就可以找出来与z相差最小的部分了。代码:
def e6(z): A = random.normal(size=100000) AD=A[A>z] AX=A[A<=z] print(AD,AX) lea,lar=0,0 if(len(AD)>0): lar=argmin(AD) else: aaa= AX[argmax(AX)] print("Closest value:", aaa) return if(len(AX)>0): lea=argmax(AX) else: aaa= AD[lar] print("Closest value:", aaa) return if(abs(AD[lar]-z)>abs(AX[lea]-z)): aaa= AX[lea] print("Closest value:", aaa) return else: aaa=AD[lar] print("Closest value:",aaa) return
全部的代码:
from numpy import * from scipy.linalg import toeplitz import time def fun_e1(A,B,lamda): return dot(A,B-lamda*eye(500)) def e1(A,B): print("A+A=",A+A) print("A*At=",dot(A,A.T)) print("At*A=",dot(A.T,A)) lamda=5 ans = fun_e1(A,B,lamda) print("A(B-lamda*I)=",ans) return ans def e2(B): b=mat(random.normal(size=(500,1))) print("b:",b) x=linalg.solve(B,b) print("Bx=b,the x is:",x) def e3(A,B): print("Frobenius norm of A",linalg.norm(A,'fro')) print("Infiinity norm of B",linalg.norm(B,inf)) print("Largest singular of B",linalg.norm(B,2)) print("Smallest singular of B",linalg.norm(B,-2)) def e4(): times=0 Z=random.standard_normal((200,200)) Z0=ones(200) Z1=zeros(200) Znorm=0 sta=time.clock() while(1): Z1=dot(Z,Z0) tem=Znorm Znorm=linalg.norm(Z1) Z0=Z1/Znorm times+=1 if(abs(Znorm-tem)<0.00005): break sto=time.clock() print("Largest eigenvalue",Znorm) print("Corresponding eigenvector",Z0) print("Iterations",times) print("Time needed",sto-sta) def e5(): p=0.5 C=random.binomial(1,p,(200,200)) print("Largest singular of C",linalg.norm(C,2)) print("Smallest singular of C",linalg.norm(C,-2)) print("The relationship between largest singular and n,p,n*p=",200*p,"≈largestsingular=",linalg.norm(C,2)) def e6(z): A = random.normal(size=100000) AD=A[A>z] AX=A[A<=z] lea,lar=0,0 if(len(AD)>0): lar=argmin(AD) else: aaa= AX[argmax(AX)] print("Closest value:", aaa) return if(len(AX)>0): lea=argmax(AX) else: aaa= AD[lar] print("Closest value:", aaa) return if(abs(AD[lar]-z)>abs(AX[lea]-z)): aaa= AX[lea] print("Closest value:", aaa) return else: aaa=AD[lar] print("Closest value:",aaa) return A=mat(random.normal(size=(200,500))) b0=mat(random.normal(size=500)) b1=mat(random.normal(size=500)) B=toeplitz(b0,b1) e1(A,B) e2(B) e3(A,B) e4() e5() z=2 e6(z)
截图: