高级编程技术 Numpy课后习题

0. Generate matrices A, with random Gaussian entries, B, a Toeplitz matrix, where A ∈Rn×m and B ∈Rm×m, for n = 200, m = 500.

import numpy as np
from scipy.linalg import toeplitz 

#0__construct A snd B 
n = 10
m = 15

A = np.random.randn(n,m)
print(A)

r = np.random.randint(0,10,size = m)
c = np.random.randint(0,10,size = m)
B = toeplitz(r,c)
print(B)

为了操作方便，所以先用n = 10，m = 15，打印了A和B，利用了numpy里的randn和randint，randn时标准正态分布，randint是离散分部的整数值。toeplitz（r，c）是scipy里的一个函数，可用来生成toeplitz矩阵。如果r[0]和c[0]不同则得到矩阵的（0，0）位使用r[0]；

得到的矩阵如下：

         

1. Exercise 9.1: Matrix operations Calculate

A + A, AAT,ATA and AB. Write a function that computes A(B−λI) for any λ.

#1__Matrix operations
A1 = A+A
A2 = np.dot(A,A.T)
A3 = np.dot(A.T,A)
AB = np.dot(A,B)

def calculateFunc(A,B,c):
	temp = B - c*np.eye(m)
	return np.dot(A,temp)
	
A4 = calculateFunc(A,B,2)

print('A+A is:\n',A1)
print('A*AT is:\n',A2)
print('AT*A is:\n',A3)
print('A*B is:\n',AB)
print('A*(B-2*I) is:\n',A4)

numpy中dot支持矩阵乘法，A.T直接就是转置

2.Exercise 9.2: Solving a linear system

Generate a vector b with m entries and solve Bx = b.

#2__ Solving a linear system
b = np.random.randint(0,10,size = m)
x = np.linalg.solve(B,b)
print('In Bx = b, the x is:\n',x)

np.linalg.solve可以直接得到解。

3.Exercise 9.3: Norms Compute

the Frobenius norm of A: ||A||F and the infinity norm of B: ||B||∞. Also find the largest and smallest singular values of B.

#3__Norms
Af = np.linalg.norm(A, ord = 'fro')
Bi = np.linalg.norm(B, ord = np.inf)

U,S,V = np.linalg.svd(B,full_matrices=True)
maxval = max(S)
minval = min(S)

print('the frobrnius norm of A is: ',Af)
print('the infinity norm of B is: ',Bi)
print('the largest singular values of B is: ',maxval)
print('the smallest singular values of B is: ',maxval)

numpy中的svd可用来求svd分解，其中返回的对角矩阵s中的每个元素都是奇异值，u和v的每个向量是奇异值相对应的向量。

4.Exercise 9.4: Power iteration

Generate a matrix Z, n × n, with Gaussian entries, and use the power iteration to find the largest eigenvalue and corresponding eigenvector of Z. How many iterations are needed till convergence?
Optional: use the time.clock() method to compare computation time when varying n.

#4__ Power iteration 
import time

def PowerItera(A):
	u = np.random.randint(0,10,size = n)
	last = 1
	lam = 0
	num = 0
	start = time.clock()
	while abs(lam - last)>0.0001:
		last = lam
		num+=1
		v = np.dot(A, u)
		lam = 0
		for vx in v:
			if abs(lam)<abs(vx):
				lam = vx
		u = v/lam		
	end = time.clock()
	return u,lam,end - start,num
	
Z = np.random.randn(n,n)
u,lam,itetime,num = PowerItera(Z)
print('the largest eigenvalue is: ', lam)
print('the corresponding eigenvector is: ', u)
print('the number of iterations is: ', num)
print('the time of iterations is: ',itetime)

代码是根据数值课本上幂法算法所得到的

当n=10，m=15时

经过63次迭代就收敛了

当n=200，m=500时

17万次就收敛了，可见迭代速度很快

5.Exercise 9.5: Singular values

Generate an n×n matrix, denoted by C, where each entry is 1 with probability 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?

#5__Singular values 
def calSigularVal(C):
	u, s, vh = np.linalg.svd(C,full_matrices=True) 
	return max(s)

for n in range(5,50,5):
	for p in [0.1,0.3,0.5,0.7,0.9]:
		C = np.random.binomial(1, p, (n,n))
		print(calSigularVal(C),end = ' ')
	print()

随着横坐标的增长p在增长，随着纵坐标的增长n在增长，由截图可见最大的奇异值也在随着n和p的增长不断地增长，与n*p成正比。

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 find this value manually. In particular, use brackets and argmin.

#6__Nearest neighbor
def nearNeighbor(A,z):
	B = A - z
	i = np.argmin(np.abs(B)) 
	return i