numpy练习

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

说明：由于这里的矩阵A，B都是随机生成的，并且都很巨大，因此不方便展示运行结果，这里的练习仅展示代码。

初始化：

import numpy as np  
import scipy.linalg  
import matplotlib.pyplot as plt  
n = 200
m = 500  
A = np.mat(np.random.normal(size = (200, 500)))  
B = np.mat(scipy.linalg.toeplitz([np.random.normal(0, 1) for i in range(m)]))  

说明：以下所有代码都是基于以上初始化

Exercise 9.1: Matrix operations Calculate A + A, AA>,A>A and AB. Write a function that computes A(B−λI) for any λ. 

def func(Lambda):  
    return A * (B - Lambda * np.eye(m))  
C1 = A + A  
C2 = A * A.T  
C3 = A.T * A  
C4 = A * B  
func(1) 

Exercise 9.2: Solving a linear system Generate a vector b with m entries and solve Bx = b.

b = np.random.normal(0, 1, m)
x = np.linalg.solve(B, b) 

Exercise 9.3: Norms Compute the Frobenius norm of A: kAkF and the infinity norm of B: kBk∞. Also find the largest and smallest singular values of B.

FNormOfA = scipy.linalg.norm(A)  
INormOfB = scipy.linalg.norm(B, np.inf)  
SmallestSingularValueOfB = min(scipy.linalg.svdvals(B))  
LargestSingularValueOfB = max(scipy.linalg.svdvals(B))  

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.

Z = np.random.normal(size = (500, 500))  
b = np.random.normal(500)  
step = 0  
while 1:  
    step += 1
    temp = np.dot(Z, b)  
    temp /= scipy.linalg.norm(temp)  
    if abs(scipy.linalg.norm(b, np.inf) - scipy.linalg.norm(temp, np.inf)) < 0.0001: break  
    b = temp  
print(step) 

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?

C = np.empty([n, n])  
for i in range(n):  
    for j in range(n):  
        x = np.random.random()  
        if x < p:  
            C[i, j] = 1  
        else:  
            C[i, j] = 0  
singulars = svdvals(C)  
print(singulars)  
print(np.linalg.norm(singulars, ord=np.inf))  

说明：发现当n比较大时，近似满足最大奇异值=np。

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(A, z):  
    B = z*np.ones([A.shape[0], A.shape[0]])  
    A -= B  
    print(A)  
    print(A.flat[abs(A).argmin()]+z)