hw11（第十一周）

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.

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

Exercise 9.2: Solving a linear system Generate a vector b with m entries and solve Bx = 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.

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.

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?

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.

代码：
ex9_x代表了第x题的代码。其中main函数中生成了题中最前面要求的两个矩阵。ex9_1中实现了求A(B - xI)的函数（以lambda表达式的形式），ex9_4中实现了幂迭代算法。

import numpy
from scipy import linalg
import time
n = 200
m = 500
# Random Gaussian Matrix
A = numpy.random.randn(n, m)
x = numpy.random.randn(m)
y = numpy.random.randn(m)
B = linalg.toeplitz(x, y)
print(A)
print(A.shape)
print(B)
print(B.shape)
# 9.1
def ex9_1(A, B):
    print("9.1")
    # A + A
    AplusA = numpy.add(A, A)
    print(AplusA)
    # AAT
    AAT = numpy.dot(A, A.T)
    print(AAT)
    # ATA
    ATA = numpy.dot(A.T, A)
    print(ATA)
    # AB
    AB = numpy.dot(A, B)
    print(AB)
    # A(B - xI) for any x
    compute = lambda A, B, x: numpy.dot(A, numpy.subtract(B, numpy.asarray(numpy.eye(B.shape[0], B.shape[1])) * x))
    print(compute(A, B, 0))

# 9.2
def ex9_2(B):
    print("9.2")
    b = numpy.random.randn(B.shape[0], 1)
    print(b)
    x = numpy.linalg.solve(B, b)
    print(x)

# 9.3
def ex9_3(A, B):
    print("9.3")
    fnA = numpy.linalg.norm(A, 'fro')
    print(fnA)
    infnB = numpy.linalg.norm(B, numpy.inf)
    print(infnB)

# 9.4    
def ex9_4(n):
    print("9.4")
    def power_iteration(A, epsilon=1e-5):
        # Ideally choose a random vector
        # To decrease the chance that our vector
        # Is orthogonal to the eigenvector
        b_k = numpy.random.rand(A.shape[0])
        b_k = b_k / numpy.linalg.norm(b_k)
        num_sim = 0
        d_k = 0
        while True:
        #for _ in range(num_simulations):
            # calculate the matrix-by-vector product Ab
            b_k1 = numpy.dot(A, b_k)

            # calculate the norm
            b_k1_norm = numpy.linalg.norm(b_k1)

            # re normalize the vector
            b_k2 = b_k1 / b_k1_norm
            d_k1 = numpy.linalg.norm(b_k2 - b_k)
            #print(abs(d_k - d_k1))
            if abs(d_k - d_k1) < epsilon:
                break
            d_k = d_k1
            num_sim += 1
            b_k = b_k2

        return (b_k, num_sim)

    Z = numpy.random.randn(n, n)
    for _ in range(5): # Make sure the eigenvalue is valid.
        print("Turn %d"%_)
        start_time = time.clock()
        (v, num) = power_iteration(Z, 1e-5)
        x = numpy.dot(numpy.dot(Z, v).T, v) / numpy.dot(v.T, v)
        end_time = time.clock()
        print(Z)
        print("eigenvector: ")
        print(v)
        print("eigenvalue: ", end="")
        print(x)
        print("iteration num: %d"%num)
        print("Time used ", end="")
        print(end_time - start_time, end="")
        print("(ms)")

def ex9_5(n, p=0.5):
    print("9.5")
    C = numpy.array([0 if numpy.random.rand() > p else 1 for _ in range(n * n)])
    C.shape = (n, n)
    print(C)
    UsVh = numpy.linalg.svd(C)
    print("Singular values: ")
    print(UsVh[1])
    print("n = %d, p = %f, sv = %f"%(n, p, UsVh[1][0]))
    print("largest sv = n * p ?")

def ex9_6(A, z=0.5):
    print("9.6")
    index = (numpy.abs(A - z)).argmin()
    #A[index] is the closest
    print(A)
    print(z)
    print(A[index])

ex9_1(A, B)
ex9_2(B)
ex9_3(A, B)
ex9_4(n)
'''
ex9_4(5)
ex9_4(10)
ex9_4(200)
ex9_4(2000)
'''
ex9_5(n, p=0.6)
ex9_6(A.reshape(A.shape[0] * A.shape[1]), z=0.5)