numpy习题答案（仅供参考）

Generate matrices A, with random Gaussian entries（高斯分布）, B, a Toeplitz matrix（Toeplitz矩阵介绍）, where A ∈ Rn×m and B ∈ Rm×m, for n = 200, m = 500.

n = 200
m = 500
A = random.randn(n,m)  #生成符合高斯分布的 n * m 矩阵
#print(A)
def toeplitz(row_vector,column_vector):  #通过一个行向量和一个列向量生成Toeplitz 矩阵
    tar = row_vector
    cur_row = row_vector
    for i in range(1,column_vector.size):
        cur_row = numpy.roll(cur_row, 1)
        cur_row[0][0] = column_vector[i, 0]
        tar = numpy.vstack((tar, cur_row))
    return tar
vec = random.randint(1,100,m)
vec1 = vec.reshape((1,m))
vec2 = vec.reshape((m,1))
B = toeplitz(vec1, vec2)
print(B)

Exercise 9.1: Matrix operations

Calculate A + A, AA⊤, A⊤A and AB. Write a function that computes A(B − λI) for any λ.

print(A+A) #向量或矩阵的加法可以直接将向量或矩阵相加，前提是大小匹配
print(numpy.dot(A,A.T))  #dot可以用来计算两个矩阵相乘， A.T表示A矩阵的转置矩阵
print(numpy.dot(A.T,A))
print(numpy.dot(A,B))
def fun1(matrix1, matrix2, cof): 
    I = numpy.eye(m,m)
    print(I)
    return numpy.dor(A, B - cof * I)

Exercise 9.2: Solving a linear system

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

from numpy import linalg
b = numpy.ones((m,1))
x = linalg.solve(B, b ) #linalg中的solve函数能够解方程组

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.

print(linalg.norm(A, 'fro'))  #norm函数的第一个参数是array_like类型的数据，第二个参数是范数类型
print(linalg.norm(B, numpy.inf))
print(linalg.norm(B, 2))
print(linalg.norm(B, -2))

norm函数的用法详见 链接 。

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.

def find_max(y):
    largest = abs(y[0, 0])
    tar = y[0, 0]
    for i in range(1, y.size):
        if abs(y[i, 0]) > largest:
            largest = abs(y[i, 0])
            ans = y[i, 0]
    return tar

def power_iteration(Z, E,n):
	start = time.process_time()
	X = random.rand(n, 1)
	count = 0
	Y = numpy.dot(Z, X)
	c_k1 = find_max(Y)
	c_k = c_k1 + E + 1
	while count == 0 or abs(c_k1 - c_k) >= E:
		X = Y / c_k1
		Y = numpy.dot(Z, X)
		c_k = c_k1
		c_k1 = find_max(Y)
		count = count + 1
		if(count >= 10000):
			break
	end = time.process_time()
	return X, c_k1, end - start, count

print(' n   eigen_value          time       num_iteration')
for i in [10,20,40]:
	Z = random.randn(i, i)
	eigen_vector, eigen_value, duration, count = power_iteration(Z, 1e-3,i)
	print('%3d'%i + '%13f'%eigen_value + '%14f'%duration + '%10d'%count)

测试结果：

 n   eigen_value          time       num_iteration
 10    -3.099009      0.000850        47
 20     5.005884      0.000702        22
 40    -6.902739      0.002538        74

根据资料，迭代次数和选取的初始向量X0有很大关系，更多详细的幂迭代资料可以戳这里 Spark2.0机器学习系列之11： 聚类(幂迭代聚类， power iteration clustering， PIC)。

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?

def max_singular_value(n, p):
	C = random.binomial(1,p,n*n)
	C = C.reshape((n,n))
	u, s, vh = linalg.svd(C)
	print('%-5s'%n, '%-5s' %p ,'%f' %numpy.max(s))

print('n     p     svd')
for i in range(50,251,50):
	max_singular_value(i,0.5)

for i in range(1,6,1):
	max_singular_value(500,i/10)

程序运行结果如图所示：

n     p     svd
50    0.5   24.729492
100   0.5   49.957149
150   0.5   75.181004
200   0.5   99.652842
250   0.5   125.230619
500   0.1   50.708960
500   0.2   100.526990
500   0.3   151.212283
500   0.4   199.728628
500   0.5   250.616401

我们可以看出 max_svd（最大奇异值） ≈ n * p

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.

def nearest_neighbor(z, A):
    D = A - z
    D = numpy.abs(D)
    index = numpy.argmin(D) #argmin返回值最小的元素的下标
    min_ = A[index // D[0].size][index % D[0].size]
    return min_

n = 5
A = random.rand(n, n)
z = 0.5
print(C)
print("The nearest neighbor is:" + str(nearest_neighbor(z, A)))