Numpy Exercise

# Numpy Exercise
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
import time

# Generate random seed
seed = numpy.int64(time.time())
numpy.random.seed(seed)
print("Random seed:\t", seed)

# Generate matrix A with random Gaussian entries
n = 200
m = 4
A = numpy.random.normal(size=[n, m])
print("A =\n", A)


# Generate Toeplitz matrix B
B = numpy.zeros([m, m])
row = numpy.random.normal(size=[m])
column = numpy.random.normal(size=[m])
for i in range(0, m):
    for j in range(i, m):
        B[j][j-i] = column[i]
    for j in range(0, m - i):
        B[j][j+i] = row[i]
print("B =\n", B)

Exercise 1: Matrix operations

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

# A + A
A_plus_A = A + A
print("A + A = ", A_plus_A)
# AT
AT = A.transpose()
print("AT = ", AT)
# A * AT
A_mul_AT = numpy.dot(A, AT)
print("A * AT = ", A_mul_AT)
# AT * A
AT_mul_A = numpy.dot(AT, A)
print("AT * A = ", AT_mul_A)
# A * B
A_mul_B = numpy.dot(A, B)
print("A * B = ", A_mul_B)


# Calculate A * (B - XI)
def A_MUL_B_Min_XI(A, B, X):
    I = numpy.eye(len(B), len(B))
    return numpy.dot(A, B - X * I)


X = 3
print("A * (B - 3I) = ", A_MUL_B_Min_XI(A, B, 3))

Exercise 2: Solving a linear system

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

# B * x = b -> return x
def solve_linear(B, b):
    return numpy.linalg.solve(B, b)

# b is generated vector and x is the solution
b = numpy.random.normal(size=[m])
print("b = ", b)
x = solve_linear(B, b)
print("x = ", x)

Exercise 3: Norms

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

# Calculate norms
A_forbenius_norm = numpy.linalg.norm(A, "fro")
B_inf_nurm = numpy.linalg.norm(B, numpy.inf)
U, S, UT = numpy.linalg.svd(B)
B_max_singular = numpy.max(S)
B_min_singular = numpy.min(S)
print("A's Forbenius norm:\t", A_forbenius_norm)
print("B's inifinity norm:\t", B_inf_nurm)
print("B's largest singular value:\t", B_max_singular)
print("B's smallest singular value:\t", B_min_singular)

Exercise 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.

# Power iteration
n = 200
Z = numpy.random.normal(size=[n, n])
end = 1e-8
eigenvector = numpy.random.normal(size=[len(Z)]).transpose()
eigenvector = eigenvector / numpy.max(eigenvector)
temp = eigenvector.copy()
deviation = numpy.inf
eigenvalue = numpy.max(temp)
iteration = 0
start_time = time.clock()
while abs(deviation) > end:
    temp = numpy.matmul(Z, eigenvector)
    max_temp = numpy.max(temp)
    deviation = eigenvalue - max_temp
    eigenvalue = max_temp
    eigenvector = temp / max_temp
    iteration += 1
finish_time = time.clock()
print("Total iterations:\t", iteration)
print("Total time:\t", format(finish_time - start_time))
print("Eigenvalue:\t", eigenvalue)
print("Eigenvector:\t", eigenvector)

Exereise 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?

# Singular values
n = 200
p_store = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
for p in p_store:
    print("------------------------------------------------------------")
    print("p = ", p)
    C = numpy.random.random(size=[n, n])
    for i in range(0, n):
        for j in range(0, n):
            if C[i][j] > p:
                C[i][j] = 1
            else:
                C[i][j] = 0
    print("C =\n", C)
    U, S, UT = numpy.linalg.svd(C)
    max_singular = numpy.max(S)
    print("S = ", S)
    print("Largest singular value:\t", max_singular)
    print("n * (1-p):\t", n * (1-p))
print("------------------------------------------------------------")

通过计算发现，最大奇异值max_singular ≈ n * (1 - p)

Exercise6: 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.

# Nearest neighbor
def nearest_neighbor(A, z):
    neighbor = numpy.argmin(numpy.abs(A - z))
    return A[neighbor // numpy.shape(A)[0]][neighbor % numpy.shape(A)[1]]

创建一个这样的矩阵，其对应(i, j)位置的值为A[i][j] - z，即A[i][j]与z的差。找出该矩阵中使该差最小的元素返回即可。